Dear Kwant developers,

I was wondering if there is a limit on the resolution of Kwant and Tkwant when calculating the density of an evolved wavefunction.

I am currently working on a code to simulate the response of a carbon nanotube spin qubit to an electric dipole spin resonance drive. I am using Kwant to define the carbon nanotube as a scattering region and then apply a Hamiltonian to the system. The lowest eigenvector of the Hamiltonian at t=0 is determined and associated with a tkwant.onebody.WaveFunction. This wavefunction is evolved over time and the density of the wavefunction calculated using kwant.operator.Density.

The relevant part of my code is given below:

class System:
    def __init__(self, hamiltonian, pertubation_type, number_of_lattices, potential_type):
        self.potential_type = potential_type
        self.hamiltonian = hamiltonian
        self.number_of_lattices = number_of_lattices
        self.pertubation_type = pertubation_type
        self.evolve_state = False
       
    def potential(self, z, time):
        if self.potential_type == 0: #infinite square well
            total_potential = 0  
        elif self.potential_type == 1: #parabolic potential
            total_potential = .5 * (z * self.omega_0_au) ** 2        
        if self.pertubation_type == "cos":
            self.pertubation = self.cosine_v_ac
        else:
            self.pertubation = self.sine_v_ac
        total_potential += self.pertubation(time, z, self.eV_0_au, self.pulse_frequency_au, self.total_length_au)
        return total_potential

    def kwant_shape(self, site):
        (z,) = site.pos
        return (-self.total_length_au / 2 <= z < self.total_length_au / 2)

    def make_system(self):
        self.template = kwant.continuum.discretize(self.hamiltonian, grid=self.lattice_size_au)  
        self.syst = kwant.Builder()
        # add the nanotube to the system
        self.syst.fill(self.template, self.kwant_shape, (0,))
        self.syst = self.syst.finalized()

        def B_constant(z):
            return  -self.g * self.mu_B_au * self.B_z_au * self.hbar_au / 2
        def C_constant(z):
            return -self.g * self.mu_B_au * self.b_sl_au * z * self.hbar_au / 2
        def D_constant(z):
            return -self.g * self.mu_B_au * self.B_y_au * self.hbar_au / 2

        # coefficient for the kinetic energy term
        A_constant = self.hbar_au ** 2 / (2 * self.m_au)

        # parameters of hamiltonian
        self.params = dict(A=A_constant, V=self.potential, B=B_constant, C=C_constant, D=D_constant)
        self.tparams = self.params.copy()  # copy the params array
        self.tparams['time'] = 0  # initial time = 0

        # compute the Hamiltonian matrix
        hamiltonian = self.syst.hamiltonian_submatrix(params=self.tparams)
        # compute the eigenvalues (energies) and eigenvectors (wavefunctions).
        eigenValues, eigenVectors = np.linalg.eig(hamiltonian)

        # sort the eigenvectors and eigenvalues according the ascending eigenvalues.
        idx = eigenValues.argsort()
        self.initial_eigenvalues = eigenValues[idx]
        eigenVectors = eigenVectors[:, idx]

        # initial wave functions unperturbed by time dependent part of hamiltonian
        self.psi_1_init = eigenVectors[:, 0]
        self.psi_2_init = eigenVectors[:, 1]

        # an object representing the spin-up state
        self.spin_up_state = tkwant.onebody.WaveFunction.from_kwant(syst=self.syst,
                                                                    psi_init=self.psi_1_init,
                                                                    energy=eigenValues[0],
                                                                    params=self.params)
        return self.syst

    def evolve(self, time_steps):  
        self.evolve_state = True
       
        total_osc_time = self.second_to_au(0.5e-8)
        times_au = np.linspace(0, total_osc_time, num=time_steps)
        times_SI = np.linspace(0, self.au_to_second(total_osc_time), num=time_steps)

        density_operator = kwant.operator.Density(self.syst)  
        psi = self.spin_up_state
       
        data = {'states': []}

        for time in times_au:
            psi.evolve(time)  # evolve the wavefunction
            density = psi.evaluate(density_operator) # compute density

            data['states'].append({'pdf': density.tolist()})
       
        json_string = json.dumps(data)
        with open(r'\{}.json'.format(timestr), 'w') as outfile:
            outfile.write(json_string)

        self.evolve_state = False

        return True
   
    def main():
        lattices = 100
        potential = 0
        system = System("((A * k_z**2) + V(z, time)) * identity(2) + B(z) * sigma_z + C(z) * sigma_x + D(z) * sigma_y",
                        pertubation_type="sin", number_of_lattices=lattices,
                        potential_type=potential)
        system.make_system()

        system.evolve(10)

if __name__ == '__main__':
    main()

Please note that I haven't included some conversion functions in the above code.

Using this code, I obtain the results shown in the animation attached. The first subplot shows the density of the wavefunction across the carbon nanotube evolving with time. The animation shows an expected smooth Gaussian response for the first few time steps, however the graph then becomes very noisy. I was wondering why this occurs and if it could be due to a limit on the resolution of Tkwant?

The second subplot shows the energy of the time-dependent perturbation over time, displaying a sin wave perturbation as a simple 'see-sawing' potential. 

In addition, I was wondering if you had any advice for making the code run faster? I ran the simulation for 10 time steps over a very short time period to obtain the attached animation and it took around 4 hours on my desktop.

Thank you for any help with this.

Best wishes,
Isobel