Probing bath-induced entanglement in a qubit pair by measuring photon correlations

Self-assembled quantum dots (QDs), quantum heterostructures in which electrons and holes are confined in all three dimensions, are artificial solid-state atoms with tailored optical and electronic coherence properties. Impressive progress on fabrication and optical manipulation techniques has enabled high fidelity preparation, control and readout of the quantum states of charge carriers confined in individual QDs [1–8] and in QD pairs [9–11]. Indeed, entanglement of photons and carrier spins has been fully characterized [12], and QDs are now recognized as superb single and entangled photon sources [13, 14]. At the same time, the coupling of QDs to their solid-state environment [15–17] provides a rich platform for the study of open system effects that may be more difficult to observe in, for instance, atomic systems.

There is considerable evidence that the decoherence effects induced by the solid state environment of QD excitons are dominated by interactions with longitudinal acoustic phonons via the deformation potential coupling [15, 16]. The QD interaction with the phonon bath leads to pure dephasing of bare excitons, or relaxation of driven excitons, in individual QDs [18–20]—but a range of not-yet-observed phenomena are predicted to appear due to the interaction of multiple QDs with a phonon bath. For example, phonon induced dephasing results in entanglement decay between two QDs at a much faster rate than the individual QD dephasing rate [21, 22]. Moreover, there has been considerable interest in phonon-assisted processes that appear when coupled QDs interact with the same phonon bath. Phonon assisted tunneling [23–25], relaxation [26, 27] and excitation transfer [10, 28, 29] have been demonstrated. However, in this paper we focus on an investigation of the phenomena that appear solely due to the collective interaction of un-coupled QDs with a common phonon reservoir, which cannot be explained by an interaction with separate reservoirs.

The properties of the phonon bath are often investigated through excitonic occupation dynamics [30–33]. Here we will study another powerful method—the analysis of emitted photon statistics. Owing to the strong optical transition dipole of semiconductor QDs, the optical properties of individual QDs [6, 34–37], as well as ensembles of QDs [38–40] have been studied extensively with this technique. For example, the single qubit second-order fluorescence intensity correlation function ${{g}^{(2)}}$ has been investigated both experimentally [41] and theoretically [42] and has been shown to yield important information about the nature of the QD solid-state environment. Furthermore, the two QD intensity correlation function has been used to characterize the coupling mechanism between two adjacent dots [10, 43, 44].

In this paper we will show that photon statistics measurements can be used to probe the QDs' immediate environment and find signatures of a common environment in a pair of un-coupled, driven QDs. We will describe how the interaction with a common phonon field results in both a coherent elastic coupling mediated by virtual phonons and an inelastic incoherent coupling with emission and absorption of real phonons. Although the interaction of excitons with phonons usually results in decoherence and entanglement decay we find that when driven dots interact with the same environment this interaction can be used to entangle the QDs, even in the steady state. This complements previous work in which environment-induced entanglement of undriven, and unseparated atoms was seen to persist ad infinitum [45]; for our QD case the assumption of zero separation is obviously unrealistic. A further study of undriven, spatially separated dots found environment-induced entanglement decays away to zero, albeit on very long timescales [46]. We find that intensity correlation measurements of emitted photons provide a signature of the common environment and can be used to measure the strength of the coherent and incoherent coupling mechanisms.

In the following section we present the Hamiltonian of the QDs coupled to the radiation and solid-state environment. In section 3 we trace out the radiation and phonon bath in a Born–Markov fashion to obtain a second-order master equation in Lindblad form for the reduced density matrix of the two QDs. We go on in section 4 to investigate the steady state solution of the master equation showing that the common phonon bath induces entanglement, before in section 5 finding signatures of a common environment in the normalized intensity correlation function ${{g}^{(2)}}$. We conclude in section 6.

The system we consider is formed of two driven QDs that do not interact with each other directly, but which are coupled to the same phonon and radiation bath. We will consistently use the parameters of a GaAs self-assembled exciton QD, although all the calculations can apply to any optical emitters in a solid state environment, so long as the approximations we used are valid. In this section we obtain the Hamiltonian for our system following closely the derivation in [47].

We model a single QD as a two level system with ground state $|0\rangle$ and excited state $|\xi \rangle$ separated by an energy difference ${{\omega }_{\xi }}$. We define the creation operator for an exciton as ${{c}^{\dagger }}=|\xi \rangle \langle 0|$ with the annihilation operator its Hermitian conjugate. We will denote an operator O acting on the Hilbert space of the first dot (second dot) as ${{O}_{A}}\equiv O\otimes I$ (${{O}_{B}}\equiv I\otimes O$) or, when notational clarity demands, as ${{O}^{A}}\equiv O\otimes I$ (${{O}^{B}}\equiv I\otimes O$). Each dot is driven by its own near-resonant laser of frequency ω_l^j which couples to each dot with a strength leading to a Rabi frequency Ω^j with $j\in \{A,B\}$.

The dots are also coupled to the electromagnetic environment, which we represent as a bath of harmonic oscillators of frequencies ${{\Theta }_{{\boldsymbol{q}} }}$ and creation operators $a_{{\boldsymbol{q}} }^{\dagger }$, where ${\boldsymbol{q}}$ denotes a mode of the radiation field with wave vector ${\boldsymbol{q}}$. The coupling Hamiltonian is:

$$\begin{eqnarray}&&{{H}_{I}}=\mathop{\sum }\limits_{j\in \{A,B\}}(c_{j}^{\dagger }+{{c}_{j}})\mathop{\sum }\limits_{{\boldsymbol{q}} }\;f_{{\boldsymbol{q}} }^{j}(a_{{\boldsymbol{q}} }^{\dagger }+{{a}_{{\boldsymbol{q}} }}).\end{eqnarray} \tag{ 1 }$$

The interaction between the dots and the radiation field is fully characterized by the photon spectral density function:

$$\begin{eqnarray}&&J_{\gamma }^{j}(\omega )\equiv 2\pi \mathop{\sum }\limits_{{\boldsymbol{q}} }|f_{{\boldsymbol{q}} }^{j}{{|}^{2}}\delta (\omega -{{\omega }_{{\boldsymbol{q}} }}).\end{eqnarray} \tag{ 2 }$$

For our calculation the spectral density will not vary significantly across the relevant frequencies and so we can consider it a constant, $J_{\gamma }^{j}(\omega )\approx 1/{{T}^{*}}$, with ${{T}^{*}}\;\approx$ 1 ns.

QDs exist in a solid state environment, and so are also coupled to a common phonon bath. This is again represented as a collection of harmonic oscillators of frequencies ${{\omega }_{{\boldsymbol{k}} }}$ and creation operators $b_{{\boldsymbol{k}} }^{\dagger }$, where ${\boldsymbol{k}}$ denotes the phonon wave vector. By assuming the strong confinement limit, we may model the electron–hole wave functions ${{\psi }_{e,h}}$ as (unentangled) product states of the corresponding single particle wave functions. The essential physics of QDs may then be captured by assuming Gaussian single particle spatial wave functions with standard deviation lengths d_e and d_h for electrons and holes respectively. The interaction of such states with longitudinal acoustic phonons coupled through the deformation potential is dominant [15], and thus we obtain the following exciton–phonon interaction Hamiltonian:

$$\begin{eqnarray}&&{{H}_{I}}=\mathop{\sum }\limits_{j\in \{A,B\}}c_{j}^{\dagger }{{c}_{j}}\mathop{\sum }\limits_{{\boldsymbol{k}} }t_{{\boldsymbol{k}} }^{j}(b_{{\boldsymbol{k}} }^{\dagger }+{{b}_{{\boldsymbol{k}} }}),\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&t_{{\boldsymbol{k}} }^{A}={\rm i}\sqrt{\frac{1}{2\mu V{{\omega }_{{\boldsymbol{k}} }}}}|{\boldsymbol{k}} |\left[ {{{\sf D} }_{e}}{\rm exp} \left( -\frac{d_{e}^{2}|{\boldsymbol{k}} {{|}^{2}}}{4} \right)-{{{\sf D} }_{h}}{\rm exp} \left( -\frac{d_{h}^{2}|{\boldsymbol{k}} {{|}^{2}}}{4} \right) \right]\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&t_{{\boldsymbol{k}} }^{B}=t_{{\boldsymbol{k}} }^{A}{{{\rm e}}^{{\rm i}{\boldsymbol{k}} \cdot {\boldsymbol{d}} }},\end{eqnarray} \tag{ 5 }$$

where μ is the GaAs mass density, V is the volume of the crystal, ${{{\sf D} }_{e,h}}$ are the electron/hole deformation potential constants. We are also able to define a phonon spectral density:

$$\begin{eqnarray}&&{{J}_{p}}(\omega )\equiv 2\pi \mathop{\sum }\limits_{{\boldsymbol{k}} }|t_{{\boldsymbol{k}} }^{j}{{|}^{2}}\delta (\omega -{{\omega }_{{\boldsymbol{k}} }})={{\Lambda }_{\mathbb{D}}}{{\left( \frac{\omega }{{{\omega }_{0}}} \right)}^{\mathbb{D}}}{{\mathcal{P}}^{2}}(\omega ),\end{eqnarray} \tag{ 6 }$$

where ω₀ is a scaling energy introduced for convenience and its value defined in table 1, along with all the other parameters we have used. $\mathbb{D}=1,2,3$ is the dimensionality of the phonons and ${{\mathcal{P}}^{2}}(\omega )$ is the form factor given by:

$$\begin{eqnarray}&&{{\mathcal{P}}^{2}}(\omega )={{\left( \frac{1}{{{{\sf D} }_{{\rm e}}}-{{{\sf D} }_{h}}} \right)}^{2}}\left( {\sf D} _{h}^{2}{{{\rm e}}^{-{{\omega }^{2}}/\omega _{h}^{2}}}+{\sf D} _{e}^{2}{{e}^{-{{\omega }^{2}}/\omega _{e}^{2}}}-2{{{\sf D} }_{h}}{{{\sf D} }_{e}}{{{\rm e}}^{-{{\omega }^{2}}/\omega _{eh}^{2}}} \right),\end{eqnarray} \tag{ 7 }$$

where ${{\omega }_{\alpha }}={{c}_{s}}\sqrt{2}/{{d}_{\alpha }}$ for $\alpha \in \{e,h\}$ and ${{\omega }_{eh}}=2{{c}_{s}}/\sqrt{d_{e}^{2}+d_{h}^{2}}$. The dimensionality dependent constant ${{\Lambda }_{\mathbb{D}}}$ is given by:

$$\begin{eqnarray}\frac{{{\Lambda }_{\mathbb{D}}}}{\omega _{0}^{\mathbb{D}}}=\left\{ \begin{array}{ccccccccccccccc} {{({{{\sf D} }_{e}}-{{{\sf D} }_{h}})}^{2}}/(2\pi {{\mu }_{3}}c_{s}^{5}), & {\rm when}\ \quad \mathbb{D}=3, \\ {{({{{\sf D} }_{e}}-{{{\sf D} }_{h}})}^{2}}/(2{{\mu }_{2}}c_{s}^{4}), & {\rm when}\ \quad \mathbb{D}=2, \\ {{({{{\sf D} }_{e}}-{{{\sf D} }_{h}})}^{2}}/({{\mu }_{1}}c_{s}^{3}), & {\rm when}\ \quad \mathbb{D}=1. \\ \end{array} \right.\end{eqnarray} \tag{ 8 }$$

${{\mu }_{1,2,3}}$ are the mass densities for different dimensionalities, and c_s is the speed of sound.

Table 1.  Properties for GaAs QDs.

Electron deformation potential ${{{\sf D} }_{e}}$	14.6 eV
Hole deformation potential ${{{\sf D} }_{h}}$	4.8 eV
Mass density in 3D μ3	$5.0\times {{10}^{3}}$ kg m−3
Mass density in 2D μ2	$3.0\times {{10}^{-6}}$ kg m−2
Mass density in 1D μ1	$1.7\times {{10}^{-15}}$ kg m−1
Velocity of sound cs	$5.11\times {{10}^{5}}$ cm s−1
Electron/hole wave function size $d_{e,h}^{j}$	5 nm
First QD bare energy WA	1.34 eV
Second QD bare energy WB	1.38 eV
Electron/hole cutoff frequency ${{\omega }_{e,h,eh}}$	1.05 meV
Electron–phonon coupling strength in 3D c3	0.46 meV
Electron–phonon coupling strength in 2D c2	8.7 meV
Electron–phonon coupling strength in 1D c1	103 meV
Scaling energy ω0	1 meV
Photon timescale T*	1 ns

The resulting Hamiltonian takes the form $(\hbar =1):$

$$\begin{eqnarray}\begin{array}{rcl} H & = & \mathop{\sum }\limits_{j\in \{A,B\}}\left[ \omega _{\xi }^{j}c_{j}^{\dagger }{{c}_{j}}+\frac{1}{2}{{\Omega }^{j}}{\rm cos} (\omega _{l}^{j}t)(c_{j}^{\dagger }+{{c}_{j}}) \right]+\mathop{\sum }\limits_{{\boldsymbol{q}} }{{\Theta }_{{\boldsymbol{q}} }}a_{{\boldsymbol{q}} }^{\dagger }{{a}_{{\boldsymbol{q}} }}+\mathop{\sum }\limits_{{\boldsymbol{k}} }{{\omega }_{{\boldsymbol{k}} }}b_{{\boldsymbol{k}} }^{\dagger }{{b}_{{\boldsymbol{k}} }} \\ {} & {} & +\mathop{\sum }\limits_{j\in \{A,B\}}(c_{j}^{\dagger }+{{c}_{j}})\mathop{\sum }\limits_{{\boldsymbol{q}} }g_{{\boldsymbol{q}} }^{j}(a_{{\boldsymbol{q}} }^{\dagger }+{{a}_{{\boldsymbol{q}} }})+\mathop{\sum }\limits_{j\in \{A,B\}}c_{j}^{\dagger }{{c}_{j}}\mathop{\sum }\limits_{{\boldsymbol{k}} }t_{{\boldsymbol{k}} }^{j}(b_{{\boldsymbol{k}} }^{\dagger }+{{b}_{{\boldsymbol{k}} }}). \\ \end{array}\end{eqnarray} \tag{ 9 }$$

In order to eliminate the time dependence in the QD part of the Hamiltonian we move to a rotating frame by applying the unitary transform $U={\rm exp} ({\rm i}{{H}_{0}}t)$ with ${{H}_{0}}={{\sum }_{j\in \{A,B\}}}\omega _{l}^{j}c_{j}^{\dagger }{{c}_{j}}$. After performing a rotating-wave approximation on the driving term and on the QD-radiative bath interaction, the transformed Hamiltonian is:

$$\begin{eqnarray}\begin{array}{rcl} \tilde{H} & = & \mathop{\sum }\limits_{j}\left[ \Delta _{e}^{j}\tilde{c}_{j}^{\dagger }{{\tilde{c}}_{j}}+\frac{1}{2}{{\Omega }^{j}}(\tilde{c}_{j}^{\dagger }+{{\tilde{c}}_{j}}) \right]+\mathop{\sum }\limits_{{\boldsymbol{q}} }{{\Theta }_{{\boldsymbol{q}} }}a_{{\boldsymbol{q}} }^{\dagger }{{a}_{{\boldsymbol{q}} }}+\mathop{\sum }\limits_{{\boldsymbol{k}} }{{\omega }_{{\boldsymbol{k}} }}b_{{\boldsymbol{k}} }^{\dagger }{{b}_{{\boldsymbol{k}} }} \\ {} & {} & +\mathop{\sum }\limits_{j}\mathop{\sum }\limits_{{\boldsymbol{q}} }g_{{\boldsymbol{q}} }^{j}({{\tilde{c}}_{j}}{{{\rm e}}^{-{\rm i}\omega _{l}^{j}}}a_{{\boldsymbol{q}} }^{\dagger }+\tilde{c}_{j}^{\dagger }{{{\rm e}}^{{\rm i}\omega _{l}^{j}t}}{{a}_{{\boldsymbol{q}} }})+\mathop{\sum }\limits_{j}\tilde{c}_{j}^{\dagger }{{\tilde{c}}_{j}}\mathop{\sum }\limits_{{\boldsymbol{k}} }t_{{\boldsymbol{k}} }^{j}(b_{{\boldsymbol{k}} }^{\dagger }+{{b}_{{\boldsymbol{k}} }}) \\ \end{array}\end{eqnarray} \tag{ 10 }$$

where we have denoted the interaction picture creation operators with $\tilde{c}_{j}^{\dagger }={{{\rm e}}^{-{\rm i}\omega _{l}^{j}}}c_{j}^{\dagger }$ and $\Delta _{e}^{j}=\omega _{\xi }^{j}-\omega _{l}^{j}$. We have now dropped the explicit $\{A,B\}$ indices over which j is summed, but this will henceforth always be assumed.

Since the two QDs do not interact with each other directly we can diagonalize the DQD part of the Hamiltonian by diagonalizing each QD part separately. For each QD the resulting eigenenergies are ${{W}_{j}}=\sqrt{{{(\Delta _{e}^{j})}^{2}}+{{({{\Omega }^{j}})}^{2}}}$; the resulting eigenstates are $|e{{\rangle }_{j}}={\rm cos} ({{\theta }_{j}}/2)|\xi {{\rangle }_{j}}\;+{\rm sin} ({{\theta }_{j}}/2)|0{{\rangle }_{j}}$ and $|g{{\rangle }_{j}}=-{\rm sin} ({{\theta }_{j}}/2)|\xi {{\rangle }_{j}}+{\rm cos} ({{\theta }_{j}}/2)|0{{\rangle }_{j}}$, and we have introduced the mixing angles ${{\theta }_{j}}={\rm arccos} (\Delta _{e}^{j}/{{W}_{j}})$. Using the Pauli spin notation (${{\sigma }_{+}}=|e\rangle \langle g|$, ${{\sigma }_{z}}=|e\rangle \langle e|-|g\rangle \langle g|$, $I=|e\rangle \langle e|+|g\rangle \langle g|$), we may now write the Hamiltonian as:

$$\begin{eqnarray}\begin{array}{rcl} \tilde{H} & = & \mathop{\sum }\limits_{j}\frac{{{W}_{j}}}{2}\sigma _{z}^{j}+\mathop{\sum }\limits_{{\boldsymbol{q}} }{{\Theta }_{{\boldsymbol{q}} }}a_{{\boldsymbol{q}} }^{\dagger }{{a}_{{\boldsymbol{q}} }}+\mathop{\sum }\limits_{{\boldsymbol{k}} }{{\omega }_{{\boldsymbol{k}} }}b_{{\boldsymbol{k}} }^{\dagger }{{b}_{{\boldsymbol{k}} }} \\ {} & {} & +\mathop{\sum }\limits_{j}\mathop{\sum }\limits_{{\boldsymbol{q}} }g_{{\boldsymbol{q}} }^{j}\left\{ \left[ \sigma _{+}^{j}{{{\rm cos} }^{2}}\frac{{{\theta }_{j}}}{2}-\sigma _{-}^{j}{{{\rm sin} }^{2}}\frac{{{\theta }_{j}}}{2}+\frac{{\rm sin} {{\theta }_{j}}}{2}\sigma _{z}^{j} \right]{{{\rm e}}^{{\rm i}\omega _{l}^{j}t}}{{a}_{{\boldsymbol{q}} }}+{\rm h}.{\rm c}. \right\} \\ {} & {} & +\mathop{\sum }\limits_{j}\left\{ -\frac{{\rm sin} {{\theta }_{j}}}{2}(\sigma _{-}^{j}+\sigma _{+}^{j})+\frac{{\rm cos} {{\theta }_{j}}}{2}\sigma _{z}^{j}+\frac{1}{2}{{I}^{j}} \right\}\mathop{\sum }\limits_{{\boldsymbol{k}} }t_{{\boldsymbol{k}} }^{j}(b_{{\boldsymbol{k}} }^{\dagger }+{{b}_{{\boldsymbol{k}} }}). \\ \end{array}\end{eqnarray} \tag{ 11 }$$

While a perturbative treatment of the radiation bath under the Born–Markov approximation is enough to accurately describe the effect of the radiation field, a range of theoretical methods have been developed to investigate the how phonon interactions affects the system dynamics. In the limit of weak coupling, one can do perturbative expansions of the QD–phonon coupling, resulting in master equation descriptions of both Markovian [15, 42, 48, 49] and non-Markovian [30, 33, 50] nature, as well as correlation expansions [31, 51, 52]. The polaron transform [53] in conjunction with a perturbative expansion in the polaron-transformed basis can account for various non-perturbative effects not captured by the weak-coupling treatment [54–56]. Non-perturbative numerically exact techniques that rely on the calculation of the path integral have also been implemented [32, 57].

Polaronic effects in the QD dynamics are smaller at lower temperatures ($T\lt 30$ K) [56] and weaker coupling, and, as we will show later in figures 5 and 6, the phonon-induced entanglement and correlations in photon emission are most pronounced at low temperatures. We are therefore able to treat both the radiation and phonon bath in a Born–Markov fashion, discussing future possible improvements on this approximation in our conclusion. The general form of the resulting second order master equation for the reduced density matrix of a system ρ is [58]:

$$\begin{eqnarray}&&\dot{\rho }(t)=-{\rm i}[{{H}_{0}},\rho (t)]-\int _{0}^{\infty }{\rm d}\tau {\rm t}{{{\rm r}}_{B}}\left\{ \left[ {{H}_{I}},\left[ {{{\rm e}}^{{\rm i}({{H}_{0}}+{{H}_{B}})\tau }}{{H}_{I}}{{{\rm e}}^{-{\rm i}({{H}_{0}}+{{H}_{B}})\tau }},\rho (t)\otimes {{\rho }_{B}} \right] \right] \right\}.\end{eqnarray} \tag{ 12 }$$

with ${{H}_{0}}({{H}_{B}})$ the system (bath) Hamiltonian and H_I the system–bath interaction; ${{\rho }_{B}}$ is the time-independent bath density matrix.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can apply this general formulation to our specific case: it is straightforward to trace out the phonon and radiation reservoirs and obtain a master equation in Lindblad form (for details see the appendix). For the case of three-dimensional (3D) phonons we obtain:

$$\begin{eqnarray}\begin{array}{rcl} {\dot{\rho }} & = & -{\rm i}\left[ H^{\prime} +{{H}_{X}},\rho \right]+\mathop{\sum }\limits_{j}{{\mathcal{L}}_{j}}[\rho ]+{{\mathcal{L}}_{X}}[\rho ], \\ H^{\prime} & = & \frac{W_{A}^{\prime }}{2}\sigma _{z}^{A}+\frac{W_{B}^{\prime }}{2}\sigma _{z}^{B}, \\ {{H}_{X}} & = & \frac{1}{2}{{t}_{X}}(\sigma _{+}^{A}\sigma _{-}^{B}+\sigma _{-}^{A}\sigma _{+}^{B})+\frac{{{\Delta }_{X}}}{2}\sigma _{z}^{A}\sigma _{z}^{B}, \\ {{\mathcal{L}}_{j}}[\rho ] & = & \left( \Gamma _{\downarrow }^{j}-|\Gamma _{\downarrow }^{X}| \right)\mathcal{D}(\sigma _{-}^{j},\rho )+\left( \Gamma _{\uparrow }^{j}-|\Gamma _{\uparrow }^{X}| \right)\mathcal{D}(\sigma _{+}^{j},\rho ), \\ {{\mathcal{L}}_{X}}[\rho ] & = & |\Gamma _{\downarrow }^{X}|\mathcal{D}(\sigma _{-}^{A}\pm \sigma _{-}^{B},\rho )+|\Gamma _{\uparrow }^{X}|\mathcal{D}(\sigma _{+}^{A}\pm \sigma _{+}^{B},\rho ), \\ \end{array}\end{eqnarray} \tag{ 13 }$$

where we define the dissipators $\mathcal{D}(x,\rho )\equiv x\rho {{x}^{\dagger }}-1/2\rho {{x}^{\dagger }}x-1/2{{x}^{\dagger }}x\rho$, and in ${{\mathcal{L}}_{X}}[\rho ]$ the plus (minus) sign is chosen for positive (negative) rates. The various new parameters are defined as:

$$\begin{eqnarray*}\begin{array}{rcl} {{W}_{j}}^{\prime} & = & {{W}_{j}}+{{\Delta }_{j}}, \\ {{\Delta }_{j}} & = & \frac{{{{\rm sin} }^{2}}{{\theta }_{j}}}{2}\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \frac{[2\bar{n}(\omega )+1]{{W}_{j}}}{{{\omega }^{2}}-W_{j}^{2}}\frac{{{J}_{p}}(\omega )}{2\pi }+{\rm cos} {{\theta }_{j}}\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \frac{1}{\omega }\frac{{{J}_{p}}(\omega )+{{J}_{X}}(\omega )}{2\pi }, \\ {{\Delta }_{X}} & = & {\rm cos} {{\theta }_{A}}{\rm cos} {{\theta }_{B}}\mathcal{P}\int _{0}^{\infty }\frac{1}{\omega }\frac{{{J}_{X}}(\omega )}{2\pi }, \\ {{t}_{X}} & = & \frac{{\rm sin} {{\theta }_{A}}{\rm sin} {{\theta }_{B}}}{2}\mathcal{P}\int _{0}^{\infty }\left( \frac{1}{{{\omega }^{2}}-W_{A}^{2}}+\frac{1}{{{\omega }^{2}}-W_{B}^{2}} \right)\frac{\omega {{J}_{X}}(\omega )}{2\pi }, \\ \Gamma _{\downarrow }^{j} & = & J_{\gamma }^{j}(\omega _{l}^{j}){{{\rm cos} }^{4}}\frac{{{\theta }_{j}}}{2}+\frac{{{{\rm sin} }^{2}}{{\theta }_{j}}}{4}{{J}_{p}}({{W}_{j}})[\bar{n}({{W}_{j}})+1], \\ \Gamma _{\uparrow }^{j} & = & J_{\gamma }^{j}(\omega _{l}^{j}){{{\rm sin} }^{4}}\frac{{{\theta }_{j}}}{2}+\frac{{{{\rm sin} }^{2}}{{\theta }_{j}}}{4}{{J}_{p}}({{W}_{j}})\bar{n}({{W}_{j}}), \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{rcl} \Gamma _{\phi }^{j} & = & J_{\gamma }^{j}(\omega _{l}^{j})\frac{{{{\rm sin} }^{2}}{{\theta }_{j}}}{4}, \\ \Gamma _{\downarrow }^{X} & = & \frac{{\rm sin} {{\theta }_{A}}{\rm sin} {{\theta }_{B}}}{8}\left\{ {{J}_{X}}({{W}_{A}})[\bar{n}({{W}_{A}})+1]+{{J}_{X}}({{W}_{B}})[\bar{n}({{W}_{B}})+1] \right\}, \\ \Gamma _{\uparrow }^{X} & = & \frac{{\rm sin} {{\theta }_{A}}{\rm sin} {{\theta }_{B}}}{8}\left\{ {{J}_{X}}({{W}_{A}})\bar{n}({{W}_{A}})+{{J}_{X}}({{W}_{B}})\bar{n}({{W}_{B}}) \right\}, \\ \end{array}\end{eqnarray} \tag{ 14 }$$

where $\bar{n}(\omega )=1/({\rm exp} [\omega /({{k}_{B}}T)]-1)$ is the average phonon number at temperature T and frequency ω. In addition to the phonon and photon spectral density function we have defined the spectral density function of the common phonon field as:

$$\begin{eqnarray}&&{{J}_{X}}(\omega )\equiv 2\pi \mathop{\sum }\limits_{{\boldsymbol{k}} }|{{t}_{{\boldsymbol{k}} }}{{|}^{2}}{{{\rm e}}^{{\rm i}{\boldsymbol{k}} \cdot {\boldsymbol{d}} }}\delta ({{\omega }_{{\boldsymbol{k}} }}-\omega )=F(\omega d/{{c}_{s}}){{J}_{p}}(\omega ),\end{eqnarray} \tag{ 15 }$$

where c_s is the speed of sound in GaAs and F(x) is a function that, in the case of 3D phonons, has the form $F(x)={\rm sinc}(x)$.

The master equation (13) is fully physical—i.e. it is in Lindblad form with all rates positive—as long as $\Gamma _{\downarrow ,\uparrow }^{X}\lt \Gamma _{\downarrow }^{j}$. This is always the case since ${{F}_{D}}(x)\leqslant 1$. In the above equations we have isolated the effects that result solely due to the interaction of the QDs with a common phonon bath from the effects that would still appear even when the dots interact with separate phonon baths. We have denoted the corresponding Hamiltonian, dissipators and rates due to this common bath with a subscript X since they only appear when there is a finite cross correlation between the local phonon environments of the two QDs.

The effects that appear when the QDs interact with independent, separate phonon and photon baths are already well-known. The master equation for this case can be obtained from equation (13) by setting ${{J}_{X}}(\omega )=0$, which corresponds physically to very distant or off resonant QDs. This leads first to a renormalization of the QD energies by Δ_j, which has two contributions. One (whose size is proportional to ${{{\rm sin} }^{2}}\;{{\theta }_{j}}$) is most effective for resonantly driven dots and disappears for far-detuned driving lasers; the other (proportional to cos θ_j) is negligible in resonantly driven dots but is significant for far-detuned driving lasers. Independent interactions also lead to relaxation and pumping of each QD, at rates $\Gamma _{\downarrow }^{j}$ and $\Gamma _{\uparrow }^{j}$ respectively. This is caused by both phonons and photons, with the former being most effective when the dots are resonantly driven. The photon contribution is more nuanced: for red-detuned driving lasers (i.e. $\theta =0$) the photons will induce relaxation, for blue-detuned driving lasers (i.e. $\theta =\pi$) they will induce pumping, while for resonantly driven QDs the photons will induce relaxation and pumping at equal rates. The photon bath will also induce pure dephasing of the QDs at rates $\Gamma _{\phi }^{j}$ that is most effective for resonantly driven dots, and which decays to zero for far-detuned driving lasers.

When terms associated with the common phonon bath are introduced, this leads most straightforwardly to modified renormalization, relaxation and pumping rates for each dot individually. However, there are additional effects embedded in H_X and ${{\mathcal{L}}_{X}}$ that lead to completely different physics. We present these additional effects schematically in figure 1 and discuss them in detail below.

Zoom In Zoom Out Reset image size

The coherent term H_X, corresponding to processes mediated by virtual phonons, contains a renormalization and an elastic interaction between the two QDs. The renormalization part causes a shifting in energy of the eigenstates $|gg\rangle$ and $|ee\rangle$ by ${{\Delta }_{X}}$, while keeping the one excitation eigenstates $|eg\rangle$ and $|ge\rangle$ fixed. This does not affect the steady-state dynamics, so it is not important in our further calculations. By contrast, the induced elastic interaction between the two QDs, with strength t_X, changes the QD coherent dynamics and steady state properties. The one exciton states $|eg\rangle$ and $|ge\rangle$ become coupled elastically such that the single exciton eigenstates of the system are entangled states split by t_X. When the two excitons are resonant the new eigenstates are the usual symmetric and antisymmetric states: $|{{\psi }^{+}}\rangle =(|eg\rangle +|ge\rangle )/\sqrt{2}$ and $|{{\psi }^{-}}\rangle =(|eg\rangle -|ge\rangle )/\sqrt{2}$.

The incoherent term ${{\mathcal{L}}_{X}}$, corresponding to processes mediated by real phonons, leads to relaxation and pumping between the two QD states. When the qubits interact separately with the phonon environment, the phonon field induces relaxation and pumping between the single QD basis states $|g{{\rangle }_{j}}$ and $|e{{\rangle }_{j}}$ as we saw previously. Since these decoherence mechanisms affect single qubit states they will tend to destroy any coherence between the two QDs. However, the interaction with a common bath induces relaxation (accompanied by phonon emission) and pumping (accompanied by phonon absorption) between the two QD entangled states $|{{\psi }^{+}}\rangle$ and $|{{\psi }^{-}}\rangle$ and the two QD states $|gg\rangle$ and $|ee\rangle$ with rates $\Gamma _{\downarrow }^{X}$ (relaxation) and $\Gamma _{\uparrow }^{X}$ (pumping)—see figure 1. If the rate constants for these processes are favourable, then entangled steady states of the two QDs can result.

The induced elastic coupling strength and the inter-QD pumping and relaxation rates are both proportional to ${\rm sin} {{\theta }_{A}}{\rm sin} {{\theta }_{B}}$ and are therefore most effective when the dots are resonantly driven. To get an idea of the strength of these couplings we plot them in figure 2, for the case of 3D phonons. The couplings exhibit an oscillatory behaviour in both plots, due to the oscillatory function ${{F}_{\mathbb{D}}}(x)$. The amplitude of the oscillations decreases with increasing distance between the dots due phonon dispersion.

Zoom In Zoom Out Reset image size

From the master equation we wish to obtain the equations of motion of the expectation values of system operators. We find that the resulting equations decouple into a set that includes all the population dynamics (i.e. the density matrix elements ${{\rho }_{ee-ee}},{{\rho }_{eg-eg}},{{\rho }_{ge-ge}},{{\rho }_{gg-gg}}$, where ${{\rho }_{x-y}}\equiv \langle x|\rho |y\rangle$), together with the one excitation coherences ρ_eg-ge and ρ_ge-eg. All other coherences decouple and decay to zero in the long time limit, and so are not relevant for steady state calculations. A convenient operator basis is therefore formed by the single qubit operators σ_z^j, the two qubit operator $\sigma _{z}^{A}\sigma _{z}^{B}$ and the two coherences, which in a rotating frame have the form $c(t)={{\chi }_{r}}+{\rm i}{{\chi }_{i}}$ where ${{\chi }_{r}}=\Re \{{{{\rm e}}^{{\rm i}(W_{B}^{\prime }-W_{A}^{\prime })t}}\langle ge|\rho |eg\rangle \}$ and ${{\chi }_{i}}=\Im \{{{{\rm e}}^{{\rm i}(W_{B}^{\prime }-W_{A}^{\prime })t}}\langle ge|\rho |eg\rangle \}$. In order to simplify the form of the resulting equations of motion we define

$$\begin{eqnarray}&&{{\gamma }_{j}}=\Gamma _{\downarrow }^{j}+\Gamma _{\uparrow }^{j}\qquad \qquad {\rm (single}\;{\rm QD}\;{\rm population}\;{\rm difference}\;{\rm relaxation}\;{\rm rate),}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\gamma _{d}^{j}=2\Gamma _{\phi }^{j}+{{\gamma }_{j}}/2\qquad \quad {\rm (single}\;{\rm QD}\;{\rm dephasing}\;{\rm rate),}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{{D}_{j}}=(\Gamma _{\uparrow }^{j}-\Gamma _{\downarrow }^{j})/{{\gamma }_{j}}\quad \ \;{\rm (single}\;{\rm QD}\;{\rm relative}\;{\rm steady}\;{\rm state}\;{\rm inversion),}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&D={{D}_{A}}+{{D}_{B}}\qquad \quad \ \;{\rm (total}\;{\rm steady}\;{\rm state}\;{\rm inversion),}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\delta D={{D}_{A}}-{{D}_{B}}\qquad \quad {\rm (inversion}\;{\rm difference}\;{\rm between}\;{\rm QDs),}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\gamma ={{\gamma }_{A}}+{{\gamma }_{B}}\qquad \qquad \ \ {\rm (total}\;{\rm relaxation}\;{\rm rate),}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{\gamma }_{d}}=\gamma _{d}^{A}+\gamma _{d}^{B}\qquad \qquad {\rm (total}\;{\rm dephasing}\;{\rm rate),}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{\Gamma }^{X}}=\Gamma _{\downarrow }^{X}+\Gamma _{\uparrow }^{X}\qquad \quad {\rm (sum}\;{\rm of}\;{\rm cross}\;{\rm relaxation}\;{\rm rates),}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\delta {{\Gamma }^{X}}=\Gamma _{\downarrow }^{X}-\Gamma _{\uparrow }^{X}\qquad \ \,({\rm difference}\;{\rm of}\;{\rm cross}\;{\rm relaxation}\;{\rm rates}).\end{eqnarray} \tag{ 24 }$$

We may then write the following closed set of equations:

$$\begin{eqnarray}&&{\boldsymbol{\dot{x}}} =M{\boldsymbol{x}} +{\boldsymbol{k}} ,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}M=\left( \begin{array}{ccccccccccccccc} -{{\gamma }_{A}} & 0 & 0 & -2\delta {{\Gamma }^{X}} & 2{{t}_{X}} \\ 0 & -{{\gamma }_{B}} & 0 & -2\delta {{\Gamma }^{X}} & -2{{t}_{X}} \\ {{\gamma }_{B}}{{D}_{B}} & {{\gamma }_{A}}{{D}_{A}} & -\gamma & 4{{\Gamma }^{X}} & 0 \\ \delta {{\Gamma }^{X}}/4 & \delta {{\Gamma }^{X}}/4 & {{\Gamma }^{X}}/2 & -{{\gamma }_{d}} & 0 \\ {{t}_{X}}/4 & -{{t}_{X}}/4 & 0 & 0 & -{{\gamma }_{d}} \\ \end{array} \right),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}{{{\boldsymbol{x}} }^{T}}=\left( \begin{array}{ccccccccccccccc} \langle \sigma _{z}^{A}\rangle & \langle \sigma _{z}^{B}\rangle & \langle \sigma _{z}^{A}\sigma _{z}^{B}\rangle & {{\chi }_{r}} & {{\chi }_{i}} \\ \end{array} \right),\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}{{{\boldsymbol{k}} }^{T}}=\left( \begin{array}{ccccccccccccccc} {{\gamma }_{A}}{{D}_{A}} & {{\gamma }_{B}}{{D}_{B}} & 0 & 0 & 0 \\ \end{array} \right).\end{eqnarray} \tag{ 28 }$$

When a steady-state is reached, the above equations of motion can be solved exactly to obtain the steady state vector ${{{\boldsymbol{x}} }_{ss}}$. Although a full exact solution is straightforward, it is cumbersome and offers little insight. Therefore, in this section we will only consider the case when the phonon induced couplings t_X and ${{\Gamma }^{X}}$ are small in comparison to the other decoherence rates in the problem (i.e. ${{t}_{X}},{{\Gamma }_{X}}\ll {{\gamma }_{d}}$). In this pertubative limit we obtain solutions for the steady-state expectation values of the system operators ${{{\boldsymbol{x}} }_{ss}}$:

$$\begin{eqnarray}&&{{\langle \sigma _{z}^{A}\rangle }_{ss}}={{D}_{A}}+\frac{2{{t}_{X}}}{{{\gamma }_{A}}}{{\chi }_{i}}-\frac{2\delta {{\Gamma }^{X}}}{{{\gamma }_{A}}}{{\chi }_{r}}\equiv {{D}_{A}}+{{\delta }_{A}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{{\langle \sigma _{z}^{B}\rangle }_{ss}}={{D}_{B}}-\frac{2{{t}_{X}}}{{{\gamma }_{B}}}{{\chi }_{i}}-\frac{2\delta {{\Gamma }^{X}}}{{{\gamma }_{B}}}{{\chi }_{r}}\equiv {{D}_{B}}+{{\delta }_{B}},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\quad {{\langle \sigma _{z}^{A}\sigma _{z}^{B}\rangle }_{ss}}={{\langle \sigma _{z}^{A}\rangle }_{ss}}{{\langle \sigma _{z}^{B}\rangle }_{ss}}-{{\delta }_{A}}{{\delta }_{B}}+\frac{2}{\gamma }\left( {{t}_{X}}\delta D{{\chi }_{i}}+\delta {{\Gamma }^{X}}D{{\chi }_{r}}+2{{\Gamma }^{X}}{{\chi }_{r}} \right),\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{{\chi }_{r}}=\frac{\delta {{\Gamma }^{X}}}{4{{\gamma }_{d}}}\left( {{\langle \sigma _{z}^{A}\rangle }_{ss}}+{{\langle \sigma _{z}^{B}\rangle }_{ss}} \right)+\frac{{{\Gamma }^{X}}}{2{{\gamma }_{d}}}{{\langle \sigma _{z}^{A}\sigma _{z}^{B}\rangle }_{ss}}\approx \frac{\delta {{\Gamma }^{X}}}{4{{\gamma }_{d}}}D+\frac{{{\Gamma }^{X}}}{2{{\gamma }_{d}}}{{D}_{A}}{{D}_{B}},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{{\chi }_{i}}=\frac{{{t}_{X}}}{4{{\gamma }_{d}}}\left( {{\langle \sigma _{z}^{B}\rangle }_{ss}}-{{\langle \sigma _{z}^{A}\rangle }_{ss}} \right)\approx -\frac{{{t}_{X}}}{4{{\gamma }_{d}}}\delta D.\end{eqnarray} \tag{ 33 }$$

In fact, the expressions obtained for ${{\chi }_{i,r}}$ are correct up to second order in $({{t}_{X}},{{\Gamma }^{X}})/{{\gamma }_{d}}$ while the expressions obtained for $\langle \sigma _{z}^{j}\rangle$ and ${{\langle \sigma _{z}^{A}\sigma _{z}^{B}\rangle }_{ss}}$ are correct to third order in $({{t}_{X}},{{\Gamma }^{X}})/{{\gamma }_{d}}$ since there is no contribution to ${{\chi }_{i,r}}$ of second order in $({{t}_{X}},{{\Gamma }^{X}})/{{\gamma }_{d}}$.

The coherences ${{\chi }_{i,r}}$ depend strongly on the values of the D_j, the qubit inversion in the absence of any QD interaction. To obtain insight into the behaviour of these coherences we consider how the D_j depend on the detuning angle θ_j. At zero temperature we find:

$$\begin{eqnarray}&&{{D}_{j}}({{\theta }_{j}})=\frac{\Gamma _{\uparrow }^{j}-\Gamma _{\downarrow }^{j}}{\Gamma _{\uparrow }^{j}+\Gamma _{\downarrow }^{j}}=\frac{-{\rm cos} {{\theta }_{j}}J_{\gamma }^{j}(\omega _{l}^{j})-\frac{{{{\rm sin} }^{2}}{{\theta }_{j}}}{4}{{J}_{p}}({{W}_{j}})}{(1-\frac{{{{\rm sin} }^{2}}{{\theta }_{j}}}{2})J_{\gamma }^{j}(\omega _{l}^{j})+\frac{{{{\rm sin} }^{2}}{{\theta }_{j}}}{4}{{J}_{p}}({{W}_{j}})}.\end{eqnarray} \tag{ 34 }$$

We can see from this that there are two potentially competing processes at work. First, the phonon coupling always tends to relax the QD to its ground state, regardless of θ. However, the photon field will invert the QD when it is driven by a far blue-detuned laser (i.e. ${{D}_{j}}(\pi )=1$), while relaxing the QD to its ground state when it is driven by a far red-detuned laser (i.e. ${{D}_{j}}(0)=-1$). The strength of the radiation coupling then determines how blue-detuned the driving lasers must be to counteract the effect of the phonon field and invert the QD. We will need to use the appearance of such an inversion to explain some of our key results below.

Up to this point we have been completely general regarding the relationship between the angles ${{\theta }_{A}}$ and ${{\theta }_{B}}$ of the two QDs. However, in order to further investigate the strength of the coherences ${{\chi }_{i,r}}$ (and therefore the effects due to the induced interaction between the QDs) we need to be more specific about the relationship between the two driving lasers. Therefore, for the rest of the paper we will focus on two particular cases:

• Case 1: similarly detuned driving lasers. Firstly we consider the case when ${{\theta }_{A}}={{\theta }_{B}}$ and the two QDs exhibit the same behaviour. In this case D_A = D_B and therefore ${{\chi }_{i}}=0$. This regime is very useful, since then all the terms involving the phonon induced coherent coupling t_X are also zero in equations (29)–(33) and so we can isolate the incoherent QD coupling mediated by real phonons. We should also note that, at zero temperature, ${{\chi }_{r}}=0$ for far red-detuned driving lasers (i.e. ${{\chi }_{r}}(\theta =0)=0$) and therefore we should use blue-detuned driving lasers to ensure a contribution of χ_r.
• Case 2: oppositely detuned driving lasers. Secondly we consider the case ${{\theta }_{A}}=\pi -{{\theta }_{B}}$ which will allow us to investigate the coherent QD coupling mediated by virtual phonons. In this case, for far off resonant driving fields one QD will be in the excited state and the other in the ground state such that ${{D}_{A}}=-{{D}_{B}}$, resulting in a finite χ_i and χ_r.

Since the coupling to the same phonon field results non-zero steady-state coherence between the two qubit states $|eg\rangle$ and $|ge\rangle$ we expect that the two QDs will actually be entangled through their interaction with the phonon bath. For a bipartite system the most common measure of entanglement between the two subsystems is the concurrence, defined as $C(\rho )={\rm max} \left\{ 0,{{\lambda }_{1}}-{{\lambda }_{2}}-{{\lambda }_{3}}-{{\lambda }_{4}} \right\}$, where λ_i are the square roots of the eigenvalues of the matrix $\rho \tilde{\rho }$ listed in decreasing order, where $\tilde{\rho }=({{\sigma }_{y}}\otimes {{\sigma }_{y}}){{\rho }^{*}}\left( {{\sigma }_{y}}\otimes {{\sigma }_{y}} \right)$ [59]. In our case the concurrence has the following expression:

$$\begin{eqnarray}\begin{array}{rcl} C & = & 2\cdot {\rm max} \left\{ 0,{\rm min} \left[ |{{\rho }_{eg-ge}}|,\sqrt{{{\rho }_{eg-eg}}{{\rho }_{ge-ge}}} \right]-\sqrt{{{\rho }_{ee-ee}}{{\rho }_{gg-gg}}} \right\}, \\ C & = & {\rm max} \left\{ 0,{\rm min} \left[ 2\sqrt{\chi _{i}^{2}+\chi _{r}^{2}},\frac{1}{2}\sqrt{{{(1-{{\langle \sigma _{z}^{A}\sigma _{z}^{B}\rangle }_{ss}})}^{2}}-{{({{\langle \sigma _{z}^{A}\rangle }_{ss}}-{{\langle \sigma _{z}^{B}\rangle }_{ss}})}^{2}}} \right] \right. \\ {} & {} & \left. \quad \quad \quad -\frac{1}{2}\sqrt{{{(1+{{\langle \sigma _{z}^{A}\sigma _{z}^{B}\rangle }_{ss}})}^{2}}-{{({{\langle \sigma _{z}^{A}\rangle }_{ss}}+{{\langle \sigma _{z}^{B}\rangle }_{ss}})}^{2}}} \right\}. \\ \end{array}\end{eqnarray} \tag{ 35 }$$

The case of oppositely detuned driving lasers is a lot more efficient in entangling the QDs since it includes contributions from both χ_r and χ_i. We plot concurrence for this case in figure 3, as a function of the driving angle θ for a fixed Rabi frequency W = 0.3 meV (left panel) and as a function of the Rabi frequency strength W for a fixed angle $\theta =0.1$ rad (right panel), both parameter regimes easily accessible experimentally. We can see that the entanglement is large far from resonance because that is where the coherences ${{\chi }_{i}},{{\chi }_{r}}$ are largest. We also notice that both the virtual phonon coupling and the real phonon coupling give rise to entanglement.

Zoom In Zoom Out Reset image size

In obtaining the master equation for the reduced density matrix we lost track of the emitted field by tracing out the photon bath operators ${{a}_{{\boldsymbol{q}} }}$. However, correlation measurements of the emitted photon fields can yield valuable information about the QD properties and the phonon bath [10, 60]. In this section we show how these correlation measurements can be used to find signatures of a common phonon bath and to measure the strength of the induced couplings t_X and ${{\Gamma }^{X}}$.

An experimentally accessible quantity which is very sensitive to QD coupling mechanisms is the (normalized) intensity cross correlation function. For two photon modes a_j and a_k the intensity cross correlation function is defined as:

$$\begin{eqnarray}&&g_{jk}^{(2)}(\tau )=\frac{\langle a_{j}^{\dagger }(0)a_{k}^{\dagger }(\tau ){{a}_{k}}(\tau ){{a}_{j}}(0)\rangle }{\langle a_{j}^{\dagger }(0){{a}_{j}}\rangle \langle a_{k}^{\dagger }(0){{a}_{k}}(0)\rangle },\end{eqnarray} \tag{ 36 }$$

$g_{jk}^{(2)}(\tau )$ is proportional to the probability of detecting a photon in mode k at time $t=\tau$ given that a photon in mode j was detected at time t = 0. If the two photons are completely uncorrelated the intensity correlation function is one; any deviation from one is a signature of correlated photons, and so also a signature of correlations in the QD states which led to the photon emission.

In order to measure $g_{jk}^{(2)}$ experimentally one can use a Hanbury Brown–Twiss interferometer. In this set-up the two photon streams emitted by the QDs are detected separately using a beam splitter and corresponding frequency filters. The photon detectors $\mathcal{A}$ and $\mathcal{B}$ used to detect the photon streams are connected to a timer which is activated when a photon is detected by detector $\mathcal{A}$ and is stopped when another photon is detected by detector $\mathcal{B}$ (negative time correlations $g_{jk}^{(2)}(-\tau )$ can be obtained by delaying the second photon stream by τ). The delay time between the two detection events is recorded and a delay-histogram can be obtained from which $g_{jk}^{(2)}(\tau )$ can be obtained. Since the timing resolution of the detectors is not perfect, the actual measured function will be:

$$\begin{eqnarray}&&\tilde{g}_{jk}^{(2)}(\tau )=\int _{-\infty }^{\infty }{\rm d}t\frac{1}{\sigma \sqrt{2\pi }}{{{\rm e}}^{-\frac{1}{2}{{\left( \frac{t-\tau }{\sigma } \right)}^{2}}}}g_{jk}^{(2)}(t),\end{eqnarray} \tag{ 37 }$$

where σ captures the timing jitter of the photon detectors. Here we will take $\sigma =150$ ps, which is achievable with state-of-the-art photon detectors [61, 62].

In order to relate the correlations in the emitted photon fields to the QD operators and thus to the phonon induced couplings we use the input–output formalism [63]. We start from the initial Hamiltonian equation (9). According to the input–output formalism, the input field driving the QDs and the output field emitted by the QDs are related through the relation ${{a}_{{\rm out}}}(t)={{a}_{{\rm in}}}(t)+\sqrt{{{\Gamma }_{A}}}{{c}_{A}}(t)+\sqrt{{{\Gamma }_{B}}}{{c}_{B}}(t)$, where ${{\Gamma }_{j}}={{J}_{\gamma }}(\omega _{l}^{j})[1+{{\bar{n}}_{\gamma }}(\omega _{l}^{j},T)]$ and for simplicityʼs sake we make the reasonable approximation ${{J}_{\gamma }}(\omega _{l}^{j}\pm {{W}_{j}})={{J}_{\gamma }}(\omega _{l}^{j})$. a_in is given by the classical driving field amplitude corresponding to the two lasers incident on the QDs: in a typical experimental set-up this contribution is eliminated. In terms of the slowly rotating operators ${{\tilde{c}}_{j}}$ the output field is then ${{a}_{{\rm out}}}(t)=\sqrt{{{\Gamma }_{A}}}{\rm exp} ({\rm i}\omega _{l}^{A}t){{\tilde{c}}_{A}}+\sqrt{{{\Gamma }_{B}}}{\rm exp} ({\rm i}\omega _{l}^{B}t){{\tilde{c}}_{B}}$. Thus we can see that the output field is composed of two separate photon fields, one emitted by the first QD with frequencies centered around ω_l^A and one emitted by the second QD with frequencies centered around ω_l^B. Therefore, so long as the dot emission spectra are well resolved, a simple grating can be used to separate the output field into ${{a}_{A}}(t)=\sqrt{{{\Gamma }_{A}}}{{c}_{A}}$ and ${{a}_{B}}(t)=\sqrt{{{\Gamma }_{B}}}{{c}_{B}}$. These two fields are then incident on detectors $\mathcal{A}$ and $\mathcal{B}$ respectively. Therefore, according to input–output formalism:

$$\begin{eqnarray}&&g_{AB}^{(2)}(\tau )=\frac{\langle c_{A}^{\dagger }(0)c_{B}^{\dagger }(\tau ){{c}_{B}}(\tau ){{c}_{A}}(0)\rangle }{\langle c_{A}^{\dagger }(0){{c}_{A}}(0)\rangle \langle c_{B}^{\dagger }(0){{c}_{B}}(0)\rangle }.\end{eqnarray} \tag{ 38 }$$

In order to evaluate this function we need to solve the time dynamics from the equations of motion, equation (25) and then use the quantum regression theorem [64] to relate the correlation functions to the QD operators. The time dynamics solution requires obtaining the eigenvalues of the matrix M and they can only be obtained numerically. Therefore, we are able to obtain $g_{AB}^{(2)}(\tau )$ only numerically. However, $g_{AB}^{(2)}(0)$ can be evaluated analytically since it only requires the evaluation of the two QD operator $c_{A}^{\dagger }{{c}_{A}}c_{B}^{\dagger }{{c}_{B}}$ in the steady state:

$$\begin{eqnarray}\begin{array}{rcl} {{\langle c_{A}^{\dagger }{{c}_{A}}c_{B}^{\dagger }{{c}_{B}}\rangle }_{ss}} & = & \langle XX|{{\rho }_{ss}}|XX\rangle =\frac{{\rm sin} ({{\theta }_{A}}){\rm sin} ({{\theta }_{B}})}{2}{{\chi }_{r}} \\ {} & {} & +\frac{1}{4}\left[ 1+{\rm cos} ({{\theta }_{A}}){{\langle \sigma _{z}^{A}\rangle }_{ss}}+{\rm cos} ({{\theta }_{B}}){{\langle \sigma _{z}^{B}\rangle }_{ss}} \right. \\ {} & {} & +\left. {\rm cos} ({{\theta }_{A}}){\rm cos} ({{\theta }_{B}}){{\langle \sigma _{z}^{A}\sigma _{z}^{B}\rangle }_{ss}} \right], \\ \end{array}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{{\langle c_{A}^{\dagger }{{c}_{A}}\rangle }_{ss}}=\langle XX|{{\rho }_{ss}}|XX\rangle +\langle X0|{{\rho }_{ss}}|X0\rangle =\frac{1}{2}\left[ 1+{\rm cos} ({{\theta }_{A}}){{\langle \sigma _{z}^{A}\rangle }_{ss}} \right].\end{eqnarray} \tag{ 40 }$$

From the above relations we obtain $g_{AB}^{(2)}(0)$:

$$\begin{eqnarray}&&g_{AB}^{(2)}(0)=1+\frac{{\rm cos} ({{\theta }_{A}}){\rm cos} ({{\theta }_{B}})\left[ {{\langle \sigma _{z}^{A}\sigma _{z}^{B}\rangle }_{ss}}-{{\langle \sigma _{z}^{A}\rangle }_{ss}}{{\langle \sigma _{z}^{B}\rangle }_{ss}} \right]+2{\rm sin} ({{\theta }_{A}}){\rm sin} ({{\theta }_{B}}){{\chi }_{r}}}{(1+{\rm cos} ({{\theta }_{A}}){{\langle \sigma _{z}^{A}\rangle }_{ss}})(1+{\rm cos} ({{\theta }_{B}})\langle \sigma _{z}^{B}\rangle )},\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&g_{AB}^{(2)}(0)\approx 1+2\frac{\frac{{\rm cos} ({{\theta }_{A}}){\rm cos} ({{\theta }_{B}})}{\gamma }\left[ {{t}_{X}}\delta D{{\chi }_{i}}+(2{{\Gamma }^{X}}+\delta {{\Gamma }^{X}}D){{\chi }_{r}} \right]+{\rm sin} ({{\theta }_{A}}){\rm sin} ({{\theta }_{B}}){{\chi }_{r}}}{(1+{\rm cos} ({{\theta }_{A}})({{D}_{A}}+{{\delta }_{A}}))(1+{\rm cos} ({{\theta }_{B}})({{D}_{B}}+{{\delta }_{B}}))}.\end{eqnarray} \tag{ 42 }$$

We can see that the incoherent coupling mediated by real phonons has both a first order and a second order (in the small parameters ${{t}_{X}}/\gamma ,\Gamma _{X}^{\uparrow ,\downarrow }/\gamma$) contribution to $g_{AB}^{(2)}(0)$ while the coherent coupling mediated by real phonons has only a second order contribution to $g_{AB}^{(2)}(0)$. However, the first order contribution is most effective at resonance while the second order contribution is most effective far from resonance. Since the denominator of equation (41) decreases at a faster rate than the coupling strengths ${{t}_{X}},{{\Gamma }^{X}}$ we expect that we will see the largest effects far from resonance, even though these will be due to the second order process. It is therefore possible to probe and characterize the common phonon environment by measuring ${{g}^{(2)}}(0)$. In the case of similarly detuned lasers, when there is no virtual phonon contribution, the real phonon coupling strength can be obtained directly from $g_{AB}^{(2)}(0)$. Then, by moving to the case of oppositely detuned lasers and inserting the known value of real phonon coupling, the virtual phonon coupling strength can also be found.

In order to establish whether this simple signature can be observed using current technology, we take into account the finite smearing of the measurement due to timing jitter. We then only have access to the smeared function $\tilde{g}_{AB}^{(2)}(\tau )$, which can only be found numerically. $\tilde{g}_{AB}^{(2)}(\tau )$ is very sensitive to the driving lasers' mixing angles ${{\theta }_{A,B}}$ and renormalized QD frequency $W_{A,B}^{\prime }$, so first, in figure 4, we take values for these parameters where the effect of the common phonon bath is most clearly seen. We show both an example of similarly detuned (right panel) and oppositely detuned lasers (left panel). We drive the first QD with a blue detuned laser (${{\theta }_{A}}=\pi -0.1$) and use Rabi frequencies $W_{A}^{\prime }=W_{B}^{\prime }=0.2$ meV and look at both $g_{AB}^{(2)}(\tau )$ (blue line) and $\tilde{g}_{AB}^{(2)}(\tau )$ (red line). The common phonon environment has a large impact on $\tilde{g}_{AB}^{(2)}(\tau )$, easily detectable experimentally. Indeed, in the case of oppositely detuned lasers the signature of the common phonon bath is extremely large. It would certainly be possible to use such data to fit for t_X and ${{\Gamma }^{X}}$.

Zoom In Zoom Out Reset image size

We next show how these signatures change as several parameters are varied for the case of similarly (oppositely) detuned lasers in figure 5 (figure 6). There are several differences between these two figures. For example, in the case of oppositely detuned lasers the largest signature is when the dot distance is around 10 nm (top left panel), a distance easily accessible experimentally. In contrast, for the case of similarly detuned lasers the signature is largest at zero dot separation, which is obviously unfeasible, but the signature is still large at 10 nm (top left panel). As expected, in the case of oppositely detuned lasers, we can see that the largest signature is at zero temperature (top middle panel) because t_X is independent of temperature while the other decoherence rates grow with increasing temperature. However, because ${{\Gamma }^{X}}$ grows with increasing temperature, in the case of similarly detuned lasers there is an ideal temperature when the signature is largest (top middle panel). We also see that, in the case of oppositely detuned lasers the largest signature is obtained when the QD frequencies $W_{A,B}^{\prime }$ are smallest while in the case of similarly detuned lasers a finite $W_{A,B}^{\prime }$ is optimal (bottom right panel).

There are also similarities between the signatures in figures 5 and 6: both are largest when the two QDs are resonant with each other $W_{A}^{\prime }=W_{B}^{\prime }$ (top right panel) and the largest signature is for off resonant driving lasers (bottom middle panel). We also see that there are ideal radiation field coupling strengths (bottom left panel) at which the signatures are largest.

Using a weak coupling master equation approach, we have shown that two nominally un-coupled driven QDs are nevertheless effectively coupled through the interaction with a common environment. Using a second order Born–Markov master equation we isolated the effects that are solely due to the interaction of the QDs with the same phonon field, and which cannot be explained by considering the QDs coupled to separate phonon fields. We showed that the interaction with the phonon bath results in two types of induced interaction between the two QDs: an elastic interaction mediated by virtual phonons and an inelastic interaction mediated by real phonons. Both types of interaction entangle the two QDs, and we have shown that photon statistics measurements can be used to obtain a signature of this common environment.

We have based our study on a weak coupling master equation. Our neglect of non-Markovian effects here will cause most deviation when the correlation time of the bath is of order the system dynamics timescale—i.e. when the QD frequencies are of the same order as the phonon bath cutoff frequencies ${{\omega }_{e,h,eh}}\approx 1$ meV. Therefore, our weak coupling predictions are quite reliable for the regime in which ${{W}_{A,B}}$ is smaller than this 1 meV scale. In future, we would like to extend this work beyond this regime by either introducing unitary transformations such as the polaron transformation [56], or by using an exact approach such as QUAPI [56, 65, 66].

We would like to thank Jonathan Keeling, Paul Eastham, Helen Cammack, Ahsan Nazir and Jacob Iles-Smith for useful discussions. BWL thanks the Royal Society for a University Research Fellowship, OC would also like to thank the Royal Society for support and the School of Physics and Astronomy in St Andrews for kind hospitality.

In this appendix we obtain a second order Born–Markov master equation for the general case of two systems coupled to the same reservoir. We start from the most general Hamiltonian:

$$\begin{eqnarray}&&H={{H}_{0}}+{{H}_{B}}+{{H}_{I}},\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}&&{{H}_{B}}={{\omega }_{{\boldsymbol{k}} }}c_{{\boldsymbol{k}} }^{\dagger }{{c}_{{\boldsymbol{k}} }},\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&{{H}_{I}}=\mathop{\sum }\limits_{{\boldsymbol{k}} }({{c}_{{\boldsymbol{k}} }}+c_{{\boldsymbol{k}} }^{\dagger })({{t}_{{\boldsymbol{k}} }}{{S}_{A}}+{{g}_{{\boldsymbol{k}} }}{{S}_{B}}+{\rm h}.{\rm c}.),\end{eqnarray} \tag{ A.3 }$$

where H₀ is the Hamiltonian of the two systems which can also contain interactions with other reservoirs and coupling between the two systems, H_B is the Hamiltonian of the reservoir, and H_I is the interaction Hamiltonian between the reservoir and the two systems. S_A and S_B can be any system operators describing the coupling to the reservoir.

A.1. Simple system-reservoir coupling operator

Before deriving a general result, we first consider a simple case when ${{S}_{A}}=a+{{a}^{\dagger }}$ and ${{S}_{B}}=b+{{b}^{\dagger }}$, where $a=|\psi _{i}^{A}\rangle \langle \psi _{j}^{A}|$ and $b=|\psi _{k}^{B}\rangle \langle \psi _{l}^{B}|$, and where $|\psi _{i/j}^{A}\rangle$ and $|\psi _{k/l}^{B}\rangle$ are eigenstates of the system Hamiltonian H₀. In this case we have the relation ${\rm exp} ({\rm i}{{H}_{0}}t)a{\rm exp} (-{\rm i}{{H}_{0}}t)={\rm exp} (-{\rm i}{{\omega }_{A}}t)a$ and ${\rm exp} ({\rm i}{{H}_{0}}t)b{\rm exp} (-{\rm i}{{H}_{0}}t)={\rm exp} (-{\rm i}{{\omega }_{B}}t)b$, with ${{\omega }_{A}}$ and ${{\omega }_{B}}$ the relevant differences in eigenstate energies for systems A and B.

In the Schrödinger picture the Born–Markov second-order master equation takes the form [58]:

$$\begin{eqnarray}&&\dot{\rho }(t)=-{\rm i}[{{H}_{0}},\rho (t)]-\int _{0}^{\infty }{\rm d}\tau {\rm t}{{{\rm r}}_{B}}\left\{ \left[ {{H}_{I}},\left[ {{{\rm e}}^{{\rm i}({{H}_{0}}+{{H}_{B}})\tau }}{{H}_{I}}{{{\rm e}}^{-{\rm i}({{H}_{0}}+{{H}_{B}})\tau }},\rho (t)\otimes {{\rho }_{B}} \right] \right] \right\}.\end{eqnarray} \tag{ A.4 }$$

Assuming a bath in thermal equilibrium we trace out the bath operators to obtain the second order master equation:

$$\begin{eqnarray}\begin{array}{rcl} {\dot{\rho }} & = & -{\rm i}[{{H}_{0}}+\Delta _{1}^{A}[{{a}^{\dagger }},a]-\Delta _{2}^{A}\{{{a}^{\dagger }},a\}+\Delta _{1}^{B}[{{b}^{\dagger }},b]-\Delta _{2}^{B}\{{{b}^{\dagger }},b\},\rho ] \\ {} & {} & +\Gamma _{\downarrow }^{A}\mathcal{D}(a,\rho )+\Gamma _{\uparrow }^{A}\mathcal{D}({{a}^{\dagger }},\rho )+\Gamma _{\downarrow }^{B}\mathcal{D}(b,\rho )+\Gamma _{\uparrow }^{B}\mathcal{D}({{b}^{\dagger }},\rho ) \\ \end{array}\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&+\mathcal{D}(\Gamma _{\downarrow X}^{A}a+\Gamma _{\uparrow X}^{A}{{a}^{\dagger }},b+{{b}^{\dagger }},\rho )+\mathcal{D}(\Gamma _{\downarrow X}^{B}b+\Gamma _{\uparrow X}^{B}{{b}^{\dagger }},a+{{a}^{\dagger }},\rho ),\end{eqnarray} \tag{ A.6 }$$

where we have defined the dissipators $\mathcal{D}$ as:

$$\begin{eqnarray}&&\mathcal{D}(x,y,\rho )\equiv \frac{1}{2}\left[ x\rho {{y}^{\dagger }}+y\rho {{x}^{\dagger }}-{{y}^{\dagger }}x\rho -\rho {{x}^{\dagger }}y \right],\end{eqnarray} \tag{ A.7 }$$

$$\begin{eqnarray}&&\mathcal{D}(x,\rho )\equiv \mathcal{D}(x,x,\rho ),\end{eqnarray} \tag{ A.8 }$$

and where

$$\begin{eqnarray}&&\Delta _{1}^{j}=\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \ \frac{[2\bar{n}({{\omega }_{j}})+1]{{\omega }_{j}}}{({{\omega }^{2}}-\omega _{j}^{2})}\frac{{{J}_{j}}(\omega )}{2\pi },\end{eqnarray} \tag{ A.9 }$$

$$\begin{eqnarray}&&\Delta _{2}^{j}=\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \ \frac{\omega }{({{\omega }^{2}}-\omega _{j}^{2})}\frac{{{J}_{j}}(\omega )}{2\pi },\end{eqnarray} \tag{ A.10 }$$

$$\begin{eqnarray}&&\Gamma _{\downarrow }^{j}={{J}_{j}}({{\omega }_{j}})[\bar{n}({{\omega }_{j}})+1],\end{eqnarray} \tag{ A.11 }$$

$$\begin{eqnarray}&&\Gamma _{\uparrow }^{j}={{J}_{j}}({{\omega }_{j}})\bar{n}({{\omega }_{j}}),\end{eqnarray} \tag{ A.12 }$$

$$\begin{eqnarray}&&\Gamma _{\downarrow X}^{j}=J_{X}^{j}({{\omega }_{j}})[\bar{n}({{\omega }_{j}})+1]+{\rm i}\cdot 2\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \left\{ \frac{\bar{n}(\omega )+1}{\omega -{{\omega }_{j}}}\frac{J_{X}^{j}(\omega )}{2\pi }-\frac{\bar{n}(\omega )}{\omega +{{\omega }_{j}}}\frac{{{[J_{X}^{j}(\omega )]}^{*}}}{2\pi } \right\},\end{eqnarray} \tag{ A.13 }$$

$$\begin{eqnarray}&&\Gamma _{\uparrow X}^{j}={{[J_{X}^{j}({{\omega }_{j}})]}^{*}}\bar{n}({{\omega }_{j}})-{\rm i}\cdot 2\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \left\{ \frac{\bar{n}(\omega )}{\omega -{{\omega }_{j}}}\frac{{{[J_{X}^{j}(\omega )]}^{*}}}{2\pi }-\frac{\bar{n}(\omega )+1}{\omega +{{\omega }_{j}}}\frac{J_{X}^{j}(\omega )}{2\pi } \right\},\end{eqnarray} \tag{ A.14 }$$

where, as in the main text, $j\in \{A,B\}$, $\bar{n}(\omega )$ is the average phonon number of frequency ω at temperature T and we have defined the following spectral density functions:

$$\begin{eqnarray}&&{{J}_{A}}(\omega )=2\pi \mathop{\sum }\limits_{{\boldsymbol{k}} }|{{t}_{{\boldsymbol{k}} }}{{|}^{2}}\delta ({{\omega }_{{\boldsymbol{k}} }}-\omega ),\end{eqnarray} \tag{ A.15 }$$

$$\begin{eqnarray}&&{{J}_{B}}(\omega )=2\pi \mathop{\sum }\limits_{{\boldsymbol{k}} }|{{g}_{{\boldsymbol{k}} }}{{|}^{2}}\delta ({{\omega }_{{\boldsymbol{k}} }}-\omega ),\end{eqnarray} \tag{ A.16 }$$

$$\begin{eqnarray}&&J_{X}^{A}(\omega )=2\pi \mathop{\sum }\limits_{{\boldsymbol{k}} }t_{{\boldsymbol{k}} }^{*}{{g}_{{\boldsymbol{k}} }}\delta ({{\omega }_{{\boldsymbol{k}} }}-\omega ),\end{eqnarray} \tag{ A.17 }$$

$$\begin{eqnarray}&&J_{X}^{B}(\omega )=2\pi \mathop{\sum }\limits_{{\boldsymbol{k}} }{{t}_{{\boldsymbol{k}} }}g_{{\boldsymbol{k}} }^{*}\delta ({{\omega }_{{\boldsymbol{k}} }}-\omega )={{[J_{X}^{A}(\omega )]}^{*}}.\end{eqnarray} \tag{ A.18 }$$

A.2. Most general system-reservoir coupling

In the previous section we have obtained the master equation in the case when the system reservoir coupling operator has a simple form ${{S}_{A}}=a+{{a}^{\dagger }}$ and ${{S}_{B}}=b+{{b}^{\dagger }}$, where a is a transition between two eigenstates of the first QD and b represents a transition between two eigenstates of the second QD. However, in general, the system reservoir coupling operator S_j can be any operator acting on QD j. Since any operator on QD j can be written as a sum of transitions between two eigenstates for that QD, then we can easily generalize the calculations in the previous question to obtain the master equation for a more general system reservoir coupling operator.

The most general system-reservoir coupling operator can be written as ${{S}_{j}}={{\sum }_{q}}{{\sigma }_{qj}}$ where ${{\sigma }_{qj}}=|\psi _{n}^{j}\rangle \langle \psi _{m}^{j}|$ where $|{{\psi }_{n,m}}\rangle$ are eigenstates of the jth system Hamiltonian; q represents any pair of indices $\{n,m\}$. Therefore ${\rm exp} ({\rm i}{{H}_{0}}t){{\sigma }_{qj}}{\rm exp} (-{\rm i}{{H}_{0}}t)={\rm exp} (-{\rm i}{{\omega }_{qj}}t)a$, and here ω_q is the relevant difference in eigenstate energies for the particular pair of indices $\{n,m\}$. Therefore, the most general Hamiltonian for the system-reservoir interaction can be written as:

$$\begin{eqnarray}&&H={{H}_{0}}+{{H}_{B}}+{{H}_{I}},\end{eqnarray} \tag{ A.19 }$$

$$\begin{eqnarray}&&{{H}_{B}}={{\omega }_{{\boldsymbol{k}} }}c_{{\boldsymbol{k}} }^{\dagger }{{c}_{{\boldsymbol{k}} }},\end{eqnarray} \tag{ A.20 }$$

$$\begin{eqnarray}&&{{H}_{I}}=\mathop{\sum }\limits_{{\boldsymbol{k}} }\mathop{\sum }\limits_{qj}({{c}_{{\boldsymbol{k}} }}+c_{{\boldsymbol{k}} }^{\dagger })[t_{{\boldsymbol{k}} }^{j}{{\sigma }_{qj}}+{\rm h}.{\rm c}.].\end{eqnarray} \tag{ A.21 }$$

After tracing out the reservoir in Born–Markov fashion we obtain the second order master equation:

$$\begin{eqnarray}\begin{array}{rcl} {\dot{\rho }} & = & -{\rm i}\left[ {{H}_{0}}+\mathop{\sum }\limits_{qj}\Delta _{1}^{qj}[\sigma _{qj}^{\dagger },{{\sigma }_{qj}}]-\Delta _{2}^{qj}\{\sigma _{qj}^{\dagger },{{\sigma }_{qj}}\},\rho \right] \\ {} & {} & +\mathop{\sum }\limits_{q,j}\left\{ \Gamma _{\downarrow }^{qj}\mathcal{D}({{\sigma }_{qj}},\rho )+\Gamma _{\uparrow }^{qj}\mathcal{D}(\sigma _{qj}^{\dagger },\rho ) \right\}+\mathop{\sum }\limits_{q,p\ne q,j}\mathcal{D}\left( \Gamma _{\downarrow x}^{qj}{{\sigma }_{qj}}+\Gamma _{\uparrow x}^{qj}\sigma _{qj}^{\dagger },{{\sigma }_{pj}}+\sigma _{pj}^{\dagger },\rho \right) \\ {} & {} & +\mathop{\sum }\limits_{q,p}\left\{ \mathcal{D}\left( \Gamma _{\downarrow X}^{qA}{{\sigma }_{qA}}+\Gamma _{\uparrow X}^{qA}\sigma _{qA}^{\dagger },{{\sigma }_{pB}}+\sigma _{pB}^{\dagger },\rho \right)+\mathcal{D}\left( \Gamma _{\downarrow X}^{pB}{{\sigma }_{pB}} \right. \right. \\ {} & {} & +\left. \left. \Gamma _{\uparrow X}^{pB}\sigma _{pB}^{\dagger },{{\sigma }_{qA}}+\sigma _{qA}^{\dagger },\rho \right) \right\}, \\ \end{array}\end{eqnarray} \tag{ A.22 }$$

where

$$\begin{eqnarray}&&\Delta _{1}^{qj}=\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \ \frac{[2\bar{n}({{\omega }_{qj}})+1]{{\omega }_{qj}}}{({{\omega }^{2}}-\omega _{qj}^{2})}\frac{{{J}_{j}}(\omega )}{2\pi },\end{eqnarray} \tag{ A.23 }$$

$$\begin{eqnarray}&&\Delta _{2}^{qj}=\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \ \frac{\omega }{({{\omega }^{2}}-\omega _{qj}^{2})}\frac{{{J}_{j}}(\omega )}{2\pi },\end{eqnarray} \tag{ A.24 }$$

$$\begin{eqnarray}&&\Gamma _{\downarrow }^{qj}={{J}_{j}}({{\omega }_{qj}})[\bar{n}({{\omega }_{qj}})+1],\end{eqnarray} \tag{ A.25 }$$

$$\begin{eqnarray}&&\Gamma _{\uparrow }^{qj}={{J}_{j}}({{\omega }_{qj}})\bar{n}({{\omega }_{qj}}),\end{eqnarray} \tag{ A.26 }$$

$$\begin{eqnarray}&&\Gamma _{\downarrow x}^{qj}={{J}^{j}}({{\omega }_{qj}})[\bar{n}({{\omega }_{qj}})+1]+{\rm i}\cdot 2\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \left\{ \frac{\bar{n}(\omega )+1}{\omega -{{\omega }_{qj}}}\frac{{{J}^{j}}(\omega )}{2\pi }-\frac{\bar{n}(\omega )}{\omega +{{\omega }_{qj}}}\frac{{{J}^{j}}(\omega )}{2\pi } \right\},\end{eqnarray} \tag{ A.27 }$$

$$\begin{eqnarray}&&\Gamma _{\uparrow x}^{qj}={{J}^{j}}({{\omega }_{qj}})\bar{n}({{\omega }_{qj}})-{\rm i}\cdot 2\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \left\{ \frac{\bar{n}(\omega )}{\omega -{{\omega }_{qj}}}\frac{{{J}^{j}}(\omega )}{2\pi }-\frac{\bar{n}(\omega )+1}{\omega +{{\omega }_{qj}}}\frac{{{J}^{j}}(\omega )}{2\pi } \right\},\end{eqnarray} \tag{ A.28 }$$

$$\begin{eqnarray}\begin{array}{rcl} \Gamma _{\downarrow X}^{qj} & = & J_{X}^{j}({{\omega }_{qj}})[\bar{n}({{\omega }_{qj}})+1] \\ {} & {} & +{\rm i}\cdot 2\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \left\{ \frac{\bar{n}(\omega )+1}{\omega -{{\omega }_{qj}}}\frac{J_{X}^{j}(\omega )}{2\pi }-\frac{\bar{n}(\omega )}{\omega +{{\omega }_{qj}}}\frac{{{[J_{X}^{j}(\omega )]}^{*}}}{2\pi } \right\}, \\ \end{array}\end{eqnarray} \tag{ A.29 }$$

$$\begin{eqnarray}\begin{array}{rcl} \Gamma _{\uparrow X}^{qj} & = & {{[J_{X}^{j}({{\omega }_{qj}})]}^{*}}\bar{n}({{\omega }_{qj}}) \\ {} & {} & -{\rm i}\cdot 2\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \left\{ \frac{\bar{n}(\omega )}{\omega -{{\omega }_{qj}}}\frac{{{[J_{X}^{j}(\omega )]}^{*}}}{2\pi }-\frac{\bar{n}(\omega )+1}{\omega +{{\omega }_{qj}}}\frac{J_{X}^{j}(\omega )}{2\pi } \right\}. \\ \end{array}\end{eqnarray} \tag{ A.30 }$$

This is the most general form of the second order master equation for two systems coupled to the same reservoir. (The master equation can be put in Lindblad form but we will only do this for the simpler case of two excitons.)

The above master equation is complete but it also contains terms that contribute very little to the dynamics and a simplified master equation can be obtained by 'secularizing' the master equation as outlined in [67]. Secularization is an approximation akin to a rotating wave approximation (RWA), which allows us to eliminate the incoherent terms in the master equation which oscillate a lot faster compared to the system timescales. Equivalently, we can say that secularization allows us to eliminate the incoherent terms in the master equation which do not conserve energy and therefore are forbidden to second order. We expect that, similarly to a RWA, this approximation holds as long as the strength of the terms ignored is small compared to their oscillation frequency.

A.3. Two un-coupled driven excitons interacting with a common phonon bath

In the case of a pair of un-coupled excitons coupled to the same phonon bath and to separate photon baths as in equations (9) and (11) we have the reservoir-system coupling operator $({{S}_{j}}+S_{j}^{\dagger })=\tilde{c}_{j}^{\dagger }{{\tilde{c}}_{j}}$ which yields ${{S}_{j}}=\frac{{\rm sin} \theta }{2}\sigma _{-}^{j}+\frac{1}{4}(\sigma _{z}^{j}{\rm cos} \theta +I)$. We also have ${{g}_{{\boldsymbol{k}} }}={{t}_{{\boldsymbol{k}} }}{{{\rm e}}^{{\rm i}{\boldsymbol{k}} \cdot {\boldsymbol{d}} }}$ [47].

The starting Hamiltonian is:

$$\begin{eqnarray}&&H=\frac{{{W}_{A}}}{2}\sigma _{z}^{A}+\frac{{{W}_{B}}}{2}\sigma _{z}^{B}+\mathop{\sum }\limits_{{\boldsymbol{k}} }{{\omega }_{{\boldsymbol{k}} }}b_{{\boldsymbol{k}} }^{\dagger }{{b}_{{\boldsymbol{k}} }}+{{H}_{\gamma }}\end{eqnarray} \tag{ A.31 }$$

$$\begin{eqnarray}&&+\mathop{\sum }\limits_{{\boldsymbol{k}} }({{b}_{{\boldsymbol{k}} }}+b_{{\boldsymbol{k}} }^{\dagger })\left[ {{t}_{{\boldsymbol{k}} }}\left( \frac{{\rm sin} \theta }{2}\sigma _{-}^{A}+\frac{1}{4}(\sigma _{z}^{A}{\rm cos} \theta +I) \right)+{\rm h}.{\rm c}. \right]\end{eqnarray} \tag{ A.32 }$$

$$\begin{eqnarray}&&+\mathop{\sum }\limits_{{\boldsymbol{k}} }({{b}_{{\boldsymbol{k}} }}+b_{{\boldsymbol{k}} }^{\dagger })\left[ {{t}_{{\boldsymbol{k}} }}{{{\rm e}}^{{\rm i}{\boldsymbol{k}} \cdot {\boldsymbol{d}} }}\left( \frac{{\rm sin} \theta }{2}\sigma _{-}^{B}+\frac{1}{4}(\sigma _{z}^{B}{\rm cos} \theta +I) \right)+{\rm h}.{\rm c}. \right],\end{eqnarray} \tag{ A.33 }$$

where we define ${{H}_{\gamma }}={{\sum }_{{\boldsymbol{q}} }}{{\Theta }_{{\boldsymbol{q}} }}a_{{\boldsymbol{q}} }^{\dagger }{{a}_{{\boldsymbol{q}} }}+{{\sum }_{j}}({{{\rm e}}^{{\rm i}\omega _{l}^{j}t}}\tilde{c}_{j}^{\dagger }+{{{\rm e}}^{-{\rm i}\omega _{l}^{j}t}}{{\tilde{c}}_{j}}){{\sum }_{{\boldsymbol{q}} }}f_{{\boldsymbol{q}} }^{j}(a_{{\boldsymbol{q}} }^{\dagger }+{{a}_{{\boldsymbol{q}} }})$; this is the Hamiltonian resulting from the coupling of the QDs to the radiation field. Since the two QDs probe the radiation field at different frequencies the photon bath does not mediate any interaction between the two QDs and therefore we can safely treat the photon bath as effectively two photon baths interacting separately and individually with each QD. We assume that the reader is familiar with the usual procedure of tracing out the photon operators ${{a}_{{\boldsymbol{q}} }}$ to obtain the contribution to the second-order Born–Markov master equation and therefore directly insert the resulting dissipators in equation (A.37) below.

We define the following spectral density functions:

$$\begin{eqnarray}&&{{J}_{p}}(\omega )=\mathop{\sum }\limits_{{\boldsymbol{k}} }|{{t}_{{\boldsymbol{k}} }}{{|}^{2}}\delta (\omega -{{\omega }_{{\boldsymbol{k}} }}),\end{eqnarray} \tag{ A.34 }$$

$$\begin{eqnarray}&&J_{\gamma }^{j}(\omega )=\mathop{\sum }\limits_{{\boldsymbol{q}} }|{{f}_{{\boldsymbol{q}} }}{{|}^{2}}\delta (\omega -{{\omega }_{{\boldsymbol{q}} }}),\end{eqnarray} \tag{ A.35 }$$

$$\begin{eqnarray}&&{{J}_{X}}(\omega )\equiv \mathop{\sum }\limits_{{\boldsymbol{k}} }|{{t}_{{\boldsymbol{k}} }}{{|}^{2}}{{{\rm e}}^{{\rm i}{\boldsymbol{k}} \cdot {\boldsymbol{d}} }}\delta (\omega -{{\omega }_{{\boldsymbol{k}} }})={{F}_{\mathbb{D}}}\left( \frac{\omega d}{{{c}_{s}}} \right){{J}_{p}}(\omega ),\end{eqnarray} \tag{ A.36 }$$

where c_s is the speed of sound in GaAs and ${{F}_{\mathbb{D}}}(x)$ is a function that depends on the dimensionality of the phonons such that ${{F}_{\mathbb{D}}}(x)={\rm sinc}(x)$ in 3D, ${{F}_{\mathbb{D}}}(x)={{J}_{0}}(x)$ (J₀ is the Bessel function of the first kind) in two-dimensional and ${{F}_{\mathbb{D}}}(x)={\rm exp} (ix)$ in one-dimensional.

Using the above notation we obtain the following secularized master equation in Lindblad form describing the interaction of two excitons with a common phonon bath but separate photon baths:

$$\begin{eqnarray}\begin{array}{rcl} {\dot{\rho }} & = & -{\rm i}\left[ H,\rho \right] \\ {} & {} & +\mathop{\sum }\limits_{j}\left\{ \left[ \Gamma _{\downarrow }^{j}-|\Re \left( \Gamma _{\downarrow }^{X} \right)|-|\Im \left( \delta \Gamma _{\downarrow }^{X} \right)| \right]\mathcal{D}(\sigma _{-}^{j},\rho ) \right. \\ {} & {} & +\left. \left[ \Gamma _{\uparrow }^{j}-|\Re \left( \Gamma _{\uparrow }^{X} \right)|-|\Im \left( \delta \Gamma _{\uparrow }^{X} \right)| \right]\mathcal{D}(\sigma _{+}^{j},\rho ) \right\} \\ {} & {} & +|\Re \left( \Gamma _{\downarrow }^{X} \right)|\mathcal{D}(\sigma _{-}^{A}\pm \sigma _{-}^{B},\rho )+|\Re \left( \Gamma _{\uparrow }^{X} \right)|\mathcal{D}(\sigma _{+}^{A}\pm \sigma _{+}^{B},\rho ) \\ {} & {} & +|\Im \left( \delta \Gamma _{\downarrow }^{X} \right)|\mathcal{D}({\rm i}\sigma _{-}^{A}\pm \sigma _{-}^{B},\rho )+|\Im \left( \delta \Gamma _{\uparrow }^{X} \right)|\mathcal{D}({\rm i}\sigma _{+}^{A}\pm \sigma _{+}^{B},\rho ), \\ \end{array}\end{eqnarray} \tag{ A.37 }$$

where in the last four terms the sign ± between the two operators forming the dissipators is the sign of their corresponding rates and where

$$\begin{eqnarray*}\begin{array}{rcl} H & = & \frac{{{W}_{A}}+{{\Delta }_{A}}}{2}\sigma _{z}^{A}+\frac{{{W}_{B}}+{{\Delta }_{B}}}{2}\sigma _{z}^{B}+\frac{{{\Delta }_{X}}}{2}\sigma _{z}^{A}\sigma _{z}^{B} \\ {} & {} & +\frac{\Im \left( \Gamma _{\downarrow }^{X}+\Gamma _{\uparrow }^{X} \right)}{2}(\sigma _{+}^{A}\sigma _{-}^{B}+\sigma _{-}^{A}\sigma _{+}^{B})+{\rm i}\frac{\Re \left( \delta \Gamma _{\downarrow }^{X}-\delta \Gamma _{\uparrow }^{X} \right)}{2}(\sigma _{+}^{A}\sigma _{-}^{B}-\sigma _{-}^{A}\sigma _{+}^{B}), \\ {{\Delta }_{j}} & = & \frac{{{{\rm sin} }^{2}}{{\theta }_{j}}}{2}\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \frac{(2\bar{n}(\omega )+1){{W}_{j}}}{{{\omega }^{2}}-W_{j}^{2}}\frac{{{J}_{p}}(\omega )}{2\pi }+{\rm cos} ({{\theta }_{j}})\mathcal{P}\int _{0}^{\infty }{\rm d}\omega \frac{1}{\omega }\frac{{{J}_{p}}(\omega )+\Re \{{{J}_{X}}(\omega )\}}{2\pi }, \\ {{\Delta }_{X}} & = & {\rm cos} {{\theta }_{A}}{\rm cos} {{\theta }_{B}}\int _{0}^{\infty }\frac{1}{\omega }\frac{\Re \{{{J}_{X}}(\omega )\}}{2\pi }, \\ \Gamma _{\downarrow }^{j} & = & {{{\rm cos} }^{4}}\left( \frac{{{\theta }_{j}}}{2} \right)J_{\gamma }^{j}(\omega _{l}^{j})+\frac{{{{\rm sin} }^{2}}({{\theta }_{j}})}{4}{{J}_{p}}({{W}_{j}})[\bar{n}({{W}_{j}})+1], \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{rcl} \Gamma _{\uparrow }^{j} & = & {{{\rm sin} }^{4}}\left( \frac{{{\theta }_{j}}}{2} \right)J_{\gamma }^{j}(\omega _{l}^{j})+\frac{{{{\rm sin} }^{2}}({{\theta }_{j}})}{4}{{J}_{p}}({{W}_{j}})\bar{n}({{W}_{j}}), \\ \Gamma _{\phi }^{j} & = & \frac{{{{\rm sin} }^{2}}{{\theta }_{j}}}{4}J_{\gamma }^{j}(\omega _{l}^{j}), \\ \Gamma _{\downarrow }^{X} & = & \frac{\Gamma _{\downarrow X}^{A}+\Gamma _{\downarrow X}^{B}}{2},\quad \Gamma _{\uparrow }^{X}=\frac{\Gamma _{\uparrow X}^{A}+\Gamma _{\uparrow X}^{B}}{2},\quad \delta \Gamma _{\downarrow }^{X}=\frac{\Gamma _{\downarrow X}^{A}-\Gamma _{\downarrow X}^{B}}{2},\quad \delta \Gamma _{\uparrow }^{X}=\frac{\Gamma _{\uparrow X}^{A}-\Gamma _{\uparrow X}^{B}}{2}, \\ \end{array}\end{eqnarray} \tag{ A.38 }$$

and where we have defined:

$$\begin{eqnarray*}\begin{array}{rcl} \Gamma _{\downarrow X}^{j} & = & \frac{{\rm sin} {{\theta }_{A}}{\rm sin} {{\theta }_{B}}}{4}\left\{ J_{X}^{j}({{W}_{j}})[\bar{n}({{W}_{j}})+1]+{\rm i}\cdot 2\mathcal{P} \right. \\ {} & {} & \times \left. \int _{0}^{\infty }{\rm d}\omega \left[ \frac{\bar{n}(\omega )+1}{\omega -{{W}_{j}}}\frac{J_{X}^{j}(\omega )}{2\pi }-\frac{\bar{n}(\omega )}{\omega +{{W}_{j}}}\frac{{{[J_{X}^{j}(\omega )]}^{*}}}{2\pi } \right] \right\}, \\ \Gamma _{\uparrow X}^{j} & = & \frac{{\rm sin} {{\theta }_{A}}{\rm sin} {{\theta }_{B}}}{4}\left\{ {{[J_{X}^{j}({{W}_{j}})]}^{*}}\bar{n}({{W}_{j}})-{\rm i}\cdot 2\mathcal{P} \right. \\ {} & {} & \times \left. \int _{0}^{\infty }{\rm d}\omega \left[ \frac{\bar{n}(\omega )}{\omega -{{W}_{j}}}\frac{{{[J_{X}^{j}(\omega )]}^{*}}}{2\pi }-\frac{\bar{n}(\omega )+1}{\omega +{{W}_{j}}}\frac{J_{X}^{j}(\omega )}{2\pi } \right] \right\}. \\ \end{array}\end{eqnarray*}$$

Since the most pronounced effects are at resonance when ${{W}_{A}}\approx {{W}_{B}}$ in analytical calculations we ignore the terms proportional to $\delta \Gamma _{\downarrow ,\uparrow }^{X}$ in A.37, although we use the most general master equation in numerical simulations.