Obtaining equations of motion for a three-level atom and plasmon system

Question

I am interested in a scenario where a three level atom (cascade/ladder type) is interacting with a spherical metal nanoparticle (MNP) with single plasmonic excitations. Let us assume that the ground, first excited and second excited states of the atom are given by $|S1\rangle$, $|S2\rangle$ and $|S3\rangle$ forming a complete basis for the atom's Hilbert space, and the creation and annihilation operators of the MNP plasmonic mode is given by $\hat{a}^\dagger = |N\rangle \langle0|$ and $\hat{a}=|0\rangle \langle N|$. An incoherent pumping field ($P_{13}$) is applied between levels $|S1\rangle$ and $|S3\rangle$ of the atom. A coherent probe field $E_{pr}=E_{0pr}(e^{i\omega_{pr}t}+e^{-i\omega_{pr}t})$ is applied between the levels $|S2\rangle$ and $|S3\rangle$. The levels $|S1\rangle$ and $|S2\rangle$ undergo dipole-dipole coupling with the plasmonic field. The interaction picture Hamiltonian of the system was found to be $$ H_{sys}=\hbar\Delta_1\hat{a}^\dagger\hat{a} + \hbar\Delta_1|S1\rangle \langle S1| +\hbar\Delta_2|S3\rangle \langle S3|+i\hbar g(\hat{\sigma}_{12}\hat{a}^\dagger - \hat{\sigma}_{12}^\dagger\hat{a})-E_{0pr}\mu_{23}(\hat{\sigma}^{\dagger}_{23}+\hat{\sigma}_{23}) $$ where $\Delta_1=\omega_{plasmon}-\omega_{(\langle S2|-\langle S1|)}$ and $\Delta_2=\omega_{(\langle S3|-\langle S2|)}-\omega_{pr}$, $g$ is a coupling constant, $\mu_{23}$ is the dipole moment between $|S2\rangle$,$|S3\rangle$, $\hat{\sigma}$ and $\hat{\sigma}^{\dagger}$ represent the atomic lowering and raising operators between the corresponding levels.

The full quantum dynamics of the coupled nano system were derived using the master equation for the interaction picture density operator as:

$$ \frac{\partial\hat{\rho}}{\partial t}=\frac{i}{\hbar}\left[\hat{\rho},H_{sys}\right]+\mathcal{L}_m+\mathcal{L}_{12}+\mathcal{L}_{23}+\mathcal{L}_p $$

where $\mathcal{L}_m, \mathcal{L}_{12},\mathcal{L}_{23}$ are the Lindblad terms for the decay of plasmon mode and atomic levels respectively and $\mathcal{L}_p$ is the Lindblad term for the incoherent pump.

$\textbf{Question:}$ I intend to find the equations of motion for the expectation values of the plasmon and atomic annihilation operators ($\langle \hat{a} \rangle$, $\langle \hat{\sigma}_{12} \rangle$, $\langle \hat{\sigma}_{23} \rangle$ ). With the atom being a three level system and the plasmon operators being two level, is it possible to use the trace as follows for this purpose?

Eg:$\frac{\partial \langle \hat{a} \rangle}{\partial t} = Tr\left[ \langle \hat{a} \rangle \frac{\partial \langle \hat{\rho} \rangle}{\partial t}\right] $, $\frac{\partial \langle \hat{\sigma}_{12} \rangle}{\partial t} = Tr\left[ \langle \hat{\sigma}_{12} \rangle \frac{\partial \langle \hat{\rho} \rangle}{\partial t}\right] $

(I intend to use the cyclic property of trace, bosonic commutator relations and orthogonality of atomic states to simplify the above.) If this doesn't seem like a valid approach, it would be much helpful if you can kindly mention the approach I should take or guide me to a similar calculation available online. Thank you

Answer 1

Some general remarks:

1. The equations of motion for the expectation values $\langle \hat{a} \rangle$, $\langle \hat{\sigma}_{12} \rangle$, $\langle \hat{\sigma}_{23} \rangle$ only have closed form at linear order in the driving fields. For higher orders the populations (e.g. $\langle \hat{\sigma}_{33} \rangle$) can change (see Why do populations only change in second order of the driving field?).
2. The equations of motions $$ \frac{\partial\hat{\rho}}{\partial t}=\frac{i}{\hbar}\left[\hat{\rho},H_{sys}\right]+\mathcal{L}_m+\mathcal{L}_{12}+\mathcal{L}_{23}+\mathcal{L}_p $$ can be expressed in a variety of matrix representations, i.e. different bases for the superspace of the density matrix. In your case it would of course be most convenient to choose a basis containing $\langle \hat{a} \rangle$, $\langle \hat{\sigma}_{12} \rangle$, $\langle \hat{\sigma}_{23} \rangle$. For a nice numerical approach see this technique.
3. I don't really understand these trace relations. I don't see how they are relevant for deriving equations of motions though. The equation above is one already, too get equations for expectation values you just have to choose a nice basis and compute a matrix form of the equation $$\frac{\partial\vec{\rho}}{\partial t} = \mathbf{K}\vec{\rho}$$ where $\vec{\rho}$ is now a vector representation of the density matrix (also called a superspace representation, see the link under 2). The step between the density matrix equations of motions and the superspace equations of motions is really just computing operator expectations with your basis states and bringing them in matrix form, so I don't see how these traces are relevant. Unless you want to do some approximations of course, but then the question as posed does not seem self-contained.