Semiconductor Bloch equations in Wannier gauge with well-behaved dephasing

TL;DR

The paper addresses numerical stability and artifacts in semiconductor Bloch equations (SBEs) for high-harmonic generation when using the Wannier gauge with a constant dephasing operator. It shows that the standard dephasing is ill-defined near band degeneracies and avoided crossings, producing spurious carrier-distribution features and unstable integrations. The authors introduce a soothed dephasing operator (SDO) whose damping depends on energy separation with width $w_S$, preserving diagonal elements while suppressing dephasing for closely spaced states; combined with a comoving Houston basis and FFT-based interpolation, this approach dramatically improves convergence and reduces computation time without materially affecting the high-harmonic spectrum. The work provides practical guidelines and parallel code to enable reliable, gauge-consistent SBEs for solid-state HHG modeling.

Abstract

The semiconductor Bloch equations (SBEs) with a dephasing operator for the microscopic polarizations are a well established approach to simulate high-harmonic spectra in solids. We discuss the impact of the dephasing operator on the stability of the numerical integration of the SBEs in the Wannier gauge. It is shown that the standard approach to apply dephasing is ill-defined in the presence of band crossings and leads to artifacts in the carrier distribution. They are caused by rapid changes of the dephasing operator matrix elements in the Wannier gauge, which render the convergence of the simulation in the stationary basis infeasible. In the comoving basis, also called Houston basis, these rapid changes can be resolved, but only at the cost of a largely increased computation time. As a remedy, we propose a modification of the dephasing operator with reduced magnitude in energetically close subspaces. This approach removes the artifacts in the carrier distribution and significantly speeds up the calculations, while affecting the high-harmonic spectrum only marginally. To foster further development, we provide our parallelized source code.

Figures (8)

• Figure 1: Band structure of CdSe (wurzite) along the A-$\Gamma$-A path and the effect of an avoided crossing on the CDO. a) Band structure along the A-$\Gamma$-A path. The different markers depict the contribution of different WFs for three selected bands. Their size is proportional to the corresponding overlap. The gray rectangle (extended by a factor of 20 in each direction) indicates the region of the avoided crossing presented in b). c) Two selected matrix elements of the CDO in Wannier gauge. The gray stripe (extended by a factor of 20) indicates the region of the avoided crossing presented in d). e) Absolute values of the density matrix in Hamiltonian gauge and f) the dephasing operator in Wannier gauge for the $\mathbf{k}$-points $\mathrm{P}$ and $\mathrm{Q}$ indicated in d). c-f) are based on a simulation of a pulse with $E_0 = 0.9\,\mathrm{V}/\mathrm{nm}$ at time $22.8\,\mathrm{fs}$ for a constant dephasing of $T_2=10\,\mathrm{fs}$.
• Figure 2: a) Time evolution of the excited carriers on a $30\times 30\times 100$ MP grid and $E_0 = 0.9\,\mathrm{V/nm}$ for the different dephasing operators. b, c) Distance of the density matrix in the SB to the reference solution in the CB for varying $N^{\mathbf{k}}_3$ and field strength $E_0$ for b) the CDO and c) the SDO.
• Figure 3: Selected occupations in a) Wannier and b) Hamiltonian gauge at the end of a pulse with $E_0=0.9\,\mathrm{V/nm}$. The density matrix was propagated in the CB along the A-$\Gamma$-A path. The colored and black lines were calculated in the CDO and the SDO, respectively. The WFs which correspond to the shown occupations participate in the avoided crossing.
• Figure 4: High-harmonic spectrum ($E_0=0.9\,\mathrm{V/nm}$) of CdSe a) for the most converged simulation in the CB for CDO and b) the difference to the simulation in the SB (purple line), as well as the simulation in SDO and either $N^{\mathbf{k}}_3=200$ (black line) or $N^{\mathbf{k}}_3=2000$ (orange dashed line). We note that the intensity of the difference spectra in b) are at least six orders of magnitude smaller than the intensity of the reference spectrum.
• Figure 5: Computation times to simulate a single pulse in the CB ($N^{\mathbf{k}}_3=100$) for all dephasing types and different maximum field strengths. a) Normalized to the overall lowest computation time. b) Normalized to the computation time of the propagation with the SDO.