Open dynamics of entanglement in mesoscopic bosonic systems

A key issue in quantum information is finding an adequate description of mesoscopic systems that is simpler than full quantum formalism yet retains crucial information about non-classical phenomena like entanglement. In particular, the study of fully bosonic systems undergoing open evolution is of great importance for the advancement of photonic quantum computing and communication. In this paper, we propose a mesoscopic description of such systems based on boson number correlations. This description allows for tracking Markovian open evolution of entanglement of both non-Gaussian and Gaussian states and their sub-Poissonian statistics. It can be viewed as a generalization of the reduced state of the field formalism (Alicki 2019 Entropy 21 705), which by itself does not contain information about entanglement. As our approach adopts the structure of the description of two particles in terms of first quantization, it allows for broad intuitive usage of known tools. Using the proposed formalism, we show the robustness of entanglement against low-temperature damping for four-mode bright squeezed vacuum state and beam-splitted single photon. We also present a generalization of the Mandel Q parameter. Building upon this, we show that the entanglement of the state obtained by beam splitting of a single occupied mode is inherited from sub-Poissonian statistics of the input state.


I. INTRODUCTION
One of the important problems that is still developing in physics is the description of systems on the mesoscopic scale.This problem arises because full quantum treatment of large systems is not feasible.Classical description allows for an accurate depiction of the macroscopic systems; however, one can be interested in systems on the scale on which a full quantum description is not feasible, but with physical phenomena based on quantum features.In such a case, the classical paradigm becomes insufficient, and some intermediate methods are needed for this mesoscopic scale.The need to develop such theoretical tools emerges even more as experiments on this scale become nowadays feasible [1,2].This need is also fundamental from the perspective of quantum information, since its protocols require complicated quantum systems, and in their principles, they are based on quantum phenomena, such as superposition and entanglement [3].
Different semi-classical methods were introduced for the simplification of the problems, for example, the widely used parametric approximation [4,5].However, such techniques still require the full quantum description of some parts of the system.A different approach is to track observables that describe collective behavior of the system, e.g., total spin.Such a mean-field approach captures well macroscopic behavior of the system; however, information about quantum features is suppressed quickly with the number of particles.One can consider another collective observable, namely fluctuations instead of average values, as it is the case for the fluctuation algebra approach to mesoscopic systems [6,7].This technique is adopted in particular to describe entanglement in mesoscopic systems [8][9][10], and also different phenomena [11,12].Other methods were also used to describe the entanglement in mesoscopic systems [13][14][15].
One of the specific classes of systems for which mesoscopic descriptions are of great use is fully bosonic systems.The reason for this is that information encoded in bosonic modes is one of the promising avenues for implementing quantum information protocols.This is because, for example, information encoded in photons can be quickly transferred over long distances [16].However, such systems are also affected by the decoherence effects coming from the environment, and their strict theoretical considerations require dealing with the whole infinite-dimensional Fock space, which can be unfeasible in multi-mode systems.The widely used method that simplifies the description is to track the evolution of the covariance matrix of position and momenta.This reduces the problem to finite-dimensional when a finite number of modes are considered [17][18][19].In fact, this approach provides a full description of Gaussian states, and for non-Gaussian states, it can serve as a mesoscopic description.However, non-Gaussian features can be lost in this description, which can result in loss of useful information about quantum features of the state such as entanglement [20].Another mesoscopic approach for bosonic systems that was recently proposed is the reduced state of the field (RSF) [21].This formalism extends the mean-field approach with additional structure by the reduction of the state One of the crucial quantum phenomena is entanglement, which is one of the core resources in many quantum information tasks.This is even more so since in the quantum field theory the entanglement of pure states is fundamentally associated with Bell non-classicality [25], which is necessary for device-independent quantum key distribution.For such tasks, particularly compelling is the encoding of information in light modes, as light can be quickly distributed over long distances [16].What is more, fully photonic quantum devices are one of the promising avenues also in different subfields of quantum information, e.g., quantum computing [26].However, as all quantum systems, light modes are also affected by the environment and undergo non-unitary open evolution.This results in decoherence of the system and thus decay of quantum phenomena such as superposition (including entanglement).
In general, the full description of a non-unitary evolution of the state in terms of infinite-dimensional Fock space can be cumbersome or even infeasible.This is especially true for complicated mesoscopic systems that involve multiple bosonic modes occupied by an undefined number of particles.A custom in the treatment of mesoscopic systems is to consider only some degrees of freedom instead of the full quantum state of the field ρF .This allows for reducing the complexity of the problem.One of such methods for bosonic systems is to track the evolution of the covariance matrix of position and momenta [19]: where ⃗ Ξ = (x 1 , p1 , ...., xN , pN ), and operators xj , pj are dimensionless position and momentum operators in j-th mode, respectively.This method is particularly important, as it allows tracking of the full Gaussian evolution of Gaussian states.The reason for this is that one can reconstruct the full Gaussian state ρF from the covariance matrix.Therefore, the covariance matrix also contains full information about the entanglement of Gaussian states, which in the case of distillable bipartite entanglement can be accessed simply by using the PPT criterion [3,27,28].However, this is not always the case for non-Gaussian states, e.g., cat-like states.For such states, one has to employ different methods of entanglement detection for which information contained in the covariance matrix is not sufficient [20].Furthermore, the covariance matrix approach, while well suited for continuous variable analysis, can be non-intuitive with its symplectic structures.This may particularly concern researchers for researchers working solely with discrete variables, for whom finite-dimensional Hilbert spaces of the first quantization equipped with the Dirac notation can be more familiar.Additionally, regaining information from the covariance matrix about discrete variables also requires extra effort, and it is not always possible for non-Gaussian states.
In fact, multiple fundamental experiments and protocols, e.g., teleportation [29], entanglement swapping [30], and quantum cryptography [31] are based on entanglement of qubits realized using light modes restricted to a single photon per party.Emblematic examples of the states used in such protocols are polarization and path-entangled singlet states.Denoting by â † j the creation operator in j-the mode and vacuum state by |Ω⟩, such states can be written as: Here, one chooses partition into two parties (subsystems), such that polarization or path modes a 1 , a 2 are associated with one party and modes a 3 , a 4 with the second one.Here we have also used the following notation for Fock states: with semicolon marking the bipartition.Then, the measurements that correspond to Pauli matrices typically considered for qubits are realized as Stokes operators Θi [32] restricted to the single-photon subspace: where i = 1, 2, 3 denotes one of the three mutually unbiased basis described by two orthogonal modes a i , a i⊥ , and Π1 stands for the projector into the subspace of single photon across these modes.Additionally, after projecting with Π1 , the total photon number operator Θ 0 = â † i âi + â † i ⊥ âi ⊥ corresponds to the identity on the qubit space denoted by σ 0 .Note that any local operation on the two-qubit system can be described as a combination of tensor products of Pauli matrices (including σ 0 ) acting on subsystems A and B. In the Fock space, this translates to , where upper indices denote the party on which the given operators act.Therefore, as Stokes operators are, in fact, differences of number operators in two orthogonal modes, all non-classical correlations in such systems are described as photon number correlations This suggests that boson number entanglement is at the core of the discrete variable bipartite entanglement of bosonic systems.What is more, such correlations also become experimentally accessible in regimes with higher photon numbers due to developments in the field of photon number resolving detectors [33,34].Let us note that the considered photonic implementation of qubits is, in fact, approximate.This is because nothing prevents the given modes from being occupied by a higher or lower number of photons.It is an issue even with event-ready preparation of the system in which a photon is confirmed to enter the input of the apparatus, as the environment can add multi-photon noise (with an undefined number of photons) or cause loss of the photon which encodes the information.Such a problem arises especially during propagation through long transmission lines.Thus, even in these highly simplified scenarios, in principle, one should consider the full Fock space, or at least consider some mesoscopic description of the system.All of this motivates us to build a mesoscopic scale description of bosonic systems with the following matrix as the main structure that contains relevant degrees of freedom of the state ρF for considerations of bipartite entanglement: In the following, we show that ρ4 , which we call extended reduced state, contains information about entanglement of both Gaussian states with high average numbers of particles, and more importantly, non-Gaussian states (for which the PPT criterion cannot reveal entanglement based on the covariance matrix).At the same time, we equip the proposed mesoscopic description with multiple intuitive and familiar tools from finite-dimensional Hilbert spaces.Furthermore, in many circumstances, the open evolution imposed on the system is approximately Markovian.Thus, further on we show the set of evolution equations for the matrix ρ4 for such scenarios.This type of evolution of the full state of the bosonic modes ρF is described by the GKLS equation [35,36].In particular, this equation can be written for N orthogonal bosonic modes a j using the independent particles approximation in such a way that it includes relevant processes of particle decay and production, random elastic scattering, and classical coherent pumping [21,24].We base our considerations upon this equation, which reads: with , where â † k (â k ) are creation (annihilation) operators associated with corresponding modes, complex parameters ξ k describe classical pumping field, Ûj are unitary transformations of annihilation operators ( Û † j âk Ûj = k ′ u j k ′ ,k âk ′ with u j being unitary matrix) corresponding to scattering processes with their probability distribution κ j , and ↓ ) describe creation (annihilation) rate.

III. SEPARABILITY IN TERMS OF EXTENDED REDUCED STATE
The key question for entanglement considerations in our mesoscopic description is how separability in the Fock space translates into the extended reduced state.The following considerations are inspired by the mappings of the entanglement indicators for Stokes-like operators [37,38].Let us choose two sets of indices I A , I B ∈ {1, N } in such a way that I A ∩ I B = ∅, which determines our bipartition of the system.We assume that the modes with indices from I A constitute one subsystem and the modes with indices from I B the second one.In addition, each set contains at least two elements, and we denote the number of these elements by d X = I X .Based on this choice, we project the operator ρ4 in the following way: where: and Let us assume that the state of the system ρF after tracing out modes not contained in I A , I B is some pure state ρF = |ψ⟩⟨ψ|.In such a case, we can see that ρ ijnm has the form: where |ψ ij ⟩ = âi âj |ψ⟩.From this it follows that ρ ijnm forms a Gramian matrix and thus is positive definite.Thus ρΠ ψ (where the subscript denotes the original state on the Fock space) is also positive definite, and after normalization it is a well-defined density matrix ρN ψ for the system of two qudits with respective dimensions d X .However, there are exceptions for which the matrix cannot be normalized, as it is traceless.The trace of this matrix is equal to ⟨ NI A NI B ⟩ where NI X is the total boson number operator in modes corresponding to I X .This can be zero only when the state is the vacuum ρF = |Ω⟩⟨Ω| or when in the superposition it consists only of states with bosons only in one subsystem.We address this problem later on.For now, assume that the matrix can be normalized.Because any mixed state can be written as a convex combination of pure states ρF = i p i |ϕ i ⟩⟨ϕ i |, it maps to a convex combination of ρΠ ϕi : After normalization, we get: Note that 1 ≤ p ′ i ≤ 0 and i p ′ i = 1.Thus, this again forms a proper density matrix on the two-qudit space.In the next step, we assume that ρF is separable in the partition into modes I A and I B .We start with the case where ρF is in some pure state |ψ⟩⟨ψ|.Due to separability, the state can be written as |ψ⟩ = F † I A F † I B |Ω⟩, where |Ω⟩ stands for vacuum state, and F † Ii is a polynomial of creation operators of modes corresponding to I X .From this follows that each ρ ijnm from ρΠ factorizes as follows: where To see this, consider a rearrangement of the creation and annihilation operators in ρ ijnm such that the operators corresponding to I B are shifted to the right.Let us then put the resolution of identity 1 = ⃗ n1,⃗ n2 |⃗ n 1 ; ⃗ n 2 ⟩⟨⃗ n 1 ; ⃗ n 2 | between the operators acting on the modes I A and I B .This results in: where |⃗ n 1 ; ⃗ n 2 ⟩ stands for the Fock basis state with occupancy of the modes I A and I B described by vectors ⃗ n 1 , ⃗ n 2 , and ⃗ 0 stands for 0 photons in given modes.This rearrangement is possible as operators acting on orthogonal modes I A and I B commute.One can notice that terms of the first sum appearing in (13) contain factor ⟨⃗ n 1 ; ⃗ 0|F I B â † m âj F † I B |Ω⟩ which is equal to 0. This is because all operators inside the expectation value act only on modes I B leaving part of the bra and ket corresponding to modes I A orthogonal.Analogously, all other sums that appear in (13) are equal to 0, and the relation is proven.Now, based on the same arguments as above, ρ X ψ,l,k := ⟨ψ X k |ψ X l ⟩ forms a proper density matrix after normalization on the one-qudit space.Thus, ρN ψ = ρ1 ψ ⊗ ρ2 ψ /N is a separable state on the two-qudit space.For the mixed separable state, we see from (11) that ρN ρ F is a convex combination of pure separable states and therefore a proper separable mixed state on the two-qudit space.Therefore, if ρN ρ F is entangled, then ρF is entangled.Let us now address the scenario in which the matrix for the pure state is not normalizable.This was not important for mapping of entanglement indicators [37,38].However, it is important if one wants to use ρN ρ F in a more general context as a density operator.All pure states that result in the not normalizable ρΠ ψ have the following form: where ζ 0 stands for amplitude of vacuum state.These states can be either separable or entangled.One can easily find that for all such states ρΠ ψ = 0.This tells us that for states of the form (14) this state reduction does not allow one to access the information about their entanglement without some additional considerations.What is more, all of them are equivalent to vacuum (which is separable) at the level of reduction to ρΠ ψ .Despite that, there is a way to gain access to the entanglement of such states in our description.We discuss this further in examples.Note that such states are effectively highly unlikely to be present in experiments due to the presence of noise, for example, thermal noise.Therefore, the important question is the impact of such states if they are part of the mixed state ρF .Without loss of generality, assume that only |ϕ 0 ⟩ is of the form (14): where p ′′ i fulfill relations for probability.Therefore, even in this case ρN ρ F is a well-defined density matrix, and clearly the previous discussion on separability holds.In fact, the presence of the states (14) in the density matrix results only in renormalization of probabilities and thus neither helps to detect entanglement nor prevents it.One can observe that all states of the form: with |ψ j ⟩ having the form ( 14) and q j , p i being independent probability distributions, are equivalent to the state ρ′ in terms of ρN ρ F .Note that separability of ρΠ ψ does not imply that ρF is separable, since we use only part of the correlations.In other words, the entanglement of the state ρΠ ψ on the reduced two-qudit Hilbert space is a sufficient condition for the entanglement of ρF .Importantly, one has to consider at least 4 modes to reveal such a kind of entanglement, as one has to have both reduced subsystems at least two-dimensional.

A. Impact of mode transformations
Note that the choice of the basis in the subsystem defined by I i has no impact on the detection of entanglement.This is because the local unitary transformation of the annihilation operators that constitutes a basis transformation translates into the local change of basis in the reduced matrix ρ Π .Let us consider a general form of such a transformation âi = k U ik ĉk , where ĉk are annihilation operators in some new basis, and U ik determines the transformation between bases.Let us consider a transformation in the first subsystem: where the superscript in ρ (a) ijnm denotes the basis used.This transformation is simply a basis transformation on the subsystem.Such transformations do not affect the separability of the state.
From this we see that by choice of partition of modes our reduction procedure maps boson number entanglement into the entanglement of two distinguishable particles with finite number of levels.The extended reduced state ρ4 contains information about all partitions I A , I B .Thus, the separability of ρN ρ F for one partition does not necessarily mean that ρ4 does not contain any information on the entanglement.In fact, ρ4 contains information about the boson number entanglement in any partition I A , I B after an arbitrary unitary basis transformation of all modes, which can be seen analogously to (17).Still, sensible partitions easily emerge in the experimental setups, where specific modes correspond to different beams clearly determining local operations.

IV. OPEN EVOLUTION AND REDUCED STATE OF THE FIELD APPROACH
In the previous section, we effectively mapped the relevant in the context of entanglement degrees of freedom of a state ρF from the infinite-dimensional Fock space into the mixed state ρN ρ F on two-particle-like finite-dimensional Hilbert space.It is important to note that while ρN ρ F possesses all mathematical features of the density matrix, it does not inherit a probabilistic interpretation of the entries on the diagonal of the matrix.Here, they describe the strength of linear correlations of boson numbers between given modes.This situation is similar to the case of the reduced state of the field approach (RSF) to the mesoscopic description of bosonic modes [21].In fact, our reduction can be seen as a higher-order extension of the RSF.
Let us briefly recall the RSF approach.The RSF formalism attempts to reduce the infinite-dimensional description of the second quantization into a single-particle-like finite-dimensional Hilbert space of the first quantization.This reduction is carried out as follows.Consider a density operator ρF on the Fock space that describes the state of N orthogonal modes {a i } N i=1 .This operator is reduced to two structures on the N -dimensional Hilbert space equipped with orthonormal basis {|k⟩} N k=1 , single-particle density matrix and averaged field: The operator ρ contains information on the occupation of modes and coherences, while |α⟩ about the local phases of the field.Note that ρ and |α⟩ after normalization have the properties of a density matrix and a pure state, respectively; however, they do not inherit the probabilistic interpretation normally associated with quantum states.One can see that the diagonal elements of ρ represent the average occupation of modes instead of probabilities.The reduction of the state is also accompanied by a reduction of observables.One can reduce a class of additive observables on the Fock space in the following way: This construction preserves the expectation values: What is important, one can find the time evolution equations for ρ, |α⟩ corresponding to the evolution equation ( 6): where ĥ stands for reduced Hamiltonian (using ( 19)), and: This formalism is also equipped with notion of entropy: where k B stands for Boltzmann constant, and ρα = ρ − |α⟩⟨α|.

A. Extended reduced state of the field description
The RSF formalism, while simple, has drawbacks limiting its applications.Although it is fully based on quantum considerations, it is a highly classical description [24].This is in the sense that RSF does not contain information about non-classical phenomena like distillable entanglement.Therefore, it was conjectured that RSF does not contain any information on entanglement, which was strictly shown for two-mode Gaussian states.What is more, the entropy of this formalism has properties similar to the semi-classical Wehrl entropy rather than the quantum von Neumann entropy.However, it turns out that the information contained in RSF is needed to track the evolution of our reduced two-particle-like state ρN ρ F .Moreover, this evolution can be put into a framework that extends RSF in a natural way.This extension can be seen as the simplest extension of RSF that adds quantum features to the formalism and is not equivalent to the covariance matrix formalism.
Let us show how our extended reduced state ρ4 can be incorporated into an RSF-like formalism.States from the single-particle space of RSF contain information about particle numbers in modes.Thus, a natural concept is to introduce a tensor product of two single-particle spaces to track correlations of particle numbers.The extended reduced state ρ4 is then an operator on the two-particle Hilbert space constructed based on two copies of the singleparticle space.Note that reduction of state ρF does not always have to result in the operator ρ4 being Hermitian.Thus, it does not have to constitute a proper state on this Hilbert space, but still ρN ρ F does.We keep the whole structure ρ4 as it provides additional useful information on quantum features of the state.We will exploit this information later.Moreover, ρ4 is necessary for deriving a closed set of evolution equations that corresponds to (6).We build an extended reduced state ρ4 (5) to resemble a single-particle density matrix ρ (18), i.e., with indices of annihilation operators responsible for labeling rows and indices of creation operators labeling columns.Analogously to (19), one can reduce the observables: with tr ρF Ô4 = tr{ρ 4 ô4 }.We also complement the reduction of the state with two additional structures, a rank 3 tensor: and an additional operator on single-particle space: These two structures are introduced solely in order to include coherent pumping in the evolution, which will become apparent from the evolution equations.The evolution of operators on two-particle space also requires tracking the evolution in a single-particle space.Thus, this two-particle reduction has to be used together with the single-particle reduction, and the total reduced state consists of ρ4 , β, r, ρ, |α⟩.The Markovian evolution equations for operators ρ4 , β, r can be obtained analogously to equations (21,22).This is achieved by multiplying both sides of ( 6) with a suitable combination of creation and annihilation operators for a given matrix element and then tracing the resulting expressions.Then, after some simplifications done with the commutation relation [â i , â † j ] = δ ij , the set of evolution equations can be put in the following form: Here, we introduced the following notation: T i stands for transposition on i-th space, τ L acts on kets as follows |n, m⟩ → |m, n⟩ and τ R acts analogously on bras, and "T." stands for the transposed expression.Also, for some operator on single-particle space ô, we denote |k 2 , k 3 ⟩⟨k One can observe that if there is no coherent pumping (i.e., ∀ i ξ i = 0), evolution of the operators ρ4 and ρ decouples from evolution of the rest of the operators.Therefore, in such circumstances, ρ4 and ρ can be considered by themselves without considering β and r, for which there is a lack of a quantum state-like interpretation.Note that the first term of the evolution equation ( 28) for ρ4 resembles the Heisenberg equation, and it is responsible for the unitary evolution generated by the Hamiltonian.The remaining terms come from decoherence and pumping terms, as in the case of ρ.Let us also comment that r together with ρ and |α⟩ allow for reconstruction of the covariance matrix.Therefore, the extension of RSF by r is the simplest extension of RSF that gives quantum features to this formalism.However, it is trivial in the sense that it reproduces the covariance matrix approach, and it is also not a natural extension of this formalism.At the same time, the extension of RSF by ρ 4 follows the (multi)particle-like state interpretation of the original RSF.Thus, when coherent pumping is omitted, it provides the simplest extension of RSF that is natural to this formalism and adds the quantum features to the description.Still, one could make analogously further extensions of RSF using three copies of single-particle space, and so on, to track higher-order correlations in photon numbers.
The important feature of the set of equations (21,22) and (28)(29)(30) is that in fact they represent a finite set of first order in time differential equations.Therefore, if the time dependence of ĥ, |ξ⟩, γ↕ is given by simple functions, one can easily find analytic solutions.Additionally, one can solve this set order by order, as solutions of tensors of lower order are independent of solutions for tensors of higher order.Simply, one solves from the lowest order, i.e., |α⟩ and use the solution as non-homogeneous part for the rest of equations.
The introduction of implicit linear coupling of different modes in the system Hamiltonian Ĥ = k,k ′ ω k,k ′ a † k âk ′ is equivalent to taking a local approach in deriving the master equation where this coupling is considered to be a perturbation.Although it allows one to work explicitly with modes in desired basis, it is not always valid, and then a global approach should be considered (see [39] for a brief review).In the global approach, one first diagonalizes the system Hamiltonian and then considers the evolution of transformed modes.In the end, one transforms the obtained results by reverse unitary transformation of modes to the original basis.As described further in the text, this can be easily done on the reduced state.Therefore, one can easily apply also global approach.In fact, the unitary transformation that diagonalizes the reduced Hamiltonian ĥ is the same as the unitary transformation of modes that diagonalizes the original Hamiltonian.This allows for consideration of the system solely from the reduced perspective.

B. Local evolution
If one is interested in entanglement in a given bipartition, the structure of interest is the reduced two-particle-like state ρN ρ F .However, based on (28), the evolution of this matrix is, in general, not self-contained, and so the full ρ4 is required.Still, if the evolution is local in terms of a given bipartition, then the evolution of ρN ρ F is closed.In this case, "local" intuitively means that the evolution does not transfer particles between modes corresponding to two parties determined by the given bipartition or does not damp or pump them in a correlated manner.This translates to ĥ, γ↕ and all ûj being block diagonal in this bipartition, i.e., they can be written as a direct sum of matrices that act only on degrees of freedom of the given party.This allows the projector Π I A ,I B to commute with tensor products of such block diagonal matrices as ĥ ⊗ 1 in equation (28).From this, after projecting (28) into the subspace corresponding to Π I A ,I B , we get: where one simply replaces ρ4 → ρΠ ρ in equation ( 28) and projects the non-homogeneous part using Π I A ,I B .This allows restricting the considerations of entanglement in the extended RSF to ρN ρ F instead of tracking the whole ρ4 .

V. EXAMPLES: PPT CRITERION FOR MULTI-MODE BOSONIC FIELDS
We have shown that the sufficient condition for the entanglement of the full state ρF is that the reduced twoparticle-like state ρN ρ F is entangled.However, this does not mean that there exists a state ρF such that ρN ρ F is entangled.Therefore, in the following, we show the existence of such states by presenting two examples.What is more, we show that ρN ρ F contains information about the entanglement for Gaussian states and for non-Gaussian states for which the PPT criterion for the covariance matrix cannot reveal entanglement.
Because one can apply any entanglement detection method for bipartitie finite-dimensional systems to the ρN ρ F , let us consider one of the most widely known entanglement criteria, i.e., the PPT criterion [3].This criterion states that it is sufficient if partial transposition of the state results in: for state ρN ρ F to be entangled.In the simplest case of a two-qubit density matrix, the PPT criterion is a sufficient and necessary condition for entanglement.Therefore, when ρN ρ F is a two-qubit density matrix and the PPT criterion does not reveal the entanglement of ρN ρ F , then our reduced two-particle-like state is not able to reveal any entanglement of ρF .This does not imply that the state ρF is necessarily a separable state.This is because the entanglement of ρN ρ F is only a sufficient condition for the entanglement of ρF .In the case of a larger number of modes, if the PPT criterion does not reveal the entanglement of ρN ρ F , the reduced state might still keep information about the entanglement of ρF .However, in such a case, it has to be accessed using a different criterion.While one needs at least four modes for considerations of entanglement using ρN ρ F , in the following, we show a method of how one can also apply this methodology to reveal the entanglement of two-mode states.

A. Bright squeezed vacuum
As our first example, we will use the 2 × 2 bright squeezed vacuum (BSV).The non-classicality of this state was considered in stationary scenarios, including the high photon number limit [40][41][42], and in the context of correlation in Bose-Einstein condensates [43].Let us consider four orthogonal modes and the bipartition given by I A = {1, 2}, I B = {3, 4}.The BSV state is obtained by unitary evolution of the vacuum generated by the interaction Hamiltonian [4]: and therefore it is a Gaussian state.Now, the BSV state can be written in the form: where the parameter Γ = γt denotes the amplification gain, which is related to the power of the coherent source impinged on the nonlinear crystal in the parametric down conversion process, and Calculation of ρΠ ψ yields: where we have restricted this matrix only to the subspace determined by Π I A ,I B , and After normalization and applying the partial transpose, we find that (ρ N ψ ) T2 has two eigenvalues: The eigenvalue λ 1 is negative for all Γ > 0. Therefore, we applied the PPT criterion to show that the BSV state is entangled for any Γ > 0 and thus for an arbitrarily high average number of bosons.As a consequence, the two-particle space can keep information about entanglement in mesoscopic scenarios.

B. Single-photon
A s the second example, let us also consider the state of a single photon symmetrically beam splitted into modes a 1 , a 3 : This state is of particular interest as it is non-Gaussian and one cannot reveal its entanglement through the PPT criterion for the covariance matrix (see A). Clearly, this is a two-mode state, and as we argued to apply the above entanglement criterion, one needs four modes.However, one can trivially extend the state to four modes by including modes a 2 , a 4 with zero photons occupying these modes: Still, this state has the form (14), for which the RSF results in ρΠ Ψ = 0, and thus does not give access to information about entanglement.One can see that such a trivial extension of the two mode state will always result in ρΠ ρ F not revealing entanglement.This is because for such an extension there can be at most one nonzero matrix element ⟨â † 1 â1 â † 3 â3 ⟩ that is necessarily positive, and so ρΠ ρ F in such a case is always a PPT state.Note that this can be seen as a manifestation of the fact that one is not able to observe entanglement of two-mode entangled state with simple photon counting without ancillary states acting as a reference frame.
To avoid reaching the form ( 14), instead of (40), one can perform a different extension of the state to four modes by, for example, introducing coherent states with real amplitude α in the additional modes: This extension resembles performing a homodyne measurement on the single-photon state.However, for example, the assumption of matching energy was not done at this step, which would be required to perform the necessary interferometry for homodyne measurement.This is also not fully equivalent to performing homodyne measurement as one is interested in measuring the number of photons instead of quadratures.Such schemes are sometimes referred to as weak homodyne measurement [44][45][46].Now, using the same partition I A = {1, 2}, I B = {3, 4}, we get the following normalized reduced state ρN Ψ : After partial transposition, there are three eigenvalues of the resulting matrix: , where the first eigenvalue is negative for any finite α > 0. Therefore, even for two-mode entangled states, one can retrieve information about entanglement from the extended RSF.Note that the lowest value of λ 1 is obtained for α = 0.This is because one needs α ̸ = 0 to allow normalization of the matrix, yet coherent states do not add entanglement to the system but rather some form of noise.Thus, optimally, one should have α ≪ 1 in this scenario.

C. Time evolution of entanglement
Let us apply the rest of the extended RSF toolbox to consider how entanglement behaves under non-unitary evolution resulting from decoherence coming from the environment.We first consider that one distributes a singlephoton state |Ψ⟩ between two distant parties that want to perform a weak homodyne measurement locally to reveal the entanglement of the state.This is an interesting example because a device-independent quantum key distribution protocol can be constructed based on such a setup [47], and maintaining entanglement during transmission is necessary for its success.However, the state during transmission of the photon necessarily undergoes decoherence.Let us assume that during the distribution of the state, the transmission line is affected by the thermal environment with temperature T , while the local oscillators (coherent states) are approximately not affected.In this scenario, the operators γ↕ are diagonal and fulfil the relation: γ↑ = e −ℏω/k b T γ↓ , where ω stands for the angular frequency of the modes, and k b is Boltzmann constant.The creation rates are given by: where δ nm stands for Kronecker delta, and coefficient γ ω determines the coupling of the bath to modes with frequency ω, and N (ω) = 1/(e ℏω/k b T − 1) is Planck distribution of the average number of photons in the bath [48].Furthermore, as all modes are separated up to the moment of the measurement, there is also no coupling between modes.Consequently, the Hamiltonian is simply a free Hamiltonian for four modes with the same frequency.Based on our discussion in section IV B, for this scenario, we can consider only the evolution of reduced single-and two-particle states: ρ, ρΠ ρΨ .Due to the trivial free evolution Hamiltonian, the master equation simplifies to: Using initial condition for the reduced state ρ: we obtain the following time evolution of the ρΠ Ψ : where λ = N (ω) (e γωt − 1).Now, in order to access information about entanglement, we apply the PPT criterion after normalization.In the result, we find that there is again one eigenvalue that can be negative: where N stands for the normalization factor of ρΠ Ψ (t).Let us first consider λ − (t) in the limit t → ∞.In this case, there are two options: Remarkably, if α ̸ = 0 and the bath temperature is T = 0 (N (ω) = 0) in this limit, one has λ − (t) = 0, which suggests that the extended RSF formalism finds the state entangled for any finite time t.To see that this is the case, first note that: Clearly, the second term of the nominator has a higher absolute value, as it is a decreasing function of time with the limit t → ∞ equal to α 4 .Therefore, the second term determines the sign of the nominator to be negative.Because the denominator is always positive (as it is the sum of expectation values containing only photon number operators that are non-negative), we find that λ − (t) is also negative in this case for any finite t.This certifies the entanglement of the state in the whole time range.Still, for any T > 0, the limit of λ − (t) as t → ∞ is positive, and therefore at some finite time the entanglement will no longer be seen within the RSF formalism.Let us stress that this does not have to imply that the state is no longer entangled, as RSF provides only a sufficient criterion.This indicates that there is no detectable entanglement of the particular kind considered by RSF, i.e., photon number entanglement based on linear correlations.Solving λ − (t) = 0 for t, one finds that this criterion does not detect entanglement after the critical time t c : which clearly goes to infinity for N (ω) → 0 and to 0 for N (ω) → ∞.Interestingly, t c is independent of α and therefore α only affects the magnitude of λ − (t) and not its qualitative behavior in terms of its sign.Let us comment that, in general, extensions of the state do not have to be equivalent.For example, one can try to engineer an optimal ancilla to maximally violate some Bell inequality [49].Here, we also see this kind of inequivalence, since the different extensions (by different coherent states) give different magnitudes of λ − (t).In addition, one can have a completely different behavior, as found, for example, with the trivial extension considered before.Therefore, it is interesting to see that different non-trivial extensions result in the same t c .However, it is not a general property of all non-trivial extensions.Assuming only that the extended part of the state is separable from the original state, one can find for the considered single-photon state that: From this, one can see that any extension for which ⟨â † 2 â4 ⟩ = 0 will not find entanglement.Thus, in this scenario, one excludes any Fock states as meaningful extensions.Performing calculations analogous to (53), one finds: For all values of α axis is reached at the same time t c .The value of λ − (0) increases with α.However, the rate of growth of α greatly decreases for higher values of α.This can be seen as the curve for α = 0.1 quickly crosses other curves while still being negative, and therefore extremely small values of α are no longer optimal as in the unitary case.
which clearly depends on the form of the extension.Then, for the single-photon state (39), the optimal extensions are the ones that maximize If the ancillary state is a separable state, then by the Cauchy-Schwarz inequality, one can find that this quantity is bounded by 1 [50].Thus, since the extension by coherent states saturates this bound, it gives an optimal value of t c .Fig. 1 presents examples of the evolution of the eigenvalue λ − (t) in different regimes.Clearly, entanglement becomes weaker with time evolution, and its decay becomes faster with increasing N (ω).In general, this shows that the single-photon entangled state is a resource that is quite robust against thermal damping when the thermal bath has a low temperature.
Let us note that one can propose one extension which is not built from a separable state.In particular, a second copy of the original state can be taken as such an extension.This still indicates that the observed entanglement is due to the entanglement of the original state.However, in such a case, it is rather sensible that the ancillary modes undergo the same evolution as the modes of the original state.This is because the state of the ancillary modes cannot be prepared locally.Thus, for the considered scenario, one should take the following creation rates: For this case, t c is exactly the same as for the case of the extension by coherent states.Thus, this simple extension by the second copy of the state is also optimal in this case.Now, let us consider an analogous scenario for the BSV state where all modes are coupled to the thermal environment with the same temperature during the distribution of the state between parties.Therefore, for this case, we consider the creation rates (56).Then, using as initial conditions reduced states (36) and we get the following solution for ρΠ ψ : where Applying the PPT criterion, we find that the negative eigenvalue evolves as follows: With the realization that the only time dependence is in λ, which for the special case of T = 0 is equal to 0, one can see that the entanglement of the BSV state is, in fact, completely stationary under the considered non-unitary evolution, i.e., This shows the exceptional robustness of the BSV state to low-temperature damping.For any non-zero temperature (T ̸ = 0), there is a point in time at which the entanglement is not detected any more: and for t → ∞ the eigenvalue λ 1 (t) converges to 1/4.Note that t c is an increasing function of the pumping parameter with the limit Compering this expression to critical time for the single-photon state (53), one can see that the critical time for the highly pumped BSV state is more robust against this type of decoherence.Still, for small values of Γ where the BSV state is dominated by the vacuum component, t c can be smaller than for the single-photon state.This happens for Γ < − log 2 − √ 2 /2 ≈ 0.27.Fig. 2 shows comparison of t c for BSV state with different values of Γ and single-photon state.

D. Relation with covariance matrix for Gaussian states
Let us stress that detection of entanglement in a case of some Gaussian states does not imply that ρN ρ F can reveal entanglement for any entangled Gaussian state.The reason for this is that while the covariance matrix keeps full information about the Gaussian state, our ρN ρ F stores information only about a particular type of correlations.
The general inequivalence of these criteria for Gaussian states can be seen based on the two-mode squeezed vacuum state.As we discussed earlier, for the trivial extension of the two-mode state, one is unable to find the entanglement in terms of ρN ρ F , and the proper extension has to be chosen.At the same time, the PPT criterion for the covariance matrix finds that the state is entangled.However, choosing a suitable extension will allow the entanglement of the considered state to be found using ρN ρ F .As in the case of the single-photon state, one can take a second copy of the state as an extension to find the entanglement.Then, the extended state is, in fact, analogous to the BSV state (with the difference that the two ancillary modes are phase-shifted in total by π), for which ρN ρ F keeps information about entanglement.To see this, note that the Hamiltonian that generates the BSV state (33) consists of two commuting parts Ĥ14 and Ĥ23 acting on different sets of modes: {a 1 , a 4 } and {a 2 , a 3 }.Thus, the unitary evolution of the modes under this Hamiltonian factorizes as e −i Ĥ14t/ℏ e −i Ĥ23t/ℏ where each exponent corresponds to two-mode squeezing.From this it is apparent that this evolution with the vacuum as the initial state creates a product of two two-mode squeezed vacuum states.However, squeezing parameters for these pairs of modes have the opposite sign based on (33).
Observe also that entanglement in Gaussian states is created only through squeezing.In particular, two-mode squeezing creates correlated particle pairs.Therefore, since this type of correlation is in the core of ρN ρ F , it should be efficient in detecting the entanglement of Gaussian states.Moreover, entanglement from single-mode squeezing appears by beam splitting or, more generally speaking, by linear couplings with other modes.From our discussion, we saw that the entanglement created by the beam splitting can be captured by ρN ρ F , and thus is similar in nature.Consider, as an example, a specific class of entangled Gaussian states.In particular, assume that one starts with the state resulting from some squeezing of different pairs of modes.Then, allow the full state of these modes to evolve based on some Hamiltonian that admits the RSF reduction and couples only modes local to a given party.Before time evolution, one should be able to find if such a state is entangled using ρN ρ F based on our disscusion (this may require taking a copy of the state as the extension).Now, the evolution equation ( 31) in our scenario is simply a Heisenberg equation.It generates a unitary evolution that factorizes with respect to subsystems: e −i ĥA t/ℏ ⊗ e −i ĥB t/ℏ , where ĥA(B) is a reduced Hamiltonian restricted to the modes of subsystem A (B). Now, one can write ρN ρ F in the product basis that diagonalizes ĥA and ĥB .Importantly, the choice of the particular local basis does not affect whether the PPT criterion finds entanglement.Then, one is left with local free evolution that does not destroy entanglement in the given bipartition.Thus, ρN ρ F detects entanglement for the considered class of Gaussian states.

E. Maximally mixed states
Let us consider the limit t → ∞ of the reduced two-particle-like state ρΠ ψ (t) in (58).In this limit, one gets: After normalization, this matrix corresponds to the maximally mixed state on the subspace of the two-particle space restricted to the given bipartition.As the bath in this example was thermal and the same for all modes, this reduced two-particle-like state corresponds to the four-mode thermal state for a given temperature T .Note that the temperature does not affect the normalized reduced state.Analogously, we have ρ = diag(N (ω), N (ω), N (ω), N (ω)), which constitutes a maximally mixed state in a single-particle space.As the inclusion of additional separable modes in the thermal state for a given T cannot add any new type of correlations seen in the reduced states, one can see that the thermal states always correspond to maximally mixed states on single and two particle spaces.

VI. PASSIVE OPTICAL ELEMENTS
Most optical devices include passive optical elements, such as beam splitters and phase shifters.The RSF formalism allows for easy introduction of such elements into the system.Therefore, this formalism allows for a straightforward analysis of experimental setups.Let us start with the phase shifter.We assume that the beam contains the mode a i at some point in time, starts the interaction with the phase shifter at t = τ 0 and ends it at t = τ e , which adds the phase to this particular mode.This can be achieved by the following modification of the Hamiltonian: where I [τ0,τe] (t) stands for the indicator function on the time interval [τ 0 , τ e ], and ϕ i determines the phase gain rate of i-th mode.The resulting final phase shift is given by: ∆ϕ = (τ e − τ 0 )ϕ i .
The beam splitting of two modes a i , a j into two modes a ′ i , a ′ j can be effectively realized as a unitary transformation of the corresponding annihilation and creation operators given by: where R and T stand for reflectivity and transmissivity, respectively.Before beam splitting, modes a ′ i , a ′ j are in the ground state, and after this transformation, modes a i , a j are left empty.Therefore, their roles interchange, and tracking modes a i , a j is of no use after the operation.Therefore, just after the operation, one can rename the modes a ′ i , a ′ j to a i , a j , keeping the structure of the reduced state.Based on equation (17), one can find that any unitary transformation U of annihilation operators reduces to the operator: where U ji stand for matrix elements of U.Then, under beam splitting operation, the full reduced state transforms as follows: Note that one can analogously incorporate to the analysis phase shifting in an instantaneous manner as it results in unitary mode transformation.This reduction of unitary transformations of modes also allows one to consider the inefficiency of the detectors, which can be modeled as a two-port beamsplitter on the path of the mode a i with transmissivity equal to efficiency T = η and second output mode traced out.Denoting the efficiency of the detector at mode a i as η i , inefficient detectors can be applied by operation: Then, the reduced state ρ transforms as ρ → ηρη, and other parts of the reduced state transform analogously.
Assuming equal efficiency η of all detectors, one can see that entanglement of the BSV state can be revealed for arbitrary efficiency, as the matrix ρΠ ψ (36) is varied only by constant coefficient η 2 , which is removed after normalization.Thus, after partial transposition, the eigenvalues are not changed.In fact, this applies to any state and any entanglement criterion by the same means.
Let us go back to the unitary transformation U d of the modes that diagonalize the original Hamiltonian Ĥ.Clearly, the diagonalized Hamiltonian Ĥd reduces to the diagonal Hamiltonian ĥd .As the unitary transformation of modes reduces by (68), one gets: where the equality stands for equality of the matrices representing operators in specific basis in which h d is diagonal.Therefore, the reduced operators ûd perform the diagonalization of the reduced Hamiltonian.On the other hand, if one finds a unitary operator ûd that fulfills (71) for some ĥd , because the reduction of the Hamiltonian is a bijection, one knows that there is a corresponding diagonal Hamiltonian in the Fock space Ĥd that must be obtained by mode transformation from the original Hamiltonian Ĥ.However, since the reduction of unitary transformations of modes is also bijective, this transformation is obtained from ûd .
A. Homodyne measurement and squeezed states pumping One of the tools used in quantum optics is homodyne measurements.In such measurements on state ρF , one interferes the coherent state through beamsplitter with state ρF and performs intensity measurement or photon counting.Such schemes can be analyzed in the RSF formalism as it requires only extending the reduced state by modes in coherent states (see the B) and performing a beam splitting operation on the obtained reduced state.However, as presented above, in some circumstances one does not even have to mathematically perform a beam splitting operation to get access to information about entanglement.
One of the drawbacks of the RSF formalism is the absence of the possibility of direct inclusion of the squeezing in the evolution of the system.However, there is a partial solution to this problem.In many experimental setups, first some modes are impulsed pumped with squeezed states coming, for example, from parametric down conversion and then combined with the rest of the system through beamsplitters.As one can easily calculate the reduced state for the squeezed states, analogously to the case of homodyne measurements, one can combine the modes in squeezed states with the reduced state of the system.Therefore, one can effectively include some form of squeezing in the system.What is more, one can also include in such a way higher-order squeezed states [5].

VII. MANDEL-Q PARAMETER AND TWO MODE ENTANGLEMENT
Our extension of RSF is not limited only to revealing entanglement in terms of non-classical phenomena.One of the non-classical phenomena in quantum optics is the sub-Poissonian photon statistics of light, which has no classical counterpart.To quantify this phenomenon, one can use the Mandel Q parameter [51]: If Q i < 0, the state of the mode a i has sub-Poissonian statistics.This non-classical behavior is still widely studied and used, for example, for source verification [52][53][54], in particular, using the time-dependent version of this parameter.
Clearly, the extended RSF formalism keeps information on photon statistics and its time dependence in the operators ρ(t), ρ4 (t) and so the Mandel Q parameter can be found as: This is one of the reasons why one would want to include ρ4 in the extension of RSF instead of only ρΠ .The Mandel Q parameter, which is a shifted and rescaled variance, allows for detecting non-classical behavior internal to the state of a single mode.Still, other parts of the covariance matrix (note that this is a different covariance matrix from the previously mentioned covariance matrix associated with quadratures) that mixes two modes a i , a j can also be used to reveal non-classicality, which is shared among multiple modes.Let us build a generalization of the Mandel Q parameter using the criterion for the two-mode entanglement from [50]: where we have used commutation relation [â j , â † j ] = 1 and the fact that ⟨â † i âi ⟩ is non negative.In fact, this parameter contains the covariance of the operators â † j âi , â † i âj , and clearly it highly resembles the Mandel Q parameter as Q ii = Q i .In terms of RSF, one can calculate the parameter Q ij as: For the maximally entangled state of a single photon shared among two modes (39), one gets Q 13 = Q 31 = 1/2.For this state, the time evolution analogous to the one considered in previous sections leads to: One can find that Q 13 = 0 at exactly the same time t c as the scenario with ancillary coherent states (53).For Q 13 , adding ancillary states to reveal the entanglement turned out to be unnecessary as this parameter uses additional information contained in ρ.
Let us also remark some additional similarities of Q ij with Q i .Consider Q 12 for a two-mode coherent state in modes a 1 , a 2 with complex amplitudes α 1 , α 2 .In this scenario, one can find that: This shows that, starting from the product input state (80), the two-mode entanglement after beam splitting is exclusively inherited from the sub-Poissonian statistics of the input state.Therefore, beam splitting redistributes the internal non-classical photon number correlation of the state to the photon number entanglement of the output modes.Note that due to this equality, Q out 12 + Q out 21 is lower bounded by -1 and minimized by the input states being Fock states ρφ = |n⟩⟨n|.

VIII. CONCLUDING REMARKS
In summary, we have proposed an extension of the reduced state of the field formalism for mesoscopic bosonic fields.Our extension allows for tracking the evolution of particle number entanglement for systems undergoing non-unitary Markovian evolution, which is not possible in the original reduced state of the field approach.The extended approach allows for studying entanglement of both Gaussian and non-Gaussian states.This was achieved by reducing the problem of entanglement of multi-mode bosonic state into the entanglement of the two-particle state on a finite-dimensional Hilbert space.Based on the proposed approach, we showed that particle number entanglement of multi-mode bosonic field is independent of the efficiency of the detectors if it is the same across all detectors.Moreover, we showed that the entanglement of the states of the beam-splitted single photon and 2 × 2 bright squeezed vacuum are robust against interaction with a low-temperature thermal environment.Furthermore, we considered the generalization of the Mandel Q parameter constructed within the proposed approach, which describes non-classical correlations both in a single mode and shared among different modes.Using this parameter, we showed that the entanglement of two-mode entangled states obtained by the beam splitting of one occupied mode is solely inherited from the sub-Poissonian statistics of the input state.
The proposed reduced state of the field description can provide a versatile and intuitive tool for the analysis of the impact of decoherence on non-classical phenomena in multiple optical experimental setups and consequently their designing.This is because our approach inherits multiple tools from standard quantum mechanics of finite-dimensional systems, with only a slight change in interpretation.As analysis of a given setup requires only specifying mesoscopic properties of the sources instead of the whole initial state, it could be done without even invoking the Fock space.Furthermore, one can easily incorporate passive optical elements and inefficiencies in the description of the evolution.In addition, ancillary modes and entanglement can also be added to the system through the use of beamsplitters.
One of the important limitations of proposed formalism is the class of accessible evolution equations.This class is restricted to master equations with terms that are at most linear in the creation and annihilation operators.Often, it is possible to transform the considered problem into the problem described by such evolution through employing, e. g., Bogoliubov transformations of modes and through simplifying interactions using the mean field.However, if one wants to transform back to the original modes, one has to consider the evolution of additional structures.This is because the information contained in ρ 4 is insufficient.Alternatively to Bogoliubov transformations, one can directly calculate the evolution of ρ 4 for quadratic Hamiltonians by including the additional structures mentioned above (see C).However, for higher-order interactions, one needs to track higher-order correlations in order to propose a sensible approximate set of closed evolution equations.Still, whenever one cannot use the convenient evolution equations for the reduced state described in this work, one can use the extended reduced state in the analysis of non-classical features of the state.Another limitation of the formalism is that it is restricted to two-party entanglement and that entanglement of the reduced state is only a sufficient criterion for entanglement.Nevertheless, one can build analogously reduced three-and more-particle-like states and upon it consider more involved types of entanglement.
Let us also note that since the extended reduced state of the field preserves a state-like interpretation of reduced states and maintains a Hilbert space structure for them, one can exploit scalar product of these Hilbert spaces for geometric-like considerations of bosonic states as, e.g., in [55].This allows for considerations of proximity of states in terms of reduced states and therefore for their classification in terms of mesoscopic properties.
As we found that for both single-photon entangled state and bright squeezed vacuum the entanglement is robust against thermal damping in low temperatures, an interesting open question arises.Is this property universal to photon number entanglement?
The proposed extension of the reduced state of the field approach could also be used to examine different quantum optical phenomena.Let us recall that the original reduced state of the field contained information about the expectation values of Stokes operators, providing a generalization to the description of polarization in terms of Jones vectors.Our extension, when restricted to two modes, additionally keeps information about second-order polarization tensors [32,56].As these tensors are associated with non-classical phenomena of hidden polarization, our approach could give additional insight into this phenomenon.
where subscript t stands for the total state and indices A and B for the corresponding parties.Note that swap operations τ L , τ R can be performed using N 2 × N 2 -dimensionl matrix τ , where N is total number of modes: Then, for some operator Ô, one gets: Appendix C: Full second-order evolution In a manner similar to the presented extended RSF approach, one could consider the full second-order evolution of the system.In such a scenario, one could consider reduction of the second-order observables to three reduced operators on single particle space: Note that reducing also first-order observables: We can now allow any Hamiltonian that can be reduced by the procedure (C1).We denote the reduced Hamiltonian by the following triple ( ĥ, ĥs , ĥ * s ).The operator ĥ governs the free evolution and exchange of photons between modes, whereas the operators ĥs , ĥ * s correspond to single-mode and two-mode squeezing.Introducing such an interaction Hamiltonian can correspond, e.g., to using a parametric approximation to some modes in the system.These modes are treated as classical coherent pumping fields which are subject to, e.g., a parametric down-conversion process in which annihilation of a single photon from the pumping field is accompanied by creation of two photons in modes which we treat in a full quantum manner.This addition to the formalism is similar to the addition of coherent pumping in (6).However, such an extended evolution requires an extension of the reduced states.If one is only interested in the evolution of the original reduced state ρ the addition of r is sufficient.Then, one can analogously obtain equations of motion for RSF, which are the following: (1 ⊗ ĥ m + m1 ⊗ ĥT ) − 2i ℏ (ρ 4 1 ⊗ ĥ † s + {1 ⊗ ĥ † s , ρ4 } T2 − qĥ s ⊗ 1) ℏ ( ĥ † s ⊗ 1 mT ) τ L + ( mĥ † s ⊗ 1) τ R + ĥ † s ⊗ 1 m + mT ĥ † s ⊗ 1 − 2i ℏ r ⊗ ĥ † s + (r ⊗ ĥ † s ) τ L + (r ⊗ ĥ † s ) τ L T2 − i ℏ ĥ † s ⊗ r + ( ĥ † s ⊗ r) τ L + ( ĥ † s ⊗ r) Note that such an extended RSF allows for calculation of any additive observable up to fourth the order in creation and annihilation operators.

FIG. 1 :
FIG.1: a) Time evolution of λ − (t) for α = 0.5 and different values of N (ω).The solid red line presents the value of λ − (t) in the case of unitary evolution without a thermal environment.For N (ω) = 0 value of λ − (t) converges to zero but is always negative, showing the entanglement of the state.With increasing N (ω) the photon number entanglement starts to decay faster, eventually disappearing.b) Time evolution of λ − (t) for N (ω) = 0.1 and different values of α.For all values of α axis is reached at the same time t c .The value of λ − (0) increases with α.However, the rate of growth of α greatly decreases for higher values of α.This can be seen as the curve for α = 0.1 quickly crosses other curves while still being negative, and therefore extremely small values of α are no longer optimal as in the unitary case.

FIG. 2 :
FIG.2: Comparison of t c for the single-photon state and the BSV state for different Γ as a function of N (ω).The critical time t c for the BSV state quickly converges to the limiting case Γ → ∞ and reaches a significant advantage in robustness over the single-photon state.Still, for small values of Γ, entanglement of the single-photon state is more robust against non-unitary evolution resulting from thermal bath.