Pedagogical Overview of Maxwell-Bloch Equations for Two Level Atoms Using Quantum Master Equation Approach

1. Introduction

In this article we will discuss about the basics of the Light Matter Interaction to be more precise, reviewing the pedagogical understanding of the Interaction of the Electromagnetic field coupled to a single atom and the extension for the ensemble of two level atoms situation. Here we will focus on the interaction of the two level atoms with the radiation Field. Before going for the full fledged quantum mechanical treatment we will discuss the semiclassical theory describing the interaction of the two-level atom with the electromagnetic field which is treated classically, a model popularly known as the Rabi Model. The dynamics of the atomic system can be described by the Master equation which Essentiality describes the time evolution of the density operator of the system and the optical Bloch equations describes the evolution of the matrix element of the density operator of the system. We will start with the Bloch Equation for the isolated Two level atom followed by the same for the interaction of the two level atom with the semiclassical Electromagnetic radiation Field before deriving the fully quantum mechanical Optical Bloch equation treating the radiation field to be Quantized. The density operator formulation plays a crucial role to describe the time evolution of the system. After that we will discuss the model describing the interaction of three level atoms with the electromagnetic field.

2. Modeling of the Hamiltonian for the Two Level Atoms

The two level atom can be visualized as a quantum system with two possible energy levels namely the ground and the excited state. Let us denote the ground state and excited state energies by $E_e$ and $E_g$. The ground and excited states represented as $|e\rangle$ and $|g\rangle$ are the eigen state of the atomic Hamiltonian with eigenvalues $E_e$ and $E_g$ respectively. So we have,

$$\widehat{H_A}|e\rangle =E_e|e\rangle ,\widehat{H_A}|g\rangle =E_g|g\rangle$$

(1)

The eigen states of the Hamiltonian forms an orthonormal basis in the 2 dimensional Hilbert space with the property of orthogonality and completeness relation listed as follows.

$$\begin{array}{c}\langle m|n\rangle ={\delta }_{mn}\mathrm{for}m,n=e,g\\ |e\rangle \langle e|+|g\rangle \langle g|={\widehat{I}}_{2\times 2}\end{array}$$

(2)

Now, using the spectral; decomposition theorem and a little bit of algebra we can write the Hamiltonian of the isolated two level atom as,

$$\widehat{H_A}=\left({\displaystyle \frac{E_e+E_g}{2}}\right)[|e\rangle \langle e|+|g\rangle \langle g|]+\left({\displaystyle \frac{E_e-E_g}{2}}\right)[|e\rangle \langle e|-|g\rangle \langle g|]$$

(3)

The first term of the equation will be proportional to the $I_{2\times 2}$ because of the completeness relation which can be dropped without loss of generality as the factor is tantamount to a shift in the energy. Now with the conventional choice of the basis vectors $|e\rangle =\left(\begin{array}{c}1\\ 0\end{array}\right)$ and $|g\rangle =\left(\begin{array}{c}0\\ 1\end{array}\right)$ followed by the introduction of the raising and lowering operators i.e $\widehat{{\sigma }^{+}}=|e\rangle \langle g|$ and the hermitian conjugate $\widehat{{\sigma }^{-}}=|g\rangle \langle e|$ we can write the Hamiltonian in the natural unit system as,

$$\widehat{H_A}={\displaystyle \frac{{\omega }_0}{2}}\widehat{{\sigma }^z}$$

(4)

Where we have defined the atomic transition frequency between the two levels as ${\omega }_0$ and used the following properties and definitions listed as follows,

$$\begin{array}{c}\widehat{{\sigma }^{+}}={\displaystyle \frac{1}{2}}\left(\widehat{{\sigma }^x}+i\widehat{{\sigma }^y}\right)=\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)\\ \widehat{{\sigma }^{-}}={\displaystyle \frac{1}{2}}\left(\widehat{{\sigma }^x}-i\widehat{{\sigma }^y}\right)=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\\ \widehat{{\sigma }^{+}}\widehat{{\sigma }^{-}}=|e\rangle \langle e|={\displaystyle \frac{1}{2}}(\widehat{I}+\widehat{{\sigma }^z})=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)\\ \widehat{{\sigma }^{-}}\widehat{{\sigma }^{+}}=|g\rangle \langle g|={\displaystyle \frac{1}{2}}(\widehat{I}-\widehat{{\sigma }^z})=\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)\\ \widehat{{\sigma }^{+}}\widehat{{\sigma }^{-}}-\widehat{{\sigma }^{-}}\widehat{{\sigma }^{+}}=[{\sigma }^{+},\widehat{{\sigma }^{-}}]=(|e\rangle \langle e|-|g\rangle \langle g|)=\widehat{{\sigma }^z}=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\\ \widehat{{\sigma }^{+}}\widehat{{\sigma }^{-}}+\widehat{{\sigma }^{-}}\widehat{{\sigma }^{+}}=\{{\sigma }^{+},\widehat{{\sigma }^{-}}\}=\widehat{I}\end{array}$$

(5)

It is important to note that, $\widehat{{\sigma }^{+}}|e\rangle =0$, $\widehat{{\sigma }^{-}}|g\rangle =0$, $\widehat{{\sigma }^{+}}|g\rangle =|e\rangle$, $\widehat{{\sigma }^{-}}|e\rangle =|g\rangle$. So in principle it is possible to introduce the fermionic creation and annihilation operators $\widehat{b}$ and ${\widehat{b}}^{†}$ with the mapping $\widehat{b}\equiv \widehat{{\sigma }^{-}}$ and ${\widehat{b}}^{†}\equiv \widehat{{\sigma }^{-}}$ with $\{\widehat{b},{\widehat{b}}^{†}\}=\widehat{I}$, $|0\rangle \equiv |g\rangle$, $|1\rangle \equiv |e\rangle$, $\{\widehat{b},\widehat{b}\}={\widehat{b}}^2=0$ and $\{{\widehat{b}}^{†},{\widehat{b}}^{†}\}={\widehat{b}}^{†2}=0$. Here the operators $\widehat{b}$ and ${\widehat{b}}^{†}$ are the fermionic creation and annihilation operators equivalent to the raising and lowering operators of a fermionic Harmonic oscillator and $|0\rangle$ , $|1\rangle$ are the two eigenstates of the Fermionic Harmonic oscillator hamiltonian with eigenvalues 0 and 1 respectiverly. This one to one correspondence between a fermionic harmonic oscillator and the two level atom is possible because of the resemblance of the anticommuation relations i.e. $\{\widehat{b},{\widehat{b}}^{†}\}=\widehat{I}$ alongwith $\{{\sigma }^{+},\widehat{{\sigma }^{-}}\}=\widehat{I}$.

2.1. A Brisk Tour of the Fermionic Harmonic Oscillator and Correspondence with the Single Isolated Two Level Atom

The hamiltonian of a fermionic harmonic oscillator as opposed to the bosonic harmonic oscillator can be modelled which is antisymmetric with respect to the operators $\widehat{b}$ and ${\widehat{b}}^{†}$ such that with,$\{\widehat{b},{\widehat{b}}^{†}\}=\widehat{I}$ and and $\{\widehat{b},\widehat{b}\}={\widehat{b}}^2=0$ and $\{{\widehat{b}}^{†},{\widehat{b}}^{†}\}={\widehat{b}}^{†2}=0$

$${\widehat{H}}_f={\displaystyle \frac{\hslash \omega }{2}}({\widehat{b}}^{†}\widehat{b}-\widehat{b}{\widehat{b}}^{†})$$

(6)

as opoosed to the completely usually known symmetric Bosonic Harmonic oscillator hamiltonian written as keeping in mind that $[\widehat{b},{\widehat{b}}^{†}]=\widehat{I}$

$$\widehat{H_B}={\displaystyle \frac{\hslash \omega }{2}}({\widehat{b}}^{†}\widehat{b}+\widehat{b}{\widehat{b}}^{†})=\hslash \omega ({\widehat{b}}^{†}\widehat{b}+{\displaystyle \frac{\widehat{I}}{2}})$$

(7)

The only difference is the commutation and anti-commutation relations used for the bosonic and fermionic systems. The former hamiltonian will have only two eigenstates denoted by $|0\rangle$ and $|1\rangle$ with eigenvalues 0 and 1 respectively after dropping the term in the fermionic harmonic oscillator hamiltonian proportional to identity matrix and using the fact that the fermionic number operator is involutory i.e.${({\widehat{b}}^{†}\widehat{b})}^2={\widehat{b}}^{†}\widehat{b}$ immediately follows from the anticommutaion relations introduced in the previous section. From now onwards the fermionic harmonic oscillator will be expressed as, ${\widehat{H}}_f=\omega {\widehat{b}}^{†}\widehat{b}$ in the natural unit system i.e. $\hslash =c=k_B=1.$ Now we can appreciate the correspondence between the harmonic oscillator and the two level atom with the mapping introduced before such that,

$${\widehat{H}}_A={\displaystyle \frac{{\omega }_0}{2}}(|e\rangle \langle e|-|g\rangle \langle g|)\equiv {\displaystyle \frac{{\omega }_0}{2}}({\widehat{b}}^{†}\widehat{b}-\widehat{b}{\widehat{b}}^{†})={\omega }_0{\widehat{b}}^{†}\widehat{b}$$

(8)

Where we have used the anti commuation relation and dropped the term proportional to the identity operator.

2.2. Discussion on the Optical Bloch Equation for the Isolated Two Level Atom

Before going to discuss about the optical bloch equation for the single isolated two levl atom we will discuss a bit more about the two level system and its connection to a single non interacting fermionic site. We start the discussion about the mapping of the Spin angular momentum operators with the Fermionic creation and annihilation operators. Because of the equivalence of the two level atom and spin-${\displaystyle \frac{1}{2}}$ we can define the operators with the mapping introduced before i.e. $\widehat{b}\equiv \widehat{{\sigma }^{-}}$ and ${\widehat{b}}^{†}\equiv \widehat{{\sigma }^{-}}$ as follows,

$$\begin{array}{c}{\widehat{R}}_1={\displaystyle \frac{1}{2}}\left({\widehat{\sigma }}^{+}+{\widehat{\sigma }}^{-}\right)={\displaystyle \frac{1}{2}}\left(\widehat{b}+{\widehat{b}}^{†}\right)\\ {\widehat{R}}_2={\displaystyle \frac{1}{2i}}\left({\widehat{\sigma }}^{+}-{\widehat{\sigma }}^{-}\right)={\displaystyle \frac{1}{2i}}\left({\widehat{b}}^{†}-\widehat{b}\right)\\ {\widehat{R}}^3={\displaystyle \frac{1}{2}}[{\widehat{\sigma }}^{+},{\widehat{\sigma }}^{-}]={\displaystyle \frac{1}{2}}\left[{\widehat{b}}^{†},\widehat{b}\right]={\displaystyle \frac{1}{2}}({\widehat{b}}^{†}\widehat{b}-\widehat{b}{\widehat{b}}^{†})\\ {\widehat{R}}_0={\displaystyle \frac{1}{2}}\widehat{I}\end{array}$$

(9)

It is evident from the above definitions that ${\widehat{R}}_i^s$$\forall \mathrm{i}=1,2,3$ are mutually non-commuting but ${\widehat{R}}_0$ commutes with each of ${\widehat{R}}_i$ with,

$$\begin{array}{c}[\widehat{R_i},\widehat{R_j}]=i{ϵ}_{ijk}{\widehat{R}}_k\\ \left[{\widehat{R}}_0,{\widehat{R}}_i\right]=0\end{array}$$

(10)

with, ${ϵ}_{ijk}$ being the completely antisymmetric Levi-civita symbol also known as the structure constant of the underlying SU(2) group. Now using the completeness relation of the eigenbasis of the atomic hamiltonian the density operator for the two level atom $\widehat{\rho }(t)$ will be a $2\times 2$ matrix which can be written as

$$\widehat{\rho }(t)={\rho }_{ee}(t)|e\rangle \langle e|+{\rho }_{gg}(t)|g\rangle \langle g|+{\rho }_{eg}(t)|e\rangle \langle g|+{\rho }_{ge}(t)|g\rangle \langle e|=\left(\begin{array}{cc}{\rho }_{ee}(t)& {\rho }_{eg}(t)\\ {\rho }_{ge}(t)& {\rho }_{gg}(t)\end{array}\right)$$

(11)

The general state of the two level atom can be described as the linear superposition of the hamiltonian’s eigensataes such that, $|\psi (t)\rangle =c_e(t)|e\rangle +c_g(t)|g\rangle$ and the density operator as $\widehat{\rho }(t)=|\psi \rangle (t)\langle \psi (t)|$, such that with the comparison we can write,

$${\rho }_{ee}(t)=|c_e{(t)|}^2,{\rho }_{gg}(t)={\left|c_g(t)\right|}^2,{\rho }_{eg}(t)=c_e(t)c_g{(t)}^{\ast },{\rho }_{ge}(t)=c_g(t)c_e{(t)}^{\ast }.$$

(12)

the conservation of the probability indicates that, $|c_e{(t)|}^2+{\left|c_g(t)\right|}^2=1=Tr[\widehat{\rho (t)}]$ With the general state of the system evolves according to the time dependent schrodinger equation i.e. $|\psi (t)\rangle =exp(-i\widehat{H}t)|\psi (0)\rangle$, with the hamiltonian being $\widehat{H}={\displaystyle \frac{{\omega }_0}{2}}\widehat{{\sigma }_z}$ the density operator of the two level atom will evolve in the schrodinger picture as,

$$\widehat{\rho }(t)=exp(-i\widehat{H}t)\widehat{\rho }(0)exp(-i\widehat{H}t)$$

(13)

The time evolution of the density operator can be described by the Von Neumann Louisville Master equation such that

$${\displaystyle \frac{d\widehat{\rho }(t)}{dt}}=i\left[\widehat{\rho (t)},\widehat{H}\right]$$

(14)

which leads to the optical bloch equations describing the time evolution of the matrix elements of the density operator of the two level atom in the eigen-basis of the matrix given by,

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{gg}(t)}{dt}}=0\\ {\displaystyle \frac{d{\rho }_{ee}(t)}{dt}}=0\\ {\displaystyle \frac{d{\rho }_{eg}(t)}{dt}}=-i{\omega }_0{\rho }_{eg}(t)\\ {\displaystyle \frac{d{\rho }_{ge}(t)}{dt}}=+i{\omega }_0{\rho }_{ge}(t)\end{array}$$

(15)

with the immediate solution leading to,

$${\rho }_{ee}(t)={\rho }_{ee}(0),{\rho }_{gg}(t)={\rho }_{gg}(0),{\rho }_{eg}(t)={\rho }_{eg}(0)exp(-i{\omega }_{0}t),{\rho }_{ge}(t)={\rho }_{ge}(0)exp(+i{\omega }_{0}t).$$

(16)

It can be easily seen that the trace of the density matrix is preserved in time which is quite obvious. The determinant of the density matrix is also preserved in time which means that if the syatem initailly starts in a pure state then it will remain in the pure state at any future time such that,

$$Det(\widehat{\rho }(t))={\rho }_{ee}(t){\rho }_{gg}(t)-{\left|{\rho }_{eg}(t)\right|}^2=Det(\widehat{\rho }(0)).$$

(17)

In general the density matrix of the two level atom can be expressed as a linear combination of the $I_{2\times 2}$ and the generators of the SU(2) group i.e. the three pauli matrices such that,

$$\widehat{\rho }(t)={\displaystyle \frac{1}{2}}\left(\widehat{I}+\overrightarrow{n}(t)·\widehat{\overrightarrow{\sigma }}\right)$$

(18)

with, a triplet $n_i(t)$ defines the three dimensional Bloch vector $\overrightarrow{n}(t)$ with, $n_i(t)=Tr[\widehat{\rho }(t){\widehat{\sigma }}_i]$ for i=1,2,3. more explicitly we can write that,

$$n_1(t)={\rho }_{eg}(t)+{\rho }_{ge}(t)=2\mathfrak{Re}({\rho }_{ge}(t)),n_2(t)=\mathfrak{Im}({\rho }_{eg}(t)),n_3(t)=({\rho }_{ee}(t)-{\rho }_{gg}(t)).$$

(19)

If the state of the system is a pure state then, with $\widehat{\rho }{(t)}^2=\widehat{\rho }(t)$ we must have $\overrightarrow{n}(t)·\overrightarrow{n}(t)=1.$ which implies that ${\sum }_{i=1}^3{n}_i{(t)}^2=1$.the magnitude of the bloch vector can be found as,

$$n{(t)}^2=\sqrt{1+4\left\{|{\rho }_{ge}{(t)|}^2-{\rho }_{ee}(t){\rho }_{gg}(t)\right\}}$$

(20)

We can see that the population at the ground state and excited state defined as the diagonal matrix elements of the density operator i.e. ${\rho }_{ee}(t)$ and ${\rho }_{gg}(t)$ does not change with time but the off diagonal matrix elements which defines the quantum coherence oscillates with time but is not going to destroy in future time. It is sometimes useful in quantum optics to define the Population Inversion defined as,

$$P(t)={\rho }_{ee}(t)-{\rho }_{gg}(t)=2{\rho }_{ee}(t)-1.$$

(21)

As the population does not change with time we can write that $P(t)=P(0)$. And with this we end our discussion of the optical Bloch equations of the Isolated Two level atom which is non interacting in nature.

3. Semiclasical and Full Quantum Mechanical Model of Light Matter Interaction

Before going for the discussion of the interaction of the atom to be more precise two level atom with a quantized radiation field we will discuss about the semiclassical model of light matter interaction which essentially describes the interaction of a two level atom with a classical radiation field i.e. the electric field associated with the incident Electromagnetic radiation field is being treated classically. Such a semi classical description is known as the Rabi model and the corresponding quantum mechanical treatment without using the Rotating Wave approximation is called the Quantum Rabi Model which we will not discuss in the present context.

3.1. Semiclassical Rabi Model

The atom considered in the present context is a two level atom and it can be either a polar or a non-polar atom. If the atom does not possess a permanent dipole moment it will be called a non-polar atom otherwise polar. If the atom is non polar then due to the incident radiation field mainly due to the incident electric field there will be a displacement of the center of gravity of the negative charges with respect to that of the positive charges and we can define the induced dipole moment as $\overrightarrow{\mu }=e\overrightarrow{r}$ where e being the electronic charge and $\overrightarrow{r}$ being the displacement vector. By definition the dipole moment vector is a polar vector under the parity operation and the corresponding dipole moment vector operator will be a polar vector operator such that,

$$\widehat{\mathsf{\Pi }}\widehat{\overrightarrow{\mu }}{\widehat{\mathsf{\Pi }}}^{†}=-\widehat{\overrightarrow{\mu }}$$

(22)

With, the parity operator defined as, $\widehat{\mathsf{\Pi }}={(-1)}^{{\widehat{\sigma }}^{+}{\widehat{\sigma }}^{-}}=exp(i\pi {\widehat{\sigma }}^{+}{\widehat{\sigma }}^{-})$ with the desirable properties, ${\widehat{\mathsf{\Pi }}}^2=\widehat{I}$ and ${\widehat{\mathsf{\Pi }}}^{†}=\widehat{\mathsf{\Pi }}$ alongside $\widehat{\mathsf{\Pi }}{\widehat{\mathsf{\Pi }}}^{†}=\widehat{I}$ and ${\widehat{\mathsf{\Pi }}}^{-1}=\widehat{\mathsf{\Pi }}$. It is also evident from the definition that the eigenvalues of the parity operators will be $\pm 1$. As things stands we can say that the Parity operator commutes with the hamiltonian of the two level atom in the absence of the interaction which means they have the simultaneous eigen-states. Then we can say that the eigensates of the atomic hamiltonian will have definite parity either even or odd. Now if the radiation field is being treated classically the interaction hamiltonian can be written in the dipolar approximation with the classical electric field $\overrightarrow{E}(t)={\overrightarrow{E}}_{0}cos\omega t$ as,

$${\widehat{H}}_{int}(t)=-\widehat{\overrightarrow{\mu }}.\overrightarrow{E(t)}=-\widehat{\overrightarrow{\mu }}.\overrightarrow{E_0}cos\omega t$$

(23)

Then the Hamiltonian of the two level atom coupled with the radiation field will now become,

$${\widehat{H}}_{AF}(t)={\widehat{H}}_A+{\widehat{H}}_{int}(t)={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z--\widehat{\overrightarrow{\mu }}.\overrightarrow{E_0}cos\omega t$$

(24)

Using the aforementioned properties of the parity operator and completeness relation of the eigenbasis of the unperturbed hamiltonian i.e. ${\widehat{H}}_A$ we can write that,

$$\widehat{\overrightarrow{\mu }}={\overrightarrow{\mu }}_{eg}|e\rangle \langle g|+{\overrightarrow{\mu }}_{ge}|g\rangle \langle e|={\overrightarrow{\mu }}_{eg}{\widehat{\sigma }}^{+}+{\overrightarrow{\mu }}_{ge}{\widehat{\sigma }}^{-}$$

(25)

with, ${\overrightarrow{\mu }}_{ij}=\langle i|\widehat{\overrightarrow{\mu }}|j\rangle$ where, $i,j=e,g$ and due to the parity argument ${\overrightarrow{\mu }}_{ee}={\overrightarrow{\mu }}_{gg}=\overrightarrow{0}$. Now with, $V_{eg}=-{\overrightarrow{\mu }}_{eg}·{\overrightarrow{E}}_0$ that the interaction hamiltonian can be written as,

$${\widehat{H}}_{int}(t)=V_{eg}{\widehat{\sigma }}^{+}cos(\omega t)+h.c$$

(26)

h.c is the acronym for hermitian conjugate. Now as the electric field is treated classically essentially there is no bath as such and the total hamiltonian of the atom coupled with the radiation field is essentially a Hamiltonian consisting of the system operators only so that, the density operator of the system will evolve with the Hamiltonian ${\widehat{H}}_{AF}(t)$ as,

$$\widehat{\rho }(t)=\widehat{U}(t)\widehat{\rho }(0){\widehat{U}}^{†}(t)$$

(27)

with, $\widehat{U}(t,0)=exp\left(-{\displaystyle \frac{i}{\hslash }}{\int }_0^t{\widehat{H}}_{AF}(t^{\prime })d{t}^{\prime }\right)$ being the unitary time evolution operator satisfying $\widehat{U}(t){\widehat{U}}^{†}(t)=\widehat{I}$. Further analysis of the dynamics of this system can be dine by using the standard time dependent perturbation theory with the general state of the system can be written as,

$$|\mathsf{\Psi }(t)\rangle =c_e(t)|e\rangle e^{-{\displaystyle \frac{i{E}_{e}t}{\hslash }}}+c_g(t)|g\rangle e^{-{\displaystyle \frac{i{E}_{g}t}{\hslash }}}$$

(28)

The normalization of the state vector demands that, $|c_e{(t)|}^2+{\left|c_g(t)\right|}^2=1$. And we can define the density operator of the system as $\widehat{\rho }(t)=|\mathsf{\Psi }(t)\rangle \langle \mathsf{\Psi }(t)|$ with $|\mathsf{\Psi }(t)\rangle$ evolves according to the time dependent Schrodinger equation given as,

$$i\hslash {\displaystyle \frac{d}{dt}}|\mathsf{\Psi }(t)\rangle ={\widehat{H}}_{AF}(t)|\mathsf{\Psi }(t)\rangle$$

(29)

Using the time dependent perturbation theory we can write in general,

$$i\hslash {\displaystyle \frac{d{c}_m(t)}{dt}}=\sum _n{c}_n(t)H_{int}^{mn}(t)e^{i{\omega }_{mn}t}$$

(30)

with, $H_{int}^{mn}(t)$ defines the matrix element of the perturbing term of the hamiltonian carrying the explicit time dependence which in our case is the ${\widehat{H}}_{int}(t)$ in the eigen basis of the unperturbed Hamiltonian which is ${\widehat{H}}_A$ in the present context such that, ${\widehat{H}}_{int}^{mn}(t)=\langle n|{\widehat{H}}_{int}(t)|m\rangle$ with, $H_{int}^{mn}{(t)}^{\ast }=H_{int}^{mn}(t)$, due to hermiticity and $n,m$ both can have the choice $e,g$. Also to note that we have defined the transition frequency as ${\omega }_{mn}={\displaystyle \frac{(E_m-E_n)}{\hslash }}$, with $E_m$’s are the unperturbed energy eigenvalues. In our case the equation for $c_e(t)$ and $c_g(t)$ will be as follows,

$$\begin{array}{c}i\hslash {\displaystyle \frac{d{c}_e(t)}{dt}}=c_e(t)H_{int}^{ee}(t)+c_g(t)H_{int}^{eg}(t)e^{i{\omega }_{0}t}=c_g(t)H_{int}^{eg}(t)e^{i{\omega }_{0}t}=c_g(t)v_{eg}cos(\omega t)e^{i{\omega }_{0}t}\end{array}$$

(31)

$$\begin{array}{c}i\hslash {\displaystyle \frac{d{c}_g(t)}{dt}}=c_g(t)H_{int}^{gg}(t)+c_e(t)H_{int}^{ge}(t)e^{-i{\omega }_{0}t}=c_e(t)H_{int}^{ge}(t)e^{-i{\omega }_{0}t}=c_e(t)v_{eg}^{\ast }cos(\omega t)e^{-i{\omega }_{0}t}\end{array}$$

(32)

with, ${\omega }_{eg}=(E_e-E_g)/\hslash ={\omega }_0$, the atomic transition frequency. with the fact that the diagonal matrix elements of ${\widehat{H}}_{int}(t)$ in the basis $\{|e\rangle ,|g\rangle \}$ will vanish such that, $H_{int}^{ee}(t)=H_{int}^{ee}(t)=0$. Now we will apply the rotating wave approximation for the exact solution of the problem which goes like this, in the bove equations there are terms like $cos(\omega t)e^{\pm i{\omega }_{0}t}$ which can be decomposed into two parts, one with $e^{\pm i(\omega +{\omega }_0)t}$ and $e^{\pm i(\omega -{\omega }_0)t}$. Near the resonance i.e. $\omega \approx {\omega }_0$the terms of the type $e^{\pm i(\omega +{\omega }_0)t}$ will be rapidly oscillating compared to the terms with terms $(\omega -{\omega }_0)$, commonly called the anti-resonant term and will be dropped for the simplicity of the calculation. Which leads to the following set of differential equations given by,

$$\begin{array}{c}i\hslash {\displaystyle \frac{d{c}_e(t)}{dt}}=c_g(t)v_{eg}{e}^{-i(\omega -{\omega }_0)t}\end{array}$$

(33)

$$\begin{array}{c}i\hslash {\displaystyle \frac{d{c}_g(t)}{dt}}=c_e(t)v_{eg}^{\ast }{e}^{i(\omega -{\omega }_0)t}\end{array}$$

(34)

The equations for $c_e(t)$ and $c_g(t)$ are coupled first-order differential equations which can be converted into uncoupled second-order differential equations for $c_e(t)$ and $c_g(t)$ respectively. If we assume that the atom was initially in the ground state then with the boundary condition $c_g(0)=1$ and $c_e(0)=0$ the solution of the differential equations can be obtained. If the atom was initially in the ground state then $|c_e{(t)|}^2$ will define the transition probability from the ground state to the excited state denoted by $P_{g\to e}(t)$. After decoupling the equations for $c_e(t)$ and $c_g(t)$ we get the following second order homogeneous differential equations as follows,

$$\begin{array}{c}{\displaystyle \frac{d^2{c}_e(t)}{d{t}^2}}+i\left({\displaystyle \frac{\omega -{\omega }_0}{2}}\right){\displaystyle \frac{d{c}_e(t)}{dt}}+c_e(t){\displaystyle \frac{|v_{eg}{|}^2}{4{\hslash }^2}}=0\\ {\displaystyle \frac{d^2{c}_g(t)}{d{t}^2}}-i\left({\displaystyle \frac{\omega -{\omega }_0}{2}}\right){\displaystyle \frac{d{c}_g(t)}{dt}}+c_g(t){\displaystyle \frac{|v_{eg}{|}^2}{4{\hslash }^2}}=0\end{array}$$

(35)

The solution of the above equations are given by,

$$\begin{array}{c}{c}_e(t)=e^{-i(\omega -{\omega }_0)t/4}\left(-{\displaystyle \frac{2i{v}_{eg}}{{\omega }_R\hslash }}\right)sin\left({\displaystyle \frac{{\omega }_{R}t}{4}}\right)\end{array}$$

(36)

$$\begin{array}{c}{c}_g(t)=e^{3i(\omega -{\omega }_0)t/4}\left[sin\left({\displaystyle \frac{{\omega }_{R}t}{4}}\right)-i\left({\displaystyle \frac{\omega -{\omega }_0}{{\omega }_R}}\right)sin\left({\displaystyle \frac{{\omega }_{R}t}{4}}\right)\right]\end{array}$$

(37)

And the probability to find the system in the ground and excited states are respectively given by,

$$\begin{array}{c}{P}_e(t)={\left|c_e(t)\right|}^2={\displaystyle \frac{4|v_{eg}{|}^2}{{\hslash }^2{\omega }_R^2}}{sin}^2\left({\displaystyle \frac{{\omega }_{R}t}{4}}\right)\end{array}$$

(38)

$$\begin{array}{c}{P}_g(t)={\left|c_g(t)\right|}^2={cos}^2\left({\displaystyle \frac{{\omega }_{R}t}{4}}\right)+{\left({\displaystyle \frac{\omega -{\omega }_0}{{\omega }_R}}\right)}^2{sin}^2\left({\displaystyle \frac{{\omega }_{R}t}{4}}\right)\end{array}$$

(39)

Where we have the Rabi oscillation frequency as,

$${\omega }_R=\sqrt{{\left(\omega -{\omega }_0\right)}^2+4{\displaystyle \frac{|v_{eg}{|}^2}{{\hslash }^2}}}$$

(40)

With the definition of the Rabi frequency we can easily check that the total probability is 1 i.e. $|c_e{(t)|}^2+{\left|c_g(t)\right|}^2=1$ and the transition probability from the ground state to the excited state will be,

$$P_{g\to e}(t)={\left|c_e(t)\right|}^2={\displaystyle \frac{4|v_{eg}{|}^2}{{\hslash }^2{\omega }_R^2}}{sin}^2\left({\displaystyle \frac{{\omega }_{R}t}{4}}\right)$$

(41)

As we can see that the probabilities oscillates with time because of the nature of the solution with a characteristics frequency defined as Rabi Frequency denoted by ${\omega }_R$.

3.2. Density Operator Approach to the Rabi Problem

Now before going to discuss about the Density operator formulation of the Semi classical Rabi problem we will discuss some basic prerequisites. In the Heisenberg picture the Time evolution of any operator is defined as,

$${\widehat{A}}^H(t)=exp(i\widehat{H}t){\widehat{A}}^{S}exp(-i\widehat{H}t)$$

(42)

with $\hslash =1$.Here, ${\widehat{A}}^S$ denotes the operator in the Schrodinger picture and $\widehat{H}$ being the hamiltonian of the quantum system. Using this we can find the equation of motion of the Dipole moment operator of the two level atom in the Heisenberg picture for the isolated two level atom. From now onwards we will take $\hslash =1$. With, $\widehat{\overrightarrow{\mu }}(t)=e\widehat{\overrightarrow{r}}(t)$, in the heisenberg picture we can write,

$${\displaystyle \frac{d\widehat{\overrightarrow{\mu }}(t)}{dt}}=e\widehat{\overrightarrow{V}}(t)=i\left[\widehat{H_A},\widehat{\overrightarrow{\mu }}(t)\right]=-i{\omega }_0[\widehat{{\sigma }^{-}}(t){\overrightarrow{\mu }}_{ge}-{\overrightarrow{\mu }}_{ge}^{\ast }{\widehat{\sigma }}^{+}(t)]$$

(43)

With, $\widehat{H_A}$ being the hamiltonian of the isolated two level atom. With the free evolution of the operators as ${\widehat{\sigma }}^{+}(t)={\widehat{\sigma }}^{+}{e}^{i{\omega }_{0}t}$ and ${\widehat{\sigma }}^{-}(t)={\widehat{\sigma }}^{-}{e}^{-i{\omega }_{0}t}$ we can directly write that,

$$\widehat{\overrightarrow{\mu }}(t)={\overrightarrow{\mu }}_{ge}(t){\widehat{\sigma }}^{-}{e}^{-i{\omega }_{0}t}+h.c$$

(44)

We have defined the Three dimensional bloch vector $\overrightarrow{n}(t)$. Every possible pure states of the two level system will lie on the surface of the three dimensional bloch sphere and the mixed state will lie in the interior of the three dimensional bloch sphere. Th bloch vector essentially defines the position vector of the point t which the state of the quantum system is mapped either on the surface or in the interior of the three dimensional sphere. Then in the spherical polar coordinate system in general it is possible to parameterize the bloch vector components as,

$$n_1=nsin\theta cos\varphi ,n_2=nsin\theta sin\varphi ,n_3=ncos\theta \mathrm{with}0\le \theta \le \pi ,0\le \varphi \le 2\pi .$$

(45)

For the pure states which are mapped on the surface of the sphere the bloch vector will be an unit vector with $n=1$. To be more precise the pure states will be mapped on the surface of the bloch sphere of radius unity. By the definition of the Bloch vector componenets introduced before we can write the matrix elements of the density operator for the two level tom in terms of $n,\theta ,\varphi$ as follows,

$$\begin{array}{c}{\rho }_{ge}={\displaystyle \frac{1}{2}}nsin\theta e^{i\varphi }\\ {\rho }_{eg}={\displaystyle \frac{1}{2}}nsin\theta e^{-i\varphi }\\ {\rho }_{ee}={\displaystyle \frac{1}{2}}\left(1+ncos\theta \right)\\ {\rho }_{gg}={\displaystyle \frac{1}{2}}(1-ncos\theta )\\ \varphi =arctan\left({\displaystyle \frac{\mathfrak{Im}({\rho }_{ge})}{\mathfrak{Re}({\rho }_{ge})}}\right)\\ cos\theta ={\displaystyle \frac{{\rho }_{ee}-{\rho }_{gg}}{n}}\end{array}$$

(46)

With, $|\mathsf{\Psi }(t)\rangle$ denoting the general state of the two level atom written as a linear combination of $|e\rangle$ and $|g\rangle$ we can write, $|\mathsf{\Psi }(t)\rangle =c_e(t)|e\rangle +c_g(t)|g\rangle$ with, $|\mathsf{\Psi }(t)\rangle =e^{-i\widehat{H_{A}t}}|\mathsf{\Psi }(0)\rangle$. With ${\rho }_{ee}={\left|c_e\right|}^2$ and ${\rho }_{gg}={\left|c_g\right|}^2$ we can write the initial state of the system as $|\mathsf{\Psi }(0)\rangle =sin{\displaystyle \frac{\theta }{2}}{e}^{i\varphi /2}|g\rangle +cos{\displaystyle \frac{\theta }{2}}{e}^{-i\varphi /2}|e\rangle$ such that the state evolved over time at any future time $t>0$ can be written as,

$$|\mathsf{\Psi }(t)\rangle =sin{\displaystyle \frac{\theta }{2}}{e}^{i(\varphi +{\omega }_0)t/2}|g\rangle +cos{\displaystyle \frac{\theta }{2}}{e}^{-i(\varphi +{\omega }_0)t/2}|e\rangle$$

(47)

The expectation value of the dipole moment operator and the squared uncertainty will be,

$$\begin{array}{c}\langle \widehat{\overrightarrow{\mu }}\rangle =n_1\mathfrak{Re}({\overrightarrow{\mu }}_{ge})+n_2\mathfrak{Im}({\overrightarrow{\mu }}_{ge}).\\ {\left(\Delta \widehat{\mu }\right)}^2={\left|\overrightarrow{{\mu }_{eg}}\right|}^2\end{array}$$

(48)

Now we are in a position to discuss the density operator formulation to study the dynamics of the two level tom interacting with the classical radiation field. The optical Bloch equations for the semiclassical description wil be,

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{gg}(t)}{dt}}={\displaystyle \frac{1}{i}}\left[{\rho }_{eg}(t)\langle g|{\widehat{H}}_{int}(t)|e\rangle -c.c]\right]\end{array}$$

(49)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{ee}(t)}{dt}}={\displaystyle \frac{1}{i}}\left[{\rho }_{eg}(t)\langle e|{\widehat{H}}_{int}(t)|g\rangle -c.c\right]\end{array}$$

(50)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{eg}(t)}{dt}}={\displaystyle \frac{1}{i}}\left[{\rho }_{eg}(t){\omega }_0+\left({\rho }_{gg}(t)-{\rho }_{ee}(t)\right)\langle e|{\widehat{H}}_{int}(t)|g\rangle \right]\end{array}$$

(51)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{ge}(t)}{dt}}={\displaystyle \frac{1}{i}}\left[-{\rho }_{ge}(t){\omega }_0+\left({\rho }_{ee}(t)-{\rho }_{gg}(t)\right)\langle g|{\widehat{H}}_{int}(t)|e\rangle \right]\end{array}$$

(52)

The following equations can also be expressed in terms of the time evolution of the bloch vector components with the later also being called the optical bloch equations which are given by,

$$\begin{array}{c}{\dot{n}}_1(t)={\displaystyle \frac{d{n}_1(t)}{dt}}=-{\omega }_0{n}_2+2n_3\mathfrak{Im}\langle g|{\widehat{H}}_{int}(t)|e\rangle \end{array}$$

(53)

$$\begin{array}{c}{\dot{n}}_2(t)={\displaystyle \frac{d{n}_2(t)}{dt}}=+{\omega }_0{n}_1-2n_3\mathfrak{Im}\langle g|{\widehat{H}}_{int}(t)|e\rangle \end{array}$$

(54)

$$\begin{array}{c}{\dot{n}}_3(t)={\displaystyle \frac{d{n}_3(t)}{dt}}=-2\mathfrak{Im}\langle g|{\widehat{H}}_{int}(t)|e\rangle n_1+2n_2\mathfrak{Re}\langle g|{\widehat{H}}_{int}(t)|e\rangle \end{array}$$

(55)

From the above equations it follows immediately that,

$${\displaystyle \frac{d}{dt}}\left(n_1^2(t)+n_2^2(t)+n_3^2(t)\right)=0$$

(56)

Which means that the length of the bloch vector does not change with time. it is sometimes important to define another three dimensional vector $\overrightarrow{Q}(t)$ with,

$$\begin{array}{c}{Q}_1(t)=2\mathfrak{Re}\langle g|{\widehat{H}}_{int}(t)|e\rangle \end{array}$$

(57)

$$\begin{array}{c}{Q}_2(t)=2\mathfrak{Im}\langle g|{\widehat{H}}_{int}(t)|e\rangle \end{array}$$

(58)

$$\begin{array}{c}{Q}_3(t)={\omega }_0\end{array}$$

(59)

Then the equations for the bloch vector components can be writen in terms of the components $Q_i^{\prime }$s as follows,

$$\begin{array}{c}{\displaystyle \frac{d{n}_1(t)}{dt}}=-Q_3(t)n_2(t)+n_3{Q}_2(t)\end{array}$$

(60)

$$\begin{array}{c}{\displaystyle \frac{d{n}_2(t)}{dt}}=Q_3(t)n_1(t)-n_3(t)Q_1(t)\end{array}$$

(61)

$$\begin{array}{c}{\displaystyle \frac{d{n}_3(t)}{dt}}=-n_1(t)Q_2(t)+Q_1(t)n_2(t).\end{array}$$

(62)

The above set of equations can be summarised as,

$${\displaystyle \frac{d\overrightarrow{n}(t)}{dt}}=\overrightarrow{Q}(t)\times \overrightarrow{n}(t)$$

(63)

Now let us look at the interaction hamiltonian defined before in the Dipolar approximation. Now it is important to note that the matrix element of the dipole moment operator can be in general complex. If we consider the atomic transition corresponding to the selection rule given by $\Delta m=0$ then the matrix element i.e. ${\overrightarrow{\mu }}_{eg}={\overrightarrow{\mu }}_{ge}={\overrightarrow{\mu }}_{eg}^{\ast }$ which means the element is real but corresponding to the transitions given by the selection rule $\Delta m=\pm 1$ the matrix element will be in general complex. First let us consider the transition in accordance with the dipole selection rule $\Delta m\pm 1$ for which ${\overrightarrow{\mu }}_{ge}$ is complex.Now with the classical electric field given by, $\overrightarrow{E}(t)=\widehat{e}\left|\epsilon (t)\right|e^{i\varphi (t)}{e}^{-i\omega t}+c.c$, which contains a time dependent amplitude and a time independent phase factor the interaction hamiltonian will be,

$${\widehat{H}}_{int}(t)=-{\widehat{\sigma }}^{-}{\overrightarrow{\mu }}_{ge}·\widehat{e}\left|\epsilon (t)\right|e^{i\varphi (t)}{e}^{-i\omega t}+h.c-{\widehat{\sigma }}^{-}{\overrightarrow{\mu }}_{ge}·{\widehat{e}}^{\ast }\left|\epsilon (t)\right|e^{-i\varphi (t)}{e}^{i\omega t}+h.c$$

(64)

It is important to note that the spatial variation of the electric field over the atomic length scale can be neglected using the long wavelength approximation such that the wavelength of the electromagnetic radiation is very large compared to the bohr radius of the atom. We choose that, ${\overrightarrow{\mu }}_{ge}=\left|{\overrightarrow{\mu }}_{ge}\right|{\displaystyle \frac{(\widehat{x}+i\widehat{y})}{\sqrt{2}}}$ and $\widehat{e}={\displaystyle \frac{(\widehat{x}+i\widehat{y})}{2}}$ where $\widehat{x}$ and $\widehat{y}$ are the unit vectors and $\widehat{e}$ being the polarization vector which describes the direction of the incident electric field associated with the radiation field. We have assumed that the wave is propagating along the z direction which justifies the choice of $\widehat{e}$ in the section. From the given it follows immediately that the real and imaginary part of the off-diagonal matrix elements of the Interaction hamiltonian are respectively,

$$\begin{array}{c}\mathfrak{Re}\langle g|{\widehat{H}}_{int}(t)|e\rangle =-|\overrightarrow{{\mu }_{ge}}\left|\right|\epsilon (t)|cos\left(\omega t-\varphi (t)\right)\\ \mathfrak{Im}\langle g|{\widehat{H}}_{int}(t)|e\rangle =|{\overrightarrow{\mu }}_{ge}\left|\right|\epsilon (t)|sin\left(\omega t-\varphi (t)\right)\end{array}$$

(65)

Now let us define the rabi frequency as, ${\mathsf{\Omega }}_R(t)=2{\overrightarrow{\mu }}_{eg}·\widehat{e}\left|\epsilon (t)\right|$ such that the equations for the bloch vector components will become,

$$\begin{array}{c}{\dot{n}}_1(t)=-{\omega }_0{n}_2(t)+n_3(t){\mathsf{\Omega }}_R(t)sin(\omega t-\varphi (t))\end{array}$$

(66)

$$\begin{array}{c}{\dot{n}}_2(t)=-{\omega }_0{n}_1(t)+n_3(t){\mathsf{\Omega }}_R(t)cos(\omega t-\varphi (t))\end{array}$$

(67)

$$\begin{array}{c}{\dot{n}}_3(t)=-n_1(t){\mathsf{\Omega }}_R(t)sin\left(\omega t-\varphi (t)\right)-n_2(t)cos\left(\omega t-\varphi (t)\right)\end{array}$$

(68)

Now, if the transition is consider under the selection rule $\Delta m=0$ then the matrix element of the dipole moment operator will be real and if we assume the Polarisation vector i.e. $\widehat{e}$ to be real then the modified equations for the time evolution of the bloch vector components will be given by,

$$\begin{array}{c}{\dot{n}}_1(t)=-{\omega }_0{n}_2(t)\end{array}$$

(69)

$$\begin{array}{c}{\dot{n}}_2(t)={\omega }_0{n}_1(t)+2{\mathsf{\Omega }}_R(t)+n_3(t)cos(\omega t-\varphi (t))\end{array}$$

(70)

$$\begin{array}{c}{\dot{n}}_3(t)=-2n_2(t){\mathsf{\Omega }}_R(t)cos(\omega t-\varphi (t))\end{array}$$

(71)

Now we summarise the componenets of the vector $\overrightarrow{Q}(t)$ for the cases with $\Delta m=0$ and $\Delta m\pm 1$ respectively as follows, For, $\Delta m\pm 1$ we have,

$$\begin{array}{c}{Q}_1(t)=-2{\mathsf{\Omega }}_R(t)cos\left(\omega t-\varphi (t)\right)\end{array}$$

(72)

$$\begin{array}{c}{Q}_2(t)=2{\mathsf{\Omega }}_R(t)sin\left(\omega t-\varphi (t)\right)\end{array}$$

(73)

$$\begin{array}{c}{Q}_3(t)={\omega }_0\end{array}$$

(74)

And for $\Delta m=0$ the components will be,

$$\begin{array}{c}{Q}_1(t)=-2{\mathsf{\Omega }}_R(t)cos\left(\omega t-\varphi (t)\right)\end{array}$$

(75)

$$\begin{array}{c}{Q}_2(t)=0\end{array}$$

(76)

$$\begin{array}{c}{Q}_3(t)={\omega }_0\end{array}$$

(77)

As we can see that the Bloch vector rotates in the three dimensional coordinate system with an angular velocity $|Q(t)|$ given by the magnitude of the $\overrightarrow{Q}(t)$. Now we will explicitly write the optical Bloch equation in this semiclassical prescription when the electromagnetic field is being treated classically. The semiclassical way of describing the dynamics of the Two level system is incomplete in the sense that it does not incorporate the quantum effects of spontaneous emission which can not be captured in the Bloch equations because the radiation field is not treated quantum mechanically. The spontaneous emission effect can be incorporated in the semi classically derived optical Bloch Equations though we will prove them in a separate section later using the full fledged Quantum Treatment later on. Still this approach plays a vital role in term of developing the understanding of the light mater interaction models. Now coming back to the discussion of Semi classical rabi Model using the time dependent perturbation theory we had obtained using the Rotating Wave approximation that,

$$\begin{array}{c}i\hslash {\displaystyle \frac{d{c}_e(t)}{dt}}=c_g(t)v_{eg}{e}^{-i(\omega -{\omega }_0)t}\end{array}$$

(78)

$$\begin{array}{c}i\hslash {\displaystyle \frac{d{c}_g(t)}{dt}}=c_e(t)v_{eg}^{\ast }{e}^{i(\omega -{\omega }_0)t}\end{array}$$

(79)

This equations along with the fact that ${\rho }_{ee}(t)={\left|c_e(t)\right|}^2$, ${\rho }_{gg}(t)={\left|c_g(t)\right|}^2$, ${\rho }_{eg}(t)=c_e(t)c_g^{\ast }(t)$, ${\rho }_{ge}(t)=c_g(t)c_e^{\ast }(t)$ and $\hslash =1$ we get,

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{ee}(t)}{dt}}={\displaystyle \frac{-i{\omega }_R}{2}}\left[{\tilde{\rho }}_{ge}(t)-{\tilde{\rho }}_{eg}(t)\right]\end{array}$$

(80)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{gg}(t)}{dt}}={\displaystyle \frac{-i{\omega }_R}{2}}\left[{\tilde{\rho }}_{eg}(t)-{\tilde{\rho }}_{ge}(t)\right]\end{array}$$

(81)

$$\begin{array}{c}{\displaystyle \frac{d{\tilde{\rho }}_{eg}(t)}{dt}}=i\Delta {\tilde{\rho }}_{eg}(t)+{\displaystyle \frac{i{\omega }_R}{2}}\left[{\rho }_{ee}(t)-{\rho }_{gg}(t)\right]\end{array}$$

(82)

$$\begin{array}{c}{\displaystyle \frac{d{\tilde{\rho }}_{ge}(t)}{dt}}=-i\Delta {\tilde{\rho }}_{eg}(t)+{\displaystyle \frac{i{\omega }_R}{2}}\left[{\rho }_{gg}(t)-{\rho }_{ee}(t)\right]\end{array}$$

(83)

Where we have defined the Rabi Frequency as, ${\omega }_R=-({\overrightarrow{\mu }}_{eg}·{\overrightarrow{E}}_0)$ with $\hslash =1$ and $\Delta$ defines the detuning frequnecy such that, $\Delta =\omega -{\omega }_0$. It is important to note that we have asumed the matrix element of the dipole moment operator is real such that,${\overrightarrow{\mu }}_{eg}={\overrightarrow{\mu }}_{ge}={\overrightarrow{\mu }}_{eg}^{\ast }$ and We have defined that,

$${\tilde{\rho }}_{ge}(t)={\rho }_{ge}(t)e^{-i\Delta t},{\tilde{\rho }}_{eg}(t)={\rho }_{eg}(t)e^{i\Delta t},{\tilde{\rho }}_{eg}(t)={\tilde{\rho }}_{ge}{(t)}^{\ast }$$

(84)

From the equations it can be easily seen that, the trace of the density matrix is preserved in time i.e. ${\displaystyle \frac{d}{dt}}Tr[\widehat{\rho }(t)]=0$. With this we end the discussion about the optical Bloch equations derived using the semiclassical Rabi Model which lacks the inclusion of the Spontaneous emission effects which will be soon incorporated in the Quantum Treatment.

3.3. Quantum Mechanical Treatment of the Interaction of Two Level Atom with External Radiation Field

So far we have assumed that the Radiation field to be clasical but now we consider the radiation field to be quantized with the Electric field for a multi mode quantized radiation field is given by,

$$\widehat{\overrightarrow{E}}(\overrightarrow{r},t)=i\sum _m\sum _{s=\pm }\left({\displaystyle \frac{{\omega }_m}{2{ϵ}_{0}V}}\right)\left[{\widehat{a}}_{ms}(t)e^{i{\overrightarrow{k}}_m·\overrightarrow{r}}-{\widehat{a}}_{ms}^{†}{e}^{-i{\overrightarrow{k}}_m·\overrightarrow{r}}\right]{\widehat{e}}_{ms}$$

(85)

with, ${\widehat{e}}_{ms}$ denotes the polarisation vector i.e. the direction of the incident electric field and the bosonic creation and annihilation operators ${\widehat{a}}_{ms}^{†}$,${\widehat{a}}_{ms}$ satisfying the standard commutation relations given by,

$$\left[{\widehat{a}}_{ms},{\widehat{a}}_{n{s}^{\prime }}^{†}\right]={\delta }_{mn}{\delta }_{s{s}^{\prime }},[{\widehat{a}}_{ms},{\widehat{a}}_{ns}]=0,[{\widehat{a}}_{ms}^{†},{\widehat{a}}_{n{s}^{\prime }}^{†}]=0.$$

(86)

With the hamiltonian of the free Electromagnetic field is given by,

$${\widehat{H}}_f=\sum _m\sum _{s=\pm }{\omega }_m{\widehat{a}}_{ms}^{†}{\widehat{a}}_{ms}$$

(87)

Here, V being the volume of the Cavity within which the electromagnetic field has been quantized such that the Electric field given above is written in the Heisenberg picture and the electric field is written with a sum over all possible modes of vibration and mode of polarisations which is denoted by $s=\pm$ generally termed as the right and left handed state of polarisations. The evolution of the operators ${\widehat{a}}_{ms}(t)$ and ${\widehat{a}}_{ms}^{†}(t)$ are given by,

$${\widehat{a}}_{ms}(t)={\widehat{a}}_{ms}{e}^{-i{\omega }_{m}t},{\widehat{a}}_{ms}^{†}(t)={\widehat{a}}_{ms}^{†}{e}^{i{\omega }_{m}t}.$$

(88)

Now let us consider the interaction of the two level atom with the quantized radiation field such that in the long wavelength approximation the spatial dependence of the incident electric field can be neglected assuming that the wavelength of the incident electromagnetic field is very large compared to the spatial extent of the two level atom i.e. a length scale which is of the order of the bohr radius such that the total hamiltonian of the Atomic system coupled to the radiation field in the schrodinger picture along with Rotating wave approximation is given by,

$$\widehat{H}={\widehat{H}}_s+{\widehat{H}}_B+{\widehat{H}}_{SB}$$

(89)

Where,

$$\begin{array}{c}{\widehat{H}}_s={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z\end{array}$$

(90)

$$\begin{array}{c}{\widehat{H}}_B=\sum _m\sum _{s=\pm }{\widehat{a}}_{ms}^{†}{\widehat{a}}_{ms}\end{array}$$

(91)

$$\begin{array}{c}{\widehat{H}}_{SB}=\sum _m\sum _{s=\pm }\left(g_{ms}{\widehat{\sigma }}^{+}{\widehat{a}}_{ms}+g_{ms}^{\ast }{\widehat{a}}_{ms}^{†}{\widehat{\sigma }}^{+}\right)\end{array}$$

(92)

Where we have identified that,

$$g_{ms}=-i\left({\displaystyle \frac{{\omega }_m}{2{ϵ}_{0}V}}\right){\overrightarrow{\mu }}_{eg}·{\widehat{e}}_{ms}{e}^{i{\overrightarrow{k}}_m·{\overrightarrow{r}}_A}$$

(93)

Just to reiterate the free evolution of the operators ${\widehat{\sigma }}^{+}(t)$ and ${\widehat{\sigma }}^{-}(t)$ are given by,

$${\widehat{\sigma }}^{+}(t)={\widehat{\sigma }}^{+}{e}^{i{\omega }_{0}t},{\widehat{\sigma }}^{-}(t)={\widehat{\sigma }}^{-}{e}^{-i{\omega }_{0}t}$$

(94)

We have used the rotating wave approximation to write the interaction hamiltonian along with dipolar approximation which is elaborated below. We will use the interaction picture here to study the dynamics of the system and any arbitrary operator in the interaction picture is defined by an unitary transformation with the free part of the total hamiltonian such that,

$${\widehat{A}}^I(t)=e^{i({\widehat{H}}_s+{\widehat{H}}_B)t}{\widehat{A}}^S{e}^{-i({\widehat{H}}_s+{\widehat{H}}_B)t}$$

(95)

with, ${\widehat{A}}^S$ denotes the corresponding operator in the schrodinger picture. In the long wave approximation we can neglect the spatial variation of the the Electric field over the atomic length scale such that with $\widehat{\overrightarrow{E}}(\overrightarrow{r},t)\approx \widehat{\overrightarrow{E}}(\overrightarrow{r_A},t)$ and the Interaction hamiltonian in the schrodinger picture will be,

$${\widehat{H}}_{SB}=-i\sum _m\sum _{s=\pm }\left[{\overrightarrow{\mu }}_{eg}{\widehat{\sigma }}^{+}+{\overrightarrow{\mu }}_{ge}{\widehat{\sigma }}^{-}\right]·\left({\displaystyle \frac{{\omega }_m}{2{ϵ}_{0}V}}\right)\left[{\widehat{a}}_{ms}(t)e^{i{\overrightarrow{k}}_m·{\overrightarrow{r}}_A}-{\widehat{a}}_{ms}^{†}{e}^{-i{\overrightarrow{k}}_m·{\overrightarrow{r}}_A}\right]{\widehat{e}}_{ms}$$

(96)

Where we have identified $|e\rangle \langle g|={\widehat{\sigma }}^{+}$ and $|g\rangle \langle e|=\widehat{{\sigma }^{-}}$. Now when the Interaction hamiltonian is written in the interaction picture there will be four different types of terms whose time evolutions are summarised as follows,

$$\begin{array}{c}{\widehat{a}}_{ms}(t){\widehat{\sigma }}^{-}(t)\sim {\widehat{a}}_{ms}{\widehat{\sigma }}^{-}{e}^{-i({\omega }_m+{\omega }_0)t}\\ {\widehat{a}}_{ms}(t){\widehat{\sigma }}^{+}(t)\sim {\widehat{a}}_{ms}{\widehat{\sigma }}^{+}{e}^{i({\omega }_0-{\omega }_m)t}\\ {\widehat{a}}_{ms}^{†}(t){\widehat{\sigma }}^{+}(t)\sim {\widehat{a}}_{ms}^{†}{\widehat{\sigma }}^{+}{e}^{i({\omega }_m+{\omega }_0)t}\\ {\widehat{a}}_{ms}^{†}(t){\widehat{\sigma }}^{-}(t)\sim {\widehat{a}}_{ms}^{†}{\widehat{\sigma }}^{-}{e}^{i({\omega }_m-{\omega }_0)t}\end{array}$$

(97)

Now we will drop the rapidly oscillating (antiresonant) terms which goes as $e^{\pm i\left({\omega }_m+{\omega }_0\right)t}$ in comparison to the terms $e^{\pm i\left({\omega }_m-{\omega }_0\right)t}$ near the point of resonance i.e.${\omega }_m\approx {\omega }_0$ such that for small detuning we can write the interaction hamiltonian in the dipolar and rotating wave approximation keeping only two possible combinations not containing the anti resonant terms given by,${\widehat{a}}_{ms}{\widehat{\sigma }}^{+}$ and ${\widehat{a}}_{ms}^{†}{\widehat{\sigma }}^{-}$, dropping he other two terms of the type ${\widehat{a}}_{ms}{\widehat{\sigma }}^{-}$ and its hermitian conjugate with the system bath interaction parameter $g_{ms}$ defined above. For large detuning the Rotating wave approximation will not be valid anymore. But we will not discuss the large detuning case here for the time being. Now, along with the assumption that the bath is thermal in nature with the density operator at time $t=0$ is described by the density operator in the cannonical ensemble i.e. ${\widehat{\rho }}_B(0)={\displaystyle \frac{e^{-\beta {\widehat{H}}_B}}{Tr[e^{-\beta {\widehat{H}}_B}]}}$ with the denominator being the cannonical partition function given by, $Z=Tr\left[e^{-\beta {\widehat{H}}_B}\right]$, this assumption immediately follows that the Bath Hamiltonian will commute with ${\widehat{\rho }}_B(0)$. So we can write,

$$[{\widehat{\rho }}_B(0),{\widehat{H}}_B]=0,Tr[{\widehat{\rho }}_B(0){\widehat{a}}_{ms}]=0=Tr[{\widehat{\rho }}_B(0){\widehat{a}}_{ms}^{†}].$$

(98)

We can also obtain the following Bath Correlation functions which appear in the master equation given by,

$$\begin{array}{c}Tr\left[{\widehat{\rho }}_B(0){\widehat{a}}_{ms}{\widehat{a}}_{n{s}^{\prime }}\right]=0\end{array}$$

(99)

$$\begin{array}{c}Tr\left[{\widehat{\rho }}_B(0){\widehat{a}}_{ms}^{†}{\widehat{a}}_{n{s}^{\prime }}^{†}\right]=0\end{array}$$

(100)

$$\begin{array}{c}Tr\left[{\widehat{\rho }}_B(0){\widehat{a}}_{ms}{\widehat{a}}_{n{s}^{\prime }}^{†}\right]=\left(1+\overline{N}({\omega }_m,T)\right){\delta }_{mn}{\delta }_{s{s}^{\prime }}\end{array}$$

(101)

$$\begin{array}{c}Tr\left[{\widehat{\rho }}_B(0){\widehat{a}}_{ms}^{†}{\widehat{a}}_{n{s}^{\prime }}\right]=\overline{N}({\omega }_m,T){\delta }_{mn}{\delta }_{s{s}^{\prime }}\end{array}$$

(102)

Where, $\overline{N}({\omega }_m,T)$ defines the avearge number of excitations present in the m th mode of vibration at absolute temperature T given by the Bose einstein distribution with zero chemical potential written as,

$$\overline{N}({\omega }_m,T)={\displaystyle \frac{1}{e^{\beta {\omega }_m}-1}}$$

(103)

The dynamics of the system can be described by the Quantum Master Equation describing the time evolution of the Density operator of the system or sometimes we call the reduced density operator of the system obtained by taking the partial trace of the global density operator of the supersystem with respect to the bath states. We can write therefore that, ${\widehat{\rho }}_s^I(t)=T{r}_B[{\widehat{\rho }}_{tot}^I(t)]$. The time evolution of th reduced density operator of the syetem can be described by the Redfield Quantum Master equation obtained with the invocation of the Born Markov Approximation with the later describing the memory less time evolution of the Quantum system given for general scenario as follows,

$${\displaystyle \frac{d{\widehat{\rho }}_s^I(t)}{dt}}=-{\int }_0^{t}T{r}_B\left[\left[{\widehat{\rho }}_s^I(t)\otimes {\widehat{\rho }}_B(0),{\widehat{H}}_{SB}(t^{\prime })\right],{\widehat{H}}_{SB}(t)\right]d{t}^{\prime }$$

(104)

along with the initial condition ${\widehat{\rho }}_{tot}(0)={\widehat{\rho }}_s(0)\otimes {\widehat{\rho }}_B(0)$ which means that there was no correlation between the system and the bath at the initial time. Without loss of generlity it is in principle possible to set the upper limit of the time integral in the Redfield QME to infinity keping in mind that the bath correlation functions diminishes to zero after a certain time known as the bath correlatiion time denoted by $\tau$ such that,

$$Tr[{\widehat{\rho }}_B(0){\widehat{a}}_{ms}(t^{\prime }-t){\widehat{a}}_{ms}^{†}(0)]\approx 0\mathrm{for}\tau =(t-t^{\prime })>>{\tau }_B$$

(105)

Now in the present scenario with the interaction Hamiltonian obtained earlier, substituting it back in the master equation and expanding the double commutator terms and performing the partial trace operations and the integrals along with using the results mentioned before we can obtain the master equation describing the time evolution for the reduced density operator of the system in the interaction picture and later we can revert it back to the Schrodinger picture, neglecting the Lamb and Sterk shift terms to finally obtain the closed form given below.

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i[{\widehat{\rho }}_s(t),{\widehat{H}}_s]+{\displaystyle \frac{\gamma \overline{n}}{2}}\widehat{\mathcal{L}}[{\widehat{\sigma }}^{+}]{\widehat{\rho }}_s(t)+{\displaystyle \frac{\gamma (\overline{n}+1)}{2}}\widehat{\mathcal{L}}[\widehat{{\sigma }^{-}}]{\widehat{\rho }}_s(t)$$

(106)

With, $\overline{N}$ being the average number of excitation calculated at the frequency of the atomic transition i.e. ${\omega }_0$ such that,

$$\overline{n}=\overline{n}({\omega }_0,T)={\displaystyle \frac{1}{e^{\beta {\omega }_0}-1}}$$

(107)

In the above equation $\gamma$ is called the spontaneous emission coefficient which is given by,

$$\gamma =2\pi \sum _{s=\pm }\int d^3\overrightarrow{k}{\left|g(\overrightarrow{k},s)\right|}^2\delta (ck-{\omega }_0)D(\overrightarrow{k})$$

(108)

the Lindblad superoperators are defined in general as,

$$\widehat{\mathcal{L}}[\widehat{c}]{\widehat{\rho }}_s(t)=2\widehat{c}\widehat{\rho }{\widehat{c}}^{†}-\left\{{\widehat{c}}^{†}\widehat{c},{\widehat{\rho }}_s(t)\right\}$$

(109)

To summarise we can say that in order to cast the master equation in the standard Lindblad form we have neglected the Lamb and Sterk shift terms while evaluating the integrals.If the thermal effects are neglected i.e. if we assume zero temperature bath then the master equation will become,

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i[{\widehat{\rho }}_s(t),{\widehat{H}}_s]+{\displaystyle \frac{\gamma }{2}}\left[2{\widehat{\sigma }}^{-}{\widehat{\rho }}_s(t){\widehat{\sigma }}^{+}-{\widehat{\sigma }}^{+}{\widehat{\sigma }}^{-}{\widehat{\rho }}_s(t)-{\widehat{\rho }}_s(t){\widehat{\sigma }}^{+}{\widehat{\sigma }}^{-}\right]$$

(110)

From the above master equation we can write the equations describing the time evolution of the matrix element of the density operator of the two level atom as ,

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{ee}(t)}{dt}}=\left[\gamma \overline{n}{\rho }_{gg}(t)-\gamma (\overline{n}+1){\rho }_{ee}(t)\right]\end{array}$$

(111)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{gg}(t)}{dt}}=\left[-\gamma \overline{n}{\rho }_{gg}(t)+\gamma (\overline{n}+1){\rho }_{ee}(t)\right]\end{array}$$

(112)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{eg}(t)}{dt}}=-\left[i{\omega }_0+{\displaystyle \frac{\gamma (2\overline{n}+1)}{2}}(t)\right]{\rho }_{eg}(t)\end{array}$$

(113)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{ge}(t)}{dt}}=\left[i{\omega }_0-{\displaystyle \frac{\gamma (2\overline{n}+1)}{2}}\right]{\rho }_{ge}(t)\end{array}$$

(114)

The following equations can also be expressed in terms of the expectation values of ${\widehat{\sigma }}_z(t)$ , ${\widehat{\sigma }}^{+}(t)$ and ${\widehat{\sigma }}^{-}(t)$ respectively with,

$$\begin{array}{c}\langle {\widehat{\sigma }}_z(t)\rangle =Tr[{\widehat{\rho }}_s{\widehat{\sigma }}_z]={\rho }_{ee}(t)-{\rho }_{gg}(t)\\ {\rho }_{ee}(t)={\displaystyle \frac{1+\langle {\widehat{\sigma }}_z(t)\rangle }{2}}\\ {\rho }_{gg}(t)={\displaystyle \frac{1-\langle {\widehat{\sigma }}_z(t)\rangle }{2}}\\ \langle {\widehat{\sigma }}^{+}(t)\rangle =Tr[{\widehat{\rho }}_s(t){\widehat{\sigma }}^{+}]={\rho }_{ge}(t)\\ \langle {\widehat{\sigma }}^{-}(t)\rangle =Tr[{\widehat{\rho }}_s(t){\widehat{\sigma }}^{-}]={\rho }_{eg}(t)\end{array}$$

(115)

such that we have the following equations,

$$\begin{array}{c}{\displaystyle \frac{d}{dt}}\langle {\widehat{\sigma }}_z(t)\rangle =-\left[2\gamma \overline{n}\langle {\widehat{\sigma }}_z(t)\rangle +\gamma (1+\langle {\widehat{\sigma }}_z(t)\rangle )\right]\end{array}$$

(116)

$$\begin{array}{c}{\displaystyle \frac{d}{dt}}\langle {\widehat{\sigma }}^{+}\rangle (t)=\left[i{\omega }_0-{\displaystyle \frac{\gamma }{2}}(2\overline{n}+1)\right]\langle {\widehat{\sigma }}^{+}(t)\rangle \end{array}$$

(117)

$$\begin{array}{c}{\displaystyle \frac{d}{dt}}\langle {\widehat{\sigma }}^{-}(t)\rangle =-\left[i{\omega }_0+{\displaystyle \frac{\gamma }{2}}(2\overline{n}+1)\right]\langle {\widehat{\sigma }}^{-}(t)\rangle \end{array}$$

(118)

3.3.1. Full Fledged Quantum Treatment for the Optical Bloch Equation

The next task is to obtain the optical bloch equations [1] for the resonantly driven Two level atom which will incorporate the effect of the spontaneous emission. If the electromagnetic radiation field is described by coherent state then using a suitably chosen Unitary transformation we can seperate out the interaction hamiltonian in two parts with the first part describing the interaction of the Two level atom with the classical radiation field and the second term describing the effect of the quantum vacuum fluctuations. The idea is to describe one of the modes of the radiation field by the coherent state with the rest of the mods of vibrations are described by the vacuum state. Let us appreciate the step to obtain the hamiltonian for the Resonantly driven two level atom by a coherent radiation field which can be a Laser Field used for the cooling of atoms! Let us consider the interaction of the two level atom with a single mode quantized Electromagnetic field described by the hamiltonian written in the Dipolar and rotating wave approximation as follows with the Electric field operator of the beam is, $\widehat{E}=i\left({\displaystyle \frac{\omega }{2{ϵ}_{0}V}}\right)\left[\widehat{a}{e}^{ik{z}_a}-{\widehat{a}}^{†}{e}^{-ik{z}_a}\right]\widehat{e}$.

$$\widehat{H}={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z+\omega {\widehat{a}}^{†}\widehat{a}+(g{\widehat{\sigma }}^{+}\widehat{a}+g^{\ast }{\widehat{\sigma }}^{-}{\widehat{a}}^{†})$$

(119)

Here the system bath coupling hamiltonian i.e. the interaction hamiltonian of atom and field is given by,

$${\widehat{H}}_{SB}=(g{\widehat{\sigma }}^{+}\widehat{a}+g^{\ast }{\widehat{\sigma }}^{-}{\widehat{a}}^{†})$$

(120)

Here, we have assumed that the matrix element of the dipole moment dipole operator of the atom is in real and we have also assumed that the beam is propagating along the z direction and $\widehat{e}$ being the direction of the Electric field. In the expression of the Interaction Hamiltonian we have identified $g=-i({\overrightarrow{\mu }}_{eg}·\widehat{e})e^{ik{z}_a}$. Now let us consider the unitary transformation of the ${\widehat{H}}_{SB}$ defined as,

$${\widehat{H}}_{SB}^{\prime }=\widehat{D(\alpha }){\widehat{H}}_{SB}{\widehat{D}}^{†}(\alpha )$$

(121)

Where, $\widehat{D}(\alpha )=exp(\alpha {\widehat{a}}^{†}-\widehat{a}{\alpha }^{\ast })$, the operator is called the displacement operator and it is unitary in nature with the following properties,

$$\begin{array}{c}\widehat{D}(\alpha ){\widehat{D}}^{†}(\alpha )=\widehat{I}\end{array}$$

(122)

$$\begin{array}{c}{\widehat{D}}^{†}(\alpha )=\widehat{D}(-\alpha )\end{array}$$

(123)

$$\begin{array}{c}\widehat{D}(\alpha )=exp(-{\displaystyle \frac{1}{2}}{\left|\alpha \right|}^2)exp(\alpha {\widehat{a}}^{†})exp(-{\alpha }^{\ast }\widehat{a})\end{array}$$

(124)

$$\begin{array}{c}\widehat{D}(\alpha )\widehat{a}{\widehat{D}}^{†}(\alpha )=\widehat{D}(\alpha )\widehat{a}\widehat{D}(-\alpha )=\widehat{a}-\alpha \end{array}$$

(125)

$$\begin{array}{c}\widehat{D}(\alpha ){\widehat{a}}^{†}\widehat{D}(-\alpha )={\widehat{a}}^{†}-{\alpha }^{\ast }\end{array}$$

(126)

Those third, fourth and fifth properties are obtained using the Baker-Hausdorff-Campbell formulae. Using the following properties we can see that,

$${\widehat{H}}_{SB}^{\prime }={\widehat{H}}_{SB}+{\widehat{H}}_{SC}$$

(127)

with, ${\widehat{H}}_{SC}=-(g\alpha {\widehat{\sigma }}^{+}+g^{\ast }{\alpha }^{\ast }{\widehat{\sigma }}^{-})$.From the above interaction hamiltonian ${\widehat{H}}_{SB}^{\prime }$ it can be inferred that, when the Electromagnetic field is described by the coherent state then the unitary transformation of the Interaction hamiltonian can be written as the sum of two hamiltonians with ${\widehat{H}}_{SB}$ describing the interaction of the Atom with the vacuum of the radiation field and the other term ${\widehat{H}}_{SC}$ describing the interaction of the Two level atom with the Classical counterpart of the Radiation field. Now, we can identify the parameter $g\alpha$ as $-({\displaystyle \frac{{\omega }_R}{2}}{e}^{-i\omega t})$, with the Classical Electric field being $\overrightarrow{E}(t)=\widehat{e}cos\omega t$ and the classical interaction hamiltonian being written in the dipolar and Rotating Wave approximation such that, ${\widehat{H}}_{SC}=-\widehat{\overrightarrow{\mu }}·\widehat{e}cos\omega t\approx -{\displaystyle \frac{1}{2}}\left[({\overrightarrow{\mu }}_{eg}·\widehat{e}){\widehat{\sigma }}^{+}{e}^{-i\omega t}+h.c\right]$. We have defined the Rabi Frequency as, ${\omega }_R=-({\overrightarrow{\mu }}_{eg}·\widehat{e})$ and as per our assumption ${\omega }_R$ is real. So, after separating the interaction hamiltonian in two parts we we can redefine the system hamiltonian as,

$${\widehat{H}}_S={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z+{\displaystyle \frac{{\omega }_R}{2}}\left[{\widehat{\sigma }}^{+}{e}^{-i\omega t}+{\widehat{\sigma }}^{-}{e}^{i\omega t}\right]$$

(128)

with, the total hamiltonian of the supersystem being,

$$\widehat{H}={\widehat{H}}_S+{\widehat{H}}_B+{\widehat{H}}_{SB}$$

(129)

with, ${\widehat{H}}_B=\omega {\widehat{a}}^{†}\widehat{a}$ and, ${\widehat{H}}_{SB}=(g{\widehat{\sigma }}^{+}\widehat{a}+g^{\ast }{\widehat{\sigma }}^{-}{\widehat{a}}^{†})$ Now, its just a matter of generalization for the multimode case where again the radiation field is described by the coherent state with one of the mode being represented by the coherent state and rest of the modes are vacuum. For such a periodically driven two level atomic system can be described by the hamiltonian,

$$\begin{array}{c}\widehat{H}={\widehat{H}}_S+{\widehat{H}}_B+{\widehat{H}}_{SB}\end{array}$$

(130)

$$\begin{array}{c}{\widehat{H}}_S={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z+{\displaystyle \frac{{\omega }_R}{2}}\left[{\widehat{\sigma }}^{+}{e}^{-i\omega t}+{\widehat{\sigma }}^{-}{e}^{i\omega t}\right]\end{array}$$

(131)

$$\begin{array}{c}{\widehat{H}}_B=\sum _m\sum _{s=\pm }{\omega }_m{\widehat{a}}_{ms}^{†}{\widehat{a}}_{ms}\end{array}$$

(132)

$$\begin{array}{c}{\widehat{H}}_{SB}=\sum _m\sum _{s=\pm }\left(g_{ms}{\widehat{\sigma }}^{+}{\widehat{a}}_{ms}+g_{ms}^{\ast }{\widehat{\sigma }}^{-}{\widehat{a}}_{ms}^{†}\right)\end{array}$$

(133)

with, $g_{ms}=-i\left({\displaystyle \frac{{\omega }_m}{2{ϵ}_{0}V}}\right){\overrightarrow{\mu }}_{eg}·{\widehat{e}}_{ms}{e}^{i{\overrightarrow{k}}_m·{\overrightarrow{r}}_A}$ Similarly the quantum master equation describing the time evolution of the reduced density operator of the system i.e. ${\widehat{\rho }}_s(t)$ for the resonantly and periodically driven two level atom with, $\omega ={\omega }_0$ and with neglecting the thermal effects i.e. $\overline{n}\to 0$ will be given by,

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i[{\widehat{\rho }}_s(t),{\widehat{H}}_s]+{\displaystyle \frac{\gamma }{2}}\left[2{\widehat{\sigma }}^{-}{\widehat{\rho }}_s(t){\widehat{\sigma }}^{+}-{\widehat{\sigma }}^{+}{\widehat{\sigma }}^{-}{\widehat{\rho }}_s(t)-{\widehat{\rho }}_s(t){\widehat{\sigma }}^{+}{\widehat{\sigma }}^{-}\right]$$

(134)

Where, ${\widehat{H}}_S={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z+{\displaystyle \frac{{\omega }_R}{2}}\left[{\widehat{\sigma }}^{+}{e}^{-i\omega t}+{\widehat{\sigma }}^{-}{e}^{i\omega t}\right]$ Now similarly from this equation we can write,

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{ee}(t)}{dt}}=\left[{\displaystyle \frac{i{\omega }_R}{2}}{\rho }_{eg}(t)e^{i{\omega }_{0}t}-{\displaystyle \frac{i{\omega }_R}{2}}{\rho }_{ge}(t)e^{-i{\omega }_{0}t}\right]-\gamma {\rho }_{ee}(t)\end{array}$$

(135)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{gg}(t)}{dt}}=\left[-{\displaystyle \frac{i{\omega }_R}{2}}{\rho }_{eg}(t)e^{i{\omega }_{0}t}+{\displaystyle \frac{i{\omega }_R}{2}}{\rho }_{ge}(t)e^{-i{\omega }_{0}t}\right]+\gamma {\rho }_{ee}(t)\end{array}$$

(136)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{eg}(t)}{dt}}=-i{\omega }_0{\rho }_{eg}(t)-{\displaystyle \frac{\gamma }{2}}{\rho }_{eg}(t)+{\displaystyle \frac{i{\omega }_R}{2}}\left({\rho }_{ee}(t)-{\rho }_{gg}(t)\right)e^{-i{\omega }_{0}t}\end{array}$$

(137)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{ge}(t)}{dt}}=i{\omega }_0{\rho }_{ge}(t)-{\displaystyle \frac{\gamma }{2}}{\rho }_{ge}(t)-{\displaystyle \frac{i{\omega }_R}{2}}\left({\rho }_{ee}(t)-{\rho }_{gg}(t)\right)e^{i{\omega }_{0}t}\end{array}$$

(138)

Alternatively we can write the following equations in terms of $\langle {\widehat{\sigma }}_z(t)\rangle$,$\langle {\widehat{\sigma }}^{-}(t)\rangle$ and $\langle {\widehat{\sigma }}^{+}(t)\rangle$ such that,

$$\begin{array}{c}{\displaystyle \frac{d}{dt}}\langle {\widehat{\sigma }}^{+}(t)\rangle =\left(i{\omega }_0-{\displaystyle \frac{\gamma }{2}}\right)\langle {\widehat{\sigma }}^{+}(t)\rangle -{\displaystyle \frac{i{\omega }_R}{2}}\langle {\widehat{\sigma }}_z(t)\rangle e^{i{\omega }_{0}t}\end{array}$$

(139)

$$\begin{array}{c}{\displaystyle \frac{d}{dt}}\langle {\widehat{\sigma }}^{-}(t)\rangle =-\left(i{\omega }_0+{\displaystyle \frac{\gamma }{2}}\right)\langle {\widehat{\sigma }}^{+}(t)\rangle +{\displaystyle \frac{i{\omega }_R}{2}}\langle {\widehat{\sigma }}_z(t)\rangle e^{-i{\omega }_{0}t}\end{array}$$

(140)

$$\begin{array}{c}{\displaystyle \frac{d}{dt}}\langle {\widehat{\sigma }}_z(t)\rangle =-\gamma (1+\langle {\widehat{\sigma }}_z(t)\rangle )+i{\omega }_R\langle {\widehat{\sigma }}^{-}(t)\rangle e^{i{\omega }_{0}t}-i{\omega }_R\langle {\widehat{\sigma }}^{+}(t)\rangle e^{-i{\omega }_{0}t}\end{array}$$

(141)

Those equations are also called the Optical Bloch equations.

3.4. Solution of the Optical Bloch Equation

The above equations can be solved algebricaly in the rotating frame of reference. As we can see that for the periodically driven two level atom the System Hamiltonian is explicitly time independent with,

$$\begin{array}{c}\widehat{H}={\widehat{H}}_s(t)+{\widehat{H}}_B+{\widehat{H}}_{SB}\\ {\widehat{H}}_S(t)={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z+{\displaystyle \frac{{\omega }_R}{2}}\left({\widehat{\sigma }}^{+}{e}^{-i\omega t}+{\widehat{\sigma }}^{-}{e}^{+i\omega t}\right)\\ {\widehat{H}}_B=\sum _m\sum _{s=\pm }{\omega }_m{\widehat{a}}_{ms}^{†}{\widehat{a}}_{ms}\\ {\widehat{H}}_{SB}=\sum _m\sum _{s=\pm }\left(g_{ms}{\widehat{\sigma }}^{+}{\widehat{a}}_{ms}+g_{ms}^{\ast }{\widehat{\sigma }}^{-}{\widehat{a}}_{ms}^{†}\right)\end{array}$$

(142)

Now, we can rewrite the master equation for the reduced density operator of the system i.e. ${\widehat{\rho }}_s(t)$ after getting rid of the explicit time dependence of the system Hamiltonian by transforming to the Rotating frame of reference. The general prescription of transformation to the rotating frame of reference is described as follows. Let $\widehat{H}(t)$ being the hamiltonian which satisfies the time dependent schrodinger equation with state vector $|\psi (t)\rangle$ such that,

$$i{\displaystyle \frac{d}{dt}}|\psi (t)\rangle =\widehat{H}(t)|\psi (t)\rangle$$

(143)

Now, let us defined an unitary transformation $|\tilde{\psi }(t)\rangle =\widehat{U}(t)|\psi (t)\rangle$ with, ${\widehat{U}}^{†}(t)\widehat{U}(t)=\widehat{I}$ demanding that, $|\tilde{\psi }(t)\rangle$ satisfies the time dependent schrodinger equation with a transformed hamiltonian ${\widehat{H}}_{Rot}$ which is not time dependent such that,

$$i{\displaystyle \frac{d}{dt}}|\tilde{\psi }(t)\rangle ={\widehat{H}}_{Rot}|\tilde{\psi }(t)\rangle$$

(144)

Now, with the substitution, $|\tilde{\psi }(t)\rangle =\widehat{U}(t)|\psi (t)\rangle$ in the above equation we find the relation between the hamiltonian in the rotating frame and that of in the "lab" frame i.e. the relation between ${\widehat{H}}_{Rot}$ and $\widehat{H}(t)$ such that,

$${\widehat{H}}_{Rot}=-\widehat{U}(t){\displaystyle \frac{d{\widehat{U}}^{†}(t)}{dt}}+\widehat{U}(t)\widehat{H}(t){\widehat{U}}^{†}(t)$$

(145)

Now, we can use the following procedure to obtain the time independent form of the system hamiltonian in the rotating frame of reference where it will be time independent such that with, $\widehat{H}(t)={\widehat{H}}_S(t)$ and $\widehat{U}(t)=e^{i\omega {\widehat{\sigma }}_{z}t/2}$ we find that,

$${\widehat{H}}_S^{Rot}\equiv {\widehat{H}}_S=-{\displaystyle \frac{\Delta }{2}}{\widehat{\sigma }}_z+{\displaystyle \frac{{\omega }_R}{2}}[{\widehat{\sigma }}^{+}+{\widehat{\sigma }}^{-}]=-{\displaystyle \frac{\Delta }{2}}{\widehat{\sigma }}_z+{\displaystyle \frac{{\omega }_R}{2}}{\widehat{\sigma }}_x$$

(146)

Where we have defined the detuning frequency as $\Delta =\omega -{\omega }_0$. Now we can once again write the Master equation for ${\widehat{\rho }}_s(t)$ with this redefined time independent system hamiltonian as follows,

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i[{\widehat{\rho }}_s(t),{\widehat{H}}_S]+{\displaystyle \frac{\gamma }{2}}\left(2{\widehat{\sigma }}^{-}{\widehat{\rho }}_s(t){\widehat{\sigma }}^{+}-{\widehat{\sigma }}^{+}{\widehat{\sigma }}^{-}{\widehat{\rho }}_s(t)-{\widehat{\rho }}_s(t){\widehat{\sigma }}^{+}{\widehat{\sigma }}^{-}\right)$$

(147)

With, ${\widehat{H}}_S=-{\displaystyle \frac{\Delta }{2}}{\widehat{\sigma }}_z+{\displaystyle \frac{{\omega }_R}{2}}\left({\widehat{\sigma }}^{+}+{\widehat{\sigma }}^{-}\right)$. Now, writting the differential equations for the matrix element of the reduced density opertor of the system we get,

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{ee}(t)}{dt}}=-\gamma {\rho }_{ee}(t)+{\displaystyle \frac{i{\omega }_R}{2}}[{\rho }_{eg}(t)-{\rho }_{ge}(t)]\end{array}$$

(148)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{gg}(t)}{dt}}=\gamma {\rho }_e(t)-{\displaystyle \frac{i{\omega }_R}{2}}[{\rho }_{eg}(t)-{\rho }_{ge}(t)]\end{array}$$

(149)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{eg}(t)}{dt}}=i\Delta {\rho }_{eg}(t)-{\displaystyle \frac{\gamma }{2}}{\rho }_{eg}(t)+{\displaystyle \frac{i{\omega }_R}{2}}\left[{\rho }_{ee}(t)-{\rho }_{gg}(t)\right]\end{array}$$

(150)

$$\begin{array}{c}{\displaystyle \frac{d{\rho }_{ge}(t)}{dt}}=-\Delta {\rho }_{eg}(t)-{\displaystyle \frac{\gamma }{2}}{\rho }_{ge}(t)-{\displaystyle \frac{i{\omega }_R}{2}}\left[{\rho }_{ee}(t)-{\rho }_{gg}(t)\right]\end{array}$$

(151)

It is important to note that the System Hamiltonian in the rotating frame can be visualized as a $2\times 2$ matrix given by,

$${\widehat{H}}_S={\displaystyle \frac{1}{2}}\left(\begin{array}{cc}-\Delta & {\omega }_R\\ {\omega }_R& \Delta \end{array}\right)$$

(152)

Now, let us denote the bloch vector components by $u,v,w$ such that,

$$\begin{array}{c}u(t)=Tr[{\widehat{\rho }}_s(t){\widehat{\sigma }}_x]=\left[{\rho }_{eg}(t)+{\rho }_{ge}(t)\right]\end{array}$$

(153)

$$\begin{array}{c}v(t)=Tr[{\widehat{\rho }}_s(t){\widehat{\sigma }}_y]=i\left[{\rho }_{eg}(t)-{\rho }_{ge}(t)\right]\end{array}$$

(154)

$$\begin{array}{c}w(t)=Tr[{\widehat{\rho }}_s(t){\widehat{\sigma }}_z]=\left[{\rho }_{ee}(t)-{\rho }_{gg}(t)\right]\end{array}$$

(155)

Such that we can write the following equations,

$$\begin{array}{c}{\displaystyle \frac{du(t)}{dt}}=\Delta v(t)-{\displaystyle \frac{\gamma }{2}}u(t)\end{array}$$

(156)

$$\begin{array}{c}{\displaystyle \frac{dv(t)}{dt}}={\displaystyle \frac{\gamma v(t)}{2}}+\Delta u(t)+{\omega }_{R}v(t)\end{array}$$

(157)

$$\begin{array}{c}{\displaystyle \frac{dw(t)}{dt}}={\omega }_{R}v(t)-\gamma (1+w(t))\end{array}$$

(158)

Now, we can solve this equations in the steady state i.e. at $t\to \infty$ for which, $\underset{t\to \infty }{lim}{\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=0$, which means at steady state,

$$\begin{array}{c}{\displaystyle \frac{du(t)}{dt}}=0\Rightarrow \Delta v^{ss}={\displaystyle \frac{\gamma }{2}}{u}^{ss}\end{array}$$

(159)

$$\begin{array}{c}{\displaystyle \frac{dv(t)}{dt}}=0\Rightarrow {\displaystyle \frac{\gamma v^{ss}}{2}}+\Delta u^{ss}+{\omega }_R{w}^{ss}=0\end{array}$$

(160)

$$\begin{array}{c}{\omega }_R{v}^{ss}=\gamma (1+w^{ss})\end{array}$$

(161)

with, $u^{ss},v^s,w^{ss}$ refers to the steady state values of the Bloch vector components i.e. at $t\to \infty$. We are particularly interested to find out the steady state population of the excited state which will be, ${\rho }_{ee}(t\to \infty )={\displaystyle \frac{1+w^{ss}}{2}}$ which will be obtained by solving the system of three linear equations mentioned above. After solving them we find,

$$\begin{array}{c}{v}^{ss}={\displaystyle \frac{2{\omega }_R\gamma }{\left[4{\Delta }^2+{\gamma }^2+2{\omega }_R^2\right]}}\end{array}$$

(162)

$$\begin{array}{c}{u}^{ss}={\displaystyle \frac{4\Delta {\omega }_R}{\left[4{\Delta }^2+{\gamma }^2+2{\omega }_R^2\right]}}\end{array}$$

(163)

$$\begin{array}{c}{w}^{ss}=-\left[{\displaystyle \frac{{\gamma }^2+4{\Delta }^2}{4{\Delta }^2+{\gamma }^2+2{\omega }_R^2}}\right]\end{array}$$

(164)

Then from the above solutions we can directly write that,

$${\rho }_{ee}(t\to \infty )={\displaystyle \frac{1+w^{ss}}{2}}={\displaystyle \frac{{\omega }_R^2}{4{\Delta }^2+{\gamma }^2+2{\omega }_R^2}}={\displaystyle \frac{\frac{{\omega }_R^2}{4}}{{\Delta }^2+{\left(\frac{\gamma }{2}\right)}^2+\frac{{\omega }_R^2}{2}}}$$

(165)

us define the steady state parameter S such that,

$$S={\displaystyle \frac{\frac{{\omega }_R^2}{2}}{{\Delta }^2+{\left(\frac{\gamma }{2}\right)}^2}}$$

(166)

Such that, we can express the steady state Bloch vector components in terms of the parameter S with,

$$\begin{array}{c}{w}^{ss}=-{\displaystyle \frac{1}{1+S}}\end{array}$$

(167)

$$\begin{array}{c}{v}^{ss}=\left({\displaystyle \frac{\gamma }{{\omega }_R}}\right){\displaystyle \frac{S}{1+S}}\end{array}$$

(168)

$$\begin{array}{c}{u}^{ss}=\left({\displaystyle \frac{2\Delta }{{\omega }_R}}\right){\displaystyle \frac{S}{1+S}}\end{array}$$

(169)

We define the scattering rate as $R_{sc}=\gamma {\rho }_{ee}(t\to \infty )$.

4. Generalization of the Model for an Ensemble of Two Level Atoms

far in the previous section we have discussed the interaction of a single two level atom with the quantized radiation field in the dipolar and rotating wave approximation and now we are in a position to generalize the formulation for the ensemble of large number of identical two level atoms [2] which are assumed to be non-interacting in nature and for this situation the total Hamiltonian of the ensemble of two level atoms coupled with the quantized radiation field can be written as,

$$\begin{array}{c}\widehat{H}={\widehat{H}}_S+{\widehat{H}}_B+{\widehat{H}}_{SB}\end{array}$$

(170)

$$\begin{array}{c}{\widehat{H}}_S={\displaystyle \frac{{\omega }_0}{2}}\sum _m{\widehat{\sigma }}_m^z\end{array}$$

(171)

$$\begin{array}{c}{\widehat{H}}_B=\sum _m\sum _{s=\pm }{\omega }_m{\widehat{a}}_{ms}^{†}{\widehat{a}}_{ms}\end{array}$$

(172)

$$\begin{array}{c}{\widehat{H}}_{SB}=\sum _a\sum _m\sum _{s=\pm }\left(g_{ms}{\widehat{\sigma }}_a^{+}{\widehat{a}}_{ms}+g_{ms}^{\ast }{\widehat{\sigma }}_a^{-}{\widehat{a}}_{ms}^{†}\right)\end{array}$$

(173)

We have defined the interaction hamiltonian in the dipolar and rotating wave approximation with,${\widehat{H}}_{SB}=-{\sum }_a{\widehat{\overrightarrow{\mu }}}_a·\widehat{\overrightarrow{E}}({\overrightarrow{r}}_{cl})$, assuming that the electric field is uniform in the region of atomic confinement such that the spatial variation can be neglected. and the Electric field is given by,

$$\widehat{\overrightarrow{E}}(\overrightarrow{r_{cl}})=i\sum _m\sum _{s=\pm }\left({\displaystyle \frac{{\omega }_m}{2{ϵ}_{0}V}}\right)\left[{\widehat{a}}_{ms}{e}^{i{\overrightarrow{k}}_m·{\overrightarrow{r}}_{cl}}-{\widehat{a}}_{ms}^{†}{e}^{-i{\overrightarrow{k}}_m·r_{cl}}\right]{\widehat{e}}_{ms}$$

(174)

Similarly we can write,${\widehat{\overrightarrow{\mu }}}_a={\overrightarrow{\mu }}_{eg}^a{\widehat{\sigma }}_a^{+}+{\overrightarrow{\mu }}_{eg}^{a\ast }{\widehat{\sigma }}_a^{-}$ but as the atoms are identical we can say that the matrix element of the dipole moment operators for every atoms will be identical such that,

$${\overrightarrow{\mu }}_{eg}^a={}^a\langle e|\widehat{\overrightarrow{\mu }}|g\rangle ^a$$

Such that, we can identify the quantity,

$$g_{ms}=-i\left({\displaystyle \frac{{\omega }_m}{2{ϵ}_{0}V}}\right){\overrightarrow{\mu }}_{eg}·{\widehat{e}}_{ms}{e}^{i{\overrightarrow{k}}_m·{\overrightarrow{r}}_{cl}}$$

The most general way to represent the interaction hamiltonian for an ensemble of two level atoms coupled to the multi-mode quantized radiation field is given by,

$${\widehat{H}}_{SB}=\sum _a\sum _m\sum _{s=\pm }\left({\kappa }_{ams}{\widehat{\sigma }}_a^{+}{\widehat{a}}_{ms}+{\kappa }_{ams}^{\ast }{\widehat{\sigma }}_a^{-}{\widehat{a}}_{ms}^{†}\right)$$

(175)

Where the coefficient ${\kappa }_{ams}$ is defined as,

$${\kappa }_{ams}=-i\left({\displaystyle \frac{{\omega }_m}{2{ϵ}_{0}V}}\right){\overrightarrow{\mu }}_{eg}^a·{\widehat{e}}_{ms}{e}^{i{\overrightarrow{k}}_m·r_a}$$

The above expression will be most general when the atoms are not identical to each other but in this scenario all the atoms are identical along with the assumption that we have neglected the spatial variation of the Electric field by assuming the fact that the electric field is uniform in the region where the atoms are confined. Now we can write the master equation for the reduced density operator of the system in the schrodinger picture neglecting the lamb and sterk [3] shift terms such that,

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i[{\widehat{\rho }}_s(t),{\widehat{H}}_S]+{\displaystyle \frac{\gamma \overline{n}}{2}}\widehat{\mathcal{L}}\left[\sum _a{\widehat{\sigma }}_a^{+}\right]{\widehat{\rho }}_s(t)+{\displaystyle \frac{\gamma (\overline{n}+1)}{2}}\widehat{\mathcal{L}}\left[\sum _a{\widehat{\sigma }}_a^{-}\right]{\widehat{\rho }}_s(t)$$

(176)

Where, $\overline{n}=\overline{n}({\omega }_0,T)$, the Bose Einstein distribution evaluated at the atomic frequency i.e. $\omega ={\omega }_0$ at absolute temperature T with chemical potential $\mu =0$ such that,

$$\overline{n}=\overline{n}({\omega }_0,T)={\displaystyle \frac{1}{e^{\beta {\omega }_0}-1}}$$

The factor $\gamma$ is called the spontaneous emission coefficient defined as,

$$\gamma =2\pi \sum _{s=\pm }\int d^3\overrightarrow{k}{\left|g(\overrightarrow{k},s)\right|}^{2}D(\overrightarrow{k})\delta (ck-{\omega }_0)$$

(177)

And the lindblad superoperators are given by,

$$\begin{array}{c}\widehat{\mathcal{L}}\left[\sum _a{\widehat{\sigma }}_a^{+}\right]\widehat{{\rho }_s(t)}=\sum _a\sum _b\left[2{\widehat{\sigma }}_a^{+}\widehat{{\rho }_s}(t){\widehat{\sigma }}_b^{-}-\left\{{\widehat{\sigma }}_a^{-}{\widehat{\sigma }}_b^{+}{\widehat{\rho }}_s(t)+{\widehat{\rho }}_s(t){\widehat{\sigma }}_a^{-}{\widehat{\sigma }}_b^{+}\right\}\right]\end{array}$$

(178)

$$\begin{array}{c}\widehat{\mathcal{L}}\left[\sum _a{\widehat{\sigma }}_a^{-}\right]{\widehat{\rho }}_s(t)=\sum _a\sum _b\left[2{\widehat{\sigma }}_a^{-}{\widehat{\rho }}_s(t){\widehat{\sigma }}_b^{+}-\left\{{\widehat{\sigma }}_a^{+}{\widehat{\sigma }}_b^{-}{\widehat{\rho }}_s(t)+{\widehat{\rho }}_s(t){\widehat{\sigma }}_a^{+}{\widehat{\sigma }}_b^{-}\right\}\right]\end{array}$$

(179)

If the atoms were non-identical then with the general form of the interaction hamiltonian [4,5] written before in dipolar an rotating wave approximation [3], the master equation [6,7] describing the collection of two levels atom in the rotating wave approximation will be given by,

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i[{\widehat{\rho }}_s(t),{\widehat{H}}_S]+\sum _a{\displaystyle \frac{{\gamma }_a{\overline{n}}_a}{2}}\widehat{\mathcal{L}}\left[{\widehat{\sigma }}_a^{+}\right]{\widehat{\rho }}_s(t)+\sum _a{\displaystyle \frac{{\gamma }_a({\overline{n}}_a+1)}{2}}\widehat{\mathcal{L}}\left[{\widehat{\sigma }}_a^{-}\right]{\widehat{\rho }}_s(t)$$

(180)

Along with,

$$\begin{array}{c}\widehat{H}={\widehat{H}}_S+{\widehat{H}}_B+{\widehat{H}}_{SB}\end{array}$$

(181)

$$\begin{array}{c}{\widehat{H}}_S=\sum _a{\displaystyle \frac{{\omega }_a}{2}}{\widehat{\sigma }}_a^z\end{array}$$

(182)

$$\begin{array}{c}{\widehat{H}}_B=\sum _m\sum _{s=\pm }{\omega }_m{\widehat{a}}_{ms}^{†}{\widehat{a}}_{ms}\end{array}$$

(183)

$$\begin{array}{c}{\widehat{H}}_{SB}=\sum _a\sum _m\sum _{s=\pm }\left(g_{ams}{\widehat{\sigma }}_a^{+}{\widehat{a}}_{ms}+g_{ams}^{\ast }{\widehat{\sigma }}_a^{-}{\widehat{a}}_{ms}^{†}\right)\end{array}$$

(184)

With, the general definition of the lindbladian [8] superoperator is given by,

$$\widehat{\mathcal{L}}[\widehat{c}]{\widehat{\rho }}_s(t)=2\widehat{c}{\widehat{\rho }}_s(t){\widehat{c}}^{†}-\left\{{\widehat{c}}^{†}\widehat{c},{\widehat{\rho }}_s(t)\right\}$$

As the atoms are not identical in nature there spontaneous emission coefficients will also be different with,

$${\gamma }_a=2\pi \sum _{s=\pm }\int d^3\overrightarrow{k}{\left|g_a(\overrightarrow{k},s)\right|}^{2}D(\overrightarrow{k})\delta (ck-{\omega }_a)$$

(185)

With, ${\omega }_a$ being the atomic transition frequency for the ath atom. Now, if we consider the Lamb and sterk shift hamiltonian as well then the master equation [6] for the single two level atom [9] interacting with the quantized radiation field will be given by,

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i\left[{\widehat{\rho }}_s(t),({\widehat{H}}_S+{\widehat{H}}_{LS})\right]+{\displaystyle \frac{\gamma \overline{n}}{2}}\widehat{\mathcal{L}}[{\widehat{\sigma }}^{+}]{\widehat{\rho }}_s(t)+{\displaystyle \frac{\gamma (\overline{n}+1)}{2}}\widehat{\mathcal{L}}[{\widehat{\sigma }}^{-}]{\widehat{\rho }}_s(t)$$

(186)

Where, the Lamb and stark shift is incorporated through ${\widehat{H}}_{LS}$ which is traditionally called the Lamb and Stark shift hamiltonian defined as,

$${\widehat{H}}_{LS}=\left(\Delta +{\Delta }^{\prime }\right){\widehat{\sigma }}^{+}{\widehat{\sigma }}^{-}-{\Delta }^{\prime }{\widehat{\sigma }}^{-}{\widehat{\sigma }}^{+}$$

(187)

where, $\Delta ,{\Delta }^{\prime }$ are defined as follows,

$$\begin{array}{c}\Delta =\sum _{s=\pm }\mathcal{P}\int d^3\overrightarrow{k}{\displaystyle \frac{D(\overrightarrow{k})|g(\overrightarrow{k},s){|}^2}{{\omega }_0-ck}}\end{array}$$

(188)

$$\begin{array}{c}{\Delta }^{\prime }=\sum _{s=\pm }\mathcal{P}\int d^3\overrightarrow{k}{\displaystyle \frac{D(\overrightarrow{k})|g(\overrightarrow{k},s){|}^2\overline{n}(ck,T)}{{\omega }_0-ck}}\end{array}$$

(189)

Here, $\mathcal{P}N$ stands for the Cauchy principle value. ow, including the Lamb and stark shift hamiltonian and with a little bit of algebra the final form of the master equation will become,

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i\left[{\widehat{\rho }}_s(t),{\displaystyle \frac{{\omega }_0^{\prime }}{2}}{\widehat{\sigma }}_z\right]+{\displaystyle \frac{\gamma \overline{n}}{2}}\widehat{\mathcal{L}}\left[{\widehat{\sigma }}^{+}\right]{\widehat{\rho }}_s(t)+{\displaystyle \frac{\gamma (\overline{n}+1)}{2}}\widehat{\mathcal{L}}\left[{\widehat{\sigma }}^{-}\right]{\widehat{\rho }}_s(t)$$

(190)

Where, the atomic transition frequency in the above master equation is being modified due to the interaction with the radiation field with, ${\omega }_0^{\prime }={\omega }_0+(\Delta +2{\Delta }^{\prime })$, with ${\omega }_0$ being the original atomic transition frequency in the absence of any interaction. Here, in order to obtain the above master equations starting from the Redfield equation or Born-Markov equation, we have to evaluate some integrals whose real and imaginary parts finally appear in the final form of the master equation casted in the Linblad form. the integrals are as follows,

$$\begin{array}{c}{\alpha }_1+i{\beta }_1=\sum _{s=\pm }\int {\int }_0^{\infty }{d}^3\overrightarrow{k}d\tau {\left|g(\overrightarrow{k},s)\right|}^{2}D(\overrightarrow{k})e^{i(ck-{\omega }_0)\tau }\left[1+\overline{n}(ck.T)\right]\end{array}$$

(191)

$$\begin{array}{c}{\alpha }_2+i{\beta }_2=\sum _{s=\pm }\int {\int }_0^{\infty }{d}^3\overrightarrow{k}d\tau {\left|g(\overrightarrow{k},s)\right|}^{2}D(\overrightarrow{k})e^{-i(ck-{\omega }_0)\tau }\overline{n}(ck.T)\end{array}$$

(192)

The following integrals has been evaluated with respect to $\tau$ by using the standard integral,

$${\int }_0^{\infty }{e}^{-i(ck-{\omega }_0)\tau }=\pi \delta (ck-{\omega }_0)+i\mathcal{P}\left({\displaystyle \frac{1}{{\omega }_0-ck}}\right)$$

And we found,

$${\alpha }_2={\displaystyle \frac{\gamma \overline{n}}{2}},{\alpha }_1={\displaystyle \frac{\gamma (\overline{n}+1)}{2}},{\beta }_1=-(\Delta +{\Delta }^{\prime }),{\beta }_2={\Delta }^{\prime }$$

. Which, leads to the final results, with the real parts of this integrals appears as the coefficients of the Lindblad super-operators ${\displaystyle \frac{\gamma \overline{n}}{2}}$ and ${\displaystyle \frac{\gamma (\overline{n}+1)}{2}}$ respectively and the imaginary parts contribute in the Lamb-Stark shift Hamiltonian ${\widehat{H}}_{LS}$. Now, this formulation can be extended for an ensemble of non-interacting two level atoms with the inclusion of the Lamb-Stark shift hamiltonian such that, the master equation becomes,

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i[{\widehat{\rho }}_s(t),\left({\widehat{H}}_S+{\widehat{H}}_{LS}\right)]+{\displaystyle \frac{\gamma \overline{n}}{2}}\widehat{\mathcal{L}}\left[\sum _a{\widehat{\sigma }}_a^{+}\right]{\widehat{\rho }}_s(t)+{\displaystyle \frac{\gamma (\overline{n}+1)}{2}}\widehat{\mathcal{L}}\left[\sum _a{\widehat{\sigma }}_a^{-}\right]{\widehat{\rho }}_s(t)$$

(193)

With, the Lamb-Stark shift Hamiltonian being,

$${\widehat{H}}_{LS}=\left(\Delta +{\Delta }^{\prime }\right)\sum _a\sum _b{\widehat{\sigma }}_a^{+}{\widehat{\sigma }}_b^{-}-{\Delta }^{\prime }\sum _a\sum _b{\widehat{\sigma }}_a^{-}{\widehat{\sigma }}_b^{+}$$

With, $\Delta ,{\Delta }^{\prime }$ has been defined above. If, the atoms were non-identical then with the Rotating wave approximation we get the exact same form of the master equation which is already mentioned before neglecting the ${\widehat{H}}_{LS}$. Now, including ${\widehat{H}}_{LS}$, we get

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i[{\widehat{\rho }}_s(t),\left({\widehat{H}}_S+{\widehat{H}}_{LS}\right)]+\sum _a{\displaystyle \frac{{\gamma }_a{\overline{n}}_a}{2}}\widehat{\mathcal{L}}\left[{\widehat{\sigma }}_a^{+}\right]{\widehat{\rho }}_s(t)+\sum _a{\displaystyle \frac{{\gamma }_a({\overline{n}}_a+1)}{2}}\widehat{\mathcal{L}}\left[{\widehat{\sigma }}_a^{-}\right]{\widehat{\rho }}_s(t)$$

(194)

With,

$$\begin{array}{c}{\overline{n}}_a={\overline{n}}_a({\omega }_a,T)={\displaystyle \frac{1}{e^{\beta {\omega }_a-1}}}\end{array}$$

(195)

$$\begin{array}{c}{\gamma }_a=2\pi \sum _{s=\pm }\int d^3\overrightarrow{k}\delta (ck-{\omega }_a){\left|g_a(\overrightarrow{k},s)\right|}^{2}D(\overrightarrow{k})\end{array}$$

(196)

$$\begin{array}{c}{\widehat{H}}_{LS}=\sum _a\left[\left({\Delta }_a+{\Delta }_a^{\prime }\right){\widehat{\sigma }}_a^{+}{\widehat{\sigma }}_a^{-}-{\Delta }_a^{\prime }{\widehat{\sigma }}_a^{-}{\widehat{\sigma }}_a^{+}\right]\end{array}$$

(197)

$$\begin{array}{c}{\Delta }_a=\sum _{s=\pm }\mathcal{P}\int d^3\overrightarrow{k}{\displaystyle \frac{D(\overrightarrow{k}){\left|g_a(\overrightarrow{k},s)\right|}^2}{{\omega }_a-ck}}\end{array}$$

(198)

$$\begin{array}{c}{\Delta }_a^{\prime }=\sum _{s=pm}\mathcal{P}\int d^3\overrightarrow{k}{\displaystyle \frac{D(\overrightarrow{k}){\left|g_a(\overrightarrow{k},s)\right|}^2{\overline{n}}_a(ck,T)}{{\omega }_a-ck}}\end{array}$$

(199)

Along with the ${\widehat{H}}_{LS}$ term and defining the modified atomic transition frequency due to the interaction with the radiation field the final form of the master equation in RWA reduces to,

$${\displaystyle \frac{d{\widehat{\rho }}_s(t)}{dt}}=i\sum _a[{\widehat{\rho }}_s(t),{\displaystyle \frac{{\omega }_a^{\prime }}{2}}{\widehat{\sigma }}_a^z]+\sum _a{\displaystyle \frac{{\gamma }_a{\overline{n}}_a}{2}}\widehat{\mathcal{L}}\left[{\widehat{\sigma }}_a^{+}\right]{\widehat{\rho }}_s(t)+\sum _a{\displaystyle \frac{{\gamma }_a({\overline{n}}_a+1)}{2}}\widehat{\mathcal{L}}\left[{\widehat{\sigma }}_a^{-}\right]{\widehat{\rho }}_s(t)$$

(200)

With,

$${\omega }_a^{\prime }={\omega }_a+\left({\Delta }_a+2{\Delta }_a^{\prime }\right)$$

Our primary focus will be the case where the atoms are non interacting and identical in nature such that, they are uniformly coupled with the radiation field which we have described earlier along with the approximation to neglect the spatial variation of the electric field, assumed to be uniform in the region of the atomic confinement of an optical trap like set up. In this case alternatively we can write the interaction hamiltonian as,

$${\widehat{H}}_{SB}=\sum _a\sum _m\left({\kappa }_m{\widehat{\sigma }}_a^{+}{\widehat{b}}_m+{\kappa }_m^{\ast }{\widehat{\sigma }}_a^{-}{\widehat{b}}_m^{†}\right)$$

With,

$$\sum _{s=\pm }{g}_{ms}{\widehat{a}}_{ms}={\kappa }_m{\widehat{b}}_m,\sum _{s=\pm }{g}_{ms}^{\ast }{\widehat{a}}_{ms}^{†}={\kappa }_m^{\ast }{\widehat{b}}_m^{†}$$

(201)

The preservation of the bosonic commutation [4,10] relation of the newly defined Hermitian conjugate pair ${\widehat{b}}_m,{\widehat{b}}_m^{†}$ we must have the following condition,

$$\begin{array}{c}\left[{\widehat{b}}_m,{\widehat{b}}_n\right]=0\mathrm{as}\left[{\widehat{a}}_{ms},{\widehat{a}}_{n{s}^{\prime }}\right]=0\end{array}$$

(202)

$$\begin{array}{c}\left[{\widehat{b}}_m^{†},{\widehat{b}}_n^{†}\right]=0\mathrm{as}\left[{\widehat{a}}_{ms}^{†},{\widehat{a}}_{n{s}^{\prime }}^{†}\right]=0\end{array}$$

(203)

$$\begin{array}{c}\left[{\widehat{b}}_m^{†},{\widehat{b}}_n^{†}\right]=\widehat{I}{\delta }_{mn}\Rightarrow {\kappa }_m{\kappa }_n^{\ast }{\delta }_{mn}=\sum _{s=\pm }{g}_{ms}{g}_{ns}^{\ast }{\delta }_{mn}\end{array}$$

(204)

$$\begin{array}{c}\left[{\widehat{b}}_m,{\widehat{b}}_m^{†}\right]=\widehat{I}\Rightarrow \sum _{s=\pm }|g_{ms}{|}^2={\left|{\kappa }_m\right|}^2\end{array}$$

(205)

$$\begin{array}{c}{\widehat{b}}_m(t)=e^{i{\widehat{H}}_{B}t}{\widehat{b}}_m{e}^{-i{\widehat{H}}_{B}t}={\widehat{b}}_m{e}^{-i{\omega }_{m}t}\end{array}$$

(206)

$$\begin{array}{c}{\widehat{b}}_m^{†}(t)=e^{i{\widehat{H}}_{B}t}{\widehat{b}}_m^{†}{e}^{-i{\widehat{H}}_{B}t}={\widehat{b}}_m^{†}{e}^{i{\omega }_{m}t}\end{array}$$

(207)

With this we end the discussion over the model describing the interaction of the single two level atom and the ensemble of non-interacting two level atoms with both the cases where they are identical and they are not identical coupled to the multi-mode electromagnetic radiation field using the semiclassical Rabi model and fully quantum mechanical description. It is worth to mention that the Jayne’s Cummings model and Tavies Cummings model can be seen as two special cases of the following discussions, where we consider the interaction of a single two level atom and the collection of identical two level atoms interacting with a single mode quantized radiation field respectively in the regime of Rotating Wave approximation. It is also possible to describe the interaction without making the rotating wave approximation for both single two level atom and the collection of large number of identical and non-interacting two level atoms which are popularly known as the Quantum Rabi Model and Dicke Model [11] of super-radiance respectively which goes beyond the RWA and there is another model which describes the interaction of the single two level atom with the single mode quantized radiation field keeping the rapidly oscillating anti-resonant terms in the interaction Hamiltonian and dropping the terms which are not rapidly oscillating for small de tuning just the opposite of the Jayne’s Cummings model known as the Anti-Jayne’s Cummings Model which will not be discussed in great details. The Hamiltonians of the above-mentioned models are given by,

$$\begin{array}{c}{\widehat{H}}_{AJC}={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z+\omega {\widehat{a}}^{†}\widehat{a}+\left(g{\widehat{\sigma }}^{+}\widehat{a}+g^{\ast }{\widehat{\sigma }}^{-}{\widehat{a}}^{†}\right)\end{array}$$

(208)

$$\begin{array}{c}{\widehat{H}}_{AJCM}={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z+\omega {\widehat{a}}^{†}\widehat{a}+(g{\widehat{\sigma }}^{-}\widehat{a}+g^{\ast }{\widehat{\sigma }}^{+}{\widehat{a}}^{†})\end{array}$$

(209)

$$\begin{array}{c}{\widehat{H}}_{QR}={\displaystyle \frac{{\omega }_0}{2}}{\widehat{\sigma }}_z+\omega {\widehat{a}}^{†}\widehat{a}+g\left({\widehat{\sigma }}^{+}+{\widehat{\sigma }}^{-}\right)\left(\widehat{a}+{\widehat{a}}^{†}\right)\end{array}$$

(210)

$$\begin{array}{c}{\widehat{H}}_{TCM}={\displaystyle \frac{{\omega }_0}{2}}\sum _m{\widehat{\sigma }}_m^z+\omega {\widehat{a}}^{†}\widehat{a}+\sum _m\left(g{\widehat{\sigma }}_m^{+}\widehat{a}+g^{\ast }{\widehat{\sigma }}_m^{-}{\widehat{a}}^{†}\right)\end{array}$$

(211)

$$\begin{array}{c}{\widehat{H}}_{DC}={\displaystyle \frac{{\omega }_0}{2}}\sum _m{\widehat{\sigma }}_m^z+\omega {\widehat{a}}^{†}\widehat{a}+{\displaystyle \frac{G}{\sqrt{N}}}\sum _m\left({\widehat{\sigma }}_m^{+}+{\widehat{\sigma }}_m^{-}\right)\left(\widehat{a}+{\widehat{a}}^{†}\right)\end{array}$$

(212)