Two-photon resonance fluorescence of two interacting non-identical quantum emitters

We study a system of two interacting, non-indentical quantum emitters driven by a coherent field. We focus on the particular condition of two-photon resonance and obtain analytical expressions for the stationary density matrix of the system and observables of the fluorescent emission. Importantly, our expressions are valid for the general situation of non-identical emitters with different transition energies. Our results allow us to determine the regime of parameters in which coherent two-photon excitation, enabled by the coherent coupling between emitters, is dominant over competing, first-order processes. Using the formalism of quantum parameter estimation, we show that the features imprinted by the two-photon dynamics into the spectrum of resonance fluorescence are particularly sensitive to changes in the distance between emitters, making two-photon excitation the optimal driving regime for the estimation of inter-emitter distance. This can be exploited for applications such as superresolution imaging of point-like sources.

Here, we focus on a particularly relevant effect that arises when interacting emitters are driven by a classical field: the coherent, non-linear excitation of the transition from the ground state to a doubly-excited state via a two-photon resonance, enabled by the emitter-emitter interaction [10,28]. This mechanism of two-photon excitation, which lies at the heart of important technological applications such as twophoton microscopy [42], has attracted great interest for the generation of squeezing [26], steady state atomic entanglement [27,43], and emission of entangled photons [20]. The emergence of an extra peak in the excitation spectrum due to this two-photon resonance has been demonstrated experimentally [28], used as a method to quantify dipole coupling * carlossmwolff@gmail.com and, indirectly, estimate the distance between quantum emitters with nanometer resolution.
Despite the high fundamental and technological relevance of this mechanism, and the apparent simplicity of the model that describes it, most theoretical studies have focused on the particular case of identical emitters, where one can find analytical expressions for the stationary density matrix of the quantum emitters by direct diagonalization [7,8,43]. A straightforward analytical solution cannot be obtained in the more complicated case of two non-identical emitters (e.g. with different transition frequencies), where only approximated solutions for steady-state populations, limited to the case of very weak coupling and driving strength have been reported [10,27]. The case of non-identical emitters is the common situation in solid-state emitters [27,28,33], and it is a relevant situation for related models describing, e.g., lightharvesting [44] and energy-transfer [29,45]. Given its importance, a full theoretical description of interacting non-identical quantum emitters under coherent driving is desirable.
In this work, we obtain general analytical expressions for the stationary density matrix of two interacting non-identical quantum emitters under coherent driving at the two-photon resonance. Our expressions are valid for an ample regime of parameters, under the only condition that the energy splitting between one-excitation eigenstates must be the largest energy scale in the system. This allows us to provide closed-form expressions in regimes that could not be described by previous analytical results [10], such as large driving strengths that saturate the two-photon transition, and large coupling between quantum emitters. Furthermore, motivated by the experimental work in [28], we establish the potential of resonance fluorescence measurements to estimate the distance between dipoles below Abe's limit of difraction. Using the framework of quantum parameter estimation [46][47][48][49][50][51][52][53], we find that the maximum precision is obtained by driving the system at the two-photon transition in the onset of saturation.
This paper is organized as follows. In Section II, we describe the model of the system and then introduce two effective models that account for second-order and first-order processes independently. In Section III, we apply the results of the previous section to analyze the steady-state observables of the light radiated by the quantum emitters. In Section IV, we describe the spectrum of two-photon resonance fluorescence. Finally, in Section V, we use quantum parameter estimation theory to analyze the potential of spectrum measurements for the estimation of the inter-emitter distance.

A. General model
Our system is composed of two quantum emitters, each of them described as a two-level system (TLS) placed at position r i (i ∈ 1, 2), with natural frequency ω i and dipole moment µ i . The system and energy levels are sketched in Fig. 1(a-c): each TLS spans a basis {|g i , |e i }, where we define the lowering operatorσ i = |g i e i |. The total basis of the composite system is {|gg , |ge , |eg , |ee }, where |gg ≡ |g 1 ⊗ |g 2 , and similarly for the other states. We consider a coherent coupling between both TLS with a coupling rate J. In this work, J will be determined by the dipole-dipole interaction between the emitters, although the same model and analysis can be applied in situations in which the interaction between quantum emitters is modified by the mediation of a photonic structure [22,23,27,31,32,[39][40][41]. We also include a coherent laser field of frequency ω L driving both emitters with a Rabi frequency Ω. In the rotating frame of the laser, the resulting time-independent Hamiltonian isĤ =Ĥ 0 +Ĥ d , whereĤ 0 is the bare Hamiltonian of the interacting quantum emitterŝ having defined ∆ ≡ (ω 1 + ω 2 )/2 − ω L , and δ ≡ (ω 2 − ω 1 )/2, and set = 1. Both quantum emitters interact with the electromagnetic field continuum that surrounds them. This field is responsible for mediating the coherent interaction between emitters in Eq. (1), and also provides a decay mechanism that de-excites the quantum emitters by spontaneous emission to free space. In the reduced Hilbert space of the emitters, this dissipative dynamics is modeled by a master equation for the density matrix [2,54,55], where γ ii is the local decay rate of spontaneous emission of the i-th emitter, and γ 12 = γ 21 is the dissipative coupling rate between emitters that emerges as a consequence of collective decay. We assume optical transitions, and therefore neglect the incoherent excitation by thermal photons. The local decay rates depend on each emitter's natural frequency and dipole moment, For the sake of simplicity and without loss of generality, we will assume that both dipole moments equal, µ 1 = µ 2 . This implies γ ii ≈ γ (assuming ω 1 , ω 2 δ) and justifies our choice of the same Rabi frequency Ω for the driving of both emitters, since, besides having the same dipole moment, their separation will be considered smaller than the resonant wavelength, kr 12 , and therefore both emitters are driven with the same amplitude, Ω(r 1 ) ≈ Ω(r 2 ) ≈ Ω. Nevertheless, all the results that we obtain can be easily generalized to the case of γ 1 = γ 2 , Ω(r 1 ) = Ω(r 2 ). We emphasize that, even if we assume equal dipole moments, we still consider the general case of non-identical emitters that may have different natural frequencies, i.e. δ = 0. The fact that our results apply to this general case is one of the main achievements of this paper.
The coherent and dissipative coupling rates depend on the dipole moments and and also on the relative separation between emitters, where r 12 = |r 12 |, k = ω 0 /c, and ω 0 = (ω 1 + ω 2 )/2, having assumed ω 0 (ω 2 −ω 1 ). We are interested in the case of two very close TLS, i.e., kr 12 1. In this regime, the collective parameters reach their maximal values Henceforth, we adopt Eqs. (6,7) in all our calculations, and also assume that the dipole moments are perpendicular to r 12 , so µ · r 12 = 0. The Hamiltonian of the undriven emitters (1) can be easily diagonalized, yielding a new basis {|gg , |A , |S , |ee } where the eigenstates in single-photon subspace are given by with β being a mixing angle defined as This diagonal basis is depicted in Fig. 1(c). In the case of identical emitters usually discussed in the literature, δ = 0 and therefore β = 0, so that |S and |A are, respectively, purely symmetrical and anti-symmetrical superpositions of the states |ge and |eg (hence the notation used). In this work, we will consider the more general case where β ∈ [0, π/2], which includes the possibility of non-identical emitters, δ = 0. The corresponding energies of these two eigenstates are E S/A = ∆ ± R, where we have defined the Rabi frequency of the dipole-dipole coupling as In this work, we will be particularly interested in the case in which R is the largest energy scale in the system, so thatĤ d can be treated perturbatively with respect to R. This approach is different to the one taken, for instance, in Ref. [10], where J was taken as a perturbative parameter. Our approach will allow us to derive analytical expressions valid in more regimes, such as the one of high dipole-dipole coupling J δ that can be relevant for closely spaced emitters. This will also imply that we focus on a regime where the two resonances ∆ = ±R are well resolved, R γ, and are therefore visible as two distinct peaks in measurements such as resonance fluorescence excitation spectra [28].
Since we are interested in using R as our energy reference [44], we will reformulate the Hamiltonian parameters in terms of R and β as When J is given by the dipole-dipole coupling, as we consider in this work, β will depend on both δ and the distance between emitters, kr 12 . This makes β range from 0, at short distances, to π/2, at long distances, as we show in Fig. 1(d).
In this figure, one can see that, as the detuning δ is decreased, the distance required to reach β = π/2 increases, until one reaches the limiting case δ = 0, where β = 0 for any value of the distance kr 12 . The distance kr 12 has thus as strong impact in the structure of the eigenstates (8), which, as we shall see, will affect the quantum optical properties of the emitted light and the response to coherent, two-photon driving. This panel also shows that Eq. (6) provides a good approximation for J even for values kr 12 ∼ 1.

B. Effective models
The full master equation in Eq. (3) does not yield analytical expressions in general case of δ = 0, which explains why the majority of literature has focused on the case δ = 0. In order to obtain analytical expressions for the system density matrix, we will make the assumption that the dynamics is governed by two different types of processes that take place independently and that can be described by two different effective models. The first process is the resonant, two-photon excitation that drives the |gg → |ee transition via a second-order process, and the subsequent incoherent decay towards |gg , passing through the single-photon subspace. This is described by an effective three-level cascade model, that we denote Model 2P, see Fig. 2(a).
The second mechanism is the excitation of the one-photon subspace {|S , |A } via first-order processes. This is described by a three-level Vee model that we denote Model 1P, see Fig. 2(b), that excludes the doubly excited state |ee .
Our key assumption will be to consider that either both types of processes take place independently and very scarcely, or that the dynamics is completely dominated by one of the processes (e.g. at the two-photon resonance condition ∆ = 0 for the second-order processes, or the one-photon resonance conditions ∆ = ±R for the first-order). In both cases, this means that the probabilities of occupation of excited states and coherences computed from each of these models contribute additively to the total density matrix. Adopting the notation ρ ij ≡ i|ρ|j , we express this as where ρ (2) ij and ρ (1) ij denote matrix elements computed from second-order processes (Model 2P) and first-order processes (Model 1P) respectively. In the following, we detail how these matrix elements are computed from each of these two approximated models.

Model 2P: second-order processes
Model 2P describes the nonlinear process of coherent, twophoton excitation by the driving laser. Therefore, it will predict sizable probabilities for the excited states only around the two-photon resonance, ∆ ∼ 0. For values of ∆ close to zero, |gg and |ee form a quasi-degenerate subspace in the rotating frame of the laser, split from the states, |S and |A by an energy difference ±R, where we assume R to be the largest energy scale in our system, so R ∆. The states |gg and |ee are not coupled to first order, ee|Ĥ|gg = 0, but they are coupled through second-order processes mediated by |S and |A , which couple to both |gg and |ee through the driving laser. We will assume Ω R so that we can describe the effective two-photon coupling rate Ω 2p between |gg and |ee by second-order perturbation theory, withĤ d as the perturbation: Ω 2p can be understood as a two-photon Rabi frequency, featuring a quadratic scaling with Ω and a strong dependence on β, or equivalently, the ratio J/δ. One obvious consequence of Eq. (13) is that, for J = 0, there is not two-photon coupling between |gg and |ee due to destructive interference between the two possible pathways that mediate the interaction. This explains why the emergence of optical features related to the two-photon process stands as a clear evidence of coherent coupling between quantum emitters [10,28]. The states |ee and |gg also experience an effective Lamb shift λ which, notably, is the same for both of them and equal to with j = gg, ee.
Beyond the coherent two-photon driving, the remaining ingredient of the dynamics is the incoherent decay from |ee to |gg , passing through one of the single-photon states, |S or |A . Since this is an incoherent process that populates in an equal manner both single-photon states, we simplify our description by considering a single intermediate single-photon state |1 , whose steady state population gives us the sum of the populations of |S and |A , see Fig. 2(a). The energy of this state is irrelevant in this picture since it is incoherently populated, and thus we set it to zero.
The state of the reduced subsystem {|gg , |ee , |1 } is described by a 3×3 density matrixχ 2p . Its dynamics is governed by the following master equation: whereĤ 2p is the effective two-photon Hamiltonian From this, we can obtain the second-order contributions to the excited-state components of the totalρ-see Eq. (12)establishing the following relations Solving for the steady-state of Eq. (15) gives direct analytical expressions for the elements ofχ 2p , and therefore forρ (2) . These read:

Model 1P: first-order processes
Our second model describes dynamics in which the singlephoton states |S and |A are directly excited by the driving field through a first-order process. In essence, Model 1P consists in removing the state |ee from our description, thus neglecting the two-photon excitation mechanisms that are described by Model 2P, and thus allowing only for first-order processes to occur. The result is a three-level Vee system comprising the basis states {|gg , |A , |S }. The Hamiltonian in this reduced model readŝ The driving rates Ω S and Ω A are simply given by S|Ĥ d |gg and A|Ĥ d |gg respectively, and are obtained directly from Eqs. (2) and (8), In the usually considered situation of identical emitters, β = 0, one finds that Ω S = √ 2Ω and Ω A = 0, so that the antisymmetric state is dark and decoupled from the driving field. Considering now spontaneous emission in the basis of |S and |A , we find the following master equation for the 3 × 3 density matrixχ 1p of Model 1P: where we have defined the following decay rates The rates γ S/A describe the standard decay of |S/A towards |gg , while γ C is the rate of incoherent coupling between |S and |A that originates from their collective decay. Describing the one-photon dynamics in this basis allows us to easily distinguish between enhanced and suppressed channels of decay. For instance, focusing on γ A , one can see that γ A = 0 for two close, identical emitters (β = 0, γ 12 = γ), meaning that the antisymmetric state |A gets completely decoupled from the dynamics and turns into a dark state whose population gets trapped [56][57][58][59]. Similarly, in that situation one finds γ S = 2γ, which clearly shows the superradiant nature of the symmetric state |S . The master equation (21) also yields analytical solutions for the stationary state. Once the master equation (21) is solved, we map the elements of the reduced density matrixχ 1p to the first-order contributions to the totalρ, which we denoteρ (1) . For elements of the one-photon subspace {|S , |A }, this map is simply given by The resulting first-order contributions to total steady-state density matrix are thus given by: and The expression provided in Eq. (25) is an approximation that we obtained by eliminating |A from Model 1P when calculating ρ S,S , and viceversa (effectively working with twolevel systems {|g , |S/A } instead of the three-level system of Model 1P). This approximation is not necessary to obtain an analytical expression from the master equation for the threelevel system in Eq. (21); however the full analytical expression of ρ (1) S,A is long and cumbersome to write, and it is very well matched by the more compact and tractable expressions in Eq. (25).
In principle, Model 1P ignores the excitation of the twophoton subspace |ee via successive, first-order excitation processes. This is justified since, focusing around the region we are most interested in, the two-photon resonance ∆ ≈ 0, the contribution of this mechanism to ρ ee is rather small compared to that of second-order processes. However, when one approaches the limit β ≈ π/2, the two-photon population originated by second-order processes, ρ (2) ee,ee , tends to zero, as can be seen from Eq. (18a), and the small contribution from first-order processes may dominate and become relevant. In that particular regime, we can infer the value of the first-order contribution ρ (1) ee,ee even if the state |ee is not included in Model 1P. This can be done by noting that, for β ∼ π/2, the two quantum emitters are essentially decoupled, and the population ρ ee,ee stems from the simultaneous but indepen-dent excitation of both emitters. This yields factorizable correlations of the type σ + 1σ + 2σ 1σ2 = σ + 1σ 1 σ + 2σ 2 . Since ρ ee,ee = σ + 1σ + 2σ 1σ2 , this factorization allows us to use the expressions of the populations obtained from Model 1P, σ + iσ i (1) , to estimate first-order contributions to the occupa- (1) , even if |ee was not explicitly included in the model. For β ∼ π/2, we have that σ + 2σ 2 σ + 1σ 1 ≈ ρ S,S ρ A,A . Therefore, we will define This expression is necessary to regularize the expected twophoton population ρ ee,ee in the limit of uncoupled emitters, β = π/2, and basically unimportant in any other case.

III. STEADY STATE OBSERVABLES OF THE EMITTED LIGHT
Following the scheme of Eq. (12), we can now combine the results provided by the effective models just discussed and obtain an estimate of the total steady-state density matrixρ. The approximations used are expected to hold particularly well around the region of interest ∆ ∼ 0, i.e. the two-photon resonance. In particular, the most relevant density matrix elements for subsequent calculations read: These density matrix elements allow us to obtain steady-state observables of the fluorescent light emitted by the system, which is done by establishing a connection between the radiated electric field operator and the annihilation operators of the quantum emitters. For instance, if we consider the emitters to be located at the origin and to have the same dipole moment µ (as we do throughout this work), the positive frequency part of far-zone electric field operator radiated by the two quantum emitters is given byÊ (+) (r, t) =Ê (+) (r, t)u x , where u x is a unit vector perpendicular to r and contained within the plane spanned by µ and r [60], and Other detection schemes can give a different spatial distribution of the field, e.g. imaging by focusing the field with a lens yields a distribution given by its point spread function. Nevertheless, the relationship will hold provided the dipole moments are equal and their separation is small compared to the resonant wavelength of the field kr 12 1 [61]. This regime is of particular interest for our work, since a good understanding of the system dynamics and quantum optical properties of the emission becomes essential to infer inter-emitter distances below the diffraction limit [28], which has important applications for microscopy and superresolution imaging [62]. Therefore, we will assume the proportionality relation (30) throughout the text.
The two main observables of the flourescent emission that we will analyze here are mean intensity and second-order correlation function.

A. Mean intensity
From the relation (30), we see that the mean intensity of the signal is given by Ê (−)Ê(+) ∝ Î , where we have defined the intensity operator The steady-state mean value of the intensity operator I ≡ Î can be expressed in terms of the density matrix elements as Equation (32), together with Eqs. (18a, 18b), Eqs. (25)(26)(27) and Eqs. (28a-28d), provides a direct analytical expression for I. Our analytical results are shown in comparison with numerical results in Fig. 3. Panel (a) depicts the intensity I versus the laser detuning ∆; (b) shows I versus both ∆ and β, and (e) represents a zoom of (a) around the two-photon excitation regime, ∆ ∼ 0. There are three characteristic highintensity peaks corresponding to the values ∆ = {−R, 0, R}. The two peaks at ±R correspond to the resonant excitation of the states |A/S and has therefore a first-order origin described by Model 1P. Their relative height depends strongly on β (i.e. on the ratio between the dipole-dipole coupling and the relative detuning between the emitters), as clearly seen in Fig. 3(b). Model 1P provides a good match with the exact, numerical calculation provided that the occupation probabilities are small enough for the simultaneous excitation probability to be negligible. When β = 0, i.e. the usually considered situation of resonant (identical) emitters, only the symmetric state |S gets significantly populated at ∆ = −R, while the resonance with the |A state at ∆ = R is suppressed. The reverse situation is found for β = π/2, where both emitters are uncoupled and thus only classical correlations between them can exist. In this situation, the curve is symmetric around ∆ = 0 and the two peaks at ∆ = ±R, which in this case correspond to the resonant driving of each of the independent QEs, have similar values.
The central peak at ∆ = 0 is arguably the most relevant feature given its implications for microscopy and imaging [28]. This peak emerges from the resonant two-photon excitation enabled by the coherent coupling between emitters, and thus is fully described by the contributions of second-order processes from Model 2P, Eqs. (18a, 18b). Our analytical solution allows us to establish that both the height and width of this peak scale as Ω 4 cos 2 β. As expected, the peak vanishes for uncoupled emitters (β = π/2) as a consequence of the destructive interference of the two-excitation pathways, which yield Ω 2p = 0. Our analytical expression of ρ ee,ee (18a) gives us the possibility to use the intensity of the two-photon peak at ∆ = 0 to infer the value of β. On the other hand, R can also be easily estimated from the position of the one-photon peaks. The knowledge of these two magnitudes can then be combined to obtain information about the internal structure of the quantum emitters, i.e., their natural energy detuning 2δ and coherent coupling rate J. In turn, this allows one to infer quantities such as the inter-emitter distance kr 12 , highly relevant for technological applications like superresolution imaging [28].
Given its importance, it is desirable to determine the set of conditions under which the two-photon peak will be visible. Following our approach of separating contributions from firstand second-order processes, the total intensity can be written as I = I (1) + I (2) . The two-photon peak arises from the resonant contribution I (2) , while the off-resonant, first-order contribution I (1) gives a featureless background at ∆ = 0 which, under certain conditions, can be brighter than the twophoton peak and hide it. In particular, since the population of second-order origin ρ peak have quartic scaling with Ω, while first-order processes yield populations that scale quadratically with Ω, there must be a value Ω v below which first-order processes dominate, see Fig. 4(a). To determine Ω v , we can define a two-photon visibility V 2p as the ratio V 2p = I (2) /I (1) , such that the twophoton peak will be visible when V 2p > 1. I (2) and I (1) can be computed from the first-and second-order contributions to ρ ee,ee , ρ S,S , ρ A,A and ρ S,A , using Eqs. (18a, 18b) and Eqs. (25)(26)(27) respectively. Introducing these into the equation V 2p = 1 and solving it, we obtain a value for the minimum necessary driving amplitude Ω v that guarantees two-photon visibility. The equation is greatly simplified in the regime that we consider in this paper, Ω, γ R, yielding where the last approximation applies provided γ R tan β. These results are summarized and confirmed by exact numerical calculations in Fig. 4. In panel (a), we show that the sum of our analytical estimations of I (1) and I (2) recovers the exact value of I computed numerically, which, as discussed above, features a transition from a ∝ Ω 2 to a ∝ Ω 4 scaling at Ω v , which marks the onset of visibility of the two-photon peak, i.e., the emergence of features characteristic of the two-photon dynamics. The two-photon visibility V 2p is shown in the full (γ, Ω) space in panel (b), where the approximated expression for Ω v provided in Eq. (33) is shown to match perfectly the condition V 2p = 1.
To end our discussion on the mean intensity of the radiated field, we observe that a destructive interference also appears for values β = π/2. This can be seen in Fig. 3(b), where a destructive interference dip manifests when the laser detuning is ∆ = R cos β. However, unlike the particular case β = π/2, this is not a destructive quantum interference between excitation pathways in the internal system dynamics, but an optical one, taking place in the radiated electric field and well described by the interference terms appearing in Eq. (32). Indeed, using our Model 1P, we find that the dip is given by point where all the one-photon subspace terms proportional to ρ S,S , ρ A,A and ρ S,A in Eq. (32) add up to zero. As we discuss below, the small amount of remaining light retains a very strong two-photon character.

B. Second-order correlation function
The zero-delay second-order correlation function is defined as g (2) (0) = Ê (−)Ê(−)Ê(+)Ê(+) / Ê (−)Ê(+) 2 , and thus it can be written as This value quantifies the probability of detecting two photons simultaneously, normalized by the probability of doing so in a classical coherent field of similar intensity. As seen in Eq. (34), in our case this is directly related to the probability of occupying the doubly excited state |ee . The numerical and analytical results are summarized in Fig. 3(c-d). Our approximated analytical methods are able to describe accurately the population of the two-photon subspace-and therefore g (2) (0)-around the two-photon resonance ∆ ∼ 0, provided Ω R. This can be seen more clearly in the zoom around the two-photon resonance depicted in Fig. 3(f). Away from the two-photon resonance, the population of the doubly excited state becomes much smaller and it is established by a more complicated mixture of one-photon and two-photon processes that our approximated models fail to capture. These small occupations of |ee , nevertheless, contribute very little to the actual intensity of radiation emitted, which away from ∆ ∼ 0 is dominated by the one-photon subspace and thus is well described by our models, c.f. Fig. 3(b). Interestingly, the maximum values of g (2) (0) are found at the dip of destructive interference discussed before, where the first-order contributions to the total emission interfere destructively and thus the small amount of remaining emission stems mainly from second-order processes, yielding a strong probability of detecting two photons.
At the two-photon resonance ∆ = 0, the value of the g (2) (0) is sharply reduced and tends to 1 from above as Ω increases, c.f. Fig. 3(f). The reason for this is that, as Ω increases, the light emitted is less coherent ( σ i → 0), and thus the emission converges to that of two incoherent quantum emitters [44]. The dip of g (2) precisely at ∆ = 0 is a typical feature of multi-photon processes at resonance [63], and it is simply a consequence of the increased intensity of emission.

IV. SPECTRUM OF TWO-PHOTON RESONANCE FLUORESCENCE
The fluorescence spectrum of the emitted radiation yields useful information about the energy transitions that can take place among the dressed states of the hybridized light-matter system. The most celebrated example of the revealing character of this type of measurement is the Mollow triplet spectrum in the emission from a two-level atom [64]. Its characteristic lineshape with three peaks provides key information about the structure of dressed energy levels [65], and it serves a source of strongly correlated non-classical light [66][67][68][69][70][71].
The resonance fluorescence spectrum of two interacting quantum emitters displays a more complex structure than the Mollow triplet, which has been reported for the case of identical emitters [8,18], and for the analogous case of a coherently driven three-level system at the two-photon resonance [72]. The resulting spectrum has a seven-peaked structure with a central peak and six sidebands. Here, we recover this same result for the case of identical emitters [β = 0, see Fig. 5(a)], but extend it for non-indentical emitters (β > 0), where we show that the spectrum develops an even more complicated structure with up to 13 peaks [see Fig. 5(b)], that end up converging to the three peaks characteristic of the Mollow triplet for completely uncoupled emitters, β = π/2 [see Fig. 5(c)]. In the perturbative regime Ω R, we are able to describe the position of these peaks from the hybridization between the quantum emitters and photon pairs using the dressed energy levels described by Model 2P. The understanding and effective Hamiltonian obtained from this effective model allows one to describe analytically the location and origin of the spectral resonances in the perturbative regime Ω R, and also to obtain insights about this structure beyond that regime, when Ω > R, which we explore numerically.
The spectrum of emission is given by the Fourier transform of the two-time correlation function of the radiated electric field Ê (−) (t)Ê (+) (t + τ ) . For convenience, we will define a general spectral function By using the quantum regression theorem to express the two-time correlation function [73] and formally integrating Eq. (35), we can write S(ω;Â,B) in a computationally con- Resonance fluorescence spectrum at the two-photon resonance ∆ = 0, for β = 0 (a), β = π/4 (b) and β = π/2 (c). In (a-c), upper plots depict the spectrum at Ω/R = 1 (solid blue lines) and the different contributions in its expansion in Eq. (50) (solid purple and dashed lines) , allowing to see cases in which peaks are not visible due to destructive interference, but visible if only the emission from a single emitter is collected. The lower panels in (a-c) show the emergence of sidebands as Ω/R increases; in the perturbative regime Ω R, these correspond to transitions between two-photon dressed states. Positive-frequency peaks are marked with dashed lines, which are identified with the corresponding transition between the dressed-atom ladder of eigenstates in panels (d,f,g,i). Panels (e,g,i) depict the structure of the eigenstates at Ω/R = 1, in the basis {|S , |A2 , |S2 , |A }. Parameters: in upper plots of (a-c), γ/R = 0.1. In the lower panels of (a-c), γ/R = 10 −3 , γ12 = 0.999γ venient form [74], whereρ is the steady-state density matrix, and L is the Liouvillian superoperator that represents the master equation [54], whereρ has a vectorial form, ∂ t |ρ = L|ρ . Considering the relation Eq. (30) between the radiated field and the raising/lowering operators of the quantum emitters, and disregarding global factors, we can express the spectrum of resonance fluorescence as S(ω) = S(ω;σ + 1 +σ + 2 ,σ 1 +σ 2 ).
The spectra computed as a function of Ω for different values of β are shown in Fig. 5, for a laser detuning at the two-photon resonance, ∆ = 0. Perturbative regime-The position of all the peaks in the spectra can be obtained from the possible transitions among energy levels in the system Hamiltonian. Within the perturbative regime Ω R, at ∆ = 0 the ground and excited states |gg and |ee are resonantly coupled by the two-photon excitation described by the Hamiltonian in Eq. (16): these two states then hybridize into a symmetric and antisymmetric states that we denote |S 2 and |A 2 , defined as The resulting set of eigenstates is given by {|S , |A 2 , |S 2 , |A }, ordered by decreasing energy, and their correspondent eigenenergies in the rotating frame of the laser are: In the same way that the eigenenergies of A 2 and S 2 take into account the Lamb shifts induced by the coupling to |A and |S -see Eq. (14)-, E A and E S also include the Lamb shifts of states |S and |A due to their coupling to |ee and |gg , given by with j = S, A. This gives λ S/A = ±2(1 − cos β)Ω 2 /R. The energy differences between these eigenvalues give us the transition frequencies that can observed as distinct peaks in the spectrum. Defining the transition energies ω i→j ≡ E i − E j , the position of the six positive-frequency sidebands with respect to the laser frequency are defined, in the perturbative regime, by the following equations: The sidebands at negative frequencies come from the reversed processes outline above, ω i→j = −ω j→i . The central peak at ω 0 = 0 is given by transitions between similar states, ω i→i . These equations provide the position of all the possible 13 peaks that can be observed in the fluorescence spectrum within the perturbative regime Ω R. The six positive sidebands in Eqs. (44)(45)(46)(47)(48)(49), however, are not visible for all values of β. The case of coupled identical emitters β = 0 features only three positive sidebands, yielding a total of seven peaks as has been noted before [8,18]. In particular, the peaks that are not visible are those that involve transitions that start or end at the one-photon antisymmetric state |A , which for the perturbative regime correspond to the peaks ω 2 , ω 5 and ω 6 in Eqs. (44)(45)(46)(47)(48)(49). The reason why these peaks are not visible is due to the destructive interference phenomena that takes place from the equal contribution from both quantum emitters to the radiated electric field, c.f. Eq. (30). This can be seen in Fig. 5(a). Indeed, expanding the expression of the total spectrum of emission, Eq. (37), we obtain: where S ij (ω) ≡ S(ω;σ + i ,σ j ) and S i (ω) ≡ S ii (ω). S 1 (ω) and S 2 (ω) describe the spectrum of emission that would be obtained by detecting only the radiation emitted by the QE 1 and 2, respectively. On the other hand, S 12 (ω) and S 21 (ω) are interference terms arising from the superposition of both fields. In Fig. 5(a), we show that the missing peaks for β = 0 (all involving the state |A ) are indeed visible if S 1 (ω) or S 2 (ω) are measured independently (e.g. if their emission is collected locally), but they interfere destructively if the fields emitted by both QEs are superimposed, explaining the absence of these peaks in the total spectrum. For β > 0, this perfect destructive interference does not occur, and all the possible 13 peaks are visible in the spectrum of emission. Finally, the opposite limit of completely decoupled nonidentical emitters, β = π/2, yields a three-peaked structure which corresponds to the Mollow triplet of emission from the two independent emitters. In this limit, many of the transitions in Eqs. (44)(45)(46)(47)(48)(49), become degenerate, so that only three peaks can be seen. The reason for having fewer peaks here is obviously different from the destructive interference seen at β = 0, since it responds to the internal structure of transition energies available within the dressed energy levels.
Away from β = π/2, the frequencies in Eqs. (44)(45)(46)(47)(48)(49) describe transitions between dressed light-matter states in which the emitters are hybridized with photon pairs. This strong twophoton character is evidenced by the quadratic scaling of these frequencies with Ω, instead of the linear scaling with Ω that one finds, e.g., for the position of the sidebands in the standard Mollow triplet . For this reason, we refer to the sidebands described by Eqs. (44)(45)(46)(47)(48)(49) as two-photon sidebands. In order to be able to observe two-photon sidebands, we need the energy separation ∆E = ω 2 − ω 3 = ω 4 − ω 5 = 2Ω 2p to be larger than the decay rate of spontaneous emission, γ. The condition 2Ω 2p > γ yields the following equation for the two-photon saturation amplitude Ω 2PS that marks the onset of the resolved two-photon sideband regime: The previous expression tells us that Ω 2P S R, meaning that the two-photon sidebands can be developed within the perturbative regime, provided that the Rabi splitting is large so that R γ, and that β is not too close to π/2. As β approaches π/2, larger values of Ω are necessary to resolve the two-photon sidebands. When values Ω 2PS ∼ R are reached, the perturbative approach used below does not apply, meaning that sidebands developed purely from twophoton hybridization are no longer observable since first-order processes dominate before the former are visible. Setting the condition Ω 2PS ≈ R, we find that this would occur approximately at a mixing angle β max ≈ arccos(γ/4R). One can see, however, that for the typical values considered in this text, i.e. γ = 10 −3 R, β max ≈ 0.9998π/2, i.e. two-photon sidebands should be visible for most of the range β ∈ [0, π/2]. As we will see, the two-photon saturation amplitude Ω 2PS has a great importance for metrological applications, since the onset of the resolved two-photon sidebands regime marks the point of maximum sensitivity for optical estimations of the interemitter distance.
Regime of strong driving-Beyond the perturbative regime, the states {|A , |S 2 , |A 2 , |S } are mixed by the driving and no longer represent the eigenstates of the system. In this situation, tractable analytical expressions for the eigenstates can only be obtained for limiting cases, e.g. β = 0 or β = π/2.
The bottom row of Fig. 5 shows the eigenstates computed numerically for Ω = R; the eigenstates are labeled |U i , with i = 1 . . . 4, in order of decreasing energy. The composition of the eigenstates depends on β, and for the particular cases β = {0, π/4, π/2}, it can be written as follows: due to the strong coherent drive. In this case, the eigenstates can be obtained analytically. To make the resulting expressions more readable, we write them here in non-normalized form: One can see that, in the limit Ω R, these equations tend to the perturbative basis used above, for β = 0. These eigenstates have the following eigenvalues: which, in principle, would yield six sidebands. As discussed in the perturbative analysis, destructive interference between the emission from both emitters make all the transitions involving the antisymmetric state |A invisible. These expressions for the eigenvalues also show that the eigenstates E 3 and E 4 cross in energy when Ω = R/ √ 2.
2. (β = π/4, ∆ = 0). This limit is not easily tractable analytically, so we limit our discussion to a description of the results from numerical calculations. As β increases, |U 1 remains mostly a superposition of the two symmetric states, with a very small component of |A , and |U 3 and |U 4 mix, meaning that |A is no longer an eigenstate of the system. This is an important observation, since it explains why the perfect destructive interference occurring for β = 0 for transitions involving |A no longer take place, and 13 peaks are visible. The eigenstates have the following structure: where the C i,j (i = 1, 3, 4, j = S, A, S 2 ) represent generic amplitudes whose numerically computed values can be seen in Fig. 6(g).
3. (β = π/2, ∆ = 0) This is the limit of detuned, uncoupled emitters, that recovers the physics of two independent, detuned Mollow triplets. The eigenstates have the following analytical form: with and eigenvalues A clear conclusion of this analysis is that the spectrum resonance fluorescence is strongly sensitive to the ratio J/δ (expressed here in terms of the mixing angle β), and therefore it can serve as a valuable source of information about the parameters that characterize the system of quantum emitters, such as the inter-emitter distance. The continuous variation of the spectral features with β is depicted in Fig. 6 for different values of ∆. These plots reproduce the emergence and disappearance of peaks with β described above, and furthermore show that similarly complex structures emerge outside the two-photon resonance, i.e. for ∆ = 0, giving distinct patterns on the excitation-emission spectra that also depend strongly on β.

V. QUANTUM PARAMETER ESTIMATION OF INTER-EMITTER DISTANCES
The results just discussed suggest that the measurement of excitation-emission spectra could provide valuable information for the estimation of internal parameters such as the value of the coherent coupling between emitters or, equivalently, their distance in real space. In this section, we address this question by establishing the metrological potential of these measurements within the formalism of quantum parameter estimation [46][47][48][49][50][51][52][53]. To do this, we consider the precision achievable in the estimation of an unknown parameter, X, from the outcome of a positive operator-valued mesaurement (POVM) Λ, consisting of a set of operators {Λ µ }, where the index µ ∈ {1, 2, . . . , M } denotes different possible measurement outcomes. The elements of the POVM add up to the identity, µ Λ µ = 1. The probability for each of the different measurement outcomes follows a distribution P (µ|X) = Tr[ρ X Λ µ ], where ρ X is the X-dependent density matrix of the system at the time of the measurement. When one infers X from P (µ|X), the best attainable sensitivity ∆ 2 X that can be achieved is given by the classical Fisher information of the probability distribution, The Cramér-Rao bound then establishes that the minimum possible variance ∆ 2 X in the estimation of X is given by the inverse of F [48], Here, we will focus on the particular example of the estimation of the inter-emitter distance kr 12 by measuring the fluorescence spectrum described in the previous section. In order to do this, we will make a series of assumptions. First, we will assume that the value of the spectrum at a frequency ω is obtained by performing a photon-counting measurement in a bosonic sensor resonant at that frequency and weakly coupled to the quantum emitters. For instance, this sensor can be understood as a tunable Fabry-Perot cavity used as a frequency filter. We assume that a discrete set of N frequencies ω n are measured. The corresponding POVM is given by the tensor product of photon-number operators of the different sensors,n 1 ⊗n 2 . . . ⊗n N , and the probability distribution that describes the measurement outcomes is of the form P (n 1 , n 2 , . . . , n N |kr 12 ), representing the joint photon-counting probability in each of the sensors. We then assume that the measurements done at different sensors are uncorrelated, so that P (n 1 , n 2 , . . . , n N |kr 12 ) = P (n 1 |kr 12 )P (n 2 |kr 12 ) . . . P (n N |kr 12 ). This ignores the possible contribution to the Fisher information of the correlations between photons emitted at different frequencies [66,71,75,76]. This assumption is justified if the lifetime of photons within the sensor is very long, so that temporal correlations are lost, or if different frequencies are measured sequentially in independent experiments, e.g. tuning the frequency of a Fabry-Perot filter. Following the same approach that is used in the theory of quantum image processing [77], we assume that the resulting state of the sensor is a coherent state, so that the corresponding measured photo-current displays Poissonian fluctuations. This means that the photon-counting distribution associated to the sensor of frequency ω is given by P (n ω |kr 12 ) = n ω e − nω n ω ! , where n ω = ηS(ω), with η a global constant that depends on the particular details of the detection scheme (e.g. detection efficiency). We set η = 1 for simplicity, since it only yields an overall factor. Given that the spectrum S(ω) is strongly dependent on the value of β, it will also vary strongly with the the inter-emitter distance kr 12 , which we emphasize by writing S(ω) = S(ω, kr 12 ). Since different sensors are uncorrelated and probability distributions factorize, the Fisher information associated to the measurement of the spectrum S(ω, kr 12 ) is given by the sum where we used Eq. (75). This quantity allows us to evaluate the metrological potential of fluorescence spectrum measurements. For this calculation, we take into account a finite detector linewidth Γ = γ in the spectrum, [75], done by the replacement ω → ω + iΓ in Eq. (36). Our results are summarized in Fig. 7. Fig. 7(a) depicts the Fisher information versus the optically tunable parameters ∆ and Ω, showing that the optimal regime of operation is at the two-photon resonance, ∆ ≈ 0, where F is found to be larger. This is explained by the fact that the mechanism of two-photon excitation is strongly dependent on the coherent coupling between emitters, as we have seen in previous sections, and this is strongly modified by the inter-emitter distance kr 12 . At ∆ = 0, we find that there is an optimum value of the driving amplitude Ω that maximizes F and consequently the precision in the estimation of kr 12 . As we show in Fig. 7(b), this maximum varies with the actual value of kr 12 , and it is well approximated by the driving amplitude of two-photon saturation Ω 2PS , that we obtained in Eq. (51), at which the two-photon sidebands begin to be resolved in the spectrum. This establishes the onset of the two-photon saturation regime as the most sensitive point of operation for the estimation of the distance between interacting quantum emitters. We note that the results shown here as a function of the optically tunable parameters ∆ and Ω represent different, independent experiments. Thus, a series of measurements over the (Ω, ∆) parameter space could provide a higher precision of estimation, with a Fisher information that would be given by the sum F = ∆,Ω F (Ω, ∆).

VI. CONCLUSIONS
In this work, we have studied a system of two interacting, non-identical quantum emitters under coherent driving, focusing particularly on the regime of two-photon excitation. We have provided, for the first time for non-identical emitters, analytic approximations of the stationary density matrix and of steady-state observables such as the intensity of fluorescent emission. These calculations provide valuable insights on how the properties of the light emitted depend on degree of the coherent coupling and the detuning between emitters, and allows us to establish the regime of parameters in which specific features of two-photon excitation, such as the characteristic two-photon resonance peak in the excitation spectrum or two-photon sidebands in resonance fluorescence, are visible. Given that this features are strongly dependent on the coupling strength between emitters, and thus on their relative distance, we have explored the potential of these optical measurements for the estimation of the inter-emitter distance in terms of their Fisher information. We have established that the onset of two-photon effects at the two-photon resonance is the most sensitive point of operation for the estimation of the inter-emitter distance, a result that can be of great relevance for the problem of imaging beyond Abe's resolution limit [28].