Mathematical background
In this page we will analyse the algorithm for simulating the outcomes of a Gaussian boson-sampling experiment in more detail. Since the formulae are already quite heavy by themselves, we will not use the braket notation (except for the outer product \(\outp{u}{v}\) of two vectors \(u\) and \(v\)), opting for a lighter mathematical style.
Gaussian states
We start by laying out some notation and fundamental facts about Gaussian states that we will need later.
Definition
Our context is the Hilbert space of \(n\) bosonic modes. For a single mode the space is the Fock space \(\fockb(\C) \iso \lseq(\C)\), and for \(n\) modes we have \(\fockb(\C)\tsp{n} \iso \fockb(\C^n)\). For \(m\in\N\), we call \(\ns{m}\in\fockb(\C)\) the eigenvector of the number operator with eigenvalue \(m\), and for a multi-index \(\multi{m}=(m\sb1, \dotsc, m\sb{n})\in\N^n\) we define \(\ns{\multi{m}} \in \fockb(\C)\tsp{n}\) as \(\ns{m\sb1} \otimes \dotsb \otimes \ns{m\sb{n}}\). These vectors form an orthonormal basis of the overall Fock space; \(\ns{0}\) is also the vacuum vector, which we will denote by \(\vacuum\).
Let \(R\) be a vector collecting the single-mode position and momentum operators as \(\xpvec\sb{2j-1} \defeq x\sb{j}\) and \(\xpvec\sb{2j} \defeq p\sb{j}\), for \(j\in\{1,\dotsc,n\}\). A Gaussian state \(\rho\) on \(\fockb(\C^n)\) is completely determined by its first and second moments, gathered in the vector \(r\) and in the covariance matrix \(\covmat\) defined as
\[\begin{gather*} \fmom_j \defeq \tr(\rho \xpvec_j),\\ \sigma_{jk} \defeq \tr\bigl(\rho \fcomm{(\xpvec-\fmom)_j, (\xpvec-\fmom)_k}\bigr) = \tr(\rho \fcomm{\xpvec_j,\xpvec_k}) - 2\fmom_j\fmom_k \end{gather*}\]
for \(j,k\in\{1,\dotsc,2n\}\). Covariance matrices are always positive-definite and satisfy \(\covmat+\iu\sympmat \geq 0\), where
\[\sympmat \defeq \imat[n] \otimes \begin{pmatrix} 0 & 1\\ -1 & 0 \end{pmatrix}\]
is the reference symplectic matrix, such that \([\xpvec\sb{j}, \xpvec\sb{k}] = \iu\sympmat\sb{jk}\). We will denote by \(\gauss{\fmom,\covmat}\) the Gaussian state given by its covariance matrix \(\covmat\) and the vector \(\fmom\) of its first moments, or just \(\gauss{\covmat}\) if \(\fmom=0\), as is often the case.
The symplectic group
The symplectic group \(\Sp{2n,\R}\) is the group \(\set{ S \in \Mat{2n,\R} }{ S \sympmat \transpose{S} = \sympmat }\) where \(\sympmat\) is a skew-symmetric \(2n \times 2n\) real matrix such that \(\sympmat^2 = -\imat[2n]\). Every symplectic matrix has determinant \(\pm 1\), is invertible and \(S^{-1} = -\sympmat \transpose{S} \sympmat\).
There are other choices for the skew-symmetric matrix \(\sympmat\) that defines the properties of the symplectic group: another popular choice is
\[\sympmat \defeq \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \otimes \imat[n],\]
with which some properties can be written in a clearer way. We can switch back and forth between these two notations with a permutation matrix \(Y\) such that \(\imat[n] \otimes \begin{psmallmatrix} 0 & 1 \\ -1 & 0 \end{psmallmatrix} = Y \bigl(\begin{psmallmatrix} 0 & 1 \\ -1 & 0 \end{psmallmatrix} \otimes \imat[n]\bigr) \transpose{Y}\).
Let \(\UtoSp \colon \Mat{n,\C} \to \Mat{2n,\R}\) be the map
\[\UtoSp(M) = Y \begin{pmatrix} M\real & -M\imag \\ M\imag & M\real \end{pmatrix} \transpose{Y},\]
where \(M\real\) and \(M\imag\) are the real and imaginary parts of \(M\), respectively. If \(U\) is an \(n\times n\) unitary matrix, then \(\UtoSp(U) \transpose{\UtoSp(U)} = \imat[2n]\) and \(\UtoSp(U) \sympmat \transpose{\UtoSp(U)} = \sympmat\), therefore \(\UtoSp(U) \in \Sp{2n,\R} \cap \Og{2n}\). Vice versa, if \(M \in \Sp{2n,\R} \cap \Og{2n}\) then there exist \(A,B \in \Mat{n,\R}\) such that \(\transpose{Y} M Y = \begin{psmallmatrix} A & -B \\ B & A \end{psmallmatrix}\) and \(A + \iu B \in \Ug{n}\). The map \(\UtoSp\) is thus actually a group isomorphism between \(\Ug{n}\) and \(\Sp{2n,\R} \cap \Og{2n}\).
Moreover, \(M \in \Mat{n, \C}\) is positive-(semi)definite if and only if \(\UtoSp(M)\) is.
Williamson decomposition
Let \(A\) be a positive-definite \(2n \times 2n\) matrix: then there exists \(S \in \Sp{2n,\R}\) and \(D = \bigoplus\sb{j=1}^n d\sb{j} \imat[2]\), with \(d\sb{j} \geq 1\), such that \(A = S D \transpose{S}\) (cfr. Eq. (3.58) in [6]).
Euler decomposition
Let \(S \in \Sp{2n,\R}\): then there exist \(P,Q \in \Sp{2n,\R} \cap \Og{2n}\) and \(d\sb{1},\dotsc,d\sb{n} > 0\) such that \(S = PDQ\) where \(D = \bigoplus\sb{j=1}^{n} \begin{psmallmatrix} d\sb{j} & 0 \\ 0 & d\sb{j}^{-1} \end{psmallmatrix}\).
We can have the outcome be unitary matrices with the aid of the \(\UtoSp\) map introduced before: given a symplectic \(2n \times 2n\) matrix \(S\), we can find \(U,V \in \Ug{n}\) such that \(\UtoSp(U) D \UtoSp(V) = S\).
Unitary representation on a symmetric Fock space
The symplectic group of order \(2n\) admits a unitary representation \(\spunitarysymbol\) on the symmetric (i.e. bosonic) Fock space with \(n\) modes, \(\fockb(\C^n)\). Given \(S \in \Sp{2n, \R}\), its action on the vector \(\xpvec\) of position and momentum operators is given by
\[\begin{equation} \adj{\spunitary{S}} \xpvec_j \spunitary{S} \defeq \sum_{k=1}^{2n} S_{jk} \xpvec_k. \end{equation}\]
such that the canonical commutation relations are preserved, since
\[\begin{equation*} [% \adj{\spunitary{S}} R_j \spunitary{S}, \adj{\spunitary{S}} R_k \spunitary{S} ] = \sum_{l=1}^{2n} \sum_{m=1}^{2n} S_{jl} S_{km} \sympmat_{lm} = \sympmat_{jk}. \end{equation*}\]
On the annihilation and creation operators of the \(j\)-th mode, \(a\sb{j} = \frac{1}{\sqrt{2}}(R\sb{2j-1} + \iu R\sb{2j})\) and \(\adj{a\sb{j}} = \frac{1}{\sqrt{2}}(R\sb{2j-1} - \iu R\sb{2j})\), it behaves as
\[\begin{multline} \adj{\spunitary{S}} a_j \spunitary{S} = \frac12 \sum_{k=1}^{n} \bigl[ (S_{2j-1,2k-1} + S_{2j,2k}) + \iu (S_{2j,2k-1} - S_{2j-1,2k}) \bigr] a_k +\\+ \frac12 \sum_{k=1}^{n} \bigl[ (S_{2j-1,2k-1} - S_{2j,2k}) + \iu (S_{2j,2k-1} + S_{2j-1,2k}) \bigr] \adj{a_k} \end{multline}\]
and similarly for \(\adj{\spunitary{S}} \adj{a\sb{j}} \spunitary{S}\). If \(S\) is also orthogonal, then with \(M = \UtoSp^{-1}(S)\) we can also write
\[\begin{equation} \begin{gathered} \adj{\spunitary{\UtoSp(M)}} a_j \spunitary{\UtoSp(M)} = \sum_{k=1}^{n} M_{jk} a_k,\\ \adj{\spunitary{\UtoSp(M)}} \adj{a_j} \spunitary{\UtoSp(M)} = \sum_{k=1}^{n} \overline{M_{jk}} \adj{a_k}. \end{gathered} \label{eq:fock_space_action_symplectic_orthogonal_matrices_shorter} \end{equation}\]
Normal-mode decomposition
The covariance matrix can always be diagonalised through a symplectic transformation: the Williamson decomposition tells us that there always exist \(S\in\Sp{2n,\R}\) and \(\{\lambda\sb{j}\}\sb{j=1}^{n}\), \(\lambda\sb{j} \geq 1\) such that
\[\begin{equation} \covmat = S \biggl( \bigoplus_{j=1}^{n} \lambda_j\imat[2] \biggr) \transpose{S} \label{eq:covmat-symplectic-diagonalisation} \end{equation}\]
The symplectic eigenvalues \(\lambda\sb{j}\) determine, in turn, the singular values of the whole state, which can be written as
\[\begin{equation} \gauss{\fmom,\covmat} = \displacement{r} \spunitary{S} \Biggl( \bigotimes_{j=1}^{n} \biggl( \sum_{m=0}^{+\infty} \nu_m(\lambda_j) \outp{\ns{m}}{\ns{m}} \biggr) \Biggr) \adj{\spunitary{S}} \adj{\displacement{r}} \label{eq:gaussian-state-normal-mode-decomposition} \end{equation}\]
where \(\spunitary{S}\) is the unitary map on \(\fockb(\C)\tsp{n}\) generated by the symplectic matrix \(S\) appearing in Eq. \eqref{eq:covmat-symplectic-diagonalisation}, \(\displacement{r}\) is the displacement operator by \(r\) and
\[\nu_m(x) \defeq \frac{2}{x+1}\biggl(\frac{x-1}{x+1}\biggr)^m.\]
We can rewrite this by exchanging tensor product and sum, obtaining
\[\begin{equation} \gauss{\fmom,\covmat} = \displacement{r} \spunitary{S} \Biggl( \sum_{\multi{m} \in \N^n} \biggl( \prod_{j=1}^{n} \nu_{m_j}(\lambda_j) \biggr) \outp{\ns{\multi{m}}}{\ns{\multi{m}}} \Biggr) \adj{\spunitary{S}} \adj{\displacement{r}} \label{eq:normal-mode-decomposition} \end{equation}\]
or, by defining \(\hat{\nu}\sb{\multi{m}}(\lambda) \defeq \prod\sb{j=1}^{n} \nu\sb{m\sb{j}}(\lambda\sb{j})\),
\[\begin{equation} \gauss{\fmom,\covmat} = \sum_{\multi{m} \in \N^n} \hat{\nu}_\multi{m}(\lambda) \, \displacement{r} \spunitary{S} \outp{\ns{\multi{m}}}{\ns{\multi{m}}} \adj{\spunitary{S}} \adj{\displacement{r}}. \label{eq:gaussian-state-normal-mode-decomposition2} \end{equation}\]
We can see that the vectors \(\displacement{r} \spunitary{S} \ns{\multi{m}}\) form an orthonormal basis of eigenvectors of \(\gauss{\fmom,\covmat}\), with eigenvalues \(\hat{\nu}\sb{\multi{m}}(\lambda)\). Note that this is both a diagonalisation and a singular value decomposition (not too obvious in general), since all the eigenvalues are positive.
There's no clear way to establish a total order on the eigenvalues: we know \(\nu\sb{m}\) is decreasing in its argument, and always less than one, for each \(m\), moreover \(\nu\sb{m}(x) \leq \nu\sb{m'}(x')\) for any \(x, x' \geq 1\) if \(m' > m\). When we consider \(\hat{\nu}\), however, it's more difficult: we can only say that \(\hat{\nu}\sb{\multi{m}}(\lambda) < \hat{\nu}\sb{\multi{m}'}(\lambda)\) if \(\multi{m} > \multi{m}'\) lexicographically, but still this leaves us with some multiindices we cannot compare.
Construction of the matrix-product state
Suppose we know the covariance matrix \(\covmat\sb{\psi}\) and the vector \(\fmom\sb{\psi}\) of first moments of a pure Gaussian state \(\psi\) over \(M\) bosonic modes. We want to find a matrix-product state representation of \(\psi\),
\[\sum_{m_1=1}^{d} \sum_{m_2=1}^{d} \dotsi \sum_{m_M=1}^{d} A^{(1,m_1)} A^{(2,m_2)} \dotsm A^{(M,m_M)} \, f_{m_1} \otimes f_{m_2} \otimes \dotsb \otimes f_{m_M}\]
in the eigenbasis \(\ns{\multi{m}}\) of the number operator. We will follow the algorithm presented in [1]. The full, “true” form of the state
\[\psi = \sum_{\multi{m}\in\N^M} c_{m_1,\dotsc,m_M} \ns{m_1} \otimes \dotsb \otimes \ns{m_M}\]
cannot usually be known, and we will use what we know from \(\covmat\sb{\psi}\) and \(\fmom\sb{\psi}\) to reconstruct, to a certain degree of approximation, the array \(c\) of coefficients.
We begin by writing the state in the basis made up by the number eigenstates \(\ns{m\sb1}\) on the first mode, and the eigenstates of the reduced state on modes \(2\) to \(M\), that is
\[\begin{equation} \rho\modepart{2\cdots M} = \tr_1 \outp{\psi}{\psi} = \sum_{\alpha_1=1}^{+\infty} (\lambda\modepart{1}_{\alpha_1})^2 \outp{b\modepart{2\cdots M}_{\alpha_1}}{b\modepart{2\cdots M}_{\alpha_1}}. \label{eq:svd-reduced-state-2-to-M} \end{equation}\]
where we assume that the terms in the sum are ordered with decreasing \((\lambda\sb{\alpha\sb1}\modepart{1})^2\). We have
\[\begin{equation} \psi = \sum_{m_1 \in \N} \sum_{\alpha_1=1}^{+\infty} A^{(1,m_1)}_{\alpha_1} f_{m_1} \otimes b\modepart{2\cdots M}_{\alpha_1}. \end{equation}\]
where \(A^{(1,m\sb1)}\sb{\alpha\sb1}\), the first tensor of the MPS, is given by
\[A^{(1,m_1)}_{\alpha_1} = \innp{\ns{m_1} \otimes b\modepart{2\cdots M}_{\alpha_1}}{\psi}.\]
We don't have access either to the full state \(\psi\) or to all eigenvectors \(b\modepart{2\cdots M}\sb{\alpha\sb1}\) from the decomposition in Eq. \eqref{eq:svd-reduced-state-2-to-M}, so we need another way to compute this inner product. Our choice is the formula given in [2].
Ref. [2] gives us a way to compute the inner product
\[\begin{equation} \innp{% \ns{\multi{m'}} }{% \displacement{\alpha} \spunitary{W'} \spunitary{\varLambda} \spunitary{W} \ns{\multi{m}} } \label{eq:franck-condon-formula} \end{equation}\]
where \(\ns{\multi{m'}}\) and \(\ns{\multi{m}}\) are number eigenstates in \(\fockb(\C^n)\), \(W, W' \in \Ug{n}\), \(\displacement{\alpha}\) is a displacement operator, and \(\varLambda = \diag(l\sb1,l\sb1^{-1},\cdots,l\sb{n},l\sb{n}^{-1})\), so that \(\spunitary{\varLambda}\) is a squeezing operator.
We need a few steps to transform this quantity into the desired form. The first step is using the Williamson decomposition of \(\psi\): since it's a pure state, there exists \(S\in\Sp{2M,\R}\) such that \(\covmat\sb\psi=S\transpose{S}\), i.e. \(\psi = \spunitary{S}\vacuum\).
At the same time, we can easily obtain the moments \(\covmat' \defeq \covmat\modepart{2\cdots M}\) and \(r' \defeq \fmom\modepart{2\cdots M}\) of the reduced state by selecting from the original matrix and vector the rows and columns relative to these modes. We compute the Williamson decomposition of the covariance matrix of the reduced state: \(\covmat' = S' Z' \transpose{(S')}\), where \(Z'\) is a diagonal matrix containing the symplectic eigenvalues \(\{z'\sb{j}\}\sb{j=2}^{M}\) (which are greater than or equal to \(1\)), and as in Eq. \eqref{eq:normal-mode-decomposition} we write
\[\begin{equation} \rho\modepart{2\cdots M} = \sum_{\multi{m}'\in\N^{M-1}}\Biggl(\,\prod_{j=2}^{M} \nu_{m'_j}(z'_j)\Biggr) \displacement{r'} \spunitary{S'} \outp{% \ns{\multi{m}'}\modepart{2\cdots M} }{% \ns{\multi{m}'}\modepart{2\cdots M} } \adj{\spunitary{S'}} \adj{\displacement{r'}} \label{eq:normal-modes-reduced-state-2-to-M} \end{equation}\]
Since the state is a positive-semidefinite operator and \((\lambda\modepart{1}\sb{\alpha\sb1})^2\) and \(\prod\sb{j=2}^{M} \nu\sb{m\sb{j}}(z'\sb{j})\) are non-negative, both Eq. \eqref{eq:svd-reduced-state-2-to-M} and Eq. \eqref{eq:normal-modes-reduced-state-2-to-M} are the diagonalisation of \(\rho\modepart{2\cdots M}\) on an orthonormal eigenbasis. The eigenvalue-eigenvector pairs are unique as long as the eigenvalue is not degenerate: therefore, in this case, we have a one-to-one correspondence
\[\begin{equation} \bigl( (\lambda\modepart{1}_{\alpha_1})^2, b\modepart{2\cdots M}_{\alpha_1} \bigr) \Longleftrightarrow \Biggl(\, \prod_{j=2}^{M} \nu_{m'_j}(z'_j), \displacement{r'} \spunitary{S'} \ns{\multi{m}'}\modepart{2\cdots M} \Biggr). \label{eq:svd-reduced-state-correspondence} \end{equation}\]
Here we need to start leaving the infinite Fock space and start truncating, since we would like to write some matrices into the computer memory. There are two truncation points:
- in Eq. \eqref{eq:svd-reduced-state-2-to-M}, we consider a finite amount \(\chi\) of terms, essentially performing a truncated singular-value decomposition of the state;
- in Eq. \eqref{eq:normal-modes-reduced-state-2-to-M}, we limit the multi-indices from \(\N^{M-1}\) to only \(\{0, \dotsc, d\}^{M-1}\), effectively truncating the Fock space, considering at most \(d\) photons on each mode.
These two parameters can potentially be different from mode to mode, but for the sake of simplicity we will consider fixed \(\chi\) and \(d\) for all modes.
Exploiting the correspondence in \eqref{eq:svd-reduced-state-correspondence}, we perform the following operations.
- Combine the symplectic eigenvalues \(z'\sb{j}\) into all possible products \(\prod\sb{j=2}^{M} \nu\sb{m'\sb{j}}(z'\sb{j})\) for \(m'\sb{2}, \dotsc, m'\sb{M} \in \{0, \dotsc, d\}\).
- Take the \(\chi\) largest values, writing down their corresponding multi-indices \(\multi{m}'\sb{\alpha\sb1}\).
- Now we know which are the \(\chi\) “most important” states making up \(\rho\modepart{2\cdots M}\) in Eq. \eqref{eq:svd-reduced-state-2-to-M}: they are \(\displacement{r'} \spunitary{S'} \ns{\multi{m}'\sb{\alpha\sb1}}\modepart{2\cdots M}\) for our chosen \(\multi{m}'\sb{\alpha\sb1}\)s.
These states are also the states we use in the MPS representation of \(\psi\). The state, taking the truncation into account, is then written as
\[\psi = \sum_{m_1=0}^{d} \sum_{\alpha_1=1}^{\chi} A^{(1,m_1)}_{\alpha_1} f_{m_1} \otimes b\modepart{2\cdots M}_{\alpha_1}\]
with
\[\begin{equation} \begin{split} A^{(1,m_1)}_{\alpha_1} &= \innp{% \ns{m_1} \otimes (\displacement{r'} \spunitary{S'} \ns{\multi{m}'_{\alpha_1}}) }{% \spunitary{S} \vacuum } =\\ &= \innp{% \ns{m_1} \otimes \ns{\multi{m}'_{\alpha_1}} }{% \adj{\bigl( \id\modepart{1} \otimes \displacement{r'}\modepart{2\cdots M} \spunitarymp[2\cdots M]{S'} \bigr)} \spunitary{S} \vacuum }. \end{split} \label{eq:first-mps-tensor} \end{equation}\]
as in Eq. (3) in [1].
We know in principle how to generate all elements on the right-hand side of \eqref{eq:svd-reduced-state-correspondence}, but they are infinite, and we can only generate them by trying all the possible multi-indices (up to a certain local dimension) and comparing all the values afterwards. Clearly our strategy allows us to explore only a finite subset of the possible values, and there is no certainty that this subset will actually contain the \(\chi\) largest values we are seeking. As we observed in Normal-mode decomposition, there are some properties that suggest that as the occupation numbers in the multi-indices increase, the eigenvalue decreases, but there is no way to compare the multi-indices apart from actually computing the associated eigenvalue.
This means that there is a small chance of error in this step of the algorithm; however, the error will consist in missing a component of the state whose coefficient is relatively small, so it shouldn't be a serious cause for concern.
From now on, let's ignore displacement operators, since we usually start from \(\fmom=0\). The symplectic matrix \(S'\) on modes \(2\) to \(M\) can be extended to all \(M\) modes as \(\imat[2] \oplus S'\), and in Eq. \eqref{eq:first-mps-tensor} we have
\[\id\modepart{1} \otimes \spunitarymp[2\cdots M]{S'} = \spunitary{\imat[2] \oplus S'}\]
The group representation \(\spunitarysymbol\) is unitary, so we can replace \(\adj{\spunitary{A}}\) by \(\spunitary{A^{-1}}\) and then merge the two \(\spunitarysymbol\)s, obtaining
\[\adj{ \bigl(\id\modepart{1} \otimes \spunitarymp[2\cdots M]{S'}\bigr) } \spunitary{S} = \spunitary{\imat[2] \oplus (S')^{-1}} \spunitary{S} = \spunitary{(\imat[2] \oplus (S')^{-1}) S};\]
Now, we apply the Euler decompostion to the symplectic matrix \((\imat[2] \oplus (S')^{-1}) S\) to find two new symplectic (and orthogonal) matrices \(V\sb1, V\sb2 \in \Sp{2M,\R} \cap \Og{2M}\) and a diagonal matrix \(D = \diag(d\sb1, d\sb1^{-1}, \dotsc, d\sb{M}, d\sb{M}^{-1})\) with \(d\sb{j} \gt 0\) such that \((\imat[2] \oplus (S')^{-1}) S = V\sb1 D V\sb2\). It follows that
\[\begin{equation} A^{(1,m_1)}_{\alpha_1} = \innp{% \ns{m_1} \otimes \ns{\multi{m}_{\alpha_1}} }{% \spunitary{V_1} \spunitary{D} \spunitary{V_2} \vacuum }. \label{eq:first-mps-tensor-simplified} \end{equation}\]
This inner product can be computed following Eq. \eqref{eq:franck-condon-formula}. We have obtained
\[\psi = \sum_{m_1=1}^{d} \sum_{\alpha_1=1}^{\chi} A^{(1,m_1)}_{\alpha_1} \ns{m_1} \otimes b\modepart{2\cdots M}_{\alpha_1}.\]
We continue the construction by further decomposing each \(b\modepart{2\cdots M}\sb{\alpha\sb1}\). The goal is once again to arrive at an expression such as Eq. \eqref{eq:first-mps-tensor-simplified}, which we know how to compute. First, we write the normal-mode decomposition and the diagonalisation of the reduced state over modes \(3\) to \(M\):
\[\begin{align} \rho\modepart{3\cdots M} &= \sum_{\multi{m}\in\N^{M-2}} \Biggl(\,\prod_{j=3}^{M} \nu_{m_j}(z''_j)\Biggr) \spunitary{S''} \outp{% \ns{\multi{m}}\modepart{3\cdots M} }{% \ns{\multi{m}}\modepart{3\cdots M} } \adj{\spunitary{S''}}, \label{eq:normal-modes-reduced-state-3-to-M} \\ \rho\modepart{3\cdots M} &= \sum_{\alpha_2=1}^{+\infty} (\lambda\modepart{2}_{\alpha_2})^2 \outp{b\modepart{3\cdots M}_{\alpha_2}}{b\modepart{3\cdots M}_{\alpha_2}}. \label{eq:svd-reduced-state-3-to-M} \end{align}\]
We express \(b\modepart{2\cdots M}\sb{\alpha\sb1}\) in the eigenbasis of \(\rho\modepart{3\cdots M}\) as
\[b\modepart{2\cdots M}_{\alpha_1} = \sum_{m_2\in\N} \sum_{\alpha_2=1}^{+\infty} A^{(2,m_2)}_{\alpha_1,\alpha_2} \ns{m_2} \otimes b\modepart{3\cdots M}_{\alpha_2},\]
where we have defined
\[A^{(2,m_2)}_{\alpha_1,\alpha_2} \defeq \innp{% b\modepart{2\cdots M}_{\alpha_1} }{% \ns{m_2} \otimes b\modepart{3\cdots M}_{\alpha_2} }.\]
We already know that \(b\modepart{2\cdots M}\sb{\alpha\sb1}\) can be written as \(\spunitarymp[2\cdots M]{S'} \ns{\multi{m}\sb{\alpha\sb1}}\). For the \(b\modepart{3\cdots M}\sb{\alpha\sb2}\) vectors, always assuming there is no degeneracy in the eigenvalues, the procedure is the following.
- We compute the symplectic eigenvalues \(z''\sb{j}\) and the symplectic matrix \(S''\) from the Williamson decomposition \(\sigma\modepart{3\cdots M} = S'' Z'' \transpose{(S'')}\).
- We find the \(\chi\) largest eigenvalues \((\lambda\modepart{2}\sb{\alpha\sb2})^2\) and the associated occupation numbers \(\multi{m}\sb{\alpha\sb2} \in \{0,\dotsc,d\}^{M-2}\) from Eq. \eqref{eq:svd-reduced-state-3-to-M}.
- We identify the \(b\modepart{3\cdots M}\sb{\alpha\sb2}\) vectors associated to the largest eigenvalues with \(\spunitarymp[3\cdots M]{S''} \ns{\multi{m}\sb{\alpha\sb2}}\).
In the end, we get
\[\begin{equation*} \begin{split} A^{(2,m_2)}_{\alpha_1,\alpha_2} &= \innp{% \spunitarymp[2\cdots M]{S'} \ns{\multi{m}\sb{\alpha\sb1}} }{ \ns{m_2} \otimes (\spunitarymp[3\cdots M]{S''} \ns{\multi{m}\sb{\alpha\sb2}}) } =\\ &= \innp{% \ns{\multi{m}\sb{\alpha\sb1}} }{ \adj{\spunitarymp[2\cdots M]{S'}} (\id\modepart{2} \otimes \spunitarymp[3\cdots M]{S''}) (\ns{m_2} \otimes \ns{\multi{m}\sb{\alpha\sb2}}) } \end{split} \end{equation*}\]
where we can decompose
\[\begin{multline} \adj{\spunitarymp[2\cdots M]{S'}} (\id\modepart{2} \otimes \spunitarymp[3\cdots M]{S''}) = \\ = \spunitarymp[2\cdots M]{(S')^{-1} (\imat[2] \oplus S'')} = \spunitarymp[2\cdots M]{V_1'} \spunitarymp[2\cdots M]{D'} \spunitarymp[2\cdots M]{V_2'}, \end{multline}\]
as we did for Eq. \eqref{eq:first-mps-tensor-simplified}, and then use Eq. \eqref{eq:franck-condon-formula}. We now have
\[\psi = \sum_{m_1=1}^{d} \sum_{m_2=1}^{d} \sum_{\alpha_1=1}^{\chi} \sum_{\alpha_2=1}^{\chi} A^{(1,m_1)}_{\alpha_1} A^{(2,m_2)}_{\alpha_1,\alpha_2} \ns{m_1} \otimes \ns{m_2} \otimes b\modepart{3\cdots M}_{\alpha_2}\]
and we repeat this process until we reach the last mode:
\[\psi = \sum_{m_1=1}^{d} \! \dotsi \!\! \sum_{m_{M-1}=1}^{d} \sum_{\alpha_1=1}^{\chi} \! \dotsi \!\! \sum_{\alpha_{M-1}=1}^{\chi} A^{(1,m_1)}_{\alpha_1} \dotsm A^{(M-1,m_{M-1})}_{\alpha_{M-2},\alpha_{M-1}} \ns{m_1} \otimes \dotsb \otimes \ns{m_{M-1}} \otimes b\modepart{M}_{\alpha_{M-1}}\]
The last step is to write
\[b\modepart{M}_{\alpha_{M-1}} = \spunitarymp[M]{S^{(M-1)}} \ns{m_{\alpha_{M-1}}},\]
where \(S^{(M-1)}\) and \(m\sb{\alpha\sb{M-1}}\) are already known from the last iteration, when we computed \(A^{(M-1,m\sb{M-1})}\), and to compute the last tensor of the MPS as
\[A^{(M,m_M)}_{\alpha_{M-1}} = \innp{\ns{m_M}}{\spunitarymp[M]{S^{(M-1)}} \ns{m_{\alpha_{M-1}}}}\]
with the usual technique. In the end, we have written the (approximated) MPS form of the initial state,
\[\psi = \sum_{m_1=1}^{d} \! \dotsi \!\! \sum_{m_M=1}^{d} \sum_{\alpha_1=1}^{\chi} \! \dotsi \!\! \sum_{\alpha_{M-1}=1}^{\chi} A^{(1,m_1)}_{\alpha_1} \dotsm A^{(M-1,m_{M-1})}_{\alpha_{M-2},\alpha_{M-1}} A^{(M,m_M)}_{\alpha_{M-1}} \ns{m_1} \otimes \dotsb \otimes \ns{m_M}.\]