Nonlinear optical response of an extended system

A short introduction

The nonlinear optical (NLO) response of an extended system can be obtained from its ground state electronic structure. Due to large computional cost, the local field effect effect is neglected. Hence, the NLO response is only frequency-dependent. The NLO response can be descibed by nonlinear susceptibility tensors. There are numerious NLO processes, depending on the number of photons involved in the process.

The second harmonic generation (SHG) is a NLO process, in which two photons at the frequency of \(\omega\) generate a photon at the frequency of \(2\omega\). The SHG response is characterized by the second-order (quadratic) susceptibility tensor, defined via

\[P_{\gamma}(t) = \sum_{\alpha\beta} \chi_{\gamma\alpha\beta}^{(2)}(\omega,\omega) E_{\alpha}E_{\beta} e^{-2i\omega t}+\textrm{c.c.}\]

where \(\alpha,\beta={x,y,z}\) denotes the Cartezain coordinates, and \(E_{\alpha}\) and \(E_{\alpha}\) are the polarization and electric fields, respectivly. For bulk systems, \(\chi_{\gamma\alpha\beta}^{(2)}\) is expressed in the units of m/V. The details of computing \(\chi_{\gamma\alpha\beta}^{(2)}\) is documented in Ref. [1].

Example 1: SHG spectrum of semiconductor: Monolayer MoS2

To compute the SHG spectrum of given structure, 3 steps are performed:

1. Ground state (GS) calculation

import numpy as np
from ase.build import mx2
from gpaw import GPAW, FermiDirac
from gpaw.nlopt.matrixel import make_nlodata
from gpaw.nlopt.shg import get_shg

# Make the structure and add the vaccum around the layer
atoms = mx2(formula='MoS2', a=3.16, thickness=3.17, vacuum=5.0)
atoms.center(vacuum=15, axis=2)

# GPAW parameters:
nk = 40
params_gs = dict(
    mode='lcao',
    symmetry={'point_group': False, 'time_reversal': True},
    nbands='nao',
    convergence={'bands': -10},
    parallel={'domain': 1},
    occupations=FermiDirac(width=0.05),
    kpts={'size': (nk, nk, 1), 'gamma': True},
    xc='PBE',
    txt='gs.txt')

# Ground state calculation:
gs_name = 'gs.gpw'
calc = GPAW(**params_gs)
atoms.calc = calc
atoms.get_potential_energy()
atoms.calc.write(gs_name, mode='all')

In this script a normal ground state calculation is performed with coarse kpoint grid. Note that LCAO basis is used here, but the PW basis set can also be used. For a smoother spectrum, a finer mesh should be employed.

2. Get the required matrix elements from the GS

   Here, the matrix elements of momentum are computed. Then, all required quantities such as energies, occupations, and momentum matrix elements are saved in a file (‘mml.npz’). The GS file cane be removed after this step.

   
   # Calculate momentum matrix:
   mml_name = 'mml.npz'
   make_nlodata(gs_name=gs_name, out_name=mml_name)
   
   

   Note that, two optional paramters are available in the gpaw.nlopt.matrixel.make_nlodata() function: ni and nf as the first and last bands used for calculations of SHG.

3. Compute the SHG spectrum

   In this step, the SHG spectrum is calculated using the saved data. There are two well-known gauges that can be used: length gauge or velocity gauge. Formally, they are equivalnent but they may generate different results. Here, SHG susceptibility is computed in both gauges and saved. The SHG susceptibility is a rank-3 symmteric tensor with at most 18 independent components. In addition, the point group symmtery reduce the number of independent tensor elements. Monolayer MoS2 has only one independent tensor element: yyy. A broadening is necessary to obtain smooth graphs, and here 50 meV has been used.

   
   # Shift parameters:
   eta = 0.05  # Broadening in eV
   w_ls = np.linspace(0, 6, 500)  # in eV
   pol = 'yyy'
   
   # LG calculation
   shg_name1 = 'shg_' + pol + '_lg.npy'
   get_shg(
       freqs=w_ls, eta=eta, pol=pol, gauge='lg',
       out_name=shg_name1, mml_name=mml_name)
   
   # VG calculation
   shg_name2 = 'shg_' + pol + '_vg.npy'
   get_shg(
       freqs=w_ls, eta=eta, pol=pol, gauge='vg',
       out_name=shg_name2, mml_name=mml_name)
   

Result

Now the calculated SHG spectra are plotted at the end. Both real and imaginary parts of the computed SHG susceptibilities, obtained from two gauges are shown. The gauge invariance is confirmed from the calculation. Note that the bulk susceptibility (with SI units of m/V) is ill-defined for 2D materials, since the volume cannot be defined without ambiguity in 2D systems. Instead, the sheet susceptibility, expressed in unit of m\(^2\)/V, is an unambiguous quantity for 2D materials. Hence, the bulk susceptibility is transformed to the unambiguous sheet susceptibility by multiplying with the width of the unit cell in the \(z\)-direction.

# web-page: shg.png
import numpy as np
import matplotlib.pyplot as plt
from ase.io import read

# Plot and save both spectra
atoms = read('gs.txt')
cell = atoms.get_cell()
cellsize = atoms.cell.cellpar()
mult = cellsize[2] * 1e-10  # make the sheet sus.
legls = []
res_name = ['shg_yyy_lg.npy', 'shg_yyy_vg.npy']
plt.figure(figsize=(6.0, 4.0), dpi=300)
for ii, name in enumerate(res_name):
    # Load the data
    shg = np.load(name)
    w_l = shg[0]
    plt.plot(np.real(w_l), np.real(mult * shg[1] * 1e18), '-')
    plt.plot(np.real(w_l), np.imag(mult * shg[1] * 1e18), '--')
    legls.append(f'{name}: Re')
    legls.append(f'{name}: Im')
    plt.xlabel(r'$\hbar\omega$ (eV)')
    plt.ylabel(r'$\chi_yyy$ (nm$^2$/V)')
    plt.legend(legls, ncol=2)
plt.tight_layout()
plt.savefig('shg.png', dpi=300)

The figure shown here is generated from scripts: shg_MoS2.py and shg_plot.py. It takes 30 minutes with 16 cpus on Intel Xeon X5570 2.93GHz.

Technical details

There are few points about the implementation that we emphasize:

• The code is parallelized over kpoints.

• The code employs only the time reversal symmtery to improve the speed.

• Refer to [1] for the details on the NLO theory.

Useful tips

It’s important to converge the results with respect to:

• nbands

• nkpt (number of kpoints in gs calc.)

• eta

• ecut

• ftol

• vacuum (if there is)

API

gpaw.nlopt.matrixel.make_nlodata(gs_name='gs.gpw', out_name='mml.npz', spin='all', ni=0, nf=0)

Get all required NLO data and store it in a file.

Writes NLO data to file: w_sk, f_skn, E_skn, p_skvnn.

Parameters:

gs_name:

Ground state file name

out_name:

Output filename

spin:

Which spin channel (‘all’, ‘s0’ , ‘s1’)

ni:

First band to compute the mml.

nf:

Last band to compute the mml (relative to number of bands for nf <= 0).

gpaw.nlopt.shg.get_shg(freqs=[1.0], eta=0.05, pol='yyy', eshift=0.0, gauge='lg', ftol=0.0001, Etol=1e-06, band_n=None, out_name='shg.npy', mml_name='mml.npz')

Calculate RPA SHG spectrum for nonmagnetic semiconductors.

Output: shg.npy file with numpy array containing the spectrum and frequencies.

Parameters:

freqs:

Excitation frequency array (a numpy array or list)

eta:

Broadening, a number or an array (default 0.05 eV)

pol:

Tensor element (default ‘yyy’)

gauge:

Choose the gauge (lg or vg)

Etol, ftol:

Tol. in energy and fermi to consider degeneracy

band_n:

List of bands in the sum (default 0 to nb)

out_name:

Output filename (default ‘shg.npy’)

mml_name:

The momentum filename (default ‘mml.npz’)

gpaw.nlopt.linear.get_chi_tensor(freqs=[1.0], eta=0.05, ftol=0.0001, Etol=1e-06, eshift=0.0, band_n=None, mml_name='mml.npz', out_name=None)

Calculate full linear susceptibility tensor for nonmagnetic semiconductors.

Input:

freqs Excitation frequency array (a numpy array or list) eta Broadening, a number or an array (default 0.05 eV) Etol, ftol Tol. in energy and fermi to consider degeneracy eshift Bandgap correction band_n List of bands in the sum (default 0 to nb) out_name If it is given: output filename mml_name The momentum filename (default ‘mml.npz’)

Output:

chi_vvl The output tensor (3, 3, nw) chi.npy If specified: array containing the spectrum and freqs

gpaw.nlopt.shift.get_shift(freqs=[1.0], eta=0.05, pol='yyy', eshift=0.0, ftol=0.0001, Etol=1e-06, band_n=None, out_name='shift.npy', mml_name='mml.npz')

Calculate RPA shift current for nonmagnetic semiconductors.

Parameters:

• freqs – Excitation frequency array (a numpy array or list).

• eta – Broadening, a number or an array (default 0.05 eV).

• pol – Tensor element (default ‘yyy’).

• Etol – Tol. in energy and fermi to consider degeneracy.

• ftol – Tol. in energy and fermi to consider degeneracy.

• band_n – List of bands in the sum (default 0 to nb).

• out_name – Output filename (default ‘shift.npy’).

• mml_name – The momentum filename (default ‘mml.npz’).

Returns:

Numpy array containing the spectrum and frequencies.

Return type:

np.ndarray

[1] (1,2)

A. Taghizadeh, K. S. Thygesen and T. G. Pedersen Arxiv, 2010.11596 (2020)