On the generalization of the exponential basis for tensor network representations of long-range interactions in two and three dimensions

TL;DR

This paper extends the exponential basis for tensor network representations of long-range interactions from 1D to higher dimensions, enabling efficient simulations of 2D and 3D quantum systems with long-range potentials.

Contribution

The authors generalize the exponential basis to 2D and 3D using Green's functions of the Helmholtz equation, providing compact tensor network representations for long-range interactions.

Findings

Tensor network representations for 2D and 3D long-range interactions are derived.

The 3D construction yields a D=4 tensor network for Coulomb-like potentials.

The method enables efficient simulation of realistic materials with long-range interactions.

Abstract

In one dimension (1D), a general decaying long-range interaction can be fit to a sum of exponential interactions $e^{-\lambda r_{ij}}$ with varying exponents $\lambda$, each of which can be represented by a simple matrix product operator (MPO) with bond dimension $D=3$. Using this technique, efficient and accurate simulations of 1D quantum systems with long-range interactions can be performed using matrix product states (MPS). However, the extension of this construction to higher dimensions is not obvious. We report how to generalize the exponential basis to 2D and 3D by defining the basis functions as the Green's functions of the discretized Helmholtz equation for different Helmholtz parameters $\lambda$, a construction which is valid for lattices of any spatial dimension. Compact tensor network representations can then be found for the discretized Green's functions, by expressing them as correlation functions of auxiliary fermionic fields with nearest neighbor interactions via Grassmann Gaussian integration. Interestingly, this analytic construction in 3D yields a $D=4$ tensor network representation of correlation functions which (asymptotically) decay as the inverse distance ($r^{-1}_{ij}$), thus generating the (screened) Coulomb potential on a cubic lattice. These techniques will be useful in tensor network simulations of realistic materials.

Tables

Table 1. Table 1: The rules for the full 3D PEPO that generates all pairwise interactions in Eq. ( 100 ) with D O = 5 subscript 𝐷 𝑂 5 D_{O}=5 . All combinations of indices not listed in this table correspond to O n k ​ n k ′ [ k ] = 0 ^ k subscript superscript 𝑂 delimited-[] 𝑘 subscript 𝑛 𝑘 subscript superscript 𝑛 ′ 𝑘 subscript ^ 0 𝑘 O^{[k]}_{n_{k}n^{\prime}_{k}}=\hat{0}_{k} , while I k subscript 𝐼 𝑘 I_{k} is simply the identity operator. Note that the local operators A ^ k subscript ^ 𝐴 𝑘 \hat{A}_{k} and B ^ k subscript ^ 𝐵 𝑘 \hat{B}_{k} do not have to be the same, although in our case from the main text they would both be the number operator n k subscript 𝑛 𝑘 n_{k} .

Rule number	Index values $(L_k,U_k,D_k,R_k,T_k,B_k)$	$O_{n_k{n}_k^{\prime }}^{\left[k\right]}$
1	(0,0,0,0,0,0)	$I_k$
2	(0,2,2,0,1,0)	$I_k$
3	(2,1,0,2,1,0)	$I_k$
4	(0,0,0,0,2,2)	$I_k$
5	(1,1,0,1,1,0)	$I_k$
6	(0,1,1,0,1,0)	$I_k$
7	(0,0,0,0,1,1)	$I_k$
8	(0,2,0,0,0,0)	${\hat{A}}_k$
9	(0,1,0,2,0,0)	${\hat{A}}_k$
10	(0,0,0,0,2,0)	${\hat{A}}_k$
11	(0,1,2,1,1,0)	${\hat{B}}_k$
12	(2,1,0,1,1,0)	${\hat{B}}_k$
13	(0,1,0,1,1,2)	${\hat{B}}_k$
14	(0,1,2,2,1,0)	$I_k$
15	(0,2,0,0,1,2)	$I_k$
16	(0,1,0,2,1,2)	$I_k$
17	(3,1,0,3,1,0)	$I_k$
18	(3,1,2,1,1,0)	$I_k$
19	(0,1,0,3,1,0)	${\hat{B}}_k$
20	(3,1,0,1,1,2)	$I_k$
21	(0,1,4,0,1,2)	$I_k$
22	(0,4,4,0,1,0)	$I_k$
23	(3,4,0,1,1,0)	$I_k$
24	(0,4,0,2,1,0)	$I_k$
25	(0,4,0,1,1,0)	${\hat{B}}_k$
${26}^{\ast }$	$P_{0,0,0,0,0,0}^{\text{upper right top corner}}$	${\hat{0}}_k$

Equations

$$\hat{V}=(\hat{W}[1]\hat{W}[2]\cdots \hat{W}[N]{)}_{11},$$

$$\displaystyle\hat{W}[i]=\left[\begin{array}[]{ccc}I&e^{-\lambda l/2}n_{i}&0\\ 0&e^{-\lambda l}I&e^{-\lambda l/2}n_{i}\\ 0&0&I\\ \end{array}\right].$$

$$(-{\nabla }^2+{\lambda }^2)G(\mathbf{r},{\mathbf{r}}^{\prime })=\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),$$

$$({\mathbf{K}}_d+{\lambda }^2{a}^2\mathbf{I})\mathbf{V}=\mathbf{I},$$

$$\displaystyle\mathbf{K}_{1}=\left[\begin{array}[]{cccccc}2&-1&0&\cdots&0&0\\ -1&2&-1&\cdots&0&0\\ 0&-1&2&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\cdots&2&-1\\ 0&0&0&\cdots&-1&2\\ \end{array}\right]_{N\times N}.$$

$$G({\mathbf{r}}_i,{\mathbf{r}}_j)$$

$$V_{ij}$$

$$Z$$

$$\displaystyle\mathbf{K}_{1}^{\prime}=\left[\begin{array}[]{cccccc}k_{b}&k_{c}&0&\cdots&0&0\\ k_{c}&k_{i}&k_{c}&\cdots&0&0\\ 0&k_{c}&k_{i}&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\cdots&k_{i}&k_{c}\\ 0&0&0&\cdots&k_{c}&k_{b}\\ \end{array}\right]_{n\times n},$$

$$k_b=\frac{1}{2}(1+\coth \xi ),k_c=-\frac{1}{2}\mathrm{csch}\xi ,k_i=\coth \xi ,\xi =\lambda l.$$

$$e^{{\bar{c}}_i{c}_j+{\bar{c}}_j{c}_i}=1+{\bar{c}}_i{c}_j+{\bar{c}}_j{c}_i+{\bar{c}}_i{c}_j{\bar{c}}_j{c}_i=\sum _{m=1}^4{\alpha }_{i,m}{\beta }_{j,m}.$$

$${\alpha }_i$$

$${\beta }_j$$

$$Z$$

$$(A_k{)}_{lr}$$

$$(B_k{)}_{lr}$$

$$(C_k{)}_{lr}$$

$$V_{ij}=\frac{1}{Z}(A_1\cdots B_i\cdots C_j\cdots A_N{)}_{11},Z=(A_1\cdots A_N{)}_{11}.$$

$$\displaystyle\hat{W}[i]=\left[\begin{array}[]{ccc}A_{i}\otimes I&B_{i}\otimes n_{i}&0\\ 0&A_{i}\otimes I&C_{i}\otimes n_{i}\\ 0&0&A_{i}\otimes I\\ \end{array}\right].$$

$$(A_k{)}_{lr}$$

$$(B_k{)}_{lr}$$

$$(C_k{)}_{lr}$$

$$V_{ij}$$

$$(O_k{)}_{dlur}=\int \mathrm{d}{\bar{c}}_k\mathrm{d}{c}_k{Q}_k{\delta }_{k,d}{\beta }_{k,l}{o}_k{\gamma }_{k,u}{\alpha }_{k,r},$$

$$(P_k{)}_{ab}={\delta }_{ab}(-1{)}^{p({\gamma }_{k,a})}=\mathrm{diag}(1,-1,-1,1),$$

$$Z$$

$$\displaystyle\begin{array}[]{llll}\alpha_{13}&\beta_{14}\alpha_{14}&\beta_{15}\alpha_{15}&\beta_{16}\\ \alpha_{9}(\gamma_{9}\delta_{13})&\beta_{10}\alpha_{10}(\gamma_{10}\delta_{14})&\beta_{11}\alpha_{11}(\gamma_{11}\delta_{15})&\beta_{12}(\gamma_{12}\delta_{16})\\ \alpha_{5}(\gamma_{5}\delta_{9})&\beta_{6}\alpha_{6}(\gamma_{6}\delta_{10})&\beta_{7}\alpha_{7}(\gamma_{7}\delta_{11})&\beta_{8}(\gamma_{8}\delta_{12})\\ \alpha_{1}(\gamma_{1}\delta_{5})&\beta_{2}\alpha_{2}(\gamma_{2}\delta_{6})&\beta_{3}\alpha_{3}(\gamma_{3}\delta_{7})&\beta_{4}(\gamma_{4}\delta_{8})\\ \end{array}$$

$$\displaystyle\begin{array}[]{llll}\alpha_{13}&\beta_{14}\alpha_{14}&\beta_{15}\alpha_{15}&\beta_{16}\\ \alpha_{9}(\gamma_{9}\delta_{13})&\beta_{10}\bar{c}_{10}\alpha_{10}(\gamma_{10}\delta_{14})&\beta_{11}\alpha_{11}(\gamma_{11}\delta_{15})&\beta_{12}(\gamma_{12}\delta_{16})\\ \alpha_{5}(\gamma_{5}\delta_{9})&\beta_{6}\alpha_{6}(\gamma_{6}\delta_{10})&\beta_{7}c_{7}\alpha_{7}(\gamma_{7}\delta_{11})&\beta_{8}(\gamma_{8}\delta_{12})\\ \alpha_{1}(\gamma_{1}\delta_{5})&\beta_{2}\alpha_{2}(\gamma_{2}\delta_{6})&\beta_{3}\alpha_{3}(\gamma_{3}\delta_{7})&\beta_{4}(\gamma_{4}\delta_{8})\\ \end{array}$$

$$(O_k{)}_{dblutr}=\int \mathrm{d}{\bar{c}}_k\mathrm{d}{c}_k{Q}_k{\delta }_{k,d}{\tau }_{k,b}{\beta }_{k,l}{o}_k{\gamma }_{k,u}{\sigma }_{k,t}{\alpha }_{k,r},$$

$$S_{yz}^{wx}={\delta }_{wz}{\delta }_{xy}(-1{)}^{p(w)p(x)},$$

$$Z$$

$$\displaystyle\begin{array}[]{lll}\alpha_{16}&\beta_{17}\alpha_{17}&\beta_{18}\\ \alpha_{13}(\gamma_{13}\delta_{16})&\beta_{14}\alpha_{14}(\gamma_{14}\delta_{17})&\beta_{15}(\gamma_{15}\delta_{18})\\ \alpha_{10}(\gamma_{10}\delta_{13})&\beta_{11}\alpha_{11}(\gamma_{11}\delta_{14})&\beta_{12}(\gamma_{12}\delta_{15})\\ \alpha_{7}(\sigma_{7}\tau_{16})&\beta_{8}\alpha_{8}(\sigma_{8}\tau_{17})&\beta_{9}(\sigma_{9}\tau_{18})\\ \alpha_{4}(\gamma_{4}(\sigma_{4}\tau_{13})\delta_{7})&\beta_{5}\alpha_{5}(\gamma_{5}(\sigma_{5}\tau_{14})\delta_{8})&\beta_{6}(\gamma_{6}(\sigma_{6}\tau_{15})\delta_{9})\\ \alpha_{1}(\gamma_{1}(\sigma_{1}\tau_{10})\delta_{4})&\beta_{2}\alpha_{2}(\gamma_{2}(\sigma_{2}\tau_{11})\delta_{5})&\beta_{3}(\gamma_{3}(\sigma_{3}\tau_{12})\delta_{6})\end{array}$$

$$\displaystyle\begin{array}[]{lll}\delta_{16}\alpha_{16}&\beta_{17}\alpha_{17}&\beta_{18}\\ \delta_{13}\gamma_{13}\alpha_{13}&\beta_{14}\alpha_{14}(\gamma_{14}\delta_{17})&\beta_{15}(\gamma_{15}\delta_{18})\\ \gamma_{10}\alpha_{10}&\beta_{11}\alpha_{11}(\gamma_{11}\delta_{14})&\beta_{12}(\gamma_{12}\delta_{15})\\ \delta_{7}\alpha_{7}(\sigma_{7}\tau_{16})&\beta_{8}\alpha_{8}(\sigma_{8}\tau_{17})&\beta_{9}(\sigma_{9}\tau_{18})\\ \delta_{4}\gamma_{4}\alpha_{4}(\sigma_{4}\tau_{13})&\beta_{5}\alpha_{5}(\gamma_{5}(\sigma_{5}\tau_{14})\delta_{8})&\beta_{6}(\gamma_{6}(\sigma_{6}\tau_{15})\delta_{9})\\ \gamma_{1}\alpha_{1}(\sigma_{1}\tau_{10})&\beta_{2}\alpha_{2}(\gamma_{2}(\sigma_{2}\tau_{11})\delta_{5})&\beta_{3}(\gamma_{3}(\sigma_{3}\tau_{12})\delta_{6})\end{array}$$

$$\displaystyle\begin{array}[]{lll}\tau_{16}\delta_{16}\alpha_{16}&\beta_{17}\alpha_{17f}&\beta_{18}\\ \tau_{13}\delta_{13}\gamma_{13}\alpha_{13}&\beta_{14}\alpha_{14}(\gamma_{14}\delta_{17})&\beta_{15}(\gamma_{15}\delta_{18})\\ \tau_{10}\gamma_{10}\alpha_{10}&\beta_{11}\alpha_{11}(\gamma_{11}\delta_{14})&\beta_{12}(\gamma_{12}\delta_{15})\\ \delta_{7}\alpha_{7}\sigma_{7}&\beta_{8}\alpha_{8}(\sigma_{8}\tau_{17})&\beta_{9}(\sigma_{9}\tau_{18})\\ \delta_{4}\gamma_{4}\alpha_{4}\sigma_{4}&\beta_{5}\alpha_{5}(\gamma_{5}(\sigma_{5}\tau_{14})\delta_{8})&\beta_{6}(\gamma_{6}(\sigma_{6}\tau_{15})\delta_{9})\\ \gamma_{1}\alpha_{1}\sigma_{1}&\beta_{2}\alpha_{2}(\gamma_{2}(\sigma_{2}\tau_{11})\delta_{5})&\beta_{3}(\gamma_{3}(\sigma_{3}\tau_{12})\delta_{6})\end{array}$$

Full text

I Introduction

To understand the electronic properties of realistic materials, it is important to account for the effects of the long-range Coulomb interaction between electrons. However, solving the many-electron Schrödinger equation (SE) including the Coulomb potential $V_{\mathrm{int}}=\sum_{i<j}1/r_{ij}$ is challenging and typically involves uncontrolled approximations. A promising approach where there is a systematic control of accuracy is provided by the density matrix renormalization group (DMRG) algorithm1, 2 and its higher-dimensional extensions such as projected entangled-pair states (PEPS) 3, 4, 5, which reduce the effective dimensionality of the SE by exploiting the locality typically found in physical systems. In fact, the DMRG has already been widely applied to solve the SE for molecules6, 7, 8 within a finite basis expansion. However, to exploit the power of higher dimensional tensor network states (TNS)9, the standard orbital (or spectral) basis10 expansion for the SE is not ideal. This is because whereas the variational freedom in the TNS parametrization scales only linearly with system size $A$, the representation of the Coulomb operator has a large number of terms that scales like $O(A^{4})$.

In our earlier work11, we proposed to combine higher-dimensional TNS with a real-space lattice discretization of the SE, which reformulates the SE as an extended Hubbard model with density-density type long-range interactions $\hat{V}_{\mathrm{int}}=\sum_{i<j}v_{ij}^{ee}n_{i}n_{j}$, where $i,j$ label lattice sites, $v_{ij}^{ee}=1/r_{ij}$, $r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|=|\vec{i}-\vec{j}|l$, $l$ is the lattice spacing, and $n_{i}$ is the number operator. In this form, the discretization error can be systematically controlled by reducing the lattice spacing $l$. The representation of the Hamiltonian is also improved for TNS simulations, as there are now $O(A^{2})$ interaction terms between the electron sites. Nonetheless, even with this reduction, a term-by-term evaluation of the Coulomb interaction energy $\langle\Psi|\hat{V}_{\mathrm{int}}|\Psi\rangle$ (in which the values of the potential are explicitly computed for each term) still leads to an undesirable quadratic computational scaling with system size.

In the one dimensional (1D) case, the above problem can be overcome by fitting the Coulomb interaction $1/r_{ij}$ to a sum of exponentials $\sum_{t=1}^{N_{t}}c_{t}e^{-\lambda_{t}r_{ij}}$12, 13, where the number of terms $N_{t}$ depends only on the target fitting accuracy rather than the system size $A$. Each exponential interaction $\hat{V}=\sum_{i<j}e^{-\lambda r_{ij}}n_{i}n_{j}$ can then be represented by a matrix product operator (MPO)12, 14, 13, 15, 16, 17 with bond dimension $D=3$,

[TABLE]

In this way, computing $\langle\Psi|\hat{V}_{\mathrm{int}}|\Psi\rangle$ with $|\Psi\rangle$ represented by a matrix product state (MPS)17 with bond dimension $D$ scales as $O(N_{t}AD^{3})$, that is, linearly with the system size. In combination with DMRG, this representation has been used to simulate interacting 1D models in the continuum limit18, 19, 20, 21 with controllable accuracy by systematically reducing the spacing $l$ and increasing the bond dimension $D$.

Generalizing such a construction for long-range interactions to higher dimensions is, however, nontrivial. In 2D or 3D, the exponential $e^{-\lambda r_{ij}}$ cannot be formed as a product of weight factors $e^{-\lambda l}$ along the path from $i$ to $j$, as it is done in 1D. In our previous work11 on the 2D case, we used the spin-spin correlation functions of the 2D classical Ising model $\langle\sigma_{i}\sigma_{j}\rangle_{\beta_{t}}$ (where $\sigma_{i},\sigma_{j}\in\{+1,-1\}$ and $\beta_{t}$ is the inverse temperature) as a basis to numerically fit the Coulomb interaction on a square lattice to a sum of the correlation functions $\langle\sigma_{i}\sigma_{j}\rangle_{\beta_{t}}$ at different temperatures, viz., $1/r_{ij}\approx\sum_{t=1}^{N_{t}}c_{t}\langle\sigma_{i}\sigma_{j}\rangle_{\beta_{t}}$. However, an important difference between the 1D and 2D cases is that $\langle\sigma_{i}\sigma_{j}\rangle_{\beta_{t}}$ at short lattice distances is not a smooth and radially isotropic function of $r_{ij}$. To control these errors in 2D, we embedded the physical lattice into a larger underlying Ising lattice, for details, see Ref. 11. [NB: In the experimental setting, a related recent proposal uses auxiliary particles on a larger underlying lattice to mediate the Coulomb interaction between physical particles for the purposes of building an analog quantum simulation of Hamiltonians with long-range interactions22.] Then, since the Ising correlation function $\langle\sigma_{i}\sigma_{j}\rangle_{\beta_{t}}$ can be represented by a classical PEPS with $D=2$5, 23, by using the finite automata construction14, 13, 24 to couple operators $n_{i}n_{j}$ with weights $\langle\sigma_{i}\sigma_{j}\rangle_{\beta_{t}}$, the long-range interaction $\sum_{i<j}\langle\sigma_{i}\sigma_{j}\rangle_{\beta_{t}}n_{i}n_{j}$ can be represented by a projected entangled-pair operator (PEPO) with maximal bond dimension $D=2\times 3=6$ or $D=2\times 4=8$, depending on the choice of finite automata rules11 to represent the operator sum $\sum_{i<j}n_{i}n_{j}$. Consequently, the long-range interaction $\hat{V}_{\mathrm{int}}$ can be approximated as a sum of PEPOs with constant bond dimension, which is analogous to the 1D case.

Nonetheless, using a numerically defined basis $\langle\sigma_{i}\sigma_{j}\rangle_{\beta_{t}}$ to expand the interaction in 2D complicates matters significantly as compared to using the analytic basis $e^{-\lambda_{t}r_{ij}}$ in 1D. For example, in Ref. 11 the performance of the fit in various limits could only be assessed numerically, and the analysis was restricted to the 2D square lattice. In this work, we define an analytic framework which contains the exponential basis construction in 1D and provides a natural generalization to lattices in any spatial dimension, although we will focus explicitly only on 2D and 3D. Importantly, this formulation produces a set of long-range basis functions with explicit tensor network (TN) representations with small, constant bond dimensions such that decaying long-range interactions can be approximated by a sum of tensor network operators (TNO) efficiently.

Specifically, in Sec. II, we introduce the framework by defining appropriate basis functions (in any dimension) as the Green’s function of the discretized Helmholtz equation. Then, using Grassmann Gaussian integration, the discretized Green’s functions can be expressed as correlation functions of auxiliary fermionic fields $\langle c_{i}\bar{c}_{j}\rangle$. This allows us to show in what sense the 1D geometry is special: there are fundamental differences between exponentials in 1D and their higher dimensional extensions. In 1D, both the discretization and finite size errors can be removed analytically, resulting in a TN representation of the continuum exponential function, while in 2D and 3D they cannot be removed in a way compatible with a simple TN representation of low bond dimension. In Sec. III, the MPO representation in Eq. (4) is re-derived within the proposed framework. In Sec. IV and V, we show that the discretized Green’s functions in 2D and 3D, respectively, also have very compact TN representations. Most interestingly, the analytic construction in 3D yields a TN representation of correlation functions decaying as $r^{-1}_{ij}$ asymptotically, which generates the (screened) Coulomb potential on the cubic lattice. The subroutines for constructing TN representations and examples for numerical contractions of the resulting TN are available online25. Finally, conclusions are drawn in Sec. VI.

II Fitting basis in any dimension

One way to view the exponential $e^{-\lambda|\mathbf{r}-\mathbf{r}^{\prime}|}$ in 1D is that it is the Green’s function of the Helmholtz equation in free space (up to a constant $\frac{1}{2\lambda}$),

[TABLE]

or in other words, it is the Fourier transform of the momentum-space kernel $(|\mathbf{k}|^{2}+\lambda^{2})^{-1}$. In 2D and 3D, the solution to Eq. (5) is given by the modified Bessel function $\frac{1}{2\pi}K_{0}(\lambda|\mathbf{r}-\mathbf{r}^{\prime}|)$ and the screened Coulomb potential $\frac{e^{-\lambda|\mathbf{r}-\mathbf{r}^{\prime}|}}{4\pi|\mathbf{r}-\mathbf{r}^{\prime}|}$, respectively, both of which decay exponentially at large distance for $\lambda>0$. Therefore, similarly to in the 1D case, if we are able to find the corresponding TN representation of these interactions, we can use them as a basis of functions with which to represent decaying long-range interactions.

To begin, we consider the discretized version of Eq. (5) on a $d$-dimensional cubic lattice with spacing $a$ (see Figures 1(a) and 1(b)),

[TABLE]

where $\mathbf{K}_{d}/a^{2}$ is the discretized version of $(-\nabla^{2})$, which in the simplest case can be represented by the central difference scheme with open-boundary conditions (OBC),

[TABLE]

Here the length of the lattice is $L=(N+1)a$, $\mathbf{K}_{2}=\mathbf{K}_{1}\otimes\mathbf{I}+\mathbf{I}\otimes\mathbf{K}_{1}$, and $\mathbf{K}_{3}=\mathbf{K}_{1}\otimes\mathbf{I}\otimes\mathbf{I}+\mathbf{I}\otimes\mathbf{K}_{1}\otimes\mathbf{I}+\mathbf{I}\otimes\mathbf{I}\otimes\mathbf{K}_{1}$. The matrix $V_{ij}$ is related to the continuum Green’s function by

[TABLE]

where the factor $a^{-d}$ in the first identity comes from scaling 111The factor $a^{-d}$ can be understood by considering the discretization of the Gaussian functional integral for $G(\mathbf{r}_{i},\mathbf{r}_{j})=\frac{\int\mathrm{D}[\phi]e^{-S}\phi(\mathbf{r}_{i})\phi(\mathbf{r}_{j})}{\int\mathrm{D}[\phi]e^{-S}}$ with $S=\frac{1}{2}\int\mathrm{d}^{d}\mathbf{r}\phi(\mathbf{r})(-\nabla^{2}+\lambda^{2})\phi(\mathbf{r})$. On a $d$-dimensional lattice with spacing $a$, $S$ becomes $S=\frac{1}{2}\phi_{i}(K_{ij}/a^{2}+\lambda^{2}\delta_{ij})\phi_{j}a^{d}$ such that by a change of variable $\tilde{\phi}_{i}=\phi_{i}a^{d/2}$, $G(\mathbf{r}_{i},\mathbf{r}_{j})=\frac{\int\prod_{k}\mathrm{d}\phi_{k}e^{-S}\phi_{i}\phi_{j}}{\int\prod_{k}\mathrm{d}\phi_{k}e^{-S}}=\frac{\int\prod_{k}\mathrm{d}\tilde{\phi}_{k}e^{-1/2\tilde{\phi}_{i}(K_{ij}/a^{2}+\lambda^{2}\delta_{ij})\tilde{\phi}_{j}}\tilde{\phi}_{i}\tilde{\phi}_{j}a^{-d}}{\int\prod_{k}\mathrm{d}\tilde{\phi}_{k}e^{-1/2\tilde{\phi}_{i}(K_{ij}/a^{2}+\lambda^{2}\delta_{ij})\tilde{\phi}_{j}}}=a^{-d}[(\mathbf{K}_{d}/a^{2}+\lambda^{2}\mathbf{I})^{-1}]_{ij}$.. In order to analytically relate $V_{ij}$ to a tensor network, we can use Gaussian integration to express $V_{ij}$ as the correlation function of an auxiliary bosonic or fermionic system with nearest neighbor couplings. We will use Grassmann variables such that $V_{ij}$ can be written as

[TABLE]

where a pair of Grassmann variables $\{\bar{c}_{i},c_{i}\}$ is associated with each lattice site. The advantage of using auxiliary fermions instead of bosons is that the resulting tensor network representations for $Z$ and $V_{ij}$ will have finite as opposed to infinite27 bond dimension (vide post).

The special nature of the 1D geometry can now be examined: As depicted in Figure 1(a), the physical lattice (red) with spacing $l$ can be placed on an infinite underlying lattice with spacing $a$ with associated Grassmann variables at each site. The infinite boundary sites can then be integrated out analytically ($N\rightarrow\infty$ while keeping $a$ fixed) to remove boundary effects, and the same can be done for the infinite interior sites ($a\rightarrow 0$ while keeping $l$ fixed) in order to remove the discretization error. We then obtain an effective matrix $\mathbf{K}_{1}^{\prime}$ defined only on the physical lattice that replaces $(\mathbf{K}_{1}+\lambda^{2}a^{2}\mathbf{I})$ in Eq. (15),

[TABLE]

It can be easily verified that the inverse of $\mathbf{K}_{1}^{\prime}$ indeed gives the exponential interaction on the $n$-site physical lattice with spacing $l$, i.e., $(\mathbf{K}_{1}^{\prime})^{-1}_{ij}=e^{-\xi|i-j|}$. However, in 2D and 3D (as shown in Figure 1(b)), integrating out the boundary and interior sites to obtain continuum limit interactions between the physical sites introduces couplings among all the physical sites. This results in a dense matrix $\mathbf{K}_{d}^{\prime}$, which does not lead to a simple exact TN representation with constant bond dimension, because every tensor would require a bond to every other tensor. Therefore, in the following discussion, we will focus on finding a TN representation for the discrete analog $V_{ij}$ (15) of the continuum Green’s function in 2D and 3D. Once this is done, the discretization error and finite size error in representing continuum interactions can be reduced by choosing a suitably large or infinite underlying Grassmann lattice and embedding the physical sites (red) in it (blue) as shown in Figure 1(b) to effectively work with a smaller lattice spacing $a$ in $a^{-d+2}V_{i=\frac{\mathbf{r}_{i}}{a},j=\frac{\mathbf{r}_{j}}{a}}$ (14), which is similar in spirit to our previous work11 where we used an underlying larger Ising lattice to mediate interactions. Interestingly, this construction in 3D yields a TN representation of correlation functions that decays as $r^{-1}_{ij}$ asymptotically for $\lambda=0$.

In order to explicitly represent $V_{ij}$ as a tensor network, we note that the partition function $Z$ introduced in Eq. (15) is similar to that of the Ising model, which is easily written as a TN 5, 23 in any dimension. Here, however, Grassmann variables are used rather than the spins $\sigma_{i}\in\{+1,-1\}$, and Grassmann integration replaces the summation over spins. Similarly to in the TN representation of the Ising $Z$, by factorizing $e^{\beta\sigma_{i}\sigma_{j}}$ into local quantities coupled by a virtual bond11, we can rewrite the nearest neighbor coupling term in Eq. (15) as

[TABLE]

The termination of the series for the exponential in the first equality is due to the nilpotency of Grassmann variables, which is the advantage of using fermionic rather than bosonic Gaussian integration in (15). The decomposition in the second equality can be performed in various ways. For simplicity, we use the following form for the local factors $\alpha_{i}$ and $\beta_{j}$,

[TABLE]

Thus, instead of $D=2$ for the bond dimension of the TN representation of $Z$ for the Ising model, we will have a $D=4$ construction here. Now the partition function $Z$ (15) can be expressed as a product of local “projectors” and terms for each bond between two sites, e.g., in 3D it reads

[TABLE]

where $Q_{k}=e^{-(\mathbf{K}_{d}+\lambda^{2}a^{2}\mathbf{I})_{kk}\bar{c}_{k}c_{k}}=1-(2d+\lambda^{2}a^{2})\bar{c}_{k}c_{k}$ with $d=3$ here, and $B^{x}_{ij}$, $B^{y}_{ij}$, and $B^{z}_{ij}$ represent the decomposed pairs in Eq. (24) in different directions. To distinguish the pairs in different directions, in the following discussion we will use different pairs of Greek letters for different directions even though they denote the same vectors as in Eq. (25): $\alpha,\beta$ for pairs in the $X$ direction; $\gamma,\delta$ for pairs in the $Y$ direction; and $\sigma,\tau$ for pairs in the $Z$ directions, see Figure 1(c). That is, $B_{ij}^{x}=\alpha_{i}\beta_{j}$, $B_{ij}^{y}=\gamma_{i}\delta_{j}$, and $B_{ij}^{z}=\sigma_{i}\tau_{j}$, where the summations over components have been omitted for simplicity.

The partition function for the Ising model can be written in the same form as Eq. (26), and by collecting the local factors belonging to the same site together, $Z$ can be represented as a tensor network, and the same strategy applies for the correlation function $\langle\sigma_{i}\sigma_{j}\rangle$. However, in our case, due to the anti-commuting property of Grassmann variables, additional sign factors will appear in moving variables in Eq. (26) to their respective local site. We note that Eq. (26) for 2D is structurally similar to the fermionic PEPS (fPEPS) 28, 29, and it can be viewed as a “classical” fPEPS without a physical index. In Sec. IV and Sec. V we will show how to express $Z$ and $V_{ij}$ as TN in 2D and 3D, and in particular, how to deal with the sign factors that appear in different dimensions by developing graphical rules similar to that for fPEPS30. However, before this, we will first show how the present construction in 1D reproduces the MPO in Eq. (4).

III Revisiting the 1D MPO representation

An important simplification in 1D is that the product $\prod_{<i,j>}B^{x}_{ij}$ in Eq. (26) is already in the desired form, viz., $\prod_{<i,j>}B^{x}_{ij}=\prod_{i=1}^{N-1}(\alpha_{i}\beta_{i+1})=\alpha_{1}\prod_{i=2}^{N-1}(\beta_{i}\alpha_{i})\beta_{N}$, where the subscripts for components in $\alpha$ and $\beta$ have been omitted for simplicity. Similarly, for correlation functions $\langle c_{i}\bar{c}_{j}\rangle$ ($i<j$), the necessary product $(\prod_{<k,l>}B_{kl}^{x})(c_{i}\bar{c}_{j})$ can be arranged into local products $\alpha_{1}(\beta_{2}\alpha_{2})\cdots(\beta_{i}c_{i}\alpha_{i})\cdots(\beta_{j}\bar{c}_{j}\alpha_{j})\cdots(\beta_{N-1}\alpha_{N-1})\beta_{N}$ without introducing any sign factors, because in moving $c_{i}$ or $\bar{c}_{j}$, the terms $(\alpha_{i}\beta_{i+1})$ that must be passed over correspond to a bond and are always of even parity, i.e., products of an even number of Grassmann variables. Therefore, by defining the following local tensors for $\mathbf{K}_{1}$ (13),

[TABLE]

the correlation function $V_{ij}$ ($i<j$) can be written as

[TABLE]

When coupled with the operators $n_{i}n_{j}$ via the finite automata construction24, 13, 14, this form of $V_{ij}$ gives an MPO representation for $\sum_{i<j}V_{ij}n_{i}n_{j}$ with bond dimension $D=3\times 4=12$, that is, the factor $\hat{W}[i]$ appearing in the analog of Eq. (4) reads

[TABLE]

This MPO construction for 1D is not optimal in the sense that it is very sparse and compressible in view of the sparse structure of $B_{k}$ and $C_{k}$. We illustrate this compression for $\mathbf{K}_{1}^{\prime}$ (23) in order to eventually recover the MPO in Eq. (4).

In this case, we can replace $\bar{c}$ in Eq. (25) by $-k_{c}\bar{c}$ to factorize $e^{-k_{c}(\bar{c}_{i}c_{j}+\bar{c}_{j}c_{i})}$. Then, the local tensors read

[TABLE]

where $(A_{k})_{11}=k_{i}$ ($k_{b}$) for interior (boundary) sites on a lattice with $n$ sites. It can be found that the products $(A_{1}\cdots A_{i-1})_{11}=(\frac{1}{2}\mathrm{csch}\xi e^{\xi})^{i-1}$ ($2\leq i\leq n$), $(A_{j+1}\cdots A_{n})_{11}=(\frac{1}{2}\mathrm{csch}\xi e^{\xi})^{n-j}$ ($1\leq j\leq n-2$), and $Z=(A_{1}\cdots A_{n})_{11}=(\frac{1}{2}\mathrm{csch}\xi e^{\xi})^{n-1}$. Due to the sparsity of the local tensors, we can observe that for $V_{ij}$ (42), only the element $(B_{i})_{12}$ can contribute and $(B_{i})_{31}$ cannot, because at the boundary $(A_{1})_{13}=0$ such that $(A_{1}\ldots A_{i-1})_{13}=0$. Similar observations apply to $C_{j}$. Thus, $V_{ij}$ in Eq. (42) can be rewritten as

[TABLE]

such that the resulting $e^{-\xi(j-i)}$ can be re-factorized into a product of factors $e^{-\xi}$ between $i$ and $j$. Therefore, the bond dimension for representing $V_{ij}$ is reduced from 4 to 1, which, when coupled with the operators $n_{i}n_{j}$, leads to the MPO (4) with $D=3$.

In 2D and 3D, such a simplification is unlikely to be possible, thus our construction will lead to a TNO with bond dimension $D=4D_{O}$. The factor 4 comes from the present TN construction for the correlation functions $V_{ij}$, while $D_{O}$ depends on the way that $V_{ij}$ is coupled with the product $n_{i}n_{j}$ to form $\sum_{i<j}V_{ij}n_{i}n_{j}$. It has been shown that in 2D11, $D_{O}=3$ for the snake MPO construction for $\sum_{i<j}n_{i}n_{j}$, and $D_{O}=4$ using a 2D finite automata construction14, 13, 24. In Sec. V, we will show in 3D, $D_{O}$ can be 3, 4, or 5, depending on whether an explicitly 1D, 2D, or 3D finite automata representation for $\sum_{i<j}n_{i}n_{j}$ is used.

IV 2D formulation

The problem of rewriting Eq. (26) and the correlation functions in 2D as products of local terms is more complicated than in 1D. To avoid immediately delving into algebraic details, we will first present the obtained results in terms of graphical rules, for which Figure 1(c) defines the local tensor configuration and Figure 2 defines the correlation functions. Then a sketch of the derivation of these rules will be given via a simple example.

IV.1 Rules for writing down TN representations

By analogy to Eq. (41) for 1D, the local tensors for 2D are defined in the following way,

[TABLE]

where $o_{k}$ can be one of $\{1,c_{k},\bar{c}_{k}\}$ for $A_{k}$, $B_{k}$, or $C_{k}$, respectively. While performing the integration manually as we did in 1D quickly becomes tedious for the different combinations of subscripts $(d,l,u,r)$, the necessary integrals (63) can be easily evaluated using a simple program25. We can give Eq. (63) a graphical representation as shown in Figure 1(c), where the factors ($\delta,\beta,o,\gamma,\alpha$) appear in a clockwise order starting from $\delta$. With this definition, the 2D partition function $Z$ (26) can be demonstrated to be given by the PEPS $Z=\mathrm{Tr}(\prod_{k}A_{k})$. However, unlike in the 1D case, the correlation function $V_{ij}$ is not simply $\mathrm{Tr}(\prod_{k\neq i,j}A_{k}B_{i}C_{j})$. Due to the anti-commutation of Grassmann variables, some additional sign factors will appear when moving $c_{i}$ or $\bar{c}_{j}$ to its local site. Assuming the site at position $(x,y)$ on an $N$-by-$N$ lattice shown in Figure 2 is indexed by $(y-1)N+x$, one finds that there are additional sign factors such as $(-1)^{p(\gamma_{k})}$ (or equivalently $(-1)^{p(\delta_{k+N})}$) appearing at the position shown by the green dot in Figure 2(a). Here, $p(\gamma_{k})$ is the parity of the bond $\gamma_{k}$ and is [math] if $\gamma_{k}$ contains an even number of Grassmanns and $1$ otherwise. It is then seen that when representing $V_{ij}$ as a Grassmann correlation function some parity tensors,

[TABLE]

need to be inserted on the upward bonds for the sites to the left of each fermionic site (red dot). The final TN representation for $V_{ij}$ is given by Figure 2(a) (omitting the red lines).

From this basic representation, we can derive various equivalent TN representations for $V_{ij}$ by using the parity conserving properties of the tensor $O_{k}$ (63), which means that if one of its virtual bonds is odd in parity, then the total parity of the other virtual bonds must also be, $p(o_{k})+1$ mod 2, otherwise the Grassmann integration vanishes. This property allows us to define a jump move similar to that in fPEPS30. Graphically, we can view the parity tensors as a result of the crossing between a fermionic line (red) connecting sites $i$ and $j$ and the bonds between local tensors, see Figure 2(a). Then, starting from this graph, one is free to deform the fermionic line freely, as long as the necessary parity factors are inserted at the crossings. This degree of freedom can be used to make the fermionic line coincide with the path used in the finite automata construction11 to derive rules for coupling with operators $n_{i}n_{j}$, thus allowing the parity factors to be inserted by the automata construction itself. Figures 2(b) and 2(c) are examples of a simple path and snake path, respectively. This eventually allows for the construction of the PEPO for $\sum_{i<j}V_{ij}n_{i}n_{j}$, as shown in our previous work11.

Finally, it should be noted that Figure 2(c) and Figure 2(d) differ by a minus sign. This is because when moving from (c) to (d), the jump through $C_{j}$ (63) introduces a minus sign as $\bar{c}_{j}$ is odd, such that $(-1)^{p(\delta_{j,d})+p(\beta_{j,l})+p(\gamma_{j,u})+p(\alpha_{j,r})}=-1$ for nonvanishing Grassmann integrations. In summary, we can express both $Z$ and $V_{ij}$ in 2D as PEPS, with the latter requiring additional parity factors on certain virtual bonds given by Figure 2(a).

IV.2 Sketch of the derivations

To illustrate how the above rules for 2D are actually derived, we consider a simple $4\times 4$ example. From Eq. (26), the partition function $Z$ can be rewritten as

[TABLE]

where a more instructive way to write the right hand side is the following 2D representation:

[TABLE]

with the even-parity factors $Q_{k}$ and the indices for components omitted for simplicity. In Eq. (70), one should read from the bottom-left factor $\alpha_{1}$ to the upper-right factor $\beta_{16}$ for $Z$ (65). In this representation, it is clear that in order to move all the factors to their local sites, we only need to move all $\delta_{k}$ one row up in Eq. (70) along the 1D sequence for $Z$. One way we found to be convenient is to move them column-by-column from left to right. That is, we first move $\delta_{13}$, $\delta_{9}$, and $\delta_{5}$ sequentially to the respective upper rows, and then consider moving $\delta_{14}$, $\delta_{10}$, and $\delta_{6}$, etc. We illustrate this explicitly for $\delta_{13}$. When moving this past $\beta_{10}$ in the second column, the factor $(-1)^{p(\delta_{13})p(\beta_{10})}$ appears due to the exchange of $\delta_{13}$ with $\beta_{10}$. Using the fact that the bond pairs (24) are always even, i.e., $p(\gamma_{9}\delta_{13})=1$ and $p(\alpha_{9}\beta_{10})=1$, this factor can be made local $(-1)^{p(\delta_{13})p(\beta_{10})}=(-1)^{p(\gamma_{9})p(\alpha_{9})}$, which can be further cancelled out by a local exchange from $\alpha_{9}\gamma_{9}$ to $\gamma_{9}\alpha_{9}$ in the product (70). This is how the ordering of factors in Figure 1(c) is derived. After exchanging $\delta_{13}$ with $\beta_{10}$, we move $\delta_{13}$ past $\alpha_{10}(\gamma_{10}\delta_{14})$, $\beta_{11}\alpha_{11}(\gamma_{11}\delta_{15})$, $\beta_{12}(\gamma_{12}\delta_{16})$, but these can be regrouped into complete bonds, $(\alpha_{10}\beta_{11})(\gamma_{10}\delta_{14})\ldots$ which are all even, thus no more signs accrue in moving $\delta_{13}$. Once the $\delta$ factors in the first column have been moved to their local sites, these sites are in their final forms as shown in Eq. (65), where the product of factors at each site is even and parity preserving. Thus, when moving the factors in the second column, we can jump over the sites in the first column without incurring any sign factor. By repeating this procedure, we can express $Z$ as a PEPS with local tensors defined in Figure 1(c).

The same process applies to the correlation functions. In this case, taking $\langle c_{7}\bar{c}_{10}\rangle$ as an example, the counterpart of Eq. (70) is

[TABLE]

The task is again to move the $\delta$ factors to their local sites, and we can apply the same procedure for Eq. (75). However, one can see that moving $\delta_{13}$ will involve an additional exchange with $\bar{c}_{10}$, which results in an additional sign factor $(-1)^{p(\delta_{13})}$. A similar situation occurs when moving $\delta_{9}$ and $\delta_{10}$ due to the exchanges with $c_{7}$. These additional sign factors give the rule for parity tensors $(P_{k})_{ab}$ (green dots) in Figure 2(a). Thus, the long-range interaction $V_{ij}=\langle c_{i}\bar{c}_{j}\rangle$ is given by the quotient of the TN diagram in Figure 2(a) and that for $Z$.

V 3D formulation

V.1 Rules for writing down TN representations

Similar to the 2D graphical representation, the local tensors in 3D can be written down according to Figure 1(c), viz.,

[TABLE]

where $o_{k}\in\{1,c_{k},\bar{c}_{k}\}$ for $A_{k}$, $B_{k}$, or $C_{k}$, respectively, and where the Grassmann integral can be conveniently evaluated via the same program25. However, the partition function $Z$ in 3D (26) is not simply given by $\mathrm{Tr}(\prod_{k}A_{k})$ as in 1D and 2D. The correct TN representation is given in Figure 3(a), where a swap tensor (black dot, see also Figure 4),

[TABLE]

needs to be introduced at each crossing between a vertical bond and a horizontal bond, when the 3D network is viewed as a projection onto 2D. The necessity for these swap tensors is explained in Sec. V.2.1. [NB: The special case of an $N\times N\times 2$ 3D network is structurally identical to the network for the overlap $\langle\Psi|\Psi\rangle$ between two fPEPS30.] For correlation functions, the rule for the TN representation can still be summarized by the fermionic line (red) in Figure 3(a). We will discuss the derivation of this rule given by Figure 5 in the next section.

Before closing this section for the 3D rules, we mention that the resulting $N\times N\times N$ 3D network can also be viewed equivalently as a $N^{2}\times N^{2}$ 2D network, see Figure 3(b), which can be readily contracted using standard algorithms for 2D PEPS. This mapping also implies that to construct the 3D TNO for $\sum_{i<j}V_{ij}n_{i}n_{j}$, we can use the same finite automata rules used in 2D11 (either the explicitly 2D rules with $D_{O}=4$ or the 1D snake MPO rules with $D_{O}=3$) to construct the tensor network representation of the operator sum $\sum_{i<j}n_{i}n_{j}$. This can be seen by indexing the physical site at the position $(x,y,z)$ by $(z-1)N^{2}+(y-1)N+x$, such that the relative ordering of physical sites is unchanged when mapped into 2D. In addition to these rules for the operators, one can also use a set of “3D” rules with $D_{O}=5$ to construct the TNO representation of $\sum_{i<j}n_{i}n_{j}$, which explicitly uses the 3D lattice structure, see Appendix. Thus, the final 3D TNO representation for $\sum_{i<j}V_{ij}n_{i}n_{j}$ will have bond dimension $D=4D_{O}$, where $D_{O}$ can be chosen to be 3, 4, or 5.

V.2 Sketch of the derivations

While the above rules for expressing $V_{ij}$ as a TN may look familiar to readers who have previously worked with fPEPS, in this section, for a more general audience, we will give a pedagogical explanation of two of the main ingredients in the derivations: (1) in the TN representation of $Z$ in 3D (26) (Figure 3(a)), how we obtain the order of factors in Eq. (76) for the local tensors, and how the swap tensors arise, (2) in the TN representation of $V_{ij}$, how the rule for the fermionic line (red) is derived.

V.2.1 Partition function, local tensors, and swap tensors

For simplicity, we consider a simple $3\times 3\times 2$ lattice shown in Figure 4. From Eq. (26), the partition function $Z$ can be rewritten as

[TABLE]

where the integrand can be written simply as

[TABLE]

which, similarly to Eq. (70), should be read from bottom-left to upper-right. Now to move all factors to local sites, apart from the need to move $\delta$ up one row as in the 2D case, the $\tau$ factors also need to be moved up one layer, which increases the complexity of finding the TN representation of $Z$ in 3D.

Similarly to in the 2D case, we found the most convenient way to move factors in 3D to be face-by-face from left to right. The $\delta$ factors ($\delta_{4}$, $\delta_{7}$, $\delta_{13}$, $\delta_{16}$) can first be moved to their respective local site in the same way as in 2D, viz.,

[TABLE]

where we have exchanged the $\gamma$ and $\alpha$ factors in the first column to compensate for the introduced sign factors. Next, we move the $\tau$ factors in the order $\tau_{16}$, $\tau_{13}$, and $\tau_{10}$, which is essential for simplifying the manipulations, to the upper layer, leading to

[TABLE]

Note that after moving $\tau_{16}$ to the upper row, the factors in site 7 are complete, such that when moving $\tau_{13}$, no additional sign factors due to the jump over this site need to be considered. Thus, the net sign factors introduced are $(-1)^{p(\tau_{16})p(\beta_{8})}$, $(-1)^{p(\tau_{13})[p(\beta_{5})+p(\beta_{8})]}$, and $(-1)^{p(\tau_{10})[p(\beta_{2})+p(\beta_{5})+p(\beta_{8})]}$, respectively. Again by noting the even parity of bonds, the factors such as $(-1)^{p(\tau_{16})p(\beta_{8})}$ can be made local $(-1)^{p(\tau_{16})p(\beta_{8})}=(-1)^{p(\sigma_{7})p(\alpha_{7})}$, and further absorbed locally by exchanging $\alpha_{7}$ and $\sigma_{7}$ in the local product $\delta_{7}\alpha_{7}\sigma_{7}$. These local exchanges to compensate the local sign factors determine the order of factors in Eq. (76) or equivalently Figure 1(c). The remaining signs that cannot be absorbed are given by $(-1)^{p(\tau_{13})p(\beta_{8})}(-1)^{p(\tau_{10})p(\beta_{5})}(-1)^{p(\tau_{10})p(\beta_{8})}$. These nonlocal terms can be exactly represented/decomposed in terms of swap tensors (black dots) shown in Figure 4. The whole process for moving $\delta$ and $\tau$ factors can be repeated for the other faces/columns such that the final TN representation of $Z$ is given by a 3D network composed of local tensors $A_{i}$ and swap tensors.

V.2.2 Correlation functions and fermionic line

For the correlation functions $\langle c_{i}\bar{c}_{j}\rangle$ in 3D, the introduced parity factors can be found in the same manner following the logic for 2D. Thus, we only describe the basic idea here, assuming $c_{i}\bar{c}_{j}$ (even parity and $i<j$) is first placed on site $i$, which means $\bar{c}_{j}$ needs to be moved to site $j$. One can show that the additional parity factors introduced by the fermionic variables $c_{i}$ ($\bar{c}_{j}$) can be classified into three groups for $c_{i}$ ($\bar{c}_{j}$), as shown in Figures 5(a,b,c), for crossings with different bonds. Specifically, the in-plane parity factors for $c_{i}$ ($\bar{c}_{j}$) in Figure 5(a) are the same as those in 2D, see Figure 2(a), while Figures 5(b,c) are new due to the existence of bonds in the $z$-direction. Summarizing all parities and swaps together leads to Figure 5(d), which can be greatly simplified into a single rule of a fermionic line (red) in Figure 5(f) by moving certain parities upwards using the exchange rule shown in Figure 5(e), viz., $\sum_{w^{\prime}}P_{ww^{\prime}}S^{w^{\prime}x}_{yz}=(-1)^{p(w)}S^{wx}_{yz}=(-1)^{p(z)}S^{wx}_{yz}=\sum_{z^{\prime}}S^{wx}_{yz^{\prime}}P_{z^{\prime}z}$ following from the definitions in (64) and (77). Therefore, the final rule shown in Figure 5(f) for half of the fermionic pair and Figure 3(a) for the whole pair $c_{i}\bar{c}_{j}$ is the same as that for 2D, see Figure 2(a).

VI Conclusions

In this work we have presented an analytic construction of the tensor network representations of the discretized Green’s function $V_{ij}$ of the Helmholtz equation in 2D and 3D using Grassmann Gaussian integration. The resulting TN representation is very compact, with bond dimension $D=4$. Interestingly, in 3D it gives an analytic TN representation of correlation functions decaying as $r_{ij}^{-1}$ asymptotically, which yields the discretized (screened) Coulomb interaction on the simple cubic lattice. The TN representation can be made compatible with the rules for finite automata by properly deforming the associated fermionic lines, such that we can construct a TNO representation for $\sum_{i<j}V_{ij}n_{i}n_{j}$. These interactions with different Helmholtz parameters can be used as basis functions to fit decaying long-range interactions in higher dimensions, as an analog of the exponential fitting procedure used in TN algorithms in 1D.

The resulting TN operators can readily be used in simulations of continuum systems, such as the uniform electron gas (UEG), via the combination of a discretized lattice representation and tensor network algorithms. Another possible direction is the simulation of the effective low-energy sectors of lattice gauge theory31 using higher dimensional TNS. Integrating out gauge degrees of freedom leads to Hamiltonians with non-local or long-range terms. In fact, this is precisely how the Coulomb interaction in Nature arises.

Acknowledgements

This work was supported by the US National Science Foundation via Grant No. 1665333. ZL is supported by the Simons Collaboration on the Many-Electron Problem. MJO is supported by a US National Science Foundation Graduate Research Fellowship under Grant No. DEG-1745301. GKC is a Simons Investigator in Physics.

Appendix: Finite automata rules for coupling with operators in 3D

In this section, we discuss how to construct a 3D PEPO representation for the distance-independent interaction operator,

[TABLE]

with a bond dimension $D_{O}=5$. In other words, we will build $\hat{V}$ as $\hat{V}=\mathrm{Tr}(\prod_{k}\hat{P}_{k})$ with $(\hat{P}_{k})_{L_{k},U_{k},D_{k},R_{k},T_{k},B_{k}}=O_{n_{k}n_{k}^{\prime}}^{[k]}$ being a local operator on site $k$. This PEPO representation can be combined11 with the 3D PEPS representation for correlation functions $V_{ij}$ described before (Sec. V), resulting in a PEPO representation of $\sum_{i<j}V_{ij}\hat{A}_{i}\hat{B}_{j}$ with $D=4D_{O}=20$.

The finite automata (also known as finite state machine) picture14, 13, 24 of a PEPO views each tensor as a node in a graph, and each virtual bond of dimension $D_{O}$ as a directed edge in the graph that can pass $D_{O}$ different signals (or has $D_{O}$ different possible states). By convention we have chosen our directed edges to point in the $+\hat{x}$, $+\hat{y}$, and $+\hat{z}$ directions, where the axes are defined in Figure 6. This allows us to impose an ordering of the sites, where we traverse the $x$ direction most quickly, then $y$, then $z$, starting from the bottom left corner. With this convention, the tensor at position $k$ has its $U_{k}$, $R_{k}$, and $T_{k}$ indices (corresponding to $+\hat{y}$, $+\hat{x}$, and $+\hat{z}$ respectively) pass “outgoing” signals while its $L_{k}$, $D_{k}$ and $B_{k}$ indices (corresponding to $-\hat{x}$, $-\hat{y}$, and $-\hat{z}$) receive “incoming” signals (see the first tensor in Fig. 6(b)). When a tensor has certain specific combinations of incoming and outgoing signals, that tensor’s two physical indices $n_{k}$ and $n^{\prime}_{k}$ (not shown in Fig. 6 for simplicity) encode a local operator $O^{[k]}_{n_{k}n^{\prime}_{k}}$, which is either a physical operator ($\hat{A}$, $\hat{B}$) or the identity operator ($I$). These special combinations of index values precisely correspond to the set of rules which generate a finite state machine (PEPO) which encodes all the terms in the sum (100). In order to avoid encoding any additional unwanted terms, the value of $O^{[k]}_{n_{k}n^{\prime}_{k}}$ is the zero operator $\hat{0}_{k}$ when the states of the six virtual indices do not match any rule which generates the desired machine. In other words, unwanted configurations of the state machine (and thus unwanted configurations of the local operators) are prevented by causing such a configuration to trigger the action of $\hat{0}_{k}$ on at least one site of the machine, rendering the entire term null.

The complete list of the rules that define the full 3D PEPO which generates all pairwise interactions in Eq. (100) with bond dimension $D_{O}=5$ is given in Table 1. The presentation is in the style of Ref. 24. In the following, we will provide an intuitive explanation for the derivation of these rules, assuming some familiarity with the simpler constructions in 1D and 2D14, 24, 11. The present 3D construction can be viewed as a generalization of the 2D rules by incorporating an additional set of rules to include interactions between sites in different layers.

VI.1 Basics

A useful way to reason about the construction of finite state machines is to assign some verbal meaning to each of the $D_{O}$ possible signals that can be passed between the nodes. In the present case we have $D_{O}=5$, meaning that each virtual bond of a tensor can take index values of (0, 1, 2, 3, and 4). The meanings that we assign to these signals are used to describe the different “messages” of information that they pass to the adjacent tensor that the bond is connected to. The “0” signal is the default signal, which generally means that the machine is in its initial state along that signal path and no physical operators have been applied yet. “1” is the “stop” signal which, when received, generally tells a tensor to avoid acting with a physical operator but instead to act with the identity operator. This is used when another tensor along that signal path has applied a physical operator and does not want an interaction to be generated along the direction that it sends the “1” message. “2” is the “start” signal, which is passed along the directed edges starting with the action of $\hat{A}$ on site $i$ and terminating with the action of $\hat{B}$ on site $j$. The path of this signal can be thought of as the “interaction path.”

With these signals, we can encode all the terms in Eq. (100) for which $j$ lies in the $+\hat{x}$/$+\hat{y}$/$+\hat{z}$ direction (or along the edges of this sector, for which $\Delta x$, $\Delta y$, or $\Delta z$ can be zero) with respect to site $i$. The rules which generate these terms are 1-16. Rules 1-7 encode the propagation of the “0”, “1”, and “2” signals along the directed edges in straight lines. Rules 8-13 encode the action of the physical operators, which begin and terminate “1” and “2” signals. Rules 14-16 allow for the “2” signal to “turn” in allowed directions. Specifically, 14 allows a “2” which is travelling in the $+\hat{y}$ direction (and is thus received by the $D_{k}$ index) to turn and propagate along the $+\hat{x}$ direction. Rule 15 encodes a turn from $+\hat{z}$ to $+\hat{y}$ and 16 allows a turn from $+\hat{z}$ to $+\hat{x}$. Note that other turns which do not violate the directions of the edges, such as $+\hat{x}$ to $+\hat{y}$ and $+\hat{y}$ to $+\hat{z}$, are not allowed in order to prevent double counting. This illuminates a more subtle convention that we have chosen: for interactions $\hat{A}_{i}\hat{B}_{j}$ in which sites $i$ and $j$ do not lie along a straight line, the “2” signal first propagates in the $+\hat{z}$ direction, then $+\hat{y}$, then $+\hat{x}$ (when $i$ and $j$ are in the same plane but not along a straight line, one of the directions in this ordering is skipped). Figure 6 provides a characteristic example of a set of tensor configurations which encodes one “basic” interaction term.

VI.2 Remaining terms

There are additional terms in the sum (100) for which site $j$ does not lie in the $+\hat{x}$/$+\hat{y}$/$+\hat{z}$ direction with respect to site $i$. Specifically, there are six additional cases:

$\{-\hat{x},+\hat{y},\Delta\hat{z}=0\}$, 2. 2.

$\{-\hat{x},\Delta\hat{y}=0,+\hat{z}\}$, 3. 3.

$\{-\hat{x},+\hat{y},+\hat{z}\}$, 4. 4.

$\{\Delta\hat{x}=0,-\hat{y},+\hat{z}\}$, 5. 5.

$\{+\hat{x}=0,-\hat{y},+\hat{z}\}$, 6. 6.

$\{-\hat{x},-\hat{y},+\hat{z}\}$.

Since $j$ lies in at least one negative direction with respect to $i$, the “2” signal that starts at site $i$ cannot propagate all the way to $j$ because at some point it will need to go against the direction of a directed edge. To account for these terms, the “3” and “4” signals can be introduced to propagate from site $j$ towards site $i$ along the $+\hat{x}$ and $+\hat{y}$ directions, respectively. To complete the interaction, these new signals can then meet up with the “2” that began propagating from site $i$ towards site $j$ along the $+\hat{x}$/$+\hat{y}$/$+\hat{z}$ directions via the introduction of new state machine rules. Below we will explain case-by-case how this is done.

Case 1: $-\hat{x}\,,+\hat{y}\,,\Delta\hat{z}=0$ direction (see Figure 7) In this case, the “2” signal starts at site $i$ and propagates in the $+\hat{y}$ direction according to some of the basic rules (8 and 2). Since $j$ lies in the $-\hat{x}$ direction, the “3” signal starts at site $j$ (rule 19) and propagates in the $+\hat{x}$ direction (rule 17). These two signals meet at their intersection point, and the interaction is completed by a new type of “turning” tensor given by rule 18.

Case 2: $-\hat{x}\,,\Delta\hat{y}=0\,,+\hat{z}$ direction (see Figure 8) In this case, the “2” signal starts at site $i$ and propagates in the $+\hat{z}$ direction according to basic rules 10 and 4. Since $j$ lies in the $-\hat{x}$ direction, the “3” signal starts at site $j$ (rule 19) and propagates in the $+\hat{x}$ direction (rule 17). These two signals meet at their intersection point, and the interaction is completed by a new type of “turning” tensor given by rule 20.

Case 3: $-\hat{x}\,,+\hat{y}\,,+\hat{z}$ direction (see Figure 9) In this case, the “2” signal starts at site $i$ and first propagates in the $+\hat{z}$ direction (rules 10 and 4). It then “turns” to the $+\hat{y}$ direction (rule 15) and propagates (rule 2). Since $j$ lies in the $-\hat{x}$ direction, the “3” signal starts at site $j$ (rule 19) and propagates in the $+\hat{x}$ direction (rule 17). These two signals meet at their intersection point, and the interaction is completed by the “turning” tensor given in rule 18.

Case 4: $\Delta\hat{x}=0\,,-\hat{y}\,,+\hat{z}$ direction (see Figure 10) In this case, the “2” signal starts at site $i$ and propagates in the $+\hat{z}$ direction according to basic rules 10 and 4. Since $j$ lies in the $-\hat{y}$ direction, the “4” signal starts at site $j$ (rule 25) and propagates in the $+\hat{y}$ direction (rule 22). These two signals meet at their intersection point, and the interaction is completed by a new type of “turning” tensor given by rule 21.

Case 5: $+\hat{x}\,,-\hat{y}\,,+\hat{z}$ direction (see Figure 11) Since our convention is to propagate the “interaction path” first in the $\hat{z}$ direction, then $\hat{y}$, then $\hat{x}$, this case is a bit less intuituve than the preceding ones. In our previous analysis, we have pictured the “4” as originating from site $j$ and propagating in the $+\hat{y}$ direction. However, the present case is more easily understood if we adopt a different (but equivalent) picture in which the “4” signal is a special component of the interaction signal propagating from site $i$ to site $j$ which is allowed to travel in the $-\hat{y}$ direction, against the directed edge.

Using this new picture, we start as usual with the “2” signal originating at site $i$ and propagating in the $+\hat{z}$ direction according to basic rules 10 and 4. Next, the signal “turns” from the $+\hat{z}$ direction to the $-\hat{y}$ direction, becoming a “4” (rule 21, as in the previous case). The “4” then propagates in the $-\hat{y}$ direction (rule 22, as above). Finally, the “interaction signal” must turn and travel in the $+\hat{x}$ direction and end at site $j$. Since the “2” can already go in the $+\hat{x}$ direction and terminate at $j$ according to basic rules 3 and 12, it can be reused instead of introducing additional rules. Thus, the “4” propagating along $-\hat{y}$ “turns” to $+\hat{x}$ and becomes a “2” again according to rule 24, and then basic rules 3 and 12 complete the interaction.

Case 6: $-\hat{x}\,,-\hat{y}\,,+\hat{z}$ direction (see Figure 12) In this final case, we combine the two pictures for the “3” and “4” signals used in previous cases. First, the “2” signal starts at site $i$ and propagates in the $+\hat{z}$ direction according to basic rules 10 and 4. Next, the signal “turns” from the $+\hat{z}$ direction to the $-\hat{y}$ direction, becoming a “4” (rule 21) and then propagates in the $-\hat{y}$ direction (rule 22). Since $j$ lies in the $-\hat{x}$ direction, the “3” signal starts at site $j$ (rule 19) and propagates in the $+\hat{x}$ direction (rule 17). The “3” and “4” then meet at their intersection point, and the interaction is completed by a new type of “turning” tensor given by rule 23.

Final rule: Up to this point, all the rules in Table 1 have been utilized except rule 26. This is a special rule that only applies to the tensor in the top right corner of the top plane of the network, where the finite state machine terminates. This rule is included to disallow the state of the machine where all tensors have virtual index values $(0,0,0,0,0,0)$ and a spurious 1 is added to Eq. (100) so that the final operator is $1+\sum_{i<j}\hat{A}_{i}\hat{B}_{j}$ instead of the target $\sum_{i<j}\hat{A}_{i}\hat{B}_{j}$.