# Regularized lattice theory for spatially dispersive nonlinear optical
conductivities

Steven Gassner

Department of Physics and Astronomy, University of Pennsylvania, Philadelphia,
PA 19104  E. J. Mele

Department of Physics and Astronomy, University of Pennsylvania, Philadelphia,
PA 19104

###### Abstract

Nonlinear optical responses are becoming increasingly relevant for
characterizing the symmetries and quantum geometry of electronic phases in
materials. Here, we develop an expanded diagrammatic scheme for calculating
spatially dispersive corrections to nonlinear optical conductivities, which we
expect to enhance or even dominate even-order responses in materials of recent
interest. Building upon previous work that enforces gauge invariance of
spatially uniform nonlinear optical responses, we review the cancellation of
diagrams required to ensure the equivalence between velocity gauge and length
gauge formulations, and provide a simple vertex rule for extending optical
responses to first order in the light wave vector $\mathbf{q}$. We then
demonstrate the method with calculations on a prototypical centrosymmetric
model where spatial dispersion admits anomalous second-harmonic generation, a
response that is symmetry-forbidden under the dipole approximation.

##  I Introduction

Optical responses are an extremely valuable tool for understanding the
electronic structure of quantum matter, from materials to molecules Boyd
(2003); Sekino and Bartlett (1986). Not only are they vital for characterizing
a host of different devices, from photovoltaics to lasers, but they are also
immensely useful as probes for symmetry breaking and quantum geometric
properties in materials Young  _et al._ (2012); Morimoto and Nagaosa (2016);
Wu  _et al._ (2017); de Juan  _et al._ (2017); Ma  _et al._ (2017); Chan  _et
al._ (2017); Holder  _et al._ (2020); Ni  _et al._ (2022); Ahn  _et al._
(2022a); Tai and Claassen (2023). Over the past several decades, the
development of a straightforward and general means of systematically
calculating these responses has been a surprisingly difficult task, revealing
a number of subtleties that are still being worked out to this day. A main
reason is that there are different choices of gauge for describing the light-
matter coupling with various advantages and disadvantages. One option is to
couple the electric field $\textbf{E}(t)$ with the position operator
$\hat{\textbf{r}}$,

$\hat{\mathcal{H}}\to\hat{\mathcal{H}}-{\rm
q}\,\hat{\textbf{r}}\cdot\textbf{E}(t)\text{,}$ (1)

which goes by the name “length gauge” (also known as “position gauge” or
”dipole gauge.”) Here, ${\rm q}=-|e|$ denotes the electron charge. Another
option, called the “velocity gauge,” employs a minimal substitution in the
$\mathbf{k}$-dependent Bloch Hamiltonian,

$\hat{\mathcal{H}}(\mathbf{k})\to\hat{\mathcal{H}}(\mathbf{k}-{\rm
q}\,\mathbf{A}(t))\text{,}$ (2)

where $\mathbf{A}(t)$ is the magnetic vector potential satisfying
$\mathbf{E}(t)=-\partial_{t}\mathbf{A}(t)$. (At first order in
$\mathbf{A}(t)$, this takes the form of a coupling $-{\rm
q}\,\hat{\mathbf{v}}\cdot\mathbf{A}(t)$ where
$\hat{\mathbf{v}}=\nabla_{\mathbf{k}}\hat{\mathcal{H}}$ is the velocity
operator; hence the name “velocity gauge.”) These two gauges are related by a
time-dependent unitary transformation, and are hence equivalent Kobe and Smirl
(1978) in a fully microscopic theory. However maintaining this equivalence for
effective Hamiltonians that are projected into a partial sector of the full
Hilbert space requires additional consideration.

Velocity gauge, while having the advantage of being diagonal in k-space and
easier for numerical calculations, has the serious disadvantage of introducing
artificial low-frequency divergences that must be regularized systematically.
Historically, this has been attempted using sum rules, first developed by
Aversa and Sipe Sipe and Ghahramani (1993); Aversa and Sipe (1995) and later
generalized by Ventura et al Ventura  _et al._ (2017); Passos  _et al._
(2018), which formally show that the weight of these divergent terms vanishes
only after a full k-space integration. In more recent work with Wannierized
tight-binding models, Schüler et al identify a sum rule encoding the
cancellation between the paramagnetic and diamagnetic currents calculated at
linear order in $\mathbf{A}(t)$, and show that enforcing this sum rule
improves the accuracy of velocity gauge calculations Schüler  _et al._ (2021).

In this work, we relate a convenient diagrammatic method Parker  _et al._
(2019) and the traditional reduced density matrix perturbation theory approach
Ventura  _et al._ (2017); Passos  _et al._ (2018). Representing terms in the
perturbation theory by their associated diagrams, we can understand in a
physically transparent manner the general cancellation between the
paramagnetic and diamagnetic contributions in the linear response, and
investigate the analogs of that cancellation that occurs for nonlinear
responses. We also introduce a scheme for extending these results to first
order in the wavevector $\mathbf{q}$ of the incident radiation. We call these
responses spatially dispersive, because they take into account at lowest order
the effects of the spatial gradient of the electric field on the optical
response. As a prototypical example, we investigate the anomalous second-
harmonic generation (SHG) response that appears in a centrosymmetric 1D model.
While SHG is typically understood to vanish in centrosymmetric systems in the
dipole approximation (and hence is often used experimentally as a probe for
inversion breaking in crystals) this selection rule is violated when one takes
into account the spatial variation of the electromagnetic field on the length
scale of the electronic states being coupled in optical transitions. This
becomes especially relevant for artificial lattices where the lattice
constants are inflated and can greatly exceed the microscopic atomic scale or
more generally for band structures in which the Wannier representation
introduces coherence over large length scales.

##  II Calculations of nonlinear optical responses

###  II.1 Generalities

To calculate nonlinear optical conductivities, we employ the density matrix
equation of motion approach detailed in Refs. Ventura  _et al._ (2017); Passos
_et al._ (2018) and connect it to the diagrammatic approach developed in Ref.
Parker  _et al._ (2019). Recall the density matrix equation of motion,

$i\hslash\,\partial_{t}\hat{\rho}(t)=[\hat{\mathcal{H}}(t),\hat{\rho}(t)]\text{,}$
(3)

Assume that $\hat{\mathcal{H}}(t)=\hat{\mathcal{H}}^{(0)}+\hat{V}(t)$, with
$\hat{V}(t)$ a weak perturbation. Transforming to Fourier representation,

$\hat{\rho}(t)=\int d\omega\,e^{-i\omega t}\,\hat{\rho}(\omega)\text{,}$ (4)
$\hat{V}(t)=\int d\omega\,e^{-i\omega t}\,\hat{V}(\omega)\text{,}$ (5)

we can re-express the equation of motion as an iterative equation in the
eigenbasis of the unperturbed Hamiltonian,
$\hat{\mathcal{H}}^{(0)}_{ab}=\epsilon_{a}\delta_{ab}$,

$(\hslash\omega^{(n)}-\epsilon_{ab})\rho^{(n)}_{ab}(\omega^{(n)})=[\hat{V}(\omega_{n}),\hat{\rho}^{(n-1)}(\omega^{(n-1)})]_{ab}$
$\implies\rho^{(n)}_{ab}(\omega^{(n)})=\frac{[\hat{V}(\omega_{n}),\hat{\rho}^{(n-1)}(\omega^{(n-1)})]_{ab}}{\hslash\omega^{(n)}-\epsilon_{ab}}\text{,}$
(6)

with $\omega^{(n)}=\sum_{j=1}^{n}\omega_{j}$ and
$\epsilon_{ab}\equiv\epsilon_{a}-\epsilon_{b}$, where $\epsilon_{a}$ are the
eigenvalues of $\hat{\mathcal{H}}^{(0)}$. Introducing a matrix
$\hat{\epsilon}$ whose entries are $\epsilon_{ab}$, (and henceforth setting
$\hslash=1$), we can compactly write the solution to the equation of motion
as,

$\hat{\rho}^{(n)}(\omega^{(n)})=\frac{1}{\omega^{(n)}-\hat{\epsilon}}\circ\left[\hat{V}(\omega_{n})\,,\,\hat{\rho}^{(n-1)}(\omega^{(n-1)})\right]\text{,}$
(7)

where $\circ$ denotes the Hadamard product, or elementwise multiplication:
$(\hat{A}\circ\hat{B})_{ab}=A_{ab}B_{ab}$. As a final bit of notation, we will
assume $\hat{V}$ is a sum of terms from which we can select a different term
$\hat{V}_{i}$ with each iteration $i$ of Eq. (7). We therefore write,

$\hat{\rho}^{(n)}_{\hat{V}_{1}\ldots\hat{V}_{n}}=\frac{1}{\omega^{(n)}-\hat{\epsilon}}\circ\left[\hat{V}_{n}\,,\,\hat{\rho}^{(n-1)}_{\hat{V}_{1}\ldots\hat{V}_{n-1}}\right]\text{,}$
(8)

dropping the implied dependence of $\hat{\rho}^{(n)}$ on $\omega^{(n)}$ as
well as the dependence of $\hat{V}_{n}$ on $\omega_{n}$. The utility of this
notation will become clear shortly.

####  II.1.1 Velocity gauge: Diagrammatic method

Let us first apply this density matrix perturbation theory to light-matter
coupling in velocity gauge. In this case,

$\hat{\mathcal{H}}_{\mathbf{k}}(t)=\hat{\mathcal{H}}^{(0)}_{\mathbf{k}+e\mathbf{A}(t)}\equiv\hat{\mathcal{H}}^{(0)}_{\mathbf{k}}+\hat{V}_{\mathbf{k}}(t)\text{,}$
(9)

which means that $\hat{V}_{\mathbf{k}}(t)$ takes the form,

$\begin{split}\hat{V}_{\mathbf{k}}(t)&=\frac{\partial\hat{\mathcal{H}}_{\mathbf{k}}^{(0)}}{\partial
k_{\alpha_{1}}}eA_{\alpha_{1}}(t)+\frac{1}{2}\frac{\partial^{2}\hat{\mathcal{H}}_{\mathbf{k}}^{(0)}}{\partial
k_{\alpha_{1}}\partial
k_{\alpha_{2}}}e^{2}A_{\alpha_{1}}(t)A_{\alpha_{2}}(t)+\ldots\\\
&=\sum_{n=1}^{\infty}\frac{e^{n}}{n!}\hat{h}^{\alpha_{1}\ldots\alpha_{n}}A_{\alpha_{1}}(t)\ldots
A_{\alpha_{n}}(t)\text{,}\end{split}$ (10)

where we introduce the shorthand
$\hat{h}^{\alpha_{1}\ldots\alpha_{n}}\equiv\partial^{\alpha_{1}}\ldots\partial^{\alpha_{n}}\hat{\mathcal{H}}^{(0)}_{\mathbf{k}}$
with $\partial^{\alpha}\equiv\frac{\partial}{\partial k_{\alpha}}$. Here and
throughout, Greek indices denote spatial directions, and repeated indices are
summed over. We can calculate the density matrix to any desired order in
$A_{\alpha}(t)$. At first order, we have,

$\hat{\rho}^{(1)}(\omega_{1})=\frac{1}{\omega_{1}-\hat{\epsilon}}\circ\left[\hat{h}^{\alpha_{1}},\hat{\rho}^{(0)}\right]A_{\alpha_{1}}(\omega_{1})\text{.}$
(11)

At second order, we have an additional term,

$\begin{split}\hat{\rho}^{(2)}(\omega^{(2)})&=\frac{1}{\omega^{(2)}-\hat{\epsilon}}\circ\left[\hat{h}^{\alpha_{2}},\hat{\rho}^{(1)}(\omega_{1})\right]A_{\alpha_{2}}(\omega_{2})\\\
&\,\,\,\,\,\,\,+\frac{1}{\omega^{(2)}-\hat{\epsilon}}\circ\left[\frac{1}{2}\hat{h}^{\alpha_{1}\alpha_{2}},\hat{\rho}^{(0)}\right]A_{\alpha_{1}}(\omega_{1})A_{\alpha_{2}}(\omega_{2})\\\
&\equiv\left(\hat{\rho}^{(2)}_{\hat{h}^{\alpha_{1}},\hat{h}^{\alpha_{2}}}(\omega^{(2)})+\hat{\rho}^{(1)}_{\frac{1}{2}\hat{h}^{\alpha_{1}\alpha_{2}}}(\omega^{(2)})\right)A_{\alpha_{1}}(\omega_{1})A_{\alpha_{2}}(\omega_{2})\text{.}\end{split}$
(12)

Hence, the utility of the notation in Eq. (8) is to split the density matrix
at $n$-th order in $A_{\alpha}(t)$ into a sum of density matrices
$\hat{\rho}^{(m)}_{\hat{V}_{1}\ldots\hat{V}_{m}}$ that are perturbed with
relatively simple operators $\hat{V}_{i}$. In fact, each of these operators
$\hat{V}_{i}$ is simply a $j$-th $\mathbf{k}$-derivative (for some $j$) of the
bare Hamiltonian $\hat{\mathcal{H}}^{(0)}_{\mathbf{k}}$, which has the
physically transparent meaning of a vertex interaction with $j$ photons. (A
word of caution: the “$(n)$” in
$\hat{\rho}^{(n)}_{\hat{V}_{1}\ldots\hat{V}_{n}}$ no longer means “at $n$-th
order in $A_{\alpha}(t)$” (as opposed to the more conventional notation
$\hat{\rho}^{(n)}$) and now simply means “perturbed with $n$ operators
$\hat{V}_{i}$, $i=1,\ldots,n$.”)

This decomposition is the basis for the connection with the diagrammatic
scheme in Ref. Parker  _et al._ (2019). At zeroth order in $A_{\alpha}(t)$,
the expectation value of the current is simply a trace of the current operator
with respect to the equilibrium density matrix,

$\left\langle\hat{j}^{\mu}\right\rangle_{0}=\text{tr}\left\\{\hat{j}_{(0)}^{\mu}\cdot\hat{\rho}^{(0)}\right\\}\text{.}$
(13)

At nonzero order in $A_{\alpha}(t)$, both the density matrix and the current
operator must be expanded in powers of $A_{\alpha}(t)$,

$\begin{split}\hat{j}^{\mu}(t)&\equiv-\frac{\partial\hat{\mathcal{H}}_{\mathbf{k}}(t)}{\partial
A_{\mu}(t)}\\\
&=-\sum_{n=0}^{\infty}\frac{e^{n+1}}{n!}\hat{h}^{\mu\alpha_{1}\ldots\alpha_{n}}A_{\alpha_{1}}(t)\ldots
A_{\alpha_{n}}(t)\text{,}\end{split}$ (14)
$\hat{\rho}(t)=\sum_{n=0}^{\infty}\hat{\rho}^{(n)}(t)\text{.}$ (15)

Keeping only terms at some desired order in $A_{\alpha}(t)$, the result is
always a sum of traces, where each trace (expressible by a diagram) is acting
on a product between some term in Eq. (14) and a density matrix of the form
$\hat{\rho}^{(n)}_{\hat{V}_{1}\ldots\hat{V}_{n}}$. We express this
diagrammatically in terms of “Feynman rules” illustrated in Figure 1. This
approach can be viewed as a hybrid approach making connection between the
density matrix equation of motion method of Refs. Ventura  _et al._ (2017);
Passos  _et al._ (2018) and the Matsubara Green’s function approach of Ref.
Parker  _et al._ (2019), which lead to the same results. That is, one can
alternatively view the calculation of
$\hat{\rho}^{(n)}_{\hat{V}_{1}\cdots\hat{V}_{n}}$ in a particular diagram
using a trace of a product of Matsubara Green’s functions multiplied by vertex
operators, summed over a single unconstrained fermionic Matsubara frequency.

Figure 1: Summary of Feynman rules employed for computing optical responses in
velocity gauge. Following the convention of Parker  _et al._ (2019), the
output vertex is marked with a “$\otimes$,” and represents some output current
operator $\hat{\mathcal{O}}$. Measuring this output current amounts to a trace
of this operator with a density matrix, which is represented by a fermion
loop. This density matrix is perturbed (as in Eq. (8)) with operators
$\hat{V}_{i}$ represented by black vertices. The vertex operators take the
form of $\mathbf{k}$-derivatives of the Hamiltonian (denoted
$\hat{h}^{\alpha_{1}\ldots\alpha_{n}}$) with the directions of the derivatives
given by the polarizations (denoted with Greek letters) of the photons
entering the vertex. The multi-photon vertex operators are symmetrized with a
factor $1/n!$, where $n$ is the number of photons excluding the photon
representing the output current (we reserve the index $\mu$ to indicate the
polarization of the output photon).

The last step is to convert the response in terms of
$A_{\alpha_{j}}(\omega_{j})$ to a response in terms of
$E_{\alpha_{j}}(\omega_{j})\equiv i\omega_{j}A_{\alpha_{j}}(\omega_{j})$. We
formally define,

$\begin{split}\left\langle
j^{\mu}(\omega^{(n)})\right\rangle=\,\,&\sigma^{\mu\alpha_{1}\ldots\alpha_{n}}(\omega^{(n)};\omega_{1},\ldots,\omega_{n})\\\
&\times E_{\alpha_{1}}(\omega_{1})\ldots
E_{\alpha_{n}}(\omega_{n})\text{.}\end{split}$ (16)

Recalling $\omega^{(n)}\equiv\sum_{j=1}^{n}\omega_{j}$, this is the fully
general nonlinear conductivity at frequency $\omega^{(n)}$ due to input
electric fields at frequencies $\\{\omega_{j}\\}_{j=1,\ldots,n}$. We calculate
this by first calculating a quantity which we will call
$\kappa^{\mu\alpha_{1}\ldots\alpha_{n}}_{\mathbf{k}}(\omega^{(n)};\\{\omega_{j}\\})$

$\begin{split}\sigma^{\mu\alpha_{1}\ldots\alpha_{n}}(\omega^{(n)};\\{\omega_{j}\\})=\,\,&-\frac{e^{n+1}}{\hslash^{n}}\prod_{j=1}^{n}\left(\frac{1}{i\omega_{j}}\right)\\\
&\times\int[d\mathbf{k}]\,\kappa_{\mathbf{k}}^{\mu\alpha_{1}\ldots\alpha_{n}}(\omega^{(n)};\\{\omega_{j}\\})\text{.}\end{split}$
(17)

The task is therefore reduced to calculating
$\kappa_{\mathbf{k}}^{\mu\alpha_{1}\ldots\alpha_{n}}$ via Feynman diagrams. We
will do this for the example cases of the linear conductivity and SHG
conductivity, and note the cancellations that occur to eliminate the apparent
$1/\omega$ and $1/\omega^{2}$ divergences, respectively.

####  II.1.2 Length gauge

In length gauge, one simply has the perturbation,

$\hat{V}(t)=e\,\hat{r}^{\alpha}E_{\alpha}(t)\text{,}$ (18)

where $\hat{r}^{\alpha}$ denotes the position operator. Care needs to be taken
when interpreting the $\mathbf{k}$-space form of the operator
$\hat{r}^{\alpha}$ Blount (1962). When calculating commutators with other
operators, $\hat{r}^{\alpha}$ acts as a covariant derivative,

$\left[\hat{r}^{\alpha},\hat{\mathcal{O}}\right]_{ab}\equiv
i\left[\hat{\mathcal{D}}^{\alpha},\hat{\mathcal{O}}\right]_{ab}=i\partial^{\alpha}\hat{\mathcal{O}}_{ab}+\left[\hat{\mathcal{A}}^{\alpha},\hat{\mathcal{O}}\right]_{ab}\text{,}$
(19)

where
$\mathcal{A}^{\alpha}_{ab}(\mathbf{k})\equiv\braket{u_{a}(\mathbf{k})}{\partial^{\alpha}u_{b}(\mathbf{k})}$
is the non-Abelian Berry connection. With this definition, nonlinear
conductivity tensors in length gauge can be computed directly as (suppressing
frequency dependence),

$\sigma^{\mu\alpha_{1}\ldots\alpha_{n}}=-\frac{e^{n+1}}{\hslash^{n}}\int[d\mathbf{k}]\,\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(n)}_{\hat{r}^{\alpha_{1}}\ldots\hat{r}^{\alpha_{n}}}\right\\}\text{,}$
(20)

with,

$\hat{\rho}^{(n)}_{\hat{r}^{\alpha_{1}}\ldots\hat{r}^{\alpha_{n}}}=\frac{1}{\omega^{(n)}-\hat{\epsilon}}\circ\left[\hat{r}^{\alpha_{n}},\hat{\rho}^{(n-1)}_{\hat{r}^{\alpha_{1}}\ldots\hat{r}^{\alpha_{n-1}}}\right]\text{,}$
(21)

without need for extra terms as in the diagrammatic scheme for velocity gauge.
While the length gauge expression is analytically more compact, it is
numerically less straightforward to use, since $\mathbf{k}$-derivatives act
iteratively on the density matrix. In velocity gauge, all the
$\mathbf{k}$-derivatives act on the Hamiltonian, so they can be computed ahead
of time to speed up calculations.

###  II.2 Examples

####  II.2.1 Linear conductivity

Figure 2: Diagrams relevant to computing the linear conductivity tensor
$\sigma^{\mu\alpha}(\omega)$ in velocity gauge.

The linear conductivity
$\sigma^{\mu\alpha}(\omega;\omega)\equiv\sigma^{\mu\alpha}(\omega)$ is defined
as follows,

$\left\langle\hat{j}^{\mu}(\omega)\right\rangle=\sigma^{\mu\alpha}(\omega)E_{\alpha}(\omega)\text{.}$
(22)

The calculation in velocity gauge can be written as the sum of two terms,
represented diagrammatically in Figure 2,

$\sigma^{\mu\alpha}(\omega)=-\frac{e^{2}}{\hslash}\left(\frac{1}{i\omega}\right)\int[d\mathbf{k}]\,\kappa^{\mu\alpha}_{\mathbf{k}}(\omega)\text{,}$
$\kappa_{\mathbf{k}}^{\mu\alpha}(\omega)=\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(1)}_{\hat{v}^{\alpha}}(\omega)\right\\}+\text{tr}\left\\{\hat{h}^{\mu\alpha}\cdot\hat{\rho}^{(0)}\right\\}\text{,}$
(23)

where each operator within the traces is understood to be evaluated at
$\mathbf{k}$. Note that $\hat{v}^{\mu}\equiv\hat{h}^{\mu}$, so we will use
these interchangeably. Physically, these two terms correspond to the well-
known paramagnetic current and diamagnetic current, respectively. The apparent
$1/\omega$ divergence in the conductivity formula is canceled when these two
terms are summed and integrated over k. It is interesting to note that this
cancellation does not occur in the case of a superconductor, where the low-
frequency divergence of the conductivity is physical and results in the
linear-in-$\mathbf{A}$ current in the London equation.

Assuming $\rho^{(0)}_{ab}=f_{a}\delta_{ab}$, we can re-express the linear
conductivity as follows,

$\sigma^{\mu\alpha}(\omega)=\frac{ie^{2}}{\hslash\omega}\int[d\mathbf{k}]\left(\sum_{ab}f_{ab}\frac{v^{\mu}_{ab}v^{\alpha}_{ba}}{\omega-\epsilon_{ba}}+\sum_{a}f_{a}h^{\mu\alpha}_{aa}\right)\text{,}$
(24)

which matches velocity gauge results from the literature Passos  _et al._
(2018).

In length gauge, recalling that
$[\hat{r}^{\alpha},\hat{\mathcal{O}}]_{ab}=i\partial^{\alpha}\mathcal{O}_{ab}+[\hat{\mathcal{A}}^{\alpha},\hat{\mathcal{O}}]_{ab}$,
where $\hat{\mathcal{A}}^{\alpha}$ is the non-Abelian Berry connection, we
write,

$\begin{split}\sigma^{\mu\alpha}_{\text{len.}}(\omega)&=-\frac{e^{2}}{\hslash}\int[d\mathbf{k}]\,\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(1)}_{\hat{r}^{\alpha}}(\omega)\right\\}\\\
&=-\frac{e^{2}}{\hslash}\int[d\mathbf{k}]\left(\sum_{ab}f_{ab}\frac{v^{\mu}_{ab}\mathcal{A}^{\alpha}_{ba}}{\omega-\epsilon_{ba}}+\sum_{a}\frac{iv^{\mu}_{aa}\partial^{\alpha}f_{a}}{\omega}\right)\text{.}\end{split}$
(25)

Using $v^{\alpha}_{ab}=i\epsilon_{ab}\mathcal{A}^{\alpha}_{ab}$ for $a\neq b$,
this allows us to write,

$\begin{split}\sigma^{\mu\alpha}_{\text{len.}}(\omega)=\frac{ie^{2}}{\hslash}\int[d\mathbf{k}]&\left(\sum_{ab}f_{ab}\frac{v^{\mu}_{ab}v^{\alpha}_{ba}}{\epsilon_{ba}(\omega-\epsilon_{ba})}\right.\\\
&\left.-\sum_{a}\frac{v^{\mu}_{aa}\partial^{\alpha}f_{a}}{\omega}\right)\text{.}\end{split}$
(26)

####  II.2.2 Second harmonic generation conductivity

Figure 3: Diagrams relevant to computing the SHG conductivity tensor
$\sigma_{\text{SHG}}^{\mu\alpha\beta}(\omega)$ in velocity gauge.

The second harmonic generation conductivity
$\sigma^{\mu\alpha\beta}(2\omega;\omega,\omega)\equiv\sigma^{\mu\alpha\beta}_{\text{SHG}}(\omega)$
is defined as follows,

$\left\langle\hat{j}^{\mu}(2\omega)\right\rangle=\sigma^{\mu\alpha\beta}_{\text{SHG}}(\omega)E_{\alpha}(\omega)E_{\beta}(\omega)$
(27)

The calculation in velocity gauge can be written as the sum of four terms,
represented diagrammatically in Figure 3,

$\sigma^{\mu\alpha\beta}_{\text{SHG}}(\omega)=-\frac{e^{3}}{\hslash^{2}}\left(\frac{-1}{\omega^{2}}\right)\int[d\mathbf{k}]\,\kappa^{\mu\alpha\beta}_{\text{SHG},\mathbf{k}}(\omega)\text{,}$
$\begin{split}\kappa^{\mu\alpha\beta}_{\text{SHG},\mathbf{k}}(\omega)=\,\,&\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(2)}_{\hat{v}^{\alpha},\hat{v}^{\beta}}(2\omega)\right\\}+\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(1)}_{\frac{1}{2}\hat{h}^{\mu\alpha}}(2\omega)\right\\}\\\
&+\text{tr}\left\\{\hat{h}^{\mu\alpha}\cdot\hat{\rho}^{(1)}_{\hat{v}^{\beta}}(\omega)\right\\}+\text{tr}\left\\{\frac{1}{2}\hat{h}^{\mu\alpha\beta}\cdot\hat{\rho}^{(0)}\right\\}\text{.}\end{split}$
(28)

Assuming $\rho^{(0)}_{ab}=f_{a}\delta_{ab}$, this becomes,

$\begin{split}\kappa^{\mu\alpha\beta}_{\text{SHG},\mathbf{k}}(\omega)=\,\,&\sum_{a,b,c}v^{\mu}_{ab}\frac{\displaystyle
v^{\beta}_{bc}\frac{v^{\alpha}_{ca}f_{ac}}{\omega-\epsilon_{ca}}-\frac{f_{cb}v^{\alpha}_{bc}}{\omega-\epsilon_{bc}}v^{\beta}_{ca}}{2\omega-\epsilon_{ba}}\\\
&+\sum_{a,b}v^{\mu}_{ab}\frac{\frac{1}{2}h^{\alpha\beta}_{ba}f_{ab}}{2\omega-\epsilon_{ab}}+\sum_{a,b}h^{\mu\alpha}_{ab}\frac{v^{\beta}_{ba}f_{ab}}{\omega-\epsilon_{ba}}\\\
&+\sum_{a}\frac{1}{2}h^{\mu\alpha\beta}_{aa}f_{a}\text{.}\end{split}$ (29)

We note that, under $\mathbf{k}$-space symmetry considerations alone, the
contributions from the first and last diagrams in Figure 3 should vanish under
$\mathbf{k}$-space integration over the full Brillouin zone (that is, their
integrands are odd in $\mathbf{k}$ while the integration region is symmetric
in $\mathbf{k}$). Therefore, the content of the cancellation required to
eliminate the $1/\omega^{2}$ prefactor in the conductivity comes solely from
the second and third diagrams, just as the cancellation in the linear
conductivity comes from only two diagrams.

##  III Calculations of spatially dispersive nonlinear optical responses

###  III.1 Generalities

A nonlinear optical conductivity tensor is in full generality a function of
the wavevectors of the input/output optical fields $\mathbf{q}_{i}$, but is
typically treated in the $\mathbf{q}\to 0$ limit. This is justified within the
so-called dipole approximation, in which it is assumed that the length scale
of the electronic states being coupled by optical transitions (for instance,
the size of electronic orbitals) is very small compared to the wavelength of
the field driving the transition. However, there are many cases where this
approximation should be expected to fail (for instance, in band systems whose
Wannier representation has poor localization). For these reasons, one may be
interested in the lowest-order contribution of the spatial variation of the
optical fields,

$E_{\alpha_{j}}(\mathbf{r},t)=E_{\alpha_{j}}(\mathbf{0},t)\left(1-i\mathbf{q}_{j}\cdot\hat{\mathbf{r}}\right)+\mathcal{O}(\mathbf{q}^{2})\text{.}$
(30)

Previous works Malashevich and Souza (2010); Ahn  _et al._ (2022b) have
identified how to expand the linear optical conductivity to first-order in
$\mathbf{q}$. Interpreting the zeroth-order contribution to the current
operator as the time-derivative of the electric dipole,
$\hat{P}^{\mu}=-e\hat{r}^{\mu}$, the first-order terms include contributions
from the electric quadrupole and magnetic dipole, Melrose and McPhedran (2009)

$\hat{j}^{\mu}=\partial_{t}\hat{P}^{\mu}-iq_{\nu}\left(\frac{1}{2}\partial_{t}\hat{Q}^{\nu\mu}+c\epsilon^{\nu\mu\rho}\hat{M}_{\rho}\right)+\mathcal{O}(\mathbf{q}^{2})\text{,}$
(31)

where $\hat{Q}^{\nu\mu}=-e\hat{r}^{\nu}\hat{r}^{\mu}$ is the electric
quadrupole operator, and
$\hat{M}_{\rho}=-\frac{e}{2c}\epsilon_{\rho\alpha\beta}\hat{r}^{\alpha}\hat{j}^{\beta}$
is the magnetic dipole operator. Eq. (31) constitutes the starting point for a
velocity gauge form for light-matter coupling up to first order in
$\mathbf{q}$. In this work, we will focus on the electric quadrupole
contribution, and determine the velocity gauge and length gauge expressions
for treating this contribution. In length gauge, the perturbation takes the
form,

$\hat{V}(\omega)=\left(\hat{r}^{\alpha}-iq_{\nu}\left(\frac{1}{2}\hat{Q}^{\nu\alpha}+\frac{c}{\omega}\epsilon^{\nu\mu\rho}\hat{M}_{\rho}\right)\right)E_{\alpha}(\omega)$
(32)

We constrain the exact form of the electric quadrupole operator by ensuring
that we recover $-\frac{\partial}{\partial
r^{\alpha}}V(\omega)=E_{\alpha}(\omega)(1-iq_{\nu}r^{\nu})+\mathcal{O}(\mathbf{q}^{2})$.
Notably, this gives a relative factor of 2 to the off-diagonal spatial
components of $\hat{Q}^{\nu\mu}$. For example, in $d=2$, we have,

$\hat{Q}^{\nu\mu}=\begin{pmatrix}\hat{Q}^{xx}&\hat{Q}^{xy}\\\
\hat{Q}^{yx}&\hat{Q}^{yy}\end{pmatrix}=\begin{pmatrix}\hat{\mathcal{A}}^{x}\cdot\hat{\mathcal{A}}^{x}&2\hat{\mathcal{A}}^{x}\cdot\hat{\mathcal{A}}^{y}\\\
2\hat{\mathcal{A}}^{y}\cdot\hat{\mathcal{A}}^{x}&\hat{\mathcal{A}}^{y}\cdot\hat{\mathcal{A}}^{y}\end{pmatrix}\text{,}$
(33)

where we recall the non-Abelian Berry connection
$\mathcal{A}^{\mu}_{ab}\equiv\left\langle
u_{a\mathbf{k}}\right.\left|\partial^{\mu}u_{b\mathbf{k}}\right\rangle$.

It has been shown Malashevich and Souza (2010) that the first-order-
in-$\mathbf{q}$ contribution to the current operator matrix element can be
cast in terms of the following object,

$\begin{split}\left\langle
u_{n,\mathbf{k}+\frac{\mathbf{q}}{2}}\right|\hat{v}^{\mu}_{\mathbf{k}}\left|u_{m,\mathbf{k}-\frac{\mathbf{q}}{2}}\right\rangle&=-\frac{i}{2}q_{\nu}\left\\{\hat{A}^{\nu}_{\mathbf{k}},\hat{v}^{\mu}_{\mathbf{k}}\right\\}_{nm}\\\
&\equiv-\frac{i}{2}q_{\nu}\,\left(\partial_{t}\hat{Q}^{\nu\mu}_{\mathbf{k}}\right)_{nm}\end{split}$
(34)

where in the last line we make the crucial identification that this should
equal the time derivative of the electric quadrupole moment. This is the basis
for a new vertex rule we introduce, illustrated in Figure 4 and summarized in
the following expression for the velocity-gauge perturbation,

$\begin{split}\hat{V}_{\mathbf{k}}(t)=&\sum_{n=1}^{\infty}\frac{e^{n}}{n!}A_{\alpha_{1}}(t)\ldots
A_{\alpha_{n}}(t)\\\
&\times\Bigg{(}\hat{\mathcal{D}}^{\alpha_{1}\ldots\alpha_{n}}_{\mathbf{k}}\left[\hat{\mathcal{H}}_{\mathbf{k}}\right]\\\
&+\frac{1}{2}\sum_{j=1}^{n}(q_{j})_{\nu}\hat{\mathcal{D}}_{\mathbf{k}}^{\alpha_{1}\ldots\alpha_{j-1}\alpha_{j+1}\ldots\alpha_{n}}\left[\hat{\mathcal{H}}_{\mathbf{k}}\right.,\left.\hat{Q}^{\nu\alpha_{j}}_{\mathbf{k}}\right]\Bigg{)}\text{,}\end{split}$
(35)

where the vectors $\mathbf{q}_{j}$ are the wavevectors from the $n$ fields
$A_{\alpha_{1}}(t),\ldots,A_{\alpha_{n}}(t)$. We can interpret this as
introducing an additional type of vertex in the diagrammatic method, one that
takes into account at first order the momentum ($\hslash\mathbf{q}_{j}$) of
each photon connected to that vertex. We will use the terms “momentum” and
“wavevector” interchangeably.

Schematically, this leads to a spatially dispersive correction to nonlinear
optical conductivities,

$\sigma^{\mu\alpha_{1}\ldots\alpha_{n}}(\mathbf{q})=\sigma^{\mu\alpha_{1}\ldots\alpha_{n}}_{(0)}+q_{\nu}\sigma^{\nu\mu\alpha_{1}\ldots\alpha_{n}}_{(1)}+\mathcal{O}(\mathbf{q}^{2})\text{,}$
(36)

where $\mathbf{q}$ is some characteristic momentum of the input fields. To
declutter notation, we have suppressed the $\omega$ dependence, which is
understood to be $(\omega_{(n)};\omega_{1}\ldots\omega_{n})$ in general. As
before, in velocity gauge, we write this as,

$\sigma^{\nu\mu\alpha_{1}\ldots\alpha_{n}}_{(1)}=-\frac{e^{n+1}}{\hslash^{n}}\prod_{j=1}^{n}\left(\frac{1}{i\omega_{j}}\right)\int[d\mathbf{k}]\,\kappa_{(1),\mathbf{k}}^{\nu\mu\alpha_{1}\ldots\alpha_{n}}\text{,}$
(37)

and calculate $\kappa^{\nu\mu\alpha_{1}\ldots\alpha_{n}}_{(1),\mathbf{k}}$
using a diagrammatic scheme. As a consequence of Eq. (35), this diagrammatic
scheme amounts to taking each diagram for the $\mathcal{O}(\mathbf{q}^{0})$
response and, for each vertex $\hat{h}^{\alpha_{1}\ldots\alpha_{n}}$, creating
a new diagram with the following replacement rule (dropping $\mathbf{k}$
dependence for brevity),

$\begin{split}\hat{h}^{\alpha_{1}\ldots\alpha_{n}}&\to\frac{1}{2}\sum_{j=1}^{n}(q_{j})_{\nu}\hat{\mathcal{D}}^{\alpha_{1}\ldots\alpha_{j-1}\alpha_{j+1}\ldots\alpha_{n}}\left[\hat{\mathcal{H}}\right.,\left.\hat{Q}^{\nu\alpha_{j}}\right]\\\
&\equiv-\frac{i}{2}\sum_{j=1}^{n}(q_{j})_{\nu}\partial_{t}\hat{Q}^{\nu\alpha_{j};\alpha_{1}\ldots\alpha_{j-1}\alpha_{j+1}\ldots\alpha_{n}}\text{,}\end{split}$
(38)

where $(q_{j})_{\nu}$ is the $\nu$ component of the wavevector for the photon
with polarization $\alpha_{j}$ (again, summation over the repeated index $\nu$
is implied). In the second line, we introduce a compact notation, whereby
indices after the semicolon indicate covariant $\mathbf{k}$-derivatives. To
make sense of Eq. (38), note that for a single-photon vertex, the replacement
rule reduces to,

$\hat{h}^{\alpha_{1}}\to-\frac{i}{2}(q_{1})_{\nu}\partial_{t}\hat{Q}^{\nu\alpha_{1}}\text{.}$
(39)

In other words, it replaces an electric dipole coupling
$\hat{h}^{\alpha_{1}}\equiv\hat{v}^{\alpha_{1}}\equiv\partial_{t}\hat{r}^{\alpha_{1}}$
with an electric quadrupole coupling $\partial_{t}\hat{Q}^{\nu\alpha_{1}}$.
The extension to multi-photon vertices can then be understood as a sum of
terms that change the coupling from dipolar to quadrupolar for one photon at a
time while having the rest of the photons contribute covariant
$\mathbf{k}$-derivatives of the resulting quadrupolar current operator. We
illustrate this new vertex rule in Figure 4.

We note that the output vertex should also be perturbed to linear-
in-$\mathbf{q}$ order in this scheme, and we denote this by a “$\boxtimes$”
symbol in the diagrams. As in the spatially uniform case, the only special
treatment that must be given to the output vertex is that the symmetry factor
$1/n!$ should only include the $n$ input photons (i.e. excluding the output
photon, for whom we reserve the index $\mu$ to denote its polarization).

Figure 4: Rule for determining the first-order-in-$\mathbf{q}$ contribution
from a vertex operator. We change a vertex from a dot to a square to indicate
that it is being perturbed to first order in $\mathbf{q}$.

###  III.2 Examples

####  III.2.1 Spatially dispersive linear conductivity

Applying the Feynman rules in Figure 4 to the linear conductivity diagrams in
Figure 2, we obtain the diagrams in Figure 5. These diagrams represent the
following two terms,

Figure 5: Diagrams relevant to computing the spatially dispersive linear
conductivity tensor to lowest order ($\sigma_{(1)}^{\nu\mu\alpha}(\omega)$) in
velocity gauge. The third diagram is shown for completeness, but has a
vanishing contribution since it only contains a vertex with net-zero momentum
flow. (We omit a common factor of $-\frac{i}{2}$ from each diagram.)

$\begin{split}\kappa_{(1),\mathbf{k}}^{\nu\mu\alpha}(\omega)=-\frac{i}{2}&\Big{(}\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(1)}_{\partial_{t}\hat{Q}^{\nu\alpha}}(\omega)\right\\}\\\
&+\text{tr}\left\\{\partial_{t}\hat{Q}^{\nu\mu}\cdot\hat{\rho}^{(1)}_{\hat{v}^{\alpha}}(\omega)\right\\}\\\
&+\text{tr}\left\\{\left(\partial_{t}\hat{Q}^{\nu\mu;\alpha}+\partial_{t}\hat{Q}^{\nu\alpha;\mu}\right)\cdot\hat{\rho}^{(0)}\right\\}\Big{)}\text{.}\end{split}$
(40)

The conductivity can therefore be written as,

$\begin{split}\sigma^{\nu\mu\alpha}_{(1)}(\omega)=\frac{e^{2}}{2\hslash\omega}\int[d\mathbf{k}]&\Bigg{\\{}\sum_{ab}f_{ab}\frac{v^{\mu}_{ab}\dot{Q}^{\nu\alpha}_{ba}+\dot{Q}^{\nu\mu}_{ab}v^{\alpha}_{ba}}{\omega-\epsilon_{ba}}\\\
&+\sum_{a}f_{a}\left(\dot{Q}^{\nu\mu;\alpha}_{aa}+\dot{Q}^{\nu\alpha;\mu}_{aa}\right)\Bigg{\\}}\text{,}\end{split}$
(41)

where we abbreviate $\left\langle
u_{a\mathbf{k}}\right|\partial_{t}\hat{Q}^{\nu\mu}\left|u_{b\mathbf{k}}\right\rangle\equiv\dot{Q}^{\nu\mu}_{ab}\equiv
i\left[\hat{\mathcal{H}}\right.,\left.\hat{Q}^{\nu\mu}\right]_{ab}$. Despite
the apparent divergence from the prefactor $1/\omega$, this response is
regularized in the low-frequency limit, as we show numerically for an example
model in Section IV. To further justify this, we can re-derive this result in
length gauge,

$\begin{split}\sigma^{\nu\mu\alpha}_{(1)}(\omega)=\frac{ie^{2}}{2\hslash}\int[d\mathbf{k}]&\left(\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(1)}_{\hat{Q}^{\nu\alpha}}(\omega)\right\\}\right.\\\
&\left.+\text{tr}\left\\{\partial_{t}\hat{Q}^{\nu\mu}\cdot\hat{\rho}^{(1)}_{\hat{r}^{\alpha}}(\omega)\right\\}\right)\text{.}\end{split}$
(42)

We show numerically in Section IV that Eqs. (41) and (42) are equivalent.

####  III.2.2 Spatially dispersive second harmonic generation conductivity

Figure 6 displays the eight diagrams for calculating the linear-
in-$\mathbf{q}$ SHG conductivity. These yield the following seven terms,

Figure 6: Diagrams for computing lowest-order spatially dispersive corrections
to the SHG conductivity tensor
($\sigma^{\nu\mu\alpha\beta}_{\text{SHG},(1)}(\omega)$) in velocity gauge.

$\begin{split}\kappa^{\nu\mu\alpha\beta}_{\text{SHG},(1),\mathbf{k}}(\omega)&=\,-\frac{i}{2}\Bigg{[}\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(2)}_{\hat{v}^{\alpha},\partial_{t}\hat{Q}^{\nu\beta}}(2\omega)\right\\}\\\
&+\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(2)}_{\partial_{t}\hat{Q}^{\nu\alpha},\hat{v}^{\beta}}(2\omega)\right\\}\\\
&+\text{tr}\left\\{\partial_{t}\hat{Q}^{\nu\mu}\cdot\hat{\rho}^{(2)}_{\hat{v}^{\alpha},\hat{v}^{\beta}}(2\omega)\right\\}\\\
&+\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(1)}_{\frac{1}{2}\partial_{t}\hat{Q}^{\nu\alpha;\beta}+\frac{1}{2}\partial_{t}\hat{Q}^{\nu\beta;\alpha}}(2\omega)\right\\}\\\
&+\text{tr}\left\\{\partial_{t}\hat{Q}^{\nu\mu}\cdot\hat{\rho}^{(1)}_{\frac{1}{2}\hat{h}^{\alpha\beta}}(2\omega)\right\\}\\\
&+\text{tr}\left\\{\hat{h}^{\mu\alpha}\cdot\hat{\rho}^{(1)}_{\partial_{t}\hat{Q}^{\nu\beta}}(\omega)\right\\}\\\
&+\text{tr}\left\\{\left(\partial_{t}\hat{Q}^{\nu\mu;\alpha}+\partial_{t}\hat{Q}^{\nu\alpha;\mu}\right)\cdot\hat{\rho}^{(1)}_{\hat{v}^{\beta}}(\omega)\right\\}\\\
&+\text{tr}\left\\{\frac{1}{2}\left(\partial_{t}\hat{Q}^{\nu\mu;\alpha\beta}+\text{[perms.]}\right)\cdot\hat{\rho}^{(0)}\right\\}\Bigg{]}\text{,}\end{split}$
(43)

where “[perms.]” indicates the two other permutations of $\mu,\alpha,\beta$
modulo exchanges of the indices after the semicolon. Assuming
$\rho^{(0)}_{ab}=f_{a}\delta_{ab}$, this becomes,

$\begin{split}\kappa^{\nu\mu\alpha\beta}_{\text{SHG},(1),\mathbf{k}}(\omega)&=-\frac{i}{2}\sum_{a,b}\Bigg{(}v^{\mu}_{ab}\frac{N_{1,ba}^{\nu\alpha\beta}(\omega)}{2\omega-\epsilon_{ba}}-\dot{Q}^{\nu\mu}_{ab}\frac{N_{2,ba}^{\alpha\beta}(\omega)}{2\omega-\epsilon_{ba}}\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\frac{1}{2}v^{\mu}_{ab}\frac{(\dot{Q}^{\nu\alpha;\beta}_{ba}+\dot{Q}^{\nu\beta;\alpha}_{ba})f_{ab}}{2\omega-\epsilon_{ab}}\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\frac{1}{2}\dot{Q}^{\nu\mu}_{ab}\frac{h^{\alpha\beta}_{ba}f_{ab}}{2\omega-\epsilon_{ab}}+h^{\mu\alpha}_{ab}\frac{\dot{Q}^{\nu\beta}_{ba}f_{ab}}{\omega-\epsilon_{ba}}\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\left(\dot{Q}^{\nu\mu;\alpha}_{ab}+\dot{Q}^{\nu\alpha;\mu}_{ab}\right)\frac{v^{\beta}_{ba}f_{ab}}{\omega-\epsilon_{ba}}\Bigg{)}\\\
&\,\,\,\,\,-\frac{i}{4}\sum_{a}\left(\dot{Q}^{\nu\mu;\alpha\beta}_{aa}+\dot{Q}^{\nu\alpha;\beta\mu}_{aa}+\dot{Q}^{\nu\beta;\mu\alpha}_{aa}\right)f_{a}\text{,}\end{split}$
(44)

where,

$\begin{split}N_{1,ba}^{\nu\alpha\beta}(\omega)&=\sum_{c}\left(\dot{Q}^{\nu\beta}_{bc}\frac{v^{\alpha}_{ca}f_{ac}}{\omega-\epsilon_{ca}}-\frac{f_{cb}v^{\alpha}_{bc}}{\omega-\epsilon_{bc}}\dot{Q}^{\nu\beta}_{ca}\right.\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\left.v^{\beta}_{bc}\frac{\dot{Q}^{\nu\alpha}_{ca}f_{ac}}{\omega-\epsilon_{ca}}-\frac{f_{cb}\dot{Q}^{\nu\alpha}_{bc}}{\omega-\epsilon_{bc}}v^{\beta}_{ca}\right)\text{,}\\\
N_{2,ba}^{\alpha\beta}(\omega)&=\sum_{c}\left(v^{\beta}_{bc}\frac{v^{\alpha}_{ca}f_{ac}}{\omega-\epsilon_{ca}}-\frac{f_{cb}v^{\alpha}_{bc}}{\omega-\epsilon_{bc}}v^{\beta}_{ca}\right)\text{.}\end{split}$
(45)

In our forthcoming results, we compare to the following compact length gauge
formula,

$\begin{split}\sigma^{\nu\mu\alpha\beta}_{(1)}(\omega)=\frac{ie^{3}}{2\hslash^{2}}\int[d\mathbf{k}]&\Bigg{(}\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(2)}_{\hat{Q}^{\nu\alpha},\hat{r}^{\beta}}\right\\}\\\
&+\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(2)}_{\hat{r}^{\alpha},\hat{Q}^{\nu\beta}}\right\\}\\\
&+\text{tr}\left\\{\partial_{t}\hat{Q}^{\nu\mu}\cdot\hat{\rho}^{(2)}_{\hat{r}^{\alpha},\hat{r}^{\beta}}\right\\}\Bigg{)}\text{.}\end{split}$
(46)

##  IV Results for minimal models

Figure 7: Results for the spatially dispersive correction to the linear
conductivity ($\sigma^{\nu\mu\alpha}_{(1)}(\omega)$) calculated on a minimal
quasi-2D model. This model takes the form of a chain of rectangular molecules
with inter-molecule couplings, with a sensible “decoupling limit” in which the
problem reduces to the well-understood problem of a multipole response in a
molecule. (Top) Agreement is demonstrated for the molecular chain between our
scheme (in both length and velocity gauge) and the expression determined in
Ref. Malashevich and Souza (2010) from the Kubo formula expanded to first
order in $\mathbf{q}$. (Bottom) Agreement is also demonstrated in the
decoupling limit between our scheme (in both length and velocity gauge) and
the electric quadrupole linear response from an isolated molecule. Note that
inversion-breaking is needed to observe this response, so we introduce a small
static on-site potential ($\Delta$) of $\pm 0.1$ on sublattices $AD$/$BC$ to
produce these results. Figure 8: Results for the spatially dispersive
correction to the second harmonic generation conductivity
($\sigma^{\nu\mu\alpha\beta}_{\text{SHG},(1)}(\omega)$) calculated on the
molecular chain model detailed in Figure 7. Agreement is shown between length
gauge and velocity gauge in both the molecular chain and in the decoupling
limit. Agreement is also shown with the analogous single-molecule calculation.
We use $t_{x}^{\prime}=0.8$

In this section, we test our scheme on a minimal tight-binding model with
three useful features: (1) inversion symmetry, (2) a quasi-2D structure
allowing for multiple different polarizations, and (3) a clear “decoupling
limit” in which the model reduces to the well-understood problem of multipole
optical responses in a molecule Cohen-Tannoudji  _et al._ (1997); Melrose and
McPhedran (2009); Jaszuński  _et al._ (2017). We will confirm that, up to
electric quadrupole order, a nonzero SHG response indeed exists in this
centrosymmetric system, and we can therefore treat $\mathbf{q}$ on the same
footing as a static inversion-breaking parameter in this model.

###  IV.1 Rectangular molecule

We first establish the response up to first order in $\mathbf{q}$ and up to
second order in the field strength for a 4-site rectangular molecule with only
nearest-neighbor hoppings between the sites. Here, the position operators are
unambiguous up to a choice of origin, and so one can straightforwardly
describe the problem in length gauge. We write the Hamiltonian as,

$\hat{\mathcal{H}}=-\begin{pmatrix}0&t_{x}&0&t_{y}\\\ t_{x}&0&t_{y}&0\\\
0&t_{y}&0&t_{x}\\\ t_{y}&0&t_{x}&0\end{pmatrix}\text{.}$ (47)

Assume $t_{x}$ and $t_{y}$ are both positive, and without loss of generality,
assume $t_{x}>t_{y}$. We use a basis in which the position operators
$\hat{r}^{x}$ and $\hat{r}^{y}$ take the form,

$\hat{r}^{x}=R_{x}\begin{pmatrix}-1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&-1\end{pmatrix},\,\,\,\hat{r}^{y}=R_{y}\begin{pmatrix}1&0&0&0\\\
0&1&0&0\\\ 0&0&-1&0\\\ 0&0&0&-1\end{pmatrix}\text{,}$ (48)

where $R_{x}$ and $R_{y}$ denote half the length of the molecule in the $x$
and $y$ directions, respectively. Note that the quadrupole operators
$\hat{r}^{x}\hat{r}^{x}$ and $\hat{r}^{y}\hat{r}^{y}$ are proportional to the
identity matrix, so they cannot mediate any nonzero responses in the density
matrix equation of motion. On the other hand,
$\hat{r}^{x}\hat{r}^{y}\equiv\hat{r}^{y}\hat{r}^{x}=R_{x}R_{y}\text{diag}(-1,1,-1,1)$,
which can mediate a nonzero response. So one should only expect a response
when the electric field is perpendicular to the incident $\mathbf{q}$ vector,
which is sensible.

This Hamiltonian is diagonalized by a unitary transformation such as,

$\hat{\mathcal{U}}=\frac{1}{2}\begin{pmatrix}1&1&1&1\\\ 1&1&-1&-1\\\
1&-1&-1&1\\\ 1&-1&1&-1\end{pmatrix}\text{,}$ (49)

yielding levels with energies $-t_{x}-t_{y}$, $-t_{x}+t_{y}$, $t_{x}-t_{y}$,
$t_{x}+t_{y}$, in ascending order. We calculate the electric quadrupole SHG
response using,

$\begin{split}\sigma^{\nu\mu\alpha\beta}_{\text{SHG},(1)}(\omega)&=\frac{ie^{3}}{2\hslash^{2}}\bigg{(}\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(2)}_{\hat{Q}^{\nu\alpha},\hat{r}^{\beta}}(2\omega)\right\\}\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\text{tr}\left\\{\hat{v}^{\mu}\cdot\hat{\rho}^{(2)}_{\hat{r}^{\alpha},\hat{Q}^{\nu\beta}}(2\omega)\right\\}\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\text{tr}\left\\{\partial_{t}\hat{Q}^{\nu\mu}\cdot\hat{\rho}^{(2)}_{\hat{r}^{\alpha},\hat{r}^{\beta}}(\omega)\right\\}\bigg{)}\text{,}\end{split}$
(50)

starting from an initial density matrix of
$\hat{\rho}^{(0)}=\text{diag}(f_{1},f_{2},f_{3},f_{4})$. We find,

$\sigma^{xxyy}_{\text{SHG},(1)}(\omega)=-\frac{4ie^{3}}{\hslash^{2}}\frac{(f_{12}-f_{34})t_{x}t_{y}\omega
R_{x}^{2}R_{y}^{2}}{(\omega^{2}-t_{x}^{2})(\omega^{2}-(2t_{y})^{2})}\text{.}$
(51)
$\sigma^{yyxx}_{\text{SHG},(1)}(\omega)=-\frac{4ie^{3}}{\hslash^{2}}\frac{(f_{12}-f_{34})t_{x}t_{y}\omega
R_{x}^{2}R_{y}^{2}}{(\omega^{2}-(2t_{x})^{2})(\omega^{2}-t_{y}^{2})}\text{.}$
(52)

This result reveals a number of interesting features. First of all, it is an
example of a nonzero second harmonic response in an inversion-symmetric
system, where the inversion breaking is allowed by the nonzero order in
$\mathbf{q}$. Secondly, the response intriguingly vanishes for half-filling,
owing to a symmetry between the Bloch states in the occupied and unoccupied
energy levels. We note that introducing diagonal hoppings breaks this
symmetry, and therefore allows a nonzero response at half-filling, at the
expense of the energetic factors in Eq. (51)-(52) becoming considerably more
complicated. In other words, the response takes the form,

$\sigma^{xxyy}_{\text{SHG},(1)}(\omega)=-\frac{ie^{3}}{\hslash^{2}}F(\omega)R_{x}^{2}R_{y}^{2}\text{,}$
(53)

where $F(\omega)$ only depends on $\omega$, the hopping amplitudes, and the
Fermi factors.

Appealingly, $\sigma^{\nu\mu\alpha\beta}_{\text{SHG},(1)}(\omega)$ can be seen
to depend on four factors of the length of a molecule, in the directions
specified by $\mathbf{q}$, $\mathbf{j}$, and the two $\mathbf{E}$-fields.
Therefore, a larger molecule exhibits a larger magnitude for this spatially
dispersive correction. The size of the molecule can be considered as a proxy
for the coherence length of electronic states. Crucially, in the case of
obstructed band systems forcing delocalized Wannier states, this
characteristic coherence length may be on the order of many unit cells, and
therefore one should expect substantial spatially dispersive corrections from
these systems.

###  IV.2 Chain of rectangular molecules

Transitioning to a band theory, we study a quasi-1D model in which two copies
of the Rice-Mele model form the legs of a ladder. We define the Hamiltonian
as,

$\hat{\mathcal{H}}^{\text{o}}_{k}=-\begin{pmatrix}\Delta&g_{\text{o}}(k)&0&t_{y}\\\
g_{\text{o}}^{*}(k)&-\Delta&t_{y}&0\\\ 0&t_{y}&-\Delta&g_{\text{o}}^{*}(k)\\\
t_{y}&0&g_{\text{o}}(k)&\Delta\end{pmatrix}\text{,}$ (54)
$g_{\text{o}}(k)=t_{x}+t_{x}^{\prime}e^{-ik}\text{.}$ (55)

This Hamiltonian reduces to the molecule Hamiltonian in Eq. (47) in the limit
$t_{x}^{\prime}\to 0$ (and $\Delta\to 0$). The purpose of the parameter
$\Delta$ (which breaks inversion if $t_{x}\neq t_{x}^{\prime}$) is to allow
calculation of the spatially dispersive correction to the linear conductivity
(which vanishes under inversion symmetry). Here, we assume translational
symmetry in the $x$-direction, with crystal momentum $k\in[0,2\pi)$.

We have written the Hamiltonian in what we will call “orbital form” (hence the
subscript “o”), in which it has the desirable property of being periodic with
respect to $k$-space translations by a reciprocal lattice vector,

$\hat{\mathcal{H}}^{\text{o}}_{k}=\hat{\mathcal{H}}^{\text{o}}_{k+2\pi}\text{.}$
(56)

However, this comes at the cost of neglecting Bloch phases $e^{ik\Delta x}$
for intracell hoppings, and hence the normal Peierls substitution $k\to
k+eA(t)$ with this Hamiltonian will neglect intracell currents. We will call
this form of the Hamiltonian “orbital form,” (hence the superscript “o”) since
it may be viewed as treating all states within the unit cell as orbitals on a
single generalized site in the unit cell. There are two equivalent options for
dealing with this. One is to define the position operator with an intercell
and intracell part,

$\hat{r}_{i}\to\hat{r}_{i}+\hat{\tau}\text{,}$ (57)

and to explicitly add the intracell part of the velocity operator
$\hat{v}_{\text{intra}}\equiv i\left[\hat{\mathcal{H}},\hat{\tau}\right]$ when
calculating currents. Alternatively, one can define a form of the $k$-space
Hamiltonian such that the Bloch phases do accurately reflect intracell
processes,

$\hat{\mathcal{H}}^{\text{s}}_{k}=-\begin{pmatrix}\Delta&g_{\text{s}}(k)&0&t_{y}\\\
g_{\text{s}}^{*}(k)&-\Delta&t_{y}&0\\\ 0&t_{y}&-\Delta&g_{\text{s}}^{*}(k)\\\
t_{y}&0&g_{\text{s}}(k)&\Delta\end{pmatrix}\text{,}$ (58)
$g_{\text{s}}(k)=t_{x}e^{ik/2}+t_{x}^{\prime}e^{-ik/2}\text{,}$ (59)

and find the unitary transformation
$\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}$ that takes this Hamiltonian back
to the orbital form,

$\hat{\mathcal{H}}^{\text{o}}_{k}=\hat{\mathcal{U}}^{\text{s}\to\text{o}\dagger}_{k}\cdot\hat{\mathcal{H}}^{\text{s}}_{k}\cdot\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}\text{.}$
(60)

We call this alternative form of the Hamiltonian “sublattice form,” since it
treats states within the unit cell as being located on different sublattice
sites. We now show that implementing the Peierls substitution with
$\hat{\mathcal{H}}^{\text{o}}_{k}$ before the unitary transformation is
equivalent to taking into account both the intercell and intracell currents
(we set the electric charge ${\rm q}=-1$ for simplicity),

$\begin{split}\hat{j}^{x}&=\frac{\partial}{\partial
A}\left(\hat{\mathcal{U}}^{\text{s}\to\text{o}\dagger}_{k}\hat{\mathcal{H}}^{\text{s}}_{k+A}\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}\right)\\\
&=\frac{\partial}{\partial
k}\left(\hat{\mathcal{U}}^{\text{s}\to\text{o}\dagger}_{k}\mathcal{H}^{\text{s}}_{k}\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}\right)\\\
&\,\,\,\,\,\,\,-\frac{\partial}{\partial
k}\hat{\mathcal{U}}^{\text{s}\to\text{o}\dagger}_{k}\hat{\mathcal{H}}^{\text{s}}_{k}\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}-\hat{\mathcal{U}}^{\text{s}\to\text{o}\dagger}_{k}\hat{\mathcal{H}}^{\text{s}}_{k}\frac{\partial}{\partial
k}\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}\\\
&=\frac{\partial\hat{\mathcal{H}}^{\text{o}}_{k}}{\partial
k}-\frac{\partial}{\partial
k}\hat{\mathcal{U}}^{\text{s}\to\text{o}\dagger}_{k}\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}\hat{\mathcal{H}}^{\text{o}}_{k}-\hat{\mathcal{H}}^{\text{o}}_{k}\hat{\mathcal{U}}^{\text{s}\to\text{o}\dagger}_{k}\frac{\partial}{\partial
k}\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}\\\
&=\frac{\partial\hat{\mathcal{H}}^{\text{o}}_{k}}{\partial
k}+i\left[\hat{\mathcal{H}}^{\text{o}}_{k}\,\,,\,\,i\,\hat{\mathcal{U}}^{\text{s}\to\text{o}\dagger}_{k}\frac{\partial}{\partial
k}\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}\right]\\\
&\equiv\hat{v}^{x}_{\text{inter}}+\hat{v}^{x}_{\text{intra}}\text{.}\end{split}$

We therefore recognize the identity
$i\,\hat{\mathcal{U}}^{\text{s}\to\text{o}\dagger}_{k}\frac{\partial}{\partial
k}\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}\equiv\hat{\tau}$, or in other
words, $\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}\equiv e^{-ik\hat{\tau}}$.
In the case at hand, we have, (up to an overall phase representing a shift in
origin)

$\hat{\mathcal{U}}^{\text{s}\to\text{o}}_{k}=\begin{pmatrix}e^{ik/4}&0&0&0\\\
0&e^{-ik/4}&0&0\\\ 0&0&e^{-ik/4}&0\\\ 0&0&0&e^{ik/4}\end{pmatrix}$ (61)
$\implies\hat{\tau}=\begin{pmatrix}-\frac{1}{4}&0&0&0\\\ 0&\frac{1}{4}&0&0\\\
0&0&\frac{1}{4}&0\\\ 0&0&0&-\frac{1}{4}\end{pmatrix}\text{.}$ (62)

That is, the position operator acquires an intracell part encoding the fact
that two of the sublattices are displaced in the $x$-direction by half the
width of the unit cell. Therefore, in what follows, we use the following form
for the Peierls-substituted Hamiltonian,

$\hat{\mathcal{H}}^{\text{s}}_{k}=-\begin{pmatrix}\Delta&g_{A}(k)&0&t_{y}\\\
g_{A}^{*}(k)&-\Delta&t_{y}&0\\\ 0&t_{y}&-\Delta&g_{A}^{*}(k)\\\
t_{y}&0&g_{A}(k)&\Delta\end{pmatrix}\text{,}$ (63)
$g_{A}(k)=e^{-ik/2}\left[(t+\delta)e^{i(k-A)/2}+(t-\delta)e^{-i(k-A)/2}\right]\text{,}$
(64)

which has both the periodicity $k\to k+2\pi$ and a dependence on $A$ that
encodes both intracell and intercell currents.

We can then define the current operators in the $x$-direction in terms of
derivatives with respect to $A$,

$\hat{v}^{x}\equiv\hat{h}^{x}=\left.\hat{\mathcal{U}}^{\dagger}_{k}\cdot\frac{\partial\hat{\mathcal{H}}^{A}_{k}}{\partial
A}\cdot\hat{\mathcal{U}}_{k}\right|_{A=0}\text{,}$ (65)
$\hat{h}^{xx}=\left.\hat{\mathcal{U}}^{\dagger}_{k}\cdot\frac{\partial^{2}\hat{\mathcal{H}}^{A}_{k}}{\partial
A^{2}}\cdot\hat{\mathcal{U}}_{k}\right|_{A=0}\text{,}$ (66)
$\hat{h}^{xxx}=\left.\hat{\mathcal{U}}^{\dagger}_{k}\cdot\frac{\partial^{3}\hat{\mathcal{H}}^{A}_{k}}{\partial
A^{3}}\cdot\hat{\mathcal{U}}_{k}\right|_{A=0}\text{,}$ (67)

where $\hat{\mathcal{U}}_{k}$ is some (non-unique) unitary transformation that
diagonalizes $\hat{\mathcal{H}}^{A=0}_{k}$. As for the current operators in
the $y$-direction, since there is no periodicity in $y$, we use,

$\hat{v}^{y}=i\,\hat{\mathcal{U}}^{\dagger}_{k}\cdot\left[\hat{\mathcal{H}}^{\text{o}}_{k}\,,\,\hat{r}^{y}\right]\cdot\hat{\mathcal{U}}_{k}\text{,}$
(68)

with $\hat{r}^{y}$ defined just as in the molecule case (Eq. (48)). For
simplicity, we set the width of the chain in the $y$-direction to be $1/2$, to
match up with the half-unit-cell width in the $x$-direction of a single
decoupled molecule. (That is, in comparing with Eq. (48), we will set
$R_{x}=R_{y}=1/4$.) One can then formally compute higher-order currents in the
$y$-direction using the covariant derivative, where in this case there is no
partial derivative, and so the covariant derivative reduces to a commutator
with the Berry connection $\hat{\mathcal{A}}^{y}\equiv\hat{r}^{y}$,

$\hat{h}^{yy}=-i\left[\hat{\mathcal{A}}^{y}\,,\,\hat{v}^{y}\right]\text{,}$
(69)

and so on for all other current operators
$\hat{h}^{\alpha_{1}\ldots\alpha_{n}}$.

###  IV.3 Numerical results

Figures 7 and 8 show numerical results for the spatially dispersive
corrections to the linear conductivity and SHG conductivity, respectively. We
benchmark our methodology in two ways. First, we compare the conductivities of
the chain in the $t_{x}^{\prime}\to 0$ limit to the conductivities calculated
for a single rectangular molecule. Second, for the linear conductivity, we
also benchmark against established results Malashevich and Souza (2010); Ahn
_et al._ (2022b) for the spatially dispersive linear conductivity calculated
by expanding the Kubo formula to first order in $\mathbf{q}$. Figures 7 and 8
show agreement in all cases. In our plots, we introduce a phenomenological
relaxation parameter $\eta$ by taking $\omega\to\omega+i\eta$, with
$\eta=0.01$. We set $t_{x}=1$, without loss of generality, and in the case for
the linear conductivity, we set $\Delta=0.1$ in both the chain and the
molecule to statically break inversion (and to therefore allow a nonzero
spatially dispersive correction). We also always assume the lowest three bands
in this model are filled,

$\hat{\rho}^{(0)}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&0\end{pmatrix}\text{,}$ (70)

because the spatially dispersive correction to SHG depends on a peculiar
combination of the filling factors (e.g. for the molecule, it is proportional
to $f_{1}-f_{2}-f_{3}+f_{4}$).

As a note on symmetries, the molecular chain we consider is inversion
symmetric when either $\Delta=0$ or $\delta=0$, where we define
$\delta=t_{x}-t_{x}^{\prime}$. This is only evident from the “sublattice form”
$\hat{\mathcal{H}}_{k}^{\text{s}}$ of the Hamiltonian in Eq. (58). When
$\delta=0$, the system has an inversion center on a vertical bond, represented
by the following matrix transformation,

$\hat{\mathcal{H}}_{k}^{\text{s}}\longrightarrow\hat{\mathcal{H}}_{-k}^{\text{s}}\text{.}$
(71)

Additionally, the system has an inversion center at the midpoint of the unit
cell whenever $\Delta=0$, represented by,

$\hat{\mathcal{H}}_{k}^{\text{s}}\longrightarrow\begin{pmatrix}0&0&1&0\\\
0&0&0&1\\\ 1&0&0&0\\\
0&1&0&0\end{pmatrix}\cdot\hat{\mathcal{H}}_{-k}^{\text{s}}\cdot\begin{pmatrix}0&0&1&0\\\
0&0&0&1\\\ 1&0&0&0\\\ 0&1&0&0\end{pmatrix}\text{.}$ (72)

So inversion is only broken when both $\delta\neq 0$ and $\Delta\neq 0$.

Figure 9: Magnitude of the total SHG response computed up to electric
quadrupole order resulting from off-axis linearly polarized light incident on
the molecular chain model, evaluated at
$\omega_{0}=\sqrt{(t_{x}-t_{x}^{\prime})^{2}+\Delta^{2}}$ (half the minimum
gap between the filled and empty bands). Light polarized in $x$ (chain
direction) leads to a conventional SHG response when inversion is broken
($\Delta\neq 0$) while light polarized in $y$ leads to a spatially dispersive
correction ($q\neq 0$). Here, $t_{x}-t_{x}^{\prime}=0.5$ and $t_{y}=0.8$.

Given these considerations, we can directly compare SHG resulting from static
inversion breaking with SHG resulting from a spatially dispersive correction.
Figure 9 plots the magnitude of the total SHG response resulting from
diagonally polarized light at a frequency of half the minimum gap in this
model. We choose diagonal polarization because $x$-polarized light triggers a
uniform SHG response ($\sigma_{\text{SHG},(0)}^{xxx}$) when $\Delta$ is
nonzero, while $y$-polarized light triggers a spatially dispersive correction
($\sigma_{\text{SHG},(1)}^{xxyy}$) at first order in $q_{x}$. In that plot, we
choose a window such that the contours are roughly circular in shape, and find
a pair of dimensionless parameters in which the two axes have roughly the same
scale. We find these relevant dimensionless parameters to be
$\Delta/t_{x}^{\prime}$ and $q_{x}R_{x}$. That is, the effects of spatial
dispersion on the length scale of half the width of an individual molecule
create a similar magnitude of SHG to that of a static inversion breaking
parameter on the energy scale of the inter-cell hopping. This heurstic has
intriguing implications for real systems where spatially dispersive effects
may be present. For instance, moiré systems famously create extremely inflated
artificial lattice constants, in which case $qR$ is very large (where $R$ is
interpreted generally as some coherence length for localized electronic
states). So long as the hopping amplitude between superlattice sites is not
too small, the spatially dispersive correction to SHG may be on the same order
as a typical static inversion-breaking mechanisms for which SHG is often used
as a probe. This has important ramifications for using SHG to detect
spontaneous breaking of inversion symmetry, since many systems of recent
interest may have sufficiently large coherence length scales $R$ to allow
spatially dispersive effects to dominate such a response.

##  V Outlook

We have introduced a new scheme for calculating spatially dispersive
corrections to nonlinear optical responses in velocity gauge. Reviewing from a
diagrammatic standpoint the subtle cancellations that make up the sum rules
guaranteeing equivalence between the length and velocity gauges, we compactly
write velocity gauge expressions for nonlinear optical conductivities and
introduce spatially dispersive corrections as an additional vertex rule.
Though we only treat the electric quadrupole contribution here, we expect the
formalism can be straightforwardly extended to include the magnetic dipole
contribution, as well as higher-multipole corrections as desired. We expect
spatially dispersive corrections to be important for understanding anomalously
strong bulk second-order responses in centrosymmetric materials where the
response is typically expected to vanish.

Our results put forward a new tool for using nonlinear optics to study quantum
materials beyond electric dipole order. We expect spatially dispersive
corrections to optical responses to reveal information about intrinsic
nonlocality in material systems, since they should be strong when the electric
field varies appreciably on the length scale of the electronic states being
coupled by optical transitions. A number of materials of recent interest fall
into this category, including Moiré materials Koshino  _et al._ (2018); Kang
and Vafek (2018) and systems whose bands exhibit nontrivial quantum geometry
Brouder  _et al._ (2007); Marzari  _et al._ (2012); Po  _et al._ (2018) for
which maximally localized Wannier representations still feature coherence over
length scales of many unit cells. Searching for measurable spatially
dispersive corrections in these systems and linking them to quantum geometry
is a promising direction for future work.

###### Acknowledgements.

The authors thank M. Claassen, C. De Beule, Z. Addison, and L. Wu for
insightful discusssions. S.G. was supported by an NSF Graduate Research
Fellowship. This work was supported by the Department of Energy under Grant
DE-FG02-84ER45118.

## References

  * Boyd (2003) R. W. Boyd,  _Nonlinear optics_ (Academic Press, London, 2003).
  * Sekino and Bartlett (1986) Hideo Sekino and Rodney J. Bartlett, “Frequency dependent nonlinear optical properties of molecules,” The Journal of Chemical Physics 85, 976–989 (1986).
  * Young  _et al._ (2012) Steve M. Young, Fan Zheng, and Andrew M. Rappe, “First-principles calculation of the bulk photovoltaic effect in bismuth ferrite,” Phys. Rev. Lett. 109, 236601 (2012).
  * Morimoto and Nagaosa (2016) Takahiro Morimoto and Naoto Nagaosa, “Topological nature of nonlinear optical effects in solids,” Science Advances 2, e1501524 (2016).
  * Wu  _et al._ (2017) Liang Wu, S. Patankar, T. Morimoto, N. L. Nair, E. Thewalt, A. Little, J. G. Analytis, J. E. Moore, and J. Orenstein, “Giant anisotropic nonlinear optical response in transition metal monopnictide weyl semimetals,” Nature Physics 13, 350–355 (2017).
  * de Juan  _et al._ (2017) Fernando de Juan, Adolfo G. Grushin, Takahiro Morimoto, and Joel E. Moore, “Quantized circular photogalvanic effect in weyl semimetals,” Nature Communications 8, 15995 (2017).
  * Ma  _et al._ (2017) Qiong Ma, Su-Yang Xu, Ching-Kit Chan, Cheng-Long Zhang, Guoqing Chang, Yuxuan Lin, Weiwei Xie, Tomás Palacios, Hsin Lin, Shuang Jia, Patrick A. Lee, Pablo Jarillo-Herrero, and Nuh Gedik, “Direct optical detection of weyl fermion chirality in a topological semimetal,” Nature Physics 13, 842–847 (2017).
  * Chan  _et al._ (2017) Ching-Kit Chan, Netanel H. Lindner, Gil Refael, and Patrick A. Lee, “Photocurrents in weyl semimetals,” Phys. Rev. B 95, 041104 (2017).
  * Holder  _et al._ (2020) Tobias Holder, Daniel Kaplan, and Binghai Yan, “Consequences of time-reversal-symmetry breaking in the light-matter interaction: Berry curvature, quantum metric, and diabatic motion,” Phys. Rev. Res. 2, 033100 (2020).
  * Ni  _et al._ (2022) Zhuoliang Ni, Nan Huang, Amanda V. Haglund, David G. Mandrus, and Liang Wu, “Observation of giant surface second-harmonic generation coupled to nematic orders in the van der waals antiferromagnet feps3,” Nano Letters 22, 3283–3288 (2022).
  * Ahn  _et al._ (2022a) Junyeong Ahn, Guang-Yu Guo, Naoto Nagaosa, and Ashvin Vishwanath, “Riemannian geometry of resonant optical responses,” Nature Physics 18, 290–295 (2022a).
  * Tai and Claassen (2023) Wai Ting Tai and Martin Claassen, “Quantum-geometric light-matter coupling in correlated quantum materials,” (2023), arXiv:2303.01597 [cond-mat.str-el] .
  * Kobe and Smirl (1978) Donald H. Kobe and Arthur L. Smirl, “Gauge invariant formulation of the interaction of electromagnetic radiation and matter,” American Journal of Physics 46, 624–633 (1978).
  * Sipe and Ghahramani (1993) J. E. Sipe and Ed Ghahramani, “Nonlinear optical response of semiconductors in the independent-particle approximation,” Phys. Rev. B 48, 11705–11722 (1993).
  * Aversa and Sipe (1995) Claudio Aversa and J. E. Sipe, “Nonlinear optical susceptibilities of semiconductors: Results with a length-gauge analysis,” Phys. Rev. B 52, 14636–14645 (1995).
  * Ventura  _et al._ (2017) G. B. Ventura, D. J. Passos, J. M. B. Lopes dos Santos, J. M. Viana Parente Lopes, and N. M. R. Peres, “Gauge covariances and nonlinear optical responses,” Phys. Rev. B 96, 035431 (2017).
  * Passos  _et al._ (2018) D. J. Passos, G. B. Ventura, J. M. Viana Parente Lopes, J. M. B. Lopes dos Santos, and N. M. R. Peres, “Nonlinear optical responses of crystalline systems: Results from a velocity gauge analysis,” Phys. Rev. B 97, 235446 (2018).
  * Schüler  _et al._ (2021) Michael Schüler, Jacob A. Marks, Yuta Murakami, Chunjing Jia, and Thomas P. Devereaux, “Gauge invariance of light-matter interactions in first-principle tight-binding models,” Phys. Rev. B 103, 155409 (2021).
  * Parker  _et al._ (2019) Daniel E. Parker, Takahiro Morimoto, Joseph Orenstein, and Joel E. Moore, “Diagrammatic approach to nonlinear optical response with application to weyl semimetals,” Phys. Rev. B 99, 045121 (2019).
  * Blount (1962) E. I. Blount, “Bloch electrons in a magnetic field,” Phys. Rev. 126, 1636–1653 (1962).
  * Malashevich and Souza (2010) Andrei Malashevich and Ivo Souza, “Band theory of spatial dispersion in magnetoelectrics,” Phys. Rev. B 82, 245118 (2010).
  * Ahn  _et al._ (2022b) Junyeong Ahn, Su-Yang Xu, and Ashvin Vishwanath, “Theory of optical axion electrodynamics and application to the kerr effect in topological antiferromagnets,” Nature Communications 13, 7615 (2022b).
  * Melrose and McPhedran (2009) D B Melrose and R C McPhedran,  _Electromagnetic Processes in Dispersive Media_ (Cambridge University Press, Cambridge, England, 2009).
  * Cohen-Tannoudji  _et al._ (1997) Claude Cohen-Tannoudji, Jacques Dupont-Roc, and Gilbert Grynberg,  _Photons and atoms_ , edited by Claude Cohen-Tannoudji, Jacques Dupont-Roc, and Gilbert Grynberg (John Wiley & Sons, Nashville, TN, 1997).
  * Jaszuński  _et al._ (2017) Michał Jaszuński, Antonio Rizzo, and Kenneth Ruud, “Molecular electric, magnetic, and optical properties,” in  _Handbook of Computational Chemistry_, edited by Jerzy Leszczynski, Anna Kaczmarek-Kedziera, Tomasz Puzyn, Manthos G. Papadopoulos, Heribert Reis, and Manoj K. Shukla (Springer International Publishing, Cham, 2017) pp. 497–592.
  * Koshino  _et al._ (2018) Mikito Koshino, Noah F. Q. Yuan, Takashi Koretsune, Masayuki Ochi, Kazuhiko Kuroki, and Liang Fu, “Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene,” Phys. Rev. X 8, 031087 (2018).
  * Kang and Vafek (2018) Jian Kang and Oskar Vafek, “Symmetry, maximally localized wannier states, and a low-energy model for twisted bilayer graphene narrow bands,” Phys. Rev. X 8, 031088 (2018).
  * Brouder  _et al._ (2007) Christian Brouder, Gianluca Panati, Matteo Calandra, Christophe Mourougane, and Nicola Marzari, “Exponential localization of wannier functions in insulators,” Phys. Rev. Lett. 98, 046402 (2007).
  * Marzari  _et al._ (2012) Nicola Marzari, Arash A. Mostofi, Jonathan R. Yates, Ivo Souza, and David Vanderbilt, “Maximally localized wannier functions: Theory and applications,” Rev. Mod. Phys. 84, 1419–1475 (2012).
  * Po  _et al._ (2018) Hoi Chun Po, Haruki Watanabe, and Ashvin Vishwanath, “Fragile topology and wannier obstructions,” Phys. Rev. Lett. 121, 126402 (2018).