# Nonlinear Quantum Electrodynamics in Dirac materials

Aydın Cem Keser  School of Physics, University of New South Wales, Sydney, NSW
2052, Australia  Australian Research Council Centre of Excellence in Low-
Energy Electronics Technologies,The University of New South Wales, Sydney
2052, Australia  Yuli Lyanda-Geller  Department of Physics and Astronomy,
Purdue University, West Lafayette, IN, 47907, USA  Oleg P. Sushkov  School of
Physics, University of New South Wales, Sydney, NSW 2052, Australia
Australian Research Council Centre of Excellence in Low-Energy Electronics
Technologies,The University of New South Wales, Sydney 2052, Australia

###### Abstract

Classical electromagnetism is linear. However, fields can polarize the vacuum
Dirac sea, causing quantum nonlinear electromagnetic phenomena, e.g.,
scattering and splitting of photons, that occur only in very strong fields
found in neutron stars or heavy ion colliders. We show that strong
nonlinearity arises in Dirac materials at much lower fields $\sim
1\>\text{T}$, allowing us to explore the nonperturbative, extremely high field
limit of quantum electrodynamics in solids. We explain recent experiments in a
unified framework and predict a new class of nonlinear magneto-electric
effects, including a magnetic enhancement of dielectric constant of insulators
and a strong electric modulation of magnetization. We propose experiments and
discuss the applications in novel materials.

Classical electromagnetism is linear and hence supports the principle of
superposition. It has been pointed out by Heisenberg and Euler in 1936 that
due to quantum mechanical effects and the presence of the Dirac sea, linearity
ceases to hold in strong fields Heisenberg and Euler (1936); *euler-
heisenberg. Quantum electrodynamics (QED) is therefore nonlinear as the
electromagnetic field polarizes the Dirac sea as though it is a material
medium. This effect becomes significant at electric and magnetic fields
$E_{\star}\simeq 1.3\times 10^{16}\>\text{V}/\text{cm}$, $B_{\star}\simeq
4.4\times 10^{9}\>\text{T}$, at which the Zeeman splitting and electric
potential over the Compton wavelength become comparable to the electron rest
energy. These are the so called Schwinger critical values Schwinger (1951) and
they are enormous on the laboratory scale. Such fields exist only in exotic
environments e.g. neutron stars Kaspi and Beloborodov (2017);
*reisenegger2001magnetic and heavy ion colliders Tuchin (2013). Nevertheless,
some low-order nonlinear QED effects, such as scattering or splitting of
photons have been observed in the laboratory Marklund and Shukla (2006);
Akhmadaliev  _et al._ (1998), and probing strong field effects is an active
area of research Fedeli  _et al._ (2021).

Dirac materials have been known for decades Cohen and Blount (1960); Lax  _et
al._ (1960); Wolff (1964). Nevertheless, their recently understood topological
properties and surface excitations have led to a surge of interest Qi  _et
al._ (2008); Zhang  _et al._ (2009); Hasan and Kane (2010); Ando (2013);
Wehling  _et al._ (2014); Culcer  _et al._ (2020). The nonlinear
electromagnetic response of Dirac materials have been studied Sodemann and Fu
(2015); Matsyshyn and Sodemann (2019); Culcer  _et al._ (2020); Rostami and
Cappelluti (2020a, b); Bhalla  _et al._ (2020); Du  _et al._ (2020); Morimoto
and Nagaosa (2016) due to their transport properties (e.g. rectification) and
possible applications in photovoltaics. In this letter, rather than transport,
we study dielectric and magnetization response of the three-dimensional (3D)
Dirac insulators and semimetals due to the Dirac vacuum, i.e. filled valence
band. We include nonlinear contributions to all orders by nonperturbatively
analyzing the Heisenberg-Euler action Heisenberg and Euler (2006); Akhiezer
and Berestestskii (1965); Berestetskii  _et al._ (1982); Peskin and Schroeder
(2005); Dunne ; *Dunne_arxiv, going both beyond known results of QED and the
general framework in condensed matter physics 111The high order nonlinear
electromagnetic effects can also emerge in two-dimensional materials
Katsnelson  _et al._ (2013), but is unrelated to the effects we discuss here,
due to dimensionality considerations. For applications to SU(N) gauge fields
see Ref. Redlich (1984).

Figure 1: Diagrammatic representation of the Heisenberg-Euler action $\delta
L_{1}+\delta L_{HE}$ (Eqs. (6) and (Nonlinear Quantum Electrodynamics in Dirac
materials)). The dashed line corresponds to the constant external
electromagnetic fields $B$,$E$ and the solid line is the Green function of a
Dirac sea electron.

In known Dirac materials, we find the typical values of Schwinger fields,
$E_{\star}\sim 10^{5}\>\text{V}/\text{cm}$, $B_{\star}\sim 1\>\text{T}$, are
easily accessible, providing a platform to explore the strong field regime of
QED and to observe quantum nonlinear electromagnetic effects in the
laboratory.

The nonlinear effects contribute to the experimentally observed high-field
magnetization in the recent work on the Weyl semimetal TaAs Zhang  _et al._
(2019) and the Dirac semimetal Bi Iwasa  _et al._ (2019), but the importance
of this observation and its origin in the Heisenberg-Euler effect has not been
recognized. In the present work we demonstrate this connection and show that
the data Zhang  _et al._ (2019); Iwasa  _et al._ (2019) agree with our
predictions. More importantly, we predict a new class of magneto-electric
effects. The most significant is the magnetic field tunable, very large
enhancement of the dielectric constant, reaching up to $\delta\epsilon_{r}\sim
10$ per every $1\>\text{T}$ of the applied magnetic field. We also predict an
electric field modulated magnetization. Both these effects are highly
anisotropic, that is, they depend on relative orientation of E and B-fields
and their crystallographic orientation.

In a material, the classical Lagrangian of the electromagnetic field is 222In
CGS units. In an anisotropic material $\bm{\epsilon},\bm{\mu}$ are symmetric
$3\times 3$ matrices. Some QED textbooks are opposite to the condensed matter
convention, see the footnote in ${\S}$ 129 of Ref. Berestetskii  _et al._
(1982). In this paper $\bm{B}=\bm{\mu}\bm{H}$ is the magnetic flux density
and, $\bm{E}=\bm{\epsilon}^{-1}\bm{D}$ is the screened electric field.

$\displaystyle
L_{cl}=\frac{1}{8\pi}(\bm{E}^{T}\bm{\epsilon}\bm{E}-\bm{B}^{T}\bm{\mu}^{-1}\bm{B}).$
(1)

One of our principal results is the quantum 1-loop, nonlinear, nonperturbative
contribution to the Lagrangian

$\delta L_{HE}\to\frac{\Delta}{24\pi^{2}\lambdabar_{D}^{3}}\left[({\bf
b}\cdot{\bf e})^{2}|{\bf b}|^{-1}+|\mathbf{b}|^{2}\ln|{\bf b}|\right],$ (2)

in a strong B, and weak E-field (See further below and also Sec. S6 of the
supplement). Here, the dimensionless vectors $\mathbf{e}$ and $\mathbf{b}$
depend on the fine structure constant $\alpha_{D}=e^{2}/\hbar v$ of the Dirac
material:

${\bf e}(\alpha_{D})=\frac{\bm{\mathcal{U}E}}{E_{\star}(\alpha_{D})},\quad{\bf
b}(\alpha_{D})=\frac{\bm{\mathcal{U}^{-1}B}}{B_{\star}(\alpha_{D})},$ (3)

and the critical ‘Schwinger’ electric $E_{\star}(\alpha_{D})$ and magnetic
$B_{\star}(\alpha_{D})$ fields in the material are defined by

$E_{\star}^{2}(\alpha_{D})=\frac{v^{2}}{c^{2}}B_{\star}^{2}(\alpha_{D})=\frac{\Delta}{\alpha_{D}\lambdabar_{D}^{3}}=\left(\frac{\Delta^{2}}{e\hbar
v}\right)^{2}.$ (4)

Eqs. (2), (3) and (4) account for the anisotropy of real materials Aronov and
Pikus (1967), for which the velocity tensor is a $3\times 3$ symmetric matrix
Wolff (1964) $\bm{\mathcal{V}}=v\bm{\mathcal{U}}$, with
$\text{det}(\bm{\mathcal{U}})=1$. The term $\bm{\mathcal{U}E}$ and
$\bm{\mathcal{U}^{-1}B}$, are linear transformations of $\bm{E}$ and $\bm{B}$
respectively (See Supplement Sec. LABEL:s-sec:aniso Sup ).

| $2\Delta$ [meV] | ${\alpha_{D}}/{\alpha}$ | $E_{\star}$ [V/cm] | $B_{\star}$ [mT]
---|---|---|---|---
QED $(\Delta=m_{e}c^{2})$ |  $10^{9}$ Thomson (1897); *Millikan; Tiesinga  _et al._ (2021) | 1 | $1.3\times 10^{16}$ | $4.4\times 10^{12}$
$\text{Pb}_{0.5}\text{Sn}_{0.5}\text{Te}$ |  $63$ Dziawa  _et al._ (2012) | 580 | $2.9\times 10^{4}$ | $5.6\times 10^{3}$
$\text{Bi}_{0.9}\text{Sb}_{0.1}$ |  $15.5$ Liu and Allen (1995); Hsieh  _et al._ (2008) | 188 | 571 | 36
TaAs | $0$ | 357 | 0 | 0
Table 1: Comparison of parameters including the band gap ($2\Delta$) edo ,
effective fine structure constant $\alpha_{D}=\frac{e^{2}}{\hbar v}$ as the
ratio $\alpha_{D}/\alpha=c/v$, and the Schwinger fields in Eq. (4).

We have defined the symbols in Eqs.(3),(4) according to convention in QED. The
“Dirac wavelength” $\lambdabar_{D}=\frac{\hbar v}{\Delta}$ and the “Dirac
magneton” $\mu_{D}=\frac{e\hbar v^{2}}{2\Delta c}$ replace the Compton
wavelength and the Bohr magneton respectively. When the fields reach the
‘Schwinger scale’, Zeeman splitting and the potential difference at
$\lambdabar_{D}$ are equal to the half of the Dirac band-gap:

$\displaystyle 2\mu_{D}B_{\star}=\lambdabar_{D}eE_{\star}=\Delta,$ (5)

and the nonlinearity becomes relevant. In Table 1 we list material parameters
considered in this work. For more details see Sec. S2 and Table S2 of the
Supplement Sup .

(a)

(b)

(c)

Figure 2: a) Nonlinear diamagnetic susceptibility $\chi_{ref}-\chi$ versus
magnetic field, $\chi_{ref}=\chi(50\>\text{T})$ in Bi and
$\chi_{ref}=\chi(30\>\text{T})$ in TaAs. In TaAs (blue) the field is along the
crystal c-direction and in Bi there are two directions, binary (red) and
bisectrix (green). The points represent numerical differentiation of TaAs and
Bi magnetization data from Refs. Zhang  _et al._ (2019) and Iwasa  _et al._
(2019) respectively. The experimental points are connected by dashed lines for
guidance. Solid lines represent our theory. The red solid line is our
prediction for $\text{Bi}_{0.9}\text{Sb}_{0.1}$. b) Predicted variation of the
dielectric constant in magnetic field, for parallel (perpendicular) field
configurations shown by solid (dashed) lines. c) Predicted electric field
modulated magnetization as a function of applied magnetic field, along the
magnetic field direction at $E=0.3E_{\star}$.

The quantum contribution to the Lagrangian can be viewed as the sum of the
infinite chain of 1-loop diagrams in Fig.1 that represent the polarization of
the Dirac sea of electrons by external electric and magnetic fields. In this
work we consider only non-magnetic crystals with inversion symmetry Fu and
Kane (2007); *Bansil2016 333TaAs is non-centrosymmetric crystal, but since it
is gapless, we take $E=0$ to avoid transport. TaAs is non-magnetic and the
time reversal symmetry requires that the effective action contains even powers
of B only. and assume the static/quasistatic approximation,
$\omega,kv\ll\Delta$, where $\omega$ and $k$ are the frequency and the wave
number of the external fields. Therefore our diagrams, Fig.1, have only even
numbers of external E-lines and B-lines. Besides diagrams in Fig.1, there are
also multi-loop diagrams suppressed by a factor of $\alpha_{D}/\epsilon\sim
0.03$ per each additional loop, where $\epsilon$ is the large dielectric
constant mainly due to the lattice and intra-ionic polarization. For the
discussion of the suppression of the multi-loop diagrams in the context of
phenomena considered here, see Sec. S3 in the supplement and also Refs. Sham
(1966); Sham and Rice (1966)

In Fig. 1, the first diagram quadratic in external fields is ultraviolet
divergent and is equal to Akhiezer and Berestestskii (1965); edo

$\delta
L_{1}=\frac{\Delta}{12\pi^{2}\lambdabar_{D}^{3}}\ln\left(\frac{\Lambda}{\Delta}\right)\left(|\bf{e}|^{2}-|\bf{b}|^{2}\right).$
(6)

Here the subscript ‘1’ indicates contribution from the first diagram in Fig.1
and $\Lambda\sim v\frac{\hbar\pi}{a}\sim 1\>\text{eV}$ is the ultraviolet cut-
off energy, where $a$ being the lattice spacing. In QED this diagram describes
the electric permittivity and magnetic permeability of vacuum and thus it is
included into the definitions of the electric charge and electromagnetic
fields. As a result, $\delta L_{1}$ does not appear explicitly in QED.
However, for Dirac materials $\delta L_{1}$ is an explicit contribution that
has to be added to the classical Lagrangian Eq. (1). Indeed, this is the
contribution of the Dirac sea (valence band) to the dielectric constant and
magnetic susceptibility.

Equating $(E^{2}-B^{2})/(8\pi)+\delta L_{1}$ to the classical Lagrangian (1),
we find the linear dielectric constant $\bm{\epsilon_{D}}$ and the linear
magnetic susceptibility $\bm{\chi_{D}}$ ($\bm{\mu}=\bm{1}+4\pi\bm{\chi}$):

$\displaystyle\bm{\epsilon_{D}}=\bm{1}+\frac{2\alpha_{D}}{3\pi}\ln\left(\frac{\Lambda}{\Delta}\right)\bm{\mathcal{U}}^{2},\quad\epsilon_{D}\sim
3,$ (7)
$\displaystyle\bm{\chi_{D}}=-\frac{\alpha_{D}}{6\pi^{2}}\frac{v^{2}}{c^{2}}\ln\left(\frac{\Lambda}{\Delta}\right)\bm{\mathcal{U}}^{-2},\quad\chi_{D}\sim-10^{-6}.$
(8)

where estimates are given for the diagonalized tensors.

Eqs. (7) and (8) define the Dirac contributions to the total dielectric and
magnetic susceptibilities. The contribution (7) is relatively small compared
to the total relative permittivity $\bm{\epsilon}$ in Eq. (1), typically
$\epsilon\sim 100$, which is primarily due to the ionic (lattice) and intra-
ionic contributions (See Supplement Table S3). The magnetic response (8)
constitutes a significant part of the diamagnetic susceptibility, which also
has contributions from lower bands and core electrons. For Bismuth, the Dirac
valence band contribution (8) has been previously considered in Ref. Fuseya
_et al._ (2015).

We describe now the nonlinear effects. The diagrams in Fig.1

beyond the first one ($n\geq 2$) are convergent at arbitrarily large
$\mathbf{|e|,|b|}$ 444$\mathbf{|e|,|b|}\ll\Lambda^{2}/\Delta^{2}$, otherwise
the effective Dirac Hamiltonian no longer applies. See Sec. S5 and are re-
summed exactly 555The exact result is obtained by directly summing the
energies of Landau levels Akhiezer and Berestestskii (1965); Berestetskii  _et
al._ (1982), or in modern language, zeta-function regularization Dunne ;
*Dunne_arxiv to yield the 1-loop, nonperturbative Heisenberg-Euler action

$\displaystyle\delta L_{HE}$ $\displaystyle=\sum_{n=2}^{\infty}\delta
L_{n}\equiv\frac{-\Delta}{8\pi^{2}\lambdabar_{D}^{3}}\int_{0}^{\infty}\frac{d\eta
e^{-\eta}}{\eta}$ $\displaystyle\times\bigg{[}A_{-}$ $\displaystyle\cot(\eta
A_{-})A_{+}\cot(\eta
A_{+})-\frac{1}{\eta^{2}}+\frac{1}{3}(A_{-}^{2}+A_{+}^{2})\bigg{]},$
$\displaystyle A_{\mp}$ $\displaystyle=-\frac{i}{2}\left[\sqrt{({\bf b}+i{\bf
e})^{2}}\mp\sqrt{({\bf b}-i{\bf e})^{2}}\right],$ (9)

which accounts for crystal anisotropy, cf., Eq. (3), as well as the strong
field behavior. The imaginary part of Eq.(Nonlinear Quantum Electrodynamics in
Dirac materials), obtained via its analytic continuation, captures the
electric breakdown, which can be avoided in weak electric fields $|{\bf e}|<1$
($E<E_{*}$ ). Then, Eq. (Nonlinear Quantum Electrodynamics in Dirac materials)
can be expanded in powers of $\mathbf{e}$. However, the magnetic field can be
much larger than $B_{\star}$, leading to the asymptotic expression Eq. (2)
(See Supplement Sec. S6). At weak magnetic fields, $|{\bf e}|,|{\bf b}|\ll 1$,
Eq. (Nonlinear Quantum Electrodynamics in Dirac materials) reduces to the 2nd
diagram in Fig.1,

$\displaystyle\delta
L_{2}=\frac{\Delta}{360\pi^{2}\lambdabar_{D}^{3}}\left[\left(|\mathbf{e}|^{2}-|\mathbf{b}|^{2}\right)^{2}+7\left({\bf
e}\cdot{\bf b}\right)^{2}\right].$ (10)

At $E=0$, the nonlinear magnetic susceptibility is

$\displaystyle\bm{\delta\chi}=\frac{\partial^{2}\delta
L_{HE}}{\partial\bm{B}\partial\bm{B}}=\bm{\mathcal{U}}^{-2}\frac{\alpha_{D}}{12\pi^{2}}\frac{v^{2}}{c^{2}}F(|{\bf
b}|);$ $\displaystyle F(|{\bf b}|)=\frac{2}{5}|{\bf b}|^{2},\ \ \|{\bf b}|\ll
1$ $\displaystyle F(|{\bf b}|)=\ln|{\bf b}|,\ \ \ |{\bf b}|\gg 1.$ (11)

The dimensionless function $F(|{\bf b}|)$ in the full range of magnetic fields
obtained by numerical integration of Eq. (Nonlinear Quantum Electrodynamics in
Dirac materials) is shown in Supplement Fig. S2. Strong and weak field limits
of $F$ follow from the actions given by Eqs.(2) and (10) respectively.

The total magnetic susceptibility of the Dirac valence band is the sum of the
linear susceptibility, Eq.(8) and the nonlinear contribution,
$\bm{\chi}=\bm{\chi_{D}}+\bm{\delta\chi}$. When $|{\bf b}|\gg 1$ we have

$\bm{\chi}=-\bm{\mathcal{U}}^{-2}\frac{\alpha_{D}}{12\pi^{2}}\frac{v^{2}}{c^{2}}\ln\left(\frac{c\Lambda^{2}}{e|\bm{\mathcal{U}^{-1}B}|\hbar
v^{2}}\right).$ (12)

Here $\bm{\chi}$ depends on $B$ but not on $\Delta$, and is well-defined in
the limit $\Delta=0$, as in the Weyl semimetal TaAs Zhang  _et al._ (2019).

According to Eqs.(Nonlinear Quantum Electrodynamics in Dirac materials),(12)
the magnetic susceptibility is nonlinear, i.e. it depends on magnetic field.
Remarkably, this Dirac nonlinearity has been recently observed, but its
connection to nonlinear electrodynamics was not identified. Here we show its
origin in the Heisenberg-Euler effect. The magnetization of Weyl semimetal
TaAs has been measured up to $B=30$ T, Ref. Zhang  _et al._ (2019), and
magnetization of Dirac semimetal Bi has been measured up to $B=60$ T, Ref.
Iwasa  _et al._ (2019). In Zhang et al. Zhang  _et al._ (2019) the valence
band contribution to magnetization at $\bm{E}=0$ was considered 666 However,
our approach from the vantage point afforded by the transformation of
Heisenberg-Euler action under dilation and contraction of space (See Sec. S1)
leads to a more accurate overall factor., and in the high magnetic field
limit, the magnetization quasi-linear in the applied B-field is investigated.
Here we study the universal nonlinear susceptibility and eliminate all
uncertainties such as the choice of ultraviolet cut-off $\Lambda$, subleading
terms and contributions from other bands or core electrons.

Both TaAs and Bi have nonzero chemical potential and hence have conduction
electrons. Therefore at weak magnetic fields both compounds show magnetic
oscillations. The conduction electrons freeze and the oscillations disappear
at $B>5$T in Bi Iwasa  _et al._ (2019) and $B>10-13$T in TaAs Zhang  _et al._
(2019). In these ranges of B, we can compare the data with our predictions. In
Fig. 2(a) the points show magnetic susceptibilities of TaAs (c-direction) and
Bi (binary and bisector directions). The points have significant spread as
they are obtained by numerical differentiation of experimental magnetizations
from Refs. Zhang  _et al._ (2019); Iwasa  _et al._ (2019). To focus on the
nonlinearity, we plot $\chi_{ref}-\chi$, where $\chi_{ref}=\chi(B=30T)$ for
TaAs and $\chi_{ref}=\chi(B=50T)$ for Bi. Solid curves present our theoretical
predictions, which is manifestly consistent with the data. For discussion of
material specific details, anisotropy etc., see Supplement Sec. S8.
Interestingly, $\text{Bi}_{0.9}\text{Sb}_{0.1}$ alloy has the band structure
very close to that of Bi, but with no conduction electrons Liu and Allen
(1995); Hsieh  _et al._ (2008), and could be an ideal test platform for our
theory. The susceptibility of this compound has not been measured yet, and the
solid red curve in Fig. 2(a) shows our theoretical prediction.

We now consider novel magneto-electric effects. The nonlinear dielectric
constant is

$\displaystyle\bm{\delta\epsilon_{D}}=4\pi\frac{\partial^{2}\delta
L_{HE}}{\partial\bm{E}\partial\bm{E}}=\bm{\mathcal{U}}^{2}\frac{\alpha_{D}}{3\pi}G_{i}(|{\bf
b}|);$ (13) $\displaystyle|{\bf b}|\ll 1:\ \ G_{||}(|{\bf
b}|)=\frac{1}{3}|{\bf b}|^{2}\ ,\ \ \ \ G_{\perp}(|{\bf
b}|)=-\frac{2}{15}|{\bf b}|^{2}$ $\displaystyle|{\bf b}|\gg 1:\ \ G_{||}(|{\bf
b}|)=|{\bf b}|\ ,\ \ \ \ \ \ \ G_{\perp}(|{\bf b}|)=-\ln(|{\bf b}|).$

Here the index $i=||,\perp$ shows the relative orientation of
$\mathbf{e},\mathbf{b}$ par . Dimensionless functions $G_{i}(|{\bf b}|)$ in
the whole range of $\bf{b}$ obtained by numerical integration of (Nonlinear
Quantum Electrodynamics in Dirac materials) are plotted in Supplement Fig. S2.
Its strong and weak field limits define actions given by Eqs.(2) and (10)
respectively. The dependence of the dielectric constant on the applied
magnetic field is a novel magneto-electric effect. For
${\mathbf{b}}\parallel{\mathbf{e}}$ the contribution $\delta\epsilon_{D}$ is
positive and can be very large, while for ${\mathbf{b}}\perp{\mathbf{e}}$ the
contribution $\delta\epsilon_{D}$ is negative. The expressions for arbitrary
angle between ${\mathbf{b}},{\mathbf{e}}$, and the relation to the angle
between applied fields $\bm{B}$ and $\bm{E}$, which is generally different due
to properties of the anisotropy transformation are given in Sec. S7.
Furthermore, according to (1), (10) there is a nonlinear contribution
quadratic in the electric field,

$\bm{\delta\epsilon_{D}}(\bm{E})=\bm{\mathcal{U}}^{2}\frac{2\alpha_{D}}{15\pi}|\mathbf{e}|^{2},\quad|\mathbf{b}|=0,$
(14)

which is suppressed by $|\mathbf{e}|^{2}/|\mathbf{b}|$ when $|\mathbf{b}|\gg
1$. Notably, at $|\mathbf{e}|,|\mathbf{b}|\ll 1$, contributions (13) and (14)
add up.

The magnetic field induced variation of the dielectric constant in Eq.(13)
scales as $\delta\epsilon_{D}\propto 1/B_{\star}\propto\Delta^{-2}$. Thus, the
effect is most significant in a small band-gap Dirac insulators. In Fig. 2(b)
we plot our predictions for $\text{Bi}_{0.9}\text{Sb}_{0.1}$. For
$\mathbf{e}||\mathbf{b}$ the effect is enormous, $\delta\epsilon_{D}\sim
10$/Tesla. For $\bm{E}\perp\bm{B}$ the effect is smaller and has the negative
sign. In the same Fig. 2(b) we also plot predictions for $\delta\epsilon_{D}$
in $\text{Pb}_{0.5}\text{Sn}_{0.5}\text{Te}$. This compound has larger gap and
therefore the effect is smaller, but still observable.

One more novel magneto-electric effect is the dependence of magnetization, on
the applied electric field. The electric field dependent magnetization
${\bm{M}^{(\mathbf{e})}}=\frac{\partial\delta L}{\partial{\bm{B}}}$, in units
of “Dirac magnetons” per “Dirac volume”, reads

$\displaystyle
4\pi{\bm{M}^{(\mathbf{e})}}=\frac{\bm{\mathcal{U}}^{-1}\mathbf{b}}{|\mathbf{b}|}\frac{\mu_{D}}{3\pi\lambdabar_{D}^{3}}|{\bf
e}|^{2}D_{i}(|{\bf b}|)$ $\displaystyle|{\bf b}|\ll 1:\ \ \ D_{||}(|{\bf
b}|)=\frac{2}{3}|{\bf b}|\ ,\ \ \ D_{\perp}(|{\bf b}|)=-\frac{4}{15}|{\bf b}|$
$\displaystyle|{\bf b}|\gg 1:\ \ \ D_{||}(|{\bf b}|)=1\ ,\ \ \ \ \ \ \
D_{\perp}(|{\bf b}|)=-\frac{1}{|{\bf b}|}.$

The direction of the magnetization (Nonlinear Quantum Electrodynamics in Dirac
materials) in a Dirac crystal is defined by the vector ${\bf b}$ and depends
on crystal anisotropy as described by Eq. (3). Dimensionless functions
$D_{i}(|{\bf b}|)$ in the whole range of $\bf{b}$ obtained by numerical
integration of (Nonlinear Quantum Electrodynamics in Dirac materials) are
plotted in Fig.S2 in Supplementary material. For
${\mathbf{b}}\parallel{\mathbf{e}}$ the magnetization is large and
paramagnetic, while for ${\mathbf{b}}\perp{\mathbf{e}}$ the magnetization is
diamagnetic par . Magnetization (Nonlinear Quantum Electrodynamics in Dirac
materials) is quadratic in the applied electric field and as a function of
magnetic field, saturates when $|{\bf b}|\gg 1$.

To enhance the magnetization in Eq. (Nonlinear Quantum Electrodynamics in
Dirac materials) one needs the electric field as strong as possible. However,
the field is limited by the dielectric strength, $E_{d}$ of the material,
beyond which dielectric breakdown occurs. The breakdown probability (rate of
Zener tunneling by electric field per unit volume) is obtained from Eq.
(Nonlinear Quantum Electrodynamics in Dirac materials) Berestetskii  _et al._
(1982) and found to be $P\propto|{\mathbf{e}}|^{2}e^{-\pi/|{\mathbf{e}}|}$
(See Sec. S4). The most important here is the exponential dependence, which
universally applies to both the Dirac spectrum and quadratic dispersion. Thus,
one expects that $E_{d}$, is proportional to $E_{\star}$. Taking two band
insulators, diamond ($2\Delta\approx 5.5$eV, $E_{d}\approx 10^{7}$ V/cm), and
silicon ($2\Delta\approx 1.14$ eV, $E_{d}\approx 3\times 10^{5}$ V/cm), as
reference materials, we observe that the dielectric strength scales as
$E_{d}\propto\Delta^{2}$. Therefore $E_{d}$ is a fixed fraction of
$E_{\star}$. Significant $E$-dependent magnetic effects Eq. (Nonlinear Quantum
Electrodynamics in Dirac materials) can then be observed for
$|\mathbf{e}|=0.1\text{-}0.3$ 777Of course, there are other factors like the
dependence of dielectric strength on impurities, the size of the sample, etc..
Furthermore, as usual in solids, setups with huge built-in electric fields in
the insulating regime can be explored Calawa  _et al._ (1960); *Rediker2.

For a fixed ${\bf e}=E/E_{\star}$, the electric field modulated magnetization
in Eq. Nonlinear Quantum Electrodynamics in Dirac materials obeys
$M^{(\mathbf{e})}\propto B_{\star}\propto\Delta^{2}$, so materials with large
gap are preferable, unlike in the dependence of dielectric constant on
magnetic field. In Fig. 2(c) we plot the predicted magnetization for
$\text{Pb}_{0.5}\text{Sn}_{0.5}\text{Te}$ versus magnetic field at
$E=10^{4}\>\text{V/cm}$, which corresponds to ${\bf e}\approx 0.3$. For the
both fields, ${\bf e}$ and ${\bf b}$, parallel to the c-axis, the electric
field driven magnetization is $4\pi M^{(\mathbf{e})}\approx 0.2\>\mu\text{T}$
at $B=1T$. When $\bm{E\perp B}$, the magnetization changes sign, see Fig.
2(c). In the same figure we also plot the magnetization in
$\text{Bi}_{0.9}\text{Sb}_{0.1}$ for ${\bf e}\approx 0.3$. Here the effect is
smaller due to the smaller Dirac gap.

The electric field driven magnetization in $\text{Bi}_{0.9}\text{Sb}_{0.1}$
($4\pi M^{(\bf e)}\sim 10^{-8}$T) and in
$\text{Pb}_{0.5}\text{Sn}_{0.5}\text{Te}$ ($4\pi M^{(\mathbf{e})}\sim 2\times
10^{-7}$T) can be feasibly detected in lock-in experiments, in an applied
electric field having a constant and an AC component (with frequency
$\omega$). The induced magnetization is then characterized by contributions
modulated at frequencies $\omega$ and $2\omega$. Of course, the condition
$\hbar\omega\ll 2\Delta$ is assumed fulfilled. Experiments on observation of
$M^{(\mathbf{e})}$ could also take advantage of the SQUID magnetometry,
sensitive to magnetization as low as $10^{-15}\>\text{T}/\sqrt{\text{Hz}}$
Gramolin  _et al._ (2021), much lower than the predicted values.

In conclusion, Based on the Heisenberg-Euler theory of the physical vacuum we
develop the theory of nonlinear electromagnetic effects in Dirac materials. We
explain the results of two recent experiments on nonlinear contribution to
magnetization of Dirac materials. We predict two novel magneto-electric
effects and discuss possible experiments and materials for their observation.

###### Acknowledgements.

Acknowledgements We thank M. O’Brien, H. Takagi, A. O. Sushkov, V. M. Shabaev,
J. Seidel, A. R. Hamilton, U. Zuelicke and Y. Ashim for useful comments. YLG
was supported by the U.S. Department of Energy, Office of Basic Energy
Sciences, Division of Materials Sciences and Engineering under Award DE-
SC0010544. He also acknowledges the Gordon Godfrey bequest for the support of
his visit to UNSW. ACK and OS acknowledge the support from the Australian
Research Council Centre of Excellence in Future Low Energy Electronics
Technologies (CE170100039).

### References

  * Heisenberg and Euler (1936) W. Heisenberg and H. Euler, Zeitschr. Phys. , 714 (1936).
  * Heisenberg and Euler (2006) W. Heisenberg and H. Euler, “Consequences of dirac theory of the positron,” (2006), arXiv:physics/0605038 [physics.hist-ph] .
  * Schwinger (1951) J. Schwinger, Phys. Rev. 82, 664 (1951).
  * Kaspi and Beloborodov (2017) V. M. Kaspi and A. M. Beloborodov, Annual Review of Astronomy and Astrophysics 55, 261 (2017), https://doi.org/10.1146/annurev-astro-081915-023329 .
  * Reisenegger (2001) A. Reisenegger, “Magnetic fields of neutron stars: an overview,” (2001), arXiv:astro-ph/0103010 [astro-ph] .
  * Tuchin (2013) K. Tuchin, Advances in High Energy Physics 2013, 490495 (2013).
  * Marklund and Shukla (2006) M. Marklund and P. K. Shukla, Rev. Mod. Phys. 78, 591 (2006).
  * Akhmadaliev  _et al._ (1998) S. Z. Akhmadaliev, G. Y. Kezerashvili, S. G. Klimenko, V. M. Malyshev, A. L. Maslennikov, A. M. Milov, A. I. Milstein, N. Y. Muchnoi, A. I. Naumenkov, V. S. Panin, S. V. Peleganchuk, V. G. Popov, G. E. Pospelov, I. Y. Protopopov, L. V. Romanov, A. G. Shamov, D. N. Shatilov, E. A. Simonov, and Y. A. Tikhonov, Phys. Rev. C 58, 2844 (1998).
  * Fedeli  _et al._ (2021) L. Fedeli, A. Sainte-Marie, N. Zaim, M. Thévenet, J. L. Vay, A. Myers, F. Quéré, and H. Vincenti, Phys. Rev. Lett. 127, 114801 (2021).
  * Cohen and Blount (1960) M. H. Cohen and E. I. Blount, The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 5, 115 (1960), https://doi.org/10.1080/14786436008243294 .
  * Lax  _et al._ (1960) B. Lax, J. G. Mavroides, H. J. Zeiger, and R. J. Keyes, Phys. Rev. Lett. 5, 241 (1960).
  * Wolff (1964) P. Wolff, Journal of Physics and Chemistry of Solids 25, 1057 (1964).
  * Qi  _et al._ (2008) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B 78, 195424 (2008).
  * Zhang  _et al._ (2009) H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nature Physics 5, 438 (2009).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Ando (2013) Y. Ando, Journal of the Physical Society of Japan 82, 102001 (2013), https://doi.org/10.7566/JPSJ.82.102001 .
  * Wehling  _et al._ (2014) T. Wehling, A. Black-Schaffer, and A. Balatsky, Advances in Physics 63, 1 (2014), https://doi.org/10.1080/00018732.2014.927109 .
  * Culcer  _et al._ (2020) D. Culcer, A. C. Keser, Y. Li, and G. Tkachov, 2D Materials 7, 022007 (2020).
  * Sodemann and Fu (2015) I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * Matsyshyn and Sodemann (2019) O. Matsyshyn and I. Sodemann, Phys. Rev. Lett. 123, 246602 (2019).
  * Rostami and Cappelluti (2020a) H. Rostami and E. Cappelluti, “Dominant role of two-photon vertex in nonlinear response of dirac materials,” (2020a), arXiv:2007.08282 [cond-mat.mes-hall] .
  * Rostami and Cappelluti (2020b) H. Rostami and E. Cappelluti, “Many-body effects in third harmonic generation of graphene,” (2020b), arXiv:2011.03824 [cond-mat.mes-hall] .
  * Bhalla  _et al._ (2020) P. Bhalla, A. H. MacDonald, and D. Culcer, Phys. Rev. Lett. 124, 087402 (2020).
  * Du  _et al._ (2020) Z. Z. Du, C. M. Wang, H.-P. Sun, H.-Z. Lu, and X. C. Xie, “Quantum theory of the nonlinear hall effect,” (2020), arXiv:2004.09742 [cond-mat.mes-hall] .
  * Morimoto and Nagaosa (2016) T. Morimoto and N. Nagaosa, Phys. Rev. B 93, 125125 (2016).
  * Akhiezer and Berestestskii (1965) A. Akhiezer and V. Berestestskii,  _Quantum electrodynamics_ , Interscience Monographs and Texts in Physics and Astronomy No. v. XI (Interscience Publishers [1965], 1965).
  * Berestetskii  _et al._ (1982) V. Berestetskii, E. Lifshitz, and L. Pitaevskii,  _Quantum Electrodynamics: Volume 4_, Course of theoretical physics (Elsevier Science, 1982).
  * Peskin and Schroeder (2005) M. Peskin and D. Schroeder,  _An Introduction to Quantum Field Theory_, Advanced book program (Levant Books, 2005).
  * (29) G. V. Dunne, “Heisenberg–euler effective lagrangians: Basics and extensions,” in  _From Fields to Strings: Circumnavigating Theoretical Physics_, pp. 445–522.
  * Dunne (2004) G. V. Dunne, “Heisenberg–euler effective lagrangians: Basics and extensions,” (2004), arXiv:arXiv:hep-th/0406216 [hep-th] .
  * Note (1) The high order nonlinear electromagnetic effects can also emerge in two-dimensional materials Katsnelson  _et al._ (2013), but is unrelated to the effects we discuss here, due to dimensionality considerations. For applications to SU(N) gauge fields see Ref. Redlich (1984).
  * Zhang  _et al._ (2019) C.-L. Zhang, C. M. Wang, Z. Yuan, X. Xu, G. Wang, C.-C. Lee, L. Pi, C. Xi, H. Lin, N. Harrison, H.-Z. Lu, J. Zhang, and S. Jia, Nature Communications 10, 1028 (2019).
  * Iwasa  _et al._ (2019) A. Iwasa, A. Kondo, S. Kawachi, K. Akiba, Y. Nakanishi, M. Yoshizawa, M. Tokunaga, and K. Kindo, Scientific Reports 9, 1672 (2019).
  * Note (2) In CGS units. In an anisotropic material $\bm{\epsilon},\bm{\mu}$ are symmetric $3\times 3$ matrices. Some QED textbooks are opposite to the condensed matter convention, see the footnote in ${\S}$ 129 of Ref. Berestetskii  _et al._ (1982). In this paper $\bm{B}=\bm{\mu}\bm{H}$ is the magnetic flux density and, $\bm{E}=\bm{\epsilon}^{-1}\bm{D}$ is the screened electric field.
  * Aronov and Pikus (1967) A. G. Aronov and G. E. Pikus, JETP 24, 188 (1967).
  * (36) See Supplemental Material at [URL will be inserted by publisher] for details of the derivation of anisotropic 1-loop effective action and the calculation of nonlinear susceptibilites and a detailed comparison of our theory in different materials and QED.
  * Thomson (1897) J. J. Thomson, Philosophical Magazine 44, 293 (1897).
  * Millikan. (1913) R. A. Millikan., Phys. Rev. 2, 109 (1913).
  * Tiesinga  _et al._ (2021) E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor, Journal of Physical and Chemical Reference Data 50, 033105 (2021), https://doi.org/10.1063/5.0064853 .
  * Dziawa  _et al._ (2012) P. Dziawa, B. J. Kowalski, K. Dybko, R. Buczko, A. Szczerbakow, M. Szot, E. Łusakowska, T. Balasubramanian, B. M. Wojek, M. H. Berntsen, O. Tjernberg, and T. Story, Nature Materials 11, 1023 (2012).
  * Liu and Allen (1995) Y. Liu and R. E. Allen, Phys. Rev. B 52, 1566 (1995).
  * Hsieh  _et al._ (2008) D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Nature 452, 970 (2008).
  * (43) When the band gap is inverted, as in topological insulators, our analysis applies to $|\Delta|$. The additional boundary term in the effective action of the electromagnetic field, proportional to $\bm{E}\cdot\bm{B}$ in topological insulators Fröhlich (2018) and Weyl semimetals Burkov (2018); Raines and Galitski (2017); Rylands  _et al._ (2021) does not affect the nonlinear response which we investigate in this letter.
  * Fu and Kane (2007) L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007).
  * Bansil  _et al._ (2016) A. Bansil, H. Lin, and T. Das, Rev. Mod. Phys. 88, 021004 (2016).
  * Note (3) TaAs is non-centrosymmetric crystal, but since it is gapless, we take $E=0$ to avoid transport. TaAs is non-magnetic and the time reversal symmetry requires that the effective action contains even powers of B only.
  * Sham (1966) L. J. Sham, Phys. Rev. 150, 720 (1966).
  * Sham and Rice (1966) L. J. Sham and T. M. Rice, Phys. Rev. 144, 708 (1966).
  * Fuseya  _et al._ (2015) Y. Fuseya, M. Ogata, and H. Fukuyama, Journal of the Physical Society of Japan 84, 012001 (2015), https://doi.org/10.7566/JPSJ.84.012001 .
  * Note (4) $\mathbf{|e|,|b|}\ll\Lambda^{2}/\Delta^{2}$, otherwise the effective Dirac Hamiltonian no longer applies. See Sec. S5.
  * Note (5) The exact result is obtained by directly summing the energies of Landau levels Akhiezer and Berestestskii (1965); Berestetskii  _et al._ (1982), or in modern language, zeta-function regularization Dunne ; *Dunne_arxiv.
  * Note (6) However, our approach from the vantage point afforded by the transformation of Heisenberg-Euler action under dilation and contraction of space (See Sec. S1) leads to a more accurate overall factor.
  * (53) The angle between $\mathbf{e}$ and $\mathbf{b}$ in general is different to the one between applied fields $\bm{B}$ and $\bm{E}$. The transformation of the angle is given in Sec. S1. The susceptibilities for arbitrary angle between ${\mathbf{b}},{\mathbf{e}}$, are given in Sec. S6 and S7.
  * Note (7) Of course, there are other factors like the dependence of dielectric strength on impurities, the size of the sample, etc.
  * Calawa  _et al._ (1960) A. R. Calawa, R. H. Rediker, B. Lax, and A. L. McWhorter, Phys. Rev. Lett. 5, 55 (1960).
  * Rediker and Calawa (1961) R. H. Rediker and A. R. Calawa, Journal of Applied Physics 32, 2189 (1961), https://doi.org/10.1063/1.1777040 .
  * Gramolin  _et al._ (2021) A. V. Gramolin, D. Aybas, D. Johnson, J. Adam, and A. O. Sushkov, Nature Physics 17, 79 (2021).
  * Katsnelson  _et al._ (2013) M. Katsnelson, G. Volovik, and M. Zubkov, Annals of Physics 331, 160 (2013).
  * Redlich (1984) A. N. Redlich, Phys. Rev. D 29, 2366 (1984).
  * Fröhlich (2018) J. Fröhlich, Reviews in Mathematical Physics 30, 1840007 (2018), https://doi.org/10.1142/S0129055X1840007X .
  * Burkov (2018) A. Burkov, Annual Review of Condensed Matter Physics 9, 359 (2018), https://doi.org/10.1146/annurev-conmatphys-033117-054129 .
  * Raines and Galitski (2017) Z. M. Raines and V. M. Galitski, Phys. Rev. B 96, 161115(R) (2017).
  * Rylands  _et al._ (2021) C. Rylands, A. Parhizkar, A. A. Burkov, and V. Galitski, Phys. Rev. Lett. 126, 185303 (2021).
  * Birrell  _et al._ (1984) N. Birrell, N. Birrell, and P. Davies,  _Quantum Fields in Curved Space_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 1984).
  * (65) L. D. Landau, A. A. Abrikosov, and I. M. Khalatnikov, Doklady Akad. Nauk S.S.S.R. .
  * Suslov (2009) I. M. Suslov, Journal of Experimental and Theoretical Physics 108, 980 (2009).
  * Gell-Mann and Brueckner (1957) M. Gell-Mann and K. A. Brueckner, Phys. Rev. 106, 364 (1957).
  * Giuliani and Vignale (2005) G. Giuliani and G. Vignale,  _Quantum Theory of the Electron Liquid_, Masters Series in Physics and Astronomy (Cambridge University Press, 2005).
  * Lippmann and Schwinger (1950) B. A. Lippmann and J. Schwinger, Phys. Rev. 79, 469 (1950).
  * Salpeter and Bethe (1951) E. E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232 (1951).
  * Wannier (1937) G. H. Wannier, Phys. Rev. 52, 191 (1937).
  * Mott (1938) N. F. Mott, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 167, 384 (1938), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1938.0137 .
  * Schrödinger (1926) E. Schrödinger, Annalen der Physik 384, 361 (1926), https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.19263840404 .
  * Abrikosov (2017) A. Abrikosov,  _Fundamentals of the Theory of Metals_ (Dover Publications, 2017).
  * Dyson (1952) F. J. Dyson, Phys. Rev. 85, 631 (1952).
  * HUET  _et al._ (2012) I. HUET, M. R. DE TRAUBENBERG, and C. SCHUBERT, International Journal of Modern Physics: Conference Series 14, 383 (2012), https://doi.org/10.1142/S2010194512007507 .
  * Lin and Kleinman (1966) P. J. Lin and L. Kleinman, Phys. Rev. 142, 478 (1966).
  * Burstein  _et al._ (1968) E. Burstein, S. Perkowitz, and M. Brodsky, Journal de Physique Colloques 29, C4 (1968).
  * Svane  _et al._ (2010) A. Svane, N. E. Christensen, M. Cardona, A. N. Chantis, M. van Schilfgaarde, and T. Kotani, Phys. Rev. B 81, 245120 (2010).
  * Kanai and Shohno (1963) Y. Kanai and K. Shohno, Japanese Journal of Applied Physics 2, 6 (1963).
  * Suzuki and Adachi (1995) N. Suzuki and S. Adachi, Japanese Journal of Applied Physics 34, 5977 (1995).
  * Herman  _et al._ (1968) F. Herman, R. L. Kortum, I. B. Ortenburger, and J. P. Van Dyke, J. Phys. Colloques C4, 62 (1968).
  * Hayasaka and Fuseya (2016) H. Hayasaka and Y. Fuseya, Journal of Physics: Condensed Matter 28, 31LT01 (2016).
  * Lefebvre  _et al._ (1987) I. Lefebvre, M. Lannoo, G. Allan, A. Ibanez, J. Fourcade, J. C. Jumas, and E. Beaurepaire, Phys. Rev. Lett. 59, 2471 (1987).
  * Spr (a) “Antimony telluride (sb2te3) dielectric constants: Datasheet from landolt-börnstein \- group iii condensed matter · volume 41c: “non-tetrahedrally bonded elements and binary compounds i” in springermaterials (https://doi.org/10.1007/10681727_1045),”  (a), copyright 1998 Springer-Verlag Berlin Heidelberg.
  * Spr (b) “Bismuth telluride (bi2te3) optical properties, dielectric constant: Datasheet from landolt-börnstein - group iii condensed matter · volume 41c: “non-tetrahedrally bonded elements and binary compounds i” in springermaterials (https://doi.org/10.1007/10681727_963),”  (b), copyright 1998 Springer-Verlag Berlin Heidelberg.
  * Spr (c) “Bismuth selenide (bi2se3) optical properties, dielectric constants: Datasheet from landolt-börnstein - group iii condensed matter · volume 41c: “non-tetrahedrally bonded elements and binary compounds i” in springermaterials (https://doi.org/10.1007/10681727_945),”  (c), copyright 1998 Springer-Verlag Berlin Heidelberg.
  * Jain  _et al._ (2013) A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. a. Persson, APL Materials 1, 011002 (2013).
  * Schwerdtfeger and Nagle (2019) P. Schwerdtfeger and J. K. Nagle, Molecular Physics 117, 1200 (2019), https://doi.org/10.1080/00268976.2018.1535143 .
  * Rudolph  _et al._ (1980) R. Rudolph, H. Krüger, B. Fellmuth, and R. Herrmann, physica status solidi (b) 102, 295 (1980), https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.2221020127 .
  * Weng  _et al._ (2015) H. Weng, C. Fang, Z. Fang, B. A. Bernevig, and X. Dai, Phys. Rev. X 5, 011029 (2015).

## Supplemental Material

###  S1 Dirac Cone Anisotropy

Here, for the sake of completeness, we provide a prescription to treat
anisotropic crystals, adopted from Ref. Aronov and Pikus (1967). In real
crystals the velocity is a tensor that is represented by a $3\times 3$ real
symmetric matrix $\bm{\mathcal{V}}$, so that the Dirac Hamiltonian reads

$\bm{H}=\bm{\beta}\Delta-\bm{1}|e|{\phi}+\sum_{i,j}\bm{\alpha}_{i}\mathcal{V}_{ij}(p_{j}+|e|A_{j}/c),\quad
i,j=1,2,3.$ (S16)

where $\bm{\alpha}_{i},\bm{\beta}$ are $4\times 4$ Dirac matrices and $\bm{1}$
is a $4\times 4$ identity matrix. We define the Dirac velocity $v$ and the
anisotropy matrix $\bm{\mathcal{U}}$ as

$v^{3}=|\text{det}(\bm{\mathcal{V}})|,\quad\bm{\mathcal{U}}=\bm{\mathcal{V}}/v.$
(S17)

####  S1.1 Dilation/contractions of space that renders the cone isotropic

The matrix $\bm{\mathcal{U}}$, being a symmetric matrix with
$\text{det}(\bm{\mathcal{U}})=1$, without loss of generality, induces a volume
preserving transformation on coordinates, that is

$\bm{\tilde{x}}=\bm{\mathcal{U}}^{-1}\bm{x},\quad\bm{\tilde{p}}=\bm{\mathcal{U}}\bm{p},\quad\bm{\tilde{A}}=\bm{\mathcal{U}}\bm{A},\quad
d^{3}\tilde{x}=d^{3}x{,}$ (S18)

with which the Dirac Hamiltonian assumes its isotropic form

$\bm{H}=v\bm{{\alpha}}\cdot(\bm{\tilde{p}}+|e|\bm{\tilde{A}}/c)+\bm{\beta}\Delta-\bm{1}|e|{\phi}{.}$
(S19)

The Lagrangian, being a scalar, is invariant under this volume preserving
transformation. Meanwhile, the scaled electric and magnetic fields are

$\bm{\tilde{E}}=\frac{\partial\phi}{\partial\bm{\tilde{x}}}=\bm{\mathcal{U}E},\quad\bm{\tilde{B}}=\bm{\mathcal{U}^{-1}B}{,}$
(S20)

where the equation for $\bm{\tilde{B}}$ follows from the identity (repeated
indices are summed)

$\mathcal{U}_{nk}\tilde{B}_{k}=\mathcal{U}_{nk}\frac{\partial\tilde{A}_{i}}{\partial\tilde{x}_{j}}\epsilon_{ijk}=\frac{\partial
A_{s}}{\partial
x_{m}}\mathcal{U}_{im}\mathcal{U}_{js}\mathcal{U}_{nk}\epsilon_{ijk}=B_{n}{.}$
(S21)

Finally, the physical susceptibilities $\chi,\chi^{e}$ in the anisotropic
crystal as a function of applied fields $\bm{E},\bm{B}$ are

$\displaystyle\bm{\chi}(\bm{E},\bm{B})=\frac{\partial^{2}L}{\partial\bm{B}\partial\bm{B}}=\bm{\mathcal{U}}^{-1}[\bm{\tilde{\chi}}(\bm{\tilde{E}},\bm{\tilde{B}})]\bm{\mathcal{U}}^{-1}$
$\displaystyle\bm{\chi^{e}}(\bm{E},\bm{B})=\frac{\partial^{2}L}{\partial\bm{E}\partial\bm{E}}=\bm{\mathcal{U}}[\bm{\tilde{\chi}^{e}}(\bm{\tilde{E}},\bm{\tilde{B}})]\bm{\mathcal{U}}{.}$
(S22)

The magnetization and polarization vectors and higher order susceptibility
tensors transform in the usual way like above. We re-iterate a point already
mentioned in the main text that even when the crystal anisotropy is not taken
into account, that is $\bm{\mathcal{U}}=\bm{1}$, the susceptibilities are
intrinsically anisotropic, since they depend on the mutual alignment of
electric and magnetic fields. Crystal anisotropy is an additional source of
directional dependence that is relevant in the experimental context.

####  S1.2 Invariance of $E\cdot B$ and the transformation of the angle
between $E$ and $B$

In account of anisotropy, the normalized fields are defined as Eq. (3). We
note that the scalar product of externally applied fields
$\mathbf{E}\cdot\mathbf{B}$ transforms as

$\bm{E}\cdot\bm{B}=\mathcal{U}^{-1}_{ik}\tilde{E}_{k}\mathcal{U}_{il}\tilde{B}_{l}=B_{\star}E_{\star}\mathbf{e}\cdot\mathbf{b},$
(S23)

i.e., is proportional to the scalar product of the transformed fields.
However, we note that $\bm{\mathcal{U}}$ is not a simple rotation. Hence,
while scalar products are proportional, it does not preserve angles or norms
taken separately. We see this when we express
$\mathbf{\hat{e}}\cdot\mathbf{\hat{b}}$ in terms of
$\bm{\hat{E}}\cdot\bm{\hat{B}}$ as

$\mathbf{\hat{e}}\cdot\mathbf{\hat{b}}=\frac{\bf{e}\cdot\bf{b}}{|\bf{e}||\bf{b}|}=\frac{\bm{E}\cdot\bm{B}}{|\bm{\mathcal{U}E}||\bm{\mathcal{U}}^{-1}\bm{B}|}=\frac{|\bm{E}||\bm{B}|}{|\bm{\mathcal{U}E}||\bm{\mathcal{U}}^{-1}\bm{B}|}\bm{\hat{E}}\cdot\bm{\hat{B}}$
(S24)

Therefore

$\bm{E}\perp\bm{B}\iff\mathbf{e}\perp\mathbf{b}.$ (S25)

We can work in a reference frame where $\bm{\mathcal{U}}$ is diagonal.
However, even when $\bm{E}\parallel\bm{B}$ are directed along a general
direction, we have

$\bm{E}\parallel\bm{B}\implies\mathbf{\hat{e}}\cdot{\mathbf{\hat{b}}}=\frac{1}{|\bm{\mathcal{U}\hat{E}}||\bm{\mathcal{U}}^{-1}\bm{\hat{E}}|}\leq
1$ (S26)

The equality is satisfied when $\bm{E}$ and $\bm{B}$ are directed along the
principle directions of the frame that diagonalizes $\bm{\mathcal{U}}$, e.g.
$\bm{\hat{E}}=\bm{\hat{B}}=(1,0,0)^{T}$. However, for general direction, the
vectors $\mathbf{e}$ and $\mathbf{b}$ are no longer parallel, for example, if

$\bm{\mathcal{U}}=\text{diag}(u,u,1/u^{2}),\quad\bm{\hat{E}}=\bm{\hat{B}}=\frac{1}{\sqrt{3}}(1,1,1)^{T}$
(S27)

we obtain

$\mathbf{\hat{e}}\cdot{\mathbf{\hat{b}}}=\left(\frac{9}{{\left(\frac{2}{u^{2}}+u^{4}\right)}\,{\left(2\,u^{2}+\frac{1}{u^{4}}\right)}}\right)^{1/2}$
(S28)

For exemplary values of $u=1.2,1.5,2$, we get the angle between
$\mathbf{\hat{e}}$ and ${\mathbf{\hat{b}}}$ as
$28^{\circ},55^{\circ},75^{\circ}$ respectively.

In terms of transformed vectors $\mathbf{e}$ and $\mathbf{b}$, Dirac materials
are described by a Hamiltonian that has the usual isotropic form which
simplifies calculations. In particular, for the perpendicular fields
$\mathbf{e}$ and $\mathbf{b}$, one can apply the Lorentz transformation that
eliminates the smaller of the fields, resulting in the purely electric field
or the purely magnetic field case. In the case when fields $\mathbf{e}$ and
$\mathbf{b}$ are not orthogonal, a particular choice of a Lorentz
transformation leads to a Hamiltonian with a new set of electric and magnetic
fields parallel to each other, for which the calculations of dielectric and
magnetic susceptibilities are simplified (Section S6). Once the results are
obtained in the frame where the fields are parallel, the inverse Lorentz
transformation allows to express them in terms of the original non-orthogonal
fields $\mathbf{e}$ and $\mathbf{b}$. Finally, inverse transformation to
initial anisotropic system will give the results in terms of applied fields
$\bm{E}$ and $\bm{B}$. In Sec. S7 we derive the susceptibilities for an
arbitrary angle between $\mathbf{e}$ and $\mathbf{b}$, see Table S6 for quick
reference.

Figure S1: Multi-loop and multi-leg diagrams for the effective action. The dashed line corresponds to the external electromagnetic fields $B$,$E$, the blue wavy lines are Coulomb interactions and the solid line is the Green function of a Dirac sea electron. Only the 1-loop, 2-leg diagram given by (S40), is ultraviolet divergent, the others are convergent. The divergence is renormalized into the dielectric constant and relative permeability. The multi-loop diagrams contain Coulomb interaction (blue) that are suppressed by the relative permittivity, hence carry the factor $\alpha_{D}/\epsilon\sim 0.03$. The 1-loop diagrams with more than 4-legs are given by (S45), when the expansion parameters $|\mathbf{e}|^{2}=\frac{\alpha_{D}E^{2}}{\Delta/\lambdabar_{D}^{3}}$ or $|\mathbf{b}|^{2}=\frac{\alpha_{D}B^{2}v^{2}/c^{2}}{\Delta/\lambdabar_{D}^{3}}$, are large, the diagrams can be exactly summed to yield Eq. (S4) Parameter | Dirac Insulator | QED
---|---|---
| symbol | expression | value$\sim$ | symbol | expression | value$\sim$
Energy gap | $2\Delta$ |  |  $100\>\text{meV}$888This estimate is based on tunable gap compound $\text{Pb}_{1-x}\text{Sn}_{x}\text{Se}$ Dziawa  _et al._ (2012). There is no gap in Weyl, Dirac semimetals (see TaAs Zhang  _et al._ (2019)), and appropriate limit has to be taken. Gap is tunable between $\pm 10\>\text{meV}$ in Bi and Sb alloys. Typical gap for materials in Table S3 is $200\>\text{meV}$. | $2m_{e}c^{2}$ |  | $1\>\text{MeV}$
Dirac velocity | $v$ | $c/400$ | $10^{6}\>\text{m}/\text{s}$ | $c$ |  | $3\times 10^{8}\>\text{m/s}$
Wavelength (reduced) | $\lambdabar_{D}$ | $\frac{\hbar v}{\Delta}$ | $10\>\text{nm}$ | $\lambdabar_{C}$ | $\frac{\hbar}{m_{e}c}$ | $3.86\times 10^{-4}\>\text{nm}$
Fine structure constant | $\alpha_{D}$ | $\frac{c\alpha}{v}$ | $3$ | $\alpha$ | $\frac{e^{2}}{\hbar c}$ | $1/137$
Magnetic moment | $\mu_{D}$ | $\frac{e\hbar v^{2}}{2\Delta c}$ | $3.7\>\text{meV/T}$ | $\mu_{B}$ | $\frac{e\hbar}{2m_{e}c}$ | $5.8\times 10^{-2}\text{meV}/\text{T}$
Ultraviolet cut-off | $\Lambda$ |  | $2\>\text{eV}$ | $\Lambda_{QED}$ |  |  $10^{28}\text{-}10^{286}\>\text{eV}$999Low estimate is based on Planck energy Birrell  _et al._ (1984), where gravitational effects come into play. Large estimate is the Landau pole, where the coupling constant diverges. Landau  _et al._ ; *Suslov2009
Dielectric screening (due to Dirac sea) | $\epsilon_{D}$ | $1+\frac{2\alpha_{D}}{3\pi}\log\frac{\Lambda}{\Delta}$ | $3$ | $\epsilon_{\Lambda}$ |  $1+\frac{2\bar{e}^{2}}{3\pi\hbar c}\log\frac{\Lambda_{QED}}{m_{e}c^{2}}$ 101010See footnote g |
Relative permittivity (renormalized) | $\epsilon$ |  $\sim\epsilon_{D}+\epsilon_{i}$111111$\epsilon_{i}$ is the lattice and intra-ionic contribution. |  $100$121212See Table S3. | $\epsilon$ |  |  $1$131313By definition, the vacuum has a permittivity of unity or $\epsilon_{0}$ in SI units.
Electric charge in use | bare: $e$ |  | $1.6\times 10^{-19}\>\text{C}$ | RG: $e_{0}=e$ |  $\bar{e}/\sqrt{\epsilon_{\Lambda}}$141414The quantity $\bar{e}$ is the running charge defined at the ultraviolet cutoff energy. Peskin and Schroeder (2005) | $1.6\times 10^{-19}\>\text{C}$
Field definition in use | screened: $E$ |  $D/\epsilon$151515$D$ is the macroscopic electric displacement that is due to free charges only. |  | RG: $\bar{E}=E$ | $\sqrt{\epsilon_{\Lambda}}E_{0}$ |
Schwinger E-Field | $E_{\star}$ | $\frac{\Delta^{2}}{e\hbar v}$ | $5\times 10^{4}\>\text{V/cm}$ | $E_{S}$ | $\frac{m_{e}^{2}c^{3}}{e\hbar}$ | $1.3\times 10^{16}\>\text{V/cm}$
Schwinger B-Field | $B_{\star}$ | $\frac{\Delta^{2}c}{e\hbar v^{2}}$ | $6.8\>\text{T}$ | $B_{S}$ | $\frac{m_{e}^{2}c^{3}}{e\hbar}$ | $4.4\times 10^{9}\>\text{T}$
1-loop expansion parameter 1 | $|\bf{e}|^{2}$ | $\frac{\alpha_{D}E^{2}}{\Delta/\lambdabar_{D}^{3}}$ |  $<0.1$161616Dielectric break down occurs beyond $|\bf{e}|\sim 0.1$ which is captured by our nonperturbative theory, see Eq. (S41) | $|\bf{e}|^{2}$ | $\frac{\alpha B^{2}}{m_{e}c^{2}/\lambdabar_{c}^{3}}$ |  $<0.1$171717Electric field accelerates free particle-hole pairs beyond this point.
1-loop expansion parameter 2 | $|\bf{b}|^{2}$ | $\frac{\alpha_{D}B^{2}}{\Delta/\lambdabar_{D}^{3}}\frac{v^{2}}{c^{2}}$ |  $<10^{3}$181818First Landau level is ejected beyond the cut-off energy past this point in external magnetic field. Therefore our theory breaks down. Note that perturbation theory is not valid beyond $|b|>1$, but the nonperturbative resummation in Eq. (S4) is still valid. | $|\mathbf{b}|^{2}$ | $\frac{\alpha B^{2}}{m_{e}c^{2}/\lambdabar_{c}^{3}}$ |  $<10^{44}$191919Magnetic field energy curves space-time beyond this point.
N-loop expansion parameter (interaction probability) |  | $\alpha_{D}/\epsilon$ |  $0.03$2020201-loop nonperturbative theory is valid as long as $\alpha_{D}<\epsilon\sim 100$, that is $v\gtrsim 10^{-4}c$. | $\alpha$ |  | $\approx 0.01$
Invariant 1: $|\bf{e}|^{2}-|\bf{b}|^{2}$ |  $\propto\mathcal{F}^{\mu\nu}\mathcal{F}_{\mu\nu}$212121$\mathcal{F}_{\mu\nu}=\tilde{\partial}_{\mu}\mathcal{A}_{\nu}-\tilde{\partial}_{\nu}\mathcal{A}_{\mu}$. This is different from true relativistic field tensor $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. See Eq. (S30). | $E^{2}-\frac{v^{2}}{c^{2}}B^{2}$ |  | $\propto{F}^{\mu\nu}{F}_{\mu\nu}$ | $E^{2}-B^{2}$ |
Invariant 2: $\bf{e}\cdot\bf{b}$ | $\propto\mathcal{F}_{\alpha\beta}\mathcal{F}_{\mu\nu}\varepsilon^{\alpha\beta\mu\nu}$ | $\frac{v}{c}\bm{E}\cdot\bm{B}$ |  | $\propto{F}_{\alpha\beta}{F}_{\mu\nu}\varepsilon^{\alpha\beta\mu\nu}$ | $\bm{E}\cdot\bm{B}$ |
Table S2: Comparison of Dirac insulator with Quantum Electrodynamics (QED).
For simplicity isotropic case is considered. For the general case see Sec. S1

###  S2 Comparison to the relativistic quantum electrodynamics

The Dirac equation that follows from Eq. (S19), after multiplying from the
left by $\bm{\beta}$, is

$i\bm{\gamma}^{0}(\partial_{t}-i|e|{\phi})-v\bm{{\gamma}}\cdot(\bm{\tilde{p}}+|e|\bm{\tilde{A}}/c)\psi-\Delta\psi=0,\quad\bm{\gamma}^{0}=\bm{\beta},\>\bm{\gamma}^{i}=\bm{\beta}\bm{\alpha}^{i}.$
(S29)

If we define the effective and true 4-position and electromagnetic 4-potential
respectively as

effective:
$\displaystyle\quad\tilde{x}^{\mu}=(vt,\bm{\tilde{x}}),\quad\mathcal{A}^{\mu}=\frac{1}{v}(\phi,v\tilde{\bm{A}}/c)$
(S30) true-relativistic:
$\displaystyle\quad{x}^{\mu}=(ct,\bm{{x}}),\quad{A}^{\mu}=\frac{1}{c}(\phi,\bm{A}),\quad\text{where
}\eta^{\mu\nu}=(+,-,-,-),$ (S31)

we can write down the Dirac action coupled to classical electromagnetism in
the (anisotropic) material as

$\text{Dirac material:}\quad S=\int
dtd^{3}\tilde{x}\>\left(\bar{\psi}[iv\gamma^{\mu}(\tilde{\partial}_{\mu}-i|e|\mathcal{A}_{\mu})-\Delta]\psi+\frac{1}{8\pi}(\bm{E}^{T}\bm{\epsilon}\bm{E}-\bm{B}^{T}\bm{\mu}^{-1}\bm{B})\right).$
(S32)

In addition, when the band gap is inverted, as in topological insulators,
there is an additional boundary term bulk-boundary term $\bm{E}\cdot\bm{B}$
Fröhlich (2018) which does not effect the nonlinear response which we
investigate in this letter. For this reason we assume $\Delta>0$ with out loss
of generality.

For comparison, the true relativistic QED action is

$\text{QED:}\quad S=\int
dtd^{3}x\>\left(\bar{\psi}[ic\gamma^{\mu}({\partial}_{\mu}-i|e|{A}_{\mu})-m_{e}c^{2}]\psi+\frac{1}{8\pi}({E}^{2}-{B}^{2})\right).$
(S33)

The first difference one can notice is the speed of the Dirac fermion, that is
$c$ in QED. In the material with a gap $2\Delta$, the Fermi-Dirac velocity is
related to the effective mass $m^{*}$ (measured in terms of electron mass) and
satisfies

$\displaystyle E=\sqrt{\Delta^{2}+p^{2}v^{2}},\quad
E\approx\Delta+\frac{v^{2}p^{2}}{2\Delta}\implies m^{*}m_{e}v^{2}=\Delta,$
$\displaystyle\Delta\approx 0.1\>\text{eV},\>m^{*}=0.01\text{-}0.5\implies
1000\gtrsim c/v\gtrsim 100{.}$ (S34)

For this reason, the effective fine structure constant of the insulator is

$\alpha_{D}=\frac{e^{2}}{\hbar v}\sim\alpha\times 400\approx 3{.}$ (S35)

The minimally coupled 4-potentials are different as outlined in Eq. (S30). The
classical action of the electromagnetic field in the two cases are obviously
different. In the material case, $\bm{\epsilon},\bm{\mu}$ are symmetric
matrices that define the relative permittivity and the relativity of the
material, respectively. In the vacuum they satisfy $\epsilon=\mu=1$ in the CGS
units.

A more subtle difference between QED and the Dirac insulator is the definition
of electric charge and electric field, when seen from a renormalization point
of view. In QED the charge is normalized at zero momentum transfer,
$e_{0}=e_{q=0}$ (referred to as the running electric charge Peskin and
Schroeder (2005) ), which is related to the bare charge $\bar{e}$ through the
ultraviolet (UV) cutoff dependent vacuum dielectric screening,
$e_{0}=\bar{e}/\sqrt{\epsilon_{\Lambda}}$. However, in condensed matter, the
charge is defined independently of scale and is given by $e\approx 1.6\times
10^{-19}\>\text{C}$. Nevertheless, due to RPA dielectric screening, the
apparent charge at the macroscopic scale is $e/\sqrt{\epsilon}$, hence
different from the bare $e$ that applies at the lattice spacing scale
$q\sim\hbar\pi/a$. Definitions of the electric field are also different. In
QED, the field is normalized to $\bm{\bar{E}}=\sqrt{\epsilon_{\Lambda}}\bm{E}$
so that the energy density of the field is fixed as $\bar{E}^{2}/(8\pi)$, and
the dielectric constant of vacuum is by definition equal to unity. However in
condensed matter, the Lagrangian of the free field, is already written in
terms of the screened field $\bm{E}$, as $\epsilon E^{2}/(8\pi)$. The coupling
term is also written in terms of the bare charge and the screened field
$e\bm{E}$, hence $\epsilon$ does not appear in the definition of the effective
fine structure constant $\alpha_{D}=e^{2}/(\hbar v)$. Importantly, the
coupling term is independent of the definition,
$e_{0}{\bar{\bm{E}}}=e{\bm{E}}$, allowing us to use the 1-loop renormalized
effective action of QED in the condensed matter setting. A detailed comparison
of parameters in Dirac insulator to their counterparts in QED is given in
Table S2.

###  S3 Effective action for Dirac insulator

The effective action can be schematically organized in terms of multi-leg and
multi-loop diagrams as in Fig. S1. The multi-loop diagrams contain at least
one internal interaction line and therefore in QED, each such line contains
the factor $\alpha\sim 0.01$.

We now justify the possibility to limit the calculation of the effective
action to the 1-loop approach and justify neglecting multi-loop diagrams in
Dirac materials, particularly insulators. In the Dirac insulator, we first
note that multi-loop diagrams are UV convergent hence the polarization is
dictated mostly by the low energy-momentum sector where the static dielectric
constant $\epsilon$ applies Sham (1966); Sham and Rice (1966). The Coulomb
interaction is therefore screened due to $\epsilon$, which is dominated by the
lattice and intra-ionic polarization contribution. In a wide range of
practical situations as seen in Table S3 $\epsilon\sim 100$ and therefore the
interaction lines are suppressed by the factor
$\frac{\alpha_{D}}{\epsilon}\sim 0.03$.

We note that, no matter how small the interaction parameter
$\alpha_{D}/\epsilon$, the multi-loop diagrams may come with large
combinatorial factors and can render the perturbation series divergent. In
both QED and condensed matter physics several examples of divergences are
known. In some of these divergences, the coupling constant $\alpha$ gets
enhanced by an additional parameter $R$, $\alpha\rightarrow R\alpha$. For
example, in an electron gas in metals, there is Random Phase Approximation
(RPA) infrared divergence, where $\alpha$ is enhanced by a large distance.
This diverging series of diagrams is eventually converging, and can be re-
summed, as was shown by Gell-Mann Brueckner. Gell-Mann and Brueckner (1957)
However, in the case of Dirac insulator, the electron gas infrared divergence
does not arise. Indeed, for example the expansion of contribution to energy
due to interactions,

$E=E_{0}+\langle 0|H_{ee}|0\rangle+\sum_{n}\frac{\langle
0|H_{ee}|n\rangle\langle n|H_{ee}|0\rangle}{E_{n}-E_{0}}+....,$ (S36)

where $E_{0}$ is the interaction-independent contribution to energy, $H_{ee}$
is the interaction Hamiltonian, index zero denotes the ground state of the
system and index n corresponds to excited states. In an electron gas, the
denominator
$E_{n}-E_{0}=\frac{\hbar^{2}}{m}(q^{2}+({\bm{k}}_{1}-{\bm{k}}_{2})\cdot{\bm{q}})$,
where ${\bm{k}}_{1}$ and ${\bm{k}}_{2}$ are momenta (wave vectors) of the two
holes inside the Fermi sphere, ${\bm{q}}$ is the transferred momentum,
${\bm{k}}_{2}+{\bm{q}}$ and ${\bm{k}}_{2}-{\bm{q}}$ are the momenta of
electrons outside the Fermi sphere. This denominator is divergent in the
second order perturbation theory, and higher order terms are even more
divergent. For the electron gas, these perturbation series can be summed.
However, for the Dirac insulator, when excitations arise only in transitions
between the full valence band and empty conduction band, every denominator in
series Eq. (S36) contains extra constant $2\Delta$, and divergences do not
appear. Furthermore, as we consider high magnetic fields, if both spin states
of the highest occupied Landau level are fully filled, excitation involves a
cyclotron energy gap. Then if the cyclotron energy is much larger than the
characteristic Coulomb interaction, the mixing of excited states is
negligible, and the non-interacting ”closed shell configuration” is
essentially the ground state Giuliani and Vignale (2005). Due to large
$\epsilon$, this happens in rather small magnetic fields in our case. Such
situation is relevant for our consideration of Bi, where we consider the limit
$\Delta=0$ and zero electric field. It becomes clear that electron gas RPA-
like contribution need not be included in our consideration.

Another type of divergence often arising due to interactions is Lippman-
Schwinger Lippmann and Schwinger (1950), Bethe-Salpeter effects Salpeter and
Bethe (1951), or, in condensed matter context, the divergence related to the
exciton bound state. Wannier (1937); Mott (1938); Schrödinger (1926) For
exciton, when $\epsilon$ is big, binding energy is small and the size is
effectively big. Small binding energy means closeness in energy to the bottom
of the conduction band. For this bound state to manifest itself in
susceptibilities, excitations have to be generated at frequencies of
electromagnetic wave close to $2\Delta$. Thus, the related contributions can
be separated from the effects at small frequencies that we mostly consider.

In susceptibilities, divergences as a result of summation of series of
interaction contributions may also indicate phase transitions, such as
ferromagnetism. In the presence of the gap $2\Delta$ for excitations in Dirac
insulators, such phase transition is not expected. Similarly, a
superconducting transition, while emerging due to binding of electrons in
Cooper pairs associated with the presence of the Fermi surface (see, e.g.
Abrikosov Abrikosov (2017)), is unlikely to arise in the presence of
excitation gap in insulators.

We note that the one-loop Heisenberg-Euler contribution at large B that we
calculated, Eq. 2 of the main text, also stems from divergence of the general
renormalization group type, when $\alpha$ is enhanced by the first power of
log.

Finally, there is potentially a divergence stemming from loop expansion that
differ from the cases discussed above. It contains the fine structure constant
$\alpha$ only, and is asymptotic, with combinatorial coefficients growing fast
with increasing the power $n$ of the fine structure constant in diagrams with
$n$ interaction lines and vertices of electromagnetic interactions. These so-
called asymptotic series were considered first by Dyson Dyson (1952) (Also see
an excellent review by Huet and co-authors HUET  _et al._ (2012)). It is now
widely accepted that for asymptotic series a summation up to an infinite order
in small $\alpha$ does not make sense, and one must terminate the series after
a few terms when doing practical calculations with QED, although a rigorous
mathematical justification for using 1-loop QED is an open problem.

Thus, we can safely restrict ourselves to 1-loop approximation since

$\frac{(K+1)\text{-loops},2n\text{-legs}}{K\text{-loops},2n\text{-legs}}\sim\frac{\alpha_{D}}{\epsilon}\approx
0.03,\quad\text{See Eq.~{}\eqref{alpha_estimate} and
Table~{}\ref{table:epsilons}}{.}$ (S37)

This applies to QED where $\alpha\approx 1/137$, as considered by Heisenberg
and Euler, and in the condensed matter contexts, when the interaction
parameter $e^{2}/(\hbar v\epsilon)$ is small at large $\epsilon$, fixed by
lattice properties. We note that our predictions pass the test of comparison
with available experimental results (See Fig. 2 of the main text), even in the
strong field regime which can not be directly probed in QED.

Compound |  $2\Delta$[eV] | $\epsilon$ | Refs.
---|---|---|---
PbTe | $0.19$ | $450$ | Lin and Kleinman (1966); Burstein  _et al._ (1968); Svane  _et al._ (2010); Kanai and Shohno (1963)
PbSe | $0.16$ | $280$ | Lin and Kleinman (1966); Burstein  _et al._ (1968); Svane  _et al._ (2010)
PbS | $0.28$ | $190$ | Lin and Kleinman (1966); Burstein  _et al._ (1968); Svane  _et al._ (2010)
SnTe | $0.27$ | $1770$ | Suzuki and Adachi (1995); Herman  _et al._ (1968); Hayasaka and Fuseya (2016)
Sb2Te3 | $0.21$ | $168,36$ | Lefebvre  _et al._ (1987); Spr (a)
Bi2Te3 | $0.17$ | $290,75$ | Ando (2013); Spr (b)
Bi2Se3 | $0.30$ | $113$ | Ando (2013); Spr (c)
GeBi${}_{4\text{-}x}$SbxTe7 | $0.10$ |  $448,30$ 222222@ $x=0$, i.e. GeBi4Te7 | Ando (2013); Jain  _et al._ (2013)
GeBi2Te4 | $0.18$ | $315,19$ | Ando (2013); Jain  _et al._ (2013)
$\text{Bi}_{1-x}\text{Sb}_{x}$ |  $0.01$ 232323gap vanishes at $x=0.04$ |  $\sim 100$ 242424@ x = 0.07, where $2\Delta\approx 0.01\>\text{eV}$ | Schwerdtfeger and Nagle (2019); Rudolph  _et al._ (1980)
Table S3: The Dirac gap and dielectric constants of Dirac materials. When
multiple numbers for $\epsilon$ are provided, anisotropy is signified (see
references for more information). The linear dielectric contribution from the
valence bands coming from Eq. (S40) is, by definition, already contained in
these measurements.

###  S4 1-Loop effective action

The 1-loop effective action of the material follows as

$S_{eff}=-i\ln\mathrm{Det}\left[v\gamma^{\mu}(\tilde{\partial}_{\mu}-i|e|\mathcal{A}_{\mu})+i\Delta\right]{.}$
(S38)

This determinant is exactly calculated in the QED case Dunne ; *Dunne_arxiv.
By comparing the QED and Dirac insulator actions in Eqs. (S32) and (S33), we
can write the 1-loop, nonperturbative, renormalized Heisenberg-Euler action in
terms of the normalized fields defined in Eqs. (3) and (4) as Akhiezer and
Berestestskii (1965); Berestetskii  _et al._ (1982)

$\displaystyle\delta L_{HE}$ $\displaystyle=$
$\displaystyle\frac{\Delta}{8\pi^{2}\lambdabar_{D}^{3}}\int_{0}^{\infty}\frac{d\eta
e^{-\eta}}{\eta^{3}}\left[(-\eta A\cot(\eta A)\eta C\coth(\eta
C)+1-\frac{1}{3}\eta^{2}(A^{2}-C^{2})\right],$ $\displaystyle A$
$\displaystyle=$ $\displaystyle-\frac{i}{2}\left[\sqrt{({\bf b}+i{\bf
e})^{2}}-\sqrt{({\bf b}-i{\bf e})^{2}}\right],$ $\displaystyle C$
$\displaystyle=$ $\displaystyle\frac{1}{2}\left[\sqrt{({\bf b}+i{\bf
e})^{2}}+\sqrt{({\bf b}-i{\bf e})^{2}}\right]{}.$ (S39)

In addition to this, there is also the cut-off dependent part of the action
that is due to the UV-divergent diagram (shown in Fig. S1 in a red box) which
evaluates to

$\delta
L_{1}=\frac{\Delta}{12\pi^{2}\lambdabar_{D}^{3}}\ln\left(\frac{\Lambda}{\Delta}\right)\left(|\mathbf{e}|^{2}-|\mathbf{b}|^{2}\right).$
(S40)

In QED, this quantity has the same form as the classical electromagnetic
Lagrangian , hence it is renormalized into the definition of electric charge.
However in the material case, it gives rise to the magnetic and electric
susceptibilities that we discuss below absorbed into $\bm{\epsilon},\bm{\mu}$
of the material. The Eq. (S4) is valid for both strong and weak fields. The
contour on the complex $\eta$ plane shall be chosen so that the poles in the
integrand are avoided. The imaginary part of the action gives the breakdown
probability (volume rate of particle hole creation in QED) . For example, when
$|\mathbf{b}|=0$, we have to leading order Berestetskii  _et al._ (1982)

$P=2\text{Im}\>\delta
L_{HE}=\frac{\Delta}{4\pi^{3}\hbar\lambdabar_{D}^{3}}|{\bf
e}|^{2}e^{-\pi/|{\bf e}|},$ (S41)

Meanwhile, another way to write Eq. (S38) is the standard form

$S_{eff}=-i\ln\mathrm{Det}\left[v\gamma^{\mu}(\tilde{\partial}_{\mu}-i|e|\mathcal{A}_{\mu})+i\Delta\right]=-i\text{Tr}\ln[1+|e|vG\gamma^{\mu}\mathcal{A}_{\mu}]-i\ln\mathrm{Det}\left[-iG^{-1}\right]{.}$
(S42)

Carrying out the formal diagrammatic expansion of the logarithm , we have

$1\text{-loop},2n\text{-legs}=\frac{i}{2n}\text{Tr}[(|e|v\gamma^{\mu}\mathcal{A}_{\mu}G)^{2n}],\quad
n=1,2,...$ (S43)

For example, the 2-leg diagram after taking the Fourier transform
$\mathcal{A}(q)=\int
d^{4}\tilde{x}\mathcal{A}(\tilde{x})e^{ik^{\mu}\tilde{x}_{\mu}}$ reads

$1\text{-loop},2\text{-legs}=\frac{ie^{2}v^{2}}{4}\int\frac{d^{4}qd^{4}k}{(2\pi)^{8}}\text{tr}\left[\not{\mathcal{A}}(q)\frac{1}{v\not{k}-\Delta}\not{\mathcal{A}}(-q)\frac{1}{v(\not{k}+\not{q})-\Delta}\right],\quad
k^{\mu}=(\omega/v,\bm{\tilde{k}}),\>\not{k}=\gamma^{\mu}k_{\mu}{,}$ (S44)

which gives the UV divergent term Eq. (S40). The higher order diagrams can be
identified by formally expanding the action Eq. (S4) in the fields. For
example , if we choose $\mathbf{e}\parallel\mathbf{B}$, we have Dunne ;
*Dunne_arxiv

$1\text{-loop},2n\geq
4\text{-legs}=-\frac{\Delta}{8\pi^{3}\lambdabar_{D}^{3}}(2n-3)!\sum_{k=0}^{n}\frac{\mathcal{B}_{2k}\mathcal{B}_{2n-2k}}{(2k)!(2n-2k)!}(2|\mathbf{b}|)^{2n-2k}(2i|\mathbf{e}|)^{2k}{,}$
(S45)

where $\mathcal{B}_{2n}$ are Bernoulli numbers. The expansion parameters are

$|\mathbf{e}|^{2}=\alpha_{D}E^{2}\frac{\lambdabar_{D}^{3}}{\Delta},\quad|\mathbf{b}|^{2}=\alpha_{D}B^{2}\frac{v^{2}}{c^{2}}\frac{\lambdabar_{D}^{3}}{\Delta}{.}$
(S46)

Therefore when the fields or the coupling constant are strong, the
diagrammatic expansion is only schematic because the series Eq. (S45) should
be summed exactly back into the original form Eq. (S4), where the non-coplanar
$\bm{E}$ and $\bm{B}$ case is also taken into account.

###  S5 Regime of validity of 1-loop effective action for a Dirac material

From previous sections, firstly, we have to ensure that the 1-loop
approximation holds by keeping the interaction strength under control

$\alpha_{D}\ll\epsilon\sim 100,\quad\text{(1-loop approximation holds)}{.}$
(S47)

When the magnetic field is too large, the Landau level separation energy
becomes comparable to the cut-off, and the Dirac description breaks down, but
at

$|\bm{B}|\ll\frac{\Lambda^{2}c}{e\hbar v^{2}},$ (S48)

the Dirac theory is valid. Finally when the E-field is strong, a formal
asymptotic expression can be easily derived from Eq. (S4). However, in our
case this becomes a somewhat an academic question, because as seen from Eq.
(S41) the electric breakdown probability increases significantly with strong
fields, which we do not consider. At small fields we have conditions

$\displaystyle|\bm{E}|$ $\displaystyle\ll\frac{\Lambda^{2}}{e\hbar
v},\quad\text{(The Dirac theory is valid)}$ (S49) $\displaystyle|\mathbf{e}|$
$\displaystyle<0.3,\quad\text{(avoids electric breakdown)}.$ (S50)

###  S6 Asymptotic expansion of the Heisenberg-Euler Action

In order to avoid electric breakdown of an insulator, we are always in the
weak $E$ limit, where the parallel field Lagrangian, being quadratic in
$|\mathbf{e}|^{2}$ and $\hat{\bf{e}}\cdot\hat{\mathbf{b}}$, can be expanded as

$f_{n}^{m}(b,x)$ | $n\rightarrow$ |  |
---|---|---|---
$m\downarrow$ | $b^{2}\left(\frac{1}{x^{2}}+\frac{1}{3}-\frac{\coth(x)}{x}\right)$ | $-\left(1+\frac{x}{2b}\right)\left(\frac{1}{x^{2}}+\frac{1}{3}-\frac{\coth(x)}{x}\right)$ | $\frac{x(x+b)}{8b^{4}}\left(\frac{1}{x^{2}}+\frac{1}{3}-\frac{\coth(x)}{x}\right)$
|  | $\frac{1}{2\sinh^{2}(x)}-\frac{3-2x^{2}}{6x}\coth(x)$ | $\left(\frac{1}{b^{2}}-\frac{x}{2b^{3}}\right)\left(\frac{1}{2\sinh^{2}(x)}-\frac{3-2x^{2}}{6x}\coth(x)\right)$
|  |  | $-\frac{6x\coth(x)+4x^{2}+9}{24b^{2}\sinh^{2}(x)}+\coth(x)\frac{225-60x^{2}+8x^{4}}{360b^{2}x}$

.

Table S4: The functions $f_{n}^{m}(x,b)$ appearing in the integrands in
$I^{m}_{n}(b)$ (Eq. (S52)) the perturbation expansion of the renormalized
Lagrangian Eq. (S51)

$I_{n}^{m}(b)$ | $n\rightarrow$ |  |
---|---|---|---
$m\downarrow$ | $\ln(b)\left(\frac{1}{3}b^{2}+b+\frac{1}{2}\right)+b+\frac{3}{4}+\frac{\pi^{2}}{72b}-\frac{\zeta(3)}{48b^{2}}$ | $-\ln(b)\left(\frac{1}{3}+\frac{1}{2b}\right)-\frac{1}{6}-\frac{1}{b}-\frac{1}{4b^{2}}$ | $\frac{1}{12b^{2}}$
|  | $\frac{b}{3}+\frac{1}{2b}\left(\ln(b)+1\right)+\frac{\zeta(3)+3}{12b^{2}}$ | $\frac{1}{6b}$
|  |  | 0

.

Table S5: Strong $B$-field expansion (up to order $1/b^{2}$) of the integrals
$I^{m}_{n}(b)$ (Eq. (S52)) in the perturbation expansion of the renormalized
Lagrangian Eq. (S51)

$\delta
L_{HE}=\frac{\Delta}{8\pi^{2}\lambdabar_{D}^{3}}\sum_{n=0}^{\infty}\sum_{m=0}^{n}|\mathbf{e}|^{2n}(\mathbf{\hat{e}}\cdot\mathbf{\hat{b}})^{2m}I_{n}^{m}(b),\quad
b=|\mathbf{b}|{.}$ (S51)

The functions $I_{n}^{m}$ are integrals of the form

$I_{n}^{m}(b)=\int_{0}^{\infty}dx\frac{e^{-x/b}}{x}f_{n}^{m}(b,x){.}$ (S52)

where the functions $f_{n}^{m}(b,x)$ are tabulated in Table S4.

Making a change of variables $x\to bx^{\prime}$, we can recheck the low field
limit of the Lagrangian in Eq. (10) of the main text

$\frac{8\pi^{2}\lambdabar_{D}^{3}}{\Delta}\delta
L_{HE}\to\frac{1}{45}\int_{0}^{\infty}dx\>xe^{-x}\left(b^{4}+7(\mathbf{e}\cdot\mathbf{b})^{2}-|\mathbf{e}|^{2}b^{2}\frac{x+2}{2}+|\mathbf{e}|^{4}\frac{x+x^{2}}{8}\right)=\frac{1}{45}[(|\mathbf{e}|^{2}-|\mathbf{b}|^{2})+7(\mathbf{e}\cdot\mathbf{b})]{.}$
(S53)

The strong $B$ limit is obtained if we use the fact that for a bounded
function $g(x)\to 0$ when $x\to 0$ the integral converges to

$\int_{0}^{\infty}\frac{dx}{x}e^{-x\delta}g(x)\to
g(\infty)\ln(1/\delta),\quad\delta\to 0{.}$ (S54)

Taking the derivative or integral with respect to $1/b$, we can generate the
subleading or super-logarithmic terms , respectively. Performing this
procedure we obtain the strong $B$ expansions of the integrals in Eqs. (S51)
and (S52), as tabulated in Table S5. Reading of the leading order terms from
Table S5 and substituting in Eq. (S51) we obtain Eq. (2) of the main text

$\frac{8\pi^{2}\lambdabar_{D}^{3}}{\Delta}\delta
L_{HE}\to\frac{1}{3}|\mathbf{b}|^{2}\log|\mathbf{b}|+\frac{1}{3}|\mathbf{e}|^{2}|\mathbf{b}|(\mathbf{\hat{e}}\cdot\mathbf{\hat{b}})^{2}.$
(S55)

Figure S2: Dimensionless functions F,G,D of Eqs. (Nonlinear Quantum Electrodynamics in Dirac materials), (13)) and (Nonlinear Quantum Electrodynamics in Dirac materials) versus the dimensionless magnetic field. Functions G and D depend on the relative orientation of electric and magnetic field. The left (right) panel corresponds to $\mathbf{b}\parallel\mathbf{e}$ ($\mathbf{b}\perp\mathbf{e}$). For arbitrary angle between $\mathbf{b},\mathbf{e}$ see Table S6. $\bm{\delta\chi}$ | $\frac{v^{2}}{c^{2}}\frac{\alpha_{D}}{12\pi^{2}}\bm{\mathcal{U}}^{-2}$ | $2|\mathbf{b}|^{2}/5$ | $|\mathbf{b}|\ll 1$
---|---|---|---
$\log(|\mathbf{b}|)$ | $|\mathbf{b}|\gg 1$
$\bm{\delta\epsilon(B)}$ | $\frac{\alpha_{D}}{3\pi}\bm{\mathcal{U}}^{2}$ | $\frac{7}{15}(\mathbf{\hat{e}}\cdot\mathbf{\hat{b}})^{2}-\frac{2}{15}|\mathbf{b}|^{2}$ | $|\mathbf{b}|\ll 1$
$-\log(|\mathbf{b}|)+(\mathbf{\hat{e}}\cdot\mathbf{\hat{b}})^{2}|\mathbf{b}|$ | $|\mathbf{b}|\gg 1$
$4\pi\bm{M^{(e)}}$ | $\frac{\mu_{D}}{3\pi\lambdabar_{D}^{3}}|\mathbf{e}|^{2}\bm{\mathcal{U}}^{-1}$ | $\frac{2}{15}\left(7\mathbf{\hat{e}}(\mathbf{\hat{e}}\cdot\mathbf{b})-2\mathbf{b}\right)$ | $|\mathbf{b}|\ll 1$
$-\hat{\mathbf{b}}(\mathbf{{\mathbf{\hat{e}}}}\cdot\mathbf{\hat{b}})^{2}+2{\mathbf{\hat{e}}}(\mathbf{{\hat{e}}}\cdot\mathbf{\hat{b}})-\frac{\mathbf{\hat{b}}}{|\mathbf{b}|}$ | $|\mathbf{b}|\gg 1$
Table S6:  Asymptotic expressions for the nonlinear magnetic susceptibility
$\bm{\chi}$, magnetically modulated dielectric constant
$\bm{\delta\epsilon(B)}$ and the electric modulated magnetization
$4\pi\bm{M^{(e)}}$ tensors for arbitrary angle between $\mathbf{\hat{e}}$ and
$\mathbf{\hat{b}}$.

###  S7 Derivation of Susceptibility tensors

Inspecting the classical Lagrangian of the electromagnetic field
$L_{0}=(E^{2}-B^{2})/(8\pi^{2})$, we can derive the (linear and nonlinear)
contributions to susceptibilities from Eq. (S40) and Eq. (S4). If the quantum
part of the action is called $\delta L$, the differential magnetic and
electric susceptibilities are

$\displaystyle\chi_{ij}$ $\displaystyle=\frac{\partial^{2}\delta L}{\partial
B_{j}ialB_{i}}=\frac{v^{2}\alpha_{D}\lambdabar_{D}^{3}}{c^{2}\Delta}\frac{\partial^{2}\delta
L}{\partial\text{b}_{l}\partial\text{b}_{k}}\mathcal{U}^{-1}_{ik}\mathcal{U}^{-1}_{jl}{,}$
(S56a) $\displaystyle\chi_{ij}^{e}$ $\displaystyle=\frac{\partial^{2}\delta
L}{\partial{E}_{j}\partial{E}_{i}}=\frac{\alpha_{D}\lambdabar_{D}^{3}}{\Delta}\frac{\partial^{2}\delta
L}{\partial\text{e}_{l}\partial\text{e}_{k}}\mathcal{U}_{ik}\mathcal{U}_{jl}{.}$
(S56b)

####  S7.1 Linear and low order nonlinear susceptibilities

The cutoff dependent UV contribution is

$\delta
L_{1}=\frac{\Delta}{12\pi^{2}\lambdabar_{D}^{3}}\ln\left(\frac{\Lambda}{\Delta}\right)\left(|\mathbf{e}|^{2}-|\mathbf{b}|^{2}\right).$
(S57)

In vacuum quantum electrodynamics , this quantity has the same form as the
classical electromagnetic Lagrangian, hence it is absorbed into the definition
of electric charge. However in the material case, it gives rise to the
isotropic valence band contribution to the magnetic and electric
susceptibilities, respectively

$(\tilde{\chi}_{1})_{ij}=-\frac{v^{2}}{c^{2}}\frac{\alpha_{D}}{6\pi^{2}}\ln\left(\frac{\Lambda}{\Delta}\right)\delta_{ij}=-\frac{v^{2}}{c^{2}}(\tilde{\chi}^{e}_{1})_{ij}{.}$
(S58)

In the anisotropic case, these become

$\bm{\chi_{D}}=-\frac{v^{2}}{c^{2}}\frac{\alpha_{D}}{6\pi^{2}}\ln\left(\frac{\Lambda}{\Delta}\right)\bm{\mathcal{U}}^{-2},\quad\bm{\epsilon_{D}}=\bm{1}+4\pi\bm{\chi_{D}^{e}}=\bm{1}+\frac{2\alpha_{D}}{3\pi^{2}}\ln\left(\frac{\Lambda}{\Delta}\right)\bm{\mathcal{U}}^{2}{.}$
(S59)

as in Eqs. (7) and (8) of the main text.

The weak field limit of Eq. (S4) contains the lowest order nonlinear terms
that are quadratic in $|\mathbf{e}|^{2}-|\mathbf{b}|^{2}$ and
$\mathbf{e}\cdot\mathbf{b}$

$\delta
L_{2}=\frac{\Delta}{360\pi^{2}\lambdabar_{D}^{3}}\left[\left(|\mathbf{e}|^{2}-|\mathbf{b}|^{2}\right)^{2}+7\left({\bf
e}\cdot{\bf b}\right)^{2}\right],$ (S60)

from which the nonlinear differential susceptibility tensors for an isotropic
system follow as

$\displaystyle\delta\tilde{\chi}_{ij}$
$\displaystyle\to\frac{v^{2}}{c^{2}}\frac{\alpha_{D}}{180\pi^{2}}\left[2|\mathbf{b}|^{2}(2\hat{b}_{i}\hat{b}_{j}+\delta_{ij})+|\mathbf{e}|^{2}(7\hat{e}_{i}\hat{e}_{j}-2\delta_{ij})\right],$
(S61a) $\displaystyle\delta\tilde{\epsilon}_{ij}$
$\displaystyle\to\frac{\alpha_{D}}{45\pi}\left[|\mathbf{b}|^{2}(7\hat{b}_{i}\hat{b}_{j}-2\delta_{ij})+2|\mathbf{e}|^{2}(2\hat{e}_{i}\hat{e}_{j}+\delta_{ij})\right],\quad\hat{e}_{i}=\frac{\text{e}_{i}}{|\mathbf{e}|},\>\hat{b}_{i}=\frac{\text{b}_{i}}{|\mathbf{b}|}$
(S61b)

These quantities can be re-expressed in the anisotropic case by using Eq.
(S1.1) as

$\displaystyle\bm{\delta\chi}$
$\displaystyle=\frac{v^{2}}{c^{2}}\frac{\alpha_{D}}{180\pi^{2}}\left(4(\bm{\mathcal{U}^{-1}}\mathbf{b})\otimes(\bm{\mathcal{U}^{-1}}\mathbf{b})+2|\mathbf{b}|^{2}\bm{\mathcal{U}}^{-2}+7(\bm{\mathcal{U}^{-1}}\mathbf{e})\otimes(\bm{\mathcal{U}^{-1}}\mathbf{e})-2\bm{\mathcal{U}}^{-2}|\mathbf{e}|^{2}\right),$
(S62a) $\displaystyle\bm{\delta\epsilon}$
$\displaystyle=\frac{\alpha_{D}}{45\pi}\left(4(\bm{\mathcal{U}}\mathbf{e})\otimes(\bm{\mathcal{U}}\mathbf{e})+2|\mathbf{e}|^{2}\bm{\mathcal{U}}^{2}+7(\bm{\mathcal{U}}\mathbf{b})\otimes(\bm{\mathcal{U}}\mathbf{b})-2\bm{\mathcal{U}}^{2}|\mathbf{b}|^{2}\right).$
(S62b)

The electric modulated magnetization then is simply
$\bm{M^{(e)}}=\bm{\delta\chi(E)B}=B_{\star}\bm{\delta\chi(E)\mathcal{U}}\mathbf{b}$
and reads

$4\pi\bm{M^{(e)}}=\frac{2\mu_{D}}{45\pi\lambdabar_{D}^{3}}\bm{\mathcal{U}}^{-1}\left(7\mathbf{e}(\mathbf{e}\cdot\mathbf{b})-2\mathbf{b}|\mathbf{e}|^{2}\right).$
(S63)

To obtain the expressions in the main text we work in coordinates where

$\bm{\mathcal{U}}=\text{diag}(U_{xx},U_{yy},U_{zz}).$ (S64)

For simplicity we assume $\bm{E}$ and $\bm{B}$ directed along principal
directions, e.g. $\bm{E}\parallel\bm{B}\parallel\mathbf{\hat{x}}$ for the
parallel field configuration and
$\bm{E}\parallel\mathbf{\hat{x}},\>\>\bm{B}\parallel\mathbf{\hat{y}}$ for the
perpendicular field configuration. When these conditions are met, we can
reduce, for example the parallel weak field magnetic susceptibility as

$(\bm{\delta\chi_{\parallel}(E)})_{xx}=\bm{\mathcal{U}}^{-2}_{xx}\frac{v^{2}}{c^{2}}\frac{\alpha_{D}}{36\pi^{2}}|\mathbf{e}|^{2},\quad
4\pi(\bm{M^{(e)}_{\parallel}})_{x}=(\bm{\mathcal{U}}^{-1}\mathbf{b})_{x}\frac{2\mu_{D}}{9\pi\lambdabar_{D}^{3}}|\mathbf{e}|^{2}.$
(S65)

For a summary of susceptibilities for arbitrary mutual orientation of E&B see
Table S6. As we discussed in section S1, using inverse anisotropy
transformation, the results for susceptibilities can be re-written in terms of
applied external fields $\bm{E}$ and $\bm{B}$. The angle between $\mathbf{e}$
and $\mathbf{b}$ as function of the angle between $\bm{E}$ and $\bm{B}$ is
given in Eq. (S24).

####  S7.2 Higher nonlinear corrections to susceptibilities

Now that we have the renormalized Lagrangian in the form Eq. (S51), with the
integrals $I^{m}_{n}$ tabulates in Table S5, we can calculate the
susceptibilities according to Eq. (S56), in notable cases. If the electric
field is zero , we have

$\frac{8\pi^{2}\lambdabar_{D}^{3}}{\Delta}\delta
L_{HE}(|\mathbf{e}|=0)\to\frac{\mathbf{b}^{2}}{3}\ln(\mathbf{b})+\mathbf{b}\ln(\mathbf{b})+...$
(S66)

consistent with Eq. (2) of the main text. Therefore the contribution to
magnetic susceptibility due to applied magnetic field is

$\bm{\delta\chi}(|\mathbf{e}|=0)\to\frac{v^{2}}{c^{2}}\frac{\alpha_{D}}{12\pi^{2}}\bm{\mathcal{U}}^{-2}\ln|\mathbf{b}|{,}$
(S67)

as in Eq. (Nonlinear Quantum Electrodynamics in Dirac materials) of the main
text.

The term in the Lagrangian that is second order in the electric field is

$\frac{8\pi^{2}\lambdabar_{D}^{3}}{\Delta}\frac{\partial\delta
L_{HE}}{\partial|\mathbf{e}|^{2}}\bigg{|}_{|\mathbf{e}|=0}=-\ln(|\mathbf{b}|)\left(\frac{1}{3}+\frac{1}{2|\mathbf{b}|}\right)+(\mathbf{\hat{\mathbf{e}}}\cdot\mathbf{\hat{b}})^{2}\left(\frac{|\mathbf{b}|}{3}+\frac{1}{2|\mathbf{b}|}\left(\ln|\mathbf{b}|+1\right)\right)+...$
(S68)

Then the leading order magnetization contribution is obtained from $M=\partial
L/\partial B$ as

$4\pi\bm{M}^{(\bm{e})}\to\frac{\mu_{D}}{3\pi\lambdabar_{D}^{3}}|\mathbf{e}|^{2}\bm{\mathcal{U}^{-1}}\left(-\hat{\mathbf{b}}(\mathbf{{\mathbf{\hat{e}}}}\cdot\mathbf{\hat{b}})^{2}+2{\mathbf{\hat{e}}}(\mathbf{{\hat{e}}}\cdot\mathbf{\hat{b}})-\frac{\mathbf{\hat{b}}}{|\mathbf{b}|}...\right),$
(S69)

as in Eq. (Nonlinear Quantum Electrodynamics in Dirac materials) of the main
text. The dielectric response, linear in the electric field, is due to the
term in the Lagrangian that is second order in the electric field , from which
we obtain the magnetic field contribution in the isotropic case as

$\bm{\delta{\epsilon}}\to\frac{\alpha_{D}}{3\pi}\left(-\ln(|\mathbf{b}|)(\bm{\mathcal{U}}\mathbf{\hat{e}})\otimes(\bm{\mathcal{U}}\mathbf{\hat{e}})+(\bm{\mathcal{U}}\mathbf{\hat{b}})\otimes(\bm{\mathcal{U}}\mathbf{\hat{b}})|\mathbf{b}|\right){.}$
(S70)

The leading order electric field contribution to the dielectric tensor comes
from the term that is $\sim|\mathbf{e}|^{4}$ in the Lagrangian

$\delta{\tilde{\epsilon}}_{ij}\to\frac{\alpha_{D}}{6\pi}\frac{|\mathbf{e}|^{2}}{\mathbf{|b|}^{2}}\left(\delta_{ij}+2\hat{e}_{i}\hat{e}_{j}+|\mathbf{b}|\delta_{ij}(\mathbf{\hat{e}}\cdot\mathbf{\hat{b}})^{2}+2|\mathbf{b}|[\hat{e}_{i}\hat{b}_{j}+\hat{e}_{j}\hat{b}_{i}](\mathbf{\hat{e}}\cdot\mathbf{\hat{b}})+\mathbf{|b|}\hat{e}_{i}\hat{e_{j}}\right){,}$
(S71)

hence small in the parameter $|\mathbf{e}|^{2}/\mathbf{|b|}\ll 1$.

###  S8 Material Applications

####  S8.1 Bismuth and $\text{Bi}_{0.9}\text{Sb}_{0.1}$

While our consideration is primarily for insulators , we, e.g., demonstrated
that the calculated magnetic susceptibility would be valid for valence band
contributions in gapless materials. However real experiments often include,
particularly for insulators with small gap, a certain density of free
electrons. Then analysis of experimental settings must take into account
contributions from free electrons, or experimental conditions must be found
when these contributions are suppressed.

Bismuth has a band gap of $2\Delta=15.5\>\text{meV}$ at the $L$-point and a
Fermi level of $\mu=35\>\text{meV}$ Fuseya  _et al._ (2015) measured from the
midgap point of the Dirac bands. The electronic Fermi surface is composed of 3
electron ellipsoids that lie on the binary-bisectrix (x-y) plane perpendicular
to the trigonal (z) axis. There is also a hole pocket along the trigonal axis.
The hole Fermi level is about $-196\>\text{meV}$ measured from the midgap
point of the hole Dirac bands. The hole band gap is about
$2\Delta=370\>\text{meV}$. Liu and Allen (1995)

The hole contribution to susceptibility in the binary-bisectrix plane is small
$\chi\sim-10^{-7}$ due to the large hole band gap and the alignment of the
hole ellipsoid.

The conduction bands are polarized in magnetic field, so that only the e1
ellipsoid is populated when a field is applied in the binary direction Iwasa
_et al._ (2019). Furthermore, the e1 conduction electrons can exhibit de Haas-
van Alphen oscillations, which are suppressed above $B=5\>\text{T}$. At higher
magnetic fields the nonlinear diamagnetism is identical to the insulating
alloy $\text{Bi}_{0.9}\text{Sb}_{0.1}$, which is identical to bismuth except
for the absence of the Fermi level.

Based on the band structure calculations Liu and Allen (1995), we write the
diagonalized velocity tensor as

$\bm{\mathcal{V}}=\text{diag}(1.7,\>1.5,\>0.4)v,\quad c/v=188,$ (S72)

The Dirac cone is located at the L-point, therefore the velocity operator in
the crystal frame is obtained by using an appropriate rotation-reflection
operator. The critical field is about $B_{\star}=36\>\text{mT}$. We find the
total susceptibility by adding the contributions from each L-point located
symmetrically in the binary axis

$\bm{\chi_{1}}+\bm{\delta\chi}(\bm{B})=\sum_{n=1}^{3}(\bm{R}^{T})^{n}\left[\bm{\chi_{1}}+\bm{\delta\chi}\left(\frac{\bm{\mathcal{U}}^{-1}\bm{R}^{n}\bm{B}}{B_{\star}}\right)\right]^{(L-point\>n)}\bm{R}^{n-1},$
(S73)

where $R$ implements rotation by $120^{\circ}$ about the trigonal axis. The
first order contribution is

$\bm{\chi_{D}}=-10^{-5}\begin{pmatrix}3&0&0\\\ 0&3&0\\\ 0&0&0.5\end{pmatrix}.$
(S74)

If we apply magnetic field in the binary ($x$) direction we have

$\delta\chi_{xx}(B_{x}=5\>\text{T})=13\times 10^{-6}.$ (S75)

For the full field range see Fig. 2. The other quantities (dielectric
enhancement, magnetization etc.) are calculated in a similar fashion. For
example the dielectric enhancement when $E\parallel B$ is in $x$-direction we
have

$\delta\epsilon_{\parallel}(B_{x})\to 10B_{x}.$ (S76)

Given that the dielectric constant of the bismuth alloy is about
$\epsilon~{}100$ (See Table S3) arising mostly from the ionic crystal, the
magnetic field induced nonlinear contribution due to the Dirac band is
enormous even at relatively weak B-fields.

To compare the nonlinear behavior of our theory with the experiment by Iwasa
et al. Iwasa  _et al._ (2019), we take the reference level of susceptibility
to be $\chi_{ref}=-\chi(50\>\text{T})$.

According to our theory

$\chi_{theory}^{binary}(50\>\text{T})=-5.5\times
10^{-6},\quad\chi_{theory}^{bisectrix}(50\>\text{T})=-8.0\times 10^{-6}.$
(S77)

The Iwasa experiment has the reference points

$\chi_{exp}^{binary}(50\>\text{T})=-18.7\times
10^{-6},\quad\chi_{exp}^{bisectrix}(50\>\text{T})=-15.4\times 10^{-6}.$ (S78)

The difference in the reference levels are due to the additional contributions
to linear magnetic susceptibility that are not contained in the Dirac theory.
Although the logarithmic Dirac contribution in Eq. (8) constitutes a
significant part of the linear diamagnetic response, non-Dirac core shell
electron contributions can be equally important. Furthermore, there can be
additional errors due to the choice of UV cut-off and subleading terms in the
summation of the Landau levels can contribute to the total. None of the above
mentioned factors affect the nonlinear response, which is dominated by the
Dirac bands. Therefore, once the overall constant offset due to linear
susceptibility is adjusted, we get an excellent agreement in the nonlinear
behavior of the measured susceptibility in bismuth, as seen in Fig. 2(a).

####  S8.2 TaAs

TaAs has 12 pairs of Weyl fermions Weng  _et al._ (2015), 4 pairs forming hole
pockets $2\>\text{meV}$ below the Fermi level and 8 forming electron pockets
$21\>\text{meV}$ above the Fermi level. The Fermi level $E_{F}=21\>\text{meV}$
contribution which is not accounted in our theory, becomes unimportant
$B\gg\frac{4E_{F}^{2}}{v_{1}v_{2}e\hbar}\approx 10\>T.$

We represent each pair by a single massless Dirac fermion with the velocity
tensor comparable to the calculations and measurements Zhang  _et al._ (2019);
Weng  _et al._ (2015)

$\bm{\mathcal{V}}=\text{diag}(1.7,\>1.7,\>0.35)v,\quad c/v=447$ (S79)

where the the z-component is the velocity in the $c$-axis. Since the system is
gapless we consider the situation where $E=0$. The nonlinear susceptibility is
the gapless limit is

$\bm{\chi_{D}+\delta\chi}\to-\frac{\alpha_{D}}{12\pi^{2}}\frac{v^{2}}{c^{2}}\ln\left(\frac{c\Lambda^{2}}{eB\hbar
v^{2}}\right)\bm{\mathcal{U}}^{-2}$ (S80)

To compare the nonlinear behavior of our theory with the experiment by Zhang
et al. Zhang  _et al._ (2019), we take the reference level of susceptibility
to be $\chi_{ref}=-\chi(30\>\text{T})=1.56\times 10^{-5}$.

####  S8.3 $\text{Pb}_{0.8}\text{Sn}_{0.2}\text{Te}$

The gap as a function of the Sn content is Hayasaka and Fuseya (2016)

$2\Delta[\text{meV}]=182-480x$ (S81)

At $x=0.235$ we have the inverse masses $m^{-1}_{\perp}=100m^{-1}_{e}$ and
$m^{-1}_{\parallel}=11.25m_{e}^{-1}$ where the gap is
$2\Delta=69.1\>\text{meV}$. At $x=0.510$ we have the same inverse massed where
the gap is similar $2\Delta=-62.8\>\text{meV}$ (TI phase). We will base our
estimates on these values . The parallel denotes the z-axis aligned in the
[001] direction. The velocity tensor is then

$\bm{\mathcal{V}}=\text{diag}(1.4,\>1.4,\>0.5)v,\quad c/v=580.$ (S82)

There are a total of 4 Dirac fermions located at the L-points symmetric about
the z-axis [001] direction of the rock salt structure. Due to the relatively
large gap the most striking property is the electric field modulated
magnetization. If we apply $\bm{B}\parallel\bm{E}$ in the $x$-direction with
$B\sim 5\>T$ and $E=E_{\star}/3$ we have

$4\pi({M^{e}_{\parallel}})_{x}(B_{x}=5\>\text{T},E_{x}=10^{4}\>\text{V/cm})\approx
0.5\>\mu\text{T}.$ (S83)