Dielectric function from Bloch-states dynamics (5.3)

Step 0: Theoretical background

The dynamics of Bloch-states is found by integrating a set of Time Dependent Effective Schrödinger Equations (TD-ESEs). Effective means that an effective Hamiltonian is considered. We consider (see also below) the many-body effective Hamiltonians including electron-electron correlation at mean-field (Hartree) and at the screened-exchange level. In addition, one can also consider effective Hamiltonians based on density-functional theory. The latter choice gives time-dependent density-functional theory.

The coupling between electrons and the external field is described by means of the Modern Theory of Polarization:

Where [math]\displaystyle{ | v_{i \mathbf{k}} \rangle }[/math] are the time-dependent Bloch states (corresponding to the filled Bloch-states at zero-field), [math]\displaystyle{ \mathcal E(t) }[/math] is the external field and the term [math]\displaystyle{ i e \partial k }[/math] is the dipole operator consistent with periodic boundary condition. For more details on this last term, see Ref. and .

[math]\displaystyle{ h_{rt} }[/math] is the Hamiltonian,

.

This contains the equilibrium Hamiltonian, that is the zero-field eigenvalues evaluated through DFT; the quasiparticle corrections, evaluated from GW or estimated from experiment, and the variation of the self-energy. The latter is a functional of the density matrix [math]\displaystyle{ \rho }[/math], which is obtained from the Bloch states.

From the solution of the equation of motion for the Bloch states, the polarization is calculated by means of Berry's phase

.

One can decide the level of theory by selecting the terms to include in the Hamiltonian:

• by including only the first term, one selects the independent-particle approximation
• by including the first and second terms, one selects the independent-quasiparticle approximation
• by including the first, second terms and the Hartree part of the self-energy, one selects the time-dependent Hartree or random-phase approximation
• by including all terms, one selects the time-dependent Screened-exchange or Bethe-Salpeter equation level of theory.

The type of "experiment" to perform, or of the response one would like to extract, depend on the time-dependent field in input and the post-processing of the Berry's phase polarization. The induced polarization is the sum of the responses of the system to all orders in the external field: [math]\displaystyle{ \bf{P}(t) = \sum_{n=-\infty}^{+\infty} \bf{P}_n (t) }[/math]

In the case we are interested in linear response only, we use a delta-like electric field that excites all frequencies of the system. We choose the intensity of the field [math]\displaystyle{ E }[/math] to be weak enough so that the response is dominated by the first order term. The response at all frequencies is obtained by taking the Fourier-transform of the polarization as

[math]\displaystyle{ \chi(\omega) = \frac{P(\omega)}{E} }[/math].

The response [math]\displaystyle{ \chi }[/math] is related to the dielectric function through the relation [math]\displaystyle{ \epsilon(\omega) = 1 + 4 \pi \chi(\omega) }[/math]. See also Ref. .

Notes:

• Length vs Velocity gauge: the equations of motion of Yambo are in the length gauge. Other codes choose instead to work in the velocity gauge. If you want to know more about the advantages and disadvantages of the two gauges read section 2.7 of Ref. .

• Advantages and disadvantages with respect to the density-matrix based formalism: Yambo implements two distinct formalisms for studying the real-time response, the one based on density-matrix and the one based on Boch states. Which formalism is best to use depends on what one wants to simulate. The dynamics formulation in terms of density matrix or non-equilibrium Green's functions allows a systematic treatment of correlation effects and electron-phonon coupling, but the price to pay is that the coupling with the external field is correct only to the linear order. This formalism is important in all those phenomena where electronic relaxation plays a major role. On the other hand, the formulation in terms of the Bloch states treats coupling with the external field in an exact way even beyond the linear regime, and for this reason, it allows the description of phenomena coherent with the external fields and to access, for example, non-linear responses.

Step 1: Prerequisites

In this example, we will consider a single layer of hexagonal boron nitride (hBN). Before running real time simulations in yambo you need to:

• download input files and Yambo databases for this tutorial here: hBN-2D-RT.tar.gz.
• perform the steps described in Prerequisites for Real Time propagation with Yambo (5.3).

Step 2: Independent-particle approximation

Let's start from the lowest level of theory. Here it is assumed you are familiar with the optical spectrum of 2D h-BN at this level of theory.

Run the simulations

You should now be inside the folder YAMBO_TUTORIALS/hBN-2D-RT/YAMBO/FixSymm. First of all, create a directory for the inputs to be run by yambo_nl:

mkdir Inputs_nl

Then, use the command:

 yambo_nl -nl n -F Inputs_nl/01_td_ip.in

to generate the input (OMP, angular and extra fields keywords are omitted for brevity):

nloptics                         # [R] Non-linear spectroscopy
% NLBands
  3 | 6  |                           # [NL] Bands range
%
NLverbosity= "high"              # [NL] Verbosity level (low | high)
NLtime= 55.000000           fs    # [NL] Simulation Time
NLintegrator= "INVINT"           # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA"             # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX")
NLLrcAlpha= 0.000000             # [NL] Long Range Correction
% NLEnRange
 -1.000000 |-1.000000 |      eV    # [NL] Energy range
%
NLEnSteps= 1                     # [NL] Energy steps
NLDamping= 0.000000        eV    # [NL] Damping (or dephasing)
RADLifeTime=-1.000000      fs    # [RT] Radiative life-time (if negative Yambo sets it equal to Phase_LifeTime in NL)
#EvalCurrent                   # [NL] Evaluate the current
#FrPolPerdic                   # [DIP] Force periodicity of polarization respect to the external field
HARRLvcs= 18475            RL    # [HA] Hartree     RL components
EXXRLvcs= 18475            RL    # [XX] Exchange    RL components
% Field1_Freq
 0.100000 | 0.100000 |         eV    # [RT Field1] Frequency
%
Field1_NFreqs= 1                 # [RT Field1] Frequency
Field1_Int=  1000.00       kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000     fs    # [RT Field1] Width
Field1_kind= "DELTA"           # [RT Field1] Kind(SIN|COS|RES|ANTIRES|GAUSS|DELTA|QSSIN)
Field1_pol= "linear"             # [RT Field1] Pol(linear|circular)
% Field1_Dir
 0.000000 | 1.000000 | 0.000000 |        # [RT Field1] Versor
%
Field1_Tstart= 0.010000    fs    # [RT Field1] Initial Time
       

Note that the standard input of yambo_nl has default settings for nonlinear response (e.g. the SOFTSIN for the Field1_kind that set the time-dependence 'kind' of the field). Set the field direction Field1_Dir along y, the field kind to DELTA, the length of the simulation NLtime to 55 fs, select the bands from 3 to 6 (NLBands) and set the NLDamping to zero. The fields you need to change with respect to the default settings are shown above in red.

Now run

 yambo_nl -F Inputs_nl/01_td_ip.in -J TD-IP_nl -C TD-IP_nl &

and then monitor the progress with

 tail -f TD-IP_nl/o-TD-IP_nl.polarization_F1
 ctrl+C

The code produces different files (in the TD-IP_nl directory, all with the suffix TD-IP_nl: o-TD-IP_nl.polarization_F1 that contains the polarization, o-TD-IP_nl.external_potential_F1 the external field we used, and finally r-TD-IP_nl_nloptics a report with all information about the simulation.

Plot the third column of o-TD-IP_nl.polarization_F1 versus the first one (time-variable) to get the time-dependent polarization along the y-direction:

gnuplot> p 'TD-IP_nl/o-TD-IP_nl.polarization_F1' u 1:3 w l

From the plot, it is apparent that the polarization oscillates at multiple frequencies. These are the eigenfrequencies of the system excited by the delta-kick.

Output postprocessing: the dielectric function

Once we obtained the polarization, we Fourier transform the polarization to obtain the dielectric function. As the polarization oscillates at all eigenfrequencies of the system, the real/imaginary part of the dielectric function obtained by the Fourier transform have poles/zeros at the eigenfrequencies of the system.

Use the command:

 ypp_nl -nl -F Inputs_nl/02_ypp_abs.in

to generate the input file:

nonlinear                    # [R] NonLinear Optics Post-Processing
Xorder= 1                    # Max order of the response functions
% TimeRange
-1.000000 |-1.000000 | fs    # Time-window where processing is done
%
ETStpsRt= 1001                # Total Energy steps
% EnRngeRt
  0.00000 | 10.00000 | eV    # Energy range
%
DampMode= "LORENTZIAN"             # Damping type ( NONE | LORENTZIAN | GAUSSIAN )
DampFactor=  0.10000   eV    # Damping parameter

where we set a Lorentzian smearing corresponding to 0.1 eV. Note that due to the finite time of our simulation, we need to introduce a finite smearing to Fourier transform the result.

To run the post-processing, use the command:

ypp_nl -F Inputs_nl/02_ypp_abs.in -J TD-IP_nl -C TD-IP_nl 

and obtain the files for the dielectric constant along with the field direction, the EELS along with the same direction, and the damped polarization.

o-TD-IP_nl.YPP-eps_along_E
o-TD-IP_nl.YPP-eels_along_E
o-TD-IP_nl.YPP-damped_polarization

inside the folder TD-IP_nl

Plot the dielectric constant and compare it with the linear response (if available - or you can generate it as an exercise):

gnuplot
plot  "TD-IP_nl/o-TD-IP_nl.YPP-eps_along_E" u 1:2 w l
rep "../CHI_IP/o-CHI_IP.eps_q1_ip" u 1:2 w l

The differences observed with the spectra obtained from the linear-response formalisms are due to how the dipoles are evaluated. Within the linear-response formalism, dipoles are computed as [math]\displaystyle{ \frac{1}{E_{c \mathbf{k}}-E_{v \mathbf{k}}} \left\langle v \mathbf{k}\left|\mathbf{p}_{\alpha}+\mathrm i\left[V^{\mathrm{NL}}, \mathbf{r}_{\alpha}\right]\right| c \mathbf{k}\right\rangle\ }[/math] (which is valid in the linear response regime). The Bloch-states formalism uses instead covariant dipoles evaluated numerically on the k-mesh (valid at all orders). The differences are due to numerical errors in evaluating the covariant dipoles. These differences disappear with a denser k-mesh.

It is a good practice to converge the dielectric function obtained from the Bloch-state dynamics with respect to the k-points using the linear-response spectrum as a reference. This gives the minimum number of k-points to be used to obtain accurate results.

Step 3: Hartree+Screened exchange approximation

Evaluating the collisions

.

K is the functional derivative of the self--energy with respect to the density matrix and corresponds to the kernel of the Bethe--Salpeter equation. In the language of the yambo code we will refer to it with the name "real time collisions" or simply "collisions". Choosing the HSEX approximation and using the GW energies for the quasi--particle corrections, the time propagation corresponds to a real-time version of the Bethe-Salpeter Equation More details can be found in Ref. (This part is common in between the yambo_nl and the yambo_rt).

Generate the input file to calculate the collisions use the command :

 yambo_nl -X s -e -v hsex -F Inputs_nl/03_coll_hsex.in

The flag -X will tell the code to calculate the dielectric constant that is required for the screened interaction. It could be even computed in an independent calculation and then loaded using the -J flags

em1s                           # [R Xs] Static Inverse Dielectric Matrix
dipoles                        # [R   ] Compute the dipoles
Chimod= "HARTREE"              # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXs
   1 |  20 |                   # [Xs] Polarization function bands
%
NGsBlkXs= 1000 mHa      # [Xs] Response block size
% DmRngeXs
  0.10000 |  0.10000 | eV      # [Xs] Damping range
%
% LongDrXs
 1.000000 | 0.000000 | 0.000000 |        # [Xs] [cc] Electric Field
%

While this is the part of the input file specific for the evaluation of the collisions

collisions                     # [R] Eval the extended Collisions
% COLLBands
  4 |  5  |                  # [COLL] Bands for the collisions
%
HXC_Potential= "SEX+HARTREE"           # [SC] SC HXC Potential
HARRLvcs= 1000 mHa      # [HA] Hartree     RL components
EXXRLvcs= 1000 mHa      # [XX] Exchange    RL components
CORRLvcs= 1000 mHa      # [GW] Correlation RL components

With this input, we calculate the HARTREE plus SEX collisions integrals. Notice that the HARTREE term in principle can be calculated on the fly, but in this way it is more efficient especially for the non-linear response. Here one has to converge the cutoff for the Hartree and the Screened Exchange. Around 5000 mHa is a reasonable value for hBN. In this example we will use 1000 mHa to speed up calculations. The collisions bands COLLBands have to be the same number of bands you want to use in the linear/nonlinear response. Then you run

 yambo_nl -F Inputs_nl/03_coll_hsex.in -J COLL_HSEX -C COLL_HSEX 

The calculation may take sone time (about 1 minute in serial). It produced many binary files ndb.COLLISIONS_HXC_fragment_* which will be needed to perform TD-HSEX simulations.

Run the simulations

To generate the input for the real-time simulation you can run

yambo_nl -nl n -V qp -F Inputs_nl/04_td-sex.in

The input file is

nloptics                       # [R NL] Non-linear optics
% NLBands
  4 |  5 |                     # [NL] Bands
%
NLverbosity= "high"             # [NL] Verbosity level (low | high)
NLtime= 20.00000       fs      # [NL] Simulation Time
NLintegrator= "CRANKNIC"         # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "SEX"           # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX")
NLLrcAlpha= 0.000000           # [NL] Long Range Correction
% NLEnRange
 -1.000000 | -1.000000 | eV      # [NL] Energy range
%
NLEnSteps= 1                   # [NL] Energy steps
NLDamping=  "0.10000    eV      # [NL] Damping
#EvalCurrent                   # [NL] Evaluate the current
#FrPolPerdic                   # [DIP] Force periodicity of polarization respect to the external field
HARRLvcs= 1000         mHa     # [HA] Hartree     RL components
EXXRLvcs= 1000         mHa     # [XX] Exchange    RL components

The next section of the input is about external quasiparticle corrections (EXTQP G). There change only

% GfnQP_E
 3.000000 | 1.000000 | 1.000000 |        # [EXTQP G] E parameters  (c/v) eV|adim|adim
%

The latter line introduces a scissor operator (a rigid shift of the conduction bands) of 3.0 eV. In principle, it is possible to perform a G₀W₀ calculation with Yambo and use the Quasi-particle band structure instead of the rigid shift.

The last section regards the applied field (RT Field1). There change these two lines:

% Field1_Dir
 0.000000 | 1.000000 | 0.000000 |        # [RT Field1] Versor
%
 Field1_kind= "DELTA"             # [RT Field1] Kind(SIN|SOFTSIN|RES|ANTIRES|GAUSS|DELTA|QSSIN)

Run the calculation:

yambo_nl -F Inputs_nl/04_td-sex.in -J "TD-SEX_nl,COLL_HSEX" -C TD-SEX_nl

Output postprocessing: dielectric function

The post-processing of the output does not depend on the level of approximation. The same steps are taken as in the independent-particle approximation.

Use the command:

ypp_nl -F Inputs_nl/ypp_abs.in -J TD-SEX_nl -C TD-SEX_nl

Plot the results:

gnuplot> set xrange [3:10]
gnuplot> set yrange [0:12]
gnuplot> plot "TD-SEX_nl/o-TD-SEX_nl.YPP-eps_along_E" u 1:2 w l title "HSEX/BSE", 
gnuplot> rep  "TD-IP_nl/o-TD-IP_nl.YPP-eps_along_E" u ($1+3):2 w l title "IPA + 3 eV scissor"

For comparison, the linear response results at the BSE level can be obtained following the BSE tutorial.

References

Links for 2025 school

• Back to Modena 2025#Tutorials
• Next step Second-harmonic generation of 2D-hBN