Fourier transform of fermionic systems and the spectral tensor network

Leveraging the decomposability of the fast Fourier transform, I propose a new class of tensor network that is efficiently contractible and able to represent many-body systems with local entanglement that is greater than the area law. Translationally invariant systems of free fermions in arbitrary dimensions as well as 1D systems solved by the Jordan-Wigner transformation are shown to be exactly represented in this class. Further, it is proposed that these tensor networks be used as generic structures to variationally describe more complicated systems, such as interacting fermions. This class shares some similarities with Evenbly&Vidal's branching MERA, but with some important differences and greatly reduced computational demands.

Leveraging the decomposability of the fast Fourier transform, I propose a new class of tensor network that is efficiently contractible and able to represent many-body systems with local entanglement that is greater than the area law. Translationally invariant systems of free fermions in arbitrary dimensions as well as 1D systems solved by the Jordan-Wigner transformation are shown to be exactly represented in this class. Further, it is proposed that these tensor networks be used as generic structures to variationally describe more complicated systems, such as interacting fermions. This class shares some similarities with Evenbly & Vidal's branching MERA, but with some important differences and greatly reduced computational demands. Solving strongly-correlated many-body systems is one of the primary challenges of modern quantum mechanics that arises in a wide variety of fields and contexts. Because direct approaches are usually intractable, an array of analytical and numerical techniques and approximations have been developed over the last century, each associated with its own region of validity, accuracy and difficulty.
One relatively modern approach is based on tensor network decompositions of many-body wavefunctions [1][2][3][4], motivated by our understanding of quantum information in this context. By decomposing the wavefunction into a particular product of tensors, one defines the structure of entanglement that is permitted (for example, the low-energy states of many local Hamiltonians obey an area law scaling for entanglement entropy [5][6][7][8]). Typically, one applies the variational principle within this class of wavefunctions to determine, for example, the tensor network state closest to the ground state. Increasing the rank of the decomposition increases the amount of allowable entanglement, and potentially leads to a more accurate description of the many-body state.
A wide variety of tensor network geometries have been proposed, beginning with the matrix-product state (MPS) [9][10][11] that is the foundation of the famous density matrix renormalization group [12] and Wilson's numerical renormalisation group [13]. This one-dimensional structure connects nearest neighbors in a chain, and is extremely accurate for 1D systems that have spatially-local correlations (proven to be exponentially close to the low-energy states of gapped, local 1D Hamiltonians [5]). Real-space renormalization group ideas lead to the development of the tree-tensor network (TTN) [14,15] and the multi-scale entanglement renormalization ansatz (MERA) [16,17]. The MERA in particular is able to reproduce larger amounts of entanglement such as that necessary to describe critical, 1D quantum systems -and has been found to correctly reproduce the features of the corresponding universality class at the renormalization fixed-point, as described by conformal field theory [18].
Extensions to higher dimensional systems are also possible. Although the MPS and TTN do not have the correct geometry to represent the area law for large 2D systems [15,19,20], the MPS can be generalized to to a two-dimensional tensor network state commonly referred to as a projected entangledpair state (PEPS) [21][22][23]. The MERA is also generalized to higher dimensions in a straightforward fashion [24]. All the above techniques have shown success in solving 2D systems (e.g. [25][26][27][28]), but are somewhat handicapped by relatively large computational demands for even slightly-entangled systems.
Furthermore, both PEPS and the 2D MERA are limited to obey the area-law for entanglement entropy, where the entanglement entropy of a region is proportional to its boundary area (or perimeter, in 2D). This excludes important classes of physical systems which violate the area law -notably any 2D system exhibiting a 1D Fermi surface (resulting in a logarithmic violation [29]), and exotic phases such as Bose metals [30][31][32].
This situation left the field of tensor networks in an awkward position, unable to accurately deal with what many outside the field consider to be the simplest models of noninteracting free-fermions, and their perturbations (i.e. Fermi liquids). More recently, a new class of tensor network that exhibits very large amounts of entanglement (up to volumetric) was proposed by Evenbly & Vidal, which they have dubbed the 'branching MERA' [33]. However, its properties and relationship to Fermi liquids has not been determined, as far as the author is aware.
In this work, I propose a new tensor-network geometry that mimics the structure of a spectral transformation and includes the quantum Fourier transform (for fermions) as a trivial example. Up to a reordering of the lattice, this 'spectral tensor network' can be seen to be a generalization of TTN, or equally a simplification of branching MERA, and local expectation values can be extracted from a wavefunction of n lattice sites in time n log n. Similarly, it is related to previously suggested circuits for use in linear-optics quantum computing [34] and its fermionic analogue [35]. Here, I will give exact wavefunctions for translationally-invariant free fermion systems in arbitrary dimensions as well as examples of 1D systems that are exactly solvable via Jordan-Wigner transformations (Heisenberg XX chains, transverse-field Ising model, etc). I will then describe a procedure for finding a variational approximation of ground state of arbitrary local Hamiltonians in an efficient manner, which may be particularly relevant to simulations of 2D systems in the Fermi liquid phase, before concluding with remarks on future directions.
Quantum fast Fourier transforms: We now discuss the unitary circuit required to transform a quantum lattice system between real-and momentum-space, and write it as a log-depth circuit using at most two-body gates. The discrete Fourier transform has long been known to decompose into a series of sparse operations [36,37], in a technique generally referred to as the fast Fourier transform (FFT). Here we will leverage this decomposition to describe the unitary circuit for the quantum Fourier transform for fermionic lattice models (similar to [35]), and in the next section show that local expectation values of wavefunctions based on this circuit can be efficiently determined.
Following the spirit of Danielson and Lanczos [36], we first demonstrate that the fermionic Fourier transform over n (even) sites can be decomposed into two parallel Fourier transforms over n/2 sites followed by simple 2-body linear gates (matchgates): (1) The fermionic creation operator for site x is denoted byĉ † x , and the LHS of the above is (proportional to) the creation operator for the momentum mode k (with wave-number 2πk/n). The RHS is the sum of two terms, each corresponding to a Fourier-transformed mode on the subset of even (or odd) sites. This sum is a linear transformation between two modes that can be implemented with a 'beam-splitter' and 'phase-delay' (graphically depicted in Fig. 1). The beam-splitter is nothing other than the two-site Fourier transform written above in the number basis {|00 , |01 , |10 , |11 }. Notice thatF 2 |11 = −|11 because of the anti-commutation relations of the creation/annihilation operators. The phase-delay takes care of the so-called twiddle-factor e 2πik/n , and is writtenω k n =Ẑ 2k/n (Ẑ =1 − 2n is the Pauli matrix,n being the occupation operator).
The decomposition in Eq. (1) can be iterated k times to fully determine the Fourier transform over 2 k sites -however one must be careful to keep track of the ordering of modes, given by a bit-reversal operation. The resulting circuit for 8 sites is given in Fig. 1 (a). The structure of the circuit is identical to the Cooley-Tukey FFT [37]. One can independently verify that the circuit produces the correct transformation because it is (a) linear in the field operators, and (b) performs the correct transformation on the single-particle subspace, in analogy with the classical FFT. For systems of higher dimensions, the Fourier transform along each dimension can be applied sequentially along parallel wires and then decomposed. Furthermore, the decomposition can be alternated along different dimensions, for instance to recover an alternating structure with gates along the x, y, etc dimensions at each length scale, as shown for 2D in Fig. 1 (b).
So far, this gives an efficient prescription for performing the Fourier transform on fermionic systems in arbitrary dimensions using a log-depth circuit and two-body gates. This is potentially useful for state preparation or measurement schemes. More generally, a linear circuit with the same structure can be applied to bosonic systems. For bosons,F 2 is the beamsplitter operation that maps two input annihilation operators (â 1 ,â 2 ) to (â 1 +â 2 ,â 1 −â 2 )/ √ 2, while the phase gatê Z a/b transformsâ 1 to (−1) a/bâ 1 . Unlike fermions, the local Hilbert-space dimension is not fixed at two and in general can grow as large as the total number of particles in the system (and thus will not lead to a numerically efficient scheme as discussed below -though this has already been suggested as a powerful tool in linear-optics quantum computation [34]). It would be interesting to determine if the structure is applicable to anyonic or spin models.
The spectral tensor network: Consider a chain of n = 2 k sites described by the second-quantization formalism, that begin in the simple Fock-state |Ψ 0 = |α 1 |α 2 . . . |α n = c †αn n . . .ĉ †α1 1 |0 . The circuit corresponding to the inverse Fourier transform can be applied to an (unentangled) productstate wavefunction in momentum space to yield a (highlyentangled) wavefunction in real space, as depicted in Fig. 2. I call this tensor network structure a 'spectral tensor network', and the wavefunction is a spectral tensor network state. In general, the connectivity of the diagram is that of a hypercube with side length 2 (even in higher spatial dimensions). Ignoring completely the permutation of sites at the bottom of Fig. 2, the state resembles a branching MERA with half of the gates removed (making it significantly cheaper to contract), or similarly a generalized TTN. However, the site reordering turns the traditional concept of real-space renormalization on its head, as the tensor network couples sites with maximal separation first. An intuitive explanation more likely comes from the decomposition of the cyclic group C n corresponding to translation invariance that motivates Eq. (1) and the FFT (where C pq = C p × (C pq /C p ), which is isomorphic to C p × C q ).
When performing calculations in the Fock basis, one must remember to respect the correct statistics of the degrees of freedom in the tensor network diagram by performing the appropriate permutation operation when two wires crossspecifically, introducing a phase of −1 whenever an odd number of fermions cross [38][39][40][41]. This phase can be accounted for by including an extra tensor at every crossing, fortunately without increasing the scaling of the computational cost.
Together, this gives us a prescription for describing translationally-invariant free-particle systems in arbitrary dimensions, because these Hamiltonians are exactly diagonalized by the Fourier transform. In particular, we have a tensornetwork wavefunction for free fermions in arbitrary dimensions. As is well known, such wavefunctions can exhibit quasi-long range correlations (polynomially decaying) and logarithmic violations to the area law. Nonetheless, we now show that extracting local expectation values from these wavefunctions can be achieved by an efficient tensor network contraction scheme, scaling (essentially) linearly with n.
First, we take advantage of the unitary nature of the tensors to vastly simplify the calculation of expectation values of local operators. Figure 3 shows the tensor network corresponding to the expectational value of a single-site operator. The number of tensors reduces from O(n log n) to O(n) and by con- For the systems we study here, both ground states and (grandcanonical) thermal states can be studied because the transformation exactly diagonalizes the Hamiltonian.
As an simple example, we examine free-fermions in 1D and 2D with nearest-neighbour tunnelinĝ where <i, j> corresponds to neighboring sites on the lattice with i < j (or similar in higher dimensions). We create the ground state wavefunction for a low filling factor and calculate the normalized first-and second-order correlation functions: The numerical results for 1D and 2D are shown in Fig. 4. These results should be considered exact, to numerical precision. The 2D simulation describes a highly-entangled system (in real-space) of more than quarter of a million sites, yet the calculation was easily completed on a laptop computer using MATLAB (of course, in this case the results could be obtained analytically or by alternative numerical techniques).
Generic free fermion systems: The quantum circuit structure above can be used to construct translationally invariant states in a very wide range of phases by a straightforward extension to more spin or spatial degrees of freedom, multiple species, and/or including a Bogoliubov transformation to treat anomalous Hamiltonian terms (such asĉ iĉj + h.c.). The first-order and second-order normalized correlation functions g (1) and g (2) for the ground state of a 1D system of free fermions (103 fermions on 1024 sites). (b) Secondorder normalized correlation function g (2) for a 2D system of free fermions (2093 fermions on a 512 × 512 lattice). In both cases the distance (∆x and ∆y) is normalized to the mean particle spacing.
In general, we can consider each site to allow s different types of fermion, described by a Hilbert space of dimension 2 s . Then the Fourier transform from a state initially unentangled in momentum space takes the same form as before, by treating each fermion type individually. The entire state can be encapsulated in a spectral tensor network with bond dimension χ = 2 s . Conversely, the Fourier transform partially diagonalizes translationally invariant free-fermion models, resulting in a Hamiltonian that is local in momentum space (but couples the different species in a momentum-dependent fashion) -meaning the spectral tensor network can describe all eigenstates of translationally invariant bilinear Hamiltonians. With s > 1 we have a richer set of possibilities in a given momentum mode than simply 'occupied' or 'unoccupied', resulting in wider range of possible phases. These include the standard multiple-band possibilities including gapped insulators, standard gapless metals, higher-dimensional systems with interesting Fermi surfaces such as the Dirac point observed in graphene, and perhaps most interestingly, topological systems with non-trivial Chern numbers, as seen in the integer quantum Hall effect. These include as a subset the chiral PEPS states given in Ref. [42].
Hamiltonians with anomalous terms are only partially diagonalized by the Fourier transform, leaving coupling between opposite momentum modes ±k. The Supplemental Material describes how to include the Bogoliubov transformation into the spectral tensor network directly in momentum space while not changing the n log n computational scaling. This addition allows us to treat systems with superconductivity, Majorana fermion modes, or symmetry-protected topological order.
Some spectral tensor network ground states, such as insulators, chiral systems, and 2D systems with Dirac points have entanglement scaling with the area law. Meanwhile, other states such as those in the metallic phase contain significantly more entanglement with a logarithmic violation. It is even possible to achieve a volumetric scaling, for example, with a high-energy eigenstate or the ground state of a long-range Hamiltonian with random couplings. As far as the author is aware, the spectral tensor network is able to exactly de-scribe a wider range of phases and dimensions than any other efficiently-scalable tensor network ansatz.
An additional result is a compact and exact tensor-network representation of strongly-correlated 1D systems that can be analytically solved using a Jordan-Wigner transformation. The Appendix provides details on the process and uses the transverse-field Ising model as an example.
Simulation of arbitrary fermionic systems: More generic systems are free to be studied by employing the variational method on the manifold of states defined by the tensor network diagram (in addition to the unitary constraints and any symmetry constraints of the model, e.g. the fermion parity). For instance, simple energy minimization may yield states very close to highly-entangled ground states, particularly in the Fermi-liquid phase where the entanglement structure should be related to that of free fermions. The bonddimension of the tensor network can be increased by combining neighboring lattice sites, which increases the fraction of Hilbert space covered by the ansatz and thus potentially the accuracy of the minimum energy state. This is an area of intense interest to the author and is the subject of further work. Furthermore, it would be interesting to see if the spectral tensor network can accurately represent systems with different boundary conditions or impurities, or be applied effectively to spin systems.

Conclusion:
The spectral tensor network has been introduced here as a way of efficiently representing the quantum fast Fourier transform and certain classes of quantum states. The Fourier transform can be used to construct free-fermion states with a logarithmic violation of the area law, and as far as the author is aware is the first efficiently contractible tensor network that has been shown to exactly capture this phase in two-and higher-dimensions.
The tensor network is shown to be a convenient form for capturing analytic solutions to 1D spin problems that are solved by the Jordan-Wigner transformation. Further, we suggest that this tensor network may be used as a variational ansatz for low-energy states of many-body systems. This may be especially effective in the case of Fermi-liquid phases, where from a local perspective the system is highly entangled and is somewhat difficult to treat with existing tensor network techniques, while the Fourier transform is able to convert the system into a relatively disentangled momentum-space basis. Work in this direction is ongoing.
The spectral tensor network is somewhat related to the branching MERA construction of Evenbly and Vidal, with two major differences. The spectral tensor network contains only half the tensors of the branching MERA and is more efficient numerically. A more striking difference is the nonlocal reordering of lattice sites placed at the bottom of the unitary circuit. In complete opposition to the spirit of real-space renormalization, the unitaries at the bottom of the circuit deal with the longest length-scales of the system, with the length scales decreasing as we ascend the diagram.
Here we give additional technical details to complement the main text. We begin by giving a fuller derivation of the quantum fast Fourier transform for fermion/bosonic systems, and discuss the inherent freedom allowed by the decomposition. We then detail the tensor contractions necessary for performing calculations with the spectral tensor network, including derivations of the computational cost of measuring one-and two-site local variables and performing ground state optimizations in arbitrary dimensions. We also consider the extended case corresponding to a Fourier transform combined with a Bogoliubov transform that couples modes with opposite momentum, which should be contracted in a slightly different order to remain efficient. The 1D transverse-field Ising model is given as an example application of how both the Jordan-Wigner transformation and Bogliubov transform can be used in conjunction with the spectral tensor network to efficiently extract properties of strongly-correlated systems.

Derivation of the quantum fast Fourier transform
The Fourier transform for fermions (and bosons) over n (assumed even) sites is given by the following linear transformation: We now define two subsystems containing the even sites: and the odd sites Each subsystem contains n/2 sites and have associated Fourier transforms: It is straightforward to verify the following relation (Eq. (1) of the main text), which decouples even and odd sites in the full n-site Fourier transform: (10) This can be written in terms of the above variables as where we note thatã † k+n/2 =ã † k and similarlyb † k+n/2 = b † k . Equation 11 is a simple two-body, linear transformation between modes {ã † k ,b † k } and {c † k ,c † k+n/2 }. Generically this transformation can be represented with the unitary This allows us to decompose the Fourier transform over n sites into smaller gates, which include the two smaller Fourier transforms for n/2 theâ † x andb † x modes (i.e. the even and odd sites), as well as n/2 one-and two-body linear gates. By iterating the decomposition recursively, the remaining Fourier transforms can decomposed into one-or two-body gates. The complete decomposition has depth ∼ log 2 (n) with a total of ∼ n log 2 n gates.
One interesting feature of the fast Fourier transform is the flexibility of in the decomposition, which comes from three sources: (a) graphical manipulations of the unitary circuit/tensor network diagram, (b) the fact that the Fourier transform is self-transpose,F n =F T n , and (c) commutability between different dimensions. Graphical manipulations include 'pushing' the permutation through to the top of the diagram, which simply rearranges the gates and means the permutation performed at the momentum space side of the diagram. The self-transpose property can be applied at different stages of the decomposition. An example using 8 sites is provided in Fig. 5. Furthermore, in two or higher dimensions the transformation layers along different directions commute -see for example Fig. 6. This flexibility is analogous to that available to numerical algorithms for the classical FFT.
As a final note, it should be possible to extend this decomposition to system sizes that are not a power of two while retaining the efficiency. For instance, the decomposition works for factoring out of n any prime factor p, but would necessitate p-body gates and therefore increase the computation scaling as a power of χ. The classical Fourier transform can calculated efficiently (i.e. with n log n scaling) even for large primes n [43], and it would be interesting to determine if a simple quantum circuit can be deduced for these situations, and furthermore whether it would be associated with an efficiently contractible tensor network. However this is beyond FIG. 7: (Color online) The expectation value of a single-body operator can be contracted efficiently by contracting indices on the outside of the diagram and working inwards. Each primitive step consists of contracting a section of the diagram like that highlighted in yellow, corresponding to the operation in Fig. 8 (a). This has to occur n times, leading to a total cost scaling as O(nχ 5 ).
the scope of the present work.

Tensor contractions for one-and two-site operators, and optimization
In the main text, we briefly showed that the expectation value of a single-site observable can be calculated with cost O(n). To reiterate this more precisely, we describe further both the process of removing unitary gates from the expectation value and the contraction scheme required to efficiently compute the final value.
First, we must recall that the tensor network diagrams here are for fermionic Hilbert spaces. Each bond (or wire) in diagrams such as Fig. 3 of the main text represent one or more fermionic modes. The tensor network language lets us graphically express the total Hilbert space as if it were a product of independent Hilbert spaces for each mode (or wire). However, we know such a decomposition is not possible for fermionic Hilbert spaces, at least in the usual sense, and thus additional rules are required to remain consistent with fermionic statistics. Specifically, when two wires cross, this is equivalent to algebraically re-ordering our fermionic modes and thus should introduce a factor of −1 when an odd number of fermions pass through each wire. One way of encoding this is to place a diagonal operator ("mini-tensor") at each of the wire crossings, and then the resulting diagram can be interpreted and computed in the usual (bosonic) sense. Furthermore, it is always possible to choose a new ordering of the fermionic modes that we use to define our Fock-space basis, in either the traditional algebraic or graphical tensor network language. In the diagrams, this is achieved by dragging the wires over each other and pulling them through tensors, each of which is allowed without changing the result (note that mini-tensors at wire crossings may be added or removed because of these manipulations). These graphical rules are a powerful tool for simplifying tensor network diagrams. As an example, consider Fig. 3 of the main text. In this equation, all of the unitaries that are not causally connected to the op-erator at the 7th site cancel. To see this cancellation clearly, we can drag around the wires in the centre (those that cross the "physical" plane) so that the cancelling unitaries multiply their conjugates without wire crossings in-between (at which point they can be removed from the diagram). In some cases (e.g. tensors further from the centre) it may even be necessary to move the unitary tensors up or down to achieve this effect. The freedom given by the allowed manipulations leads to the conclusion that wire crossings do not impair the application of the unitarity identityÛÛ † =Î.
Next we consider the contraction scheme for efficiently computing this expectation value, as depicted in Fig. 7. The contraction proceeds by working from the outside inwards, taking advantage of the very small amount of connectivity in the diagram (consider that the diagram folded around a horizontal line through the center is a simple tree). Each of these basic steps follows the diagrams in Fig. 8 (a) and have cost O(χ 5 ), and occur n times, leading to a total cost of O(χ 5 n).
The expectation value of few-site observables can also be calculated efficiently with cost O(n), although with a cost which grows exponentially with the size of the operator in question. In Fig. 9, the necessary contraction for a widelyseparated, two-site operator is shown. Note here that we are really working in the fermionic Hilbert space and it is unnecessary to include a string of operators between the two observables, as would be typically seen in the Jordan-Wigner transformation (this is a demonstration of the power and simplicity of the graphical, fermionic tensor network language). Like we saw earlier, many of the unitaries cancel each other resulting in a sparser tensor network diagram.
As before, this diagram can be contracted using simple steps starting from the outside of the diagram and working towards the physical indices, layer-by-layer. These steps each cost up to O(χ 8 ) for two-site operators, and are shown in Fig. 8. At the outer-layers, the single-site contractions ( Fig. 8 (a)) suffice until the two sites 'fuse' in the tree, described by Fig. 8 (b). After this point, a effective two-site 'density operators' must be carried through the diagram as in Fig. 8 (c). For a given observable, n such steps are necessary, making the total cost O(χ 8 n) for a two-site observable. The exponent of χ would grow for observables over more sites.
Fortunately, when computing a range of observables many computations can be reused to reduce the total cost. For instance, in a non-translationally invariant model, the entire set of one-body observables Â i can be calculated with cost O(χ 5 n log n). This is achieved by noting that, over all values of i, only two versions of the outside contraction is needed and can be reused (corresponding to the left and right versions of the diagrams in Fig. 8 (a)). However, there are n/2 two-body unitaries in that layer of the circuit, leading to n unique calculations. At the next layer, there are now four possible 'directions' to account for but only n/4 unitaries, so there are again n unique calculations to perform -and so on for every layer. Because there are log 2 n layers, the total number of unique contraction steps required is reduced to just n log 2 n , while a naïve scheme would perform n 2 contraction steps. As another example of this speed-up, in a translationallyinvariant system, the entire set of two-body observables Â 0Bi can be calculated with cost O(χ 8 n log n) (as was done to generate Fig. 5 in the main paper). The same scaling applies to the expectation value of local, two-body Hamiltonians (with no assumption on translational invariance). With a simple modification to the above procedure, one can additionally recover the derivative of the expectation value with respect to the tensors themselves (the so-called 'environment' of each tensor) with the same O(χ 8 n log n) scaling. This allows us to follow a derivative-based optimization procedure for determining the state with minimum energy, such as steepest-descent within the unitary subspace [44], singular-  value decomposition updates [45] or more complicated approaches [46].

Jordan-Wigner transformation for spin chains and Bogoliubov transformations
The spectral tensor network is also useful in treating nonfermionic systems. For example, the Jordan-Wigner transformation can be used to transform some 1D spin systems into non-interacting fermionic models. A straight-forward example is the fully-anisotropic Heisenberg spin-1/2 chain (XX  Fig. 10. Efficient contraction relies on unitarity (first equality) followed by contracting the diagram from the outside in. The typical step for large systems is highlighted in yellow, corresponding to first diagram in Fig. 12. model), which transforms into the free-fermion system studied in the main paper. Other models, such as the quantum transverse-field Ising chain on n sites with Hamiltonian are slightly more complicated because the require a Bogoliubov transformation between opposite momentum modes. This system corresponds to a free fermion system by the transformationĉ i = 1 2 (X i + iŶ i ) j<iẐ j whereĉ i are fermion annihilation operators, resulting in This Hamiltonian is bilinear in the fermionic operators and translationally invariant, and so can be solved using a Fourier transform. In the case of a quadratic Hamiltonian with anomalous pair terms, like above, a Bogoliubov transformation coupling modes with opposite momentum will diagonalize the Hamiltonian. That is, we coupleĉ k withĉ † −k , whereĉ k = j exp(2πijk)ĉ j / √ n, by applying an additional transformation at the 'top' of the tensor network diagram, as shown in Fig. 10.
However, contracting the tensor network naïvely in the previous order will result in growing bond-dimensions in the intermediates steps. To perform the contraction efficiently, we simply perform the layer-by-layer 'renormalization' on the coupled pairs instead of the individual sites, as shown for a single-site observable in Fig. 11. At each step, we combine two pairs of sites with a unitary transform and partial trace, resulting in a density matrix for a single pair of sites. We can see the steps required for calculating single-site observables in Fig. 12 O(n log n) steps for a full set on every site by reusing calculations). A similar procedure applies to the case of observables over two or more sites.
We have implemented and tested this procedure for the transverse-field Ising model. In Fig. 13 we present results for the magnetic susceptibility in the field direction, based on a simple numerical derivative evaluating We observe a sharp peak at the critical point at h = 1, as expected.
The scheme presented here might be improved on in several ways. For instance, it has not yet been investigated if the spectral tensor network can represent the Bogoliubov transformation natively, that is, without the extra tensors at the top linking ±k modes, which is the subject of future work.