# Supplementary Material: Fundamental Limits to Coherent Photon Generation
with Solid-State Atomlike Transitions

Z. X. Koong<EMAIL_ADDRESS>SUPA, Institute of Photonics and Quantum Sciences,
Heriot-Watt University, Edinburgh EH14 4AS, Scotland, United Kingdom  D.
Scerri  SUPA, Institute of Photonics and Quantum Sciences, Heriot-Watt
University, Edinburgh EH14 4AS, Scotland, United Kingdom  M. Rambach  SUPA,
Institute of Photonics and Quantum Sciences, Heriot-Watt University, Edinburgh
EH14 4AS, Scotland, United Kingdom  T. S. Santana  Departamento de Física,
Universidade Federal de Sergipe, Sergipe, 49100-000, Brazil  S. I. Park
Center for Opto-Electronic Materials and Devices Research, Korea Institute of
Science and Technology, Seoul 02792, Republic of Korea  J. D. Song  Center for
Opto-Electronic Materials and Devices Research, Korea Institute of Science and
Technology, Seoul 02792, Republic of Korea  E. M. Gauger  SUPA, Institute of
Photonics and Quantum Sciences, Heriot-Watt University, Edinburgh EH14 4AS,
Scotland, United Kingdom  B. D. Gerardot<EMAIL_ADDRESS>SUPA,
Institute of Photonics and Quantum Sciences, Heriot-Watt University, Edinburgh
EH14 4AS, Scotland, United Kingdom

(October 16, 2019)

##  I Experimental Setup and Lifetime Measurements

Fig. S1A shows the schematic of the experimental setup used to perform
spectroscopy measurements on the emission from a semiconductor InGaAs Quantum
dot (QD) mounted in the Voigt configuration. A tunable continuous wave (CW)
laser is used to excite the QD which is kept cold at $\mathrm{T}=4\,K$ in a
closed-cycle helium flow cryostat. In resonance fluorescence (RF), a cross-
polarization scheme using a pair of orthogonally oriented linear polarizers
(LP) on the excitation and the collection arms of the confocal microscope is
used to suppress the background scattering laser up to $10^{7}$. A quarter-
wave plate ($\lambda/4$) is used to correct any birefringence due to optics in
the propagation path. The excitation power is measured by the photodiode (PD).
The photons scattered from the QD are coupled into a single mode (SM) fibre
which is directed to (i) a spectrometer, (ii) a $27.5\,\mathrm{MHz}$
resolution, $5.5\,\mathrm{GHz}$ free spectral range fiber Fabry–Pérot
interferometer (FPI), (iii) a Hanbury-Brown and Twiss interferometer (HBT) or
(iv) an unbalanced Mach-Zehnder interferometer (MZI) setup (with
interferometric delay of $\Delta\mathrm{T}=49.7\,\mathrm{ns}$) to measure
intensity correlations and two-photon interference. An example of the RF
spectrum as the QD is detuned across the laser resonance for the neutral
exciton $\mathrm{X^{0}}$ is shown in Fig. S1B. To perform a lifetime
measurement on the emission from a neutral exciton, instead of a CW laser, the
QD is excited with a mode-locked resonant pulse laser (with ps pulse width and
$80.3\,\mathrm{MHz}$ repetition rate). The scattered photons are detected
using superconducting nanowire single-photon detectors (SNSPD) (timing jitter
$\sim 100\,\mathrm{ps}$). The signal from the SNSPD is then sent to a time-
correlated single photon counting module (PicoHarp 300) with the start channel
synced to the pulsed laser. The results from the lifetime measurement is shown
in Fig. S1C for both QDs studied in the paper. The coincidence histogram shows
a single exponential decay term (corresponding to the lifetime of the emitter,
$T_{1}$) and a fast oscillation term (corresponding to the fine-structure
splitting, $\Delta$ of the neutral exciton). A simple exponential decay with a
sinusoidal function gives $T_{1}=625\,(2)\,\mathrm{ps}$,
$\Delta=18.3\,(1)\,\mu\mathrm{eV}$ and $T_{1}=679\,(1)\,\mathrm{ps}$,
$\Delta=21.3\,(2)\,\mu\mathrm{eV}$ for QD1 and QD2, respectively.

Figure S1: Measurement setup and lifetime of emitters. (A) Schematic of the
experimental setup for the spectroscopy measurement on an InGaAs QD device.
(B) Example of the RF voltage sweeps taken from QD1 and QD2 showing the
neutral exciton $\mathrm{X^{0}}$ lines used in the main text. The excitation
wavelength for QD1 and QD2 are $\lambda=963.45\,\mathrm{nm}$ and
$\lambda=965.9\,\mathrm{nm}$ respectively. (C) Lifetime measurement of the
scattered photons from $\mathrm{X^{0}}$ exhibit a single exponential decay and
a quantum beating. For QD1 (QD2), the lifetime is
$T_{1}=625\,(2)\,\mathrm{ps}$ ($T_{1}=679\,(1)\,\mathrm{ps}$) and the fine
structure splitting (given by the frequency of the beating) is
$\Delta=18.3\,(1)\,\mu\mathrm{eV}$ ($\Delta=21.3\,(2)\,\mu\mathrm{eV}$).
Figure S2: Input ($E_{1,2}(t)$) and output ($E_{3,4}(t)$) field modes at the
second beam splitter of the unbalanced Mach-Zehnder interferometer, with a
temporal delay $\Delta T$ in the upper arm.

##  II Post-Select Two-Photon Interference

Figure S3: Autocorrelation measurement and Two-photon Interference. (A, B)
Correlation measurement for the two-photon interference of indistinguishable
$g^{(2)}_{\parallel}(\tau)$ and distinguishable photons
$g^{(2)}_{\perp}(\tau)$ at saturation, $\Omega\approx\,\Omega_{\mathrm{sat}}$
with their respective fit functions (orange) within a correlation window of
$\pm 100\,\mathrm{ns}$ (A) and $\pm 10\,\mathrm{ns}$ (B). (C) Visibility of
the resonantly scattered photons from QD1 extracted using Eq. S1 at
$\Omega\approx\,\Omega_{\mathrm{sat}}$. The shaded region in
$V_{\mathrm{HOM}}(\tau)$ gives the coalescence time window (CTW) of
$1.24\,(10)\,\mathrm{ns}$. (D) Autocorrelation measurement, $g^{(2)}(\tau)$ at
$\Omega\approx 0.1\,\Omega_{\mathrm{sat}}$ on the resonantly scattered photons
exhibits a suppressed multi-photon emission probability of
$g^{(2)}(0)=0.046\,(13)$.

We perform Hong-Ou-Mandel type two-photon interference measurements Hong  _et
al._ (1987) by sending the resonantly scattered photons to an unbalanced Mach-
Zehnder interferometer (as shown in Fig. S1A(iv)). Here, the photons are split
equally (transmission, $T\sim 53\,\%$) along two paths with one of the paths
being delayed by $\Delta T=49.7\,\mathrm{ns}$ before they interfere on a fibre
beam splitter ($T\sim 51\,\%$), shown in Fig. S2. Here, we perform two sets of
measurements: $g^{(2)}_{\perp}(\tau)$ and $g^{(2)}_{\parallel}(\tau)$. The
visibility is obtained as follows:

$V_{\mathrm{HOM}}(\tau)=(g^{(2)}_{\perp}(\tau)-g^{(2)}_{\parallel}(\tau))/g^{(2)}_{\perp}(\tau).$
(S1)

Fig. S3(A-C) show the results of the two-photon interference measurement from
a neutral exciton, $\mathrm{X^{0}}$ from QD1. The bias applied to the QD
device is selected such that only one of the fine structure of the neutral
exciton is addressed. The solid lines on $g^{(2)}_{\parallel}(\tau)$ and
$g^{(2)}_{\perp}(\tau)$ are produced using the equations from Ref. Proux  _et
al._ (2015), convolved with the setup instrument response function (IRF),
which is determined by a Gaussian peak with a full-width at half maximum
(FWHM) of $160\,\mathrm{ps}$. We take into account the timescale for the
additional bunching (caused by fluctuations in the electronic environment) in
$g^{(2)}_{\perp,\parallel}(\tau)$ by modifying the autocorrelation function
$g^{(2)}(\tau)$ shown in Fig. S3D. By introducing this additional bunching
timescale $\tau_{R}$, the modified autocorrelation function is then given by

$g^{(2)}(\tau)=g^{(2)}_{RF}(\tau)\times(1+b\,e^{-|\tau|/\tau_{R}}),$ (S2)

where $b$ is the amplitude of the bunching and $g^{(2)}_{RF}(\tau)$ is the
autocorrelation function for resonant excitation of a two-level system Scully
and Zubairy (1997). Fitting Eq. (S2) to the experimental data, we obtain
$b=0.245\,(1)$ and $\tau_{R}=137\,(1)\,\mathrm{ns}$.

##  III Coalescence Time Window

The coalescence time window (CTW) is defined as

$\mathrm{CTW}=\int\mathrm{d}\tau\,V_{\mathrm{HOM}}(\tau),$ (S3)

which gives the integrated area of the HOM visibility
$V_{\mathrm{HOM}}(\tau)$. The CTW was first introduced in Ref. (2) as a figure
of merit to quantify photon indistinguishability. The motivation to use CTW
rather than the two photon visibility at zero time delay is that, for
continuous wave excitation, the latter is solely determined by the detector
response time and is not indicative of indistinguishability of the photon
wavepackets Legero  _et al._ (2004). The advantage of using the CTW to
characterise the indistinguishability is that it is not affected by the timing
jitter of the detectors: averaging over temporal delay, it gives the
characteristic time beyond which no two photon interference can be observed.
Indeed, one can see how Eq. (S3) does not depend on detector jitter by noting
that the latter can be accurately described by a convolution of the signal
$V_{\mathrm{HOM}}(\tau)$ with a normal distribution having width proportional
to the detector resolution:

$\displaystyle\begin{split}CTW&=\int^{\infty}_{-\infty}\mathrm{d}\tau\leavevmode\nobreak\
V_{\mathrm{HOM}}(\tau)\star\phi_{\tau_{d}}(\tau)\\\
&=\int^{\infty}_{-\infty}\mathrm{d}\tau\leavevmode\nobreak\
V_{\mathrm{HOM}}(\tau)\int^{\infty}_{-\infty}\mathrm{d}\tau\leavevmode\nobreak\
\phi_{\tau_{d}}(\tau)\\\
&=\int^{\infty}_{-\infty}\mathrm{d}\tau\leavevmode\nobreak\
V_{\mathrm{HOM}}(\tau)\leavevmode\nobreak\ ,\end{split}$ (S4)

where $\tau_{d}$ is the detector resolution, and $\phi_{\tau_{d}}(\tau)$ is
the normal distribution with FWHM $\tau_{d}$ centred at $\tau=0$. The CTW
remains sensitive to experimental imperfections other than the detector
jitter: for example, it will still be affected by the presence of background
scattered light, spatial overlap on the beam splitter and imperfection in
polarization Baudin  _et al._ (2019); Giesz  _et al._ (2015).

Having motivated the use of the CTW over the visibility at zero delay, we now
briefly outline the derivation of our theoretical model, as well as an
approximation leading to the simplified visibility model given in Ref. Iles-
Smith  _et al._ (2017). We denote the operators of the field inputs prior to
the second beam splitter by $E_{1}(t)$ and $E_{2}(t)$, respectively (c.f. Fig.
S2). Since for an unbalanced Mach-Zehnder setup the two input operators differ
only by the temporal delay in one of the arms, we may write these as
$E_{1}(t)=E(t+\Delta T)$ and $E_{2}(t)=E(t)$. Furthermore, we can relate the
input field operators to the QD emission by using $E(t)=E_{0}\sigma_{-}(t)$,
where $E_{0}$ is the vacuum field, and $\sigma_{-}(t)$ is the Heisenberg
picture dipole lowering operator Ficek  _et al._ (2005a). Thus, we may write
the output field modes after the secondary beam splitter as

$\displaystyle\begin{split}E_{3}(t)&=\frac{1}{\sqrt{2}}(E_{1}(t)+E_{2}(t))=\frac{1}{\sqrt{2}}(\sigma_{-}(t+\Delta
T)+\sigma_{-}(t))\leavevmode\nobreak\ ,\\\
E_{4}(t)&=\frac{1}{\sqrt{2}}(E_{1}(t)-E_{2}(t))=\frac{1}{\sqrt{2}}(\sigma_{-}(t+\Delta
T)-\sigma_{-}(t))\leavevmode\nobreak\ .\end{split}$ (S5)

We are now in a position to fully expand the second order correlation function
for the output field modes $E_{3}(\tau)$ and $E_{4}(\tau)$, given by

$g^{(2)}(\tau)=\frac{\Braket{E^{*}_{3}(t)E^{*}_{4}(t+\tau)E_{4}(t+\tau)E_{3}(t)}}{\Braket{E^{*}_{3}(t+\tau)E_{3}(t+\tau)}\Braket{E^{*}_{4}(t+\tau)E_{4}(t+\tau)}}\leavevmode\nobreak\
.$ (S6)

Expanding $g^{(2)}(\tau)$, we obtain 16 terms involving the dipole operators.
Whilst we refrain from writing the full expression here, we note that each
term has one of the three forms:

(1) $\Braket{\sigma_{+}(t)\sigma_{+}(t+\tau+k\Delta
T)\sigma_{-}(t+\tau+k\Delta T)\sigma_{-}(t)},\leavevmode\nobreak\ k\in\\{0,\pm
1\\}$: which is essentially the standard second order correlation function
evaluated at different delays.

(2) Correlation operators involving an odd number of dipole operator
contributions from the two interferometer arms (such as
$\Braket{\sigma_{+}(t+\Delta
T)\sigma_{+}(t+\tau)\sigma_{-}(t+\tau)\sigma_{-}(t)}$). Whilst these terms
vanish for low laser coherence $T_{L}$ compared to $\Delta T$ Lebreton  _et
al._ (2013); Baudin  _et al._ (2019), we are in the opposite regime, where
$T_{L}\gg\Delta T$, and thus we keep these terms in our calculations.

(3) $\Braket{\sigma_{+}(t+\Delta T)\sigma_{+}(t+\tau)\sigma_{-}(t+\tau+\Delta
T)\sigma_{-}(t)}$ and its complex conjugate. If we treat the output from both
arms as independent, we can write this term as

$\displaystyle\begin{split}&\Braket{\sigma_{+}(t+\Delta
T)\sigma_{+}(t+\tau)\sigma_{-}(t+\tau+\Delta T)\sigma_{-}(t)}\\\
&\approx\Braket{\sigma_{+}(t+\Delta T)\sigma_{-}(t+\tau+\Delta
T)}\Braket{\sigma_{+}(t+\tau)\sigma_{-}(t)}\\\
&=|g^{(1)}(\tau)|^{2}\leavevmode\nobreak\ .\end{split}$ (S7)

We find that this approximation is valid for a delay $\Delta T=49.7$ ns. In
the absence of the vibrational environment, it has been shown that the CTW
increases beyond $2\,T_{1}$ as $\Omega\rightarrow 0$ due to the increasing
elastic fraction of scattered photons Proux  _et al._ (2015); Baudin  _et al._
(2019).

##  IV Spin-$\Lambda$ System Coupled to Nuclear Spin Bath

Figure S4: Resonance fluorescence spectra from the same negatively-charged
exciton, $\mathrm{X^{1-}}$ transition at below saturation, $\Omega\approx
0.1\,\Omega_{sat}$. (A) At $B_{\perp}=0\,\mathrm{T}$, the fluorescence spectra
of the photons scattered from the negatively-charged exciton
($\mathrm{X^{1-}}$) shows an elastic peak (red dashed line) and two displaced
inelastic Raman peaks (green solid line). The origin of the splitting is due
to the coupling of the electron spin to the nuclear spin bath, whereby it
induces a Zeeman splitting proportional to the mean nuclear spin field,
$\left<B_{N}\right>$. The linewidth of the Raman peaks gives
$\Gamma=421\,(11)\,\mathrm{MHz}$ and is proportional to the nuclear spin
fluctuation, $\delta B$. (B) At $B_{\perp}=2\,\mathrm{T}$, the emission
spectra from the Raman transition shows a single Gaussian line-shape. The
linewidth $\Gamma$, extracted from the fit (orange) is
${\Gamma}=439\,(20)\,\mathrm{MHz}$, which shows agreement with the linewidth
obtained in the zero field.

Fig. S4 shows the emission spectra from a resonantly driven negatively charged
exciton, $\mathrm{X^{1-}}$ at zero field (A) and Voigt field (4 T) (B) in the
weak excitation regime, $\Omega\approx 0.1\,\Omega_{\mathrm{sat}}$. At zero
field, the emission spectra of $\mathrm{X^{1-}}$ consists of an elastic peak
(given by the red dashed line) and two Raman peaks (given by the green dashed-
dotted line). The origin of the presence of the Raman peaks is due to the non-
zero coupling of the ground state electron spins to the nuclear spin bath
Urbaszek  _et al._ (2013); Malein  _et al._ (2016). This gives rise to non-
zero nuclear-spin induced magnetic field (Overhauser field),
$\left<B_{N}\right>\approx 346\,(11)\,\mathrm{MHz}\approx 25\,\mathrm{mT}$,
which splits the ground state of $\mathrm{X^{1-}}$. The broken symmetry forms
a four-level double-$\Lambda$ system (see Fig. 1(D) in main text). This allows
the previously forbidden diagonal transition (in the spin-z
($\sigma_{\mathrm{z}}$) basis) which gives the Raman side peaks Malein  _et
al._ (2016). The ground state spin flip rate is proportional to the Overhauser
field fluctuations, $\delta B_{N}=188\,(5)\,\mathrm{MHz}$. In the absence of
the finite exciton-phonon coupling, since the coherence of the Raman photons
is determined by the ground state dephasing rate Fernandez  _et al._ (2009);
He  _et al._ (2013); Sun  _et al._ (2016), the linewidth of the Raman side
peaks is given by the spin flip rate, which is $2\sqrt{2\log 2}\,\delta
B_{N}=421\,(11)\,\mathrm{MHz}$. This result matches the linewidth of the
scattered Raman photons at $B_{\mathrm{ext}}=B_{\perp}=2\,\mathrm{T}$,
$\Gamma=439\,(20)\,\mathrm{MHz}$, within one standard deviation.

Alternatively, one could solve the rate equations of a three level
spin-$\Lambda$ system analytically for a strong external magnetic field,
$B_{\mathrm{ext}}\gg\delta B_{N}$ such that the Overhauser field fluctuations
are partially suppressed Malein  _et al._ (2016). In this regime, the
component of $\delta B_{N}$ perpendicular to the applied field,
$B_{\mathrm{ext}}$ can be ignored to a good approximation. Assuming both
ground states are equally populated initially, the resulting emission
spectrum, $S(\omega)$, contains two delta functions, one being the elastic
peak with the same frequency as the laser, $\omega$, and the other being a
shifted inelastic peak with frequency given by $\omega-B$,

$S(\omega)=\frac{\Omega^{2}}{2\gamma}[\delta(\omega)+\delta(\omega-B)],$ (S8)

where $B=B_{\mathrm{ext}}+B_{N}$, $\Omega$ is the driving Rabi frequency and
$\gamma$ is the spontaneous emission rate. Averaging the emission spectrum
$S(\omega)$ over the Overhauser field fluctuations, the resulting spectrum,
$\langle S(\omega)\rangle$ is

$\langle
S(\omega)\rangle=\frac{\Omega^{2}}{2\gamma}\left[\delta(\omega)+\frac{e^{\left(\omega-
B_{\mathrm{ext}}\right)^{2}/2\delta B_{N}^{2}}}{\delta
B_{N}\sqrt{2\pi}}\right],$ (S9)

where the first term represents the elastically scattered Rayleigh photons,
given by the delta function and the second term describes the inelastically
scattered Raman photons which is given by a Gaussian with a full-width at half
maximum (FWHM) of $2\sqrt{2\log 2}\,\delta B_{N}=439\,(20)\,\mathrm{MHz}$. In
the experiment, since a cross polarization technique is used to collect only
the orthogonally polarized Raman photons, we could neglect the first term in
Eq. S9. Hence we are left with the Gaussian term which describes the spectrum
of the scattered Raman photons. Using the linewidth of the Gaussian fit on the
spectrum in Fig. S4B, we get $\delta B_{N}=186\,(8)\,\mathrm{MHz}$ which shows
agreement with the results obtained from fitting the zero field spectra with
the model described in Ref. Malein  _et al._ (2016).

Disregarding phonon interactions, this agrees with the previous predictions
that the coherence of the Raman zero phonon line is given by the ground (spin)
state dephasing rate, proportional to the Overhauser field fluctuations (Refs.
Fernandez  _et al._ (2009); Malein  _et al._ (2016); Sun  _et al._ (2016)).
However, the presence of the phonon sideband that accompanied by the Raman
zero phonon line (as shown in Fig. 4A) modifies this picture and renders the
scattered Raman photons partially incoherent.

##  V X1- spin-$\Lambda$ system dynamics

In this section we present the model used for the theoretical spectra using
the polaron frame master equation formalism. We focus on the $\mathrm{X^{1-}}$
$\Lambda$ system formed by driving the spin-conserving
$\Ket{\uparrow}\rightarrow\Ket{T_{\Uparrow}}$ transition and collect emission
from the spin-flipping $\Ket{T_{\Uparrow}}\rightarrow\Ket{\downarrow}$
transition. The full Hamiltonian for the four level system discussed in the
main text can be written as:

$\displaystyle\begin{split}H_{4LS}=&H_{0}+H_{\Omega}\\\
H_{0}=&\frac{\delta_{e}}{2}\left(\Ket{\downarrow}\Bra{\downarrow}-\Ket{\uparrow}\Bra{\uparrow}\right)+\frac{\delta_{h}}{2}\left(\Ket{T_{\Uparrow}}\Bra{T_{\Uparrow}}-\Ket{T_{\Downarrow}}\Bra{T_{\Downarrow}}\right)\\\
&-\frac{\omega_{0}}{2}\left(\Ket{\uparrow}\Bra{\uparrow}+\Ket{\uparrow}\Bra{\uparrow}\right)\\\
&+\frac{\omega_{0}}{2}\left(\Ket{T_{\Uparrow}}\Bra{T_{\Uparrow}}+\Ket{T_{\Downarrow}}\Bra{T_{\Downarrow}}\right)\\\
H_{\Omega}=&\leavevmode\nobreak\
\Omega\cos(\omega_{L}t)\left(\Ket{T_{\Uparrow}}\Bra{\uparrow}+\Ket{\uparrow}\Bra{T_{\Uparrow}}\right)\leavevmode\nobreak\
,\end{split}$ (S10)

where $\Ket{\uparrow,\downarrow}$ and
$\Ket{T_{\Uparrow,\Downarrow}}=\Ket{\uparrow\downarrow\Uparrow,\uparrow\downarrow\Downarrow}$
are the negatively charged ground and trion states, respectively. $\omega_{0}$
and $\omega_{L}$ are the exciton transition and laser frequencies,
respectively.

Starting with the system initialized in the $\Ket{\uparrow}$ ground state, we
drive using vertically polarized light (driving the
$\Ket{\uparrow}\rightarrow\Ket{T_{\Uparrow}}$ transition), and collect cross
polarized light ($\Ket{T_{\Uparrow}}\rightarrow\Ket{\downarrow}$). After
moving to the rotating frame with respect to the driven transition, and
subtracting a term proportional the identity, we focus on the spin-$\Lambda$
system formed by the states
$\\{\Ket{\uparrow},\Ket{\downarrow},\Ket{T_{\Uparrow}}\\}$ in order to obtain
the Hamiltonian

$\displaystyle\begin{split}H=&\frac{\delta}{2}\Ket{\downarrow}\Bra{\downarrow}+\frac{\Omega}{2}\left(\Ket{T_{\Uparrow}}\Bra{\uparrow}+\Ket{\uparrow}\Bra{T_{\Uparrow}}\right)\leavevmode\nobreak\
,\\\ \end{split}$ (S11)

where $\delta=2\delta_{e}$.

###  V.1 Exciton–photon interaction

The photonic environment can be modelled by the Hamiltonian

$H^{pt}_{E}=\sum_{\mathbf{q},\,\lambda}\nu_{\mathbf{q}}a^{\dagger}_{\mathbf{q}\lambda}a_{\mathbf{q}\lambda}\leavevmode\nobreak\
,$ (S12)

where $a^{\dagger}_{\mathbf{q}\lambda}$ ($a_{\mathbf{q}\lambda}$) is the
creation (annihilation) operator for a photon of momentum $\mathbf{q}$ and
polarization $\lambda$. In the dipole approximation, the photon interaction
Hamiltonian is of the form

$\displaystyle\begin{split}H^{pt}_{I}=&-\mathbf{d}_{\Uparrow\downarrow}\cdot\mathbf{E}(\mathbf{r}_{d})(\Ket{\downarrow}\Bra{T_{\Uparrow}}+\Ket{T_{\Uparrow}}\Bra{\downarrow})\\\
&-\mathbf{d}_{\Uparrow\uparrow}\cdot\mathbf{E}(\mathbf{r}_{d})(\Ket{\uparrow}\Bra{T_{\Uparrow}}+\Ket{T_{\Uparrow}}\Bra{\uparrow})\leavevmode\nobreak\
,\end{split}$ (S13)

with $\mathbf{E}(\mathbf{r}_{d})$ being the Schrödinger picture electric field
in free space at the location $\mathbf{r}_{d}$ of the dipole, and
$\mathbf{d}_{\Uparrow\downarrow}=\Braket{\downarrow}{\mathbf{D}}{T_{\Uparrow}}$
and
$\mathbf{d}_{\Uparrow\uparrow}=\Braket{\uparrow}{\mathbf{D}}{T_{\Uparrow}}$
are the transition matrix element of the dipole operator for the
$\Ket{T_{\Uparrow}}\rightarrow\Ket{\downarrow}$ and
$\Ket{T_{\Uparrow}}\rightarrow\Ket{\uparrow}$ transitions respectively.

###  V.2 Exciton–phonon interaction

Unlike for atomic systems, the vibrational environment plays a key role in the
dynamics of a confined electron in a QD. As an electron is excited from the
valence to the conduction band, the charge configuration of the semiconductor
is modified accordingly. This results in a shift of the lattice ions’
equilibrium positions, giving rise to an exciton–phonon coupling depending on
the excitons’s state. The phonon bath can be described via the Hamiltonian

$H^{pn}_{E}=\sum_{\mathbf{k}}\omega_{\mathbf{k}}b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}\leavevmode\nobreak\
,$ (S14)

where $b^{\dagger}_{\mathbf{k}}$ ($b_{\mathbf{k}}$) is the creation
(annihilation) operator for a phonon of momentum $\mathbf{k}$. The
corresponding Hamiltonian governing the exciton–phonon interaction dynamics is
then given by

$H^{pn}_{I}=\Ket{T_{\Uparrow}}\Bra{T_{\Uparrow}}\sum_{\mathbf{k}}g_{\mathbf{k}}(b^{\dagger}_{\mathbf{k}}+b_{\mathbf{k}})\leavevmode\nobreak\
,$ (S15)

where $g_{\mathbf{k}}$ is the coupling strength of the excited electronic
configuration with phonon mode $\mathbf{k}$, given by Mahan (2000)

$g_{\mathbf{k}}=\left(\frac{\hbar}{2\rho\nu\omega_{\mathbf{k}}}\right)^{1/2}\left[\tilde{M}_{X}(\mathbf{k})\mathrm{e}^{-d^{2}_{X}|\mathbf{k}|^{2}/4}-\tilde{M}_{0}(\mathbf{k})\mathrm{e}^{-d^{2}_{0}|\mathbf{k}|^{2}/4}\right]\leavevmode\nobreak\
,$ (S16)

where $\rho$ is the mass density of the solid, $\nu$ is the lattice volume,
$\tilde{M}_{X}(\mathbf{k})$ and $\tilde{M}_{0}(\mathbf{k})$ are the long-
wavelength phonon coupling matrix elements for the excited and ground state,
respectively, and $d_{X}$ and $d_{0}$ characterize the size of the electron
and hole wavefunctions, respectively.

###  V.3 Master equation dynamics: strong vibrational coupling

To summarize, our system consists of the following Hamiltonian:

$\displaystyle\begin{split}H_{tot}=&H+H^{pt}_{E}+H^{pn}_{E}+H^{pt}_{I}+H^{pn}_{I}\leavevmode\nobreak\
;\\\
\vspace{5mm}H=&\delta\Ket{\downarrow}\Bra{\downarrow}+\frac{\Omega}{2}\left(\Ket{T_{\Uparrow}}\Bra{\uparrow}+\Ket{\uparrow}\Bra{T_{\Uparrow}}\right)\leavevmode\nobreak\
,\\\
H^{pt}_{E}=&\sum_{\mathbf{q},\,\lambda}\nu_{\mathbf{q}}a^{\dagger}_{\mathbf{q}\lambda}a_{\mathbf{q}\lambda}\leavevmode\nobreak\
,\\\
H^{pn}_{E}=&\sum_{\mathbf{k}}\omega_{\mathbf{k}}b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}\leavevmode\nobreak\
,\\\
H^{pt}_{I}=&-\mathbf{d}_{\Uparrow\downarrow}\cdot\mathbf{E}(\mathbf{r}_{d})(\Ket{\downarrow}\Bra{T_{\Uparrow}}+\Ket{T_{\Uparrow}}\Bra{\downarrow})\\\
&-\mathbf{d}_{\Uparrow\uparrow}\cdot\mathbf{E}(\mathbf{r}_{d})(\Ket{\uparrow}\Bra{T_{\Uparrow}}+\Ket{T_{\Uparrow}}\Bra{\uparrow})\leavevmode\nobreak\
,\\\
H^{pn}_{I}=&\Ket{T_{\Uparrow}}\Bra{T_{\Uparrow}}\sum_{\mathbf{k}}g_{\mathbf{k}}(b^{\dagger}_{\mathbf{k}}+b_{\mathbf{k}})\leavevmode\nobreak\
.\end{split}$ (S17)

Having described full exciton, vibrational and photonic environment systems,
we are now in a position to briefly go through the details of the master
equation describing the excitonic dynamics. For a solid state system, we
expect a significantly strong coupling to the vibrational environment. Hence,
the polaron, or Lang-Firsov Lang and Firsov (1962), frame adequately captures
the dynamics of our system Nazir and McCutcheon (2016); Scerri  _et al._
(2017). The two transitions $\Ket{T_{\Uparrow}}\rightarrow\Ket{\downarrow}$
and $\Ket{T_{\Uparrow}}\rightarrow\Ket{\uparrow}$ couple to the same
vibrational environment, and hence the polaron transformation can be written
as $U_{p}=\mathrm{e}^{S_{\Uparrow}}$, where

$S_{\Uparrow}=\Ket{T_{\Uparrow}}\Bra{T_{\Uparrow}}\sum_{\mathbf{k}}\frac{g_{\mathbf{k}}}{\omega_{\mathbf{k}}}(b^{\dagger}_{\mathbf{k}}-b_{\mathbf{k}})\leavevmode\nobreak\
.$ (S18)

It can be easily shown that $U_{p}$ can be simplified to

$\displaystyle\begin{split}U_{P}&=\Ket{\uparrow}\Bra{\uparrow}+\Ket{\downarrow}\Bra{\downarrow}+\Ket{T_{\Uparrow}}\Bra{T_{\Uparrow}}B_{+}\leavevmode\nobreak\
;\\\
B_{+}&=\prod_{k}\mathrm{e}^{\frac{g_{\mathbf{k}}}{\omega_{\mathbf{k}}}(b^{\dagger}_{\mathbf{k}}-b_{\mathbf{k}})}\leavevmode\nobreak\
,\end{split}$ (S19)

and that, in the absence of driving, diagonalizes the phonon interaction
$H^{pn}_{I}$. Using this transformation, we obtain the Hamiltonian in the
polaron frame, indexed by the subscript $P$:

$\displaystyle\begin{split}H_{tot\hskip
2.84526ptP}=&H_{P}+H^{pt}_{E}+H^{pn}_{EP}+H^{pt}_{IP}+H^{pn}_{IP}\leavevmode\nobreak\
;\\\
\vspace{5mm}H_{P}=&\delta_{P}\Ket{\downarrow}\Bra{\downarrow}+\frac{\Omega_{P}}{2}\left(\Ket{T_{\Uparrow}}\Bra{\uparrow}+\Ket{\uparrow}\Bra{T_{\Uparrow}}\right)\leavevmode\nobreak\
,\\\
H^{pt}_{EP}=&\sum_{\mathbf{q},\,\lambda}\nu_{\mathbf{q}}a^{\dagger}_{\mathbf{q}\lambda}a_{\mathbf{q}\lambda}\leavevmode\nobreak\
,\\\
H^{pn}_{EP}=&\sum_{\mathbf{k}}\omega_{\mathbf{k}}b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}\leavevmode\nobreak\
,\\\
H^{pt}_{IP}=&-\mathbf{d}_{\Uparrow\downarrow}\cdot\mathbf{E}(\mathbf{r}_{d})(B_{-}\Ket{\downarrow}\Bra{T_{\Uparrow}}+B_{+}\Ket{T_{\Uparrow}}\Bra{\downarrow})\\\
&-\mathbf{d}_{\Uparrow\uparrow}\cdot\mathbf{E}(\mathbf{r}_{d})(B_{-}\Ket{\uparrow}\Bra{T_{\Uparrow}}+B_{+}\Ket{T_{\Uparrow}}\Bra{\uparrow})\leavevmode\nobreak\
,\\\
H^{pn}_{IP}=&\frac{\Omega}{2}\Big{[}\Ket{\uparrow}\Bra{T_{\Uparrow}}(B_{-}-\Braket{B})+\Ket{T_{\Uparrow}}\Bra{\uparrow}(B_{+}-\Braket{B})\Big{]}\leavevmode\nobreak\
.\end{split}$ (S20)

where $B_{-}=B^{\dagger}_{+}$, $\Braket{B}=\Braket{B_{+}}=\Braket{B_{-}}$,
$\Omega_{P}=\Braket{B}\Omega$, and
$\delta_{P}=\delta-\sum_{\mathbf{k}}g^{2}_{\mathbf{k}}/\omega_{\mathbf{k}}\rightarrow\delta-\int^{\infty}_{0}\mathrm{d}\omega
J_{pn}(\omega)/\omega$, with $J_{pn}(\omega)$ being the phonon spectral
density, i.e.
$J_{pn}(\omega)=\alpha\omega^{3}\mathrm{e}^{-\frac{\omega^{2}}{\omega^{2}_{c}}}$,
where $\alpha$ is a measure of the coupling strength and $\omega_{c}$ is the
phonon frequency cut-off, which depends on the size and confinement of the
quantum dot. The parameters were chosen to be $\alpha=0.03\,\mathrm{ps}^{2}$
and $\omega_{c}=2.2\,\mathrm{ps}^{-1}$, which agree with the standard values
for self-assembled GaAs quantum dots Ramsay  _et al._ (2010). The new
vibrational interaction Hamiltonian $H^{pn}_{IP}$ appears in the total
Hamiltonian due to the driving term. The details of the derivation can be
found in Refs. (17; 18; 20; 21), the main result being that the new
interaction can now be treated perturbatively in the master equation (ME)
derivation.

In order to derive the corresponding emission rates, it would be beneficial to
write the two interaction Hamiltonians in a more compact form. Thus, with the
definitions $A_{1,\uparrow}^{pt}=\Ket{\uparrow}\Bra{T_{\Uparrow}}$
$\left(A_{1,\downarrow}^{pt}=\Ket{\downarrow}\Bra{T_{\Uparrow}}\right)$,
$A_{2,\uparrow}^{pt}=A_{1,\uparrow}^{pt\dagger}$
$\left(A_{2,\downarrow}^{pt}=A_{1,\downarrow}^{pt\dagger}\right)$,
$B^{pt}_{1/2}\equiv B_{\mp}$,
$C_{1}=i\sum_{\mathbf{q},\lambda}\mathbf{d}\cdot\mathbf{u}^{*}_{\mathbf{q}\lambda}(\mathbf{r}_{d})a^{\dagger}_{\mathbf{q}\lambda}$,
and $C_{2}=C^{\dagger}_{1}$, and expanding $\mathbf{E}(\mathbf{r}_{d})$ into
the corresponding photonic creation and annihilation operators Breuer and
Petruccione (2007), we can write the Hamiltonian $H^{pt}_{IP}$ Scerri  _et
al._ (2017) as

$H^{pt}_{IP}=\sum_{j\in\\{\uparrow,\downarrow\\}}\sum^{2}_{i=1}A^{pt}_{i,j}\otimes
B^{pt}_{i}\otimes C_{i}\leavevmode\nobreak\ .$ (S21)

Similarly we can write a more compact phonon interaction Hamiltonian

$H^{pn}_{IP}=\sum_{i=1}^{2}A^{pn}_{i,\uparrow}\otimes
B^{pn}_{i}\leavevmode\nobreak\ ,$ (S22)

where $B^{pn}_{1/2}=\mathcal{B}_{\mp}=B_{\mp}-\Braket{B}$,
$A^{pn}_{1,\uparrow}=\Omega/2\,\Ket{\uparrow}\Bra{T_{\Uparrow}}$ and
$A^{pn}_{2,\uparrow}=A_{1,\uparrow}^{pn\dagger}$. Following the steps of Refs.
(17; 18), we can now obtain the following ME in the polaron frame:

$\displaystyle\begin{split}\diff{}{t}&\rho_{SP}(t)=\\\
&-\int_{0}^{\infty}\mathrm{d}\tau\mathrm{Tr}^{pn}_{E}[H^{pn}_{IP}(t),[H^{pn}_{IP}(t-\tau),\rho_{SP}(t)\otimes\rho^{pn}_{E}(0)]]\leavevmode\nobreak\
\\\
&-\int_{0}^{\infty}\mathrm{d}\tau\mathrm{Tr}_{E}[H^{pt}_{IP}(t),[H^{pt}_{IP}(t-\tau),\rho_{SP}(t)\otimes\rho_{E}(0)]]\leavevmode\nobreak\
,\end{split}$ (S23)

where $\rho_{SP}(t)$ is the density matrix of the $\Lambda$-system in the
polaron frame, and $\rho_{E}$ is the joint optical-vibrational density matrix.
The above was derived under the assumption that the (initial) environmental
state is thermal, hence $\rho_{E}(0)$ factorizes as
$\rho_{E}(0)=\rho^{pn}_{E}(0)\otimes\rho^{pt}_{E}(0)$.

In the ME formalism, the rate $\gamma(\omega)$ of a dissipative process is
given by
$\gamma(\omega)=2\mathrm{Re}\left[\int_{0}^{\infty}\mathrm{d}sK(s)\right]$,
where $K(s)$ is the relevant correlation function [c.f. Eq. (3.137) in Ref.
Breuer and Petruccione (2007)]. For our phonon dissipator (first term in Eq.
S23), these functions are given by

$\displaystyle\begin{split}C^{pn}_{ii}(\tau)&=\mathrm{Tr}^{pn}_{E}\left[\mathcal{B}^{\dagger}_{\pm}(\tau)\mathcal{B}_{\pm}(0)\rho^{pn}_{E}(0)\right]\\\
&=\langle B\rangle^{2}(\mathrm{e}^{\phi(\tau)}-1)\leavevmode\nobreak\
,\end{split}$ (S24)
$\displaystyle\begin{split}C^{pn}_{ij}(\tau)&=\mathrm{Tr}^{pn}_{E}\left[\mathcal{B}^{\dagger}_{\pm}(\tau)\mathcal{B}_{\mp}(0)\rho^{pn}_{E}(0)\right]\\\
&=\langle B\rangle^{2}(\mathrm{e}^{-\phi(\tau)}-1)\leavevmode\nobreak\
,\end{split}$ (S25)

where $i,j\in\\{1,2\\}$, $i\neq j$. The temperature-dependent phonon
propagator $\phi(\tau)$ is given by

$\phi(\tau)=\int^{\infty}_{0}\frac{J_{pn}(\omega)}{\omega^{2}}\left(\cos(\omega\tau)\coth(\beta\omega/2)-i\sin(\omega\tau)\right)\leavevmode\nobreak\
,$ (S26)

where $\beta=1/k_{B}T$, with $k_{B}$ being the Boltzmann constant and $T$ the
temperature. This allows us to write the operator expectation $\langle
B\rangle$ as

$\displaystyle\begin{split}\langle
B\rangle&=\exp\left[-\frac{1}{2}\int^{\infty}_{0}\frac{J_{pn}(\omega)}{\omega^{2}}\coth(\beta\omega/2)\right]\\\
&=\exp\left[-\frac{1}{2}\phi(0)\right]\leavevmode\nobreak\ .\end{split}$ (S27)

After some algebra, we obtain a phonon dissipator of the form

$\displaystyle\begin{split}&\gamma^{pn}(\omega^{\prime})(\mathcal{L}[\sigma^{\uparrow}_{-}]+\mathcal{L}[\sigma^{\downarrow}_{-}])+\gamma^{pn}(-\omega^{\prime})(\mathcal{L}[\sigma^{\uparrow}_{+}]+\mathcal{L}[\sigma^{\downarrow}_{+}])\\\\[10.0pt]
&-\gamma^{pn}_{cd}(\omega^{\prime})(\mathcal{L}_{cd}[\sigma^{\uparrow}_{-}]+\mathcal{L}_{cd}[\sigma^{\downarrow}_{-}])\\\\[10.0pt]
&-\gamma^{pn}_{cd}(-\omega^{\prime})(\mathcal{L}_{cd}[\sigma^{\uparrow}_{+}]+\mathcal{L}_{cd}[\sigma^{\downarrow}_{+}])\leavevmode\nobreak\
,\end{split}$ (S28)

where
$\sigma^{\uparrow,\downarrow}_{-}=\Ket{\uparrow,\downarrow}\Bra{T_{\Uparrow}}$,
$\sigma^{\uparrow,\downarrow}_{+}=\Ket{T_{\Uparrow}}\Bra{\uparrow,\downarrow}$,
$\mathcal{L}[C]=C\rho_{SP}C^{\dagger}-\frac{1}{2}\\{C^{\dagger}C,\rho_{SP}\\}$
and $\mathcal{L}_{cd}[C]=C\rho_{SP}C-\frac{1}{2}\\{C^{2},\rho_{SP}\\}$. The
phonon absorption ($\gamma^{pn}(\omega^{\prime})$), emission
($\gamma^{pn}(-\omega^{\prime})$) and cross-dephasing rates
($\gamma^{pn}_{cd}(-\omega^{\prime})$) Scerri  _et al._ (2017); Roy-Choudhury
and Hughes (2015a, b) are given by

$\displaystyle\gamma^{pn}(\pm\omega^{\prime})$
$\displaystyle=\frac{\left|\Omega^{pn}\right|^{2}}{4}\int_{-\infty}^{\infty}\mathrm{d}\tau\;\mathrm{e}^{\pm
i\omega^{\prime}\tau}\left(\mathrm{e}^{\phi(\tau)}-1\right)\leavevmode\nobreak\
,$ $\displaystyle\gamma^{pn}_{cd}(\omega^{\prime})$
$\displaystyle=\frac{\left(\Omega^{pn}\right)^{2}}{4}\int_{-\infty}^{\infty}\mathrm{d}\tau\;\cos(\omega^{\prime}\tau)\left(1-\mathrm{e}^{-\phi(\tau)}\right)\leavevmode\nobreak\
,$ $\displaystyle\gamma^{pn}_{cd}(-\omega^{\prime})$
$\displaystyle=\frac{\left(\Omega^{pn}\right)^{2}}{4}\int_{-\infty}^{\infty}\mathrm{d}\tau\;\cos(\omega^{\prime}\tau)\left(1-\mathrm{e}^{-\phi(\tau)}\right)\leavevmode\nobreak\
,$

Having derived the corresponding phonon absorption/emission rates, we now
briefly turn our attention to the second term in Eq. (S23). This term governs
the system’s interaction with the electromagnetic environment in the polaron
frame, and thus gives the corresponding correlation functions

$\displaystyle
C^{pt}_{ij,\uparrow}(\tau)=C^{pt}_{ij,\downarrow}(\tau)\coloneqq
C^{pt}_{ij}(\tau)$ (S29)
$\displaystyle=\mathrm{Tr}_{E}\left[\left(B^{pt\dagger}_{i}(\tau)\otimes
C^{\dagger}_{i}(\tau)\right)\left(B^{pt}_{j}(0)\otimes
C_{j}(0)\right)\rho_{E}(0)\right]\leavevmode\nobreak\ $
$\displaystyle=\mathrm{Tr}^{pn}_{E}\left[B^{pt\dagger}_{i}(\tau)B^{pt}_{j}(0)\rho^{pn}_{E}(0)\right]\mathrm{Tr}^{pt}_{E}\left[C^{\dagger}_{i}(\tau)C_{j}(0)\rho^{pt}_{E}(0)\right]\leavevmode\nobreak\
,$

where $i,j\in\\{1,2\\}$. As the two transitions
$\Ket{\uparrow}\leftrightarrow\Ket{T_{\Uparrow}}$ and
$\Ket{\downarrow}\leftrightarrow\Ket{T_{\Uparrow}}$ are coupled to both the
vibrational and electromagnetic baths, the above correlation function is the
same for both transitions of the $\Lambda$ system. We also note that the cross
correlation terms involving both transitions will vanish due to the orthogonal
dipole moments of the two transitions Ficek  _et al._ (2005b). After
substituting for the bath operators, and making use of the
creation/annihilation commutation relations, we get that the only the
$C^{pt}_{11}$ correlation function is non zero, giving the standard photon
emission rates $\gamma_{11,j}=\omega^{3}_{j}d^{2}_{j}/(3\pi\epsilon_{0}\hbar
c^{3})$, $j\in\\{\uparrow,\downarrow\\}$ Ficek  _et al._ (2005b).

##  VI X0 two-level system dynamics

In this section, we briefly sketch the derivation for the X0 exciton master
equation, analogous to the charged exciton master equation from the previous
section. Denoting the ground and excited states for this two-level system as
$\Ket{0}$ and $\Ket{X}$, respectively, the Hamiltonian in the frame rotating
at the exciton transition frequency $\omega_{0}$ is given by

$\displaystyle\begin{split}H_{tot}=&H+H^{pt}_{E}+H^{pn}_{E}+H^{pt}_{I}+H^{pn}_{I}\leavevmode\nobreak\
;\\\
\vspace{5mm}H=&\frac{\Omega}{2}\left(\Ket{X}\Bra{0}+\Ket{0}\Bra{X}\right)\leavevmode\nobreak\
,\\\
H^{pt}_{E}=&\sum_{\mathbf{q},\,\lambda}\nu_{\mathbf{q}}a^{\dagger}_{\mathbf{q}\lambda}a_{\mathbf{q}\lambda}\leavevmode\nobreak\
,\\\
H^{pn}_{E}=&\sum_{\mathbf{k}}\omega_{\mathbf{k}}b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}\leavevmode\nobreak\
,\\\
H^{pt}_{I}=&-\mathbf{d}_{X0}\cdot\mathbf{E}(\mathbf{r}_{d})(\Ket{0}\Bra{X}+\Ket{X}\Bra{0})\\\
&-\mathbf{d}_{X0}\cdot\mathbf{E}(\mathbf{r}_{d})(\Ket{0}\Bra{X}+\Ket{X}\Bra{0})\leavevmode\nobreak\
,\\\
H^{pn}_{I}=&\Ket{X}\Bra{X}\sum_{\mathbf{k}}g_{\mathbf{k}}(b^{\dagger}_{\mathbf{k}}+b_{\mathbf{k}})\leavevmode\nobreak\
,\end{split}$ (S30)

where $\mathbf{d}_{X0}$ is the transition matrix element for the exciton
transition. Thus, using the polaron transformation
$S_{X}=\Ket{X}\Bra{X}\sum_{\mathbf{k}}\frac{g_{\mathbf{k}}}{\omega_{\mathbf{k}}}(b^{\dagger}_{\mathbf{k}}-b_{\mathbf{k}})$
and following the derivation in the previous section and Ref. (Nazir and
McCutcheon, 2016), we arrive at a similar master equation for the neutral
exciton, with the phonon dissipator given by

$\displaystyle\begin{split}&\gamma^{pn}(\omega^{\prime})\mathcal{L}[\sigma_{-}]+\gamma^{pn}(-\omega^{\prime})\mathcal{L}[\sigma_{+}]\\\\[10.0pt]
&-\gamma^{pn}_{cd}(\omega^{\prime})\mathcal{L}_{cd}[\sigma_{-}]-\gamma^{pn}_{cd}(-\omega^{\prime})\mathcal{L}_{cd}[\sigma_{+}]\leavevmode\nobreak\
,\end{split}$ (S31)

where $\sigma_{-}=\Ket{0}\Bra{X}$, whilst the electromagnetic environment
dissipation can be described by a Lindblad superoperator
$\mathcal{L}[\sigma_{-}]$ and a dissipation rate
$\gamma_{0}=\omega^{3}_{0}d^{2}_{j}/(3\pi\epsilon_{0}\hbar c^{3})$, that is,
$\gamma_{0}\mathcal{L}[\sigma_{-}]$.

##  VII Polaron frame correlation operators

Without going into the explicit form of every term of the HOM $g^{(2)}(\tau)$
in the polaron frame, we note that a general correlation operator of the form
$\Braket{\sigma_{+}(t_{1})\sigma_{+}(t_{2})\sigma_{-}(t_{3})\sigma_{-}(t_{4})}$
in the polaron frame is given by

$\displaystyle\begin{split}&\Braket{\sigma_{+}(t_{1})B_{+}(t_{1})\sigma_{+}(t_{2})B_{+}(t_{2})\sigma_{-}(t_{3})B_{-}(t_{3})\sigma_{-}(t_{4})B_{-}(t_{4})}\\\
&\hskip
28.45274pt\approx\Braket{B_{+}(t_{1})B_{+}(t_{2})B_{-}(t_{3})B_{-}(t_{4})}\\\
&\hskip
42.67912pt\times\Braket{\sigma_{+}(t_{1})\sigma_{+}(t_{2})\sigma_{-}(t_{3})\sigma_{-}(t_{4})}\leavevmode\nobreak\
,\end{split}$ (S32)

where the approximation in the second holds due to the difference in
timescales associated with the exciton lifetime (nanoseconds) and the phonon
bath relaxation time (picoseconds) Iles-Smith  _et al._ (2017). The phonon
correlation operator can be easily expanded as

$\displaystyle\begin{split}&\Braket{B_{+}(t_{1})B_{+}(t_{2})B_{-}(t_{3})B_{-}(t_{4})}\\\
&=\Braket{\mathcal{B}^{t_{1},t_{2}}_{t_{3},t_{4}}}\\\
&\quad\times\exp\left[i\Im\\{\phi(t_{2}-t_{1})\\}\right]\exp\left[i\Im\\{\phi(t_{4}-t_{3})\\}\right]\\\
&\quad\times\exp\left[i\Im\\{\phi(t_{1}+t_{2}-t_{3}-t_{4})\\})\\}\right]\leavevmode\nobreak\
,\end{split}$ (S33)

where $\Im$ denotes the imaginary part, and
$\Braket{\mathcal{B}^{t_{1},t_{2}}_{t_{3},t_{4}}}$ is given by

$\displaystyle\begin{split}&\Braket{\mathcal{B}^{t_{1},t_{2}}_{t_{3},t_{4}}}\\\
&=\exp\left[-\frac{1}{2}\int^{\infty}_{0}\frac{J_{pn}(\omega)}{\omega^{2}}|\mathcal{K}^{t_{1},t_{2}}_{t_{3},t_{4}}|^{2}\coth(\beta\omega/2)\right]\leavevmode\nobreak\
,\end{split}$ (S34)

where $\mathcal{K}^{t_{1},t_{2}}_{t_{3},t_{4}}=\mathrm{e}^{i\omega
t_{1}}+\mathrm{e}^{i\omega t_{2}}-\mathrm{e}^{i\omega
t_{3}}-\mathrm{e}^{i\omega t_{4}}$. These expressions can be simplified
further by noting that the phonon propagator vanishes on a picosecond
timescale, resulting in an easier numerical evaluation of Eq. (S33).

##  VIII Resonance fluorescence spectrum

In this section, we outline the derivation of the RF spectral function in the
presence of strong vibrational coupling. In the polaron frame, this quantity
is simply given by the Fourier transform of the (steady-state) first order
correlation function
$\mathrm{lim}_{t\rightarrow\infty}\langle\mathbf{E}^{(-)}(\mathbf{R},t)\mathbf{E}^{(+)}(\mathbf{R},t+\tau)\rangle$,
where $\mathbf{E}^{(-)}(\mathbf{R},t)$ and $\mathbf{E}^{(+)}(\mathbf{R},t)$
are, respectively, the negative and positive components of the electric field
operator evaluated at the position $\mathbf{R}$ of the detector Ficek  _et
al._ (2005b).

$\displaystyle\begin{split}S(\omega)\propto\int_{-\infty}^{\infty}\mathrm{d}\tau&\sum_{j\in{\uparrow,\downarrow}}\mathrm{e}^{-i(\omega-\omega^{\prime})\tau}\times\\\
&\langle\sigma^{j}_{+}(\tau)B_{+}(\tau)\sigma^{j}_{-}(0)B_{-}(0)\rangle_{s}\leavevmode\nobreak\
,\\\ \end{split}$ (S35)

where we have exploited the temporal homogeneity of the stationary correlation
function, and where the subscript ‘s’ denotes the trace taken with respect the
steady-state density matrix. Due to the different timescales Scerri  _et al._
(2017) for the phonon and photon processes, Eq. (S35) simplifies to

$\displaystyle\begin{split}S(\omega)\propto\langle
B\rangle^{2}\int_{-\infty}^{\infty}\mathrm{d}\tau&\sum_{j\in{\uparrow,\downarrow}}\mathrm{e}^{-i(\omega-\omega^{\prime})\tau}\mathrm{e}^{\phi(\tau)}\times\\\
&\langle\sigma^{j}_{+}(\tau)\sigma^{j}_{-}(0)\rangle_{s}\leavevmode\nobreak\
.\\\ \end{split}\leavevmode\nobreak\ .$ (S36)

Similarly, for the neutral exciton, the spectral function is given by

$\displaystyle\begin{split}S(\omega)\propto\langle
B\rangle^{2}\int_{-\infty}^{\infty}\mathrm{d}\tau&\mathrm{e}^{-i(\omega-\omega^{\prime})\tau}\mathrm{e}^{\phi(\tau)}\times\\\
&\langle\sigma_{+}(\tau)\sigma_{-}(0)\rangle_{s}\leavevmode\nobreak\ .\\\
\end{split}\leavevmode\nobreak\ $ (S37)

## References

  * Hong  _et al._ (1987) C. K. Hong, Z. Y. Ou, and L. Mandel, Physical Review Letters 59, 2044 (1987).
  * Proux  _et al._ (2015) R. Proux, M. Maragkou, E. Baudin, C. Voisin, P. Roussignol, and C. Diederichs, Phys. Rev. Lett. 114, 067401 (2015).
  * Scully and Zubairy (1997) M. O. Scully and M. S. Zubairy,  _Quantum optics_ (Cambridge University Press, Cambridge, 1997).
  * Legero  _et al._ (2004) T. Legero, T. Wilk, M. Hennrich, G. Rempe, and A. Kuhn, Physical Review Letters 93, 070503 (2004).
  * Baudin  _et al._ (2019) E. Baudin, R. Proux, M. Maragkou, P. Roussignol, and C. Diederichs, Physical Review A 99, 013842 (2019).
  * Giesz  _et al._ (2015) V. Giesz, S. L. Portalupi, T. Grange, C. Antón, L. De Santis, J. Demory, N. Somaschi, I. Sagnes, A. Lemaître, L. Lanco, A. Auffèves, and P. Senellart, Phys. Rev. B 92, 161302 (2015).
  * Iles-Smith  _et al._ (2017) J. Iles-Smith, D. P. S. McCutcheon, J. Mørk, and A. Nazir, Physical Review B 95, 201305 (2017).
  * Ficek  _et al._ (2005a) Z. Ficek, T. Asakura, S. Swain, K. Brenner, T. Hänsch, T. Kamiya, F. Krausz, B. Monemar, W. Rhodes, H. Venghaus, _et al._ ,  _Quantum Interference and Coherence: Theory and Experiments_, Springer Series in Optical Sciences (Springer, 2005).
  * Lebreton  _et al._ (2013) A. Lebreton, I. Abram, R. Braive, I. Sagnes, I. Robert-Philip, and A. Beveratos, Phys. Rev. A 88, 013801 (2013).
  * Urbaszek  _et al._ (2013) B. Urbaszek, X. Marie, T. Amand, O. Krebs, P. Voisin, P. Maletinsky, A. Högele, and A. Imamoğlu, Reviews of Modern Physics 85, 79 (2013).
  * Malein  _et al._ (2016) R. Malein, T. Santana, J. Zajac, A. Dada, E. Gauger, P. Petroff, J. Lim, J. Song, and B. Gerardot, Physical Review Letters 116, 257401 (2016).
  * Fernandez  _et al._ (2009) G. Fernandez, T. Volz, R. Desbuquois, A. Badolato, and A. Imamoğlu, Physical Review Letters 103, 087406 (2009).
  * He  _et al._ (2013) Y. He, Y.-M. He, Y.-J. Wei, X. Jiang, M.-C. Chen, F.-L. Xiong, Y. Zhao, C. Schneider, M. Kamp, S. Höfling, C.-Y. Lu, and J.-W. Pan, Physical Review Letters 111, 237403 (2013).
  * Sun  _et al._ (2016) Z. Sun, A. Delteil, S. Faelt, and A. Imamoğlu, Physical Review B 93, 241302 (2016).
  * Mahan (2000) G. Mahan,  _Many-Particle Physics_, Physics of Solids and Liquids (Springer, 2000).
  * Lang and Firsov (1962) I. Lang and Y. A. Firsov, Journal of Experimental and Theoretical Physics 43, 1301 (1962).
  * Nazir and McCutcheon (2016) A. Nazir and D. P. S. McCutcheon, Journal of Physics: Condensed Matter 28, 103002 (2016).
  * Scerri  _et al._ (2017) D. Scerri, T. S. Santana, B. D. Gerardot, and E. M. Gauger, Phys. Rev. B 95, 165403 (2017).
  * Ramsay  _et al._ (2010) A. J. Ramsay, A. V. Gopal, E. M. Gauger, A. Nazir, B. W. Lovett, A. M. Fox, and M. S. Skolnick, Physical Review Letters 104, 017402 (2010).
  * Roy-Choudhury and Hughes (2015a) K. Roy-Choudhury and S. Hughes, Opt. Lett. 40, 1838 (2015a).
  * Roy-Choudhury and Hughes (2015b) K. Roy-Choudhury and S. Hughes, Phys. Rev. B 92, 205406 (2015b).
  * Breuer and Petruccione (2007) H. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (OUP Oxford, 2007).
  * Ficek  _et al._ (2005b) Z. Ficek, T. Asakura, S. Swain, K. Brenner, T. Hänsch, T. Kamiya, F. Krausz, B. Monemar, W. Rhodes, H. Venghaus, _et al._ ,  _Quantum Interference and Coherence: Theory and Experiments_, Springer Series in Optical Sciences (Springer, 2005).