<jupyter_start><jupyter_text>Nonlinear time evolution and the GPEThis notebook shows how to use QuSpin to study time evolution. Below we show three examples:* imaginary time evolution to find the lowest energy state of the GPE* real time evolution with a user-defined function for a non-linear ODE* unitary quantum evolutionThe 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<jupyter_code>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]<jupyter_output><empty_output><jupyter_text>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`.<jupyter_code>##### 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')<jupyter_output>Hermiticity check passed!
Symmetry checks passed!
Particle conservation check passed!<jupyter_text>Imaginary Time EvolutionLet 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:<jupyter_code>#########################################################
##### 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) )<jupyter_output>$J\tau=0.00,\ E_\mathrm{GS}(\tau)=2095.4510J$
$J\tau=0.50,\ E_\mathrm{GS}(\tau)=864.1360J$
$J\tau=1.00,\ E_\mathrm{GS}(\tau)=551.3869J$
$J\tau=1.50,\ E_\mathrm{GS}(\tau)=395.5782J$
$J\tau=2.00,\ E_\mathrm{GS}(\tau)=298.5153J$
$J\tau=2.50,\ E_\mathrm{GS}(\tau)=231.0927J$
$J\tau=3.00,\ E_\mathrm{GS}(\tau)=181.1400J$
$J\tau=3.50,\ E_\mathrm{GS}(\tau)=142.5232J$
$J\tau=4.00,\ E_\mathrm{GS}(\tau)=111.7613J$
$J\tau=4.50,\ E_\mathrm{GS}(\tau)=86.7082J$
$J\tau=5.00,\ E_\mathrm{GS}(\tau)=65.9578J$
$J\tau=5.50,\ E_\mathrm{GS}(\tau)=48.5440J$
$J\tau=6.00,\ E_\mathrm{GS}(\tau)=33.7774J$
$J\tau=6.50,\ E_\mathrm{GS}(\tau)=21.1510J$
$J\tau=7.00,\ E_\mathrm{GS}(\tau)=10.2816J$
$J\tau=7.50,\ E_\mathrm{GS}(\tau)=0.8739J$
$J\tau=8.00,\ E_\mathrm{GS}(\tau)=-7.3046J$
$J\tau=8.50,\ E_\mathrm{GS}(\tau)=-14.4398J$
$J\tau=9.00,\ E_\mathrm{GS}(\tau)=-20.6824J$
$J\tau=9.50,\ E_\mathrm{GS}(\tau)=-26.1563J$
$J\tau=10.00,\ E_\mathrm{GS}(\tau)=-30.9646J$
$J\tau=10.50,\ E_\mathrm{GS}(\tau)=-35.1936J$
$J\tau=11.00,\ E[...]<jupyter_text>Real Time Evolution with User-Specified ODENext, 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.<jupyter_code>#########################################################
############## 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) )<jupyter_output>$Jt=0.00,\ E(t)-E_\mathrm{GS}=-0.0000J$
$Jt=0.40,\ E(t)-E_\mathrm{GS}=-1.2925J$
$Jt=0.80,\ E(t)-E_\mathrm{GS}=-2.5853J$
$Jt=1.20,\ E(t)-E_\mathrm{GS}=-3.8787J$
$Jt=1.60,\ E(t)-E_\mathrm{GS}=-5.1731J$
$Jt=2.00,\ E(t)-E_\mathrm{GS}=-6.4687J$
$Jt=2.40,\ E(t)-E_\mathrm{GS}=-7.7658J$
$Jt=2.80,\ E(t)-E_\mathrm{GS}=-9.0649J$
$Jt=3.20,\ E(t)-E_\mathrm{GS}=-10.3663J$
$Jt=3.60,\ E(t)-E_\mathrm{GS}=-11.6702J$
$Jt=4.00,\ E(t)-E_\mathrm{GS}=-12.9771J$
$Jt=4.40,\ E(t)-E_\mathrm{GS}=-14.2872J$
$Jt=4.80,\ E(t)-E_\mathrm{GS}=-15.6011J$
$Jt=5.20,\ E(t)-E_\mathrm{GS}=-16.9190J$
$Jt=5.60,\ E(t)-E_\mathrm{GS}=-18.2412J$
$Jt=6.00,\ E(t)-E_\mathrm{GS}=-19.5683J$
$Jt=6.40,\ E(t)-E_\mathrm{GS}=-20.9005J$
$Jt=6.80,\ E(t)-E_\mathrm{GS}=-22.2383J$
$Jt=7.20,\ E(t)-E_\mathrm{GS}=-23.5820J$
$Jt=7.60,\ E(t)-E_\mathrm{GS}=-24.9320J$
$Jt=8.00,\ E(t)-E_\mathrm{GS}=-26.2887J$
$Jt=8.40,\ E(t)-E_\mathrm{GS}=-27.6525J$
$Jt=8.80,\ E(t)-E_\mathrm{GS}=-29.0238J$
$Jt=9.20,\ E(t)-E_\mathrm{GS}=-30.4030J$
$Jt=9.60,\ E(t)-E_\mathr[...]<jupyter_text>Real Time Evolution under a HamiltonianThe 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.<jupyter_code>#######################################################################################
##### 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) )<jupyter_output>$Jt=0.00,\ E(t)-E_\mathrm{GS}=-390.2907J$
$Jt=0.80,\ E(t)-E_\mathrm{GS}=-392.8745J$
$Jt=1.60,\ E(t)-E_\mathrm{GS}=-395.4518J$
$Jt=2.40,\ E(t)-E_\mathrm{GS}=-398.0162J$
$Jt=3.20,\ E(t)-E_\mathrm{GS}=-400.5615J$
$Jt=4.00,\ E(t)-E_\mathrm{GS}=-403.0818J$
$Jt=4.80,\ E(t)-E_\mathrm{GS}=-405.5714J$
$Jt=5.60,\ E(t)-E_\mathrm{GS}=-408.0248J$
$Jt=6.40,\ E(t)-E_\mathrm{GS}=-410.4371J$
$Jt=7.20,\ E(t)-E_\mathrm{GS}=-412.8033J$
$Jt=8.00,\ E(t)-E_\mathrm{GS}=-415.1192J$
$Jt=8.80,\ E(t)-E_\mathrm{GS}=-417.3806J$
$Jt=9.60,\ E(t)-E_\mathrm{GS}=-419.5838J$
$Jt=10.40,\ E(t)-E_\mathrm{GS}=-421.7255J$
$Jt=11.20,\ E(t)-E_\mathrm{GS}=-423.8026J$
$Jt=12.00,\ E(t)-E_\mathrm{GS}=-425.8125J$
$Jt=12.80,\ E(t)-E_\mathrm{GS}=-427.7530J$
$Jt=13.60,\ E(t)-E_\mathrm{GS}=-429.6222J$
$Jt=14.40,\ E(t)-E_\mathrm{GS}=-431.4184J$
$Jt=15.20,\ E(t)-E_\mathrm{GS}=-433.1405J$
$Jt=16.00,\ E(t)-E_\mathrm{GS}=-434.7875J$
$Jt=16.80,\ E(t)-E_\mathrm{GS}=-436.3588J$
$Jt=17.60,\ E(t)-E_\mathrm{GS}=-437.8541J$
$Jt=18.40,\ E(t)-E_\math[...]
{'fraction_non_alphanumeric': 0.18768347160178686, 'fraction_numerical': 0.052010210593490744, 'mean_word_length': 4.889139421270199, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '39928', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6049, 'n_tokens_neox': 5517, 'n_words': 1617, 'path': 'examples/notebooks/GPE.ipynb', 'repo_name': 'nelsond/QuSpin', 'size': '15670'}