You can calculate the particle current from site "j" to site "i"
due to a state "psi":

    I_ij = -2 Im(psi_i* H_ij psi_j)

Where "Im" is the imaginary part, "H_ij" is the Hamiltonian matrix
element between sites "i" and "j" and "*" is complex conjugation.

H_ij can be obtained from your finalized system, "sys", with
`sys.hamiltonian(i, j, args=args)`. The wavefunctions
"psi" can be obtained from `kwant.wave_function(sys, args)`.

The following Python snippet should give you what you want::

    lat = kwant.lattice.square()
    B = 1.0  # your magnetic field
    E = 1.0  # energy @ which you want the current
    sys = make_system()

    wf = kwant.wave_function(sys, energy=E, arg=(B,))

    site_i = lat(4, 5)  # whatever sites you want
    site_j = lat(4, 6)
    # get the indices of the sites in the wavefunction
    i = sys.sites.index(site_i)
    j = sys.sites.index(site_j)

    def current(psi):
        H_ij = sys.hamiltonian(i, j, args=(B,))
        return -2 * (psi[i].conjugate() * H_ij * psi[j]).imag

    I_ij = 0
    # add the contributions from all open modes in all leads
    for l, lead in enumerate(sys.leads):
        for psi in wf(l):
            I_ij += current(psi)

There are several caveats to the above:

1. It will not work if we have sites with more than 1 orbital
   as we have assumed that there is a 1-1 mapping between sites
   and elements in "psi".

2. No statistical physics has been included in the above --
   the quantity calculated is *not* strictly the particle current,
   but the "particle current at energy E". To get
   the actual particle current you need to do the statistical physics
   right. This will involve integrating over energy and weighting
   the contributions from each lead by some distribution function
   (e.g. Fermi-Dirac) to get the full particle current. Depending on
   what kind of system you are studying you may be able to justify
   foregoing the energy integration on physical grounds, but this is a
   question of the physics you are looking at and is independent of

Hope that helps,


P.S. I have not tested the above code snippet, and it's entirely
possible that I've got a sign wrong, missed a few syntax errors etc.
The point is just to give you the gist of how you would go about
implementing such a calculation.

Attachment: signature.asc
Description: OpenPGP digital signature

Reply via email to