Design: GreenFunctions and EffectiveHamiltonians

[pablosanjose]

Closes #1, #2

This is a concrete design of scattering systems and Green functions that combines ideas from #1 and #2, and then some. There are still some implementation details to clarify that could affect some parts of the proposal.

Problems

The basic problems to solve are:

• GreenFunction: Compute Green's function (dense) matrices, or derived quantities such as spectral functions, within a given unit cell of any finite or unbounded AbstractHamiltonians.
  • Allow the above also for semi-infinite systems (in any dimension) by the method of images. A particular case is a semi-infinite quasi-1D lead in a scattering problem
  • Allow to compute Green function matrices between arbitrary points, even in different cells. In periodic systems (even semibounded), this should take into account periodicity to compute them efficiently
  • Allow different algorithms to be used interchangingly. For example, in 1D one could use Schur decompositions or bandstructure integration.
• EffectiveHamiltonian: Allow to compute effective, frequency dependent (sparse) AbstractHamiltonians that incorporate the effect of coupling to other infinite/finite system through open boundaries (i.e. through self-energies or ). Case on point: effective Hamiltonian for the first cell of a semi-infinite lead.
  • Allow for the coupling of several effective Hamiltonians in arbitrary configurations, and produce an effective Hamiltonian out of them. Case on point: effective Hamiltonian of a central region coupled to several leads (independent or coupled between them)
• Compute the GreenFunction between arbitrary sites in any EffectiveHamiltonian
  • Compute conductance between a subset of sites in an EffectiveHamiltonian
  • Compute the corresponding scattering matrices
  • Compute the corresponding scattering states
  • Compute local density of states...

API

As suggested by the list of problems to solve, the whole API revolves around two core types: GreenFunction and EffectiveHamiltonian <: AbstractHamiltonian. Given an h::AbstractHamiltonian{B} (possibly a ParametricHamiltonian) we build each of these over a set of sites and cells, and in the presence of certain boundaries, with

g = green(h, solver; boundaries, cells, siteselectors...)         -> GreenFunction
heff = effective(h, solver; boundaries, cells, siteselectors...)  -> EffectiveHamiltonian

solver is a GreenSolvers.AbstractSolver, one of several possible, each with their own kwargs.

gω = g(ω; params...) is then the Green function matrix gω::Matrix{B} at frequency ω.

hω = heff(ω; params...) is a 0-dimensional hω::Hamiltonian, and hω[] its SparseMatrixCSC{B}.

Note: EffectiveHamiltonian is essentially a ParametricHamiltonian under the hood (with an extra compulsory ω parameter). It must therefore remain possible to parametrically modify one with

parametric(heff, modifiers...)  -> EffectiveHamiltonian

Open and external sites

Both g and h have an underlying 0-dimensional Lattice on which the matrices are defined. We do not consider (in this proposal) any periodic GreenFunctions or EffectiveHamiltonians. In this lattice there are two subsets of special sites: open sites (coupled to the environment) and, optionally, external sites (which represent the environment). The latter can be absent, depending on what the solver provides, in which case the environment has been fully integrated out, inducing a self-energy on open sites.

Only EffectiveHamiltonian's can have external sites; they are always integrated out in GreenFunctions. If hω has no external sites, then the following holds

gω = inv(ω-hω[])

Otherwise, the matrix hω[] will be larger than gω.

The purpose of keeping external sites as part of the EffectiveHamiltonian is that it can be numerically more stable to compute (ω-hω[]) \ sources with them included in hω than to first integrate them out and then do the ldiv.

The GreenFunction of a given EffectiveHamiltonian is obtained with

g = green(heff)

green(heff)(ω; params...) is guaranteed to be equal to gω within numerical precision.

Open sites are obtained with

opensites(heff) == opensites(g) -> Vector{Int}

The self-energy of the environment on said sites is computed with

selfenergy(g, ω; params...) == selfenergy(heff, ω; params...) -> Matrix{B}

Note that this matrix will be square and of size equal to the number of open sites (external sites are again integrated out).

Solvers

Solvers for green and effective are subtypes of AbstractGreenSolver, and are all bundled in the submodule GreenSolvers, aliased as GS. Two solvers are currently envisioned

s = GS.Schur(; deflate = true, inversefree = true)   # for 1D `AbstractHamiltonians`
s = GS.Bands(mesh; bands_kwargs...)                  # For arbitrary `Hamiltonian`'s

To make use of them, these need to be applied first to an h::AbstractHamiltonian with

solver = s(h)

which yields an AbstractAppliedGreenSolver.

The internal API of an AbstractAppliedGreenSolver is the following

green(solver, ω; params...)       -> Matrix{B}
effective(solver, ω; params...)   -> SparseMatrixCSC{B}
ldos(solver, ω; params...)        -> Vector{diag{B}}
selfenergy(solver, ω; params...)  -> Matrix{B}

The functionality of other methods in the external API rely on these internal ones.

Integrating out sites

Apart from the option to select sites of interest upon constructing a GreenFunction or an EffectiveHamiltonian as above, we also have the option to integrate out a subset of sites after construction. This is done with

g´ = green(g; siteselectors...)
heff´ = effective(heff; siteselectors...)

These work by bundling the internal solver of g or heff into a special type of solver that selects a view of the original solver,

GS.GreenView(solver, indices) <: AbstractAppliedGreenSolver

Note that GreenView also needs to know which sites are open after projecting out, namely any site that was originally open or that was coupled to a site that is not included in indices.

Coupling systems

We need a general way to couple a number of EffectiveHamiltonians into another EffectiveHamiltonian, e.g. to be able to define arbitrary ScatteringSystem for transport (where the leads are all instances of EffectiveHamiltonians), or to compute local density of states composite unbounded systems. To this end we need to generalize the couple functionality in master.

The coupling API will now rely on a Coupler object with the following API

c = coupler(hs::AbstractHamiltonian...)
c[1,2] = c[2,1] = hopping_model

# to finalize c:
hamiltonian(c)  -> Hamiltonian
parametric(c)   -> ParametricHamiltonian
effective(c)    -> EffectiveHamiltonian

The setindex! method applies hopping_model to a given pair of hs, populating the corresponding blocks in the final Hamiltonian harmonics upon finalizing the coupler.

Note that one can succintly apply an all-to-all coupling model using broadcasting

c = coupler(hs::AbstractHamiltonian...) .= hopping_model

Note: Any coupled pair should be "compatible", which essentially means both h should have the same Bravais vectors (up to ±1 factors).

Scattering systems

A Scattering system is little more than an EffectiveHamiltonian, where we group different open sites into a number of contacts. We build one simply with

scat = scattering(heff, ss::SiteSelector...)  -> ScatteringSystem

where the ss select the different contact sites among open sites.

A ScatteringSystem supports additional methods to compute the conductance and scattering matrix between contacts i and j

conductance(scat, i, j, ω; Nambu = false, params...)  ->  Real
smatrix(scat, [i, j,] ω; params...) -> Matrix{<:Matrix}

Implementation

A possible sketch of core types:

struct EffectiveHamiltonian{T,E,B,S<:AbstractAppliedGreenSolver} <: AbstractHamiltonian{T,E,0,B}
    h::ParametricHamiltonian{T,E,0,B}
    solver::S
    open::Vector{Int}
    external::Vector{Int}
end

struct GreenFunction{T,E,B,S<:AbstractAppliedGreenSolver}
    solver::S
    lattice::Lattice{T,E,0}
    orbstruct::OrbitalStructure{B}
    open::Vector{Int}
end

struct ScatteringSystem{T,E,B,S<:AbstractAppliedGreenSolver}
    h::EffectiveHamiltonian{T,E,B,S}
    contacts::Vector{Vector{Int}}
end

struct Coupler{T,E,L,B,N<:NTuple{<:Any,AbstractHamiltonian{T,E,L,B}}}
    hs::N
    lattice::Lattice{T,E,L}
    couplings::Vector{Hamiltonian{T,E,L,B}}
    pairs::Vector{Pair{Int,Int}}
    blockoffsets::Vector{Int}  # not sure this is needed?
end

The Coupler type can perhaps be restructured to allow for more efficient coupling insertions. For example, instead of building a whole coupling Hamiltonian when applying a given hopping_model, one could start by decomposing and merging all the non-zero elements of the different hs harmonics (with findnz), and then appending to that list the ones from each of the couplings when finalizing.

[pablosanjose]

Variation 1

Given the fact that ScatteringSystem is essentially the same as EffectiveHamiltonian, it would be possible to drop the former by storing open sites in the latter in groups from the start.

struct EffectiveHamiltonian{T,E,B,S<:AbstractAppliedGreenSolver} <: AbstractHamiltonian{T,E,0,B}
    h::ParametricHamiltonian{T,E,0,B}
    solver::S
    open::Vector{Vector{Int}}
    external::Vector{Int}
end

Whenever we couple two effective Hamiltonians we could just keep their open sites (and external sites) in different groups. The we would just do

conductance(heff, i, j, ω; params...)
smatrix(heff, i, j, ω; params...)

We could also introduce some way to redefine contacts by regrouping open sites into different groups. I can think of one use case. If you integrate out a collection of sites in an heff with heff with effective(heff; siteselectors...), the result will loose its open site grouping. Then you would need to do something like

contacts!(heff, ss::Selectors...)

to regroup them. This should check that no open site is repeated or left out.

[pablosanjose]

An important observation that came up while implementing Coupler. When coupling two AbstractHamiltonians, sublattices with the same names are not merged, but renamed and kept separate. Doing the opposite introduces a surprising degree of complexity in two fronts: (1) the resulting Hamiltonian matrices do not have a block structure, so coupling blocks with c[1,2] = model implies keeping track of disjoint index ranges in the matrix in general; and (2) when updating the modifier pointers to the coupled Hamiltonian sparse matrices, the code becomes gnarly.

For this reason, we defer to a future PR the functionality of merging sublattices, including the corresponding update to modifier pointers in a ParametricHamiltonian. We reserve the name Base.merge(lattices...) for that, and adopt combine(lattices...) for the present functionality.

[pablosanjose]

EDIT: Removed due to revision

[pablosanjose]

Variation 2

(Now implemented in the greens branch, to be soon merged into greatrefactor)

• No coupler approach. Instead we recover tree-like connection of subsystems without the possibility of loops via a syntax h |> attach(args...; site_selection...), where args can be different things
• The result of attach-ing something to an AbstractHamiltonian creates an OpenHamiltonian which wraps one or more SelfEnergy. One can attach several systems to the OpenHamiltonian, each time adding a SelfEnergy to it
• A SelfEnergy wraps a latslice::LatticeSlice (portion of h's lattice selected by site_selection, can span one or more unit cells) and a solver::SelfEnergySolver, which produces an Σ::AbstractMatrix over the orbitals in latslice when we do call!(solver, ω; params...).
• Advanced: call!(solver, ω; params...) can also produce three AbstractMatrix, Vᵣₑ, gₑₑ⁻¹, Vₑᵣ when solver is an ExtendedSelfEnergySolver instead of a RegularSelfEnergySolver, so that Σ = Vᵣₑ gₑₑ Vₑᵣ.
• We build a GreenFunction with g = greenfunction(o::OpenHamiltonian, GS.SomeGreenSolver(options...)), or g = o |> greenfunction(GS.SomeGreenSolver(options...))
• This call performs any work that can be done towards computing the green function at some points and energy of o, without knowing the point or the energy
• One can fix the energy by doing gω = g(ω; params...), which produces a gω::GreenSolution, that wraps, in particular, a GreenSlicer that will produce a flat matrix when we pass the points.
• gω[(; site_selection...), (; site_selection´...)] or gω[cellsites(cell, sites), cellsites(cell´, sites´)] or gω[contact(i), contact(j)] produces the final matrix at energy ω and the selected points. contact(i) denotes all sites selected when we attached the i-th self-energy to o

The above machinery relies on a number of new types, such as Contacts or ContactBlockStructure that are invisible to the user.

The possible syntaxes for attach will be expanded as we go. For the moment we have attach(h, model; ...) to define a selfenergy through a model (like we define a Hamiltonian), and attach(h, glead; ...) to attach a 1D lead with Green function glead. We will also add attach(h, glead, model;...) to use a specific coupling between lead and h, and in general attach(h, g;...) for an arbitrary g::GreenFunction. GreenFunction solvers live in the GreenSolvers === GS submodule. For the moment we have GS.Schur for 1D systems and GS.SparseLU for 0D systems.