Recent Changes - Search:

HomePage

PmWiki

pmwiki.org

IBC


Figure 1: Tensor network diagrams of (a) finite, (b) infinite and (c) infinite boundary condition (IBC) MPS wavefunctions.

Introduction

In this tutorial, we will introduce the tools for handling infinite boundary condition (IBC) wavefunctions [1–4].

A typical problem that we are interested in is to take a translation-invariant state (usually the ground state of some Hamiltonian), apply some local perturbation on it, and look at how that perturbation evolves with time. Now, while a translation-invariant state can be represented directly in the thermodynamic limit by using an infinite MPS (iMPS), when we apply this local perturbation to the state, it is no longer translation invariant. The usual way people go about simulating this problem is to just use a finite-size MPS, where we find the ground state for some finite system, apply the perturbation somewhere in the middle, and evolve it in time.

However, this is inefficient in a couple of ways. Firstly, we are limited by the finite size of the system: once the perturbation spreads out and reaches the boundaries, the results will start to be spoiled by the boundary effects. Hence, we must anticipate how long we want the simulation to take when deciding the system size. Secondly, as we will often have to use a very large finite system in order to reach the desired simulation times, a significant amount of time will have to be spent evolving sites far away from the perturbation wavefront, which will not change much until the perturbation reaches that site. Moreover, while not directly part of the time evolution simulation, if we prepare the initial state using finite DMRG, the time taken will scale linearly with system size, and if we are using a system size of hundreds of sites, this step could potentially be hundreds of times slower than finding the initial state using iDMRG. (There are some other more subtle inefficiencies, which we will mention later.)

A more effective solution is to use an infinite boundary condition wavefunction, which is essentially a finite MPS window coupled to two semi-infinite MPSs on its left and right boundaries. We can represent the initial perturbation as a small window only containing a single site or a few sites, and when we perform time evolution, we can incorporate more sites into the window from the boundaries as the wavefronts spreads throughout the system. As such, as long as we continue to increase the size of the window appropriately, we overcome the first inefficiency mentioned about, and since we only evolve sites contained inside of the window, we can overcome the second inefficiency.


Figure 2: As the wavefront spreads out during time evolution, we incorporate more sites into the IBC window.

Example 1: XX chain domain wall

In this example, we extend the example used in the tutorial for finite TDVP, where we evolve the domain wall initial state {$$|\Psi(t=0)\rangle = |\cdots\uparrow\uparrow\uparrow\uparrow\uparrow\downarrow\downarrow\downarrow\downarrow\downarrow\cdots\rangle.$$} under the XX Hamiltonian {$$H_\text{XX} = \sum_j (S^x_j S^x_{j+1} + S^y_j S^y_{j+1}) = \frac{1}{2} \sum_j (S^-_j S^+_{j+1} + S^+_j S^-_{j+1}).$$} In this case, we have an IBC wavefunction where the two semi-infinite boundaries are different states: one with all spins pointing up, and the other with all spins pointing down. To generate the lattice and the left and right boundary iMPS wavefunction, we use

spinchain-u1 -o lattice
mp-construct -l lattice -o psi+ 0.5
mp-construct -l lattice -o psi- " -0.5"

(We need to use " -0.5" since -0.5 would be interpreted as an option for mp-construct.) To generate the initial IBC wavefunction representing this domain wall, we use mp-ibc-create

mp-ibc-create psi+ psi- -o psi.t0.0 -H lattice:H_J1t

Note that we need to specify the Hamiltonian when specifying two different left and right boundary wavefunctions: this is since the singular values in the bond between the two boundaries need to be determined, which we do in this case by minimising the Hamiltonian (although, since the boundary wavefunctions are both product states, the only singular value is 1, so this optimisation doesn’t need to be done, but this will be required if both boundaries are nontrivial MPSs). We can then perform a time evolution simulation of this initial state using mp-ibc-tdvp

mp-ibc-tdvp -H lattice:H_J1t -w psi.t0.0 -o psi -t 0.1 -n 100 -d 1e-6

As we can see, as the simulation progresses, the (evolution) window size gets larger: the tolerance for whether new sites are added can be controlled by the --fidtol (-f) option.

We can then calculate expectation values in a similar manner as for finite wavefunctions

parallel -k 'echo -n "{1} {2} "; mp-expectation psi.t{1} lattice:"Sz({2})" -r' ::: $(seq 0.0 0.1 10.0) ::: $(seq -20 19)

In the initial state, the domain wall bond between the two boundary states is between sites -1 and 0. Note that we can calculate expectation values for any site in the system, whether or not it is in the IBC window or in one of the semi-infinite boundaries.

Example 2: Spin-1 Heisenberg chain spin correlation function

For this example, we will consider the spin-1 Heisenberg model {$$H_\text{Heisenberg} = \sum_j \mathbf{S}_j \cdot \mathbf{S}_{j+1}.$$} In this case, we need to optimise the ground state wavefunction first: for this example, a bond dimension of 50 will be sufficient,

spinchain-u1 -o lattice
mp-idmrg-s3e -H lattice:H_J1 -w psi -m 10..50x100,50x100 --create -q 0

We can create an empty IBC wavefunction using just this background with

mp-ibc-create psi -o psi-empty

Since we only have a single background wavefunction here, we do not need to optimise the singular values in the central bond. We can then create a spin excitation on top of this background with mp-ibc-apply

mp-ibc-apply -w lattice:"Sp(0)" psi-empty psi.t0.0 --normalize

Note that unlike mp-apply and mp-iapply, we use the --window (-w) flag to specify a finite operator acting on the window (we could also use --left (@-l@) to specify a string operator acting on the left boundary).

We can then perform the time evolution simulation with

mp-ibc-tdvp -H lattice:"H_J1-sum_unit(-1.4014839481575*I(0))"" -w psi.t0.0 -o psi -t 0.1 -n 100 -d 1e-3 -f 1e-8

We subtract the ground state energy density from the Hamiltonian here, otherwise the energy would decrease whenever the window expands, creating a spurious phase factor when calculating the two-point correlation function. Note that we use a larger value for the fidelity tolerance for window expansion (-f 1e-8): this is since the ground state we are using is only a finite bond dimension approximation, so it is nonstationary under time evolution with an asymptotic tolerance of around {$10^{-10}$} far from the window.

To calculate the unequal-time, two-point correlation function {$C(t, x) = \langle \Psi(0) | S^-_0 S^+_x | \Psi(t) \rangle$}, we use mp-ibc-overlap,

parallel -k 'echo -n "{1} {2} "; mp-ibc-overlap psi.t{1} psi.t0.0 --quiet --rotate {2}' ::: $(seq 0.0 0.1 10.0) ::: $(seq 0 25) > corr.dat

Here, --rotate x shifts the bra wavefunction x sites to the left before calculating the overlap (for a finite-sized system, we would have to calculate the overlap with a state with the initial excitation on site x). By spatial reflection and time reversal symmetries, {$C(t, -x) = C(t, x)$} and ${C(-t, x) = C(t, x)^*$}, so we only need to calculate the correlation function for {$t, x \geq 0$}.

We can rewrite the correlation function into the form {$C(t, x) = \langle \Psi(-t/2) | S^-_0 S^+_x | \Psi(t/2) \rangle$}, and using time reversal symmetry, we can get {$|\Psi(-t/2)\rangle$} from the complex conjugate of {$\Psi(t/2)\rangle$}: using this form, we can get the correlation function up to time {$2T$} from a time evolution simulation to time {$T$}.

parallel -k 'echo -n $(perl -e "print 2*{1}")" {2} "; mp-ibc-overlap psi.t{1} psi.t{1} --conj --quiet --rotate {2}' ::: $(seq 0.0 0.1 10.0) ::: $(seq 0 50) > corr.dat

To do this for a finite-sized system, if we wanted the correlation function for up to {$x=N_x$}, we would need to perform time evolution simulations for {$N_x$} wavefunctions with the initial excitation on site {$0, \ldots, N_x-1$}. To plot this correlation function in gnuplot (including the {$-x$} part), we can use

plot "corr.dat" u 2:1:(sqrt($3**2+$4**2)) w image, "corr.dat" u (-$2):1:(sqrt($3**2+$4**2)) w image

We can calculate the spectral function for this correlation function using a discrete Fourier transform (convoluted with a Gaussian function {$\exp(-\alpha (t/t_\text{max})^2)$} to remove the oscillations from the finite-time effects) with a short python script.

import numpy as np

##############
# USER INPUT #
##############
t_steps = 101
x_steps = 51
t_step = 0.2

# Strength of the Gaussian window.
alpha = 5
# The multiplier of the samples of the correlation function to use to create the spectral function.
sample_multiplier = 5

# Load data.
data = np.genfromtxt("corr.dat", delimiter="")

t = data[:, 0].reshape((t_steps, x_steps))
x = data[:, 1].reshape((t_steps, x_steps))
corr = (data[:, 2] + 1j*data[:, 3]).reshape((t_steps, x_steps))

# Multiply by Gaussian window.
corr *= np.exp(-alpha * (t/np.max(t))**2)

# Add -x axis.
t = np.concatenate((np.flip(t[:,1:], axis=1), t), axis=1)
x = np.concatenate((-np.flip(x[:,1:], axis=1), x), axis=1)
corr = np.concatenate((np.flip(corr[:,1:], axis=1), corr), axis=1)

# Add -t axis.
corr = np.concatenate((np.flip(corr[1:,:], axis=0).conj(), corr), axis=0)

# Calculate spectral function.
spec = np.fft.fftshift(np.fft.fft2(corr, ((2*t_steps-1)*sample_multiplier, (2*x_steps-1)*sample_multiplier)), axes=(0,))

nw, nk = np.shape(spec)
k = np.outer(np.ones(nw), np.linspace(0, 2, nk))
w = np.outer(np.linspace(-np.pi/t_step, np.pi/t_step, nw), np.ones(nk))

# Save spectral function to spec.dat.
np.savetxt("spec.dat", np.transpose([w.reshape(w.size), k.reshape(k.size), np.real(spec.reshape(spec.size)), np.imag(spec.reshape(spec.size))]))

We can then plot this spectral function in gnuplot with

set xrange [0:2]
set yrange [-2:5]
plot "spec.dat" u 2:($1/2):(sqrt($3**2+$4**2)) w image

Extra technical considerations

The output for mp-ibc-tdvp shows two different window sizes: WindowSize and EvWindowSize. The former is the total number of sites in the IBC window, but we only perform time evolution on a subset of this window, which we call the evolution window. The reason is that we want to be able to expand the left environment of the leftmost site in the evolution window, which would involve modifying the site to its left (and similarly for the RHS of the window). However, we may not want to incorporate this environment site into the evolution window just yet (there are situations where this would lead to the window to expand indefinitely), so we add it to the IBC window, while leaving it out of the evolution window.

Controlling the evolution window independently from the IBC window allows for some interesting setup. For instance, we may choose to use an evolution window of a fixed width which moves forward in line with the wavefront (this can be done with the --comoving option): this is an efficient way to get the state of the wavefront with high accuracy while neglecting the rest of the wavefunction [2,4].

An important consideration is that the size of the window does not have to be as large as the perturbed part of the wavefunction. This is since modifying a single MPS tensor does not just affect the wavefunction at that site alone, but also affects the sites around it in a region proportional to the correlation length. For instance, the initial state for the spin-1 Heisenberg simulation above has a single site window, but if you look at the correlation function at {$t=0$}, the perturbation spreads out far beyond this window. This means we can usually get away with using quite a small window, but the --fidtol parameter is often quite conservative in adding a lot of sites.

(On the other hand, for a product state background with correlation length zero, we do need the window to be a large as the perturbed part of the wavefunction.)

A final consideration is that the bond expansion shouldn’t be neglected. For these sort of local perturbations, the entanglement will be greater around the perturbation than in the background, so we will usually need a larger bond dimension within the window than in the background state. And as is usually the case with time evolution, this window bond dimension will generally need to continue expanding throughout the time evolution as the wavefunction becomes more entangled here.

References

[1] H. N. Phien, G. Vidal and I. P. McCulloch, Infinite boundary conditions for matrix product state calculations, Phys. Rev. B 86, 245107 (2012), doi:10.1103/PhysRevB.86.245107.
[2] H. N. Phien, G. Vidal and I. P. McCulloch, Dynamical windows for real-time evolution with matrix product states, Phys. Rev. B 88, 035103 (2013), doi:10.1103/PhysRevB.88.035103.
[3] A. Milsted, J. Haegeman, T. J. Osborne and F. Verstraete, Variational matrix product ansatz for nonuniform dynamics in the thermodynamic limit, Phys. Rev. B 88, 155116 (2013), doi:10.1103/PhysRevB.88.155116.
[4] V. Zauner, M. Ganahl, H. G. Evertz and T. Nishino, Time evolution within a comoving window: scaling of signal fronts and magnetization plateaus after a local quench in quantum spin chains, J. Phys.: Condens. Matter 27, 425602 (2015), doi:10.1088/0953-8984/27/42/425602.
Edit - History - Print - Recent Changes - Search
Page last modified on March 29, 2024, at 06:23 AM