Examples¶

Expectation values of the Hamiltonian and kinetic operators¶

The following code block gives a simple example of initializing a state and calculating the expectation values of the Hamiltonian and kinetic operators and the norm of the state after the evolution.

import numpy as np
from trottersuzuki import *

grid = Lattice2D(256, 15)  # create a 2D lattice

potential = HarmonicPotential(grid, 1, 1)  # define an symmetric harmonic potential with unit frequecy
particle_mass = 1.
hamiltonian = Hamiltonian(grid, potential, particle_mass)      # define the Hamiltonian:

frequency = 1
state = GaussianState(grid, frequency)  # define gaussian wave function state: we choose the ground state of the Hamiltonian

time_of_single_iteration = 1.e-4
solver = Solver(grid, state, hamiltonian, time_of_single_iteration)  # define the solver

# get some expected values from the initial state
print("norm: ", solver.get_squared_norm())
print("Total energy: ", solver.get_total_energy())
print("Kinetic energy: ", solver.get_kinetic_energy())

number_of_iterations = 1000
solver.evolve(number_of_iterations)  # evolve the state of 1000 iterations

# get some expected values from the evolved state
print("norm: ", solver.get_squared_norm())
print("Total energy: ", solver.get_total_energy())
print("Kinetic energy: ", solver.get_kinetic_energy())

Imaginary time evolution to approximate the ground-state energy¶

import numpy as np
from trottersuzuki import *

grid = Lattice2D(256, 15)  # create a 2D lattice

potential = HarmonicPotential(grid, 1, 1)  # define an symmetric harmonic potential with unit frequecy
particle_mass = 1.
hamiltonian = Hamiltonian(grid, potential, particle_mass)  # define the Hamiltonian:

frequency = 3
state = GaussianState(grid, frequency)  # define gaussian wave function state: we choose the ground state of the Hamiltonian

time_of_single_iteration = 1.e-4
solver = Solver(grid, state, hamiltonian, time_of_single_iteration)  # define the solver

# get some expected values from the initial state
print("norm: ", solver.get_squared_norm())
print("Total energy: ", solver.get_total_energy())
print("Kinetic energy: ", solver.get_kinetic_energy())

number_of_iterations = 40000
imaginary_evolution = true
solver.evolve(number_of_iterations, imaginary_evolution)      # evolve the state of 40000 iterations

# get some expected values from the evolved state
print("norm: ", solver.get_squared_norm())
print("Total energy: ", solver.get_total_energy())
print("Kinetic energy: ", solver.get_kinetic_energy())

Imprinting of a vortex in a Bose-Einstein Condensate¶

import numpy as np
import trottersuzuki as ts

grid = ts.Lattice2D(256, 15)  # create a 2D lattice

potential = HarmonicPotential(grid, 1, 1)  # define an symmetric harmonic potential with unit frequecy
particle_mass = 1.
coupling_intra_particle_interaction = 100.
hamiltonian = Hamiltonian(grid, potential, particle_mass, coupling_intra_particle_interaction)  # define the Hamiltonian:

frequency = 1
state = GaussianState(grid, frequency)  # define gaussian wave function state: we choose the ground state of the Hamiltonian
def vortex(x, y):  # vortex to be imprinted
    z = x + 1j*y
    angle = np.angle(z)
    return np.exp(1j * angle)

state.imprint(vortex)  # imprint the vortex on the condensate

time_of_single_iteration = 1.e-4
solver = Solver(grid, state, hamiltonian, time_of_single_iteration)  # define the solver

Dark Soliton Generation in Bose-Einstein Condensate using Phase Imprinting¶

This example simulates the evolution of a dark soliton in a Bose-Einstein Condensate. For a more detailed description, refer to this notebook.

from __future__ import print_function
import numpy as np
import trottersuzuki as ts
from matplotlib import pyplot as plt

grid = ts.Lattice2D(300, 50.)  # # create a 2D lattice

potential = ts.HarmonicPotential(grid, 1., 1./np.sqrt(2.))  # create an harmonic potential
coupling = 1.2097e3
hamiltonian = ts.Hamiltonian(grid, potential, 1., coupling)  # create the Hamiltonian

state = ts.GaussianState(grid, 0.05)  # create the initial state
solver = ts.Solver(grid, state, hamiltonian, 1.e-4)  # initialize the solver
solver.evolve(10000, True)  # evolve the state towards the ground state

density = state.get_particle_density()
plt.pcolor(density)  # plot the particle denisity
plt.show()

def dark_soliton(x,y):  # define phase imprinting that will create the dark soliton
    a = 1.98128
    theta = 1.5*np.pi
    return np.exp(1j* (theta * 0.5 * (1. + np.tanh(-a * x))))

state.imprint(dark_soliton)  # phase imprinting
solver.evolve(1000)  # perform a real time evolution

density = state.get_particle_density()
plt.pcolor(density)  # plot the particle denisity
plt.show()

The results are the following plots:

Imaginary time evolution in a 1D lattice using radial coordinate¶

import matplotlib.pyplot as plt
import numpy as np
import trottersuzuki as ts

angular_momentum = 1 # quantum number
radius = 100 # Physical radius
dim = 50 # Lattice points
time_step = 1e-1
fontsize = 16

# Set up the system
grid = ts.Lattice1D(dim, radius, False, "cylindrical")

def const_state(r):
    return 1./np.sqrt(radius)

state = ts.State(grid, angular_momentum)
state.init_state(const_state)

def pot_func(r,z):
    return 0.

potential = ts.Potential(grid)
potential.init_potential(pot_func)

hamiltonian = ts.Hamiltonian(grid, potential)
solver = ts.Solver(grid, state, hamiltonian, time_step)

# Evolve the system
solver.evolve(30000, True)

# Compare the calculated wave functions with respect to the groundstate function
psi = np.sqrt(state.get_particle_density()[0])
psi = psi / np.linalg.norm(psi)

groundstate = ts.BesselState(grid, angular_momentum)
groundstate_psi = np.sqrt(groundstate.get_particle_density()[0])

# Plot wave functions
plt.plot(grid.get_x_axis(), psi, 'o')
plt.plot(grid.get_x_axis(), groundstate_psi)
plt.grid()
plt.xlim(0,radius)
plt.xlabel('r', fontsize = fontsize)
plt.ylabel(r'$\psi$', fontsize = fontsize)
plt.show()