## 1 Introduction
If we take the laser frequency of the applied field at exact resonance with the \(|1\rangle\leftrightarrow|4\rangle\) transition, the density matrix equations, under the dipole and rotating wave approximations, take the form:  
\[\dot{\rho}_{11} = \gamma_{41}\rho_{44}+\gamma_{31}\rho_{33}-2\Omega\rho_{41}^{I}, \tag{1}\] \[\dot{\rho}_{22} = \gamma_{32}\rho_{33}+\gamma_{42}\rho_{44}+2\Omega\dot{\rho}_{23}^ {I},\] (2) \[\dot{\rho}_{33} = -(\gamma_{31}+\gamma_{32})\rho_{33}-2\Omega\rho_{23}^{I},\] (3) \[\dot{\rho}_{44} = -(\gamma_{41}+\gamma_{42})\rho_{44}+2\Omega\rho_{41}^{I},\] (4) \[\dot{\rho}_{41}^{I} = -\frac{\gamma_{41}+\gamma_{42}}{2}\rho_{41}^{I}+\Omega(\rho_{11} -\rho_{44}),\] (5) \[\dot{\rho}_{23}^{R} = -\frac{\gamma_{31}+\gamma_{32}}{2}\bar{\rho}_{23}^{R}+(\omega_{21 }+\omega_{43})\bar{\rho}_{23}^{I},\] (6) \[\dot{\rho}_{23}^{I} = -\frac{\gamma_{31}+\gamma_{32}}{2}\bar{\rho}_{23}^{I}-(\omega_{21 }+\omega_{43})\bar{\rho}_{23}^{R}\] (7) \[+\Omega(\rho_{33}-\rho_{22}),\]  
where \(\rho_{nm}=\rho_{mn}^{*}\), and we defined \(\dot{\rho}_{23}=\rho_{23}e^{-i(\omega_{41}+\omega_{43})t}\). Furthermore, with the superscripts \(R,I\), we denote, respectively, the real and imaginary part. The energy of level \(|n\rangle\) is \(\hbar\omega_{n},n=1,2,3,4\), and the time-dependent, in general, Rabi frequency, which determines the interaction of the QD with the applied laser field, is \(\Omega(t)\). As we took the laser field to be at exact resonance with the transition \(|1\rangle\leftrightarrow|4\rangle\), we chose its frequency as \(\omega_{0}=\omega_{4}-\omega_{1}=\omega_{41}\). We note that the dynamics that will be studied below will be from a few hundred ps to some ns. The spin-flip trion relaxation that may connect directly states \(|4\rangle\) and \(|3\rangle\) had a characteristic time of 1 \(\upmu\)s [10], so it did not contribute in the system dynamics in the time scales of interest, and it was omitted. This approximation was done in all relevant studies; see, e.g., [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. The density matrix equations also include the effect of spontaneous decay from the trion to the spin states, where \(\gamma_{nm}\) denotes the spontaneous decay rate for the transition \(|n\rangle\rightarrow|m\rangle\) with:  
\[\gamma_{41}=\gamma_{32}=\gamma_{x}=F_{x}\gamma,\quad\gamma_{42}=\gamma_{31}= \gamma_{z}=F_{z}\gamma, \tag{8}\]  
and \(\gamma\) is the free-space decay rate, for which we use the value \(\hbar\gamma=1.2\)\(\upmu\)eV [10, 21] in the calculations. Furthermore, \(F_{l}\), with \(l=x,z\), is the Purcell factor for the influence of the MoS\({}_{2}\) monolayer on the spontaneous decay rates. The geometry we considered contains a homogeneous host dielectric environment with permittivity \(\varepsilon_{1}\) with a MoS\({}_{2}\) monolayer in it; see Figure 1b. The MoS\({}_{2}\) monolayer lied on the \(x-y\) plane at \(z=0\). The QD was at distance \(z=R\) above the monolayer. The Purcell factors at \(z=R\) were calculated from first-order perturbation theory, with the application of [PERSON]'s golden rule, summing over all final states, and calculation of the relevant electromagnetic Green's tensor at the position of the QD, with [PERSON]'s tensor representing the response of the geometry under consideration to a point-like excitation [34]. For details of the calculation, we refer the reader to [28]:  
\[\dot{\rho}_{11} = \gamma_{41}\rho_{44}+\gamma_{31}\rho_{33}-2\Omega\rho_{41}^{I}, \tag{1}\]0 \[\dot{\rho}_{11} = \gamma_{41}\rho_{44}+\gamma_{31}\rho_{33}-2\Omega\rho_{41}^{I}, \tag{1}\]1  
Here, \(k_{s}\), \(k_{z1}=\sqrt{k_{1}^{2}-k_{s}^{2}}\) denote the in-plane and perpendicular components of the wavevector, respectively, and \(k_{1}=\frac{\omega}{c}\sqrt{\varepsilon_{1}}\) denotes the wavenumber in the host dielectric environment.