# Maximal entanglement distillation with single-photon subtraction

Kun Zhang1 Jietai Jing1,3,4,5<EMAIL_ADDRESS>Nicolas Treps2 Mattia
Walschaers2<EMAIL_ADDRESS>1State Key Laboratory of Precision
Spectroscopy, Joint Institute of Advanced Science and Technology, School of
Physics and Electronic Science, East China Normal University, Shanghai 200062,
China  2Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-Université
PSL, Collège de France, 4 place Jussieu, F-75252 Paris, France  3CAS Center
for Excellence in Ultra-intense Laser Science, Shanghai 201800, China
4Department of Physics, Zhejiang University, Hangzhou 310027, China
5Collaborative Innovation Center of Extreme Optics, Shanxi University,
Taiyuan, Shanxi 030006, China

(August 29, 2024)

###### Abstract

Entanglement is an indispensable quantum resource for quantum information
technology. In continuous-variable quantum optics, photon subtraction can
increase the entanglement between Gaussian states of light, but for mixed
states the extent of this entanglement increase is poorly understood. In this
work, we use an entanglement measure based the Rényi-2 entropy to prove that
single-photon subtraction increases bipartite entanglement by no more than
$\log 2$. This value coincides with the maximal amount of bipartite
entanglement that can be achieved with one photon.

Quantum entanglement A. Einstein  is commonly seen as an important resource
for quantum metrology V. Giovannetti , quantum communication C. H. Bennett2
and quantum computation P. W. Shor . The increase of entanglement is an
important task, since it can be used to strengthen the precision of quantum
metrology S. Steinlechner , the fidelity of quantum information processing B.
Schumacher ; S. Lloyd  and the efficiency of quantum computing R. Raussendorf
; E. Knill . However, in many platforms entanglement is fragile and rapidly
vanishes due to environment-induced decoherence. In this regard, light is a
particularly interesting platform due to its intrinsic resilience against its
environment.

Continuous-variable (CV) quantum optics S. L. Braunstein  allows for the
deterministic generation of Gaussian entanglement between modes of light. In
particular, it was used to produce large-scale entangled states M. Chen ; S.
Armstrong ; X. Su ; J. Roslund ; S. Gerke ; J.-i. Yoshikawa ; K. Zhang ; W.
Wang . However, because the entangled modes are Gaussian, they do not offer a
genuine quantum advantage. To implement protocols that cannot be efficiently
simulated by classical computers non-Gaussian states are required S. D.
Bartlett ; A. Mari ; S. Rahimi-Keshari .

A common method to produce non-Gaussian modes of light is adding A. Zavatta
or subtracting A. Ourjoumtsev ; V. Parigi ; Y.-S. Ra ; J. Wenger  a photon.
Furthermore, this operation can also increase the entanglement between optical
modes. For example, by adding/subtracting a photon to/from a two-mode squeezed
vacuum state, entanglement of the resulting states is greater than that of the
original squeezed state Yang Yang . The entanglement of two-mode squeezed
vacuum state was shown to be progressively enhanced with the number of photons
added (subtracted) when acting on one mode only Carlos Navarrete-Benlloch .
Similar results are obtained for single-photon addition and subtraction for
pure multimode Gaussian systems M. Walschaers2 ; M. Walschaers1 ; M.
Walschaers3 .

Entanglement distillation is much harder to study for mixed states. Typical
entanglement measures for multimode mixed states such as entanglement of
formation require convex roof constructions R. Horodecki ; C. H. Bennett ; G.
Toth  and the evaluation of the von Neumann entropy. Both are problematic for
mixed non-Gaussian states. Here, we focus on the entanglement measure based on
Rényi-2 entropy as a viable alternative, following Kim_2010 ; G. Adesso1 .
Because the Rényi-2 entropy is related to the purity of the state, it can be
obtained directly from the Wigner function. As such, we will quantify the
increase of entanglement due to photon addition or subtraction.

In this work, we will restrict ourselves to discussing bipartite entanglement.
As shown in Fig. 1(a), for an arbitrary global pure state
$\hat{\rho}=|\psi\rangle\langle\psi|$, the entanglement between two subsystems
$\mathcal{A}$ and $\mathcal{B}$ can be measured through the Rényi-$2$ entropy
as

$\mathcal{E}_{R}(|\psi\rangle)=-\mathrm{log}\mu(\hat{\rho}_{\mathcal{A}}),$
(1)

where $\hat{\rho}_{\mathcal{A}}=\mathrm{tr}_{\mathcal{B}}\hat{\rho}$ is the
reduced density operator for subsystem $\mathcal{A}$. The measure is based on
the purity $\mu_{\mathcal{A}}=\mathrm{tr}(\hat{\rho}_{\mathcal{A}}^{2})$,
which can directly be calculated from the Wigner function by

$\mu_{\mathcal{A}}=(4\pi)^{{m_{\mathcal{A}}}}\int_{\mathbb{R}^{2{m_{\mathcal{A}}}}}|W_{\mathcal{A}}(\beta)|^{2}\mathrm{d}^{2{m_{\mathcal{A}}}}\beta,$
(2)

where ${m_{A}}$ is the total number of modes in subsystem $\mathcal{A}$.
Importantly, this method works both for Gaussian and non-Gaussian states. Its
extension to mixed states can be obtained by the following convex
construction:

${\cal
E}_{R}(\hat{\rho})=\inf_{\\{p(\lambda),\left\lvert\psi_{\lambda}\right\rangle\\}}\int{\rm
d}\lambda\,p(\lambda){\cal E}_{R}(\left\lvert\psi_{\lambda}\right\rangle),$
(3)

where we minimize over all decompositions $\hat{\rho}=\int{\rm
d}\lambda\,p(\lambda)\left\lvert\psi_{\lambda}\right\rangle\left\langle\psi_{\lambda}\right\rvert$.
We use this measure to focus on the upper limit of the entanglement increase.
We will show that for any Gaussian state –mixed or pure– the maximal increase
of bipartite “Rényi-2 entanglement” through single-photon addition/subtraction
is $\log 2$.

In this work, we focus on entanglement increase via photon subtraction. This
operation is theoretically described by the action of an annihilation operator
on the (in our case Gaussian) state. We assume that this operation is
implemented on subsystem $\mathcal{A}$, as shown in Fig. 1. To study the
increase of entanglement, we start by analysing the purity of the photon
subtracted state in subsystem $\mathcal{A}$:

$\displaystyle\hat{\rho}^{-}_{\mathcal{A}}=\frac{\hat{a}_{g}\hat{\rho}_{\mathcal{A}}\hat{a}_{g}^{{\dagger}}}{\mathrm{tr}(\hat{a}_{g}^{{\dagger}}\hat{a}_{g}\hat{\rho}_{\mathcal{A}})}.$
(4)

The label $g$ indicates the mode in which the photon is subtracted. Our main
tool to evaluate the purity of the state (4) are Gaussian transformations C.
Weedbrook . As shown in Fig. 1(b), we can write a thermal decomposition form
for an arbitrary subsystem $\mathcal{A}$ as follows

$\displaystyle\hat{\rho}_{{\mathcal{A}}}=\mathcal{\hat{D}}\hat{U}\mathop{\otimes}\limits_{i=1}^{{m_{A}}}\hat{\rho}_{i}\hat{U}^{{\dagger}}\mathcal{\hat{D}}^{{\dagger}},$
(5)

where $\hat{U}$ is a Bogoliubov transformation and
$\mathcal{\hat{D}}=\prod\limits_{i=1}^{m_{\mathrm{A}}}\hat{D}(\alpha_{i})$ is
a displacement operator. Furthermore, single-mode thermal states
$\hat{\rho}_{i}$ can be fully characterised by a covariance matrix
$V_{i}=\mathrm{diag}(n_{i},n_{i})$, where $n_{i}$ is the fraction of thermal
noise compared to shot noise. Since the Gaussian unitary
$\mathcal{\hat{D}}\hat{U}$ does not affect the purity, the purity of subsystem
$\hat{\rho}_{{\mathcal{A}}}$ can be directly expressed as
$\mu_{{\mathcal{A}}}=1/\prod\limits_{i=1}^{m_{\mathcal{A}}}n_{i}$.

As shown by the green circles in Fig. 1(b) and Fig. 1(c), if we subtract a
photon from the mode $g$ of subsystem ${\mathcal{A}}$, we can use Bogoliubov
transform to convert the photon-subtraction operation at mode $g$ into a
combined operator $\hat{b}$, which adds and subtracts a photon on each thermal
mode

$\displaystyle\hat{\rho}_{{\mathcal{A}}}^{-}=\frac{\mathcal{\hat{D}}\hat{U}\hat{b}\mathop{\otimes}\limits_{i=1}^{{m_{\mathcal{A}}}}\hat{\rho}_{i}\hat{b}^{{\dagger}}\hat{U}^{{\dagger}}\mathcal{\hat{D}}^{{\dagger}}}{\mathrm{tr}(\hat{b}^{{\dagger}}\hat{b}\mathop{\otimes}\limits_{i=1}^{{m_{\cal
A}}}\hat{\rho}_{i})},\text{ with
}\quad\hat{b}=\hat{U}^{{\dagger}}\mathcal{\hat{D}}^{{\dagger}}\hat{a}_{g}\mathcal{\hat{D}}\hat{U}.$
(6)

When $m_{\mathcal{A}}$ is the number of modes in subsystem $\mathcal{A}$, we
define
$\vec{\hat{a}}^{{\dagger}}=(\hat{a}_{1}^{{\dagger}},\cdots,\hat{a}_{{m_{\mathcal{A}}}}^{{\dagger}})^{T}$
and $\vec{\hat{a}}=(\hat{a}_{1},\cdots,\hat{a}_{{m_{\mathcal{A}}}})^{T}$. From
the properties of Gaussian unitary transformations, we find that

$\hat{U}^{{\dagger}}\mathcal{\hat{D}}^{{\dagger}}\vec{\hat{a}}\mathcal{\hat{D}}\hat{U}=K\vec{\hat{a}}^{{\dagger}}+L\vec{\hat{a}}+\vec{\alpha},$
(7)

with $L^{{\dagger}}L-K^{{\dagger}}K=\mathds{1}$ and
$K^{{\dagger}}L=-L^{{\dagger}}K$. We can thus explicitly write that
$\hat{b}=\vec{k}\cdot\vec{\hat{a}}^{{\dagger}}+\vec{l}\cdot\vec{\hat{a}}+\alpha_{g},$
where $\vec{k}=(k_{1},\cdots,k_{{m_{\mathcal{A}}}})$ and
$\vec{l}=(l_{1},\cdots,l_{{m_{\mathcal{A}}}})$ are
${m_{\mathcal{A}}}$-dimensional complex vectors, and $\alpha_{g}$ is a complex
number.

Based on both the Williamson decomposition (5) and Bogoliubov transformation
(6), we can derive (see Supplementary Material) a general expression of
relative purity

$\displaystyle\frac{\mu_{{\mathcal{A}}}^{-}}{\mu_{{\mathcal{A}}}}=$
$\displaystyle\frac{1}{2}+\Bigg{[}\frac{1}{2}(\sum_{i=1}^{{m_{\mathcal{A}}}}\frac{\tilde{{N}}_{i}}{n_{i}})^{2}+\frac{1}{2}|\alpha_{g}|^{4}+|\sum_{i=1}^{{m_{\mathcal{A}}}}k_{i}l_{i}\frac{n_{i}^{2}-1}{2n_{i}}|^{2}$
(8)
$\displaystyle+{\alpha_{g}^{\ast}}^{2}(\sum_{i=1}^{{m_{\mathcal{A}}}}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}})+\alpha_{g}^{2}(\sum_{i=1}^{{m_{\mathcal{A}}}}\frac{k_{i}^{\ast}l_{i}^{\ast}(n_{i}^{2}-1)}{2n_{i}})$
$\displaystyle+|\alpha_{g}|^{2}\sum_{i=1}^{{m_{\mathcal{A}}}}N_{i}\Bigg{]}/{(\sum_{i=1}^{{m_{\mathcal{A}}}}{N}_{i}+|\alpha_{g}|^{2})^{2}},$

where ${N}_{i}=|k_{i}|^{2}\frac{n_{i}+1}{2}+|l_{i}|^{2}\frac{n_{i}-1}{2}$, and
${\tilde{N}}_{i}=|k_{i}|^{2}\frac{n_{i}+1}{2}-|l_{i}|^{2}\frac{n_{i}-1}{2}.$
By minimising this expression, we find that

$\frac{\mu_{{\mathcal{A}}}^{-}}{\mu_{{\mathcal{A}}}}\geqslant{1/2}.$ (9)

We thus show that photon subtraction decreases the purity of an arbitrary
state by at most a factor two. In the examples at the end of the manuscript,
we show that this lower bound is tight in the sense that there are states for
which ${\mu_{{\mathcal{A}}}^{-}}={\mu_{{\mathcal{A}}}}/2$. Furthermore, when
there is no entanglement in the initial Gaussian state, i.e. when $\mu_{\cal
A}=1$, we find that $\mu^{-}_{\cal A}=1$.

When the state for the global system is pure, we can immediately combine the
definition (1) with (9) to obtain that the entanglement increase is governed
by

$\Delta{\cal E}_{R}=\log\mu_{\cal A}-\log\mu^{-}_{\cal A}\leqslant\log 2.$
(10)

For mixed states, on the other hand, we use the definition of (3), which
implies that we must consider all possible decompositions of the state. For
each state in the decomposition we can use (2) to calculate
$\mathcal{E}_{R}(|\psi_{\lambda}\rangle)$ based on the Wigner function. This
construction naturally implies that
$\mathcal{E}_{R}(\hat{\rho}^{G})\leqslant\int{\rm d}\lambda
p^{\prime}(\lambda)\mathcal{E}_{R}(|\psi_{\lambda}\rangle)$, for any possible
decomposition of $\hat{\rho}$ in pure states. For Gaussian states this idea
was explicitly used to define Gaussian entanglement of formation M. M.Wolf  as
an upper bound to the entanglement of formation. At least for some classes of
Gaussian states it was later proven that Gaussian entanglement of formation is
the same as the general entanglement of formation P. Marian ; Y. Akbari-
Kourbolagh .

Figure 1: (Color online) (a) Schematic diagram of assigning all modes in a
global pure state to subsystems ${\cal A}$ and ${\cal B}$. The green dots
refer to non-Gaussian modes, and the small black circles refer to Gaussian
modes. (b) Thermal decomposition of subsystems $\mathcal{A}$, in which a
photon is subtracted (surrounded by a green circle) in the mode-$g$ after the
displacement operator $\mathcal{\hat{D}}$. (c) A combined operator $\hat{b}$
(surrounded by a green ellipses), which adds and subtracts a photon to and
from each thermal mode before the unitary operator $\hat{U}$. The operations
in (b) and (c) are different, but the Bogoliubov transform tells us that the
ultimate output state is strictly the same.

The Wigner function C. Weedbrook  of a Gaussian state $\hat{\rho}^{G}$ can be
written as

$\displaystyle W^{G}(\vec{\beta})$ $\displaystyle=$
$\displaystyle\frac{\mathrm{exp}[{-\frac{1}{2}(\vec{\beta}-\vec{\alpha}_{0})^{t}V^{-1}(\vec{\beta}-\vec{\alpha}_{0})}]}{{(2\pi)^{{m}}\sqrt{\mathrm{det}(V)}}},$
(11)

where $\vec{\beta}=({x}_{1},\cdots,{x}_{m},$
${p}_{1},\cdots,{p}_{m})^{t}\in\mathbb{R}^{2{m}}$ is a set of amplitude and
phase quadrature and $\xi$ is the corresponding displacement vector. And the
covariance matrix $V$ can be decomposed as $V=V_{p}+V_{c}$, where $V_{p}$ is
the covariance matrix of a pure state, and $V_{c}$ is some positive-definite
matrix that describes additional classical noise. This decomposition is
generally not unique, but we can use a constructive technique to prove that
this decomposition always exists M. Walschaers1 . We can think of the state
$\hat{\rho}^{G}$ as being generated by injecting a pure state $\hat{\rho}_{p}$
into a noisy Gaussian channel, so that the Gaussian state $\hat{\rho}^{G}$ can
be decomposed M. M.Wolf  as

$\displaystyle\hat{\rho}^{G}$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{2m}}\mathrm{d}^{2m}\vec{\alpha}\,\mathcal{\hat{D}}(\vec{\alpha})\hat{\rho}^{G}_{p}\mathcal{\hat{D}}^{{\dagger}}(\vec{\alpha})p_{c}^{G}(\vec{\alpha}),$
(12)

where $\hat{\rho}_{p}$ is a pure squeezed vacuum state with covariance matrix
$V_{p}$. We also introduce the probability distribution

$p_{c}^{G}(\vec{\alpha})=\frac{e^{{-\frac{1}{2}(\vec{\alpha}-\vec{\alpha}_{0})^{T}V_{c}^{-1}(\vec{\alpha}-\vec{\alpha}_{0})}}}{{(2\pi)^{{m}}\sqrt{\mathrm{det}(V_{c})}}},$
(13)

with $\vec{\alpha}_{0}$ the mean field of $\hat{\rho}^{G}$. This gives us an
explicit way to decompose any Gaussian state as a mixture of squeezed vacuum
states. This decomposition of the Gaussian state can be used to decompose the
associated photon-subtracted state as

$\displaystyle\hat{\rho}^{-}$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{2m}}\mathrm{d}^{2m}\vec{\alpha}\hat{\rho}_{p,\vec{\alpha}}^{-}p_{c}^{-}(\vec{\alpha}),$
(14)

where $\hat{\rho}_{p,\alpha}^{-}$ is the state that is obtained by subtracting
a single photon from
$\mathcal{\hat{D}}(\alpha)\hat{\rho}_{p}\mathcal{\hat{D}}^{{\dagger}}(\alpha)$,
and $p^{-}_{c}(\alpha)$ a probability distribution for which the explicit
expression is found in M. Walschaers1 .

The decomposition (14) can be used to derive an upper bound for the amount of
distilled entanglement, optimised over all possible decompositions of $V$.
First, we find that

${\cal E}_{R}(\hat{\rho}^{-})\leqslant\inf_{V_{p}\leqslant
V}\int_{\mathbb{R}^{2m}}\mathrm{d}^{2m}\vec{\alpha}p_{c}^{-}(\vec{\alpha})\mathcal{E}_{R}(\hat{\rho}_{p,\vec{\alpha}}^{-}).$
(15)

Because $\hat{\rho}_{p,\alpha}^{-}$ is a single-photon subtracted pure state,
we can use (10) to obtain that ${\cal
E}_{R}(\rho^{-}_{p,\vec{\alpha}})\leqslant\log 2+{\cal
E}_{R}(\mathcal{\hat{D}}(\vec{\alpha})\hat{\rho}_{p}\mathcal{\hat{D}}(\vec{\alpha})^{{\dagger}})$.
Furthermore, because entanglement is unchanged under local unitary
transformations, we find that ${\cal
E}_{R}(\mathcal{\hat{D}}(\vec{\alpha})\hat{\rho}_{p}\mathcal{\hat{D}}(\vec{\alpha})^{{\dagger}})={\cal
E}_{R}(\hat{\rho}_{p})$. Inserting all these elements in (15), we find that

${\cal E}_{R}(\hat{\rho}^{-})\leqslant\log 2+\inf_{V_{p}\leqslant
V}[\mathcal{E}_{R}(\hat{\rho}_{p})].$ (16)

The second term of the right side is by definition the Gaussian Rényi-2
entanglement of the mixed state $\hat{\rho}^{G}$ G. Adesso1 . This
entanglement measure is defined by minimising over all possible Gaussian
decompositions, but it is a priori not clear that there cannot be any
decomposition in non-Gaussian states that yields a lower entanglement. In
other words, we know that $\inf\limits_{V_{p}\leqslant
V}\mathcal{E}_{R}(\hat{\rho}_{p})\leqslant\mathcal{E}_{R}(\hat{\rho}^{G})$,
but it is not clear whether or not this is actually an equality.

Hence, we have shown that photon subtraction on any arbitrary Gaussian state
$\hat{\rho}^{G}$ can increase the Rényi-2 entanglement by at most $\log 2$ as
compared to its Gaussian Rényi-2 entanglement. A completely analogous
calculation for photon-added state results in the same bound even though
photon subtraction and photon addition tend to distill different amounts of
entanglement depending on the state.

Figure 2: (Color online) Two examples of multimode entangled pure states. The
optical modes are indicated by black circles, in which the color code
indicates the increase (${\Delta\cal E}_{R}$) in Rényi-2 entanglement between
the corresponding mode and the rest of the system. The mode $g$ with one
photon subtracted is indicated by a downward green arrow and its neighboring
mode, mode $g^{\prime}$, are marked with a green circle. (a) is a parallel
optical network with 10 modes, and all these modes are in the initial vacuum
sates. $S$ indicates the two-mode squeezer operator with parameter $r=1$. $D$
is the displacement operator, which makes mode $g$ produce a displacement of
$\alpha_{g}=0.5$. (b) is a $3\times 3$ graph state. All these modes are in the
initial squeezed vacuum with $10$ dB of squeezing. And the mode $g$ also has a
displacement of $\alpha_{g}=0.5$.

To get a better understanding of the quality of the bound and its physical
interpretation, we will consider a range of examples. Because there is no
exact method to calculate the entanglement for non-Gaussian mixed states, we
restrict ourselves to pure states where we can use (10) to calculate the
entanglement increase. We consider the following two typical entangled pure
states as examples:

  1. 1.

We subtract a photon from a linear entangled network, which is composed of
several two-mode squeezer gates, as shown in Fig. 2(a). All the modes are
initially in vacuum states $\hat{\rho}_{\mathrm{o}}$, then the output state
can be expressed as
$\hat{\rho}=\hat{\mathcal{D}}\hat{U}\hat{\rho}_{\mathrm{o}}\hat{U}^{{\dagger}}{\hat{\mathcal{D}}^{{\dagger}}}$,
where $\hat{U}=\hat{S}_{m-1}\cdots\hat{S}_{1}$ with
$\hat{S}_{i}=\mathrm{exp}[r(\hat{a}_{i}\hat{a}_{i+1}-\hat{a}_{i}^{{\dagger}}\hat{a}_{i+1}^{{\dagger}})/2]$.
To test the effect of a mean field, we add the operator $\hat{\mathcal{D}}$
that adds a displacement $(\alpha_{g})$ to the mode-g. Thus the covariance
matrix of the state $\hat{\rho}$ is $V_{\mathrm{network}}=C^{t}C$, where $C$
is a symplectic matrix that implements the multimode squeezing operation on
the quadratures $\hat{U}^{{\dagger}}\hat{\beta}\hat{U}=C\hat{\beta}$.

  2. 2.

We subtract a photon from a graph state, as shown in Fig. 2(b), which is the
backbone of measurement-based CV quantum computation. We start out from a set
of independent squeezed modes, with a joint covariance matrix given by
$V_{0}=\mathrm{diag}(s,\cdots,s,s^{-1},\cdots,s^{-1})$. Every edge of the
graph corresponds to a $C_{Z}$ gate, which turn those non-correlated modes
into a graph state. On the level of covariance matrices, the action of the set
of $C_{Z}$ gates is represented by a symplectic transformation $C$, applied to
$V_{0}$ M. Walschaers3 : $V_{\mathrm{graph}}=C^{t}V_{0}C$, with
$C=\left(\begin{array}[]{cc}\mathds{1}&\mathbb{A}\\\ 0&\mathds{1}\\\
\end{array}\right)$, where $\mathbb{A}$ is $m\times m$ adjacency matrix.

We can use methods purely based on phase space representations M. Walschaers1
; M. Walschaers2  to describe the Wigner function of a reduced photon
subtracted state for any subsystem ${\cal A}$, which leads to the general
expression

$\displaystyle W^{-}_{\cal A}(\vec{\beta}_{\cal A})=$
$\displaystyle\frac{\lVert XV^{-1}_{\cal A}(\vec{\beta}_{\cal
A}-\vec{\alpha}_{\cal A})-\vec{\alpha}_{g}\rVert^{2}+{\rm
tr}(V_{g}-XV^{-1}_{\cal A}X^{t})-2}{\lVert\vec{\alpha}_{g}\rVert^{2}+{\rm
tr}(V_{g})-2}$ (17) $\displaystyle\times W^{G}_{\cal A}(\vec{\beta}_{\cal
A}),$

where $\vec{\beta}_{\cal A}\in\mathbb{R}^{2{m_{\cal A}}}$ and
$X=G^{t}(V-\mathds{1})A$. Here, $G$ is a $2m\times 2$ matrix where the two
columns are the basis vectors $\vec{g}^{(x)}$ and $\vec{g}^{(p)}$ associated
with the two phase space axes of mode $g$:

$G=\begin{pmatrix}\mid&\mid\\\ \vec{g}^{(x)}&\vec{g}^{(p)}\\\
\mid&\mid\end{pmatrix}.$ (18)

Analogously, $A$ is a $2m\times 2m_{\cal A}$ matrix where the columns are the
symplectic basis vectors $\vec{a}^{(x)}_{1},\dots,\vec{a}^{(x)}_{m_{\cal
A}},\vec{a}^{(p)}_{1},\dots,\vec{a}^{(p)}_{m_{\cal A}}$ that generate the
phase space of subsystem $\cal A$:

$A=\begin{pmatrix}\mid&&\mid&\mid&&\mid\\\
\vec{a}_{1}^{(x)}&\cdots&\vec{a}_{m_{\cal
A}}^{(x)}&\vec{a}_{1}^{(p)}&\cdots&\vec{a}_{m_{\cal A}}^{(p)}\\\
\mid&&\mid&\mid&&\mid\end{pmatrix}.$ (19)

In (17), we then obtain the matrices $V_{\cal A}=A^{t}VA$ and $V_{g}=G^{t}VG$,
and the displacement vectors $\vec{\alpha}_{\cal A}=A^{t}\vec{\alpha}$ and
$\vec{\alpha}_{g}=G^{t}\vec{\alpha}$. We can now use standard techniques of
Gaussian integrals to evaluate the purity (2) from $W^{-}_{\cal
A}(\vec{\beta}_{\cal A})$ for an arbitrary bipartition ${\cal A}$. This in
turn allows to directly evaluate the entanglement increase (10).

In Fig. 2, each circle represents a mode in a multimode quantum state, and the
color code inside the circle indicate the increase in entanglement between the
mode (subsystem $\mathcal{A}$) and the rest system (subsystem $\mathcal{B}$).
When a photon is subtracted from the mode $g$ ( marked by the green arrow) in
the pure multimode systems, non-Gaussian features will spread from the mode
$g$ to other modes, whose entanglement can also increase. In Fig. 3(a) and
Fig. 3(b), we respectively show the increase of entanglement (${\Delta\cal
E}_{R}$) related to mode $g$ and mode $g^{\prime}$ (representing the
neighboring mode of $g$ and marked by a green circle in Fig. 2). As shown by
the blue and pink curves, ${\Delta\cal E}_{R}$ increases as the degree of
squeezing increases until it approaches a maximum. The maximum values are
usually different for different initial Gaussian states, but they are all less
than $\log 2$. In Fig. 3(c), each dots represents a different bi-partition for
the example of Fig. 2(a), which we group depending on the number of modes in
subsystem $\mathcal{A}$ (similar data for the example of Fig. 2(b) are shown
in Supplementary Material). The entanglement increase varies depending on the
modes that are comprised in each subsystem, but it is never greater than $\log
2$.

Finally, we aim to provide an intuition for the value $\log 2$ of entanglement
increase. To this goal, we restrict to a two-mode system prepared in a non-
displaced pure state $\left\lvert\psi\right\rangle=\hat{U}(\left\lvert
0\right\rangle\otimes\left\lvert 0\right\rangle)$, where $\hat{U}$ is a
Bogoliubov transformation. When we subtract a photon in one of the two modes,
we can use (6) to write the photon-subtracted state
$\left\lvert\psi^{-}\right\rangle=\hat{U}\hat{b}(\left\lvert
0\right\rangle\otimes\left\lvert 0\right\rangle)/{\cal N}$ (with ${\cal N}$ a
normalisation factor). We recall that
$\hat{b}=\hat{U}^{{\dagger}}\hat{a}_{g}\hat{U}=\vec{k}\cdot\vec{\hat{a}}^{{\dagger}}+\vec{l}\cdot\vec{\hat{a}}$,
which implies that

$\left\lvert\psi^{-}\right\rangle=\frac{1}{\cal
N}\hat{U}\left(k_{1}\left\lvert 1\right\rangle\otimes\left\lvert
0\right\rangle+k_{2}\left\lvert 0\right\rangle\otimes\left\lvert
1\right\rangle\right).$ (20)

In particular in the limit where the squeezing induced by $\hat{U}$ is small,
we find that $\left\lvert\psi\right\rangle$ resembles the vacuum state with
negligible entanglement. On the other hand, $\left\lvert\psi^{-}\right\rangle$
is close to a single photon that passed through a beam splitter. The maximal
entanglement in this low squeezing limit is achieved for the Bell state
$\left\lvert\psi^{-}\right\rangle\approx(\left\lvert
1\right\rangle\otimes\left\lvert 0\right\rangle\pm\left\lvert
0\right\rangle\otimes\left\lvert 1\right\rangle)/\sqrt{2}$, which has a
Rényi-2 entanglement ${\cal E}_{R}(\left\lvert\varphi\right\rangle)=\log 2$.
This corresponds exactly to the maximal entanglement increase and thus we can
interpret our findings as a generalisation of this phenomenon to mixed states
of arbitrary many modes with arbitrary amounts of squeezing.

Figure 3: (Color online) (a) and (b) are the increase of entanglement
${\Delta\cal E}_{R}$ as the degree of squeezing changes, corresponding to the
network and graph states in Fig. 2, respectively. The solid and dash curves
correspond to the photon-subtracted mode (mode-g) and its neighboring mode
(mode-$g^{\prime}$), respectively. The blue (cyan) and the red (pink) curves
correspond to the small displacement $\alpha_{g}=0.5$ and $\alpha_{g}=0$,
respectively. The single-mode Wigner functions are shown for the curves
pointed by the green thin arrow. (c) Each scattered dot represents a different
bi-partition with $m_{\cal A}$ modes in subsystem $\cal A$. The blue and the
red scattered dots correspond to the first example in Fig. 2(a).

In conclusion, we strictly prove that single-photon subtraction can reduce the
purity of any Gaussian state by at most a factor $1/2$. This means that for
any Gaussian state the amount of Rényi-2 entanglement that can be distilled by
single-photon subtraction is at most $\log 2$. This result holds for both pure
and mixed states. A fully equivalent treatment of single-photon addition
yields the same result (see Supplementary Material).

We evaluate the actual Rényi-2 entanglement of two examples based on the
Wigner function and show that the bound is tight. The results of this work
mark the limits of the non-Gaussian features induced by photon subtraction (or
addition), which are in many ways complementary to other studies on induced
Wigner negativity M. Walschaers PRL ; M. Walschaers PRXQ .

Acknowledgements— K.Z. and J.J. acknowledge financial support through the
National Natural Science Foundation of China (11874155, 91436211); Basic
Research Project of Shanghai Science and Technology Commission (20JC1416100);
Natural Science Foundation of Shanghai (17ZR1442900); Minhang Leading Talents
(201971); the 111 project (B12024). N.T. and M.W. received funding from the
European Union’s Horizon 2020 research and innovation programme under grant
agreement No 899587.

## References

  * (1) A. Einstein, B. Podolsky, and N. Rosen, Phys. Rev. A 47, 777 (1935).
  * (2) V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev. Lett. 96, 010401 (2006).
  * (3) C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, Phys. Rev. Lett. 70, 1895 (1993).
  * (4) P. W. Shor, Phys. Rev. A 52, R2493 (1995).
  * (5) S. Steinlechner, J. Bauchrowitz, M. Meinders, H. Müller-Ebhardt, K. Danzmann, R. Schnabel, Nat. Photonics 7, 626 (2013).
  * (6) B. Schumacher, and M. A. Nielsen, Phys. Rev. A 54, 2629 (1996).
  * (7) S. Lloyd, Phys. Rev. A 55, 1613 (1997).
  * (8) R. Raussendorf, and H. J. Briegel, Phys. Rev. Lett. 86, 5188 (2001).
  * (9) E. Knill, R. Laflamme, and G. J. Milburn, Nature 409, 46 (2001).
  * (10) S. L. Braunstein and P. van Loock, Rev. Mod. Phys. 77, 513 (2005).
  * (11) X. Su, Y. Zhao, S. Hao, X. Jia, C. Xie, and K. Peng, Opt. Lett. 37, 5178 (2012).
  * (12) S. Armstrong, J.-F. Morizur, J. Janousek, B. Hage, N. Treps, P. K. Lam, and H.-A. Bachor, Nat. Commun. 3, 1206 (2012).
  * (13) M. Chen, N. C. Menicucci, and O. Pfister, Phys. Rev. Lett. 112, 120505 (2014).
  * (14) J. Roslund, R. M. de Araújo, S. Jiangl, C. Fabre, and N. Treps, Nat. Photonics 8, 109 (2014).
  * (15) S. Gerke, J. Sperling, W. Vogel, Y. Cai, J. Roslund, N. Treps, and C. Fabre, Phys. Rev. Lett. 114, 050501 (2015).
  * (16) J.-i. Yoshikawa, S. Yokoyama, T. Kaji, C. Sornphiphatphong, Y. Shiozawa, K. Makino, and A. Furusawa, APL Photonics 1, 060801 (2016).
  * (17) K. Zhang, W. Wang, S. Liu, X. Pan, J. Du, Y. Lou, S. Yu, S. Lv, N. Treps, C. Fabre, and J. Jing, Phys. Rev. Lett. 125, 140501 (2020).
  * (18) W. Wang, K. Zhang, and J. Jing, Phys. Rev. Lett. 125, 140501 (2020).
  * (19) S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto, Phys. Rev. Lett. 88, 097904 (2002).
  * (20) A. Mari, and J. Eisert, Phys. Rev. Lett. 109, 230503 (2012
  * (21) S. Rahimi-Keshari, T. C. Ralph, and C. M. Caves, Phys. Rev. X 6, 021039 (2016).
  * (22) A. Zavatta, V. Parigi, and M. Bellini, Phys. Rev. A 75, 052106 (2007).
  * (23) A. Ourjoumtsev, A. Dantan, R. Tualle-Brouri, and P. Grangier, Phys. Rev. Lett. 98, 030502 (2007).
  * (24) V. Parigi, A. Zavatta, M. Kim, and M. Bellini, Science 317, 1890 (2007).
  * (25) Y.-S. Ra, A. Dufour, M. Walschaers, C. Jacquard, T. Michel, C. Fabre, and N. Treps, Nat. Phys. 16, 144 - 147 (2020).
  * (26) J. Wenger, R. Tualle-Brouri, and P. Grangier, Phys. Rev. Lett. 92, 153601 (2004).
  * (27) Yang Yang, and F.-L. Li, Phys. Rev. A 80, 022315 (2009).
  * (28) C. Navarrete-Benlloch, R. García-Patrón, J. H. Shapiro, and N. J. Cerf, Phys. Rev. A 86, 012328 (2012).
  * (29) M. Walschaers, S. Sarkar, V. Parigi, and N. Treps, Phys. Rev. Lett. 121, 220501 (2018).
  * (30) M. Walschaers, C. Fabre, V. Parigi, and N. Treps, Phys. Rev. Lett. 119, 183601 (2017).
  * (31) M. Walschaers, C. Fabre, V. Parigi, and N. Treps, Phys. Rev. A. 96, 053835 (2017).
  * (32) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * (33) C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, Phys. Rev. A 54, 3824 (1996).
  * (34) G. Tóth, T. Moroder, and O. Gühne, Phys. Rev. Lett. 114, 160501 (2015).
  * (35) J. S. Kim, and B. C. Sanders, J. Phys. A: Math. Theor. 43, 445305 (2010).
  * (36) G. Adesso, D. Girolami, and A. Serafini, Phys. Rev. Lett. 109, 190502 (2012).
  * (37) C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Rev. Mod. Phys 84, 621 (2012).
  * (38) M. M.Wolf, G. Giedke, O. Krüger, R. F.Werner, and J. I. Cirac, Phys. Rev. A 69, 052320 (2004).
  * (39) P. Marian, and T. A. Marian, Phys. Rev. Lett. 101, 220403 (2008).
  * (40) Y. Akbari-Kourbolagh, and H. Alijanzadeh-Boura, Quantum Information Processing 14, 4179 (2015).
  * (41) M. Walschaers, N. Treps, Phys. Rev. Lett. 124, 150501 (2020).
  * (42) M. Walschaers, V. Parigi, N. Treps, PRX Quantum 1, 020305 (2020).

## Supplementary Material

##  I $\mathbf{A}$. Decomposition of single-photon subtraction operation

In CV quantum optics, for an arbitrary $m$-mode Gaussian state
$\hat{\rho}_{0}$, we can write a thermal decomposition as

$\displaystyle\hat{\rho}_{0}=\mathcal{\hat{D}}\hat{U}\hat{\rho}_{M}\hat{U}^{{\dagger}}\mathcal{\hat{D}}^{{\dagger}},$
(21)

where $\hat{U}$ is an arbitrary canonical unitary and
$\mathcal{\hat{D}}=\prod\limits_{i=1}^{m}\hat{D}(\alpha_{i})$ is displacement
operator, which displace the vacuum state to generate coherent states
$|\alpha_{i}\rangle=\hat{D}(\alpha_{i})|0\rangle$. Furthermore,
$\hat{\rho}_{M}=\bigotimes\limits_{i=1}^{m}\hat{\rho}_{i}$, where single- mode
thermal states $\hat{\rho}_{i}=(0,V_{i})$ can be fully characterised by a
covariance matrix $V_{i}=\mathrm{diag}(n_{i},n_{i})$, , where $n_{i}$ is the
fraction of thermal noise compared to shot noise. Then the purity of
$\hat{\rho}_{0}(\bar{\mathrm{x}},V)$ can be expressed as
$\mu_{0}=\prod\limits_{i=1}^{m}\mathrm{tr}(\hat{\rho}_{i}^{2})=\prod\limits_{i=1}^{m}\frac{1}{n_{i}}$.

When subtracting a photon to the mode-$\mathrm{g}$, it corresponds to a
unitary Bogoliubov transformation

$\displaystyle\hat{U}^{{\dagger}}\mathcal{\hat{D}}^{{\dagger}}\hat{a}_{\mathrm{g}}\mathcal{\hat{D}}\hat{U}$
$\displaystyle=$
$\displaystyle\hat{b}+\alpha_{\mathrm{g}}^{\ast},\quad\mathrm{with}\quad\hat{b}=\vec{k}\cdot\vec{\hat{a}}^{{\dagger}}+\vec{l}\cdot\vec{\hat{a}},$
(22)

where
$\vec{\hat{a}}^{{\dagger}}=(\hat{a}_{1}^{{\dagger}},\cdots,\hat{a}_{m}^{{\dagger}})^{T}$
and $\vec{\hat{a}}=(\hat{a}_{1},\cdots,\hat{a}_{m})^{T}$ . A set of $1\times
m$ complex matrices $\vec{k}=(k_{1},\cdots,k_{m})$ and
$\vec{l}=(l_{1},\cdots,l_{m})$ corresponds an canonical unitary $\hat{U}$.
Thus, we convert the single-photon subtraction operation on the
mode-$\mathrm{g}$ into the addition and subtraction operation on each thermal
mode. Thus the single-photon subtracted state can then be written as

$\displaystyle\frac{\hat{a}_{\mathrm{g}}\hat{\rho}_{0}\hat{a}^{{\dagger}}_{\mathrm{g}}}{\mathrm{tr}(\hat{a}_{\mathrm{g}}\hat{\rho}_{0}\hat{a}^{{\dagger}}_{\mathrm{g}})}=\frac{\mathcal{\hat{D}}\hat{U}(\hat{b}+\alpha_{\mathrm{g}})\hat{\rho}_{M}(\hat{b}^{{\dagger}}+\alpha^{\ast}_{\mathrm{g}})\hat{U}^{{\dagger}}\mathcal{\hat{D}}^{{\dagger}}}{\mathcal{N}},$
(23)

where the normalization factor
$\mathcal{N}=\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}})+|\alpha_{\mathrm{g}}|^{2}$.
Due to $\mathcal{\hat{D}}\hat{U}$ will not change the purity, so the purity of
an arbitrary single-photon subtracted state can be expressed as

$\displaystyle\mu{\mathcal{N}^{2}}$ $\displaystyle=$
$\displaystyle\mathrm{tr}([\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}}+|\alpha_{\mathrm{g}}|^{2}\hat{\rho}_{M}+\alpha_{\mathrm{g}}\hat{\rho}_{M}\hat{b}^{{\dagger}}+\alpha_{\mathrm{g}}^{\ast}\hat{b}\hat{\rho}_{M}]^{2})$
(24) $\displaystyle=$
$\displaystyle\mathrm{tr}((\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}}+|\alpha_{\mathrm{g}}|^{2}\hat{\rho}_{M})^{2}+(\alpha_{\mathrm{g}}\hat{\rho}_{M}\hat{b}^{{\dagger}}+\alpha^{\ast}_{\mathrm{g}}\hat{b}\hat{\rho}_{M})^{2})$
$\displaystyle=$
$\displaystyle\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}})+|\alpha_{\mathrm{g}}|^{4}\mathrm{tr}(\hat{\rho}_{M}^{2})+{\alpha_{\mathrm{g}}}^{2}\mathrm{tr}(\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{\rho}_{M}\hat{b}^{{\dagger}})+{\alpha^{\ast}_{\mathrm{g}}}^{2}\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{b}\hat{\rho}_{M})+$
$\displaystyle|\alpha_{\mathrm{g}}|^{2}[\mathrm{tr}(\hat{\rho}_{M}\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}})+\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{\rho}_{M})]+|\alpha_{\mathrm{g}}|^{2}[\mathrm{tr}(\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{b}\hat{\rho}_{M})+\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{\rho}_{M}\hat{b}^{{\dagger}})].$

where we have used that
$\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{\rho}_{M}\hat{b}^{{\dagger}})=\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{b}\hat{\rho}_{M})=\mathrm{tr}(\hat{\rho}_{M}\hat{\rho}_{M}\hat{b}^{{\dagger}})=\mathrm{tr}(\hat{\rho}_{M}\hat{b}\hat{\rho}_{M})=0$,
since $\hat{\rho}_{i}$ is a diagonal matrix, and $\hat{a}_{i}^{{\dagger}}$
$(\hat{a}_{i})$ is off-diagonal matrices. Similarly, for single-photon
addition, the calculation is equivalent the same expression just by swapping
the values of $\vec{k}$ and $\vec{l}$ and the values of
$\alpha_{\mathrm{g}}^{\ast}$ with $\alpha_{\mathrm{g}}$.

For the same Gaussian state, different form Eq.(21), we can also take the
following decomposition

$\displaystyle\hat{\rho}_{0}=\hat{U}\mathcal{\hat{D}}(\bar{\mathrm{x}})\hat{\rho}_{M}\mathcal{\hat{D}}(\bar{\mathrm{x}})^{{\dagger}}\hat{U}^{{\dagger}}.$
(25)

In this case, we just need to to replace $\alpha_{\mathrm{g}}^{\ast}$ with
$\sum_{i=1}^{m}k_{i}\alpha^{\ast}+l_{i}\alpha$ in Eq.(24), which is still
applicable, regardless of the decomposition form.

##  II $\mathbf{B}$. Generalized expression of relative purity

###  II.1 1\. When $\alpha_{\mathrm{g}}=0$

When $\alpha_{\mathrm{g}}=0$, only the first term of (24) is retained.

Let us take $m=1$, so that (24) becomes

$\displaystyle\mu_{1}{N}_{1}^{2}=$
$\displaystyle\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}})$
(26) $\displaystyle=$
$\displaystyle|k_{1}|^{4}\mathrm{tr}(\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1})+|l_{1}|^{4}\mathrm{tr}(\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}}\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}})+k_{1}^{2}l_{1}^{\ast
2}\mathrm{tr}(\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}}\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}})+k_{1}^{\ast
2}l_{1}^{2}\mathrm{tr}(\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1})$
(27)
$\displaystyle+|k_{1}|^{2}|l_{1}|^{2}[\mathrm{tr}(\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}})+\mathrm{tr}(\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}}\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1})]+|k_{1}|^{2}|l_{1}|^{2}[\mathrm{tr}(\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}}\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1})+\mathrm{tr}(\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}})]$
$\displaystyle+|k_{1}|^{2}k_{1}l_{1}^{\ast}[\mathrm{tr}(\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}})+\mathrm{tr}(\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}}\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1})]+|k_{1}|^{2}k_{1}^{\ast}l_{1}[\mathrm{tr}(\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1})+\mathrm{tr}(\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1})]$
$\displaystyle+|l_{1}|^{2}k_{1}l_{1}^{\ast}[\mathrm{tr}(\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}}\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}})+\mathrm{tr}(\hat{a}_{1}^{{\dagger}}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}}\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}})]+|l_{1}|^{2}k_{1}^{\ast}l_{1}[\mathrm{tr}(\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}}\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1})+\mathrm{tr}(\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}\hat{a}_{1}\hat{\rho}_{1}\hat{a}_{1}^{{\dagger}})].$

Remove the zero traces and replace the subscript 1 with $i$, we get any
single-mode purity

$\displaystyle\mu_{i}{N}_{i}^{2}=$
$\displaystyle|k_{i}|^{4}\mathrm{tr}(\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i}\hat{a}_{i}\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i}\hat{a}_{i})+|l_{i}|^{4}\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}}\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}})+2|k_{i}|^{2}|l_{i}|^{2}\mathrm{tr}(\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i}\hat{a}_{i}\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}})+2|k_{i}|^{2}|l_{i}|^{2}\mathrm{tr}(\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}}\hat{a}_{i}\hat{\rho}_{1}\hat{a}_{i}).$
(28)

When we generalise this to a state with an arbitrary number of modes $m$, we
find that the purity can be expressed as

$\displaystyle\mu_{m}(\sum_{i=1}^{m}{N}_{i})^{2}=$
$\displaystyle\sum_{i=1}^{m}\mu_{i}{N}_{i}^{2}\cdot
n_{i}\cdot\prod_{i=1}^{m}\frac{1}{n_{i}}+4\cdot[|\sum_{i=1}^{m}[k_{i}l_{i}\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i})\cdot
n_{i}|^{2}-\sum_{i=1}^{m}|k_{i}l_{i}\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i})\cdot
n_{i}|^{2}]\cdot\prod_{i=1}^{m}\frac{1}{n_{i}}+$ (29)
$\displaystyle[(\sum_{i=1}^{m}[|k_{i}|^{2}\mathrm{tr}(\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i}^{2}\hat{a}_{i})+|l_{i}|^{2}\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}^{2}\hat{a}_{i}^{{\dagger}})]\cdot
n_{i})^{2}-\sum_{i=1}^{m}[|k_{i}|^{2}\mathrm{tr}(\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i}^{2}\hat{a}_{i})+|l_{i}|^{2}\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}^{2}\hat{a}_{i}^{{\dagger}})]^{2}\cdot
n_{i}^{2}]\cdot\prod_{i=1}^{m}\frac{1}{n_{i}}+$
$\displaystyle[(\sum_{i=1}^{m}[(|k_{i}|^{2}+|l_{i}|^{2})\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i})]\cdot
n_{i})^{2}-\sum_{i=1}^{m}[(|k_{i}|^{2}+|l_{i}|^{2})\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i})]^{2}\cdot
n_{i}^{2}]\cdot\prod_{i=1}^{m}\frac{1}{n_{i}},$

where
$\mathcal{N}_{i}=|k_{i}|^{2}\mathrm{tr}(\hat{a}_{i}^{{\dagger}}\rho_{i}\hat{a}_{i})+|l_{i}|^{2}\mathrm{tr}(\hat{a}_{i}\rho_{i}\hat{a}_{i}^{{\dagger}})$.
Therefore, the traces in (29) are related to the single-mode photon-added and
-subtracted thermal states, which are the key to deriving the expression of
final purity. Fortunately, Winger function provides us with an easy way to
calculate these traces. For a single-mode Gaussian thermal state
$\hat{\rho}_{i}$, the Wigner function is

$\displaystyle w_{i}(q_{i},p_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi
n_{i}}\mathrm{exp}[-\frac{q_{i}^{2}+p_{i}^{2}}{2n_{i}}],$ (30)

where $q_{i}$ and $p_{i}$ are amplitude and phase quadrature. Then the Wigner
functions due to photon addition or subtraction are

$\displaystyle w_{(i,+1)}(q_{i},p_{i})$ $\displaystyle=$
$\displaystyle(\frac{n_{i}+1}{2n_{i}}(q_{i}^{2}+p_{i}^{2})-\frac{1}{n_{i}})w_{i}(q_{i},p_{i}),\quad\mathrm{and}\quad
w_{(i,-1)}(q_{i},p_{i})=(\frac{n_{i}-1}{2n_{i}}(q_{i}^{2}+p_{i}^{2})+\frac{1}{n_{i}})w_{i}(q_{i},p_{i}).$
(31)

Thus, we can calculate the traces by the above Wigner functions:

$\displaystyle\mathrm{tr}(\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}})$
$\displaystyle=$
$\displaystyle\int\frac{q_{i}^{2}+p_{i}^{2}}{4}w_{i}\mathrm{d}q_{i}\mathrm{d}p_{i}-1=\frac{n_{i}-1}{2},\quad\mathrm{and}\quad\mathrm{tr}(\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a})=\int\frac{q_{i}^{2}+p_{i}^{2}}{4}w_{i}\mathrm{d}q_{i}\mathrm{d}p_{i}-1+1=\frac{n_{i}+1}{2},$
(32)

which are the quantities related to the average-photon in the thermal state.
Similarly, we can get that

$\displaystyle\mathrm{tr}(\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}}\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}})$
$\displaystyle=$
$\displaystyle\frac{\mathrm{tr}(\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}}\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}})}{\mathrm{tr}(\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}})^{2}}\cdot\mathrm{tr}(\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}})^{2}=\frac{(1+n_{i}^{2})}{2n_{i}^{3}}(\frac{n_{i}-1}{2})^{2},\quad\mathrm{and}\quad\mathrm{tr}(\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a}\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a})=\frac{(1+n_{i}^{2})}{2n_{i}^{3}}(\frac{n_{i}+1}{2})^{2},$
(33)

which are quantities related to the purity of the photon-added and -subtracted
states.

$\displaystyle\mathrm{tr}(\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a}\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}})$
$\displaystyle=$
$\displaystyle\mathrm{tr}(\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}}\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a})=\mathrm{tr}(\hat{\rho}_{i}\hat{a}^{{\dagger}}\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a}\hat{a})=\frac{(n_{i}^{2}-1)^{2}}{8n_{i}^{3}},$
(34)

which can be obtained by photon-added and -subtracted Wigner functions through
the following integration

$\displaystyle\mathrm{tr}(\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a}\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}})$
$\displaystyle=$ $\displaystyle(4\pi\int w_{(i,+1)}\cdot
w_{(i,-1)}\mathrm{d}q_{i}\mathrm{d}p_{i})\cdot(\int\frac{q_{i}^{2}+p_{i}^{2}}{4}w_{i}\mathrm{d}q_{i}\mathrm{d}p_{i}-1)\cdot(\int\frac{q_{i}^{2}+p_{i}^{2}}{4}w_{i}\mathrm{d}q_{i}\mathrm{d}p_{i})=\frac{(n_{i}^{2}-1)^{2}}{8n_{i}^{3}}.$
(35)

With the help of Wigner function of single-photon added or subtracted and
Gaussian thermal state Wigner function, we can get

$\displaystyle\mathrm{tr}(\hat{\rho}_{i}\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a})$
$\displaystyle=$
$\displaystyle\mathrm{tr}(\hat{\rho}_{i}\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}})=(4\pi\int
w_{i}\cdot
w_{(i,-1)}\mathrm{d}q_{i}\mathrm{d}p_{i})\cdot(\int\frac{q_{i}^{2}+p_{i}^{2}}{4}w_{i}\mathrm{d}q_{i}\mathrm{d}p_{i}-1)=\frac{n_{i}+1}{2n_{i}^{2}}(\frac{n_{i}-1}{2}).$

Due to the noncommutative nature of creation and annihilation operators, we
have

$\displaystyle\mathrm{tr}(\hat{a}\hat{\rho}_{i}\hat{a}\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a}^{{\dagger}})=\mathrm{tr}(\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}}\hat{a}\hat{\rho}_{i}\hat{a}^{{\dagger}})+\mathrm{tr}(\hat{\rho}_{i}^{2}\hat{a}^{{\dagger}}\hat{a})=\mathrm{tr}(\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a}\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a})-\mathrm{tr}(\hat{\rho}_{i}^{2}\hat{a}^{{\dagger}}\hat{a})-\mathrm{tr}(\hat{\rho}_{i}^{2}),$

then we can deduce the following traces

$\displaystyle\mathrm{tr}(\hat{\rho}_{i}^{2}\hat{a}^{{\dagger}}\hat{a})$
$\displaystyle=$
$\displaystyle\frac{(n_{i}-1)^{2}}{4n_{i}^{2}},\quad\mathrm{tr}(\hat{\rho}_{i}^{2}\hat{a}\hat{a}^{{\dagger}})=\frac{(n_{i}+1)^{2}}{4n_{i}^{2}}\quad\mathrm{and}\quad\mathrm{tr}(\hat{a}\hat{\rho}_{i}\hat{a}\hat{a}^{{\dagger}}\hat{\rho}_{i}\hat{a}^{{\dagger}})=\frac{(n_{i}^{2}-1)^{2}}{8n_{i}^{3}}.$
(36)

By bringing the results of these traces into (29), we get

$\displaystyle{\mu_{m}^{-}}=$
$\displaystyle[\frac{1}{2}(\sum_{i=1}^{m}{N}_{i})^{2}+\frac{1}{2}(\sum_{i=1}^{m}\frac{\tilde{{N}}_{i}}{n_{i}})^{2}+|\sum_{i=1}^{m}k_{i}l_{i}\frac{n_{i}^{2}-1}{2n_{i}}|^{2}]{\mu_{0}}/{(\sum_{i=1}^{m}{N}_{i})^{2}}.$

where ${N}_{i}=|k_{i}|^{2}\frac{n_{i}+1}{2}+|l_{i}|^{2}\frac{n_{i}-1}{2}$ and
${\tilde{N}}_{i}=|k_{i}|^{2}\frac{n_{i}+1}{2}-|l_{i}|^{2}\frac{n_{i}-1}{2}$.
Hence, we can immediately know the relative purity is

$\displaystyle\frac{\mu_{m}^{-}}{\mu_{0}}=$
$\displaystyle\frac{1}{2}+[\frac{1}{2}(\sum_{i=1}^{m}\frac{\tilde{{N}}_{i}}{n_{i}})^{2}+|\sum_{i=1}^{m}k_{i}l_{i}\frac{n_{i}^{2}-1}{2n_{i}}|^{2}]/{(\sum_{i=1}^{m}{N}_{i})^{2}}\geqslant
1/2.$

###  II.2 2\. When $\alpha_{\mathrm{g}}\neq 0$

According to similar approach above, we can derive the following traces

$\displaystyle|\alpha_{\mathrm{g}}|^{4}\mathrm{tr}(\hat{\rho}_{M}^{2})$
$\displaystyle=$ $\displaystyle|\alpha_{\mathrm{g}}|^{4}\mu_{0},$ (37)
$\displaystyle{\alpha_{\mathrm{g}}}^{2}\mathrm{tr}(\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{\rho}_{M}\hat{b}^{{\dagger}})$
$\displaystyle=$ $\displaystyle
2{\alpha_{\mathrm{g}}}^{2}(\sum_{i=1}^{m}k_{i}^{\ast}l_{i}^{\ast}\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i})n_{i})\mu_{0}={\alpha_{\mathrm{g}}}^{2}(\sum_{i=1}^{m}\frac{k_{i}^{\ast}l_{i}^{\ast}(n_{i}^{2}-1)}{2n_{i}})\mu_{0},$
(38)
$\displaystyle{\alpha^{\ast}_{\mathrm{g}}}^{2}\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{b}\hat{\rho}_{M})$
$\displaystyle=$ $\displaystyle
2{\alpha^{\ast}_{\mathrm{g}}}^{2}(\sum_{i=1}^{m}k_{i}l_{i}\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i})n_{i})\mu_{0}={\alpha^{\ast}_{\mathrm{g}}}^{2}(\sum_{i=1}^{m}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}})\mu_{0}$
(39)

and

$\displaystyle|\alpha_{\mathrm{g}}|^{2}[\mathrm{tr}(\hat{\rho}_{M}\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}})+\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{\rho}_{M})+\mathrm{tr}(\hat{\rho}_{M}\hat{b}^{{\dagger}}\hat{b}\hat{\rho}_{M})+\mathrm{tr}(\hat{b}\hat{\rho}_{M}\hat{\rho}_{M}\hat{b}^{{\dagger}})]$
$\displaystyle=2|\alpha_{\mathrm{g}}|^{2}[\sum_{i}^{m}(|k_{i}|^{2}+|l_{i}|^{2})\mathrm{tr}(\hat{a}_{i}\hat{\rho}_{i}\hat{a}_{i}^{{\dagger}}\hat{\rho}_{i})n_{i}+\sum_{i}^{m}(|k_{i}|^{2}\mathrm{tr}(\hat{\rho}_{i}^{2}\hat{a}_{i}\hat{a}_{i}^{{\dagger}})n_{i}+|l_{i}|^{2}\mathrm{tr}(\hat{\rho}_{i}^{2}\hat{a}_{i}^{{\dagger}}\hat{a}_{i})n_{i})]\mu_{0}$
$\displaystyle=2|\alpha_{\mathrm{g}}|^{2}\sum_{i}^{m}N_{i}\mu_{0}.$ (40)

Thus the relative purity is

$\displaystyle\frac{\mu_{m}^{-}}{\mu_{0}}=$
$\displaystyle[\frac{1}{2}(\sum_{i=1}^{m}{N}_{i})^{2}+\frac{1}{2}(\sum_{i=1}^{m}\frac{\tilde{{N}}_{i}}{n_{i}})^{2}+2|\alpha_{\mathrm{g}}|^{2}\sum_{i=1}^{m}N_{i}+|\sum_{i=1}^{m}k_{i}l_{i}\frac{n_{i}^{2}-1}{2n_{i}}|^{2}+|\alpha_{\mathrm{g}}|^{4}+$
$\displaystyle{\alpha_{\mathrm{g}}}^{2}(\sum_{i=1}^{m}\frac{k_{i}^{\ast}l_{i}^{\ast}(n_{i}^{2}-1)}{2n_{i}})+{\alpha^{\ast}_{\mathrm{g}}}^{2}(\sum_{i=1}^{m}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}})]/{(\sum_{i=1}^{m}{N}_{i}+|\alpha_{\mathrm{g}}|^{2})^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}+[\frac{1}{2}(\sum_{i=1}^{m}\frac{\tilde{{N}}_{i}}{n_{i}})^{2}+\frac{1}{2}|\alpha_{\mathrm{g}}|^{4}+|\sum_{i=1}^{m}k_{i}l_{i}\frac{n_{i}^{2}-1}{2n_{i}}|^{2}$
(41)
$\displaystyle|\alpha_{\mathrm{g}}|^{2}\sum_{i=1}^{m}N_{i}+{\alpha_{\mathrm{g}}}^{2}(\sum_{i=1}^{m}\frac{k_{i}^{\ast}l_{i}^{\ast}(n_{i}^{2}-1)}{2n_{i}})+{\alpha^{\ast}_{\mathrm{g}}}^{2}(\sum_{i=1}^{m}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}})]/{(\sum_{i=1}^{m}{N}_{i}+|\alpha_{\mathrm{g}}|^{2})^{2}}.$

Since $(\alpha_{\mathrm{g}}+\alpha_{\mathrm{g}}^{\ast})^{2}\geqslant 0$ and
$(\alpha_{\mathrm{g}}-\alpha_{\mathrm{g}}^{\ast})^{2}\leqslant 0$, thus we
have
$-2|\alpha_{\mathrm{g}}|\leqslant\alpha_{\mathrm{g}}^{2}+{\alpha_{\mathrm{g}}^{\ast}}^{2}\leqslant
2|\alpha_{\mathrm{g}}|^{2}$. Then

$\displaystyle
B^{2}(\sum_{i=1}^{m}\frac{k_{i}^{\ast}l_{i}^{\ast}(n_{i}^{2}-1)}{2n_{i}})+{B^{\ast}}^{2}(\sum_{i=1}^{m}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}})\geqslant-2|B|^{2}\sqrt{\sum_{i=1}^{m}\frac{k_{i}^{\ast}l_{i}^{\ast}(n_{i}^{2}-1)}{2n_{i}})}\cdot\sqrt{\sum_{i=1}^{m}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}}},$
(42)

where $B=\alpha_{\mathrm{g}}^{\ast}$. Furthermore, The left side of the
equation has

$\displaystyle|\sum_{i=1}^{m}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}}|^{2}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\frac{|k_{i}|^{2}|l_{i}|^{2}(n_{i}^{2}-1)^{2}}{(2n_{i})^{2}}+\sum_{i=1}^{m}\sum_{j=1}^{m}(k_{i}l_{i}{k_{j}}^{\ast}{l_{j}}^{\ast}+{k_{i}}^{\ast}{l_{i}}^{\ast}k_{j}l_{j})\frac{(n_{i}^{2}-1)}{2n_{i}}\frac{(n_{j}^{2}-1)}{2n_{j}}$
(43) $\displaystyle\leqslant$
$\displaystyle\sum_{i=1}^{m}\frac{|k_{i}|^{2}|l_{i}|^{2}(n_{i}^{2}-1)^{2}}{(2n_{i})^{2}}+\sum_{i=1}^{m}\sum_{j=1}^{m}(2\sqrt{|k_{i}|^{2}|l_{i}|^{2}|k_{j}|^{2}|l_{j}|^{2}})\frac{(n_{i}^{2}-1)}{2n_{i}}\frac{(n_{j}^{2}-1)}{2n_{j}}$
$\displaystyle=$
$\displaystyle(\sum_{i=1}^{m}\sqrt{|k_{i}|^{2}|l_{i}|^{2}}\frac{(n_{i}^{2}-1)}{2n_{i}})^{2}.$

So that

$\displaystyle
B^{2}(\sum_{i=1}^{m}\frac{k_{i}^{\ast}l_{i}^{\ast}(n_{i}^{2}-1)}{2n_{i}})+{B^{\ast}}^{2}(\sum_{i=1}^{m}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}})\geqslant-|B|^{2}\sum_{i=1}^{m}\sqrt{|k_{i}|^{2}|l_{i}|^{2}}\frac{(n_{i}^{2}-1)}{n_{i}},$
(44)

and we add $|B|^{2}\sum_{i=1}^{m}N_{i}$ to both sides of the equation

$\displaystyle|B|^{2}\sum_{i=1}^{m}N_{i}+B^{2}(\sum_{i=1}^{m}\frac{k_{i}^{\ast}l_{i}^{\ast}(n_{i}^{2}-1)}{2n_{i}})+{B^{\ast}}^{2}(\sum_{i=1}^{m}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}})\geqslant|B|^{2}(\sum_{i=1}^{m}N_{i}-\sum_{i=1}^{m}\sqrt{|k_{i}|^{2}|l_{i}|^{2}}\frac{(n_{i}^{2}-1)}{n_{i}}).$
(45)

Due to $\frac{n_{i}-1}{2}\geq 0$, thus

$\displaystyle(\sum_{i=1}^{m}N_{i}-\sum_{i=1}^{m}\sqrt{|k_{i}|^{2}|l_{i}|^{2}}\frac{(n_{i}^{2}-1)}{n_{i}})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{m}[(\sqrt{|k_{i}|^{2}\frac{n_{i}+1}{2}}-\sqrt{|l_{i}|^{2}\frac{n_{i}-1}{2}})^{2}+\sqrt{|k_{i}|^{2}|l_{i}|^{2}(n_{i}^{2}-1)}-\sqrt{|k_{i}|^{2}|l_{i}|^{2}}\frac{(n_{i}^{2}-1)}{n_{i}})]$
(46) $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}[(\sqrt{|k_{i}|^{2}\frac{n_{i}+1}{2}}-\sqrt{|l_{i}|^{2}\frac{n_{i}-1}{2}})^{2}+\sqrt{{|k_{i}|^{2}|l_{i}|^{2}(n_{i}^{2}-1)}}(\frac{\sqrt{n_{i}^{2}}-\sqrt{n_{i}^{2}-1}}{n_{i}})]$
$\displaystyle\geqslant$ $\displaystyle 0,$

so that

$\displaystyle|B|^{2}\sum_{i=1}^{m}N_{i}+B^{2}(\sum_{i=1}^{m}\frac{k_{i}^{\ast}l_{i}^{\ast}(n_{i}^{2}-1)}{2n_{i}})+{B^{\ast}}^{2}(\sum_{i=1}^{m}\frac{k_{i}l_{i}(n_{i}^{2}-1)}{2n_{i}})\geqslant
0.$ (47)

Therefore, the relative purity is always $1/2+$ a positive number, so that

$\displaystyle\frac{\mu_{m}^{-}}{\mu_{0}}\geqslant 1/2.$ (48)

It is a proof that for any Gaussian state, single-photon subtraction (or
addition) can reduce the purity with at most a factor 1/2.

The case of photon addition yields the same expressions with $\alpha_{g}$ and
$\alpha_{g}^{\ast}$ interchanged. In practice, the case of photon addition
will also lead to different values for $k_{i}$ and $l_{i}$. However, this does
not change any of the proofs and as a consequence, we also find
$\mu_{m}^{+}/\mu_{0}\geqslant 1/2$ for photon addition.

##  III C. Entanglement increase in all possible partitions for graph states

In analogy to Fig. 3(c) of the main text, we present in Fig. 4 the
entanglement increase for all the possible bi-partitions in the graph state
example of Fig. 2(b) in the main text. Every dot represents the entanglement
increase of a different bi-partition and the bi-partitions are grouped in
terms of the number of modes $m_{\cal A}$ in the subsystem ${\cal A}$ where
the photon is subtracted.

Figure 4: Each scattered dot represents a different bi-partition with $m_{\cal
A}$ modes in subsystem $\cal A$. The cyan and the pink scattered dots
correspond to the second example in Fig. 2 of the main text.