Gaussian boson sampling at finite temperature

Gaussian boson sampling (GBS) is a promising candidate for an experimental demonstration of quantum advantage using photons. However, sufficiently large noise might hinder a GBS implementation from entering the regime where quantum speedup is achievable. Here, we investigate how thermal noise affects the classical intractability of generic quantum optical sampling experiments, GBS being a particular instance of the latter. We do so by establishing sufficient conditions for an efficient simulation to be feasible, expressed in the form of inequalities between the relevant parameters that characterize the system and its imperfections. We demonstrate that the addition of thermal noise has the effect of tightening the constraints on the remaining noise parameters, required to show quantum advantage. Furthermore, we show that there exist a threshold temperature at which any quantum sampling experiment becomes classically simulable, and provide an intuitive physical interpretation by relating this occurrence with the disappearance of the quantum state's non-classical properties.


I. INTRODUCTION
Gaussian boson sampling (GBS) is a computational task that, under widely accepted complexity-theoretic conjectures, is believed to be intractable using classical machines [1][2][3].The original formulation of the problem consists of sampling from a non-classical m−mode Gaussian state obtained by sending identical squeezed vacuum states through a passive linear optical network (LON), with photon-number-resolving (PNR) detectors.More experimentally-feasible variants of the task employing threshold detectors [4] and click-counting detectors [5] have been proposed since.The advancements of photonic platforms in recent years made an experimental GBS demonstration of quantum advantage feasible with current technological capabilities, with state-of-the-art experiments consisting of 216 modes and up to 125 measured photons [6].Besides the quest for quantum advantage, Gaussian boson samplers can also be used to tackle problems of practical interest such as simulating molecular vibronic spectra [7], measuring graph similarity [8], perfect matching counting [9], identifying dense subgraphs [10,11], and predicting stable molecular docking configurations for drug design [12,13].
As near term quantum devices do not benefit from error correction, it is well known that sufficient noise can prevent GBS experiments from entering the regime where quantum advantage is, in principle, attainable.Extensive research has been conducted to investigate the impact of various sources of noise and imperfections on the classical simulability of the sampling task.These include losses [14], partial distinguishability [15,16], and detectors' inefficiencies and random counts [17,18].In particular, Ref. [17] provides sufficient conditions for efficient classical simulability of generic quantum optical experiments − GBS being a specific instance − expressed in the form of inequalities that involve the noise parameters.The method is based on expressing the output probability distribution in terms of phase-space quasi-probability distributions (PQDs) of the input state, the measurement operators and the transition function associated with the specific quantum process, and further identifies their negativity as a necessary con-dition to achieve quantum advantage [19,20].
It is well known that a sampling task with thermal state inputs can be efficiently simulated (i.e., in polynomial time) on a classical computer [21].This fact suggests that there should exist a transition in the computational complexity of GBS as temperature grows.Nevertheless, finite-temperature effects have received limited attention in this setting thus far.This paper investigates in this direction by assessing how the addition of thermal noise affects the classical intractability of simulating quantum optical sampling experiments.We model thermal noise by introducing additional environmental modes initialized in thermal states that interact with the system's modes via a (passive) linear bosonic channel [22].Significant attention is then dedicated to predicting how much thermal noise can a GBS experiment tolerate before it becomes classically efficiently simulable.
In addition to fundamental interest, thermal noise effects are particularly relevant for experiments conducted in the MHz domain, e.g.those involving the phononic modes of motion of trapped ions [23,24], or those in the GHz domain that employ superconducting architectures.These platforms offer highly efficient photon-number-resolving detection [25], enabled by quantum-non-demolition measurements [26] that allow for repeated detection, ultimately leading to higher measurement fidelities.Moreover, superconducting circuits provide an excellent degree of control over the required interactions − namely beam splitter and phase shifter operations − to build an arbitrary passive LON [27].Squeezing and displacement operations may also be achieved in circuit QED, thus allowing for the implementation of GBS experiments [28].The non-linearities provided by Josephson junctions enable efficient preparation of non-classical states of light [29], including multi-photon Fock states [30,31], and make it possible to engineer non-linear operations that would otherwise be challenging to implement in optical systems.A notable example is given by Kerr-type unitaries, which have recently been introduced into the boson sampling framework as a mean to enhance the system's robustness against noise [32].
The rest of this paper is structured as follows.In Sec.II we revise key aspects of the phase-space formulation of quantum optics and generalize the formalism introduced in Ref. [17] by including finite-temperature effects.This allows us to obtain a general sufficient condition for the classical simulability of a generic quantum optical experiment, that we then apply to a noisy GBS instance employing threshold detectors.We find that, as one might expect, the additional thermal noise has the effect of reducing the detection imperfections sufficient to efficiently simulate the sampling task.In Sec.III we show that, under the assumption of a uniform loss rate, there exists a threshold temperature − which depends on the system's losses − at which any sampling experiment becomes classically simulable, even in the presence of ideal detection.We provide a physical interpretation of this phenomenon and show that it is linked to the disappearance of the state's genuine quantum properties.In Sec.IV we build on the approach of Ref. [18] and present a sufficient condition for the classical simulability of noisy GBS experiments that takes into account both thermal noise and approximate sampling, overcoming the main limitation of Ref. [17].Lastly, in Sec.V we summarize our findings and provide concluding remarks.

II. SUFFICIENT CONDITIONS FOR EFFICIENT CLASSICAL SIMULATION OF QUANTUM OPTICS AT FINITE TEMPERATURE
In this section we investigate how thermal noise affects the noise thresholds in Ref. [17] that are sufficient for an efficient classical simulation of a generic bosonic experiment.Our main result resides in Eq. ( 20), the only assumption to derive such classical simulability condition being the model of noisy evolution we employ, outlined later in this section.We then consider the special case of a noisy GBS device and show that the latter may be efficiently simulated by classical means if the following inequality is satisfied Here, r is the squeezing parameter of the input states, n is the mean number of environmental thermal photons, η L denotes the transmission of the LON, and η D and p D are the threshold detectors' efficiency and their dark count rate, respectively.A generic quantum optical experiment is described in terms of an m-mode initial state ρ, an m-mode quantum process represented by a completely positive (CP) map E and a measurement performed on the final state E(ρ).A quantum measurement is characterized by a positive operator-valued measure (POVM), i.e., a collection of operators {Π n } satisfying the conditions Π n ≥ 0 and n Π n = I, where I denotes the identity operator on the Hilbert space.The probability of obtaining a specific measurement outcome n is given by the Born rule p(n) = Tr{E(ρ)Π n }.The outcome probability can alternatively be expressed in terms of ordered phase-space quasi-probability distributions as Here, W Πn and ρ denote the PQD of the POVM element Π n and that of the input state ρ, respectively.In what follows, the dagger transposes a vector of complex numbers to a column vector and takes a complex conjugate.The s−ordered PQD (s−PQD) of a generic m-mode Hermitian operator O is defined as Here, s = diag(s 1 , . . ., s m ) is the diagonal matrix of the ordering parameters s j ∈ R and D(ξ) is the m-mode displacement operator a = (a 1 , . . ., a m ) being the vector of bosonic operators.
The well known Husimi Q-function, Wigner function and Glauber-Sudarshan P -function are retrieved for s = −I m , s = 0 and s = I m respectively, where I m denotes the m−dimensional identity matrix.It is worth noting that the s−PQD of a quantum state is normalized to one, however it can in general attain negative values and diverge more severely than a delta function.The remaining function appearing in Eq. ( 2) is the transition function associated with the quantum process E, defined as One can also prove that meaning that the action of E on a multimode coherent state |γ⟩ is enough to completely characterize the transition function in Eq. ( 5).If there exist values of the ordering parameters t and s such that W ρ (α) are all non-negative and well-behaved, then it is possible to sample from p(n) efficiently.We emphasize that this condition is only sufficient, and there might exist other efficient simulation strategies that succeed even in regimes where the PQDs exhibit negativities.It should also be noted that this framework lets us address the feasibility of efficient exact simulations only, i.e., sampling from p(n), rather than sampling from an approximation of the latter.Despite this shortcoming, the strength of this formalism resides in the wide range of applicability enabled by its modular nature, which allows us to investigate the classical simulability of generic quantum optical experiments.
Let us now assume that E is a CP map describing a LON subject to both photon loss and thermal noise.We adopt a simple model where, alongside the system's m modes, we consider m additional environmental modes, each initialized in the thermal state Here n is the mean photon number for the given temperature and a is the annihilation operator of the corresponding environmental mode.These 2m modes interact by means of a fictitious ideal interferometer described by the block unitary matrix The unitarity of U implies that i.e., L is a subunitary matrix when losses are present in the system.Hence, the action of the noisy LON on an m-mode coherent state |γ⟩ reads where U is the unitary operator associated with the larger 2mmode interferometer and the trace is taken over the environmental degrees of freedom.In Fig. (1) we display a schematic representation of the noise model employed to describe the LON.We can expand the m-mode thermal state over the coherent state basis, making use of its P -function representation with where k = 2n + 1.We can thus write Eq. ( 10) as The action of the larger LON on a multi-mode coherent state is easily computed thus leading to (15) By substituting this expression into Eq.( 6), we can compute the trace appearing in Eq. ( 5) ) where we have used the identity as well as the unitarity constraint Eq. ( 9) and standard multidimensional Gaussian integration.We now plug Eq. ( 16) into Eq.( 5) and finally obtain the transition function The latter is non-negative, well-behaved and has multi-variate Gaussian form iff Furthermore, we can always find t, s ∈ R m such that the input state t−PQD is non-negative for t ≤ t and the (−s)−PQD associated with the quantum measurement is non-negative for s ≥ s.Hence, the following inequality constitutes our sufficient classicality condition for an efficient simulation of the sampling task described above to be feasible.We emphasize once more that this condition is only sufficient.On the other hand, the modular nature of the formalism enables its wide applicability, our sole assumption being the noise model of the linear optical evolution given by Eq. (10).
As expected, the results of Ref. [17] are retrieved in the zerotemperature limit k = 1.
We can now apply Eq. ( 20) to a noisy GBS experiment.In particular, let us consider an initial state comprising of m identical squeezed vacuum states m j=1 S(r) |0⟩, where is the single-mode squeezing operator and r > 0 is the squeezing parameter.One can show that the t−PQD of a generic m-mode Gaussian state ρ reads Here, σ is the covariance matrix and α is the displacement vector that fully characterizes ρ.The conventions used are such that the covariance matrix of the single-mode thermal state Eq. ( 7) is proportional to the identity matrix and reads σ = kI 2 = (2n + 1)I 2 .Furthermore, t is a diagonal matrix defined as Consequently, the t−PQD of a Gaussian state ρ is nonnegative iff If this condition is satisfied we will say that ρ belongs to the set of t−classical Gaussian states, which we denote with G .For the input state m j=1 S(r) |0⟩, the above condition simplifies to Let us also consider noisy threshold photo-detectors characterized by sub-unit efficiency 0 ≤ η D ≤ 1 and by a dark count rate 0 ≤ p D ≤ 1.The "off" and "on" elements of the POVM associated with this measurement respectively read A close inspection of Eq. ( 26) reveals that Π 0 is an unnormalized thermal state, hence one can analytically compute the (−s)−PQDs of both POVM elements using Eq. ( 22) and prove that they are non-negative for s ≥ 1 − 2p D /η D .If we consider m identical noisy threshold detectors as described above at the output ports of our LON, then the (−s)−PQD of the m-mode measurement is simply the product of the (−s j )−PQD of the single-mode measurements, and it is clearly non-negative for If one further assumes that losses are uniform across all possible paths across the network − usually a good approximation for integrated setups − then the linear transformation of the input modes is described by L = √ η L W , where W is a unitary matrix and 0 ≤ η L ≤ 1 denotes the transmission of the interferometer.By substituting the threshold values t and s into Eq.( 20) we obtain the classical simulability condition for the noisy GBS experiment described above, i.e., The term −n(1 − η L ) on the right-hand side (r.h.s.) of the inequality above represents a finite-temperature correction that accounts for thermal effects.Notice how the latter is always negative, meaning that thermal noise has the effect of reducing the detection noise needed for a classical simulation of the sampling task to be feasible, with respect to the zerotemperature scenario.Evidently, Eq. ( 29) is automatically satisfied whenever the r.h.s.becomes negative, i.e., if We also remind the reader that the mean photon number of a thermal state is related to the environment's temperature T via where ω is the mode's frequency and k B is the Boltzmann constant.Therefore, Eq. ( 30) predicts the existence of a threshold temperature above which the sampling task becomes classically efficiently simulable even when ideal detectors are employed.This is somewhat expected, as it is well established that a noiseless boson sampling task with thermal input states leads to a classically simulable problem.Consequently, we envision a transition in the computational complexity of the task as the environment's temperature increases.
In the remaining part of this section we formalize these ideas and provide a physical interpretation of this phenomenon.In Appendix A we show that the system's losses and thermal noise effects can be fully absorbed into the initial state, while retaining a unitary evolution via an effective ideal LON.In particular, the noisy evolution given by E with L = √ η L W is equivalent to a loss-less LON described by the unitary matrix W , preceded by m identical single-mode maps F that mix each input mode with a thermal state by means of a beam splitter with transmissivity equal to η L , and finally taking the trace over the environmental degrees of freedom (See Fig. (2) for a schematic representation of the channel decomposition).Hence, the action of the map F on a generic single-mode state ρ reads Here, U BS represents the beam splitter unitary operator acting on the system's mode and the corresponding ancillary environmental mode, and we have explicitly displayed the parameter k that completely identifies the thermal state.Focusing on GBS, each input mode of this loss-less LON is fed with F(S(r) |0⟩⟨0| S † (r)), a Gaussian state whose covariance matrix reads diag{a + , a − }, with a ± = η L e ±2r + k(1 − η L ).
As the covariance matrix of the vacuum state is simply the identity matrix, it is clear that quadrature squeezing vanishes if a − ≥ 1.One then easily proves that this happens when condition Eq. ( 30) is satisfied.As a result, the threshold temperature described in Eq. ( 30) can be physically interpreted as the temperature at which genuine quantum features of the input state completely disappear, so that efficient sampling on a classical machine becomes feasible regardless of the presence of noise in the detectors.In the following section we extend this argument to a generic quantum optical experiment.

III. THERMAL CLASSICALIZATION
Let us consider a generic m-mode input state ρ undergoing noisy linear optical evolution via the quantum map E defined in the previous section, followed by a generic quantum measurement.If we further assume that losses are uniform within the interferometer, then the classicality condition Eq. ( 20) may be recast as The r.h.s. of this inequality (i.e., the threshold temperature) diverges at η L = 1, which can be understood as the loss parameter also measures the ability of the LON to couple the system with the environment.Eq. ( 33) allows us to compute the temperature for an efficient classical simulation once the input state and the POVM have been specified.Since s ≤ I m and t ≥ −I m , it is clear that the most restrictive scenario is obtained by substituting s = I m and t = −I m into Eq.( 33), resulting in or equivalently To summarize, when the condition condition Eq. ( 35) is satisfied, it enables the efficient simulation of a general noisy quantum optical experiment, our sole assumption pertaining to the model employed to describe the linear optical evolution.
We can also find a physical meaning to the threshold temperature given by Eq. (35).We have already pointed out that the system's losses and thermal noise effects can be absorbed into the initial state by applying the map F to each input mode.Furthermore, we introduce the notation F m ≡ F ⊗m for the corresponding m−mode map.
where we have exploited the Glauber P −function representation of ρ Our objective is to compute the P -function of F m (ρ), as it captures the non-classical properties of the noisy input state at study [33,34].In particular, it is well known that having a well-behaved and non-negative P −function is a necessary and sufficient condition for a state to admit a classical description, i.e., it can be expressed as a statistical mixture of coherent states.This property is also referred to as P −classicality.
As previously mentioned, the P -function is obtained by substituting s = I m into Eq.( 3), namely We then substitute Eq. ( 36) into the previous expression to obtain The trace above may be computed using standard Gaussian calculation techniques, yielding where λ = k(1 − η L ) + η L .Putting everything together, we obtain We then perform straightforward Gaussian integration and obtain namely the convolution of ρ's P −function and a Gaussian distribution.At the threshold temperature − i.e., k = 1+η L 1−η L or equivalently λ = 1 + 2η L − we have meaning that the P -function of the noiseless state ρ and that of F m (ρ) are related by a Weierstrass transform (Gaussian filter).The latter has a smoothing effect on P ρ that removes its negativities and divergencies, resulting in a full suppression of the state's genuine quantum features, as we prove in the remainder of this section.To this end, we recall that the P −function and the Q−function are related via the identity Comparing Eq. ( 43) and Eq. ( 44) allows us to establish a connection between the P −function of ρ and that of its noisy counterpart F m (ρ) It is well known that any Q−function is positive definite and well behaved, thus proving our thesis: any input state ρ becomes P −classical after interacting with a thermal state via a beam splitter with trasmissivity η L , at the threshold temperature given by Eq. (35).In particular, we find that the P −function of the noisy state F m (ρ) is proportional to the Q−function of the input state ρ, properly rescaled.It is worth noting the close connection between the temperature bounds derived here and the definition of nonclassicality depth as introduced in Ref. [35].

IV. APPROXIMATE CLASSICAL SIMULATION OF GAUSSIAN BOSON SAMPLING UNDER THERMAL NOISE
In Ref. [18] the authors investigate the classical simulability of a noisy GBS experiment whose lossy linear optical evolution is given by Eq. (A2) and whose imperfect threshold detection is described by the POVM elements Eq. ( 26) and Eq.(27).As it will be clear in the following, their approach can account for approximate sampling, thus overcoming the principal limitation of the formalism outlined in Section II.However, this advancement is achieved at the expense of generality, limiting the applicability of this methodology to a noisy sampling task as described above.In particular, it is showed that the latter may be efficiently simulated up to error ε if the following (sufficient) condition is satisfied where q D = p D η D and Θ(x) = max (x, 0) is the ramp function.If the previous inequality does not admit a solution for any 0 ≤ ε ≤ 1, then the classical simulation algorithm fails and we say that the Gaussian boson sampler has passed the non-classicality test.Building on the approach of Ref. [18], we establish the following classical simulability condition that accounts for the influence of thermal noise where k = 2n + 1 and n denotes the mean number of environmental thermal photons.Furthermore, we show that at zero temperature Eq. (47) constitutes a tighter bound than Eq. ( 46).Additionally, note that it always exists ε such that Eq. ( 47) admits a solution for values of ε such that 0 ≤ ε ≤ ε ≤ 1.
Hence, we will say that the non-classicality test is passed if ε is larger than some fixed (small) threshold.
As stated earlier, the assumption of uniform losses enables us to absorb all imperfections and thermal noise effects into the initial state state, while retaining an ideal linear optical evolution described by the unitary operator W. Each input port of this loss-less interferometer is fed with the single-mode Gaussian state which has a diagonal covariance matrix that reads diag{a + , a − } with a ± = η L e ±2r + k(1 − η L ).The probability p(n) of observing a specific measurement outcome n = (n 1 , . . ., n m ) with n i ∈ {0, 1} is thus given by where Π n = m i=1 Π ni is the POVM associated with m−mode noisy threshold photo-detection.Let us now consider a related sampling problem, where the same loss-less LON is fed with m identical t−classical Gaussian states τ .The resulting outcome probability distribution reads We can then use the classical simulability condition Eq. ( 20) − setting L = W unitary matrix and s = (1 − 2q D )I m − to deduce that sampling (exactly) from the output state W τ ⊗n W † can be efficiently performed if The idea is that when the noisy input state τ is similar enough to τ ∈ C G for some t ∈ [1 − 2q D , 1] the corresponding sampling problem can be efficiently simulated up to a small error.In the remainder of this section, we formalize this intuition.The total variational distance (TVD) between the two output probability distributions p and p is upper bounded as follows Here, F (ρ, τ ) = (Tr{ √ ρτ √ ρ}) 2 denotes the quantum fidelity and is the trace norm, which is invariant under unitary transformations acting on ρ and τ .We emphasize that the bound on the TVD in Eq. ( 52) is more stringent compared to the one presented in Ref. [18].There, the authors exploited a generalization of the Pinsker's inequality to bound the total variational distance with the Rényi relative entropy.However, the latter is an unbounded quantity, which in turn leads to a looser bound that also explains why Eq. ( 46) does not always admit a solution.On the other hand, in Eq. ( 52) we bound the TVD with a quantity that is always less or equal than one.
As any t−classical state with t ∈ [1 − 2q D , 1] leads to an efficiently simulable instance of GBS, we further minimize the TVD (i.e., maximize the fidelity) over all possible choices of τ , namely where In Appendix B we show how to analytically carry out the optimization above, following the techniques outlined in Ref. [18].Notice that, if τ is itself t−classical, then clearly fidelity is maximised (and equal to one) for τ = τ , and it is possible to efficiently simulate the sampling task exactly.In general, we obtain which we then substitute into Eq.(54) to retrieve Eq. (47), i.e. our final sufficient condition for classical simulability of noisy GBS at finite temperature .As expected, we observe that in the zero temperature limit k = 1 we obtain a bound that is more restrictive than that presented in Eq. (46).Lastly, notice that by fixing ε = 0 (i.e., sampling from the exact probability distribution of the noisy experiment), we retrieve Eq. ( 29).

V. CONCLUSIONS
Using a phase space method based on the negativity of the relevant quasi-probability distributions, we have established a sufficient condition for the efficient classical simulation of generic quantum (linear) optical experiments affected by loss and thermal noise.Our results show how finite temperature effects reduce the threshold of the system's imperfections that are sufficient for an efficient simulation of the experiment to be feasible.We then turned our attention to a GBS task employing threshold detectors, and provided a non-classicality condition in the form of an inequality involving the squeezing and noise parameters (photon loss, mean thermal photon number, detectors inefficiencies and dark count rates), that any potential candidate for an experimental demonstration of quantum advantage must satisfy.Furthermore, we showed that there exist a threshold temperature at which any sampling experiment becomes efficiently simulable, even in the presence of ideal detectors.We presented a physical interpretation of this phenomenon by establishing a connection with the vanishing of the genuine quantum features of the state.We hope that this work inspires rigorous studies on the transition occurring in the computational complexity of GBS subject to increasing levels of thermal noise.We first optimize F (τ, τ ) over all possible t−classical states τ , i.e. we find the point (n * τ , r * τ ) that maximizes the fidelity at fixed t, subject to the constraint Eq. (B4).Clearly, if τ is itself a t−classical state − i.e. r τ ≤ r 0 (τ ) − it follows that fidelity is maximised for (n * τ , r * τ ) = (n τ , r τ ) and F max = 1 for every value of t.On the other hand, when r τ > r 0 (τ ) we find that the maximum fidelity at fixed t reads sech (r τ − 1/2 ln ((2n τ + 1)/t)).One then easily no-tices that the further optimization of this expression over t ∈ [1 − 2q D , 1] is achieved by setting t = 1 − 2q D .Finally, using Eq.(B2) and Eq.(B3), and combining the two optimization regimes using the ramp function Θ, after some algebra we obtain the desired result, namely . (B8)

FIG. 1 .
FIG.1.Schematics of the noise model used throughout this work, described by the CP map E. The system's initial m-mode state ρ interacts with the environment − initialized in a thermal state − through a loss-less 2m-mode interferometer represented by a unitary operation U.At the output ports, the system's modes are measured, while the crosses represent tracing over the environmental degrees of freedom.

FIG. 2 .
FIG.2.Decomposition of the CP map E describing the noisy linear optical evolution.Under the assumption of uniform losses we can absorb the latter and thermal noise effects into the initial state by means of the quantum channel F, while retaining an ideal evolution described by the unitary matrix W .
. is part of the AppQInfo MSCA ITN which received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 956071.H.K. is supported by the KIAS Individual Grant No. CG085301 at Korea Institute for Advanced Study.MSK acknowledges the KIST Open Research Programme, Samsung GRC programme and the KIAS visiting professorship.The project was supported by the UK EPSRC through EP/Y004752/1 and EP/W032643/1.