HOWTO /
IBCIntroductionIn 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 translationinvariant 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 translationinvariant 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 finitesize 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 semiinfinite 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. Example 1: XX chain domain wallIn 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 semiinfinite 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 spinchainu1 o lattice mpconstruct l lattice o psi+ 0.5 mpconstruct l lattice o psi " 0.5" (We need to use mpibccreate 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 mpibctdvp H lattice:H_J1t w psi.t0.0 o psi t 0.1 n 100 d 1e6 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 We can then calculate expectation values in a similar manner as for finite wavefunctions parallel k 'echo n "{1} {2} "; mpexpectation 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 semiinfinite boundaries. Example 2: Spin1 Heisenberg chain spin correlation functionFor this example, we will consider the spin1 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, spinchainu1 o lattice mpidmrgs3e 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 mpibccreate psi o psiempty 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 mpibcapply w lattice:"Sp(0)" psiempty psi.t0.0 normalize Note that unlike We can then perform the time evolution simulation with mpibctdvp H lattice:"H_J1sum_unit(1.4014839481575*I(0))"" w psi.t0.0 o psi t 0.1 n 100 d 1e3 f 1e8 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 twopoint correlation function.
Note that we use a larger value for the fidelity tolerance for window expansion ( To calculate the unequaltime, twopoint correlation function {$C(t, x) = \langle \Psi(0)  S^_0 S^+_x  \Psi(t) \rangle$}, we use parallel k 'echo n "{1} {2} "; mpibcoverlap psi.t{1} psi.t0.0 quiet rotate {2}' ::: $(seq 0.0 0.1 10.0) ::: $(seq 0 25) > corr.dat Here, 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} "; mpibcoverlap 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 finitesized 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_x1$}. 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 finitetime 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_steps1)*sample_multiplier, (2*x_steps1)*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 considerationsThe output for 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 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 spin1 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 (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 realtime 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/09538984/27/42/425602.
