⇦Chinese Physics Letters, 2019, Vol. 36, No. 9, Article code 090301 Methods for Derivation of Density Matrix of Arbitrary Multi-Mode Gaussian States from Its Phase Space Representation * Sheng-Li Zhang (张胜利)1**, Song Yang (杨颂)2,1 Affiliations 1Beijing Key Laboratory of Nanophotonics & Ultrafine Optoelectronic Systems, School of Physics, Beijing Institute of Technology, Beijing 100081 2Beijing Institute of Space Mechanics & Electricity, Beijing 100076 Received 9 April 2019, online 23 August 2019 *Supported by the National Natural Science Foundation of China under Grant Nos 11574400 and 11204379, the Beijing Institute of Technology Research Fund Program for Young Scholars, and the NSFC-ICTP Proposal under Grant No 11981240356.
**Corresponding author. Email: zhangshengli@bit.edu.cn Citation Text: Zhang S L and Yang S 2019 Chin. Phys. Lett. 36 090301    Abstract We present a method for derivation of the density matrix of an arbitrary multi-mode continuous variable Gaussian entangled state from its phase space representation. An explicit computer algorithm is given to reconstruct the density matrix from Gaussian covariance matrix and quadrature average values. As an example, we apply our method to the derivation of three-mode symmetric continuous variable entangled state. Our method can be used to analyze the entanglement and correlation in continuous variable quantum network with multi-mode quantum entanglement states. DOI:10.1088/0256-307X/36/9/090301 PACS:03.67.Mn, 03.67.Hk, 42.50.Dv © 2019 Chinese Physics Society Article Text Phase space method is an efficient tool to describe entanglement and dynamics evolution of continuous variable (CV) Gaussian entangled state.[1,2] However, the drawback of phase space is that it applies only to the Gaussian state and Gaussian channels. In practice, non-Gaussian state and non-Gaussian processing are more general. Actually, non-Gaussian operations are the necessary building blocks for continuous variable entanglement distillation,[3–5] error correction[6] and many other information processing protocols.[7–9] After non-Gaussian operation, the states in the output are non-Gaussian and the density matrix in the Fock state basis could provide all the information of the non-Gaussian state. The most straightforward way to obtain the density matrix is to represent every interaction and every state with photon number basis. Thus, the number of the density matrix entries scales exponentially with the optical mode number $N$, which imposes a great limit to the optical modes that we can process and analyze in realistic CV quantum networks. An alternative method is to write the non-Gaussian state as a linear combination of a few Gaussian states in phase space. Then, by deriving the density matrix of each Gaussian state one by one and again using the linear combination, one can finally obtain the density matrix of the non-Gaussian state.[10–12] This method has been successfully applied to describe the quantum states generated in continuous variable entanglement distillation,[13] photon catalysis[14] and noiseless linear amplification.[15] Deriving the density matrix from the phase space is non-trivial and requires large amounts of calculation. However, research in this area is rare and the already-known methods can be divided into two categories. First, the Hermite polynomials method: for a Gaussian state, the Husimi $Q$ function of the quantum state has a Gaussian form and the density matrix entries in the Hilbert space can be expressed by terms of multivariate Hermite polynomials.[16] Second, the characteristic functions method: the density matrix entries can be expressed by the multi-variant differential of the characteristic function of the Gaussian state.[17,18] In Ref.  [17] the density matrix elements in the Hilbert Space are obtained assuming that the Gaussian entangled state resides in two-mode Hilbert space and has zero-displacement (i.e., quadrature average values are zero). However, their density matrix is not normalized and later in Ref.  [12] we improve Chen's idea and develop an auto-normalized state transfer algorithm for zero-displaced two-mode Gaussian entanglement. Up till now, most of the studies in this area have been restricted in two-mode entangled states, which imposes a general limit on the number of modes that we can deal with in the continuous variable information processing. In this Letter, we present an algorithm to derive the density matrix of an arbitrary multi-mode Gaussian entangled state from its phase space representation. This result is of great value in the study of the exact photon number population in the Fock state space and may find widespread applications in the study of the entanglement dynamics of multi-mode CV entanglement.[19,20] For an $N$-mode CV quantum system, each mode is an infinite dimensional subsystem and can be conveniently expressed by field quadrature position operator $\hat{x}_k=(\hat{a}_k+\hat{a}_k^†)/\sqrt{2}$ and momentum operator $\hat{p}_k=(\hat{a}_k-\hat{a}_k^†)/(i\sqrt{2})$, where $k=1, 2,\ldots N$ denotes the index of mode, $\hat{a}_k$ and $\hat{a}_k^†$ denote the photon annihilation and creation operators for the $k$th mode. For any density matrix ${\rho}$ in the Hilbert space, there is a corresponding characteristic function defined in phase space $\chi(\zeta)={\rm Tr}[\exp[i\hat{R}\zeta^{\rm T} {\rho}]]$, with $\hat{R}=(\hat{x}_1, \hat{p}_1, \hat{x}_2, \hat{p}_2,\ldots, \hat{x}_N, \hat{p}_N)$ and $\zeta\in \mathbb{R}^{2N}$ being a $2N$-dimensional real row vector. The characteristic function of a Gaussian state can be completely described by the covariance matrix and quadrature average values, $\chi(\zeta)\equiv\chi_G(\zeta,V,\overline{\boldsymbol R})=\exp[-\frac{1}{2}\zeta V \zeta^ {\rm T} +i \overline{\boldsymbol R}\zeta^{\rm T} ]$, with $V$ being the $2N\times 2N$ covariance matrix (CM) $$\begin{alignat}{1} \!\!\!\!V_{lm}=\frac{1}{2}{\rm Tr}[(\hat{R}_l\hat{R}_m+\hat{R}_m\hat{R}_l){\rho}]-{\rm Tr}[\hat{R}_l{\rho}]{\rm Tr}[\hat{R}_m {\rho}],~~ \tag {1} \end{alignat} $$ with $l, m =1, 2, 3,\ldots,2N$ and $\overline{\boldsymbol R}$ being the average values of quadrature operators $\hat{R}$. To derive the density matrix for characteristic function $\chi$, we note that the relation follows $$\begin{align} \rho=\int \frac{d\boldsymbol{\mu}}{\pi^N} \,{\rm Tr} [\rho D(\boldsymbol \mu)]D(-\boldsymbol\mu),~~ \tag {2} \end{align} $$ with $\boldsymbol{\mu}=(\mu_1,\mu_2,\ldots,\mu_N)\in\mathbb{C}^N, D(\boldsymbol{\mu})=\exp[\boldsymbol{\mu}\boldsymbol{a}^†-\boldsymbol{\mu}^*\boldsymbol{a}] =\exp[(\boldsymbol{\mu},{{\boldsymbol \mu}^*})(\boldsymbol{a}^†,-\boldsymbol{a})^{{\rm T}}]$ being the $N$-mode displacement operator. The matrix elements of $\rho$ can be conveniently obtained through[12] $$\begin{align} &\langle k_1,k_2,\ldots,k_N|D(-\boldsymbol{\mu})|m_1,m_2,\ldots, m_N\rangle\\ =\,&\frac{\exp[-\frac{|\boldsymbol{\mu}|^2}{2}]}{\sqrt{\prod_{i=1}^N k_i! \prod_{i=1}^N m_i!}}\prod_{i=1}^N \partial_{t_i}^{k_i}\prod_{i=1}^N \partial_{t_i^\prime}^{m_i} \\ &\times\exp[{\boldsymbol t}{{\boldsymbol t}^\prime}^{\rm T}-{\boldsymbol t}\boldsymbol{\mu}^{\rm T}+{{\boldsymbol t}^\prime}{{\boldsymbol \mu}^*}^{\rm T}]|_{{\boldsymbol t}={{\boldsymbol t}^\prime}=0},~~ \tag {3} \end{align} $$ where ${\boldsymbol t}=(t_1,t_2,\ldots,t_N)$ and ${{\boldsymbol t}^\prime}=(t_1^\prime,t_2^\prime,\ldots,t_N^\prime)$ are $N$-dimensional real vectors. It should be noted that $(\boldsymbol{a}^†,-\boldsymbol{a})^{\rm T}=i L_N(\hat{x}_1,\hat{p}_1,\hat{x}_2,\hat{p}_2,\ldots,\hat{x}_N,\hat{p}_N)^{\rm T}$, $$ L_N=\Big(\left(\begin{array}{cc} \frac{-i}{\sqrt{2}} &\frac{-1}{\sqrt{2}}\\ \frac{i}{\sqrt{2}} & \frac{-1}{\sqrt{2}} \end{array}\right)\otimes I_N\Big)P, $$ $$ P_{kl}=\begin{cases} \delta_{2k-1,l} & k\le N,\\ \delta_{2(k-N),l} & k > N, \end{cases}~~ \tag {4} $$ with $I_N$ the identity matrix in $N$ dimensions. We obtain (after integrating over $\boldsymbol{\mu}$) the following theorem. Theorem 1: $$\begin{align} &\langle k_1,k_2,\ldots,k_N|\rho|m_1,m_2,\ldots, m_N\rangle\\ =\,&\frac{\prod_{i=1}^N \partial_{t_i}^{k_i}\prod_{i=1}^N \partial_{t_i^\prime}^{m_i} F |_{{\boldsymbol t}={{\boldsymbol t}^\prime}=0}} {\sqrt{\prod_{i=1}^N k_i! \prod_{i=1}^N m_i!}},~~ \tag {5} \end{align} $$ with $F=\exp[\frac{1}{2}({\boldsymbol t},{{\boldsymbol t}^\prime})R({\boldsymbol t},{{\boldsymbol t}^\prime})^{\rm T}]\Big /\sqrt{\det(V+I_{2N}/2)}$, $R =\sigma_x\otimes I_N+(\sigma_z\otimes I_N) L_N^{*}(V+I_{2N}/2)^{-1}L_N^†(\sigma_z\otimes I_N)$, and $\sigma_x$, $\sigma_z$ the usual Pauli matrices. Then, one can check that the state $\rho$ as given by Eq. (5) is now automatically normalized; i.e., ${\rm Tr}[\rho]=1$. In the following, we further simplify this theorem and finally present a computer algorithm that enables us to perform a numerical calculation. We emphasize that the theorem presented here is quite general and versatile, making it very useful for further applications, especially when entanglement has to be quantitatively evaluated with density matrix in the Fock state basis. In what follows, we give a detailed derivation of the density matrix from a Gaussian-state covariance matrix with $N=3$. For simplicity, we use the notation $t_{N+i}=t_i^\prime (i=1,2,\ldots N)$. For a given three-mode Gaussian state with CM $V$ and zero displacements, the density matrix is given by $$\begin{align} &\langle k_1,k_2,k_3|\rho|m_1,m_2,m_3\rangle\\ =\,&\frac{\partial_{t_1}^{k_1}\partial_{t_2}^{k_2}\partial_{t_3}^{k_3} \partial_{t_4}^{m_1}\partial_{t_5}^ {m_2}\partial_{t_6}^{m_3}\exp[\tau]|_{t_1=t_2=t_3=t_4=t_5=t_6=0}} {\sqrt{k_1!k_2!k_3!m_1!m_2!m_3!\det(V+I_6/2)}},~~ \tag {6} \end{align} $$ with $R=\sigma_x\otimes I_3+\sigma_z\otimes I_3 L_3^{*}(V+I_6/2)^{-1}L_3^†\sigma_z\otimes I_3, \tau=[\frac{1}{2}(t_1,t_2,t_3,t_4,t_5,t_6) R(t_1,t_2,t_3,t_4,t_5,t_6)^{\rm T}]=\sum_{i=1}^6 \sum_{j=i}^6 c_{ij} t_it_j$. In this equation, $c_{ij}$ can be expressed by a $6\times 6$ upper-triangular matrix $$\begin{align} c_{ij}=\begin{cases} R_{ii}/2& i=j,\\ (R_{ij}+R_{ji})/2& i < j,\\ 0 & i>j. \end{cases}~~ \tag {7} \end{align} $$ The exponential can be calculated as $$\begin{align} \exp[\tau]=\,&\sum_{n=0}^\infty \frac{1}{n!}(\sum_{i=1}^6 \sum_{j=i}^6 c_{ij} t_it_j)^n~~ \tag {8} \end{align} $$ $$\begin{align} =\,&\sum_{n=0}^\infty \sum_{{\hat{\boldsymbol n}} \in {\it \Omega}(n)}\frac{1}{\prod_{i=1}^6\prod_{j=i}^6 {\hat{\boldsymbol n}}_{ij}!}\prod_{i=1}^6\prod_{j=i}^6(c_{ij}t_it_j)^{{\hat{\boldsymbol n}}_{ij}}\\ =\,&\sum_{n=0}^\infty \sum_{{\hat{\boldsymbol n}} \in {\it \Omega}(n) }\frac{\prod_{i=1}^6\prod_{j=i}^6(c_{ij})^{{\hat{\boldsymbol n}}_{ij}}} {\prod_{i=1}^6\prod_{j=i}^6 {\hat{\boldsymbol n}}_{ij}!}\\ & \times t_1^{2{\hat{\boldsymbol n}}_{11}+{\hat{\boldsymbol n}}_{12}+{\hat{\boldsymbol n}}_{13} +{\hat{\boldsymbol n}}_{14} +{\hat{\boldsymbol n}}_{15}+{\hat{\boldsymbol n}}_{16}}\\ & \times t_2^{2{\hat{\boldsymbol n}}_{22}+{\hat{\boldsymbol n}}_{12}+{\hat{\boldsymbol n}}_{23} +{\hat{\boldsymbol n}}_{24}+{\hat{\boldsymbol n}}_{25}+{\hat{\boldsymbol n}}_{26}}\\ &\times t_3^{2{\hat{\boldsymbol n}}_{33}+{\hat{\boldsymbol n}}_{13}+{\hat{\boldsymbol n}}_{23} +{\hat{\boldsymbol n}}_{34}+{\hat{\boldsymbol n}}_{35}+{\hat{\boldsymbol n}}_{36}}\\ & \times t_4^{2{\hat{\boldsymbol n}}_{44}+{\hat{\boldsymbol n}}_{14}+{\hat{\boldsymbol n}}_{24} +{\hat{\boldsymbol n}}_{34}+{\hat{\boldsymbol n}}_{45}+{\hat{\boldsymbol n}}_{46}}\\ &\times t_5^{2{\hat{\boldsymbol n}}_{55}+{\hat{\boldsymbol n}}_{15}+{\hat{\boldsymbol n}}_{25} +{\hat{\boldsymbol n}}_{35}+{\hat{\boldsymbol n}}_{45}+{\hat{\boldsymbol n}}_{56}}\\ &\times t_6^{2{\hat{\boldsymbol n}}_{66}+{\hat{\boldsymbol n}}_{16}+{\hat{\boldsymbol n}}_{26} +{\hat{\boldsymbol n}}_{36}+{\hat{\boldsymbol n}}_{46}+{\hat{\boldsymbol n}}_{56}},~~ \tag {9} \end{align} $$ where ${\hat{\boldsymbol n}}$ is again a $6\times 6$ upper-triangular matrix, and $\boldsymbol{\hat{n}_{ij}}$ is the corresponding matrix element for the $i$th row and $j$th column $$ {\hat{\boldsymbol n}}=\left(\begin{array}{ccccc} {\hat{\boldsymbol n}}_{11} {\hat{\boldsymbol n}}_{12}&\ldots,{\hat{\boldsymbol n}}_{15} &{\hat{\boldsymbol n}}_{16}\\ 0 {\hat{\boldsymbol n}}_{22}& {\hat{\boldsymbol n}}_{23} \ldots,& {\hat{\boldsymbol n}}_{26}\\ \vdots \vdots & \vdots \vdots &\vdots\\ 0 0& 0 {\hat{\boldsymbol n}}_{55}& {\hat{\boldsymbol n}}_{56}\\ 00&00&{\hat{\boldsymbol n}}_{66} \end{array}\right),~~ \tag {10} $$ with ${\it \Omega}(n)$ being the set of all possible matrices ${\hat{\boldsymbol n}}$ such that $$\begin{align} {\it \Omega}(n)=\{{\hat{\boldsymbol n}}| \sum_{i=1}^6\sum_{j=i}^6{\hat{\boldsymbol n}}_{ij}=n \}.~~ \tag {11} \end{align} $$ Under the condition $t_i=0$ in Eq. (6), the density matrix can be explicitly written as $$\begin{align} &\langle k_1,k_2,k_3|\rho|m_1,m_2,m_3\rangle\\ =\,&\frac{1}{\sqrt{\det(V+I_6/2)}}\sum_{n=0}^\infty\sum_{\boldsymbol{\hat{n}\in {\it \Omega}(n)}}\\ &\cdot\prod_{i=1}^6\prod_{j=i}^6(c_{ij})^{{\hat{\boldsymbol n}}_{ij}} f(k_1,k_2,k_3,{\hat{\boldsymbol n}}),~~ \tag {12} \end{align} $$ with $$\begin{align} &f(k_1,k_2,k_3,{\hat{\boldsymbol n}})\\ =\,&\sqrt{{k_1 \choose {\hat{\boldsymbol n}}_{11},{\hat{\boldsymbol n}}_{11},{\hat{\boldsymbol n}}_{12}, {\hat{\boldsymbol n}}_{13},{\hat{\boldsymbol n}}_{14},{\hat{\boldsymbol n}}_{15},{\hat{\boldsymbol n}}_{16}}}\\ &\times\sqrt{{k_2\choose {\hat{\boldsymbol n}}_{22},{\hat{\boldsymbol n}}_{22},{\hat{\boldsymbol n}}_{12}, {\hat{\boldsymbol n}}_{23},{\hat{\boldsymbol n}}_{24},{\hat{\boldsymbol n}}_{25},{\hat{\boldsymbol n}}_{26}}}\\ &\times\sqrt{{k_3\choose{\hat{\boldsymbol n}}_{33},{\hat{\boldsymbol n}}_{33},{\hat{\boldsymbol n}}_{13}, {\hat{\boldsymbol n}}_{23},{\hat{\boldsymbol n}}_{34},{\hat{\boldsymbol n}}_{35},{\hat{\boldsymbol n}}_{36}}}\\ &\times\sqrt{{m_1\choose {\hat{\boldsymbol n}}_{44},{\hat{\boldsymbol n}}_{44},{\hat{\boldsymbol n}}_{14}, {\hat{\boldsymbol n}}_{24}, {\hat{\boldsymbol n}}_{34},{\hat{\boldsymbol n}}_{45},{\hat{\boldsymbol n}}_{46}}}\\ &\times\sqrt{{m_2\choose {\hat{\boldsymbol n}}_{55},{\hat{\boldsymbol n}}_{55},{\hat{\boldsymbol n}}_{15}, {\hat{\boldsymbol n}}_{25},{\hat{\boldsymbol n}}_{35},{\hat{\boldsymbol n}}_{45},{\hat{\boldsymbol n}}_{56}}}\\ &\times\sqrt{{m_3\choose {\hat{\boldsymbol n}}_{66},{\hat{\boldsymbol n}}_{66},{\hat{\boldsymbol n}}_{16}, {\hat{\boldsymbol n}}_{26},{\hat{\boldsymbol n}}_{36},{\hat{\boldsymbol n}}_{46},{\hat{\boldsymbol n}}_{56}}},~~ \tag {13} \end{align} $$ and ${k\choose n_1,n_2,n_3,n_4,n_5,n_6,n_7}=k!\delta_{k,\sum_{i=1}^7 n_i}/(\prod_{i=1}^7 n_i!)$ being the multinomial coefficient. Moreover, from the definition of multinomials, we know $$\begin{alignat}{1} &2n=2\sum_{i=1}^6\sum_{j=i}^6 n_{ij}=\sum_{i=1}^3k_i+\sum_{i=1}^3 m_i,~~ \tag {14} \end{alignat} $$ $$\begin{alignat}{1} &2{\hat{\boldsymbol n}}_{11}+{\hat{\boldsymbol n}}_{12}+{\hat{\boldsymbol n}}_{13}+{\hat{\boldsymbol n}}_{14} +{\hat{\boldsymbol n}}_{15}+{\hat{\boldsymbol n}}_{16}=k_1,~~ \tag {15} \end{alignat} $$ $$\begin{alignat}{1} &2{\hat{\boldsymbol n}}_{22}+{\hat{\boldsymbol n}}_{12}+{\hat{\boldsymbol n}}_{23}+{\hat{\boldsymbol n}}_{24} +{\hat{\boldsymbol n}}_{25}+{\hat{\boldsymbol n}}_{26}=k_2,~~ \tag {16} \end{alignat} $$ $$\begin{alignat}{1} &2{\hat{\boldsymbol n}}_{33}+{\hat{\boldsymbol n}}_{13}+{\hat{\boldsymbol n}}_{23}+{\hat{\boldsymbol n}}_{34} +{\hat{\boldsymbol n}}_{35}+{\hat{\boldsymbol n}}_{36}=k_3,~~ \tag {17} \end{alignat} $$ $$\begin{alignat}{1} &2{\hat{\boldsymbol n}}_{44}+{\hat{\boldsymbol n}}_{14}+{\hat{\boldsymbol n}}_{24}+{\hat{\boldsymbol n}}_{34} +{\hat{\boldsymbol n}}_{45}+{\hat{\boldsymbol n}}_{46}=m_1,~~ \tag {18} \end{alignat} $$ $$\begin{alignat}{1} &2{\hat{\boldsymbol n}}_{55}+{\hat{\boldsymbol n}}_{15}+{\hat{\boldsymbol n}}_{25}+{\hat{\boldsymbol n}}_{35} +{\hat{\boldsymbol n}}_{45}+{\hat{\boldsymbol n}}_{56}=m_2,~~ \tag {19} \end{alignat} $$ $$\begin{alignat}{1} &2{\hat{\boldsymbol n}}_{66}+{\hat{\boldsymbol n}}_{16}+{\hat{\boldsymbol n}}_{26}+{\hat{\boldsymbol n}}_{36} +{\hat{\boldsymbol n}}_{46}+{\hat{\boldsymbol n}}_{56}=m_3.~~ \tag {20} \end{alignat} $$ Finally, the matrix elements $\langle k_1 k_2 k_3|\rho|m_1 m_2 m_3\rangle$ are given by $$\begin{align} \!\!&\langle k_1,k_2,k_3|\rho|m_1,m_2,m_3\rangle\\ \!\!\!\!\!\!\!\!\!\!=\,&\sum_{\boldsymbol{\hat{n}\in {\it \Omega}(\sum_{i=1}^3\frac{k_i+ m_i}{2})}} \prod_{i=1}^6\prod_{j=i}^6(c_{ij})^{{\hat{\boldsymbol n}}_{ij}} \frac{f(k_1,k_2,k_3,{\hat{\boldsymbol n}})}{\sqrt{\det(V\!+\!I_6/2)}}.~~ \tag {21} \end{align} $$ To save computer time and storage, every mode is truncated at $\overline{D}$; i.e., only the contributions of the terms $|0\rangle,|1\rangle,\ldots,|D-1\rangle$ are finally kept. The value of $\overline{D}$ is chosen to be as large as possible such that our state transfer is performed with very high precision, corresponding to an error less than $5\times 10^{-5}$. A possible computer algorithm that realizes this state transfer looks as follows: Algorithm START (1) Set $\overline{D}=6$. (2) Generate all the possible $({\hat{\boldsymbol n}}_{11},{\hat{\boldsymbol n}}_{12},\ldots,{\hat{\boldsymbol n}}_{66})$ satisfying Eqs. (14)-(20). (3) FOR each $(k_1,k_2,k_3,m_1,m_2,m_3)$ such that $0\le k_1,k_2,k_3,m_1,m_2,m_3 \le \overline{D}-1$     IF $k_1+k_2+k_3+m_1+m_2+m_3$ is even               (i) Generate the set ${\it \Omega}(\sum_{i=1}^3 \frac{k_i+m_i}{2})$, i.e., all possible ${\hat{\boldsymbol n}}$ such that $$\begin{align} \sum_{i=1}^6\sum_{j=i}^6{\hat{\boldsymbol n}}_{ij}=\sum_{i=1}^3 \frac{k_i+m_i}{2}.~~ \tag {22} \end{align} $$               (ii) Calculate Eq. (21).         ELSE                          $$\begin{align} \langle k_1 k_2 k_3|\rho| m_1 m_2 m_3\rangle=0.~~ \tag {23} \end{align} $$         END IF    END FOR Algorithm END With the set ${\it \Omega}$ and the definition of $c_{ij}$ in Eq. (7), we can calculate the density matrix entry $\langle k_1,k_2,k_3|{\rho}|m_1,m_2,m_3\rangle$. Furthermore, by varying the values of $(k_1, k_2,k_3, m_1, m_2,m_3)$, the total density matrix ${\rho}$ can be numerically obtained. To check the correctness of our algorithm, here we use the genuinely multipartite entangled 3-mode symmetric Gaussian states as an example. These states can be experimentally prepared using a sequence of $N-1$ phase-free beam splitters and $N$ squeezed input states,[21–23] with the corresponding CM can be given by $$\begin{align} V_N=|\mathcal{I}\rangle\langle \mathcal{I}| \otimes \boldsymbol{\epsilon}+(\boldsymbol{\alpha}-\boldsymbol{\epsilon})\otimes I_N,~~ \tag {24} \end{align} $$ in which $|\mathcal{I}\rangle$ is the unnormalized $N$-dimensional real vector $|\mathcal{I}\rangle=(1,1,1\ldots, 1)^{\rm T}\in \mathbb{R}^N$ and $\boldsymbol{\alpha}=\mathrm{diag}(a,b), \boldsymbol{\epsilon}=\mathrm{diag}(c,d)$ are $2\times 2$ diagonal matrices, with (Eq. (32) in Ref.  [22] multiplied by a factor 2 due to the different definitions of CV) $$\begin{alignat}{1} a=\,&\frac{1}{2N}[e^{2r_1}+(N-1)e^{-2r_2}],\\ c=\,&\frac{1}{2N}(e^{2r_1}-e^{-2r_2}),\\ b=\,&\frac{1}{2N}[e^{-2r_1}+(N-1)e^{2r_2}],\\ d=\,&\frac{1}{2N}(e^{-2r_1}-e^{2r_2})~~ \tag {25} \end{alignat} $$ and $r_1$ being the squeezing parameter of one squeezed input and $r_2$ being the squeezing parameter of the other two equally squeezed inputs. For simplicity, we consider the so-called unbiased states for which $r_1=r_1(r_2)=\frac{1}{2}\ln[(N-1)\sinh(2r_2)]+\frac{1}{2} \ln[\sqrt{1+[(N-1)\sinh(2r_2)]^{-2}}+1]$. With our computer algorithm, we are able to calculate the density matrix $\rho$ for arbitrary three-mode Gaussian state. Assuming that the density matrix obtained from our algorithm is ${\rho}_{\mathrm{algo}}$. One good way to evaluate how well our algorithm works is to calculate the covariance matrix $V_{\mathrm{algo}}$ of our output matrix ${\rho}_{\mathrm{algo}}$. Then, we can use the maximal norm $||A||_\infty=\max\{a_{ij}\}$ to measure the error during our state derivation $$\begin{align} V_{\mathrm{err}}=||V_\mathrm{alog}-V_{N}||_\infty.~~ \tag {26} \end{align} $$ We can use the logarithmic negativity[24] to quantify the entanglement, which is more important in practical quantum communications.[25,26] The logarithmic negativity of a multipartite state $\rho$ is defined as[24] $$\begin{align} E_N(\rho)=\log_2(1+2\mathcal{N}(\rho)),~~ \tag {27} \end{align} $$ in which $\mathcal{N}(\rho)$ is given by the absolute value of the sum of negative eigenvalues of the partially transposed (e.g., with respect to the first mode) density operator $\rho^{{\it \Gamma}_{1}}$. For Gaussian entangled states, we can evaluate the logarithmic negativity from the CM $V_{N}$, which follows $$\begin{align} E_\mathrm{th}=\mathrm{max}[0,-\mathrm{log}_2 \widetilde{\nu}],~~ \tag {28} \end{align} $$ with $\nu$ being the minimal positive symplectic eigenvalues of $TV_{N}T$ for $T=(\begin{array}{cc} 1,0 0,-1 \end{array})\oplus I_4$.[27] Because this algorithm also gives us a method for calculating the density matrix of $\rho$ in the Fock state space, we can numerically find $E_N$ by a straightforward diagonalization of the partially transposed density matrix. The comparison between the theoretic and numerical values $|E_N-E_\mathrm{th}|$ gives us another figure of merit to evaluate the performance of our algorithm.
cpl-36-9-090301-fig1.png
Fig. 1. Numerical check of correctness of our numerical algorithm for $r_2=0.2$, $N=3$. (a) Error of trace of density matrix: $1-{\rm Tr}[{\rho}_{\rm algo}]$. (b) Error of covariance matrix. (c) Error of logarithmic negativity.
cpl-36-9-090301-fig2.png
Fig. 2. Performance of our numerical algorithm for increasing squeezing $r$. Other parameters are chosen as $N=3$, $\overline{D}=8$. (a) Error of trace of density matrix. (b) Error of covariance matrix. (c) Error of logarithmic negativity.
In Fig. 1, we take $r_2=0.2$ as an example and give a plot of the performance of our numerical algorithm for different values of truncation number $\overline{D}$. As shown in Fig. 1, the figure of merit of the truncated state ${\rho}_{\mathrm{algo}}$ quickly converges for a larger $\overline{D}$. When $\overline{D}=6$, the error due to our truncation is almost negligible. In particular, Fig. 1(a) shows the error of the trace of our derived state ${\rho}_{\rm algo}$, which once again proves the auto-normalization of our algorithm. In Fig. 2, we fix $\overline{D}=8$ and evaluate the performance of state transferring, for variations of squeezing $r$. The error increases with $r$ and when $r=0.6$ the average photon number in each modes of symmetric state is about 0.6342. Although this is much smaller than $\overline{D}=8$, we can see that this is still strong enough to influence the entanglement $E_N$. Finally, we give a remark on the storage complexity and computational complexity. In fact, for $\overline{D}=3,4,5,6,7,8$, there are $785, 11874, 1.3\times 10^5, 1.1\times 10^6, 7.18 \times 10^6, 4.0\times 10^7$ groups of possible $({\hat{\boldsymbol n}}_{11},{\hat{\boldsymbol n}}_{12},\ldots,{\hat{\boldsymbol n}}_{66})$. Especially for $\overline{D}=8$, it costs 806 Megabyte storage and 590 min to complete the numerical evaluation with a single-core 2.0 Ghz-CPU-equipped personal computer (PC). Large $\overline{D}$ requires more storage and more computation, which is necessary for multi-mode entanglement state with larger average photon number. Fortunately, we note that the computational complexity in the large squeezing scenario can be further improved thanks to the fact that our algorithm deals with the density matrix entries one by one, which can be certainly updated with a parallel computation technique and run with the multi-core environments. In conclusion, we have provided an auto-normalized computer method for evaluating arbitrary multi-mode continuous variable Gaussian entangled stats. Finally, we apply our method into the symmetric three-mode Gaussian entangled state. The performance of our algorithm is numerically evaluated, in terms of errors of trace of density matrix, covariance matrix, and logarithmic negativity. We envisage that our state transfer method will be widely used in characterizing CV entanglement dynamics evolution and multipartite nonlocality test[28] in infinite-dimensional quantum systems. References Quantum information with continuous variablesGaussian quantum informationDistilling Gaussian States with Gaussian Operations is ImpossibleGaussian Transformations and Distillation of Entangled Gaussian StatesCharacterization of Gaussian operations and distillation of Gaussian statesError Correction for Continuous Quantum VariablesRole of Initial Entanglement and Non-Gaussianity in the Decoherence of Photon-Number Entangled States Evolving in a Noisy ChannelOptimal estimation of losses at the ultimate quantum limit with non-Gaussian statesContinuous-variable quantum teleportation with non-Gaussian resourcesProposal for a Loophole-Free Bell Test Using Homodyne DetectionProposed Test of Quantum Nonlocality for Continuous VariablesLocal Gaussian operations can enhance continuous-variable entanglement distillationContinuous-variable-entanglement distillation with photon additionPhoton catalysis acting as noiseless linear amplification and its application in coherence enhancementImproving noiseless linear amplification for optical quantum communication with quadrature squeezingImproving entanglement concentration of Gaussian states by local displacementsFock-space inseparability criteria of bipartite continuous-variable quantum statesSeparability Criterion for One-Sided Gaussian ChannelsNonclassical correlations in continuous-variable non-Gaussian Werner statesMultipartite Entanglement for Continuous Variables: A Quantum Teleportation NetworkDetecting genuine multipartite continuous-variable entanglementExperimental Creation of a Fully Inseparable Tripartite Continuous-Variable StateComputable measure of entanglementEquivalence between Entanglement and the Optimal Fidelity of Continuous Variable TeleportationDistillation of Atmospherically Disturbed Continuous Variable Quantum Entanglement with Photon SubtractionEntanglement, Purity, and Information Entropies in Continuous Variable SystemsConnecting Quantum Contextuality and Genuine Multipartite Nonlocality with the Quantumness Witness
[1] Braunstein S and van Loock P 2005 Rev. Mod. Phys. 77 513
[2] Weedbrook C, Pirandola S, García-Patrón R, Cerf N, Ralph T, Shapiro J and Lloyd S 2012 Rev. Mod. Phys. 84 621
[3] Eisert J, Scheel S and Plenio M 2002 Phys. Rev. Lett. 89 137903
[4] Fiurášek J 2002 Phys. Rev. Lett. 89 137904
[5] Giedke G and Cirac J 2002 Phys. Rev. A 66 032316
[6] Braunstein S 1998 Phys. Rev. Lett. 80 4084
[7] Allegra M, Giorda P and Paris M 2010 Phys. Rev. Lett. 105 100503
[8] Adesso G, Dell'Anno F, Siena S, Illuminati F and Souza L 2009 Phys. Rev. A 79 040305
[9] Dell'Anno F, De Siena S, Albano L and Illuminati F 2007 Phys. Rev. A 76 022301
[10] Garcá-Patrón R, Fiuráek J, Cerf N, Wenger J, Tualle-Brouri R and Grangier P 2004 Phys. Rev. Lett. 93 130409
[11] Nha H and Carmichael H 2004 Phys. Rev. Lett. 93 020401
[12] Zhang S L and van Loock P 2011 Phys. Rev. A 84 062309
[13] Zhang S L, Dong Y L, Zou X B, Shi B S and Guo G C 2013 Phys. Rev. A 88 032324
[14] Zhang S L and Zhang X D 2018 Phys. Rev. A 97 043830
[15] Yang S, Zhang S L, Zou X B, Bi S W and Lin X L 2013 Phys. Rev. A 87 024302
[16] Fiurášek J 2011 Phys. Rev. A 84 012335
[17] Chen X 2007 Phys. Rev. A 76 022309
[18]Perelomov A 1986 Generalized Coherent States (Berlin: Springer Verlag)
[19] Hoelscher-Obermaier J and van Loock P 2010 arXiv:1001.2225
[20] Tatham R, Mišta J, Adesso G and Korolkova N 2012 Phys. Rev. A 85 022326
[21] van Loock P and Braunstein S 2000 Phys. Rev. Lett. 84 3482
[22] van Loock P and Furusawa A 2003 Phys. Rev. A 67 052315
[23] Aoki T et al 2003 Phys. Rev. Lett. 91 080404
[24] Vidal G and Werner R 2002 Phys. Rev. A 65 032314
[25] Adesso G and Illuminati F 2005 Phys. Rev. Lett. 95 150503
[26] Zhang S L, Guo J S, Shi J H and Zou X B 2016 Chin. Phys. Lett. 33 070303
[27] Adesso G, Serafini A and Illuminati F 2005 Open Syst. Inf. Dyn. 12 189
[28] Chen X, Su H Y and Chen J L 2016 Chin. Phys. Lett. 33 010302