The Quantum Regression Formula

The quantum regression formula (QRF) found by Lax (1963), (1967) provides a method to compute any two-time correlator from a master equation of the form $d\rho/dt=\mathcal{L}\rho$ [Eq. (2.71)] (for system interacting with Markovian reservoirs). As demonstrated in the book by Carmichael (2002), once one has found a set of operators $C_{\{\eta\}}$ that satisfy

$\displaystyle \Tr(C_{\{\eta\}}\mathcal{L}\Omega_1)=\sum_{\{\lambda\}}M_{\{\eta\lambda\}}\Tr(C_{\{\lambda\}}\Omega_1)$

(2.98)

for the general operator $\Omega_1$ , and the corresponding matrix elements $M_{\{\alpha\beta\}}$ , then, the equations of motion for the two-time correlators follow as (for $\tau\geq 0$ ):

$\displaystyle \frac{d}{d\tau}\langle\Omega_1(t)C_{\{\eta\}}(t+\tau)\rangle =\su... ...ambda\}}M_{\{\eta\lambda\}}\langle\Omega_1(t)C_{\{\lambda\}}(t+\tau)\rangle \,.$

(2.99)

The Hilbert space of correlators is separated in manifolds, just as the Hilbert space of states is separated in manifolds of excitation. The order of a manifold is the minimum number of particles that should be in the system (regardless of the regression matrix) so that the correlator is nonzero. Equivalently, it is the minimum manifold of excitations that should be probed in the dynamics. We will refer to the two-time and one-time correlator manifolds as $\mathcal{N}_k$ and $\tilde{\mathcal{N}}_k$ , respectively.

The first step is to find the set of operators that are needed to compute the correlator of interest. For example, in the case of bosons, in order to set the equations to obtain $\langle\ud{a}(t)a(t+\tau)\rangle$ , we consider $\Omega_1=\ud{a}$ . If this is the only field involved in the dynamics, the most general set of operators in normal order can be written as $\{C_{\{m,n\}}=\ud{a}^ma^n\}$ . For the simple problem of a thermal bath, it is enough to consider only . The only matrix element is $M_{1,1}=-i\omega_a-\frac{\gamma_a-P_a}{2}$ and the correlator can be trivially integrated taking as the initial state the SS value $n_a^\mathrm{SS}$ : $\langle\ud{a}(t)a(t+\tau)\rangle =e^{-i\omega_a \tau}e^{-(\gamma_a-P_a)/2 \tau}n_a^\mathrm{SS}$ . The spectra is again a Lorentzian

$\displaystyle S(\omega)=\frac{1}{\pi}\frac{\frac{\Gamma_a}2}{\big(\frac{\Gamma_a}2\big)^2+(\omega-\omega_a)^2}\,,$

(2.100)

this time broadened by a renormalized bosonic decay rate:

$\displaystyle \Gamma_a=\gamma_a-P_a=\frac{\gamma_a}{1+n_a^{\mathrm{SS}}}\,.$

(2.101)

In this case, in order to have a physical SS (finite populations and correlations which decay with time), we know that the parameters $\gamma_a$ and

correspond to a given thermal bath. This implies that $n_a^{\mathrm{SS}}=\bar n_T$ and $\Gamma_a=\kappa_a$ . On the other hand, if we deal with a free 2LS and we want to compute $\langle\ud{\sigma}(t)\sigma(t+\tau)\rangle$ , the equivalent procedure yields an effective broadening:

$\displaystyle \Gamma_\sigma=\gamma_\sigma+P_\sigma=\frac{\gamma_\sigma}{1-n_\sigma^{\mathrm{SS}}}$

(2.102)

which is different from the decay rate in vacuum $\kappa_\sigma$ . We can also distinguish the nature of the particles in the appearance of the bosonic (

) or the fermionic ( $1-n_\sigma$ ) factors in the dependence of the effective decays with the average occupation. The narrowing (broadening) of the linewidth with pump, as the number of particles increases, is a bosonic (fermionic) effect.

In more complicated systems, the correlator of interest will depend on other correlators, giving rise to a set of coupled equations of the form of Eq. (2.99). The initial values of these equations at $\tau=0$ must also be found, either in the SS ( $t\rightarrow \infty$ ) or for a general time in the SE case. The equations for $\langle\Omega_1(t)C_{\{\eta\}}(t)\rangle$ can be equally found applying the same QRF with $\tilde{\Omega}_1=1$ and a new set of operators $C_{\{\tilde\eta\}}$ where the previous operator $\Omega_1C_{\{\eta\}}$ is included:

$\displaystyle \frac{d}{d t}\langle C_{\{\tilde\eta\}}(t)\rangle = \sum_{\{\tild... ...da\}}M_{\{\tilde\eta\tilde\lambda\}}\langle C_{\{\tilde\lambda\}}(t)\rangle \,.$

(2.103)

In the SE case, the initial state of these differential equations is the initial value of the correlators $\langle C_{\{\tilde\eta\}}(0)\rangle$ . In the SS case, the equations are set to zero and solved as a linear system to find $\langle C_{\{\tilde\eta\}}\rangle^{\mathrm{SS}}$ . The one-time mean values can also be obtained from the density matrix resulting from the master equation by applying the general relation $\langle C_{\{\tilde\eta\}}(t)\rangle =\Tr{[C_{\{\tilde\eta\}}\rho(t)]}$ .

In the general problem of two coupled modes, , , we refer with the label $\{\eta\}=(m,n,\mu,\nu)$ to the two-time correlator $\langle\Omega_1(t) C_{\{\eta\}}(t+\tau)\rangle$ with $C_{\{\eta\}}=\ud{a}^m a^n \ud{b}^\mu b^\nu$ , regardless $\Omega_1$ . This is the most general form for the correlators, grouped in manifolds $\mathcal{N}_k$ . The emission of particles () corresponds to $\Omega_1=\ud{a}$ ( $=\ud{b}$ ). Each two-time correlator will have as initial condition ( $\tau=0$ ) the one-time correlator $\langle C_{\{m+1,n,\mu,\nu\}}\rangle$ ( $\langle C_{\{m,n,\mu+1,\nu\}}\rangle$ ) that belongs to the corresponding manifold of the same order $\tilde{\mathcal{N}}_k$ . The QRF for them, with $\tilde\Omega_1=1$ , applies in a new set of operators $\{\tilde\eta\}=(m+1,n,\mu,\nu)$ ( $\{\tilde\eta\}=(m,n,\mu+1,\nu)$ ). The final result for the correlator of interest, $\langle C_\mathrm{i}(t,\tau)\rangle$ , will be, as we will see in the following Chapters, always of the form

$\displaystyle \langle C_\mathrm{i}(t,\tau)\rangle =\sum_p (l_p(t)+i k_p(t)) e^{-i\omega_p\tau}e^{-\frac{\gamma_p}{2}\tau}\,,$

(2.104)

where all the new parameters, the weights

, the frequencies $\omega_p$ and the effective decay rates $\gamma_p$ , are real. The values

and

depend on the one-time mean-values $\langle C_{\{\tilde\eta\}}(t)\rangle$ , either in the SS or in time

, depending on the case. In the SS, they are already constants and can be normalized $L_p^{\mathrm{SS}}=l_p^{\mathrm{SS}}/\langle C_\mathrm{i}\rangle ^{\mathrm{SS}}$ , $K_p^{\mathrm{SS}}=k_p^{\mathrm{SS}}/\langle C_\mathrm{i}\rangle ^{\mathrm{SS}}$ . In the SE, they have to be integrated and normalized accordingly: $L_p^{\mathrm{SE}}=\int_0^\infty l_p^{\mathrm{SE}}(t) dt/\int_0^\infty\langle C_\mathrm{i}(t,0)\rangle dt$ and $K_p^{\mathrm{SE}}=\int_0^\infty k_p^{\mathrm{SE}}(t) dt/\int_0^\infty\langle C_\mathrm{i}(t,0)\rangle dt$ . From here, Eq. (2.94) leads to the spectrum,

$\displaystyle S_\mathrm{i}(\omega)=\frac1{\pi}\sum_{p}\left[L_p\frac{\frac{\gam... ...{\omega-\omega_p}{\big(\frac{\gamma_p}{2}\big)^2+(\omega-\omega_p)^2}\right]\,,$

(2.105)

which is in general a sum of many peaks labelled with

, with two type of contributions. The first one is the already introduced Lorentzian shape, which is the only contribution for free fields. The second term is the dispersive part that breaks the symmetry of the Lorentzian around the frequency $\omega_p$ . This is due to the interference between nearby resonances. They appear whenever the particles that are emitted are not exactly the eigenstates of the system, due to decoherence or due to the fact that the detectors measure in the bare state basis (photons). They disappear in weak coupling, where bare modes are predominant (giving rise to other manifestations of interferences), but also in the limit of very strong coupling, where polaritons rule the dynamics completely.

This way of computing the spectra, in which we take the Fourier transform explicitly, gives us the structure of the lines in a transparent way. $\omega_p$ and $\gamma_p$ are the line positions and broadenings. They originate from the energy levels structure and uncertainties, whose skeleton is the Hamiltonian eigenstates, but that can be greatly distorted by decoherence. As such, they are independent of the channel of detection (cavity or direct exciton emission) and independent of time. Coefficients and depend on the one time correlators and, therefore, they are different in the SS or the SE cases. They determine which lines actually appear in the spectra, and with which intensity depending on the channel of emission and the quantum state in the system.

Most of the authors, like Savage (1989), Clemens et al. (2004), Porras & Tejedor (2003) or Perea et al. (2004), compute the spectrum with completely numerical methods from the density matrix and master equation. Their results are blind to the underlying individual lines and, therefore, miss all the information on the manifold structure that the spectra contains. This is a very dramatic loss if one is interested in quantum features or the crossover from quantum to classical regime, like in the case of this thesis. However, the lack of this information is not so important when the system is essentially classical or in the classical regime, where there is no quantized manifold structure. We should/will prefer this kind of ``blind'' methods then. In this direction and taking advantage of the SS properties, it is possible to go from the density matrix to the spectra without any need to compute the correlator. This method is described by Mølmer (1996) in his notes:

First, we choose a basis of states $\{\ket{i}\}$ , ordered in a given way and labelled with the index $i=1,\dots , n_\mathrm{max}$ . Then, we obtain the density matrix, in its $n_\mathrm{max}\times n_\mathrm{max}$ matrix form $\rho_{ij}=\bra{i}\rho\ket{j}$ , in the SS. That is, we solve the master equation

$\displaystyle \frac{d}{dt}\rho_{ij}=\sum_{k,l} \mathcal{M}_{ij,kl} \rho_{kl}=0$

(2.106)

where $\mathcal{M}_{ij,kl}$ is a supermatrix $n_\mathrm{max}^2\times n_\mathrm{max}^2$ . The spectral function can be defined in terms of the two operators

and

which constitute the two-time correlator of interest, $C_\mathrm{i}(t,\tau)=A(t)B(t+\tau)$ , as

$\displaystyle s_{AB}(\omega)=\frac{1}{\pi}\Re\int_{0}^\infty \langle A(t)B(t+\tau)\rangle e^{i\omega\tau} d\tau \,.$

(2.107)

By means of the QRF, the spectrum can be put in terms of the matrix form of the operators $A_{ij}$ ( $=\bra{i}A\ket{j}$ ) and $B_{ij}$ and the SS density matrix:

$\displaystyle s_{AB}(\omega)=-\frac{1}{\pi}\Re\sum_{i,j,k,l,m}[(\mathcal{M}+i\omega I)^{-1}]_{ij,kl}\rho_{km}^{\mathrm{SS}} A_{ml}B_{ji}\,.$

(2.108)

is the $n_\mathrm{max}^2\times n_\mathrm{max}^2$ identity matrix. This formula is valid as long as $\langle A\rangle =0$ in the SS. It is useful only when the Hilbert space is not very large, as it requires the inversion of a $n_\mathrm{max}^2\times n_\mathrm{max}^2$ matrix for each point $\omega$ of the spectra.