Numerical continuum tensor networks in two dimensions

We describe the use of tensor networks to numerically determine wave functions of interacting two-dimensional fermionic models in the continuum limit. We use two different tensor network states: one based on the numerical continuum limit of fermionic projected entangled pair states obtained via a tensor network formulation of multi-grid, and another based on the combination of the fermionic projected entangled pair state with layers of isometric coarse-graining transformations. We first benchmark our approach on the two-dimensional free Fermi gas then proceed to study the two-dimensional interacting Fermi gas with an attractive interaction in the unitary limit, using tensor networks on grids with up to 1000 sites.


I. INTRODUCTION
Understanding the collective behavior of quantum manybody systems is a central theme in physics. While it is often discussed using lattice models, there are systems where a continuum description is essential. One such case is found in superfluids, 1 where recent progress in precise experiments on ultracold atomic Fermi gases has opened up new opportunities to probe key aspects of the phases. 2 On the theoretical side, this requires solving a continuum fermionic quantum many-body problem. For example, various quantum Monte Carlo (QMC) methods have been applied to study the cross-over from Bardeen-Cooper-Schrieffer superfluidity to Bose-Einstein condensation in two-dimensional Fermi gases. [3][4][5][6] However, the applicability of (unbiased) quantum Monte Carlo is restricted to special parameter regimes due to the fermion sign problem. 7 Thus, devising numerical methods that can address general continuum quantum many-body physics remains an important objective.
Tensor network states (TNS) are classes of variational states that have become widely used in quantum lattice models. They are complementary to QMC methods as TNS algorithms are typically formulated without incurring a sign problem. In 1D, matrix product states (MPS) now provide almost exact numerical results via the DMRG algorithm. 8 In 2D, reaching a similar level of success has been harder, but much progress has been made using projected entangled pair states (PEPS), 9 which generalize MPS to higher dimensions in a natural fashion. PEPS calculations now provide accurate results for a broad range of quantum lattice problems [10][11][12][13][14] and there have been many developments to extend the range of the techniques, for example to long-range Hamiltonians, [15][16][17] thermal states, 18 and real-time dynamics. 19 Also, much work has been devoted to improving the numerical efficiency and stability of PEPS computations. [20][21][22][23][24] Formulating tensor network states and the associated algorithms in the continuum remains a challenge. In 1D, so-called continuous MPS 25 provide an analytical ansatz in the continuum and have been applied to several problems, including 1D interacting bosons/fermions and quantum field theories. [26][27][28][29] Alternatively, the continuum description can be reached by taking the numerical limit of a set of tensor network states formulated on lattices with a discretization parameter , for → 0. This kind of numerical continuum MPS calculation has also been demonstrated in conjunction with a variety of optimization algorithms. [30][31][32] In two dimensions, despite several proposals, 25,31,33,34 the appropriate analytical form of the continuum PEPS ansatz remains unclear. In this work, we carry out continuum tensor network calculations in 2D by taking the numerical limit of a lattice discretization parameter. We explore two types of ansatz to approach the continuum. The first uses the numerical continuum limit of the lattice fermionic PEPS. Here, to connect the lattice PEPS at different scales when taking this limit and to ensure an efficient optimization on finer scales we use a multi-grid like algorithm (a generalization of the MPS multigrid algorithm). The second is based on a combination of fermionic PEPS with a tree of isometries that successively coarse grains the continuum into discrete lattices. Using these 2D numerical continuum tensor network states, we demonstrate how we can study fermionic physics in the continuum limit, applying the ansatz both to the challenging (for tensor networks) case of the free Fermi gas, as well as the attractive interacting Fermi gas in the unitary limit that can be realized in ultracold atom experiments. In the first case, we can benchmark against exact results, while in the second we can perform direct comparisons of the tensor network results to recent QMC calculations at half-filling (where there is no sign problem) in the continuum and thermodynamic limits.
The remainder of the paper is organized as follows. We first introduce the fermionic continuum Hamiltonian of interest in Sec. II and describe how to discretize it in a manner consistent with open boundary conditions in Sec. III. The two types of fermionic tensor network states are discussed in Sec. IV along with the optimization algorithms used for them. We then present our numerical benchmarks for the free and interacting Fermi gases in Sec. V. We summarize our work in Sec. VI and discuss future research directions.

II. INTERACTING FERMI GAS
For concreteness, it is useful to define a particular continuum model. In this work we will consider the free and interacting Fermi gases in two dimensions. The Hamiltonian is 1. (Color online) (a) A two-dimensional square L × L box containing spinful fermions. The fermions are confined to the box by an infinite potential at the walls of the box. (b) Lattice discretization of the L×L box into 8×8 grid points. Using Dirichlet boundary conditions, the wave function is zero at the gray points (boundary points). The active points where the wave function takes non-trivial values are defined on the 6 × 6 lattice, with a lattice spacing of = L 7 . The grid points with blue and orange colors require special treatment in the fourth-order finite difference approximation, in order to provide an accurate discretization.
given by where ψ † σ (r) and ψ σ (r) are fermionic field operators, creating and annihilating a fermion with spin σ at position r, respectively. The fermionic field operators satisfy the anticommutation relation {ψ † σ (r), ψ σ (r )} = δ(r − r )δ σσ . The coupling parameters µ and g denote the chemical potential (which controls the number of particles) and strength of interaction in the system. We assume the system is confined in a L × L square box (0 < r x , r y < L), so that the potential outside the box is V (r) = ∞. When g = 0, we have a free fermion gas confined to the box.

III. HAMILTONIAN LATTICE DISCRETIZATION
Because we define the continuum properties as a numerical limit, we need to first discretize the continuum Hamiltonian H. To do so, we replace the continuum space L × L box by a lattice containing (N + 1) × (N + 1) grid points with lattice spacing = L N . Due to the Dirichlet (open) boundary conditions, the wave function is zero on the boundary of the grid. Thus the non-trivial part of the quantum state is defined on (N − 1) × (N − 1) grid points (we refer to these as the "active points"). The lattice discretization is shown in Figs. 1(a, b).
The kinetic energy operator can be represented using a finite difference stencil on the grid. Using 2nd and 4th order finite difference approximations, the Laplacian operator where c, c † are fermionic lattice operators and the symbols ij and ij denote nearest neighbor and third nearest neighbor pairs.
Because the wave function is not smooth at the edge of the box, some care must be taken in applying the stencils to ensure that boundary errors of lower order in than implied by the stencil formula do not appear. In the 2nd order approximation, the second-derivative takes the form ∼ −4ψ0+ψ1+ψ−1 2 where the indices on ψ denote the x coordinate. Taking the index −1 to refer to the left boundary, we see that representing the wave function on only the interior N − 1 active points, while using the spacing = 1/(N + 1) (rather than the naive spacing of = 1/N ) is consistent with choosing the boundary condition ψ −1 = 0 (and similarly for the right boundary). Using this definition of = 1/(N + 1) we obtain the full 2nd order lattice discretized Hamiltonian as whereḡ is a regularized δ function interaction parameter, whose regularization procedure is described in detail in Secs. V B and A. The same continuum wave function |Ψ can also be represented by a single fPEPS on a coarse lattice connected to layers of isometric tensors (here two layers are shown). We refer to this as fPEPS-tree.
For the 4th order approximation, the second-derivative takes the form . Again, taking index −1 to refer to the left boundary, we see that the value of ψ −2 is left unspecified. To maintain the accuracy of the finite difference expression, we should choose ψ −2 to smoothly continue the wavefunction past the boundary. In our case, we choose ψ −2 = −ψ 0 (i.e. the wavefunction is antisymmetric around the boundary). This means that at the boundary, the second-derivative should be replaced by . Continuing this argument, the coefficient 60 12 2 in Eq. 2 is replaced with different values depending on the nature of the boundary points: at the red, blue and green points shown in Fig. 1(b), the coefficients become 58 12 2 , 59 12 2 and 60 12 2 . The final form of the 4th order discretized lattice Hamiltonian thus becomes with = 1/(N +1) and where the colors red, blue and orange correspond to the colored grid points in Fig. 1 To illustrate the importance of the correct representation of the Laplacian for Dirichlet boundary conditions, in Fig. 2 we show the relative error in the ground-state energy of two free fermions in the box, using different treatments of the Laplacian, as a function of . In Fig. 2(a), we compare the lattice spacing = 1/(N + 1) that is consistent with the boundary conditions to the naive spacing = 1/N for the Hamiltonian H (2) , showing that the quadratic convergence in is achieved only for the former spacing. In Fig. 2(b) we show the effect of using an antisymmetric continuation of the wavefunction in H (4) compared to simply setting the value of the wave function outside of the box to zero; much faster convergence is obtained using the antisymmetric continuation.

IV. CONTINUUM FERMIONIC TENSOR NETWORK ANSATZ
We explore two different fermionic tensor networks to approach the continuum limit ground state |Ψ , depicted in Fig. 3(a, b). One ( Fig. 3(a)) is a standard fermionic PEPS (fPEPS) defined by a set of local tensors connected by virtual bonds corresponding to the geometry of the lattice discretization. The associated bond dimension of the virtual bonds is denoted by D and controls the accuracy of the fPEPS ansatz. To enforce fermion statistics (i) all fPEPS local tensors are set to be symmetric under the action of Z 2 or U (1) symmetry groups, and (ii) each line crossing in the network is replaced by a fermionic swap gate. Such an fPEPS captures fermionic states obeying an entanglement area law. As the lattice spacing goes to zero, the fPEPS then provides a numerical representation of the continuum limit. The primary numerical challenge is ensuring that the fPEPS tensors on the finest scales are properly optimized. This can be done by taking the numerical limit, i.e. connecting fPEPS representations at different discretizations, using a multi-grid algorithm discussed further below. We refer to this numerical continuum representation by fPEPS as fPEPS-fine.
The second consists of several layers of isometric tensors, with a fPEPS placed at the top, see Fig. 3(b). We denote this ansatz an fPEPS-tree. The isometric tensors are chosen to map four fermion sites onto one effective fermion site with bond dimension χ o . The isometries describe a coarsegraining transformation, where the parameter χ o controls the accuracy of the transformation. The amount of entanglement in the ansatz is controlled by the bond dimension of the topmost fPEPS, i.e. D, which encodes quantum entanglement between effective fermions on the coarsest level. The flexibility of the ansatz is thus controlled by both {χ o , D}. We expect this representation to work well in the dilute regime, where the effective area occupied by each fermion represents a coarse-grained length-scale and the isometries connect the finest (continuum) scale to that length-scale; in general, we expect χ o , D ≈ e ρ , where ρ is the particle density. Although it is formally desirable, during coarse-graining, to decouple entanglement at each length-scale (as is the basis of fermionic multi-scale entanglement renormalization (fMERA)) 35,36 the computational cost to work with fMERA is much higher than that of the fPEPS-tree. Thus, in fPEPS-tree we account for the short-range entanglement entirely within the topmost fPEPS tensors.

A. Tensor optimization and contraction techniques
We now briefly summarize some of the techniques used to optimize and contract the different types of tensors appearing in the two ansatzes above. The fPEPS tensors are optimized towards a representation of the ground state using imaginary-time evolution, using the "full-update" method to perform bond truncations. [37][38][39] The full-update builds the environment from the entire wave function around each bond before truncation. We use the single-layer boundary contraction The isometric tensors in the fPEPS-tree are optimized using techniques similar to those used for the fMERA as described in Ref. 36. These are based on linearizing the respective cost functions with respect to the isometric tensors. In the fPEPS-tree, contracting a single layer of isometric tensors costs O(χ 9 o ). One advantage of the fPEPS-tree is that after a single isometric layer contraction, the fourth-order discretized Hamiltonian (H (4) ) is renormalized into a nearest-neighbour Hamiltonian, which then retains its nearest neighbour form through subsequent isometric layers. This simplifies the optimization of the topmost fPEPS layer, which can be performed using standard nearest-neighbour imaginary-time evolution.
To improve efficiency, we exploit Z 2 and U (1) symmetry in all tensors. We adopt the techniques developed in Refs. 41 and 42 to implement U (1) symmetry, choosing relevant symmetric sectors during the optimization. We use a simple-update strategy (based on a direct SVD decomposition) to obtain an initial guess for the symmetry sectors, and those are then further dynamically updated during the full-update optimization by using a similar strategy.

B. Multigrid fPEPS-fine optimization
Although the fPEPS-fine wave function on the finest lattice is a straightforward representation of the (near)-continuum wave function, direct optimization of such an fPEPS leads to numerical difficulties, such as slow convergence and be- ing stuck in local minima. This is analogous to what is seen in MPS simulations on very fine lattices 31 and also what is seen in solutions of partial differential equations on fine grids. Consequently, it is necessary to construct the fPEPS on finer lattices from those on coarser lattices, which can be seen as taking the continuum limit on the fPEPS tensors in an algorithmic sense. This we achieve using a multigrid-inspired algorithm.
The main idea in the multigrid approach is to interleave optimization and interpolation steps for the fPEPS tensors that are determined on lattices with different discretizations. In our version of the multigrid algorithm (i) we first approximate the ground state on the coarsest level by using a N ×N fPEPS ansatz where N ∈ {2, 3}, (ii) we attach a layer of isometric tensors to the N × N fPEPS to create an fPEPS-tree with a single layer of isometries for the 2N ×2N lattice, and we subsequently perform energy optimization of the isometries and fPEPS tensors, (iii) we use a splitting map to map the fPEPStree to a 2N × 2N fPEPS, then we relax the energy again on this finer lattice, (iv) we repeat steps (ii) and (iii) until the de- sired discretization level is reached, yielding the final fPEPSfine wavefunction. A schematic of the multigrid algorithm is depicted in Fig. 4. Note that the above is only one realization of a multigrid algorithm and many alternative choices can be made.
Two steps need to be further specified, namely (i) initializing the isometric tensors, and (ii) constructing the splitting map. We initialize the isometric tensors by diagonalizing the local Hamiltonian, defined on the 2 × 2 finer lattice, and picking the χ o lowest-energy eigenvectors. This provides an approximate initial guess for the isometric tensors. The key steps in the splitting map are shown in Fig. 4. We first add a resolution of the identity U † U = I onto the virtual bonds of the fPEPS, where U is an isometric matrix with dimension D 2 × D. Such an isometric tensor U (shown in red) splits one virtual bond into two virtual bonds, without changing the overall tensor network state. Then, one approximately solves the equation: The parameter controlling the accuracy of this transformation is the bond dimension of the virtual bonds connecting the tensors, denoted χ r ; as χ r → ∞ the transformation becomes exact. The tensors U can in principle be considered to be variational parameters in order to best satisfy the above equation. However, in this paper, we fix the form of U ; some of the diagonal elements are set to 1, and the rest are set to 0, i.e. U ij,m = δ i×D+j,m . This appears sufficient to obtain our desired accuracy (see Sec. V). To solve Eq. 3, we carry out sequential SVD to obtain guesses for four resulting tensors, and then direct optimize the fidelity using alternating least squares to improve the accuracy of the splitting. In Fig. 5, we provide the details of this procedure.
To summarize, the accuracy of calculations with fPEPSfine using the multi-grid algorithm is controlled by four parameters: the initial bond dimension D; the boundary bond dimension χ b controlling the accuracy of the environment contractions; and the accuracy of the splitting map controlled by parameters χ o and χ r , denoting the bond dimension of the isometries and the bond dimension of the resulting fPEPS on the finer lattice. In practice, we can reasonably set parameters χ b ∼ D 2 and χ r ∼ D, thus the essential controlling parameters are only D, χ o . As an example of the bond dimensions used in this paper, for the 16 × 16 lattice fPEPS-fine simulation we used (D, χ b ) = (9, 200) and (D, χ b ) = (12, 250) for Z 2 and U (1) symmetries, respectively.

A. Free Fermi gas benchmarks
We first assess the accuracy of the wave functions and algorithms discussed above using the free Fermi gas (g = 0) as a benchmark system. This system is exactly solvable by reduction to single-particle quantities but is a challenging problem for tensor networks in the continuum and thermodynamic limit as the ground-state violates the entanglement area law by logarithmic terms, i.e. ∼ ρ 1 2 A log A, where ρ is the particle density and A is the boundary length. 43,44 Consequently, to obtain accurate results, large bond dimensions in all steps of the algorithms are required and approximation errors are magnified. In the calculations below, we shall use a fixed box side-length of L = 6.
We first compute numerical continuum results using the fPEPS-tree ansatz using the 4th order spinless Hamiltonian discretization. Here we use an fPEPS-tree with two layers of isometries. In Fig. 6(a), we show the relative error of the ground-state energy ∆E as a function of particle density ρ. As expected, the relative error increases sharply when we increase the particle density using fixed bond dimensions (D, χ o ) = (8,16). However, the inset shows that going to a finer lattice does not affect the accuracy significantly (the error increases only slightly in the continuum limit). In Fig. 6(b), we plot the ground-state energy E versus the lattice spacing , where we use the fourth-order polynomial function E( ) = E →∞ + b −4 to extract the numerical continuum limit of the ground-state energy. Note that due to the use of coarse-graining, we no longer observe a simple 4th order discretization error. For each data point for a specific lattice spacing , we also perform a bond dimension extrapolation (χ o → ∞). A second-order polynomial function E(χ o ) = E χo→∞ + aχ −1 o + bχ −2 o is used to estimate the extrapolated results, as shown in the inset in Fig. 6(b). Because the energy is a function of the bond dimensions (D, χ o , χ b ), to obtain a good extrapolated estimation, we need to make sure the results are converged with respect to the bond dimensions (D, χ b ). Thus at each χ o , we use as large (D, χ b ) as possible (up to (12, 250)). The extrapolated results are shown in Table I. We find accuracies of roughly 1% or better are obtained with these bond dimensions for the lowest 3 densities where ρ < 0.4.
We next explore the behaviour of the fPEPS-fine ansatz in the continuum free fermion problem. As the multi-grid optimization involves several steps, we first illustrate the accuracy and errors associated with the individual steps. We start with the accuracy of the splitting map between coarser and finer lattices. This incurs an error which can be measured as and E f are the energies per lattice site before and after the splitting map, respectively.
As χ r increases we expect this error to go to zero, and for a practical method, it is important that a moderate χ r ∼ D is sufficient for good accuracy. In Fig. 7, we plot ∆ versus the in-ρ fPEPS-tree, H (4) I. Ground-state energy of the free Fermi gas for different particle densities in the continuum limit. Data with the symbol are for spin-1 2 fermions, and data with no such symbol are for spinless fermions. Data is better converged (and extrapolations are more accurate) at low density.
ternal bond dimension χ r . We see that the relative error drops when increasing χ r so that when χ r ∼ D it is ∼ 10 −2 −10 −3 . One might expect that if we proceed to finer lattice spacings , (i.e. where the fine lattice has more sites), the relative error might increase due to the accumulation of individual tensor splitting errors. This is in fact seen in Fig. 7, although the relative error increases quite slowly with . We also find that the accuracy of the splitting map does not depend on the particle density, as the relative error remains the same in the left and right panels. We further illustrate the numerical behavior of the multigrid algorithm in Fig. 8(a). As discussed above, applying isometries to the fPEPS to obtain a single-layer fPEPS-tree provides the initial guess for the splitting map (fPEPS-fine) with subsequent optimization being carried out after the splitting map is performed. We see that the ground-state energy jumps due to the infidelity of the splitting map, however, further optimization of the fPEPS on the finer level rapidly improves the energy. To show the importance of the multi-grid algorithm, in Fig. 8(b), we compare the fPEPS energies initialized by a simple-update method 45 and the ones initialized from coarser lattices by the multi-grid algorithm. We observe that the fPEPS energies obtained by the multi-grid algorithm are much more accurate.
Using the above multigrid calculation with fPEPS-fine with up to four layers (a fine lattice of 16 × 16 sites) and the 2nd order Hamiltonian discretization, we can estimate the energy of the fPEPS-fine ansatz in the numerical continuum limit. We use a linear extrapolation in 1 D (using a few largest values) to estimate the large-D limit for each lattice spacing. 46 A polynomial function E( ) = E →∞ + b −2 + c −4 is used to estimate the numerical continuum limit for fPEPS-tree. The estimated energies are shown in Table I. We see that in the very dilute regime, the fPEPS-tree is slightly more accurate than the multigrid algorithm, likely because it uses the 4th order Hamiltonian discretization. However, in the non-dilute regime, the multigrid algorithm performs significantly better. Indeed for densities ρ < 0.4, the fPEPS-fine results are accurate to better than 0.4%.

B. Interacting Fermi gas
We next present calculations using the fPEPS-fine ansatz for the spin-balanced interacting (g < 0) Fermi gas. We fo-   cus on the spin-balanced regime because at this point auxiliary field quantum Monte Carlo has no sign problem, and thus provides a reliable comparison; however, away from this point, sign problems manifest, for which the methods developed here remain suitable. Rather than using the interaction g appearing in the continuum Hamiltonian, the interaction strength is commonly parametrized using the dimensionless coupling parameter η = , where e f is the non-interacting Fermi energy and e b is the two-particle binding energy. We will be interested in the simultaneous continuum and thermodynamic limit, which is obtained in the simultaneous limit of infinite particle number and lattice size N e , N → ∞. For each finite lattice size N , an effective interactionḡ can be specified to be consistent with a given η, which defines the numerical fPEPS-fine lattice Hamiltonian. The procedure to determineḡ is given in Appendix B, and follows that in Ref. 6. The physics of the system is governed by η and in the limits η 0, η 0, the system is in the Bose-Einstein condensate (BEC) and Bardeen-Cooper-Schrieffer (BCS) phases, respectively. Recent QMC studies have suggested that the BCS-BEC crossover occurs near η ∼ 1.
We present fPEPS-fine results for the ground-state energy (for various N and N e ) in units of the non-interacting Fermi energy E = H Nee f /2 . E = lim Ne,N →∞ E is the desired continuum and thermodynamic limit result. We compare to results from auxiliary field QMC (AFQMC) which is exact up to statistical error (for given N e , N ). The AFQMC data is computed using periodic boundary conditions whereas the fPEPS-fine energies are computed for open boundary conditions, however, in the thermodynamic limit they should approach the same value. In Fig. 9(a, b), the ground-state TN results as a function of N e , N are compared against the AFQMC data (shown as lines, extrapolated to the thermodynamic limit for each N e ). To extrapolate the TN data to the thermodynamic limit, we use the second-order polynomial function E = E N →∞ + bN −1 + cN −2 . For the point η = −0.5 ( Fig. 9(b)), the thermodynamic and continuum limits are rapidly approached as seen in both the TN and AFQMC data; however, this is more challenging for the crossover point η = 1.0 ( Fig. 9(b)) where there are sizable finite-size effects in both the TN and AFQMC results. However, by using up to N e = 40 particles, and lattices with up to 32 × 32 sites, the continuum TN data and AFQMC extrapolations provide consistent estimates.
We further show the extrapolated continuum and thermodynamic limit fPEPS-fine and AFQMC ground-state energies across a range of renormalized coupling parameters, in Fig. 9(c). We find good agreement between the fPEPS-fine and AFQMC energies, with the largest errors around the transition point η ∼ 1.0, which again comes mainly from the uncertainty in the N e , N → ∞ extrapolations required in both methods.

VI. CONCLUSIONS
We have described the numerical continuum limit of twodimensional tensor network states based on two variants of projected entangled pair states, as well as the numerical algorithms used to work with them. Using continuum grids with up to approximately 1000 sites, our initial calculations show promising results for two fermionic continuum systems in two dimensions: the entanglement law violating free Fermi gas, as well as the interacting unitary Fermi gas of much interest in cold atom experiments. In the latter case, our results compare well to auxiliary field quantum Monte Carlo calculations that are feasible at the spin-balanced point. However, the strength of the continuum tensor network approach is that it is not limited to special points in the phase diagram. This opens up the use of tensor networks to address unresolved questions in spin-polarized Fermi gases, [47][48][49] as well as in other problem areas, such as the realistic description of electronic structure with tensor networks, and the numerical study of field theories in two-dimensions and higher.