Frequency-resolved photon correlations in cavity optomechanics

Frequency-resolved photon correlations have proven to be a useful resource to unveil nonlinearities hidden in standard observables such as the spectrum or the standard (color-blind) photon correlations. In this manuscript, we analyze the frequency-resolved correlations of the photons being emitted from an optomechanical system where light is nonlinearly coupled to the quantized motion of a mechanical mode of a resonator, but where the quantum nonlinear response is typically hard to evidence. We present and unravel a rich landscape of frequency-resolved correlations, and discuss how the time-delayed correlations can reveal information about the dynamics of the system. We also study the dependence of correlations on relevant parameters such as the single-photon coupling strength, the filtering linewidth, or the thermal noise in the environment. This enriched understanding of the system can trigger new experiments to probe nonlinear phenomena in optomechanics, and provide insights into dynamics of generic nonlinear systems.


I. INTRODUCTION
As light emerges from an open system, it carries a lot of information about the system and its dynamics. It is up to our ingenuity to learn how to extract that information. For example, by counting the number of photons at a given frequency ω [1] using a photodetector with spectral resolution Γ, we can obtain the emission spectrum S Γ (ω), and extract information about the underlying level structure of the system. If instead we perform a Hanbury Brown-Twiss (HBT) experiment [2], splitting the emitted light into two beams and measuring their intensity correlation, we can measure its second-order coherence g (2) (τ ) [3,4] which informs us about the statistical nature of the emitted light, e.g., whether it is of a quantum or classical character. These two observables (S Γ (ω) and g (2) (τ )) are arguably the most fundamental ones to characterize open quantum optical setups. However, sometimes the information they carry is not sufficient to unravel the dynamics of complex quantum systems -most notably, when several processes lead to multiple emission lines with competing statistics.
An intriguing system that is known to exhibit very rich physics, but whose TPS has not yet been considered, is single-mode cavity optomechanics (OM) [26], in which the optical and motional degrees of freedom of a resonator are nonlinearly coupled (see schematic in Fig. 1(a)). Cavity OM systems are particularly interesting as a platform for studying frequency-resolved intensity correlations, since the typical OM emission spectrum includes several lines from competing processes involving the creation or annihilation of vibrational quanta -phonons ( Fig. 1(b)). Intensity correlations between such processes have been used for the heralded generation of single phonons [27][28][29][30][31][32]. However, the theoretical descriptions of these correlations are based on simplified models [33], and provide a limited picture of the complex landscape of frequency-resolved correlations. Besides, these experiments have been performed in the linear regime, where the nonlinearity of the coupling is removed by strongly driving the cavity. Thus, finding signatures of the nonlinear OM couplings [34][35][36] is still an open challenge of the field, that could open the path to many OM quantum applications.
In this manuscript, we present the first complete analysis of the emission statistics from a generic OM system by studying its frequency-resolved photon correlations g (2) Γ (ω 1 , ω 2 ; τ ), and identifying its features with the arXiv:2009.06216v2 [quant-ph] 2 Oct 2020 system's underlying processes: (i) the effective Kerr cavity nonlinearity induced by the photon-phonon interaction, and (ii) a family of higher-order terms defining multi-phonon transitions, linear in the optical degree of freedom. We study the evolution of these features in terms of both the parameters of the OM system, e.g., the optomechanical single-photon coupling g 0 /κ, as well as the characteristics of the external measurement setup, e.g., the frequency filter linewidth Γ. We also calculate the temporal dynamics of frequency-resolved correlations g (2) Γ (ω 1 , ω 2 ; τ ), demonstrate how they encode information about the nature of emission processes, and discuss the relationship between the spectral and temporal resolution of the measurement setup. Finally, we also show how some of these frequency regions can be associated with the emission of non-classical light by studying the violation of the Cauchy-Schwarz inequality [17,37].
The manuscript is structured as follows. In Section II, we introduce the theoretical foundations of the paper, introducing the single-mode cavity OM Hamiltonian in Section II A, and defining the spectra S Γ (ω) and frequency-resolved correlations g (2) Γ (ω 1 , ω 2 ; τ ) in Section II B. In Section III we study the TPS of two elementary nonlinear Hamiltonians, namely, a coherently driven Kerr cavity and cavity-multi-phonon interactions, which are instrumental for understanding the frequency-resolved correlations of the cavity OM Hamiltonian characterized in Section IV. Finally, in Section V we demonstrate the violation of Cauchy-Schwarz inequalities in some frequency regions of the TPS, and summarize our findings in Section VI.

II. THEORETICAL FRAMEWORK: SETUP AND OBSERVABLES
A. Single-mode cavity optomechanics Single-mode cavity optomechanics studies the interaction of a single quantized mode of an optical cavity, with frequency ω a , with a quantized vibrational, or phonon, mode of frequency ω b , as schematically depicted in Fig. 1(a). The photon-phonon coupling can be implemented in various physical systems using radiation pressure [26], the photoelastic effect [38], or Raman scattering in molecular systems [39][40][41]. Irrespective of the physical mechanism inducing such interaction, the OM Hamiltonian can be written as (using = 1 throughout the manuscript): where a † (a)/b † (b) are the bosonic creation (annihilation) operators of the photon/phonon mode, and g 0 denotes the single-photon coupling parameter. The system is typically coherently driven with a laser exciting the cavity field, described by the following Hamiltonian: FIG. 1. (a) Schematic of a generic optomechanical (OM) system in which an optical cavity (with photon creation/annihilation operators a † /a) is dispersively coupled to a mechanical mode (with phonon creation/annihilation operators b † /b). The spectrum of light emitted from the cavity mode a, SΓ[a](ω1), is measured by frequency-blind detection of the photons passing through a spectral filter tuned to frequency ω1 with resolution Γ. Directly extending this detection setup, we can measure the frequency-resolved correlations g Γ [a](ω1, ω2) by splitting the light emitted from the cavity, filtering each of the beams separately, and measuring intensity correlations of the photocurrents from the frequency-blind detectors. (b) Illustrative one-photon spectrum SΓ[a](ω) of an OM system illuminated by a laser with frequency ωL, and operating in the single-photon strong coupling regime. It includes a dominant elastic scattering peak at ω = ωL, broadened by the filter resolution Γ, and lower-and higher-frequency peaks corresponding to the phonon creation (Stokes) and annihilation (anti-Stokes) lines shifted from ωL by multiples of the phonon frequency ω b . (c) Schematic two photon spectrum (TPS; Eq. (11)) of the cavity mode a of an OM system, with red and blue regions denoting frequency bunched g (2) Γ [a](ω1, ω2) > 1 and antibunched g (2) Γ [a](ω1, ω2) < 1 emission regions, respectively. Throughout this work, the spectra and color maps in TPSs are given in logarithmic scale.
where Ω is the driving amplitude and ω L the laser frequency. In a frame rotating with ω L the total Hamiltonian, i.e., H(t) = H OM + H L (t), becomes timeindependent: where ∆ a = ω a − ω L is the detuning between the cavity and laser frequencies.
Importantly, neither the optical cavity nor the mechanical mode are isolated from the photonic and phononic environments, inducing dissipation into them at rates κ and γ, respectively. To formally account for such losses, we describe the state of the OM system using a density matrix ρ. Assuming that the environmental timescales are much faster than the system ones (Born-Markov approximation), the dynamics of the system is then described by the following master equation [42]: Kossakowski-Sudarshan-Lindblad (GKSL) terms, and n th b is the thermal population of the phonon bath, which both increases the decay rate of phonons, and governs the rate of incoherent pumping of mechanical vibrations by the thermal environment. Importantly, this thermal phonon bath population depends on the ratio between the phonon energy ω b and the thermal energy k B T , such that it will be very small for high-frequency optical phonon modes, e.g., in molecular systems [39,40,43,44]. For the sake of illustration we will assume n th b (T ) ≈ 0 throughout most of this manuscript to keep the discussion of the physics simpler, and only briefly consider the effect of thermal population on frequency-resolved correlations in Section IV. We will also ensure that the optical cavity is always weakly populated a † a 1. Throughout this manuscript, we will always consider the sideband-resolved regime, where ω b > κ, now commonplace in many cavity OM systems [26]. The naming of this limit refers to the ability of the cavity to resolve the main relevant scattering processes identified in the system (see Fig. 1(b)): elastic/Rayleigh scattering corresponding to emission processes that involve no exchange of excitations with the mechanical mode, so that the frequency of scattered photons matches that of the incident laser ω R = ω L ; Stokes and anti-Stokes processes in which the cavity photon either loses or absorbs the energy of one phonon (or an integer number n of phonons) so that light is emitted at ω S/aS = ω L ∓ nω b . Furthermore, we will consider systems operating from the more accessible weak coupling limit (g 0 κ) to the more demanding strong single-photon coupling regime (g 0 κ) approached in physical setups based on cold atoms [45], molecular optomechanics [39,40,43,44], microwave micromechanics [46], and phoxonic cystals [47] [48]. To gain insight into the underlying dynamics of the OM system, and the statistics of the emitted light, it is instructive to perform the polaron transformation U = exp (g 0 /ω b )a † a(b † − b) [34,35] on the OM Hamiltonian. This transformation decouples the photon and phonon modes in H OM , at the expense of transforming the harmonic cavity Hamiltonian into an anharmonic one: with Kerr nonlinearity parameter ∆ g = g 2 0 /ω b . The eigenstates of the transformed Hamiltonian are then defined by the photon/phonon number states |n a , n b and have the following energies: In the regime where the nonlinearity ∆ g becomes comparable to the cavity losses κ, the system cannot readily absorb two photons with the same frequency, resulting in the well-known OM photon blockade effect [34,35].
In the transformed picture with non-vanishing pumping Ω = 0, the explicit interaction between photons and phonons appears in the transformed coherent pumping term which, in the frame rotating with ω L , takes on the form: as well as in the transformed GKSL terms, which we will for now omit. The transformed HamiltonianH L can be broken into several contributions by expanding the exponential: where f n = (g 0 /ω b ) n /n!. When g 0 /ω b 1, the first term of the expansion -the standard coherent driving term -dominates.
We therefore see that the transformed driven OM Hamiltonian,H = U † HU includes two terms which can induce nonlinear dynamics: (i) the anharmonicity of the optical mode described by HamiltonianH OM , and governed by ∆ g = g 2 0 /ω b , and (ii) a series of cavityphonon coupling terms inH L , determined by parameters f n ∝ (g 0 /ω b ) n , and describing n-phonon mediated processes. Most of the literature discussing optomechanics in the single-photon strong coupling regime is focused on exploring the photon blockade effect induced by the Kerr nonlinearity ∆ g > κ, and identified in frequency-blind correlations in the limit of Ω → 0 [34,35]. Here instead, we study the frequency-resolved photon correlations and explicitly consider both sources of nonlinearities. In fact, for pedagogical purposes, we will analyze in Section III separately the TPS of these two contributions, namely, of a coherently driven Kerr-cavity Hamiltonian (in Subsection III A), and of a higher-order multiphonon terms of Eq. (8) (in Subsection III B), which will help us to understand the TPS features of the complete cavity-driven OM Hamiltonian studied in Section IV.

B. Optical observables: spectrum and correlations
As mentioned in the introduction, light emitted by the cavity carries information that can be used to characterize its underlying dynamics. The simplest measurement that one can perform is to count the number of photons emitted around a given frequency ω, within a frequency window Γ determined by the resolution of the photodetector, or the spectral width of the filter set up before the color-blind detector (see schematic in Fig. 1(a)). We label this magnitude as the one-photon spectrum, and calculate it as follows [1]: (9) Throughout this manuscript, we always use bracket notation, i.e., [a], to denote the field operator that is being measured by the detector, e.g., here the cavity mode a. In the above equation, lim t→∞ indicates that the dynamics of the one-photon correlator a † (t+τ )a(t) should be calculated in the steady state. This definition also assumes a Lorentzian filter profile with linewidth Γ, which naturally broadens the emission lines. For example, elastic scattering from the cavity will no longer appear as a Dirac delta δ(ω L ), but rather as a Lorentzian with linewidth Γ centered at ω L which, as we will see, sometimes masks other features in the OM spectrum. Theoretically, this elastic contribution can be removed by observing that the optical mode a can be represented as the sum of its mean value and fluctuations a = a + δa, and calculating only the spectrum of the operator δa, which we will denote as S Γ [δa]. In an experiment, this removal can be achieved by self-homodyning the emitted light with the one of the driving laser [49,50].
Another widely used quantity in the characterization of quantum optical setups is the second-order coherence function, labelled as g (2) [a](τ ) [3,4] (denoting again in the bracket the field operator being measured), and defined as: Experimentally, g (2) [a](τ ) is measured with a HBT setup by dividing the light emitted from the cavity with a beam splitter, and then measuring intensity correlations between the photon detection in each of the beams [2]. This quantity allows one to distinguish between the classical and quantum nature of the emission. For example, the detection of g (2) [a](0) < 1 (subpoissonian) or g (2) [a](0) < g (2) [a](τ ) (antibunched) can both only be obtained with quantum light fields [3,4]. Importantly, g (2) [a](τ ) as defined in Eq. (10) accounts for all the different emission processes occurring in the system regardless of their frequencies. In complex quantum systems, however, where several emission processes with different frequencies and statistics simultaneously occur, this results in a loss of information which can be recovered by placing frequency filters in each of the paths of the HBT configuration (see Fig. 1(a)). This upgrade results in the measurement of the frequencyresolved two-photon correlations, which can be calculated as follows [5][6][7][8][9]: Here are the operators describing light passing through the Lorentzian frequency filters, and T and : : enforce the time-and normalordering of the cavity mode operators a. For τ = 0, this magnitude defines the two-photon spectrum (TPS) [14][15][16] Γ [a](ω 1 , ω 2 ), which carries information about the correlations of the photons emitted at frequencies ω 1 and ω 2 , given a filter linewidth Γ (see Fig. 1(c) for an example of a TPS from an OM system). Generalizing the notation inherited from the standard photon correlations, we will refer to the emission with g Γ [a](ω 1 , ω 2 ) < 1 as frequency bunched and antibunched, respectively [51]. As shown in other works [10,[14][15][16], the TPS is typically characterized by a grid of horizontal and vertical features, crossed by antidiagonal ones. The latter corresponds to filtering frequencies corresponding to two-photon processes in which the intermediate state is virtual, i.e., not an eigenstate of the system, dubbed in other works as leapfrog processes [14]. The former (vertical/horizontal structure) corresponds to fixing one of the filters at a transition frequency between the eigenfrequency of the system. In the limit of very large filter linewidths, as expected, we recover the standard colorblind correlation measurements, that is g Similarly to what it occurs for S Γ [a](ω), the frequencyresolved photon correlations near the elastic scattering frequencies might also be dominated by those of the laser light g Thus, in order to unveil the intrinsic dynamics of the OM interaction Hamiltonian, we will -when explicitly noted -consider both the TPS of the mode a, and of the fluctuations δa. The latter will be denoted as g (2) Γ [δa](ω 1 , ω 2 ), and calculated from Eq. (11) by replacing the operator a with δa. In an experiment, this measurement could be performed though the extension of the self-homodyning setup described in Ref. [49], where the light emitted from the cavity is mixed with the driving laser before splitting it in the HBT setup. Another option to remove the elastic spectral components is to use notch filters. However, to model the effect of these filters properly, one would need to extend the formalism used here, suited for Lorentzian filters, e.g., by adopting the approach developed by Kamide et al. [52], which lies beyond the scope of this work.
The numerical framework for calculating the TPS is based on the contributions from del Valle et al. [10] and Holdaway et al. [13], and described in more detail in Appendix A.

III. CORRELATIONS OF UNDERLYING NONLINEAR PROCESSES
In Section II we have shown how, using the polarontransformed picture, the OM Hamiltonian can be mapped to that describing a pair of decoupled harmonic (b) and anharmonic (a) oscillators (see Eq. (5)), whereas the cavity driving Hamiltonian can be expanded into a series of terms which include the standard coherent drive, but also higher-order nonlinear interaction terms (see Eq. (8)). Since all of these processes contribute to the frequency-resolved correlations of the full OM Hamiltonian, in this section we consider the frequency correlations induced by each of these individual nonlinear processes separately. This will help us to understand the complete picture when we analyze it in Section IV.
Before we continue, we should note that throughout this section we discuss the spectra and correlations of mode a for the respective Hamiltonians expressed through this operator, as if these Hamiltonians were given in an untransformed picture. Similarly, we will consider the GSKL term in the original form, given by the second term on the right-hand side of Eq. (4).

A. Coherently driven Kerr system
The first ingredient of the transformed optomechanical Hamiltonian that gives rise to non-trivial correlations is the phonon-mediated Kerr cavity interaction. Its Hamiltonian reads: which, to account for losses, has to be complemented with the aforementioned GKSL term κ/2L a [ρ]. Furthermore, we assume that the cavity is driven by a coherent laser described by the Hamiltonian H (0) = −iΩ(a † − a), and define the complete Hamiltonian H Kerr,Ω = H Kerr +H (0) . We note that the TPS of Kerr cavities was considered before in Ref. [53], but in a very different scenarioin the regime of very large Kerr nonlinearity, where one must redefine the observables to obtain physical results, and under incoherent cavity pumping. Before considering directly the TPS of the coherently driven Kerr Hamiltonian, H Kerr,Ω , we will first develop some intuition of its emergent features by expanding the Kerr nonlinearity around the fluctuations δa of mode a using a → δa + α and α = a [54], arriving at: The terms in the first line of this expansion represent energy shifts and an additional coherent driving of δa. The second line describes one-mode squeezing, or degenerate parametric amplification, corresponding to the simultaneous creation or annihilation of two cavity photons by the incident laser, which results in a strong correlation of field components oscillating at frequencies ω 1 and ω 2 with ω 1 + ω 2 = 2ω L . The third line describes cubic processes, and the fourth the Kerr nonlinearity of the fluctuations δa.
In Fig. 2(b), we plot the spectra S Γ [δa] (upper panel) and TPS g (2) Γ [δa] (lower panel) of the fluctuations δa including only the terms describing the squeezing of δa for a system exhibiting a strong Kerr nonlinearity ∆ g /κ = 0.5 and fixing the laser detuning to ∆ a = ∆ g [55]. While the one-photon spectrum shows a single Lorentzian peak, the TPS exhibits a clear leapfrog bunching behavior along the antidiagonal ω 1 + ω 2 = 2ω L , as we predicted to occur due to the degenerate parametric amplification. We also identify indistinguishability bunching along the diagonal ω 1 = ω 2 [14,15], which emerges because the photons emitted within the finite-time response time of the filters (Γ −1 ), appear as if they arrived simultaneously at the detector, resulting in an increase of the frequency-resolved bunching along this line.
In Fig. 2(c) we plot the spectra and TPS corresponding to a nonlinear Kerr interaction between the δa operators, including also the driving term ∝ α 3 (δa † + δa) for the same parameters as panel (b). In this case, we find that besides the leapfrog and indistinguishability bunching features, a blue region of weak frequency antibunching emerges around the central region (see the zoom in the inset panel). Its origin can be attributed to the nonlinear energy shift of the excited states due to the Kerr nonlinearity, which prevents the driving term from populating efficiently the state |2 (see Fig. 2(a)). The system  then becomes reminiscent of a two-level system, whose TPS was discussed in Ref. [14], exhibiting particularly strong antibunching near the resonance ω i ∼ ω L along the antidiagonal ω 1 + ω 2 = 2ω L (see inset in Fig, 2(c)). In the case of frequency-blind correlations, this is equivalent to the photon blockade effect.
Finally, we plot in Fig. 2(d) the spectra and the TPS for the coherently driven Kerr Hamiltonian H Kerr,Ω , for both the cavity mode a, and its fluctuations δa (plotting each TPS separately by cutting the correlation map along the diagonal defined by ω 1 = ω 2 ). Both the spectra and TPS of mode a are nearly identical to those exhibited by the fluctuations δa in Fig. 2(c). The main difference between the spectra of a and δa in Fig. 2(d) is that the latter does not include the filter-broadened elastic scattering peak. Further, the TPS of δa for the H Kerr,Ω Hamiltonian does not show the blue dip on the antidiagonal (see inset in Fig, 2(d)). This indicates that the admixing of the coherent field a to the light emerging from a Kerr system can switch the statistical properties of the leapfrog-dominated emission from strong frequency bunching to antibunching. A similar effect was recently analyzed by Casalengua et al. in the context of frequency-blind correlations [56].

B. Multi-phonon processes
The second nonlinear ingredient of the polarontransformed OM Hamiltonian is described by the series of terms, written in Eq. (8), linear in the optical degree of freedom a, and increasingly nonlinear in the mechanical mode operator b. In Fig. 3 we analyze how they contribute to the emergence of additional features in the spectra and TPS of the system: a. 1st termH Since this Hamiltonian does not include a laser driving term, we find that a = 0. Thus, the one-photon spectra and TPS of the a mode will be the same as the one of δa, which is what is plotted in Fig. 3(b). In that figure, we observe how the spectra S Γ [a](ω) develops in this case three peaks corresponding to the elastic (ω ≈ ω L ) and Stokes/anti-Stokes emission processes (ω ≈ ω L ± ω b ), denoted as processes (iii) and (iv) in Fig. 3(a). Regarding the corresponding TPS, its most prominent feature corresponds to the strong bunching antidiagonal around ω 1 + ω 2 ≈ 2ω L . This feature can be understood by separating the passive (beamsplitter-like) ab † +a † b and active (two-mode-squeezing) ab + a † b † interaction terms. In a perturbation picture, the latter will drive the system into states |i a , i b , generating strong a ↔ b intensity correlations, and the former will transfer these correlations into intensity autocorrelations of mode a (simultaneously, b),   populating the state |2 a , 0 b . This is effectively the same mechanism (quadratic in both g 0 and Ω) as the degenerate squeezing identified in the expansion of the Kerr interaction, and thus results in the strong antidiagonal bunching line in the TPS (we discuss these features in more detail in Appendix B).
On top of this strong bunching antidiagonal, the TPS develops a vertical/horizontal grid of correlations ≈ 1 (white regions) when fixing one of the filters to the Stokes/anti-Stokes frequencies, marked as (iii) and (iv) in Fig. 3(a-b). These correlations correspond to emission processes between the real energy levels of the system and will be discussed in more detail in the subsequent section.
Finally, let us note that the HamiltonianH L is identical to the linearized OM Hamiltonian, and describes the dynamics of OM systems operating in the single-photon weak-coupling limit (g 0 κ), which is where the majority of OM systems currently work.
b. 2nd termH In Fig. 3(c) we plot the one-photon spectra and TPS of this Hamiltonian for both the cavity mode a and its fluctuations δa. The one-photon spectra are both dominated by two peaks: one at the laser frequency ω ≈ ω L , and another one displaced twice the phonon frequency ω ≈ ω L − 2ω b . This translates into a TPS developing two strong bunching antidiagonals at both ω 1 + ω 2 ≈ 2ω L , 2ω L − 2ω b . Again, these features can be intuitively understood by expandingH (2) L into two terms: (i) In the absence of the b-squeezing terms (sq) or additional driving terms, the b mode population vanishes b † b = 0, and the dispersive term dis turns into a coherent drive for mode a, decoupled from b, generating TPS features ∼ 1 at ω 1/2 = 0. Alternatively, in the absence of the dis terms, this Hamiltonian becomes similar tõ H L with two-phonon operators replacing single-phonon operators. This explains the onset of the strong antidiagonal ω 1 + ω 2 = 2ω L and indistinguishability bunching lines along the diagonal of the TPS of HamiltonianH  Fig. 3(c). (ii) The two-photon Stokes emission peak in the one-photon spectra and the shifted antidiagonal ω 1 + ω 2 = 2ω L − 2ω b can be explained by tracking the Hilbert space accessed by the HamiltonianH (2) L when both the sq and dis terms are present: |i a , i b such that i b is even. This space includes both states |2 a , 0 b and |0 a , 2 b , that can be connected by a two-photon processes with energy given by ω 1 + ω 2 = 2ω L − 2ω b , which defines the off-centre diagonal observed in Fig. 3(c). Furthermore, we note that the non-vanishing coherent amplitude contribution α to the a operator generates a strong elastic, filter-broadened emission line in the spectrum S Γ [a], which also homogeneously lowers the values of the observed TPS g (2) Γ (ω 1 , ω 2 ), by flooding the detection with coherent light (with values ∼ 1). This effect can be seen by comparing the TPS correlations obtained for a and for the δa fluctuations (regions below and above the diagonal in Fig. 3(c)) c. Full HamiltonianH L . Finally, we show in Fig. 3(d) the spectra and TPS of the entire series of linear drive, single-and multi-phonon processesH L given in Eq. (8) (see caption of Fig. 3 for parameters), for both operator a and its fluctuations δa. The spectrum is made up of multiple emission peaks at ω i = ω L ±nω b , which are also visible in the TPS as vertical and horizontal lines. We also identify a series of antidiagonal bunching features ω 1 + ω 2 = 2ω L + nω b resulting from the multiphonon leapfrog transitions, including one (n = −1) depicted in panels (a) and (d) as (v), which describes emission into the excited vibrational state |0 a , 1 b . Interestingly, the TPS of the cavity mode g (2) Γ [a](ω 1 , ω 2 ) and its fluctuations g Fig. 3(d) appear to be qualitatively very different. In particular, one observes that the absence of the elastic scattering contribution seemingly removes all the traces of frequency antibunched (blue) emission. This is particularly striking along the central antidiagonal line, where the strong frequency antibunching in g (2) Γ [a](ω 1 , ω 2 ) near the elastic emission peak is entirely replaced by a bunched characteristic of fluctuations g (2) Γ [δa](ω 1 , ω 2 ), similarly to what occurs in the TPS of the coherently driven Kerr cavity that we show in Fig. 2(d). We note that this strong contrast, and more generally, the presence of large white and blue areas in the TPS g Fig. 3(d), is largely due to the linear drive included inH L , but absent inH (1,2) L . In the presence of coherent driving, in Fig. 3(b,c) we would also find mostly white g (2) Γ [a], with some significantly frequency antibunched (blue) areas.

IV. FREQUENCY-RESOLVED CORRELATIONS IN SINGLE-CAVITY OPTOMECHANICS CAVITIES
In the previous section we studied the features of the TPS of the two terms that form the H OM in the polaron picture. With this knowledge as a basis, in this section we finally consider the emission from single-cavity OM systems as described in Eqs. (3) and (4), and analyze their TPS in detail. We relate to the knowledge developed in the previous section when possible, and expand it when new features emerge in the correlation maps of the OM setup.
In particular, we study the dependence of the TPS on relevant parameters of the system such as optomechanical coupling and thermal phonon population, in subsection IV A. Then, in subsection IV B we consider the dynamics of the correlations with time-delay τ of some of the features of the TPS: the Stokes-anti-Stokes, and leapfrog correlations.
For the weakest coupling g 0 /κ = 0.1, the one-photon spectrum of mode a, shown in solid gray line in Fig. 4(a), is almost entirely dominated by the elastic peak, as it also occurs for the corresponding color-blind photon correlations, which display a value very close to 1, i.e., g (2) (0) = 0.99. Interestingly, the corresponding TPS already reveals a non-trivial correlation structure with some regions of frequency bunching and antibunching. For example, the TPS features a vertical and horizontal grid of uncorrelated transitions characterized by g (2) Γ [a](ω 1 , ω 2 ) ≈ 1 (marked by a grid of white lines, and including Stokes and anti-Stokes processes denoted schematically as (iii) and (iv) in Fig. 3(a)). Conversely, the spectrum and TPS of the fluctuations are similar to those shown and discussed in Fig. 3(b) corresponding to the Hamiltoniañ H (1) L : in the spectrum (dashed black line in Fig. 4(a)) the elastic scattering peak becomes suppressed, which favours the observation of the phonon sideband peaks; in the TPS, the weak frequency-antibunched regions disappear favoring the appearance of frequency-bunched regions. The most prominent features are still the horizontal/vertical grid of uncorrelated (white) Stokes/anti-Stokes transitions, and the bunched antidiagonal line at ω 1 + ω 2 = 2ω L . This qualitative resemblance with the features ofH (1) L corroborates our expectation that in the regime of weak single-photon coupling g 0 /κ 1 analyzed in Fig. 4(a), the TPS of an OM system will be nearly identical to that described by linearizing the OM Hamiltonian.
As we increase the coupling g 0 (Figs. 4(b,c)) towards the single-photon strong coupling regime, where the linearized Hamiltonian fails to describe the response of the optomechanical system [34,35], we begin to observe features originating from the nonlinear form of the coupling. In the spectra S Γ , we find additional Stokes and anti-Stokes peaks, corresponding to the multi-phonon processes. In the TPS, we also identify antidiagonal bunch- Γ [δa](ω1, ω2) of a cavity OM system (Eqs. (3,4)) in the sideband-resolved regime (ω b /κ = 2), for three values of g0/κ = (0.1, 0.5, 1). For these parameters, the photon nonlinear shift induced by the phonons is ∆g/κ = (0.005, 0.125, 0.5), and the values of the frequency-blind twophoton correlations g (2) given above the spectra indicate the onset of a strong photon blockade. (d) Guidelines of different processes (i-v) appearing in the TPS of OM, shown schematically in panels in Fig. 2(a) and in Fig. 3(a). We also mark as A, B, C and D the specific pairs of filter wavelengths that are discussed in the main text. (e-g) Magnified regions of the TPS maps near the Stokes-anti-Stokes emission correlation ((e,g), ω 1/2 = ωL ± ω b ), and near the elastic scattering peaks (f), for both the fluctuations δa and dressed cavity field a. For all the systems we set the filter linewidth as Γ/κ = 0.05, decay rate of phonons as γ/κ = 0.1, and consider a laser with amplitude Ω/κ = 0.1 tuned to the first excited eigenstate ωL = ωa − ∆g.
ing lines (at ω 1 + ω 2 = 2ω L − nω b for n = 1, 2, 3, ...) originating from leapfrog emission driving the system to the excited vibrational states (including processes (v) marked in Fig. 3(a) for n = 1), and a much weaker antidiagonal line at ω 1 +ω 2 = 2ω L +ω b , describing emission induced by the laser driving the system from the excited vibrational states. For the largest coupling (panel (c)), we find -as in Fig. 3 -a weak and narrow antibunching along the central antidiagonal (see panel (f)), resulting from the interference of elastic scattering peak and leapfrog emission. We also find that some of the previously highlighted features become stronger, e.g., the antidiagonal bunching line at ω 1 + ω 2 = 2ω L , as well as the frequency-antibunched (blue) regions for operator a.
Finally, the TPS for both the cavity operator a and the fluctuations δa include a grid of vertical and horizontal antibunching features coming from single-photon transitions between the relevant energy levels of the system, including ω i = ω L ± ω b (processes (iii, iv) in Fig. 3(a)). Let us now highlight how the TPS look like around some relevant crossings of this grid: • The Stokes-Stokes correlation (point C in Fig. 4(d)) features a structure, in both a and δa, similar to one identified in Kerr cavities (see Fig. 2(d)) or two-level systems [14], with vertical (ω 1 = ω S ) and horizontal (ω 2 = ω S ) lines depicting strong antibunching, and crossing with the bunching features associated with the leapfrog ω 1 + ω 2 = 2ω L − 2ω b and indistinguishability. This two-level-like TPS results from the effective Kerr detuning between the states connected via Stokes transitions: |0 a , 0 b , |1 a , 1 b and |2 a , 2 b .
• The anti-Stokes-anti-Stokes (point D) crossing is mostly dominated by the elastic response since the probability for these anti-Stokes processes to occur is very low with these system parameters.
significantly suppressed (note the logarithmic color scale) when the filters are tuned to the Stokes and anti-Stokes lines [33]. This change of behavior occurs because different photon-emission processes start to dominate the correlations (strongly affecting both the numerator and denominator in Eq. 10) : a) Stokes-anti-Stokes transitions mediated by a real, one-phonon state |0 a , 1 b , described by sequential transitions (iii) and (iv) depicted in Fig. 3(a). As described in Ref. [39], in the absence of thermal phonon population, and for sufficiently weak coherent pumping, an anti-Stokes photon is necessarily accompanied by a previous Stokes process, yielding strong bunching. b) Transitions from a two-photon state, similar to the leapfrog process (i), but where the intermediate state is a proper energy level of the system, c) Stokes and anti-Stokes transitions mediated by the thermal phonons (when present). The contribution from these three mechanisms and, when present, interference with the elastic scattering, can give rise to a nontrivial dependence of Stokes-anti-Stokes correlations on the coupling parameter g 0 and the time delay τ , which we discuss in detail in the following subsections.

B. Stokes-anti-Stokes and leapfrog correlationsdependence on parameters
In this section, we will explore in more detail the dependence on the parameters of the OM system of some of the features of the TPS discussed in the previous section. In particular, we choose two points of the TPS of Fig. 4, namely, a point that corresponds to a leapfrog correlation line (B in Fig. 4(d)) and the Stokes-anti-Stokes correlations (point A in Fig. 4(d)). We plot in green/orange lines, respectively, in Fig. 5(a) the evolution of their correlations as a function of the granularity parameter g 0 /κ (the importance of the width of the filter is discussed in Appendix C) . We also plot together with them the evolution of the colorblind photon correlation g (2) Γ [a](0) (solid blue). The correlations g (2) Γ [a](ω 1 , ω 2 ) in this panel are obtained for the cavity operator a. We observe that the leapfrog correlation grows from the elastic-field induced g (2) Γ [a](ω 1 , ω 2 ; τ = 0) ≈ 1, to become very strongly bunched for g 0 /κ ∼ 0.7, and then relaxing back to 1 for g 0 /κ ∼ 1. On the contrary, the Stokes-anti-Stokes correlations exhibit a strong frequency antibunching for coupling parameters g 0 /κ ∼ 0.05, and then grows to develop a strongly bunched signal until it relaxes back to one, similarly to the leapfrog point correlations. In both cases, it is important to note the higher sensitivity to g 0 /κ of the frequency-resolved correlations compared to its colorblind counterpart, which only starts deviating from the elastically dominated correlations to show the expected OM blockade [57,58] for g 0 /κ 0.3. Note that all these observations still hold if we consider a nonvanishing incoherent thermal pumping of the mechanical degree of freedom (n th b = 0.1, depicted by dashed lines of the same color).
A particularly striking feature of the frequencyresolved correlations depicted in Fig. 5(a) is the strong frequency antibunching displayed by the Stokes-anti Stokes correlations for small g 0 /κ parameters. To learn more about its behaviour, as we did in previous Sections, we will compare the results of Fig. 5(a) with the correlations of the fluctuations δa for the same parameters, which is shown in Fig. 5(b). There, we find no trace of frequency antibunching for any of those points, as the Stokes-anti-Stokes, leapfrog, and even the frequencyblind correlations exhibit bunching for the entire range of the coupling parameter. This calculation suggests that the frequency antibunching of the Stokes-anti-Stokes point could be related to the effective interference between the Stokes emission, and the elastically scattered laser light, in an effect which we dub as interference antibunching. We will confirm this intuition in the next subsection by studying in detail the time dynamics of such correlation points.

C. Dynamics of frequency-resolved correlations
In this section, we will consider the dynamics of the frequency-resolved correlations in the Stokes-anti-Stokes and leapfrog emission processes discussed in the previous section. Furthermore, we will also study in more detail the mechanism of interference antibunching previously suggested as an explanation for the Stokes-anti-Stokes correlation behaviour, and briefly comment on the con- nection between the temporal and spectral resolution of our setup.

Stokes-anti-Stokes correlation dynamics
To illustrate the origin of the interference antibunching of the Stokes-anti-Stokes correlations found for g 0 /κ < 0.1, we plot in Fig. 6(a) the time-delayed correlations g (2) Γ [a](ω S , ω aS ; τ ) (solid orange line) for the parameters that maximize the antibunching found in Fig. 5(a): (g 0 , Γ)/κ = (0.08, 0.05). We observe in these time-delayed correlations that the antibunching g (2) Γ [a](ω S , ω aS ; τ = 0) < 1 results from the largeamplitude oscillations with frequency ω b , which are not present when considering correlations of the fluctuations g (2) Γ [δa](ω S , ω aS ; τ ) (dashed red line). This corroborates our previous statement that the interference antibunching results from the elastic component dominating over the anti-Stokes emission line, rather than originating from the Stokes-anti-Stokes emission mediated by a real one-phonon state |0 a , 1 b .
Notably, the strong bunching of the fluctuations in Fig. 6(a) is asymmetric near τ = 0, and decays as exp(−Γ|τ |) (with filter response time) for larger delays. Since this emission pathway is mediated by an excited phonon state, we would expect the correlations to decay over time with the characteristic rate of phonon decay (γ) instead. This discrepancy can be attributed to the choice of very narrow filters Γ/γ = 0.5, Γ/κ = 0.05, made to prioritize the spectral resolution of the TPS. This narrow filter linewidth imposes a poor temporal response of the setup, which effectively masks the intrinsic dynamics of the Stokes-anti-Stokes emission pathway. This issue can be solved by choosing a range of parameters such that Γ ∼ ω b , κ γ, which would reduce the spectral selectivity of the detection but could offer insights into the dynamics of such processes. We consider such an arrangement in Fig. 6(b), setting wide filters Γ/κ = 0.5, and simultaneously increasing the coupling g 0 /κ = 1 to amplify the anti-Stokes emission. We plot the correlations with (g (2) Γ [a], solid orange line), and without the elastic component (g (2) Γ [δa], dashed red line). As previously observed [59], g Γ [a](ω S , ω aS ; τ ) is now strongly asymmetric with respect to delay τ . This reflects the fact that in the regime of near-zero steady-state phonon populations, a phonon-annihilating anti-Stokes emission has to be preceded by a Stokes emission in which that phonon is generated. In other words, the two-photon process occurs with a given temporal order. Furthermore, the correlations decay approximately as exp(−γτ ) for τ > 0, which demonstrates that the increased time resolution, offered by a broad Γ γ filter, allows us access to the vibrational dynamics of the system. In Appendix C we provide additional calculations that illustrate the impact of the filtering linewidth on the frequency-resolved correlations of these systems.

Leapfrog dynamics
We can perform a similar analysis for the frequencyresolved correlations between the strongly bunched photons emitted in a leapfrog process. The temporal dynamics of this process (Fig. 6(c)) again shows large oscillations resulting from the interference with the elastic peak (solid green line), that are removed if we instead consider the correlations of the fluctuations δa (dashed line). The resulting correlations are clearly symmetric with respect to the delay time τ , highlighting the fact that there is no  (16)) in an optomechanical system in the (a) weak g0/κ = 0.1 and (b) strong g0/κ = 1 single-photon coupling regime. Besides the parameters given in the plots, for all the systems we set the linewidth of the filter as Γ = 0.05κ, the energy and decay rate of phonons as ω b /κ = 2 and γ/κ = 0.1, and consider laser with amplitude Ω/κ = 0.1 tuned to the first excited eigenstate ωL = ωa − ∆g. particular time order of the emission, and decay exponentially over time. In principle, since this two-photon transition is not mediated by the emission or absorption of phonons, the decay of correlations over time should be governed by the characteristics of the cavity -κ -and the temporal resolution of the detection setup determined by Γ. For the regime of parameter chosen, i.e., Γ/κ = 0.05, it is the temporal response of the filter which provides the decay time of the correlations. As with the Stokesanti-Stokes correlations, to resolve the intrinsic dynamics of the leapfrog correlations, governed by κ, one would require much broader filters Γ > κ, and a system where the leapfrog processes would be separated far enough from Stokes and elastic emission lines.
It is, however, possible to test the non-classical character of the correlations between different frequency channels ω 1 and ω 2 by using the information encoded in the TPS g Γ [a](ω 1 , ω 2 ) (or in the correlations of any other mode, for example δa). The idea consists in harnessing one of the features of classical correlations between two random variables, that is, the Cauchy-Schwarz inequality (CSI) [37]. Identifying the intensities at the two frequencies with such variables, we can write down the CSI as Γ [a](ω 2 , ω 2 ) . (15) Since the CSI should hold for any classically correlated variables, one can define a parameter (introduced in Ref. [19]) that indicates the degree of CSI violation: yielding a sufficient condition for the observation of nonclassical correlations that is R Γ [a](ω 1 , ω 2 ) > 1. Furthermore, we can also measure the degree of non-classicality of the inelastically scattered light, by calculating R Γ from the fluctuations δa, denoted as R Γ [δa](ω 1 , ω 2 ). In Figs. 7(a,b) we plot the frequency maps of R Γ [a] and R Γ [δa] for the same parameters as in Figs. 4(a,c), showing how the strongly bunched regions of the TPS often display a large CSI violation (in green), being therefore a source of non-classical correlations. The results in (a) correspond to g 0 /κ = 0.1 and reveal a strong violation of CSI over a broad region near the ω 1 + ω 2 = 2ω L antidiagonal. For g 0 /κ = 1 (panel (b)) the CSI is violated for narrower regions around multiple ω 1 + ω 2 = 2ω L + nω b antidiagonals. As in Fig. 4(c), these additional features originate form the intrinsic nonlinearity of the optomechanical Hamiltonian, and are therefore absent in the CSI maps of its linearized form (not shown here). We also find that the violation of CSI can be enhanced if we remove the coherent elastic field which dominates the emission in (a). One of the differences with the already observed CSI violation in OM systems [30] is that here one does not rely on a heralding preparation step.

VI. CONCLUSIONS & OUTLOOK
In summary, we present a systematic study of frequency-resolved correlations in cavity optomechanical systems in the sideband-resolved regime and for moderate to strong single-photon coupling strengths. We show how the two-photon correlation spectra unveil a rich landscape of correlations hidden in other observables. We also provide an intuitive picture that explains these correlations based on the anharmonic level structure and multi-mode squeezing Hamiltonians appearing through the nonlinear optomechanical coupling, and test their non-classical nature based on the Cauchy-Schwarz inequality violation. Importantly, non-trivial frequencyresolved correlations appear already for smaller coupling strengths than for the frequency blind correlations, thus opening new avenues to observe nonlinear phenomena in optomechanical systems operating far from the singlephoton strong coupling limit. Furthermore, we believe this work opens many research directions that one can follow. For example, although we focused on a particular range of parameters involving low optical quality factors, Q's, and high-frequency optical phonons, which best describes the novel implementation of optomechanics in molecular systems [39,40], there are many other relevant questions to be answered, e.g., what will be the role of phonon population in systems with low-frequency phonon modes where this population will be non-negligible?; how will the correlation maps change with incoherent pumping?; how will the balance of effects of Kerr and multi-phonon nonlinearities change in larger-Q, weakly pumped cavities? Another interesting direction could be to harness the knowledge acquired through these frequency-resolved correlations to connect it with recent findings in photonic Cooper pairs [60], to design non-classical photon [17,18] or, as recently proposed in [61], phonon sources. AGT acknowledges support from CSIC Research Platform on Quantum Technologies PTI-001 and from Spanish project PGC2018-094792-B-100 (MCIU/AEI/FEDER, EU). RE and JA acknowledge project PID2019-107432GB-I00 from the Spanish Ministry of Science and Innovation, project H2020-FET Open "THOR" Nr. 829067 from the European Commission, and grant IT1164-19 from the Basque Government for consolidated groups of the Basque University.
for (ς † i ) 2 ς 2 i , quadratic in the cavity-sensor coupling parameter i , is found by following the recipe detailed in Section IIIA of Ref. [13] up to Eq. (17), and -adopting the notation from that contribution -defining the following auxiliary conditional states: ρ 2 0 (in vectorized form denoted as |ρ 2 0 ), ρ 2 1 (|ρ 2 1 ) and ρ 2 2 (|ρ 2 2 ) as and identifying For the details of this method and clarification of the notation used, we direct the reader to Ref. [13].
We should note that this autocorrelation could equivalently be calculated by considering the cross-correlation between two identical sensors (ω 1 = ω 2 and ε 1 = ε 2 ) using the algorithm proposed in Ref. [13]. However, we found that, for each point on the map of CSI, this method requires 9 calls to the scipy.sparse.linalg.spsolve procedure to calculate the vectorized density matrices (one for each of Eqs. (25a-h) in Ref. [13]), compared to 5 required for the implementation described above (two for |ρ 0 1 and |ρ 1 1 and 3 for Eqs. (A5-A7)). The Python implementation of this method is available upon request from the corresponding authors.