Dear Felix, I am not sure to follow precisely your concern but I can make some comments. You are claiming that you are eliminating negative current by choosing hoppings from x to x+1. This is not true. You do not have to do it in both directions. The current is well defined after you choose your hoppings and it can be positive or negative (Local current). I also noticed that you are looking at the local current density through the wire and then plot the result on the site where the hopping starts. It is one way of looking at the current density. Kwant usually plots the current on hoppings and not on sites. If you want to look at the current density on the cut, you can do the following. 1) Create a new system that contains the slice where you do the cut. The order of the sites in this new system will be the same as in the original (Tested on a square lattice. Needs to be verified on other complex lattices). 2) Pick up the wave function in the original system and choose the values on the chosen sites and keep the order. 3) Calculate the current using these chosen wave functions and the new system. 4) Plot your current density as usual with kwant.plotter.current You can check the enclosed file for an example. You can do the same for a longitudinal cut (instead of the transverse cut.) I hope this helps, Adel On Thu, Jun 10, 2021 at 1:35 PM wrote:

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

-- Abbout Adel

Dear Felix, I understood the way you did in method one and it seems fine. For the second method, although the general idea seems correct, I doubt the accuracy of the way you handle it. You wrote : * field[x_cut,ny,nz,0]) # field at n_x = x_cut, iterating through ny, nz and adding just Jx* The shape of "Field" is not very clear to me. I do not understand exactly how it depends on the parameter n (the default value is n=9) in the documentation. So, it is not good to put x_cut as a first index since the shape of "field" changes with n. The last parameter should not refer to Jx or any other component. It seems that it depends on the dimension of the system. (You can check this function with a spinless model, in 2D and 3D) Probably @Christoph Groth (or another person) can better explain this since the documentation is poor about this function. I hope this helps, Adel On Mon, Jun 14, 2021 at 12:57 PM wrote:

Dear Adel,

Thank you for your quick reply and also for providing some example code!

You created a 3D system, then a 2D system with same in plane hoppings (2d plane cut trough the wire) and the projected the in plane current density onto these hoppings via kwant.plotter.current.

What I want is not displaying the in plane current, but the current density penetrating the 2D cut. Following the Kwant logic I would have to project onto the out of plane hoppings. Therefore I need to create a system with out of plane hoppings, which automatically is 3D. But in this case kwant.plotter.current does not work any more.

My idea in my first method was, to specify "where" in kwant.operator.Current with where= cut and cut return site0.pos[0]==x_cut and site1.pos[0]==x_cut+1 (so specifying the out of plane x hoppings).

I created a 2D system with just onsite terms and no hoppings, so that the current density in x direction is maped on the onsites (I cannot map on the x hoppings because yz 2D system has no x hoppings) Then using kwant.plotter.map to display this.

But I suspect something is wrong with this.

In my method 2 I plotted the J_z component instead of the J_x component. I updated the minimal example accordingly

https://github.com/Quaki96/KwantQuestion_2D-current-density-profile-through-...

My Question in particular is, whether my second method of doing things is correct, because in the case of no magnetic field, the current density does not vanish, which makes me somewhat sceptical.

PS. If a reader of this thread is interested in an small extension of Adels example with norbs=2 the following thread would be useful: https://mail.python.org/archives/list/kwant-discuss@python.org/thread/LJFBAX...

Best,

Felix

-- Abbout Adel

Dear Felix, The current will not vanish when you put B=0. You are calculating the contribution of the left lead alone. The one which vanishes is the measured one at equilibrium. It corresponds to the current you obtain with kwant when you sum the contributions of all the leads and then integrate them over the whole Fermi sea. If one lead is set to a slightly higher potential, then the sum of all the contributions will vanish but a small part of the order ev of that lead, and thus, the current will be proportional to this ev. You can have more details about this in the documentation of Tkwant or the thesis of Joseph Weston for example. For the sign you are mentioning (between method 1 and 2), it is a global sign, that comes from the way plotter.interpolate returns the result. (usually, the current is calculated on all the hoppings in syst.graph which contains the hoppings (a, b) as well as (b, a)) I hope this helps, Adel On Tue, Jun 15, 2021 at 11:44 AM wrote:

Dear Adel,

Thank you for your quick reply!

I am sceptical about the results my method 1 yields, because depending on how I define the "cut" function I obtain either a positiv or negativ current density profil. Furthermore I would expect for zero magnetic field no current density at all and even more critical I would expect vanishing total current (no bias applied). Because the model is totally time reversal symmetric. I made a minimal example highlighting these differences:

https://github.com/Quaki96/KwantQuestion_2D-current-density-profile-through-... the only difference here is the definition of the "cut" function either defined as def cut(site0,site1): return site0.pos[0]==x_cut and site1.pos[0]==x_cut +1 or as def cut(site0,site1): return site0.pos[0]==x_cut+1 and site1.pos[0]==x_cut

For the second method. It seems like that the interpolation increases somehow the number of sites, (probably for some interpolation reasons) but I think the labelling modulo some shift is #field[ x integer lattice position index, y integerlattice position index, z integer lattice position index, vec field]. I agree some explanation or an improved documentation would be really desirable.

Happy Kwanting

Felix

-- Abbout Adel

Hi Felix, Adel, Abbout Adel wrote:

The shape of "Field" is not very clear to me. I do not understand exactly how it depends on the parameter n (the default value is n=9) in the documentation. So, it is not good to put x_cut as a first index since the shape of "field" changes with n. The last parameter should not refer to Jx or any other component. It seems that it depends on the dimension of the system. (You can check this function with a spinless model, in 2D and 3D) Probably @Christoph Groth (or another person) can better explain this since the documentation is poor about this function.

I would be happy to improve the documentation of kwant.plotter.interpolate_current, but I am not sure that I understand what exactly is not clear to you. The function interpolate_current interpolates the current field of a n-dimensional Kwant system in real space and than samples this continuous field over a n-dimensional hyper-square lattice (i.e. cubic in 3d). This is what the docstring means when it says: “This routine samples the smoothed field on a regular (square or cubic) grid.” So, the points of the returned field do not in general correspond to Kwant sites. After all, the Kwant system does not have to live on a square lattice (or any regular lattice). Moreover, there are, depending on the ‘n’ parameter, typically more field points than system sites. Thanks to this, running the routine with the default n=9 yields a smooth field that looks nicely. The default of n=9 was chosen such that a “bump” is resolved reasonably well. Sampling the continuous interpolated current density field with a coarser grid is of course possible, but the plot will not look smooth, and while it is certainly possible to smooth it out while plotting, this smoothing will not be physically correct. Therefore it is best to display interpolated current density fields without further smoothing. I hope that now it is clear to you what interpolate_current does. Please do not hesitate to suggest clarification of the docstring. ---------------- Regarding the original problem, I expect that if there are indeed hoppings that are perpendicular to the plane, and all the other hoppings are in-plane, it should be possible to use a Kwant current operator and specify the piercing hoppings as “where”. Then, plotting of the resulting *scalar* field could be done with Kwant’s “interpolate_density”. In the general case of an arbitrary cutting plane and an arbitrary system, if “interpolate_current” is too slow for the whole system, a generalization of Adel’s suggestion should work. One creates a copy of the original system that retains only the sites and hoppings that have an influence on the cutting plane. (The influence of each hopping is restricted to “Capsule” [1] along the hopping with radius “r” that is computed as using the “relwidth” and “abswidth” parameters.) Please do let me know if anything remains unclear. Christoph [1] https://en.wikipedia.org/wiki/Capsule_(geometry)