Mailman 3 June 2021 - Kwant-discuss

Inquiry about Kwant boundstate algorithm (unphysical spectra at some sets of parameters)
by Zhang, Chi 25 Jun '21

25 Jun '21

Dear Colleagues, Hi, we are 2 students at Imperial College looking to incorporate your boundstate algorithm (https://scipost.org/submissions/1711.08250v1<https://scipost.org/submissions/1711.08250v1/>) in the context of a lateral Josephson junction. The code we used was copied and pasted from here: ( https://gitlab.kwant-project.org/jbweston/kwant/-/blob/feature/boundstate/k…) As an exercise we are trying to reproduce fig. 9 in this paper. Implementing the same parameters, we managed to obtain similar graphs shown below for a value of t=0.7. [cid:d6350847-a84d-413c-966e-9f76a7e16b14] However, we notice that at some values of t, the boundstate algorithm starts yielding non-physical results. Examples include: 1. 1. At t=-1 (which is what we assumed t to be initially from the captions in fig. 7), the algorithm returned energies above the band gap in the JJ in topological regime for some values of phi, the phase difference, and returned true values for some values of phi but not for others in the trivial regime. 2. [cid:ab6b6f7d-c949-4867-a1de-fc904e9d617f] 1. In the B scan for t=-1 boundstates start disappearing from B=0.325 onwards, well before crossing the degeneracy point, and in the topological regime they do not line up with the band edge. 2. [cid:1390ff07-aa3a-4545-ae1b-11b041706548] 1. We also noticed that the algorithm results are dependent on the system setup. To be able to independently replicate fig. 9 my partner and I set up a 2-site (left lead - super - normal - right lead) and 3-site (left lead - super - normal - super - right lead) system respectively to implement the phase from crossing the S-N interface. A schematic of this is shown below. 2. [cid:ccd40751-9e75-4c9f-9dc1-66343648e586] 3. While most standard kwant functionalities seem to recognise the 2 systems are equivalent to each other, we obtained drastically different results from the boundstate algorithm. Examples include a t-scan at B=1 to uncover possible t values where physical ABS spectra could be returned. 4. [cid:d9ca0c7e-5d9f-4546-aadb-091d36aa3542] 5. 2 site 3 site 6. (We checked that the thing causing the drastically different behaviour is the number of sites in the system. We were able to replicate each other's results only by changing the number of sites in the scattering region.) 1. We plotted f(E) in the find_boundstate function for t=-1, B=0.3 (where the algorithm succeeds in detecting the bound state) and B=0.35 (where the algorithm fails to do so), respectively. We noticed that the 2 roots from B=0.3 which are validated to be the true bound states are omitted in B=0.35, not because the validation operator rejects them, but because min_eigenvalue() (or lowest value of sparse_diag()) at these 2 energies indeed returned a value much larger than tol=1e-8 (0.003). We incremented the B field to see when the boundstate to no-boundstate transition occurs and it's at B=0.325 that the min_eigenvalue suddenly shoots up to 1e-4 (for B=0.324 we have 1e-9). We do not understand how this comes to be. 2. [cid:b7e35a85-363c-47d1-b756-52827f73d3f5] 3. The above disappearance of boundstates past B=0.325 is encountered when we set sparse=False too. We would like to know if similar problems with the boundstate algorithm have been encountered before, and if so, what do these signatures tell us is wrong, or does this indeed suggest something is incomplete in the implementation? Finally, we would also like to ask about your specific system setup in Kwant to obtain figure 9, as well as the parameter t that was used, as this was not specified explicitly. Many thanks in advance for answering our query! Sincerely, Chi Zhang & Ryan Tiew Blackett Laboratory Imperial College London

4 3

Plotting band structure of 2D continuum system in 1D strip with domain wall
by thauwa＠gmail.com 20 Jun '21

20 Jun '21

Hi all, I am somewhat new to Kwant and am trying to achieve the following: I have a 2D 3-band continuum model that I want to see edge states of using band and density diagrams via a 1D real-space strip. I make the system finite in the x direction (implying translational symmetry in the y direction, making k_y a good quantum number and so constant). I setup the domain wall by choosing between two values for a potential, dependent on whether the site is before length/2 or after. I use kwant.continuum.discretize() for builder(). However, I am having trouble with the setup in general. For one, I cannot seem to apply translational symmetry to builder() because I am not sure what vector to use (I tried (0,1) and (1,0)), but to no avail. Second, if I proceed without translational symmetry, kwant.plotter.bands(syst, show=False) gives me "TypeError: Expecting an instance of InfiniteSystem." When I browsed the email archives, people suggested syst.leads[0], etc: but the issue is that I don't have these leads explicitly and that gives me an error. My understanding is that I can get a proper band diagram showing edge states if I plot eigenvalues from looping through different values of k_y that I fixed earlier. However, I don't think I saw examples of the Kwant examples doing it this way. So, it seems as if I have confused some basic applications in Kwant. Would someone mind clarifying my confusions for me? All I want to see is the a typical band structure diagram for this 3-band continuum model in a domain wall setup. Following is my code (I know that the given Hamiltonian will not give edge states, but I just want to get the code working for starters). Thank you for your time! # Minimal example import kwant.continuum import kwant import numpy as np import matplotlib.pyplot as plt import tinyarray import scipy lat_const = 3.19 # lattice constant hami = """ [[2* t0 * ((1-(2*((1/2)*k_x*a))**2/2)+2*(1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2))+e1, M1(x,y)-2*sqrt(3)*t2*(((1/2)*k_x*a))*((momentumY*a*sqrt(3)/2)) + 2*1j*t1*((2*((1/2)*k_x*a)) + (((1/2)*k_x*a))*(1-((momentumY*a*sqrt(3)/2))**2/2)), 2*t2*((1-(2*((1/2)*k_x*a))**2/2) - (1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2)) + 2*sqrt(3)*1j*t1*(1-(((1/2)*k_x*a))**2/2)*((momentumY*a*sqrt(3)/2))], [-M1(x,y)-2*sqrt(3)*t2*(((1/2)*k_x*a))*((momentumY*a*sqrt(3)/2)) - 2*1j*t1*((2*((1/2)*k_x*a)) + (((1/2)*k_x*a))*(1-((momentumY*a*sqrt(3)/2))**2/2)), 2*t11*(1-(2*((1/2)*k_x*a))**2/2) + (t11+3*t22)*(1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2)+e2, -M2(x,y)+sqrt(3)*(t22-t11)*(((1/2)*k_x*a))*((momentumY*a*sqrt(3)/2))+4*1j*t12*(((1/2)*k_x*a))*((1-(((1/2)*k_x*a))**2/2)-(1-((momentumY*a*sqrt(3)/2))**2/2))], [2*t2*((1-(2*((1/2)*k_x*a))**2/2) - (1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2)) - 2*sqrt(3)*1j*t1*(1-(((1/2)*k_x*a))**2/2)*((momentumY*a*sqrt(3)/2)), M2(x,y)+sqrt(3)*(t22-t11)*(((1/2)*k_x*a))*((momentumY*a*sqrt(3)/2))-4*1j*t12*(((1/2)*k_x*a))*((1-(((1/2)*k_x*a))**2/2)-(1-((momentumY*a*sqrt(3)/2))**2/2)), 2*t22*(1-(2*((1/2)*k_x*a))**2/2) + (3*t11+ t22)*(1-(((1/2)*k_x*a))**2/2)*(1-((momentumY*a*sqrt(3)/2))**2/2)+e2]] """ template = kwant.continuum.discretize(hami, grid = lat_const) #print(template) L = 100 def wire_shape(site): (x, ) = site.pos return (0<=x<L) val1 = 10.07822923*1j val2 = 12.42088476*1j def pot1(x,y=0): return[0,val1][x<L/2] def pot2(x,y=0): return[0,val2][x<L/2] syst = kwant.Builder(); syst.fill(template, wire_shape,(0,)); syst = syst.finalized(); def get_ham(momY): ham = syst.hamiltonian_submatrix(params=dict(M1=pot1,M2=pot2, momentumY = momY, y=0, a = 3.19, e1 = 1.046, e2 = 2.104, t0 = -0.184, t1 = 0.401, t2 = 0.507, t11 = 0.218, t12 = 0.338, t22 = 0.057),sparse = False) return ham #plt.spy(ham) #block off-diagonal #plt.show() # Band structure diagram??? # different momentumY's to iterate over v1=-2*np.pi; v2=2*np.pi; v_segs=10; vary_range = np.linspace(v1, v2, num=v_segs) # np.linspace(kx0_1, kx0_2, num=kx_segs, retstep=True) plt.figure() for i in range(v_segs): ham = get_ham(vary_range[i]) ev,evecs = scipy.linalg.eigh(ham) plt.plot(ev) # plot eigenvalues plt.show() # Wavefunction density diagram??? Should show peak at edge states, if there are any. plt.figure() for i in range(v_segs): ham = get_ham(vary_range[i]) ev,evecs = scipy.linalg.eigh(ham) wavefunction = evecs[:, len(ham)-1] up,middle,down = wavefunction[::3],wavefunction[1::3],wavefunction[2::3] density = np.abs(up)**2 + np.abs(middle)**2 + np.abs(down)**2 plt.plot(density) ##kwant.plotter.map(syst, density, show=False) plt.show() Best, Tharindu

1 0

2D current density profile through 3D nanowire cut
by felixmende2＠web.de 18 Jun '21

18 Jun '21

Hello Kwant Community, I created a 3D nanowire grown along x direction with leads attached in the x direction and want to plot the current density profile through a cut at a x = const plane. I think this Problem was never really solved in the public community communication channels, and personally I think this is quiet useful, because people might want to show surface currents of TI's etc. I created a minimal example displaying two methods to calculate and display the current density through a plane x = const: see my jupyter notebook on GitHub: https://github.com/Quaki96/KwantQuestion_2D-current-density-profile-through… There at "# current density through x=const plane METHOD 1": I build a 2D slice of the 3D nanowire and use kwant.operator.Current(fsyst, where=cut) and specifying for which hoppings in should calculate the current density: def cut(site0,site1): return site0.pos[0]==x_cut and site1.pos[0]==x_cut+1 # calculating current just for this x hopping This yields a current density profil with only non negativ current densities even if I have no magnetic field, hence this somehow breaks time reversal symmetry. If I instead define def cut2(site0,site1): return site0.pos[0]==x_cut+1 and site1.pos[0]==x_cut # calculating current just for this x hopping I get the same picture except that the sign changed completely. "# current density through plane x = const METHOD 2 Here I use kwant.operator.Current(fsyst), and calculate the current density over all hoppings of the system. Then I use kwant.plotter.interpolate_current(fsyst, current_density) to interpolate the current density vectorfield which returns it as field, box. In a simple inefficient for loop I then create an array of y_vals and z_vals and their corresponding current density z (2) component values = values.append(field[x_cut,ny,nz,2]) at x = x_cut. This yields a current density profile which seems to respect time reversal symmetry (vanishes when B = 0) is non zero when B != 0 but because no bias is applied integration over the plain yields zero total current. I think this might be a correct implementation. But I am really not sure because I don't know what I might did wrong in the first implementation. Happy Kwanting Felix

3 8

RunTime error due to different number of left and right going modes in honeycomb lattice with SO
by 170070051＠iitb.ac.in 15 Jun '21

15 Jun '21

raise RuntimeError("Numbers of left- and right-propagating " RuntimeError: Numbers of left- and right-propagating modes differ, possibly due to a numerical instability. Hi all, I am getting the above error message while doing transport calculations in Silicene(next nearest neighbour hopping due to intrinsic SO). But it disappears when I set the the intrinsic SO to zero. What could be the possible reason. Thanks and Regards, Koustav Jana

3 7

SdH oscillation in kwant
by abd.gm90＠yandex.ru 13 Jun '21

13 Jun '21

Hi! I'm new at kwant, can anyone help me with my problem. Is it possible to get SdH oscillation for layered system with kwant? There is two type of layers in the system, normal insulators and topological insulators with edge states and we build a system by alternating different types of layers, and there is hopping between edge states through normal insulator layer (Burkov model - https://journals.aps.org/prb/abstract/10.1103/PhysRevB.85.165110 ). It will be very nice if someone give me advice

2 1