This notebook shows how to use QuSpin to study time evolution. Below we show three examples:
The Gross-Pitaevskii equation (GPE) describes the physics of weakly-interacting bosonic systems, and is given by $$ i\partial_t\psi_j(t) = -J\left[ \psi_{j-1}(t) + \psi_{j+1}(t)\right] + \frac{1}{2}\kappa_i(j-j_0)^2\psi_j(t) + U|\psi_j(t)|^2\psi_j(t) $$ where $J$ is the hopping matrix element, $\kappa_i$ is the harmonic trap 'frequency' [we use the subindex $i$ to indicate an initial value which will play a role later on], and $U$ -- the interaction strength. The lattice sites are labelled by $j=0,\dots,L-1$, and $j_0$ is the centre of the 1d chain. We set the lattice constant to unity, and use open boundary conditions.
It will prove useful to define the GPE in vectorised form. Let $\vec \psi$ be the vector whose elements are the magnitude of the function $\psi_j$ on every site. The GPE then reduces to $$ i\partial_t\vec{\psi}(t) = H_\mathrm{sp}\vec{\psi}(t) + U \vec{\psi}^*(t)\circ \vec{\psi}(t)\circ \vec{\psi}(t)$$ where $H_\mathrm{sp}$ is a single-particle Hailtonian which contains the hopping term and the armonic potential, and the simbol $\circ$ defines element-wise multiplication: $(\vec\psi\circ\vec\phi)_j = \psi_j\phi_j$.
We start by constructing the single-particle Hamiltonian $H_\mathrm{sp}$. For the sake of saving code, it would be advantageous to view this Hamiltonian as the $t=0$ limit of a more-generic time-dependent Hamiltonian $H_\mathrm{sp}(t)$, which is defined by
$$ H_\mathrm{sp}(t) = -J\sum_{j=0}^{L-2} (a^\dagger_{j+1}a_j + \mathrm{h.c.}) + \frac{1}{2}\kappa_\mathrm{trap}(t)\sum_{j=0}^{L-1}(j-j_0)^2n_j $$$$\kappa_\mathrm{trap}(t)=(\kappa_f-\kappa_i)t/t_\mathrm{ramp}+ \kappa_i $$In the limit $t=0$, we have $\kappa_\mathrm{trap}(0) = \kappa_i $.
First, we load the required libraries and define the model parameters
In [1]:
from quspin.operators import hamiltonian # Hamiltonians and operators
from quspin.basis import boson_basis_1d # Hilbert space boson basis
from quspin.tools.evolution import evolve
import numpy as np # generic math functions
from six import iteritems # loop over elements of dictionary
import matplotlib.pyplot as plt
#
##### define model parameters #####
L=300 # system size
# calculate centre of chain
if L%2==0:
j0 = L//2-0.5 # centre of chain
else:
j0 = L//2 # centre of chain
sites=np.arange(L)-j0 # centred chain sites around zero
# static parameters
J=1.0 # hopping
U=1.0 # Bose-Hubbard interaction strength
# dynamic parameters
kappa_trap_i=0.001 # initial chemical potential
kappa_trap_f=0.0001 # final chemical potential
t_ramp=40.0/J # set total ramp time
# ramp protocol
def ramp(t,kappa_trap_i,kappa_trap_f,t_ramp):
return (kappa_trap_f - kappa_trap_i)*t/t_ramp + kappa_trap_i
# ramp protocol parameters
ramp_args=[kappa_trap_i,kappa_trap_f,t_ramp]
where we defined the function ramp() as the protocol $\kappa_\mathrm{trap}(t)$. Pay special attention to its arguments: the first argument must necessarily be the time $t$, followerd by all optional arguments (the parameters). These same parameters are then stored in the variable ramp_args.
Defining the static part of the Hamiltonian is straightforward and proceeds as before. However, due to the time dependence of $H_\mathrm{sp}(t)$, we now use a non-empty dynamic list. The structure of dynamic lists is similar to that of static lists: first comes the operator string, then the corresponding site-coupling list. The new thing is that we also need to parse the time-dependent function ramp, and its arguments ramp_args.
In [2]:
##### construct single-particle Hamiltonian #####
# define site-coupling lists
hopping=[[-J,i,(i+1)%L] for i in range(L-1)]
trap=[[0.5*(i-j0)**2,i] for i in range(L)]
# define static and dynamic lists
static=[["+-",hopping],["-+",hopping]]
dynamic=[['n',trap,ramp,ramp_args]] # <-- definition of dyamic lists
# define basis
basis = boson_basis_1d(L,Nb=1,sps=2)
# build Hamiltonian
Hsp=hamiltonian(static,dynamic,basis=basis,dtype=np.float64)
E,V=Hsp.eigsh(time=0.0,k=1,which='SA')
Let us set $t=0$ and consider the GPE for the initial trap parameter $\kappa_i=\kappa_\mathrm{trap}(0)$. The ground state of the single-particle Hamiltonian $H_\mathrm{sp}(t=0)$ is squeezed due to the trap. However, repulsive interactions have an opposite effect trying to push the particles apart. Thus, the true state of the system is a compromise of hte two effects.
Our first goal is to find the GS of the GPE, which is formally defined as the state of minimal energy: $$\vec\psi_\mathrm{GS} = \inf_{\vec{\psi}} \bigg( \vec{\psi}^t H_\mathrm{sp}(t=0)\vec{\psi} + \frac{U}{2}\sum_{j=0}^{L-1}|\psi_j|^4\bigg)$$
One way to find the configuration $\vec\psi_\mathrm{GS}$, is to solve the GPE in imaginary time ($it\to \tau$), which induces exponential decay in all modes of the system, except for the lowest-energy state. In doing so, we keep the norm of the solution fixed: $$\partial_{\tau}\vec\varphi(\tau) = -\bigg[H_\mathrm{sp}(0)\vec\varphi(\tau) + U \vec\varphi^*(\tau)\circ \vec\varphi(\tau)\circ \vec\varphi(\tau)\bigg],\qquad ||\vec\varphi(\tau)||=\mathrm{const.}$$ $$\vec{\psi}_\mathrm{GS} = \lim_{\tau\to\infty}\vec\varphi(\tau)$$
Any initial value problem requires us to pick an initial state. In the case of imaginary evolution, this state can often be arbitrary, but needs to possess the same symmetries as the true GPE ground state. Here, we choose the ground state of the single-particle Hamiltonian for an initial state, and normalise it to one particle per site. We also define the imaginary time vector tau. This array has to contain sufficiently long times so that we make sure we obtain the long imaginary time limit $\tau\to\infty$. Since imaginary time evolution is not unitary, QuSpin will be normalising the vector every $\tau$-step. Thus, one also needs to make sure these steps are small enough to avoid convergence problems of the ODE solver.
Performing imaginary time evolution is done using the evolve() method of the measurements tool. This function accepts an initial state phi0, initial time tau[0], and a time vector tau and solves the user-defined ODE GPE_imag_time. The first two arguments of this user-defined ODE function must be the time variable and the state. The parameters of the ODE are passed using the keyword argument f_params=GPE_params. To ensure the normalisation of the state at each $\tau$-step we use the flag imag_time=True. Real-valued output can be specified by real=True. Last, we request evolve() to create a generator object using the keyword argument iterate=True. Many of the keyword arguments of evolve() are the same as in the H.evolve() method of the hamiltonian class: for instance, one can choose a specific SciPy solver and its arguments, or the solver's absolute and relative tolerance.
Last, looping over the generator phi_tau we have access to the solution, which we display in a form of a movie:
In [3]:
#########################################################
##### imaginary-time evolution to compute GS of GPE #####
################################################### ######
def GPE_imag_time(tau,phi,Hsp,U):
"""
This function solves the real-valued GPE in imaginary time:
$$ -\dot\phi(\tau) = Hsp(t=0)\phi(\tau) + U |\phi(\tau)|^2 \phi(\tau) $$
"""
return -( Hsp.dot(phi,time=0) + U*np.abs(phi)**2*phi )
# define ODE parameters
GPE_params = (Hsp,U)
# define initial state to flow to GS from
phi0=V[:,0]*np.sqrt(L) # initial state normalised to 1 particle per site
# define imaginary time vector
tau=np.linspace(0.0,35.0,71)
# evolve state in imaginary time
psi_tau = evolve(phi0,tau[0],tau,GPE_imag_time,f_params=GPE_params,imag_time=True,real=True,iterate=True)
#
# display state evolution
for i,psi0 in enumerate(psi_tau):
# compute energy
E_GS=(Hsp.matrix_ele(psi0,psi0,time=0) + 0.5*U*np.sum(np.abs(psi0)**4) ).real
print('$J\\tau=%0.2f,\\ E_\\mathrm{GS}(\\tau)=%0.4fJ$'%(tau[i],E_GS) )
Next, we use our GPE ground state, to time-evolve it in real time according to the trap widening protocol $\kappa_\mathrm{trap}(t)$ hard-coded into the single-particle Hamiltonian $H_\mathrm{sp}(t)$. We proceed analogously -- first we define the real-time GPE and the time vector. In defining the GPE function, we split the ODE into a time-independent static part and a time-dependent dynamic part. The single-particle Hamiltonian for the former is accessed using the hamiltonian attribute Hsp.static which returns a sparse matrix. We can then manually add the non-linear cubic mean-field interaction term. In order to access the time-dependent part of the Hamiltonian, and evaluate it, we loop over the dynamic list Hsp.dynamic, reading off the corresponding sparse matrix Hd together with the time-dependent function f which multiplies it, and its arguments f_args. Last, we multiply the final output vector by the Schr\"odinger $-i$, which ensures the unitarity of for real-time evolution.
To perform real-time evolution we once again use the evolve() function. This time, however, since the solution of the GPE is anticipated to be complex-valued, and because we do not do imaginary time, we do not need to pass the flags real and imag_time. Instead, we decided to show the flags for the relative and absolute tolerance of the solver.
In [4]:
#########################################################
############## real-time evolution of GPE ###############
#########################################################
def GPE(time,psi):
"""
This function solves the complex-valued time-dependent GPE:
$$ i\dot\psi(t) = Hsp(t)\psi(t) + U |\psi(t)|^2 \psi(t) $$
"""
# solve static part of GPE
psi_dot = Hsp.static.dot(psi) + U*np.abs(psi)**2*psi
# solve dynamic part of GPE
for f,Hd in iteritems(Hsp.dynamic):
psi_dot += f(time)*Hd.dot(psi)
return -1j*psi_dot
# define real time vector
t=np.linspace(0.0,t_ramp,101)
# time-evolve state according to GPE
psi_t = evolve(psi0,t[0],t,GPE,iterate=True,atol=1E-12,rtol=1E-12)
#
# display state evolution
for i,psi in enumerate(psi_t):
# compute energy
E=(Hsp.matrix_ele(psi,psi,time=t[i]) + 0.5*U*np.sum(np.abs(psi)**4) ).real
print('$Jt=%0.2f,\\ E(t)-E_\\mathrm{GS}=%0.4fJ$'%(t[i],E-E_GS) )
The last example we show demonstrates how to use the hamiltonian class method evolve(), which is almost the same as the measurement function evolve(). The idea behind it is that any Hamiltonian, defins a unique unitary evolution through the Schroedinger equation.
Suppose we put aside the GPE altogether, setting $U=0$. We ask the question how does the GPE GS (obtained via imaginary time evolution) evolve under the free Hamiltonian $H_\mathrm{sp}(t)$ which widens thr trap.
Below, we show how to evolve the GPE ground state under the single-particle Hamiltonian, which des not know about the interactions. This can be though of as quenching the interaction strength $U$ to zero and observing the time evolution of the state in a slowly changing harmonic trap. More precisely, we want to solve the linear initial value problem
$$ i\partial_t\vec{\psi}(t) = H_\mathrm{sp}(t)\vec{\psi}(t),\ \ \ \vec \psi(0) = \vec\psi_\mathrm{GS} $$This time, there is no need for a user-defined function for the ODE -- Schroedinger's equation (in real and imaginary time) is provided in QuSpin by default.
In [5]:
#######################################################################################
##### quantum real time evolution from GS of GPE with single-particle Hamiltonian #####
#######################################################################################
# define real time vector
t=np.linspace(0.0,2*t_ramp,101)
# time-evolve state according to linear Hamiltonian Hsp (no need to define a GPE)
psi_sp_t = Hsp.evolve(psi0,t[0],t,iterate=True,atol=1E-12,rtol=1E-12)
#
# display state evolution
for i,psi in enumerate(psi_sp_t):
# compute energy
E=Hsp.matrix_ele(psi,psi,time=t[i]).real
print('$Jt=%0.2f,\\ E(t)-E_\\mathrm{GS}=%0.4fJ$'%(t[i],E-E_GS) )
In [ ]: