Generalize to non-Hermitian models

Author: thchr

CLAUDE.md

@@ -22,7 +22,7 @@ hopping ranges, the package:
 ## Key types and API
 
 ### Core pipeline
-- `tb_hamiltonian(cbr, Rs; antihermitian)` — top-level entry point; returns `TightBindingModel`
+- `tb_hamiltonian(cbr, Rs, Val{hermiticity})` — top-level entry point; returns `TightBindingModel`
 - `obtain_symmetry_related_hoppings(Rs, br_a, br_b)` — enumerates hopping orbits
 - `sgrep_induced_by_siteir(br, op)` — site-symmetry induced space group representation
 

docs/make.jl

@@ -25,6 +25,7 @@ makedocs(;
         "Band symmetry" => "band-symmetry.md",
         "Berry curvature & Chern numbers" => "berry.md",
         "Symmetry breaking" => "symmetry-breaking.md",
+        "Nonhermitian models" => "nonhermitian.md",
         "API" => "api.md",
         "Internal API" => "internal-api.md",
         "Theory" => "theory.md",

docs/src/nonhermitian.md

@@ -0,0 +1,60 @@
+# [Non-Hermitian tight-binding models](@id nonhermitian)
+
+By default, the models returned by `tb_hamiltonian` are Hermitian. In addition to Hermitian models, however, it is also possible to return anti-Hermitian or generically non-Hermitian (neither Hermitian or anti-Hermitian) models. Here, we showcase the latter.
+
+## Hatano--Nelson model
+
+The architypical non-Hermitian model is the 1D Hatano--Nelson model, consisting of a single site and nearest-neigbor hoppings, in a setting with only time-reversal symmetry.
+It is simple to build this model with SymmetricTightBinding.jl:
+
+```@example hatano-nelson
+brs = calc_bandreps(1, 1)  # EBRs in plane group 1, with time-reversal symmetry
+pin_free!(brs, [1=>[0]]) # the 1a Wyckoff position in plane group 1 has a free parameter: set to 0 for definiteness
+cbr = @composite brs[1]    # single-site model
+tbm = tb_hamiltonian(cbr, [[1]], NONHERMITIAN) # nearest neigbor hoppings
+```
+
+The model is very simple: two different hopping terms, corresponding to right- and left-directed hopping terms. The absence of hermiticity allows the hopping amplitudes in either direction to differ, contrasting the Hermitian case:
+
+```@example hatano-nelson
+tb_hamiltonian(cbr, [[1]], HERMITIAN)
+```
+The non-Hermitian model reduces to the Hermitian model when the left- and right-directed hopping amplitudes are equal. When the two are _not_ equal, the Hatano-Nelson model features exceptional points and spontaneous symmetry breaking of the real spectrum, as we can verify by example (using Brillouin.jl and GLMakie.jl for dispersion plotting):
+
+```@example hatano-nelson¨
+ptbm = tbm([0.9, 1.1]) # a model with 0.9 hopping amplitude to right, 1.1 to the left
+
+using Brillouin, GLMakie
+kp = irrfbz_path(1, directbasis(1, 1)) # a k-path in plane group 1
+kpi = interpolate(kp, 500) # interpolated over 500 points
+Es = spectrum(ptbm, kpi)
+Es_re = sort(real.(Es); dims=2)
+Es_im = sort(imag.(Es); dims=2)
+
+using GLMakie
+plot(kpi, Es_re, Es_im; color=[:blue, :red])
+```
+
+We can also explore the consequences of breaking time-reversal symmetry:
+```
+brs_notr = calc_bandreps(1, 1; timereversal=false)  # EBRs in plane group 1, without time-reversal symmetry
+pin_free!(brs_notr, [1=>[0]])
+tbm_notr = tb_hamiltonian((@composite brs_notr[1]), [[0], [1]], NONHERMITIAN) # on-site terms _and_ nearest-neighbor hoppings
+```
+
+## A more complicated example: exceptional surfaces in p4
+
+We can also construct more complicated examples where symmetry plays a role. Consider for example the following simple extension of the Hatano-Nelson model to a 2D lattice with p4 symmetry:
+```
+brs = calc_bandreps(10, Val(2))
+cbr = @composite brs[1]
+tbm_H  = tb_hamiltonian(cbr, [[0,0], [1,0]], Val(HERMITIAN))
+tbm_NH = tb_hamiltonian(cbr, [[0,0], [1,0]], Val(NONHERMITIAN))
+```
+
+It is instructive to visualize both the Hermitian and non-Hermitian models and compare the involved hopping terms:
+
+```
+plot(tbm_H)
+plot(tbm_NH)
+```

docs/src/symmetry-breaking.md

@@ -47,7 +47,7 @@ This allows three additional terms. Conversely, we could have also tried to brea
 Δtbm_tr = subduced_complement(tbm, 11; timereversal = false) # break TR
 ```
 
-[^1]: For mirror symmetry-breaking, the absence of new terms is a result of looking only at a limited set of hopping orbits (in the original model `tb_hamiltonian(cbr, [[0,0], [1,0]])`): by including longer-range hopping orbits, we would eventually find new mirror-symmetry-broken terms. This is not so for time-reversal breaking, however: in *p*4mm, mirror symmetry and hermicity jointly impose an effective time-reversal symmetry.
+[^1]: For mirror symmetry-breaking, the absence of new terms is a result of looking only at a limited set of hopping orbits (in the original model `tb_hamiltonian(cbr, [[0,0], [1,0]])`): by including longer-range hopping orbits, we would eventually find new mirror-symmetry-broken terms. This is not so for time-reversal breaking, however: in *p*4mm, mirror symmetry and hermiticity jointly impose an effective time-reversal symmetry.
 
 However, by breaking both mirror and time-reversal symmetry simultaneously, additional terms do appear:
 

examples/sg213.jl

@@ -10,19 +10,22 @@ using Makie
 ## --------------------------------------------------------------------------------------- #
 
 const MaybeCoefficient = Union{Nothing, <:AbstractVector{<:Number}}
-
-@recipe(HoppingOrbitPlot, h, Rs, t, offdiag) do Scene
-    Attributes(;
-        origins = Attributes(; color = :firebrick2, label = "Origins (a)"),
-        destinations = Attributes(; color = :royalblue1, label = "Destinations (b+R)"),
-        markersize = 0.05,
-        bonds = Attributes(; color = :gray37, linewidth = 2.0, label = "Bonds"),
-        unitcell = Attributes(; color = :gray65, linewidth = 2.0, label = "Unit cell",
-                                patchcolor = :gray97),
-        context = Attributes(; include = true, color = :gray85, 
-                               linecolor = :gray90, linewidth= 1.25,
-                               limits = nothing),
-    )
+const default_context_attributes = Attributes(;
+    include = true, color = :gray85, linecolor = :gray90, linewidth = 1.25, limits = nothing
+) # TODO: we define this so we can manually merge: remove this hack once Makie v0.25 is out
+
+@recipe HoppingOrbitPlot (h, Rs, t, offdiag) begin
+    origins = Attributes(; color = :firebrick2, label = "Origins (a)")
+    destinations = Attributes(; color = :royalblue1, label = "Destinations (b+R)")
+    markersize = 0.05
+    bonds = Attributes(; color = :gray37, linewidth = 2.0, label = "Bonds")
+    unitcell = Attributes(; color = :gray65, linewidth = 2.0, label = "Unit cell",
+                            patchcolor = :gray97)
+    context = default_context_attributes
+    # TODO: Broken as-is: caller must give _full_ attributes for any field that is itself
+    #       an `Attributes(...)`: no automatic merging of default attributes with partial
+    #       caller attributes.
+    #       Will apparently be fixed in Makie v0.25+ (cf. ffreyer, Slack).
 end
 
 """
@@ -50,7 +53,9 @@ function Makie.plot!(
     h = p.h[]   # TODO: actually do the Observables updates; just so tedious...
     Rs = p.Rs[]
     t = p.t[]
-    offdiag = p.offdiag[] # implies that (anti-)hermicity requires adding reversed hoppings
+    offdiag = p.offdiag[] # implies that (anti-)hermiticity requires adding reversed hoppings
+    context = p.context[]
+    context = merge(context, default_context_attributes) # TODO: remove manual merging once Makie v0.25+ is out
 
     P, V = Point{D, Float32}, Vec{D, Float32}
     Rm = stack(Rs)
@@ -135,7 +140,7 @@ function Makie.plot!(
     )
 
     # set square axis limits, centered around unit cell center
-    bbox = if isnothing(p.context[].limits[])
+    bbox = if isnothing(context.limits[])
         bbox_coords = Makie.GeometryBasics.coordinates(data_limits(p))
         cntr = sum(Rs) ./ 2
         max_dist = maximum(abs,
@@ -145,11 +150,11 @@ function Makie.plot!(
         width = max_dist + pad*0.1
         Rect{D, Float32}(cntr-V(width), V(2*width))
     else
-        p.context[].limits[] :: Rect{D, Float32}
+        context.limits[] :: Rect{D, Float32}
     end
 
     # include symmetry-related sites, that we didn't hop to
-    if p.context[].include[]
+    if context.include[]
         i_lower = Rm \ (bbox.origin)
         i_lower = ntuple(d->floor(Int, i_lower[d])-1, Val(D))
         i_upper = Rm \ (bbox.origin .+ bbox.widths)
@@ -169,14 +174,14 @@ function Makie.plot!(
             D == 3 && continue # doesn't look nice for 3D
             unitcell′ = unitcell .+ Ref(R)
             any(in(bbox), unitcell′) || continue
-            lines!(unitcell′; color=p.context[].linecolor, linewidth=0.5, depth_shift=0)
+            lines!(unitcell′; color=context.linecolor, linewidth=0.5, depth_shift=0)
         end
         scatter!(
             p,
             rs,
             marker = :circle,
             markersize = 11,
-            color = p.context[].color,
+            color = context.color,
         )
     end
     limits!(bbox)
@@ -321,7 +326,10 @@ _orbit(h::HoppingOrbit) = h
 _coefficients(tbt::TightBindingTerm) = _coefficients(tbt.block)
 _coefficients(tbb::TightBindingBlock) = tbb.t
 _offdiag(tbt::TightBindingTerm) = _offdiag(tbt.block)
-_offdiag(tbb::TightBindingBlock) = !tbb.diagonal_block
+function _offdiag(tbb::TightBindingBlock{D, S}) where {D, S}
+    S === NONHERMITIAN && return false
+    return !tbb.diagonal_block
+end
 
 ## --------------------------------------------------------------------------------------- #
 # Plotting entire tight-binding models by tiling their hopping terms
@@ -347,8 +355,12 @@ function Makie.convert_arguments(
     for idx in LinearIndices(axs)
         if idx ≤ length(tbm)
             tbt = tbm[idx]
-            plots = [S.HoppingOrbitPlot(_orbit(tbt), Rs, _coefficients(tbt), _offdiag(tbt);
-                                        context = (; limits = bbox))]
+            plots = [
+                S.HoppingOrbitPlot(
+                    _orbit(tbt), Rs, _coefficients(tbt), _offdiag(tbt);
+                    context = (; limits = bbox)
+                )
+            ]
         else
             plots = PlotSpec[]
         end
@@ -417,7 +429,8 @@ function layout_limits(hs::AbstractVector{HoppingOrbit{D}}, Rs::DirectBasis{D})
 
     # add points from unit cell
     rect = Rect{D, Float32}(P(0), V(1)) # unit cube at origin
-    sites = P.(Ref(Rm) .* Makie.GeometryBasics.coordinates(rect))
+    rect_coords = D == 1 ? [P(0), P(1)] : Makie.GeometryBasics.coordinates(rect)
+    sites = P.(Ref(Rm) .* rect_coords)
     filter!(!isnan, sites)
     pts = @view sites[1:end] # for later referencing only the unit cell
 
@@ -437,13 +450,14 @@ function layout_limits(hs::AbstractVector{HoppingOrbit{D}}, Rs::DirectBasis{D})
     lower_bound = ntuple(d -> minimum(v -> getindex(v, d), sites), Val(D))
     upper_bound = ntuple(d -> maximum(v -> getindex(v, d), sites), Val(D))
     data_bbox = Rect{D, Float32}(P(lower_bound), V(upper_bound .- lower_bound))
-    bbox_coords = Makie.GeometryBasics.coordinates(data_bbox)
+    bbox_coords = D == 1 ? [P(lower_bound), P(upper_bound)] : Makie.GeometryBasics.coordinates(data_bbox)
     cntr = sum(Rs) ./ 2
     max_dist = maximum(abs,
         ntuple(d->maximum(v->abs(cntr[d]-getindex(v, d)), bbox_coords), Val(D)))
     pad = maximum(abs, 
         ntuple(d -> splat(-)(extrema(v -> getindex(v, d), filter(!isnan, pts))), Val(D)))
     width = max_dist + pad*0.1
+    # TODO: the padding and added width should surely not be the same in all dimensions
     return Rect{D, Float32}(cntr-V(width), V(2*width))
 end
 function layout_limits(tbm::TightBindingModel{D}, Rs::DirectBasis{D}) where D

ext/SymmetricTightBindingMakieExt.jl

@@ -10,9 +10,10 @@ import SymmetricTightBinding: fit
 # Define loss as sum of absolute squared error (MSE, up to scaling)
 
 function fg!(
-    F, G, cs, tbm::TightBindingModel, Em_r, ks;
+    F, G, cs, tbm::TightBindingModel{D, S}, Em_r, ks;
     lasso::Union{Nothing,Real} = nothing
-)
+) where {D, S}
+    S ≠ HERMITIAN && error("loss function can currently only handle HERMITIAN models")
     ptbm = tbm(cs)
     if !isnothing(G)
         fill!(G, zero(eltype(G)))

ext/SymmetricTightBindingOptimExt.jl

@@ -23,6 +23,8 @@ const ZASSENHAUS_ATOL_DEFAULT = NULLSPACE_ATOL_DEFAULT
 include("types.jl")
 export HoppingOrbit
 export TightBindingBlock
+export Hermiticity, HERMITIAN, ANTIHERMITIAN, NONHERMITIAN
+export hermiticity
 export TightBindingModel
 export ParameterizedTightBindingModel
 include("show.jl")

src/SymmetricTightBinding.jl

@@ -48,13 +48,19 @@ The gradient is returned as column vectors, one for each band, with each column
 the gradient of the corresponding energy with respect to the hopping coefficients of `ptbm`.
 """
 function energy_gradient_wrt_hopping(
-    ptbm::ParameterizedTightBindingModel{D},
+    ptbm::ParameterizedTightBindingModel{D, S},
     k::ReciprocalPointLike{D},
     (Es, us) = solve(ptbm, k; bloch_phase=Val(false)) # "unperturbed" energies & eigenstates
     ;
     degen_rtol::Float64 = 1e-12,
     degen_atol::Float64 = 1e-12
-) where D
+) where {D, S}
+    if S === NONHERMITIAN
+        # TODO: requires left/right version of Feynman-Hellmann theorem + possibly handling
+        #       of defective case
+        error("energy gradient with respect to hopping is not currently implemented for \
+               NONHERMITIAN models")
+    end
     Nᶜ = length(ptbm.tbm) # number of hopping terms
     Nᵇ = ptbm.tbm.N       # number of bands
 
@@ -82,19 +88,21 @@ function energy_gradient_wrt_hopping(
     end
 
     # apply Feynman-Hellmann theorem, either in degenerate or non-degenerate variants
-    ∇ᶜEs = Matrix{Float64}(undef, Nᶜ, Nᵇ)
+    ResultType = S == HERMITIAN ? Float64 : ComplexF64
+    ∇ᶜEs = Matrix{ResultType}(undef, Nᶜ, Nᵇ)
     for i in 1:Nᶜ
         ∂ᵢH = tbmg(k, i)
         for ns in bands
             if length(ns) == 1 # non-degenerate band
                 n = @inbounds ns[1]
                 uₙ = @view us[:, n]
                 ∂ᵢEₙ = dot(uₙ, ∂ᵢH, uₙ)
-                ∇ᶜEs[i, n] = real(∂ᵢEₙ)
+                ∇ᶜEs[i, n] = S == HERMITIAN ? real(∂ᵢEₙ) : ∂ᵢEₙ
             else               # degenerate bands
                 us′ = @view us[:, ns]
-                M = us′' * ∂ᵢH * us′ # Mₙₘ = ⟨uₙ|∂ᵢH|uₘ⟩ for n,m ∈ `ns`
-                ∂ᵢEs′ = eigvals!(Hermitian(M)) # ∂ᵢEₙ for n in `ns`
+                _M = us′' * ∂ᵢH * us′ # Mₙₘ = ⟨uₙ|∂ᵢH|uₘ⟩ for n,m ∈ `ns`
+                M = S == HERMITIAN ? Hermitian(_M) : _M
+                ∂ᵢEs′ = eigvals!(M) # ∂ᵢEₙ for n in `ns`
                 ∇ᶜEs[i, ns] = ∂ᵢEs′
             end
         end
@@ -215,12 +223,12 @@ The function is analogous to `evaluate_tight_binding_term!`, but computes moment
 gradient components rather than the Hamiltonian matrix itself.
 """
 function evaluate_tight_binding_momentum_gradient_term!(
-    tbt::TightBindingTerm{D},
+    tbt::TightBindingTerm{D, S},
     k::ReciprocalPointLike{D},
     components::NTuple{C, Int},
     c::Union{Nothing, <:Number} = nothing,
     ∇Hs::NTuple{C, Matrix{ComplexF64}} = ntuple(_ -> zeros(ComplexF64, size(tbt)), Val(C)),
-) where {D, C}
+) where {D, S, C}
     block = tbt.block
     block_i, block_j = tbt.block_ij
     is = tbt.axis[Block(block_i)] # global row indices
@@ -246,8 +254,9 @@ function evaluate_tight_binding_momentum_gradient_term!(
                 ∇Hᵢⱼ *= 2im * π # (-2πi) factor from ∂/∂kᵢ of e^{-2πik·δ}
                 isnothing(c) || (∇Hᵢⱼ *= c) # multiply by coefficient if provided
                 ∇H[i, j] += ∇Hᵢⱼ
-                i == j && continue # don't add diagonal elements twice
-                ∇H[j, i] += tbt.hermiticity == ANTIHERMITIAN ? -conj(∇Hᵢⱼ) : conj(∇Hᵢⱼ)
+                if S !== NONHERMITIAN && i ≠ j # off-diagonal contribution (& don't double-add diagonal)
+                    ∇H[j, i] += S == ANTIHERMITIAN ? -conj(∇Hᵢⱼ) : conj(∇Hᵢⱼ)
+                end
             end
         end
     end

src/gradients.jl

@@ -4,10 +4,10 @@
         h_orbit::HoppingOrbit{D},
         brₐ::NewBandRep{D},
         brᵦ::NewBandRep{D},
+        antihermitian::Bool,
         orderingₐ::OrbitalOrdering{D} = OrbitalOrdering(brₐ),
         orderingᵦ::OrbitalOrdering{D} = OrbitalOrdering(brᵦ),
-        Mm::AbstractArray{Int, 4} = construct_M_matrix(h_orbit, brₐ, brᵦ, orderingₐ, orderingᵦ);
-        antihermitian::Bool = false,
+        Mm::AbstractArray{Int, 4} = construct_M_matrix(h_orbit, brₐ, brᵦ, orderingₐ, orderingᵦ)
     ) where {D}
 
 Constructs a basis for the coefficient vectors `t⁽ⁿ⁾` that span the space of Hermitian (or
@@ -27,10 +27,10 @@ function obtain_basis_free_parameters_hermiticity(
     h_orbit::HoppingOrbit{D},
     brₐ::NewBandRep{D},
     brᵦ::NewBandRep{D},
+    antihermitian::Bool,
     orderingₐ::OrbitalOrdering{D} = OrbitalOrdering(brₐ),
     orderingᵦ::OrbitalOrdering{D} = OrbitalOrdering(brᵦ),
-    Mm::AbstractArray{Int, 4} = construct_M_matrix(h_orbit, brₐ, brᵦ, orderingₐ, orderingᵦ);
-    antihermitian::Bool = false,
+    Mm::AbstractArray{Int, 4} = construct_M_matrix(h_orbit, brₐ, brᵦ, orderingₐ, orderingᵦ)
 ) where {D}
     # Determine a basis for coefficient vectors t⁽ⁿ⁾ that span the space of Hermitian (or
     # `antihermitian` if `true`) TB Hamiltonians Hₛₜ(k) = vᵢ(k) Mᵢⱼₛₜ tⱼ = vᵀ(k) M⁽ˢᵗ⁾ t. 
@@ -65,8 +65,8 @@ function constraint_hermiticity(
 ) where {D}
     Q = Array{Int}(undef, size(Mm))
     opI = inversion(Val(D)) # inversion operation
-    Pᵀ =
-        transpose(_permute_symmetry_related_hoppings_under_symmetry_operation(h_orbit, opI))
+    P = _permute_symmetry_related_hoppings_under_symmetry_operation(h_orbit, opI)
+    Pᵀ = transpose(P)
     for s in axes(Mm, 3)
         for t in axes(Mm, 4)
             Q[:, :, s, t] .= Pᵀ * Mm[:, :, t, s] # Pᵀ M⁽ᵗˢ⁾