Spectrum not working for OpenHamiltonian

[seadra]

Hi,

I have a composite system made of leads attached to a scattering region, where leads is defined as a self energy term. After attaching the leads, the result is an OpenHamiltonian object. Calling spectrum on it gives an error, which says spectrum doesn't do OpenHamiltonian objects.

Is this planned for future? For a simple scattering region in 1D, I can obtain eigenenergies analytically for this system, but of course, it doesn't work for more complicated situations.

[seadra]

Maybe somewhat relevant: tkwant seems to has an API for bound states https://tkwant.kwant-project.org/doc/dev/tutorial/boundstates.html

[pablosanjose]

Hi @seadra, thanks for posting! Interestingly, a similar question was recently posted in discourse. As you see from the answer there the question of finding the "spectrum" in an open system is somewhat ill-defined. Finding the spectral density is well defined, though. You can do that with ldos, instead of spectrum

If what you are really after is the discrete part of the spectrum (which may or may not exist for any given problem), ldos can still help, but you will need to broaden states using complex ω's so you can see them as a density peak of finite width. Using the tkwant example you link to, we could do it like this:

julia> using GLMakie, Quantica

julia> h_inf = LP.linear() |> onsite(0) - hopping(1);

julia> g_lead = h_inf |> greenfunction(GS.Schur(boundary = 0));

julia> g_open = h_inf |> supercell |> @onsite!(o -> o - 2.5) |> attach(g_lead) |> attach(g_lead, reverse = true) |> greenfunction;

julia> ρ = ldos(g_open[]);

julia> ωs = subdiv(-4,2,200); ρs = [sum(ρ(ω+0.0001im)) for ω in ωs];

julia> lines(ωs, ρs)

As you see there is a peak below the continuum that corresponds to a bound state.

An alternative way is the one suggested in the tkwant documentation. One can do a finite system by incorporating a large portion of the leads into the central region, which can be seen as an approximation to the open Hamiltonian. Since this is however finite, its spectrum can be computed. It exhibits a bound state at the energy shown by ldos.

julia> h_finite = h_inf |> supercell(region = r -> -40 <= r[1] <= 40) |> @onsite!(o -> o - 2.5, region = r -> iszero(r));

julia> es, _ = spectrum(h_finite);

julia> qplot(bands(h_inf))

julia> foreach(real.(es)) do e
           lines!([-1,1], [e, e], color = :blue)
       end

[pablosanjose]

There is also the interesting option of using GreenSolvers.KPM over h_finite, which can efficiently handle really large systems. This allows to compute the spectral density, including the discrete spectrum, without any artificial broadening (the broadening is optimally produced by the order of the Chebychev expansion used by the solver)

julia> h_finite = h_inf |> supercell(region = r -> -1000 <= r[1] <= 1000) |> @onsite!(o -> o - 2.5, region = iszero);

julia> g_finite = h_finite() |> attach(nothing, region = iszero) |> greenfunction(GS.KPM(order = 1000));
┌ Warning: Computing spectrum bounds... Consider using the `bandrange` option for faster performance.
└ @ Quantica ~/.julia/dev/Quantica/src/solvers/green/kpm.jl:167
┌ Warning: Computed real bandrange = (-3.201562118716358, 1.9999766493443067)
└ @ Quantica ~/.julia/dev/Quantica/src/solvers/green/kpm.jl:170

julia> ρ = ldos(g_finite[1]);

julia> ωs = subdiv(-4,2,300); ρs = [sum(ρ(ω)) for ω in ωs];

julia> lines(ωs, ρs)

julia> ylims!(-0.001, 0.1)

[seadra]

Thank you so much, this is quite amazing!

That makes sense, I was only interested in the discrete bound state spectrum, but indeed, for infinite leads, there is a continuous spectrum as well which complicates things.

The speed of KPM solver is quite impressive, and since it gives the exact bound state eigenenergies, I think I'll go with a finite system.

None what you wrote is this trivial, and I imagine similar questions will pop up every now and then (like the Discord one). I think it would be very helpful to have a section on bound states in the documentation.

Lastly, thank you so much for making Quantica.jl available for everyone!

[pablosanjose]

Closing now, with a mental note of possibly adding some example of bounds states to the docs in the future, if we can make it useful and somewhat general