Gauge covariances and nonlinear optical responses

TL;DR

This paper develops a covariant formalism for calculating nonlinear optical responses in crystals using the reduced density matrix approach, comparing length and velocity gauges, and analyzing the effects of band truncation.

Contribution

It introduces a covariant derivative representation of the position operator, providing compact expressions for nonlinear responses in both gauges and analyzing gauge equivalence under band truncation.

Findings

Derived compact expressions for first and third order nonlinear responses in graphene.

Showed that gauge equivalence is broken by band truncation and how to properly implement truncation.

Compared length and velocity gauge formalisms, highlighting their formal similarities and differences.

Abstract

The formalism of the reduced density matrix is pursued in both length and velocity gauges of the perturbation to the crystal Hamiltonian. The covariant derivative is introduced as a convenient representation of the position operator. This allow us to write compact expressions for the reduced density matrix in any order of the perturbation which simplifies the calculations of nonlinear optical responses; as an example, we compute the first and third order contributions of the monolayer graphene. Expressions obtained in both gauges share the same formal structure, allowing a comparison of the effects of truncation to a finite set of bands. This truncation breaks the equivalence between the two approaches: its proper implementation can be done directly in the expressions derived in the length gauge, but require a revision of the equations of motion of the reduced density matrix in the velocity gauge.

Equations

$$V(\mathbf{r})=-e\left(V_{\mathbf{q}}{e}^{i\mathbf{q}\cdot \mathbf{r}}+V_{\mathbf{q}}^{\ast }{e}^{-i\mathbf{q}\cdot \mathbf{r}}\right),$$

$${\psi }_{\mathbf{k}s}\to e^{i{\theta }_s(\mathbf{k})}{\psi }_{\mathbf{k}s},$$

$${\mathcal{A}}_{\mathbf{k}{\mathbf{k}}^{\prime }s{s}^{\prime }}\to e^{-i({\theta }_s(\mathbf{k})-{\theta }_{s^{\prime }}({\mathbf{k}}^{\prime }))}{\mathcal{A}}_{\mathbf{k}{\mathbf{k}}^{\prime }s{s}^{\prime }}.$$

$$H_E(t)=\int d^d\mathbf{r}{\Psi }^{†}(\mathbf{r})\left[{\mathcal{H}}_0\left(\mathbf{r},\frac{\nabla }{i}\right)+e\mathbf{E}(t)\cdot \mathbf{r}\right]\Psi (\mathbf{r}),$$

$$H_A(t)=\int d^d\mathbf{r}{\Psi }^{†}(\mathbf{r})\left[{\mathcal{H}}_0\left(\mathbf{r},\frac{\nabla }{i}+\frac{e}{\hslash }\mathbf{A}(t)\right)\right]\Psi (\mathbf{r}).$$

$$i\hslash \frac{\partial |\psi (t)⟩}{\partial t}$$

$$|\bar{\psi }(t)⟩$$

$$i\hslash \frac{\partial |\bar{\psi }(t)⟩}{\partial t}=\left[\mathcal{U}(t)H_A(t){\mathcal{U}}^{†}(t)+i\hslash \frac{d\mathcal{U}(t)}{dt}{\mathcal{U}}^{†}(t)\right]|\bar{\psi }(t)⟩.$$

$$\mathcal{U}(t)=\exp \left[i\frac{e}{\hslash }\int d^d\mathbf{r}\mathbf{A}(t)\cdot \mathbf{r}\rho (\mathbf{r})\right],$$

$$\mathcal{U}(t)H_A(t){\mathcal{U}}^{†}(t)+i\hslash \frac{d\mathcal{U}(t)}{dt}{\mathcal{U}}^{†}(t)=H_E(t),$$

$$\mathcal{H}{\psi }_{\mathbf{k}s}(\mathbf{r})={ϵ}_{\mathbf{k}s}{\psi }_{\mathbf{k}s}(\mathbf{r}),$$

$$\mathcal{H}=\frac{{\hslash }^2}{2m}{\left(\frac{\nabla }{i}\right)}^2+V(\mathbf{r}),$$

$${\psi }_{\mathbf{k}s}(\mathbf{r})=e^{i\mathbf{k}\cdot \mathbf{r}}{u}_{\mathbf{k}s}(\mathbf{r}),$$

$$\mathcal{H}(\mathbf{k})u_{\mathbf{k}s}(\mathbf{r})$$

$$u_{\mathbf{k}s}(\mathbf{r})=u_{\mathbf{k}s}(\mathbf{r}+\mathbf{R}).$$

$$⟨u_{\mathbf{k}s}|u_{\mathbf{k}{s}^{\prime }}⟩=\frac{1}{v_c}{\int }_{uc}{d}^d\mathbf{r}{u}_{\mathbf{k}s}^{\ast }\left(\mathbf{r}\right)u_{\mathbf{k}{s}^{\prime }}\left(\mathbf{r}\right)={\delta }_{s{s}^{\prime }}.$$

$$H_0$$

$$c_{\mathbf{k}s}^{†}$$

$$⟨{\psi }_{\mathbf{k}s}|{\psi }_{{\mathbf{k}}^{\prime }{s}^{\prime }}⟩$$

$$\bigl{\{}c_{\mathbf{k}s},c_{\mathbf{k^{\prime}}s^{\prime}}^{\dagger}\bigr{\}}=(2\pi)^{d}\delta_{ss^{\prime}}\delta(\mathbf{k}-\mathbf{k}^{\prime}),$$

$$\sum _{\mathbf{R}}{e}^{i(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{R}}=\frac{(2\pi {)}^3}{v_C}\delta (\mathbf{k}-{\mathbf{k}}^{\prime }),$$

$$\mathbf{r}{}_{\mathbf{k}{\mathbf{k}}^{\prime },s{s}^{\prime }}$$

$${\mathbf{\xi }}_{\mathbf{k}s{s}^{\prime }}$$

$$=$$

$$\Psi (\mathbf{r})=\sum _s\int \frac{d^d\mathbf{k}}{(2\pi {)}^d}{\Phi }_s\left(\mathbf{k}\right){\psi }_{\mathbf{k}s}(\mathbf{r}),$$

$$⟨{\psi }_{\mathbf{k}s}|\mathbf{r}|\Psi ⟩$$

$$=i\sum _{s^{\prime }}$$

$${\mathbf{D}}_{\mathbf{k}s{s}^{\prime }}:={\delta }_{s{s}^{\prime }}{\nabla }_{\mathbf{k}}-i{\mathbf{\xi }}_{\mathbf{k}s{s}^{\prime }},$$

$$⟨{\psi }_{\mathbf{k}s}|\mathbf{r}|\Psi ⟩=\sum _{s^{\prime }}i{\mathbf{D}}_{\mathbf{k}s{s}^{\prime }}{\Phi }_{s^{\prime }}\left(\mathbf{k}\right),$$

$$u_{\mathbf{k}s}\to e^{i{\theta }_s(\mathbf{k})}{u}_{\mathbf{k}s},$$

$${\mathbf{D}}_{\mathbf{k}s{s}^{\prime }}$$

Full text

I Introduction

The calculation of nonlinear optical (NLO) coefficients in crystals has seen a renewed impetus, spurred by the strong nonlinear properties of layered materials like graphene Mikhailov and Ziegler (2008); Bonaccorso et al. (2010); Glazov and Ganichev (2014); Cheng et al. (2014); Peres et al. (2014); Cheng et al. (2015); Pedersen (2015); Mikhailov (2016).

Perturbative calculations of NLO coefficients in bulk semiconductors, with a full quantum treatment of matter, date back to early nineties of the previous century, and have not been entirely trouble free Moss et al. (1990); Sipe and Ghahramani (1993). In the long wavelength limit—in which the spatial dependence of the radiation electric field is neglected—, there are two representations of the radiation field: by a time dependent vector potential $\mathbf{A}(t)$, with the electric field given by $\mathbf{E}(t)=-\partial\mathbf{A}/\partial t$; by the electric dipole scalar potential, $V(\mathbf{r})=e\mathbf{E}(t)\cdot\mathbf{r}$. The advantage of the first method, known as the velocity gauge, is that the perturbation introduces no extra spatial dependence to the crystal Hamiltonian, thus preserving the crystal’s translational symmetry. This leads to a decoupling of the system’s response in momentum space: it becomes a sum of independent contributions of each $\mathbf{k}$ value in the Brillouin zone. Early attempts to calculate NLO coefficients using this approach were, however, plagued by unphysical contributions, diverging at low frequencies Moss et al. (1990). Several authors addressed this issue by separating the treatment of inter and intra band contributions, using time-dependent basis sets Genkin and Mednis (1968); Sipe and Ghahramani (1993). Later Aversa and Sipe Aversa and Sipe (1995) revisited the problem, emphasizing the gauge freedom that allows you to choose either form of the coupling to the radiation field. They recognized that the unphysical divergences mentioned above actually have coefficients that are exactly zero, expressing sum rules that they derived explicitly in first order response, and claimed to hold in all orders. Because these sums rules are easily violated in approximations, they end up advocating using the scalar potential method, also referred to as the length gauge, in actual calculations.

Similar problems were found in earlier calculations of NLO response of atoms Lamb et al. (1987). As far back as 1951, ref. Lamb (1952), W. Lamb Jr. recommended as more convenient the scalar potential gauge in perturbative calculations of the Hydrogen atom fine structure, and, for a while, the view that the two choices of gauge lead to different results was widely held Lamb et al. (1987).

The obvious advantage of the scalar potential gauge is that it is written in terms of a gauge invariant entity, the electric field, even though it expresses a specific choice of gauge for the electromagnetic field. But because the scalar potential contains the position variable, $\mathbf{r}$, the perturbation is no longer diagonal in Bloch momentum space, and couples different $\mathbf{k}$ values. Furthermore, the position operator is highly singular in momentum space, and its matrix elements can only be properly defined in the infinite crystal limit.

This choice of representation was used on the recent reduced density matrix (RDM) calculations of the nonlinear optical response of graphene Cheng et al. (2014, 2015); Mikhailov (2016). In ref. Cheng et al. (2014), the derivative term in the RDM equations of motion was removed by means of a $\mathbf{k}$-space translation thus decoupling them in crystal momentum space; this is not without cost, as the system’s response is now expressed in terms of both the $\mathbf{A}$ and the $\mathbf{E}$** **fields. In their subsequent work Cheng et al. (2015), the authors retained this derivative term as well as introduced relaxation terms to the RDM equations of motion. A different approach was proposed by Mikhailov Mikhailov (2016), who avoided the problem of the singular intra band term of the position operator by using a finite wavelength perturbation that satisfies periodic boundary conditions,

[TABLE]

and taking the limit $\mathbf{q}\to 0$ at the end of the calculation.

In addition to the freedom in expressing the external electric field, we consider the freedom of choice of the phase of the Bloch functions, in each point in the Brillouin zone,

[TABLE]

since any expectation value is necessarily independent from the choice of $\theta_{s}(\mathbf{k})$. This sets the transformation law of any observable’s matrix elements to

[TABLE]

Striving to make this property explicit leads to the concept of the covariant derivative, which will be used to derive a considerably simpler form of the RDM formalism in the length gauge.

This paper is organized as follows. In Section II, we present an overview of some concepts on gauge invariance and the key ideas regarding the crystal Hamiltonian. We then introduce the covariant derivative in momentum space, which captures the phase freedom mentioned in the previous paragraph. Section III is dedicated to the reduced density matrix. We shall derive the RDM equations of motions and, more importantly, the relation between these two objects. The rest of the section is dedicated to the equivalence between observables in the two formalisms. In Section IV, we write the solutions to the RDM equations of motion. As a proof of concept, Section V is dedicated to the study of graphene’s current response by means of the scalar potential formalism Cheng et al. (2014, 2015); Mikhailov (2016). This is followed by Section VI, where we discuss the breakdown of the scalar potential/vector potential equivalence, upon truncation of the expressions for the current for a finite set of bands Rzążewski and Boyd (2004). The last section is dedicated to a summary of our results.

II A GENERAL DESCRIPTION

The two possible representations of the uniform electric field entail two different, but equivalent, ways to write the many-body Hamiltonian. One can either add the dipole interaction to the single particle Hamiltonian of the unperturbed system, $\mathcal{H}_{0}$, 111The charge $q$ has already been replaced by the charge of an electron, $q=-e$.

[TABLE]

or use the minimum coupling procedure in $\mathcal{H}_{0}$,

[TABLE]

The electron field, $\Psi(\mathbf{r})$, and its Hermitian conjugate, $\Psi^{\dagger}(\mathbf{r})$, satisfy the usual anti-commutation relations.

The many-body state vector in the vector potential approach, $\left|\psi(t)\right\rangle$, evolves in time according to the Hamiltonian $H_{A}(t)$

[TABLE]

A second state vector $\left|\bar{\psi}(t)\right\rangle$, obtained via a time-dependent unitary transformation of $\left|\psi(t)\right\rangle$,

[TABLE]

has an equation of motion

[TABLE]

If the unitary transformation is chosen as

[TABLE]

where $\rho(\mathbf{r}):=\Psi^{\dagger}(\mathbf{r})\Psi(\mathbf{r})$ is the density operator, it is straightforward to show that

[TABLE]

implying that $\left|\bar{\psi}(t)\right\rangle$ is the state vector in scalar potential gauge.

Observables in the two gauges are also related by a unitary transformation, $O_{E}:=\mathcal{U}(t)\,O_{A}(t)\,\mathcal{U}^{\dagger}(t)$. The exception is the Hamiltonian; because the unitary transformation, $\mathcal{U}(t)$, is time dependent, the time evolution operator in the length gauge, $H_{E}(t)$, is not simply $\mathcal{U}(t)\,H_{A}(t)\,\mathcal{U}^{\dagger}(t)$, but has an additional term involving the time derivative of $\mathcal{U}(t)$, Eq. (7). The existence of this transformation between the two descriptions of the radiation field establishes their complete equivalence Lamb et al. (1987); Aversa and Sipe (1995).

Next we recall some important results of electron eigenstates in an unperturbed crystal. The single particle Schrödinger equation is Ashcroft and Mermin (1976)

[TABLE]

with

[TABLE]

and $V(\mathbf{r})=V(\mathbf{r}+\mathbf{R})$, for $\mathbf{R}$ any Bravais lattice vector. According to Bloch’s theorem, the eigenfunctions have the form of a plane wave times a periodic function,

[TABLE]

allowing the eigenvalue problem to be expressed in terms of the $\mathbf{k}$-dependent Hamiltonian, $\mathcal{H}(\mathbf{k}):=e^{-i\mathbf{k}\cdot\mathbf{r}}\,\mathcal{H}\,e^{i\mathbf{k}\cdot\mathbf{r}}$,

[TABLE]

The function $u_{\mathbf{k}s}(\mathbf{r})$ is a periodic function in the real space unit cell,

[TABLE]

Each $\mathbf{k}$-point in the First Brillouin Zone (FBZ) defines an Hamiltonian operator, $\mathcal{H}(\mathbf{k})$, that acts on functions whose domain is the real space unit cell, and that satisfy the boundary condition of Eq. (12). The eigenfunctions of $\mathcal{H}_{\mathbf{k}}$, $\left\{\left|u_{\mathbf{k}s}\right\rangle,\,s=0,1,\dots\right\}$ are a basis of such functions. We can assume this basis to be orthonormal with an inner product defined as the integral over the real space unit cell (volume $v_{c}$),

[TABLE]

Different values of $\mathbf{k}$ have different basis, for they are eigenfunctions of different Hamiltonians. Here, the Bloch wave vector is a continuous parameter, even in the finite volume crystal, as the eigenvalues in Eq. (11) are well defined for every $\mathbf{k}$ in the FBZ. The $\mathbf{k}$ value selection by periodic boundary conditions only involves the plane wave factor of the Bloch function, and has no bearing on the periodic part, $u_{\mathbf{k}s}\left(\mathbf{r}\right)$ Ashcroft and Mermin (1976). As such, derivatives with respect to $\mathbf{k}$ of $u_{\mathbf{k}s}\left(\mathbf{r}\right)$ are always well defined whereas derivatives of plane wave factors require the infinite volume limit. We shall work in this limit from the start, due to the difficulties of properly defining the position operator, $\mathbf{r}$, in a finite system with periodic boundary conditions.

In the $\Omega\to\infty$ limit, ($\Omega$, the volume of the crystal) the momentum sums are replaced by $d$-dimensional integrals over the FBZ. The many-body crystalline Hamiltonian then reads as

[TABLE]

for $c_{\mathbf{k}s}^{\dagger}$, $c_{\mathbf{k}s}$ the creation and destruction operators of Bloch states,

[TABLE]

The Bloch state orthogonality relation, the anti-commutation relations, and the lattice sum rule are suitably modified,

[TABLE]

[TABLE]

[TABLE]

so that the operator $c_{\mathbf{k}s}^{\dagger}c_{\mathbf{k}s}$ is a density in momentum space and not a dimensionless number operator as in the finite volume case.

In the scalar potential approach, the perturbation is written in terms of the position operator, $\mathbf{r}$. Its matrix elements are ill-defined in the finite volume system, but can be computed for $\Omega\rightarrow\infty$. In that limit they read Blount (1962),

[TABLE]

where the Berry connection, $\boldsymbol{\xi}_{\mathbf{k}ss^{\prime}}$, is defined as a scalar product in the real space unit cell, independent of the crystal’s volume,

[TABLE]

The somewhat awkward looking expression of Eq. (19), can be cast in a more transparent form if we bear in mind that, for a continuous non-normalizable basis, the matrix elements of an operator are a kernel of an integral transform. In the Bloch representation, a general single particle state is represented as

[TABLE]

for $\Phi_{s}\left(\mathbf{k}\right)=\langle\psi_{\mathbf{k}s}|\Psi\rangle$. The wave function for the state $\mathbf{r}\left|\Psi\right\rangle$ is

[TABLE]

where, to reach the final expression, we have integrated by parts. The integration is over the FBZ, and, given the periodicity of any function of $\mathbf{k}$, $\phi(\mathbf{k})=\phi(\mathbf{k}+\mathbf{G})$, there are no surface terms. This prompts us do define the covariant derivative operator Cheng et al. (2015),

[TABLE]

such that

[TABLE]

*i.e., *the position operator is $i\mathbf{D}_{\mathbf{k}ss^{\prime}}$ in the Bloch representation Pitaevskii and Lifshitz (1980). The designation of *covariant *refers to its behavior under a local gauge transformation in momentum space,

[TABLE]

for which

[TABLE]

The gradient term of the phase $\theta_{s^{\prime}}(\mathbf{k})$ is canceled by the transformation of the Berry connection; this is in complete parallel to the definition of covariant derivative in gauge theories in real space.

In the vector potential approach, the perturbation is written in terms of the velocity matrix elements, $\mathbf{v}{}_{\mathbf{k}\mathbf{k}^{\prime}ss^{\prime}}:=\mathbf{p}{}_{\mathbf{k}\mathbf{k}^{\prime}ss^{\prime}}/m_{e}$, which are diagonal in** $\mathbf{k}$**-space, and expressible as matrix elements in the basis of periodic functions,

[TABLE]

Since the velocity operator is quite generally $(i\hbar)^{-1}\left[\mathbf{r},H\right]$, it is not surprising to find that

[TABLE]

This turns out to be the expression of a more general result, which will prove useful later and which we now discuss.

Any matrix in the space generated by the basis $\left\{\left|u_{\mathbf{k}s}\right\rangle,\,s=1,2,\dots\right\}$ , parametrized by $\mathbf{k}$, defines an operator in band space by

[TABLE]

Two examples of this are the Hamiltonian, $\mathcal{H}(\mathbf{k}):=\epsilon_{\mathbf{k}s}\delta_{ss^{\prime}}$ and the velocity operator. Consider the matrix elements of the $\mathbf{k}$-derivative of $\mathcal{\mathcal{O}}(\mathbf{k})$,

[TABLE]

The first term on the right hand side is simply the $\mathbf{k}$-gradient of the matrix element, $\nabla_{\mathbf{k}}\mathcal{\mathcal{O}}_{\mathbf{k}ss^{\prime}}$. The other two can be expressed in terms of the Berry connection, by application of the completeness relations \sum_{r}$$\left|u_{\mathbf{k}r}\right\rangle\left\langle u_{\mathbf{k}r}\right|=\hat{1},

[TABLE]

The matrix element, Eq. (31), then reads as a commutator with the covariant derivative,

[TABLE]

which can be alternatively represented in operator form,

[TABLE]

The relation between observables in the velocity and length gauges can be cast in this language. First, we note that

[TABLE]

Using the form of the Bloch functions, Eq. (10), this integral over all space can reduced to one over a unit cell, summed over all of them, giving,

[TABLE]

or,

[TABLE]

The relationship between operators in the two descriptions can also be expressed in terms of these objects defined in the same k-point. To do so, we expand the RHS of that equation in powers of $\mathbf{A}(t)$. It follows from Eq. (36), that

[TABLE]

where a sum over the repeated cartesian indexes $\alpha_{j}$ is left implied. To conclude this brief account of the use of the covariant derivative, we point out that the canonical commutation relation,

[TABLE]

is expressed in the Bloch basis as

[TABLE]

since $\hat{r}^{\alpha}=iD^{\alpha}$ and $\hat{p}^{\alpha}=m_{e}\mathcal{V}^{\alpha}$. Eq. (43) can also be explicitly derived from the form of the $\mathbf{D}_{\mathbf{k}}$ and $\mathbf{v}_{\mathbf{k}}$ matrices.

We can now write the general many-body Hamiltonians, (2) and (3), in the Bloch description,

[TABLE]

as well as their respective current operators,

[TABLE]

The expression of the current in the velocity gauge is a consequence of Eqs. (40) and (43). Any component of velocity in operator form satisfies the general relation between $\mathbf{k}$ diagonal observables in both gauges,

[TABLE]

An expansion of the right hand side in powers of A produces only two terms of orders $A^{(0)}$ and $A^{(1)}$. The remaining terms involve two or more derivatives with respect to k, which, following Eq. (36), can be expressed in terms of commutators of the covariant derivative and \bigl{[}D_{\mathbf{k}}^{\beta},\mathcal{V}^{\alpha,E}(\mathbf{k})\bigr{]}. As the latter is proportional to the identity operator, Eq. (43), the commutators are exactly zero and the higher order terms vanish.

[TABLE]

In matrix form this is

[TABLE]

where $\mathbf{v}{}_{\mathbf{k}ss^{\prime}}$ is the velocity matrix in the length gauge.

III the reduced density matrix

The reduced density matrix (RDM), is a matrix in band space, defined by the average of momentum conserving inter band transitions ($s\neq s^{\prime}$) and the intra band ($s=s^{\prime}$) transitions 222Note the switch in the band indexes of the RDM..

[TABLE]

The superscript $\alpha$ now denotes the gauge in which this object is computed, $\alpha=A,E$. The average of any operator $\hat{A}$ is the trace

[TABLE]

where, $\rho^{\alpha}$, the full many-body density matrix is

[TABLE]

and \bigl{\{}\left|\psi_{n}^{\alpha}\right\rangle\bigr{\}} is a complete set of state vectors. In the Schrödinger picture, the time evolution of $\rho^{\alpha}(t)$ is governed by the time-evolution of the state vectors, $\left|\psi_{n}(t)\right\rangle$,

[TABLE]

The equation of motion of the RDM takes the form

[TABLE]

The gauge freedom to express the uniform electric field is thereby carried into the calculation of the system’s dynamics, which can be described either in terms of $\rho_{\mathbf{k}ss^{\prime}}^{E}(t)$ or $\rho_{\mathbf{k}ss^{\prime}}^{A}(t)$.

Computing the commutators on the right hand side of Eq. (54), using the Hamiltonians, Eqs. (LABEL:eq:HAMB1) and (45), we obtain closed equations of motion for the RDM 333In order to simplify notation we will drop the crystal momentum label from the objects inside the commutators and write it alongisde the band indexes. ,

[TABLE]

Here we have defined $\epsilon_{\mathbf{k}ss^{\prime}}:=\epsilon_{\mathbf{k}s}-\epsilon_{\mathbf{k}s^{\prime}}$. The equation for the scalar potential RDM is found in references Aversa and Sipe (1995); Cheng et al. (2014, 2015), but not, as here, cast in terms of the covariant derivative. The presence of the derivative with respect to $\mathbf{k}$ in Eq. (55) couples the response at different values of $\mathbf{k}$, whereas, its counterpart for the vector potential gauge, Eq. (56), is completely decoupled in crystal momentum, $\mathbf{k}$, and can thus be solved independently for each point of the FBZ. Averages of single particle observables, diagonal in momentum space, such as the currents, (46) and (47), can be obtained from the RDM’s as traces over band space

[TABLE]

for $\hat{1}$ the identity in band space. Given that these two alternative formulations are related by a unitary transformation, Eq. (57) and Eq. (LABEL:eq:CUR2) have to yield the same results, although this is far from obvious at this point. To relate these two objects, we must first establish the relation between RDMs.

The full many-body density matrix $\rho^{A}$,

[TABLE]

can be expressed in the state vectors of its counterpart description by means of a unitary transformation, Eq. (5),

[TABLE]

The vector potential RDM is therefore expressible as averages with $\rho^{E}$, of suitably modified operators

[TABLE]

The new creation operator $\tilde{c}_{\mathbf{k}s^{\prime}}^{\dagger}$ is obtained from $c_{\mathbf{k}s^{\prime}}^{\dagger}$ by the same unitary transformation

[TABLE]

and creates an electron with wave function

[TABLE]

While this is clearly a Bloch state with wave vector $\mathbf{q}=\mathbf{k}+e\mathbf{A}(t)/\hbar$, it is not $\psi_{\mathbf{k}+e\mathbf{A}(t)/\hbar\,s^{\prime}}$, but can be expanded as a linear combination of the shifted Bloch states,

[TABLE]

This allows us to relate $\tilde{c}_{\mathbf{k}s^{\prime}}^{\dagger}$ and $\tilde{c}_{\mathbf{k}s^{\prime}}$ to the original creation and destruction operators

[TABLE]

and express the RHS of Eq. (61) in terms of the shifted Bloch operators. Since

[TABLE]

using this and Eq. (65) in Eq. (61) we obtain the relation between $\rho_{ss^{\prime}}^{A}(\mathbf{k},t)$ and $\rho_{ss^{\prime}}^{E}(\mathbf{k},t)$ as

[TABLE]

Following our definition of operators in band space, Eq. (29), this equality can be cast in operator form444The hats will be henceforth dropped.,

[TABLE]

The simplicity of this relation between the density matrices in the velocity and length gauges, which amounts to a simple shift in the value of $\mathbf{k}$, is only true in the operator representation; as can be seen in Eq. (67), the matrix elements of these two operators between states of a single basis $\left\{\left|u_{\mathbf{k}s}\right\rangle,\,s=0,1,\dots\right\}$ with the same** **$\mathbf{k}$, do not satisfy this simple relation.

We can now show that the expectation values of observables in the two descriptions are exactly the same. Consider,

[TABLE]

Both the RDM and the observable operators are equal to their scalar potential counterparts when their argument is adequately translated, Eqs. (40) and (68),

[TABLE]

This shifts the FBZ by a constant, which is irrelevant as the integrand is periodic. This means that,

[TABLE]

The sum rules in reference Aversa and Sipe (1995) can be traced back to this equivalence (see Appendix A).

IV Solutions to the rdm equations of motion

A system’s current response (in either scalar or vector potential formalism) is obtained from the solution to its respective RDM equation of motion. To write these solutions explicitly, we must employ the general procedure from nonlinear physics: break the RDM into contributions of different powers on the external field, $\rho=\rho^{(0)}+\rho^{(1)}+(...)$ and then proceed to iteratively solve the equations of motion for each order. For $\rho^{(n)}$, these read as

[TABLE]

[TABLE]

In expressing the time-dependent objects in terms of their Fourier decomposition, we assume adiabatic switching of the perturbation ($\eta\rightarrow 0^{+}$),

[TABLE]

The time derivative in the equations of motion is replaced by a frequency factor that is collected into an energy denominator,

[TABLE]

From this point on, the frequency argument of these energy denominators is understood to have an infinitesimal imaginary part. Using the Hadamard product of two matrices,

[TABLE]

we write recursion relations for the RDM solutions,

[TABLE]

[TABLE]

The successive application of these expressions, brings the $n$-th order solution, $\rho_{\mathbf{k}ss^{\prime}}^{(n)}(\omega)$, to the form of nested commutators of the zeroth order one, $\rho_{\mathbf{k}ss^{\prime}}^{(0)}$, which is the Fermi-Dirac distribution function times the unit matrix in band space, $\rho_{\mathbf{k}ss^{\prime}}^{(0)}=f_{\mathbf{k}s}\delta_{ss^{\prime}}$. With one last bit of notation,

[TABLE]

we write $\rho_{\mathbf{k}ss^{\prime}}^{(n)}$ in each formalism as

[TABLE]

[TABLE]

The use of the nested commutator and the factor $d(\omega)$ allow for a compact form of these solutions (albeit hiding their considerable complexity). These solutions are entirely written in terms of three objects (and their derivatives, in the case of $\rho^{E}$): the band energies, $\epsilon_{\mathbf{k}s}$; the Berry connection, $\xi_{\mathbf{k}ss^{\prime}}$; and the structure of filled/empty bands, $f_{\mathbf{k}s}$. In principle, by determining these, one determines a system’s nonlinear response, provided that one can compute the FBZ integrals in Eqs.(57) and (LABEL:eq:CUR2). In section VI, we use these expressions to discuss whether it is possible to truncate the sums over bands, implicit in these matrix products, to a reduced set. But before, we study the first two non-trivial orders of the monolayer graphene response to a uniform electric field, to illustrate that these reproduce the results of references Cheng et al. (2014, 2015); Mikhailov (2016).

V Linear and third order response in the monolayer graphene

Monolayer graphene is (usually) described by a tight-binding model that considers only nearest neighbour hoppingCastro Neto et al. (2009) and for which the expressions for the two bands, $\epsilon_{\mathbf{k}s}$, and the periodic functions, $u_{\mathbf{k}s}$, can be explicitly computed. This allows one to derive some useful properties.

First, a double band-index sum, such as the ones in $\langle\mathbf{J}_{E}\rangle$ and $\langle\mathbf{J}_{A}\rangle$, over an antisymmetric object $\theta_{ss^{\prime}}=-\theta_{s^{\prime}s}$, can be reduced to a single sum, where $\bar{s}$ reads as the band opposite to $s$,

[TABLE]

Second, the Berry connection for the monolayer is described by an intra band and an inter band term. With adequate choice of gauge (Eq. 24) the following properties can be obtained Cheng et al. (2014, 2015); Mikhailov (2016),

[TABLE]

In addition, we can choose these to be even under $\mathbf{k}\rightarrow-\mathbf{k}$.

Linear order response

The first order term of the RDM is

[TABLE]

Writing the current, Eq. (57), in terms of its Fourier components, we obtain two contributions Cheng et al. (2014, 2015); Mikhailov (2016),

[TABLE]

The intra band contribution is a generalized Drude term,

[TABLE]

and the inter band contribution, $\Pi_{\textrm{inter}}^{\alpha\beta,(1)}(\mathbf{k},\omega)$, which involves the Berry connection,

[TABLE]

Third order response

Because graphene has inversion symmetry, its second order response is zero (see Appendix B). The first nonlinear contribution to the current is thus the third order one, which is to be computed here. From the general RDM solution, Eq. (LABEL:eq:RECURSIVE1), we obtain for $n=3$,

[TABLE]

Expanding this by means of Eq. (34), produces a plethora of terms which we organize, following Mikhailov, by the number of intra band (derivatives) and inter band (Berry connections) factors. These can be manipulated independently and give contributions to the current, that can be labeled by $i=1,2,3$, the number of intra band factors.

The field factors and the FBZ integrations are the same in all these contributions, and so

[TABLE]

The $\Pi_{3}$ contribution is the single term in Eq. (LABEL:eq:RHO3) with three $\mathbf{k}$-derivatives, while $\Pi_{2}$ gathers three terms with two $\mathbf{k}$-derivatives, to which we can apply the property Eq. (84). Contributions $\Pi_{1}$ and $\Pi_{0}$ require additional manipulations, presented in Appendix C. Consistently with our notation $\epsilon_{\mathbf{k}s\bar{s}}:=\epsilon_{\mathbf{k}s}-\epsilon_{\mathbf{k}\bar{s}},$ we abbreviate the factor $f_{\mathbf{k}s}-f_{\mathbf{k}\bar{s}}$ as $f_{\mathbf{k}s\bar{s}}$. The expressions for the $\Pi_{i}$ are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Apart from differences in the Cartesian and frequency indexes, which can always be relabeled, these are the expressions found in Cheng et al. (2014, 2015); Mikhailov (2016), in the limit where $\Gamma_{(i)}$, $\Gamma_{(e)}$, the phenomenological scattering rates, are set to zero.

VI Effective hamiltonians

As we have seen in the previous sections, the current response in the Schrödinger problem can be computed by two different but equivalent procedures. Although this is conceptually important, one is ultimately interested in computing the current in materials described by effective Hamiltonians, that account only for a finite number of bands.

We will now show that, by truncating the band space, the scalar and vector potential currents are no longer the same; the latter contains relevant contributions from bands that are left out of the effective Hamiltonian.

In Fig. (1), we present a conceptual picture of a spectrum which has a cluster of bands close to the Fermi level, which we deem relevant, well separated in energy by the bands below (filled) and above (empty). We denote the energy scale in the subspace $\mathcal{E}_{0}$ of relevant bands by $\delta$ and the energy separation to other bands by $\Delta\gg\delta$; we assume the frequency of the external field, $\omega$, to be of the order, $\omega\sim\delta/\hbar$. The question we wish to answer is whether the bands outside $\mathcal{E}_{0}$ can be ignored in the calculation of the current.

Energy denominators, $d(\omega)$, involving bands inside $\mathcal{E}_{0}$,

[TABLE]

will be larger then those involving transition to and from bands in $\mathcal{E}_{0}$ and those outside

[TABLE]

In the $n$-th order contribution to the current in the scalar potential gauge, each term is a trace of a product $n+1$ matrices in band space.

[TABLE]

One such term is the fully intra band one, in which we pick only the diagonal part of each covariant derivative operator (see an example in Eq. (94)). This term is non zero only for the band that contains the Fermi level; it has no contributions from bands outside $\mathcal{E}_{0}.$

All inter band contributions contain at least one difference of occupation factors, thus allowing us to discard any terms that involve only filled, or only empty bands. What remains are contributions of three types: (a) terms that involve transition between filled and empty bands, both outside $\mathcal{E}_{0}$; (b) terms that involve transition between bands in $\mathcal{E}_{0}$ and bands outside $\mathcal{E}_{0}$; (c) inter band terms among the bands in $\mathcal{E}_{0}$.

In Eq. (99), if $s^{\prime}$ belongs to a filled band outside $\mathcal{E}_{0}$, and $s$ to a empty band, also outside $\mathcal{E}_{0}$ (or vice-versa), $\epsilon_{\mathbf{k}ss^{\prime}}\sim\Delta\gg\hbar\omega$ and

[TABLE]

the matrix that follows this term in Eq. (99), must include at least an energy denominator $\sim\mathcal{O}(1/\Delta)$ because the last band index is the same as the first, $s^{\prime}$. If $s^{\prime}$ and $s$ both refer to a filled band (or an empty one) outside $\mathcal{E}_{0}$ the inter band terms between filled and empty bands must have at least two such energy denominators. In other words, terms of type (a) have a least and energy denominator of order $\sim\mathcal{O}(1/\Delta).$ An identical argument can be made for transitions between bands inside and and outside $\mathcal{E}_{0}$, i.e. for transitions of type (b). We conclude that, in the limit of $\Delta\gg\delta$, the dominant terms come from bands in $\mathcal{E}_{0}$ (terms of type (c)) and truncation to this subspace is a valid approximation.

The same argument does not carry to the corresponding contribution to the current in the velocity gauge:

[TABLE]

In this case, every energy denominator is associated with a velocity matrix element. In any transition involving energies $\epsilon_{rr^{\prime}}\sim\mathcal{O}(\Delta)$, the $d_{rr^{\prime}}(\omega)\sim\mathcal{O}(1/\Delta)$, as before, but the corresponding velocity matrix element has an off-diagonal contribution $v_{\mathbf{k}r^{\prime}r}^{\beta}\sim\xi_{\mathbf{k}rr^{\prime}}^{\beta}\mathcal{O}(\Delta)$, and such terms give relevant contributions no matter how large $\Delta$ is. In other words, bands away from the Fermi surface contribute just as much as those in $\mathcal{E}_{0}$ for the expression of the current in the velocity gauge.

VII Summary and conclusions

The concept of covariant derivative in $\mathbf{k}$-space has been shown here to be of considerable value in the calculation of the nonlinear current response. It is a very convenient representation of the position operator [Eq. (23)]; it clarifies the structure of the velocity matrix [Eq. (27)]; it allows a complete parallel development of the structure of the reduced density matrices (RDM) in the length (scalar potential) and velocity (vector potential) gauges [Eqs. (55) and (56)] and it provides compact expressions for the perturbative solutions of the equations of motion of the RDM [Eqs. (LABEL:eq:RECURSIVE1) and (LABEL:eq:RECURSIVE2)]. It also allowed us to see how the equivalence between objects in the two gauges breaks down when the band space is truncated.

Furthermore that the approximation is legitimate in the length gauge, but fails in the velocity one is made clear by considering the truncation of the covariant derivative and the velocity operators to a restricted set of bands. The commutator of these two quantities, Eq. (43), which is constant for the case of the infinite bands of the Schrödinger Hamiltonian, no longer holds for a truncated subset of these. In fact, the statement that the commutator is constant is equivalent to the Hamiltonian being linear or quadratic in k (the former is the case of the Dirac Hamiltonian, where the commutator gives zero). As a result, Eq. (48) is no longer valid for a general effective Hamiltonian and there is no equality between currents in the two gauges, at least as they are written in Eqs. (46) and (47). In order to use the velocity gauge in actual calculations, one must start from the beginning with the effective Hamiltonian and perform the minimum coupling then. Naturally, this means modifying both the equation of motion of the RDM, Eq. (56), and the current operator in the velocity gauge. This will be the subject of a future paper.

VIII Acknowledgments

The work of G.B.V and D.J.P is supported by Fundação para a Ciência e Tecnologia (FCT) under the grants PD/BI/129220/2017 and PD/BD/135019/2017 respectively. N.M.R.P. acknowledges funding from the European Commission within the project “Graphene-Driven Revolutions in ICT and Beyond” (ref. no. 696656) and the Portuguese Foundation for Science and Technology (FCT) in the framework of the Strategic Financing UID/FIS/04650/2013.

Appendix A The Aversa and Sipe Sum rules

In reference Aversa and Sipe (1995), Aversa and Sipe sketch how problems might arise in the velocity gauge as it involves additional terms (with respect to the length gauge) that are zero only when they are treated exactly. They show this explicitly in the linear response.

In our formalism this can be seen in full generality, starting from Eqs. (41) and (70). Expanding in powers of $\mathbf{A}(t)$, we can write $\mathcal{\langle O}^{A}\left(\mathbf{k},t\right)\rangle$ in terms of $\langle\mathcal{O}^{E}\left(\mathbf{k},t\right)\rangle$ and contributions that depend explicitly on $\mathbf{A}(t)$, that is, the $n\geq 1$ terms of the sum,

[TABLE]

The equivalence of the two gauges requires that everything other than $\langle\mathcal{O}^{E}\left(\mathbf{k},t\right)\rangle$ in the right hand side to be zero. Each contribution contains a factor which, for arbitrary order $n$, reads as

[TABLE]

for $\mathrm{G}(\mathbf{k},t)$ is some matrix in band space,

[TABLE]

The two terms in the expression (102) amount to

[TABLE]

The first term integrates to zero due to the periodicity of $\mathrm{G}(\mathbf{k},t)$ in reciprocal space. As for the second one, the trace of any commutator of two matrices is zero, which follows from the cyclic invariance of the trace.

The cancellation of all terms in the sum on the RHS of Eq. (101) constitute the sum rules referred to by Aversa and Sipe. They are an order by order formulation of the result that in an integral over the FBZ of a function with the period of the reciprocal lattice, the argument of the integrand can be shifted without changing the integral.

Nevertheless, it is clear from this formulation that these extra terms exist only if, starting from the velocity gauge, we try to reduce our expressions to the ones in the length gauge.

Appendix B Current response in a centrosymmetric material

For a centrosymmetric crystal, the spatial inversion operator $\mathcal{P}$ commutes with the crystal Hamiltonian, $\mathcal{H}$,

[TABLE]

and the solutions to the $\mathbf{k}$-dependent Hamiltonian, $u_{\mathbf{k}s}(-\mathbf{r})$ and $u_{-\mathbf{k}s}(\mathbf{r})$, are related by a phase phase factor,

[TABLE]

If we take $\mathbf{r}\rightarrow-\mathbf{r}$ in the integral that defines the Berry connection, Eq. (21), we can determine how it behaves when we exchange the sign of the crystal momentum,

[TABLE]

and,

[TABLE]

This determines the transformation law for the covariant derivative:

[TABLE]

We can now determine how the $\mathbf{k}$-dependent factor of $\rho^{(1)}$,

[TABLE]

transforms upon $\mathbf{k}\rightarrow-\mathbf{k}$, by recalling that $f_{\mathbf{k}s}$ and $d_{\mathbf{k}ss^{\prime}}(\omega)$ are even, since they depend only on the band energies,

[TABLE]

This can be extended to the higher order contributions of the RDM, in particular to the second order one, easily extracted from Eq. (LABEL:eq:RECURSIVE1),

[TABLE]

It follows from Eqs. (107) and (109) that

[TABLE]

This object picks up the same $\mathbf{k}$-space phase factors, but unlike its first order counterpart, the sign does not change. Combining this with the transformation law for the velocity matrix element,

[TABLE]

we see that the integrand in the FBZ integral, $v_{\mathbf{k}s^{\prime}s}^{\beta}\Pi_{\mathbf{k}ss^{\prime}}^{(2),\alpha_{1}\alpha_{2}}(\omega,\omega_{1})$ is an odd function of $\mathbf{k}$, the second order current vanishes.

This argument also carries for an arbitrary order $n$. The $\mathbf{k}$-parity of $v_{\mathbf{k}s^{\prime}s}^{\beta}\Pi_{\mathbf{k}ss^{\prime}}^{(n),\alpha_{1}(...)\alpha_{n}}$ is determined by its number of covariant derivatives in $\Pi_{\mathbf{k}ss^{\prime}}^{(n),\alpha_{1}(...)\alpha_{n}}$. For $n$ even, the integrand is odd under $\mathbf{k}\rightarrow-\mathbf{k}$, so even order contributions to the current vanish in a centrosymmetric material.

Appendix C Deriving the expressions: (97)-(98)

Consider $\Pi_{1}$, the collection of terms with only one intra band factor, where we have used the two-band character of the monolayer graphene, Eq. (84),

[TABLE]

For $s^{\prime}=\bar{s}$, the two terms in each line cancel out by application of the Berry connection properties, Eq. (85) and (86). This fixes $s^{\prime}=s$, and the $\Pi_{1}$ contribution reads as

[TABLE]

The $\Pi_{0}$ portion of the current for the two band material reduces to,

[TABLE]

Since the Berry connection is even under $\mathbf{k}\rightarrow-\mathbf{k}$, the velocity matrix element in band space, Eq. (28), is written as the sum of two contributions of opposite parity: an odd intra band term and an even inter band term. The integration over the FBZ carries cancels all odd terms, and so the intra band part of $\mathbf{v}_{\mathbf{k}ss^{\prime}}$ can be ignored: this fixes $s^{\prime}=\bar{s}$, and

[TABLE]

Setting $r=s$ and using (85) and (86), the first two terms of this expression cancel out; the last two cancel when $r=\bar{s}$.

Finally, $\Pi_{0}$ reduces to,

[TABLE]