Adiabatic corrections for velocity-gauge simulations of electron dynamics in periodic potentials

TL;DR

This paper introduces analytical corrections for velocity-gauge simulations that reduce computational complexity and resolve divergence issues in modeling electron dynamics in crystalline solids.

Contribution

The authors derive band-energy-dependent corrections that improve convergence and eliminate divergences in velocity-gauge electron dynamics simulations.

Findings

Significantly reduces the number of energy bands needed.

Overcomes divergence problems in susceptibility calculations.

Enhances accuracy and efficiency of light-matter interaction modeling.

Abstract

We show how to significantly reduce the number of energy bands required to model the interaction of light with crystalline solids in the velocity gauge. We achieve this by deriving analytical corrections to the electric current density. These corrections depend only on band energies, the matrix elements of the momentum operator, and the macroscopic vector potential. Thus, the corrections can be evaluated independently from modeling the interaction with light. In addition to improving the convergence of velocity-gauge calculations, our analytical approach overcomes the long-standing problem of divergences in expressions for linear and nonlinear susceptibilities.

Equations

$$\hat{H}(t)={\hat{H}}^{(0)}+{\hat{H}}_{\mathrm{int}}={\hat{H}}^{(0)}+\frac{e}{m_0}\mathbf{A}(t)\cdot \hat{\mathbf{p}}.$$

$${\mathbf{p}}_{nm}(\mathbf{k})=\frac{1}{\Omega }{\int }_{\mathrm{cell}}{\mathrm{d}}^3(\mathbf{r}){\psi }_{n\mathbf{k}}^{\ast }(\mathbf{r})\hat{\mathbf{p}}{\psi }_{m\mathbf{k}}(\mathbf{r}),$$

$$|{\psi }_{\mathbf{k}}(t)⟩=\sum _m{\alpha }_{m\mathbf{k}}(t)|m\mathbf{k}⟩.$$

$$\mathrm{i}\hslash \frac{\partial {\alpha }_{q\mathbf{k}}}{\partial t}={ϵ}_q(\mathbf{k}){\alpha }_{q\mathbf{k}}+\frac{e}{m_0}\mathbf{A}(t)\cdot \sum _m{\mathbf{p}}_{qm}(\mathbf{k}){\alpha }_{m\mathbf{k}}.$$

$$\mathbf{J}(t)=\sum _{n\in \mathrm{occ}}{\int }_{\mathrm{BZ}}\frac{{\mathrm{d}}^3\mathbf{k}}{(2\pi {)}^3}{f}_n(\mathbf{k}){\mathbf{j}}_{n\mathbf{k}}(t).$$

$$\mathbf{j}_{n\mathbf{k}}(t)=-\frac{e}{m_{0}}\Biggl{(}e\mathbf{A}(t)+\\ \operatorname{Re}\biggl{[}\sum_{ij}\left(\alpha_{i\mathbf{k}}^{(n)}(t)\right)^{*}\alpha_{j\mathbf{k}}^{(n)}(t)\mathbf{p}_{ij}(\mathbf{k})\biggr{]}\Biggr{)}.$$

$$\mathbf{P}(t)={\int }_{-\infty }^t\mathrm{d}{t}^{\prime }\mathbf{J}(t^{\prime }).$$

$$\mathbf{A}(t)=\mathrm{Re}\left[\mathbf{a}{\mathrm{e}}^{\gamma t-\mathrm{i}{\omega }_{0}t}\right],$$

$$\lim_{\gamma\to 0+}\lim_{\omega_{0}\to 0}\mathbf{J}(t)\equiv\mathbf{J}_{\mathrm{ad}}\bigl{(}\mathbf{A}(t)\bigr{)}$$

$${\mathbf{J}}_{\mathrm{ad}}(\mathbf{A})=-\frac{e}{m_0}\sum _{n\in \mathrm{occ}}{\int }_{\mathrm{BZ}}\frac{{\mathrm{d}}^3\mathbf{k}}{(2\pi {)}^3}{f}_n(\mathbf{k}){\mathbf{p}}_{nn}\left(\mathbf{k}+\frac{e}{\hslash }\mathbf{A}\right).$$

$$\Delta\mathbf{J}_{N}=\lim_{\gamma\to 0+}\lim_{\omega_{0}\to 0}\bigl{[}\mathbf{J}(t)-\mathbf{J}_{N}(t)\bigr{]}=\mathbf{J}_{\mathrm{ad}}-\lim_{\gamma\to 0+}\lim_{\omega_{0}\to 0}\mathbf{J}_{N}.$$

$${\mathbf{J}}_{\mathrm{ad}}(\mathbf{A})={\mathbf{J}}_{\mathrm{ad}}^{(0)}+{\mathbf{J}}_{\mathrm{ad}}^{(1)}(\mathbf{A})+{\mathbf{J}}_{\mathrm{ad}}^{(2)}(\mathbf{A})+{\mathbf{J}}_{\mathrm{ad}}^{(3)}(\mathbf{A})+\dots$$

$${\mathbf{J}}_{\mathrm{ad}}^{(0)}$$

$${\mathbf{J}}_{\mathrm{ad}}^{(1)}$$

$$\alpha_{q\mathbf{k}}^{(n)}(t)\approx\exp\left[-\frac{\mathrm{i}}{\hbar}\epsilon_{n}(\mathbf{k})t\right]\delta_{qn}-\\ \frac{e}{2\hbar m_{0}}\exp\left[\left(\gamma-\frac{\mathrm{i}}{\hbar}\epsilon_{n}(\mathbf{k})\right)t\right]\times\\ \Biggl{(}\frac{\mathrm{e}^{-\mathrm{i}\omega_{0}t}\bigl{(}\mathbf{p}_{qn}(\mathbf{k})\cdot\mathbf{a}\bigr{)}}{\omega_{qn}(\mathbf{k})-\omega_{0}-\mathrm{i}\gamma}+\frac{\mathrm{e}^{\mathrm{i}\omega_{0}t}\bigl{(}\mathbf{p}_{qn}(\mathbf{k})\cdot\mathbf{a}^{*}\bigr{)}}{\omega_{qn}(\mathbf{k})+\omega_{0}-\mathrm{i}\gamma}\Biggr{)},$$

$${\omega }_{qn}(\mathbf{k})=\left({ϵ}_q(\mathbf{k})-{ϵ}_n(\mathbf{k})\right)/\hslash .$$

$$\lim_{\omega_{0}\to 0}\mathbf{j}_{n\mathbf{k}}=-\frac{e}{m_{0}}\mathbf{p}_{nn}+\frac{e^{2}}{m_{0}}\mathbf{A}-\frac{2e^{2}}{\hbar m_{0}^{2}}\times\\ \sum_{i}\frac{\operatorname{Re}\left[(\omega_{in}+\mathrm{i}\gamma)\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\mathbf{p}_{ni}\right]}{\gamma^{2}+\omega_{in}^{2}}.$$

$$\lim_{\gamma\to 0+}\lim_{\omega_{0}\to 0}\mathbf{J}\approx-\frac{e}{m_{0}}\sum_{n\in\mathrm{occ}}\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\mathbf{k}}{(2\pi)^{3}}\,f_{n}(\mathbf{k})\mathbf{p}_{nn}-\\ \frac{e^{2}}{m_{0}}\sum_{n\in\mathrm{occ}}\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\mathbf{k}}{(2\pi)^{3}}\,f_{n}(\mathbf{k})\Biggl{\{}\mathbf{A}-\frac{2}{m_{0}}\sum_{i\neq n}\frac{\operatorname{Re}\Bigl{[}\bigl{(}\mathbf{p}_{ni}\cdot\mathbf{A}\bigr{)}\mathbf{p}_{in}\Bigr{]}}{\hbar\omega_{in}}\Biggr{\}}.$$

$$\frac{2}{m_{0}}\sum_{i\neq n}\frac{\bigl{|}\bigl{(}\mathbf{e}_{\beta}\cdot\mathbf{p}_{in}\bigr{)}\bigr{|}^{2}}{\hbar\omega_{in}}=1-\hbar\frac{\partial\,\bigl{(}\mathbf{e}_{\beta}\cdot\mathbf{p}_{nn}\bigr{)}}{\partial k_{\beta}},$$

$$\Delta\mathbf{J}^{(1)}=\mathbf{J}_{\mathrm{ad}}^{(1)}(\mathbf{A})+\frac{e^{2}}{m_{0}}\sum_{n\in\mathrm{occ}}\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\mathbf{k}}{(2\pi)^{3}}\,f_{n}(\mathbf{k})\Biggl{\{}\mathbf{A}-\\ \frac{2}{m_{0}}\sum_{i\neq n}^{N}\frac{\operatorname{Re}\Bigl{[}\bigl{(}\mathbf{p}_{ni}\cdot\mathbf{A}\bigr{)}\mathbf{p}_{in}\Bigr{]}}{\hbar\omega_{in}}\Biggr{\}}\\ =\frac{e^{2}}{m_{0}}\sum_{n\in\mathrm{occ}}\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\mathbf{k}}{(2\pi)^{3}}\,f_{n}(\mathbf{k})\Biggl{\{}\mathbf{A}-\frac{1}{\hbar}\bigl{(}\mathbf{A}\cdot\nabla_{\mathbf{k}}\bigr{)}\mathbf{p}_{nn}-\\ \frac{2}{m_{0}}\sum_{i\neq n}^{N}\frac{\operatorname{Re}\Bigl{[}\bigl{(}\mathbf{p}_{ni}\cdot\mathbf{A}\bigr{)}\mathbf{p}_{in}\Bigr{]}}{\hbar\omega_{in}}\Biggr{\}}.$$

$$\Delta\mathbf{J}_{\mathrm{dielectric}}^{(1)}=\frac{e^{2}}{m_{0}}\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\mathbf{k}}{(2\pi)^{3}}\,\Biggl{\{}N_{\mathrm{VB}}\mathbf{A}-\\ \frac{2}{m_{0}}\sum_{n\in\mathrm{VB}}\sum_{i\neq n}^{N}\frac{\operatorname{Re}\Bigl{[}\bigl{(}\mathbf{p}_{ni}\cdot\mathbf{A}\bigr{)}\mathbf{p}_{in}\Bigr{]}}{\hbar\omega_{in}}\Biggr{\}}.$$

$$\Delta\mathbf{J}_{\mathrm{dielectric}}^{(1)}=\mathbf{A}(t)\frac{e^{2}}{m_{0}}\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\mathbf{k}}{(2\pi)^{3}}\,\Biggl{\{}N_{\mathrm{VB}}-\\ \frac{2}{m_{0}}\sum_{n\in\mathrm{occ}}\sum_{i\neq n}^{N}\frac{\Bigl{|}\mathbf{e}_{\mathrm{L}}\cdot\mathbf{p}_{in}\Bigr{|}^{2}}{\hbar\omega_{in}}\Biggr{\}},$$

$$N_{\mathrm{VB}}^{\mathrm{eff}}=\frac{2\Omega}{m_{0}}\sum_{n\in\mathrm{VB}}\sum_{i\neq n}^{N}\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\mathbf{k}}{(2\pi)^{3}}\,\frac{\Bigl{|}\mathbf{e}_{\mathrm{L}}\cdot\mathbf{p}_{in}\Bigr{|}^{2}}{\hbar\omega_{in}}.$$

$$\Delta\mathbf{J}^{(2)}=\Delta\mathbf{J}_{\mathrm{ad}}^{(2)}+\\ \frac{e^{3}}{\hbar^{2}m_{0}^{3}}\sum_{n\in\mathrm{occ}}\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\mathbf{k}}{(2\pi)^{3}}\,f_{n}(\mathbf{k})\Biggl{\{}\Biggl{(}\sum_{i\neq n}^{N}\sum_{j\neq n}^{N}\frac{1}{\omega_{in}\omega_{jn}}\\ \times\Biggl{[}\mathbf{p}_{jn}\bigl{(}\mathbf{p}_{ij}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{ni}\cdot\mathbf{A}\bigr{)}+\\ \bigl{(}\mathbf{p}_{jn}\cdot\mathbf{A}\bigr{)}\biggl{(}\mathbf{p}_{ij}\bigl{(}\mathbf{p}_{ni}\cdot\mathbf{A}\bigr{)}+\mathbf{p}_{ni}\bigl{(}\mathbf{p}_{ij}\cdot\mathbf{A}\bigr{)}\biggr{)}\Biggr{]}\Biggr{)}-\\ \sum_{i\neq n}^{N}\frac{2\bigl{(}\mathbf{p}_{nn}\cdot\mathbf{A}\bigr{)}\operatorname{Re}\left[\mathbf{p}_{ni}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\right]+\mathbf{p}_{nn}\left|\mathbf{p}_{ni}\cdot\mathbf{A}\right|^{2}}{\omega_{in}^{2}}\Biggr{\}}.$$

$$\Delta\mathbf{J}^{(3)}=\Delta\mathbf{J}_{\mathrm{ad}}^{(3)}+\frac{e^{4}}{\hbar^{3}m_{0}^{4}}\sum_{n\in\mathrm{occ}}\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\mathbf{k}}{(2\pi)^{3}}\,f_{n}(\mathbf{k})\Biggl{\{}-\sum_{i\neq n}^{N}\sum_{j\neq n}^{N}\sum_{k\neq n}^{N}\Biggl{(}\frac{1}{\omega_{in}\omega_{jn}\omega_{kn}}\operatorname{Re}\Biggl{[}\bigl{(}\mathbf{p}_{ji}\cdot\mathbf{A}\bigr{)}\times\\ \biggl{(}\bigl{(}\mathbf{p}_{nk}\cdot\mathbf{A}\bigr{)}\Bigl{[}\mathbf{p}_{in}\bigl{(}\mathbf{p}_{kj}\cdot\mathbf{A}\bigr{)}+\mathbf{p}_{kj}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\Bigr{]}+\mathbf{p}_{nk}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{kj}\cdot\mathbf{A}\bigr{)}\biggr{)}+\\ \mathbf{p}_{ji}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{kj}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{nk}\cdot\mathbf{A}\bigr{)}\Biggr{]}\Biggr{)}+\sum_{i\neq n}^{N}\sum_{j\neq n}^{N}\Biggl{(}\frac{\omega_{in}+\omega_{jn}}{2\omega_{in}^{2}\omega_{jn}^{2}}\Biggl{[}\mathbf{p}_{ji}\bigl{(}\mathbf{p}_{nn}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{nj}\cdot\mathbf{A}\bigr{)}+\\ \mathbf{p}_{jn}\bigl{(}\mathbf{p}_{ni}\cdot\mathbf{A}\bigr{)}\Bigl{[}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{nj}\cdot\mathbf{A}\bigr{)}+\bigl{(}\mathbf{p}_{nn}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{ij}\cdot\mathbf{A}\bigr{)}\Bigr{]}+\\ \bigl{(}\mathbf{p}_{ji}\cdot\mathbf{A}\bigr{)}\Bigl{[}\bigl{(}\mathbf{p}_{nn}\cdot\mathbf{A}\bigr{)}\Bigl{(}\mathbf{p}_{in}\bigl{(}\mathbf{p}_{nj}\cdot\mathbf{A}\bigr{)}+\mathbf{p}_{nj}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\Bigr{)}+\mathbf{p}_{nn}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{nj}\cdot\mathbf{A}\bigr{)}\Bigr{]}+\\ \bigl{(}\mathbf{p}_{jn}\cdot\mathbf{A}\bigr{)}\Biggl{(}\mathbf{p}_{ni}\Bigl{(}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{nj}\cdot\mathbf{A}\bigr{)}+\bigl{(}\mathbf{p}_{nn}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{ij}\cdot\mathbf{A}\bigr{)}\Bigr{)}+\\ \bigl{(}\mathbf{p}_{ni}\cdot\mathbf{A}\bigr{)}\Bigl{[}\mathbf{p}_{ij}\bigl{(}\mathbf{p}_{nn}\cdot\mathbf{A}\bigr{)}+\mathbf{p}_{in}\bigl{(}\mathbf{p}_{nj}\cdot\mathbf{A}\bigr{)}+\mathbf{p}_{nj}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}+\mathbf{p}_{nn}\bigl{(}\mathbf{p}_{ij}\cdot\mathbf{A}\bigr{)}\Bigr{]}\Biggr{)}\Biggr{]}\Biggr{)}-\\ \bigl{(}\mathbf{p}_{nn}\cdot\mathbf{A}\bigr{)}\sum_{i\neq n}^{N}\frac{\mathbf{p}_{ni}\bigl{(}\mathbf{p}_{nn}\cdot\mathbf{A}\bigr{)}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}+\bigl{(}\mathbf{p}_{ni}\cdot\mathbf{A}\bigr{)}\Bigl{(}\mathbf{p}_{in}\bigl{(}\mathbf{p}_{nn}\cdot\mathbf{A}\bigr{)}+2\mathbf{p}_{nn}\bigl{(}\mathbf{p}_{in}\cdot\mathbf{A}\bigr{)}\Bigr{)}}{\omega_{in}^{3}}\Biggr{\}}.$$

$$\Delta J(t)=\sum _q{c}_q{A}^q(t).$$

$$c_{1}=\frac{e^{2}}{m_{0}}\int_{\mathrm{BZ}}\frac{\mathrm{d}k}{2\pi}\,\Biggl{\{}N_{\mathrm{VB}}-\frac{2}{\hbar m_{0}}\sum_{n\in\mathrm{VB}}\sum_{i\neq n}^{N}\frac{|p_{in}|^{2}}{\omega_{in}}\Biggr{\}},$$

$$c_{2}=\frac{3e^{3}}{\hbar^{2}m_{0}^{3}}\sum_{n\in\mathrm{VB}}\int_{\mathrm{BZ}}\frac{\mathrm{d}k}{2\pi}\,\Biggl{\{}\sum_{i\neq n}^{N}\sum_{j\neq n}^{N}\frac{p_{ij}p_{ni}p_{jn}}{\omega_{in}\omega_{jn}}-\\ p_{nn}\sum_{i\neq n}^{N}\frac{\left|p_{in}\right|^{2}}{\omega_{in}^{2}}\Biggr{\}},$$

$$c_{3}=-\frac{4e^{4}}{\hbar^{3}m_{0}^{4}}\times\\ \sum_{n\in\mathrm{VB}}\int_{\mathrm{BZ}}\frac{\mathrm{d}k}{2\pi}\,\Biggl{\{}\sum_{i\neq n}^{N}\sum_{j\neq n}^{N}\sum_{\ell\neq n}^{N}\frac{\operatorname{Re}[p_{ji}p_{in}p_{\ell j}p_{n\ell}]}{\omega_{in}\omega_{jn}\omega_{\ell n}}-\\ \sum_{i\neq n}^{N}\sum_{j\neq n}^{N}\frac{(\omega_{in}+\omega_{jn})\left(\frac{\left|p_{in}p_{jn}\right|^{2}}{2}+p_{nn}\operatorname{Re}[p_{ij}p_{ni}p_{jn}]\right)}{\omega_{in}^{2}\omega_{jn}^{2}}+\\ p_{nn}^{2}\sum_{i\neq n}^{N}\frac{\left|p_{in}\right|^{2}}{\omega_{in}^{3}}\Biggr{\}}.$$

$$V(x)=\sum_{q}\Biggl{\{}-2.2\,\text{sech}^{2}\bigl{(}0.9(x-qa)\bigr{)}+\\ 0.01\sin\left(\frac{2\pi(x-qa)}{a}\right)\Biggr{\}}.$$

$$A(t)=-\theta ({\tau }_{\mathrm{L}}-|t|)\frac{E_0}{{\omega }_0}{\cos }^4\left(\frac{\pi t}{2{\tau }_{\mathrm{L}}}\right)\sin ({\omega }_{0}t).$$

Full text

I Introduction

The velocity and length gauges are two most frequent choices for a theoretical description of the interaction between electromagnetic radiation and matter. Even though, in principle, all physical observables must be gauge-independent, the gauge choice is very important for numerical simulations because approximations and discretization errors violate the gauge invariance in its strict sense Scully and Zubairy (1997). In atomic and molecular physics, where quantum systems are finite, both gauges are easily implemented. The situation is different for modeling the interaction of light with bulk crystals, where an infinite periodic lattice potential and the degenerate energy states present difficulties for length-gauge implementations Aversa and Sipe (1995); Virk and Sipe (2007). These difficulties also pertain to the closely related approach of modeling quantum dynamics in the basis of Houston states Houston (1940), also known as accelerated Bloch states Krieger and Iafrate (1986). The fundamental problem here is that these methods require differentiation with respect to the crystal momentum, $\mathbf{k}$, while numerically evaluated Bloch states are, in general, not smooth functions of $\mathbf{k}$ Marzari et al. (2012); Attaccalite and Grüning (2013). The problem becomes particularly severe for modeling strong-field phenomena, where electrons traverse a significant part of the Brillouin zone during a laser cycle Ghimire et al. (2011); Schubert et al. (2014); Luu et al. (2015); Wismer et al. (2016). Solutions to these problems are known Virk and Sipe (2007); Grüning et al. (2016), but they either do not ensure the periodicity with respect to $\mathbf{k}$ Lindefelt et al. (2004) or require the evaluation of so-called covariant derivatives Nunes and Gonze (2001); Souza et al. (2004); Virk and Sipe (2007). Evaluating the covariant derivatives at each propagation step slows down computations. In contrast, the velocity gauge does not require differentiation with respect to the crystal momentum, which often makes it advantageous for modeling light-driven electron dynamics Yabana et al. (2012); Goncharov (2013); Simonsen et al. (2014); Wachter et al. (2014); Yakovlev et al. (2015); Chizhova et al. (2016). However, performing velocity-gauge simulations in a basis of Bloch states has a serious drawback: for convergence, they usually require many more energy bands than length-gauge or Houston-basis simulations Wu et al. (2015). In this paper, we propose a solution to this problem. We do so by deriving analytical corrections to the polarization response evaluated with a relatively small number of bands (section III). Although we derive our corrections in the limit of a weak external field, we show in section IV that they work surprisingly well even for strong optical fields.

II Velocity-gauge description of light-solid interaction

Let $\left|n\mathbf{k}\right\rangle$ denote a Bloch state with band index $n$ and crystal momentum $\mathbf{k}$. The Bloch states are eigenstates of an unperturbed Hamiltonian $\hat{H}^{(0)}$. That is, $\hat{H}^{(0)}\left|n\mathbf{k}\right\rangle=\epsilon_{n}(\mathbf{k})\left|n\mathbf{k}\right\rangle$, where $\epsilon_{n}(\mathbf{k})$ is energy. Even though $\hat{H}^{(0)}$ is a single-particle operator, it may represent an outcome of self-consistent many-body calculations. We describe the interaction with a classical electromagnetic field by adding an interaction operator, $\hat{H}_{\mathrm{int}}$, to the Hamiltonian. Unitary transformations allow one to write $\hat{H}_{\mathrm{int}}$ in different forms. One of them is $\hat{H}_{\mathrm{int}}=e(\mathbf{A}\cdot\hat{\mathbf{p}}+\hat{\mathbf{p}}\cdot\mathbf{A})/(2m_{0})+e^{2}A^{2}/m_{0}$. Here, $\mathbf{A}$ is the vector potential of the electromagnetic field, $\hat{\mathbf{p}}$ is the momentum operator, $e>0$ is the elementary charge, and $m_{0}$ is the electron mass. In SI units, the relation between $\mathbf{A}$ and the electric field is $\mathbf{E}=-\partial\mathbf{A}/\partial t$. In the following, we assume the dipole approximation, where the dependence of $\mathbf{A}$ on coordinate, $\mathbf{r}$, is neglected. The term $e^{2}A^{2}/m_{0}$ can be eliminated by the unitary transformation $\left|\psi\right\rangle\to\exp[-\mathrm{i}e^{2}/(2\hbar m_{0})\int_{-\infty}^{t}A^{2}(t^{\prime})\mathrm{d}t^{\prime}]\left|\psi\right\rangle$. The results is the velocity-gauge form of the Hamiltonian:

[TABLE]

In this equation, the external field preserves the spatial periodicity of the lattice potential (which may or may not be local). Therefore, the Bloch theorem applies not only to $\hat{H}^{(0)}$, but also to $\hat{H}(t)$. Consequently, transitions induced by the external field formally preserve the crystal momentum, which is one of the most important reasons why velocity-gauge simulations are so attractive for numerical calculations.

However, the convenience of preserving the crystal momentum comes at a price. This becomes apparent if we consider the interaction of a dielectric with a constant electric field ($\mathbf{E}=\text{const}$) that is sufficiently weak to neglect interband tunneling. Even though, in the absence of free charge carriers, all physical observables are expected to take constant values, the velocity-gauge wave functions have a nontrivial dependence on time because the vector potential, $\mathbf{A}(t)=-\int_{-\infty}^{t}\mathbf{E}(t^{\prime})\mathrm{d}t^{\prime}$, makes the Hamiltonian time-dependent. In this case, numerical approximations, such as truncating the basis, lead to a spurious divergence of the polarization response in the low-frequency limit Aspnes (1972); Moss et al. (1990); Sheik-Bahae et al. (1991); Sipe and Ghahramani (1993); Aversa and Sipe (1995). This is one of the underlying reasons why velocity-gauge calculations require many more bands than length-gauge or Houston-basis ones.

Another, closely related, problem with the velocity gauge can be understood by considering a particular Bloch state exposed to a weak long half-cycle of the electric field. Because the field performs only a half of the oscillation, the vector potential has a nonzero value at the end of the pulse. In the adiabatic limit, the pulse transforms an initial state $\left|n\mathbf{k}_{0}\right\rangle$ into another Bloch state, $\left|n\mathbf{k}\right\rangle$, provided that the $n$-th band is nondegenerate along the reciprocal-space path prescribed by the acceleration theorem: $\mathbf{k}(t)=\mathbf{k}_{0}+e\hbar^{-1}\mathbf{A}(t)$. This adiabatic intraband motion is guaranteed in length-gauge and Houston-basis calculations, making these approaches a good choice when only ad hoc transition matrix elements are available Apalkov and Stockman (2012); Schiffrin et al. (2013). In contrast, the correct adiabatic behavior is not explicitly encoded in the velocity-gauge equations of motion—it is implicit in the parameters of the model, demanding accurate energies and matrix elements.

The only matrix elements required in the velocity-gauge simulations are those of the momentum operator:

[TABLE]

where the integration is performed over a unit cell, $\Omega$ is the cell volume, and $\psi_{n\mathbf{k}}(\mathbf{r})=\left\langle\mathbf{r}|n\mathbf{k}\right\rangle$ are Bloch functions in the coordinate representation. Indeed, let us search for time-dependent wave functions using the following ansatz:

[TABLE]

Inserting this ansatz into the time-dependent Schrödinger equation (TDSE), $\mathrm{i}\hbar\partial_{t}\left|\psi\right\rangle=\hat{H}(t)\left|\psi\right\rangle$, and assuming that $\left\langle m\mathbf{k^{\prime}}|n\mathbf{k}\right\rangle=\delta_{mn}\delta(\mathbf{k}-\mathbf{k^{\prime}})$, we obtain a system of coupled differential equations for the probability amplitudes:

[TABLE]

In these equations, the medium properties are fully defined by $\epsilon_{q}(\mathbf{k})$ and $\mathbf{p}_{qm}(\mathbf{k})$. Once the probability amplitudes are evaluated, the electric current density can be obtained according to

[TABLE]

Here, $f_{n}(\mathbf{k})$ is the Fermi factor of band $n$ in the initial (stationary) state. The summation in Eq. (5) is carried out over all the bands that contain electrons in the unperturbed state (these would be valence bands for a dielectric), the integration is performed over the first Brillouin zone (BZ), and the contribution from crystal momentum $\mathbf{k}$ in initial band $n$ is given by

[TABLE]

In this equation, we have added a superscript to the probability amplitudes in order to denote the initial band: $\alpha_{i\mathbf{k}}^{(n)}(t_{\mathrm{min}})=\delta_{in}$. That is, in each initial state, only one band is occupied; by adding the contributions from all the initial states, as prescribed by Eq. (5), we model a solid with occupied valence bands. Once $\mathbf{J}(t)$ is evaluated, the polarization response, which is the key quantity in linear and nonlinear optics, is readily given by King-Smith and Vanderbilt (1993); Resta (1998)

[TABLE]

III Velocity-gauge corrections

In the previous section, we argued that velocity-gauge simulations tend to violate adiabaticity where adiabaticity is expected, which is why they demand high-quality band energies and matrix elements. In this section, we derive analytical corrections that alleviate this issue in the case where a finite number of energy bands is the main source of discretization errors. In the next section, we will demonstrate that our corrections work well even if the optical response is far from being adiabatic (that is, the polarization $\mathbf{P}$(t) is not a function of the electric field at time $t$). In spite of that, we derive the corrections by considering an adiabatic limit. It is worth clarifying this point before we start the derivation. If velocity-gauge simulations require $N_{\mathrm{VG}}$ bands to produce physically meaningful results, while length-gauge simulations require only $N_{\mathrm{LG}}<N_{\mathrm{VG}}$ bands, then the bands that are unnecessary in the length gauge serve a rather numerical than physical purpose. They are required in the velocity gauge to ensure that the basis set is reasonably complete. After the interaction with an external field, the occupation of these highly excited states is negligibly small. Even if the polarization response due to “physically relevant” bands is very nontrivial, the contributions due to very high conduction bands (and very low valence ones) are expected to be, in some sense, “simple”. We argue that this simplicity consists in the adiabaticity with respect to the vector potential: the difference between the exact current density, $\mathbf{J}(t)$, and that evaluated with a finite number of bands, $\mathbf{J}_{N}(t)$ is, approximately, a function of the vector potential, $\mathbf{A}$, at time $t$. In this case, one may use any function $\mathbf{A}(t)$ that is particularly well suited for analytical calculations. We chose

[TABLE]

where $\mathbf{a}$ is a constant vector (it is complex-valued unless $\mathbf{A}$ is linearly polarized), $\omega_{0}$ is a carrier-wave frequency, and $\gamma>0$ is a small parameter that controls how slowly the field is turned on. The next step will be the evaluation of the current density using time-dependent perturbation theory, where we only consider a time interval where the external field may be viewed as a perturbation. Once this step is accomplished, we will take the adiabatic limit: $\omega_{0}\to 0$. Here comes our key idea. When the amplitude of the vector potential, $|\mathbf{a}|$, is fixed and its frequency is decreased, then the electric field decreases as well; consequently, we expect $\lim_{\gamma\to 0+}\lim_{\omega_{0}\to 0}\mathbf{J}(t)\equiv 0$ for a dielectric and, more generally,

[TABLE]

if there are partially occupied states in the unperturbed solid. The adiabatic intraband current density, $\mathbf{J}_{\mathrm{ad}}$, is produced by charge carriers that change their crystal momentum according to the acceleration theorem, without leaving their initial bands:

[TABLE]

Here, it is sufficient to carry out the sum over partially occupied bands. Indeed, the contribution from empty bands is obviously zero, while for a fully occupied band, where $f_{n}(\mathbf{k})\equiv 1$, the integral can be shown to be zero: $\int_{\mathrm{BZ}}\mathrm{d}^{3}\mathbf{k}\,\mathbf{p}_{nn}(\mathbf{k}+e\hbar^{-1}\mathbf{A})=m_{0}\hbar^{-1}\int_{\mathrm{BZ}}\mathrm{d}^{3}\mathbf{k}\,\nabla_{\mathbf{k}}\epsilon_{n}(\mathbf{k}+e\hbar^{-1}\mathbf{A})=0$.

In the velocity gauge, the equations of motion do not automatically satisfy Eq. (9). We will soon see that this equation is, in general, violated by truncating the basis set, which is an unavoidable operation for numerical calculations. However, knowing that Eq. (9) holds in the ideal case, we arrive at a simple recipe for adiabatic corrections to the electric current density evaluated with a finite number of bands:

[TABLE]

In the adiabatic limit, $\Delta\mathbf{J}_{N}$ is a function of $\mathbf{A}$. In the rest of this section, we will be searching for a decomposition of $\Delta\mathbf{J}_{N}$ into the powers of the vector potential. As a first step, we decompose the intraband adiabatic current:

[TABLE]

The first two terms in this decomposition are

[TABLE]

The other terms are straightforward to obtain by continuing the Taylor expansion.

We proceed with a detailed derivation of the first-order correction to the current density. In first-order perturbation theory, the solution of Eq. (4) with an initial condition $\lim_{t\to-\infty}\alpha_{q\mathbf{k}}^{(n)}(t)\exp[\frac{\mathrm{i}}{\hbar}\epsilon_{q}(\mathbf{k})t]=\delta_{qn}$ is

[TABLE]

where we have introduced transition frequencies:

[TABLE]

From now on, we will omit, for brevity, $\mathbf{k}$ in the arguments of the transition frequencies and momentum matrix elements. Substituting Eq. (15) into Eq. (6), neglecting the terms quadratic with respect to $\mathbf{a}$, reintroducing $\mathbf{A}(t)$ by substituting $\mathbf{a}^{*}$ with $2\mathbf{A}(t)\mathrm{e}^{-\gamma t-\mathrm{i}\omega_{0}t}-\mathbf{a}\mathrm{e}^{-2\mathrm{i}\omega_{0}t}$, and taking the limit $\omega_{0}\to 0$, we obtain

[TABLE]

Before we take the limit $\gamma\to 0+$, we notice that those terms under the sum where $\omega_{in}=0$ vanish in this limit. In the rest of this paper, the notation $\sum_{i\neq n}$ will imply that not only the term $i=n$ but also all other terms where $\omega_{in}=0$, if present, must be omitted (otherwise, a zero would appear in the denominator). With this in mind, we take the limit $\gamma\to 0+$, use Eq. (5), and write the sum of the zero- and first-order approximations to the adiabatic current density as

[TABLE]

The first term on the right-hand side of this expression, representing the electric current that flows without external field, is identical with the zero-order expansion term of the intraband adiabatic current, Eq. (13). The two terms cancel when we evaluate $\Delta\mathbf{J}_{N}$ according to Eq. (11). In the rest of Eq. (17), we recognize the Thomas-Reiche-Kuhn sum rule in the form applicable to solids Luttinger and Kohn (1955). For our purposes, it is convenient to write this rule as

[TABLE]

where $\beta\in\{x,y,z\}$ is a Cartesian coordinate, and $\mathbf{e}_{\beta}$ is a unit vector along the corresponding coordinate axis. Remembering Eq. (14), we see that if the summation in Eq. (17) is performed over the entire infinite manifold of bands, then the requirement encoded in Eq. (9) leads to the Thomas-Reiche-Kuhn sum rule. If we truncate the basis set, we violate the sum rule.

Following Eq. (11), we obtain the following first-order correction to the electric current density obtained by solving the TDSE numerically with a finite number of bands:

[TABLE]

Here and in the following, summation is performed only over those bands that $\mathbf{J}_{N}(t)$ was evaluated with. As a reminder that the basis is truncated, we have replaced $\sum_{i\neq n}$ with $\sum_{i\neq n}^{N}$.

For a dielectric, where $\mathbf{J}_{\mathrm{ad}}(\mathbf{A})\equiv 0$, $f_{n}(\mathbf{k})=1$ for valence bands, and $f_{n}(\mathbf{k})=0$ for conduction bands, we can write Eq. (19) as

[TABLE]

Here, $N_{\mathrm{VB}}=\sum_{n\in\mathrm{occ}}f_{n}$ is the number of valence bands in a simulation to which the correction is applied; we denote this set of bands as “VB”. To emphasize that we consider the case where all the valence bands are initially occupied, while the conduction bands are empty, we have replaced $n\in\mathrm{occ}$ with $n\in\mathrm{VB}$. The summation over $i$ goes over all the valence and conduction bands, except band $n$ and those terms where $\omega_{in}(\mathbf{k})=0$; however, it is easy to see that the terms with $i\in\mathrm{VB}$ cancel.

For dielectrics with a crystal symmetry that demands that the electric current must flow along the laser polarization (e.g. cubic symmetry), Eq. (20) is equivalent to substituting the actual number of valence bands, $N_{\mathrm{VB}}$, with an effective one, $N_{\mathrm{VB}}^{\mathrm{eff}}$. Indeed, by multiplying both sides of Eq. (20) with $\mathbf{A}$ and relying on $\Delta\mathbf{J}^{(1)}\parallel\mathbf{A}$, we obtain

[TABLE]

where $\mathbf{e}_{\mathrm{L}}$ is a unit vector pointing along $\mathbf{A}$. Comparing this expression with Eq. (6), we see that our correction gives the same effect as substituting, in Eq. (6), $\mathbf{A}(t)$ with $\mathbf{A}(t)N_{\mathrm{VB}}^{\mathrm{eff}}/N_{\mathrm{VB}}$, where

[TABLE]

The concept of an effective number of valence electrons, $N_{\mathrm{VB}}^{\mathrm{eff}}$, is well known Philipp and Ehrenreich (1963); Smith and Shiles (1978); Wooten (2013), but it is not as widely known that $N_{\mathrm{VB}}^{\mathrm{eff}}$ can be used to adjust the polarization response obtained with an insufficient number of bands. Also, in the next section, we will demonstrate that the applicability of this concept is not limited to the linear response—Eq. (19) works surprisingly well even if the polarization response is very nonlinear.

Note that $N_{\mathrm{VB}}$ does not depend on the field magnitude; consequently, our corrections retain their importance in the weak-field limit. In other words, no matter how weak the electric field is, one has to use a sufficient number of bands to either satisfy the Thomas-Reiche-Kuhn sum rule or be able to compensate for its violation.

The derivation of higher-order corrections is similar to the one presented above, but it is much more laborious. Therefore, we only present the final results. For the second-order correction, we get

[TABLE]

This expression necessarily equals zero if $\hat{H}^{(0)}$ is symmetric with respect to time reversal and $f_{n}(-\mathbf{k})=f_{n}(\mathbf{k})$.

Finally, for the third-order correction, we obtain

[TABLE]

Even though it is a rather lengthy expression, we will demonstrate that it may be very useful.

IV A one-dimensional demonstration of the method

We derived our adiabatic corrections in three spatial dimensions; nevertheless, we demonstrate their power by means of one-dimensional simulations. On the one hand, it is much easier to control and achieve numerical convergence in one-dimensional calculations; on the other hand, the equations derived in the previous section become more manageable once we replace all vector quantities with scalar ones. Also, we further simplify the equations by considering a dielectric and eliminating $f_{n}(\mathbf{k})$. We write our corrections to the current density as an expansion in powers of the vector potential:

[TABLE]

In Eqs. (5), (19), (22), and (23), we replace $\mathrm{d}^{3}\mathbf{k}/(2\pi)^{3}$ with $\mathrm{d}k/(2\pi)$, substitute scalar quantities for the vector ones, and, after some simplifications, obtain

[TABLE]

[TABLE]

and

[TABLE]

These coefficients depend only on the stationary electronic structure of the solid, and they do not depend on the laser pulse. Also, evaluating the coefficients, we have dropped the adiabatic intraband current density, which is supposed to be zero for dielectrics but may acquire nonzero values when the integral in Eq. (10) is evaluated numerically. We found that these values were so small that accounting for them had a negligible effect on the results presented in this section.

In the following, we use atomic units ($\hbar=e=m_{0}=1$) unless specified otherwise. For our numerical demonstration, we employ a stationary Hamiltonian, $\hat{H}^{(0)}=-\frac{1}{2}\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}+V(x)$, with a noncentrosymmetric lattice potential:

[TABLE]

Here, $a=9.45\ \text{at.\,u.}=5$ Å is the lattice constant. Evaluating the band structure, we Fourier-expanded $V(x)$ into 81 plane waves. We regard the lowest two energy bands of this potential as the valence ones. The energy gap between the second and the third bands is equal to $\epsilon_{3}(0)-\epsilon_{2}(0)=9$ eV, which is close to the band gap of quartz. $\hat{H}^{(0)}$ is symmetric with respect to time reversal; therefore, up to round-off errors, $c_{2}$ in Eq. (26) evaluates to zero.

We numerically solved the TDSE for 61 uniformly spaced crystal momenta using only those stationary states that had an energy below a certain level. We measure this level from the bottom of the lowest conduction band and refer to it as the cut-off energy, $\epsilon_{\mathrm{cut}}$. In other words, the subset of eigenstates of $\hat{H}^{(0)}$ that we use to solve the TDSE consists of all the Bloch states with energies $\epsilon_{q}(k)\leq\epsilon_{\mathrm{cut}}+\epsilon_{3}(0)$.

We used the following expression for the vector potential of the external field:

[TABLE]

Here, $E_{0}$ is the amplitude of the electric field, $\omega_{0}$ is the central frequency of the laser pulse, $\theta(x)$ is the Heaviside step function, and $\tau_{\mathrm{L}}$ is related to the full width at half maximum (FWHM) of $A^{2}(t)$ by

[TABLE]

We chose the central wavelength of the laser pulse to be 750 nm ($\omega_{0}=2.51\ \text{fs}^{-1}=0.06$ at. u.; $\hbar\omega_{0}=1.65$ eV) and set the FWHM to 4 fs ($\tau_{\mathrm{L}}=7.66\ \text{fs}=317$ at. u.).

In Fig. 1, we compare electric current densities evaluated in different approximations for two peak values of the electric field: $E_{0}=0.1$ V/Å for panel (a) and $E_{0}=1$ V/Å for panel (b). The spectra of these current densities are compared in Fig. 2. We multiplied the current density with a constant factor (0.96333) that made the refractive index at the central laser wavelength equal to 1.47, which is the refractive index of thin-film fused silica at 750 nm Gao et al. (2013). This procedure allows us to label the $y$-axis of Fig. 1 with units of $e/(\text{nm}^{2}\,\text{fs})$.

The current density labeled as “accurate” was obtained with $\epsilon_{\mathrm{cut}}=2391.4$ eV, which corresponds to using 40 energy bands. We verified that, for this large value of $\epsilon_{\mathrm{cut}}$, corrections are unnecessary, so we did not apply them. The other results in Fig. 1 were evaluated for $\epsilon_{\mathrm{cut}}=25$ eV, which corresponds to using just the lowest three conduction bands of our one-dimensional model. Without corrections, the velocity gauge fails dramatically for this low value of $\epsilon_{\mathrm{cut}}$. However, already the first-order correction, specified by Eqs. (24) and (25), greatly improves the accuracy. In Fig. 1(a), the curves representing the accurate and corrected current densities are virtually indistinguishable. When we increase the peak electric field to $E_{0}=1$ V/Å, we see significant deviations between the accurate $J(t)$ (thick gray curve) and that obtained in the simulation with $\epsilon_{\mathrm{cut}}=25$ eV, followed by applying the first-order correction (solid blue curve). These deviations disappear as soon as we add the third-order correction to the first-order one (red dashed curve).

In the frequency domain, the first-order correction only affects the spectral range at the central laser frequency, while the third-order correction affects the spectral components at $\omega_{0}$ and $3\omega_{0}$. This is clearly seen in Fig. 2, where we also observe that high frequencies ($\omega\gtrsim 4\omega_{0}$) are much less demanding to the number of bands than the low frequencies are. Note that the current density in Figs. 1(b) and Figs. 2(b) shows signatures of highly nonlinear processes. In particular, the fast oscillations in Figs. 1(b) emerge due to the multiphoton absorption Korbman et al. (2013) (in our case, the ratio of the band gap to the photon energy is equal to 5.4). On average, $6.6\times 10^{-4}$ electrons per unit cell are excited to conduction bands by the pulse with $E_{0}=1$ V/Å; this value grows up to $1.3\times 10^{-2}$ for $E_{0}=1.5$ V/Å. For $t>\tau_{\mathrm{L}}=7.66\ \text{fs}$, the three curves obtained with $\epsilon_{\mathrm{cut}}=25$ eV coincide because the vector potential is zero, and so are the correction terms. The fact that these curves are very close to our most accurate current density means that a negligible number of electrons is excited to energies above $\epsilon_{\mathrm{cut}}$. This is expected because bands that carry residual population are responsible for processes more complex than those that our adiabatic corrections account for.

We quantify the discrepancy between an accurate current density, $J(t)$, and an approximate one, $J_{N}(t)$, using the following measure:

[TABLE]

Figure 3 shows how $\delta$ depends on $\epsilon_{\mathrm{cut}}$.

For the moderately strong peak field equal to $E_{0}=0.1$ V/Å, already the first-order correction decreases $\delta$ by two orders of magnitude in the range of cut-off energies between 25 eV and 200 eV. The third-order correction decreases the discrepancy even further. The improvement due to the analytical corrections is not as dramatic when the laser pulse has a peak field as large as $E_{0}=1$ V/Å, but we see that the third-order correction is particularly useful in this case. For example, without corrections, we get $\delta=0.0022$ at $\epsilon_{\mathrm{cut}}=200$ eV; with the first-order corrections, this level of discrepancy is first reached at $\epsilon_{\mathrm{cut}}=176$ eV, while the third-order correction provides this value of $\delta$ already at $\epsilon_{\mathrm{cut}}=83$ eV.

Figure 4 gives more details on how the discrepancy depends on the peak electric field. In this figure, we represent $\delta$ using false colors.

Panel (a) shows the result of applying the first-order correction. For panel (b), we corrected the current density using both $c_{1}$ and $c_{3}$. This figure confirms the usefulness of the third-order correction. In both panels, we observe the same general trend: a stronger laser field requires a larger value of $\epsilon_{\mathrm{cut}}$. However, when the third-order correction is used, the dependence is not monotonous. We also notice that the local minima of $\delta(\epsilon_{\mathrm{cut}})$ in Fig. 4(b) tend to maintain their positions as $E_{0}$ is increased.

V Conclusions and outlook

We have shown how to significantly reduce the number of bands required by velocity-gauge simulations performed in a stationary basis of Bloch states. The original motivation for this work was to speed up modeling of highly nonlinear ultrafast phenomena that take place when intense few-cycle laser pulses interact with crystalline solids. However, our method is very general, and it may find other applications. In particular, it can improve real-time methods of evaluating linear and nonlinear optical susceptibilities (the results presented here are only relevant to $\chi^{(1)}$, $\chi^{(2)}$, and $\chi^{(3)}$, but higher-order corrections can be derived as well).

Equations (19), (22), and (23) present the main results of this paper. The key idea of our approach was to consider the limit where the frequency of the field, $\omega_{0}$, approaches zero, while the amplitude of the vector potential is fixed. In this limit, the electric field becomes infinitesimally weak, which is a significant difference between our approach and that used to analyze the polarization divergence Aspnes (1972); Sipe and Ghahramani (1993), where the low-frequency limit was applied with a fixed amplitude of the electric field, which made the vector potential grow indefinitely.

Our approach thrives on the fact that the adiabatic approximation works well for those contributions to the polarization response that involve highly excited stationary states, where the excitation energy is much larger than $\hbar\omega_{0}$. Hence, our method works best when $\omega_{0}$ is small. We have demonstrated excellent performance for a near-infrared laser pulse interacting with a model dielectric.

Acknowledgements

The authors thank S. Yu. Kruchinin for his helpful remarks. Supported by the DFG Cluster of Excellence: Munich-Centre for Advanced Photonics. M. S. W. was supported by the International Max Planck Research School of Advanced Photon Science (IMPRS-APS).