# Tailoring the statistics of light emitted from two interacting quantum emitters

Adrián Juan-Delgado, Ruben Esteban, Álvaro Nodar, Jean-Baptiste Trebbia, Brahim Lounis, Javier Aizpurua

► To cite this version:

Adrián Juan-Delgado, Ruben Esteban, Álvaro Nodar, Jean-Baptiste Trebbia, Brahim Lounis, et al. Tailoring the statistics of light emitted from two interacting quantum emitters. Physical Review Research, 2024, 6 (2), pp.023207. 10.1103/PhysRevResearch.6.023207. hal-04610147

HAL Id: hal-04610147

https://hal.science/hal-04610147v1

Submitted on 17 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Tailoring the statistics of light emitted from two interacting quantum emitters

Adrián Juan-Delgado,¹,*, Ruben Esteban,¹,² Álvaro Nodar,¹ Jean-Baptiste Trebbia,³,⁴
Brahim Lounis,³,⁴ and Javier Aizpurua,²,⁵,⁶,†

¹Materials Physics Center, CSIC-UPV/EHU, 20018 Donostia-San Sebastián, Spain

²Donostia International Physics Center (DIPC), 20018 Donostia-San Sebastián, Spain

³Université de Bordeaux, LP2N, F-33405 Talence, France

⁴Institut d'Optique & CNRS, LP2N, F-33405 Talence, France

⁵IKERBASQUE, Basque Foundation for Science, 48009 Bilbao, Spain

⁶Department of Electricity and Electronics, University of the Basque Country, 48940 Leioa, Spain

<figure><img src="image_2.png" /></figure>

(Received 26 October 2023; accepted 18 April 2024; published 24 May 2024)

The interaction between quantum emitters leads to the formation of superradiant and subradiant states with possible applications in quantum technologies. To improve the characterization of light emission from these systems, we present here a systematic theoretical analysis of the intensity correlation from two strongly interacting quantum emitters at cryogenic temperatures as a function of the frequency and intensity of the excitation laser. This analysis effectively accounts for the effect of vibrational modes of the emitters and of phonons of the environment through the combined Debye-Waller/Franck-Condon factor. First, we analyze the color-blind intensity correlation and show that it can be tailored from strong antibunching to strong bunching by tuning the laser from the two-photon resonance to the transition frequency of the superradiant state. We also find a particularly complex behavior of the intensity correlation when the laser frequency is tuned to that of the transition of the subradiant state, giving raise to the possibility of emitting bunched and antibunched light depending on the laser intensity and the detuning between the two emitters. The numerical results are supported by analytical equations that can be used for the experimental characterization of the interacting emitters. Additionally, by selecting photons of particular frequencies, we analyze the rich landscape of frequency-resolved intensity correlations, which also depend on the laser detuning and intensity. The analysis of the frequency-resolved correlations provides further information about the different relaxation processes underlying the photon emission, unveiling one-photon and two-photon emission processes that cannot be resolved neither in the emission spectrum nor in the color-blind intensity correlation. These results show that two interacting emitters are a versatile and practical source of quantum light and highlight the usefulness of the intensity correlation to unveil complex dynamics in this system.

DOI: 10.1103/PhysRevResearch.6.023207

I. INTRODUCTION

Quantum emitters that behave as two-level systems (TLSs) constitute efficient single-photon sources [1], which are key ingredients in quantum cryptography [2–4], quantum teleportation [5,6], and quantum computing [7,8]. Some examples of this kind of emitters are organic molecules at cryogenic temperatures [9–12] and solid-state emitters [13–17], among others [1]. The control and manipulation of these emitters is thus of great interest in quantum information and quantum communication. Additionally, the interaction between quantum emitters leads to cooperative effects that modify the resonant frequencies and the coherence times of the in-

dependent emitters, with the emergence of subradiant and superradiant states [18]. In this context, recent experiments have demonstrated the capability to manipulate the degree of delocalization between two organic molecules [19,20], as well as the capacity to couple two quantum emitters via a waveguide or a microresonator [21–23]. These interacting systems can be further characterized by analyzing their light emission.

A fundamental tool to interrogate the nature of light emitted from two interacting emitters is the (color-blind) intensity correlation $g^{(2)}(\tau)$ [24,25]. This quantity can be measured with a Hanbury-Brown Twiss interferometer [26], as schematically illustrated in Fig. 1(a) (without including the optical filters in the scheme). $g^{(2)}(\tau)$ provides information on the statistical properties of the light emitted, allowing for distinguishing between antibunched, coherent, and bunched sources of light [24]. For example, intensity correlation has been used to demonstrate that single organic molecules and single quantum dots can be efficient single-photon sources, since they exhibit a strongly antibunched emission [9–17]. On the other hand, if two quantum emitters are coherently and strongly coupled, experiments have shown that the light emitted can be either bunched or antibunched depending on the

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*adrianjuand1996@gmail.com</footnote>

<footnote>†aizpurua@ehu.eus</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Hanbury-Brown Twiss interferometer and electronic eigenstates of two TLS-like emitters. (a) Scheme of the Hanbury-Brown Twiss interferometer, where two optical filters (surrounded by dashed boxes) are included to obtain the frequency-resolved intensity correlation $g^{(2)}(\omega_1, \omega_2; \tau)$ of light emitted from two interacting emitters, such as polycyclic aromatic hydrocarbon molecules embedded in a naphthalene film. The color-blind intensity correlation $g^{(2)}(\tau)$ can be calculated by removing the optical filters. [(b),(c)] Schematic representation of the energy levels of the undriven emitters (b) when they are not interacting and (c) when they do interact. The emitters have transition dipole moments $\mu_i$, which are assumed to be identical, and transition frequencies $\omega_{e,i}$ from the ground state $|g_i\rangle$ to the excited state $|e_i\rangle$. The transition frequencies are detuned by $\delta = \omega_{e,1} - \omega_{e,2}$ from each other for the noninteracting system, and the spontaneous decay rate $\gamma_0$ is identical for both emitters. When the dipole-dipole coupling $V$ is included, the energy splitting between the delocalized states $|\Lambda_+\rangle$ (higher energy) and $|\Lambda_-\rangle$ (lower energy) becomes $2\Lambda = \sqrt{4V^2 + \delta^2}$ and the delocalized states in general have different decay rates $\gamma_+$ and $\gamma_-$. For the configuration considered in the calculations (molecular J-aggregate like, with the dipoles aligned and parallel), $|\Lambda_-\rangle$ is the superradiant state and $|\Lambda_+\rangle$ is the subradiant state.</figcaption></figure>

laser frequency and intensity [19,27]. However, the theoretical analyses of the intensity correlation from two coupled emitters have focused on the case of a laser tuned to the two-photon resonance, and usually the effect of emitter vibrations and of phonons of the environment are not included in the description [28–33].

Further insights into the emission of the system can be gained with the frequency-resolved intensity correlation (FRIC) $g^{(2)}(\omega_1, \omega_2; \tau)$ [34]. This quantity measures the correlation between pairs of photons of particular frequencies $\omega_1$ and $\omega_2$ by incorporating optical filters in the Hanbury-Brown Twiss interferometer [see Fig. 1(a)]. The direct calculation of this quantity is difficult since it involves complicated

four-dimensional time integrals [35,36], but it can be simplified assuming a Lorentzian transmission of the optical filters [36,37], as obtained with Fabry-Pérot filters. In this way, $g^{(2)}(\omega_1, \omega_2; \tau)$ has been investigated for different systems such as a single quantum dot [38], a biexciton quantum dot in a cavity [39], and cavity optomechanics systems [40]. FRIC has been calculated for two interacting TLSs, for the particular case where they are identical and the laser is tuned to the two-photon resonance [41].

Here, we present a detailed analysis of the color-blind and the frequency-resolved intensity correlation of light emitted from two strongly interacting TLS-like emitters. Notably, we include in this theoretical analysis the combined Debye-Waller/Franck-Condon factor $\alpha$, which accounts for the proportion of photons that are emitted in the purely electronic transition [19,42], i.e., it accounts in an effective manner for the internal vibrations of the emitter as well as for the phonons of the environment. We demonstrate in Sec. III that a very rich landscape of color-blind intensity correlations can be obtained by considering different detunings between the emitters, as well as different laser frequencies and intensities. In particular, we analyze the dependence of the intensity correlation on the emitters detuning and laser intensity for three different values of laser frequency corresponding to (i) the laser tuned to the two-photon resonance, (ii) the transition frequency of the superradiant state, and (iii) the transition frequency of the subradiant state. We show that the intensity correlation is drastically different for the three laser detunings, which allows for tailoring the statistics of the emission from extreme anti-bunching to extreme bunching. The results are analyzed with the help of analytical equations, which could be applied to experimentally estimate parameters such as the distance and the detuning between the two emitters. We then investigate in Sec. IV the FRIC for these three laser detunings and for weak, moderate, and strong laser intensities. We discuss how the FRIC unveils emission processes that are hidden both in the emission spectrum and in the color-blind intensity correlation.

## II. THEORETICAL MODEL

In this section, we introduce the theoretical framework used to describe the dynamics of two interacting TLS-like quantum emitters under laser illumination and at cryogenic temperatures. We consider the emitters at fixed positions $r_i$, where $i \in \{1, 2\}$ labels each of them. The emitters are described as simple TLS with transition frequencies $\omega_{e,1}$ and $\omega_{e,2} = \omega_{e,1} - \delta$, where $\delta$ is the detuning between them. These transition frequencies contain the self-energy shifts induced by the vacuum field (see Refs. [43,44]). The ground and the excited states of emitter $i$ are denoted as $|g_i\rangle$ and $|e_i\rangle$, respectively, and the transition dipole moment as $\mu_i$. The Hamiltonian of the uncoupled emitters can be written as

$$ \hat{H}_0 = \hbar \frac{\omega_{e,1}}{2} \hat{\sigma}_{z,1} + \hbar \frac{\omega_{e,2}}{2} \hat{\sigma}_{z,2}, \quad (1) $$

where $\hat{\sigma}_{z,i} = |e_i\rangle\langle e_i| - |g_i\rangle\langle g_i|$ is the z-Pauli matrix in the Hilbert space of the emitter $i$. The corresponding energy scheme is shown in Fig. 1(b).

The coherent interaction of the quantum emitters through the vacuum field is incorporated via the interaction
Hamiltonian [45]

$$ \hat{H}_{\text{int}} = \hbar V (\hat{\sigma}_1^\dagger \hat{\sigma}_2 + \hat{\sigma}_1 \hat{\sigma}_2^\dagger), \quad (2) $$

here $\hat{\sigma}_i^\dagger = |e_i\rangle \langle g_i|$ and $\hat{\sigma}_i$ are the raising and lowering Pauli operators of the emitter $i$, respectively. The coupling strength between the two emitters is given by the dipole-dipole coupling as

$$ V = 3\alpha\gamma_0 \frac{(\mu_1 \cdot \mu_2) - 3(\mu_1 \cdot \mathbf{r}_{12})(\mu_2 \cdot \mathbf{r}_{12})/r_{12}^2}{4\mu_1\mu_2(kr_{12})^3}, \quad (3) $$

where the spontaneous decay rate $\gamma_0$ of both emitters is assumed to be identical, $\mathbf{r}_{12} = \mathbf{r}_2 - \mathbf{r}_1$ is the relative position vector between both emitters and $k = n\omega_0/c$ is the wavenumber corresponding to the mean frequency $\omega_0 = (\omega_{e,1} + \omega_{e,2})/2$, with $n$ the refraction index of the medium and $c$ the speed of light in vacuum. The dipole-dipole coupling in Eq. (3) is valid for molecules at short distances (in comparison to the wavelength of light), i.e., $kr_{12} \ll 1$ [44,45]. Figure 1(c) shows a schematic representation of the energy levels of the coupled emitters for no illumination. The energy difference between the ground state $|G\rangle = |g_1g_2\rangle$ and the doubly excited state $|E\rangle = |e_1e_2\rangle$ is given by the sum of the frequencies of the two excited states $2\omega_0$. The energy difference between the two single-excitation states $|\Lambda_+\rangle$ (higher energy) and $|\Lambda_-\rangle$ (lower energy) is given by $2\Lambda = \sqrt{4V^2 + \delta^2}$. These single-excitation states correspond to delocalized states that result from the dipole-dipole coupling and they can be written as $|\Lambda_+\rangle = \cos\theta|g_1e_2\rangle + \sin\theta|e_1g_2\rangle$ and $|\Lambda_-\rangle = -\sin\theta|g_1e_2\rangle + \cos\theta|e_1g_2\rangle$, where the angle $\theta$ fulfills $\sin 2\theta = V/\Lambda$ and $\cos 2\theta = \delta/(2\Lambda)$. These $|\Lambda_-\rangle$ and $|\Lambda_+\rangle$ delocalized states correspond to the superradiant and subradiant states (which one is the superradiant state depends on $\sin 2\theta$) and are discussed below.

Additionally, we consider that the two quantum emitters are illuminated by a laser of frequency $\omega_{las}$ and intensity $I$. The coherent illumination is introduced through the Hamiltonian

$$ \hat{H}_{\text{las}} = \hbar \sum_{i \in \{1,2\}} \frac{\mu_i \cdot \mathcal{E}}{2\hbar} (\hat{\sigma}_i e^{i\omega_{\text{las}}t} + \hat{\sigma}_i^\dagger e^{-i\omega_{\text{las}}t}), \quad (4) $$

where $\mathcal{E}$ is the amplitude of the electric field, which is assumed to be the same at the position of the two emitters. Equation (4) does not include any spatial phase because we consider that both emitters are located in the same plane and that the illumination impinges normally to this plane. Additionally, we assume that both emitters have identical transition dipole moments ($\mu_1 = \mu_2$) and that the electric field is linearly polarized in the same direction than the transition dipole moments. In this configuration, both Rabi frequencies become equal $\Omega = -\mu_1 \cdot \mathcal{E}/\hbar = -\mu_2 \cdot \mathcal{E}/\hbar$. The laser intensity $I$ and the Rabi frequency $\Omega$ are related as $I = \Omega^2 \hbar \omega_0^3 n^2 / (6\pi c^2 \gamma_0) = 2I_{\text{sat}} \Omega^2 / \gamma_0^2$ as we do not consider incoherent dephasing processes [46]. Here, $I_{\text{sat}} = \hbar \omega_0^3 \gamma_0 n^2 / (12\pi c^2)$ is the saturation intensity of a single TLS-like emitter with transition frequency $\omega_0$ and decay rate $\gamma_0$.

The total Hamiltonian $\hat{H} = \hat{H}_0 + \hat{H}_{\text{int}} + \hat{H}_{\text{las}}$ can be written in the interacting basis $\{|E\rangle, |\Lambda_+\rangle, |\Lambda_-\rangle, |G\rangle\}$ and in the

rotating frame of the laser as

$$ \begin{aligned} \hat{H} = \hbar \Delta_0 (|E\rangle \langle E| - |G\rangle \langle G|) + \hbar \Lambda (|\Lambda_+\rangle \langle \Lambda_+| - |\Lambda_-\rangle \langle \Lambda_-|) \\ - \hbar \sum_{k \in \{+,-\}} \frac{\Omega_k}{2} (\hat{\sigma}_{Gk} + \hat{\sigma}_{kE} + \hat{\sigma}_{kG} + \hat{\sigma}_{Ek}), \end{aligned} \quad (5) $$

where we have introduced the laser detuning $\Delta_0 = \omega_0 - \omega_{las}$ with respect to the mean frequency $\omega_0$ and the effective pumping rates $\Omega_\pm = \Omega(\cos\theta \pm \sin\theta)$. The pumping rate $\Omega_+$ corresponds to the driving strength of the transitions $|G\rangle \to |\Lambda_+\rangle$ and $|\Lambda_-\rangle \to |E\rangle$, whereas $\Omega_-$ to the driving strength of the transitions $|G\rangle \to |\Lambda_-\rangle$ and $|\Lambda_+\rangle \to |E\rangle$. We have also introduced the operator notation $\hat{\sigma}_{G\pm} = |G\rangle \langle \Lambda_\pm|$, $\hat{\sigma}_{E\pm} = |E\rangle \langle \Lambda_\pm|$, $\hat{\sigma}_{G\pm} = |\Lambda_\pm\rangle \langle G|$, and $\hat{\sigma}_{E\pm} = |\Lambda_\pm\rangle \langle E|$, for convenience. The terms in the first line of Eq. (5) represent the undriven Hamiltonian, including the dipole-dipole interaction between the two emitters, while the terms in the second line account for the driving of the system.

We incorporate the incoherent losses of the system using the standard Lindblad master equation [47]. The dynamics of the density matrix $\hat{\rho}$ that describes the state of the composite system is given by [45]

$$ \frac{d}{dt} \hat{\rho} = -\frac{i}{\hbar} [\hat{H}, \hat{\rho}] + \mathcal{L}_{\text{inc}} \hat{\rho}, \quad (6) $$

where the superoperator

$$ \mathcal{L}_{\text{inc}} = \frac{\gamma_0}{2} \left( \sum_{i \in \{1,2\}} \mathcal{D}[\hat{\sigma}_i] + \frac{\tilde{\gamma}}{2} (\mathcal{D}[\hat{\sigma}_1, \hat{\sigma}_2] + \mathcal{D}[\hat{\sigma}_2, \hat{\sigma}_1]) \right) \quad (7) $$

accounts for the incoherent dynamics, with $\mathcal{D}[\hat{a}, \hat{b}]\hat{\rho} = 2\hat{a}\hat{\rho}\hat{b}^\dagger - \{\hat{b}^\dagger\hat{a}, \hat{\rho}\}$ and $\mathcal{D}[\hat{a}]\hat{\rho} = \mathcal{D}[\hat{a}, \hat{a}]\hat{\rho}$ the crossed and the standard Lindblad dissipators, respectively. The first term (the summation) on the right-hand side of Eq. (7) corresponds to the spontaneous emission of photons of each emitter and the decay rate is given by $\gamma_0 = k^3|\mu_i|^2 / (3\pi\epsilon_0\hbar n^2)$ [45,48], with $\epsilon_0$ the vacuum permittivity. The second term (proportional to $\tilde{\gamma}$) stands for the dissipative coupling between the emitters induced by the vacuum field [45,48]. Assuming $kr_{12} \ll 1$, the strength of this crossed decay is given by $\tilde{\gamma} = \alpha\gamma_0\mu_1 \cdot \mu_2 / (\mu_1\mu_2)$ [19,45], which for parallel transition dipole moments becomes $\tilde{\gamma} = \alpha\gamma_0$. Equation (7) neglects the thermal population of the electronic excited states $|e_i\rangle$ and, additionally, does not contain dephasing Lindblad operators, as we focus on quantum emitters at cryogenic temperatures.

The decay rates $\gamma_\pm$ of the delocalized states $|\Lambda_+\rangle$ and $|\Lambda_-\rangle$ are obtained by writing the superoperator $\mathcal{L}_{\text{inc}}$ in the interacting basis (see Appendix A), which gives $\gamma_\pm = \gamma_0(1 \pm \alpha \sin 2\theta)$. Thus, the delocalized states have different decay rates if the interaction between the emitters is nonnegligible ($\sin 2\theta \neq 0$). If this condition is fulfilled, one of the delocalized states has a larger decay rate than the spontaneous decay rate of the noninteracting emitters $\gamma_0$, i.e., it couples more efficiently with light. This state is usually called superradiant state. In contrast, the other delocalized state has a smaller decay rate than $\gamma_0$ and is known as subradiant state. As the subradiant state has a longer coherence time, it could turn into a good candidate for applications in quantum information storage [49–51]. Notably, the combined Debye-Waller/Franck-Condon factor has a strong impact on
the decay rates of the superradiant and subradiant states, e.g., in the case of identical emitters ($\delta = 0$) the decay rates of the superradiant and subradiant states are no longer $2\gamma_0$ and 0, respectively, but they are $\gamma_0(1 \pm \alpha)$.

Hereafter, as an example of application of our system, we focus on two dibenzanthanthrene (DBATT) molecules immersed in a naphthalene film with $n = 1.5$ in a J-aggregate configuration, i.e., we assume that $r_{12}$ is parallel to the dipole moments. In this configuration, the lower-energy delocalized state $|\Lambda_{-}\rangle$ is the superradiant state and the higher-energy delocalized state $|\Lambda_{+}\rangle$ is the subradiant state. We also assume the following parameters, based on the experiments in Ref. [19]: (i) the mean frequency is $\omega_0/(2\pi) \simeq 485$ THz, corresponding to a wavelength of 618 nm; (ii) the spontaneous decay rate is $\gamma_0/(2\pi) = 21.5$ MHz; (iii) the combined Debye-Waller/Franck-Condon factor is $\alpha = 0.3$; and (iv) the dipole-dipole coupling is $V = -20\gamma_0$, i.e., $r_{12} \approx 18.5$ nm. We emphasize that the resulting intensity correlations are also valid for any other quantum emitters having the same coupling strength $V$ and combined Debye-Waller/Franck-Condon factor $\alpha$. Additionally, equivalent results are obtained for the case of a molecular H-aggregate configuration ($r_{12}$ is perpendicular to the transition dipole moments), with the only difference that in this case the superradiant state is the higher-energy delocalized state $|\Lambda_{+}\rangle$ and the subradiant state is the lower-energy delocalized state $|\Lambda_{-}\rangle$.

### III. COLOR-BLIND INTENSITY CORRELATION

In this section, we examine the color-blind intensity correlation, which considers all the photons emitted from the interacting molecules and is defined as [24,52,53]

$$g^{(2)}(\tau) = \frac{\langle \hat{E}^{(-)}(0) \hat{E}^{(-)}(\tau) \hat{E}^{(+)}(\tau) \hat{E}^{(+)}(0) \rangle_{ss}}{\langle E^{(-)}(0) E^{(+)}(0) \rangle_{ss} \langle E^{(-)}(\tau) E^{(+)}(\tau) \rangle_{ss}} \quad (8)$$

In this expression, $\langle \hat{O} \rangle_{ss} = \text{Tr}(\hat{O}\hat{\rho}_{ss})$ denotes the expected steady-state value (with $\hat{\rho}_{ss}$ the steady-state density matrix) and $\hat{E}^{(+)}(t)$ and $\hat{E}^{(-)}(t)$ represent the standard positive-frequency and negative-frequency electric field operators in the Heisenberg picture at time $t$, respectively. We consider that light is detected in a direction normal to the orientation of the transition dipole moments and, thus, $\hat{E}^{(+)}(t) \propto \hat{\sigma}_1(t) + \hat{\sigma}_2(t)$ [43,54,55]. Hereafter, we focus on the intensity correlation at zero delay $g^{(2)}(0)$. Let us remind that a coherent source of light has Poissonian statistics and fulfills $g^{(2)}(0) = 1$. When $g^{(2)}(0) < 1$ the statistics of light is sub-Poissonian, but following a usual convention we call it antibunched. In contrast, if $g^{(2)}(0) > 1$ the statistics of light is super-Poissonian and the emission is said to be bunched.

We plot in Figs. 2(a)–2(c) the dependence of $g^{(2)}(0)$ on the molecular detuning $\delta$ (which can be tuned experimentally via the Stark effect [19,20,22,23]) and the laser intensity $I$ for different laser frequencies. We find that $g^{(2)}(0)$ exhibits a broad range of values, ranging from almost perfect antibunching [$g^{(2)}(0) \approx 0$] to extreme bunching [$g^{(2)}(0) \gg 1$].

First, in Fig. 2(a) the laser frequency is chosen to enable the resonant excitation of the doubly excited state $|E\rangle$ through a two-photon process ($2\omega_{las} = 2\omega_0$). For this laser detuning and under weak illumination, $g^{(2)}(0)$ is strongly bunched, which reveals that the emission of photons in cascade is strongly

enhanced with respect to light emission obeying Poissonian statistics. Additionally, the difference between the laser frequency $\omega_{las} = \omega_0$ and the resonance frequencies $\omega_0 \pm \Lambda$ of the single-excitation states $|\Lambda_{\pm}\rangle$ increases with the molecular detuning. Thus, the ratio between two-photon emission processes and single-photon emission processes grows with the molecular detuning. As a consequence, $g^{(2)}(0)$ also increases with the molecular detuning $\delta$, as shown in Fig. 2(a). In Appendix B, we provide a more extensive analysis of the light statistics for this laser detuning, supported by analytical expressions valid for arbitrary laser intensity.

Furthermore, Fig. 2(b) indicates that the emission is strongly antibunched when the laser is tuned to the transition frequency of the superradiant state $|\Lambda_{-}\rangle$ and under weak illumination. Consequently, the two coupled molecules act as a single-photon source for this laser frequency, as occurs for a single molecule, but with a larger decay rate $\gamma_0(1 + \alpha)$ and larger transition dipole moment (i.e., the emission of single photons is accelerated by a factor $1 + \alpha$ with respect to the case of a single emitter). In Appendix C, we analytically quantify the small deviation of the light statistics from that of an ideal single photon source [with $g^{(2)}(0) = 0$] for this laser detuning, as well as the increase of $g^{(2)}(0)$ observed in Fig. 2(b) for larger laser intensities.

Moreover, a more complex dependence of $g^{(2)}(0)$ on the laser intensity and on the molecular detuning is found for a laser tuned to the transition frequency of the subradiant state $|\Lambda_{+}\rangle$ in Fig. 2(c). In this case, when the two molecules are very far of resonance ($|\delta| \gg 30\gamma_0$), the system behaves again similar to a single TLS and the emission is strongly antibunched under weak illumination. In the opposite case, if the molecules are identical ($\delta = 0$) the light emitted can have approximately Poissonian statistics [$g^{(2)}(0) \approx 1$]. Interestingly, if the molecules are slightly off-resonance, then both bunched and antibunched emission can be obtained for resonant excitation of the subradiant state, depending on the laser intensity $I$, as discussed below.

We further note that, for any laser frequency and molecular detuning, the intensity correlation converges to $g^{(2)}(0) = 1$ for strong enough laser intensity because the molecules saturate and thus become uncorrelated [56–59]. Additionally, we investigate in Appendix D the impact of the combined Debye-Waller/Franck-Condon factor $\alpha$ on the intensity correlation $g^{(2)}(0)$ for the three different laser frequencies analyzed in Fig. 2. We find that $g^{(2)}(0)$ is particularly sensitive to $\alpha$ when the laser is tuned to the transition frequency of the subradiant state.

To summarize the broad range of color-blind intensity correlations that can be obtained by tuning the laser frequency and intensity, we fix $\delta = 10\gamma_0$ in Fig. 2(d) and plot the dependence of $g^{(2)}(0)$ on the laser intensity for a laser tuned to (red line) the two-photon resonance, (blue line) the transition frequency of the superradiant state and (green line) the transition frequency of the subradiant state. We have checked that the intensity correlations shown in Fig. 2 are robust against misalignment in the laser polarization and in the orientation of the transition dipole moments.

Furthermore, Fig. 2 indicates different ways to tune the emission from antibunching to bunching. For example, for two slightly detuned emitters and a laser tuned to the transition
<figure><img src="image_4.png" /><figcaption>FIG. 2. Dependence of the color-blind intensity correlation $g^{(2)}(0)$ on the laser intensity $I$ and molecular detuning $\delta$. In (a) the laser is tuned to the two-photon resonance ($2\omega_{las} = 2\omega_0$), which leads to bunched emission. In (b) the laser excites resonantly the superradiant state $|\Lambda_{-}\rangle$ ($\omega_{las} = \omega_0 - \Lambda$) giving raise to antibunched emission. In (c) the laser is tuned to the transition frequency of the subradiant state $|\Lambda_{+}\rangle$ ($\omega_{las} = \omega_0 + \Lambda$) and $g^{(2)}(0)$ shows a more complex dependence on the laser intensity and molecular detuning. In (a)–(c), color blue represents antibunched emission [$g^{(2)}(0) < 1$], red bunched emission [$g^{(2)}(0) > 1$], and white Poissonian statistics [$g^{(2)}(0) = 1$]. Cuts of the results in (a)–(c) for $\delta = 10\gamma_0$ (along the dashed-black lines) are plotted in (d), where the red line corresponds to a laser tuned to the two-photon resonance, the blue line is obtained for a laser tuned to the superradiant state and the green line corresponds to a laser tuned to the subradiant state. In (d), red and blue shaded background correspond to bunched and antibunched statistics, respectively. The insets illustrate the different laser detunings considered. The dipole-dipole coupling is fixed as $V = -20\gamma_0$ and the combined Debye-Waller/Franck-Condon factor as $\alpha = 0.3$.</figcaption></figure>

frequency of the subradiant state, we can obtain both types of light statistics by modifying the laser intensity. Additionally, we can use weak illumination and tune the laser either to the two-photon resonance or to the superradiant state in order to obtain extreme values of bunching and antibunching, respectively. This level of controlled crossover between bunched and antibunched emission depending on a reasonable variation of experimental parameters can be of potential technological interest in engineering photon sources. Alternative physical configurations have been already proposed to achieve this crossover, such as tuning the relative phase between the light scattered coherently and incoherently by an ensemble of

trapped ions in a nanofiber [60], manipulating the position of the detector for the case of two trapped ions [61], controlling the temperature of a thermal atomic vapor in a cell [62] or tailoring the laser polarization and the material gain of a plasmonic nanosphere close to two quantum emitters [63].

### A. Analytical insights

A better understanding of the intensity correlations in Fig. 2 can be obtained from the steady-state populations and coherences of the interacting system. Substituting $\hat{E}^{(+)}(0) \propto \hat{\sigma}_1 + \hat{\sigma}_2$ into Eq. (8) yields

$$g^{(2)}(0) = \frac{4\rho_E}{[2\rho_E + (1 + \sin 2\theta)\rho_{\Lambda_+} + (1 - \sin 2\theta)\rho_{\Lambda_-} - \cos 2\theta \text{Re}\rho_{\Lambda_+}]^2}, \quad (9)$$
<figure><img src="image_5.png" /><figcaption>FIG. 3. Color-blind intensity correlation and steady-state populations and coherences when the laser is tuned resonantly to the subradiant state |Λ+⟩. (a) Dependence on the laser intensity of the color-blind intensity correlation g^{(2)}(0). The solid-black line is obtained numerically using the exact expression in Eq. (9), whereas the dashed-blue line is calculated numerically using the approximate Eq. (10). The inset illustrates the laser frequency used. (b) Dependence on I of the steady-state populations of (blue) the superradiant state |Λ−⟩, (green) the subradiant state |Λ+⟩, (red) the doubly excited state |E⟩, and (pink) the absolute value of the real part of the steady-state coherence ρ+− = ⟨Λ+|ρ̂ss|Λ−⟩. Solid lines correspond to the exact numerical calculations, whereas dots represent the analytical results obtained using Eqs. (E2), (E5), and (E7) derived in Appendix E. Black-dashed lines represent the linear and quadratic scaling of ρ_E with laser intensity. (c) Dependence on the laser intensity of the coherences (dashed-orange line) −Ω±/2γ₀ Imρ_E- and (solid-purple line) Ω±/2γ₀ Imρ_E+, whose difference corresponds to (red-solid line) the population of the doubly excited state according to Eq. (12). The inset corresponds to a zoom (within the interval from I/I_{sat} = 7 × 10⁻² to I/I_{sat} = 10⁻¹) showing that −Ω±/2γ₀ Imρ_E- is slightly larger than Ω±/2γ₀ Imρ_E+. We have fixed the values V = −20γ₀, α = 0.3, and δ = 10γ₀ in all panels.</figcaption></figure>

where we have introduced the steady-state populations $\rho_k = \langle k | \hat{\rho}_{ss} | k \rangle$, with $|k\rangle$ a state of the interacting basis $\{|E\rangle, |\Lambda_+\rangle, |\Lambda_-\rangle, |G\rangle\}$, and the coherence $\rho_{+−} = \langle \Lambda_+ | \hat{\rho}_{ss} | \Lambda_- \rangle$.

For weak laser intensities ($I \lesssim I_{sat}$) and for the three laser detunings considered in this paper, we find numerically that the population $\rho_E$ of the doubly excited state $|E\rangle$ and the coherence $\rho_{+−}$ are much smaller than the populations $\rho_{\Lambda_{\pm}}$ of the delocalized states $|\Lambda_{\pm}\rangle$. Thus, the intensity correlation can be approximated under weak illumination as

$$g^{(2)}(0)\Big|_{I \lesssim I_{\text{sat}}} = \frac{4\rho_E}{[(1 + \sin 2\theta)\rho_{\Lambda_+} + (1 - \sin 2\theta)\rho_{\Lambda_-}]^2}. \quad (10)$$

The numerator in Eq. (10) accounts for the probability of emitting two photons, which is proportional to the steady-state population of the doubly excited state. Similarly, the denominator is related the probability of single-photon emission and is given by the steady-state populations of the single-excitation delocalized states. Equation (10) indicates that if the population of the doubly excited state is much smaller than the square of the steady-state populations of the delocalized states, then the intensity correlation becomes strongly antibunched. In the opposite limit, strong bunching occurs.

In the following, we discuss in more detail the light emission processes behind the complex behavior of the intensity correlation when the laser is tuned to the transition frequency of the subradiant state $|\Lambda_+\rangle$ ($\omega_{las} = \omega_0 + \Lambda$) and the molecular detuning is fixed at $\delta = 10\gamma_0$. For this value of molecular detuning, $\sin 2\theta = -0.97$ and $\cos 2\theta = 0.24$ in Eqs. (9) and (10). The intensity correlation for the other two laser detunings at the same molecular detuning is analyzed in detail in Appendices B and C, where we provide analytical expressions of g^{(2)}(0) that are valid for arbitrary laser intensity.

## B. Laser tuned to the subradiant state

When the laser is tuned to the transition frequency of the subradiant state there is an important competition between two excitation mechanisms: (i) the resonant excitation of the subradiant state $|\Lambda_+\rangle$, a state that couples very inefficiently to light; and (ii) the nonresonant excitation of the superradiant state $|\Lambda_-\rangle$, which couples very well to light. As a consequence, we obtain a complex dependence of the intensity correlation g^{(2)}(0) on the laser intensity, as shown in Fig. 3(a) [corresponding to the green line in Fig. 2(d)].

The light emitted is antibunched for weak laser intensities ($I \lesssim I_{sat}$), but g^{(2)}(0) remains larger than 0.25, much larger than the minimum value of g^{(2)}(0) when a laser of the same intensity is tuned to the transition frequency of the superradiant state $|\Lambda_-\rangle$, see Fig. 2(d). This larger value indicates that the coupled molecules cannot be approximated as a single TLS for this laser frequency and molecular detuning ($\delta = 10\gamma_0$).

Moreover, for weak enough laser intensities ($I \ll I_{sat}$), the steady-state populations in Fig. 3(b) can be obtained from the populations of two independent three-level systems, and g^{(2)}(0) becomes independent of the laser intensity (see Appendix E for the detailed derivation of the subsequent equations). One of these three-level systems consists of the states $|G\rangle, |\Lambda_+\rangle, |E\rangle$ and the other one of $|G\rangle, |\Lambda_-\rangle, |E\rangle$, which accounts for the two competing processes to excite the doubly excited state $|E\rangle$. In this manner, the populations $\rho_{\Lambda_+}$ and $\rho_{\Lambda_-}$ of the single-excitation delocalized states $|\Lambda_+\rangle$ and $|\Lambda_-\rangle$ scale linearly with the laser intensity as $\rho_{\Lambda_+}|_{I \ll I_{sat}}^{\Lambda_+} = \frac{1}{I_{sat}}\frac{\gamma_0^2(1+\sin 2\theta)}{\gamma_+^2}$ and $\rho_{\Lambda_-}|_{I \ll I_{sat}}^{\Lambda_-} = \frac{1}{I_{sat}}\frac{1-\sin 2\theta}{8\eta^2}$, where $\eta = (2\Lambda)/\gamma_0$ and we have assumed $\Lambda \gg \gamma_0$ to further simplify these expressions. The population of the subradiant state $|\Lambda_+\rangle$ [green line in Fig. 3(b)] is two orders of magnitude larger than that of the superradiant state $|\Lambda_-\rangle$ (blue line), since the former is driven resonantly and the later nonresonantly (with laser detuning equal to $2\Lambda$). On the other hand, in the denominator of
$g^{(2)}(0)$ in Eq. (10) the population $\rho_{\Lambda_{-}}$ of the superradiant state $|\Lambda_{-}\rangle$ is multiplied by a factor $1 - \sin 2\theta \approx 1.97$, which is two orders of magnitude larger than the factor $1 + \sin 2\theta \approx 0.03$ multiplying the population $\rho_{\Lambda_{+}}$ of the subradiant state $|\Lambda_{+}\rangle$, since the superradiant state couples much more efficiently to light. Both delocalized states thus contribute comparably to the denominator of Eq. (10) and need to be taken into account in the calculation of $g^{(2)}(0)$.

Further, for these laser intensities ($I \ll I_{\text{sat}}$), the doubly excited state $|E\rangle$ is populated through two independent two-photon processes (one via the subradiant state and another one via the superradiant state) and thus its population $\rho_E$ [red line in Fig. 3(b)] scales quadratically with the laser intensity as $\rho_E|_{I \ll I_{\text{sat}}} = \frac{I^2}{I_{\text{sat}}^2} \left[ \frac{(1+\sin 2\theta)\gamma_0^2}{8\Lambda\gamma_+} \right]^2 + \frac{I^2}{I_{\text{sat}}^2} \left[ \frac{1-\sin 2\theta}{8\eta^2} \right]^2$, where the first term of the right-hand side accounts for the population excited via the subradiant state $|\Lambda_{+}\rangle$ and the second term via the superradiant state $|\Lambda_{-}\rangle$. According to Eq. (10), $g^{(2)}(0)$ therefore becomes independent of the laser intensity under weak illumination and can be obtained as

$$g^{(2)}(0)|_{I \ll I_{\text{sat}}} = 4 \frac{\left(\frac{1-\sin 2\theta}{(4\Lambda)^2}\right)^2 + \left(\frac{1+\sin 2\theta}{4\Lambda\gamma_+}\right)^2}{\left[\left(\frac{1-\sin 2\theta}{4\Lambda}\right)^2 + \left(\frac{1+\sin 2\theta}{\gamma_+}\right)^2\right]^2}. \quad (11)$$

However, for this laser detuning the intensity correlation has a larger value [$g^{(2)}(0)|_{I \ll I_{\text{sat}}} \approx 0.25$] than when the laser is tuned to the transition frequency of the superradiant state [$g^{(2)}(0)|_{I \ll I_{\text{sat}}} \approx 2.5 \times 10^{-4}$], which emphasizes that the dynamics of the coupled molecules is now more complex than that of a single TLS-like emitter.

Additionally, Eq. (11) can be used to estimate molecular parameters in experiments. For example, the distance between the two emitters $r_{12}$ and the molecular detuning $\delta$ could be simultaneously estimated through the measurement of $g^{(2)}(0)|_{I \ll I_{\text{sat}}}^{|\Lambda_{+}|}$ and the splitting $2\Lambda = \sqrt{4V^2 + \delta^2}$ between the transition frequencies of the superradiant and subradiant states. Alternatives to extract $r_{12}$ are based on superresolution imaging techniques [19,27], which are time demanding.

The value of the intensity correlation $g^{(2)}(0)$ increases for stronger laser intensities ($10^{-1} \lesssim I/I_{\text{sat}} \lesssim 10^1$). In this range of laser intensity, the slope of the population $\rho_E$ of the doubly excited state $|E\rangle$ scales superquadratically, which indicates that the excitation of the doubly excited state $|E\rangle$ can no longer be understood as due to two independent two-photon processes (i.e., its population $\rho_E$ does not correspond to the population given by two independent three-level systems). In contrast, the excitation of $|E\rangle$ via the resonantly driven subradiant state $|\Lambda_{+}\rangle$ and via the nonresonantly driven superradiant state $|\Lambda_{-}\rangle$ coherently interfere. This effect is revealed, on the one hand, in the increase in the slope of the real part of the coherence $\rho_{+-} = \langle \Lambda_{+} | \hat{\rho}_{ss} | \Lambda_{-} \rangle$ between the superradiant and subradiant states [represented with a pink-solid line in Fig. 3(b)]. On the other hand, the interference between the excitation of the doubly excited state from the superradiant state and from the subradiant state also has an effect in the complex interplay between the imaginary part of the coherence $\rho_{E+} = \langle E | \hat{\rho}_{ss} | \Lambda_{+} \rangle$ (related to the excitation of the doubly excited state via the subradiant state $|\Lambda_{+}\rangle$) and the imaginary part of the coherence $\rho_{E-} = \langle E | \hat{\rho}_{ss} | \Lambda_{-} \rangle$ (related to the excitation via the superradiant state $|\Lambda_{-}\rangle$). These latter

coherences contribute to the population of the doubly excited state as

$$\rho_E = -\frac{\Omega_+}{2\gamma_0} \text{Im}\rho_{E+} - \frac{\Omega_-}{2\gamma_0} \text{Im}\rho_{E-}. \quad (12)$$

Figure 3(c) shows that the two terms on the right-hand side of Eq. (12) are of opposite sign and are comparable in magnitude. Under very weak illumination, both terms scale quadratically with the intensity and thus so does $\rho_E$, as the negative term is slightly smaller in magnitude than the positive one [see inset in Fig. 3(c)]. For increasing laser intensities the slope of the negative term starts to decrease and thus its cancellation with the larger positive term is reduced, which leads to the superquadratic behavior of $\rho_E$ with the laser intensity. Additionally, the populations of the single-excitation delocalized states $|\Lambda_{+}\rangle$ and $|\Lambda_{-}\rangle$ still depend linearly on the laser intensity [see Fig. 3(b)] and thus the superquadratic scaling of the doubly excited state leads to an increase in the intensity correlation [see Fig. 3(a)], according to Eq. (10). Importantly, $g^{(2)}(0)$ eventually crosses from antibunching to bunching, which indicates that the correlation can be of either one type or the other for this laser detuning and molecular detuning ($\delta = 10\gamma_0$), depending on the intensity of the laser. This crossover is also captured by Eq. (10) [dashed-blue line in Fig. 3(a)], which shows a good agreement with the results obtained with the exact Eq. (9) for $I \lesssim I_{\text{sat}}$.

For laser intensities larger than $10I_{\text{sat}}$, the population $\rho_{\Lambda_{+}}$ of the subradiant state $|\Lambda_{+}\rangle$ reaches a maximum value of 0.33. Additionally, the population of the doubly excited state $|E\rangle$ becomes linearly dependent on the laser intensity due to the large population of the subradiant state for these laser intensities. The intensity correlation then reaches a maximum of value 3.71. For even larger laser intensities ($I \gtrsim 10^3 I_{\text{sat}}$), all populations converge to 0.25 (i.e., both molecules become saturated and thus uncorrelated) and the intensity correlation converges to $g^{(2)}(0) = 1$. Notably, the evolution of $g^{(2)}(0)$ from its maximum value to 1 is not a simple decay, but it reaches a minimum before growing again, which highlights the complexity of the excitation and emission processes between states for this laser detuning.

## IV. FREQUENCY-RESOLVED INTENSITY CORRELATION

The color-blind intensity correlation $g^{(2)}(\tau)$ analyzed in the previous section measures the correlation between all the photons emitted from the molecules, without any discrimination on their frequency. In contrast, the frequency-resolved intensity correlation (FRIC) $g^{(2)}(\omega_1, \omega_2; \tau)$ addresses the correlation between pairs of photons emitted at particular frequencies $\omega_1$ and $\omega_2$. $g^{(2)}(\omega_1, \omega_2; \tau)$ provides information about the different emission processes of the molecules and can unveil processes that are not revealed in the emission spectrum. For example, it is known that FRIC can unveil two-photon transitions through virtual states [34,38–40]. Additionally, we show in this section that FRIC can also unveil one-photon transitions that are not revealed in the one-photon spectrum.

The FRIC of light emitted from two interacting emitters was investigated in Ref. [34] for the case of incoherent driving. We consider here coherent illumination instead, and
show that the FRIC maps exhibit significant differences depending on the laser frequency. In particular, we investigate the FRIC of the same two-molecule system as in the previous section (i.e., with molecular detuning $\delta = 10\gamma_0$, dipole-dipole coupling $V = -20\gamma_0$ and combined Debye-Waller/Franck-Condon factor $\alpha = 0.3$). We analyze the cases in which the laser is tuned to the transition frequency of the superradiant state and to the transition frequency of the subradiant state. The FRIC for a laser tuned to the two-photon resonance is discussed in Appendix H and the results are similar to those reported in Ref. [41].

### A. Methodology

The FRIC can be measured experimentally by incorporating two Fabry-Pérot filters centered at frequencies $\omega_1$ and $\omega_2$ in the Hanbury-Brown Twiss interferometer [see Fig. 1(a)], with each filter placed in the path of one of the light beams that emerge from the beam splitter towards the detectors. In this paper, Lorentzian filters are assumed and $g^{(2)}(\omega_1, \omega_2; \tau)$ is calculated theoretically using the positive-frequency and negative-frequency electric field operators filtered at frequency $\omega$ [64], $\hat{E}_{\omega}^{(\pm)}(t) = (\Gamma/2) \int_{0}^{\infty} e^{-(i\omega+\Gamma/2)t'} \hat{E}^{(\pm)}(t-t')dt'$, where $\Gamma$ represents the linewidth of the filters. In terms of these operators, the FRIC is obtained as

$$g^{(2)}(\omega_1, \omega_2; \tau) = \frac{\langle T: [\hat{E}_{\omega_1}^{(-)}(0) \hat{E}_{\omega_2}^{(-)}(\tau) \hat{E}_{\omega_2}^{(+)}(\tau) \hat{E}_{\omega_1}^{(+)}(0)] \rangle_{\text{ss}}}{\langle E_{\omega_1}^{(-)}(0) E_{\omega_1}^{(+)}(0) \rangle_{\text{ss}} \langle E_{\omega_2}^{(-)}(\tau) E_{\omega_2}^{(+)}(\tau) \rangle_{\text{ss}}}, \quad (13)$$

where $T$ denotes the product of the time-ordered and normal-ordered operators.

The direct calculation of the FRIC through Eq. (13) involves complicated time integrals [35,36,64–67]. However, this computation is simplified by coupling weakly two TLSs to the system of interest [36,37,68]. These additional TLSs act as the two Fabry-Pérot filters in the experiment (at frequencies $\omega_1$ and $\omega_2$) and the corresponding Hamiltonian is $\hat{H}^{fl} = \sum_{i=1}^{2} \hbar\omega_i \hat{\zeta}_i^\dagger \hat{\zeta}_i$, where $fl$ stands for filter, $\hat{\zeta}_i^\dagger$ and $\hat{\zeta}_i$ represent the raising and lowering operators of filter $i$, and $\omega_i$ the corresponding central frequency. The interaction of the filters with the system is modeled by the Hamiltonian $\hat{H}_{\text{int}}^{fl} = \sum_{i=1}^{2} \hbar\varepsilon[\hat{\zeta}_i(\hat{\sigma}_1^\dagger + \hat{\sigma}_2^\dagger) + \text{H.c.}]$, with the coupling strength $\varepsilon$ weak enough not to disturb the dynamics of the system. The linewidth of the filters is incorporated via the additional Lindblad operators $\mathcal{L}_{fl} \hat{\rho} = \frac{\Gamma}{2} \mathcal{D}[\hat{\zeta}_1] \hat{\rho} + \frac{\Gamma}{2} \mathcal{D}[\hat{\zeta}_2] \hat{\rho}$. FRIC can then be calculated with the normalized two-time correlations of the steady-state populations of the two filters [36], i.e.,

$$g^{(2)}(\omega_1, \omega_2; \tau) = \lim_{\varepsilon \to 0} \frac{\langle \hat{\zeta}_1^\dagger \hat{\zeta}_1(0) \hat{\zeta}_2^\dagger \hat{\zeta}_2(\tau) \rangle_{\text{ss}}}{\langle \hat{\zeta}_1^\dagger \hat{\zeta}_1(0) \rangle_{\text{ss}} \langle \hat{\zeta}_2^\dagger \hat{\zeta}_2(\tau) \rangle_{\text{ss}}}. \quad (14)$$

We focus here on the FRIC at $\tau = 0$ and follow the numerical methodology proposed in Ref. [37] to calculate Eq. (14). This methodology is based on perturbation theory, valid for small values of $\varepsilon$, and leads to expressions that ultimately depend only on the dynamics of the system of interest, i.e., the Hilbert space does not increase in the numerical calculations due to the incorporation of the filters. Additionally, within

this methodology, the results do not depend on the particular choice of the coupling strength $\varepsilon$.

### B. Laser tuned to the superradiant state

We analyze first the FRIC when the laser frequency corresponds to the transition frequency of the superradiant state $|\Lambda_-\rangle$ ($\omega_{\text{las}} = \omega_0 - \Lambda$), for weak, moderate, and strong laser intensities.

#### 1. Weak laser intensity

We fix $I = 0.1 I_{\text{sat}}$ and show the normalized emission spectrum (top) and the FRIC map (bottom) in Fig. 4(a). In the FRIC map, the $x$ and $y$ axis correspond to the normalized laser detuning of the photons arriving at detector 1, ($\omega_1 - \omega_{\text{las}}$)/$\Lambda$, and at detector 2, ($\omega_2 - \omega_{\text{las}}$)/$\Lambda$, respectively. We use the standard color convention [34,38–41], where the blue color represents antibunched emission [$g^{(2)}(\omega_1, \omega_2; 0) < 1$], red bunched emission [$g^{(2)}(\omega_1, \omega_2; 0) > 1$], and white Poissonian statistics [$g^{(2)}(\omega_1, \omega_2; 0) = 1$]. To improve the visibility of the antibunched emission as well as of the bunched emission, we additionally make use of a linear scale in the range $0 \le g^{(2)}(\omega_1, \omega_2; 0) \le 1$ and a logarithmic scale in the range $1 \le g^{(2)}(\omega_1, \omega_2; 0) \le 100$. Regarding the emission spectrum $S(\omega) \propto \int_{-\infty}^{\infty} d\tau \langle \delta\hat{E}^{(-)}(\tau) \delta\hat{E}^{(+)}(0) \rangle_{\text{ss}} e^{-i\omega t}$, we calculate it using the fluctuation operators [e.g., $\delta\hat{E}^{(-)}(\tau) = \hat{E}^{(-)}(\tau) - (\hat{E}^{(-)}(\tau))_{\text{ss}}$] to avoid the Dirac delta peak at $\omega = \omega_{\text{las}}$ and we normalize it in such a way that the integration over the complete range of frequencies gives the intensity of the incoherently emitted light, i.e., $\int_{0}^{\infty} d\omega S(\omega) = \gamma_0 \langle \delta\hat{\sigma}_1^\dagger \delta\hat{\sigma}_1 + \delta\hat{\sigma}_2^\dagger \delta\hat{\sigma}_2 \rangle_{\text{ss}} + \tilde{\gamma} \langle \delta\hat{\sigma}_1^\dagger \delta\hat{\sigma}_2 + \delta\hat{\sigma}_2^\dagger \delta\hat{\sigma}_1 \rangle_{\text{ss}}$ [48].

The FRIC map in Fig. 4(a) is characterized by a strongly antibunched background, which can be understood from the very small value of the color-blind intensity correlation for these laser parameters [$g^{(2)}(0) \approx 10^{-3}$]. On top of this strongly antibunched background we observe a vertical line at $\omega_1 = \omega_{\text{las}}$, a horizontal line at $\omega_2 = \omega_{\text{las}}$ (marked by solid-gray arrows) and a bunched antidiagonal line (dashed-gray arrow) at $\omega_1 + \omega_2 = 2\omega_{\text{las}}$. These lines constitute a signature of the mechanisms of photon emission and to analyze them is convenient to account for the quantum nature of the illumination [34]. The eigenstates of the full quantized system are distributed in infinite rungs, with the eigenstates that belongs to the same rung having equal number of total excitations (i.e., molecular excitations plus photons in the laser field). For example, the eigenstates in the rung with $n$ excitations (rung $n$) under weak illumination are the hybrid states $\{|G, n\rangle, |\Lambda_-, n-1\rangle, |\Lambda_+, n-1\rangle, |E, n-2\rangle\}$, with $|n\rangle$ the $n$-Fock state of the laser field. Notably, for this laser detuning and intensity, the states $|\Lambda_-, n-1\rangle$ and $|G, n\rangle$ become degenerate in energy. In the same way, the states $|E, n-2\rangle$ and $|\Lambda_+, n-1\rangle$ are also degenerate in energy. Thus, for these laser parameters, each rung is composed by two doubly degenerate states. The eigenstates of three subsequent rungs are depicted in the level scheme in Fig. 4(b).

On the one hand, the frequencies at which the horizontal and vertical lines emerge in the FRIC map correspond to the frequencies of the one-photon transitions between eigenstates of the system, e.g., from the superradiant state $|\Lambda_-, n\rangle$ in
<figure><img src="image_6.png" /><figcaption>FIG. 4. FRIC for a laser tuned to the transition frequency of the superradiant state $|\Lambda_{-}\rangle$, $\omega_{las} = \omega_{0} - \Lambda$. [(a),(c),(e)] Normalized emission spectra (top panel) and FRIC maps (bottom panel) for laser intensity (a) $I = 0.1I_{sat}$, (c) $I = 50I_{sat}$, and (e) $I = 1000I_{sat}$. The inset in (a) illustrates the molecular energy levels $|G\rangle$ and $|\Lambda_{-}\rangle$ and the resonant excitation of the superradiant state. In the FRIC maps we use the standard color convention, where blue color represents antibunched emission, while red color bunched emission and white Poissonian statistics. In the color bar, we use a linear scale in the interval $0 \le g^{(2)}(\omega_1, \omega_2; 0) \le 1$ and a logarithmic scale in the range $1 \le g^{(2)}(\omega_1, \omega_2; 0) \le 100$. We consider the filter linewidth $\Gamma = \gamma_0/10$. [(b),(d),(f)] Schematic representation of the eigenstates (accounting for the quantum nature of the laser field) in runs $n-1, n$, and $n+1$ for laser intensity (b) $I = 0.1I_{sat}$, (d) $I = 50I_{sat}$, and (f) $I = 1000I_{sat}$. Each run contains four eigenstates with equal number of excitations. Solid arrows are used to mark the horizontal and vertical lines in the FRIC maps and the one-photon transitions in the diagrams of the eigenstates, whereas dashed arrows are used for the antidiagonal lines and the two-photon processes through virtual states. The dipole-dipole coupling is fixed as $V = -20\gamma_0$, the combined Debye-Waller/Franck-Condon factor as $\alpha = 0.3$ and the molecular detuning as $\delta = 10\gamma_0$.</figcaption></figure>

the runn $n + 1$ to the ground state $|G, n\rangle$ in the runn $n$. We use solid arrows in Figs. 4(a) and 4(b) to depict these one-photon transitions in the level scheme as well as to mark the corresponding horizontal and vertical lines (the same arrow convention is used in the rest of the FRIC maps). The emission of one photon of energy $\omega_{las}$ (here, $\omega_{las} = \omega_0 - \Lambda$) due to these transitions does not influence the frequency of the next emitted photon. As a consequence, the correlation $g^{(2)}(\omega_1 = \omega_{las}, \omega_2 = \omega'; 0)$ and $g^{(2)}(\omega_1 = \omega', \omega_2 = \omega_{las}; 0)$ is mostly independent of the frequency $\omega'$ of the second photon, which gives raise to vertical and horizontal lines at $\omega_{las} = \omega_0 - \Lambda$ of almost constant value (except at the positions that cross the antidiagonal lines discussed below). These one-photon transitions can also be identified by the peaks in the emission spectrum of the system, shown on top of Fig. 4(a).

On the other hand, the antidiagonal line in the FRIC map at $\omega_1 + \omega_2 = 2\omega_{las}$ is a signature of relaxation processes known as leapfrog [34]. These leapfrog processes consist in the emission of a photon of frequency $\omega_1$ ($\omega_2$) due to the transition from an initial eigenstate in the runn $n + 1$ to a virtual state and the emission of a photon of frequency $\omega_2$ ($\omega_1$) due to the

relaxation from the virtual state to another eigenstate in the runn $n - 1$. The energy difference between the initial and final eigenstates determines the value of $\omega_1 + \omega_2$ satisfied along the antidiagonal line. This antidiagonal line as well as the corresponding two-photon transition are marked with dashed arrows in Figs. 4(a) and 4(b) and hereafter. Additionally, as these antidiagonal lines correspond to the correlated emission of two photons, they are usually characterized by strong bunching. Importantly, these two-photon processes cannot be unveiled in the emission spectrum, and neither in the color-blind intensity correlation, which stresses the advantages of the FRIC to gain additional information about the transitions in the system.

Moreover, the FRIC map in Fig. 4(a) resembles that of a TLS under weak illumination [34,69], which we review in Appendix F. This similarity confirms that for these laser parameters the system behaves almost as a TLS. The impact of the subradiant state and the doubly excited state on the FRIC is thus very weak for this laser intensity and detuning. For example, the decay from the weakly populated subradiant state $|\Lambda_+, n\rangle$ to the ground state $|G, n\rangle$ is reflected in
the emission spectrum with the emergence of a very small peak at frequency $\omega_{las} + 2\Lambda = \omega_0 + \Lambda$ (solid-pink arrow). The corresponding lines are not appreciated in the FRIC map in Fig. 4(a) because the associated increase of $g^{(2)}(\omega_1, \omega_2; 0)$ is too small.

## 2. Moderate laser intensity

Figure 4(c) shows the FRIC map for $I = 50I_{sat}$. This laser intensity is strong enough to modify the eigenstates of the system, which are now referred to as dressed eigenstates [46,53]. For example, the rung $n$ consists of four dressed eigenstates $|E_i^{(n)}\rangle$ (with $i \in \{1, 2, 3, 4\}$), which can be written as different linear combinations of the bare states $\{|G, n\rangle, |\Lambda_-, n-1\rangle, |\Lambda_+, n-1\rangle, |E, n-2\rangle\}$. In particular, for the laser parameters used here, the two lowest-energy dressed eigenstates $|E_1^{(n)}\rangle$ and $|E_2^{(n)}\rangle$ in the rung $n$ can be accurately described as antisymmetric $|A_{G-}^{(n)}\rangle$ and symmetric $|S_{G-}^{(n)}\rangle$ combinations of the states $|G, n\rangle$ and $|\Lambda_-, n-1\rangle$ (which are degenerate in energy under weak illumination), i.e., $|E_1^{(n)}\rangle \approx |A_{G-}^{(n)}\rangle = (|G, n\rangle - |\Lambda_-, n-1\rangle)/\sqrt{2}$ and $|E_2^{(n)}\rangle \approx |S_{G-}^{(n)}\rangle = (|G, n\rangle + |\Lambda_-, n-1\rangle)/\sqrt{2}$. Similarly, $|E_3^{(n)}\rangle \approx |S_{E+}^{(n)}\rangle = (|E, n-2\rangle + |\Lambda_+, n-1\rangle)/\sqrt{2}$ and $|E_4^{(n)}\rangle \approx |A_{E+}^{(n)}\rangle = (|E, n-2\rangle - |\Lambda_+, n-1\rangle)/\sqrt{2}$ are symmetric and antisymmetric combinations of $|E, n-2\rangle$ and $|\Lambda_+, n-1\rangle$. We remark that the photon field occupation is assumed large enough so that the energy splitting between the different eigenstates within the same rung are identical for the rungs $n-1, n$, and $n+1$ in Fig. 4(d).

The dressed states $|S_{G-}^{(n)}\rangle$ and $|A_{G-}^{(n)}\rangle$ have a significant energy splitting equal to $|(\cos\theta - \sin\theta)\Omega| \approx 0.11\Lambda$ (see Appendix G), in contrast to the case of weak illumination where they are degenerate in energy. The one-photon transitions $|S_{G-}^{(n)}\rangle \rightarrow |S_{G-}^{(n-1)}\rangle$ and $|A_{G-}^{(n)}\rangle \rightarrow |A_{G-}^{(n-1)}\rangle$ [marked with solid gray arrows in Fig. 4(d)] lead to the emission of photons of frequency $\omega_{las} = \omega_0 - \Lambda$, whereas the one-photon transitions $|S_{G-}^{(n)}\rangle \rightarrow |A_{G-}^{(n-1)}\rangle$ and $|A_{G-}^{(n)}\rangle \rightarrow |S_{G-}^{(n-1)}\rangle$ give raise to the emission of photons of frequencies $\omega_{las} + |(\cos\theta - \sin\theta)\Omega|$ and $\omega_{las} - |(\cos\theta - \sin\theta)\Omega|$, respectively. As a consequence, three peaks emerge in the emission spectrum on top of Fig. 4(c), which are analogous to the Mollow triplet in the emission spectrum of a single TLS under strong enough illumination (see Appendix F). In the FRIC map, these one-photon transitions lead to a broad vertical line and a broad horizontal line, which emerge due to the overlapping of three vertical and horizontal lines at frequencies $\omega_{las}$ and $\omega_{las} \pm |(\cos\theta - \sin\theta)\Omega|$.

Moreover, new antidiagonal lines emerge in the FRIC map satisfying $\omega_1 + \omega_2 = 2\omega_{las} \pm |(\cos\theta - \sin\theta)\Omega|$ [see the light green and dark green dashed arrows in Fig. 4(c)]. These new lines unveil additional two-photon transitions through virtual states and emerge due to the energy splitting between the dressed states $|S_{G-}^{(n)}\rangle$ and $|A_{G-}^{(n)}\rangle$, which are linear combinations of the superradiant and ground molecular states. These lines are also analogous to the ones emerging in a single TLS under the same laser intensity (see Appendix F).

On the other hand, we observe additional lines in the FRIC map that are due to the non-negligible effect of the subradiant state and of the doubly excited state for this laser intensity. First, we observe two horizontal and vertical lines centered

at frequency $\omega_{las} + 2\Lambda$, which corresponds to the one-photon transitions from one of the dressed states $|A_{E+}^{(n)}\rangle$ and $|S_{E+}^{(n)}\rangle$ to one of the dressed states $|S_{G-}^{(n-1)}\rangle$ and $|A_{G-}^{(n-1)}\rangle$ in a lower rung. The FRIC map also exhibits two additional antidiagonal lines corresponding to the two-photon transitions from $|A_{E+}^{(n+1)}\rangle$ and $|S_{E+}^{(n+1)}\rangle$ to one of the dressed states $|S_{G-}^{(n-1)}\rangle$ and $|A_{G-}^{(n-1)}\rangle$ through an intermediate virtual state. In summary, the FRIC map in Fig. 4(c) again mostly resembles the FRIC map of a TLS consisting of the ground and the superradiant state and driven resonantly under the same laser intensity, but with additional lines accounting for the impact of the subradiant state and the doubly excited state on the light emission.

Remarkably, not all the transitions between eigenstates are reflected in the FRIC map. In particular, we do not observe horizontal and vertical lines centered at frequency $\omega_{las} - 2\Lambda = \omega_0 - 3\Lambda$ corresponding to the one-photon transitions from $|S_{G-}^{(n)}\rangle$ or $|A_{G-}^{(n)}\rangle$ to $|A_{E+}^{(n-1)}\rangle$ or $|S_{E+}^{(n-1)}\rangle$. Similarly, no peak appears at frequency $\omega_{las} - 2\Lambda = \omega_0 - 3\Lambda$ in the emission spectrum on top of Fig. 4(c). To understand this effect we calculate the probability of these transitions using Fermi's golden rule, which states that the transition probability between a initial eigenstate $|i\rangle$ and a final eigenstate $|f\rangle$ is proportional to $|(f)\langle(\sigma_1 + \sigma_2)|i\rangle|^2$. The emission operator can be written in the interacting basis as

$$ \sigma_1 + \sigma_2 = (\cos\theta + \sin\theta)(|G\rangle\langle\Lambda_+| + |\Lambda_+\rangle\langle E|) + (\cos\theta - \sin\theta)(|G\rangle\langle\Lambda_-| + |\Lambda_-\rangle\langle E|). \quad (15) $$

As a consequence, $\langle S_{E+}^{(n-1)}| (\sigma_1 + \sigma_2) |S_{G-}^{(n)}\rangle = \langle A_{E+}^{(n-1)}| (\sigma_1 + \sigma_2) |S_{G-}^{(n)}\rangle = \langle S_{E+}^{(n-1)}| (\sigma_1 + \sigma_2) |A_{G-}^{(n)}\rangle = \langle A_{E+}^{(n-1)}| (\sigma_1 + \sigma_2) |A_{G-}^{(n)}\rangle = 0$, which means that all the one-photon transitions with energy difference $\omega_{las} - 2\Lambda = \omega_0 - 3\Lambda$ are forbidden. Finally, we note that the FRIC map does not show antidiagonal lines corresponding to two-photon transitions from the states $|S_{G-}^{(n+1)}\rangle$ or $|A_{G-}^{(n+1)}\rangle$ to any of the states $|S_{E+}^{(n-1)}\rangle$ or $|A_{E+}^{(n-1)}\rangle$, which indicates that these transitions are also forbidden.

## 3. Strong laser intensity

Last, we plot in Fig. 4(e) the FRIC map for intense laser illumination, $I = 10^3 I_{sat}$. For this laser intensity, we obtain the dressed eigenstates $|E_i^{(n)}\rangle$ numerically.

One of the main differences with the case of moderate laser intensity is that the three Mollow-like vertical, horizontal, and antidiagonal lines (marked with light green, gray, and dark green arrows) in the FRIC map do not overlap and can be clearly distinguished. The reason is that the energy splitting between the dressed states $|E_1^{(n)}\rangle$ and $|E_2^{(n)}\rangle$ is much larger than for the case of moderate laser intensity.

Another noticeable effect of the increase of laser intensity is that we observe a wider variety of vertical, horizontal, and antidiagonal lines in the FRIC map. There are two main factors involved in the emergence of new lines. The first one is the larger energy splitting between $|E_3^{(n)}\rangle$ and $|E_4^{(n)}\rangle$, which was negligible in the case of moderate intensity, and allows for resolving new transitions. The second factor responsible of the emergence of new lines in the FRIC map is that the forbidden transitions for moderate laser intensities now become allowed due to the more complex nature of the dressed eigenstates,
<figure><img src="image_7.png" /><figcaption>FIG. 5. FRIC for a laser tuned to the transition frequency of the subradiant state $|\Lambda_+\rangle$, $\omega_{las} = \omega_0 + \Lambda$. [(a),(c),(e)] Normalized emission spectra (top panel) and FRIC maps (bottom panel) for laser intensity (a) $I = 0.1I_{sat}$, (c) $I = 50I_{sat}$, and (e) $I = 1000I_{sat}$. We depict the molecular energy levels $|G\rangle$ and $|\Lambda_+\rangle$ and the resonant excitation of the subradiant state in the inset in (a). We follow the standard color convention in the FRIC maps and use blue color to represent antibunched emission, red for bunched emission, and white for Poissonian statistics. We make use of a linear scale in the range $0 \le g^{(2)}(\omega_1, \omega_2; 0) \le 1$ and a logarithmic scale in the interval $1 \le g^{(2)}(\omega_1, \omega_2; 0) \le 100$ in the color bar. We assume a filter with linewidth $\Gamma = \gamma_0/10$. [(b),(d),(f)] Diagrams of the eigenstates, with the laser field quantized, in the rungs $n-1, n$, and $n+1$ for laser intensity (b) $I = 0.1I_{sat}$, (d) $I = 50I_{sat}$, and (f) $I = 1000I_{sat}$. Each of the rungs have four eigenstates with equal number of excitations. We use solid arrows to mark the horizontal and vertical lines in the FRIC maps and to indicate the one-photon transitions in the diagrams of the eigenstates, while we use dashed arrows for antidiagonal lines related to leapfrog processes. We have fixed the values $V = -20\gamma_0$, $\alpha = 0.3$, and $\delta = 10\gamma_0$.</figcaption></figure>

which under strong illumination are not simple symmetric and antisymmetric combinations of two states of the interacting basis $\{|G, n\rangle, |\Lambda_-, n-1\rangle, |\Lambda_+, n-1\rangle, |E, n-2\rangle\}$.

Finally, we remark that the diagonal bunched line at $\omega_1 = \omega_2$ corresponds to the simultaneous detection of two identical photons. The intensity correlation along this diagonal is twice that of at the nearby points in the map [34,36,70]. The reason is that a coincidence is counted no matter which photon arrives to detector 1 and which one does it to the detector 2, in contrast to the case in which $\omega_1 \ne \omega_2$, which doubles the probability of a coincidence in the detectors.

### C. Laser tuned to the subradiant state

Next, we analyze the FRIC for a laser tuned to the transition frequency of the subradiant state $|\Lambda_+\rangle$ (i.e., $\omega_{las} = \omega_0 + \Lambda$). For this laser frequency, the FRIC maps become more complex than for the previous detuning due to the emergence of additional lines and also to the possibility of having FRIC maps with both antibunched or bunched background depending on the laser intensity. However, we emphasize that the general principles governing the eigenstates of the system

are very similar to the case of the previous section, where the laser was tuned to the transition frequency of the superradiant state (see Appendix G).

#### 1. Weak laser intensity

For $I = 0.1I_{sat}$, we plot the FRIC map in Fig. 5(a) and a scheme of the energy levels in Fig. 5(b). As in the case of weak illumination in the previous section, the eigenstates of the full quantized system are $\{|G, n\rangle, |\Lambda_-, n-1\rangle, |\Lambda_+, n-1\rangle, |E, n-2\rangle\}$. Now, the states $|G, n\rangle$ and $|\Lambda_+, n-1\rangle$ in the rung $n$ are degenerated in energy and, in the same way, so are the states $|\Lambda_-, n-1\rangle$ and $|E, n-2\rangle$.

The background of the FRIC map is antibunched, which is consistent with the value of the color-blind intensity correlation $g^{(2)}(0) \approx 0.3$ for these laser parameters. Over this background, three horizontal, three vertical, and three antidiagonal lines emerge. Thus, in contrast with the previous section, this FRIC map does not resemble the one of a TLS under weak illumination. The reason is the competition between the resonant excitation of the subradiant state and the nonresonant excitation of the superradiant state (which
couples very efficiently to light), as discussed for the color-blind intensity correlation in Sec. III B. Both excitation pathways lead to the emission of photons of frequency $\omega_{las} = \omega_0 + \Lambda$ giving raise to the central horizontal and vertical lines in the FRIC map. Additionally, the central antidiagonal line in the FRIC map corresponds to the emission of two photons of frequencies $\omega_1$ and $\omega_2$, via an intermediate virtual state, fulfilling $\omega_1 + \omega_2 = 2\omega_{las}$, so that the final state is the same as the initial one but two rungs down. Further, photons of frequency $\omega_{las} - 2\Lambda = \omega_0 - \Lambda$ can be also emitted [lines marked with solid red arrows in the FRIC map and also in Fig. 5(b)], as well as leapfrog processes fulfilling $\omega_1 + \omega_2 = 2\omega_{las} - 2\Lambda$ can take place (red-dashed arrows).

Remarkably, the excitation of the subradiant state $|\Lambda_+, n\rangle$ in rung $n$ can lead to the emission of photons of frequency $\omega_{las} + 2\Lambda = \omega_0 + 3\Lambda$ due to the transition to any of states $|E, n - 2\rangle$ or $|\Lambda_-, n - 1\rangle$ in rung $n - 1$, which gives raise to the horizontal and vertical lines marked with purple solid arrows in the FRIC map in Fig. 5(a). However, the emission of these photons is not reflected in the emission spectrum on top of Fig. 5(a) and, additionally, these lines should not appear in the FRIC map according to the Fermi's golden rule (following similar arguments than in the previous section). We attribute the emergence of these lines in the map as a consequence of the very small modification of the eigenstates $\{|G, n\rangle, |\Lambda_-, n - 1\rangle, |\Lambda_+, n - 1\rangle, |E, n - 2\rangle\}$ due to the weak illumination (see Appendix G), so that the transition is no longer strictly forbidden. Therefore, for this laser detuning, the FRIC map reveals one-photon transitions that occur with very low probability and which are hidden in the emission spectrum. Similarly, two-photon transitions via intermediate virtual states satisfying $\omega_1 + \omega_2 = 2\omega_{las} + 2\Lambda$ are also allowed due to the very small dressing of the eigenstates (marked wit the dashed-purple arrows).

Furthermore, we analyze the impact of the combined Debye-Waller/Franck-Condon factor $\alpha$ on the FRIC map at this laser frequency and intensity in Appendix D. We find that the FRIC map can be drastically different when $\alpha = 1$, as typically assumed in the literature, which emphasizes the importance of accounting for $\alpha$ in the description of the emitters when intensity correlations are investigated.

## 2. Moderate and strong laser intensities

We present in Fig. 5(c) the FRIC map for a laser of intensity $I = 50I_{sat}$. For these laser parameters, the dressed eigenstates can be approximated again as symmetric and antisymmetric combinations of the states that are degenerate under weak illumination, i.e., $|E_1^{(n)}\rangle \approx |A_{E-}^{(n)}\rangle = (|E, n - 2\rangle - |\Lambda_-, n - 1\rangle)/\sqrt{2}$, $|E_2^{(n)}\rangle \approx |S_{E-}^{(n)}\rangle = (|E, n - 2\rangle + |\Lambda_-, n - 1\rangle)/\sqrt{2}$, $|E_3^{(n)}\rangle \approx |S_{G-}^{(n)}\rangle = (|G, n\rangle + |\Lambda_+, n - 1\rangle)/\sqrt{2}$, and $|E_4^{(n)}\rangle \approx |A_{G+}^{(n)}\rangle = (|G, n\rangle - |\Lambda_+, n - 1\rangle)/\sqrt{2}$.

In contrast to the case in which we tune the laser to the transition frequency of the superradiant state with the same laser intensity, we observe that the FRIC map in Fig. 5(c) is bunched $g^{(2)}(\omega_1, \omega_2; 0) > 1$ for most pairs of filtered frequencies $\omega_1$ and $\omega_2$. This general bunching can be expected from the bunched color-blind intensity correlation for these laser parameters [see Fig. 3(b)].

Over this general bunched background, a variety of horizontal, vertical, and antidiagonal lines are observed. Similar to the case of moderate laser intensity for the previous detuning, this wide variety of lines emerge again due to the significant energy splitting $|\Omega(\cos\theta - \sin\theta)| = 0.11\Lambda$ between the dressed states $|S_{E-}^{(n)}\rangle$ and $|A_{E-}^{(n)}\rangle$. This energy splitting increases the number of different transitions that lead to emission of photons [see Fig. 5(d)], with respect to the case of weak laser intensity.

Finally, we present in Fig. 5(e) the FRIC map for $I = 1000I_{sat}$, which also exhibits a significant bunching for almost every pair of filtered frequencies. For this intensity, the numerical calculation of the dressed eigenstates $|E_i^{(n)}\rangle$ in rung $n$ is again necessary. The energy splittings between all the dressed eigenstates $|E_i^{(n)}\rangle$ become significant, as schematically represented in Fig. 5(f). As a consequence, in the FRIC map we observe more horizontal, vertical, and antidiagonal lines than in the case of moderate laser intensity, with the corresponding one-photon and two-photon transitions depicted in Fig. 5(f).

## V. CONCLUSIONS

We have investigated light emission from two strongly interacting TLS-like quantum emitters. In particular, we have analyzed the dependence of color-blind and frequency-resolved intensity correlations on the laser intensity for a laser tuned to (i) the two-photon resonance, (ii) the transition frequency of the superradiant state, and (iii) the transition frequency of the subradiant state.

We have first demonstrated that the color-blind intensity correlation can be tailored from extreme bunching to strong antibunching by modifying the laser parameters. We have also shown how this dependence can be altered by changing the frequency detuning between the two emitters. For weak enough laser intensities, the color-blind intensity correlation is characterized by an extreme bunching if the laser is tuned to the two-photon resonance. In contrast, if the laser is tuned to the transition frequency of the superradiant state, the intensity correlation becomes strongly antibunched, so that the system behaves as an efficient single-photon source (with larger decay rate than that of a single emitter). For these two laser detunings, we have derived analytical expressions (in Appendix A), valid for nonidentical molecules, $\alpha \neq 1$ and arbitrary laser intensity, which accurately reproduce the intensity correlation.

Additionally, we have found that the color-blind intensity correlation shows a more complex dependence on the laser intensity and on the detuning between the emitters when the laser is tuned to the transition frequency of the subradiant state. For this laser detuning, light can be either bunched, antibunched or Poissonian, which is a signature of the complexity of the excitation and emission processes at this laser frequency. These emission and excitation processes may become even more complex in the presence of dephasing effects, which have been neglected in this work and require a more sophisticated description [71,72]. Besides, we have derived an analytical expression of $g^{(2)}(0)$ when the laser is tuned resonantly to the subradiant state at weak laser intensities. This expression, as well as analogous expressions for the
other laser detunings, can help to determine experimentally the distance and detuning between the emitters.

Furthermore, we have also explored the frequency-resolved intensity correlation (FRIC) for the three laser detunings and shown that they provide further information onto a variety of relaxation processes of the driven coupled molecules. More specifically, the FRIC accounts for the photons emitted via one-photon transitions between eigenstates of the system and also via two-photon transitions through intermediate virtual states.

For a laser tuned to the transition frequency of the super-radiant state, we have demonstrated that the FRIC generally resembles that of a resonantly driven TLS, but it also exhibits some differences due to the influence of the subradiant state and the doubly excited state. Interestingly, for weak and moderate laser intensities, the FRIC maps indicate that some of the transitions are forbidden (absence of lines at certain frequencies), which is consistent with the peaks in the emission spectrum. On the other hand, when the laser is tuned to the transition frequency of the subradiant state the FRIC becomes more complex due to the important competition between the resonant excitation of the subradiant state and the nonresonant excitation of the superradiant state. We find that the FRIC can unveil one-photon transitions that have such low probability to occur that are not resolved by the emission spectrum. Thus, the FRIC manifests a stronger sensitivity to reveal transitions that occur with very low probability. These results stress the interest of the intensity correlations to better characterize and understand the emission from coupled quantum emitters.

The data used for the figures in this article can be found at [73].

## ACKNOWLEDGMENTS

We thank Alejandro Gonzalez-Tudela and Mikolaj K. Schmidt for fruitful discussions. A.J.-D., R.E., A.N., and J.A. acknowledge financial support through the Grant No.

PID2022-139579NB-I00 funded by the Spanish Government Agency of Research MCIN/AEI/10.13039/501100011033 and by "ERDF A way of making Europe" and also through the Grant IT 1526-22 funded by the Department of Education, Research and Universities of the Basque Government. A.J.-D. acknowledges financial support through the Grant No. PRE2020-095013 funded by the Spanish Government Agency of Research MCIN/AEI/10.13039/501100011033 and by "ESF Investing in your future". J.-B.T. and B.L. acknowledge the financial support from the French National Agency for Research (ANR-22-CE47-0015), Région Nouvelle-Aquitaine, Idex Bordeaux (Research Program GPR Light), and the EUR Light S&T (PIA3 Program, ANR-17-EURE-0027). We acknowledge financial support from the Laboratory for Trans-border Cooperation LTC TRANS-LIGHT from University of Bordeaux and University of the Basque Country.

## APPENDIX A: STEADY-STATE ANALYTICAL RESOLUTION

In this Appendix we derive the equations of motion of all the elements of the density matrix in the interacting basis {|G⟩, |Λ−⟩, |Λ+⟩, |E⟩}. We use these equations of motion to obtain analytical expressions of the steady-state populations and color-blind intensity correlation g(2)(0) when the laser is tuned to the two-photon resonance and also to the transition frequency of the superradiant state |Λ−⟩. According to the general expression in Eq. (9), the intensity correlation is obtained from the steady-state coherence ρ+− = ⟨Λ+|ρss|Λ−⟩ and from the steady-state populations of the superradiant state |Λ−⟩, of the subradiant state |Λ+⟩, and of the doubly excited state |E⟩. To obtain these quantities, we write the master equation [see Eq. (6)] in the interacting basis {|G⟩, |Λ+⟩, |Λ−⟩, |G⟩} and derive the equations of motion of the populations and coherences in this basis.

First, the incoherent Liouvillian given in Eq. (7) is given in the interacting basis as

$$ \mathcal{L}_{\text{inc}} \hat{\rho} = \frac{\gamma_+}{2} (\mathcal{D}[\hat{\sigma}_{G+}]\hat{\rho} + \mathcal{D}[\hat{\sigma}_{+E}]\hat{\rho}) + \frac{\gamma_-}{2} (\mathcal{D}[\hat{\sigma}_{G-}]\hat{\rho} + \mathcal{D}[\hat{\sigma}_{-E}]\hat{\rho}) + \frac{\tilde{\gamma}_+}{2} (\mathcal{D}[\hat{\sigma}_{G+}, \hat{\sigma}_{+E}]\hat{\rho} + \mathcal{D}[\hat{\sigma}_{+E}, \hat{\sigma}_{G+}]\hat{\rho}) \\ + \frac{\tilde{\gamma}_-}{2} (\mathcal{D}[\hat{\sigma}_{G-}, \hat{\sigma}_{-E}]\hat{\rho} + \mathcal{D}[\hat{\sigma}_{-E}, \hat{\sigma}_{G-}]\hat{\rho}) + \frac{\tilde{\gamma}_0}{2} (\mathcal{D}[\hat{\sigma}_{G-}, \hat{\sigma}_{+E}]\hat{\rho} + \mathcal{D}[\hat{\sigma}_{+E}, \hat{\sigma}_{G-}]\hat{\rho}) \\ + \frac{\tilde{\gamma}_0}{2} (\mathcal{D}[\hat{\sigma}_{G+}, \hat{\sigma}_{-E}]\hat{\rho} + \mathcal{D}[\hat{\sigma}_{-E}, \hat{\sigma}_{G+}]\hat{\rho}) + \frac{\tilde{\gamma}_{12}}{2} (\mathcal{D}[\hat{\sigma}_{G+}, \hat{\sigma}_{G-}]\hat{\rho} + \mathcal{D}[\hat{\sigma}_{G-}, \hat{\sigma}_{G+}]\hat{\rho}) \\ + \frac{\tilde{\gamma}_{12}}{2} (\mathcal{D}[\hat{\sigma}_{+E}, \hat{\sigma}_{-E}]\hat{\rho} + \mathcal{D}[\hat{\sigma}_{-E}, \hat{\sigma}_{+E}]\hat{\rho}), \qquad (\text{A1}) $$

where $\tilde{\gamma}_{\pm} = \tilde{\gamma} \pm \gamma_0 \sin 2\theta$, $\tilde{\gamma}_0 = \gamma_0 \cos 2\theta$, and $\tilde{\gamma}_{12} = \tilde{\gamma} \cos 2\theta$. We have also introduced the operators $\hat{\sigma}_{EE} = |E\rangle\langle E|$, $\hat{\sigma}_{++} = |\Lambda_+\rangle\langle \Lambda_+|$, $\hat{\sigma}_{--} = |\Lambda_-\rangle\langle \Lambda_-|$, $\hat{\sigma}_{EG} = |E\rangle\langle G|$, and $\hat{\sigma}_{+-} = |\Lambda_+\rangle\langle \Lambda_-|$. The terms in the first line of Eq. (A1) correspond to one-photon transitions between states of the interacting basis, whereas the rest of the terms correspond to crossed incoherent processes.

Next, we compute from Eqs. (5) and (A1) the set of equations of motion of the mean value 〈·〉 of the operators $\hat{\sigma}_{ii}$, whose expected values correspond to the populations of the states of the interacting basis, and the mean values of the coherences $\hat{\sigma}_{ij}$. The resulting equations are

$$ \frac{d}{dt} \langle \hat{\sigma}_{EE} \rangle = -2\gamma_0 \langle \hat{\sigma}_{EE} \rangle + i \frac{\Omega_+}{2} \langle \hat{\sigma}_{E+} - \hat{\sigma}_{+E} \rangle + i \frac{\Omega_-}{2} \langle \hat{\sigma}_{E-} - \hat{\sigma}_{-E} \rangle, \qquad (\text{A2a}) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{++} \rangle = -\gamma_+ \langle \hat{\sigma}_{++} - \hat{\sigma}_{EE} \rangle - \frac{\tilde{\gamma}_{12}}{2} \langle \hat{\sigma}_{-+} + \hat{\sigma}_{+-} \rangle + i \frac{\Omega_+}{2} \langle \hat{\sigma}_{+E} - \hat{\sigma}_{E+} + \hat{\sigma}_{+G} - \hat{\sigma}_{G+} \rangle, \qquad (\text{A2b}) $$
<figure><img src="image_8.png" /><figcaption>FIG. 6. Dependence of the steady-state coherences on the laser intensity for different laser detunings. Absolute value of the coherences when the laser is tuned to (a) the two-photon resonance, (b) the transition frequency of the superradiant state $|\Lambda_-\rangle$, and (c) the transition frequency of the subradiant state $|\Lambda_+\rangle$. The insets in (a)–(c) illustrate the different resonant excitation processes promoted by the different laser frequencies. Solid lines represent the coherences that involve the subradiant state $|\Lambda_+\rangle$, namely (red) $\rho_{+-} = \langle\Lambda_+|\hat{\rho}_{ss}|\Lambda_-\rangle$, (purple) $\rho_{+G} = \langle\Lambda_+|\hat{\rho}_{ss}|G\rangle$, and (brown) $\rho_{E+} = \langle E|\hat{\rho}_{ss}|\Lambda_+\rangle$. Dashed lines are used to represent the coherences that do not involve the subradiant state, namely (orange) $\rho_{-G} = \langle\Lambda_-|\hat{\rho}_{ss}|G\rangle$, (green) $\rho_{E-} = \langle E|\hat{\rho}_{ss}|\Lambda_-\rangle$, and (blue) $\rho_{EG} = \langle E|\hat{\rho}_{ss}|G\rangle$. The dipole-dipole coupling is fixed as $V = -20\gamma_0$, the combined Debye-Waller/Franck-Condon factor as $\alpha = 0.3$ and the molecular detuning as $\delta = 10\gamma_0$.</figcaption></figure>

$$ \frac{d}{dt} \langle \hat{\sigma}_{--} \rangle = -\gamma_{-} \langle \hat{\sigma}_{--} - \hat{\sigma}_{EE} \rangle - \frac{\tilde{\gamma}_{12}}{2} \langle \hat{\sigma}_{+-} + \hat{\sigma}_{+-} \rangle + i \frac{\Omega_{-}}{2} \langle \hat{\sigma}_{-E} - \hat{\sigma}_{-E} + \hat{\sigma}_{-G} - \hat{\sigma}_{-G} \rangle, \quad (A2c) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{EG} \rangle = \langle \hat{\sigma}_{EG} \rangle (i2\Delta_0 - \gamma_0) + i \frac{\Omega_+}{2} \langle \hat{\sigma}_{E+} - \hat{\sigma}_{+G} \rangle + i \frac{\Omega_-}{2} \langle \hat{\sigma}_{E-} - \hat{\sigma}_{-G} \rangle, \quad (A2d) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{E+} \rangle = -i \left( \frac{\Omega_+}{2} (-\hat{\sigma}_{EE} + \hat{\sigma}_{++} - \hat{\sigma}_{EG}) + \frac{\Omega_-}{2} \hat{\sigma}_{+-} \right) + \langle \hat{\sigma}_{E+} \rangle \left( i(\Delta_0 - \Lambda) - \gamma_0 - \gamma_+/2 \right) - \frac{\tilde{\gamma}_{12}}{2} \langle \hat{\sigma}_{E-} \rangle, \quad (A2e) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{E-} \rangle = -i \left( \frac{\Omega_-}{2} (-\hat{\sigma}_{EE} + \hat{\sigma}_{--} - \hat{\sigma}_{EG}) + \frac{\Omega_+}{2} \hat{\sigma}_{+-} \right) + \langle \hat{\sigma}_{E-} \rangle \left( i(\Delta_0 + \Lambda) - \gamma_0 - \gamma_-/2 \right) - \frac{\tilde{\gamma}_{12}}{2} \langle \hat{\sigma}_{E+} \rangle, \quad (A2f) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{+G} \rangle = -i \left( \frac{\Omega_+}{2} (-\hat{\sigma}_{++} + \hat{\sigma}_{GG} + \hat{\sigma}_{EG}) - \frac{\Omega_-}{2} \hat{\sigma}_{+-} \right) + \langle \hat{\sigma}_{+G} \rangle \left( i(\Delta_0 + \Lambda) - \gamma_+/2 \right) \langle \hat{\sigma}_{+G} \rangle - \frac{\tilde{\gamma}_{12}}{2} \langle \hat{\sigma}_{-G} \rangle + \tilde{\gamma}_+ \langle \hat{\sigma}_{E+} \rangle + \tilde{\gamma}_0 \langle \hat{\sigma}_{E-} \rangle, \quad (A2g) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{-G} \rangle = -i \left( \frac{\Omega_-}{2} (-\hat{\sigma}_{--} + \hat{\sigma}_{GG} + \hat{\sigma}_{EG}) - \frac{\Omega_+}{2} \hat{\sigma}_{+-} \right) + \langle \hat{\sigma}_{-G} \rangle \left( i(\Delta_0 - \Lambda) - \gamma_-/2 \right) \langle \hat{\sigma}_{-G} \rangle - \frac{\tilde{\gamma}_{12}}{2} \langle \hat{\sigma}_{+G} \rangle + \tilde{\gamma}_- \langle \hat{\sigma}_{E-} \rangle + \tilde{\gamma}_0 \langle \hat{\sigma}_{E+} \rangle, \quad (A2h) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{+-} \rangle = -i \frac{\Omega_+}{2} \langle \hat{\sigma}_{E-} + \hat{\sigma}_{G-} \rangle + i \frac{\Omega_-}{2} \langle \hat{\sigma}_{+E} + \hat{\sigma}_{+G} \rangle + (i2\Lambda - \gamma_0) \langle \hat{\sigma}_{+-} \rangle + \frac{\tilde{\gamma}_{12}}{2} \langle 2\hat{\sigma}_{EE} - \hat{\sigma}_{++} - \hat{\sigma}_{--} \rangle, \quad (A2i) $$

The above system of equations is valid for arbitrary laser intensity and detuning, but its exact analytical resolution leads to very complex expressions. Hence, it becomes convenient to perform some approximations when looking for simple analytical expressions of the intensity correlation. For example, when the laser is either tuned to the two-photon resonance ($\Delta_0 = 0$) or to the transition frequency of the superradiant state $|\Lambda_-\rangle$ ($\Delta_0 = \Lambda$), the steady-state coherences between the subradiant state $|\Lambda_+\rangle$ and any other state are much smaller than the largest coherence of the system (in the interacting basis) for any value of laser intensity (see Fig. 6). Thus, we can neglect the steady-state coherences related to the subradiant state for these two laser detunings. As a consequence, the set of Eqs. (A2) is reduced in the steady state ($\frac{d}{dt}\hat{\rho} = 0$) to

$$ 0 = -2\gamma_0\rho_E - \Omega_- \mathrm{Im}\rho_{E-}, \quad (A3a) $$

$$ 0 = -\gamma_+(\rho_{\Lambda+} - \rho_E), \quad (A3b) $$

$$ 0 = -\gamma_{-}(\rho_{\Lambda_{-}} - \rho_{E}) + \Omega_{-}(\mathrm{Im}\rho_{E_{-}} - \mathrm{Im}\rho_{-G}), \quad (A3c) $$

$$ 0 = (i2\Delta_0 - \gamma_0)\rho_{EG} + i \frac{\Omega_-}{2} (\rho_{E-} - \rho_{-G}), \quad (A3d) $$

$$ 0 = -i \frac{\Omega_-}{2} (-\rho_E + \rho_{\Lambda_-} - \rho_{EG}) + (i(\Delta_0 + \Lambda) - \gamma_0 - \gamma_-/2)\rho_{E-}, \quad (A3e) $$

$$ 0 = -i \frac{\Omega_-}{2} (-\rho_{\Lambda_-} + \rho_G + \rho_{EG}) + (i(\Delta_0 - \Lambda) - \gamma_-/2)\rho_{-G} + \tilde{\gamma}_- \rho_{E-}, \quad (A3f) $$

where $\rho_{ij} = \mathrm{Tr}(\hat{\sigma}_{ij}\hat{\rho}_{ss})$ and $\rho_i = \mathrm{Tr}(\hat{\sigma}_{ii}\hat{\rho}_{ss})$. This simplified set of equations gives raise to manageable expressions for the steady-state populations and the color-blind intensity correlation $g^{(2)}(0)$. We note that, as a consequence of Eq. (A3b), the analytical population of the doubly excited state $\rho_E$ is identical to the analytical population of the subradiant state

$\rho_{\Lambda+}.$
For the case in which the laser is tuned to the two-photon resonance ($2\omega_{las} = 2\omega_0$), we substitute $\Delta_0 = 0$ into Eq. (A3) and solve the set of equations. The population of $|E\rangle$ and $|\Lambda_+\rangle$ are equal $\rho_E|2^{PR} = \rho_{\Lambda_+}|2^{PR}$ (with $2PR$ standing for two-photon resonance) and given by

$$ \rho_E \Big|^{2^{PR}} = \frac{\Omega_-^4}{4\Omega_-^4 + 2\Omega_-^2\gamma_0(4\gamma_0 + 3\gamma_- - 4\tilde{\gamma}_-) + (4\gamma_0\Lambda)^2}, \quad (A4) $$

which converges to $1/4$ for strong laser intensities. The population of the superradiant state $|\Lambda_-\rangle$ becomes related to the above population as

$$ \rho_{\Lambda_{-}} |^{2^{PR}} = \rho_{E} |^{2^{PR}} \left( 1 + \left( \frac{2\gamma_{0}}{\Omega_{-}} \right)^{2} \right), \quad (A5) $$

which leads to $\rho_{\Lambda_-}|_{I \gg I_{sat}}^{2^{PR}} = \rho_E|_{I \gg I_{sat}}^{2^{PR}} = \rho_{\Lambda_+}|_{I \gg I_{sat}}^{2^{PR}} = 1/4$ for strong enough laser intensities, but $\rho_{\Lambda_-}|_{I \ll I_{sat}}^{2^{PR}} \gg \rho_E|_{I \ll I_{sat}}^{2^{PR}} = \rho_{\Lambda_+}|_{I \ll I_{sat}}^{2^{PR}}$ for weak intensities. Substituting these analyti-

cal populations into Eq. (9), with $Re\rho_{+-} = 0$, the intensity correlation results

$$ g^{(2)}(0) \Big|^{2^{PR}} = \frac{\Omega_-^4 + \Omega_-^2\gamma_0(2\gamma_0 + 3\gamma_-/2 - 2\tilde{\gamma}_-) + (2\gamma_0\Lambda)^2}{[\Omega_-^2 + \gamma_0^2(1 - \sin 2\theta)]^2}. \quad (A6) $$

Similarly, we solve the set of Eqs. (A3) with $\Delta_0 = \Lambda$ to obtain the populations and intensity correlation for a laser tuned to the transition frequency of the superradiant state. The resulting steady-state populations of the doubly excited state and the subradiant state are given as

$$ \rho_E \Big|^{|\Lambda_->} = \frac{\Omega_-^4}{4\Omega_-^4 + 2\Omega_-^2[(2\xi)^2 + \gamma_0(\gamma_- - 4\tilde{\gamma}_-) + (2\xi\gamma_-)^2]}, \quad (A7) $$

with $\xi = \sqrt{(2\Lambda)^2 + (\gamma_0 + \gamma_-/2)^2}$ and which also converges to $1/4$ for strong laser intensities. The population of the superradiant state now is

$$ \rho_{\Lambda_{-}} |^{|\Lambda_{-}>} = \rho_{E} |^{|\Lambda_{-}>} + \frac{\Omega_{-}^{2}(2\xi)^{2}}{4\Omega_{-}^{4} + 2\Omega_{-}^{2}[(2\xi)^{2} + \gamma_{0}(\gamma_{-} - 4\tilde{\gamma}_{-})] + (2\xi\gamma_{-})^{2}}, \quad (A8) $$

and the intensity correlation thus becomes

$$ g^{(2)}(0) \Big|^{|\Lambda_{-}>} = \frac{4\Omega_{-}^{4} + 2\Omega_{-}^{2}[(2\xi)^{2} + \gamma_{0}(\gamma_{-} - 4\tilde{\gamma}_{-})] + (2\xi\gamma_{-})^{2}}{4[\Omega_{-}^{2} + (1 - \sin 2\theta)\xi^{2}]^{2}}, \quad (A9) $$

## APPENDIX B: DEPENDENCE OF $g^{(2)}(0)$ ON THE LASER INTENSITY WHEN THE LASER IS TUNED TO THE TWO-PHOTON RESONANCE

In this Appendix, we consider that the molecules are driven by a laser tuned to the two-photon resonance ($2\omega_{las} = 2\omega_0$) and discuss in detail the dependence of $g^{(2)}(0)$ on the laser intensity at $\delta = 10\gamma_0$ [red line in Fig. 2(d) and black-solid line in Fig. 7(a)]. To gain analytical insights, we use the expressions of the steady-state populations in Eqs. (A4) and (A5) [plotted with dots in Fig. 7(b)], which agree well with the numerical results [solid lines in Fig. 7(b)].

In the laser intensity range $I \ll I_{sat}$, the superradiant state $|\Lambda_-\rangle$ is excited through single-photon processes and its population $\rho_{\Lambda_-}$ increases linearly with the laser intensity according to $\rho_{\Lambda_-}|_{I \ll I_{sat}}^{2^{PR}} = \frac{I}{I_{sat}} \frac{1 - \sin 2\theta}{2\eta^2}$ [obtained from Eq. (A5) and where $\eta = \frac{2\Lambda}{\gamma_0}$]. Further, the doubly excited state $|E\rangle$ is excited through two-photon processes and its population depends quadratically on the intensity $\rho_E|_{I \ll I_{sat}}^{2^{PR}} = \frac{I^2}{I_{sat}^2} \frac{(1 - \sin 2\theta)^2}{16\eta^2}$ [obtained from Eq. (A4)].

Additionally, within this intensity range ($I \ll I_{sat}$), we observe in Fig. 7(b) that the population $\rho_{\Lambda_+}$ of the subradiant state $|\Lambda_+\rangle$ is two orders of magnitude smaller than the population $\rho_{\Lambda_-}$ of the superradiant state $|\Lambda_-\rangle$ and thus $\rho_{\Lambda_+}$ can

be neglected in the calculation of the intensity correlation in Eq. (10). As a consequence, both the numerator ($4\rho_E$) and the denominator ($\rho_{\Lambda_-}^2(1 - \sin 2\theta)^2$) in Eq. (10) scale quadratically with the intensity and $g^{(2)}(0)$ becomes independent of the particular value of $I/I_{sat}$ in this range of laser intensity. The resulting expression is given by

$$ g^{(2)}(0) \Big|_{I \ll I_{sat}}^{2^{PR}} = \left( \frac{\eta}{1 - \sin 2\theta} \right)^2, \quad (B1) $$

which agrees with the numerical results in Fig. 7(a) and becomes equal to 439.5 in our configuration. Similar to the analytical expression in Eq. (11) of $g^{(2)}(0)$ at weak enough illumination and for a laser tuned resonantly to the subradiant state, the simple expression of $g^{(2)}(0)$ in Eq. (B1) can also be used for the experimental estimation of molecular parameters. Moreover, Eq. (B1) shows that the emission is strongly bunched [$g^{(2)}(0)|_{I \ll I_{sat}}^{2^{PR}} \gg 1$] if the splitting $2\Lambda$ is much larger than the spontaneous decay $\gamma_0$ (as it occurs when $V \gg \gamma_0$).

For increasing laser intensities ($I \gtrsim 10^{-1}I_{sat}$), the intensity correlation strongly decreases. Figure 7(b) shows that the population $\rho_E$ of the doubly excited state $|E\rangle$ becomes comparable to the population $\rho_{\Lambda_-}$ of the superradiant state $|\Lambda_-\rangle$ in this case and, as a consequence, Eq. (10) is no longer a good approximation. By relating $\Omega$ with $I$, we can rewrite
Eq. (A6) as

$$g^{(2)}(0) = \frac{(4\Lambda)^2 + 2\gamma_0(1 - \sin 2\theta)(2\gamma_0 + 3\gamma_-/2 - 2\tilde{\gamma}_-)(I/I_{\text{sat}}) + (1 - \sin 2\theta)^2\gamma_0^2(I/I_{\text{sat}})^2}{[(2 + I/I_{\text{sat}})(1 - \sin 2\theta)\gamma_0]^2} \quad (B2)$$

with $\tilde{\gamma}_- = \tilde{\gamma} - \gamma_0 \sin 2\theta$. Equation (B2) is consistent with the expression reported in Ref. [74] for identical emitters ($\delta = 0$) and where the combined Debye-Waller/Franck-Condon factor is not taken into account (i.e., $\alpha = 1$). Figure 7(a) shows the excellent agreement between Eq. (B2) (blue dots) and the numerical results (black-solid line), including the convergence of both to $g^{(2)}(0) = 1$ for strong enough laser intensities.

<figure><img src="image_9.png" /><figcaption>FIG. 7. Dependence on the laser intensity of (a) $g^{(2)}(0)$ and (b) the steady-state populations when the laser is tuned resonantly to the two-photon resonance. In (a), the solid-black line is obtained numerically using Eq. (9), whereas blue dots are calculated with the full analytical expression in Eq. (B2). The inset in (a) represents the excitation at the two-photon resonance ($2\omega_{\text{las}} = 2\omega_0$). In (b), blue color corresponds to the steady-state populations of the superradiant state $|\Lambda_-\rangle$, green to that of the subradiant state $|\Lambda_+\rangle$, red to that of the doubly excited state $|E\rangle$ and pink to absolute value of the real part of the steady-state coherence $\rho_{+-} = \langle \Lambda_+ | \hat{\rho}_{ss} | \Lambda_- \rangle$. Solid lines correspond to the numerical calculations, whereas dots represent the analytical results obtained from the expressions in Eqs. (A4) and (A5). We have fixed the values $V = -20\gamma_0$, $\alpha = 0.3$, and $\delta = 10\gamma_0$.</figcaption></figure>

## APPENDIX C: DEPENDENCE OF $g^{(2)}(0)$ ON THE LASER INTENSITY WHEN THE LASER IS TUNED RESONANTLY TO THE SUPERRADIANT STATE

We examine next the color-blind intensity correlation when the laser is tuned to the transition frequency of the

<figure><img src="image_10.png" /><figcaption>FIG. 8. Dependence on the laser intensity of (a) the color-blind intensity correlation and (b) steady-state populations when the laser is tuned resonantly to the superradiant state $|\Lambda_+\rangle$. In (a), the solid-black line is obtained numerically using Eq. (9), while orange dots are calculated with the analytical expression in Eq. (C2) and blue dots with the analytical expression in Eq. (C3). The resonant excitation at the transition frequency of the superradiant state is illustrated in the inset in (a). In (b), color blue corresponds to the steady-state population of the superradiant state $|\Lambda_-\rangle$, green to that of the subradiant state $|\Lambda_+\rangle$, red to that of the doubly excited state $|E\rangle$, and pink to the absolute value of the real part of the steady-state coherence $\rho_{+-} = \langle \Lambda_+ | \hat{\rho}_{ss} | \Lambda_- \rangle$. Solid lines correspond to the numerical calculations, whereas dots represent the analytical results obtained from the expressions in Eqs. (A7) and (A8). We have fixed the values $V = -20\gamma_0$, $\alpha = 0.3$, and $\delta = 10\gamma_0$.</figcaption></figure>
superradiant state $|\Lambda_{-}\rangle$, i.e., $\omega_{las} = \omega_{0} - \Lambda$. In this case, the single-excitation delocalized state $|\Lambda_{-}\rangle$, which couples more efficiently to light than $|\Lambda_{+}\rangle$, is excited resonantly so that under weak enough illumination ($I \ll I_{sat}$), it dominates the response and thus the effect of the doubly excited state $|E\rangle$ and of the subradiant state $|\Lambda_{+}\rangle$ is negligible. The system then behaves in close analogy to a resonantly driven TLS, with strongly antibunched emission, as shown by the solid-black line in Fig. 8(a) [corresponding to the blue line in Fig. 2(d)]. This TLS-like behavior is consistent with the populations in Fig. 8(b), because (i) the population $\rho_{\Lambda_{-}}$ of the superradiant state $|\Lambda_{-}\rangle$ is much larger than the population $\rho_{\Lambda_{+}}$ of the subradiant state $|\Lambda_{+}\rangle$, and (ii) the square of the population $\rho_{\Lambda_{-}}$ of the superradiant state $|\Lambda_{-}\rangle$ is also much larger than the population $\rho_{E}$ of the doubly excited state $|E\rangle$. Thus, Eq. (10) is again a good approximation and indicates that $g^{(2)}(0) \ll 1$.

To obtain further analytical insights of the intensity correlation, we employ the expressions of the steady-state populations in Eqs. (A7) and (A8) [plotted with dots in Fig. 8(b)]. In the laser intensity range $I \ll I_{sat}$, these analytical expressions reduce to $\rho_{\Lambda_{-}}|_{I \ll I_{sat}}^{|\Lambda_{-}|} = \frac{I}{I_{sat}}\frac{\gamma_{0}^{2}(1-\sin 2\theta)}{2\gamma_{-}^{2}}$ and $\rho_{E}|_{I \ll I_{sat}}^{|\Lambda_{-}|} = \rho_{\Lambda_{+}}|_{I \ll I_{sat}}^{|\Lambda_{-}|} = \frac{I^{2}}{I_{sat}^{2}}(\frac{\gamma_{0}^{2}(1-\sin 2\theta)}{4\xi\gamma_{-}})^{2}$, with $\xi^{2} = (2\Lambda)^{2} + (\gamma_{0} + \gamma_{-}/2)^{2}$. Thus, the population $\rho_{\Lambda_{-}}$ of the superradiant state $|\Lambda_{-}\rangle$ depends again linearly on the laser intensity, and the population $\rho_{E}$ of the doubly excited $|E\rangle$ does it quadratically. The population of the subradiant state $|\Lambda_{+}\rangle$, which can be safely neglected in the denominator of Eq. (10), also depends quadratically on the laser intensity, because it is populated through the relaxation of the doubly excited state $|E\rangle$ for these laser intensities. Thus, both the numerator and the denominator scale quadratically in Eq. (10) and the intensity correlation becomes again independent of the laser intensity for weak enough intensities

$$g^{(2)}(0)|_{I \ll I_{sat}}^{|\Lambda_{-}|} = \left[ \frac{\gamma_{-}}{\xi(1 - \sin 2\theta)} \right]^{2}. \quad (C1)$$

which reproduces the numerical results for arbitrary laser intensity [compare blue dots and black line in Fig. 8(a)], including the saturation to $g^{(2)}(0) = 1$ for very large intensities.

## APPENDIX D: IMPACT OF THE COMBINED DEBYE-WALLER/FRANCK-CONDON FACTOR $\alpha$ ON THE LIGHT EMISSION STATISTICS

In this Appendix we analyze the impact of the combined Debye-Waller/Franck-Condon factor $\alpha$ on the light statistics. First, we investigate the influence of $\alpha$ on the color-blind intensity correlation $g^{(2)}(0)$. Figure 9 shows the dependence of $g^{(2)}(0)$ on the laser intensity $I$ and $\alpha$, calculated at $\delta = 10\gamma_{0}$. We fix the dipole-dipole coupling at $V = -20\gamma_{0}$ [i.e., $r_{12}$ is modified simultaneously with $\alpha$ according to Eq. (3)] to focus on the impact of $\alpha$ on the decay rates $\gamma_{-}$ and $\gamma_{+}$ of the superradiant state $|\Lambda_{-}\rangle$ and the subradiant state $|\Lambda_{+}\rangle$, respectively.

Equation (C1) becomes equal to $2.5 \times 10^{-4}$ for our molecular configuration. Crucially, this equation shows that as long as the losses are much smaller than $\xi$ (e.g., strong dipole-dipole coupling), then the intensity correlation is characterized by strong antibunching (in contrast to the bunching obtained when the laser is tuned to the two-photon resonance) and thus the coupled emitters behave as a single-photon source with larger emission rate and transition dipole moment than a single molecule.

The value of $g^{(2)}(0)$ increases for stronger laser intensities ($I \gtrsim 10^{-1}I_{sat}$). This increase is due to the fact that the population $\rho_{\Lambda_{-}}$ of the superradiant state $|\Lambda_{-}\rangle$ reaches a maximum value 0.49 (which is very close to the excited-state population of a saturated TLS) and remains almost constant, while the population of the doubly excited state $|E\rangle$ now increases linearly with the laser intensity. As a consequence, $g^{(2)}(0)$ also increases linearly with the laser intensity and is given by

$$g^{(2)}(0)|_{I \lesssim I_{sat}}^{|\Lambda_{-}|} = \frac{\gamma_{-}^{2} + (1 - \sin 2\theta)\gamma_{0}^{2}I/I_{sat}}{\xi^{2}(1 - \sin 2\theta)^{2}}. \quad (C2)$$

We plot the dependence of $g^{(2)}(0)$ on the laser intensity according to this equation (orange dots) in Fig. 8(a), which shows an excellent agreement with the numerical results (solid-black line) up to $I \approx 10^{3}I_{sat}$. We remark that the population of the doubly excited state $|E\rangle$ increases linearly with the laser intensity due to the large and constant population of the superradiant state.

For even stronger laser intensities, Eq. (10) is no longer a good approximation because the contributions of the population of the doubly excited state and subradiant state to the denominator of Eq. (9) need to be accounted for. By relating $\Omega$ with $I$, Eq. (A9) can be rewritten as

$$g^{(2)}(0)|_{I \ll I_{sat}}^{|\Lambda_{-}|} = \frac{(2\xi\gamma_{-})^{2} + (I/I_{sat})\gamma_{0}^{2}(1 - \sin 2\theta)[(2\xi)^{2} + \gamma_{0}(\gamma_{-} - 4\tilde{\gamma}_{-})] + (I/I_{sat})^{2}\gamma_{0}^{4}(1 - \sin 2\theta)^{2}}{(1 - \sin 2\theta)^{2}[(I/I_{sat})\gamma_{0}^{2} + 2\xi^{2}]^{2}}, \quad (C3)$$

radiant state $|\Lambda_{-}\rangle$ and the subradiant state $|\Lambda_{+}\rangle$, respectively. For identical emitters ($\delta = 0$) and $\alpha = 1$, as usually assumed in the literature, the decay rate of the superradiant state yields $\gamma_{-} = 2\gamma_{0}$ at short distances, while the decay rate of the subradiant state vanishes. In contrast, smaller values of $\alpha$ reduces considerably the difference between $\gamma_{-}$ and $\gamma_{+}$ according to $\gamma_{\pm} = \gamma_{0}(1 \pm \alpha \sin 2\theta)$ (here, $\sin 2\theta = V/\Lambda = -0.97$).

We find that the modification of $\gamma_{\pm}$ through $\alpha$ does not significantly impact on the color-blind intensity correlation $g^{(2)}(0)$ when the laser is tuned to the two-photon resonance [Fig. 9(a)] or to the transition frequency of the superradiant state $|\Lambda_{-}\rangle$ [Fig. 9(b)]. In contrast, when the laser is tuned resonantly to the subradiant state $|\Lambda_{+}\rangle$, Fig. 9(c) shows that the minimum and maximum values of $g^{(2)}(0)$ can be notably modified by $\alpha$. However, the qualitative behavior, including the presence of a crossover from antibunched emission (for
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

FIG. 9. Dependence of the color-blind intensity correlation $g^{(2)}(0)$ on the laser intensity $I$ and on the combined Debye-Waller/Franck-Condon factor $\alpha$. The laser is tuned to (a) the two-photon resonance ($2\omega_{las} = 2\omega_0$), (b) the transition frequency of the superradiant state $|\Lambda_\rangle$ ($\omega_{las} = \omega_0 - \Lambda$), and (c) the transition frequency of the subradiant state $|\Lambda_+\rangle$ ($\omega_{las} = \omega_0 + \Lambda$). In (a)–(c), color blue represents antibunched emission [$g^{(2)}(0) < 1$], red bunched emission [$g^{(2)}(0) > 1$] and white Poissonian statistics [$g^{(2)}(0) = 1$]. Cuts of the numerical results in (a)–(c) for $I = 10^{-3}I_{sat}$ are plotted in (d), where the red line corresponds to a laser tuned to the two-photon resonance, the blue line is obtained for a laser tuned to the superradiant state and the green line corresponds to a laser tuned to the subradiant state. Dots in (d) correspond to analytical results obtained with (red) Eq. (B1), (blue) Eq. (C1), and (green) Eq. (11). In (d), red and blue shaded background represent bunched and antibunched statistics, respectively. The insets illustrate the different laser detunings considered. The molecular detuning is fixed at $\delta = 10\gamma_0$ and the dipole-dipole coupling at $V = -20\gamma_0$, i.e., the modification of $\alpha$ is accompanied by a modification of the distance between the emitters according to Eq. (3).

weak enough laser intensities) to bunched emission (for increasing laser intensities), is maintained when $\alpha$ is changed.

Moreover, we plot in Fig. 9(d) cuts of the results in Figs. 9(a)–9(c) at $I = 10^{-3}I_{sat}$. $g^{(2)}(0)$ is almost constant when the laser is tuned to the two-photon resonance (red-solid line) and increases slightly with $\alpha$ [from $g^{(2)}(0) \approx 1.5 \times 10^{-4}$ when $\alpha \to 0$ to $g^{(2)}(0) \approx 6 \times 10^{-4}$ at $\alpha = 1$] if the laser is tuned resonantly to the superradiant state $|\Lambda_\rangle$ (blue-solid line). More significantly, when the laser is tuned to the transition frequency of the subradiant state $|\Lambda_+\rangle$ (green-solid line), $g^{(2)}(0)$ decreases considerably with $\alpha$ [from $\approx 0.4$ when $\alpha \to 0$ up to $\approx \gamma_0^2/(2\Lambda)^2 \approx 6 \times 10^{-4}$ at $\alpha = 1$]. The analytical

results [dots in Fig. 9(d)] obtained with (red) Eq. (B1), (blue) Eq. (C1), and (green) Eq. (11) show an excellent agreement with the numerical results for all values of $\alpha$.

Furthermore, we analyze in Fig. 10 the influence of $\alpha$ on the FRIC maps at $I = 0.1I_{sat}$ and for a laser tuned the transition frequency of the subradiant state $|\Lambda_\rangle$. We choose this laser frequency because of the large sensitivity of the color-blind intensity correlation to $\alpha$ found in Fig. 9. Here, we fix again $V = -20\gamma_0$ and, thus, changes in $\alpha$ only modify the decay rates $\gamma_\pm$. Figure 10(a) shows the FRIC map at $\alpha = 0.05$, which is very similar to the FRIC map obtained at $\alpha = 0.3$ in Fig. 10(b),
<figure><img src="image_15.png" /><figcaption>FIG. 10. Influence of α on the FRIC maps when the laser is tuned to the transition frequency of the subradiant state |Λ+⟩. The inset in (a) illustrates the laser frequency. The combined Debye-Waller/Franck-Condon factor is fixed at (a) α = 0.05, (b) α = 0.3, (c) α = 0.75, and (d) α = 1. We follow the standard color convention in the FRIC maps. According to this convention, blue color means antibunched emission, red bunched emission, and white Poissonian statistics. We use a logarithmic scale in the color bar in the complete range of values of g^{(2)}(ω₁, ω₂; 0). We fix V = -20γ₀, I = 0.1Iₐₐₐ, and Γ = 0.1γ₀.</figcaption></figure>

corresponding to the FRIC map in Fig. 5(a), discussed in Sec. IV C. This large similarity is attributed to the small difference between the γ−/γ+ ratio in both cases (γ−/γ+ ≈ 1.1 when α = 0.05, while γ−/γ+ ≈ 1.8 when α = 0.3). However, Fig. 10(c) shows that the FRIC map is substantially modified at α = 0.75 (γ−/γ+ ≈ 6.3), with a more strongly antibunched background and larger values of intensity correlation in the bunched antidiagonal lines. Most notably, neglecting the combined Debye-Waller/Franck-Condon factor (α = 1), as usual in the literature, leads to a drastically different FRIC map [Fig. 10(d)]. In this case γ−/γ+ ≈ 66 and the FRIC map exhibits a single antidiagonal line (which is not bunched) and a single horizontal/vertical line.

## APPENDIX E: WEAK ILLUMINATION APPROXIMATION WHEN THE LASER IS TUNED TO THE TRANSITION FREQUENCY OF THE SUBRADIANT STATE

In this Appendix we obtain the steady-state populations under weak illumination and for a laser tuned to the transition

frequency of the subradiant state |Λ+⟩. With this aim, we calculate the populations of two independent three-level systems, one of them consisting of the states {|G⟩, |Λ+⟩, |E⟩}, corresponding to the excitation of the doubly excited state via the resonantly pumped subradiant state. The second three-level system is composed by the states {|G⟩, |Λ−⟩, |E⟩}, corresponding to the excitation of the doubly excited state via the nonresonantly driven superradiant state |Λ−⟩. The pumping rates for the transitions |G⟩ → |Λ±⟩ and |Λ±⟩ → |E⟩ are Ω± and the decay rates of the transitions |E⟩ → |Λ±⟩ and |Λ±⟩ → |G⟩ are γ±.

The dynamics equations of the elements of the density matrix for the first three level-system {|G⟩, |Λ+⟩, |E⟩} are given by

$$ \frac{d}{dt} \langle \hat{\sigma}_{EE} \rangle = -\gamma_{+} \langle \hat{\sigma}_{EE} \rangle + i \frac{\Omega_{+}}{2} \langle \hat{\sigma}_{E+} - \hat{\sigma}_{+E} \rangle, \quad (E1a) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{++} \rangle = -\gamma_{+} \langle \hat{\sigma}_{++} - \hat{\sigma}_{EE} \rangle + i \frac{\Omega_{+}}{2} \langle \hat{\sigma}_{+E} - \hat{\sigma}_{E+} + \hat{\sigma}_{+G} - \hat{\sigma}_{G+} \rangle, \quad (E1b) $$
$$ \frac{d}{dt} \langle \hat{\sigma}_{+G} \rangle = -i \frac{\Omega_+}{2} \langle -\hat{\sigma}_{++} + \hat{\sigma}_{GG} + \hat{\sigma}_{EG} \rangle - \frac{\gamma_+}{2} \langle \hat{\sigma}_{+G} \rangle, \qquad (E1c) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{E+} \rangle = i \frac{\Omega_+}{2} \langle \hat{\sigma}_{EE} - \hat{\sigma}_{++} + \hat{\sigma}_{EG} \rangle - (i2\Lambda + \gamma_+/2) \langle \hat{\sigma}_{E+} \rangle, \qquad (E1d) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{EG} \rangle = i \frac{\Omega_+}{2} \langle \hat{\sigma}_{E+} - \hat{\sigma}_{+G} \rangle - (i2\Lambda + \gamma_+/2) \langle \hat{\sigma}_{EG} \rangle. \qquad (E1e) $$

Next, by neglecting the steady-state coherence $\rho_{EG} = \langle E|\hat{\rho}_{ss}|G \rangle$ between the ground state $|G\rangle$ and the doubly excited state $|E\rangle$ we solve the steady-state density matrix. The resulting steady-state populations can be additionally simplified assuming $\Lambda \gg \gamma_0$, which leads to the steady-state population of the subradiant state $|\Lambda_+\rangle$,

$$ \rho_{\Lambda_+}^{|\Lambda_+\rangle} = \frac{1 + (\Omega_+/4\Lambda)^2}{2 + 3(\Omega_+/4\Lambda)^2 + (\gamma_+/\Omega_+)^2}, \qquad (E2) $$

and to the steady-state population of the doubly excited state $|E\rangle$,

$$ \rho_E^{(+)}|^{\Lambda_+} = \frac{(\Omega_+/4\Lambda)^2}{2 + 3(\Omega_+/4\Lambda)^2 + (\gamma_+/\Omega_+)^2}. \qquad (E3) $$

We follow the same procedure with the three-level system $\{|G\rangle, |\Lambda_-\rangle, |E\rangle\}$. The transition $|G\rangle \to |\Lambda_-\rangle$ is nonresonantly pumped (with detuning $2\Lambda$), whereas the transition $|\Lambda_+\rangle \to |E\rangle$ is resonantly driven. The dynamics equations are then given by

$$ \frac{d}{dt} \langle \hat{\sigma}_{EE} \rangle = -\gamma_- \langle \hat{\sigma}_{EE} \rangle + i \frac{\Omega_-}{2} \langle \hat{\sigma}_{E-} - \hat{\sigma}_{-E} \rangle, \qquad (E4a) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{--} \rangle = -\gamma_- \langle \hat{\sigma}_{--} - \hat{\sigma}_{EE} \rangle + i \frac{\Omega_-}{2} \langle \hat{\sigma}_{-E} - \hat{\sigma}_{E-} + \hat{\sigma}_{-G} - \hat{\sigma}_{G-} \rangle, \qquad (E4b) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{-G} \rangle = -i \frac{\Omega_-}{2} \langle -\hat{\sigma}_{--} + \hat{\sigma}_{GG} + \hat{\sigma}_{EG} \rangle - (i2\Lambda + \gamma_-/2) \langle \hat{\sigma}_{-G} \rangle, \qquad (E4c) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{E-} \rangle = i \frac{\Omega_-}{2} \langle \hat{\sigma}_{EE} - \hat{\sigma}_{--} + \hat{\sigma}_{EG} \rangle - \frac{\gamma_-}{2} \langle \hat{\sigma}_{E-} \rangle, \qquad (E4d) $$

$$ \frac{d}{dt} \langle \hat{\sigma}_{EG} \rangle = i \frac{\Omega_-}{2} \langle \hat{\sigma}_{E-} - \hat{\sigma}_{-G} \rangle - (i2\Lambda + \gamma_-/2) \langle \hat{\sigma}_{EG} \rangle. \qquad (E4e) $$

In this case, we do not neglect any coherence and only assume that $\Lambda \gg \gamma_0$. The steady-state population of the subradiant state $|\Lambda_-\rangle$ is thus given by

$$ \rho_{\Lambda-}^{|\Lambda_+\rangle} = \frac{1}{1 + (4\Lambda/\Omega_-)^2}, \qquad (E5) $$

and the population of the doubly excited state $|E\rangle$ by

$$ \rho_E^{(-)}|^{\Lambda_+} = \frac{2(\Omega_-/4\Lambda)^4}{2 + 3(\Omega_-/4\Lambda)^2 + (\Omega_-/4\Lambda)^4}. \qquad (E6) $$

Therefore, the total population $\rho_E$ of the doubly excited state $|E\rangle$, understood as the sum of the population $\rho_E^{(+)}$ excited

via the resonantly driven subradiant state $|\Lambda_+\rangle$ and the population $\rho_E^{(-)}$ excited via the nonresonantly driven superradiant state $|\Lambda_-\rangle$ becomes

$$ \rho_E^{|\Lambda_+\rangle} = \rho_E^{(+)}|^{\Lambda_+\rangle} + \rho_E^{(-)}|^{\Lambda_+\rangle}. \qquad (E7) $$

Finally, according to the expression of the intensity correlation $g^{(2)}(0)|_{I \le I_{sat}}$ in Eq. (10), which is valid under weak illumination, we obtain the intensity correlation for this laser detuning given by Eq. (11).

## APPENDIX F: FRIC OF A SINGLE TLS

In this Appendix we summarize the FRIC of a single TLS that is coherently and resonantly driven by a laser [34,69]. The Hamiltonian in this case is given as $\hat{H} = \hbar \frac{\omega_0}{2} \hat{\sigma}_z - \hbar \frac{\Omega_2}{2} (\hat{\sigma}^\dagger e^{-i\omega_0 t} + \hat{\sigma} e^{i\omega_0 t})$ and the spontaneous emission rate $\gamma_0$ of the emitter is incorporated through the Lindblad dissipator $\frac{\gamma_0}{2} \mathcal{D}[\hat{\sigma}]$.

We present in Fig. 11(a) the FRIC map for the case of laser intensity $I = 0.1 I_{sat}$. For such weak driving the eigenstates of the system are the ground and excited states $|g\rangle$ and $|e\rangle$ and thus the emission spectrum [shown on top of Fig. 11(a)] shows a single resonance at frequency $\omega_0$. Here, we normalize the emission spectrum $S(\omega) \propto \int_{-\infty}^{\infty} d\tau \langle \delta\hat{E}^{(-)}(\tau) \delta\hat{E}^{(+)}(0) \rangle_{ss} e^{-i\omega_0 t}$ [with $\delta\hat{E}^{(-)}(\tau) = \sigma(t) - \langle\sigma\rangle_{ss}$], in such a way that $\int_0^{\infty} d\omega S(\omega) = \gamma_0 \langle \delta\hat{\sigma}^\dagger \delta\hat{\sigma} \rangle_{ss}$. Within the picture where the laser field is quantized, the two eigenstates $|e, n\rangle$ and $|g, n+1\rangle$ are degenerated in the rung $n$ under weak illumination. Thus the one-photon transitions between eigenstates in rung $n+1$ and rung $n$ give raise to the single emission peak at $\omega_{las} = \omega_0$, as shown by $S(\omega)$. This is translated into the FRIC with the emergence of a single vertical line at $\omega_1 = \omega_{las}$ and a horizontal line at $\omega_2 = \omega_{las}$. Additionally, the leapfrog processes necessarily fulfill $\omega_1 + \omega_2 = 2\omega_{las}$ and give raise to the single antidiagonal line that emerges in the FRIC map.

For increasing laser intensities, the two eigenstates in the same rung are no longer degenerate. The two eigenstates in the rung $n$ become the dressed states $|S_{eg}^{(n)}\rangle = (|e, n-1\rangle + |g, n\rangle)/\sqrt{2}$ and $|A_{eg}^{(n)}\rangle = (|e, n-1\rangle - |g, n\rangle)/\sqrt{2}$, with energy splitting $\hbar\Omega$. We plot in Fig. 11(b) the dressed states in rungs $n-1, n$, and $n+1$ and the different one-photon and two-photon transitions. In addition to the photons emitted at $\omega_{las} = \omega_0$ [gray-solid arrows in Fig. 11(b)], photons of frequency $\omega_{las} - \Omega = \omega_0 - \Omega$ and $\omega_{las} + \Omega = \omega_0 + \Omega$ can also be emitted due to the one-photon transitions $|S_{eg}^{(n+1)}\rangle \to |A_{eg}^{(n)}\rangle$ [blue-solid arrow in Fig. 11(b)] and $|A_{eg}^{(n+1)}\rangle \to |S_{eg}^{(n)}\rangle$ [green-solid arrow], respectively. These new one-photon transitions lead to the emergence of additional vertical and horizontal lines in the FRIC map and to the emergence of two additional peaks in the emission spectrum (Mollow triplet) as the intensity of the incident laser increases. Similarly, there are now two additional leapfrog processes, one with $|S_{eg}^{(n+1)}\rangle$ as initial state in rung $n+1$ and with $|A_{eg}^{(n+1)}\rangle$ as final state in rung $n-1$ [which fulfills $\omega_1 + \omega_2 = 2\omega_{las} + |\Omega|$ and is represented with a dashed blue arrow in Fig. 11(b)] and another one with $|A_{eg}^{(n+1)}\rangle$ as initial state in rung $n+1$ and with $|S_{eg}^{(n+1)}\rangle$ as final state in rung $n-1$.
<figure><img src="image_16.png" /><figcaption>FIG. 11. FRIC for a resonantly driven TLS. [(a),(c),(d)] Normalized emission spectra (top panel) and FRIC maps (bottom panel) for laser intensity (a) $I = 0.1I_{sat}$, (c) $I = 50I_{sat}$, and (d) $I = 1000I_{sat}$. We illustrate in (b) the dressed states in rungs $n-1$, $n$, and $n+1$, together with the corresponding one-photon transitions between eigenstates (solid arrows) and two-photon transitions through intermediate virtual states (dashed arrows). We follow the standard color convention in the FRIC maps, i.e., blue color means antibunched emission, while red color bunched emission, and white Poissonian statistics. In the color bar, we utilize a linear scale in the interval $0 \le g^{(2)}(\omega_1, \omega_2; 0) \le 1$ and a logarithmic scale in the range $1 \le g^{(2)}(\omega_1, \omega_2; 0) \le 100$. Solid arrows are used to mark the horizontal and vertical lines in the FRIC maps, while dashed arrows are used for the antidiagonal lines related to the two-photon processes through intermediate virtual states. We use a filter linewidth $\Gamma = \gamma_0/10$ and fix the values $V = -20\gamma_0$, $\alpha = 0.3$, and $\delta = 10\gamma_0$.</figcaption></figure>

state in rung $n-1$ [which fulfills $\omega_1 + \omega_2 = 2\omega_{las} - |\Omega|$ and is represented with a dashed-green arrow in Fig. 11(b)].

We represent in Fig. 11(c) the FRIC map for laser intensity $I = 50I_{sat}$. For this intensity, we observe the Mollow triplet in the emission spectrum on top of Fig. 11(c); however, in the FRIC map we observe a broad horizontal line and a broad vertical line due to the overlapping of the three horizontal and vertical lines. Additionally, we also observe that the antidiagonal line is much broader than in Fig. 11(a), due to the actual presence of three antidiagonal lines, which are very close to each other. Finally, we present in Fig. 11(d) the FRIC

map for $I = 1000I_{sat}$. For this laser intensity we observe three well-differentiated vertical, horizontal, and antidiagonal lines, since the laser intensity is strong enough such that the energy splitting between the dressed states $\hbar\Omega$ is much larger than $\hbar\gamma_0$.

## APPENDIX G: HAMILTONIAN DIAGONALIZATION WHEN THE LASER IS TUNED TO THE SUPERRADIANT OR TO THE SUBRADIANT STATE

In this Appendix we show how to perform an approximate diagonalization of the Hamiltonian in Eq. (5) when the laser
<figure><img src="image_17.png" /><figcaption>FIG. 12. (a) Numerical eigenenergies $E_i$ (normalized by $\Lambda$) of the Hamiltonian in Eq. (5) when the laser is tuned either to the transition frequency of the superradiant state or to the transition frequency of the subradiant state ($E_i$ are identical for both laser frequencies). (b) Numerical energy splitting (normalized by $|\Omega|)$ between the more energetic eigenstates $|E_3\rangle$ and $|E_4\rangle$ (purple line) and energy splitting between the less energetic eigenstates $|E_1\rangle$ and $|E_2\rangle$ (brown line). Horizontal-black-dashed lines represent the same energy splittings under the analytical approximation, see Eq. (G1). (c) Deviation of the analytical eigenstates $|A_{Y-}\rangle$, $|S_{Y-}\rangle$, $|S_{X+}\rangle$, $|A_{X+}\rangle$ from the exact eigenstates $|E_1\rangle$, $|E_2\rangle$, $|E_3\rangle$, $|E_4\rangle$ obtained numerically. The deviation is estimated as the difference between 1 and the squared overlap of the analytical eigenstate and the corresponding numerical eigenstate. For the case in which the laser is tuned to the transition frequency of the superradiant state then $X = E$ and $Y = G$, whereas if the laser is tuned to the transition frequency of the superradiant state then $X = G$ and $Y = E$. The molecular detuning is fixed as $\delta = 10\gamma_0$, the dipole-dipole coupling as $V = -20\gamma_0$ and the combined Debye-Waller/Franck-Condon factor $\alpha = 0.3$.</figcaption></figure>

is tuned either to the transition frequency of the superradiant state or to the transition frequency of the subradiant state. For both laser detunings and very weak or no driving, there are two pairs of degenerate eigenstates at energies $\pm\Lambda$ in the rotating frame at the laser frequency. Namely, if the laser is tuned to the transition frequency of the superradiant state $|\Lambda_-\rangle$ ($\Delta_0 = \Lambda$), the states $|G\rangle$ and $|\Lambda_-\rangle$ become degenerate (with energy $-\Lambda$) and the states $|E\rangle$ and $|\Lambda_+\rangle$ become degenerate (with energy $\Lambda$). In contrast, if the laser is tuned to the transition frequency of the subradiant state $|\Lambda_+\rangle$ ($\Delta_0 = -\Lambda$), the states $|E\rangle$ and $|\Lambda_-\rangle$ become degenerate (with energy $-\Lambda$) and the states $|G\rangle$ and $|\Lambda_+\rangle$ become degenerate (with energy $\Lambda$).

The approximation that we perform consists in neglecting the driving terms in the second line of Eq. (5) that connect states that are not degenerated for no driving. As a consequence, for a laser tuned to the transition frequency of the superradiant state $|\Lambda_-\rangle$, the driving Hamiltonian becomes $\hat{H}_{las} = -\hbar\frac{\Omega_+}{2}\hat{\sigma}_{E_+} - \hbar\frac{\Omega_-}{2}\hat{\sigma}_{G_-} + \text{H.c.}$, whereas if the laser is tuned to the transition frequency of the subradiant state $|\Lambda_+\rangle$ then the laser Hamiltonian is approximated as $\hat{H}_{las} = -\hbar\frac{\Omega_+}{2}\hat{\sigma}_{G_+} - \hbar\frac{\Omega_-}{2}\hat{\sigma}_{E_-} + \text{H.c.}$. The resulting total Hamiltonians can be analytically diagonalized.

For a laser tuned to the transition frequency of the superradiant state, the approximation leads to the eigenstates $|A_{E+}\rangle = (|E\rangle - |\Lambda_+\rangle)/\sqrt{2}$, $|S_{E+}\rangle = (|E\rangle + |\Lambda_+\rangle)/\sqrt{2}$, $|S_{G-}\rangle = (|G\rangle + |\Lambda_-\rangle)/\sqrt{2}$, and $|A_{G-}\rangle = (|G\rangle - |\Lambda_-\rangle)/\sqrt{2}$, while for a laser tuned to the transition frequency of the subradiant state, the approximation leads to the eigenstates $|A_{G+}\rangle = (|G\rangle - |\Lambda_+\rangle)/\sqrt{2}$, $|S_{G+}\rangle = (|G\rangle + |\Lambda_+\rangle)/\sqrt{2}$, $|S_{E-}\rangle = (|E\rangle + |\Lambda_-\rangle)/\sqrt{2}$, and $|A_{E-}\rangle = (|E\rangle - |\Lambda_-\rangle)/\sqrt{2}$. The corresponding eigenenergies are equal to

$$E(|A_{X+}\rangle) = \Lambda + \frac{\Omega}{2}(\cos \theta + \sin \theta), \quad (G1a)$$

$$E(|S_{X+}\rangle) = \Lambda - \frac{\Omega}{2}(\cos \theta + \sin \theta), \quad (G1b)$$

$$E(|S_{Y-}\rangle) = -\Lambda - \frac{\Omega}{2}(\cos \theta - \sin \theta), \quad (G1c)$$

$$E(|A_{Y-}\rangle) = -\Lambda + \frac{\Omega}{2}(\cos \theta - \sin \theta), \quad (G1d)$$

where $X = E$ and $Y = G$ if the laser is tuned to the transition frequency of the superradiant state $|\Lambda_-\rangle$ and $X = G$ and $Y = E$ if the laser is tuned to the transition frequency of the subradiant state $|\Lambda_+\rangle$. We note that the energy splitting between the two higher-energy states $|A_{X+}\rangle$ and $|S_{X+}\rangle$ is equal to $|\Omega(\cos \theta + \sin \theta)|$, while the energy splitting between the lower-energy states $|S_{Y-}\rangle$ and $|A_{Y-}\rangle$ is equal to $|\Omega(\cos \theta - \sin \theta)|$. We next compare the eigenenergies obtained analytically in Eqs. (G1), with the numerical eigenenergies $E_i$ ($i \in \{1, 2, 3, 4\}$) of the exact numerical eigenstates $|E_i\rangle$. We plot in Fig. 12(a) the exact numerical eigenenergies $E_i$ (normalized by $\Lambda$) and in Fig. 12(b) the numerical energy splitting (normalized by the absolute value of the Rabi frequency $|\Omega|)$ between the two higher-energy states (purple line), as well as between the two lower-energy states (brown line). The numerical energy splittings in Fig. 12(b) show an excellent agreement with the analytical energy splittings calculated from Eqs. (G1) (horizontal-dashed-black lines) up to a laser intensity $I \approx 100I_{sat}$.

Finally, we analyze the accuracy of the eigenstates obtained analytically with the approximation presented in this Appendix. With this aim, we plot in Fig. 12(c) the difference between 1 and the squared overlap of the numerical and analytical eigenstates, which serves as a measure of the error of the analytical eigenstates. For example, for the lowest energetic state, the deviation of the eigenstate $|A_{Y-}\rangle$ (obtained analytically) from the exact eigenstate $|E_1\rangle$ (calculated numerically) is given by $1 - |\langle A_{Y-}|E_1\rangle|^2$. For weak and moderate laser intensities the analytical eigenstates are found to be a very accurate approximation of the exact eigenstates of the system.
<figure><img src="image_18.png" /><figcaption>FIG. 13. Dependence of the eigenenergies on the laser intensity for a laser tuned to the two-photon resonance ($2\omega_{las} = 2\omega_0$). Solid lines represent the numerical eigenenergies (normalized by $\Lambda$) of the full Hamiltonian in Eq. (5), whereas dots represent the analytical approximation of the eigenenergies given by Eq. (H1). The molecular detuning is fixed as $\delta = 10\gamma_0$, the dipole-dipole coupling as $V = -20\gamma_0$ and the combined Debye-Waller/Franck-Condon factor $\alpha = 0.3$.</figcaption></figure>

## APPENDIX H: FRIC WHEN THE LASER IS TUNED TO THE TWO-PHOTON RESONANCE

In this Appendix we show the analytical approximation that we perform to diagonalize the Hamiltonian in Eq. (5) for a laser tuned to the two-photon resonance ($2\omega_{las} = 2\omega_0$) and arbitrary laser intensity. We then present the FRIC maps for this detuning and the same laser intensities than for the other laser detunings in Sec. IV, again for a fixed value of molecular detuning $\delta = 10\gamma_0$.

In order to diagonalize the Hamiltonian for a laser tuned to the two-photon resonance, we neglect in Eq. (5) the terms that are proportional to $\Omega_+$. The reason is that both the superradiant state $|\Lambda_-\rangle$ and the subradiant state $|\Lambda_+\rangle$ are non-resonantly driven and, additionally, the effective pumping $\Omega_+$ through the subradiant state is one order of magnitude smaller than the effective pumping $\Omega_-$ through the superradiant state. Neglecting these terms we can analytically diagonalize the Hamiltonian in Eq. (5), which leads to the eigenenergies

$$E_1 = -\frac{\Lambda}{2} - \sqrt{\left(\frac{\Lambda}{2}\right)^2 + \left(\frac{\Omega_-}{\sqrt{2}}\right)^2}, \quad (H1a)$$

$$E_2 = 0, \quad (H1b)$$

$$E_3 = -\frac{\Lambda}{2} + \sqrt{\left(\frac{\Lambda}{2}\right)^2 + \left(\frac{\Omega_-}{\sqrt{2}}\right)^2}, \quad (H1c)$$

$$E_4 = \Lambda. \quad (H1d)$$

We plot in Fig. 13 the analytical (dots) and numerical (solid lines) eigenenergies for $\delta = 10\gamma_0$ showing that they are in

very good agreement. The corresponding analytical eigenstates result in

$$|E_1\rangle = \frac{1}{\sqrt{\left(\frac{\Omega_-}{\sqrt{2}}\right)^2 + E_3^2}} \left( \frac{\Omega_-}{\sqrt{2}} |\Lambda_-\rangle + E_3 |S_{EG}\rangle \right), \quad (H2a)$$

$$|E_2\rangle = |A_{EG}\rangle, \quad (H2b)$$

$$|E_3\rangle = \frac{1}{\sqrt{\left(\frac{\Omega_-}{\sqrt{2}}\right)^2 + E_1^2}} \left( -\frac{\Omega_-}{\sqrt{2}} |\Lambda_-\rangle - E_1 |S_{EG}\rangle \right), \quad (H2c)$$

$$|E_4\rangle = |\Lambda_+\rangle, \quad (H2d)$$

where $|S_{EG}\rangle = (|E\rangle + |G\rangle)/\sqrt{2}$ and $|A_{EG}\rangle = (|E\rangle - |G\rangle)/\sqrt{2}$. Equation (H2) shows that $|A_{EG}\rangle$ and $|\Lambda_+\rangle$ are eigenstates of the system for any laser intensity, whereas the other two dressed eigenstates ($|E_1\rangle$ and $|E_3\rangle$) are superpositions of $|S_{EG}\rangle$ and $|\Lambda_-\rangle$, with coefficients that depend on the laser intensity. Additionally, we note that under weak enough illumination, the eigenstates $|E_2\rangle = |A_{EG}\rangle$ and $|E_3\rangle = |S_{EG}\rangle$ become degenerate in energy.

Next, we plot in Fig. 14(a) the emission spectrum and the FRIC map for $I = 0.1I_{sat}$. This map is in general strongly bunched, which can be expected from the strongly bunched color-blind intensity correlation for these laser parameters [see Fig. 3(a)]. Over this strongly bunched background, we can distinguish two antibunched vertical lines at $\omega_1 = \omega_{las}$ (here, $\omega_{las} = \omega_0$) and $\omega_1 = \omega_{las} - \Lambda$ and two horizontal lines at $\omega_2 = \omega_{las}$ and $\omega_2 = \omega_{las} - \Lambda$ (marked by grey and light green solid arrows) in the FRIC map. These two lines are related to the one-photon transitions depicted in the level scheme in Fig. 14(b). Further, we observe in the emission spectrum on top of Fig. 14(a) a third peak at frequency $\omega_1 = \omega_{las} + \Lambda$, related to the one-photon transitions from the subradiant state $|\Lambda_+, n\rangle$ in rung $n+1$ to the ground state $|G, n\rangle$ in rung $n$, whose corresponding horizontal and vertical lines cannot be distinguished in the FRIC map due to the strong background bunching. In contrast, three antidiagonal lines with stronger bunching than the background emerge in the FRIC map in Fig. 14(a) (see the broad line marked by a dashed gray arrow and the subtle lines marked by a dashed-light-green arrow and a dashed-dark-green arrow) corresponding to three different type of leapfrog processes that can take place for these laser parameters.

As the laser intensity is increased up to $I = 50I_{sat}$, the eigenstates of the system become slightly dressed, according to Eq. (H2). Figure 14(c) shows the corresponding spectrum and the FRIC map for $I = 50I_{sat}$, and Fig. 14(d) the energy level structure. Notably, as the value of the color-blind intensity correlation $g^{(2)}(0)$ decreases, then the value of the frequency-resolved intensity correlation $g^{(2)}(\omega_1, \omega_2; 0)$ of the background of the FRIC map also decreases and the three horizontal and three vertical lines can now be clearly distinguished. This map is very similar to the one discussed in Ref. [41] for identical emitters driven by a laser tuned to the two-photon resonance. For even stronger laser intensities, the eigenstates of the system become strongly dressed and a variety of horizontal and vertical lines emerge in the FRIC
<figure><img src="image_19.png" /><figcaption>FIG. 14. FRIC for a laser tuned to the two-photon resonance, $2\omega_{las} = 2\omega_0$. [(a),(c),(e)] Normalized emission spectra (top panel) and FRIC maps (bottom panel) for laser intensity (a) $I = 0.1I_{sat}$, (c) $I = 50I_{sat}$ and (e) $I = 1000I_{sat}$. The inset in (a) illustrates the molecular energy levels and the resonant excitation of the doubly excited state $|E\rangle$ with two photons. In the FRIC maps we use the standard color convention, where blue color represents antibunched emission, while red color bunched emission and white Poissonian statistics. In the color bar, we use a linear scale in the interval $0 \le g^{(2)}(\omega_1, \omega_2; 0) \le 1$ and a logarithmic scale in the range $g^{(2)}(\omega_1, \omega_2; 0) \ge 1$. We consider a filter linewidth $\Gamma = \gamma_0/10$. [(b),(d),(f)] Schematic representation of the eigenstates (accounting for the quantum nature of the laser field) in rungs $n-1, n$ and $n+1$ for laser intensity (b) $I = 0.1I_{sat}$, (d) $I = 50I_{sat}$, and (f) $I = 1000I_{sat}$. Each rung contains four eigenstates with equal number of excitations. Solid arrows are used to mark the horizontal and vertical lines in the FRIC maps and the one-photon transitions in the diagrams of the eigenstates, whereas dashed arrows are used for the antidiagonal lines related to the two-photon processes through virtual states. The dipole-dipole coupling is fixed as $V = -20\gamma_0$, the combined Debye-Waller/Franck-Condon factor as $\alpha = 0.3$ and the molecular detuning as $\delta = 10\gamma_0$.</figcaption></figure>

map [41]. We plot in Fig. 14(e) the emission spectrum and the FRIC map for $I = 1000I_{sat}$, with the corresponding one-

photon transitions as well as leapfrog processes illustrated in Fig. 14(f).

[1] B. Lounis and M. Orrit, Single-photon sources, Rep. Prog. Phys. **68**, 1129 (2005).

[2] C. H. Bennett and G. Brassard, Quantum cryptography: Public key distribution and coin tossing, in *Proceedings of the IEEE International Conference on Computers, Systems and Signal Processing, 1984* (IEEE, Piscataway, NJ, 1984), pp. 175-179.

[3] C. H. Bennett and G. Brassard, *Advances in Cryptology* (Springer, Berlin, 1985), pp. 475-480.

[4] C. H. Bennett, F. Bessette, G. Brassard, L. Salvail, and J. Smolin, Experimental quantum cryptography, *J. Cryptology* **5**, 3 (1992).

[5] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, Teleporting an unknown quantum state via dual classical and Einstein-Podolsky-Rosen channels, *Phys. Rev. Lett.* **70**, 1895 (1993).

[6] D. Fattal, E. Diamanti, K. Inoue, and Y. Yamamoto, Quantum teleportation with a quantum dot single photon source, *Phys. Rev. Lett.* **92**, 037904 (2004).

[7] E. Knill, R. Laflamme, and G. J. Milburn, A scheme for efficient quantum computation with linear optics, *Nature (London)* **409**, 46 (2001).

[8] T. C. Ralph, A. G. White, W. J. Munro, and G. J. Milburn, Simple scheme for efficient linear optics quantum gates, *Phys. Rev. A* **65**, 012314 (2001).

[9] C. Brunel, B. Lounis, P. Tamarat, and M. Orrit, Triggered source of single photons based on controlled single molecule fluorescence, *Phys. Rev. Lett.* **83**, 2722 (1999).

[10] C. Toninelli, I. Gerhardt, A. S. Clark, A. Reserbat-Plantey, S. Götzinger, Z. Ristanović, M. Colautti, P. Lombardi, K. D. Major, I. Deperasińska *et al.*, Single organic molecules for photonic quantum technologies, *Nat. Mater.* **20**, 1615 (2021).
[11] J. Wrachtrup, C. Von Borczyskowski, J. Bernard, M. Orrit, and R. Brown, Optical detection of magnetic resonance in a single molecule, *Nature* (London) **363**, 244 (1993).

[12] S. Pazzagli, P. Lombardi, D. Martella, M. Colautti, B. Tiribilli, F. S. Cataliotti, and C. Toninelli, Self-assembled nanocrystals of polycyclic aromatic hydrocarbons show photostable single-photon emission, *ACS Nano* **12**, 4295 (2018).

[13] B. Lounis, H. A. Bechtel, D. Gerion, P. Alivisatos, and W. E. Moerner, Photon antibunching in single CdSe/ZnS quantum dot fluorescence, *Chem. Phys. Lett.* **329**, 399 (2000).

[14] A. Zrenner, E. Beham, S. Stufler, F. Findeis, M. Bichler, and G. Abstreiter, Coherent properties of a two-level system based on a quantum-dot photodiode, *Nature* (London) **418**, 612 (2002).

[15] S. Stufler, P. Ester, A. Zrenner, and M. Bichler, Quantum optical properties of a single $In_xGa_{1-x}As-GaAs$ quantum dot two-level system, *Phys. Rev. B* **72**, 121301(R) (2005).

[16] P. Senellart, G. Solomon, and A. White, High-performance semiconductor quantum-dot single-photon sources, *Nat. Nanotechnol.* **12**, 1026 (2017).

[17] Y. Arakawa and M. J. Holmes, Progress in quantum-dot single photon sources for quantum information technologies: A broad spectrum overview, *Appl. Phys. Rev.* **7**, 021309 (2020).

[18] M. Gross and S. Haroche, Superradiance: An essay on the theory of collective spontaneous emission, *Phys. Rep.* **93**, 301 (1982).

[19] J. B. Trebbia, Q. Deplano, P. Tamarat, and B. Lounis, Tailoring the superradiant and subradiant nature of two coherently coupled quantum emitters, *Nat. Commun.* **13**, 2962 (2022).

[20] C. M. Lange, E. Daggett, V. Walther, L. Huang, and J. D. Hood, Superradiant and subradiant states in lifetime-limited organic molecules through laser-induced tuning, *Nat. Phys.* (2024).

[21] A. Tiranov, V. Angelopoulou, C. J. van Diepen, B. Schrinski, O. A. D. Sandberg, Y. Wang, L. Midolo, S. Scholz, A. D. Wieck, A. Ludwig et al., Collective super-and subradiant dynamics between distant optical quantum emitters, *Science* **379**, 389 (2023).

[22] X. L. Chu, C. Papon, N. Bart, A. D. Wieck, A. Ludwig, L. Midolo, N. Rotenberg, and P. Lodahl, Independent electrical control of two quantum dots coupled through a photonic-crystal waveguide, *Phys. Rev. Lett.* **131**, 033606 (2023).

[23] D. Rattenbacher, A. Shkarin, J. Renger, T. Utikal, S. Götzinger, and V. Sandoghdar, On-chip interference of scattering from two individual molecules, *Optica* **10**, 1595 (2023).

[24] L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge University Press, Cambridge, 1995).

[25] R. E. K. Fishman, R. N. Patel, D. A. Hopper, T. Y. Huang, and L. C. Bassett, Photon-emission-correlation spectroscopy as an analytical tool for solid-state quantum defects, *PRX Quantum* **4**, 010202 (2023).

[26] R. Hanbury brown and R. Q. Twiss, A test of a new type of stellar interferometer on Sirius, *Nature* (London) **178**, 1046 (1956).

[27] C. Hettich, C. Schmitt, J. Zitzmann, S. Kuhn, I. Gerhardt, and V. Sandoghdar, Nanometer resolution and coherent optical dipole coupling of two individual molecules, *Science* **298**, 385 (2002).

[28] R. D. Griffin and S. M. Harris, Two-atom resonance fluorescence including the dipole-dipole interaction, *Phys. Rev. A* **25**, 1528 (1982).

[29] Z. Ficek, R. Tanas, and S. Kielich, Effect of interatomic interactions on resonance fluorescence of two atoms coherently

driven by a strong resonant laser field, *Opt. Acta* **30**, 713 (1983).

[30] S. V. Lawande, B. N. Jagatap, and Q. V. Lawande, Inhibition of fluorescence in a system of two interacting two-level atoms: A quantum jump like behaviour, *Opt. Commun.* **73**, 126 (1989).

[31] A. Beige and G. C. Hegerfeldt, Transition from antibunching to bunching for two dipole-interacting atoms, *Phys. Rev. A* **58**, 4133 (1998).

[32] J. Gillet, G. S. Agarwal, and T. Bastin, Tunable entanglement, antibunching, and saturation effects in dipole blockade, *Phys. Rev. A* **81**, 013837 (2010).

[33] A. Vivas-Viaña and C. Sánchez Muñoz, Two-photon resonance fluorescence of two interacting nonidentical quantum emitters, *Phys. Rev. Res.* **3**, 033136 (2021).

[34] A. Gonzalez-Tudela, F. P. Laussy, C. Tejedor, M. J. Hartmann, and E. Del Valle, Two-photon spectra of quantum emitters, *New J. Phys.* **15**, 033036 (2013).

[35] J. D. Cresser, Intensity correlations of frequency-filtered light fields, *J. Phys. B: At. Mol. Phys.* **20**, 4915 (1987).

[36] E. del Valle, A. Gonzalez-Tudela, F. P. Laussy, C. Tejedor, and M. J. Hartmann, Theory of frequency-filtered and time-resolved n-photon correlations, *Phys. Rev. Lett.* **109**, 183601 (2012).

[37] D. I. H. Holdaway, V. Notararigo, and A. Olaya-Castro, Perturbation approach for computing frequency-and time-resolved photon correlation functions, *Phys. Rev. A* **98**, 063828 (2018).

[38] M. Peiris, B. Petrak, K. Konthasinghe, Y. Yu, Z. C. Niu, and A. Muller, Two-color photon correlations of the light scattered by a quantum dot, *Phys. Rev. B* **91**, 195125 (2015).

[39] C. Sánchez Muñoz, F. P. Laussy, C. Tejedor, and E. Del Valle, Enhanced two-photon emission from a dressed biexciton, *New J. Phys.* **17**, 123021 (2015).

[40] M. K. Schmidt, R. Esteban, G. Giedke, J. Aizpurua, and A. González-Tudela, Frequency-resolved photon correlations in cavity optomechanics, *Quantum Sci. Technol.* **6**, 034005 (2021).

[41] E. Darsheshdar, M. Hugbart, R. Bachelard, and C. J. Villas-Boas, Photon-photon correlations from a pair of strongly coupled two-level emitters, *Phys. Rev. A* **103**, 053702 (2021).

[42] T. Basché, W. E. Moerner, M. Orrit, and U. P. e. Wild, *Single-Molecule Optical Detection, Imaging and Spectroscopy* (John Wiley & Sons, New York, 2008).

[43] G. S. Agarwal, Quantum statistical theories of spontaneous emission and their relation to other approaches, in *Quantum Optics*, edited by G. Höhler, Springer Tracts in Modern Physics, Vol. 70 (Springer, Berlin, 1974), pp. 1–128.

[44] A. Stokes and A. Nazir, A master equation for strongly interacting dipoles, *New J. Phys.* **20**, 043022 (2018).

[45] U. Akram, Z. Ficek, and S. Swain, Decoherence and coherent population transfer between two coupled systems, *Phys. Rev. A* **62**, 013413 (2000).

[46] D. A. Steck, Quantum and Atom Optics (2007), available online at http://steck.us/teaching (revision 0.8.3).

[47] H. P. Breuer and F. Petruccione, *The Theory of Open Quantum Systems* (Oxford University Press on Demand, New York, 2002).

[48] Z. Ficek and S. Swain, *Quantum Interference and Coherence: Theory and Experiments*, Vol. 100 (Springer Science & Business Media, New York, 2005).

[49] A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, and D. E. Chang, Exponential improvement in photon
storage fidelities using subradiance and “selective radiance” in atomic arrays, *Phys. Rev. X* **7**, 031024 (2017).

[50] P. O. Guimond, A. Grankin, D. V. Vasilyev, B. Vermersch, and P. Zoller, Subradiant bell states in distant atomic arrays, *Phys. Rev. Lett.* **122**, 093601 (2019).

[51] K. E. Ballantine and J. Ruostekoski, Quantum single-photon control, storage, and entanglement generation with planar atomic arrays, *PRX Quantum* **2**, 040362 (2021).

[52] R. Loudon, *The Quantum Theory of Light* (Oxford University Press, Oxford, 2000).

[53] C. Gerry, P. Knight, and P. L. Knight, *Introductory Quantum Optics* (Cambridge University Press, Cambridge, 2005).

[54] T. G. Rudolph, Z. Ficek, and B. J. Dalton, Two-atom resonance fluorescence in running-and standing-wave laser fields, *Phys. Rev. A* **52**, 636 (1995).

[55] M. O. Scully and M. S. Zubairy, *Quantum Optics* (Cambridge University Press, Cambridge, 1999).

[56] R. Loudon, Non-classical effects in the statistical properties of light, *Rep. Prog. Phys.* **43**, 913 (1980).

[57] D. Meiser and M. J. Holland, Intensity fluctuations in steady-state superradiance, *Phys. Rev. A* **81**, 063827 (2010).

[58] A. Auffèves, D. Gerace, S. Portolan, A. Drezet, and M. F. Santos, Few emitters in a cavity: From cooperative emission to individualization, *New J. Phys.* **13**, 093020 (2011).

[59] R. Jones, R. Saint, and B. Olmos, Far-field resonance fluorescence from a dipole-interacting laser-driven cold atomic gas, *J. Phys. B: At., Mol. Opt. Phys.* **50**, 014004 (2017).

[60] M. Cordier, M. Schemmer, P. Schneeweiss, J. Volz, and A. Rauschenbeutel, Tailoring photon statistics with an atom-based two-photon interferometer, *Phys. Rev. Lett.* **131**, 183601 (2023).

[61] S. Wolf, S. Richter, J. von Zanthier, and F. Schmidt-Kaler, Light of two atoms in free space: Bunching or antibunching? *Phys. Rev. Lett.* **124**, 063603 (2020).

[62] S. Ribeiro, T. F. Cutler, C. S. Adams, and S. A. Gardiner, Collective effects in the photon statistics of thermal atomic ensembles, *Phys. Rev. A* **104**, 013719 (2021).

[63] T. J. Arruda, R. Bachelard, J. Weiner, S. Slama, and P. W. Courteille, Controlling photon bunching and antibunching of two quantum emitters near a core-shell sphere, *Phys. Rev. A* **101**, 023828 (2020).

[64] C. Sánchez Muñoz and F. Schlawin, Photon correlation spectroscopy as a witness for quantum coherence, *Phys. Rev. Lett.* **124**, 203601 (2020).

[65] L. Knoll and G. Weber, Theory of n-fold time-resolved correlation spectroscopy and its application to resonance fluorescence radiation, *J. Phys. B: At. Mol. Phys.* **19**, 2817 (1986).

[66] L. Knöll, W. Vogel, and D. G. Welsch, Quantum noise in spectral filtering of light, *J. Opt. Soc. Am. B* **3**, 1315 (1986).

[67] G. Nienhuis, Spectral correlations in resonance fluorescence, *Phys. Rev. A* **47**, 510 (1993).

[68] J. C. López Carreño, E. del Valle, and F. P. Laussy, Frequency-resolved Monte Carlo, *Sci. Rep.* **8**, 1 (2018).

[69] J. C. López Carreño, E. del Valle, and F. P. Laussy, Photon correlations from the Mollow triplet, *Laser Photonics Rev.* **11**, 1700090 (2017).

[70] E. del Valle, Distilling one, two and entangled pairs of photons from a quantum dot with cavity qed effects and spectral filtering, *New J. Phys.* **15**, 025019 (2013).

[71] C. Clear, R. C. Schofield, K. D. Major, J. Iles-Smith, A. S. Clark, and D. P. S. McCutcheon, Phonon-induced optical dephasing in single organic molecules, *Phys. Rev. Lett.* **124**, 153602 (2020).

[72] M. Reitz, C. Sommer, B. Gurlek, V. Sandoghdar, D. Martin-Cano, and C. Genes, Molecule-photon interactions in phononic environments, *Phys. Rev. Res.* **2**, 033270 (2020).

[73] A. Juan-Delgado, R. Esteban, Á. Nodar, J.-B. Trebbia, B. Lounis, and J. Aizpurua, Research data supporting “Tailoring the statistics of light emitted from two interacting quantum emitter”, http://hdl.handle.net/10261/357088 (2024).

[74] Z. Ficek, R. Tanaś, and S. Kielich, Photon antibunching and squeezing in resonance fluorescence of two interacting atoms, *Phys. Rev. A* **29**, 2004 (1984).