View source for Second-harmonic generation of 2D-hBN

= Step 0: Theoretical framework =

In this tutorial, we compute the Second-harmonic generation (SHG) from the dynamics of the Bloch-states. The equation-of-motion of the Bloch-states is explained in the tutorial on [[Linear response from Bloch-states dynamics]]. Independently of the 'experiment' we simulate, the part of the integration of motion stays the same. This means that any experiment, including non-linear optics, can be performed at the level of the theory listed for the computation of the dielectric function:
* independent (quasi)particle
* time-dependent Hartree (RPA level)
* time-dependent DFT (and DPFT)
* time-dependent Hartree+Screened exchange (BSE level)
 
 
What changes when we want to calculate the SHG (or higher harmonics) is:
* the time-dependence of the input electric field and
* the post-processing of the signal  


'''Time dependence of the electric Field''': we choose a sinusoidal electric field, thus a monochromatic laser field. This allows one to expand the polarization in the form <math>\bf{P}(t) = \sum_{n=-\infty}^{+\infty} \bf{p}_n e^{-i\omega_n t} </math> where the coefficient <math>\bf{p}_1,...,\bf{p}_n </math> are related to <math>\chi^{(1)},...,\chi^{(n)} </math> through a coefficient depending on a power of the strength of the field (with the power depending on the order of the response).

At difference with a delta-like perturbation, a real-time simulation gives the response at the laser-field only. Then, to obtain the spectrum for the desired range of frequency, we have to perform so many simulations as the frequencies in the desired range.  

'''Post-processing of the signal''': The switch-on of the electric field corresponds to a weak delta-like kick. Thus, even if it may not be noticeable, the polarization results from both the sinusoidal field and the delta-like kick. The latter excites all eigenfrequencies of the system. Though weak, the resulting signal is comparable with the second-harmonic signal. To eliminate the signal from the eigenfrequencies, we apply a dephasing (<math>\gamma_{deph}</math>). To have a signal useful to sample the second-harmonic signal, we need to wait a time <math>\bar t</math>  much larger than the dephasing-time <math>1/\gamma_{deph} </math>. Note that we cannot choose a too short dephasing time (to shorten the simulation time), as this would result in a too large broadening of the spectrum. After <math>\bar t</math>, we sample $2N+1$ times in a period and use discrete Fourier transform to extract <math>{p}_n</math> for <math> n = 0,...,N</math>, where n=2 gives the SHG. This is schematically illustrated in figure:
     
[[File:Pt analysis.png|600px|center|Non-linear response analysis]]


The scheme from Ref. <ref name="Attaccalite2013">C. Attaccalite and M. Gruning [https://arxiv.org/abs/1309.4012v2 Rev. B, '''88''', 235113 (2013)]</ref> summarised the workflow for computing the SHG (or higher-harmonics) in a energy range <math>[\Omega_1,\Omega_2]</math>.

[[File:Scheme nl.png|350px|center|Schematic representation of real-time calculations]]

= Step 1: Prerequisites =

In this example, we will consider a single layer of hexagonal boron nitride (hBN). Before running real time simulations to compute the SHG in <code>yambo</code> you need to download input files and Yambo databases for this tutorial here: [http://www.yambo-code.org/educational/tutorials/files/hBN-2D-RT.tar.gz hBN-2D-RT.tar.gz].
Unzip the tarball and change directory to <code>YAMBO</code> and run the setup (<code>yambo_nl -nompi</code>). 
Execute the command 
 ypp_nl -fixsym -F 00_removesym.in  
to create the input to remove the symmetries which are not compatible with the field: 
 fixsyms                          # [R] Remove symmetries not consistent with an external perturbation
 % Efield1
  <span style="color:red">1.000000 | 1.000000 | 0.000000 |</span>        # First external Electric Field
 %
 % Efield2
  0.000000 | 0.000000 | 0.000000 |        # Additional external Electric Field
 %
 BField= 0.000000           T     # [MAG] Magnetic field modulus
 Bpsi= 0.000000             deg   # [MAG] Magnetic field psi angle [degree]
 Btheta= 0.000000           deg   # [MAG] Magnetic field theta angle [degree]
 #RmAllSymm                     # Remove all symmetries
 <span style="color:red">RmTimeRev</span>                     # Remove Time Reversal
 #RmSpaceInv                    # Remove Spatial Inversion
 #KeepKGrid                     # Do not expand the k-grid
where we chose the field along the xy direction and we remove time-reversal symmetry.


Then run the pre-processing job:
 ypp_nl -F 00_removesym.in J 00_removesym
This creates

= Step 2: Independent-particle approximation =

Use the command:
  yambo_nl -u -F Inputs_shg/01_td_ip.in
to generate the input (NOTE THIS IS A PLACEHOLDER):

 nloptics                      # [R NL] Non-linear optics
 NL_Threads= 1                # [OPENMP/NL] Number of threads for nl-optics
 % NLBands
   <span style="color:red">3</span> |  <span style="color:red">6 </span>|                   # [NL] Bands
 %
 NLstep=   0.0100       fs    # [NL] Real Time step length
 NLtime=55.000000      fs    # [NL] Simulation Time
 NLverbosity= "<span style="color:red">high</span>"           # [NL] Verbosity level (low | high)
 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
  0.200000 | 8.000000 | eV    # [NL] Energy range 
 %
 NLEnSteps=<span style="color:red"> 1  </span>               # [NL] Energy steps
 NLDamping=<span style="color:red"> 0.000000 </span>   eV    # [NL] Damping
 % ExtF_Dir
  0.000000 | <span style="color:red">1.000000</span>  | 0.000000 |        # [NL ExtF] Versor
 %
 ExtF_kind= <span style="color:red"> "DELTA" </span>         # [NL ExtF] Kind(SIN|SOFTSIN|RES|ANTIRES|GAUSS|DELTA|QSSIN)

ADD A DESCRIPTION OF THE MAIN PARAMETERS, LINKING THEM WITH THE THEORETICAL PART 

Run the simulation:
  yambo_nl -F Inputs_shg/01_td_ip.in -J TD-IP_nl -C TD-IP_nl

TO DO: The run now is performing a real-time simulation per frequency, and therefore it takes approximately X times more time than the linear response simulation (NB check as well length of the simulation and amend)

TO DO: Add a plot of the polarization. Note that it is sinusoidal etc. and the nonlinear part is not visible since every order is as a rule of a thumb 6-7 orders of magnitude smaller than the previous.

== Output post-processing: the dielectric function ==

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

Use the command:

  ypp_nl -u -F Inputs_shg/ypp_shg.in

to generate the input file (THIS IS A PLACEHOLDER):

 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 | <span style="color:red">10.00000</span> | eV    # Energy range
 %
 DampMode= "<span style="color:red">LORENTZIAN</span>"             # Damping type ( NONE | LORENTZIAN | GAUSSIAN )
 DampFactor=  <span style="color:red">0.10000</span>   eV    # Damping parameter

         
TO-DO: describe the main keys and parameters.

To run the post-processing, use the command:
 ypp_nl -F Inputs_shg/ypp_shg.in -J TD-IP_nl -C TD-IP_nl 
TO-DO describe which files are produced

TO-DO: Plot the SHG and comments


= Step 3: Hartree+Screened exchange approximation (BSE level) =

= References = 
<references/>