# Real time approach to non-linear response

## Second Harmonic Generation in AlAs

In this tutorial, we will calculate the second harmonic generation of bulk AlAs, the Yambo databases can be downloaded here: AlAs.tar.gz (10 MB). The first steps of this tutorials are the same as the one on Prerequisites for Real-Time propagation with Yambo with the only difference that we will work on AlAs and we will consider an external field in the direction ` 1.000000 | 1.000000 | 0.000000 |`

. The DFT input are available here: ABINIT or QuantumEspresso.
You can run the DFT calculation with ABINIT with the command:

abinit < AlAs.in > output_AlAs

or using QuantumEspresso:

pw.x -inp AlAs.scf.in > output_scf pw.x -inp AlAs.nscf.in > output_scf

and then import the ABINIT wave-function with the command:

a2y -F alas-o_DS2_KSS

and the QuantumEspresso one, in the folder AlAs.save with the command:

p2y

Now we consider an external field in the `[1,1,0]`

direction and remove symmetries not compatible with this field, as explained in the tutorial Prerequisites for Real-Time propagation with Yambo.

## Real-time simulation for the SHG

You can generate the input file with the command `yambo_nl -u n`

:

nloptics # [R NL] Non-linear optics % NLBands 3 | 6 | # [NL] Bands % NLstep= 0.0100 fs # [NL] Real Time step length NLverbosity= "high" # [NL] Verbosity level (low | high) NLtime=-1.000000 fs # [NL] Simulation Time NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK4/RK2EXP/HEUN/INVINT/CRANKNIC") NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/JGM/SEX/HF") NLLrcAlpha= 0.000000 # [NL] Long Range Correction % NLEnRange 1.000000 | 5.000000 | eV # [NL] Energy range % NLEnSteps= 5 # [NL] Energy steps NLDamping= 0.150000 eV # [NL] Damping % Field1_Dir 1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor Field1_kind= "SOFTSIN" # [NL ExtF] Kind(SIN|SOFTSIN|RES|ANTIRES|GAUSS|DELTA|QSSIN)

The line `|1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor`

referees to the direction of the external field (x,y,0). At present only an external field is present but the code can be easily generalized to deal with more fields.
The default parameters of Yambo/Lumen are already tuned for second-harmonic generation, so the only thing you have to change is the band range, between 3 and 6 and the energy range between 1.0-5.0 eV and the number of energy steps in this interval that we set to 10. Notice that you cannot set to zero the lowest value of the energy range because this will requires a simulation that lasts infinite time, see below. Finally, consider that Yambo performs a separate calculation for each frequency, so if you set many energy steps the computational time grows linearly with this number.
Notice that we set `NLverbosity= "high"`

in this way the code will produce a file for each laser frequency containing the time dependent polarization.
Run `yambo_nl`

. Calculations will take about fifteen minutes on a single processor PC, you have time to study the next section that explains how non-linear response is extracted from the real-time simulations.
In order to speed up calculations, you can run them in parallel.

## Non-linear response with Yambo

In order to calculate the non-linear response, the system is excited with different laser fields with a sinusoidal shape at frequencies [math]\displaystyle{ \omega_1, \omega_2, .... , \omega_n }[/math] determined by the parameters NLEnRange and NLEnSteps. A dephasing term is added to the Hamiltonian [math]\displaystyle{ \gamma = }[/math]NLDamping to simulate a finite broadening and to remove the eigenmodes that are excited by the sudden turning on of the external field. After the dephasing time the outgoing signal is analyzed to extract the non-linear coefficients as shown in the figure below:

In the report file `r_nlinear`

you will find the length of the dephasing part and of the sampling one:

Dephasing Time [fs]: 43.880793 Sampling Time [fs]: 4.185667 Total simulation time [fs]: 48.066460

The length of the dephasing interval is inversely proportional to the damping term [math]\displaystyle{ T_{depth} \simeq 1/\gamma }[/math] while the length of the sampling is dictated by the smallest frequency we are interested in: [math]\displaystyle{ T_{samp} \simeq 1/\omega_1 }[/math]. For this reason if [math]\displaystyle{ \omega_1=0 }[/math] or [math]\displaystyle{ \gamma = 0 }[/math] simulation time goes to infinity. The response at zero frequency can be calculated as limit of small frequency perturbation.
You can have a look to the file `src/nloptics/NL_initialize.F`

to see how simulation lengths are defined.
Calculations can take some time, run them in parallel, the best number of processors = number of frequencies or reduce the number of frequencies step in the SHG.

## Analysis of the results

In the sampling region we suppose that the polarization can be written as: [math]\displaystyle{ \bf{P}(t) = \sum_{n=-\infty}^{+\infty} \bf{p}_n e^{-i\omega_n t} }[/math] where the coefficient [math]\displaystyle{ \bf{p}_1,...,\bf{p}_n }[/math] are related to [math]\displaystyle{ \chi^{(1)},...,\chi^{(n)} }[/math]. We sample the polarization signal at different times and invert the previous equation by truncating the sum at a finite order ^{[1]}.
This is done with the command `ypp_nl -u`

that automatically produce an input with the correct values:

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

where `Xorder`

is the order where the previous sum is truncated, and the TimeRange specifies the sampling region. Notice that differently from the first tutorial, in this case, we do not need Damping in `ypp_nl`

because we already included it in the real-time dynamics. Run `ypp`

and it will produce a file called `o.YPP-X_probe_order_2`

. This file contains the different components of [math]\displaystyle{ \chi^2_{xy*} }[/math] and in particular the columns 6 and 7 correspond to the imaginary and real part of [math]\displaystyle{ \chi^2_{xyz} }[/math]. You can plot the result as:

and compare it with the [math]\displaystyle{ \chi^2_{zxy} }[/math] calculated with more frequencies and with the converged result on a [math]\displaystyle{ 18x18x18 }[/math] k-point grid and bands between 2 and 10. In the figure you can find also the comparison with the results of Ref. ^{[2]}.
Notice that the QuantumEspresso results are slightly different from the Abinit ones. This is due to the different pseudo-potential employed in the calculations if pseudopotentials were the same calculation would have been identical. You can download the script to generate this plot and the converged results here.

**Notice that** : the following parameters are not used in the non-linear response analysis: ` EnRngeRt`

, `ETStpsRt`

, `DampMode`

, `DampFactor`

. The damping factor is set directly in the real-time simulation.

System of units in Non-linear Optics: notice that in Yambo the non-linear response functions are in the gaussian system of units, for example the [math]\displaystyle{ \chi^2 (\omega) }[/math] is in cm/statvolt. In order to convert the non-linear coefficients between different system of units, you can have a look at the Appendix C of the book "Non-linear optics" by Robert W. Boyd.

More tutorials on non-linear response (Third harmonic generation, parallelization, spin-orbit etc..) with Yambo can be found here: More Tutorials.

## How to choose the external field direction for the different X^{2}

In the X^{2}_{abc} response you have three field directions. The first two *a* and *b* are determined by the external field and the third *c*
is the direction where you measure the response.

The external field in Yambo is set using the variable

% ExtF_Dir 1.000000 | 1.000000 | 0.000000 | # [NL ExtF] Versor .

For example if you set in Yambo field direction as:

1.0 | 0.0 | 0.0 | ---> you calculate all components X^{2}_{*xx}0.0 | 1.0 | 0.0 | ---> you calculate all components X^{2}_{*yy}0.0 | 0.0 | 1.0 | ---> you calculate all components X^{2}_{*zz}1.0 | 1.0 | 0.0 | ---> you calculate all components X^{2}_{*xy}and X^{2}_{*yx}1.0 | 0.0 | 1.0 | ---> you calculate all components X^{2}_{*xz}and X^{2}_{*zx}0.0 | 1.0 | 1.0 | ---> you calculate all components X^{2}_{*yz}and X^{2}_{*zy}

the * correspond to the different direction in the response file `o.YPP-X_probe_order_2`

.
In the `o.YPP-X_probe_order_2`

you have seven columns:

The first column is the energy omega, the 2nd and 3rd columns are the imaginary and real part of X^{2}_{x**}the 4th and 5th columns are the imaginary and real part of X^{2}_{y**}the 6th and 7th columns are the imaginary and real part of X^{2}_{z**}

where the ** are the direction of the incoming field, see above.

For example if you want to calculate the X^{2}_{zxx} you put the external field in the direction ` 1.0 | 0.0 | 0.0 |`

and then plot the column 6 and 7 of the `o.YPP-X_probe_order_2`

file.

**Nota bene:** if you change the direction of the external field you have to remove the corresponding symmetries as explained in the tutorial
on real-time linear response tutorial.

## Non-linear response of low-dimensional structures: 2D or 1D

The Yambo code is a 3D periodic code, so when you want to study a low-dimensional system 2D or 1D, you have to use a supercell approach. For example to simulate the non-linear response in a 2D-crystal one, can put the crystal in the xy plane and choose a larger distance in the z-direction in order to reduce the interaction between the periodic replica. Clearly in this case k-point sampling will be only in the x,y directions. For calculations on this type of systems, these two warnings must be taken into account:

- if you calculate GW corrections or include electron-hole interaction in the linear/non-linear response it is a good idea to use a
**Coulomb cutoff**, similar to the 2D BSE case How to treat low dimensional systems. You can add the cutoff just adding the "-r" in the input generation.

- the calculation of SHG and THG are always performed respect to the supercell, therefore you have to
**re-scale**the result to the effective thickness of the layer, usually the inter-layer distance of the corresponding bulk material for 2D systems.

## References

- ↑ C. Attaccalite and M. Grüning, Phys. Rev. B 88, 235113(2013)
- ↑ Luppi et al. PRB B,
**82**, 235201(2010)

Prev: Linear response from real time simulations | Now: ICTP Tutorials --> Non Linear Response | Next: Correlation effects in the non-linear response |