Hello everyone, I am new to KWANT but have been trying to implement a tight-binding model over a (2D) square lattice under the influence of a constant perpendicular magnetic field. This example is explored in section IV B of [Graf, M., & Vogl, P. Electromagnetic fields and dielectric response in empirical tight-binding theory. Physical Review B, 51(8), 4940 (1995)]. My implementation is inspired by the KWANT tutorials. I start with the following function to set up the model: def make_system(a: float=1, eps0: float=4.0, hop: float=1.0, r: float=20): """ Creates a tight-binding model Args: a (float): Lattice constant eps0 (float): On-site energy hop (float): Hopping energy r (float): Radius of the quantum dot Returns: syst (kwant.Builder): The system lat (kwant.lattice): The lattice """ lat = kwant.lattice.square(a, norbs=1) syst = kwant.Builder() # Define the quantum dot def circle(pos): (x, y) = pos rsq = x ** 2 + y ** 2 return rsq <= r ** 2 def hopx(tosite, fromsite, peierls): phi = peierls(tosite, fromsite) return -hop * t * phi def hopy(tosite, fromsite, peierls): phi = peierls(tosite, fromsite) return -hop * t * phi syst[lat.shape(circle, (0, 0))] = eps0 * t # hoppings in x-direction syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx # hoppings in y-directions syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy # Closed system, no leads return syst, lat And then I diagonalize the resulting Hamiltonian for different values of the magnetic field (in units of the flux quantum hc/e in cgs units): def sorted_eigs(ev): evals, evecs = ev # Sort eigenvalues in ascending order and get the corresponding indices sorted_indices = np.argsort(evals) # Use the sorted indices to sort evals sorted_evals = evals[sorted_indices] # Use the sorted indices to sort eigenvectors accordingly sorted_eigenvectors = evecs[:, sorted_indices] return sorted_evals, sorted_eigenvectors for iB in np.linspace(0, 0.03, 30): B = iB def B_syst(pos): return B peierls_syst = gauge(B_syst) params = dict(peierls=peierls_syst) ham_mat = syst.hamiltonian_submatrix(params=params, sparse=True) # Diaglonalize the Hamiltonian evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=10, sigma=0.05)) If I compare the results to Fig. 2 (a) of [Graf, M., & Vogl, P. Electromagnetic fields and dielectric response in empirical tight-binding theory. Physical Review B, 51(8), 4940 (1995)], I get results that are those for a magnetic field B that is ~2 times as small. Looking at the source code, specifically line 791 of https://gitlab.kwant-project.org/kwant/kwant/blob/v1.5.0/kwant/physics/gauge... I find that the phases are written with a factor of pi in the exponent. However, looking at [Graf, M., & Vogl, P. Electromagnetic fields and dielectric response in empirical tight-binding theory. Physical Review B, 51(8), 4940 (1995)], if I express the magnetic field in units of the flux quantum / l^2 (as is suggested in the docstring at line 1005) I expect to get a factor of h / hbar = 2*pi in the exponent, consistent with the discrepancy. Is my reasoning correct? Is the KWANT definition for the flux quantum perhaps hc/(2e) (ie the superconducting flux quantum)? Any help would be greatly appreciated. Thanks.