Evaluation of time-dependent correlators after a local quench in iPEPS: hole motion in the t− J model

Infinite projected entangled pair states (iPEPS) provide a convenient variational description of infinite, translationally-invariant two-dimensional quantum states. However, the simulation of local excitations is not directly possible due to the translationally-invariant ansatz. Furthermore, as iPEPS are either identical or orthogonal, expectation values between different states as required during the evaluation of non-equal-time correlators are ill-defined. Here, we show that by introducing auxiliary states on each site, it becomes possible to simulate both local excitations and evaluate non-equal-time correlators in an iPEPS setting under real-time evolution. We showcase the method by simulating the t − J model after a single hole has been placed in the half-filled antiferromagnetic background and evaluating both return probabilities and spin correlation functions, as accessible in quantum gas microscopes.


Introduction
While tensor network methods in the form of matrix-product states have become the method of choice for the simulation of one-dimensional quantum systems and provide both excellent ground-state data [1] and good accuracy for time-dependent quantities [2], the study of twodimensional systems remains more difficult. The limited system size of methods such as exact diagonalisation or matrix-product states on a cylinder [3] becomes particularly relevant when studying time-dependent correlators after local excitations, as the system must be able to accommodate the spread of those correlations over time and avoid their interaction with any boundaries. Infinite projected entangled pair states [4][5][6] (iPEPS) on the other hand allow for the simulation of ground-state properties of infinite two-dimensional systems with high accuracy by repeating a finite unit cell of tensors infinitely in both directions. iPEPS were also recently shown to allow for the simulation of global quenches [7][8][9] at least for short times. This simulation of a real-time evolution following a global quantum quench is relatively straightforward: evolution methods exist [10][11][12], the quench can be enacted by a change of the Hamiltonian governing this evolution and translational invariance is retained. Equal-time correlators can also be evaluated as usual for each of the computed time-evolved post-quench states.
However, when attempting to simulate a local quench and evaluate non-equal-time correlators, one encounters two problems: First, it is not possible to simply apply an operator (such asĉ † ) to a single site of the quantum state to create the local excitation: To follow this route, one would have to apply this operator to a specific site, repeated on each unit cell. While making the unit cell itself relatively large is feasible, in this case one merely recovers the case of a finite PEPS calculation and loses the inherent infinity of the iPEPS ansatz. The handling of fermionic commutation rules further complicates this approach.
Second, when pursuing this avenue to simulate the evolution of many excitations -one per unit cell -over time, it is then still not possible to evaluate non-equal-time correlators: These correlators are calculated as expectation values between two different quantum states. However, evaluating the norms of those states will yield either 0 or 1 in the thermodynamic limit and the scale of the correlator is hence not known. In comparison, equal-time correlators are evaluated as Ô (t) = ψ(t)|Ô|ψ(t) ψ(t)|ψ(t) , but the denumerator is clearly ill-defined for a correlator Ô (t , t) between two different infinite quantum states |ψ(t ) and |ψ(t) .
Here, we avoid both problems by adding one auxiliary site to each of the physical sites of our system while preserving translational invariance. We demonstrate the method by evaluating the return probability and diagonal-spin-correlators of a single hole in the twodimensional antiferromagnetic background of the t − J model [13][14][15][16][17][18].

Local excitations and non-equal-time correlators
Consider a system composed of physical local state spaces H p i repeated on each site i of an infinite lattice. We will later focus on the case of a square two-dimensional lattice, but the method likewise applies to other lattice geometries. The total Hilbert space is the tensor product of the local spaces, We can represent a translationally invariant quantum state |ψ p ∈ H p using a tensor network ansatz if it has low entanglement, which is typically true for ground states of local Hamiltonians. If |ψ p is only invariant under translation by multiple sites (such as e.g. an antiferromagnetic state under translation by two instead of one site), we can also capture this by using a sufficiently large unit cell of tensors in the ansatz.
To simulate a local excitation without breaking translational invariance, we now create a translationally invariant superposition of excitations on top of our initial state, simulate the time evolution of this superposition under some HamiltonianĤ and then select the part of the superposition which contains an excitation at a specific local site [19,20].
To create the superposition of local excitations, one could apply e.g. 1 + x p i with some creation or annihilation operatorx p i and a small prefactor governing the density of excitations on each site asŶ = If we let this operator act on our quantum state, we obtain a superposition By including a suitable operator (e.g. the particle number operator) in expectation values later, we can select one of the states with an excitation (e.g. a hole at a particular site), which is most likely one of the summands in the second term if is small. Crucially, we can also do so after a real-time evolution ofŶ |ψ p , in this way post-selecting the evolution of a single excitation out of the translationally invariant background. This approach usingŶ has two downsides: First, the operatorx p i alone typically breaks some symmetry of the system such as spin projection, particle conservation or fermionic parity. While the former two merely lead to a less efficient simulation (as those symmetries then cannot be used in the tensor network ansatz), the breaking of fermionic parity is a serious problem which makes the simulation of fermionic systems impossible. Furthermore, while it is possible to post-select a quantum state with an excitation present at a particular site after the time evolution, we cannot post-select for a state where the excitation was created at a particular site initially.
To circumvent both problems, we add an auxiliary state space H a i of the same dimension as H p i to each site of our lattice. The total Hilbert space H is then defined as the tensor product of the auxiliary and physical tensor product spaces on each lattice site The initial quantum state |ψ p is extended by a suitably-chosen empty quantum state |0 a to form a state in the full Hilbert space |ψ = |ψ p ⊗ |0 a . In the case of the t − J model, for example, |0 a is the state with zero particles on each site in the auxiliary system. The HamiltonianĤ used for the time evolution still only acts on the physical system. We then replace the excitation operatorŶ by a form which conserves all symmetries of the system, namelyX = where for convenience with existing implementations, we then instead use the local exponential Instead of creating excitations from nothing asŶ did,X now moves (e.g.) particles from the physical to the auxiliary system and thereby creates an excitation in the physical sector. The density of particles moved and hence the density of local excitations is given by , ideally we want to consider the case → 0. No symmetry is broken during this process if we account for auxiliary particles in the same way as we account for physical particles andX hence leaves the fermionic parity of the state well-defined. Additionally, it is now possible to not only post-select based on the physical state of some particular site (to select an excitation present there after the evolution), but also to post-select based on the auxiliary state of some particular site. Because there are no dynamics in the auxiliary layer, the auxiliary state at time t is equal to the auxiliary state at time 0 and hence allows for the selection of an excitation which was created at a particular site initially.

Application to the t − J model
Specifically, we consider the two-dimensional t − J model on the square lattice with a local physical three-dimensional state space Taking a second such space H a i increases the local physical dimension of the iPEPS tensor from three to nine, but iPEPS methods scale favourably in this dimension, so this is not a concern. Letĉ p( †) iσ annihilate (create) a physical fermion on site i with spin σ, letŝ iσ annihilate (create) an auxiliary fermion on site i with spin σ. Finally, letn p i (n a i ) denote the particle number operator (0 or 1) on the physical (auxiliary) site i.
The Hamiltonian acts on the physical sector only and is the standard t − J Hamiltonian linking all nearestneighbour sites i, j . Here, we fix t = 1 and J = 1 /3. Now take |GS to be an approximation of the infinite ground state ofĤ at a given iPEPS bond dimension D and half-filling (one fermion per site) in the physical sector, with the auxiliary sector being entirely empty: The physical ground state |GS p is simply the ground-state of the Heisenberg Hamiltonian, which can be reasonably well approximated by a D = 4 or D = 5 iPEPS (other states may of course require a larger bond dimension). This state breaks translational invariance, so we use a 2 × 2 unit cell. It preserves both U(1) N particle number and U(1) S z spin-projection symmetry and we make use of both [21]. Fermionic commutation relations are ensured using the fermionic tensor network ansatz [22,23] as implemented in SyTen's STensor class [24,25].
Given |GS as described above, we create the initial excitation with the operator This operator will move particles from the occupied physical sector to the empty auxiliary sector and results in new state |ψ(0) with a finite hole density on each physical site. Evolving this state under the physical HamiltonianĤ is straightforward and for a given time t results in a state |ψ(t) = e −itĤ |ψ(0) .
In the following, we are particularly interested in (a) the return probability p R (t) of a hole to its creation site and (b) the diagonal spin-spin correlator z diag (t) at time t with a hole present at time t between the two spins. The return probability p R (t) is given by where the numerator evaluates the joint probability of a hole created at site i (via the density on the auxiliary site,n a i ) present there at a later time (via the density on the physical site, n p i ) with the denumerator conditioning on the initial creation of a hole at this site. As the hole density is low, we neglect the case of the hole created at site i moving away and another hole created at some neighbouring site j taking its place.
These correlators are sketched in Fig. 1. Note that, if desired and with larger computational effort, it would be conceivable to repeat the same calculation at different values of and subsequently extrapolate → 0.

Results
In the following, we apply the method described above to evaluate the return probability and diagonal-nearest-neighbour spin correlators in the t − J model after the effective introduction of a single hole. We also simulate this system using time-dependent matrix-product states [2] on cylinders of width 4 and 6 to obtain comparison data for short times.
Time-dependent matrix-product states on cylindrical geometries are used to provide comparison data, assumed to be valid at least for short times when the finite circumference of the cylinders is not yet relevant. We compute the ground-states of the t − J model at halffilling and apply an excitationĉ 0,↑ +ĉ 0,↓ in the centre of the system. The resulting excited state is then time-evolved with either the 2TDVP [26] or the MPO W II method [27][28][29] using the SyTen [24,25] and TeNPy toolkits [30] respectively. The return probability is given simply as 1 −n 0 (t) . On cylinders of width W = 4, convergence is easy to achieve at modest bond dimensions m = 1000, increasing the bond dimension further (up to m = 5000) does not lead to different results. As the MPS bond dimension scales exponentially with the circumference of the cylinder, convergence is more difficult on W = 6 cylinders. Running the time evolution at the same fixed bond dimension as the initial ground state does not converge well. Preparing the initial ground state at a smaller bond dimension 200 and then running the time evolution at bond dimension m = 1000 leads to results at least on short times very similar to the W = 4 cylinder (cf. Fig. 2), which is expected as the short-time dynamics are independent of the spin background and hence governed by the hole motion only. Departing from the short-time regime, however, the results become uncontrolled. Increasing the bond dimension further or evolving with the same bond dimension as the initial state does not lead to good convergence. Additionally, while the hole spreads isotropically along the x-and y-direction on the W = 4 cylinder, this is not the case on the W = 6 cylinder (not shown).
Overall, we only obtain reliable data for the return probability on cylinders of width W = 4 and qualitative data for cylinders of width W = 6.
In the iPEPS simulation, we use the fast full update (FFU, [11,12]) to obtain the initial ground state and perform the subsequent evolution with the simple update (SU). While the (fast) full update would be able to make better use of the bond dimension of our state, we have encountered some stability issues [8] resulting from this update method which lead to very limited time scales. The simple update may not make perfect use of the iPEPS bond dimension but, given a sufficiently large bond dimension, still provides good results without any of the stability issues observed with the FFU. We prepare the initial (ground) state at an initial bond dimension D = 4 and create an excitation density of 10 −2 . During the subsequent real-time evolution, we allow a range of bond dimensions D = 4, . . . , 16. We focus on even bond dimensions D, as odd bond dimensions show slightly worse convergence behaviour due to truncation within spin multiplets. Future computational and algorithmic advances may make bond dimensions D > 17 possible. We use a time step size δt = 0.01 together with a second-order Trotter decomposition of the time-evolution operator.
Exploratory calculations at D = 5 and/or hole density ≈ 10 −4 result in decreased hole mobility at a given evolution bond dimension D as the competition between spin and hole entanglement during the iPEPS state truncation favour the spin sector disproportionally when it is initially more strongly entanglend (D = 5) or there are fewer holes. Hole mobility still increases when increasing the evolution bond dimension D, but convergence is much slower than when starting with D = 4.
Expectation values are calculated using the corner transfer matrix at increasing bond dimensions χ until the difference between results of two successive dimensions χ and 2χ are sufficiently small; error bars are smaller than symbol sizes in all cases. Fig. 3 and Fig. 4 show the short-time dynamics of the return probability p R (t) and diagonal spin-spin correlator z diag (t) calculated with iPEPS. We observe good convergence in the bond dimension starting from D ≥ 8 for short times. There, the td-MPS results are reproduced. In particular, the motion of the hole away from its initial site on times of the order of the nearest-neighbour hopping is captured well. At the same time, z diag (t) becomes negative because the moving hole distorts the original antiferromagnetic background. Hence, spin   Fig. 3 for longer times t ≥ 1. The return probability shows qualitative features common to all calculations at large bond dimensions, but quantitative convergence is difficult. The revival around t ≈ 3.5 is not expected and likely due to limited entanglement in our ansatz.
correlators between both originally nearest-neighbour and originally next-nearest-neighbour fermions contribute to z diag (t). The stronger nearest-neighbour correlators then dominate the sum and cause the observed sign change. Because the SU(2)-spin symmetry is spontaneously broken along the preferred z-axis in the iPEPS calculation but still present in the finite td-MPS calculations, a comparison of numerical values is not meaningful in this case.
For longer times, convergence is very difficult, as our ansatz is inherently limited in entanglement and -due to the simple update -does not make optimal use of the available bond dimension. However, the first revival of the return probability observed in the td-MPS data is still reproduced well by the iPEPS results around t ≈ 1.5, cf. Fig. 5. The iPEPS data also contains a second, much larger revival at later times t ≈ 3.5 which is not observed in the td-MPS data and not physically expected either (instead we expect the hole to move away from its creation point with frustrated spins left behind healed by spin flips [13]). At the moment, it is unclear whether this revival is due to limited entanglement in the iPEPS ansatz which hinders healing of frustrated spins through spin-exchange interactions and hence increases the cost of moving the hole further from its origin or a side-effect of the typically overestimated magnetisation in the iPEPS ground state which may lead to more Ising-like physics.

Conclusion
We have shown that both the simulation of local excitations and the evaluation of timedependent correlators is possible within the iPEPS formalism. Our predictions, such as the sign-change of diagonal correlators around the hole in Fig. 4, can already be tested in state-ofthe-art quantum-gas microscopes [31][32][33][34]. Future work using an environment-based truncation scheme such as the FFU together with a stabilised environment (e.g. as introduced in Ref. [35]) will be in a position to make much better use of the available bond dimension than the simple update employed here and hence will be able to analyse the physics of the system for longer times and potentially open an alternative avenue [36] to obtaining spectral functions of twodimensional systems.