Implementing the sine transform of fermionic modes as a tensor network

TL;DR

This paper develops a tensor network implementation of the discrete sine transform of the first kind (DST-I) for fermionic modes, enabling efficient lattice computations with complex boundary conditions.

Contribution

It introduces a recursive, algebraic, and diagrammatic method to second-quantize the DST-I as a tensor network, extending Ferris' spectral tensor network to non-trivial boundary conditions.

Findings

Network complexity scales as 1.25 n log n without swap gates

Provides a systematic approach for generalizing spectral tensor networks

Enables efficient fermionic lattice computations with complex boundaries

Abstract

Based on the algebraic theory of signal processing, we recursively decompose the discrete sine transform of first kind (DST-I) into small orthogonal block operations. Using a diagrammatic language, we then second-quantize this decomposition to construct a tensor network implementing the DST-I for fermionic modes on a lattice. The complexity of the resulting network is shown to scale as $\frac 54 n \log n$ (not considering swap gates), where $n$ is the number of lattice sites. Our method provides a systematic approach of generalizing Ferris' spectral tensor network for non-trivial boundary conditions.

Figures

Equations

$${\tilde{x}}_a:=\sum _{b=0}^{n-1}{\mathrm{e}}^{-2\pi \mathrm{i}ab/n}{x}_b.$$

$$\tilde{x}={\mathrm{DFT}}_{n}x\text{with}{\mathrm{DFT}}_n^{ab}={\omega }_n^{a\cdot b}.$$

$$\begin{gathered} {\mathrm{DFT}}_{2m}=L_{2m}({\mathrm{DFT}}_m\oplus {\mathrm{DFT}}_m) \\ \times \left({\mathbb{I}}_m\oplus \mathrm{diag}({\omega }_{2m}^0,\dots ,{\omega }_{2m}^{m-1})\right)({\mathrm{DFT}}_2\otimes {\mathbb{I}}_m) \end{gathered}$$

$${\mathbb{I}}_m^{ab}:={\delta }_{ab},L_{2m}:={\bar{L}}_{2m-1}\oplus 1,\text{and}$$

$${\bar{L}}_{2m-1}^{ab}:=\{\begin{array}{ll}1& \text{iff}b\equiv am\mathrm{mod}(2m-1)\\ 0& \text{otherwise}\end{array}$$

$$F:={\mathrm{DFT}}_2=\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)$$

$${\mathrm{DFT}}_m\oplus {\mathrm{DFT}}_m\mathrm{diag}({\omega }_{2m}^0,\dots ,{\omega }_{2m}^{m-1}),$$

$${\operatorname{DFT}}_{8}\;=\begin{gathered}\includegraphics[scale={1.0}]{DFT_8.pdf}\end{gathered}$$

$$\begin{gathered}\includegraphics[scale={1.0}]{single_particle_rules.pdf}\end{gathered}$$

$$\begin{gathered}\includegraphics[scale={1.0}]{shorthand_hline.pdf}\end{gathered}$$

$${\widehat{\mathrm{DFT}}}_n:=\frac{1}{\sqrt{n}}{\mathrm{DFT}}_n$$

$$F\to \hat{F}:=\frac{1}{\sqrt{2}}F.$$

$${DST-I}_n^{ab}$$

$${DST-III}_n^{ab}$$

$${\operatorname{DST-I}}_{2m-1}=\\ \bar{L}_{2m-1}\big{(}{\operatorname{DST-III}}_{m}\oplus{\operatorname{DST-I}}_{m-1}\big{)}B_{2m-1}.$$

$$B_{2m-1}:=\left(\begin{array}{ccc}{\mathbb{I}}_{m-1}& & {\mathbb{J}}_{m-1}\\ & 1& \\ {\mathbb{I}}_{m-1}& & -{\mathbb{J}}_{m-1}\end{array}\right)$$

$${\mathbb{J}}_{m-1}:=\begin{pmatrix}&&1\\ &{\mathchoice{\reflectbox{$\ddots$}\displaystyle{}}{\reflectbox{$\ddots$}\textstyle{}}{\reflectbox{$\ddots$}\scriptstyle{}}{\reflectbox{$\ddots$}\scriptscriptstyle{}}}&\\ 1&&\end{pmatrix}.$$

$$B_{2m-1}=\left(\begin{array}{ccc}{\mathbb{I}}_{m-1}& & {\mathbb{I}}_{m-1}\\ & 1& \\ {\mathbb{I}}_{m-1}& & -{\mathbb{I}}_{m-1}\end{array}\right)\cdot \left(\begin{array}{ccc}{\mathbb{I}}_{m-1}& & \\ & 1& \\ & & {\mathbb{J}}_{m-1}\end{array}\right),$$

$${\operatorname{DST-I}}_{7}=\begin{gathered}\includegraphics[scale={1.0}]{DST_1_7.pdf}\end{gathered}.$$

$${\operatorname{DST-III}}_{2m}=K_{2m}\big{(}{\operatorname{DST-III}}_{m}\oplus{\operatorname{DST-III}}_{m}\big{)}\\ \times\big{(}X_{m}^{-}\oplus X_{m}^{+}\big{)}\big{(}{\operatorname{\overline{DST}-III}}_{2}\otimes\,{\mathbb{I}}_{m}\big{)}\bar{B}_{2m},$$

$${\bar{B}}_{2m}:=\left(\begin{array}{ccc}{\mathbb{I}}_{m-1}& & {\mathbb{J}}_{m-1}\\ & 1& \\ & & {\mathbb{I}}_{m-1}\end{array}\right)\oplus 1,$$

$$A:=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$$

$$C:={\overline{\mathrm{DST}}-III}_2=F\mathrm{diag}(1,\sqrt{2}).$$

$$K_{2m}:=({\mathbb{I}}_2\oplus {\mathbb{J}}_2\oplus {\mathbb{I}}_2\oplus {\mathbb{J}}_2\oplus \dots )L_{2m}.$$

$${\operatorname{DST-III}}_{8}=\begin{gathered}\includegraphics[scale={1.0}]{DST_3_8.pdf}\end{gathered}.$$

$$X_{m}^{\pm}:=\begin{pmatrix}c_{m}^{1}&&&&s_{m}^{\pm(m-1)}&0\\ &\ddots&&{\mathchoice{\reflectbox{$\ddots$}\displaystyle{}}{\reflectbox{$\ddots$}\textstyle{}}{\reflectbox{$\ddots$}\scriptstyle{}}{\reflectbox{$\ddots$}\scriptscriptstyle{}}}&&\\ &&c_{m}^{m/2}+s_{m}^{\pm m/2}&&&\vdots\\ &{\mathchoice{\reflectbox{$\ddots$}\displaystyle{}}{\reflectbox{$\ddots$}\textstyle{}}{\reflectbox{$\ddots$}\scriptstyle{}}{\reflectbox{$\ddots$}\scriptscriptstyle{}}}&&\ddots&&\\ s_{m}^{\pm 1}&&&&c_{m}^{m-1}&0\\ 0&&\cdots&&0&c_{m}^{m}\\ \end{pmatrix}$$

$$c_m^{\ell }:=\cos \frac{\ell \pi }{4m}\text{and}{s}_m^{\ell }:=\sin \frac{\ell \pi }{4m}.$$

$$Y_{\ell ,m}^{\pm }:=\left(\begin{array}{cc}{c}_m^{\ell }& s_m^{\pm (m-\ell )}\\ s_m^{\pm \ell }& c_m^{m-\ell }\end{array}\right).$$

$$y^{\pm }:=c_m^{m/2}+s_m^{\pm m/2}=\sqrt{1\pm \frac{1}{\sqrt{2}}}\text{and}{c}_m^m=\frac{1}{\sqrt{2}},$$

$$X_{8}^{\pm}=\begin{gathered}\includegraphics[scale={1.0}]{X_8.pdf}\end{gathered},$$

$${\operatorname{DST-I}}_{1}=1\quad\text{and}\quad{\operatorname{DST-III}}_{2}=F\operatorname{diag}\Big{(}\frac{1}{\sqrt{2}},1\Big{)}$$

$${\widehat{\mathrm{DST}}-I}_n:={DST-I}_n\sqrt{\frac{2}{n+1}},$$

$$\begin{gathered}\includegraphics[scale={1.0}]{DST_1o_def.pdf}\end{gathered}$$

$${\operatorname{\widehat{DST}-III}}_{n}:={\operatorname{DST-III}}_{n}\sqrt{\frac{2}{n}}\operatorname{diag}\Big{(}1,\dots,1,\frac{1}{\sqrt{2}}\Big{)}$$

$$\begin{gathered}\includegraphics[scale={1.0}]{DST_3o_def.pdf}\end{gathered},$$

$${X^{\prime }}_m^{\pm }:=\mathrm{diag}(1,\dots ,1,\sqrt{2})X_m^{\pm }$$

$$\text{and}{A}^{\prime }:=\mathrm{diag}(1,\sqrt{2})A,$$

$$Y_{\ell ,m}^{\pm }={\hat{R}}_{\pm \ell ,m}{Z}^{\pm },$$

Full text

I Introduction

The study of tensor networks is currently a topic of growing interest both in condensed matter physics and quantum computation. In condensed matter physics, tensor networks can be used to model the entanglement structure of quantum states and are therefore well suited for the study of ground states of strongly correlated systems Eisert (2013); Orús (2014a, b). Specifically, the formulation of the multiscale entanglement renormalization ansatz (MERA) in this framework has shown to be very fruitful Vidal (2008); Evenbly and Vidal (2009); Corboz and Vidal (2009); Orús (2014b) and can be understood as a lattice realization of the holographic principle Susskind (1995); Swingle (2012), contributing to a better understanding of gauge-gravity type dualities Maldacena (1998). In quantum computation, on the other hand, tensor networks are known as quantum circuits Nielsen and Chuang (2000) and provide a natural framework for the decomposition of a large, complicated unitary operation into a sequence of small local building blocks – a simple example being the factorization of the quantum Fourier transform Nielsen and Chuang (2000) into a sequence of Hadamard and phase gates, which causes the efficiency of Shor’s algorithm Shor (1997).

A unification of both subjects—namely, the simulation of strongly correlated quantum systems by a quantum circuit—was already proposed by Feynman in 1982 Feynman (1982) and has recently been realized by a circuit Verstraete et al. (2009) for implementing the exact dynamics of free fermions on a circle. More recently, Ferris refined this idea Ferris (2014), giving it the interpretation of a spectral tensor network, which implements the Fourier transform of fermionic modes, hence diagonalizes the Hamiltonian of free fermions Lieb et al. (1961). Additionally, the geometry of this network generalizes to higher dimensions, while always staying easily contractible, such that it can be used for the classical simulation of phase transitions in more than one dimension – a feature that is absent in the MERA Evenbly and Vidal (2014a, b).

The starting point for the construction of Ferris’ spectral tensor network is a recursive algorithm for the discrete Fourier transformation (DFT), widely known as Fast Fourier transform (FFT) Cooley and Tukey (1965). The network is then chosen such that it implements the FFT in the one-particle sector of the theory Ferris (2014). This means that every particle is transformed separately, subject to the constraint of their respective indistinguishability. Note that this is entirely different from the quantum Fourier transform, which calculates a single DFT on the space spanned by all particles.

A closely related transformation, on which we will focus in the present work, is the discrete sine transformation of first kind (${\operatorname{DST-I}}$). This is a variant of the DFT for vanishing Dirichlet boundary conditions and thus indispensable for diagonalization procedures of systems on an open chain Lieb et al. (1961). However, the boundary conditions break the convenient cyclic translational symmetry vital to the decomposition of the DFT, meaning that the original FFT algorithm can no longer be used.

The algebraic theory of signal processing Püschel and Moura (2008a, b) tackles this problem by constructing an algebraic framework for spectral transformations. Based on these notions, decompositions of whole classes of transformations were derived in Püschel and Moura (2003, 2008c), including the ${\operatorname{DST-I}}$ as a special case. However, in contrast to the simple case of the standard FFT, the resulting decomposition almost entirely consists of non-unitary elementary transformations. While this is acceptable for numerical recipes restricted to a single particle, we run into problems when translating the decomposition into a quantum circuit, implementing the transformation for indistinguishable particles.

In this paper we unitarize a recursive algorithm for the ${\operatorname{DST-I}}$, originally given in Püschel and Moura (2008c), and show how to second-quantize the resulting algorithm, obtaining a spectral tensor network for a fermionic chain with open ends. To this end we reorganize the network in a non-trivial way. To keep the formal ballast as small as possible, the present paper explains most of the steps in terms of diagrams, showing selected parts of the decomposition and an explicit example of a complete algorithm.

II The fast Fourier transform and diagram notation

The conventional discrete Fourier transform (DFT) converts a sequence of $n$ complex numbers $x_{0},x_{1},\ldots,x_{n-1}$ into another sequence of complex numbers $\tilde{x}_{0},\tilde{x}_{1},\ldots,\tilde{x}_{n-1}$ by means of the linear transformation

[TABLE]

Defining the phase factor $\omega_{n}:=\mathrm{e}^{-2\pi\mathrm{i}/n}$, this transformation can be written as

[TABLE]

Here and in the following, we use the convention that the lower index $n$ denotes the dimension of the respective vector space, while upper indices $0\leq a,b<n$ label components of the matrix.

For even $n=2m$, the recursion formula

[TABLE]

Diagrammatically, the recursion relation (1a) for $n=8$ can be represented as

[TABLE]

running from bottom to top. Here, blocks represent matrices, the ingoing lines are columns, while the outgoing lines lines are rows. The composition rules for such diagrams are

[TABLE]

i.e., vertical composition couples rows of the lower block to columns of the upper block and therefore represents ordinary matrix multiplication, while horizontal composition gives rise to the direct sum of matrices 111Note that this is in strict contrast to the usual diagrammatics of tensor networks and stems from the fact that we are not considering multiple particles yet. This will change in section V, where horizontal composition will give rise to the tensor product. We thank the referee for suggesting to emphasize this point.. Since the matrix $F$ is applied to non-neighboring lines, we do not draw it as a box but rather use the shorthand notation

[TABLE]

with two bullets to indicate on which lines the matrix acts. Note that lines can be crossed arbitrarily, allowing us to move the remaining unaffected lines freely.

Now, if $n$ is a power of $2$, we can use eq. 1a, anchored at ${\operatorname{DFT}}_{2}=F$, to implement the ${\operatorname{DFT}}_{n}$ using only $2\times 2$-matrices. In the above diagrammatic language, this amounts to having no blocks act on more than two lines. By construction, the number of blocks in such a diagram then scales as $n\log_{2}n$, which therefore serves as an upper bound for the computational complexity of the DFT.

In order to interpret the DFT as a change of orthonormal bases in a single-particle Hilbert space, it needs to be unitary. This is achieved by the normalization

[TABLE]

and, as can be easily checked, eqs. 1a and 1b remain valid under the replacement ${\operatorname{DFT}}\to{\widehat{\operatorname{DFT}}}$, provided we also normalize $F$ by

[TABLE]

We thus obtain a decomposition of the large unitary $\smash{{\widehat{\operatorname{DFT}}}_{n}}$ into small unitary building blocks, which again carry the interpretation of basis transformations in the single particle Hilbert space.

The process of second quantization Berezin (1966) then maps the $\smash{{\widehat{\operatorname{DFT}}}_{n}}$ to a basis change in the many-particle Hilbert space. Remarkably, the FFT scheme (1b) is still valid in the many-particle case, provided that we slightly change the composition rules (2) and replace the blocks by their respective second quantizations. We will discuss second quantization of Fourier transforms in more detail later in section V.

III Decomposing the discrete sine transform

The discrete sine transform (DST) is a real linear transformation which captures essentially the imaginary part of the DFT. As it expands the data in sinusoids, this transformation is particularly useful for discrete systems with Dirichlet boundary conditions. However, depending on the implementation of the boundary conditions and the respective periodic continuation, one distinguishes various types of DSTs, which are usually labeled by Roman numbers from I to VIII Püschel and Moura (2008b).

We focus on the ${\operatorname{DST-I}}$ here, which corresponds to a periodic continuation that is odd around both $x_{-1}$ and $x_{n}$. Thus, we have $x_{-1}=x_{n}=0$, making the ${\operatorname{DST-I}}$ suitable for systems with open boundary conditions. As we will see below, for a suitable recursion scheme we also need the ${\operatorname{DST-III}}$. The two DSTs are defined by the matrices

[TABLE]

with $0\leq a,b<n$, which are non-orthogonal. An advantage of considering the ${\operatorname{DST-I}}$ is that we only need the ${\operatorname{DST-III}}$ and the ${\operatorname{DST-I}}$ itself in the corresponding recursion. In principle, it is possible to consider other types of sine (or cosine) transforms and orthogonalize them in a similar way as described in the next section, though the corresponding recursions may be more complicated. Specifically, the ${\operatorname{DST-I}}$ of odd size $n=2m-1$ can be expressed recursively in terms of ${\operatorname{DST-I}}$ and ${\operatorname{DST-III}}$ as Püschel and Moura (2008c)

[TABLE]

In the diagrammatic notation established before, the recursion relation (7a) for $n=7$ can be represented as

[TABLE]

The ${\operatorname{DST-III}}$ appearing in this recursion can be further reduced by means of another recursion relation. For the ${\operatorname{DST-III}}$ of size $n=2m$ we use

[TABLE]

The remaining parts to consider are the matrices $X_{m}^{\pm}$. For even size $m$, they are given by

[TABLE]

with

[TABLE]

Clearly, these matrices can also be decomposed into blocks acting only on pairs of components $(\ell-1,m-1-\ell)$, $\ell=1,\dots,m/2-1$ via

[TABLE]

Further, in the center and the right lower corner of the matrix (9), we have the factors

[TABLE]

acting on the components $m/2-1$ and $m-1$. For example, the matrix $X_{m}^{\pm}$ of size $m=8$ can be drawn in our diagrammatic notation as

[TABLE]

where the indices $\ell$ and $m$ of the $2\times 2$ matrices $Y_{\ell,m}^{\pm}$ are given on the left and right of the corresponding symbol, respectively.

Putting all this together, the two mutually dependent recursions relations (7a) and (8a) together with the closing conditions

[TABLE]

allow us to calculate the ${\operatorname{DST-I}}$ of size $n=2^{k}-1$ using only $2\times 2$ matrices. However, in the existing formulation, all the occurring matrices, except for permutations, are still non-orthogonal. This is a problem in the corresponding quantum version, since non-orthogonal building blocks cannot be interpreted as elementary changes of orthonormal bases. Fortunately, it is possible to reformulate the recursion relations in an orthogonal way, as will be shown in the next section.

IV Orthogonalization of the recursion relations

We now show step by step how to obtain an orthogonal version of the recursion relations (7a) and (8a) discussed in the previous section. For convenience, we label all orthogonal matrices by a hat symbol.

As stated before, the DSTs defined by eqs. 5 and 6 are non-orthogonal. However, all DSTs can be made orthogonal by a suitable scaling of rows and columns. Let us start with the ${\operatorname{DST-I}}$. Multiplying the corresponding matrix (5) by its transpose, we find that the correct scaling is given by

[TABLE]

in our graphical notation. Now we can try to recast the recursion (7a) in terms of orthogonal matrices. The procedure is shown in fig. 1 for a ${\operatorname{DST-I}}$ of size ${n=2^{k+1}-1}$.

The diagram on the left is just eq. 7b with the proper scaling factors according to eq. 11 added at the bottom. In the first step, we split up the factors for each component and move factors of $1/\sqrt{2^{k-1}}$ above the matrices $F$. This is possible since we have the same factor at every component. As no matrix $F$ is acting on the component in the middle, we can also move the remaining factor of $1/\sqrt{2}$ for this component further up, as indicated by the arrow. Now, all we have to do is to absorb all factors into the matrices directly above them. We obtain orthogonal matrices $\hat{F}$ as defined in eq. 4. Further, the factors absorbed into the ${\operatorname{DST-I}}$ are just the ones from eq. 11, so we obtain an orthogonal version of this transform. Since we know that the whole transform in fig. 1 is orthogonal and all other building blocks are orthogonal, too, we can conclude that the ${\operatorname{DST-III}}$ together with the factors in the dotted box must also be orthogonal. Defining the orthogonalized ${\operatorname{DST-III}}$ as

[TABLE]

we obtain the right diagram of fig. 1, where all occurring matrices are orthogonal.

Knowing the orthogonal version of the ${\operatorname{DST-III}}$, we now turn to the corresponding recursion relation (8a). In a first step, we replace the ${\operatorname{DST-III}}$ by its orthogonal version, as it is shown in fig. 2 for a ${\operatorname{DST-III}}$ of size $n=2^{k+1}$. We start with the recursion relation from eq. 7b with scaling factors according to eq. 12 added at the bottom of the diagram. In the first step, we use eq. 12 to replace the two smaller DSTs by their orthogonal versions and the corresponding inverse scaling factors. Further, we express $C=F\operatorname{diag}(1,\sqrt{2})$ in terms of $\hat{F}$ using eq. 4. All factors except the ones in the dotted boxes cancel out. In the second step, we absorb those into the matrices below, obtaining the matrices

[TABLE]

which unfortunately are still non-orthogonal. The part of the recursion relation that remains to be orthogonalized is indicated by a dashed box in the right diagram of fig. 2.

Let us have a look at ${X^{\prime}}_{\!\!m}^{\pm}$ first. The diagonal matrix in eq. 13 just cancels the factor $c_{m}^{m}=1/\sqrt{2}$ in the lower right corner of $X_{m}^{\pm}$. Further, the occurring matrices $Y_{\ell,m}^{\pm}$, defined in eq. 10, can be decomposed as

[TABLE]

with the rotation matrix

[TABLE]

and the non-orthogonal matrix

[TABLE]

Thus, the diagrammatic representation of ${X^{\prime}}_{\!\!m}^{\pm}$, drawn below for $m=8$, is

[TABLE]

Using this representation, the operations in the dashed box in the right diagram of fig. 2 can be redrawn as shown in the left diagram of fig. 3. We have doubled the number of components in this diagram to show all relevant structures and used the abbreviation $m=2^{k}$. Again, the non-orthogonal part is highlighted by a dotted box.

In order to obtain the expression made up from orthogonal matrices given by the right diagram in fig. 3, we observe that the part in the dotted box can be decomposed into three sorts of blocks for any size $n=2^{k+1}$. On the pair of components $(2^{k}-1,2^{k+1}-1)$, we have a trivial block

[TABLE]

which is already orthogonal. Further, the pair of components $(2^{k-1}-1,2^{k}+2^{k-1}-1)$ is coupled by

[TABLE]

where the reformulation on the right hand side only contains matrices which are orthogonal. All other operations decompose into blocks acting on four components

[TABLE]

with $\ell=0,\dots,2^{k-1}-2$. These blocks may be orthogonalized by the relation

[TABLE]

where we introduced the matrix

[TABLE]

Replacing the operations in the dashed box in the right diagram of fig. 2 by the right hand side of fig. 3 (in the appropriate size), we obtain a recursion relation for the ${\operatorname{DST-III}}$ that only contains orthogonal operations.

This completes all steps that are required to obtain a completely orthogonal recursion relation for the ${\operatorname{DST-I}}$ of size $n=2m-1=2^{k+1}-1$. Let us summarize our final results: For the ${\operatorname{DST-I}}_{n}$, we arrive at a binary recursion, which can be represented diagrammatically, e.g. for ${n=7}$, as

[TABLE]

Here, we defined the matrix

[TABLE]

which is just the orthogonalized version of the matrix $B_{2m-1}$ from the non-orthogonal relation (7a).

For the ${\operatorname{DST-III}}_{n}$, we found the diagram

[TABLE]

The newly introduced matrix $\hat{Q}_{m}^{\pm}$ is acting on pairs of components $(\ell-1,m-1-\ell)$, $\ell=1,\dots,m/2-1$ via $\hat{R}_{\pm\ell,m}$ and therefore has a similar structure as $X_{m}^{\pm}$ defined in eq. 9. The other new matrix $\hat{N}_{2m}$ couples the pair of components $(m/2-1,3m/2-1)$ by $\hat{R}_{-m,2m}$ and the pairs of components $(m/2-1+\ell,3m/2-1-\ell)$, $\ell=1,\dots,m/2-1$ via $\hat{G}$.

The corresponding closing conditions for eqs. 14 and 15 are

[TABLE]

such that the complete recursion leads to a well-defined network of operations. As an example, we drew the complete network for the $\smash{{\operatorname{\widehat{DST}-I}}_{31}}$ in fig. 4.

To calculate the computational complexity of the derived algorithm, denote by $\mathcal{C}^{\mathrm{I}}_{n}$ and $\mathcal{C}^{\mathrm{III}}_{n}$ the number of elementary operations in the ${\operatorname{\widehat{DST}-I}}_{n}$ and ${\operatorname{\widehat{DST}-III}}_{n}$, respectively, neglecting permutations. From eq. 14, we find

[TABLE]

where the last summand is the number of $\hat{F}$ matrices in each recursion step. On the other hand, eq. 15 implies

[TABLE]

the additional summands being the number of $\hat{R}_{\ell,m}$, $\hat{F}$, and $\hat{G}$ matrices per recursion step, in that order. Evaluating these formulae, anchored at $\mathcal{C}^{\mathrm{I}}_{1}=0$ and $\mathcal{C}^{\mathrm{III}}_{2}=1$, we obtain

[TABLE]

V Second quantization of diagrams

We will now outline a general method for performing second quantization of diagrammatic algorithms. We shall discuss fermions only – all results, however, extend naturally to the bosonic case Ferris (2014); Dereziński and Gérard (2013).

Let us start with a basis transformation $U$ on the single-particle Hilbert space. Its matrix representation is given by

[TABLE]

We can extend $U$ to a basis transformation $\Gamma_{\!U}$ of the multi-particle Hilbert space, by having it leave the vacuum $\left|\Omega\right\rangle$ invariant and act linearly on creation operators $c^{\dagger}$ Bogolyubov (1958); Dereziński and Gérard (2013). This means that in the occupation number basis

[TABLE]

we define the second quantization $\Gamma_{\!U}$ of $U$ by

[TABLE]

where the transformed mode

[TABLE]

creates a fermion in the state $U\left|a\right\rangle$.

Similarly to eq. 18, we then have

[TABLE]

where eqs. 19 and 20 can be used to obtain

[TABLE]

Denote now by

[TABLE]

the (ordered) lists of occupied modes in $\left|l\right\rangle$ and $\left|k\right\rangle$, respectively, where $p_{l}$ and $p_{k}$ are their total numbers. Obviously, we have $\Gamma_{\!U}^{lk}=0$ for $p_{l}\neq p_{k}$, since then the modes in eq. 21 do not match up in pairs. Let us thus consider the case where $p_{k}=p_{l}=:p$. Here, we have

[TABLE]

where the expectation value on the right hand side again vanishes if $(b_{0},\ldots,b_{p-1})\neq L$ as sets. Since the summand also changes sign under odd permutations of $(b_{0},\ldots,b_{p-1})$, we obtain

[TABLE]

which shall serve as a recipe for the calculation of $\Gamma_{\!U}$.

The Slater determinant expression (22) enables us to check the functorial relations Dereziński and Gérard (2013)

[TABLE]

Also, we see that second quantization preserves unitarity and orthogonality, since

[TABLE]

We can now use eq. 23 to second-quantize the diagrams from the preceding sections. This amounts to replacing all single particle matrices $U$ by their respective second quantizations $\Gamma_{U}$. Each vertical line then represents an occupation number, hence the vertical composition of blocks as on the left hand side of eq. 2 turns into a tensor contraction

[TABLE]

while the right hand side turns into the tensor product

[TABLE]

Note that this does not affect the structure of the diagrams at all, but just amounts to a reinterpretation of what they are supposed to mean.

Since only scalars $\alpha$ and $2\times 2$-matrices $U$ appear in the diagrams we use, we can explicitly evaluate eq. 22 to obtain

[TABLE]

in the Kronecker basis ${\left\{\left|00\right\rangle,\left|01\right\rangle,\left|10\right\rangle,\left|11\right\rangle\right\}}$.

Finally, we have to give a meaning to the crossings of vertical lines, as in eq. 3. For a single particle, these are represented by the swap matrix

[TABLE]

hence in the multi-particle case, we can apply eq. 25 to find

[TABLE]

which correctly picks up a negative sign if two fermions switch places.

As an example, consider the normalized version of the FFT diagram in eq. 1b. We can recursively break it down, so that we only need the second quantization of $\hat{F}$ from eq. 4, given by the matrix

[TABLE]

and local terms as in eq. 24. The entire process can be found in fig. 5, resulting in a diagram that exactly reproduces the spectral tensor network by Ferris (2014); Verstraete et al. (2009). Correspondingly, we can use the same rules on the unitary decompositions (14) and (15) derived in the preceding section, to obtain a quantum circuit implementing the ${\operatorname{DST-I}}$ for fermions and thus generalizing the spectral tensor network for open boundary conditions. Therefore, eq. 17 gives an upper bound for the quantum computational complexity Nielsen and Chuang (2000) of the ${\operatorname{DST-I}}$ and ${\operatorname{DST-III}}$ of fermionic modes. Note, however, that we omitted any permutations in the calculation leading to eq. 17. While this is fine if we deal with just a single particle, exchange statistics need to be incorporated for multiple fermions, hence permutations need to be decomposed into gates of type (26). This leads to additional $\sim\frac{7}{6}n^{2}$ gates, where the coefficient arises from the most naive decomposition and can likely be dropped to a smaller value.

VI Discussion

The present work is based on previous work by Püschel and Moura, who gave a purely algebraic framework for the study of spectral transformations for various types of boundary conditions Püschel and Moura (2008a, b) and their recursive decomposition into sparse matrices Püschel and Moura (2003, 2008c). Although this work is remarkable and very general, its practical use for quantum mechanics is limited in so far as the building blocks of the resulting decomposition are not unitary – a crucial property of any transformation which is to be interpreted as a change of orthonormal bases in a Hilbert space. Since unitary variants of these Fourier transforms exist, it is of course near at hand to expect that also the corresponding recursion relations can be reformulated in terms of unitary building blocks, but carrying out such a unitarization could be quite cumbersome.

In this paper, we have used a diagrammatic language to explicitly demonstrate said unitarization at the example of a discrete sine transformation. This led to a decomposition where the sparse matrices are direct sums of unitary elementary operations, hence well suited for parallelization. We expect that other generalized Fourier transformations can be made unitary in a similar way, although then the technicalities are probably even more involved.

The fact that the resulting recursion relations consist of block-diagonal unitaries is particularly important when turning to many particles in the context of second quantization. We have shown that such a second-quantized version of a sine transformation for fermions can be obtained by replacing the unitary building blocks of the diagram by appropriate second-quantized counterparts and suitable modifications for line crossings, to implement the particle statistics Orús (2014b). Doing so, we arrive at a tensor network, whose structure is essentially the same as that of the original diagram, a circumstance that has already been noted in the context of the ordinary DFT on a circle by Ferris Ferris (2014).

Another great advantage of unitary building blocks becomes apparent when calculating local expectation values with the resulting tensor network: as in the case of the MERA, a causal structure emerges Vidal (2008) and parts of the network that are not causally connected to the considered region can be contracted to trivial operations. Therefore, the effective complexity of the computation of one- and two-point functions drops even further, scaling just linearly with the system size.

The network thus makes it possible to numerically study boundary effects in one dimensional free fermion models, which are exactly solvable by means of a spectral transformation. The Jordan-Wigner transformation extends the applicability to spin-$\frac{1}{2}$ models, such as the XY model Lieb et al. (1961). Imposing the variational method on the tensor coefficients, while fixing the overall topology, the proposed network could also provide a starting point for the simulation of a wider class of models with non-cyclic boundary conditions. Furthermore, since Ferris’ similar construction for the DFT naturally generalizes to higher dimensions, we expect that this also holds for the ${\operatorname{DST-I}}$.

Finally, the observation that second quantization preserves the structure of diagrams, which can be seen as some kind of Bohr correspondence principle on a higher level, seems to be very fundamental and is linked to the underlying category theory, as we shall discuss in a forthcoming paper.