Basis-independent spectral methods for non-linear optical response in arbitrary tight-binding models

2.1.1. Green's functions.

To use the Keldysh formalism for fermions, we need the definitions of the time-ordered, lesser, greater and anti-time-ordered Green's functions. Respectively,

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\rm i}G_{ab}^{T}(t,t') = \left\langle T\left[c_{a}(t)c_{b}^{\dagger}(t')\right]\right\rangle \label{eq:Green's functions definition} \nonumber \end{align} \tag{ 1 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\rm i}G_{ab}^{<}(t,t') = -\left\langle c_{b}^{\dagger}(t')c_{a}(t)\right\rangle \nonumber \end{align} \tag{ 2 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\rm i}G_{ab}^{>}(t,t') = \left\langle c_{a}(t)c_{b}^{\dagger}(t')\right\rangle \nonumber \end{align} \tag{ 3 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\rm i}G_{ab}^{\tilde{T}}(t,t') = \left\langle \widetilde{T}\left[c_{a}(t)c_{b}^{\dagger}(t')\right]\right\rangle . \nonumber \end{align} \tag{ 4 }$$

All the creation and annihilation operators are in the Heisenberg picture and the labels a and b denote states belonging to a complete single-particle basis. T is the time-ordering operator and $\tilde{T}$ the anti-time-ordering operator. The average $\left\langle \cdots\right\rangle$ stands for ${{\rm Tr}}\left[\rho(t_{0})\cdots\right]/{{\rm Tr}}\left[\rho(t_{0})\right]$ in the grand canonical ensemble, $\rho$ is the density matrix and t₀ denotes the time at which the external perturbation has been switched on. These are the building blocks of the Keldysh formalism. The advanced and retarded Green's functions are a simple combination of the previous objects:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G^{R} = G^{T}-G^{<}\label{eq:Relation RT}\nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G^{A} = -G^{\tilde{T}}+G^{<}.\label{eq:relation AT} \nonumber \end{align} \tag{ 6 }$$

The non-perturbed versions of these Green's functions are denoted by a lowercase g.

2.1.2. Expected value of an operator.

The expected value of the current $\boldsymbol{J}(t)$ (or any one-particle operator) may be evaluated with resort to these Green's functions by tracing over its product with the perturbed lesser Green's function:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{J}(t)=\left\langle \hat{\boldsymbol{J}}(t)\right\rangle =-{{\rm Tr}}\left[\hat{\boldsymbol{J}}(t){\rm i}G^{<}(t,t)\right]. \label{current in terms of G<} \nonumber \end{align} \tag{ 7 }$$

The Fourier transform¹ of $\boldsymbol{J}(t)$ is shown diagrammatically in figure 1. The circles stand for the full, perturbed operators in the presence of an external field.

Zoom In Zoom Out Reset image size

2.1.3. Conductivity.

We use the same definition for the nonlinear optical conductivity as in [8, 10]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & J^{\alpha}\left(\omega\right)=\sigma^{\alpha\beta}\left(\omega\right)E^{\beta}\left(\omega\right)+\int\frac{{{\rm d}}\omega_{1}}{2\pi}\int\frac{{{\rm d}}\omega_{2}}{2\pi} \label{conductivity definition}\nonumber \\ & \times \sigma^{\alpha\beta\gamma}\left(\omega_{1},\omega_{2}\right)E^{\beta}\left(\omega_{1}\right)E^{\gamma}\left(\omega_{2}\right)2\pi\delta\left(\omega_{1}+\omega_{2}-\omega\right)+\cdots\nonumber \end{align} \tag{ 8 }$$

where $E^{\alpha}$ is the component of the electric field along the $\alpha$ direction and the repeated indices are assumed to be summed over. The coefficients of this expansion are the conductivities at each order in the expansion. The next section is devoted to finding the perturbation expansion of G^<. In this paper we are dealing with tight-binding models, in which case the current operator will itself be a power series of the external field.

Our system is described by the many-particle time-dependent Hamiltonian

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{H}(t)=\mathcal{H}_{0}+\mathcal{H}_{{\rm ext}}(t) \nonumber \end{align} \tag{ 9 }$$

where $\mathcal{H}_{0}$ is an Hamiltonian that we can solve exactly and $\mathcal{H}_{{\rm ext}}(t)$ is the time-dependent external perturbation. Here we restrict ourselves to non-interacting Hamiltonians since we're dealing with non-interacting electrons. These operators are expressed in terms of their single-particle counterparts as

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mathcal{H}_{{\rm ext}}(t)=\sum_{ab}\left[H_{{\rm ext}}(t)\right]_{ab}c_{a}^{\dagger}(t)c_{b}(t) \nonumber \end{align} \tag{ 10 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mathcal{H}_{0}=\sum_{ab}\left[H_{0}\right]_{ab}c_{a}^{\dagger}(t)c_{b}(t). \nonumber \end{align} \tag{ 11 }$$

The expansion of the perturbed lesser Green's function G^< will be expressed in terms of the unperturbed Green's functions g^>, g^<, g^R and g^A in Fourier space:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm i}\tilde{g}^{<}\left(\omega\right) = -2\pi f\left(\hbar\omega\right)\delta\left(\omega-H_{0}/\hbar\right) \label{eq: lesser green's function} \nonumber \end{align} \tag{ 12 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm i}\tilde{g}^{>}\left(\omega\right) = 2\pi\left[1-f\left(\hbar\omega\right)\right]\delta\left(\omega-H_{0}/\hbar\right) \nonumber \end{align} \tag{ 13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm i}\tilde{g}^{R}\left(\omega\right) = \frac{{\rm i}}{\omega-H_{0}/\hbar+{\rm i}0^{+}} \nonumber \end{align} \tag{ 14 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm i}\tilde{g}^{A}\left(\omega\right) = \frac{{\rm i}}{\omega-H_{0}/\hbar-{\rm i}0^{+}}, \nonumber \end{align} \tag{ 15 }$$

where $\newcommand{\e}{{\rm e}} f\left(\epsilon\right)=\left(1+{\rm e}^{\beta(\epsilon-\mu)}\right){}^{-1}$ is the Fermi–Dirac distribution, $\beta$ is the inverse temperature and $\mu$ is the chemical potential. The Keldysh formalism and Langreth's rules provide the perturbation expansion of G^< [15]. Defining $V(t)=\left({\rm i}\hbar\right){}^{-1}H_{{\rm ext}}(t)$, the zeroth-order term in the expansion is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm i}\tilde{G}^{<(0)}(\omega)=\int{{\rm d}}\omega_{1}{\rm i}\tilde{g}^{<}(\omega_{1})\delta(\omega)\label{eq:green zeroth order frequencies} \nonumber \end{align} \tag{ 16 }$$

and the first-order one is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & {\rm i}\tilde{G}^{<(1)}(\omega)=\int\frac{{{\rm d}}^{3}\omega_{123}}{(2\pi)^{3}}(2\pi)^{2}\delta(\omega_{1}-\omega_{2}-\omega_{3}) \label{eq:green first order frequencies}\nonumber \\ & \times\delta(\omega+\omega_{3}-\omega_{1})\left[{\rm i}\tilde{g}^{R}(\omega_{1})\tilde{V}(\omega_{2}){\rm i}\tilde{g}^{<}(\omega_{3})+\right.\nonumber \\ & \left.+{\rm i}\tilde{g}^{<}(\omega_{1})\tilde{V}(\omega_{2}){\rm i}\tilde{g}^{A}(\omega_{3})\right].\nonumber \end{align} \tag{ 17 }$$

$\int{{\rm d}}^{n}\omega_{1\cdots n}$ is a shorthand for $\int\cdots\int{{\rm d}}\omega_{1}\cdots{{\rm d}}\omega_{n}$. The second-order term is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & {\rm i}\tilde{G}^{<(2)}(\omega)=\int\frac{{{\rm d}}^{5}\omega_{1\cdots5}}{(2\pi)^{5}}(2\pi)^{3}\delta(\omega_{5}+\omega-\omega_{1})\nonumber \\ & \times\delta(\omega_{1}-\omega_{2}-\omega_{3})\delta(\omega_{3}-\omega_{4}-\omega_{5})\nonumber \\ & \times \left[{\rm i}\tilde{g}^{R}(\omega_{1})\tilde{V}(\omega_{2}){\rm i}\tilde{g}^{R}(\omega_{3})\tilde{V}(\omega_{4}){\rm i}\tilde{g}^{<}(\omega_{5})\right.\nonumber \\ & +{\rm i}\tilde{g}^{R}(\omega_{1})\tilde{V}(\omega_{2}){\rm i}\tilde{g}^{<}(\omega_{3})\tilde{V}(\omega_{4}){\rm i}\tilde{g}^{A}(\omega_{5})\label{eq:second_order_green_frequencies-1}\nonumber \\ & \left.+{\rm i}\tilde{g}^{<}(\omega_{1})\tilde{V}(\omega_{2}){\rm i}\tilde{g}^{A}(\omega_{3})\tilde{V}(\omega_{4}){\rm i}\tilde{g}^{A}(\omega_{5})\right].\nonumber \end{align} \tag{ 18 }$$

Diagrammatically, the expansion of ${\rm i}G^{<}\left(\omega\right)$ is represented by figure 2. Each wavy line ending in a circle represents an external perturbation $\tilde{V}$. There are three different types of Green's functions that may appear in these expansions, with a certain regularity: a lesser Green's function $\tilde{g}^{<}$, which is always present, retarded Green's functions $\tilde{g}^{R}$ and advanced Green's functions $\tilde{g}^{A}$. Diagrammatically, $\tilde{g}^{<}$ is represented by a dashed line while the solid lines represent retarded or advanced Green's functions. To identify whether a line represents a retarded or advanced Green's function, one needs to read the diagram and identify the position of the lesser Green's function and the outgoing line. Reading clockwise (anti-clockwise) until finding the outgoing line, there can only be advanced (retarded) Green's functions. In each intersection, the corresponding external perturbation $\tilde{V}$ is inserted. An exception is made for the intersection with the line representing $\omega$, as it still needs to be contracted.

Zoom In Zoom Out Reset image size

If the external perturbation were a simple external field $\boldsymbol{E}(t)$, then the coupling would be $H_{{\rm ext}}(t)=e\boldsymbol{E}(t)\cdot\boldsymbol{r}$ and the previous expressions coupled with equation (7) would suffice. Now we will turn to tight-binding Hamiltonians, for which the external coupling is actually an infinite series of operators due to the way the electromagnetic field is introduced. This affects not only the $V$ operators but also the expression for the current operator.

Let us consider the following tight-binding Hamiltonian:

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mathcal{H}_{0}=\sum_{\boldsymbol{R}_{i},\boldsymbol{R}_{j}}\sum_{\sigma_{1},\sigma_{2}}t_{\sigma_{1}\sigma_{2}}\left(\boldsymbol{R}_{i},\boldsymbol{R}_{j}\right)c_{\sigma_{1}}^{\dagger}\left(\boldsymbol{R}_{i}\right)c_{\sigma_{2}}\left(\boldsymbol{R}_{j}\right).\label{eq:TB hamiltonian} \nonumber \end{align} \tag{ 19 }$$

The $\boldsymbol{R}_{i}$ represent the lattice sites and the $\sigma_{i}$ the other degrees of freedom unrelated to the position, such as the orbitals and spin. The electromagnetic field is introduced through Peierls' substitution:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle t_{\sigma_{1}\sigma_{2}}\left(\boldsymbol{R}_{i},\boldsymbol{R}_{j}\right)\rightarrow {\rm e}^{\frac{-{\rm i}e}{\hbar}\int_{\boldsymbol{R}_{j}}^{\boldsymbol{R}_{i}}\boldsymbol{A}(\boldsymbol{r}',t)\cdot {\rm d}\boldsymbol{r}'}t_{\sigma_{1}\sigma_{2}}\left(\boldsymbol{R}_{i},\boldsymbol{R}_{j}\right).\label{eq:peierls substitution} \nonumber \end{align} \tag{ 20 }$$

To introduce both a static magnetic field and a uniform electric field, we use the following vector potential:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{A}(\boldsymbol{r},t)=\boldsymbol{A}_{1}(\boldsymbol{r})+\boldsymbol{A}_{2}(t). \nonumber \end{align} \tag{ 21 }$$

The electric and magnetic fields are obtained from $\boldsymbol{E}(t)=-\partial_{t}\boldsymbol{A}_{2}(t)$ and $\boldsymbol{B}(\boldsymbol{r})=\boldsymbol{\nabla}\times\boldsymbol{A}_{1}(\boldsymbol{r})$. The introduction of the magnetic field only changes the $t_{\sigma_{1}\sigma_{2}}\left(\boldsymbol{R}_{i},\boldsymbol{R}_{j}\right)$ without introducing a time dependency. Therefore, we may assume that a magnetic field is always present without any loss of generality for the following discussion while keeping in mind that its introduction broke translation invariance. Since the magnetic field only affects the hopping parameters, from now on, the term in the vector potential that provides the electric field will be denoted by $\boldsymbol{A}(t)$. The external perturbation is obtained by expanding the exponential in equation (20) and identifying the original Hamiltonian.

3.1.1. Expansion of the external perturbation.

Expanding the exponential in equation (20) yields an infinite series of operators for the full Hamiltonian:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{\boldsymbol{A}}(t)=H_{0}+H_{{\rm ext}}(t) \nonumber \end{align} \tag{ 22 }$$

from which we identify, after a Fourier transform,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \tilde{V}\left(\omega\right)=\frac{e}{{\rm i}\hbar}h^{\alpha}\tilde{A}^{\alpha}\left(\omega\right)+\frac{{\rm e}^{2}}{{\rm i}\hbar}\frac{h^{\alpha\beta}}{2!}\int\frac{{{\rm d}}\omega'}{2\pi}\int\frac{{{\rm d}}\omega''}{2\pi} \quad\label{external perturbation}\nonumber \\ & \times \tilde{A}^{\alpha}\left(\omega'\right)\tilde{A}^{\beta}\left(\omega''\right)2\pi\delta\left(\omega'+\omega''-\omega\right)+\cdots.\nonumber \end{align} \tag{ 23 }$$

Repeated spatial indices are understood to be summed over. We have defined

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{h}^{\alpha_{1}\cdots\alpha_{n}}=\frac{1}{\left({\rm i}\hbar\right)^{n}}\left[\hat{r}^{\alpha_{1}},\left[\cdots\left[\hat{r}^{\alpha_{n}},H_{0}\right]\right]\right] \label{eq:generalized velocity} \nonumber \end{align} \tag{ 24 }$$

where $\hat{\boldsymbol{r}}$ is the position operator. In first order, $\hat{h}^{\alpha}$ is just the single-particle velocity operator. Under PBC, the position operator $\hat{\boldsymbol{r}}$ is ill-defined but its commutator with the Hamiltonian is not. In real space, this commutator is simply the Hamiltonian matrix element connecting the two sites i and j  multiplied by the distance vector $\boldsymbol{d}_{ij}$ between them. If we define this distance vector as the distance between neighbors instead of the difference of the two positions, it will be well defined in PBC. Using this strategy, all the $\hat{h}$ operators may be evaluated in position space by assigning to each bond the Hamiltonian matrix element multiplied by the required product of difference vectors $h_{ij}^{\alpha_{1}\cdots\alpha_{n}}=\left({\rm i}\hbar\right){}^{-n}H_{ij}d_{ij}^{\alpha_{1}}\cdots d_{ij}^{\alpha_{n}}$.

In figure 3, we see how the diagrammatic representation of the external perturbation unfolds into an infinite series of external fields. The wavy line represents $\tilde{\boldsymbol{A}}$ and the number of external fields connected to the same point is the number of commutators in equation (24).

Zoom In Zoom Out Reset image size

3.1.2. Expansion of the current.

The current operator is calculated directly from the Hamiltonian, using $\hat{J}^{\alpha}=-\Omega^{-1}\partial H/\partial A^{\alpha}$ ($\Omega$ is the volume of the sample), which also follows a series expansion due to the presence of an infinite number of $\boldsymbol{A}(t)$ in $H_{{\rm ext}}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \hat{J}^{\alpha}(t)=-\frac{e}{\Omega}\left(\hat{h}^{\alpha}+e\hat{h}^{\alpha\beta}A^{\beta}(t)\right.\nonumber \\ &\qquad \left.\quad\quad+\frac{{\rm e}^{2}}{2!}\hat{h}^{\alpha\beta\gamma}A^{\beta}(t)A^{\gamma}(t)+\cdots\right).\label{eq: current expansion} \nonumber \end{align} \tag{ 25 }$$

Figure 4 depicts the diagrammatic representation of this operator in Fourier space.

Zoom In Zoom Out Reset image size

The complexity of this expansion becomes clear. In equation (7), both the current operator and the Green's functions follow a perturbation expansion. Furthermore, each interaction operator in every one of the terms in the Green's function expansion also follows a similar expansion. We now have all the objects needed for the perturbative expansion of the conductivity.

In the previous sections we laid out the expressions for each individual operator in our expansion and represented their corresponding diagrammatic depictions. In this subsection, we put together all the elements of the previous sections to provide the full diagrammatic representation of the first and second-order conductivities. This expansion closely resembles that of [16] but has several differences due to the usage of these specific Green's functions. The only thing left to do is to replace the perturbed objects in the diagrammatic representation of the expected value of the current operator by their expansions. It is straightforward to see how the diagrams fit together in figure 5, which shows all the contributing diagrams up to second order.

Zoom In Zoom Out Reset image size

Obtaining the conductivity from the current is a matter of expressing the frequencies $\omega'$, $\omega''$ and $\omega'''$ in terms of $\omega_{1}$, $\omega_{2}$ and $\omega$ and using $\boldsymbol{E}\left(\omega\right)={\rm i}\omega\boldsymbol{A}\left(\omega\right)$. The Dirac delta in equation (8) simply means that $\omega$ is to be replaced by the sum of external frequencies entering the diagram. Thus, the nth order conductivity may be found using the following rules:

• 1.  
  Draw all the diagrams with n wavy lines coming in the diagram, one going out and one dashed interconnecting line. Integrate over the internal frequencies and ignore the conservation of momentum in the vertex containing $\omega$, as that is already taken into account by the Dirac delta in the definition of the conductivity.
• 2.  
  Reading clockwise starting from the vertex containing $\omega$, insert, by order, a generalized velocity operator $h^{\alpha_{1}\cdots\alpha_{k}}$ at each vertex and a Green's function at each edge. Each $\alpha_{i}$ is the label of a frequency line connecting to the vertex. If the edge is a dashed line, the Green's function is ig^<. All the edges before that correspond to ig^R and the ones after it to ig^A. Trace over the resulting operator.
• 3.  
  Multiply by $\Omega^{-1}{\rm e}^{n+1}\prod\nolimits_{k=1}^{n}\left({\rm i}\omega_{k}\right){}^{-1}\left({\rm i}\hbar\right){}^{1-N}$, where n is the number of dashed lines and N is the number of interconnecting lines. For each vertex, divide by the factorial of the number of outgoing lines.

Following these rules and replacing ig^< by equation (12), the first-order conductivity is found:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \sigma^{\alpha\beta}\left(\omega\right)=\frac{{\rm i}e^{2}}{\Omega\omega}\int_{-\infty}^{\infty}{\rm d}\epsilon f(\epsilon){{\rm Tr}}\left[\hat{h}^{\alpha\beta}\delta\left(\epsilon-H_{0}\right)+\frac{1}{\hbar}\hat{h}^{\alpha}g^{R}\right. \nonumber \\ & \left.\left(\epsilon/\hbar+\omega\right)\hat{h}^{\beta}\delta\left(\epsilon-H_{0}\right)+\frac{1}{\hbar}\hat{h}^{\alpha}\delta\left(\epsilon-H_{0}\right)\hat{h}^{\beta}g^{A}\left(\epsilon/\hbar-\omega\right)\right]. \label{cond 1st order} \nonumber \end{align} \tag{ 26 }$$

Similarly, for the second-order conductivity:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \sigma^{\alpha\beta\gamma}\left(\omega_{1},\omega_{2}\right)=\frac{1}{\Omega}\frac{{\rm e}^{3}}{\omega_{1}\omega_{2}}\int_{-\infty}^{\infty}{\rm d}\epsilon f(\epsilon){{\rm Tr}}\left[\frac{1}{2}\hat{h}^{\alpha\beta\gamma}\delta(\epsilon-H)\right.\nonumber \\ &\left.+\frac{1}{\hbar}\hat{h}^{\alpha\beta}g^{R}\left(\epsilon/\hbar+\omega_{2}\right)\hat{h}^{\gamma}\delta(\epsilon-H)\right. \label{eq: cond 2nd order}\nonumber \\ & +\frac{1}{\hbar}\hat{h}^{\alpha\beta}\delta(\epsilon-H)\hat{h}^{\gamma}g^{A}\left(\epsilon/\hbar-\omega_{2}\right)+\frac{1}{2\hbar}\hat{h}^{\alpha}g^{R}\left(\epsilon/\hbar+\omega_{1}+\omega_{2}\right)\hat{h}^{\beta\gamma}\nonumber \\ &\delta(\epsilon-H)+\frac{1}{2\hbar}\hat{h}^{\alpha}\delta(\epsilon-H)\hat{h}^{\beta\gamma}g^{A}\left(\epsilon/\hbar-\omega_{1}-\omega_{2}\right)\nonumber \\ & +\frac{1}{\hbar^{2}}\hat{h}^{\alpha}g^{R}\left(\epsilon/\hbar+\omega_{1}+\omega_{2}\right)\hat{h}^{\beta}g^{R}\left(\epsilon/\hbar+\omega_{2}\right)\hat{h}^{\gamma}\delta(\epsilon-H)\nonumber \\ &+\frac{1}{\hbar^{2}}\hat{h}^{\alpha}g^{R}\left(\epsilon/\hbar+\omega_{1}\right)\hat{h}^{\beta}\delta(\epsilon-H)\hat{h}^{\gamma}g^{A}\left(\epsilon/\hbar-\omega_{2}\right)\nonumber \\ & \left.+\frac{1}{\hbar^{2}}\hat{h}^{\alpha}\delta(\epsilon-H)\hat{h}^{\beta}g^{A}\left(\epsilon/\hbar-\omega_{1}\right)\hat{h}^{\gamma}g^{A}\left(\epsilon/\hbar-\omega_{1}-\omega_{2}\right)\right].\nonumber \end{align} \tag{ 27 }$$

The procedure is exactly the same for the nth order conductivity, which will have $2^{n-1}(n+2)$ diagrams. The higher-order expansions will not be obtained because a realistic computation of physical quantities with those formulas using spectral methods would require tremendous computational power. This point will be further explained in the next section. Despite these limitations, we have provided the required framework for obtaining those formulas, which might be useful in the future.

The Chebyshev polynomials of the first kind are a set of orthogonal polynomials defined in the range $\left[-1,1\right]$ by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_{n}(x)=\cos\left(n\arccos(x)\right). \nonumber \end{align} \tag{ 28 }$$

They satisfy a recursion relation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_{0}(x) = 1 \nonumber \end{align} \tag{ 29 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_{1}(x) = x \nonumber \end{align} \tag{ 30 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_{n+1}(x) = 2xT_{n}(x)-T_{n-1}(x) \nonumber \end{align} \tag{ 31 }$$

and the following orthogonality relation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{-1}^{1}T_{n}(x)T_{m}(x)\frac{{{\rm d}}x}{\sqrt{1-x^{2}}}=\delta_{nm}\frac{1+\delta_{n0}}{2}. \nonumber \end{align} \tag{ 32 }$$

These polynomials may be used to expand functions of the Hamiltonian provided its spectrum has been scaled by a factor $\Delta_{E}$ to fit in the range $\left[-1,1\right]$. The only functions of the Hamiltonian appearing in the expansion are Dirac deltas and Green's functions, which have the following expansions in terms of Chebyshev polynomials:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta(\epsilon-H_{0})=\sum_{n=0}^{\infty}\Delta_{n}(\epsilon)\frac{T_{n}(H_{0})}{1+\delta_{n0}}\label{eq: delta function cheb expansion} \nonumber \end{align} \tag{ 33 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g^{\sigma,\lambda}(\epsilon,H_{0})=\frac{\hbar}{\epsilon-H_{0}+{\rm i}\sigma\lambda}=\hbar\sum_{n=0}^{\infty}g_{n}^{\sigma,\lambda}(\epsilon)\frac{T_{n}(H_{0})}{1+\delta_{n0}}\label{eq: green function cheb expansion} \nonumber \end{align} \tag{ 34 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta_{n}(\epsilon)=\frac{2T_{n}(\epsilon)}{\pi\sqrt{1-\epsilon^{2}}} \nonumber \end{align} \tag{ 35 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g_{n}^{\sigma,\lambda}\left(\epsilon\right)=-2\sigma {\rm i}\frac{{\rm e}^{-n{\rm i}\sigma\arccos\left(\epsilon+{\rm i}\sigma\lambda\right)}}{\sqrt{1-\left(\epsilon+{\rm i}\sigma\lambda\right)^{2}}}. \nonumber \end{align} \tag{ 36 }$$

The function $g^{\sigma,\lambda}$ encompasses both retarded and advanced Green's functions in the limit $\lambda\rightarrow0^{+}$: $g^{+,0^{+}}$ is the retarded Green's function and $g^{-,0^{+}}$ the advanced one. The operator part has been completely separated from its other arguments. All the Dirac deltas and Green's functions may therefore be separated into a polynomial of H₀ and a coefficient which encapsulates the frequency and energy parameters. The trace in the conductivity now becomes a trace over a product of polynomials and $\hat{h}$ operators, which can be encapsulated in a new object, the $\Gamma$ matrix:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Gamma_{n_{1}\cdots n_{m}}^{\boldsymbol{\alpha}_{1},\cdots,\boldsymbol{\alpha}_{m}}=\frac{{{\rm Tr}}}{N}\left[\tilde{h}^{\boldsymbol{\alpha}_{1}}\frac{T_{n_{1}}(H_{0})}{1+\delta_{n_{1}0}}\cdots\tilde{h}^{\boldsymbol{\alpha}_{m}}\frac{T_{n_{m}}(H_{0})}{1+\delta_{n_{m}0}}\right]. \nonumber \end{align} \tag{ 37 }$$

The upper indices in bold stand for any number of indices: $\boldsymbol{\alpha}_{1}=\alpha_{1}^{1}\alpha_{1}^{2}\cdots\alpha_{1}^{N_{1}}$. Here we have used $\tilde{h}^{\boldsymbol{\alpha}_{1}}=\left({\rm i}\hbar\right){}^{N_{1}}\hat{h}^{\boldsymbol{\alpha}_{1}}$ rather than $\hat{h}$ to avoid using complex numbers when the Hamiltonian matrix is purely real in our numerical simulations. It's very important to keep in mind that these new operators are no longer hermitian. The commas in $\Gamma$ separate the various $\tilde{h}$ operators. N is the number of unit cells in the sample being studied and ensures that $\Gamma$ is an intensive quantity. Some examples:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Gamma_{nm}^{\alpha,\beta\gamma} = \frac{{{\rm Tr}}}{N}\left[\tilde{h}^{\alpha}\frac{T_{n}(H_{0})}{1+\delta_{n0}}\tilde{h}^{\beta\gamma}\frac{T_{m}(H_{0})}{1+\delta_{m0}}\right] \nonumber \end{align} \tag{ 38 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Gamma_{n}^{\alpha\beta} = \frac{{{\rm Tr}}}{N}\left[\tilde{h}^{\alpha\beta}\frac{T_{n}(H_{0})}{1+\delta_{n0}}\right] \nonumber \end{align} \tag{ 39 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Gamma_{nmp}^{\alpha,\beta,\gamma} = \frac{{{\rm Tr}}}{N}\left[\tilde{h}^{\alpha}\frac{T_{n}(H_{0})}{1+\delta_{n0}}\tilde{h}^{\beta}\frac{T_{m}(H_{0})}{1+\delta_{m0}}\tilde{h}^{\gamma}\frac{T_{p}(H_{0})}{1+\delta_{p0}}\right]. \nonumber \end{align} \tag{ 40 }$$

The $\Gamma$ matrix only depends on the physical system itself as it is merely a function of the Hamiltonian and the $\tilde{h}$ operators. The coefficients of the Chebyshev expansion may similarly be aggregated into a matrix, which we denote by $\Lambda$. Some examples:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Lambda_{n} = \int_{-\infty}^{\infty}{\rm d}\epsilon f(\epsilon)\Delta_{n}(\epsilon) \nonumber \end{align} \tag{ 41 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Lambda_{nm}\left(\omega\right) &= \hbar\int_{-\infty}^{\infty}{\rm d}\epsilon f(\epsilon)\left[g_{n}^{R}\left(\epsilon/\hbar+\omega\right)\Delta_{m}\left(\epsilon\right)\right.\nonumber \\ & \left.\quad\quad\quad\quad+\Delta_{n}\left(\epsilon\right)g_{m}^{A}\left(\epsilon/\hbar-\omega\right)\right] \nonumber \end{align} \tag{ 42 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \Lambda_{nmp}\left(\omega_{1},\omega_{2}\right)\nonumber \\ &= \hbar^{2}\int_{-\infty}^{\infty}{\rm d}\epsilon f(\epsilon)\left[g_{n}^{R}\left(\epsilon/\hbar+\omega_{1}+\omega_{2}\right)g_{m}^{R}\left(\epsilon/\hbar+\omega_{2}\right)\Delta_{p}\left(\epsilon\right)\right.\nonumber \\ & \quad\quad\quad\,\,\,\,\,\quad+g_{n}^{R}\left(\epsilon/\hbar+\omega_{1}\right)\Delta_{m}\left(\epsilon\right)g_{p}^{A}\left(\epsilon/\hbar-\omega_{2}\right)\nonumber \\ & \left.\quad\quad\quad\,\,\,\,\,\quad+\Delta_{n}\left(\epsilon\right)g_{m}^{A}\left(\epsilon/\hbar-\omega_{1}\right)g_{p}^{A}\left(\epsilon/\hbar-\omega_{1}-\omega_{2}\right)\right]\nonumber .\nonumber \end{align} \tag{ 43 }$$

In terms of these new objects, the conductivities become

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma^{\alpha\beta}\left(\omega\right)=\frac{-{\rm i}e^{2}}{\Omega_{c}\hbar^{2}\omega}\left[\sum_{n}\Gamma_{n}^{\alpha\beta}\Lambda_{n}+\sum_{nm}\Lambda_{nm}\left(\omega\right)\Gamma_{nm}^{\alpha,\beta}\right]\label{eq:conductivity_1st_order} \nonumber \end{align} \tag{ 44 }$$

in first order and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \sigma^{\alpha\beta\gamma}\left(\omega_{1},\omega_{2}\right)=\frac{{\rm i}e^{3}}{\Omega_{c}\omega_{1}\omega_{2}\hbar^{3}}\left[\frac{1}{2}\sum_{n}\Lambda_{n}\Gamma_{n}^{\alpha\beta\gamma}\right.\nonumber \\ & +\sum_{nm}\Lambda_{nm}\left(\omega_{2}\right)\Gamma_{nm}^{\alpha\beta,\gamma}+\frac{1}{2}\sum_{nm}\Lambda_{nm}\left(\omega_{1}+\omega_{2}\right)\Gamma_{nm}^{\alpha,\beta\gamma}\nonumber \\ & \left.+\sum_{nmp}\Lambda_{nmp}\left(\omega_{1},\omega_{2}\right)\Gamma_{nmp}^{\alpha,\beta,\gamma}\right] \nonumber \end{align} \tag{ 45 }$$

in second order. $\Omega_{c}$ is the volume of the unit cell.

Naturally, one cannot expect to sum the entire Chebyshev series, so it has to be truncated at a certain number of polynomials $N_{{\rm max}}$. Each of the entries in a $\Gamma$ matrix represents a complex number. Numerically, this is represented as two double-precision floating-point numbers, each taking up 8 bytes of storage. The amount of storage needed to store a $\Gamma$ matrix of dimension n is 16$N_{{\rm max}}^{n}$. The number of Chebyshev polynomials needed to obtain a decent resolution depends heavily on the problem at hand, but a typical number may be $N_{{\rm max}}=1024$. A one-dimensional $\Gamma$ matrix would take up 16 KiB of storage, a two-dimensional matrix 16 MiB and a three-dimensional matrix 16 GiB. Three-dimensional matrices appear in the second-order conductivity. The third-order conductivity would require a four-dimensional matrix and as such, 16 TiB of storage. Numbers like these make it unrealistic to go beyond second order conductivity.

Let a be the distance between consecutive atoms in the honeycomb lattice. Then, the primitive vectors between unit cells are (see figure 6)

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{a}_{1} &= a\left(\sqrt{3},0\right)\nonumber \\ \boldsymbol{a}_{2} &= a\left(\frac{\sqrt{3}}{2},\frac{3}{2}\right) \nonumber \end{align*}$$

and the distance vectors between nearest neighbors are

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\delta}_{1} &= \frac{a}{2}\left(\sqrt{3},-1\right)\nonumber \\ \boldsymbol{\delta}_{2} &= a(0,1)\nonumber \\ \boldsymbol{\delta}_{3} &= \frac{a}{2}\left(-\sqrt{3},-1\right). \nonumber \end{align*}$$

Zoom In Zoom Out Reset image size

The area of the unit cell is $\Omega_{c}=\frac{3\sqrt{3}}{2}a^{2}$. Starting from equation (19), the graphene Hamiltonian is obtained by invoking translational invariance of the unit cell $t_{\mu\nu}\left(\boldsymbol{R}_{m},\boldsymbol{R}_{n}\right)\,=$$t_{\mu\nu}\left(\boldsymbol{R}_{m}-\boldsymbol{R}_{n}\right)$ and

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle t_{AB}\left(\boldsymbol{\delta}_{1}\right)=t_{AB}\left(\boldsymbol{\delta}_{2}\right)=t_{AB}\left(\boldsymbol{\delta}_{3}\right)=-t. \nonumber \end{align*}$$

The remaining non-zero hopping integrals are found by using $t_{AB}=t_{BA}$. The on-site energies $t_{AA}\left(\boldsymbol{0}\right)$ and $t_{BB}\left(\boldsymbol{0}\right)$ are taken to be zero. A factor of two is included due to spin degeneracy.

These parameters were used to obtain the first-order optical conductivity for graphene, as seen in figure 7. We used a lattice with 4096 unit cells in each direction and 2048 Chebyshev moments in the expansion. The resulting plot is compared to the results obtained in [8] through $\boldsymbol{k}$-space integration of a translation-invariant system. The curves are indistinguishable.

Zoom In Zoom Out Reset image size

The only difference relative to regular graphene is found in the on-site energies. Let $t_{AA}\left(\boldsymbol{0}\right)=\Delta/2$ and $t_{BB}\left(\boldsymbol{0}\right)=-\Delta/2$. With the opening of a sizeable gap, the effects of excitons become relevant. In this simplified model, we ignore these effects to focus on the intrinsic second-order respose of the material. In these conditions, the one- and three-dimensional $\Gamma$ matrices are identically zero², so the second-order conductivity may be calculated resorting only to two-dimensional $\Gamma$ matrices. The calculation is thus simplified tremendously because the second-order conductivity reduces to

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle & \sigma^{\alpha\beta\gamma}\left(\omega_{1},\omega_{2}\right)=\frac{{\rm i}e^{3}}{\Omega_{c}\omega_{1}\omega_{2}\hbar^{3}}\nonumber \\ & \times \left[\sum_{nm}\Lambda_{nm}\left(\omega_{2}\right)\Gamma_{nm}^{\alpha\beta,\gamma}+\frac{1}{2}\sum_{nm}\Lambda_{nm}\left(\omega_{1}+\omega_{2}\right)\Gamma_{nm}^{\alpha,\beta\gamma}\right]. \nonumber \end{align*}$$

The indices $n,m$ are understood to be summed over. The photogalvanic effect may be reproduced from this formula by setting $\omega_{1}=\omega=-\omega_{2}$ and the numerical results are shown in figure 8. Again, we used 4096 unit cells in each lattice direction and 2048 Chebyshev moments and compare the results with the ones obtained by integrating in $\boldsymbol{k}$-space, just like in the previous subsection. For convenience, we define the constant $\sigma_{2}={\rm e}^{3}a/4t\hbar$ [21] in terms of the hopping integral t and the lattice parameter a.

Zoom In Zoom Out Reset image size

This particular example benefits considerably from the cancellation of the most complicated objects that needed to be calculated. In appendix A.1 we present an example with less symmetry that confirms the complete agreement between our method and the $\boldsymbol{k}$-space formalism.

Our formalism does not rely on translation invariance, and so may be used to study disordered systems. To show this, we now introduce to gapped graphene a simple model for disorder by letting each atomic site have a random local energy taken from a uniform distribution $\left[-W/2,W/2\right]$ (Anderson disorder [22]):

$$\begin{align*} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mathcal{H}_{W}=\sum_{\boldsymbol{R}}\sum_{\sigma}W_{\sigma}\left(\boldsymbol{R}\right)c_{\sigma}^{\dagger}\left(\boldsymbol{R}\right)c_{\sigma}\left(\boldsymbol{R}\right) \nonumber \end{align*}$$

where $\boldsymbol{R}$ is the position of the unit cell and $\sigma$ labels the atoms inside each unit cell. The presence of disorder is expected to smooth out the sharp features of the optical response. As disorder increases, we should see a decrease in conductivity due to Anderson localization. This is the exact behavior that is seen in figure 9 where we plot the photogalvanic effect in gapped graphene in the presence of Anderson disorder of varying strength. Some fluctuations exist at the features, which are expected to disappear as the system size approaches the thermodynamic limit.

Zoom In Zoom Out Reset image size

As expected, the introduction of Anderson disorder produces a broadening of the sharp features of the nonlinear optical conductivity. This broadening also means that there will be a larger response to the external electric field at frequencies smaller than the gap.

The large oscillations near the origin reveal something interesting about the numerical details of our formalism. Equation (27) is comprised of a complicated sum of several terms. Individually, some of these terms may be very large, but there may be cancellations among them. For each of these terms, the Chebyshev expansion is exact in the limit of infinite polynomials. For a finite number of polynomials, there will be slight differences between the exact result and the expansion, and if the exact result is very large, this difference will be considerable. It is highly unlikely that this difference will be the same for each term, and so their sum may not cancel out in the end. This is the typical behavior at the lower-frequency regime that is related to the singularities that plagued the velocity gauge approach. This has been discussed since the early work of Sipe and challenged for a long time the equivalence between the velocity and length gauges [4, 8, 10]. This effect could be fully mitigated by greatly increasing the number of polynomials, but here we are interested in the finite frequency behavior.

In realistic samples, vacancies and impurities may exist due to imperfections in the fabrication process, as well as other more complex structural defects. In this section, we show that our method allows us to obtain the second-harmonic generation of a system with structural disorder. Using equation (27), we show in figure 10 the effect of vacancies of varying concentration in the SHG of gapped graphene. Unlike Anderson disorder, the addition of vacancies to the system does not change the gap. Their most noticeable effect is to flatten the features of the second-harmonic generation. As discussed in the previous section, the lower frequencies are dominated by oscillations and would require many more polynomials to fully converge. Therefore, we omit that region and only represent the remaining regions, which have already converged within the desired accuracy.

Zoom In Zoom Out Reset image size

In this section, we briefly discuss some convergence properties of our method. For a more thorough discussion, see [17]. The convergence to the exact value depends on several factors:

• 1.  
  Spectral methods rely on the self-averaging properties of random vectors, yielding an associated variance. The error bar decreases as $1/\sqrt{N_{R}N}$, where N_R is the number of random vectors and N is the size of the sample.
• 2.  
  In the thermodynamic limit of an infinite lattice, the spectrum becomes continuous and so we expect the conductivity curve to be smooth. However, the systems used in simulations are finite and so have a typical energy level spacing, which we denote by $\delta\varepsilon$. This has important consequences for the resolution. Details characterized by a smaller energy scale than that of $\delta\varepsilon$ are meaningless because they cannot be distinguished from the contribution of individual energy levels. The maximum resolution is therefore limited by the energy level spacing. For our concrete examples with the honeycomb lattice, we use $\delta\varepsilon=3\pi t/L$, the energy level spacing at the Dirac point in graphene for a system of linear dimension L.
• 3.  
  The resolution may be controlled through $\lambda$, the broadening parameter of the Green's functions. Energy differences smaller than $\lambda$ become indistinguishable from one another. On the one hand, a small $\lambda$ is required in order to resolve the sharp features of the curve accurately. On the other hand, when $\lambda\lesssim \delta\varepsilon$, the discrete nature of the spectrum starts to become visible through the roughness of the curve. For sufficiently small $\lambda$, the expected sharp features of the curve become indistinguishable from the contributions of the individual energy levels. If these issues are not solved, they become a major source of systematic error in the final results. Therefore, if we want to see the expected thermodynamic limit, we have to ensure $\lambda\gtrsim\delta\varepsilon$.

In figure 11, the yy optical conductivity of graphene is represented for several values of $\lambda$. In this example, $\delta\varepsilon=5.3~{\rm meV}$. As $\lambda$ is decreased, the curve becomes sharper, but when $\lambda=2.3~{\rm meV}$ the discreteness of the spectrum starts to become noticeable through the roughness of the curve. It is starting to diverge from the expected smooth curve of the thermodynamic limit.

Zoom In Zoom Out Reset image size

In the lower inset, we study the convergence as a function of the number of polynomials at $\hbar\omega=4.66~{\rm eV}$, a region of rapidly changing conductivity. The smaller the $\lambda$, the more polynomials are required in order to have a fully converged result. Within the accuracy $\delta\sigma/\sigma_{0}\simeq0.1$, all the curves have already converged at $1.6\times10^{4}$ polynomials. These calculations were repeated for several different initial random vectors. In the plot we show only one of these calculations. The error bar associated with the random vectors is too small to be distinguished from the curves themselves.

In the upper inset, we do the same thing, but now in a very small region around $\hbar\omega=2.33~{\rm eV}$, a region of slowly increasing conductivity. The plot shows three sets of curves with different colors. Inside each set, we represent a collection of frequencies, ranging from $\hbar\omega=2.3300~{\rm eV}$ to $\hbar\omega=2.3316~{\rm eV}$. The darker curves correspond to higher frequencies. The main graph shows that all these curves have converged to the same value in a region of slowly increasing conductivity. The inset, however, shows a different picture. The red ($\lambda=23~{\rm meV}$) and black ($\lambda=230~{\rm meV}$) sets of curves show a variation consistent with the expected increasing conductivity. If one zooms in to those sets of curves, it is possible to check that they are indeed increasing in value as $\omega$ increases. The green curve ($\lambda=2.3~{\rm meV}$) is not only changing in a scale much larger than expected, but it is also decreasing. This variation comes from the individual contribution of the energy levels, not from features of the conductivity and is therefore artificial. Within the accuracy $\delta\sigma/\sigma_{0}\simeq10^{-3}$ each of these curves has completely converged at $1.6\times10^{4}$ polynomials but this level of accuracy is meaningless for $\lambda=2.3~{\rm meV}$. The error bars are not shown for clarity, but their values are the following: at $\lambda=230~{\rm meV}$, $\delta\sigma/\sigma_{0}=10^{-3}$; at $\lambda=23~{\rm meV}$, $\delta\sigma/\sigma_{0}=3\times10^{-3}$; at $\lambda=2.3~{\rm meV}$, $\delta\sigma/\sigma_{0}=5\times10^{-3}$. At this scale, the error bars are comparable to the variation due to the number of polynomials and to the value of $\lambda$.

These frequencies were chosen to compare the conductivity in a place where it is expected to converge quickly and another where it is expected to converge slowly. Looking at these graphs, it is possible to estimate how many polynomials are required to converge to the final value of the conductivity for the specified parameter $\lambda$ within a given accuracy. A rough estimate of the scaling is given by $N\sim\lambda^{-1}$.

In sum, given a fixed resolution $\lambda$, the number of polynomials should be large enough to ensure that the curves have converged, and the system size L should be large enough to ensure that the discreteness of the spectrum cannot be seen.

A similar analysis may be done for the second-order conductivity. We will not present it here for two reasons. Firstly, the main points of the previous paragraphs remain the same. Secondly, we cannot do such an analysis because the computational cost would be tremendously higher.

Our formalism provides a very general framework with which to compute the nonlinear optical response up to any order. Once the formulas were obtained, we chose to use spectral methods to perform the computation. This is not always the most efficient approach: for systems with translation invariance and periodic boundary conditions, we can specify the formulas for $\boldsymbol{k}$-space and then perform the explicit integration. Then, for a given set of parameters (temperature, broadening, Fermi energy) the computation time will scale as $L^{D}N_{\omega}$ where L^D is the number of points in the Brillouin zone (which is also the number of lattice sites), D is the dimensionality and $N_{\omega}$ is the number of frequencies we want to compute. For each $\boldsymbol{k}$ and each frequency, this method comes down to diagonalizing the Bloch Hamiltonian $H_{\boldsymbol{k}}$, and then summing over the whole set of $\boldsymbol{k}$ points. This method is extremely efficient at computing the optical conductivity at any order using the velocity gauge.

Using spectral methods, the computation is split into the calculation of the Chebyshev moments and the final matrix product of the $\Gamma$ matrices with the $\Lambda$ matrices. The first part is the most demanding and is independent of the parameters mentioned above. Its computation time scales as $L^{D}N^{n+1}$, where n is the order of the conductivity and N is the number of Chebyshev polynomials. More concretely, if we want to calculate the conductivity for a certain $\lambda$, using $N\sim\lambda^{-1}\sim N_{\omega}$, we find that the $\boldsymbol{k}$-space calculation scales much more favorably.

If the system has no translation invariance, $\boldsymbol{k}$-space integration is no longer useful and we would need to numerically diagonalize the full Hamiltonian. This method scales as $L^{3D}N_{\omega}$ which is highly unfavorable and because of that we would be limited to very small systems. In this context, spectral methods become the preferred choice.

For the examples used in this paper, the computation of the second-order conductivity with the $\boldsymbol{k}$-space formalism in a system with L  =  2048 took around 2 min on a Xeon E5-2650 with 16 threads. In comparison, the same computation took 3 h for translation-invariant gapped graphene with 2048 polynomials, and 70 h for gapped graphene with Anderson disorder/vacancies and 512 polynomials. Despite the discrepancy in computational efficiency, we know of no other more efficient way to compute the nonlinear optical conductivity for disordered systems.

The calculation of the photogalvanic effect for gapped graphene was very efficient due to the cancellation of the three-dimensional $\Gamma$ matrices. In this appendix, we provide an extra example, which does not benefit from that property. By changing the relative position of the two sublattices, we are able to obtain non-zero values in all the $\Gamma$ matrices, which enables us to test the remainder of the formula. All the hopping parameters in this system are exactly the same as in regular gapped graphene. The only difference is in the distance between atoms, which changes the velocity operators while keeping the Hamiltonian intact (See figure A1).

Zoom In Zoom Out Reset image size

The primitive vectors are identical, but the nearest-neighbor vectors are different:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\delta}_{1} &= a\left(\frac{\sqrt{3}}{2},-1\right)\nonumber \\ \boldsymbol{\delta}_{2} &= a\left(0,\frac{1}{2}\right)\nonumber \\ \boldsymbol{\delta}_{3} &= a\left(-\frac{\sqrt{3}}{2},-1\right). \nonumber \end{align*}$$

One of the sublattices was translated in the y  direction by a/2. The second-order xxx conductivity remains zero, but now the xxy photogalvanic effect is no longer zero and can be seen in figure A2. The lattice size and number of polynomials used was reduced to 1024 and 512 respectively, due to the greatly increased computational cost. At lower frequencies, the results start to diverge because there are not enough polynomials to resolve this region. The results are in great agreement with the ones obtained by $\boldsymbol{k}$-space integration. The small oscillations in the imaginary part are expected to disappear as the number of polynomials is increased.

Zoom In Zoom Out Reset image size