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 <felixmende2@web.de> 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-3D-nanowire-cut/blob/4dde2854731c3f787c294779b9d13907002acba5/current_density_minimal_example.ipynb

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