# Resonance-enhanced optical nonlinearity in the Weyl semimetal TaAs

Shreyas Patankar  Department of Physics, University of California, Berkeley,
California 94720, USA  Materials Science Division, Lawrence Berkeley National
Laboratory, Berkeley, California 94720, USA  Liang Wu  Department of Physics,
University of California, Berkeley, California 94720, USA  Materials Science
Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720,
USA  Department of Physics & Astronomy, University of Pennsylvania,
Philadelphia, PA 19104  Baozhu Lu  Department of Physics, Temple University,
Philadelphia, Pennsylvania 19122, USA  Temple Materials Institute, Temple
University, Philadelphia, PA 19122, USA  Manita Rai  Department of Physics,
Temple University, Philadelphia, Pennsylvania 19122, USA  Temple Materials
Institute, Temple University, Philadelphia, PA 19122, USA  Jason D. Tran
Department of Physics, Temple University, Philadelphia, Pennsylvania 19122,
USA  Temple Materials Institute, Temple University, Philadelphia, PA 19122,
USA  T. Morimoto  Department of Physics, University of California, Berkeley,
California 94720, USA  D. Parker  Department of Physics, University of
California, Berkeley, California 94720, USA  Adolfo G. Grushin  Department of
Physics, University of California, Berkeley, California 94720, USA  Institut
Néel, CNRS and Université Grenoble Alpes, F-38042 Grenoble, France  N. L. Nair
Department of Physics, University of California, Berkeley, California 94720,
USA  J. G. Analytis  Department of Physics, University of California,
Berkeley, California 94720, USA  Materials Science Division, Lawrence Berkeley
National Laboratory, Berkeley, California 94720, USA  J. E. Moore  Department
of Physics, University of California, Berkeley, California 94720, USA
Materials Science Division, Lawrence Berkeley National Laboratory, Berkeley,
California 94720, USA  J. Orenstein<EMAIL_ADDRESS>Department of
Physics, University of California, Berkeley, California 94720, USA  Materials
Science Division, Lawrence Berkeley National Laboratory, Berkeley, California
94720, USA  Darius H. Torchinsky  Department of Physics, Temple University,
Philadelphia, Pennsylvania 19122, USA  Temple Materials Institute, Temple
University, Philadelphia, PA 19122, USA

While all media can exhibit first-order conductivity describing current
linearly proportional to electric field, $\displaystyle E$, the second-order
conductivity, $\displaystyle\sigma^{(2)}$, relating current to $\displaystyle
E^{2}$, is nonzero only when inversion symmetry is broken. Second-order
nonlinear optical responses are powerful tools in basic research, as probes of
symmetry breaking shen1986surface ; shen1989surface ; dahn1996symmetry ;
petersen2006nonlinear ; harter2017parity , and in optical technology as the
basis for generating currents from far-infrared robinsonFIR  to X-ray
wavelengths chang1997generation . The recent surge of interest in Weyl
semimetals with acentric crystal structures has led to the discovery of a host
of $\displaystyle\sigma^{(2)}$-related phenomena in this class of materials,
such as polarization-selective conversion of light to dc current
(photogalvanic effects) ma2017direct ; sun2017circular ;
osterhoudt2017colossal ; ji2018spatially  and the observation of giant second-
harmonic generation (SHG) efficiency in TaAs at photon energy 1.5 eV
wu2017giant ; li2018 . Here, we present measurements of the SHG spectrum of
TaAs revealing that the response at 1.5 eV corresponds to the high-energy tail
of a resonance at 0.7 eV, at which point the second harmonic conductivity is
approximately 200 times larger than seen in the standard candle nonlinear
crystal, GaAs. This remarkably large SHG response provokes the question of
ultimate limits on $\displaystyle\sigma^{(2)}$, which we address by a new
theorem relating frequency-integrated nonlinear response functions to the
third cumulant (or “skewness”) of the polarization distribution function in
the ground state. This theorem provides considerable insight into the factors
that lead to the largest possible second-order nonlinear response,
specifically showing that the spectral weight is unbounded and potentially
divergent when the possibility of next-neighbor hopping is included.

Second-harmonic generation (SHG) is a subset of
$\displaystyle\sigma^{(2)}$-related phenomena in which light of frequency
$\displaystyle 2\omega$ is generated by excitation at frequency
$\displaystyle\omega$. Figure 1 shows a schematic of the apparatus used to
measure the spectrum of the SHG response tensor, which is formally written as
$\displaystyle\sigma_{ijk}(2\omega;\omega,\omega)$ and hereafter shortened to
$\displaystyle\sigma^{\text{shg}}_{ijk}(\omega)$, or simply
$\displaystyle\sigma_{ijk}$. The measurements are enabled by a laser/optical
parametric amplifier system that generates pulses of duration
$\displaystyle\sim$50 fs at a 5 kHz repetition rate in the range of photon
energies 0.5 - 1.5 eV. As the crystals are opaque (optical penetration depth
of approximately 200 nm xuEllipsometry ) the SHG intensity is measured in
reflection. As shown in the Figure 1a., a polarizer $\displaystyle P$, and a
halfwave plate in front of the sample select the orientation of the electric
field (referred to by $\displaystyle\hat{\mathbf{{e}}}$) of the fundamental
light at frequency $\displaystyle\omega$, and a second analyzing polarizer
$\displaystyle A$ selects the polarization of the light at frequency
$\displaystyle 2\omega$ that reaches the detector.

Figure 1: Schematic of the experimental apparatus used to study the spectrum
of the nonlinear optical response in TaAs. Laser pulses of pulse width 50 fs,
at a repetition rate of 5 kHz with photon energy between 0.5 eV - 1.6 eV are
generated by an optical parametric amplifier coupled to an regeneratively
amplified laser. The incident radiation passes through a neutral density
filter (ND), polarizer ($\displaystyle P$), and a low pass filter (LF), to
lower the incident intensity, and to remove spurious polarizations, and
wavelengths, respectively. A halfwave plate (HW) is used to rotate the
incident linear polarization over 360o. The incident laser radiation is then
focused on to the TaAs sample (S) using a reflective objective (RO), which is
insensitive to the photon energy, in contrast with conventional refractive
microscope objectives. The SHG radiation from the sample is collimated by the
RO, and picked off by a D-mirror (D), and passed through an analyzing
polarizer ($\displaystyle A$), whose polarization state is determined by the
scan being measured (see text). The radiation then passes through a stack of
short pass filters (SF) to remove any incident light, and then focused by a
lens (L) onto a photomultiplier tube (PMT) for detection. Inset: A single
crystal of TaAs with the $\displaystyle(112)$ facet visible.

The structure of TaAs (point group $\displaystyle 4mm$) contains two
perpendicular glide planes, which are equivalent to mirror planes for optical
response functions TaAsWeng ; li2016weyl . The normal directions of these two
planes determine the two equivalent directions of the tetragonal unit cell,
which we label as $\displaystyle x$ and $\displaystyle y$. The direction
perpendicular to the $\displaystyle xy$ plane is the polar, or $\displaystyle
z$ axis. The $\displaystyle 4mm$ point group allows three independent
nonvanishing components of the conductivity tensor,
$\displaystyle\sigma_{zzz}$, $\displaystyle\sigma_{zxx}=\sigma_{zyy}$, and
$\displaystyle\sigma_{xzx}=\sigma_{yzy}$. To avoid using a surface obtained by
cutting and subsequent polishing, which are known to affect the nonlinear
optical response, we perform nonlinear reflectance measurements using the
(112) surface, as it forms a naturally flat and smooth growth facet. Figure 1b
shows a photograph of the (112) surface used for the measurement. The two high
symmetry directions in the plane of this surface are [1,-1,0], which is
perpendicular to the polar axis, and [1,1,-1].

We performed two pairs of polarization scans at each wavelength. In the first
pair, $\displaystyle\hat{\mathbf{{e}}}$ and $\displaystyle A$ are co-rotated
with their relative angle fixed at 0 and 90∘. This is equivalent to rotating
the sample, but with the advantage that the position of the laser focus on the
sample surface remains constant. These scans are primarily used to determine
the high symmetry directions. However, as will be discussed below, they
provide additional evidence for the existence of a sharp resonance in the
nonlinear response function in the infrared. In the second pair of scans
$\displaystyle A$ is fixed parallel to either of the [1,-1,0] or [1,1,-1] high
symmetry directions and $\displaystyle\hat{\mathbf{{e}}}$ is rotated over
360∘. From this pair of scans we obtain the three combinations of tensor
components that are available from reflection measurements performed on the
(112) surface: $\displaystyle|\sigma_{zxx}|$, $\displaystyle|\sigma_{xzx}|$,
and
$\displaystyle|\sigma_{\mathrm{eff}}|\equiv|\sigma_{zzz}+2\sigma_{zxx}+4\sigma_{xzx}|$
seeApp . Analysis of the polarization scans at a given wavelength determines
the relative magnitudes of these tensor components, i.e.,
$\displaystyle|\sigma_{zxx}|/|\sigma_{\mathrm{eff}}|$ and
$\displaystyle|\sigma_{xzx}|/|\sigma_{\mathrm{eff}}|$.

Measurements of absolute, as opposed to relative, magnitudes of
$\displaystyle\sigma_{ijk}$ components over a broad range of frequency were
accomplished by using GaAs to calibrate the response at 1.5 eV. Determining
the effective $\displaystyle\sigma_{ijk}(\omega)$ at lower photon energies
required characterization of the wavelength dependence of all components of
the optical setup, such as wavelength-sorting filters, attenuators, and
photodetectors, as well as spectral variation of the laser focal spot size and
pulse duration. An extensive discussion of the calibration procedure is
provided in the Supplementary Information.

Polarization scans were performed at incident wavelengths in the interval from
800 to 2500 nm and a complete library of these data is included in the
Supplementary Information. Figure 2(a-c) shows representative scans at
wavelengths 800, 1560, and 2200 nm (1.55, 0.80, and 0.56 eV, respectively)
with co-rotating parallel polarizations. The co-rotation plot at 800 nm
illustrates the extreme anisotropy of the nonlinear optical response of TaAs.
The solid curve through the data points is proportional to
$\displaystyle|\sigma_{zzz}|^{2}\cos^{6}\theta$, which is the dependence
predicted for a crystal that responds only to $\displaystyle E$ parallel to
its polar axis and radiates second harmonic light that is likewise fully
$\displaystyle z$-polarized. The polar plots at longer wavelength indicate
that the amplitudes of the off-diagonal components,
$\displaystyle|\sigma_{zxx}|$ and $\displaystyle|\sigma_{xzx}|$, begin to
increase relative to $\displaystyle|\sigma_{\text{eff}}|$ as the fundamental
photon energy approaches 0.7 eV from above, although they remain approximately
an order of magnitude smaller, as is shown below.

Figure 2: Second harmonic generation polarimetry. SHG intensity measured in
the “parallel” scan as a function of the incident polarization angle, plotted
on a normalized scale, for three different incident photon wavelengths. The
polar pattern changes from a two-fold pattern at shorter wavelengths (800
nm=1.55 eV) to a six-fold pattern at longer wavelengths (2200 nm=0.56 eV)

Figure 3 illustrates the spectral dependence of the SHG response function of a
TaAs crystal at room temperature. Figure 3a shows
$\displaystyle|\sigma_{\mathrm{eff}}|$, $\displaystyle|\sigma_{zxx}|$, and
$\displaystyle|\sigma_{xzx}|$ vs. fundamental photon energy as solid circles
(the solid lines are guides to the eye). The $\displaystyle y$-axis scale
expresses $\displaystyle\sigma^{\text{shg}}$ in units of the conductance
quantum per Volt, to facilitate comparison with theory. In SI units, the
corresponding peak value of $\displaystyle\sigma^{\text{shg}}$ at 0.7 eV is
approximately $\displaystyle 5\times 10^{-3}\;(\Omega\text{ V})^{-1}$. The
dashed line is the spectrum of $\displaystyle\sigma_{xyz}$ in GaAs as reported
in Ref. bergfeld2003second , multiplied by a factor 10 so that it can be
compared with the TaAs spectra. As reported previously, the SH response
measured at 1.5 eV in TaAs exceeds the peak response of GaAs by a factor
$\displaystyle\sim{}10$. The measurements reported here show that
$\displaystyle|\sigma_{\mathrm{eff}}|$ becomes even larger with decreasing
frequency, reaching a peak near 0.7 eV where it is $\displaystyle\sim 2\times
10^{2}$ larger than the maximum response of GaAs.

Although the diagonal response function $\displaystyle|\sigma_{zzz}|$ is not
determined directly from the SHG intensity, its range of values is highly
constrained through analysis of the polarization scans seeApp . The shaded
region in Figure 3b illustrates the upper and lower bounds on
$\displaystyle|\sigma_{zzz}|$. The dashed line that lies within the shaded
region is a best fit to a phenomenological model described below, which is
seen to capture the basic features of the resonance enhancement of
$\displaystyle|\sigma_{zzz}|$.

The key features of the nonlinear conductivity in TaAs –– strong anisotropy
and a single, sharp resonant peak –– suggest a description in terms of a one-
dimensional (1D) two-band model, modified for a 3D system. Recently it was
shown that the Rice-Mele (RM) tight-binding model of a polar semiconductor
RM82  provides an excellent description of the spectrum of the shift current
response in single-layer monochalcogenide semiconductors rangel2017large . As
we discuss below, for such models there is precise correspondence between
shift current and SHG spectra MN16 . Moreover, the energies probed in our
experiment are much larger than either the Fermi energy xu2015discovery ;
TaAsLv  or plasma frequency xu2016optical  of TaAs, suggesting that the
semimetallic part of its bandstructure is not relevant to the nonlinear
response at photon energies of order $\displaystyle 1$ eV.

These considerations lead us to compare the measured SHG spectra with the
predictions for a system described by the RM Hamiltonian,

$H_{RM}=t\cos(k_{z}a/2)\sigma_{z}+\delta\sin(k_{z}a/2)\sigma_{y}+\Delta\sigma_{z},$
(1)

where the $\displaystyle\sigma_{i}$ are the Pauli matrices. As shown in Figure
4, $\displaystyle H_{RM}$ describes a 1D chain of atoms along the
$\displaystyle z$ direction, in which inversion symmetry is broken by
alternating on-site energies ($\displaystyle\pm\Delta$) and hopping amplitudes
($\displaystyle t\pm\delta$). Despite its relative simplicity, the RM model
has played a key role in the development of the modern theory of polarization
resta2007theory  and recently in strategies to enhance nonlinear optical
response functions tan2017upper .

The optical response of the RM Hamiltonian captures qualitative features of
the nonlinear optical response in TaAs, as
$\displaystyle\sigma^{\text{shg}}_{ijk}(\omega)$ is polarized parallel to the
polar axis and exhibits a sharp peak at the threshold for absorption. However,
to allow a quantitative comparison with experiment the RM model must be
modified to include an interchain coupling to smooth the 1D van Hove
singularity that is obtained for an isolated chain.

To extend the model to 3D, we consider an array of RM chains parallel to
$\displaystyle z$ and allow electrons to hop laterally in the $\displaystyle
A$ and $\displaystyle B$ atom layers (see Figure 4). The interchain hopping
adds an additional term $\displaystyle t_{AB}(\cos k_{x}a_{x}+\cos
k_{y}a_{y})$ to the coefficient of $\displaystyle\sigma_{z}$ in $\displaystyle
H_{RM}$, where $\displaystyle t_{AB}\equiv t_{\parallel,A}-t_{\parallel,B}$.
The shift current spectrum in this model can be calculated following
SipeShkrebtii  (see seeApp ),

$\sigma_{zzz}^{\text{shift}}(\omega)=\frac{2e^{3}}{\hbar\Delta}\left(\frac{1}{4\pi}\right)^{3}\frac{c^{2}}{ab}F(\tilde{\omega};\tilde{\delta},\tilde{t},\tilde{t}_{AB}),$
(2)

where $\displaystyle F$ is a dimensionless function of the parameters of
$\displaystyle H_{RM}$, and of frequency $\displaystyle\omega$ (the tildes
indicate normalization by $\displaystyle\Delta$); $\displaystyle a$,
$\displaystyle b$, and $\displaystyle c$ are lattice constants in the
$\displaystyle x$, $\displaystyle y$, and $\displaystyle z$ directions,
respectively seeApp . The real part of the SHG conductivity can be expressed
in terms of the shift conductivity as follows MN16 ,

$\text{Re}\left\\{\sigma^{\text{shg}}_{zzz}(\omega)\right\\}=-\frac{1}{2}\sigma_{zzz}^{\text{shift}}(\omega)+\sigma_{zzz}^{\text{shift}}(2\omega)$
(3)

Equation (3) shows that the SHG response is the sum of two terms that are
opposite in sign but have the same functional form when the frequency argument
is scaled by a factor of two. The first (“one-photon”) term corresponds to the
resonant condition $\displaystyle\hbar\omega=E_{\text{gap}}(k)$ and the second
(“two-photon”) term corresponds to
$\displaystyle\hbar\omega=E_{\text{gap}}(\textbf{k})/2$, where $\displaystyle
E_{\text{gap}}(\textbf{k})$ is the energy gap at wavevector k. Clearly the
total $\displaystyle\sigma^{\text{shg}}_{zzz}$ response integrates to zero,
but we shall show momentarily that for some models each term has a geometrical
interpretation in terms of the skewness or intrinsic asymmetry of the ground-
state polarization distribution.

Figure 3: Measured nonlinear optical conductivity as a function of incident
photon energy. (a) Spectra of the conductivity components
$\displaystyle|\sigma_{zxx}|$, $\displaystyle|\sigma_{xzx}|$, and
$\displaystyle|\sigma_{\text{eff}}|=|\sigma_{zzz}+4\sigma_{xzx}+2\sigma_{zxx}|$.
The solid lines are guides to the eye. The nonlinear optical conductivity of
GaAs, $\displaystyle|\sigma_{xyz}|$ is multiplied by 100 and plotted for
comparison. (b) The (expt) lower and upper bounds of
$\displaystyle|\sigma_{zzz}|$ are estimated using measured values of
$\displaystyle|\sigma_{\text{eff}}|$, $\displaystyle|\sigma_{zxx}|$, and
$\displaystyle|\sigma_{xzx}|$ from experiment. The dashed black line depicts
the best fit to $\displaystyle|\sigma_{zzz}|$ using the phenomenological model
described in the text.

The dash-dotted curve in Figure 3b is a fit to
$\displaystyle|\sigma^{\text{shg}}_{zzz}|$ based on the two-photon term in
equation (3). Given the simplicity of the model relative to the complexity of
the full 3D bandstructure of TaAs, the 3D version of the RM model describes
the data remarkably well. The peak position, low, and high energy tails of the
spectrum are best described using parameters $\displaystyle\tilde{t}=1.5$,
$\displaystyle\tilde{\delta}=1.4$, $\displaystyle\tilde{t}_{AB}=0.02$, and
$\displaystyle\Delta=0.428$. The lattice parameters for TaAs are
$\displaystyle c=1.165$ nm, and $\displaystyle b=a=0.344$ nm.

The theory discussed above helps to identify the factors that contribute to a
large SHG response. The first is a polar crystal structure such as
$\displaystyle 4mm$ (TaAs, NbAs, etc.) or $\displaystyle mm2$ (single layer
monochalcogenides) consistent with a phenomenological description in terms of
coupled ferroelectric chains. However, even when this description is possible
there is a global bound on the frequency integrated second-order response,
$\displaystyle\Sigma_{z}^{\text{shift}}$, of a system described by
$\displaystyle H_{RM}$ that is given by,

$\Sigma_{z}^{\text{shift}}\equiv\int\sigma_{zzz}^{\text{shift}}(\omega)d\omega=\frac{e^{3}}{\hbar^{2}}\frac{c^{2}}{ab}G(\tilde{t},\tilde{\delta}).$
(4)

where $\displaystyle G$ is also a dimensionless function of the dimensionless
parameters, with a global maximum 0.604 seeApp . Given this bound on
$\displaystyle\Sigma_{z}^{\text{shift}}$, it is clear that the single most
important factor in optimizing the second-order response of a 3D crystal is
the effective aspect ratio parameter $\displaystyle c^{2}/ab$.

We conclude by addressing the question of ultimate bounds on
$\displaystyle\Sigma^{\text{shift}}_{z}$, which is relevant not only to SHG,
but to sensitivity of photodetectors and efficiency of solar cells based on
Berry curvature generated intrinsic photogalvanic effects MN16 ;
MooreOrenstein  as well. We have discovered a new connection between
$\displaystyle\Sigma^{\text{shift}}_{z}$ and the modern theory of polarization
that generalizes the formulae for spectral weight beyond the RM Hamiltonian.
Specifically we show that $\displaystyle\Sigma^{\text{shift}}_{z}$ is
unbounded and potentially divergent when the possibility of next-neighbor
hopping is included.

As demonstrated in the Supplementary Information,
$\displaystyle\Sigma_{z}^{\text{shift}}$ in two band systems, regardless of
dimensionality $\displaystyle d$ and range of electron hopping amplitude is
given by,

$\Sigma_{z}^{\text{shift}}=\frac{2\pi e^{3}}{V\hbar^{2}}C_{3},$ (5)

where,

$C_{3}\ =\ -\frac{V}{(2\pi)^{d}}\int
d^{d}k\;\text{Im}\left[c_{3}-3c_{2}c_{1}+2c^{3}_{1}\right],$ (6)

$\displaystyle c_{n}\equiv\langle
u_{0}(k)|(i\partial_{k})^{n}|u_{0}(k)\rangle$, and the periodic gauge is
assumed for the valence-band Bloch wavefunction
$\displaystyle\psi_{0k}(r)=u_{0k}(r)e^{ikr}$. The quantity $\displaystyle
C_{3}$ is a member of a set of gauge invariant quantities, $\displaystyle
C_{n}$, that are cumulants of the electronic polarization hetenyi_cumulants ;
souza2000polarization . The quantity $\displaystyle C_{1}$ is exactly the
average macroscopic polarization, which coincides with the first moment of the
polarization distribution resta2007theory . Accordingly, $\displaystyle C_{3}$
is the third cumulant, or “skewness” of the distribution, which vanishes in
the presence of inversion symmetry. $\displaystyle C_{1}$ is controlled by the
center-of-mass location of the polarization, while $\displaystyle C_{3}$
describes the intrinsic asymmetry in the shape of the polarization
distribution, independent of its center-of-mass location.

In addition to providing a very satisfying connection between the nonlinear
response and ground state fluctuations of the polarization, equation (5)
speeds the search for effective Hamiltonians with large nonlinear response.
For instance, we find that the bound on
$\displaystyle\Sigma_{z}^{\text{shift}}$ that is obtained in the 3D RM model
is exceeded when next-neighbor hopping is introduced, and can diverge when the
fundamental gap is driven to zero. This result is consistent with the recent
analysis of Tan and Rappe tan2017upper , who suggested that
$\displaystyle\Sigma_{z}^{\text{shift}}$ can be enhanced when the range of
inter-site hopping becomes larger than the lattice parameter.

Figure 4: The Rice-Mele model. Depiction of the model of coupled Rice-Mele
chains (center) used to reproduce the experimentally observed second harmonic
generation in TaAs (left). A zoom on the shaded unit-cell is depicted on the
right specifying the different hopping terms present in the model.

In this work we reported SHG spectra of TaAs demonstrating a ten-fold resonant
enhancement of the response at 0.7 eV as compared with the previously
reported, and already quite large, response at 1.5 eV. We showed that the
spectral dependence and absolute magnitude of the SHG response function is
well captured by a model based on coupled ferroelectric Rice-Mele chains,
suggesting that the giant nonlinear response in TaAs is related to its strong
polar nature. The relationship between the RM phenomenology and the actual
electronic wavefunctions and bandstructure of TaAs remains to be understood.
Finally, we presented a new theorem linking the spectral weight of the
nonlinear conductivity to the third gauge-invariant cumulant of the ground
state. The cumulant $\displaystyle C_{3}$ is likely to be more computationally
tractable, and we anticipate that this theorem will facilitate algorithmic
searches for materials with large nonlinear response functions.

##  I Acknowledgements

Measurements and modeling were performed at the Lawrence Berkeley National
Laboratory in the Quantum Materials program supported by the Director, Office
of Science, Office of Basic Energy Sciences, Materials Sciences and
Engineering Division, of the U.S. Department of Energy under Contract No. DE-
AC02-05CH11231. Spectroscopy measurements were performed at the Department of
Physics, Temple University. The authors would like to thank B. Xu for sharing
ellipsometry data of TaAs. We would also like to thank Balasz Hetényi for
helpful communication. J.O. and L.W. received support for performing and
analyzing optical measurements from the Gordon and Betty Moore Foundation’s
EPiQS Initiative through Grant GBMF4537 to J.O. at UC Berkeley. Sample
synthesis by N.N. and J.A. was supported by the Gordon and Betty Moore
Foundation’s EPiQS Initiative through Grant GBMF4374. Work by N.N. and J.A.
was supported by the Office of Naval Research under the Electrical Sensors and
Network Research Division, Award No. N00014-15-1-2674. T.M. and J.E.M. were
supported by the Quantum Materials program under the U.S. Department of
Energy. T.M was supported by the Gordon and Betty Moore Foundation’s EPiQS
Initiative Theory Center Grant to UC Berkeley. A.G.G. was supported by the
Marie Curie Program under EC Grant agreement No. 653846. J.E.M. received
support for travel from the Simons Foundation. D.H.T. acknowledges startup
funds from Temple University. D.P. received support from the NSF GRFP, DGE
1752814. The authors would like to thank Nobumichi Tamura for his help in
performing crystal diffraction and orientation on beamline 12.3.2 at the
Advanced Light Source. N. Tamura and the ALS were supported by the Director,
Office of Science, Office of Basic Energy Sciences, of the U.S. Department of
Energy under Contract No. DE-AC02-05CH11231.

Competing Interests:  The authors declare that they have no competing
financial interests.

Correspondence:  Correspondence and requests for materials should be addressed
to D.H.T<EMAIL_ADDRESS>and J.O. (jworenstein@lbl.gov).

##  II Author Contribution

J.O., L.W., S.P., and D.H.T. conceived the project. D.H.T., S.P., L.W., B.L.,
M.R., and J.D.T. performed the SHG spectroscopy and calibration measurements.
J.O., S.P., L.W., and D.H.T. analyzed the data. J.O., A.G.G., T.M., D.P., and
J.E.M. developed the one-dimensional model and calculated theoretical
parameters. N.N. and J.A. grew the crystals and characterized the crystal
structure. J.O., S.P., L.W., D.P., T.M., A.G., and D.H.T. prepared the
manuscript with assistance from all authors.

##  Appendix A Data acquisition and processing

###  A.1 Extracting susceptibility components from SHG polar patterns

To calculate the components of the third rank susceptibility tensor
$\displaystyle\overleftrightarrow{\bm{\chi}}^{(2)}$ from the measured data, we
start with the relation defining the second harmonic susceptibility in terms
of the electric fields of the incident $\displaystyle\mathbf{E}^{1\omega}$ and
radiated $\displaystyle\mathbf{E}^{2\omega}$ fields:

$E^{2\omega}_{i}=\chi_{ijk}^{\text{shg}}E^{1\omega}_{j}E^{1\omega}_{k},$ (7)

where a sum over repeated indices is implied. We use the notation
$\displaystyle\hat{x},\hat{y},\hat{z}$, to correspond to the crystalline axes
$\displaystyle[1,0,0]$, $\displaystyle[0,1,0]$, $\displaystyle[0,0,1]$, and
for brevity, we shall drop the superscript “shg” in the following. The crystal
structure of TaAs corresponds to the point group $\displaystyle 4mm$, which
allows only three distinct components of the $\displaystyle\chi_{ijk}$ tensor:
$\displaystyle\chi_{zzz}$; $\displaystyle\chi_{zxx}=\chi_{zyy}$;
$\displaystyle\chi_{xxz}=\chi_{xzx}=\chi_{yyz}=\chi_{yzy}$. As noted in the
main text, the experiments were performed by measuring the second harmonic
intensity with light incident on the naturally occuring $\displaystyle(112)$
facet of a TaAs single crystal. The SHG intensity was measured in three
polarization “channels”:

  1. 1.

parallel: the polarization state of the analyzer is set to be parallel to the
polarization state of the incident radiation, while both rotate from
$\displaystyle 0^{\circ}$ to $\displaystyle 360^{\circ}$

  2. 2.

vertical: the polarization state of the analyzer is set to be parallel to the
in-plane polar direction $\displaystyle(1,1,-1)$ of the TaAs crystal, while
the polarization state of the incident radiation rotates from $\displaystyle
0^{\circ}$ to $\displaystyle 360^{\circ}$

  3. 3.

horizontal: the polarization state of the analyzer is set to be perpendicular
to the in-plane polar direction $\displaystyle(1,1,-1)$ of the TaAs crystal,
while the polarization state of the incident radiation rotates from
$\displaystyle 0^{\circ}$ to $\displaystyle 360^{\circ}$

It is useful to label the in-plane component of the polar axis
$\displaystyle[1,1,-1]$ as the axis “$\displaystyle\gamma$”, and
correspondingly label the perpendicular in-plane direction
$\displaystyle[-1,1,0]$ as the axis “$\displaystyle\alpha$”. If the incident
light has intensity $\displaystyle I_{0}$, with linearly polarized electric
field expressed by
$\displaystyle\mathbf{E}^{1\omega}=\sqrt{I_{0}}\left(\hat{\alpha}\sin\theta+\hat{\gamma}\cos\theta\right)$,
in the parallel channel, the intensity can be expressed in terms of the angle
$\displaystyle\theta$ of the incident polarization state with respect to the
“$\displaystyle\gamma$” axis:

$\displaystyle\displaystyle I_{\text{para}}(\theta)$
$\displaystyle\displaystyle=\left|E^{2\omega}_{\alpha}\sin\theta+E^{2\omega}_{\gamma}\cos\theta\right|^{2}$
$\displaystyle\displaystyle=\left|\frac{1}{\sqrt{2}}\left(-E^{2\omega}_{x}+E^{2\omega}_{y}\right)\sin\theta+\frac{1}{\sqrt{3}}\left({E^{2\omega}_{x}+E^{2\omega}_{y}-E^{2\omega}_{z}}\right)\cos\theta\right|^{2}.$

We wish to express this in terms of the tunable quantities $\displaystyle
E^{1\omega}_{z}=-\tfrac{E^{1\omega}_{\gamma}}{\sqrt{3}}$ and $\displaystyle
E^{1\omega}_{\gamma}=\sqrt{I_{0}}\cos\theta$. To do this, we first expand
$\displaystyle\bm{E}^{2\omega}$ in terms of $\displaystyle\bm{E}^{1\omega}$
via

$E^{2\omega}_{x}=2\chi_{xxz}E^{1\omega}_{x}E^{1\omega}_{z},\qquad
E^{2\omega}_{y}=2\chi_{yyz}E^{1\omega}_{y}E^{1\omega}_{z},\qquad
E^{2\omega}_{z}=\chi_{zxx}\left(E^{1\omega}_{x}\right)^{2}+\chi_{zyy}\left(E^{1\omega}_{y}\right)^{2}+\chi_{zzz}\left(E^{1\omega}_{z}\right)^{2},$
(8)

and transform back from the $\displaystyle xyz$ to the
$\displaystyle\alpha\gamma$ basis via

$E^{1\omega}_{x}=-\frac{E^{1\omega}_{\alpha}}{\sqrt{2}}+\frac{E^{1\omega}_{\gamma}}{\sqrt{3}},\qquad
E^{1\omega}_{y}=\frac{E^{1\omega}_{\alpha}}{\sqrt{2}}+\frac{E^{1\omega}_{\gamma}}{\sqrt{3}},\qquad
E^{1\omega}_{z}=-\frac{E^{1\omega}_{\gamma}}{\sqrt{3}}.$ (9)

This yields

$\displaystyle\displaystyle I_{\text{para}}(\theta)$
$\displaystyle\displaystyle=\frac{I_{0}^{2}}{3}\left|\chi_{xxz}\cos\theta\sin^{2}\theta+\left(\frac{4}{3}\chi_{xxz}\cos^{2}\theta+\chi_{zxx}\left(\sin^{2}\theta+\frac{2}{3}\cos^{2}\theta\right)+\frac{1}{3}\chi_{zzz}\cos^{2}\theta\right)\cos\theta\right|^{2}$
$\displaystyle\displaystyle=\frac{I_{0}^{2}}{27}\left|3\chi_{xxz}\cos\theta\sin^{2}\theta+\left(4\chi_{xxz}\cos^{2}\theta+\chi_{zxx}\left(3\sin^{2}\theta+2\cos^{2}\theta\right)+\chi_{zzz}\cos^{2}\theta\right)\cos\theta\right|^{2}$
$\displaystyle\displaystyle=\frac{I_{0}^{2}}{27}\left|(6\chi_{xxz}+3\chi_{zxx})\cos\theta\sin^{2}\theta+\left(4\chi_{xxz}+2\chi_{zxx}+\chi_{zzz}\right)\cos^{3}\theta\right|^{2}.$

We define $\displaystyle\chi_{\text{eff}}=\chi_{zzz}+2\chi_{zxx}+4\chi_{xxz}$,
which gives us the expression,

$I_{\text{para}}=\frac{I_{0}^{2}}{27}\left|3(2\chi_{xxz}+\chi_{zxx})\cos\theta\sin^{2}\theta+\chi_{\text{eff}}\cos^{3}\theta\right|^{2}.$
(10)

We can similarly obtain expressions for the other two channels:

$\displaystyle\displaystyle
I_{\text{vertical}}(\theta)=\left|E^{2\omega}_{3}\right|^{2}=\left|\frac{1}{\sqrt{3}}\left({E^{2\omega}_{x}+E^{2\omega}_{y}-E^{2\omega}_{z}}\right)\right|^{2}=\frac{1}{27}\left(3\chi_{zxx}\sin^{2}\theta+\chi_{\text{eff}}\cos^{2}\theta\right),$
(11)

and

$\displaystyle\displaystyle
I_{\text{horiz}}(\theta)=\left|E^{2\omega}_{\alpha}\right|^{2}=\left|\frac{1}{\sqrt{2}}\left(-E^{2\omega}_{x}+E^{2\omega}_{y}\right)\right|^{2}=\frac{1}{3}I_{0}^{2}|\chi_{xxz}|^{2}\;\sin
2\theta.$ (12)

Figure 5: Minimized mean-squared deviation fits (solid lines) to the SHG
intensity data (points) as a function of the incident polarization angle
$\displaystyle\theta$, in three polarization channels (a) parallel, (b)
vertical, (c) horizontal, for incident photon wavelength
$\displaystyle\lambda=1400$nm, plotted on a normalized scale. (d) Fits to the
SHG intensity in the vertical channel assuming the relative sign between
$\displaystyle\chi_{zxx}$ and $\displaystyle\chi_{zzz}$ to be
$\displaystyle(+)$ (dashed line), and $\displaystyle(-)$ (solid line).

Figure 5d shows the SHG intensity data as a function of the polarization of
incident photons in the parallel, vertical, and horizontal channel, with fits
to the expressions in equation (10), (11), (12) respectively. The amplitudes
of the fitting parameters
$\displaystyle|\chi_{xxz}|,|\chi_{zxx}|,|\chi_{\text{eff}}|$ thus obtained are
then multiplied by various wavelength-dependent correction factors described
below.

Although the fitting scheme described above is not sufficient to obtain the
complete amplitude and phase information for the three components
$\displaystyle|\chi_{xxz}|,|\chi_{zxx}|,|\chi_{\text{eff}}|$, we are
nevertheless able to definitively say that for all incident photon energies,
$\displaystyle\chi_{zxx}$ and $\displaystyle\chi_{zzz}$ have a relative phase
of $\displaystyle\pi$. Figure 5 shows the best fits to the data at incident
photon wavelength $\displaystyle\lambda=$\displaystyle
1400\text{\,}\mathrm{nm}$$ in the “vertical” channel, using two contrasting
assumptions for the relative signs of the fitting parameters,
$\displaystyle(+1)$ or $\displaystyle(-1)$. The fits with relative
$\displaystyle(-)$ signs are are observed to be better approximations to the
data.

With the additional information about the relative signs, the calculated
values of
$\displaystyle|\chi_{\text{eff}}|=|\chi_{zzz}+2\chi_{zxx}+4\chi_{xxz}|$ can be
used to place bounds on the values of $\displaystyle|\chi_{zzz}|$ with an
upper bound (UB) and lower bound (LB) expressed as (see figure 3 in the main
text):

$\displaystyle\displaystyle|\chi_{zzz}^{UB}|$
$\displaystyle\displaystyle=|\chi_{\text{eff}}|+2|\chi_{zxx}|+4|\chi_{xxz}|$
(13) $\displaystyle\displaystyle|\chi_{zzz}^{LB}|$
$\displaystyle\displaystyle=|\chi_{\text{eff}}|+2|\chi_{zxx}|-4|\chi_{xxz}|.$
(14)

###  A.2 Wavelength dependent correction factors

We provide a brief summary of the factors that must necessarily be taken into
consideration to accurately determine the second harmonic susceptibility
tensor elements $\displaystyle\chi_{ijk}^{\text{shg}}$ and the procedure we
followed to experimentally measure them.

The electric field of the incident laser pulse can be well approximated by a
simple Gaussian field as

$\mathbf{E}_{\omega}(\mathbf{r},t)=E_{0}e^{-i(k_{0}z-\omega
t)}\exp\left(-\frac{x^{2}+y^{2}}{w_{0}^{2}}\right)\exp{\left[-\left(\frac{z-ct}{\tau}\right)^{2}\right]}\hat{x}$
(15)

where $\displaystyle E_{0}$ is the electric field amplitude, $\displaystyle
k_{0}$ is the field wavevector magnitude, $\displaystyle\omega$ is the
fundamental light frequency, $\displaystyle c$ is the speed of light, and
$\displaystyle w_{0}$ is the electric field beam waist Quesnel1998theory ,
here taken as the point from which the SHG field arises. We have assumed that
the beam is transversely polarized purely in the $\displaystyle\hat{x}$
direction for convenience. As a consequence of Poynting’s Theorem, the
relationship between the integrated intensity of the incident, fundamental
frequency field $\displaystyle I_{\omega}$ and the associated electric field
can be given by

$I_{\omega}=\frac{1}{2}\epsilon_{0}\int\left\langle|\mathbf{E}_{\omega}(\mathbf{r},t)|^{2}\right\rangle
d\mathbf{r}$ (16)

where the brackets denote time averaging. For the field as given by equation
(15), the integration is easily performed to yield

$I_{\omega}\propto E_{0}^{2}w_{0}\tau$ (17)

The penetration depth of the wavelengths used in this study is relatively
short in TaAs. Thus we may consider the limit of non-depleted incident fields.
The radiated second harmonic electric field is then proportional to the
second-order induced polarization as
$\displaystyle\mathbf{E}_{2\omega}(\mathbf{r},t)\propto\mathbf{P}_{2\omega}(\mathbf{r},t)=\epsilon_{0}\overleftrightarrow{\bm{\chi}}^{\text{shg}}\mathbf{E}_{\omega}(\mathbf{r},t)\mathbf{E}_{\omega}(\mathbf{r},t)$.
Hence, the intensity of the second harmonic field $\displaystyle I_{2\omega}$
is given by

$I_{2\omega}=\frac{1}{2}\epsilon_{0}\int\left\langle|\mathbf{E}_{2\omega}(\mathbf{r},t)|^{2}\right\rangle
d\mathbf{r}\propto\int\left\langle|\overleftrightarrow{\bm{\chi}}^{\text{shg}}|^{2}|\mathbf{E}_{\omega}(\mathbf{r},t)|^{4}\right\rangle
d\mathbf{r}\propto(E_{0}^{2}\chi_{ijk}^{\text{shg}})^{2}w_{0}\tau.$ (18)

When this latter result is combined with equation (17) to eliminate the
electric field $\displaystyle E_{0}$, we determine that $\displaystyle
I_{2\omega}\propto(\chi_{ijk}^{\text{shg}})^{2}I_{\omega}^{2}/w_{0}\tau$.
However, for Gaussian fields such as in equation (15), the beam waist at the
focus $\displaystyle w_{0}$ can be described in terms of the incident beam
diameter $\displaystyle d$, lens focal length $\displaystyle f$ and beam
wavelength $\displaystyle\lambda$ as $\displaystyle w_{0}\propto f\lambda/d$,
yielding

$I_{2\omega}\propto(\chi_{ijk}^{\text{shg}})^{2}\frac{I_{\omega}^{2}d}{f\lambda\tau}.$
(19)

Therefore, accurate measurement of the elements of
$\displaystyle\chi^{\text{shg}}$ must be controlled for incident integrated
pulse intensity (equivalently, energy), diameter, wavelength and duration as
well as the extrinsic factor of the lens focal length. Since the value for
$\displaystyle\chi_{ijk}^{\text{shg}}$ at $\displaystyle
800\text{\,}\mathrm{nm}$ is known from previous study wu2017giant , all
measurements were normalized to our measured response at $\displaystyle
800\text{\,}\mathrm{nm}$. In addition, the relative enhancement of second
harmonic response between TaAs and GaAs was checked independently using a two
color Er:doped fiber laser with source photon wavelengths of 1550nm and 780nm
respectively.

###  A.3 Experimental measurement of correction factors

####  A.3.1 Intrinsic factors

Figure 6: Measured FROG spectrum for input fields $\displaystyle\lambda=1100$
nm (nominally) showing the projection of the data along the temporal (above)
and spectral (right) axes. Fits of the projected data to Guassian models, as
described in the text, produced the observed pulse duration
$\displaystyle\tau$, the center SHG wavelength $\displaystyle\lambda$ and the
spectral width of the pulse $\displaystyle\Delta\lambda$.

Even though a commercially supplied wavelength separator was employed, a
number of parasitic wavelengths were present due to a combination of leakage
of other responses from the Optical Parametric Amplifier, as well as their
sum-frequency, difference-frequency, and second harmonic responses among the
signal, idler and second stage pump fields. At each individual wavelength, we
computed these parasitic wavelengths and assembled a combination of longpass
and shortpass filters of minimum ND8 to remove the parasitics before measuring
the incident power. This power was recorded after the last longpass filter (LF
in figure 1 of the main text), and a set of reflective neutral density filters
of value ND chosen to maintain a calculated incident fluence of roughly
$\displaystyle\sim$ 20 mJ/cm2. As the measured average power $\displaystyle P$
is simply proportional to the energy per pulse through a factor of the
repetition rate, the value $\displaystyle P/10^{ND}$ was divided from the
measured power to provide a measure of the incident integrated intensity
$\displaystyle I_{\omega}$.

The spectral and temporal characteristics of the OPA beam were measured using
a homebuilt Frequency Resolved Optical Gating (FROG) apparatus over the range
from incident wavelength $\displaystyle\lambda=800$ nm to
$\displaystyle\lambda=2000$ nm. In all cases, the FROG measurement was
performed on the beam before its entrance into the reflective objective (i.e.,
in between LF and RO in figure 1 of the main text) in order to account for the
filtering and dispersive characteristics of the preceding optics. While we did
not use the same ND filters for the FROG measurement as we did in the
experiments, we did use the same number of filters of identical material (UV
grade fused silica) and thickness as those used in data collection. A
representative FROG trace is shown in figure 6. The FROG apparatus used a
$\displaystyle 100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick
$\displaystyle\beta$-Barium Oxide (BBO) crystal, the radiated second harmonic
spectrum was determined via simulation with lab2.de softwareLab2  to be
virtually identical to that which would derive from a thin layer of a bulk
response reflection experiment as occured here on TaAs.

The pulse duration of the incident second harmonic field was determined from
numerically summing the FROG trace along the spectral dimension and fitting
the resulting data to a Gaussian of the form $\displaystyle
I_{2\omega}(t)=I_{0}\exp\left[-\left((t-t_{0})/\tau\right)^{2}\right]$, as
shown in the top panel of figure 6. We note that from evaluation of the
convolution of a Gaussian with itself, the pulse duration at $\displaystyle
2\omega$ is related to that at $\displaystyle\omega$ through a numerical
factor as $\displaystyle\tau_{2\omega}=\tau_{\omega}/\sqrt{2}$. In order to
determine the radiated central $\displaystyle 2\omega$ wavelength and
bandwidth, we numerically summed the FROG trace along the temporal dimension
and fit the result to a Gaussian of the form $\displaystyle
I_{2\omega}(\lambda)=I_{0}\exp\left[-\left((\lambda-\lambda_{0})/\Delta\lambda\right)^{2}\right]$,
as shown on the righthand panel. These data allowed us to determine the
fundamental wavelength as $\displaystyle 2\lambda_{0}$ and the spectral
bandwidth of the SHG pulse $\displaystyle\Delta\lambda$.

Figure 7: (a) Representative data from knife-edge measurement of the beam
diameter and corresponding fit for $\displaystyle\lambda_{0}=1100$ nm incident
light. (b) Experimentally measured transmittance of the 540 nm bandpass
filter. (c) The linear optical correction factor $\displaystyle BP(\omega)$ as
a function of incident photon energy.

The knife-edge method was used to measure the beam diameter $\displaystyle
w_{0}$ immediately before the reflective objective as a function of incident
wavelength $\displaystyle\lambda$ for all wavelengths in the study.
Representative data at $\displaystyle\lambda=1100$ nm and the corresponding
fit to the complementary error function $\displaystyle
I(x)=I_{0}\operatorname{erfc}{\left(\sqrt{2}(x-x_{0})/w_{0}\right)}$ are shown
in figure 7. A basic Gaussian beam propagation could be used to determine the
beam at the focus. However, this was unnecessary since the identical, all-
reflective optical setup was used for every data point, necessitating only
measurement of the beam diameter $\displaystyle w_{0}$.

####  A.3.2 Extrinsic factors

Experimentally, the measured second harmonic intensity is also proportional to
the spectral response of the detector $\displaystyle\mathcal{D}(\lambda)$ and
the transmittance $\displaystyle\mathcal{T}(\lambda)$ of the filters in front
of the detector.

The detector spectral response $\displaystyle\mathcal{D}(\lambda)$ was
determined from the manufacturer specifications for each detector and divided
out from the measured intensity. In the range from $\displaystyle\lambda=800$
nm to $\displaystyle\lambda=1600$ nm, we used a Hamamatsu R12829 in a
Hamamatsu C12597-01 socket for power supply. The photomultiplier output was
directly connected to the current input of a Zurich Instruments MFLI
digitizing lockin which yielded a current value proportional to the incident
SHG power. The range from $\displaystyle\lambda=1560$ nm to
$\displaystyle\lambda=1700$ nm was measured with an unbiased Thorlabs FDS010
photodiode, while measurements from $\displaystyle\lambda=1600$ nm to
$\displaystyle\lambda=2600$ nm performed with a Thorlabs FGA01 InGaAs
photodiode. For their respective wavelength ranges, the photodiodes were
connected to the input of a Cremat CR-Z-110 charge preamplifier followed by a
CR-S-8$\displaystyle\mu$s shaping amplifier whose output were connected to the
voltage input of the lockin. The overlap between the wavelength ranges of the
three detectors allowed us to account for the current to voltage conversion of
the charge integrating electronics in relation to the PMT, permitting the Si
photodiode to be a bridge between the PMT and the InGaAs detector.

In all cases, we used a set of two shortpass filters and a single bandpass
filter centered at the nominal second harmonic wavelength, which determined
the spectral transmittance $\displaystyle\mathcal{T}(\lambda)$. These filters
were chosen in order to provide attenuation of $\displaystyle>$ND10 at the
fundamental wavelength while not attenuating the SHG by more than a factor of
2. We verified the transmittance of the filters to match the supplier’s
specifications through individual UV-Vis measurements, as shown in figure 7.
It was necessary to account for the fact that the bandwidth of the filter was
more narrow than that of the second harmonic pulse, in particular at the
longest wavelengths measured. Thus, after dividing by
$\displaystyle\mathcal{D}(\lambda)$, the true incident SHG power was
determined by dividing out the numerically integrated product of the measured
filter data and a normalized Gaussian representing the fitted parameters of
the measured second harmonic spectrum the FROG traces as

$\frac{1}{\sqrt{\pi}\Delta\lambda}\int_{-\infty}^{\infty}\mathcal{T}(\lambda)\exp\left[-((\lambda-\lambda_{0})/\Delta\lambda)^{2}\right]d\lambda.$
(20)

###  A.4 Translating reflection geometry parameters to bulk nonlinear
conductivity

To calculate the bulk nonlinear optical parameters
$\displaystyle\chi^{\text{shg}}$, and $\displaystyle\sigma^{\text{shg}}$, we
also need to correct for the linear optical parameters of TaAs. We use the
formula derived by Bloembergen and Pershanbloembergen1962light ; wu2017giant
to account for the linear optical properties such as the refractive index and
transmission coefficients:

$\chi_{R}^{\text{shg}}=\frac{\chi^{\text{shg}}}{\left(\sqrt{\epsilon(2\omega)}+\sqrt{\epsilon(\omega)}\right)\left(\sqrt{\epsilon(2\omega)}+1\right)}T(\omega)^{2}$
(21)

where $\displaystyle\chi_{R}^{\text{shg}}$ denotes the nonlinear
susceptibility in reflection geometry, as measured in the experiment,
$\displaystyle\epsilon(\omega)$, and $\displaystyle T(\omega)$, are
respectively the dielectric constant and the transmission coefficients of
TaAs.

We define the Bloembergen-Pershan correction coefficient $\displaystyle
BP(\omega)\equiv\frac{T(\omega)^{2}}{\left(\sqrt{\epsilon(2\omega)}+\sqrt{\epsilon(\omega)}\right)\left(\sqrt{\epsilon(2\omega)}+1\right)}$,
which can be calculated as a function of incident photon energy based on
values of the dielectric constant as measured by ellipsometry.111The
dielectric constant was measured using ellipsometry of TaAs by Bing Xu at
Université de Fribourg, and obtained in private communication. (See figure 7.)
Optical conductivity is more conveniently modeled than susceptibility, and can
be obtained using the relation:
$\displaystyle\sigma^{\text{shg}}_{ijk}(\omega)=-i\epsilon_{0}\,\omega\,\chi^{\text{shg}}_{ijk}(\omega)$

##  Appendix B Quasi-One Dimensional Model for Second Harmonic Generation

###  B.1 Model details

This supplement describes the details of the phenomenological model of Rice-
Mele chains used in the main text, and details the calculation of its second
harmonic response. The model is inspired by the polar character of the TaAs
lattice as well as the observation that $\displaystyle\\\ {\sigma_{zzz}}$ is
the dominant component of the SHG response. This is reminiscent of a 1D
system, in which $\displaystyle\sigma_{zzz}$ is the only non-zero component.
We consider a quasi-one dimensional model composed of 1D Rice Mele chains in
the $\displaystyle z$-direction RM82  with only weak couplings in the
$\displaystyle xy$ plane. In momentum space it takes the form of a two band
model $\displaystyle H_{\bm{k}}=\bm{d}_{\bm{k}}\cdot\bm{\sigma}$ with

$\displaystyle d_{x}$ $\displaystyle=t\sin(ck_{z}/2)$ (22) $\displaystyle
d_{y}$ $\displaystyle=\delta\cos(ck_{z}/2)$ $\displaystyle d_{z}$
$\displaystyle=\Delta+t_{AB}(\cos(bk_{x})+\cos(bk_{y})).$

Here $\displaystyle b$ and $\displaystyle c$ are the dimension of the unit
cell in the $\displaystyle x$ and $\displaystyle z$ directions. Along lines in
the $\displaystyle z$-direction there is a staggered onsite potential
$\displaystyle\pm\Delta$ and staggered hopping $\displaystyle t\pm\delta$. The
parameter $\displaystyle t_{AB}=t_{\parallel,A}-t_{\parallel,B}\ll\Delta$
represents the difference between the inter-chain hopping strengths. When
$\displaystyle t_{AB}=0$, this reduces to an ensemble of independent Rice-Mele
models. The overall energy scale is set by fixing $\displaystyle\Delta$; there
are three independent parameters of the model:
$\displaystyle\tilde{t}=t/\Delta$,
$\displaystyle\tilde{\delta}=\delta/\Delta$, and
$\displaystyle\tilde{t}_{AB}=t_{AB}/\Delta$.

###  B.2 Calculation of SHG Response

The second-order conductivity tensor is defined via

$J_{a}(\omega_{0})=\sum_{b,c}\sigma_{abc}(\omega_{0};\omega_{1},\omega_{1})E_{b}(\omega_{1})E_{c}(\omega_{2})$
(23)

where $\displaystyle\omega_{0}=\omega_{1}+\omega_{2}$ is the frequency of the
emitted photon The second harmonic response is
$\displaystyle\sigma_{abc}^{\text{shg}}(\omega)=\sigma_{abc}(2\omega;\omega,\omega)$.
We will also frequently invoke the shift current,
$\displaystyle\sigma_{abc}^{\text{shift}}(\omega)=\sigma_{abc}(0;-\omega,\omega)$.
The shift current may be thought of as the “solar panel response”: the amount
of DC current generated under illumination.

In two band models, the diagonal components of $\displaystyle\sigma_{abc}$ are
particularly simple. For $\displaystyle\omega>0$,MN16 ; yang2017divergent

$\displaystyle\displaystyle\operatorname{Re}\sigma_{aaa}^{\text{shift}}(\omega)\
$ $\displaystyle\displaystyle=\ \frac{2\pi
e^{3}}{\hbar^{2}\omega^{2}}\int[d\bm{k}]f_{10}v_{01}^{a}w_{10}^{aa}\delta(\omega_{10}-\omega)$
(24)
$\displaystyle\displaystyle\operatorname{Re}\sigma_{aaa}^{\text{shg}}(\omega)\
$ $\displaystyle\displaystyle=\ \frac{\pi
e^{3}}{2\hbar^{2}\omega^{2}}\int[d\bm{k}]\Big{[}f_{10}v_{01}^{a}w_{10}^{aa}\delta(\omega_{10}-2\omega)-2f_{10}v_{01}^{a}w_{10}^{aa}\delta(\omega_{10}-\omega)\Big{]}$
(25)

where $\displaystyle\int[d\bm{k}]=\int\frac{d^{d}\bm{k}}{(2\pi)^{d}}$ is the
normalized integral over the Brillouin zone, $\displaystyle 0$ and
$\displaystyle 1$ refer to the valence and conduction bands respectively,
$\displaystyle f_{01}$ is the difference in Fermi factors,
$\displaystyle\omega_{01}$ is the difference in band frequencies,
$\displaystyle v^{a}=\partial_{k_{a}}H$ is the velocity operator, and
$\displaystyle w^{aa}=\partial_{k_{a}}\partial_{k_{a}}H$. Taken together,
these lead to the convenient identity (see equation (3) of the main text)

$\operatorname{Re}\sigma_{aaa}^{\text{shg}}(\omega)=-\frac{1}{2}\operatorname{Re}\sigma_{aaa}^{\text{shift}}(\omega)+\operatorname{Re}\sigma_{aaa}^{\text{shift}}(2\omega).$
(26)

Again, these equations only hold for two band models.

We now evaluate these equations for our specific model. We focus on the two-
photon term of the SHG. Differentiating equation (22),

$\hat{v}^{z}=\frac{c}{\hbar
2}\left[-t\sin\left(\frac{k_{z}c}{2}\right)\sigma^{x}+\delta\cos\left(\frac{k_{z}c}{2}\right)\sigma^{y}\right]$
(27)

and

$\hat{w}^{zz}=\frac{c^{2}}{4\hbar}\left[-t\cos\left(\frac{k_{z}c}{2}\right)\sigma^{x}-\delta\sin\left(\frac{k_{z}c}{2}\right)\sigma^{y}\right].$
(28)

Hence, using
$\displaystyle\sigma^{i}\sigma^{j}=i\varepsilon^{ij}_{k}\sigma^{k}$,

$\hat{v}^{z}\hat{w}^{zz}=\frac{1}{\hbar^{2}}\left(\frac{c}{2}\right)^{3}\left[\frac{t^{2}-\delta^{2}}{2}\sin\left(\frac{2k_{z}c}{2}\right)+it\delta\right]\sigma^{z}.$
(29)

We must now evaluate the matrix element $\displaystyle\langle
0|\sigma^{z}|0\rangle$. Adopting the Bloch sphere representation,

$\bm{d}=d\left(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta\right)$
(30)

which leads to the eigenvector equation $\displaystyle
H|0\rangle=\varepsilon_{0}|0\rangle$ with

$|0\rangle=\begin{pmatrix}\cos\frac{\theta}{2}\\\
\sin\frac{\theta}{2}e^{i\varphi}\end{pmatrix},\qquad$ (31)

and eigenvalue $\displaystyle\varepsilon_{v}=-d$. Therefore

$\langle
0|\sigma^{z}|0\rangle=\left(\cos^{2}\frac{\theta}{2}-\sin^{2}\frac{\theta}{2}\right)=\cos\theta=\frac{d_{z}}{d}.$
(32)

Combining (29) and (32) shows that the two-photon contribution is

$\operatorname{Re}\sigma_{zzz}^{\text{shg2p}}(\omega)=\frac{\pi
e^{e}}{2\hbar^{2}\omega^{2}}\int[d\bm{k}]f_{10}v_{01}^{z}w_{10}^{zz}\delta(\omega_{01}-2\omega)=i\frac{\pi
e^{3}}{2\hbar^{2}\omega^{2}}\int[d\bm{k}](-1)\frac{1}{\hbar^{2}}\left(\frac{c}{2}\right)^{3}it\delta\frac{d_{z}}{d}\delta(\omega_{10}-2\omega).$
(33)

Note that the $\displaystyle\sin(k_{z}c)$ term integrates to zero and can be
ignored. What remains is to de-dimensionalize the integral and simplify the
expression somewhat. Define
$\displaystyle\varepsilon(k)=d(k)=\hbar\omega_{10}/2$. Then

$\operatorname{Re}\sigma_{zzz}^{\text{shg2p}}(\omega)=\frac{\pi
e^{3}}{2\hbar^{2}}\left(\frac{c}{2}\right)^{3}t\delta\int\frac{d^{3}\bm{k}}{\left(2\pi\right)^{3}}\frac{d_{z}}{\varepsilon(k)\omega^{2}}\frac{\hbar}{2}\delta(\varepsilon(k)-\hbar\omega).$
(34)

To account for physical broadening of the peak, we use a Lorentzian
$\displaystyle\delta(\varepsilon-\omega)\to\frac{1}{\pi}\frac{\gamma}{\left(\varepsilon-\omega\right)^{2}+\gamma^{2}}$.
Furtheremore, we convert the remaining factors of $\displaystyle\omega$ into
$\displaystyle\varepsilon$ to remove an unphysical diverges at
$\displaystyle\omega\to 0$ in numerical evaluation. So

$\operatorname{Re}\sigma_{zzz}^{\text{shg2p}}(\omega)=\frac{\pi
e^{3}}{\hbar}\left(\frac{c}{2}\right)^{3}t\delta\frac{1}{4\pi}\int[d\bm{k}]\frac{d_{z}}{\varepsilon(k)^{3}}\frac{\gamma}{\left(\varepsilon(k)-\omega\right)^{2}+\gamma^{2}}.$
(35)

Now we de-dimensionalize by normalizing energies to $\displaystyle\Delta$, and
henceforth represent $\displaystyle\Delta$-normalizes quantities with a tilde.
For example, $\displaystyle\tilde{t}\equiv t/\Delta$. Furthermore, let us
define $\displaystyle X\equiv k_{x}b,Y\equiv k_{y}b,Z\equiv k_{z}c$. This
gives

$\displaystyle\displaystyle\widetilde{\varepsilon}=\sqrt{\tilde{\delta}^{2}\cos(Z/2)^{2}+\tilde{t}^{2}\sin(Z/2)^{2}+\widetilde{d_{z}}^{2}}\text{
and }\widetilde{d_{z}}=1+\tilde{t}_{AB}\left[\cos(X)+\cos(Y)\right].$ (36)

Thus

$\displaystyle\displaystyle\sigma_{zzz}^{\text{shg2p}}(\omega)$
$\displaystyle\displaystyle=\frac{e^{3}}{\hbar}\left(\frac{c}{2}\right)^{3}\frac{\widetilde{t}\widetilde{\delta}}{\Delta}\frac{1}{4}\frac{1}{(2\pi)^{3}abc}\;\int_{-\pi}^{\pi}dX\,dY\,dZ\,\frac{\widetilde{d}_{z}}{\tilde{\varepsilon}(k)^{3}}\frac{\tilde{\gamma}}{\left(\tilde{\omega}-\tilde{\varepsilon}(k)\right)^{2}+\tilde{\gamma}^{2}}.$
(37)

To match with experimental data, one should multiply by degeneracy factors
that account for the spin degeneracy $\displaystyle g_{s}=2$ and orbital
degeneracy $\displaystyle g_{O}=4$ for the number of Rice-Mele chains per unit
cell. Rearranging, we arrive at the two-photon contribution to the SHG

$\sigma_{zzz}^{\text{shg2p}}(\omega)=\left[\frac{e^{2}}{\hbar}\frac{e}{\Delta}\right]\frac{c^{2}}{b^{2}}\left(\frac{1}{4\pi}\right)^{3}F(\tilde{\omega};\tilde{\delta},\tilde{t},\tilde{t}_{AB},\tilde{\gamma})$
(38)

where the term in brackets has units of conductance per Volt, and the integral
has been re-packaged as

$\displaystyle\displaystyle
F(\tilde{\omega};\tilde{\delta},\tilde{t},\tilde{t}_{AB},\tilde{\gamma})$
$\displaystyle\displaystyle=g_{s}g_{O}\frac{\tilde{\delta}\tilde{t}}{4}\int_{-\pi}^{\pi}dX\,\int_{-\pi}^{\pi}dY\,\int_{-\pi}^{\pi}dZ\;\frac{1+\tilde{t}_{AB}\left[\cos(X)+\cos(Y)\right]}{\left(\tilde{\delta}^{2}\cos(Z/2)^{2}+\tilde{t}^{2}\sin(Z/2)^{2}+\left(1+\tilde{t}_{AB}\left[\cos(X)+\cos(Y)\right]\right)^{2}\right)^{3/2}}$
$\displaystyle\displaystyle\hskip 20.00003pt\times\
\frac{\tilde{\gamma}}{\left(\tilde{\omega}-\left[\tilde{\delta}^{2}\cos(Z/2)^{2}+\tilde{t}^{2}\sin(Z/2)^{2}+\left(1+\tilde{t}_{AB}\left[\cos(X)+\cos(Y)\right]\right)^{2}\right]^{1/2}\right)^{2}+\tilde{\gamma}^{2}}.$
(39)

By Equation (26),
$\displaystyle\operatorname{Re}\sigma^{\text{shift}}_{zzz}(\omega)=-2\sigma_{zzz}^{\text{shg2p}}(\omega)$
(see equation (2) of the main text).

It is worth examining the special case where $\displaystyle t_{AB}=0$, i.e. an
ensemble of coupled Rice-Mele chains. Here the integral becomes analytically
tractable. Converting the Lorentzian back to a $\displaystyle\delta$ function,

$\displaystyle\displaystyle\int_{0}^{\infty}d\omega\;F\left(\widetilde{\omega};\tilde{\delta},\tilde{t},\tilde{t}_{AB}\to
0,\tilde{\gamma}\to 0\right)$ (40) $\displaystyle\displaystyle=\
\frac{\pi\Delta
g_{s}g_{O}}{\hbar}\left(2\pi\right)^{2}\frac{\tilde{\delta}\tilde{t}}{4}\int_{-\pi}^{\pi}\frac{dZ}{\left(\tilde{\delta}^{2}\cos(Z/2)^{2}+\tilde{t}^{2}\sin(Z/2)^{2}+1\right)^{3/2}}$
(41) $\displaystyle\displaystyle=\
g_{s}g_{0}\frac{\Delta}{\hbar}\frac{2(2\pi)^{3}\tilde{\delta}\tilde{t}}{4}\frac{E_{2}\left(\frac{\tilde{\delta}^{2}-\tilde{t}^{2}}{1+\tilde{\delta}^{2}}\right)}{\left(1+\tilde{t}^{2}\right)\sqrt{1+\tilde{\delta}^{2}}},$
(42)

where $\displaystyle E_{2}$ is the complete elliptic integral of the second
kind. Therefore, using the fact that the spectral weights of the shift current
and two-photon contribution to the SHG are equal,

$\Sigma^{\text{shift}}_{z}=2\Sigma^{\text{shg
2p}}_{z}=\left[\frac{e^{3}}{\hbar^{2}}\right]\frac{c^{2}}{b^{2}}G(\tilde{\delta},\tilde{t})\text{
where }G(\tilde{\delta},\tilde{t})\equiv
g_{s}g_{O}\frac{\tilde{\delta}\tilde{t}}{8}\frac{E_{2}\left(\frac{\tilde{\delta}^{2}-\tilde{t}^{2}}{1+\tilde{\delta}^{2}}\right)}{\left(1+\tilde{t}^{2}\right)\sqrt{1+\tilde{\delta}^{2}}}.$
(43)

The function $\displaystyle G$ is bounded and attains its global maximum at
$\displaystyle G(\sqrt{2},\sqrt{2})=\pi 3^{-3/2}\approx 0.604$ (see equation
(4) of the main text). This shows that the total spectral weight in the two-
photon contribution is bounded, i.e. there is a maximum amount of SHG for
Rice-Mele models.

##  Appendix C Gauge-Invariant Cumulants and a new sum rule

This appendix reviews the theoretical machinery of gauge-invariant cumulants
(GICs) that underlies the relationship between ground-state polarization
distributions and the frequency-integrated nonlinear response described in the
last part of the main text. We restrict the analysis to two-band tight-binding
models, relevant for this work, to show that this connection can be made
exact. It generalizes the known relationships between the Berry connection and
polarization, as well as the relation between the spread of Wannier functions
with linear response. The last part of the supplement shows how this
connection gives a guide to construct of models whose spectral weight exceeds
the maximum spectral weight possible in the Rice-Mele model RM82 .

###  C.1 The polarization distribution

Let us start with the macroscopic polarization of a solid
souza2000polarization ; zak1989berry ; king1993theory ; resta1994macroscopic ;
resta2007theory ; PhysRevB.49.14202  and Kohn’s theory of the insulating state
kohn1964theory . Consider a solid with $\displaystyle N$ electrons and
$\displaystyle M$ nuclei in any dimension $\displaystyle d$. The macroscopic
polarization operator is

$\displaystyle\displaystyle\widehat{\bm{P}}=\frac{1}{V}\left(e\widehat{\bm{X}}+q_{\text{nuc}}\widehat{\bm{X}}_{\text{nuc}}\right),$
(44)

where $\displaystyle\widehat{\bm{X}}=\sum_{i=1}^{N}\hat{\bm{x}}_{i}$ is the
center-of-mass position of the electrons (resp.
$\displaystyle\widehat{X}_{\text{nuc}}$ of the nuclei). Within the clamped
nuclei approximation, the nuclei do not move and we can set
$\displaystyle\widehat{X}_{\text{nuc}}=0$ by a choice of coordinates. The
macroscopic polarization is the expectation
$\displaystyle\langle\widehat{\bm{P}}\rangle$.

The GICs are a systematic way to extract further information from
$\displaystyle\widehat{\bm{P}}$ that we now exploit. Computing
$\displaystyle\langle\widehat{\bm{P}}\rangle=\langle\Psi|\widehat{\bm{P}}|\Psi\rangle$,
averages over many different centers of charge. Following Souza et al
souza2000polarization  we define the distribution of those centers of charge,
i.e. the spatial distribution of polarization, via

$p(\bm{X})=\Braket{\Psi}{\delta(\widehat{\bm{X}}-\bm{X})}{\Psi}.$ (45)

This should be interpreted as the probability that the center of charge is
exactly at the position $\displaystyle\bm{X}$. Different electronic properties
of the solid have already been directly linked to $\displaystyle p(\bm{X})$,
including the Berry connection, localization length of Wannier functions, and
the f-sum rulesouza2000polarization .We can now add the second-harmonic to
that list. As a cautionary note, the polarization is distinct from other real-
space distributions in solids; $\displaystyle p(\bm{X})$ is an different
object from the electron density, or the Wannier functions. The subtle
relationship between these is discussed carefully by Souza et. al.
souza2000polarization .

The great utility of this distribution is tempered by the indirect way it must
be computed, a complication due to the fact that the position operator
$\displaystyle\widehat{\bm{X}}$ is ill-defined resta1998quantum . This implies
that, even knowing the exact eigenvectors of the system, $\displaystyle p(X)$
cannot be straightforwardly computed and, moreover, it is not immediately
clear if $\displaystyle p(X)$ is even a well-defined distribution. A key
result of Souza et al is that $\displaystyle p(X)$ can be carefully defined,
rendering its cumulants computable.

###  C.2 Definition of gauge-invariant cumulants

Recall that the first several cumulants $\displaystyle C_{n}$ of a
distribution $\displaystyle p(X)$ can be written in terms of the moments as

$\displaystyle\displaystyle C_{1}^{i}$ $\displaystyle\displaystyle=\langle
X^{i}\rangle$ (46) $\displaystyle\displaystyle C_{2}^{ij}$
$\displaystyle\displaystyle=\langle X^{i}X^{j}\rangle-\langle
X^{i}\rangle\langle X^{j}\rangle$ (47) $\displaystyle\displaystyle
C_{3}^{ijk}$ $\displaystyle\displaystyle=\langle
X^{i}X^{j}X^{k}\rangle+2\langle X^{i}\rangle\langle X^{j}\rangle\langle
X^{k}\rangle$ (48) $\displaystyle\displaystyle-\langle
X^{i}X^{j}\rangle\langle X^{k}\rangle-\langle X^{i}X^{k}\rangle\langle
X^{j}\rangle-\langle X^{j}X^{k}\rangle\langle X^{i}\rangle$

where $\displaystyle i,j,k\in\left\\{x,y,z\right\\}$ are spatial indices. The
first cumulant is the same as the first moment, the second is the variance,
and the third cumulant is often called the skew of the distribution. A
schematic of the first three cumulants is shown in figure (8).

Figure 8: A schematic of the first three cumulants of the polarization
distribution. In each panel, one cumulant is varied while the others are held
fixed. Note that the true polarization distribution is periodic in
$\displaystyle L$.

A convenient way to compute the moments is in terms of the characteristic
function

$\mathcal{C}(\bm{\alpha})=\langle e^{-i\bm{\alpha}\cdot\bm{X}}\rangle,$ (49)

so that $\displaystyle\langle X^{i}X^{j}\cdots
X^{k}\rangle=i\partial_{\alpha_{i}}i\partial_{\alpha_{j}}\cdots
i\partial_{\alpha_{k}}C(\bm{\alpha})\big{|}_{\bm{\alpha}=0}$. The cumulants
are likewise obtained by differentiating the log of the characteristic
function

$C_{n}^{ij\dots k}=i\partial_{\alpha_{i}}i\partial_{\alpha_{j}}\cdots
i\partial_{\alpha_{k}}\ln\mathcal{C}(\bm{\alpha})\big{|}_{\bm{\alpha}=0}.$
(50)

Following souza2000polarization , we define

$\ln\mathcal{C}(\bm{\alpha})=\frac{V}{\left(2\pi\right)^{d}}\int[d\bm{k}]\;\ln\langle\Psi_{\bm{k}}|e^{-i\bm{\alpha}\cdot\widehat{\bm{X}}}|\Psi_{\bm{k}+\bm{\alpha}}\rangle$
(51)

where $\displaystyle V$ is the volume of the system,
$\displaystyle\int[d\bm{k}]$ denotes the normalized integral over the
Brillouin zone, and $\displaystyle|\Psi_{k}\rangle$ is the many-body
wavefunction with boundary conditions twisted by $\displaystyle e^{ikL}$. In a
single-particle description,
$\displaystyle|\Psi_{k}\rangle=\prod_{n}u_{\bm{k}n}c^{\dagger}_{\bm{k}n}|0\rangle$
where the Bloch wavefunctions
$\displaystyle\psi_{\bm{k}n}(\bm{r})=u_{\bm{k}n}(\bm{r})e^{i\bm{k}\cdot\bm{r}}$
are chosen to satisfy a smooth, periodic gauge:
$\displaystyle\psi_{\bm{k}n}\equiv\psi_{\bm{k}+\bm{G},n}$ for reciprocal
lattice vectors $\displaystyle\bm{G}$. The cumulants of the polarization
distribution can be computed via differentiating (51). The first three are

$\displaystyle\displaystyle C_{1}^{i}\ $ $\displaystyle\displaystyle=\
iV\int[d\bm{k}]\;\operatorname{Tr}\left[c_{1}^{i}\right]$ (52)
$\displaystyle\displaystyle C_{2}^{ij}\ $ $\displaystyle\displaystyle=\
i^{2}V\int[d\bm{k}]\;\operatorname{Tr}\left[c^{ij}_{2}-c_{1}^{i}c_{1}^{j}\right]$
(53) $\displaystyle\displaystyle C_{3}^{ijk}\ $ $\displaystyle\displaystyle=\
i^{3}V\int[d\bm{k}]\;\operatorname{Tr}\left[c_{3}^{ijk}-c_{2}^{ij}c_{1}^{k}-c_{2}^{jk}c_{1}^{i}-c_{2}^{ki}c_{1}^{j}+2c_{1}^{i}c_{1}^{j}c_{1}^{k}\right]$
(54)

where the trace is over occupied bands, and the $\displaystyle c$’s are
matrices in the band space whose matrix elements are

$\displaystyle\displaystyle c_{1}^{i}$ $\displaystyle\displaystyle=i\langle
u_{kn}|\partial_{k_{i}}u_{km}\rangle,$ (55) $\displaystyle\displaystyle
c_{2}^{ij}$ $\displaystyle\displaystyle=i^{2}\langle
u_{kn}|\partial_{k_{i}}\partial_{k_{j}}u_{km}\rangle,$ (56)
$\displaystyle\displaystyle c_{3}^{ijk}$
$\displaystyle\displaystyle=i^{3}\langle
u_{kn}|\partial_{k_{i}}\partial_{k_{j}}\partial_{k_{j}}u_{km}\rangle.$ (57)

One can show that the cumulants are gauge-invariant— hence the name — and
purely real in the presence of time-reversal symmetry.

To compute the cumulants numerically, it is advantageous to use an alternative
formulation (written here for a single occupied band in 1D)hetenyi_cumulants

$\Delta k\prod_{i=0}^{M-1}\ln\langle
u_{k_{i}}|u_{k_{i+1}}\rangle=\sum_{n=1}\frac{(i\Delta k)^{n}}{n!}C_{n}$ (58)

where $\displaystyle k_{i+1}-k_{i}=\Delta k=2\pi/M$, and again using a
periodic gauge for the Bloch wavefunctions $\displaystyle\psi_{k}$. The
cumulants can be measured by computing the left-hand side on successively
finer meshes of $\displaystyle k$-points and fitting the result to a power
series in $\displaystyle\Delta k$, the mesh size. The first several gauge-
invariant cumulants in the Rice-Mele model have been computed recently,
providing intuitive visualizations of the polarization distribution
hetenyi_cumulants .

###  C.3 Interpretation and physical intuition

Each of the GICs has a dual interpretation: as a measure of the spatial
properties of the electrons, or as an electromagnetic response. The first
cumulant tells us the mean of the center of charge polarization, i.e. the
total deviation of electrons from their companion nuclei resta1994macroscopic
. Via equation (44), the polarization of the solid has a straightforward
formula in terms of the first cumulant

$\langle\widehat{\bm{P}}\rangle=\frac{e}{V}\bm{C}_{1}=e\sum_{n\text{
occ}}\int[d\bm{k}]\bm{A}_{\bm{k}n},$ (59)

where $\displaystyle\bm{A}_{\bm{k}n}$ is the Berry connection. The last
expression is familiar from the modern theory of polarization. Equation (59)
shows that, up to dimensionful prefactors, the first cumulant is the
electrical response of the system in the absence of an external field, a
purely static quantity.

The second cumulant measures the covariance of the polarization distribution.
In an insulator, electrons are exponentially localized. Their localization is
in the $\displaystyle i\;(=x,y,z)$ direction is given by (no
summation)souza2000polarization

$\xi_{i}=\lim_{N\to\infty}\sqrt{\frac{1}{N}C_{2}^{ii}}.$ (60)

The second cumulant then gives a great deal of information about the character
of the material. In a metal-insulator transition, for instance, the
localization length diverges, and hence so does the second cumulant.
Physically, the electrons in a metal are almost entirely delocalized, so the
polarization distribution should be nearly flat, with variance on the order of
the size of the system.

It is useful to contrast the information in $\displaystyle C_{2}$ with the
Wannier functions. In 1d, the maximally localized Wannier functions have the
property that the their centers are governed by the centers of polarization,
and their variance is proportional to the squared localization length. We
stress, however, that this does not imply the Wannier density is the same as
the polarization distribution. In dimensions greater than one, there are no
unique maximally-localized single-particle Wannier functions, and all Wannier
functions have variance strictly larger than the squared localization length.

As an optical response, the second cumulant encodes the all-frequency linear
response of the system. A standard application of the flucutation-dissipation
theorem shows that the fluctuations in the polarization distribution
($\displaystyle C_{2}$) is related to the total current response of the
system:souza2000polarization

$\frac{\pi
e^{2}}{V^{2}\hbar}C_{2}^{\alpha\beta}=\frac{d\omega}{\omega}\operatorname{Re}\sigma^{\alpha\beta}(\omega)$
(61)

where $\displaystyle\sigma^{\alpha\beta}$ is the linear conductivity.

###  C.4 A New Sum Rule

Now that we have seen the dual nature — spatial and optical — of the first two
GICs, it is perhaps not too surprising that the third cumulant can give a non-
linear sum rule. The third cumulant measures the skewness of the polarization
distribution which, in one dimension, says if the left or right “shoulder” of
the distribution is larger.

On the optical side, we find that the second cumulant is related to the
second-harmonic response of the system. Given that the first cumulant
determines the polarization and the second cumulant gives a linear sum rule,
it is perhaps not too surprising that the third cumulant can give a non-linear
sum rule. More precisely, we show that for a generic two-band model (see
equation (5) of the main text)

$\Sigma^{\text{shift}}_{a}=\frac{2\pi e^{3}}{V\hbar^{2}}C_{3}^{a}$ (62)

where

$\Sigma_{a}^{\text{shift}}=\int_{0}^{\infty}d\omega\;\operatorname{Re}\sigma^{aaa}(0;-\omega,\omega)$
(63)

and, specializing to the case of two bands,

$C_{3}^{a}=-V\int\frac{d^{d}\bm{k}}{(2\pi)^{d}}\operatorname{Im}\left[c_{3}-3c_{1}c_{2}+2c_{1}^{3}\right]$
(64)

where $\displaystyle c_{n}=\langle
0|\left(i\partial_{k_{a}}\right)^{n}|0\rangle$ for the valence band Bloch
wavefunction $\displaystyle|0\rangle=|u_{0}(\bm{k})\rangle$.

Integrating (25), the spectral weight of the shift current in a two-band model
is

$\displaystyle\displaystyle\Sigma^{\text{shift}}_{a}$
$\displaystyle\displaystyle=\frac{2\pi
e^{3}}{\hbar^{2}}\int[d\bm{k}]\frac{\operatorname{Im}[\langle
0|\partial_{k}H|1\rangle\langle
0|\partial_{k}^{2}H|0\rangle]}{(\epsilon_{1}-\epsilon_{0})^{2}},$ (65)

where $\displaystyle|0\rangle$ and $\displaystyle|1\rangle$ denote valence and
conduction bands, respectively, and
$\displaystyle\varepsilon_{n}=\varepsilon_{n}(k)$ are the band energies.

The Schrödinger equation and its $\displaystyle k$ derivatives give

$\displaystyle\displaystyle H|n\rangle$
$\displaystyle\displaystyle=\epsilon_{n}|n\rangle,$ (66)
$\displaystyle\displaystyle\partial_{k}H|n\rangle+H|\partial_{k}n\rangle$
$\displaystyle\displaystyle=\partial_{k}\epsilon_{n}|n\rangle+\epsilon_{n}|\partial_{k}n\rangle,$
(67)
$\displaystyle\displaystyle\partial_{k}^{2}H|n\rangle+\partial_{k}H|\partial_{k}n\rangle+H|\partial_{k}^{2}n\rangle$
$\displaystyle\displaystyle=\partial_{k}^{2}\epsilon_{n}|n\rangle+\partial_{k}\epsilon_{n}|\partial_{k}n\rangle+\epsilon_{n}|\partial_{k}^{2}n\rangle.$
(68)

Taking inner products and applying $\displaystyle|0\rangle\langle
0|+|1\rangle\langle 1|=I$ implies

$\displaystyle\displaystyle\langle 0|\partial_{k}H|1\rangle$
$\displaystyle\displaystyle=(\epsilon_{1}-\epsilon_{0})\langle
0|\partial_{k}1\rangle,$ (69) $\displaystyle\displaystyle\langle
1|\partial_{k}^{2}H|0\rangle$
$\displaystyle\displaystyle=(\epsilon_{0}-\epsilon_{1})[\langle
1|\partial_{k}^{2}0\rangle-2\langle 1|\partial_{k}0\rangle\langle
0|\partial_{k}0\rangle]+2(\partial_{k}\epsilon_{0}-\partial_{k}\epsilon_{1})\langle
1|\partial_{k}0\rangle.$

Substituting these into the integral in the spectral weight yields

$\displaystyle\frac{1}{(\epsilon_{1}-\epsilon_{0})^{2}}\operatorname{Im}[\langle
0|\partial_{k}H|1\rangle\langle 0|\partial_{k}^{2}H|0\rangle]$ (70)
$\displaystyle=\operatorname{Im}\left\\{\langle
0|\partial_{k}1\rangle\left[-\langle 1|\partial_{k}^{2}0\rangle+2\langle
1|\partial_{k}0\rangle\langle
0|\partial_{k}0\rangle]\\}-\frac{\partial_{k}\epsilon_{0}-\partial_{k}\epsilon_{1}}{\epsilon_{0}-\epsilon_{1}}\operatorname{Im}[\langle
0|\partial_{k}1\rangle\langle 1|\partial_{k}0\rangle\right]\right\\}.$

We can drop the last term since $\displaystyle\langle
0|\partial_{k}1\rangle\langle 1|\partial_{k}0\rangle$ is real. Applying the
resolution of the identity, the first term is

$\displaystyle\langle 0|\partial_{k}1\rangle\langle
1|\partial_{k}^{2}0\rangle$
$\displaystyle=-\langle\partial_{k}0|\partial_{k}^{2}0\rangle+\langle\partial_{k}0|0\rangle\langle
0|\partial_{k}^{2}0\rangle$ (71) $\displaystyle=-\partial_{k}(\langle
0|\partial_{k}^{2}0\rangle)+\langle 0|\partial_{k}^{3}0\rangle-\langle
0|\partial_{k}0\rangle\langle 0|\partial_{k}^{2}0\rangle$
$\displaystyle=\partial_{k}c_{2}+ic_{3}-ic_{1}c_{2},$

while the second term becomes

$\displaystyle\langle 0|\partial_{k}1\rangle\langle
1|\partial_{k}0\rangle\langle 0|\partial_{k}0\rangle$
$\displaystyle=-\langle\partial_{k}0|\partial_{k}0\rangle\langle
0|\partial_{k}0\rangle+\langle\partial_{k}0|0\rangle\langle
0|\partial_{k}0\rangle\langle 0|\partial_{k}0\rangle$ (72)
$\displaystyle=-\partial_{k}(\langle 0|\partial_{k}0\rangle)\langle
0|\partial_{k}0\rangle+\langle 0|\partial_{k}^{2}0\rangle\langle
0|\partial_{k}0\rangle-c_{1}^{3}$
$\displaystyle=+\frac{1}{2}\partial_{k}(c_{1}^{2})+ic_{1}c_{2}-ic_{1}^{3}.$

The total derivatives vanish after integration, and we obtain equation (62).

This sum rule leads to intriguing conclusions. For two-band models, non-linear
optical responses can be understood as a facet of the spatial distribution of
polarization. This provides a clear physical picture that may be more
intuitive than the normal expressions for SHG, which involve k-space sums over
virtual transitions. Moreover, the sum rule can predict the shift current —
and hence the SHG response — of a material as a ground state property.

###  C.5 Upper bounds in the Rice-Mele model and beyond

In light of this relation between the SHG and $\displaystyle C_{3}$, it is
worth revisiting the above fact, equation (43), that there is a maximum SHG in
Rice-Mele models. While this maximum is a good benchmark as to how much
second-harmonic can be produced by a system, it not an absolute bound. By
designing a Hamiltonian with a large $\displaystyle C_{3}$, we will see this
bound may be exceeded.

Figure 9: (Left) Reconstructed polarization distributions in the extended
Rice-Mele model via the maximum entropy methodhetenyi_cumulants . (Right) The
gauge-invariant cumulants and spectral weight in the next-nearest neighbor
extension of the Rice-Mele model. The maximum in the normal Rice-Mele model,
equation (43), is the black line, and the spectral weight as a function of
$\displaystyle\gamma$ is computed via equation 65. All parameters are given in
the text.

For concreteness, consider a generalization of the Rice-Mele model with a
next-nearest neighbor hopping $\displaystyle H=H_{\text{RM}}+H_{\text{NNN}}$
where

$\displaystyle\displaystyle H_{\text{RM}}\ $ $\displaystyle\displaystyle=\
\sum_{n}\Delta(-1)^{n}c_{n}^{\dagger}c_{n}+\left(\frac{t}{2}+(-1)^{n}\frac{\delta}{2}\right)c^{\dagger}_{n}c_{n+1}+\text{h.c.}$
(73) $\displaystyle\displaystyle H_{\text{NNN}}\ $
$\displaystyle\displaystyle=\
\sum_{n}\left(\frac{t^{\prime}}{2}+\frac{\delta^{\prime}}{2}\right)c^{\dagger}_{An}c_{A,n+1}$
(74) $\displaystyle\displaystyle\hskip
30.00005pt+\left(\frac{t^{\prime}}{2}-\frac{\delta^{\prime}}{2}\right)c^{\dagger}_{Bn}c_{B,n+1}+\text{h.c.}$
(75)

In $\displaystyle k$-space this becomes

$H=\sum_{k}\bm{c}_{k}^{\dagger}\left[h_{\text{RM}})+h_{\text{NNN}}\right]\bm{c}_{k}$
(76)

where

$\displaystyle\displaystyle h_{\text{RM}}(\bm{k})$
$\displaystyle\displaystyle=t\cos\left(\frac{kc}{2}\right)\sigma_{x}+\delta\sin\left(\frac{kc}{2}\right)\sigma_{y}+\Delta\sigma_{z}$
(77) $\displaystyle\displaystyle h_{\text{NNN}}(\bm{k})$
$\displaystyle\displaystyle=t^{\prime}\cos(kc)\operatorname{Id}_{2}+\delta^{\prime}\cos(ka)\sigma_{z}.$
(78)

By applying the maximum entropy method to approximately solve the moment
problem hetenyi_cumulants , we can reconstruct the polarization distribution.
Figure 9 shows the reconstructed polarization distributions in the next-
nearest neighbor model with parameters $\displaystyle\Delta=1$,
$\displaystyle\delta=\sqrt{2}$, $\displaystyle t=\sqrt{2}-2\gamma$,
$\displaystyle\delta^{\prime}=\gamma$ and $\displaystyle t^{\prime}=0$. As
$\displaystyle\gamma$ is tuned past zero, $\displaystyle C_{3}$ increases, and
exceeds the maximum possible in the Rice-Mele model. We can thus tune the
model to achieve an arbitrary large $\displaystyle C_{3}$ and, at a metal-
insulator transition, cause it to diverge. Finding materials realizing a large
third cumulant would be excellent candidates for giant shift current or second
harmonic generation.

## References

  * (1) Shen, Y.  Surface second harmonic generation: a new technique for surface studies.  _Annual Review of Materials Science_ 16, 69–86 (1986).
  * (2) Shen, Y.  Surface properties probed by second-harmonic and sum-frequency generation.  _Nature_ 337, 519 (1989).
  * (3) Dähn, A., Hübner, W. & Bennemann, K.  Symmetry analysis of the nonlinear optical response: Second harmonic generation at surfaces of antiferromagnets.  _Physical Review Letters_ 77, 3929 (1996).
  * (4) Petersen, J. C. _et al._ Nonlinear optical signatures of the tensor order in Cd2Re2O7.  _Nature Physics_ 2, 605 (2006).
  * (5) Harter, J., Zhao, Z., Yan, J.-Q., Mandrus, D. & Hsieh, D.  A parity-breaking electronic nematic phase transition in the spin-orbit coupled metal Cd2Re2O7.  _Science_ 356, 295–299 (2017).
  * (6) Far-infrared wave generation.  In Robinson, L. (ed.) _Physical Principles of Far-Infrared Radiation_ , vol. 10 of _Methods in Experimental Physics_ , 10 – 81 (Academic Press, 1973).
  * (7) Chang, Z., Rundquist, A., Wang, H., Murnane, M. M. & Kapteyn, H. C.  Generation of coherent soft x rays at 2.7 nm using high harmonics.  _Physical Review Letters_ 79, 2967 (1997).
  * (8) Ma, Q. _et al._ Direct optical detection of Weyl fermion chirality in a topological semimetal.  _Nature Physics_ 13, 842 (2017).
  * (9) Sun, K. _et al._ Circular photogalvanic effect in the Weyl semimetal TaAs.  _Chinese Physics Letters_ 34, 117203 (2017).
  * (10) Osterhoudt, G. B. _et al._ Colossal photovoltaic effect driven by the singular Berry curvature in a Weyl semimetal.  _arXiv preprint arXiv:1712.04951_ (2017).
  * (11) Ji, Z. _et al._ Spatially dispersive circular photogalvanic effect in a Weyl semimetal.  _arXiv preprint arXiv:1802.04387_ (2018).
  * (12) Wu, L. _et al._ Giant anisotropic nonlinear optical response in transition metal monopnictide Weyl semimetals.  _Nature Physics_ 13, 350 (2017).
  * (13) Li, Z. _et al._ Second harmonic generation in the Weyl semimetal TaAs from a quantum kinetic equation.  _Physical Review B_ 97, 085201 (2018).
  * (14) Penetration depth calculated using refractive index, as measured by ellipsometry, obtained in private communication from Bing Xu, Institute of Phyiscs, Beijing.
  * (15) Weng, H., Fang, C., Fang, Z., Bernevig, B. A. & Dai, X.  Weyl semimetal phase in noncentrosymmetric transition-metal monophosphides.  _Physical Review X_ 5, 011029 (2015).
  * (16) Li, Z. _et al._ Weyl semimetal TaAs: Crystal growth, morphology, and thermodynamics.  _Crystal Growth & Design_ 16, 1172–1175 (2016).
  * (17) See appendices for additional details.
  * (18) Bergfeld, S. & Daum, W.  Second-harmonic generation in gaas: Experiment versus theoretical predictions of $\displaystyle\chi_{xyz}^{(2)}$.  _Physical Review Letters_ 90, 036801 (2003).
  * (19) Rice, M. J. & Mele, E. J.  Elementary excitations of a linearly conjugated diatomic polymer.  _Physical Review Letters_ 49, 1455–1459 (1982).
  * (20) Rangel, T. _et al._ Large bulk photovoltaic effect and spontaneous polarization of single-layer monochalcogenides.  _Physical Review Letters_ 119, 067402 (2017).
  * (21) Morimoto, T. & Nagaosa, N.  Topological nature of nonlinear optical effects in solids.  _Science Advances_ 2, e1501524 (2016).
  * (22) Xu, S.-Y. _et al._ Discovery of a Weyl fermion semimetal and topological Fermi arcs.  _Science_ 349, 613–617 (2015).
  * (23) Lv, B. Q. _et al._ Experimental discovery of Weyl semimetal TaAs.  _Physical Review X_ 5, 031013 (2015).
  * (24) Xu, B. _et al._ Optical spectroscopy of the Weyl semimetal TaAs.  _Physical Review B_ 93, 121110 (2016).
  * (25) Resta, R. & Vanderbilt, D.  Theory of polarization: a modern approach.  In _Physics of Ferroelectrics_ , 31–68 (Springer, 2007).
  * (26) Tan, L. Z. & Rappe, A. M.  Upper limit on nonlinear optical processes: shift current and second harmonic generation in extended systems.  _arXiv preprint arXiv:1708.05433_ (2017).
  * (27) Sipe, J. E. & Shkrebtii, A. I.  Second-order optical response in semiconductors.  _Physical Review B_ 61, 5337–5352 (2000).
  * (28) Moore, J. E. & Orenstein, J.  Confinement-induced Berry phase and helicity-dependent photocurrents.  _Physical Review Letters_ 105, 026805 (2010).
  * (29) Yahyavi, M. & Hetényi, B.  Reconstruction of the polarization distribution of the Rice-Mele model.  _Physical Review A_ 95, 062104 (2017).
  * (30) Souza, I., Wilkens, T. & Martin, R. M.  Polarization and localization in insulators: Generating function approach.  _Physical Review B_ 62, 1666 (2000).
  * (31) Quesnel, B. & Mora, P.  Theory and simulation of the interaction of ultraintense laser pulses with electrons in vacuum.  _Physical Review E_ 58, 3719 (1998).
  * (32) Lab2 add on to LabVIEW® software (1999).  URL http://www.lab2.de.
  * (33) Bloembergen, N. & Pershan, P.  Light waves at the boundary of nonlinear media.  _Physical Review_ 128, 606 (1962).
  * (34) Yang, X., Burch, K. & Ran, Y.  Divergent bulk photovoltaic effect in weyl semimetals.  _arXiv preprint arXiv:1712.09363_ (2017).
  * (35) Zak, J.  Berry’s phase for energy bands in solids.  _Physical Review Letters_ 62, 2747 (1989).
  * (36) King-Smith, R. & Vanderbilt, D.  Theory of polarization of crystalline solids.  _Physical Review B_ 47, 1651 (1993).
  * (37) Resta, R.  Macroscopic polarization in crystalline dielectrics: the geometric phase approach.  _Reviews of modern physics_ 66, 899 (1994).
  * (38) Ortiz, G. & Martin, R. M.  Macroscopic polarization as a geometric quantum phase: Many-body formulation.  _Physical Review B_ 49, 14202–14210 (1994).
  * (39) Kohn, W.  Theory of the insulating state.  _Physical Review_ 133, A171 (1964).
  * (40) Resta, R.  Quantum-mechanical position operator in extended systems.  _Physical Review Letters_ 80, 1800 (1998).