Dielectric function from Bloch-states dynamics

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 the many-body effective Hamiltonians seen in the Linear response from real time simulations. In addition, one can also consider effective Hamiltonians based on density-functional theory. The latter choice gives time-dependent density-functional theory.

At difference with what seen in Linear response from real time simulations, 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.

Step 2: Independent-particle approximation

Similar to what you have seen in the tutorial for the density matrix formalism (Real time approach to linear response) we start from the lowest level of theory. Here it is assumed you ran the tutorial on the density-matrix formalism and so you are familiar with the optical spectrum of 2D h-BN at this level of theory. You will recognise that the input and output for the Bloch-state formalism are very similar to those for the density-matrix formalism.

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 -nompi -nl n -F Inputs_nl/01_td_ip.in

to generate the input:

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

 nohup yambo_nl -nompi -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.polarization_F1 that contains the polarization, o.external_potential_F1 the external field we used, and finally r_optics_nloptics a report with all information about the simulation. The simulation takes longer than the density-matrix counterpart. This is the price to pay for computing the polarization as a Berry Phase. While there is no gain for calculating the linear response, this allows computing nonlinear coefficients, which cannot be obtained within the density based approach. Have a look at the time profile of the run towards the end of the r_optics_nloptics (absolute time may differ):

      DIPOLE_buil_covariants :      2.9848s CPU
               NL Integrator :      3.9656s CPU ( 5501 calls,   0.001 sec avg)
                     Dipoles :      5.3532s CPU
              NL Hamiltonian :     23.7586s CPU (11002 calls,   0.002 sec avg)
            NL Berry Pol NEQ :     24.1255s CPU ( 5502 calls,   0.004 sec avg)

Most of the time is spent computing the covariant dipoles and the non-equilibrium Berry Polarization, which involves the inversion of small matrices via the Lapack inversion routine. Some time is also required by the integrator NL_Integrator of the equation of motion inversion, which is coded via the solution of a linear system of equations as implemented in the Lapack library. This integrator is more accurate than the Runge-Kutta 2nd order employed for the corresponding calculation within the density matrix formalism.

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
gnuplot> rep 'TD-IP_rt/o-TD-IP_rt.polarization' u 1:3 w p

Real-time polarization in the y-direction

For comparison, in the figure above, the polarization is obtained either through the Bloch-states or the density-matrix. There are few differences which could be due to the following reasons:

• Coupling with the external field: the yambo_rt evaluates the coupling with the external field through the dipoles (by default computing the commutator [x,Hnl]; the yambo_nl evaluates the coupling using an expression based on the Berry Phase. For the latter, a numerical differentiation in reciprocal space is performed. The 10x10 2D k-mesh may not be fine enough. For denser k-mesh, the Berry phase approach converges to the one based on the dipoles. To avoid this one can use DipApproach="Covariant" in the scheme used by yambo_rt.
• Integrator: the integrator used for the yambo_nl is more accurate. The fact that the differences are more pronounced for long times may hint to inaccuracies in the integration of the equation of motions in yambo_rt. One can use a similar integrator with yambo_rt, by setting in input Integrator= "INV RK2"

Output postprocessing: the dielectric function

Once we obtained the polarization, we Fourier transform the polarization to obtain the dielectric function.

Use the command:

 ypp_nl -nl -F Inputs_nl/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/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 and the density-matrix formalism (available from the Real time approach to linear response tutorial):

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

Imaginary part of the dielectric constant

The differences observed with the spectra obtained from either the density-matrix and linear-response formalisms are due to how the dipoles are evaluated. Within the density-matrix and linear-response formalisms, dipoles are computed in the same way. The Bloch-states formalism uses instead covariant dipoles evaluated numerically on the k-mesh. 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

The evaluation of the collisions. (see Introduction to Real Time propagation in Yambo to remember what are the collisions) This part is common in between the yambo_nl and the yambo_rt. If you have computed the collisions for the corresponding yambo_rt, you can re-used those results. If not, follow the related steps in the real-time BSE tutorial.

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:

nohup 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 (compare with results obtained from the density-matrix formalism):

gnuplot> set xrange [3:10]
gnuplot> set yrange [0:12]
gnuplot> plot "TD-SEX_rt/o-TD-SEX_rt.YPP-eps_along_E" u 1:2 w l
gnuplot> rep "TD-SEX_nl/o-TD-SEX_nl.YPP-eps_along_E" u 1:2 w l
gnuplot> rep "TD-IP_rt/o-TD-IP_rt.YPP-eps_along_E" u ($1+3):2 w l 

In plane absorption of hBN at the TD-HSEX level compared with IP absorption

Similarly to what we have seen in the independent-particle case, the differences between the two approaches are due to the accuracy in evaluating numerically the covariant dipoles in the Bloch-state dynamics approach.

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