<jupyter_start><jupyter_text>TLS description of optical absorption In this notebook we present a two level system (TLS) based description of optical absorption. We discuss the equationof 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$$H_0 = -\frac{1}{2}\omega_0 \sigma_z \, , $$where $\sigma_z$ is the Pauli matrix. $H_0$ is diagonal in the chosen representation with a energy shift between the two levels given by  $\omega_0$. The eigenstates are represented as$$|1\rangle = \left(\begin{array}{c}1\\0\\\end{array} \right) \, , \quad|2\rangle = \left(\begin{array}{c}0\\1\\\end{array} \right) \, .$$Note that in this notation $|1\rangle$ describes the GS of the system, while $|2\rangle$ is the excited state. 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$$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 usingthe structure of both the Hamitonian defined above. We obtain:$$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})$$where $H^I_{12} = \langle 1|H^I|2\rangle$ and $H^I_{21} = \langle 2|H^I|1\rangle$, furthermore$H^I_{12} = H^{I*}_{21}$ since $H^I$ is hermitian. These equations imply that $\rho_{21}$ is the complexconjugate of $\rho_{12}$ and that $\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).It is convenient to describe the system in terms of _polarization_ and of the _inversion_ variables defined as$$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$$\dot{p} = i\omega_0 p +i H_I I \\\dot{I} = -2i\left(H_Ip^*-H_I^*p\right)$$  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$$H^I = - \mathbf{\mu}\cdot\mathbf{E}(t)$$where $\mathbf{\mu}$ is the dipole operator. We choose a linear polarized electric field (say along the $x$ axis).The time dependece of the field is described by a monocromatic wave of energy $\omega$ times a real slow enevelope function $E_0(t)$.Due to this assumption the matrix element $H^I_{12}$ can be expressed as$$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 $\mu^{x}_{12}$ is matrix element of the component of the dipole parallel to $\mathbf{E}$ evaluated between the states 1 and 2 andwe have introduced the (time-dependent) _Rabi frequency_ $\Omega(t)$:$$\Omega(t) = \Omega_0e^{-\frac{(t-t_0)^2}{2w^2}} \, , \, \mathrm{where} \,\,\Omega_0 = \mu^{x}_{12}E_0$$here $\Omega_0$ is the usual (constant) Rabi frequency (we have set $\hbar=1$ in this analysis). Note that since $\mu^{x}_{12}$ is a complexquantity so it is the Rabi frequency and we introduce its real and imaginary parts $\Omega_R(t)$ and $\Omega_I(t)$ (as well as the associated constant factors $\Omega^0_R$ and $\Omega^0_I$) for the susequent analysis. 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$$\omega = \omega_0 + \delta$$where the condition $\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 $p$ can be written as$$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 theaddend $e^{i\omega t'}$ gives relevant contributions since the fast oscillating terms cancel. So the interaction matrix element in the RWAreads$$H_I = \frac{i}{2} \Omega(t)e^{i\omega t}$$Accordingly the EQM are expressed as$$\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$$Note that this argument is usually discussed and justified for __constant__ Rabi frequency, however if the envelope function $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. 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$$p' = e^{-i\omega t}p \, , \quad I' = I$$In this frame the EQM are$$\dot{p}' = -i\delta p' -\frac{1}{2}\Omega I' \\\dot{I}' = \Omega p'^* +\Omega^* p'$$  Bloch vector representation The Bloch vector representation is a change of variables from the polarization and inversionto the 3-dimensional vector $u$ defined as$$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 theBloch vecotr are real quantities and can be easily managed for numerical integration of the EQM.The EQM in this basis read$$\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$$  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$$\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$$  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$$P = p x_{12}^* + p^* x_{12}$$that can be expressed in the terms of the components of the Bloch vector as$$P = u_1 Re(x_{21}) + u_2 Im(x_{21})$$So we observe that the elemnts $u_1$ and $u_2$ are responsible for the behavior of the polarization.  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 theRabi coupling, since $\Omega_0 = Ad/\hbar$. This procedure preserves the value of the observable(the mean on the k point index is understood)$$\mu = Tr(\hat{\rho}\hat{d}) = pd^*+p^*d$$Indeed we can express the complex dipole as $e^{i\alpha}|d|$ and if we define complex rotated polarization as$$\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 $$\dot{\hat{p}} = -i\delta\hat{p} -\frac{1}{2}\Omega I \\\dot{I} = \Omega\left(\hat{p} +\hat{p}^*\right)$$where in these equations the envelope function $\Omega(t)$ is real and depends on the module of theRabi coupling $|\Omega_0|$.The observable polarization in this basis is expressed as$$\mu = |d|(\hat{p}+\hat{p}^*)$$So in this basis we deal with EQM with real parameters that provide equivalent results for the observable.The equations for the Bloch vector in this basis read$$\dot{u}_1 = -\delta u_2 - \Omega u_3\\\dot{u}_2 = \delta u_1  \\\dot{u}_3 = \Omega u_1$$while the observable polarization is computed as $\mu=|d|u_1$.  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$$x(t) = Asin(\Delta t) + Bcos(\Delta t) + C \,\,\, , with \quad\Delta = \sqrt{\delta^2+|\Omega|^2}$$the constants $A,B,C$ will be determined imposing the boundary conditions. Plugging this guess in the EQM for $x$ we obtain$$y(t) = \frac{B\Delta}{\delta}sin(\Delta t) - \frac{A\Delta}{\delta}cos(\Delta t) $$and pluggin this expression in the EQM for $y$ allows us to derive the expression of $z$ $$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 $z$ we obtain $-|\Omega| y(t)$ as expected, so the guess is correct. The constants are determined through the boundary conditions$$x_0 = B+C \, , \quady_0 = - \frac{A\Delta}{\delta} \, , \quadz_0 = \frac{B|\Omega|}{\delta} -\frac{C\delta}{|\Omega|}$$and their solution reads$$A = -\frac{\delta}{\Delta}y_0 \, , \quadB = \frac{\delta^2}{\Delta^2}x_0 + \frac{\delta|\Omega|}{\Delta^2}z_0  \, , \quadC = \frac{|\Omega|^2}{\Delta^2}x_0 - \frac{\delta|\Omega|}{\Delta^2}z_0$$ The general solution reads$$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)$$ 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 modificationis that the argument of the sine and cosine function becomes $$\theta(t) = \int_0^t dt' |\Omega(t')|$$We observe that $\dot{\theta}(t) = |\Omega(t)|$ and so the $x$ vector written in terms of $\theta(t)$ satisfy the correct EQM witha time-dependent Rabi frequency$$x(t) = x_0 \\y(t) = y_0cos(\theta(t))  + z_0sin(\theta(t)) \\z(t) = - y_0sin(\theta(t)) + z_0cos(\theta(t))$$ 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$$\Omega(t) = \Omega_0 e^{-\frac{1}{2}\left(\frac{t-t_0}{w}\right)^2}$$the $n\pi$ condition becomes$$\sqrt{2}w|\Omega_0| = n\sqrt{\pi}$$  Numerical solutions for arbitrary field and detuning (to be updated...)<jupyter_code># 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<jupyter_output><empty_output><jupyter_text>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<jupyter_code>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
pulse = G.gaussianPulse(time,energy=energy,amplitude=pars['field_amplitude'],fwhm=fwhm)
plt.plot(time,pulse)
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()
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()<jupyter_output><empty_output><jupyter_text>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 energyrange.<jupyter_code>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))<jupyter_output><empty_output><jupyter_text>We set the properties of the system and of the pulse. Here we assume for symplicity that all k the points havethe same value of the transition dipole<jupyter_code>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<jupyter_output><empty_output><jupyter_text>We choose a field that realizes the $\pi/$2 condition (for the points with zero detuning)<jupyter_code>theta1= 0.5*np.pi # set the area of the pulse
pars1 = TLS.pulseParametersFromTheta(dipole,theta1,fwhm=fwhm)
pars1
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()
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()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>theta2= np.pi # set the area of the pulse
pars2 = TLS.pulseParametersFromTheta(dipole,theta2,fwhm=fwhm)
pars2
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))
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()
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()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14797078893810606, 'fraction_numerical': 0.02657727762480546, 'mean_word_length': 4.934991119005328, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '438111', 'license': 'MIT', 'n_tokens_mistral': 5913, 'n_tokens_neox': 5591, 'n_words': 1902, 'path': 'sphinx_source/tutorials/Model_TLS_optical_absorption.ipynb', 'repo_name': 'marcodalessandro76/MPPI', 'size': '16706'}