TLS description of optical absorption

In this notebook we present a two level system (TLS) based description of optical absorption. We discuss the equation of motion of the system (EQM) in the rotating frame (RF) and we derive the Bloch vector representation of the model.

The Hamiltonian of the system is written as

\[ \begin{align}\begin{aligned} H_0 = -\frac{1}{2}\omega_0 \sigma_z \, ,\\where :math:`\sigma_z` is the Pauli matrix. :math:`H_0` is diagonal in the chosen representation with a energy shift between the two levels given by :math:`\omega_0`. The eigenstates are represented as\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned}\begin{split} |1\rangle = \left(\begin{array}{c} 1\\ 0\\ \end{array} \right) \, , \quad |2\rangle = \left(\begin{array}{c} 0\\ 1\\ \end{array} \right) \, .\end{split}\\Note that in this notation :math:`|1\rangle` describes the GS of the system, while :math:`|2\rangle` is the excited state.\end{aligned}\end{align} \]

The interaction of the optical pump with the system is described by an Hamiltonian \(H^I\) that has only non-diagonal non vanishing matrix elements.

The equation of motion (EQM) of the density matrix (DM) are expressed by the Liouville equation

\[ \begin{align}\begin{aligned} i \dot{\rho} = \left[H_0+H^I,\rho\right]\\and we can explicitly write the associated equations for the matrix elements of the DM in the chosen basis using the structure of both the Hamitonian defined above. We obtain:\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned}\begin{split} i \dot{\rho}_{11} = H^I_{12}\rho_{21} - H^I_{21}\rho_{12} \\ i \dot{\rho}_{22} = H^I_{21}\rho_{12} - H^I_{12}\rho_{21} \\ i \dot{\rho}_{12} = -\omega_0\rho_{12} + H^I_{12}(\rho_{22}-\rho_{11}) \\ i \dot{\rho}_{21} = \omega_0\rho_{21} - H^I_{21}(\rho_{22}-\rho_{11})\end{split}\\where :math:`H^I_{12} = \langle 1|H^I|2\rangle` and :math:`H^I_{21} = \langle 2|H^I|1\rangle`, furthermore :math:`H^I_{12} = H^{I*}_{21}` since :math:`H^I` is hermitian. These equations imply that :math:`\rho_{21}` is the complex conjugate of :math:`\rho_{12}` and that :math:`\rho_{11}+\rho_{22}` is constant, that is the conservation of the charge (in what follows the total charge is assumed to be equal to 1).\end{aligned}\end{align} \]

It is convenient to describe the system in terms of polarization and of the inversion variables defined as

\[ \begin{align}\begin{aligned} p = \rho_{12} \, , \quad I=\rho_{11}-\rho_{22} \, , \quad H_I = H^{I}_{12}\\In terms of these variables the EQM for the density matrix read\end{aligned}\end{align} \]

\[\begin{split}\dot{p} = i\omega_0 p +i H_I I \\ \dot{I} = -2i\left(H_Ip^*-H_I^*p\right)\end{split}\]

Interaction Hamiltonian and rotating wave approximation (RWA)

We consider an interaction between the system and the optical pump described by a dipole interaction, that is

\[ \begin{align}\begin{aligned} H^I = - \mathbf{\mu}\cdot\mathbf{E}(t)\\where :math:`\mathbf{\mu}` is the dipole operator. We choose a linear polarized electric field (say along the :math:`x` axis). The time dependece of the field is described by a monocromatic wave of energy :math:`\omega` times a real slow enevelope function :math:`E_0(t)`. Due to this assumption the matrix element :math:`H^I_{12}` can be expressed as\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} H^I_{12} = - \mu^{x}_{12} E_0(t) e^{-\frac{(t-t_0)^2}{2w^2}}sin(\omega t) \doteq - \Omega(t) sin(\omega t)\\where :math:`\mu^{x}_{12}` is matrix element of the component of the dipole parallel to :math:`\mathbf{E}` evaluated between the states 1 and 2 and we have introduced the (time-dependent) *Rabi frequency* :math:`\Omega(t)`:\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \Omega(t) = \Omega_0e^{-\frac{(t-t_0)^2}{2w^2}} \, , \, \mathrm{where} \,\, \Omega_0 = \mu^{x}_{12}E_0\\here :math:`\Omega_0` is the usual (constant) Rabi frequency (we have set :math:`\hbar=1` in this analysis). Note that since :math:`\mu^{x}_{12}` is a complex quantity so it is the Rabi frequency and we introduce its real and imaginary parts :math:`\Omega_R(t)` and :math:`\Omega_I(t)` (as well as the associated constant factors :math:`\Omega^0_R` and :math:`\Omega^0_I`) for the susequent analysis.\end{aligned}\end{align} \]

The RWA can be used if we are interested in probing the system with energy of the pump not very different from the energy gap of the system, that is if we introduce the energy shift \(\delta\) as

\[ \begin{align}\begin{aligned} \omega = \omega_0 + \delta\\where the condition :math:`\delta \ll \omega_0` is satisfied. This fact allows us to introduce the rotating wave approximation. To understand this point consider that, for instance, the generic solution for :math:`p` can be written as\end{aligned}\end{align} \]

\[p = p_0e^{i\omega_0t} -i e^{i\omega_0t}\int_{0}^{t}dt'e^{-i\omega_0t'}\Omega(t')sin(\omega t') I(t')\]

Due to the presence of the factor \(e^{-i\omega_0 t'}\) in the integral the sine can be splitted into complex exponentials and only the addend \(e^{i\omega t'}\) gives relevant contributions since the fast oscillating terms cancel. So the interaction matrix element in the RWA reads

\[ \begin{align}\begin{aligned} H_I = \frac{i}{2} \Omega(t)e^{i\omega t}\\Accordingly the EQM are expressed as\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned}\begin{split} \dot{p} = i\omega_0 p -\frac{1}{2}\Omega e^{i\omega t} I \\ \dot{I} = \Omega e^{i\omega t}p^* +\Omega^* e^{-i\omega t}p\end{split}\\Note that this argument is usually discussed and justified for **constant** Rabi frequency, however if the envelope function :math:`E_0(t)` is slow we can split the previous integral into various terms in which the frequency is constant and apply the same argument to all the addends.\end{aligned}\end{align} \]

If the RWA is adopted it is possible to discuss the physics of the system in the rotating frame (RF), where the fast oscillating terms are absent. The change of variables to the RF is

\[ \begin{align}\begin{aligned} p' = e^{-i\omega t}p \, , \quad I' = I\\In this frame the EQM are\end{aligned}\end{align} \]

\[\begin{split}\dot{p}' = -i\delta p' -\frac{1}{2}\Omega I' \\ \dot{I}' = \Omega p'^* +\Omega^* p'\end{split}\]

Bloch vector representation

The Bloch vector representation is a change of variables from the polarization and inversion to the 3-dimensional vector \(u\) defined as

\[ \begin{align}\begin{aligned} u_1 = 2Re(p) \, , \quad u_2 = -2Im(p) \, \quad I = u_3\\This formulation is particularly useful for computational reason since the components of the Bloch vector are real quantities and can be easily managed for numerical integration of the EQM.\end{aligned}\end{align} \]

The EQM in this basis read

\[\begin{split}\dot{u}_1 = -\delta u_2 - \Omega_R u_3\\ \dot{u}_2 = \delta u_1 + \Omega_I u_3 \\ \dot{u}_3 = \Omega_R u_1-\Omega_I u_2\end{split}\]

Inclusion of the coherence dephasing time

The polarization-inversion or equivalently the Bloch equation description allows us to add an explicit dephasing time for the coherence.

This is done by adding the term \(-\eta p\) in the EQM for the polarization. In the Bloch vector notation the complete equations read

\[\begin{split}\dot{u}_1 = -\delta u_2 - \Omega_R u_3 -\eta u_1\\ \dot{u}_2 = \delta u_1 + \Omega_I u_3 -\eta u_2\\ \dot{u}_3 = \Omega_R u_1-\Omega_I u_2\end{split}\]

Representation of the polarization and of the number of carriers

It is useful to express some observables like the polarization or the density of charge in the excited states in this formalism. These quantities can be compared with the results provided by the abinitio real time simulations.

The density of charge in the excited state is given by \(\rho_{22}\) that, using the condition \(\rho_{11}+\rho_{22} = 1\) can be expressed as

\[\rho_{11} = \frac{1-I}{2}\]

For what concerns the polarization, assuming that the only non vanishing matrix elements of the dipole operator (say \(x\) component) are the non diagonal ones we have

\[ \begin{align}\begin{aligned} P = p x_{12}^* + p^* x_{12}\\that can be expressed in the terms of the components of the Bloch vector as\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} P = u_1 Re(x_{21}) + u_2 Im(x_{21})\\So we observe that the elemnts :math:`u_1` and :math:`u_2` are responsible for the behavior of the polarization.\end{aligned}\end{align} \]

Complex rotation to real dipoles and Rabi coupling

It is possible to reabsorb the complex phase of the dipole \(d=d_{12}\), and consequently of the Rabi coupling, since \(\Omega_0 = Ad/\hbar\). This procedure preserves the value of the observable (the mean on the k point index is understood)

\[ \begin{align}\begin{aligned} \mu = Tr(\hat{\rho}\hat{d}) = pd^*+p^*d\\Indeed we can express the complex dipole as :math:`e^{i\alpha}|d|` and if we define complex rotated polarization as\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \hat{p} = e^{-i\alpha}p \, , \,\, \mathrm{where} \,\, \alpha = atan\left(\frac{Im(d)}{Re(d)}\right)\\The EQM in terms of these variables read\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned}\begin{split} \dot{\hat{p}} = -i\delta\hat{p} -\frac{1}{2}\Omega I \\ \dot{I} = \Omega\left(\hat{p} +\hat{p}^*\right)\end{split}\\where in these equations the envelope function :math:`\Omega(t)` is real and depends on the module of the Rabi coupling :math:`|\Omega_0|`. The observable polarization in this basis is expressed as\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \mu = |d|(\hat{p}+\hat{p}^*)\\So in this basis we deal with EQM with real parameters that provide equivalent results for the observable.\end{aligned}\end{align} \]

The equations for the Bloch vector in this basis read

\[ \begin{align}\begin{aligned}\begin{split} \dot{u}_1 = -\delta u_2 - \Omega u_3\\ \dot{u}_2 = \delta u_1 \\ \dot{u}_3 = \Omega u_1\end{split}\\while the observable polarization is computed as :math:`\mu=|d|u_1`.\end{aligned}\end{align} \]

Analytic solution for constant field (to be updated…)

We seek for the general solution of the EQM in the RF for a constant field.

We start from a guess for \(x\) of the form

\[ \begin{align}\begin{aligned} x(t) = Asin(\Delta t) + Bcos(\Delta t) + C \,\,\, , with \quad \Delta = \sqrt{\delta^2+|\Omega|^2}\\the constants :math:`A,B,C` will be determined imposing the boundary conditions. Plugging this guess in the EQM for :math:`x` we obtain\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} y(t) = \frac{B\Delta}{\delta}sin(\Delta t) - \frac{A\Delta}{\delta}cos(\Delta t)\\and pluggin this expression in the EQM for :math:`y` allows us to derive the expression of :math:`z`\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} z(t) = \frac{A\Omega_r}{\delta}sin(\Delta t) +\frac{B|\Omega|}{\delta}cos(\Delta t) -\frac{C\delta}{|\Omega|}\\pluggin this expression in the EQM for :math:`z` we obtain :math:`-|\Omega| y(t)` as expected, so the guess is correct.\end{aligned}\end{align} \]

The constants are determined through the boundary conditions

\[ \begin{align}\begin{aligned} x_0 = B+C \, , \quad y_0 = - \frac{A\Delta}{\delta} \, , \quad z_0 = \frac{B|\Omega|}{\delta} -\frac{C\delta}{|\Omega|}\\and their solution reads\end{aligned}\end{align} \]

\[A = -\frac{\delta}{\Delta}y_0 \, , \quad B = \frac{\delta^2}{\Delta^2}x_0 + \frac{\delta|\Omega|}{\Delta^2}z_0 \, , \quad C = \frac{|\Omega|^2}{\Delta^2}x_0 - \frac{\delta|\Omega|}{\Delta^2}z_0\]

The general solution reads

\[\begin{split}x(t) = -\frac{\delta}{\Delta}y_0 sin(\Delta t) + \left(\frac{\delta^2}{\Delta^2}x_0 + \frac{\delta|\Omega|}{\Delta^2}z_0\right)cos(\Delta t) + \left( \frac{|\Omega|^2}{\Delta^2}x_0 - \frac{\delta|\Omega|}{\Delta^2}z_0 \right)\\ y(t) = y_0cos(\Delta t) + \left( \frac{\delta}{\Delta}x_0 + \frac{|\Omega|}{\Delta}z_0\right) sin(\Delta t) \\ z(t) = - \frac{|\Omega|}{\Delta}y_0sin(\Delta t) + \left( \frac{\delta|\Omega|}{\Delta^2}x_0 + \frac{|\Omega|^2}{\Delta^2}z_0\right)cos(\Delta t) + \left( \frac{\delta^2}{\Delta^2}z_0 - \frac{\delta|\Omega|}{\Delta^2}x_0 \right)\end{split}\]

The solution can be expressed in terms of the Bloch vector in the rotating frame by using the linear relations between the \(u'_i\) and the \(x_i\) components defined above.

Analytic solution for time-dependent field (at \(\delta=0)\) (to be updated…)

For \(\delta=0\) the solution can be parametrized like in the time-independent case. The only modification is that the argument of the sine and cosine function becomes

\[ \begin{align}\begin{aligned} \theta(t) = \int_0^t dt' |\Omega(t')|\\We observe that :math:`\dot{\theta}(t) = |\Omega(t)|` and so the :math:`x` vector written in terms of :math:`\theta(t)` satisfy the correct EQM with a time-dependent Rabi frequency\end{aligned}\end{align} \]

\[\begin{split}x(t) = x_0 \\ y(t) = y_0cos(\theta(t)) + z_0sin(\theta(t)) \\ z(t) = - y_0sin(\theta(t)) + z_0cos(\theta(t))\end{split}\]

On the basis of these equations we see that the area of the module of the pulse is the object that realizes the \(\pi/2\) or the \(\pi\) conditions: * \(\theta=\pi/2\) : the system initially in the ground state \(u^0=(0,0,-1)\) reach the state with \(u_2=-1\) and \(u_3=0\), that correspond to an occupation level of the excited state equal to \(1/2\). * \(\theta=\pi\) : if \(x_0=0\) the Bloch vector changes sign. If the system starts from the ground state the occupation level of the excited state is maximum.

Since \(\Omega(t)\) is parametrized as

\[ \begin{align}\begin{aligned} \Omega(t) = \Omega_0 e^{-\frac{1}{2}\left(\frac{t-t_0}{w}\right)^2}\\the :math:`n\pi` condition becomes\end{aligned}\end{align} \]

\[\sqrt{2}w|\Omega_0| = n\sqrt{\pi}\]

Numerical solutions for arbitrary field and detuning (to be updated…)

# useful to autoreload the module without restarting the kernel
%load_ext autoreload
%autoreload 2

import numpy as np
from mppi.Utilities import Constants as C
from mppi.Models import TwoLevelSystems as TLS
from mppi.Models import GaussianPulse as G
import matplotlib.pyplot as plt

We present various example to discuss the features of the dynamics of the TLS.

Dynamics for a single gaussian pulse

We define the parameters for the numerical analysis and show the results

T = 600 #fs
tstep = int(1e4)
dipole = 1 +1j # the weight of the real and imaginary part can be changed
time = np.linspace(0,T,tstep)

fwhm =  100 # fs
energy = 1.5 # pulse energy in eV
h_red = C.Planck_reduced_ev_ps*1e3 # hbar in eV*fs
omega = energy/h_red # angular frequency of the pulse

theta = 0.5*np.pi # set the area of the pulse
pars = TLS.pulseParametersFromTheta(dipole,theta,fwhm=fwhm)
pars

time unit: fs
the width parameter of the pulse is 42.466090014400955 fs
Rabi coupling (fs^-1): (0.010434524642511517+0.010434524642511517j)
Rabi coupling (module) (fs^-1): 0.014756646266356059
field amplitude (V/m): 129788828.49284193
field intensity (kW/cm^2) : 4471405.517491325

{'Omega0': (0.010434524642511517+0.010434524642511517j),
 'Omega0_abs': 0.014756646266356059,
 'field_amplitude': 129788828.49284193,
 'intensity': 4471405.517491325}

pulse = G.gaussianPulse(time,energy=energy,amplitude=pars['field_amplitude'],fwhm=fwhm)

time unit: fs - energy unit: eV
period of the oscillations 2.757111798154881 fs
width of the pulse 42.466090014400955 fs
fwhm of the pulse 100 fs

plt.plot(time,pulse)

[<matplotlib.lines.Line2D at 0x7fee646c5130>]

Omega = lambda t : G.gaussianPulse(t,energy=energy,amplitude=pars['Omega0_abs'],fwhm=fwhm,envelope_only=True,verbose=False)

uprime0 = np.array([0.,0.,1.])

delta = 30 # in meV
delta_fsm1 = 1e-3*delta/h_red # in fs^-1

uprime = TLS.solveBlochEq(uprime0,time,Omega,mu12=dipole,delta=delta_fsm1)

plt.figure(figsize=(10,7))
plt.title('Carriers')
plt.plot(time,(1+uprime[2])/2,label='carriers')
plt.legend()
plt.show()

plt.figure(figsize=(10,7))
plt.plot(time,uprime[0],label="u'1")
plt.plot(time,uprime[1],label="u'2")
plt.title("u'1 and u'2")
plt.legend()

<matplotlib.legend.Legend at 0x7fee5c558b50>

u = TLS.convertToRotatingFrame(omega,time,uprime,invert=True)

plt.figure(figsize=(10,7))
plt.plot(time,u[0],label='u1')
plt.plot(time,u[1],label='u2')
plt.title('u1 and u2 in the original frame')
plt.legend()

<matplotlib.legend.Legend at 0x7fee5c4e9550>

Global response of many \(k\)-points to a gaussian pulse. Emergence of the FID.

We consider the global response of the system with a large samplings of \(k\)-points.

To each point there is an associated detuning and we assume a uniform distribution in a given energy range.

import random as rand

numk = 1000
delta_range = 40 # meV
deltas = [0.] # the first kpoint has zero detuning
for i in range(numk-1):
    deltas.append(delta_range*(rand.random()-0.5))

We set the properties of the system and of the pulse. Here we assume for symplicity that all k the points have the same value of the transition dipole

T = 4000 #fs
tstep = int(4e3)
dipole = 1 +1j # the weight of the real and imaginary part can be changed
time = np.linspace(0,T,tstep)

fwhm =  100 # fs
energy = 1.5 # pulse energy in eV
omega = energy/(C.Planck_reduced_ev_ps*1e3) # angular frequency of the pulse

h_red = C.Planck_reduced_ev_ps*1e3 # hbar in eV*fs

We choose a field that realizes the \(\pi/2\) condition (for the points with zero detuning)

theta1= 0.5*np.pi # set the area of the pulse
pars1 = TLS.pulseParametersFromTheta(dipole,theta1,fwhm=fwhm)
pars1

time unit: fs
the width parameter of the pulse is 42.466090014400955 fs
Rabi coupling (fs^-1): (0.010434524642511517+0.010434524642511517j)
Rabi coupling (module) (fs^-1): 0.014756646266356059
field amplitude (V/m): 129788828.49284193
field intensity (kW/cm^2) : 4471405.517491325

{'Omega0': (0.010434524642511517+0.010434524642511517j),
 'Omega0_abs': 0.014756646266356059,
 'field_amplitude': 129788828.49284193,
 'intensity': 4471405.517491325}

Omega = lambda t : G.gaussianPulse(t,energy=energy,amplitude=pars1['Omega0_abs'],fwhm=fwhm,envelope_only=True,verbose=False)

#plt.plot(time,Omega(time))

uprime0 = np.array([0.,0.,1.]) # set the system in the GS before the pulse

uprime = np.zeros([3,len(time),numk])
for k in range(numk):
    delta_fsm1 = 1e-3*deltas[k]/h_red
    uprime[:,:,k] = TLS.solveBlochEq(uprime0,time,Omega,mu12=dipole,delta=delta_fsm1)
uprime_global = np.mean(uprime,axis=2)

k = 3
plt.plot(time,uprime[0,:,k],label="u'1")
plt.plot(time,uprime[1,:,k],label="u'2")
plt.plot(time,0.5*(1+uprime[2,:,k]),label="carriers")
plt.legend()

<matplotlib.legend.Legend at 0x7fee5c391ac0>

plt.figure(figsize=(10,7))
plt.plot(time,uprime_global[0],label="u'1")
plt.plot(time,uprime_global[1],label="u'2")
plt.plot(time,0.5*(1-uprime_global[2]),label="carriers")
plt.legend()

<matplotlib.legend.Legend at 0x7fee5c558b80>

The emergence of the Free Induction Decay (FID) mechanism is clearly visible.

Echo mechanism in the \(\pi/2\) - \(\pi\) configuration

We add a second \(\pi\) pulse to realize the echo mechanism.

theta2= np.pi # set the area of the pulse
pars2 = TLS.pulseParametersFromTheta(dipole,theta2,fwhm=fwhm)
pars2

time unit: fs
the width parameter of the pulse is 42.466090014400955 fs
Rabi coupling (fs^-1): (0.020869049285023034+0.020869049285023034j)
Rabi coupling (module) (fs^-1): 0.029513292532712117
field amplitude (V/m): 259577656.98568386
field intensity (kW/cm^2) : 17885622.0699653

{'Omega0': (0.020869049285023034+0.020869049285023034j),
 'Omega0_abs': 0.029513292532712117,
 'field_amplitude': 259577656.98568386,
 'intensity': 17885622.0699653}

Omega = lambda t : G.doubleGaussianPulse(t,energy=energy,amplitude1=pars1['Omega0_abs'],amplitude2=pars2['Omega0_abs'],
                    fwhm1=fwhm,fwhm2=fwhm,t_start2=800,envelope_only=True,verbose=False)

plt.plot(time,Omega(time))

[<matplotlib.lines.Line2D at 0x7fee5c1fa850>]

uprime0 = np.array([0.,0.,1.]) # set the system in the GS before the pulse

uprime = np.zeros([3,len(time),numk])
for k in range(numk):
    delta_fsm1 = 1e-3*deltas[k]/h_red
    uprime[:,:,k] = TLS.solveBlochEq(uprime0,time,Omega,mu12=dipole,delta=delta_fsm1)
uprime_global = np.mean(uprime,axis=2)

k = 0
plt.plot(time,uprime[0,:,k],label="u'1")
plt.plot(time,uprime[1,:,k],label="u'2")
plt.plot(time,0.5*(1-uprime[2,:,k]),label="carriers")
plt.legend()

<matplotlib.legend.Legend at 0x7fee5c1b07f0>

plt.figure(figsize=(10,7))
plt.plot(time,uprime_global[0],label="u'1")
plt.plot(time,uprime_global[1],label="u'2")
plt.plot(time,0.5*(1-uprime_global[2]),label="carriers")
plt.legend()

<matplotlib.legend.Legend at 0x7fee5c122bb0>

The echo signal emerges.