How does one get this steady-state solution?

Question

I am studying the concept of a "gain medium" inside a laser cavity. The active atoms are taken to have the two levels $a$ and $b$, separated by energy $\hbar \omega$ and represented by a density matrix $\rho$. The atoms are stationary.

The equation of motion of the density matrix is

$$\dot{\rho} = -i[H, \rho] - \dfrac{1}{2}\left(\Gamma \rho + \rho \Gamma \right) + \lambda,$$

where

$$\rho = \begin{bmatrix} \rho_{aa} & \rho_{ab} \\ \rho_{ba} & \rho_{bb} \end{bmatrix}, \ \ \ \ \ H = \begin{bmatrix} W_{a} & V \\ V & W_{b} \end{bmatrix}, \\ \Gamma = \begin{bmatrix} \gamma_{a} & 0 \\ 0 & \gamma_{b} \end{bmatrix}, \ \ \ \ \ \lambda = \begin{bmatrix} \lambda_{a} & 0 \\ 0 & \lambda_{b} \end{bmatrix}$$

The perturbation Hamiltonian is $\hbar V$, and the unperturbed energies of the levels are $\hbar W_a$ and $\hbar W_b$. Furthermore, the two levels decay with damping constants $\gamma_a$ and $\gamma_b$, and are populated by pumping at rates $\lambda_a$ and $\lambda_b$.

Therefore, using the Fourier expansion of the electric field $E(z, t) = \sum\limits_n A_n(t) u_n(z)$, where $u_n(z) = \sin(k_n z)$ and $k_n = \dfrac{n \pi}{L}$, my notes claim that the electric dipole approximation for the perturbation becomes

$$V(t) = -A(t) \wp u(z) / \hbar$$

Note that we are now considering a single cavity mode, so we drop the subscript $n$.

[See this related question.]

We now write the time dependence of the electric field $A(t)$ in terms of the amplitude $E(t)$ and phase $\phi(t)$, which vary slowly in an optical period $2\pi/\nu$:

$$A(t) = E(t) \cos[\nu t + \phi(t)]$$

[See this related question.]

Now take

$$\delta(t) = E(t) e^{-i[\nu t + \phi(t)]}$$

to be the positive frequency part of the electric field.

[See this related question.]

So we can now write the off-diagonal terms in the density-matrix equation as

$$\dot{\rho}_{ab} = -(i\omega + \gamma)\rho_{ab} - i(\rho_{aa} - \rho_{bb})A(t) \wp u(z)/\hbar,$$

where $\gamma = \dfrac{1}{2} (\gamma_a + \gamma_b)$

[See this related question.]

The steady-state solution of this equation gives

$$\rho_{ab} = -\dfrac{1}{2} \dfrac{i(\rho_{aa} - \rho_{bb}) \delta(t) \wp u(z)}{\hbar[\gamma + i(\omega - \nu)]},$$

where it is assumed that the population inversion varies slowly compared with $\rho_{ab}$, and the "rotating-wave approximation" has been used.

How does one get the steady-state solution

$$\rho_{ab} = -\dfrac{1}{2} \dfrac{i(\rho_{aa} - \rho_{bb}) \delta(t) \wp u(z)}{\hbar[\gamma + i(\omega - \nu)]},$$

and how does the assumption that the population inversion varies slowly compared with $\rho_{ab}$, and the use of the "rotating-wave approximation," play a part in this?

---

EDIT

The furthest I can seem to go with this is as follows:

$$\begin{align} \dot{\rho}_{ab} = \dfrac{-\frac{1}{2} i (\rho_{aa} - \rho_{bb})(-i \nu) e^{-i \nu t} \wp u(z)}{\hbar [\gamma + i(\omega - \nu)]} = \dfrac{-(\rho_{aa} - \rho_{bb}) A(t) \wp u(z)}{\hbar [\gamma + i(\omega - \nu)]}, \end{align}$$

since $\cos(\nu t) = \dfrac{e^{i \nu t} + e^{-i \nu t}}{2} \ \Rightarrow 2\cos(\nu t) - e^{i \nu t} = e^{-i \nu t} = 2 A(t)$, since, for some reason, we seem to be assuming that $\delta(t) = E(t) e^{-i[\nu t + \phi(t)]} \approx e^{-i \nu t}$ and $A(t) = E(t) \cos[\nu t + \phi(t)] \approx \cos[\nu t]$, and, if I'm not mistaken, $e^{i \nu t}$ is the negative frequency part of the electric field, so we discard it (we only want the positive frequency part, as mentioned above).

The confusing thing is the presence of the $\rho_{ab}$ in

$$\dot{\rho}_{ab} = -(i\omega + \gamma)\rho_{ab} - i(\rho_{aa} - \rho_{bb})A(t) \wp u(z)/\hbar$$

I don't actually understand how this could be in $\dot{\rho}_{ab}$, since $\dot{\rho}_{ab}$ itself is obviously differentiated, meaning that $\rho_{ab}$ should have also been differentiated.

I'd appreciate an answer that takes the time to explain all of this, so that I can clearly understand what's going on here.

Answer 1

Here is my solution:

Start from the equation

$$\dot{\rho}_{ab} = -(i\omega + \gamma)\rho_{ab} - i(\rho_{aa} - \rho_{bb})A(t) \wp u(z) / \hbar.$$

The solution to this for $\rho_{ab}$ can be found by first assuming $\rho_{aa}$, $\rho_{bb}$ are constants, as this is the steady state. To keep things simple, lets set $\eta = i(\rho_{aa} - \rho_{bb}) \wp u(z) / \hbar$.

We then have a differential equation of the form

$$\dot{\rho}_{ab} + (i\omega + \gamma)\rho_{ab} = -A(t)\eta.$$

To solve this, you use the method of using an integration factor. Here we simply multiply through by $u(t) = e^{\int (i\omega + \gamma) \ dt}$ to get

$$\dot{\rho}_{ab} u(t) + \dot{u}(t)\rho_{ab} = -A(t) \eta u(t)$$

Notice that the LHS is the derivative $\frac{d(\rho_{ab}u(t))}{dt}$, so we now write

$$\frac{d(\rho_{ab}u(t))}{dt} = -A(t) \eta u(t).$$

We now integrate both sides wrt $t$ and then divide by $u(t)$. This gives

$$\rho_{ab} = \frac{- \eta \int A(t) u(t) \ dt}{u(t)}.$$

If we calculate $A(t) u(t)$, we find

$$A(t)u(t) = \frac{E(t)}{2}\left(\exp{[i(\omega + \nu)t + \gamma t + \phi(t) + c]} + \exp{[i(\omega - \nu)t + \gamma t -\phi(t) + c]}\right)$$

This is where we use the rotating wave approximation: $(\omega+\nu)$ will be a rapidly oscillating term, and we can say that this averages out over all reasonable timescales. This then gives the expression for $\rho_{ab}$:

$$\rho_{ab} = -\frac{\eta\int E(t)\exp{[i(\omega - \nu)t + \gamma t - \phi(t) + c]} \ dt}{2\exp{[i\omega t + \gamma t + c}]}$$

The next step is to do integration by parts on the numerator:

$$\int E(t)\exp{[i(\omega - \nu)t + \gamma t + \phi(t) + c]} \ dt = \frac{1}{{(i(\omega - \nu) + \gamma - \dot{\phi}(t))}}\left( E(t)\exp{[i(\omega - \nu)t + \gamma t - \phi(t) + c]} - \int\dot{E}(t)\exp{[i(\omega - \nu)t + \gamma t - \phi(t) + c]} \ dt \right)$$

Then we take advantage of the fact that $E(t)$ and $\phi(t)$ are slowly varying, so we set their derivatives $\dot{E}(t) = \dot{\phi}(t) = 0$. We can then write our expression for $\rho_{ab}$ as

$$\rho_{ab} = -\frac{\eta E(t)\exp{[i(\omega - \nu)t + \gamma t - \phi(t) + c]}}{2(i(\omega - \nu) + \gamma)\exp{[i\omega t + \gamma t + c]}} = -\frac{\eta E(t)\exp{[-i\nu t - \phi(t)]}}{2(i(\omega - \nu) + \gamma)}$$

Where $\delta(t) = E(t)\exp{[-i\nu t-\phi(t)]}$ and inserting $\eta$ back in:

$$\rho_{ab} = -\frac{i(\rho_{aa} - \rho_{bb})\wp u(z)\delta (t)}{2\hbar (i(\omega - \nu) + \gamma)}.$$

Note: You asked why $\dot{\rho}_{ab}$ depended on $\rho_{ab}$. This just means that the rate of change depends on it's value. Think of the rate of photons emitted from an upper state depends on the population of that upper state.

Answer 2

Apparently, it is assumed that the time dependence of $\rho_{ab}$ is mostly $\sim e^{-i\nu t}$, so $\dot{\rho_{ab}}\approx (-i\nu)\rho_{ab}$. One gets the factor $\frac{1}{2}$ from the expression for cosine as half sum of exponential functions. Only one exponential function being left in the expression is, I believe, the rotating wave approximation.