Finding four-wave-mixing terms in photon-atom interaction (by adding for fields to EIT)

Question

The Problem

I'm interested in solving for the four-wave-mixing terms in Electromagnetically Induced Transparency (EIT) when additional fields are added (while maintaining 3 atomic levels).

Normally when working with an N-level system we use (N-1) fields, (EIT has 3 atomic levels and 2 fields). Adding additional fields are known to cause "wave-mixing" terms in which light of one frequency is converted to another frequency.

I am interested in the following level scheme as shown in the diagram:

The level scheme here is identical to normal EIT, except that we have an additional EIT system that is detuned above the normal EIT.

We expect that when one of the weak fields ($\Omega_s, \Omega_p$) is off, light of the first field is created in the frequency of the second field. To see this visually we expect that when $\Omega_p, \Omega_c,\Omega_d $ are on, light is created in the frequency $\omega_s$, and we get four-wave-mixing (FWM) as indicated by the figure below:

I expect to see terms of the form $E_{out} \propto e^{i \omega_p t}\Omega_s \Omega_c \Omega_d$ and $E_{out} \propto e^{i \omega_s t}\Omega_p \Omega_c \Omega_d$

And a similar (four-level, instead of three-level, system) solution has been solved here. In their four-level system they have terms $\rho_{31}= \frac{\Omega_p |\Omega_d|^2-\Omega_s \Omega_c \Omega_d^*}{D}$ and $\rho_{41}= \frac{\Omega_s |\Omega_c|^2-\Omega_s \Omega_c^* \Omega_d}{D}$.

My Attempt

I can get an analytical solution for $p_{12}$ and $p_{13}$, but I do not see the expected four-wave-mixing terms.I've written below a summary of how I obtained this (in hopes someone can help me figure out what I'm doing wrong).

The Hamiltonian I start out with is:

\begin{equation*} H_{EIT} = \\ \begin{bmatrix} \omega_a & 0 & \mu_{13}\\ 0 & \omega_b & \mu_{23}\\ \mu_{13}^* & \mu_{23}^* & \omega_c \\ \end{bmatrix} \left( E_{p} Cos(\omega_p t+\phi_p)+ E_{s} Cos(\omega_s t+\phi_s) + E_{c}Cos(\omega_c+\phi_c)+E_{d}Cos(\omega_d+\phi_d)\right) \end{equation*}

Applying the RWA-approximation seems pretty straightforward. It doesn't seem as though having these additional fields break anything about the RWA, and we can remove the counterpropagating terms as usual. $E_s$ and $E_p$ are similar so the same rules for $\omega_p$ apply for $\omega_s$ and likewise for $E_c$ and $E_d$ : $(\omega_s-\omega_{ac})<<(\omega_c-\omega_{ac})$ )

Removing the counterpropagating terms, this can be reduced to

\begin{equation*} \left( \begin{array}{ccc} \text{$\omega_1$}-\text{$\omega_3 $} & 0 & e^{i t \text{$\Delta_p $}} \text{$\Omega_p $}+e^{i t (\text{$\Delta_p $}-\text{$\Delta_1 $})} \text{$\Omega_s $} \\ 0 & \text{$\omega_2 $}-\text{$\omega_3 $} & e^{i t \text{$\Delta $c}} \text{$\Omega $c}+e^{i t (\text{$\Delta_c $}-\text{$\Delta_2 $})} \text{$\Omega_2 $} \\ e^{-i t \text{$\Delta_p $}} \text{$\Omega_p $}+e^{-i t (\text{$\Delta_p $}-\text{$\Delta_1 $})} \text{$\Omega_s $} & e^{-i t \text{$\Delta_c $}} \text{$\Omega_c $}+e^{-i t (\text{$\Delta_c $}-\text{$\Delta_2 $})} \text{$\Omega_d $} & 0 \\ \end{array} \right) \end{equation*} Or in the original reference frame:

\begin{equation*} \left( \begin{array}{ccc} \text{$\omega_1$}-\text{$\omega_3 $} & 0 & e^{i t \text{$\omega_p $}} \text{$\Omega_p $}+e^{i t (\text{$\omega_p $}-\text{$\Delta_1 $})} \text{$\Omega_s $} \\ 0 & \text{$\omega_2 $}-\text{$\omega_3 $} & e^{i t \text{$\omega $c}} \text{$\Omega $c}+e^{i t (\text{$\omega_c $}-\text{$\Delta_2 $})} \text{$\Omega_2 $} \\ e^{-i t \text{$\omega_p $}} \text{$\Omega_p $}+e^{-i t (\text{$\omega_p $}-\text{$\Delta $})} \text{$\Omega_s $} & e^{-i t \text{$\omega_c $}} \text{$\Omega_c $}+e^{-i t (\text{$\omega_c $}-\text{$\Delta_2 $})} \text{$\Omega_d $} & 0 \\ \end{array} \right) \end{equation*}

For EIT, we look rotate our system in the corotating frame, look at the Von-Neumann equation with some additional decay terms ($\frac{dP}{dt} = [H,P]-LP$), solve the system at steady-state ($\frac{dP}{dt} \implies 0$) using the adiabatic approximation ($\rho_{11} \approx 1, \rho_{22} \approx \rho_{33} \approx 0$). Then when having a solution for the coherence terms we see that our output $\propto \rho_{12}$ and we have our answer.

For this system, it's not obvious to me is there exists a corotating frame that get's rid of the time dependencies of the equation.

Now since my goal is to have an analytic solution for the behavior of the polarization at steady-state (so that I can find what kind of output light this interaction produces), my thought is that maybe I can solve for the density matrix elements without rotating out their time dependence. Additionally, I've decided to stay in the corotating frame, simply to put the diagonal in terms of the detuning:

Using the standard corotating rotation: \begin{equation} U = \left( \begin{array}{ccc} e^{-i t \text{$\omega $p}} & 0 & 0 \\ 0 & e^{-i t \text{$\omega $c}} & 0 \\ 0 & 0 & 1 \\ \end{array} \right) \end{equation}

We can write writing our effective hamiltonian $(H_{eff} = i\frac{dU}{dt}.U^\dagger + U.H.U^\dagger$):

\begin{equation} H_{eff} = \left( \begin{array}{ccc} 0 & 0 & \text{$\Omega $p}+e^{-i t (\text{$\omega $p}-\text{$\omega $s})} \text{$\Omega $s} \\ 0 & \text{$\Delta $p} & \text{$\Omega $c} \\ \text{$\Omega_p^*$}+e^{i t (\text{$\omega $p}-\text{$\omega $s})} \text{$\Omega_s^*$} & \text{$\Omega^*_c$} & \text{$\Delta $c} \\ \end{array} \right) \end{equation}

Note that my Hamiltonian still has time-dependence. Now I'm not sure if this step is correct, but I'm not going to change anything about the standard procedure for solving for the coherence terms. That is, I'm going to simply say that $\frac{dP}{dt}\rightarrow 0$ and solve for (which seems a little bit crazy...should steadystate not be time dependent?) Proceeding with these naive steps:

\begin{equation*} \frac{dP}{dt} = [H(t), P] + \mathcal{L} P\\ \frac{dP}{dt} \rightarrow 0 \\ [H(t), P] + \mathcal{L} P\rightarrow 0 \\ \end{equation*}

We get a system of equations of the form M.P = b. The matrix M is:

\begin{equation} \left( \begin{array}{cccccc} 0 & \text{$\Omega $c}+e^{i t \text{$\Delta $2}} \text{$\Omega $d} & 0 & 0 & 0 & -\text{$\Omega $p}-e^{-i t \text{$\Delta $1}} \text{$\Omega $s} \\ \text{$\Omega $c}+e^{-i t \text{$\Delta $2}} \text{$\Omega $d} & -\Delta & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\text{$\Omega $p}-e^{i t \text{$\Delta $1}} \text{$\Omega $s} & \text{$\Omega $c}+e^{-i t \text{$\Delta $2}} \text{$\Omega $d} & 0 \\ 0 & 0 & \text{$\Omega $p}+e^{-i t \text{$\Delta $1}} \text{$\Omega $s} & \frac{\Gamma }{2}-\Delta & 0 & 0 \\ 0 & 0 & -\text{$\Omega $c}-e^{i t \text{$\Delta $2}} \text{$\Omega $d} & 0 & \frac{\Gamma }{2}+\Delta & 0 \\ -\text{$\Omega $p}-e^{i t \text{$\Delta $1}} \text{$\Omega $s} & 0 & 0 & 0 & 0 & \frac{\Gamma }{2}+\Delta \\ \end{array} \right)\end{equation}

and P and b are

\begin{equation} P = \left( \begin{array}{c} \text{$\rho $12} \\ \text{$\rho $13} \\ \text{$\rho $21} \\ \text{$\rho $23} \\ \text{$\rho $31} \\ \text{$\rho $32} \\ \end{array} \right), b = \left( \begin{array}{c} 0 \\ \text{$\Omega $p}+\text{$\Omega $s} e^{-i \text{$\Delta $1} t} \\ 0 \\ 0 \\ -\text{$\Omega $p}-\text{$\Omega $s} e^{i \text{$\Delta $1} t} \\ 0 \\ \end{array} \right) \end{equation}

Now I have a solution for my coherence terms $\rho_{13}(t), \rho_{12}(t)$:

\begin{align*} \rho_{13}(t) &= \frac{2 e^{-2 i \text{$\Delta $1} t} \left(\text{$\Omega $s}+\text{$\Omega $p} e^{i \text{$\Delta $1} t}\right)^2 \left(\text{$\Omega $p}+\text{$\Omega $s} e^{i \text{$\Delta $1} t}\right)}{D_1} \\ \rho_{23}(t) &= \frac{2 \left(\text{$\Omega $s}+\text{$\Omega $p} e^{i \text{$\Delta $1} t}\right) \left(\text{$\Omega $p}+\text{$\Omega $s} e^{i \text{$\Delta $1} t}\right) \left(\text{$\Omega $d}+\text{$\Omega $c} e^{i \text{$\Delta $2} t}\right)}{D_2} \\ \end{align*} $D_1 = \Gamma \left(\text{$\Omega $c}^2+\text{$\Omega $d}^2\right)+2 \Delta \left(\text{$\Omega $c}^2+\text{$\Omega $d}^2-\text{$\Omega $p}^2-\text{$\Omega $s}^2\right)+2 \text{$\Omega $c} \text{$\Omega $d} (\Gamma +2 \Delta ) \cos (\text{$\Delta $2} t)-4 \Delta \text{$\Omega $p} \text{$\Omega $s} \cos (\text{$\Delta $1} t)$ $D_2 = e^{i t (\text{$\Delta $1}+\text{$\Delta $2})} \left(\Gamma \left(\text{$\Omega $c}^2+\text{$\Omega $d}^2+\text{$\Omega $p}^2+\text{$\Omega $s}^2\right)+2 \Delta \left(-\text{$\Omega $c}^2-\text{$\Omega $d}^2+\text{$\Omega $p}^2+\text{$\Omega $s}^2\right)\right)+\text{$\Omega $c} \text{$\Omega $d} (\Gamma -2 \Delta ) e^{i t (\text{$\Delta $1}+2 \text{$\Delta $2})}+\text{$\Omega $p} \text{$\Omega $s} (\Gamma +2 \Delta ) e^{i t (2 \text{$\Delta $1}+\text{$\Delta $2})}+\text{$\Omega $c} \text{$\Omega $d} (\Gamma -2 \Delta ) e^{i \text{$\Delta $1} t}+\text{$\Omega $p} \text{$\Omega $s} (\Gamma +2 \Delta ) e^{i \text{$\Delta $2} t}$

Now we can use these terms to find FWM terms that I'm interested in. Using the definition of polarization: \begin{equation*} \textbf{P} = N Tr\left( \begin{bmatrix} \rho_{aa} & \rho_{ab} & \rho_{ac}\\ \rho_{ba} & \rho_{bb} & \rho_{bc}\\ \rho_{ca} & \rho_{cb} & \rho_{cc} \\ \end{bmatrix} \begin{bmatrix} 0 & 0 & \mu_{13}\\ 0 & 0 & \mu_{23}\\ \mu_{13} & \mu_{23} & 0 \\ \end{bmatrix}\right) \end{equation*}

We can see our output polarizability is: $P_{out} = \mu_{31} \rho_{13} + \mu_{32} \rho_{23} + c.c.$

Now when I extract out the terms that are proportional to $e^{i \omega_p t}$ and $e^{i \omega_st}$, I don't see ANY terms that are linear in the probe and control field (after doing a taylor expansion). This is very strange to me as there are a few sources that this should act like a frequency beamsplitter, meaning that I should have terms proportional to $\Omega_s \Omega_c^* \Omega_d$ and $\Omega_p \Omega_c^* \Omega_d$.

Maybe someone can point out what I'm doing wrong? Maybe I'm messing up the time-dependent part of the solution? Maybe there is a different framework for working with four-wave mixing?

Answer 1

I hope it's okay if instead of trying to find errors in your approach, I just describe how I'd solve the problem in a similar way. By the way, I think you're actually going in the right direction but the occasional typos and inconsistencies made it a bit hard for me to single out the problematic bits.

First of all, as pointed out in one of the comments, there seem to be some errors in the original Hamiltonian. I believe what you wanted to write is something like

$$ H(t) = \left(\begin{array}{ccc} \omega_1 & 0 & 0\\ 0 & \omega_2 & 0\\ 0 & 0 & \omega_3\\ \end{array}\right) + \left(\begin{array}{ccc} 0 & 0 & \mu_{13}\\ 0 & 0 & \mu_{23}\\ \mu_{13} & \mu_{23} & 0\\ \end{array}\right) E(t), $$ where $$ E(t) = E_p \cos(\omega_p t + \phi_p) + E_s \cos(\omega_s t + \phi_s) + E_c \cos(\omega_c t + \phi_c) + E_d \cos(\omega_d t + \phi_d). $$ Here I'm assuming that $\mu_{13}$ and $\mu_{23}$ are real which can be done without loss of generality.

I think it also makes sense to simplify the problem as much as possible before trying to understand it. So instead of the general case where the conversion between frequencies $\omega_p$ and $\omega_s$ can happen in both directions, I will look at the simpler case where we apply only the $\omega_p$ pump and look for output at frequency $\omega_s$ (i.e. $E_s=0$ in the Hamiltonian). The opposite conversion can be analyzed in a completely analogous way.

Before we do any complicated math, here is my intuitive understanding of how wave mixing would happen here: in the absence of $E_s$, we can look at this as almost the standard EIT problem with pumps $E_p$ and $E_c$ except now there is an additional signal with amplitude $E_d$ which is slightly detuned (by $\Delta_1$) from $\omega_c$. If this detuning is small, we can look at the combination of the $E_c$ and $E_d$ drives as something with frequency $\omega_c$ but with a slowly varying amplitude and phase. The slow variation modulates the transmission coefficient $T$ of the EIT system for the weak probe tone $E_p$. So in principle you should expect in the output field a term $T(t)E_p\cos(\omega_p t+\phi_p)$, where $T(t)$ is periodic with frequency $\Delta_1$. By multiplying terms oscillating at $\omega_p$ and $\Delta_1$, you of course get some component at $\omega_p+\Delta_1 = \omega_s$. And voila, that's your frequency mixing process.

Now to actually derive this. First we do the frame transformation $$ U(t) = \left(\begin{array}{ccc} \mathrm{e}^{\mathrm{i}\omega_1 t} & 0 & 0\\ 0 & \mathrm{e}^{\mathrm{i}\omega_2 t} & 0\\ 0 & 0 & \mathrm{e}^{\mathrm{i}(\omega_1+\omega_p) t}\\ \end{array}\right) $$

This way, the effective Hamiltonian is $$ H_{\mathrm{eff}}(t) = \left(\begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & \omega_3-\omega_1-\omega_p\\ \end{array}\right) + \left(\begin{array}{ccc} 0 & 0 & \mu_{13}\mathrm{e}^{-\mathrm{i}\omega_p t}\\ 0 & 0 & \mu_{23}\mathrm{e}^{\mathrm{i}(\omega_2-\omega_1-\omega_p) t}\\ \mu_{13}\mathrm{e}^{\mathrm{i}\omega_p t} & \mu_{23}\mathrm{e}^{-\mathrm{i}(\omega_2-\omega_1-\omega_p) t} & 0\\ \end{array}\right) E(t), $$

We know that $\omega_3-\omega_1-\omega_p=-\Delta_2$ and $\omega_2-\omega_1-\omega_p=-\omega_c$. The only terms which will have slowly varying components that will survive the rotating wave approximation are those proportional to $\mu_{13}E_p$, $\mu_{23}E_c$ and $\mu_{23}E_d$ (we assume the detuning between $\omega_c$ and $\omega_d$ is not large enough to drop the $\mu_{23}E_d$ term). After the approximation, we get \begin{align} H_{\mathrm{eff}}(t) =& \left(\begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & -\Delta_2\\ \end{array}\right) + \\ &\frac{1}{2} \left(\begin{array}{ccc} 0 & 0 & \mu_{13}E_p\mathrm{e}^{\mathrm{i}\phi_p}\\ 0 & 0 & \mu_{23}(E_c\mathrm{e}^{\mathrm{i}\phi_c}+E_d\mathrm{e}^{\mathrm{i}(\Delta_1 t+\phi_d)})\\ \mu_{13}E_p\mathrm{e}^{-\mathrm{i}\phi_p} & \mu_{23}(E_c\mathrm{e}^{-\mathrm{i}\phi_c}+E_d\mathrm{e}^{-\mathrm{i}(\Delta_1 t+\phi_d)}) & 0\\ \end{array}\right). \end{align}

For simplicity, we will denote the coupling terms by $\Omega_p$ and $\Omega_c(t)$, which gives us $$ H_{\mathrm{eff}}(t) = \left(\begin{array}{ccc} 0 & 0 & \frac{1}{2}\Omega_p\\ 0 & 0 & \frac{1}{2}\Omega_c(t)\\ \frac{1}{2}\Omega_p^{*} & \frac{1}{2}\Omega_c^{*}(t) & -\Delta_2\\ \end{array}\right). $$

Now I believe an important step is to assume that $\Delta_1$ is much slower than the timescale on which the system approaches equilibrium. This way, you can assume the density matrix $\rho(t)$ of the system at each time $t$ is approximately the steady state corresponding to the Hamiltonian $H_{\mathrm{eff}}(t)$. Just like you, I'll make the assumption $\rho_{11}\approx 1$, $\rho_{22}\approx\rho_{33}\approx 0$ and write the system of steady-state equations for $\rho_{12}$, $\rho_{13}$ and $\rho_{32}$ (the other three off-diagonal element are always complex conjugates of the first three, so we might as well save ourselves the trouble and not worry about them for now): \begin{align} \mathrm{i}\rho_{13}\Omega_{c}^{*}/2 - \gamma\rho_{12}/2 - \mathrm{i}\rho_{32}\Omega_{p}/2 &= 0\\ \mathrm{i}\rho_{12}\Omega_{c}/2 - (\mathrm{i}\Delta_2+\Gamma/2) \rho_{13} &= -\mathrm{i}\Omega_{p}/2\\ -\mathrm{i}\rho_{12}\Omega_{p}^{*}/2 + (\mathrm{i}\Delta_2-\Gamma/2) \rho_{32} &= 0\\ \end{align}

Let's pause now and note that these equations are quite similar to the ones you get, except for the term with $\gamma$ in my first equation which describes loss of coherence between states 1 and 2 and which you don't have. Why did I add this? It turns out that without this decoherence mechanism, the solution of these equations is somewhat ill-behaved. For instance, if all the drives are turned off ($\Omega_p = \Omega_c = 0$), $\rho_{12}$ is indeterminate - it completely drops out of the equations. Moreover, look what happens if you fix the ratio $\eta:=\Omega_p/\Omega_c$ and let $\Omega_c\to 0$. You can show that $\rho_{12}$ has a finite limit which depends on $\eta$. This means that the density matrix elements as functions of the two variables $\Omega_p$ and $\Omega_c$ show a rather nasty behavior around $\Omega_p=\Omega_c=0$. Most importantly, they can't be expanded into a Taylor series in $\Omega_p$ and $\Omega_c$. But this is exactly what you'd like to do in your quest to find the wave mixing terms.

Another way to look at this is that to identify the wave mixing process, you'd like to treat the pumps as a perturbation. But the generator of the free dynamics of the system is in a certain sense degenerate (it has multiple steady states), and therefore perturbation theory fails.

At any rate, the decoherence term I added should fix this as it ensures that the free non-driven system has only a single steady state with vanishing off-diagonal density matrix elements. Solving the equations then gives us \begin{align} \rho_{13} &= -\frac{\mathrm{i}\Omega_p}{2\mathrm{i}\Delta_2+\Gamma} -\frac{\mathrm{i}\Omega_p|\Omega_c|^2}{\gamma(2\mathrm{i}\Delta_2+\Gamma)^2} + \ldots\\ \rho_{12} &= \frac{\Omega_p \Omega_c^{*}}{\gamma(2\mathrm{i}\Delta_2+\Gamma)} + \ldots, \end{align} where we are skipping terms of fourth and higher order in the drive strengths.

Now we recall that to transform back into the lab frame, we need to undo $U(t)$. This means multiplying $\rho_{12}$ by $\exp(\mathrm{i}(\omega_2-\omega_1)t)$ and $\rho_{13}$ by $\exp(\mathrm{i}\omega_p t)$. We also need to express $\Omega_p$ and $\Omega_c(t)$ in terms of $E_p$, $E_c$ and $E_d$. A simple calculation then tells us that $\rho_{12}$ consists of terms which oscillate at frequencies $\omega_2-\omega_1$ and $\omega_2-\omega_1+\Delta_1$. None of these is equal to the desired $\omega_s$. The first-order term in the matrix element $\rho_{13}$ oscillates at $\omega_p$ and is not interesting. The second term on the other hand has a component at the frequency $\Delta_1+\omega_p = \omega_s$. This term represents the mixing process that you are after. Explicitly, we can write it as $$ \rho_{13} = -\frac{ \mathrm{i}\mu_{13}\mu_{23}^2 E_p E_c E_d\mathrm{e}^{\mathrm{i}(\omega_s t+\phi_p+\phi_d-\phi_c)} }{\gamma(2\mathrm{i}\Delta_2+\Gamma)^2} + \ldots $$

Disclaimer: I am pretty sure minus signs, factors of 2 and other pesky things got lost in the process due to my carelessness. But I believe the gist of the whole business should still hold. I hope it will be useful.