Hi Anton, Thank you for the quick reply. I hope it is OK that I ask two related follow-up questions. I have updated my code (see .txt attachment), to include 4 orbitals, the conservation law, and express my matrices via Kronecker products (mentioned previously). I'm still trying to reproduce Fig. 4 in Phys. Rev. Lett. 97, 067007 (2006). Firstly, I'm confused about how to write down the expression for the conductance N - R_{e,\uparrow} - R_{e,\downarrow} + R_{h,\uparrow} + R_{h,\downarrow} in KWANT. Here \uparrow and \downarrow refers to spin up and spin down, e and h refers to electron and hole. As far as I understand smatrix.transmission((i, a), (j, b)) computes the transmission probability from block b of lead j to block a of lead i. In my conservation law the eigenvalues are -1, -1, 1, and 1. Therefore I would expect block 0, 1, 2, and 3 to refer to electron spin up, electron spin down, hole spin up, and hole spin down respectively. However, I suspect this is not true because I obtain an out of bounds error. In the attached code I have naively used the code provided in tutorial 2.6, but I find that the conductance vanishes (which it should not do). My question is loosely related to a question asked in https://mail.python.org/archives/list/kwant-discuss@python.org/message/NLPII... In short, how can I express the formula for the conductance using the notation in KWANT in order to obtain the correct conductance? Secondly, I'm wondering how to express the particle hole symmetry. In my problem (see .PDF attachment) the particle-hole symmetry reads P.H^*.P = -H. In the simplest case where H = H^* this is not a problem. However, eventually I want to generalize my code to a case where H != H^*, but H is still hermitian H = H^{\dagger}. As the complex conjugate operator does not have a matrix representation (that I'm aware of) is there a nice way to implement this symmetry in KWANT? To be more precise let me give my interpretation. If I write lead0 = kwant.Builder(... particle_hole=np.kron(tau_x,s_0)) I think that KWANT interpretes this as P.H.P^-1 = -H, where P = np.kron(tau_x,s_0). In this case the necessary complex conjugation is ignored. Is this a correct interpretation or does KWANT include the complex conjugation? Yours sincerely, Martin F Jakobsen ________________________________ Fra: Anton Akhmerov Sendt: onsdag 19. august 2020 14:07:40 Til: Martin Fonnum Jakobsen Kopi: kwant-discuss@python.org Emne: Re: [Kwant] Andreev reflection in Graphene Hi Martin, Your interpretation of the conservation law is correct, and so does your implementation of it. However I see that you create the graphene lattice as having 2 orbitals, whereas if you have spin and particle-hole, you should have 4 orbitals instead. That is also what the error message says. I usually define all matrix-valued functions via np.kron, so that I remember which Pauli matrix refers to which degree of freedom. Best, Anton On Wed, 19 Aug 2020 at 14:03, Martin Fonnum Jakobsen wrote:

Hello,

I have recently started using KWANT. I hope this is the correct way to ask questions.

I’m interested in reproducing Fig. 4 in PRL 97, 067007 (2006) – Specular Andreev Reflection in Graphene. This has been done in KWANT in the PhD thesis of Tibor Sekera – “Quantum transport of fermions in honeycomb lattices and cold atomic systems”, chapter 4. The problem is similar to tutorial 2.6.

In my problem, I have a tight-binding Hamiltonian living on a hexagonal lattice, from which I can construct a discrete BdG equation (See attachment). In addition, to electron-hole degrees of freedom I’m also interested in including spin because later I will add magnetic fields and such.

I consider a normal metal lead attached to a superconducing scattering region with zero barrier at the interface. Since I have both electron-hole and spin degrees of freedom my BdG equation is 4x4. To account for electron-hole degrees of freedom I have specified that at each site there should be two orbitals. In the tutorial 2.6 it is emphasized that it is important to specify the correct conservation law in the N lead when we want to calculate the conductance. My question is how should the conservation law (and particle hole symmetry) in the normal-metal lead look like for my case?

I have tried the standard expression L_0 = kwant.Builder(sym,conservation_law=s_z) where s_z is the Pauli matrix [[1,0],[0,-1]]. In this case I obtain a dimension mismatch error.

I have also tried a generalized 4x4 matrix version of this

L_0 = kwant.Builder(sym,conservation_law=tinyarray.array([[1,0,0,0],[0,-1,0,0],[0,0,1,0],[0,0,0,-1]]))

Or

L_0 = kwant.Builder(sym,conservation_law=tinyarray.array([[1,0,0,0],[0,1,0,0],[0,0,-1,0],[0,0,0,-1]])).

In these cases I obtain the error “Single `onsite` matrix of shape (4, 4) provided but there are 2 orbitals per site in the system”.

I think my issue is that I do not fully understand what the Conservation law actually does. Initially I thought the conservation law was a matrix C such that C.H.C^-1 should return a block diagonal matrix, where the eigenvalues of C refers to each block. Is this not an accurate interpretation?

I have attached a PDF with the BdG equation, and my python code as a .py and .txt file. The content of the .py and .txt file are identical, which format do you prefer in the future?

I hope this is OK, and not too overwhelming.

I also tried to subscribe (with this e-mail) to the mailing list, but I’m not sure that it worked.

Yours sincerely,

Martin F Jakobsen

Dear Sergey, In terms of scaling, the most convenient is the KPM method, since it scales linearly with the inverse energy resolution (that is, the number of moments), and linearly with the system size. Could you clarify what is the issue with the approach 3)? I am not sure what is the energy resolution that you want to achieve. You can conveniently set this value when creating a `kwant.kpm.SpectralDensity` instance via the arguments `energy_resolution` or (exclusive) `num_moments`. The latter is related to the energy resoution by `\delta \pi / N`, where `\delta` is the full width of the spectrum, and `N` the size of your system. In the case of a 2D, that will be similar to `L^2`, where `L` is the linear size. Furthermore, you could progressively increase the energy resolution by adding moments: ``` spectrum = kwant.kpm.SpectralDensity(...) spectrum.add_moments(...) ``` If you want to resolve individual energy levels then the approach 1) and 3) will have the same scaling, however, the KPM method will not give you eigenvalues. Finally, since you want to get the LDOS at a single site or a small set of sites near the STS tip, you should use the `kwant.kpm.LocalVectors` vector factory, setting the `where` argument to the desired spot. ``` def center(site): return site.tag == (0, 0) vector_factory = kwant.kpm.LocalVectors(syst, where=center) spectrum = kwant.kpm.SpectralDensity(syst, num_moments=100, num_vectors=None, vector_factory=vector_factory) ``` Note that `num_vectors=None` will let the `kwant.kpm` module to figure out how many vectors does the vector factory produce. Otherwise, it should be at most the *total number of orbitals* defined by the `where` function, that is, sites times orbitals per site. As you can see, the `kwant.kpm` module is the most suited for this problem, specially when the sample is very large. I hope this helps. Best regards, Pablo

Thank you, Pablo and Joe,   Yes, KPM method did help, but I want to mention that energy_resolution parameter can be a bit misleading and for sharp DOS features one really needs to set it to notably lower value than naively expected. Here is a bit of my code: def compute_ldos4(sys, phi, depth, width, gap, epoints): print("computing LDOS for fi=" +str(phi)) siteA = sys.id_by_site[a(0,0)] siteB = sys.id_by_site[b(0,0)] where =lambda s: (s.tag == numpy.array([0,0]))and ((s.family == a)or (s.family == b)) vector_factory = kwant.kpm.LocalVectors(sys, where) local_dos = kwant.kpm.SpectralDensity(sys,params={'phi': phi,'depth': depth,'width': width,'gap': gap}, num_vectors=None, vector_factory=vector_factory, mean=False,num_moments=10000) densities = local_dos(epoints) return densities Best wishes, Sergey On 20/08/2020 16:58, Joseph Weston (Aquent LLC - Canada) wrote:

Hi Sergey

Pablo already said a few words about using KPM

...

Does kwant ldos  function work correctly if all the leads are insulating for a given energy?  Does it correctly catch the localized states in the scattering region? kwant.ldos can only give the contribution to the ldos from the continuous part of the spectrum (i.e. the contribution from the *scattering* states). If, as you indicate, none of your leads have open propagating modes at the energy that you are looking at then kwant.ldos will give 0.

Calculating the bound state spectrum for an open system (i.e. in the presence of leads) is not completely trivial. This paper [1] outlines an algorithm to do so, but it's not foolproof (especially in the presence of almost-degeneracies). I attempted to implement this algorithm in Kwant [2], but didn't finish due to time constraints.

If I understand your problem correctly you would actually like to study an *infinite* 2D sheet? As you have discovered this is not really possible in current Kwant and you have to resort to using approximations (e.g. attaching many leads). The good news is that there is work afoot to enable these kinds of calculations in Kwant, based on the approach outlined in this paper [3]. The bad news is that these things tend to move slowly and implementing [3] in a way that is robust remains a challenge.

On balance I would suggest a large finite sample and use of KPM may be your best bet for now.

Happy Kwanting,

Joe

[1]: https://arxiv.org/abs/1711.08250 [2]: https://gitlab.kwant-project.org/kwant/kwant/-/merge_requests/320 [3]: https://arxiv.org/abs/1906.09210