Categories: Electromagnetism, Physics, Quantum mechanics.

Maxwell-Bloch equations

For an electron in a two-level system with time-independent states \(\ket{g}\) (ground) and \(\ket{e}\) (excited), consider the following general solution to the full Schrödinger equation:

\[\begin{aligned} \ket{\Psi} &= c_g \: \ket{g} \exp\!(-i E_g t / \hbar) + c_e \: \ket{e} \exp\!(-i E_e t / \hbar) \end{aligned}\]

Perturbing this system with an electromagnetic wave introduces a time-dependent sinusoidal term \(\hat{H}_1\) to the Hamiltonian. In the electric dipole approximation, \(\hat{H}_1\) is given by:

\[\begin{aligned} \hat{H}_1(t) = - \hat{\vb{p}} \cdot \vb{E}(t) \qquad \quad \hat{\vb{p}} \equiv q \hat{\vb{x}} \qquad \quad \vb{E}(t) = \vb{E}_0 \cos\!(\omega t) \end{aligned}\]

Where \(\vb{E}\) is an electric field, and \(\hat{\vb{p}}\) is the dipole moment operator. From Rabi oscillation, we know that the time-varying coefficients \(c_g\) and \(c_e\) can then be described by:

\[\begin{aligned} \dv{c_g}{t} &= i \frac{q \matrixel{g}{\hat{\vb{x}}}{e} \cdot \vb{E}_0}{2 \hbar} \exp\!\big( i \omega t \!-\! i \omega_0 t \big) \: c_e \\ \dv{c_e}{t} &= i \frac{q \matrixel{e}{\hat{\vb{x}}}{g} \cdot \vb{E}_0}{2 \hbar} \exp\!\big(\!-\! i \omega t \!+\! i \omega_0 t \big) \: c_g \end{aligned}\]

We want to rearrange these equations a bit. Therefore, we split the electric field \(\vb{E}\) like so, where the amplitudes \(\vb{E}_0^{-}\) and \(\vb{E}_0^{+}\) may be slowly varying:

\[\begin{aligned} \vb{E}(t) = \vb{E}^{-}(t) + \vb{E}^{+}(t) = \vb{E}_0^{-} \exp\!(i \omega t) + \vb{E}_0^{+} \exp\!(-i \omega t) \end{aligned}\]

Since \(\vb{E}\) is real, \(\vb{E}_0^{+} = (\vb{E}_0^{-})^*\). Similarly, we define the transition dipole moment \(\vb{p}_0^{-}\):

\[\begin{aligned} \vb{p}_0^{-} \equiv q \matrixel{e}{\vb{x}}{g} \qquad \quad \vb{p}_0^{+} \equiv (\vb{p}_0^{-})^* = q \matrixel{g}{\vb{x}}{e} \end{aligned}\]

With these, the equations for \(c_g\) and \(c_e\) can be rewritten as shown below. Note that \(\vb{E}^{-}\) and \(\vb{E}^{+}\) include the driving plane wave, and the rotating wave approximation is still made:

\[\begin{aligned} \dv{c_g}{t} &= \frac{i}{\hbar} \vb{p}_0^{+} \cdot \vb{E}^{-} \exp\!(- i \omega_0 t) \: c_e \\ \dv{c_e}{t} &= \frac{i}{\hbar} \vb{p}_0^{-} \cdot \vb{E}^{+} \exp\!(i \omega_0 t) \: c_g \end{aligned}\]

Optical Bloch equations

For \(\ket{\Psi}\) as defined above, the corresponding pure density operator \(\hat{\rho}\) is as follows:

\[\begin{aligned} \hat{\rho} = \ket{\Psi} \bra{\Psi} = \begin{bmatrix} c_e c_e^* & c_e c_g^* \exp\!(-i \omega_0 t) \\ c_g c_e^* \exp\!(i \omega_0 t) & c_g c_g^* \end{bmatrix} \equiv \begin{bmatrix} \rho_{ee} & \rho_{eg} \\ \rho_{ge} & \rho_{gg} \end{bmatrix} \end{aligned}\]

Where \(\omega_0 \equiv (E_e \!-\! E_g) / \hbar\) is the resonance frequency. We take the \(t\)-derivative of the matrix elements, and insert the equations for \(c_g\) and \(c_e\):

\[\begin{aligned} \dv{\rho_{gg}}{t} &= \dv{c_g}{t} c_g^* + c_g \dv{c_g^*}{t} \\ &= \frac{i}{\hbar} \vb{p}_0^{+} \cdot \vb{E}^{-} \exp\!(- i \omega_0 t) \: c_e c_g^* - \frac{i}{\hbar} \vb{p}_0^{-} \cdot \vb{E}^{+} \exp\!(i \omega_0 t) \: c_g c_e^* \\ \dv{\rho_{ee}}{t} &= \dv{c_e}{t} c_e^* + c_e \dv{c_e^*}{t} \\ &= \frac{i}{\hbar} \vb{p}_0^{-} \cdot \vb{E}^{+} \exp\!(i \omega_0 t) \: c_g c_e^* - \frac{i}{\hbar} \vb{p}_0^{+} \cdot \vb{E}^{-} \exp\!(- i \omega_0 t) \: c_e c_g^* \\ \dv{\rho_{ge}}{t} &= \dv{c_g}{t} c_e^* \exp\!(i \omega_0 t) + c_g \dv{c_e^*}{t} \exp\!(i \omega_0 t) + i \omega_0 c_g c_e^* \exp\!(i \omega_0 t) \\ &= \frac{i}{\hbar} \vb{p}_0^{+} \cdot \vb{E}^{-} \: c_e c_e^* - \frac{i}{\hbar} \vb{p}_0^{+} \cdot \vb{E}^{-} \: c_g c_g^* + i \omega_0 c_g c_e^* \exp\!(i \omega_0 t) \\ \dv{\rho_{eg}}{t} &= \dv{c_e}{t} c_g^* \exp\!(-i \omega_0 t) + c_e \dv{c_g^*}{t} \exp\!(-i \omega_0 t) - i \omega_0 c_e c_g^* \exp\!(- i \omega_0 t) \\ &= \frac{i}{\hbar} \vb{p}_0^{-} \cdot \vb{E}^{+} \: c_g c_g^* - \frac{i}{\hbar} \vb{p}_0^{-} \cdot \vb{E}^{+} \: c_e c_e^* - i \omega_0 c_e c_g^* \: \exp\!(- i \omega_0 t) \end{aligned}\]

Recognizing the density matrix elements allows us to reduce these equations to:

\[\begin{aligned} \dv{\rho_{gg}}{t} &= \frac{i}{\hbar} \Big( \vb{p}_0^{+} \cdot \vb{E}^{-} \rho_{eg} - \vb{p}_0^{-} \cdot \vb{E}^{+} \rho_{ge} \Big) \\ \dv{\rho_{ee}}{t} &= \frac{i}{\hbar} \Big( \vb{p}_0^{-} \cdot \vb{E}^{+} \rho_{ge} - \vb{p}_0^{+} \cdot \vb{E}^{-} \rho_{eg} \Big) \\ \dv{\rho_{ge}}{t} &= i \omega_0 \rho_{ge} + \frac{i}{\hbar} \vb{p}_0^{+} \cdot \vb{E}^{-} \big( \rho_{ee} - \rho_{gg} \big) \\ \dv{\rho_{eg}}{t} &= - i \omega_0 \rho_{eg} + \frac{i}{\hbar} \vb{p}_0^{-} \cdot \vb{E}^{+} \big( \rho_{gg} - \rho_{ee} \big) \end{aligned}\]

These equations are correct if nothing else is affecting \(\hat{\rho}\). But in practice, these quantities decay due to various processes, e.g. spontaneous emission (see Einstein coefficients).

Let \(\rho_{ee}\) decays with rate \(\gamma_e\). Since the total probability \(\rho_{ee} + \rho_{gg} = 1\), we thus have:

\[\begin{aligned} \Big( \dv{\rho_{ee}}{t} \Big)_{e} = - \gamma_e \rho_{ee} \quad \implies \quad \Big( \dv{\rho_{gg}}{t} \Big)_{e} = \gamma_e \rho_{ee} \end{aligned}\]

Meanwhile, for whatever reason, let \(\rho_{gg}\) decay into \(\rho_{ee}\) with rate \(\gamma_g\):

\[\begin{aligned} \Big( \dv{\rho_{gg}}{t} \Big)_{g} = - \gamma_g \rho_{gg} \quad \implies \quad \Big( \dv{\rho_{gg}}{t} \Big)_{g} = \gamma_g \rho_{gg} \end{aligned}\]

And finally, let the diagonal (perpendicular) matrix elements both decay with rate \(\gamma_\perp\):

\[\begin{aligned} \Big( \dv{\rho_{eg}}{t} \Big)_{\perp} = - \gamma_\perp \rho_{eg} \qquad \quad \Big( \dv{\rho_{ge}}{t} \Big)_{\perp} = - \gamma_\perp \rho_{ge} \end{aligned}\]

Putting everything together, we arrive at the optical Bloch equations governing \(\hat{\rho}\):

\[\begin{aligned} \boxed{ \begin{aligned} \dv{\rho_{gg}}{t} &= \gamma_e \rho_{ee} - \gamma_g \rho_{gg} + \frac{i}{\hbar} \Big( \vb{p}_0^{+} \cdot \vb{E}^{-} \rho_{eg} - \vb{p}_0^{-} \cdot \vb{E}^{+} \rho_{ge} \Big) \\ \dv{\rho_{ee}}{t} &= \gamma_g \rho_{gg} - \gamma_e \rho_{ee} + \frac{i}{\hbar} \Big( \vb{p}_0^{-} \cdot \vb{E}^{+} \rho_{ge} - \vb{p}_0^{+} \cdot \vb{E}^{-} \rho_{eg} \Big) \\ \dv{\rho_{ge}}{t} &= - \Big( \gamma_\perp - i \omega_0 \Big) \rho_{ge} + \frac{i}{\hbar} \vb{p}_0^{+} \cdot \vb{E}^{-} \Big( \rho_{ee} - \rho_{gg} \Big) \\ \dv{\rho_{eg}}{t} &= - \Big( \gamma_\perp + i \omega_0 \Big) \rho_{eg} + \frac{i}{\hbar} \vb{p}_0^{-} \cdot \vb{E}^{+} \Big( \rho_{gg} - \rho_{ee} \Big) \end{aligned} } \end{aligned}\]

Many authors simplify these equations a bit by choosing \(\gamma_g = 0\) and \(\gamma_\perp = \gamma_e / 2\).

Including Maxwell’s equations

This two-level system has a dipole moment \(\vb{p}\) as follows, where we use Laporte’s selection rule to remove diagonal terms, by assuming that the electron’s orbitals are odd or even:

\[\begin{aligned} \vb{p} &= \matrixel{\Psi}{\hat{\vb{p}}}{\Psi} \\ &= q \Big( c_g c_g^* \matrixel{g}{\hat{\vb{x}}}{g} + c_e c_e^* \matrixel{e}{\hat{\vb{x}}}{e} + c_g c_e^* \matrixel{e}{\hat{\vb{x}}}{g} \exp\!(i \omega_0 t) + c_e c_g^* \matrixel{g}{\hat{\vb{x}}}{e} \exp\!(-i \omega_0 t) \Big) \\ &= q \Big( \rho_{ge} \matrixel{e}{\hat{\vb{x}}}{g} + \rho_{eg} \matrixel{g}{\hat{\vb{x}}}{e} \Big) = \vb{p}_0^{-} \rho_{ge}(t) + \vb{p}_0^{+} \rho_{eg}(t) \equiv \vb{p}^{-}(t) + \vb{p}^{+}(t) \end{aligned}\]

Where we have split \(\vb{p}\) analogously to \(\vb{E}\) by defining \(\vb{p}^{+} \equiv \vb{p}_0^{+} \rho_{eg}\). Its equation of motion can then be found from the optical Bloch equations:

\[\begin{aligned} \dv{\vb{p}^{+}}{t} = \vb{p}_0^{+} \dv{\rho_{eg}}{t} = - \vb{p}_0^{+} \Big( \gamma_\perp + i \omega_0 \Big) \rho_{eg} + \frac{i}{\hbar} \vb{p}_0^{+} \Big( \vb{p}_0^{-} \cdot \vb{E}^{+} \Big) \Big( \rho_{gg} - \rho_{ee} \Big) \end{aligned}\]

Some authors do not bother multiplying \(\rho_{ge}\) by \(\vb{p}_0^{+}\). In any case, we arrive at:

\[\begin{aligned} \boxed{ \dv{\vb{p}^{+}}{t} = - \Big( \gamma_\perp + i \omega_0 \Big) \vb{p}^{+} - \frac{i}{\hbar} \Big( \vb{p}_0^{-} \cdot \vb{E}^{+} \Big) \vb{p}_0^{+} d } \end{aligned}\]

Where we have defined the population inversion \(d \in [-1, 1]\) as follows, which quantifies the electron’s excitedness:

\[\begin{aligned} d \equiv \rho_{ee} - \rho_{gg} \end{aligned}\]

From the optical Bloch equations, we find its equation of motion to be:

\[\begin{aligned} \dv{d}{t} &= \dv{\rho_{ee}}{t} - \dv{\rho_{gg}}{t} = 2 \gamma_g \rho_{gg} - 2 \gamma_e \rho_{ee} + \frac{i 2}{\hbar} \Big( \vb{p}^{-} \cdot \vb{E}^{+} - \vb{p}^{+} \cdot \vb{E}^{-} \Big) \end{aligned}\]

We can rewrite the first two terms in the following intuitive form, which describes a decay with rate \(\gamma_\parallel \equiv \gamma_g + \gamma_e\) towards an equilbrium \(d_0\):

\[\begin{aligned} 2 \gamma_g \rho_{gg} - 2 \gamma_e \rho_{ee} = \gamma_\parallel (d_0 - d) \qquad \quad d_0 \equiv \frac{\gamma_g - \gamma_e}{\gamma_g + \gamma_e} \end{aligned}\]

With this, the equation for the population inversion \(d\) takes the following final form:

\[\begin{aligned} \boxed{ \dv{d}{t} = \gamma_\parallel (d_0 - d) + \frac{i 2}{\hbar} \Big( \vb{p}^{-} \cdot \vb{E}^{+} - \vb{p}^{+} \cdot \vb{E}^{-} \Big) } \end{aligned}\]

Finally, we would like a relation between the polarization and the electric field \(\vb{E}\), for which we turn to Maxwell’s equations. We start from Faraday’s law, and split \(\vb{B} = \mu_0 (\vb{H} + \vb{M})\):

\[\begin{aligned} \nabla \cross \vb{E} = - \pdv{\vb{B}}{t} = - \mu_0 \pdv{\vb{H}}{t} - \mu_0 \pdv{\vb{M}}{t} \end{aligned}\]

We assume that there is no magnetization \(\vb{M} = 0\). Then we we take the curl of both sides, and replace \(\nabla \cross \vb{H}\) with Ampère’s circuital law:

\[\begin{aligned} \nabla \cross \big( \nabla \cross \vb{E} \big) = - \mu_0 \pdv{}{t} \big( \nabla \cross \vb{H} \big) = - \mu_0 \pdv{}{t} \Big( \vb{J}_\mathrm{free} + \pdv{\vb{D}}{t} \Big) \end{aligned}\]

Inserting the definition \(\vb{D} = \varepsilon_0 \vb{E} + \vb{P}\) together with Ohm’s law \(\vb{J}_\mathrm{free} = \sigma \vb{E}\) yields:

\[\begin{aligned} \nabla \cross \big( \nabla \cross \vb{E} \big) = - \mu_0 \sigma \pdv{\vb{E}}{t} - \mu_0 \varepsilon_0 \pdv[2]{\vb{E}}{t} - \mu_0 \pdv[2]{\vb{P}}{t} \end{aligned}\]

Where \(\sigma\) is the medium’s conductivity, if any; many authors assume \(\sigma = 0\). Next, we rewrite the left side using a vector identity, and assume no net charge \(\nabla \cdot \vb{E} = 0\):

\[\begin{aligned} \nabla^2 \vb{E} - \nabla \big( \nabla \cdot \vb{E} \big) = \nabla^2 \vb{E} = \mu_0 \sigma \pdv{\vb{E}}{t} + \mu_0 \varepsilon_0 \pdv[2]{\vb{E}}{t} + \mu_0 \pdv[2]{\vb{P}}{t} \end{aligned}\]

After some rearranging, we arrive at a variant of the electromagnetic wave equation:

\[\begin{aligned} \nabla^2 \vb{E} - \mu_0 \sigma \pdv{\vb{E}}{t} - \mu_0 \varepsilon_0 \pdv[2]{\vb{E}}{t} &= \mu_0 \pdv[2]{\vb{P}}{t} \end{aligned}\]

It is trivial to show that \(\vb{E}\) and \(\vb{P}\) can be replaced by \(\vb{E}^{+}\) and \(\vb{P}^{+}\). It is also simple to convert the dipole \(\vb{p}^{+}\) and inversion \(d\) into their macroscopic versions \(\vb{P}^{+}\) and \(D\):

\[\begin{aligned} \vb{P}^{+}(\vb{r}, t) = \sum_{n} \vb{p}^{+}_n \: \delta(\vb{r} \!-\! \vb{r}_n) \qquad \quad D(\vb{r}, t) = \sum_{n} d_n \: \delta(\vb{r} \!-\! \vb{r}_n) \end{aligned}\]

We thus arrive at the Maxwell-Bloch equations, which are relevant for laser theory:

\[\begin{aligned} \boxed{ \begin{aligned} \mu_0 \pdv[2]{\vb{P}^{+}}{t} &= \nabla^2 \vb{E}^{+} - \mu_0 \sigma \pdv{\vb{E}^{+}}{t} - \mu_0 \varepsilon_0 \pdv[2]{\vb{E}^{+}}{t} \\ \pdv{\vb{P}^{+}}{t} &= - \Big( \gamma_\perp + i \omega_0 \Big) \vb{P}^{+} - \frac{i}{\hbar} \Big( \vb{p}_0^{-} \cdot \vb{E}^{+} \Big) \vb{p}_0^{+} D \\ \pdv{D}{t} &= \gamma_\parallel (D_0 - D) + \frac{i 2}{\hbar} \Big( \vb{P}^{-} \cdot \vb{E}^{+} - \vb{P}^{+} \cdot \vb{E}^{-} \Big) \end{aligned} } \end{aligned}\]


  1. F. Kärtner, Ultrafast optics: lecture notes, 2005, MIT.
  2. H. Haken, Light: volume 2: laser light dynamics, 1985, North-Holland.
  3. H.J. Metcalf, P. van der Straten, Laser cooling and trapping, 1999, Springer.

© Marcus R.A. Newman, a.k.a. "Prefetch". Available under CC BY-SA 4.0.