# Experimental validation of a micromechanically-based compaction law for
soft/hard grain mixtures

Manuel Cárdenas-Barrantes<EMAIL_ADDRESS>Laboratoire de Mécanique et Génie Civil, UMR 5508 CNRS-University Montpellier,
34095 Montpellier, France  Jonathan Barés<EMAIL_ADDRESS>Laboratoire de Mécanique et Génie Civil, UMR 5508 CNRS-University Montpellier,
34095 Montpellier, France  Mathieu Renouf<EMAIL_ADDRESS>Laboratoire de Mécanique et Génie Civil, UMR 5508 CNRS-University Montpellier,
34095 Montpellier, France  Émilien Azéma<EMAIL_ADDRESS>Laboratoire de Mécanique et Génie Civil, UMR 5508 CNRS-University Montpellier,
34095 Montpellier, France  Institut Universitaire de France (IUF), Paris,
France

###### Abstract

In this letter, we report on an experimental study which analyzes the
compressive behavior of 2D bidisperse granular assemblies made of soft
(hyperelastic) and hard grains in varying proportions ($\kappa$). By means of
a recently developed uniaxial compression set-up [1] and using advanced
Digital Image Correlation (DIC) method, we follow, beyond the jamming point,
the evolution of the main mechanical observables, from the global scale down
to the strain field inside each deformable grain. First, we experimentally
validate and extend to the uni-axial case a recently proposed micro-mechanical
compaction model linking the evolution of the applied pressure $P$ to the
packing fraction $\phi$ [2]. Second, we reveal two different linear regimes
depending on whether the system is above or below a cross-over strain
unraveling a transition from a discrete to a continuous-like system. Third,
the evolution of these linear laws are found to vary linearly with $\kappa$,
up to a saturation point around $\kappa=80$% of hard particles. These results
provide a comprehensive experimental and theoretical framework that can now be
extended to a more general class of polydisperse soft granular systems.

Granular materials are ubiquitous in nature and human activities. In their
most generic form, they are composed of grains of shape, dimension and more
importantly bulk properties that are very different, even in the same packing.
Among them, granular systems composed of squeezable compounds are the ones
with the most singular behaviors. This is the case of biological tissues
composed of soft cells [3, 4], liquid foams [5, 6], emulsions [7, 8], or
sintered material [9, 10] to name a few. The specificity of materials composed
of both soft and hard particles lies in the fact that under any loading, and
even deep in the jammed state [11], two mechanisms compete in the system:
particle rearrangement and deformation. When loaded, hard particles have only
the ability to rearrange abruptly [12, 13] while it is easier for soft
particles to deform in order to sustain a given strain [1]. In these
materials, both mechanisms cohabit, prioritizing either one or the other
depending on the soft to hard particle ratio and global loading.

Compressed systems made of purely hard or slightly deformable particles have
been extensively studied these last decades from a theoretical [11, 14, 15],
numerical [16, 17] and experimental [18, 19, 20] point of view. Most of the
fundamental aspects of their behavior under many solicitations are now well-
known [11, 21, 15]. The behavior of systems made of highly deformable grains
at high packing fraction is also more and more studied [8, 1, 7, 22, 23] and
at least in compaction, reliable models exist [2]. In between, and even if
they are very common materials, systems composed of grain with very different
– soft and hard – rheologies, driven far above the jamming transition, have a
bewildering but fascinating behavior that stays mostly misunderstood. They
have been experimentally studied in very specific applications like for stress
releasing [24, 25], seismic isolation [26, 27] or foundation damping [28, 29].
However, to our best knowledge, no local measurements have been performed to
understand micro-processes leading to their very specific macroscopic
behavior. Only very recently, through numerical approach, a micro-mechanical-
based compaction model has been proposed [22, 23] to describe the evolution of
these systems in compression. But an experimental validation is still lacking
as well as a clear understanding of local processes.

In this Letter, we aim to fill this gap by extending the recently proposed
compaction law and experimentally validating it, and by stating constitutive
laws linking both local structural parameters and local strain field with
global observables such as pressure and packing fraction. The evolution of
these laws as a function of the fraction of soft grain is also investigated,
and saturation limits observed.

Figure 1: (color online) a: Experimental set-up. A bidimensional bidisperse
granular system, composed of soft and rigid particles, lays on the top glass
of a flatbed scanner. A uniaxial compression device stresses it stepwisely
while it is imaged from below by the scanner. b: Composite view of measured
fields. Rigid particles are colored in blue. Raw image, in gray level, is
shown on the left, particle boundaries are in red. Von Mises strain field,
$\mathcal{C}$, is shown on the right with a color scale going from dark red
(low value) to yellow (high value). Contacts are shown in white. c: Evolution
of the measured global force, $f$, as a function of the global strain,
$\varepsilon$, for different softness ratios, $\kappa$, given in the inset.

Experiment – The experiments were carried out using a set-up partially
introduced by [1, 30]. As shown in fig. 1a, it consists of a bidimensional
piston of initial dimension $270\times 202$ mm2 which compresses uniaxially a
bidisperse collection of rigid and soft cylinders laying on top of a flatbed
scanner. The loading piston is composed of a stepper motor rotating a screw
which translates the moving edge of the piston in the inward direction. Two
force sensors are attached to this latter. They record the global force ($f$)
evolution (see fig. 1c) when compressing the granular media. The induced
global stress ($P$) is directly deduced and measured continuously at a
frequency of $100$ Hz, while the system is compressed stepwisely. For each
loading step the piston moves of $0.5$ mm at a speed of $2$ mm/min to stay in
the quasistatic regime. Then, the system relaxes during $1$ min and is imaged
from below with the scanner [31] at a resolution of $2400$ dpi ($10.6$
$\mu$m/px).

The compressed granular systems are composed of both hyperelastic silicone
cylinders [32, 33] and rigid PVC cylinders. Their Young moduli are
$E_{0}=0.45$ MPa [34] and $1.2$ GPa, respectively and the Poisson ratio of
silicone is $\nu=0.5$. Both soft and rigid granular packings are bidisperse
cylinders of diameters $20$ mm and $30$ mm, and height $15$ mm. For each
experiment, about $n=100$ particles are used, the ratio between rigid and
small ones is the softness ratio, $\kappa=\rm{n}_{\rm{soft}}/\rm{n}$. This
ratio is equally spread among soft and rigid particles, and is varied from
$0.2$ to $1$. To avoid basal friction the scanner glass is coated with oil
[30]. The bottom of each particle is coated with thin metallic glitter, which
creates a random pattern with a correlation length of about $50$ $\mu$m [1].

For each experiment a set of about $90$ pictures of $\sim 500$ Mpx displays
the evolution of the bottom face of the granular system. These pictures are
post-processed with an algorithm modified from [1, 35]. First, particles are
detected individually by thresholding the undeformed picture. Each particle is
then tracked along the full set of pictures. Then, from the particle’s solid
rigid motion measurement, subsets of images are extracted following each
particle and correcting its translation and rotation. For soft particles a DIC
algorithm aimed for large deformations already presented in [34] is used to
obtain the displacement field ($\vec{u}$) inside each particle. For the rigid
particles a more classical DIC algorithm is used. This latter correlates all
images with the initial one as classically done in the small deformation
assumption.

As shown in fig. 1b, the particle boundaries are obtained from the
displacement fields $\vec{u}$ [1]. The system packing fraction is directly
deduced from it as well as the particle aspherisities, $a=p^{2}/(4\pi s)$,
where $p$ and $s$ are the perimeter and surface of the particles,
respectively. For the soft particles, the right Cauchy-Green strain tensor
field, $\vec{\vec{C}}$ is computed [1, 36] and its von Mises measure,
$\mathcal{C}$, is deduced. For the rigid particles the local stresses are too
low to induce any significant deformation of the PVC [30]. Contacts between
particles and their length ($l$) are measured from the proximity of the
boundaries and the von Mises strain [1]. The jamming transition [11] at
packing fraction $\phi_{J}$ is detected for each system for $P=1$ kPa and
correspond with $Z_{J}=3.77\pm 0.05$ (see SM) in agreement with the
slipperiness of the particles [15].

Figure 2: (color online) a: Evolution of the coordination, $Z$, relative to
the coordination at jamming, $Z_{J}\approx 3.8$, as a function of the distance
to the jamming transition for the packing fraction, $\phi-\phi_{J}$. The
dashed line corresponds to a power-law with exponent $0.7$. b: Evolution of
the global stress applied to the system, $P$, as a function of
$\phi-\phi_{J}$. Plane line curves correspond with eq.1. In both panels, the
different curves correspond with different softness ratios, $\kappa$, given in
inset of (b).

Results – In the fig. 2a, we show how the increase of the coordination from
the jamming point, $Z-Z_{J}$, evolves with the relative packing fraction,
$\phi-\phi_{J}$. Whatever the softness ratio, curves collapse on a single one
fitted by a power-law with exponent $\alpha=0.7\pm 0.15$ and prefactor
$k=4.8$. Within the error-bar, this value is close from what has been already
observed [37, 38, 33].

Following a recent micro-mechanical compaction model numerically tested on
isotropically compressed soft grains [2, 22, 22], we extended it to the case
of uniaxial compression (See supplemental material). In fig. 2b, compaction
curves ($P$ vs. $\phi-\phi_{J}$) are plotted with :

$P\simeq\frac{E^{*}}{4\Gamma}(1+\mu_{M})\frac{\phi_{\rm{max}}-\phi_{J}}{\phi_{J}}(Z_{0}+k(\phi-\phi_{J})^{0.7})\phi\ln\left(\frac{\phi_{\rm{max}}-\phi_{0}}{\phi_{\rm{max}}-\phi}\right)$
(1)

where $E^{*}=E_{0}/(2\kappa(1-\nu^{2}))$ is the mean effective material Young
modulus and $\phi_{\rm{max}}$ is the asymptotic maximum packing fraction. For
$\kappa=1$, $\phi_{\rm{max}}=0.97$ is close to $1$ while for $\kappa=0.2$,
$\phi_{\rm{max}}=0.88$ closer to the random close packing (See supplemental
material). $\Gamma$ is a geometrical parameter directly measured on the local
vs. global strain curve (See supplemental material). $\mu_{M}$ is the
macroscopic friction known to be close to 0.25 for an assembly of 2D particles
[39, 40]. Although not based on fitting parameters, the prediction given by
this model is excellent for all $\kappa$-values, we only note small
disagreements with experimental data arising when the grains sharply
rearrange. This effect is also exacerbated by the fact that the system is
relatively small.

Figure 3: (color online) a: Linear evolution of the soft particle average
asphericity, $\left<a\right>$, as a function of the global stress applied to
the system, $P$. Their slope is given in the inset as a function of the
softness ratio, $kappa$. b: Linear evolution of the average value of the
relative contact length, $\left<l\right>$, as a function of the distance to
the jamming transition for the packing fraction, $\phi-\phi_{J}$. Their slope
is given in the inset as a function of $\kappa$. In both panels the different
curves correspond with different $\kappa$, given in inset of (b).

Fig. 3a shows how the average soft particles asphericity, $\left<a\right>$,
increases with global stress, $P$. Apart from few sharp grain rearrangements,
this increase is linear whatever the softness ratio, $\kappa$. The slope
increases linearly with $\kappa$ up to $0.8$ and then plateaus around $5\times
10^{-7}$. This means the fewer particles are soft, the more they are deformed
up to a point where this deformation saturates. In Fig. 3b we show the
evolution of the average relative contact length, $\left<l\right>$. It
increases linearly with the packing fraction for any $\kappa$ and the
corresponding slope also increases linearly with this softness ratio before
saturating around $2.8$ from $\kappa\approx 0.8$.

Figure 4: (color online) a-b: Evolution of the global stress, $P$, and of the
relative packing fraction $\phi-\phi_{J}$, as a function of the average value
of the von Mises strain, $\left<\mathcal{C}\right>$ in soft particles, for
different $\kappa$ values, given in inset of (b) (for the three panels). In
panels (a) and (b) straight black dashed lines show a slope of $2.2$ kPa and
$5.4$ respectively as a guide for the eyes. The vertical gray line and the
shaded area show the cross-over strain $\mathcal{C}_{c}=0.991\pm 1.2\times
10^{-3}$ and its errorbar, respectively. It splits the space horizontally
between small and large deformation levels. c: Evolution of the standard
deviation of the von Mises strain in soft particles as a function of the
relative packing fraction $\phi-\phi_{J}$. The straight dashed line with slope
$4.0\times 10^{-4}$ is a fit of the curves collapsed in their linear regime.
These curves leave their linear regime at certain cut-off values given in the
inset as a function of the softness ratio, $\kappa$.

Under the small deformation assumption the infinitesimal strain tensor,
$\vec{\vec{\varepsilon}}$, is related to $\vec{\vec{C}}$ following:
$\vec{\vec{\varepsilon}}=1/2(\vec{\vec{C}}-\vec{\vec{I}})$, with
$\vec{\vec{I}}$ the second order identity tensor. Hence, as shown in fig. 4a
and b, when the system is compressed, $\mathcal{C}$ decreases from $1$ while
the pressure and the packing fraction increase. At low compression level, for
$\left<\mathcal{C}\right>>\mathcal{C}_{c}\approx 0.991$, the average von Mises
strain, $\left<\mathcal{C}\right>$, decreases linearly with the packing
fraction. In this regime, whatever the softness ratio, curves collapse fairly
well. On the contrary at high compression level, for
$\left<\mathcal{C}\right><\mathcal{C}_{c}$, it decreases linearly with the
global pressure. Not only the average strain evolves during compression, but
its distribution also gets wider [1] (See supplemental material). In fig. 4c,
we show how the standard deviation of $\mathcal{C}$ increases with the packing
fraction. When shifted from $\phi_{J}$, curves collapse on a linear law
whatever $\kappa$, up to a cut-off which increases linearly with $\kappa$ up
to a saturation point near $\kappa=0.8$.

_Discussion & conclusion_ – First, we point out that, even if the system seems
small, macroscopic measurements stay repeatable. For example, the coordination
at the jamming point is detected in a narrow band of $2.5\%$ for all the
experiments. The exponent $\alpha$ of the $Z-Z_{J}$ vs. $\phi-\phi_{J}$ power-
law is slightly overestimated compared to previous experimental and numerical
studies [37, 33, 2, 22, 23]. Still, considering the large error bars, it stays
in the same domain: $\alpha\approx 0.5$. As explained in [2, 22] this law can
be plugged to existing models [41] to obtain an isotropic compaction equation.
In the specific case of grain mixture, $\kappa$ is introduced in the effective
Young Modulus to consider the softness variation. Also, to consider the
uniaxiality of the compression, the induced effective friction is added in the
law prefactor to obtain eq.1, and describe the compaction curve with no
fitting parameters. In this latter equation, $\phi_{\rm{max}}$ goes ideally
from $1$ for purely soft systems to $\phi_{J}$ for purely rigid systems, which
is what we observe in our experimental data. So it is here remarkable that
this model matches very well the experimental data whatever the grain mixture.
This constitutes a validation of the model.

On top of extending and validating this compaction model, we introduce new
relations between global observables and deformed local structures. Hence
linear relations are observed between the particle asphericity and the global
pressure ($\left<a\right>\sim P$) as well as for the relative contact length
and the packing fraction ($\left<l\right>\sim(\phi-\phi_{J})$). This first
observation permits indirect measurement of pressure in a granular system just
by measuring the grain boundary deformations. For this latter linear relation,
in the small deformation regime, a scaling consideration of the Hertz contact
law as well as $Z$ vs. $\phi-\phi_{0}$ equation permits to predict this linear
relation (See supplemental material). It is important to note that this
linearity persists far beyond the small deformation regime and even for any
particle mix, which makes it useful to indirectly deduce the packing fraction
evolution of a compressed system from the contact geometry observation. In
both linear relations, it is also remarkable that linear coefficients increase
linearly with the softness ratio up to a point where it saturates near
$\kappa=0.8$. This saturation comes from the fact that when $\kappa$ gets
higher, for a given load or a given packing fraction, most of the deformation
is carried out by a smaller number of grains which highly deform, leaving more
rapidly the small deformation assumption and entering in a different,
saturated, regime.

So far, never observed relations between global observables and the local
deformation field are also revealed in our experimental results. For a small
level of deformation, the average strain inside particles evolves linearly
with the packing fraction ($\left<\mathcal{C}\right>\sim\phi-\phi_{J}$) while
in the case of an important loading, it scales with the global pressure
applied to the system ($\left<\mathcal{C}\right>\sim P$). This evidences two
distinct regimes deep in the jammed state as already observed in [1] and
separated by a cross-over value $\mathcal{C}_{c}\approx 0.991$ independent of
$\kappa$. In the first regime, linearity is reminiscent of the fact that the
global strain scale with the packing fraction (See supplemental material), so
this seems to be extrapolated to the local strain. In the second regime,
linearity between strain and stress suggests that the material behaves like a
bulk one (except for small rearrangements) this is consistent with the fact
that the material is very dense with almost no interstitial porosity. It is
also consistent with the fact that, at large compression levels, the pressure
seems to decrease faster with the strain for lower softness ratio, the
material gets stiffer for lower $\kappa$. On the contrary, linear relations
collapse in the case of small deformations and do not significantly depend on
the softness ratio. Finally, a scaling between strain standard deviation and
packing fraction is evidenced ($\sigma_{\mathcal{C}}\sim\phi-\phi_{J}$) and
does not depend on the softness ratio. Only the linearity cut-off increases
with $kappa$. This is explained by the fact that for low $\kappa$ soft
particles deform more and more rapidly, entering faster in a non-linear
regime.

Granular systems composed of purely soft or mixed rigidity particles are
ubiquitous, so we believe that the results obtained in this letter, more
particularly the compaction law and scaling relations, can release blockages
in domains as different as biology, geoscience or in the industry. Using our
experimental set-up, these results could be confronted in the case of non
circular particles or the case mixture of particles with different rheologies.

###### Acknowledgements.

We thank Gilles Camp and Stéphan Devic for their technical help with the
experimental set-up and the particle making process.

## References

  * Vu and Barés [2019] T.-L. Vu and J. Barés, Soft-grain compression: Beyond the jamming point, Physical Review E 100, 042907 (2019).
  * Cantor  _et al._ [2020] D. Cantor, M. Cárdenas-Barrantes, I. Preechawuttipong, M. Renouf, and É. Azéma, Compaction model for highly deformable particle assemblies, Physical Review Letters 124, 208003 (2020).
  * Wyatt  _et al._ [2015] T. P. J. Wyatt, A. R. Harris, M. Lam, Q. Cheng, J. Bellis, A. Dimitracopoulos, A. J. Kabla, G. T. Charras, and B. Baum, Emergence of homeostatic epithelial packing and stress dissipation through divisions oriented along the long cell axis, Proceedings of the National Academy of Sciences 112, 5726 (2015).
  * Mauer  _et al._ [2018] J. Mauer, S. Mendez, L. Lanotte, F. Nicoud, M. Abkarian, G. Gompper, and D. A. Fedosov, Flow-induced transitions of red blood cell shapes under shear, Physical review letters 121, 118103 (2018).
  * Bolton and Weaire [1990] F. Bolton and D. Weaire, Rigidity loss transition in a disordered 2d froth, Physical review letters 65, 3449 (1990).
  * Katgert and van Hecke [2010] G. Katgert and M. van Hecke, Jamming and geometry of two-dimensional foams, EPL (Europhysics Letters) 92, 34002 (2010).
  * Brujić  _et al._ [2003] J. Brujić, S. F. Edwards, D. V. Grinev, I. Hopkinson, D. Brujić, and H. A. Makse, 3d bulk measurements of the force distribution in a compressed emulsion system, Faraday discussions 123, 207 (2003).
  * Zhou  _et al._ [2006] J. Zhou, S. Long, Q. Wang, and A. D. Dinsmore, Measurement of forces inside a three-dimensional pile of frictionless droplets, Science 312, 1631 (2006).
  * Cooper and EATON [1962] J. A. Cooper and L. E. EATON, Compaction behavior of several ceramic powders, Journal of the American ceramic society 45, 97 (1962).
  * Kawakita and Lüdde [1971] K. Kawakita and K.-H. Lüdde, Some considerations on powder compression equations, Powder technology 4, 61 (1971).
  * Liu and Nagel [1998] A. J. Liu and S. R. Nagel, Jamming is not just cool any more, Nature 396, 21 (1998).
  * Hayman  _et al._ [2011] N. W. Hayman, L. Ducloué, K. L. Foco, and K. E. Daniels, Granular controls on periodicity of stick-slip events: kinematics and force-chains in an experimental fault, Pure and applied geophysics 168, 2239 (2011).
  * Barés  _et al._ [2017] J. Barés, D. Wang, D. Wang, T. Bertrand, C. S. O’Hern, and R. P. Behringer, Local and global avalanches in a two-dimensional sheared granular medium, Physical Review E 96, 052902 (2017).
  * Jacquin  _et al._ [2011] H. Jacquin, L. Berthier, and F. Zamponi, Microscopic mean-field theory of the jamming transition, Physical review letters 106, 135702 (2011).
  * van Hecke [2009] M. van Hecke, Jamming of soft particles: geometry, mechanics, scaling and isostaticity, Journal of Physics: Condensed Matter 22, 033101 (2009).
  * Radjai  _et al._ [1996] F. Radjai, M. Jean, J.-J. Moreau, and S. Roux, Force distributions in dense two-dimensional granular systems, Physical review letters 77, 274 (1996).
  * O’hern  _et al._ [2003] C. S. O’hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Jamming at zero temperature and zero applied stress: The epitome of disorder, Physical Review E 68, 011306 (2003).
  * Majmudar  _et al._ [2007] T. S. Majmudar, M. Sperl, S. Luding, and R. P. Behringer, Jamming transition in granular systems, Physical review letters 98, 058001 (2007).
  * Cox  _et al._ [2016] M. Cox, D. Wang, J. Barés, and R. P. Behringer, Self-organized magnetic particles to tune the mechanical behavior of a granular system, EPL (Europhysics Letters) 115, 64003 (2016).
  * Papadopoulos  _et al._ [2016] L. Papadopoulos, J. G. Puckett, K. E. Daniels, and D. S. Bassett, Evolution of network architecture in a granular material under compression, Physical Review E 94, 032908 (2016).
  * Jop  _et al._ [2006] P. Jop, Y. Forterre, and O. Pouliquen, A constitutive law for dense granular flows, Nature 441, 727 (2006).
  * Cárdenas-Barrantes  _et al._ [2020] M. Cárdenas-Barrantes, D. Cantor, J. Barés, M. Renouf, and E. Azéma, Compaction of mixtures of rigid and highly deformable particles: A micromechanical model, Physical Review E 102, 032904 (2020).
  * Cárdenas-Barrantes  _et al._ [2021] M. Cárdenas-Barrantes, D. Cantor, J. Barés, M. Renouf, and E. Azéma, Micromechanical description of the compaction of soft pentagon assemblies, Physical Review E 103, 062902 (2021).
  * Indraratna  _et al._ [2019] B. Indraratna, Y. Qi, T. N. Ngo, C. Rujikiatkamjorn, T. Neville, F. B. Ferreira, and A. Shahkolahi, Use of geogrids and recycled rubber in railroad infrastructure for enhanced performance, Geosciences 9, 30 (2019).
  * Khatami  _et al._ [2020] H. Khatami, A. Deng, and M. Jaksa, The arching effect in rubber–sand mixtures, Geosynthetics International 27, 432 (2020).
  * H.-H.Tsang [2008] H.-H. H.-H.Tsang, Seismic isolation by rubber–soil mixtures for developing countries, Earthquake engineering & structural dynamics 37, 283 (2008).
  * Tsiavos  _et al._ [2019] A. Tsiavos, N. A. Alexander, A. Diambra, E. Ibraim, P. J. Vardanega, A. Gonzalez-Buelga, and A. Sextos, A sand-rubber deformable granular layer as a low-cost seismic isolation strategy in developing countries: Experimental investigation, Soil Dynamics and Earthquake Engineering 125, 105731 (2019).
  * Anastasiadis  _et al._ [2012] A. Anastasiadis, K. Senetakis, and K. Pitilakis, Small-strain shear modulus and damping ratio of sand-rubber and gravel-rubber mixtures, Geotechnical and Geological Engineering 30, 363 (2012).
  * Mashiri  _et al._ [2015] M. S. Mashiri, J. S. Vinod, M. N. Sheikh, and H. H. Tsang, Shear strength and dilatancy behaviour of sand–tyre chip mixtures, Soils and Foundations 55, 517 (2015).
  * Barés  _et al._ [2021] J. Barés, M. Cárdenas-Barrantes, D. Cantor, É. Azéma, and M. Renouf, Highly strained mixtures of bidimensional soft and rigid grains: an experimental approach from the local scale, in  _EPJ Web of Conferences_ , Vol. 249 (EDP Sciences, 2021) p. 05004.
  * [31] Canon™ canoscan 9000f mark ii.
  * [32] Moldstar™ 15 slow, castable silicone from smooth-on, https://www.smooth-on.com/products/mold-star-15-slow/.
  * Vu  _et al._ [2019a] T.-L. Vu, J. Barés, S. Mora, and S. Nezamabadi, Numerical simulations of the compaction of assemblies of rubberlike particles: A quantitative comparison with experiments, Physical Review E 99, 062903 (2019a).
  * Vu  _et al._ [2019b] T.-L. Vu, J. Barés, S. Mora, and S. Nezamabadi, Deformation field in diametrically loaded soft cylinders, Experimental Mechanics 59, 453 (2019b).
  * [35] Python codes to post-process experimental images (particle tracking, dic, fields computation, etc.), https://git-xen.lmgc.univ-montp2.fr/Image/granular_dic.
  * Taber [2004] L. A. Taber,  _Nonlinear theory of elasticity: applications in biomechanics_ (World Scientific, 2004).
  * Andreotti  _et al._ [2013] B. Andreotti, Y. Forterre, and O. Pouliquen,  _Granular media: between fluid and solid_ (Cambridge University Press, 2013).
  * Nezamabadi  _et al._ [2019] S. Nezamabadi, X. Frank, J.-Y. Delenne, J. Averseng, and F. Radjai, Parallel implicit contact algorithm for soft particle systems, Computer Physics Communications 237, 17 (2019).
  * Estrada  _et al._ [2008] N. Estrada, A. Taboada, and F. Radjai, Shear strength and force transmission in granular media with rolling resistance, Physical Review E 78, 021301 (2008).
  * Skinner [1969] A. E. Skinner, A note on the influence of interparticle friction on the shearing strength of a random assembly of spherical particles, Geotechnique 19, 150 (1969).
  * Kim and Carroll [1987] K. T. Kim and M. M. Carroll, Compaction equations for strain hardening porous materials, International journal of plasticity 3, 63 (1987).


# N2, NO and O2 molecules in LiGaO2 in both Ga and Li sites and their relation
to the vacancies.

Klichchupong Dabsamut  Department of Physics, Faculty of Science, Kasetsart
University, Bangkok 10900 Thailand  Department of Physics, Case Western
Reserve University, 10900 Euclid Avenue, Cleveland, Ohio 44106-7079, USA
Adisak Boonchun<EMAIL_ADDRESS>Department of Physics, Faculty of Science,
Kasetsart University, Bangkok 10900 Thailand  Walter R. L. Lambrecht
<EMAIL_ADDRESS>Department of Physics, Case Western Reserve
University, 10900 Euclid Avenue, Cleveland, Ohio 44106-7079, USA

###### Abstract

Doping of the ultrawide band gap semicodnctor LiGaO2 ($E_{g}=5.6$ eV) with N2,
NO and O2 molecules placed in either Ga or Li-vacancies is studied using
first-principles calculations. These molecular dopants are considered as
potential acceptors for $p$-type doping. Their optimal placement and
orientation relative to the lattice is studied as well as their transition
levels and energy of formation. Unfortunately, they are found to have deep
acceptor level transition states. We discuss the relation of the transition
levels to the one-electron levels, their spin state and defect wave functions.
They are found to be closely related to those of the corresponding vacancies.

##  I Introduction

LiGaO2 was recently proposed to be a promising ultrawide band gap
semiconductor. Its band gap is 5.6 eV as verified by optical absorption
measurements [1, 2, 3] and QS$GW$ (quasiparticle self-consistent $GW$)
calculatons.[4, 5] Similar to $\beta$-Ga2O3 density functional theory (DFT)
calculations calculations using a hybrid functional HSE06 predict that this
material could be $n$-type doped by Si or Ge, which form shallow donors.[6, 7]
On the other hand, comparing to $\beta$-Ga2O3, it has a larger gap, somewhat
lower valence band masses and a simpler wurtzite related crystal structure,
the $Pna2_{1}$ structure,[8] in which all atoms are tetraehedrally
coordinated. It can be grown in bulk form.[9]

In the only study thus far of the doping opportunities,[7] $p$-type doping was
still found to be problematic, as is also the case for its parent compound ZnO
and for many other oxides, including $\beta$-Ga2O3. Standard candidate
acceptor dopants like NO were found to give too deep levels in the gap and in
fact in this case amphoteric character. Zn in LiGaO2 on the other hand suffers
from site competition between the acceptor ZnGa and donor ZnLi. Thus, more
involved schemes need to be pursued for establishing $p$-type doping. Here we
take some inspiration from previous work in ZnO. [10] In ZnO, it is known that
a shallow acceptor exists[11] that is related to N but is not the simple
substitutional defect. Several candidates have been suggested. Lautenschlaeger
$et$ $al.$[12] proposed that H is involved in a
$\mathrm{N}_{\mathrm{O}}-\mathrm{H}-\mathrm{N}_{\mathrm{O}}$ complex. Bang et
al. [13] proposed doping with NH3 and a NH3 on the Zn-site. Boonchun and
Lambrecht proposed a N2 molecule on a Zn-site,[10]. Another model is a
N${}_{\mathrm{O}}-V_{\mathrm{Zn}}$ complex.[14]

The latter was shown to be possible to incorporate making use of the surface
chemistry of adsorption of N on a ZnO O-terminated Zn-polar surface. In that
case, N incorporates preferentially on the Zn sites but then leads to a
O-vacancy in the next layer. After this N${}_{\mathrm{Zn}}-V_{\mathrm{O}}$
complex is burried inside the sample, the N can easily overcome a barrier and
interchange places and hop into the adjacent $V_{\mathrm{O}}$ thereby creating
the desired complex, which was found to be shallow. This scheme was made to
work in practice by a careful sequence of growth and annealing treatments by
Reynolds et al.[15, 16] and is perhaps one of the most promising routes to
$p$-type doping in ZnO.

The alternative proposed scheme with N2 has not yet been pursued
experimentally. It is based on the idea of aligning the molecular levels of N2
with those in ZnO. A N2 molecule has 10 valence electrons (2$s$2$p$) in its
neutral charge state and is extremely stable because all the $\sigma$ and
$\pi$ bonding states are filled. When placed in a Zn site, it needs to act as
a $2+$ ion and thus compared with Zn including the filled $3d$-band which lies
at about the same energy as the deeper N2 molecule $\sigma$-bonds, it lacks 2
electrons compared to Zn2+. This implies it is a double acceptor. In the
single negative charge state of the defect, (N2)${}_{\mathrm{Zn}}^{-1}$ state
it would correspond to a N${}_{2}^{+1}$ molecule. This state with an unpaired
spin has been observed clearly in Electron Paramagnetic Resonance (EPR) [17]
although these authors assumed it was located on a O-site. However, Boonchun
and Lambrecht[10] showed that on the O-site it would behave as a donor and
demonstrated that the Zn-site has qualitatively the correct $g$-tensor and
hyperfine splitting. While they proposed it is a shallow acceptor and thus
suitable for $p$-type doping, that depends somewhat on the choice of
functional and later calculations by Petretto and Bruneval[18] found a
somewhat different location of the N2 molecule more strongly bonded to the
surrounding O of the $V_{\mathrm{Zn}}$ which behaved as a deep acceptor. An
experimental study of the recharging behavior of the N2 in ZnO also proposed
that it is a deep center.[19] Nonetheless, the verdict on the shallow or deep
nature of this defect center is still out and there are alternative
explanations for the observed behavior. For example, there might be both a
deep and a shallow configuration of the defect depending on how strongly the
N2 bonds to the O surrounding it. Rotational entropy might favor the less
strongly bound form at higher temperatures and the charging/recharging
behavior of the defect taking into account the higher charge states of the
defect has not yet been studied in detail by computational means.

Here we extend the idea of N2 incorporation on a cation site to LiGaO2. There
are some differences with the corresponding system in ZnO. In LiGaO2, the N2
on a Ga site would be a triple acceptor instead of a double acceptor. The
Ga-$3d$ level lies much deeper than the Zn-$3d$ so the alignment with the N2
molecular levels might be different. The $V_{\mathrm{Ga}}$ is a defect of high
energy of formation, so it might be less easy to obtain N2 in this site, but
the $V_{\mathrm{Li}}$ has low energy of formation. The alignment of the
molecular levels with the host levels is not a-priori clear and is
investigated here. Here we study N2, NO and O2 as molecular doping species on
both cation vacancy sites and evaluate their formation energy, structural
optimization for different orientations of the molecule and their transition
levels.

##  II Computational Methods

Our study is based on density functional calculations (DFT) using the Heyd-
Scuseria-Ernzerhof (HSE) hybrid functional[20, 21] with the standard fraction
of exact screened Hartree-Fock type exchange $\alpha$ = 0.25 and range
parameter of $\mu$ = 10 Å, giving a gap of $E_{g}$ = 5.1 eV slightly smaller
than the experimental value.[22, 23, 24, 25] Since this work extends our
previous published results, we keep the computational approach used here the
same as in [6] and [7]. The calculations are performed within the Vienna Ab-
Initio Simulation Package (VASP).[26, 27] The electron ion interactions are
described by means of the Projector Augmented Wave (PAW) method.[28, 29] We
use a well-converged energy cutoff of 500 eV for the projector augmented plane
waves. The dopants are modeled in a 128-atom supercell, for which a shifted
single k-point Brillouin zone sampling was found to be sufficiently converged.

The analysis of the defect levels follows the standard defect approach as
outlined in e.g. Freysoldt et al. [30] The energy of formation of the defect
$D^{q}$ in charge state $q$ is given by

$\begin{split}E_{f}(D^{q})=E_{tot}(C:D^{q})-E_{tot}(C)-\sum_{i}\Delta
n_{i}\mu_{i}\\\ +q(\epsilon_{v}+\epsilon_{F}+V_{align})+E_{cor}\end{split}$
(1)

where $E_{tot}(C:D^{q})$ is the total energy of the supercell containing the
defect and $E_{tot}(C)$ is the total energy of the perfect crystal supercell.
The energy for adding or removing atoms from the crystal to a reservoir in the
process of producing the defect is represented by the chemical potentials
$\mu_{i}$. The chemical potentials of the host elements that we chose are
presented in our previous paper.[6] The $\Delta n_{i}$ value represents the
change in the number of atoms in the species $i$. For the impurities, the
chemical potential of N is calculated from the total energy per atom of an
isolated N2 molecule, where the chemical potential of O is fixed which respect
to the extreme condition as reported in the previous study.[6] The alignment
and image charge corrections are calculated using the Freysoldt approach.[31]
These energies of formation are used to study the transition levels between
different charge states, each with fully relaxed structures but kept at the
same volume.

The formation energies of $\rm(N_{2})_{Ga}$, $\rm(O_{2})_{Ga}$ and
$\rm(NO)_{Ga}$ are plotted as function of Fermi level $\varepsilon_{F}$ in
Fig.2 for chemical potential conditions E and F. Note that chemical potential
conditions E and F (see figure 2 in Ref. 6) correspond respectively to more
Li-rich and more Ga-rich conditions but restricted by the formation of
competing binary compounds. Realistic conditions on the O chemical potential
corresponding to the typical growth temperature are used in both cases. The
absolute chemical potentials at point E are $\mu_{O}=-8.65$ eV,
$\mu_{Li}=-4.01$ eV and $\mu_{Ga}=-6.58$ eV, while the absolute chemical
potentials at point F are $\mu_{O}=-8.65$ eV, $\mu_{Li}=-6.50$ eV and
$\mu_{Ga}=-4.10$ eV.

The absolute chemical potentials used for the dopants are $\mu_{N}=-10.24$ eV.
However, the formation energies of all three impurities are high with this
choice of chemical potentials E and F, implying low solubility in LiGaO2.

##  III Results

###  III.1 Structural models and relaxation.

For each molecule we study various configurations of the relative orientation
of the molecule relative to the crystal and fully relax them. For the N2 and
O2 cases, we have first positioned N2 and O2 along the original bond of Ga and
O. Both of $\rm(N_{2})_{Ga}$ and $\rm(O_{2})_{Ga}$, they can be categorized
into planar configurations (P1, P2 and P3) and vertical configuration (V1) as
shown in Fig.S1 in supplementary information (SI). Similarly to the NO case,
the initial position of NO was set along the original Ga-O bond. In this case,
we have 6 initial planar structures (planar (P1-a, P1-b, P2-a, P2-b, P3-a, and
P3-b) and 2 vertical structures (V1-a and V1-b) as shown in Fig.S2 in SI. The
reason why there are more structural models is that for the $\rm(NO)_{Ga}$
case, we need to distinguish how the dipole of the the NO molecule is
directed.

Figure 1: The optimized structure of (a) top view of $\rm(N_{2})_{Ga}$, (b)
side view of $\rm(N_{2})_{Ga}$, (c) top view of $\rm(O_{2})_{Ga}$, (d) side
view of $\rm(O_{2})_{Ga}$, (e) top view of $\rm(NO)_{Ga}$ and (f) side view of
$\rm(NO)_{Ga}$. We note that we removed other layers of LiGaO2 to make it
easier to spot the molecules. Yellow, brown and red balls represent Li, Ga and
O atoms, respectively.

Next, we found that the most energetically favourable is the P3 arrangement in
the case of N2 and O2. The optimized structure of $\rm(N_{2})_{Ga}$ and
$\rm(O_{2})_{Ga}$ are presented in Fig.1(a)-(d). In the top view, we can see
only a slight shift of the positions, but in the side view we may notice a
larger shift of the N2 or O2 molecule from their initial positions. For the
NO-case, the most energetically favorable is the P2-a arrangement. The
optimized structure is presented in Fig.1(e)-(f). Compared to the N2 and O2
cases, we here see a larger distortion of the O atoms surrounding the vacancy
both in the top and side views. One of the O below has significantly moved
closer to the N atom of the NO to the expense of weakening its bond with the
next nearest neighbor Ga.

Figure 2: The formation energies of N2, NO and O2 in (a) Ga site for point E,
(b) Ga site for point F, (c) Li site for point E and (d) Li site for point F.

###  III.2 Energies of formation and transition levels

The energies of formation for different charge states, choosing the
configuration with lowest structural energy are given in Fig. 2 for two cases
of chemical potentials, more Li-rich (E) and more Ga-rich (F) cases. Both
points are chosen to pick a realistic O chemical potential corresponding to
the growth conditions as explained in Ref. 6.

They show first that the energies of formation on the Ga site (Fig.2(a)-(b))
are relatively high. The reason for this is simply that the Ga vacancy itself
is a high-energy of formation defect as reported in our previous study.[6] One
could in principle still realize these defect cases by first creating
vacancies using for example high-energy irradiation and subsequently introduce
the dopants by reactions from the surface inwards, by exposure to N2, NO or O2
gas.

Second, they show that N2 occurs in three charge states depending on the
position of the Fermi level, neutral, $-1$ and $-3$. This indicates a negative
$U$ behavior between $-1$ and $-3$ charge states. Interestingly, in the
corresponding Ga-vacancy, the $-2$ charge state also has already a narrow
energy range of stability. Likewise, NO occurs in neutral and $-2$ charge
states only. However, O2 occurs only in neutral and $-1$ charge states and
exhibits no negative-$U$ behavior.

Third, the transition levels all are quite deep with the lowest $0/-1$
transition level occurring at about 2.85 eV above the valence band maximum
(VBM) for the N2 case.

The energy of formation results for the various molecules on the Li site are
given in Fig.2(c)-(d). These have slightly lower energies of formation
consistent with the lower energy of formation of the $V_{\mathrm{Li}}$ as
reported in our previous study.[6] Structurally, there is ample room for these
molecules in both types of vacancy sites. The transition levels are slightly
shallower for the Li site occupation. For example, the $0/-1$ transition level
for (N2)Li is 1.57 eV instead of (N2)Ga where it is 2.85 eV. The (N2)Li exists
only in 0 and $-1$ states and the same holds for the O2 case, while the NO
molecule exists in three charge states, 0, $-1$, $-2$.

###  III.3 One-electron levels and defect wave functions.

Figure 3: PDOS of N2 in the Ga-vacancy in the neutral charge state of the
defect: (a) isolated N2, (b) interacting system, and (c) Ga-vacancy. We note
that the dash lines indicate the highest occupied state.  Figure 4: Net spin
density of (N2)Ga in the neutral charge state.

In this section, we examine the one-electron levels, their defect wave
functions, and partial densities of states (PDOS) to clarify the spin states
of the defects and their relation to the transition levels.

In Fig. 3 we show the PDOS of N2 placed in the Ga-vacancy in panel (b), while
panels (a) and (c) show the isolated molecule and the Ga-vacancy PDOS
respectively. In panel (a), we have labeled the peaks according to the
standared nomenclature for molecular eigenstates. First, we recognize the
lowest empty molecular orbital (LUMO) $\pi_{g}$ antibonding state lying above
10 eV. Next, the $\sigma_{g}^{+}$ state forms the highest occupied molecular
orbital (HOMO). We see that it consists mostly of N-$2p$ orbitals with a small
N-$2s$ contribution. The HOMO-1 level is the $\pi_{u}$ bonding state, which
has pure $p$ character and is doubly degenerate and thus has twice the peak
height of the HOMO. The next level is the $\sigma_{u}^{-}$ state which is a
bonding state between N-$2s$ and N-$2p_{z}$. Its antibonding counter part lies
at higher energy above the CBM and is not shown. The $\sigma_{g}^{-}$ state at
$-20$ eV is the bonding counterpart to the $\sigma_{g}^{+}$ state and has
mostly N-$2s$ character with a little bit of N-$2p$. Here, the labels $g,u$
refer to even or odd with respect to inversion and $\pm$ refers to
(anti)bonding character. We need this extra label only for $\sigma$ states
because there is more than one $\sigma_{g}$ and $\sigma_{u}$. From panel (b)
we can see that there is a net spin-polarization in the neutral charge state,
which is related to its odd number of electrons. We can still clearly
recognize the $\pi_{g}$ antibonding state of the N2 molecule, lying at the
conduction band minimum (CBM). A minority spin state is seen at mid gap and
states of both majority and minority spin character occur just above the VBM.
We here aligned the deep O-$2s$ states with those of the host system
containing the Ga-vacancy shown in panel (c) and indicate the Fermi level at
zero temperature (highest occupied state) by the dashed line. We can see
several defect levels of minority spin also in panel (c) for the neutral Ga-
vacancy with only the lowest one of these occupied. These levels look dangling
bond like (see Fig.S6 in SI) but are spread over several neighboring O sites.

The N-$2p$ projected PDOS of the N2 molecule inside the Ga-vacancy shows a
peak about 1 eV below the VBM with a net difference between spin-up and spin
down, indicating that part of the spin density must be localized on the
molecule. However, the bonding $\pi_{u}$ state of the molecule are more
perturbed by the interaction with the host. It appears that this state forms
both bonding and antibonding states with the O-$2p$ valence states with the
bonding states lying just below the O-$2p$ bands and the antibonding ones
slightly below the VBM. In fact, we find three peaks localized on the molecule
lying just below the O-$2p$ valence band, two with predominantly N-$2p$
character and one with predominantly N-$2s$ character. The $\pi_{u}$ state in
the molecule has pure $p$ character but when placed in the Ga-vacancy it
acquires a small amount of $s$ character and the peak height is now the same
as for the $\sigma_{g}^{+}$ state because its weight is also partially found
in the states just below the VBM and distributed to some extent throughout the
valence band. The states near $-14$ eV are seen to be the Ga-$3d$ states and
some hybridization between the deep lying $\sigma_{g}^{-}$ states with the
O-$2s$ band can also be seen. All of this indicates a significant interaction
between the molecular states and the vacancy states. The net spin density in
the neutral charge state of (N2)Ga is shown in Fig. 4 and confirms this
picture.

Figure 5: The partial density of state (PDOS) and the partial charge density
of (N2)Ga. We note that the dash line separates the occupied and unoccupied
state.

In Fig.5, we show the partial densities of states (PDOS) and defect
wavefunctions modulo squared for several energy ranges corresponding to defect
levels for the N2 on Ga site. We note several defect levels or resonances in
the bands introduced by the molecule insertion into the vacancy. First, the
state (labeled II) corresponds clearly to the molecular $\pi^{*}$ antibonding
state and lies just below the CBM of the host crystal. A resonance also occurs
just below the VBM and corresponds to the highest occupied molecular orbital
(HOMO) of the N2 molecule and has $\sigma$-bonding character. In addition to
these molecular levels, we recognize a defect level (labeled I) which is
localized on one of the O surrounding the defects. It is accompanied by a
polaronic distortion localizing the wave function on only one of the O nearest
neighbors of the Ga vacancy in which the N2 is inserted. A defect level closer
to the VBM, (labeled III) still shows also localization on one of the nearest
neighbor oxygens but spreads somewhat further to further neighbor O. All of
these were obtained for the neutral defect.

For the N2 on the Ga site, which has 3 holes in the defect or molecule levels
compared to the perfect crystal since it can exist also in -1, -2 and -3
charge states by successively adding electrons, we have also examined the
possibility of a higher spin $S=3/2$ state with three parallel spins. This
configuration, however, was found to have 1.53 eV higher energy and can thus
safely be excluded from further consideration. We find from the corresponding
PDOS that the N-$p$ related states occur deeper below the VBM, as shown in
Fig.S6 in SI.

Figure 6: PDOS of N2 on Li-site in neutral charge state. (a) isolated
molecule, (b) molecule in the Li-site and (c) Li-vacancy. We remark that the
dash lines separate the occupied and unoccupied states.

For N2 placed on the Li-site we similarly show the comparison between the N2
PDOS in the isolated molecule and in the Li site and with the Li-vacancy case,
again, all for the neutral charge state. The $\pi_{g}$ state is again most
readily identified as lying near the CBM at about 6 eV. Defect states related
to the Li-vacancy are also found in the gap. However, in this case the N2 HOMO
state is located near the bottom of the O-$2p$ states and much less
interaction is observed than for the Ga-vacancy. This indicates a different
electrostatic potential in which the N2 molecule finds itself. In fact, it is
a less repulsive potential because removing a Li+ is less repulsive than
removing a Ga3+, thereby placing these states deeper relative to the O-$2p$
bands. The molecular states are thus differently aligned with respect to the
LiGaO2 host bands when the N2 is located in the Ga than in the Li vacancy. A
weaker interaction is also observed for the deep $\sigma_{g}^{-}$ state with
the O-$2s$ band. There is is still a net spin-density as shown in Fig. 7.

Figure 7: Net spin density for (N2)Li in the neutral charge state.

Until now we have focused on the neutral charge state and on the N2 molecule
related states and their interaction with their surroundings. However, it is
also clear that several states exist in the gap and inspection of their wave
functions shows that these states are primarily localized on a single oxygen
and are closely related to the corresponding vacancy states. Adding an
electron in the negative charge states makes the system non-spin polarized as
there is now an even number of electrons. We checked that the high-spin $S=1$
configuration for the $q=-1$ charge state has higher energy by 1.05 eV. Again,
we find that the occupied one-electron states of the N2 molecule occur farther
below the VBM. We can study the addition of an electron in two steps, first
keeping the structure frozen in that of the neutral charge state and
subsequently relaxing the atomic positions. We find that the levels in the gap
are spread deeper into the conduction band in the frozen geometry of the
neutral state and move closer to the VBM upon relaxation. This is shown in
Fig.S8 in the SI. Likewise in the $-2$ charge state, we have again a net spin
polarization and so on. After filling the empty states of the molecule, the
next electrons essentially are located in the vacancy related dangling bond
states and correspond closely to the corresponding vacancy case. However for
the $V_{\mathrm{Ga}}$ we find several states with localization on more than
one O-neighbor to the vacancy.

Figure 8: Energy levels in the supercell as function of charge state at the
shifted ${\bf k}$-point, labeled by the eigenvalue number. Solid lines connect
spin up states and dashed lines spin down states.

The one-electron energy levels lying inside the band gap are shown in Fig.8 as
function of charge state. They are shown for both spin states. One can see
that in the neutral charge state only one defect level (No. 415 $\uparrow$)
close to the VBM (levels 414 for $\uparrow$ and $\downarrow$ spin) is occupied
while three remain empty. In the $q=-1$ charge state, the lower one of these,
which is the minority spin level 415, is lowered significantly while it
becomes occupied. Likewise, in the $q=-2$ and $q=-3$ charge states, additional
levels come down in energy as they become occupied. This reflects the fact
that in hybrid functional, due to the partial inclusion of exact exchange, the
potential is orbital-dependent and occupied levels shift down relative to
empty levels.

We should here caution the reader that these one-electron levels, also seen as
peaks in the PDOS, reflect the band positions in the 128 cell model, which is
not truly reflective of the dilute limit. An electrostatic correction
$-2q\alpha/\varepsilon L$ applies to the single particle levels for charged
defect states due to their interaction with the image charges and the
neutralizing background, similar to the well-known corrections to the energies
of formation. These corrections allow one to extrapolate to the dilute limit
and would shift the one-electron levels of the isolated defect deeper into the
gap with a shift proportional to their charge and inversely proportional the
size of the supercell $L$. Here $\alpha$ is a Madelung constant and
$\varepsilon$ is the dielectric constant of the host. Because these defect
one-electron levels are pushed deeper in the gap than indicated by the PDOS,
filling them only occurs when the Fermi level is also rather deep in the gap
and this essentially explains the deep nature of the transition levels even
though rather shallow states are seen in the PDOS figures or for the occupied
one-electron levels.

We can further examine the PDOS of these different charge states and the
associated wave functions. These details are given in Fig.S9 in the SI for
charge states -1,-2, -3. The main conclusion is that the peaks in the PDOS can
be associated with the levels given in Fig. 8 and their associated states show
mostly vacancy type localized states with in some cases a partial contribution
on the molecule and for states closer to the VBM a wider spread farther away
from the defect site. We thus find that the states in the gap are closely
related to those of the corresponding vacancy and have generally the character
of dangling bond states on the surrounding oxygen atoms. The relaxation tends
to localize them on a single oxygen rather than being distributed over several
oxygens. This indicates the polaronic nature of these defect states.

Moving on to the NO and O2 as inserted molecules, the changes can basically be
understood in terms of the added valence electron in each step. The PDOS of
the NO and O2 on the Ga and Li site in the neutral charge state are shown in
SI. They indicate that the NO molecule in the neutral charge state is similar
to the N2 in the $-1$ state and the O2 molecule in the neutral state is
similar to the N2 in the $-2$ state. The defect levels in the gap in all these
cases are mostly localized on a single oxygen neighbor of the vacancy in which
the molecule is placed.

##  IV Discussion

In the present calculations, we find deep levels and polaronic acceptor
behavior for each of these molecules when inserted in either a Li or Ga
vacancy. This is manifested in the localization of the defect wave functions
each time on a single O. In part, this behavior is emphasized by the hybrid
functional choice. All of these results were obtained with the standard HSE
setting of 25 % exact exchange. As a final check on our results, we therefore
investigated the satisfaction of the generalized Koopmans theorem (gKT) by
calculating the non-Koopman’s energy[32],

$\begin{split}E_{NK}=\varepsilon_{N}-E_{I}=\varepsilon_{N}-(E_{N}-E_{N-1})\end{split}$
(2)

where $\varepsilon_{N}$ is a localized state’s Kohn-Sham quasiparticle energy
in an N electron system, $E_{I}$ is the system’s ionization energy, which is
equivalent to the difference between the total energies $E_{N}$ and $E_{N-1}$
of the N and N-1 electron systems with defect, respectively. In the Table 1,
we calculated the non-Koopman’s energy ($E_{NK}$) for 20%, 25% and 30% HF
mixing parameters from our previous SnGa calculation.[7] In the case of the
gKT, $E_{NK}$ should equal be zero. As a result, we found that $E_{NK}$ for
30% HF mixing parameter is almost zero (-0.012 eV), indicating that gKT is
fulfilled. Moreover, the energy gap from this calculation is very close to the
experiment value as presented in Table 1. Although this indicates that a 30 %
HF would have been preferable, we have kept the 25 % results to be consistent
with our previous papers. Even for the 25% HF mixing parameter, $E_{NK}$ is
already relatively small (-0.121 eV), thus the self-interaction error may be
considered small in all the present calculations.

Table 1: The calculated energy gap ($E_{g}$), the localized state’s Kohn-Sham quasiparticle energy ($\varepsilon_{N}$), the system’s ionization energy ($E_{I}$) and the non-Koopman’s energy ($E_{NK}$) for 20%, 25% and 30% HF mixing parameters in SnGa system. | HF mixing parameter |
---|---|---
| 20% | 25% | 30% | Exp. value
$E_{g}$ (eV) | 4.70 | 5.10 | 5.50 | 5.3-5.6[22, 23, 24, 25]
$\varepsilon_{N}$ (eV) | 3.818 | 4.010 | 4.070 |
$E_{I}$ (eV) | 4.057 | 4.132 | 4.082 |
$E_{NK}$ (eV) | -0.240 | -0.121 | -0.012 |

##  V Conclusions

In this paper we have examined the behavior of inserted diatomic molecules N2,
NO and O2 in both the Li and Ga site. We found that these systems have high
energy of formation, especially for the Ga site, which is related to the
corresponding vacancy formation energy. Even though these systems all behave
as acceptor systems, we found deep transition levels in the gap, which are not
compatible with $p$-type doping. Analysis of the PDOS and one-electron levels
shows that as successive levels become occupied in the negative charge states,
they move closer to the VBM. This is a direct result of the hybrid functional
used which lowers occupied vs.  empty levels. While the occupied ones appear
much shallower than the transition levels, this is an artifact of the
supercell size. Corrections for the one-electron levels in charged systems
will tend to move these levels deeper in the gap proportional to their charge
state, so that the transition levels and one-electron levels become consistent
with each other in labeling these systems as deep acceptors. The behavior of
NO and O2 in charge state $q$ correspond to the corresponding $q-1$, $q-2$
states of N2. We find a stronger interaction of the N2 molecular states with
the O-states of the host for the Ga than for the Li position in terms of the
occupied states of the molecule. The defect levels in the gap are closely
related to the corresponding vacancy levels and are dangling bond like and
characterized by polaronic distortions localizing the defect wave functions on
a single O neighbor to the vacancy. We find that the Koopmans theorem is
already fairly well satisfied by the 25 % exact exchange hybrid functional
used here and in our previous work on native defects, although in future work
a 30 % model would both improve the band gap and the satisfaction of Koopmans
theorem.

###### Acknowledgements.

The work at CWRU was supported by the U.S. National Science Foundation under
grant No. 1755479. K.D. was supported by the National Research Council of
Thailand (NRCT), No. NRCT5-RGJ63002-028. A.B. has been funded by the Office of
National Higher Education Science Research and Innovation Policy Council
(NXPO), Thailand, through Program Management Unit for Competitiveness (PMU C),
Contract Number C10F630073. Calculations were in part performed at the Ohio
Supercomputer Center.

## References

  * Wolan and Hoflund [1998a] J. T. Wolan and G. B. Hoflund, “Chemical alteration of the native oxide layer on LiGaO2(001) by exposure to hyperthermal atomic hydrogen,” J. Vac. Sci. Tech. A 16, 3414–3419 (1998a).
  * Johnson, McLeod, and Moewes [2011a] N. W. Johnson, J. A. McLeod, and A. Moewes, “The electronic structure of lithium metagallate,” Journal of Physics: Condensed Matter 23, 445501 (2011a).
  * Chen  _et al._ [2014a] C. Chen, C.-A. Li, S.-H. Yu, and M. M. Chou, “Growth and characterization of $\beta$-LiGaO2 single crystal,” Journal of Crystal Growth 402, 325 – 329 (2014a).
  * Radha, Ratnaparkhe, and Lambrecht [2021] S. K. Radha, A. Ratnaparkhe, and W. R. L. Lambrecht, “Quasiparticle self-consistent $GW$ band structures and high-pressure phase transitions of $\mathrm{Li}\mathrm{Ga}{\mathrm{O}}_{2}$ and $\mathrm{Na}\mathrm{Ga}{\mathrm{O}}_{2}$,” Phys. Rev. B 103, 045201 (2021).
  * Boonchun and Lambrecht [2011] A. Boonchun and W. R. L. Lambrecht, “Electronic structure, doping, and lattice dynamics of LiGaO2,” in  _Oxide-based Materials and Devices II_, Proceedings of SPIE, Vol. 7940, edited by F. H. Terani, D. C. Look, and D. J. Rogers, International Society for Optics and Photonics (SPIE, 2011) pp. 129–134.
  * Boonchun, Dabsamut, and Lambrecht [2019] A. Boonchun, K. Dabsamut, and W. R. L. Lambrecht, “First-principles study of point defects in LiGaO2,” Journal of Applied Physics 126, 155703 (2019).
  * Dabsamut, Boonchun, and Lambrecht [2020] K. Dabsamut, A. Boonchun, and W. R. L. Lambrecht, “First-principles study of n- and p-type doping opportunities in LiGaO2,” Journal of Physics D: Applied Physics (2020).
  * Marezio [1965] M. Marezio, “The crystal structure of LiGaO2,” Acta Crystallographica 18, 481–484 (1965).
  * Ishii, Tazoh, and Miyazawa [1998] T. Ishii, Y. Tazoh, and S. Miyazawa, “Single-crystal growth of LiGaO2 for a substrate of GaN thin films,” J. Crystal Growth 186, 409 – 419 (1998).
  * Lambrecht and Boonchun [2013] W. R. L. Lambrecht and A. Boonchun, “Identification of a N-related shallow acceptor and electron paramagnetic resonance center in ZnO: N${}_{2}{}^{+}$ on the Zn site,” Phys. Rev. B 87, 195207 (2013).
  * Zeuner  _et al._ [2002] A. Zeuner, H. Alves, D. Hofmann, B. Meyer, A. Hoffmann, U. Haboeck, M. Strassburg, and M. Dworzak, “Optical Properties of the Nitrogen Acceptor in Epitaxial ZnO,” physica status solidi (b) 234, R7–R9 (2002).
  * Lautenschlaeger  _et al._ [2011] S. Lautenschlaeger, M. Hofmann, S. Eisermann, G. Haas, M. Pinnisch, A. Laufer, and B. K. Meyer, “A model for acceptor doping in ZnO based on nitrogen pair formation,” physica status solidi (b) 248, 1217–1221 (2011).
  * Bang  _et al._ [2015] J. Bang, Y.-Y. Sun, D. West, B. K. Meyer, and S. Zhang, “Molecular doping of ZnO by ammonia: a possible shallow acceptor,” J. Mater. Chem. C 3, 339–344 (2015).
  * Liu  _et al._ [2012] L. Liu, J. Xu, D. Wang, M. Jiang, S. Wang, B. Li, Z. Zhang, D. Zhao, C.-X. Shan, B. Yao, and D. Z. Shen, “$p$-Type Conductivity in N-Doped ZnO: The Role of the ${\mathrm{N}}_{\mathrm{Zn}}\mathrm{\text{$-$}}{V}_{\mathrm{O}}$ Complex,” Phys. Rev. Lett. 108, 215501 (2012).
  * Reynolds  _et al._ [2013] J. G. Reynolds, C. L. Reynolds, A. Mohanta, J. F. Muth, J. E. Rowe, H. O. Everitt, and D. E. Aspnes, “Shallow acceptor complexes in p-type ZnO,” Applied Physics Letters 102, 152114 (2013).
  * Reynolds and Reynolds [2014] J. G. Reynolds and C. L. Reynolds, “Progress in ZnO Acceptor Doping: What Is the Best Strategy?” Advances in Condensed Matter Physics 2014, 457058 (2014).
  * Garces  _et al._ [2003] N. Y. Garces, L. Wang, N. C. Giles, L. E. Halliburton, G. Cantwell, and D. B. Eason, “Molecular nitrogen (N2-) acceptors and isolated nitrogen (N-) acceptors in ZnO crystals,” Journal of Applied Physics 94, 519–524 (2003).
  * Petretto and Bruneval [2014] G. Petretto and F. Bruneval, “Comprehensive Ab Initio Study of Doping in Bulk ZnO with Group-V Elements,” Phys. Rev. Applied 1, 024005 (2014).
  * Philipps  _et al._ [2014] J. M. Philipps, J. E. Stehr, I. Buyanova, M. C. Tarun, M. D. McCluskey, B. K. Meyer, and D. M. Hofmann, “Recharging behavior of nitrogen-centers in ZnO,” Journal of Applied Physics 116, 063701 (2014).
  * Heyd, Scuseria, and Ernzerhof [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, “Hybrid functionals based on a screened Coulomb potential,” J. Chem. Phys. 118, 8207–8215 (2003).
  * Heyd, Scuseria, and Ernzerhof [2006] J. Heyd, G. E. Scuseria, and M. Ernzerhof, “Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)],” J. Chem. Phys. 124, 219906 (2006).
  * Wolan and Hoflund [1998b] J. T. Wolan and G. B. Hoflund, “Chemical alteration of the native oxide layer on ligao 2 (001) by exposure to hyperthermal atomic hydrogen,” Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films 16, 3414–3419 (1998b).
  * Johnson, McLeod, and Moewes [2011b] N. Johnson, J. McLeod, and A. Moewes, “The electronic structure of lithium metagallate,” Journal of Physics: Condensed Matter 23, 445501 (2011b).
  * Ohkubo  _et al._ [2002] I. Ohkubo, C. Hirose, K. Tamura, J. Nishii, H. Saito, H. Koinuma, P. Ahemt, T. Chikyow, T. Ishii, S. Miyazawa, _et al._ , “Heteroepitaxial growth of $\beta$-ligao 2 thin films on zno,” Journal of applied physics 92, 5587–5589 (2002).
  * Chen  _et al._ [2014b] C. Chen, C.-A. Li, S.-H. Yu, and M. M. Chou, “Growth and characterization of $\beta$-ligao2 single crystal,” Journal of crystal growth 402, 325–329 (2014b).
  * [26] https://www.vasp.at/.
  * Kresse and Furthmiiller [1996] G. Kresse and J. Furthmiiller, “Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set,” Computational Materials Science 6, 15–50 (1996).
  * Blöchl [1994] P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B 50, 17953–17979 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, “From ultrasoft pseudopotentials to the projector augmented-wave method,” Phys. Rev. B 59, 1758–1775 (1999).
  * Freysoldt  _et al._ [2014] C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, and C. G. Van de Walle, “First-principles calculations for point defects in solids,” Rev. Mod. Phys. 86, 253–305 (2014).
  * Freysoldt, Neugebauer, and Van de Walle [2009] C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, “Fully Ab Initio finite-size corrections for charged-defect supercell calculations,” Phys. Rev. Lett. 102, 016402 (2009).
  * Ivády  _et al._ [2013] V. Ivády, I. Abrikosov, E. Janzén, and A. Gali, “Role of screening in the density functional applied to transition-metal defects in semiconductors,” Physical Review B 87, 205201 (2013).


# Supplementary information of “Absence of quantization in the circular
photogalvanic effect in disordered chiral Weyl semimetals”

Ang-Kun Wu<EMAIL_ADDRESS>Department of Physics and Astronomy, Center for
Materials Theory, Rutgers University, Piscataway, New Jersey 08854, USA
Theoretical Division, T-4 and CNLS, Los Alamos National Laboratory, Los
Alamos, New Mexico 87545,USA  Daniele Guerci  Center for Computational Quantum
Physics, Flatiron Institute, New York, NY 10010, USA  Yixing Fu  Department of
Physics and Astronomy, Center for Materials Theory, Rutgers University,
Piscataway, New Jersey 08854, USA  Justin H. Wilson  Department of Physics and
Astronomy, Louisiana State University, Baton Rouge, Louisiana 70803, USA
Center for Computation and Technology, Louisiana State University, Baton
Rouge, Louisiana 70803, USA  J. H. Pixley<EMAIL_ADDRESS>Department of Physics and Astronomy, Center for Materials Theory, Rutgers
University, Piscataway, New Jersey 08854, USA  Center for Computational
Quantum Physics, Flatiron Institute, New York, NY 10010, USA

##  I KPM formalism for the CPGE

To reach large system sizes, we employ the kernel polynomial method (KPM) and
expand the delta function and Green’s functions in the response
$\chi^{\alpha\beta\gamma}$ in terms of Chebyshev polynomials. Based on the
work from [1, 2], we find the expression is part of the second-order
conductivity

$\begin{split}&\chi^{\alpha\beta\gamma}(\omega_{1},\omega_{2})=\frac{e^{3}}{\hbar^{3}}\sum_{nmp}\Lambda_{nmp}(\omega_{1},\omega_{2})\Gamma^{\alpha,\beta,\gamma}_{nmp},\\\
&\Gamma^{\alpha,\beta,\gamma}_{nmp}=\frac{\mathrm{Tr}}{N}\bigg{[}\hat{j}^{\alpha}\frac{T_{n}(\tilde{H})}{1+\delta_{n0}}\hat{j}^{\beta}\frac{T_{m}(\tilde{H})}{1+\delta_{m0}}\hat{j}^{\gamma}\frac{T_{p}(\tilde{H})}{1+\delta_{p0}}\bigg{]}.\end{split}$
(1)

where $\tilde{H}=(H-b)/D$ the rescaled Hamiltonian for the KPM calculation,
and $\delta_{ij}$ here is the Kronecker delta (different from the Dirac deltas
of the system defined below). The function $\Lambda_{nmp}$ is a numerically
tensor

$\begin{split}\Lambda_{nmp}(\omega_{1},\omega_{2})&=\hbar^{2}\int_{-\infty}^{\infty}\mathrm{d}\epsilon
f(\epsilon)\bigg{[}g_{n}^{R}(\epsilon/\hbar+\omega_{1}+\omega_{2})g_{m}^{R}(\epsilon/\hbar+\omega_{2})\Delta_{p}(\epsilon)\\\
&+g_{n}^{R}(\epsilon/\hbar+\omega_{1})\Delta_{m}(\epsilon)g_{p}^{A}(\epsilon/\hbar-\omega_{2})+\Delta_{n}(\epsilon)g_{m}^{A}(\epsilon/\hbar-\omega_{1})g_{p}^{A}(\epsilon/\hbar-\omega_{1}-\omega_{2})\bigg{]}\end{split}$
(2)

where the Jackson kernel $K_{n}^{J}$ is absorbed into the delta and Green’s
function

$\begin{split}&\Delta_{n}(\epsilon)=\frac{2T_{n}(\epsilon)}{\pi\sqrt{1-\epsilon^{2}}}K_{n}^{J},\quad
g_{n}^{\pm,0^{+}}(\epsilon)=\mp 2i\frac{e^{\mp in\arccos(\epsilon\pm
i0^{+})}}{\sqrt{1-(\epsilon\pm i0^{+})^{2}}}K_{n}^{J}\\\
&K_{n}^{J}=\frac{1}{N_{C}+1}\bigg{[}(N-n+1)\cos\frac{\pi
n}{N_{C}+1}+\sin\frac{\pi
n}{N_{C}+1}\cot\frac{\pi}{N_{C}+1}\bigg{]},\end{split}$ (3)

and $T_{n}(x)=\cos(n\arccos(x))$ is the first kind of the Chebyshev
polynomials. Here, $g_{n}^{+,0^{+}}$ represents the retarted Green’s function
while $g_{n}^{-,0^{+}}$ is the advanced one. In deriving the Eq. 1, the Dirac
deltas and Green’s functions are implicitly expanded in terms of Chebyshev
polynomials

$\begin{split}\delta(\epsilon-\tilde{H})&=\sum_{n=0}^{N_{C}-1}K_{n}^{J}\Delta_{n}(\epsilon)\frac{T_{n}(\tilde{H})}{1+\delta_{n0}},\\\
g^{\pm,0^{+}}(\epsilon,\tilde{H})&=\frac{\hbar}{\epsilon-\tilde{H}\pm
i0^{+}}=\hbar\sum_{n=0}^{N_{C}-1}K_{n}^{J}g_{n}^{\pm,0^{+}}(\epsilon)\frac{T_{n}(\tilde{H})}{1+\delta_{n0}}.\end{split}$
(4)

From the main text, Eq. (5), we know that the CPGE response is diverging in
the DC $\Omega\to 0$ and thermodynamic limit $L\to\infty$. Taking $\Omega=0$
we have the KPM estimate $\beta(E_{F})=y(E_{F})/\Omega_{B}$ where the plateau
value of $\beta$, denoted $\beta_{P}$ will be strongly dependent on $N_{C}$.
For sufficiently large system sizes $\Omega_{B}=\delta_{\mathrm{KPM}}/N_{C}$
is $L$ independent and we are able to find
$\delta_{\mathrm{KPM}}=31.74\beta_{0}/D^{3}$ where $D$ is the full bandwidth,
as shown in Fig. 1. As a result, the height of the plateau in the clean limit
with KPM goes as $\beta_{P}=N_{C}\beta_{0}/\delta_{\mathrm{KPM}}$. Note that
for convenience, the default energy scale is all in terms of $t$, the hopping
strength of the two band model, i.e., $E_{F}\equiv E_{F}/t,W\equiv W/t$ in
this supplementary.

Figure 1: The bare response of the CPGE with various KPM expansions $N_{C}$
and system sizes $L$. Clean Weyl systems are simulated with the periodic
boundary condition. The bare value of $\Omega$ is taken to be exactly 0. We
take $\omega_{1}=\gamma=0.8t$ to have relatively large plateau. Upper panel:
Bare response as a function of $E_{F}$ with $L=100$. The quantization is
manifested by the clear plateau while the heights increase with $N_{C}$. Lower
panel: The height of plateau as a function of system sizes. For fixed $N_{C}$,
quantization is converged with large $L$.

Note that in the main text, results are shown without error bars due to the
averaging over the KPM tensor $\Gamma_{nmp}^{\alpha,\beta,\gamma}$

$\langle\beta\rangle\sim\langle\chi^{\alpha\beta\gamma}\rangle\sim\langle\Gamma_{nmp}^{\alpha,\beta,\gamma}\rangle,$
(5)

where $\langle\cdots\rangle$ represents disorder average, $\sim$ means
linearly dependent and $\Lambda_{nmp}(\omega_{1},\omega_{2})$ is system-
independent integrals as we fixed $\omega_{1},\omega_{2}$. By averaging
$\Gamma_{nmp}^{\alpha,\beta,\gamma}$, we do not need to compute the triple sum
$\mathcal{O}(N_{C}^{3})$ for $\chi^{\alpha\beta\gamma}$, and thus save a lot
of computational effort. The disorder-sample-dependent effect is small and
discussed with the rare region effects in the last section. Also, in order to
conduct the KPM tensor averaging, we fix the half-bandwidth $D$ for all
disorder samples (instead of varying $D$ for every sample, and $b=0$ in
$\tilde{H}=(H-b)/D$) to keep the energy in the same scale, where we sacrifice
a little KPM resolution $\mathcal{O}(D/N_{C})$ by slightly larger $D$ for
computational advantage.

##  II Perturbation Theory for the CPGE from effective low-energy model

Consider the effective low-energy continuum theory,

$H_{0}=\sum_{\mathbf{k}}\bigg{[}\psi_{1\mathbf{k}}^{\dagger}(v_{F}\mathbf{k}\cdot\boldsymbol{\sigma}-\mu_{1})\psi_{1\mathbf{k}}+\psi_{2\mathbf{k}}^{\dagger}(-v_{F}\mathbf{k}\cdot\boldsymbol{\sigma}-\mu_{2})\psi_{2\mathbf{k}}\bigg{]},$
(6)

with $v_{F}=t,\mu_{1}<0,\mu_{2}>0$ (the [3] construction). Indices $1,2$
denote the two nodes. The Matsubara Green’s function for node 1 is

$G_{1}(i\epsilon_{n},\mathbf{k})=\frac{1}{i\epsilon_{n}-v_{F}\mathbf{k}\cdot\boldsymbol{\sigma}+\mu_{1}}=\frac{P_{+}(\mathbf{k})}{i\epsilon_{n}-v_{F}k+\mu_{1}}+\frac{P_{-}(\mathbf{k})}{i\epsilon_{n}+v_{F}k+\mu_{1}},$
(7)

where
$P_{\pm}(\mathbf{k})\equiv(I\pm\hat{\mathbf{k}}\cdot\boldsymbol{\sigma})/2$
with $\hat{\mathbf{k}}=\mathbf{k}/k$ projects on the eigenstate of the Dirac
operator with helicity $\pm$. Note also that
$P^{2}_{\pm}(\mathbf{k})=P_{\pm}(\mathbf{k})$. Similarly, for node 2, we find
the following

$G_{2}(i\epsilon_{n},\mathbf{k})=\frac{P_{+}(\mathbf{k})}{i\epsilon_{n}+v_{F}k+\mu_{2}}+\frac{P_{-}(\mathbf{k})}{i\epsilon_{n}-v_{F}k+\mu_{2}}.$
(8)

Note that due to the minus sign in the Hamiltonian
$(-v_{F}\mathbf{k}\cdot\boldsymbol{\sigma}-\mu_{2})$, we have a reverse sign
for the upper and lower bands for node 2. This effective low-energy theory for
the Weyl semi-metal is obtained by setting a momentum cutoff
$|\mathbf{k}-\mathbf{k}_{a}|<\Lambda$ around each Weyl node $a$. Adding the
onsite random scalar potential

$H=H_{0}+H_{\mathrm{dis}}=\sum_{\mathbf{k}}H_{0}(\mathbf{k})+\sum_{\mathbf{r}}\Psi^{\dagger}(\mathbf{r})V(\mathbf{r})\Psi(\mathbf{r})$
(9)

where $V(\mathbf{r})$ is drawn from a Gaussian distribution
$\mathcal{N}(0,W^{2})$ with zero average and variance $W^{2}$:

$\overline{V(\mathbf{r})}=0,\quad\overline{V(\mathbf{r})V(\mathbf{r}^{\prime})}=W^{2}\delta(\mathbf{r}-\mathbf{r}^{\prime}).$
(10)

Figure 2: The Feynman diagrams for the bare second order conductivity. Lines
with arrows represent Fermion fields and Green’s functions with certain
frequency. Curly curves represent external fields.

Our aim is to compute the second order in $W$ correction to the CPGE which is
obtained from the second-order conductivity

$\begin{split}j^{\gamma}(\Omega)&=\frac{\chi^{\alpha\beta\gamma}(\omega_{1},\omega_{2})+\chi^{\beta\alpha\gamma}(\omega_{2},\omega_{1})}{\omega_{1}\omega_{2}}E^{\alpha}(\omega_{1})E^{\beta}(\omega_{2}),\\\
\chi^{\alpha\beta\gamma}(i\omega_{1},i\omega_{2})&=T\sum_{\varepsilon_{n}}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\mathrm{Tr}[\hat{j}^{\alpha}G(i\varepsilon_{n}-i\omega_{1},\mathbf{k})\hat{j}^{\beta}G(i\varepsilon_{n}-i\Omega,\mathbf{k})\hat{j}^{\gamma}G(i\varepsilon_{n},\mathbf{k})],\end{split}$
(11)

where $\Omega=\omega_{1}+\omega_{2}\to 0$ and the diagram for $\chi$ is shown
in Fig. 2(b). Note that in the basis $(\psi_{1},\psi_{2})^{T}$ the bare
current operator for Weyl fermions is given by

$\hat{j}^{\alpha}=e\frac{\delta\hat{H}_{0}(\mathbf{k})}{\delta
k^{\alpha}}=ev_{F}\begin{pmatrix}\sigma^{\alpha}&0\\\
0&-\sigma^{\alpha}\end{pmatrix}.$ (12)

Before moving on, we observe that the numerical results from the KPM indicate
that the bare plateau value $y_{q}$ of the function

$y\equiv\mathrm{Im}\bigg{[}\frac{\chi^{\alpha\beta\gamma}(\omega_{1},\omega_{2})+\chi^{\beta\alpha\gamma}(\omega_{2},\omega_{1})}{\omega_{1}\omega_{2}}\bigg{]},$
(13)

scales linearly with the expansion order $N_{C}$ as shown in Fig. 1. For a
fixed $N_{C}$, the results $y_{q}$ are stabilized for large system sizes $L$.
As we already take the bare value $\Omega=0,$ or $\omega_{2}=-\omega_{1}$, the
analytical divergence of the quantization $\beta_{0}=y_{q}/\Omega$ is now
governed by the finite broadening $O(1/N_{C})$ in the Green’s function (finite
life time $O(N_{C})$) introduced via the KPM. It is known that for the
disordered electron gas (metal), the disorder could also induce a finite
lifetime via self-energy correction, which can also change the value of
$y_{q}$. Numerically in Fig. 3, we observed that there is a clear linear
scaling of $1/y_{q}$ versus the variance of disorder $W^{2}$. In order to
understand this results, we study the first-order perturbation theory with
disorder upon the $\chi^{\alpha\beta\gamma}$ below.

Figure 3: The scaling of the plateau value $1/y_{q}$ versus disorder strength
$W^{2}$, where the $W^{2}$ is the variance of the Gaussian disorder. We take
$\omega_{1}=\gamma=0.8t$ to have relatively large plateau. $y_{q}$ value is
taken at approximately the maximal value of the plateau at around
$E_{F}=-0.1$. The disordered results are averaged over 50 realizations.

Figure 4: The Feynman diagrams for the first-order correction to the
$\chi^{\alpha\beta\gamma}$. (a-b) second-order corrections (bubble diagrams)
to the three single-particle Green’s functions. (c-d) Vertex corrections to
the current.

The lowest order correction to $\chi^{\alpha\beta\gamma}$ gives 12
contributions (diagrams), where 6 of them are corrections to the single
Green’s function with three frequencies as shown in Fig. 4 (a-d) and the rest
are vertex corrections to the current as shown in Fig. 4 (e-h). We now compute
the self-energy correction given by disorder scattering processes. Within the
lowest order Born approximation, we find:

$\Sigma(i\epsilon)=\frac{W^{2}}{N}\sum_{\mathbf{k}}G(\mathbf{k},i\epsilon)=\frac{W^{2}}{N}\sum^{\Lambda}_{\mathbf{k}}\sum_{l}G_{l}(\mathbf{k},i\epsilon),$
(14)

where we consider the integral over momentum is performed in a circle of
radius $\Lambda$ around the two low energy Weyl nodes. We solve the Dyson
equation (14) to lowest order in $W^{2}$, which consists of replacing $G$ in
Eq. (14) with the bare Green’s function given in Eqs. (7) and (8). To make
further progress, we compute the integral:

$T_{l}(i\epsilon)=\frac{1}{N}\sum_{\mathbf{k}}G_{l}(\mathbf{k},i\epsilon)=\frac{1}{2}\sum_{\eta=\pm}\int\frac{d^{3}\mathbf{k}}{\Omega_{\rm
BZ}}\left[\mathbf{1}+\eta\frac{\mathbf{k}}{|\mathbf{k}|}\cdot{\pmb{\sigma}}\right]\frac{1}{i\epsilon-\eta
v_{Fl}|\mathbf{k}|+\mu_{l}},$ (15)

where $v_{F1}=-v_{F2}=v_{F}$ and $T_{l}(i\epsilon)$ is a $2\times 2$ matrix in
the spin degrees of freedom. Performing the angular integral, we find:

$T_{l}(i\epsilon)=\sum_{\eta=\pm}\frac{2\pi}{\Omega_{\rm
BZ}}\int^{\Lambda}_{0}dk\frac{k^{2}}{i\epsilon-\eta v_{Fl}k+\mu_{l}}.$ (16)

Notice that the UV cutoff must be introduced since the integrand is not
bounded and diverges for large $k$. The integral over $k$ can be performed
analytically. It is convenient to perform the integral employing the
approximation:

$\frac{1}{i\epsilon-\eta v_{Fl}k+\mu_{l}}=\frac{1}{\mu_{l}-\eta
v_{Fl}k}-i\pi\text{sign}\epsilon\delta(\eta v_{Fl}k-\mu_{l}).$ (17)

Employing the latter expression, we get:

$T_{l}(i\epsilon)=\frac{4\mu_{l}\pi}{\Omega_{\rm
BZ}}P\int^{\Lambda}_{0}\frac{k^{2}}{\mu^{2}_{l}-v^{2}_{Fl}k^{2}}-i\pi\frac{2\pi}{\Omega_{\rm
BZ}}\text{sign}(\epsilon)\frac{\mu^{2}_{l}}{v^{3}_{Fl}}.$ (18)

In the limit of small chemical potential, we are left with the simple
contribution:

$T(i\epsilon)=\rho(\mu_{l})\left[-\frac{\Lambda
v_{Fl}}{\mu_{l}}+\text{Atanh}\frac{v_{Fl}\Lambda}{\mu_{l}}\right]-\frac{i\pi}{2}\text{sign}(\epsilon)\rho(\mu_{l}),$
(19)

which gives the self-energy:

$\Sigma(i\epsilon)=W^{2}\sum_{a=1}^{2}T_{a}(i\epsilon).$ (20)

In the previous expression we have introduced the self-energy

$\rho(\omega)=\frac{3\epsilon^{2}}{2(v_{F}\Lambda)^{2}},$ (21)

with $\epsilon\in[-v_{F}\Lambda,v_{F}\Lambda]$.

The vertex corrections can be similarly carried out. For a single current
$\hat{j}^{\alpha}$, the vertex correction for each node is

$\begin{split}[\delta
j^{\alpha}]_{1}&=\frac{W^{2}}{2N}\sum_{\mathbf{k}}G^{(1)}(i\epsilon_{n},\mathbf{k})ev_{F}\sigma^{\alpha}G^{(1)}(i\epsilon_{n}-i\omega_{1},\mathbf{k})\\\
[\delta
j^{\alpha}]_{2}&=\frac{W^{2}}{2N}\sum_{\mathbf{k}}G^{(2)}(i\epsilon_{n},\mathbf{k})(-1)ev_{F}\sigma^{\alpha}G^{(2)}(i\epsilon_{n}-i\omega_{1},\mathbf{k}).\end{split}$
(22)

We notice that the cross term vanishes to $W^{2}$ order in the perturbation
theory. With the prototype integral
$\int\frac{x^{2}dx}{(a+bx)(c+dx)}=\frac{x}{bd}+\frac{a^{2}\ln(a+bx)}{b^{2}(bc-
ad)}-\frac{c^{2}\ln(c+dx)}{d^{2}(bc-ad)}$, consider the 4 integrals

$\begin{split}&\int_{0}^{\Lambda}\frac{4\pi
k^{2}dk}{(2\pi)^{3}}\frac{1}{i\epsilon_{n}-v_{F}k+\mu_{1}}\frac{1}{i\epsilon_{n}-v_{F}k+\mu_{1}-\omega_{1}}\\\
&=\frac{1}{2\pi^{2}}\bigg{[}\frac{\Lambda}{v_{F}^{2}}+\frac{(i\epsilon_{n}+\mu_{1})^{2}}{v_{F}^{3}\omega_{1}}\ln\frac{i\epsilon_{n}+\mu_{1}-v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}}-\frac{(i\epsilon_{n}+\mu_{1}-\omega_{1})^{2}}{v_{F}^{3}\omega_{1}}\ln\frac{i\epsilon_{n}+\mu_{1}-\omega_{1}-v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-\omega_{1}}\bigg{]}\\\
&\int_{0}^{\Lambda}\frac{4\pi
k^{2}dk}{(2\pi)^{3}}\frac{1}{i\epsilon_{n}-v_{F}k+\mu_{1}}\frac{1}{i\epsilon_{n}+v_{F}k+\mu_{1}-\omega_{1}}\\\
&=\frac{1}{2\pi^{2}}\bigg{[}-\frac{\Lambda}{v_{F}^{2}}-\frac{(i\epsilon_{n}+\mu_{1})^{2}}{v_{F}^{3}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}-v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}}+\frac{(i\epsilon_{n}+\mu_{1}-\omega_{1})^{2}}{v_{F}^{3}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}-\omega_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-\omega_{1}}\bigg{]}\\\
&\int_{0}^{\Lambda}\frac{4\pi
k^{2}dk}{(2\pi)^{3}}\frac{1}{i\epsilon_{n}+v_{F}k+\mu_{1}}\frac{1}{i\epsilon_{n}-v_{F}k+\mu_{1}-\omega_{1}}\\\
&=\frac{1}{2\pi^{2}}\bigg{[}-\frac{\Lambda}{v_{F}^{2}}+\frac{(i\epsilon_{n}+\mu_{1})^{2}}{v_{F}^{3}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}}-\frac{(i\epsilon_{n}+\mu_{1}-\omega_{1})^{2}}{v_{F}^{3}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}-\omega_{1}-v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-\omega_{1}}\bigg{]}\\\
&\int_{0}^{\Lambda}\frac{4\pi
k^{2}dk}{(2\pi)^{3}}\frac{1}{i\epsilon_{n}+v_{F}k+\mu_{1}}\frac{1}{i\epsilon_{n}+v_{F}k+\mu_{1}-\omega_{1}}\\\
&=\frac{1}{2\pi^{2}}\bigg{[}\frac{\Lambda}{v_{F}^{2}}-\frac{(i\epsilon_{n}+\mu_{1})^{2}}{v_{F}^{3}\omega_{1}}\ln\frac{i\epsilon_{n}+\mu_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}}+\frac{(i\epsilon_{n}+\mu_{1}-\omega_{1})^{2}}{v_{F}^{3}\omega_{1}}\ln\frac{i\epsilon_{n}+\mu_{1}-\omega_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-\omega_{1}}\bigg{]}\end{split}$
(23)

Summing over these integrals gives

$\begin{split}&[\delta
j^{\alpha}]_{1}=\frac{W^{2}ev_{F}\sigma^{\alpha}}{2\pi^{2}}\bigg{[}\frac{(i\epsilon_{n}+\mu_{1})^{2}}{v_{F}^{3}\omega_{1}}\ln\frac{i\epsilon_{n}+\mu_{1}-v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}+v_{F}\Lambda}+\frac{(i\epsilon_{n}+\mu_{1}-\omega_{1})^{2}}{v_{F}^{3}\omega_{1}}\ln\frac{i\epsilon_{n}+\mu_{1}-\omega_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-\omega_{1}-v_{F}\Lambda}\\\
&+\frac{(i\epsilon_{n}+\mu_{1})^{2}}{v_{F}^{3}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-v_{F}\Lambda}+\frac{(i\epsilon_{n}+\mu_{1}-\omega_{1})^{2}}{v_{F}^{3}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}-\omega_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-\omega_{1}-v_{F}\Lambda}\bigg{]}\\\
&=\frac{ev_{F}\sigma^{\alpha}W^{2}}{2\pi^{2}}\bigg{[}-\frac{2(i\epsilon_{n}+\mu_{1})^{2}(i\epsilon_{n}+\mu_{1}-\omega_{1})}{v_{F}^{3}\omega_{1}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-v_{F}\Lambda}\\\
&+\frac{2(i\epsilon_{n}+\mu_{1}-\omega_{1})^{2}(i\epsilon_{n}+\mu_{1})}{v_{F}^{3}\omega_{1}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}-\omega_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-\omega_{1}-v_{F}\Lambda}\bigg{]}\end{split}$
(24)

Sending $\mu_{1}\to\mu_{1}-\omega_{1},\omega_{1}\to-\omega_{1}$, we have

$\begin{split}&[\delta
j^{\beta}]_{1}=\frac{ev_{F}\sigma^{\beta}W^{2}}{2\pi^{2}}\bigg{[}-\frac{2(i\epsilon_{n}+\mu_{1}-\omega_{1})^{2}(i\epsilon_{n}+\mu_{1})}{-v_{F}^{3}\omega_{1}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}-\omega_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-\omega_{1}-v_{F}\Lambda}\\\
&+\frac{2(i\epsilon_{n}+\mu_{1})^{2}(i\epsilon_{n}+\mu_{1}-\omega_{1})}{-v_{F}^{3}\omega_{1}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-v_{F}\Lambda}\bigg{]}\end{split}$
(25)

$\omega_{1}\to 0$ gives $[j^{\gamma}]^{(f)}=0$. Thus, combine the contribution
from each current, the correction to the CPGE response due to node $1$

$\begin{split}&[\delta\beta(\omega_{1})]_{1}=\frac{\beta_{0}W^{2}}{2\pi^{2}}\bigg{[}\frac{4(i\epsilon_{n}+\mu_{1}-\omega_{1})^{2}(i\epsilon_{n}+\mu_{1})}{v_{F}^{3}\omega_{1}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}-\omega_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-\omega_{1}-v_{F}\Lambda}\\\
&-\frac{4(i\epsilon_{n}+\mu_{1})^{2}(i\epsilon_{n}+\mu_{1}-\omega_{1})}{v_{F}^{3}\omega_{1}(2i\epsilon_{n}+2\mu_{1}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{1}-v_{F}\Lambda}\bigg{]}\end{split}$
(26)

And similarly for node $2$, the opposite monopole charge has opposite
correction.

$\begin{split}&[\delta\beta(\omega_{1})]_{2}=\frac{-\beta_{0}W^{2}}{2\pi^{2}}\bigg{[}\frac{4(i\epsilon_{n}+\mu_{2}-\omega_{1})^{2}(i\epsilon_{n}+\mu_{2})}{v_{F}^{3}\omega_{1}(2i\epsilon_{n}+2\mu_{2}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{2}-\omega_{1}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{2}-\omega_{1}-v_{F}\Lambda}\\\
&-\frac{4(i\epsilon_{n}+\mu_{2})^{2}(i\epsilon_{n}+\mu_{2}-\omega_{1})}{v_{F}^{3}\omega_{1}(2i\epsilon_{n}+2\mu_{2}-\omega_{1})}\ln\frac{i\epsilon_{n}+\mu_{2}+v_{F}\Lambda}{i\epsilon_{n}+\mu_{2}-v_{F}\Lambda}\bigg{]}\end{split}$
(27)

In the limit of $i\epsilon_{n}\to 0,v_{F}\Lambda\gg\mu_{1},\omega_{1}$,

$\begin{split}&[\delta\beta(\omega_{1})]_{1}\approx\beta_{0}\frac{2W^{2}}{\pi^{2}}\bigg{[}\frac{(\mu_{1}-\omega_{1})^{2}\mu_{1}}{v_{F}^{3}\omega_{1}(2\mu_{1}-\omega_{1})}(2\frac{\mu_{1}-\omega_{1}}{v_{F}\Lambda})-\frac{\mu_{1}^{2}(\mu_{1}-\omega_{1})}{v_{F}^{3}\omega_{1}(2\mu_{1}-\omega_{1})}(2\frac{\mu_{1}}{v_{F}\Lambda})\bigg{]}\\\
&=-\beta_{0}\frac{4W^{2}}{\pi^{2}}\frac{\mu_{1}(\mu_{1}-\omega_{1})}{v_{F}^{4}\Lambda}\\\
&\delta\beta(\omega_{1})=-\beta_{0}\frac{4W^{2}}{\pi^{2}v_{F}^{4}\Lambda}[\mu_{1}(\mu_{1}-\omega_{1})-\mu_{2}(\mu_{2}-\omega_{1})]=-\beta_{0}\frac{4W^{2}}{\pi^{2}v_{F}^{4}\Lambda}[\mu_{1}^{2}-\mu_{2}^{2}+\omega_{1}(\mu_{2}-\mu_{1})]\end{split}$
(28)

##  III Numerical Disorder dependence of the CPGE

The disorder dependence of the CPGE behaves differently at various Fermi
energies $E_{F}$ as shown in Fig. 5. At the edges of the quantization plateau,
$E_{F}=-1.3,0.3$, the CPGE response linearly correlates with disorder

$\beta\approx\beta^{(0)}+\frac{\beta^{(2)}}{2!}W^{2}+\cdots,$ (29)

where $\beta^{(n)}$ is the $n$th derivative of $\beta$ with respective to $W$
and $\beta^{(2)}>0$. In contrast, at the quantized plateau,
$E_{F}=-1.0,-0.8,-0.6$, the CPGE decreases linearly, $\beta^{(2)}<0$. From the
perspective of quasi-particle lifetime, the plateau value of the CPGE scales
inversely with disorder $W^{2}$

$\frac{\beta}{\beta_{0}}=\frac{1}{1+CW^{2}}\to\frac{1}{\beta}=\frac{1}{\beta_{0}}+\frac{C}{\beta_{0}}W^{2}\approx\frac{1}{\beta^{(0)}}+\frac{1}{2!}\bigg{[}\frac{2\beta^{(1)}}{(\beta^{(0)})^{3}}-\frac{\beta^{(2)}}{(\beta^{(0)})^{2}}\bigg{]}W^{2}+\cdots,$
(30)

where the inverse dependence is confirmed in Fig. 5 $\mathbf{b}$ at small
$W^{2}$ limit.

aleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
bleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]

Figure 5:  (a) The CPGE quantization in the clean limit and its suppression
with increasing disorder strength with $L=100,N_{C}=2^{9}$. In KPM, the half-
bandwidth $D=6t$ to fit with various disorder strength. Each disorder result
$W$ is averaged over $50$ disorder realizations and twisted boundary
conditions. The normalization factor is the same as in the clean case
$C_{\mathrm{KPM}}=31.74\beta_{0}$. Dashed curves mark perturbative corrections
with quasiparticle lifetime, frequency shift, and vertex correction. Vertical
dotted lines mark locations $E_{F}=[-1.0,-0.8,-0.6]t$, while vertical dashed
lines mark $E_{F}=[-1.3,0.3]t$. (b) The normalized CPGE response $1/\beta$ as
a function of disorder strength $W^{2}$ at various Fermi levels with
$N_{C}=512$ from marked $E_{F}$. The dashed curve shows a linear fit to the
first $5$ small $W$ of the middle curve, where $1/\beta\approx
0.999+7.39W^{2}$. The dash-dotted curves are the second-order polynomial fit
with first $5$ points. Inset: The disorder broadening $\beta$ at two different
$E_{F}$ (plateau edges) from (a), which have vanished responses from the clean
systems. The dashed curves are linear fits from the first $5$ data points and
dash-dotted curves are second-order polynomial fits.

##  IV The CPGE response to rare states in Weyl semimetals

It has been shown that rare states in disordered Weyl semimetals give non-
perturbative contribution to the DOS and transport [4], when the two Weyl
nodes are at the same energy level $\mu_{1}=\mu_{2}=0$ (Fig. 6). However,
shifting the two Weyl nodes to different energies by introducing finite
$\gamma$ in the two-band model makes it hard to identify rare states from
squared energy (Fig. 6 $\mathbf{b}$) as the low energy states are burried in
the continuum (Fig. 6 $\mathbf{a}$) due to the shift. To identify rare states
under finite $\gamma$, we first find rare disorder realization from $\gamma=0$
from the squared energy, then slowly breaking the inversion via increasing
$\gamma$ at the rare disorder realization. Using the spatial profile of the
rare state (eigenstate with power-law decay from it maximal, Fig. 6
$\mathbf{c}$) and the inverse participation ratio (IPR, Fig. 6 $\mathbf{d}$)

$\mathrm{IPR}(n)=\sum_{i}|\phi_{n}(i)|^{4},$ (31)

we find that the rare state survive up to $\gamma=0.3$. Note that we need the
full profile of eigenstates at a finite $\gamma$ and the exact diagonalization
is required for the lattice systems, thus we can only use a smaller system
size, i.e., $L=18$.

aleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
bleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
cleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]
dleft3pt+/2right-3pt-/2up-3pt-/2down3pt+/2[left,up][vc,hc]

Figure 6: Rare states from increasing $\gamma$ for lattice systems $L=18$. (a)
The density of states with and without finite $\gamma$. Inset: the
corresponding dispersion for the two $\gamma$ cases. (b) The squared
eigenvalues across the disorder realizations. Vertical dashed lines mark the
rare state sample (red) and the perturbative sample (black). (c) The power-law
decay of the wavefunction as a function of the real space distance $r$ from
the maximum values of the eigenstates with various $\gamma$. The red line
marks the power law fitting. (d) The IPR of eigenstates for various $\gamma$.

Focusing on the finite $\gamma=0.3$, where the rare state survives, we look at
the CPGE response for rare and perturbative disorder realizations in Fig. 7
and compare them with the clean and averaged results over disorder. From Fig.
7, there is no qualitative difference between the rare, perturbative, and
averaged CPGE response, indicating that the rare states do not have a
prominent influence on the CPGE. Note that we are limited to a small system
size $L=18$ in order to keep track of rare states using exact diagonalization.
Thus, the quantization of the CPGE are less steep than the larger system
results in the main text. Also, due to the larger finite size effect, the
quantization value is larger in small systems, as shown in the inset of Fig.
7.

Figure 7: The CPGE responses to different disorder realizations at $W=0.5$
with $L=18,N_{C}=2^{9},\omega_{1}=\gamma=0.3t$. Dashed vertical lines mark the
locations $E_{F}=\pm\omega_{1}$. The clean, rare and perturbative disorder
samples are averaged over $100$ twisted boundary conditions while the averaged
result has $100$ disorder realizations. Inset: the CPGE plateau value at
various system sizes at $\gamma=0.8$ from Fig. 1.

## References

  * João and Lopes [2019] S. M. João and J. M. V. P. Lopes, J. Phys. Condens. Matter 32, 125901 (2019).
  * García  _et al._ [2015] J. H. García, L. Covaci, and T. G. Rappoport, Phys. Rev. Lett. 114, 116602 (2015).
  * Avdoshkin  _et al._ [2020] A. Avdoshkin, V. Kozii, and J. E. Moore, Phys. Rev. Lett. 124, 196603 (2020).
  * Pixley and Wilson [2021] J. Pixley and J. H. Wilson, Ann. Phys. 435, 168455 (2021).