The Bethe-Salpeter equation and Excitons¶
For a brief introduction to the Bethe-Salpeter equation and the details of its implementation in GPAW, see Bethe-Salpeter Equation - Theory.
Absorption spectrum of bulk silicon¶
We start by calculating the ground state density and diagonalizing the resulting Hamiltonian.
Below we will set up the Bethe-Salpeter Hamiltonian in a basis of the 4 valence bands and 4 conduction bands. However, the screened interaction that enters the Hamiltonian needs to be converged with respect the number of unoccupied bands. The calculaton is erfored with the following script gs_Si.py. It takes a few minutes on a single CPU. The last line in the script creates a .gpw file which contains all the informations of the system, including the wavefunctions.
Next we calculate the dynamical dielectric function using the Bethe-Salpeter equation. The imaginary part is proportional to the absorption spectrum. The calculation can be done with the script eps_Si.py, which also calculates the dielectric function within the Random Phase Approximation (see Linear dielectric response of an extended system). It takes about ~12 hours on a single CPU but parallelizes very well. Note the .csv output files that contains the spectre. The script also generates a .dat file that contains the eigenvalues of the BSE eigenvalues for easy application. The spectrum essentially consists of a number of peaks centered on the eigenvalues. It can be plottet with plot_Si.py and the result is shown below
The parameters that need to be converged in the calculation are the k-points in the initial ground state calculation. In addition the following keywords in the BSE object should be converged: the plane wave cutoff ecut, the numbers of bands used to calculate the screened interaction nbands, the list of valence bands valence_bands and the list of conduction bands conduction_bands included in the Hamiltonian. It is also possible to provide an array gw_skn, with GW eigenvalues to be used in the non-interacting part of th Hamiltonian. Here, the indices denote spin, k-points and bands, which has to match the spin, k-point sampling and the number of specified valence and conduction bands in the ground state calculation.
Excitons in monolayer MoS2 with Spin-orbit Coupling¶
Spectrum from the Bethe-Salpeter equation¶
The screening plays a fundamental role in the Bethe-Salpeter equation and for 2D systems the screening requires a special treatment. In particular we use a truncated Coulomb interaction inorder to decouple the screening between periodic images. We refer to Ref. [1] for details on the truncated Coulomb interaction in GPAW. As before, we calculate the ground state of \(MoS_2\) with the script gs_MoS2.py, which takes a few minutes. Note the large density of k-points, which are required to converge the BSE spectrum of two-dimensional systems.
The macroscopic dielectric function is calculated as an average of the microscopic screening over the unit cell. Clearly, for a 2D system this will depend on the unit cell size in the direction orthogonal to the slab and in the converged limit the dielectric function becomes unity. Instead we may calculate the longitudinal part of 2D polarizability which is independent of unit cell size. This is done in RPA as well as BSE with the scripts pol_MoS2.py, which takes ~20 hours on 16 CPUs. Note that the BSE polarizability is calculated with and without Coulomb truncation for comparison. In both case spin-orbit coupling is included through the spinors keyword. We refer to Ref. [2] for details on the spin-orbit implementation. The results can be plottet with plot_MoS2.py and is shown below.
The excitonic effects are much stronger than in the case of Si due to the reduced screening in 2D. In particular, we can identify a distinct spin-orbit split exciton well below the band edge. Note that without Coulomb truncation, the BSE spectrum is shifted upward in energy due the screening of electron-hole interactions from periodic images.
2D screening with and without Coulomb truncation¶
To see the effect of the Coulomb truncation, which eliminates the screening from layers in periodic images, we will now calculate the dielectric constant evaluated at the center of the layer \(z_0\) and averaged in the plane. This is accomplished with
The script get_2d_eps.py carries out this calculations with and without Coulomb truncation and the result is shown below plot_2d_eps.py. Note that the truncated screening is bound to become one at \(\Gamma\) due to the different behavior of Coulomb interaction (in \(q\)-space) in 2D systems. For small values of \(q\) the screening is linear, which makes convergence tricky in standard Brillouin zone sampling schemes. Since the \(\Gamma\)-point is always sampled, the screening is typically underestimated and the exciton binding energy is too high at finite \(k\)-point samplings.
Mott-Wannier model for excitons in 2D materials¶
In 3D materials the Mott-Wannier model of excitons has been highly succesful and simply regards the exciton as a “hydrogen atom” with bindings energies that has been rescaled by the exciton effective mass and dielectric screening. Thus in atomic units the binding energy is
where \(\mu^{-1}=m_v^{-1}+m_c^{-1}\) and \(m_v\) and \(m_c\) are the masses of valence and conduction electrons respectively. The 3D expression relies on the fact that the screening is local in real space and thus approximately independent of \(q\). This is clearly not the case in 2D where we always have
for small \(q\). It is thus expected that the hydrogenic binding energy in 2D becomes renormalized by the slope \(\alpha\) in addition to the effective mass. Indeed in Ref. [3] it was shown that the binding energy in 2D can be approximated by
From the band structure of MoS2 it is straigtforward to obtain \(\mu=0.27\) and all we need now is \(\alpha\). In principle we could read of the slope from the figure above, but there is a more direct an accurate way to do it. As it turns out, the slope is needed for any calculation of the response function in the optical limit and it is simply obtained with the script alpha_MoS2.py. This runs on a single CPU in a minute or so. It should produce a value of \(\alpha=5.27\) Å. Transforming to atomic units and inserting into the formula above yields
which is in good agreement with the BSE computation above