Entanglement and squeezing of continuous-wave stationary light

Spectral components of continuous squeezed fields are entangled. In this article we review and clarify this phenomenon by analyzing systematically the relations between the correlations of modes filtered from stationary continuous fields and the cross power spectrum between the operators of the corresponding spectral components. Moreover, we study the specific spectral components that are filtered in homodyne or heterodyne detections and their entanglement properties. In particular, we establish the equivalence between two-mode squeezing variance and logarithmic negativity for the spectral components of continuous stationary fields, thereby demonstrating that the measurement of the homodyne or heterodyne spectrum is, in fact, a direct measurement of the logarithmic negativity between specific spectral modes. As an illustrative example, we apply these concepts to the analysis of entanglement in ponderomotive squeezing.


I. INTRODUCTION
Quantum optical fields are exploited in the development of a large class of new technologies which make use of quantum mechanics to push their efficiency to the limit [1]. In particular, squeezed light plays a pivotal role in the continuous variable domain [2][3][4][5]. After the first experimental demonstrations of optical squeezing, both in the continuouswave [6,7] and in the pulsed regime [8], and of the corresponding EPR entanglement [9][10][11], nowadays squeezed optical fields are routinely produced and employed in many experiments aimed at investigating the potentiality of quantum based technologies. They range, for example, from the demonstration of quantum information tasks such as quantum teleportation [12,13] and other essential elements of scalable universal quantum computation [14][15][16][17][18], to the design of high-resolution metrology applications [19][20][21][22][23], of novel spectroscopic methods [24,25], and of enhanced optical communication schemes [26].
Squeezing and entanglement are two very related concepts. In practice squeezed fields are, for example, used to produce two-mode entangled resources by, simply, mixing them on beam splitters [10,12,16]. From a more theoretical point of view squeezing variance can be used to construct entanglement criteria [27][28][29][30]. It is also well known that the spectral components of continuous-wave squeezed light are endowed with non-trivial correlations [31][32][33][34][35]. In particular, specific spectral modes of continuous squeezed fields are entangled [36][37][38][39] realizing EPR spectral beams that have been proposed as convenient quantum communication channels [37,40].
In this article we study squeezed continuous fields in the stationary regime, and we analyze the entanglement properties of the corresponding spectral components. We aim at establishing a direct connection between the entanglement theory of continuous-variable systems and the spectral properties of squeezed light fields in the stationary continuous-wave regime. The spectral modes of a continuous field can be operatively defined as the temporal modes filtered from the total field, with a long time filter. Their entanglement and squeezing properties are therefore readily defined as the long time limit of those that are found for finite temporal modes. By employing this approach, we derive general conditions for entanglement and squeezing between two spectral components of stationary continuous fields, and we show that their twomode squeezing variance can be expressed in terms of the corresponding logarithmic negativity. We also discuss the properties of the specific spectral modes that are probed with homodyne and heterodyne detection [41,42], and we establish that the squeezing spectrum that can be measured with these techniques can be interpreted as a direct measurement of the logarithmic negativity between specific spectral modes.
Finally, we apply these ideas to the analysis of ponderomotive squeezing [43][44][45][46], namely the squeezing that is obtained as a result of the optomechanical interaction between a mechanical resonator and the light in an optical cavity. Optomechanics provides a novel approach to quantum non-linear optics which, recently, has attracted much attention for its potential applications in quantum enhanced technologies [47,48]. In this work, we investigate a two-sided cavity with a membrane in the middle and we identify the spectral components of the output fields that exhibit larger entanglement and that are experimentally accessible with homodyne and heterodyne techniques.
Part of this article comprises a review of already known results, however, rephrased with the intent of providing a clear and complete introduction to the scope of our research. In details, the article is organized as follow. In Sec. II we review the basic properties of continuous fields and of their spectral properties. We introduce the filtered temporal modes and study the correlations of the corresponding field operators in the stationary regime. In Sec. III we demonstrate the equivalence between logarithmic negativity and two-mode squeezing variance of the spectral components of stationary continuous fields. Then, in Sec. IV we review homodyne and heterodyne detection techniques and we study how they can be used to directly measure the logarithmic negativity between spectral modes. Finally, in Sec. V, we apply the concepts developed in the preceding sections to ponderomotive squeezing. The three appendices provide additional informations regarding, respectively, the basic properties of entanglement and squeezing of discrete bosonic modes, the homodyne and heterodyne techniques, and the input-output theory applied to the investigation of an optomechanical system.

II. CONTINUOUS QUANTUM OPTICAL FIELDS
In this section we introduce the objects of our investigation, namely continuous fields, and we discuss the properties of temporal modes that can be filtered from them [49][50][51]. In particular we define the operators for the spectral components in terms of the modes that are filtered with a long-time filter and that describe narrow bands of frequencies. These operators are particularly suited for the study of the entanglement properties of the spectral components of continuous field using standard techniques of entanglement theory.
In detail, we investigate the freely propagating continuous field E(t) at the output of a quantum optical system. It can be decomposed into the positive and negative frequency components , where σ is the cross section of the propagating field, and a(ω) is the annihilation operator for the spectral component at frequency ω, which satisfy the standard commutation relation a(ω), {a(ω ′ )} † = δ(ω − ω ′ ). In general, in an optical system, only a relatively narrow band of frequencies ∆ ω is relevant. This band is centred around the carrier frequency ω L of the signal field E(t), which is typically defined by the frequency of a laser driving the system, and fulfils the relation ∆ ω ≪ ω L . In practice the relevant bandwidth is set by the typical linewidth Γ of the system under investigation, and eventually by the response time, T , of the detector such that {Γ, 1/T } ≪ ∆ ω . Under these conditions the range of frequency integration in the expression for the field can be extended from −∞ to ∞, and the relevant wave numbers k can be approximated with the central value k ∼ ω L /c = k L . By this means the quantum optical continuous field can be expressed as where we have introduced the continuous field annihilation operator a(t). It is related to the operators for the spectral components by the Fourier transform a(t) = where, here, ω is the frequency relative to the carrier. It is also useful to introduce the operators a(ω) = a(ω L + ω) e −iω L t relative to the carrier frequency, which are equal to the Fourier transform of a(t) and { a(ω)} † = a † (−ω). These operators satisfy the standard commutation relation for continuous fields

A. Filtered modes and spectral components of the field
In reality one has access only to a finite time interval (and correspondingly to a finite band of frequencies) of the total field a(t). These detectable intervals of the total field correspond to specific temporal modes. They can be physically defined, for example, by the temporal profile of the pumping field in pulsed experiments [8,10,[52][53][54], they can also be extracted by post-processing the previously recorded time signal [55,56], or they can be selected by the measurement apparatus as a result of the corresponding response and detection times [56][57][58][59]. In particular, in the case of experiments involving stationary fields, the detection time can be so long to select a well defined spectral component of the total signal (as realized, for example, with an electronic spectrum analyzer) [6,7,42].
In general a temporal filtered mode can be introduced in terms of a filter function φ τ (t), which defines the time profile of the mode with a duration of order τ. Correspondingly, it defines a band of spectral components, of width 1/τ, that are combined into the filtered signal [60,61]. The generic form for the operators of a filtered mode can be expressed as and where the symbol indicates filtered quantities. The parameter Ω defines the central frequency of the filter, and the filter function φ τ (t) is real and normalized according to Consequently, the filtered operators are discrete bosonic operators which satisfy the standard commutation relation The corresponding equivalent form of the filtered operators, in terms of the spectral components of the field, is where φ τ (ω) is the Fourier transformed filter function φ τ (ω) = 1 √ 2π ∞ −∞ ds e iωs φ τ (s), that is peaked at ω = 0 and has a width of the order of 1/τ. Two particular cases are worth mentioning. The exponential filter with φ exp which has been used, for example, in Ref. [57] to introduce the physical spectrum of light, and the step filter function which we will connect to homodyne and heterodyne detection techniques in the following. In this form the time t in the filtered operator a τ (Ω, t) corresponds to the final time of the filtering process, that is a τ (Ω, t) = t −∞ ds e iΩs φ τ (t − s) a(s) .
Although not strictly relevant for the results presented in this article, this choice is physically motivated by the fact that in this way we define, at time t, a causal operator a τ (Ω, t) which depends only on the past of the continuous field a(t) [60].
In the limit of long filtering times, τ → ∞, the filter selects a single spectral component of the field. In the following we will focus on the spectral components of stationary fields for which we will use the following simplified notation where we drop the label τ, the limit symbol, and the time argument t. In particular, the time t in this operators, is irrelevant for stationary fields, in the sense that (as shown in the next section) the correlations of operators of this form, in the limit of large τ, are independent from the time arguments. We finally note that, we are describing the field in a reference frame rotating at the carrier frequency ω L , therefore a (Ω) is, in fact, the operator for the spectral component of the field at the sideband frequency ω L + Ω.

B. Correlations of filtered spectral modes of stationary fields
In this article we are interested in the squeezing properties of the electromagnetic field, which refers to reduced fluctuations or reduced variance of specific quadratures below the vacuum noise level, and in the corresponding entanglement features. Squeezing can be revealed form the analysis of the second order correlations of fields operators. Therefore in this section we analyze the basic properties of the correlations of filtered modes. In particular we focus on the spectral properties of stationary fields. The corresponding field operators, a(t) and a(ω), have diverging correlation functions, as a consequence of the commutation relations in Eq. (3). In this case is therefore instructive to analyze the fluctuations of the spectral modes in terms of the filtered spectral operators defined in Eq. (6), whose correlations, on the contrary, are always finite, also in the limit of long integration time τ → ∞. This approach is particularly useful for the study of the corresponding entanglement properties, and it has the advantage to provide a clear physical definition of discrete modes corresponding to the specific spectral components, hence allowing for a transparent application of the techniques developed in entanglement theory, which indeed deals with discrete modes (see App. A).
Let us study the correlations between the spectral components of two stationary continuous fields with annihilation operators a 1 (t) and a 2 (t) respectively, which fulfil the commutation relation a j (t), a k (t ′ ) † = δ j,k δ(t − t ′ ) (the same results that we discuss below for two continuous fields can be applied with minor modifications to different spectral components belonging to a single field). By straightforward application of the Wiener-Khintchine theorem, we find that all the informations about the correlations between the spectral components are contained into the power spectrum matrix P(ω), defined as the Fourier transform of the two-time stationary correlation matrix, which can be expressed in terms of the elements of the column vector of operators a(t) = a 1 (t), a 2 (t), a † 1 (t), a † 2 (t) T , as the matrix A(t) = a(t) a(0) T , whose elements are {A(t)} j,k = {a(t)} j {a(0)} k , where j, k ∈ {1, 2, 3, 4} are vector indices, not to be confused with the indices of the modes. To be specific where we use the fact that two-time correlation functions of stationary signals depends only on the difference of the time arguments. In particular the correlations between the Fourier transformed operators a(ω) = 1 √ 2π ∞ −∞ dt e iωt a(t) are diverging and are related to the power spectrum matrix by In order to gain insight into the physical meaning of these diverging quantities we employ the narrow filtered modes introduced in the previous section. We construct the vector of filtered spectral components a(Ω) = lim τ→∞ ∞ −∞ ds e iΩs φ τ (t − s) a(t), that is given by a(Ω) = a 1 (Ω), a 2 (Ω), a † 1 (Ω), a † 2 (Ω) T , and we compute the corresponding matrix of correlations These quantities can be evaluated by noting that for large τ, the square modulus of the filter function approaches a delta function, lim τ→∞ φ τ (ω) 2 = δ(ω), while its integral goes to zero, lim τ→∞ dω φ τ (ω) = 0. And, correspondingly, given a generic finite function f (ω), the relation holds. Consequently Eqs. (10) and (12) can be used in Eq. (11) to find that shows that the power spectrum is directly related to the correlations of narrow filtered modes. In other terms, in the limit of large integration time τ, i.e. when the bandwidth selected by the filter is sufficiently small, the correlation functions reduce to the power spectrum of the continuous field [57]. This result is valid when τ is much larger than the decay time of the signals correlations τ C (the memory time of the signals), τ C ≪ τ. We note in particular that this relation implies the stationarity of the signal which is reached on a time scale of the order of τ C . The correlations between two modes are conveniently analyzed in terms of the corresponding correlation matrix, from which the corresponding squeezing and entanglement properties can be readily derived (see App. A for a short review). We remark, however, that the matrix P(Ω) is not a correlation matrix for two modes. In fact, it contains the correlations between the four spectral modes corresponding to the two pairs of sidebands at frequency ±Ω of the two continuous fields. On the other hand, the elements of the power spectrum matrix can be used to construct the correlation matrix for two narrow modes filtered from the two stationary continuous fields as follows. We consider two modes, at frequencies Ω and Ω ′ , respectively described by the filtered operators where the narrow bandwidth limit (τ → ∞) is implicit in their definition [see Eq. (8)]. The corresponding correlation matrix is defined using the vector a(Ω,

III. EQUIVALENCE BETWEEN TWO-MODE SQUEEZING VARIANCE AND LOGARITHMIC NEGATIVITY OF THE SPECTRAL COMPONENTS OF STATIONARY CONTINUOUS FIELDS
Having, in the previous section, introduced our notation, and reviewed the basic properties of stationary continuous fields, we are now in the position to study the general conditions for the squeezing and the entanglement between two spectral components, that can be inferred from Eq. (15). In particular, in the case of Gaussian states, we establish the equivalence between the logarithmic negativity and the twomode squeezing variance of two spectral modes.
In general, given two modes described by the operators a 1 and a 2 , two-mode squeezing is characterized by non vanishing correlations of the form a 1 a 2 and a † 1 a † 2 . According to Eq (15), the correlation between the annihilation operators for two filtered spectral components of stationary fields, a 1 (Ω) and a 2 (Ω ′ ), can be non-vanishing only for opposite frequencies, that is when Ω = −Ω ′ . In this case the matrix in Eq. (15) reduces to the form where with n ± real and positive. Here we have used the general properties of the power spectrum matrix P(Ω) − P(−Ω) T = , where 1 2 is the 2 × 2 identity matrix and the missing blocks are null matrices, and P(Ω) * Eq. (16) represent the general form for the correlation matrix between two spectral components at opposite sideband frequencies of stationary continuous fields. This matrix can be exploited to derive general results regarding the corresponding squeezing and entanglement properties.
In general squeezing refers to the reduced fluctuations of field quadratures. Let us therefore define the quadrature operators for a spectral mode Thereby, we can formalize the condition for two-mode squeezing of the two spectral components as follows. The two components are two-mode squeezed when the variance of a generic composite quadrature of the form is below the shot noise level for some value of θ ± (see App. A for further remarks). We first note that, in the case of a stationary field, for which a(t) = α is constant, the average of a corresponding filtered mode is given by a τ (Ω, t) = √ 2π α e i Ω t φ τ (−Ω), and it approaches zero for large τ. Therefore, according to our definitions, the fields are two-mode squeezed if for some values of θ ± the autocorrelation function of the combined quadrature is smaller then one, ∆X (θ ± ) (Ω) < 1. More in general one can construct composite quadratures with different weights of the two components It has been shown [29] that the variance of quadratures of this form can be used to define entanglement criteria. Specifically when the relation is satisfied for some values of θ ± and ξ ± , then the two modes are entangled. In general this is a sufficient condition for entanglement, but it is also necessary in the case of Gaussian fields and for an appropriate choice of ξ ± . The calculation of the autocorrelation function of these composite quadratures is straightforward using the matrix of correlations in Eq. (16). The result is In particular, we find that (Ω), hence, in our case, the condition for entanglement reduces to ∆X (θ + ,θ − ) ξ + ,ξ − (Ω) < 1. The corresponding optimized squeezing spectrum can be defined as the minimum of this quantity over the quadrature of the field. Specifically we can identify two different minimization strategy. If we restrict to composite quadrature that are symmetric superposition of the two components (ξ + = ξ − ) as in Eq. (19) then the minimization runs only over the phases θ ± , and the corresponding phaseoptimized squeezing spectrum takes the general form This is the quantity that is obtained, for example, by the homodyne measurement of a continuous field, where the phases θ ± , in Eq. (22), are directly related to the phase of the local oscillator (see Sec. IV for further details). If, on the other hand, we consider the more general quadratures with the two components scaled by factors ξ ± as in Eq. (20), then the minimization can be performed both over the phases and over the parameters ξ ± , and the corresponding globally-optimized squeezing spectrum reduces to In general S min (Ω) ≤ S (Ω), and they are equal in the case of symmetric spectral components, for which n + = n − . In the next section we will describe how to measure both quantities in few specific cases with homodyne and heterodyne techniques. Here we emphasize that the occurrence of S min (Ω) < 1 implies that the entanglement criterion in Eq. (21) is satisfied and, in turn, it entails that squeezing spectrum smaller the one is always a signature of entanglement. We further note that Eq. (24) is, indeed, smaller then one if and only if Consequently this relation can be interpreted as a sufficient condition for the entanglement between the spectral components at opposite sideband frequencies of stationary continuous fields. And, as already noted, it is also a necessary condition in the case of Gaussian fields.
Let us now focus to the Gaussian regime. In this case the logarithmic negativity, that is a measure of bipartite entanglement [62], can be expressed as where the parameter ν is equal to the smallest symplectic eigenvalues of the covariance matrix corresponding to the partially transposed state of the two modes (see App. A). Using Eq. (16) we find that the parameter ν evaluated for each pair of spectral components at the sideband frequencies ±Ω is equal to the squeezing spectrum in Eq. (24), This is a general result valid for the spectral components of stationary continuous fields (a specific example has been discussed in Ref. [63]). In particular, this relation implies that, if the two stationary fields are Gaussian, then a measure of the minimum variance of a composite quadrature of two spectral components, of the form of Eq. (20), is a direct measurement of the corresponding logarithmic negativity.

IV. HOMODYNE AND HETERODYNE DETECTION OF THE SPECTRAL COMPONENTS OF STATIONARY CONTINUOUS FIELDS
The quadratures of continuous electromagnetic fields are routinely measured in experiments with homodyne and heterodyne techniques [5,41,[64][65][66][67]. The photocurrents resulting from homodyne and heterodyne detections are, in fact, proportional to specific quadratures of the detected field. In turn, the power spectrum of the photocurrent, namely the homodyne or heterodyne spectrum, measures the fluctuations of the quadratures at specific frequencies. Such spectra are therefore directly related to the squeezing and entanglement properties of the spectral components of the electromagnetic field, and in particular to the squeezing spectra defined in Eqs. (23) and (24). Specifically, we will show that the autocorrelation function of the photocurrent minimized over experimentally accessible parameters as for example the phase of the local oscillator can be always cast in the form of Eqs. (23) or (24), with corresponding parameters n ± and m evaluated for specific spectral modes. This justifies the interpretation of the optimized homodyne and heterodyne spectra as a direct measurement of the logarithmic negativity of these modes.

A. Single-mode homodyne spectrum and entangled spectral components
In homodyne detection the signal field is mixed on a 50:50 beam splitter with a strong monochromatic field (the local oscillator) at the same frequency of the carrier signal. The fields at the two output ports of the beam splitter are detected and the corresponding photocurrents are subtracted resulting in a signal which contains informations about a field quadrature [65], and that can be described by a photocurrent operator of the form (see App. B) where θ is the phase of the local oscillator. The power spectrum of the photocurrent contains informations about the spectral components of the detected field, and in particular it quantifies the strength of the fluctuations at specific frequencies. We will refer to it as the homodyne spectrum, and it can be expressed as the autocorrelation function of the filtered photocurrent, integrated over a long time τ, of the form In detail, the homodyne spectrum can be written as We note that, for stationary processes, this quantity is independent from the phase of the filter ϕ (see App. B). However, this phase is relevant and can be useful when considering combinations of filtered photocurrents at different phases which, as discussed below, can be exploited to probe arbitrary superpositions of spectral modes. Moreover the same results for the power spectrum in Eq. (29) is obtained when, in the filtered photocurrent J (θ,ϕ) τ (ǫ, t), one uses an exponential oscillating function, as in the filters of Sec. II A, in place of the sinusoidal function introduced above. The use of a sinusoidal function is, however, more convenient because, in this case, the filtered photocurrent is an hermitian operator thereby making the relation between the photocurrent and the field observables more transparent. Specifically the filtered photocurrent in the limit of long filtering time can be expressed as the sum of two filtered quadrature operators for the two spectral modes at frequencies ±ǫ (see App. B) Similarly to Eq. (19), it is a symmetric superposition of two quadratures, defined as in Eq. (18), corresponding to the two spectral components, which in this case are filtered from the same field, and whose annihilation and creation operators are a(ǫ) and a(−ǫ) , as illustrated in Fig. 1. The corresponding single-mode homodyne spectrum (where single-mode indicates that it results form the detection of a single continuous field) is equal to Eq. (22) with ξ + = ξ − and θ + + θ − = 2θ, i.e. G (θ) (ǫ) = 1 + n (I) Thus, the single-mode phase-optimized squeezing spectrum, that is experimentally accessible by tuning the phase of the local oscillator, S (I) (ǫ) = min θ G (θ) (ǫ), is equal to Eq. (23) evaluated for the parameters in Eq. (32). When it is smaller then one, it indicates that the two sideband modes a(±ǫ) are entangled [37][38][39][40]. Their logarithmic negativity, that as discussed in Sec. III is directly related to the squeezing spectrum in Eq. (24), could be measurable if one could construct a filtered photocurrent similar to Eq. (20), which is a nonsymmetric superposition of the quadratures of the two modes. This photocurrent is, in fact, realizable by combining two filtered photocurrents evaluated at appropriately tuned phases of the local oscillator θ, and of the filter ϕ. Specifically one should first detect the filtered photocurrent J (θ,ϕ) (ǫ) for some value of θ and ϕ, and then a second one J (θ ′ ,ϕ ′ ) (ǫ), with the phases tuned to different values θ ′ and ϕ ′ . The two photocurrents are then summed resulting in the total photocurrent that we have appropriately normalized, and where We note that Eq. (33) has, indeed, the form of the composite quadrature defined in Eq. (20). This can be experimentally realized, for example, by recording the photocurrent for a sufficiently long time and then post processing the recorded signal.
The corresponding single-mode globally-optimized squeezing spectrum S (I) min (ǫ) = min θ ± ,ξ ± J (θ + ,θ − ) ξ + ,ξ − (ǫ) 2 is then equal to Eq. (24) evaluated for the parameters in Eq. (32), and it can be measured by minimizing the homodyne spectrum over the phases of both the local oscillator and of the filter.

B. Two-mode homodyne spectrum and entangled spectral components
In the preceding section we have seen that the single-mode homodyne spectrum, which is measurable from a single stationary continuous field, provides informations about the entanglement between two spectral components. Let us now study the two-mode squeezing spectrum obtained from the combination of two homodyne photocurrents which result from the measurement of two signal fields a 1 (t) and a 2 (t), as depicted in Fig. 2. This strategy detects the correlations between four spectral components [68]. Specifically, we consider the situation in which the photocurrents corresponding to the two detected fields, I (θ 1 ) 1 (t) and I (θ 2 ) 2 (t) have the form of Eq. (28), and are combined to construct the total photocurrent where we have introduced the scaling parameters µ 1 and µ 2 , which weight differently the two quadratures, and hence provide a means to select arbitrary collective modes as discussed below. These parameters are controllable experimentally including asymmetrical amplification and/or attenuation of the two photocurrents. The total photocurrent is then analyzed in frequency. Similarly to the previous case, the filtered pho- c (s) can be decomposed into two spectral components at the frequencies ±ǫ where θ ± = θ 1 +θ 2 2 ± ϕ, and where we have introduced the quadrature operators for the collective spectral modes defined as FIG. 2: Two-modes homodyne detection: Two stationary continuous fields are detected by homodyne techniques. The two photocurrents are summed and then analyzed in frequency. As in the single-mode homodyne detection, the resulting total filtered photocurrent is a superposition of spectral modes at opposite sideband frequencies ±ǫ. However, here, each spectral mode [c(±ǫ)] can be decomposed as the superposition of two spectral components [a 1 (±ǫ) and a 2 (±ǫ)], at the same sideband frequency, each filtered from one of the two fields.
with the collective annihilation and creation operators given by where c(±ǫ), c † (∓ǫ) = 1 and θ c = (θ 1 − θ 2 ) /2. Also in this case, the filtered photocurrent is a composite quadrature of the form of Eq. (30). Thus, although it is constructed form the collective operators in Eq. (38), we can still apply the results of Sec. III to conclude that the two-mode squeezing spectrum, S (II) (ǫ), obtained as the minimum over the phases θ ± of the autocorrelation function of the filtered current, has the form of Eq. (23) but evaluated with the parameters µ 2 1 e i(θ 1 +ϕ) w (11) + µ 2 2 e −i(θ 1 +ϕ) w (22) +µ 1 µ 2 e i(θ 2 +ϕ) w (12) + e −i(θ 2 +ϕ) w (21) , where v ( jk) for j, k = 1, 2. We remark that S (II) (ǫ) is found by minimizing the autocorrelation function of the photocurrent in Eq. (36) over θ ± only, while θ c and µ j define the specific collective modes that are being probed and are fixed. Therefore, in this case the condition S (II) (ǫ) < 1 indicates the entanglement of the collective modes described by the operators in Eq. (38), each of which is a superpositions of the spectral components at the same sideband frequencies of the two fields. Moreover similarly to the discussions of the single-mode homodyne spectrum, the corresponding logarithmic negativity can be, in turn, measured by summing two filtered photocurrent, at different phases, of the form of Eq. (36), and then calculating the corresponding autocorrelation function. The two-mode squeezing spectrum S (II) min (ǫ) is then found by minimizing this quantity over both the phases of the local oscillator and of the filter, and the result, in full similarity with the single-mode squeezing spectrum, is equal to Eq. (24), but now evaluated for the parameters in Eq. (39).

C. Detecting single spectral components with homodyne techniques
As discussed in the previous sections it is possible to construct arbitrary superposition of spectral modes at opposite sideband frequencies by the superposition of two homodyne photocurrents. The corresponding total photocurrent is then given by Eq. (33). Similarly, when the phases in Eq. (34) are set to some values for which one of the two parameters ξ ± is equal to zero, then a single spectral mode is detected.
When applied to two distinct fields, this approach would permit the investigation of the correlations between two distinct spectral modes each belonging to a different field. Let us, for example, assume that we perform the type of measurement resulting in the photocurrent in Eq. (33) on two different continuous fields a 1 (t) and a 2 (t). The two resulting filtered photocurrents are then given by J (θ +, j ,θ −, j ) ξ +, j ,ξ −, j (ǫ) = ξ +, j x (θ +, j ) (ǫ) + ξ −, j x (θ −, j ) (−ǫ) / ξ 2 +, j + ξ 2 −, j with the parameters defined as in Eq. (34), and where, here, j = 1, 2 distinguish the parameters corresponding to the measurement of the first and of the second field respectively.
The two photocurrents are then summed together, after being multiplied by appropriately chosen scaling factors ζ ± , so that the resulting total photocurrent is where the single mode quadratures x (θ) j (±ǫ) are defined in Eq. (18). Thus, this protocol detects the combined quadrature defined in Eq. (20). The corresponding squeezing spectrum, S (III) (ǫ), obtained, for ζ + = ζ − , and from the minimization of the autocorrelation function of the total photocurrent over θ +,1 and θ −,2 is equal to Eq. (23). Similarly, S (III) min (ǫ), obtained as the minimization of the power spectrum of the total photocurrent over θ +,1 , θ −,2 and ζ ± is equal to Eq. (24). In particular, also in this case we can conclude that this quantity is equivalent to the logarithmic negativity between a 1 (ǫ) and a 2 (−ǫ) when the fields are gaussian.

D. Two-modes heterodyne spectrum and entangled spectral components
An alternative strategy to probe single spectral modes and hence to measure the squeezing spectrum, and the logarithmic negativity between two distinct spectral components of two distinct fields, as in Sec. IV C, is based on heterodyne measurements.
In heterodyne techniques, the local oscillator is detuned from the carrier frequency of the signal field by a quantity ∆ = ω LO − ω L , and the corresponding operator for the photocurrent reads Hence, the corresponding filtered photocurrent, J (θ,ϕ) , is the superposition of the quadratures for the spectral components at the frequencies ∆ ± ǫ, namely J (θ,ϕ) Thus, while homodyne detection probes spectral components at opposite sideband frequencies, heterodyne techniques measure two components asymmetrically located with respect to the carrier frequency of the signal field, but symmetric with respect to the local oscillator frequency, that, in our description (where all frequencies are relative to the carrier signal), is equal to ∆ (see Fig. 3).
In particular heterodyne techniques can be used to detect a single spectral component, as discussed in the following. Say we want to detect the component at frequency Ω, then we set the detuning ∆ at a value much larger then the typical bandwidth of the signal |∆| ≫ ∆ signal , that is the band of frequencies that are populated by the signal photons. We also assume that the frequency Ω is a relevant frequency for the field |Ω| ≤ ∆ signal . Then the corresponding heterodyne photocurrent is filtered at the frequency ǫ = ∆ − Ω, so that, as depicted in Fig. 3, J (θ,ϕ) ∆ (ǫ) is the superposition of the field quadratures at the frequencies Ω and 2∆ − Ω Since the signal covers a bandwidth much smaller then ∆ then the mode at 2∆−Ω is basically in vacuum and only the photons of the sideband Ω are detected. However in doing this the vacuum fluctuations of the empty component at 2∆ − Ω are added to the signal resulting in higher noise. This approach can be exploited to measure the correlations between two spectral modes belonging to two separate fields. Specifically the two-modes heterodyne spectrum is obtained when detecting two signals with two heterodyne measurement (with ∆ ≫ ∆ signal ). The corresponding photocurrents are filtered independently at the frequencies ǫ 1 and ǫ 2 respectively, and then combined, with appropriately chosen scaling factors ξ j , in order to construct the total filtered photocurrent J ∆,tot ∝ ξ 1 J (θ 1 ,ϕ 1 ) . If, in particular, we are interested in the spectral components at frequency Ω of the first field, and at frequency −Ω of the second, such that |Ω| ≤ ∆ signal , we consider the filtered photocurrents at the frequencies ǫ 1 = ∆ − Ω and ǫ 2 = ∆ + Ω. Similarly to Eq. (42) they are equal, respectively, to the superpositions of the two quadratures at frequencies Ω and 2 ∆ − Ω of the first field, and of the two quadratures at frequencies −Ω and 2 ∆ + Ω of the second. Correspondingly, the total photocurrent is given by where the quadrature operators for a single component are defined in Eq. (18). Its autocorrelation function is therefore given by where the term 1 2 in the right hand side is due to the vacuum fluctuations of the spectral modes at 2∆ ± Ω, the collective quadrature X (θ 1 ,θ 2 ) ξ 1 ,ξ 2 (Ω) has the same form of the one defined in Eq. (20), and its autocorrelation function, that is given in Eq. (22), quantifies the correlations between the modes a 1 (Ω) and a 2 (−Ω).
Also in this case we define two kinds of optimized squeezing spectra. One is obtained by minimizing the autocorrelation function in Eq. (44), with ξ 1 = ξ 2 , only over the phases of the local oscillators T (Ω) = min θ j J ∆,tot 2 ; the other is obtined when the minimization runs also over the scaling parameters, T min (Ω) = min θ j ,ξ j J ∆,tot 2 . In both cases they can be expressed in terms of the squeezing spectra resulting form the protocol described in Sec. IV C, as Thereby, according to Eq. (27), in the gaussian case, T min (Ω) measures the entanglement between the modes whose operators are a 1 (Ω) and a 2 (−Ω).

V. APPLICATION TO PONDEROMOTIVE SQUEEZING IN A TWO SIDED CAVITY
Here we study ponderomotive squeezing and the conditions under which the spectral components of the field emitted by an optomechanical system are squeezed and entangled. Moreover we determine the squeezing spectra, and we identify the detectable spectral modes that exhibit maximum entanglement.
In detail, we investigate a Fabry-Perot cavity with a membrane in the middle [70]. A single optical mode is relevant in the system dynamics. It loses photons at rates κ j from the mirrors j = 1, 2, and is driven by a laser at a frequency detuned by δ from the relevant cavity resonance. Only one mechanical mode of the membrane at frequency ω m interacts significantly with the cavity field with a linearized coupling strength g. The decay rate of the membrane is γ, and the number of thermal mechanical excitations n T . The corresponding linearized optomechanical dynamics is Gaussian [47,48] and is efficiently analyzed in terms of the standard input-output theory [69]. Here we describe the results for the field emitted through the two cavity mirrors at the steady state of the system dynamics in the regime of optomechanical stability, referring the reader to the App. C for further details and derivations.
The photons lost through the two cavity mirrors are described by the output field operators a out 1 , a out 1 † , a out 2 , and a out 2 † . The corresponding power spectrum matrix, defined in Eq. (9), can be evaluated for the vector of operators a out = a out 1 , a out 2 , a out , and the result is given by where Q out is a diagonal matrix whose diagonal elements are √ 2κ 1 , √ 2κ 2 , √ 2κ 1 , √ 2κ 2 , Y is a matrix whose only non zero elements are {Y} 1,3 = {Y} 3,4 = 1, and the matrix W is given in terms of the parameters (with α complex even function of ω, and β ±ω real and positive) as This matrix contains all the informations about the spectral properties of the output fields and can be used to construct the correlation matrix for two spectral modes as in Eq. (15). In the following we will use this matrix and the results of sections III and IV to study the corresponding entanglement properties.

A. Homodyne and heterodyne spectra and entangled components of the emitted field
In Sec. IV we have described three different strategies for the experimental investigation of the spectral properties of stationary continuous fields, which are based on homodyne and heterodyne techinques. They probe different pairs of spectral modes given, respectively, by Eqs. (31), (38) and (41). When applied to the investigation of the optomechanical system, these techniques allow the detection of the corresponding spectral components of the two output fields a out j (t) as depicted in Fig. 4. In particular, using these techniques it is possible to study composite quadratures of these pairs of modes and their squeezing and entanglement properties. We have identified two different kinds of optimized squeezing spectra, corresponding to different experimental approaches for the measurement and the minimization of the homodyne photocurrent fluctuations. Specifically, in order to probe symmetric superpositions of quadratures of the two modes it is sufficient to apply standard homodyne techniques, thereby the minimization is realized tuning the relative phase of the two quadratures, which is controlled experimentally by the phase of the local oscillator. We indicate this phase-optimized spectrum with the symbol S (ℓ) (ω) where the label ℓ = I, II, III is used to distinguish the three detection strategies (see Fig. 4). On the other hand, non-symmetric superpositions, with different weights of the two quadratures, can be probed by combining different filtered photocurrent detected with appropriately selected phases of both the filter and the local oscillator. The globallyoptimized squeezing spectrum S (ℓ) min (ω) is then obtained minimizing the corresponding fluctuations over both the phases of the local oscillator and of the filter.
In all cases the squeezing spectra S (ℓ) (ω) and S (ℓ) min (ω) are equal to Eq. (23) and (24), evaluated in each case for the specific parameters n ± and m which correspond to the detected spectral modes. Specifically, S (I) (ω) and S (I) min (ω) are obtained by single-mode homodyne detection of a single field (either one of the two output fields), as discussed in Sec. IV A, and if applied to the output from the first mirror then n (I) ± = a out 1 † (∓ω) a out 1 (±ω) and m (I) = a out 1 (ω) a out 1 (−ω) . The second strategy is based on the two-mode homodyne detection of the two output fields (see Sec. IV B) and the corresponding spectra, S (II) (ω) and S (II) min (ω), are evaluated for , where c (out) (±ω) have the same form of Eq. (38) but, in this case, they are constructed as the superpositions of the two filtered output fields a out j (±ω). Finally, S (III) (ω) and S (III) min (ω), correspond to the two-modes heterodyne detection of the two output fields as discussed in Sec. IV D. If it is applied to the spectral component at frequency ω of the first output and at −ω of the second then, S (III) (ω) and S (III) min (ω) are evaluated for n (III) ± = a out 1 † (∓ω) a out 2 (±ω) and m (III) = a out 1 (ω) a out 2 (−ω) . We remark that these last two spectra can also be retrieved by combining various homodyne photocurrents as discussed in Sec. IV C.
The power spectrum matrix in Eq. (46) can be used to find In the case of the strategy (II), the parameters µ j and θ c , in the expression for q (II) ± , determine the specific detected composite modes, that are defined as in Eq. (38). We observe that when θ c = 0 and µ 1 /µ 2 = κ 1 /κ 2 (51) then q (II) ± = κ 1 + κ 2 , hence we can conclude that S (II) min (ω) evaluated for a two-sided configuration with decay rates κ 1 and κ 2 , is equal to S (I) min (ω) when it is evaluated for a single-sided cavity with the decay rate equal to the sum of the decay rates κ 1 + κ 2 of the two-sided configuration. It indicates that the entanglement between a out 1 (ω) and a out 1 (−ω), in a single-sided cavity, is redistributed between the four spectral components a out 1 (±ω) and a out 2 (±ω) in the case of a two-sided cavity. For this reason, the same amount of squeezing found in the case of a single-sided cavity, can be recovered when the informations from the two decay channels of a two-sided cavity are properly combined. In particular q (II) ± is maximum for the parameters of Eq. (51), and consequently the corresponding modes are the collective modes that are maximally squeezed (and entangled).
We also note that according to Eq. (25) we find that, in all cases, the pairs of spectral components are entangled (and squeezed), although possibly with different degree of entanglement (and squeezing), when Moreover we remark that the logarithmic negativity between the pair of spectral components detected with each strategy is obtained applying the definition in Eq. (26) to the parameter ν (ℓ) (ω) = S (ℓ) min (ω) , which is equal to the minimum symplectic eigenvalue of the corresponding partially transposed covariance matrix.
In Fig. 5 we compare the results for the spectra evaluated for realistic parameters and corresponding to the three detection strategies, and hence to different pair of spectral components. Each plot in Fig. 5 is evaluated for different values of the relative decay rate κ 2 /κ 1 of the two mirrors, while the total decay rate κ 1 + κ 2 and all the other parameters are kept fixed. In plot (a) we study a single sided cavity with κ 2 = 0. In plot (b) the mirrors are lossy and non-symmetric, and finally in (c) the two mirrors are symmetric κ 1 = κ 2 . When κ 2 = 0, in plot (a), the curves for S (II) , ν (II) , S (III) and ν (III) correspond to the situation in which the detector on the second mirror detects only vacuum fluctuations. It is therefore clear that the curves for S (III) and ν (III) that measures the correlations between the single spectral component a out 1 (ω) of the field lost form the first mirror and the single spectral component a out 2 (−ω) of the field lost from the second show no squeezing. On the other hand maximum two-mode squeezing and entanglement, is observed for the single-mode homodyne spectra corresponding to the curves S (I) and ν (I) that measure the correlations between a out 1 (ω) and a out 1 (−ω). The curves S (II) , ν (II) , that measure the correlations between the composite modes in Eq. (38), are at an intermediate value as a result of the vacuum fluctuations of the modes a out 2 (±ω) which reduces the visibility of the two-mode squeezing between a out 1 (ω) and a out 1 (−ω). Plot (c) corresponds to a symmetric two-sided cavity with κ 1 = κ 2 . In this case the maximum squeezing and entanglement is obtained for the curves S (II) and ν (II) , indicating that in this configuration the maximally entangled spectral components correspond to the composite modes in Eq. (38). In particular S (II) and ν (II) are equal to the curves for S (I) and ν (I) in the single-sided cavity reported in plot (a). In the two cases, the correlations of the cavity field built by the optomechanical dynamics are equal, as a result of the equal total decay rate. The only difference is that, in one case all the photons are lost through a single mirror, while in the other the they are split in the two decay channels and only their combined detection can reveals the corresponding total degree of squeezing and entanglement. Moreover in plot (c) we observe that S (III) and ν (III) are equal to S (I) and ν (I) . In this case, the two output fields are symmetric, hence the correlations between a (out) 1 (ω) and a (out) 1 (−ω) are equal to that between a (out) 1 (ω) and a (out) 2 (−ω). Finally, plot (b) corresponds to an intermediate situation between the two described in (a) and (c), and shows that the three detection strategies can display different degrees of squeezing. In general, the values of all the squeezing spectra can lay at any value between the extremes, set by the corresponding curves in plot (a) and (c), depending on the actual value of the ratio κ 2 /κ 1 ∈ [0, 1].
We also emphasize that the values of S (II) and ν (II) are reported, in the three plots, for the same values of the parameters µ j and θ c which define the specific superposition of spectral components that are being probed as defined in Eq. (38). However for each value of κ 2 /κ 1 the values of µ j and θ c can be appropriately tuned in order to find the composite modes that are characterized by the same maximum amount of squeezing in as in plot (c). This is shown for the parameters of plot (b) and at ω = 0 in Fig. 6 where maximum of squeezing and entanglement (recovering the maximum value found in Fig. 5 (c)) is found when µ 2 /µ 1 = √ κ 2 /κ 1 . We further observe that, in Fig. 5, the results for S (I) and S (II) are very close to ν (I) and ν (II) respectively. They are sensibly different only close to the mechanical resonances (ω = ±ω m ), where although the two spectral modes are entangled (ν (ℓ) < 1), this feature is not reflected in the corresponding squeezing spectrum S (ℓ) , which is sensibly larger and very close to one. This happens when the two corresponding spectral modes are significantly asymmetric so that n (ℓ) − . This situation is realized very close to the condition ω = ±ω m and for a bandwidth of the order of the mechanical dissipation rate γ which, in typical optomecanical system, can be very small as shown in the insets of Fig. 5. On the other hand, the discrepancy between S (III) and ν (III) can be considerably larger, covering, for example, the full spectrum in plots (a) and (b). This is due to the fact that, in this case the difference between the corresponding n (III) + and n (III) − is proportional to the difference between κ 1 and κ 2 , which is relatively large in plots (a) and (b). In plot (c), on the contrary, the two mirrors are symmetric and correspondingly S (III) and ν (III) are very close, reproducing the results for S (I) and S (I) min .

VI. CONCLUSIONS
In conclusion, we have studied the relations between entanglement and squeezing in the continuous output field of a quantum optical system. In the stationary regime, the output field is conveniently analyzed in frequency. And correspondingly, the squeezing properties of the field can be investigated with homodyne or heterodyne techniques. In particular the long time integration of the homodyne or heterodyne signal, provides informations about specific spectral component of the field and of the corresponding squeezing. We have analyzed the discrete bosonic operators describing such spectral modes and we have studied the corresponding entanglement properties, thereby demonstrating that the measurable squeezing spectrum resulting from the spectral homodyne or heterodyne analysis of the field is, indeed a direct measurement of the corresponding logarithmic negativity.
When applied to an optomechanical system comprising a two-sided Fabriy-Perot cavity with a membrane in the middle, these findings help in identifying the specific spectral components of the output fields that are maximally entangled, showing, in particular, that maximum squeezing and entanglement is found between specific modes constituted of the superposition of carefully selected spectral components of the two outputs. We have also described how to probe such specific modes with homodyne and heterodyne techniques.
In general, a continuous-wave squeezed field combines, in a single spatial mode, a large number of spectral entangled sideband modes. It is, therefore, logical to ask, whether and how one could exploit such rich entanglement structure for real quantum-enhanced applications. Such question has been addressed, for example in Ref. [37][38][39][40], where it is discussed how to spatially separate the spectral sidebands of a continuous squeezed field in order to create N spatially independent entangled pairs and, hence, to prepare N quantum communication channels, whose actual number is limited only by the spectral resolution of the experimental apparatus and by the bandwidth of the squeezed signal. On a similar perspective, it is intriguing to ask if such large amount of entangled pairs, could be exploited as a resource for frequency encoded multimode entangled networks, in the stationary continuous-wave regime, alternative to that realized with pulsed frequency combs [18].
We finally remark that although, here we have focused on spectral components, the approach based on the filtered modes that we have described in Sec. II A is sufficiently general to be applicable to a wider area of experimental situations in which finite-time filtered mode are relevant [61,[71][72][73][74].
where ξ j are real and positive. A sufficient condition for entanglement can be defined in terms of the quantity Specifically, when E S < 2, for some values of ξ j , and φ j , then the two modes are entangled. In the case of Gaussian fields this criterion becomes also a necessary condition for entanglement (for appropriate values of ξ j ) [2]. In the analysis of the entanglement properties of Gaussian systems, for which all the informations are contained in the first and second moments of the field operators, it is useful to introduce the following matrix notation. We consider the column vector of operators b = (b 1 , b 2 , b † 1 , b † 2 ) T and the corresponding correlation matrix which is given by whose elements are {A} j,k = {b} j {b} k . The corresponding covariance matrix, C, namely the symmetric matrix of correlations of the quadrature x (0) j and x (π/2) j , can be used to compute entanglement measures, such as the logarithmic negativity. It is given by C = T A+A T 2 T T , where we have introduced the matrix The logarithmic negativity for Gaussian states is then computed as E N = max 0, − log 2 (ν) , where ν is the smallest symplectic eigenvalue of the covariance matrix of the partially transposed state, that can be expressed as C ′ = Π C Π, where Π is a diagonal matrix whose diagonal elements are (1, −1, 1, 1) [62]. In particular, this relation implies that the state is entangled when ν < 1.
A generic covariance matrix can always be transformed, using only local symplectic transformation, into the standard form [29] where a, b, c and c ′ are reals. In this case the corresponding matrix of correlations for the field operators reads where m ± = (c ± c ′ )/4, n 1 = (a − 1)/2, n 2 = (b − 1)/2.
The symplectic eigenvalue ν in the definition of the logarithmic negativity can be expressed in terms of the elements of these matrices as Correspondingly, these matrices can be used to determine an explicit expression for E S E S = 2       1 + 2 n 1 ξ 2 1 + 2 n 2 ξ 2 2 + 4 ξ 1 ξ 2 m − cos (φ 1 + φ 2 ) ξ 2 1 + ξ 2 2       (A6) that is minimized for and the corresponding minimum is min φ j ,ξ j E S = 2 1 + n 1 + n 2 − 4 m 2 − + (n 1 − n 2 ) 2 . (A8) We note that if m + = 0 (i.e. c = −c ′ ), then 2ν is equal to Eq. (A8), that is This result is important because joins directly an entanglement measure, namely the logarithmic negativity, to the field observables, namely the variances of the field quadratures.
As we have seen this is true only for the specific class of states for which m + = 0, that correspond to the condition b j b † k = 0 with j k = 1, 2 (or equivalently . A related result has been previously discussed in Ref. [77] where the symplectic eigenvalue has been shown to be equal to the EPR correlations in the case of symmetric states, for which a = b. As discussed in the main text, the condition m + = 0 is relevant for the study of entanglement between the spectral components of stationary continuous fields. The general, corresponding correlation matrix takes the form of Eq. (16), for which the correlations of the form b j b † k are zero. We note that in this case Eq. (A6) is equal to twice Eq. (22) when m − = |m| and φ 1 + φ 2 = θ 1 + θ 2 + arg[m], where arg(m) is the phase of the complex parameter m, that is introduced in Eq. (16). The minimization of Eq. (22) is therefore similar to Eqs. (A7) and (A8) (see Eq. (24) in the main text). The corresponding covariance matrix is given by We note that this matrix and the matrix in Eq. (16) are not, in general, in the standard form described by Eqs. (A3) and (A4). They can be cast in standard form by means, for example, of the single mode rotation that perform the transformation a(Ω) → e −i arg(m) a(Ω). Thereby the resulting matrices are equal to Eqs. (A3) and (A4) with c ′ = −c and m + = 0. Therefore the corresponding minimum symplectic eigenvalue of the partially transposed covariance matrix has the form of Eq. (A5) with m − = |m| and m + = 0, and it is explicitly given by Eqs. (27) and (24).
of Eq. (B10) is physically motivated by the fact that it results in a real filtered photocurrent that corresponds to an hermitian quantum operator, and hence it makes transparent the relation between the spectral properties of the detected photocurrent and the corresponding quantum observables of the stationary field. In particular the filtered photocurrent can be described by the hermitian operator ds cos (ǫ s + ϕ) I (θ) (s) where the normalization factor N τ is appropriately chosen in order to satisfy the commutation relation for quadrature operators J (θ) ∆,τ (ǫ, t), J (θ+ π 2 ) ∆,τ (ǫ, t) = 2 i, namely N τ = √ (1 + τ ǫ)/(2 + τ ǫ). Thus, the filtered photocurrent can be expressed as the sum of two filtered quadrature operators for the two frequency bands of width 1/τ each, centred at the frequencies ∆ ± ǫ, with the annihilation operator, of a single band filtered mode, defined as in Eq. (4). In the limit of large τ, Eq. (B14) reduces to Eq. (30) and Eq. (42), when, respectively, ∆ = 0 and ∆ 0.
and the corresponding power spectrum matrix is where C in is the correlation matrix of the input noise operators defined as δ(ω + ω ′ ) C in = ã in (ω)ã in (ω) T , and it is given by where n T is the number of thermal excitations of the mechanical oscillator. The explicit result for P out (ω) is given in Eq. (46).