# Massive graviton from diffeomorphism invariance

João M. L. de Freitasa, Iberê Kuntza

aDepartamento de Física, Universidade Federal do Paraná
PO Box 19044, Curitiba – PR, 81531-980, Brazil

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We describe a mechanism in which the graviton acquires a mass from the
functional measure without violating the diffeomorphism symmetry nor including
Stückelberg fields. Since gauge invariance is not violated, the number of
degrees of freedom goes as in general relativity. For the same reason,
Boulware-Deser ghosts and the vDVZ disconinuity do not show up. The graviton
thus becomes massive at the quantum level while avoiding the usual issues of
massive gravity.

##  1 Introduction

General relativity and the standard model have both been extremely successful
in describing fundamental phenomena. Nonetheless, many problems in the
interface of gravity and high-energy physics, such as the dark sector and
singularities, cannot be explained by either of these theories. This has led
to a plethora of attempts to modify these models.

One interesting modification of general relativity regards massive gravity.
Attempts of giving the graviton a mass dates back to the 30s [1], but it was
only recently that consistent theories of massive gravity have been found [2,
3, 4, 5, 6] (see also Ref. [7] for an in-depth review). A massive graviton
indeed introduces many issues, such as the violation of diffeomorphism (gauge)
invariance, the van Dam-Veltman-Zakharov (vDVZ) discontinuity and the presence
of Boulware-Deser ghosts [8]. Gauge invariance can be reinstated via
Stückelberg fields, which is legitimate but does require new degrees of
freedom. The vDVZ discontinuity [9, 10], namely the disagreement with general
relativity in the massless limit, is usually conjectured to be solved at the
non-linear regime by the Vainshtein screening mechanism [11, 12]. Only a few
models have successfully implemented the Vainshtein mechanism [12, 13, 14, 15,
16], which can also bring along the possibility of superluminal velocities [3,
15, 17, 18, 19, 20].

In this paper, we propose a novel procedure to give the graviton a mass
without introducing any of the aforementioned issues and without modifying
classical general relativity. Our finding is based on a non-trivial path-
integral measure, which is required for obtaining gauge invariant correlation
functions. The functional measure introduces non-linear loop corrections,
which act as a gravitational potential and results in a (quantum) mass for the
graviton in the linear regime. The most important point is the preservation of
the diffeomorphism invariance, which is responsible for keeping the theory
free of ghosts and of the zDVZ discontinuity. The so obtained mass is,
however, pure imaginary, thus precluding the existence of gravitons as
asymptotic states.

This paper is organized as follows. In Sec. 2, we review some aspects of the
functional measure in quantum field theory. The one-loop correction induced by
the functional measure is then studied in Sec. 3, where the graviton mass is
calculated. Comparison with experimental data yields stringent bounds on the
model. In Sec. 4, we obtain Newton’s potential and discuss its consequences.
Finally, we draw our conclusions in Sec. 5.

##  2 Functional measure in quantum field theory

A central object in quantum field theory is the generating functional:

$Z[J]=\int\mathrm{d}\mu[\varphi]e^{i\left(S[\varphi^{i}]+J_{i}\varphi^{i}\right)},$
(2.1)

where $S[\varphi]$ is the classical action for some generic set of fields
$\varphi^{i}=(\phi(x),A_{\mu}(x),g_{\mu\nu}(x),\ldots)$. All information
regarding any physical system is contained in its corresponding generating
functional $Z[J]$. In particular, correlation functions are obtained by simple
functional differentiations, from which scattering amplitudes can be obtained
using the LSZ (Lehmann–Symanzik–Zimmermann) formula. Despite its importance in
field theory, a rigorous mathematical foundation remains unknown, particularly
with respect to the functional measure $\mathrm{d}\mu[\varphi]$.
Operationally, one can however define the aforementioned object as [21, 22,
23, 24]

$\mathrm{d}\mu[\varphi]=\mathcal{D}\varphi^{i}\sqrt{\text{Det}\,G_{ij}}\ ,$
(2.2)

where $\mathcal{D}\varphi^{i}=\prod_{i}\mathrm{d}\varphi^{i}$ and
$\text{Det}\,G_{ij}$ denotes the functional determinant of the configuration-
space metric $G_{ij}$. The factor $\sqrt{\text{Det}\,G_{ij}}$ is required to
account for a non-trivial configuration space, whose typical example is that
of a non-linear sigma model.

The presence of a non-trivial functional measure is usually sidestepped by
writing $\delta^{(4)}(0)=0$ in dimensional regularization, in which case one
finds

$\displaystyle\text{Det}\,G_{ij}$
$\displaystyle=e^{\delta^{(4)}(0)\int\mathrm{d}^{4}x\sqrt{-g}\,\tr\log
G_{ij}}$ $\displaystyle\stackrel{{\scriptstyle dim.\,reg.}}{{=}}\ 1\ .$ (2.3)

One should, however, be extremelly careful when setting infinities to zero. In
fact, a Gaussian regularization

$\delta^{(4)}(x)=\frac{\Lambda^{4}}{(2\pi)^{2}}e^{\frac{-x^{2}\Lambda^{2}}{2}}\
,$ (2.4)

for some (soft) cutoff $\Lambda$, tells a completely different story. For
$\delta^{(4)}(0)=\tfrac{\Lambda^{4}}{(2\pi)^{2}}$, we find [25]:

$Z[J]=\int\mathcal{D}\varphi^{i}e^{i\left(S_{\text{eff}}[\varphi^{i}]+J_{i}\varphi^{i}\right)},$
(2.5)

with the Wilsonian effective action

$S_{\text{eff}}=\int\mathrm{d}^{4}x\sqrt{-g}\left(\mathcal{L}-i\zeta\,\tr\log
G_{IJ}\right)\ ,$ (2.6)

for some bare Lagrangian $\mathcal{L}$. Here $\zeta=\zeta(\Lambda)$ is a
Wilsonian coefficient whose running is such that

$\Lambda\frac{dZ[J]}{d\Lambda}=0\ .$ (2.7)

Note that the Dirac delta divergence is polynomial. The coefficient $\zeta$ is
thus expected to be UV sensitive.

One should note that the configuration-space metric $G_{ij}$ must be seen as
part of the definition of the theory. Although it is typically identified with
the bilinear form in the kinetic term [27, 28, 26, 29, 30], this
identification is not based on physical reasonings. Nonetheless, the
functional measure must be invariant under the underlying symmetries, which
allows one to determine $G_{ij}$ in the same spirit as effective field
theories. In particular, for pure gravity this procedure results in the well-
known DeWitt metric [31]. To lowest order, the most general configuration-
space metric for arbitrary fields yields [32, 33, 34]

$S_{\text{eff}}=\int\mathrm{d}^{4}x\sqrt{-g}\left(\mathcal{L}-i\gamma\,\tr\log|g_{\mu\nu}|\right)\
,$ (2.8)

where $\gamma$ is obtained from $\zeta$ via a finite renormalization.

In spite of the form of the correction, we stress that Eq. (2.8) does not
violate diffeomorphism invariance. The apparent violation results from the
fact that $\sqrt{\text{Det}\,G_{ij}}$ transforms as a (functional) scalar
density, thus so does the last term in Eq. (2.8). However, the
$\mathcal{D}\varphi^{i}$ also transforms as a scalar density in such a way
that the full measure $\mathcal{D}\varphi^{i}\sqrt{\text{Det}\,G_{ij}}$ is
invariant 111Note that $S_{\text{eff}}$ is the Wilsonian action before
performing the path integral. Because the functional measure and the classical
action are invariant under diffeomorphisms, the 1PI effective action
$\Gamma[g]$ naturally reflects such symmetry. At the one-loop level of
$\Gamma[g]$, the configuration-space metric enters the usual correction
$\log\text{Det}\,\mathcal{H}_{ij}$, transforming the Hessian
$\mathcal{H}_{ij}$, whose determinant is basis-dependent, into
$\mathcal{H}^{i}_{\ j}$, whose determinant is invariant [23, 35].. Therefore,
variations of the apparent symmetry-breaking term under spacetime
diffeomorphisms are canceled by the functional Jacobian that shows up from
$\mathcal{D}\varphi^{i}$, keeping the quantum theory and all observables
invariant 222Strictly speaking, the invariance of off-shell quantities, such
as the effective action, requires a connection in configuration space [26,
36]. Because we are mainly interested in the phenomenology of the functional
measure, we shall not dwell on this topic.. This is, in fact, the reason one
can generate a mass for the graviton without violating the gauge symmetry.

Finally, the above findings were obtained in the Lorentzian path integral. The
imaginary factor in Eq. (2.6) (hence in Eq. (2.8)) shows up when the
$i=\sqrt{-1}$ in the argument of the exponential in the Lorentzian path
integral is pulled out to write the measure as a correction to the classical
action. Defining the functional measure in the Euclidean formalism

$Z_{E}[J]=\int\mathrm{d}\mu[\varphi]e^{-\left(S_{\text{eff}}^{E}[\varphi^{i}]+J_{i}\varphi^{i}\right)}$
(2.9)

yields a real one-loop correction:

$S_{\text{eff}}^{E}=\int\mathrm{d}^{4}x\sqrt{-g}\left(\mathcal{L}-\gamma\,\tr\log|g_{\mu\nu}|\right)\
.$ (2.10)

One thus faces the problem of whether the path-integral measure should be
defined in the Euclidean space (and rotated back to real time) or straight in
the Lorentzian space. The former is required for a better mathematical
construction of the path integral, albeit still largely formal. On the other
hand, as far as our current experiments are concerned, Nature is fundamentally
Lorentzian. For this reason, we shall perform our calculations by defining the
measure on the latter. Naturally, predictions shall be different in different
schemes. At this level of formality, only time will tell which one, if any, is
correct.

##  3 Massive graviton

When the bare Lagrangian $\mathcal{L}$ is the Einstein-Hilbert term, the
functional measure changes the dynamics of space-time as follows 333We adopt
the metric signature $(-+++)$.:

$S_{\text{eff}}=\int\mathrm{d}^{4}x\sqrt{-g}\left(\frac{M_{p}^{2}}{2}R-i\gamma\,\tr\log|g_{\mu\nu}|+\mathcal{L}_{m}\right)\
,$ (3.11)

where $M_{p}$ is the reduced Planck mass, $R$ is the Ricci scalar and
$\mathcal{L}_{m}$ is the Lagrangian for matter fields. The corresponding
equations of motion read 444One should note that $T_{\mu\nu}$ is covariantly
conserved in the sense of the Slavnov-Taylor identities, i.e.
$\langle\nabla^{\mu}T_{\mu\nu}\rangle=0$. In particular, such conservation
holds at every loop order.:

$G_{\mu\nu}+i\frac{2\gamma}{M_{p}^{2}}\left[1+\frac{1}{2}\log(-g)\right]=\frac{1}{M_{p}^{2}}T_{\mu\nu}\
,$ (3.12)

where $T_{\mu\nu}$ is the energy-momentum tensor for $\mathcal{L}_{m}$. One
should note that general relativity is smoothly recovered in the limit
$\gamma\to 0$. The parameter $\gamma$ is proportional to the graviton mass
(see Eq. (3.15)), thus there is no tension with the experimental tests of
general relativity. As we shall see in the next section, gauge invariance
prevents the vDVZ discontinuity from appearing.

Perturbing the metric around Minkowski

$g_{\mu\nu}=\eta_{\mu\nu}+\frac{2}{M_{p}}h_{\mu\nu}$ (3.13)

leads to the Fierz-Pauli action

$\displaystyle S_{\text{eff}}=\int\mathrm{d}^{4}x\bigg{[}$
$\displaystyle-\frac{1}{2}\partial_{\lambda}h_{\mu\nu}\partial^{\lambda}h^{\mu\nu}+\frac{1}{2}\partial_{\lambda}h\partial^{\lambda}h-\partial_{\mu}h^{\mu\nu}\partial_{\nu}h+\partial_{\mu}h_{\nu\lambda}\partial^{\nu}h^{\mu\lambda}$
$\displaystyle-\frac{i}{2}m^{2}_{g}(h^{2}-h_{\mu\nu}h^{\mu\nu})+M_{p}^{-1}h_{\mu\nu}T^{\mu\nu}\bigg{]}\
.$ (3.14)

The graviton mass is thus given by

$m^{2}\equiv im^{2}_{g}=\frac{4i\gamma}{M_{p}^{2}}\ .$ (3.15)

Notice that the mass term in Eq. (3.14) comes with the correct relative
coefficient between $h^{2}$ and $h_{\mu\nu}h^{\mu\nu}$ for a ghost-free
theory. We stress that such a coefficient is not finely tuned by hand, it
follows directly from the quantum correction due to the functional measure.
Comparing the modulus of Eq. (3.15) to the bound found by LIGO [37]:

$m_{g}<1.2\times 10^{-22}\,\text{eV}$ (3.16)

translates into a bound on $\gamma$:

$\gamma<2\times 10^{-26}\,\text{GeV}^{4}\ .$ (3.17)

A few comments are in order. Because of (2.7), the graviton mass (3.15) runs
with the energy scale. In particular, in the classical limit $\hbar\to 0$, the
functional measure correction vanishes and so does the graviton mass $m$. We
stress that the massless limit $m\to 0$ (or, equivalently, $\gamma\to 0$) is
smooth as can be seen from the non-linear theory Eq. (3.12). Therefore, all
tests of general relativity are automatically satisfied. Secondly, since the
gauge symmetry is not broken, the counting of degrees of freedom goes as usual
for general relativity. Namely, a symmetric second rank tensor contains 10
degrees of freedom, 8 of which can be eliminated by gauge transformations,
yielding only 2 propagating modes rather than 5. This follows because one has
not started with massive gravity ab initio. Indeed, since the functional
measure does not contain derivatives, no new degrees of freedom show up and
the spectrum continues to be determined from the classical Einstein-Hilbert
action. The mass is generated only after quantization for the propagating
modes that had already been present in the classical theory. As a result,
Boulware-Deser ghosts do not appear in the theory as the action Eq. (3.11)
does not contain higher derivatives.

The resulting mass squared (3.15) is, however, pure imaginary. The imaginary
part of the mass corresponds to its width, which is a measure of its lifetime.
When the width is much smaller than the real part of the mass, the particle is
interpreted as a resonance, which can still be approximately treated as an
asymptotic state. A pure imaginary mass, on the other hand, is highly
unstable. They can only exist as virtual particles at intermediate stages of
physical processes, mediating the gravitational interaction, but rapidly
decaying into lighter by-products. Our result thus suggests that gravitons can
never appear as asymptotic states. Finally, we stress that our proposal has
also the advantage of being a top-down approach. We indeed know the non-linear
theory of massive gravity from the onset.

##  4 Newtonian potential

Although the above findings suggest that gravitons cannot be measured at
asymptotic states, virtual gravitons with complex mass do affect interactions.
An immediate consequence of such a massive graviton regards the modification
on the Newtonian potential. As a result of the presence of mass, one usually
expects a Yukawa potential. However, because the graviton mass is complex, the
resulting potential shall develop an oscillating behavior modulated by the
Yukawa decay, as we shall now see.

From Eq. (3.14), we obtain the effective equations of motion for the graviton:

$\Box h_{\mu\nu}-\Box
h\eta_{\mu\nu}+\partial_{\mu}\partial_{\nu}h+\partial_{\alpha}\partial_{\beta}h^{\alpha\beta}\eta_{\mu\nu}-\partial^{\lambda}\partial_{\mu}h_{\lambda\nu}-\partial^{\lambda}\partial_{\nu}h_{\lambda\mu}=-M_{p}^{-1}T_{\mu\nu}-m^{2}(h_{\mu\nu}-h\eta_{\mu\nu})$
(4.18)

The divergence and the trace of Eq. (4.18) read

$\displaystyle\partial^{\mu}h_{\mu\nu}$ $\displaystyle=\partial_{\nu}h$ (4.19)
$\displaystyle h$ $\displaystyle=\frac{M_{p}^{-1}}{3m^{2}}T\ .$ (4.20)

Using Eqs. (4.19)–(4.20) in Eq. (4.18) gives

$(\Box+m^{2})h_{\mu\nu}=-M_{p}^{-1}\left(T_{\mu\nu}-\frac{1}{2}T\eta_{\mu\nu}\right)+\left(\partial_{\mu}\partial_{\nu}h-\frac{1}{2}\eta_{\mu\nu}m^{2}h\right)\
.$ (4.21)

The first term on the RHS of Eq. (4.21) is the usual general relativistic
result (apart from the mass term on the LHS), which is then modified by the
second term on the RHS, leading to the vDVZ discontinuity. In our case,
because gauge invariance is not broken, the second term can be eliminated by a
choice of gauge. In fact, under a diffeomorphism such a term transforms as

$\partial_{\mu}\partial_{\nu}h-\frac{1}{2}\eta_{\mu\nu}m^{2}h\to\partial_{\mu}\partial_{\nu}h-\frac{1}{2}\eta_{\mu\nu}m^{2}h+2\partial_{\mu}\partial_{\nu}\partial_{\alpha}\xi^{\alpha}-\eta_{\mu\nu}m^{2}\partial_{\alpha}\xi^{\alpha}\
.$ (4.22)

Thus choosing $\partial_{\alpha}\xi^{\alpha}=-h/2$ results in:

$(\Box+m^{2})h_{\mu\nu}=-M_{p}^{-1}\left(T_{\mu\nu}-\frac{1}{2}T\eta_{\mu\nu}\right)\
.$ (4.23)

One should notice the apperance of the factor of $1/2$ instead of the infamous
$1/3$ of gauge-violating massive theories. Eq. (4.23) shows that the mass
produced by the functional measure is perfectly consistent with all general
relativistic tests as no vDVZ discontinuity takes place.

One can easily solve Eq. (4.23) in momentum space:

$h_{\mu\nu}=M_{p}^{-1}\int\frac{\differential[4]{p}}{(2\pi)^{4}}\frac{e^{ip_{\alpha}x^{\alpha}}}{p^{2}-m^{2}}\left[\tilde{T}_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}\tilde{T}\right]\
.$ (4.24)

For a static point source of mass $M$ at the origin:

$T_{\mu\nu}=M\delta_{\mu 0}\delta_{0\nu}\delta(\vec{x})\ ,$ (4.25)

one finds

$h_{00}=\frac{1}{2}\frac{M}{M_{p}}\frac{1}{4\pi r}e^{imr}\ .$ (4.26)

Despite Eq. (4.26) having a Yukawa-like functional form, the complex
exponential leads to novel predictions. Indeed, its real part provides the
Newtonian potential 555We recall that $h_{00}=-2V$.

$\displaystyle V(r)=-\frac{M}{16\pi
M_{p}}\frac{e^{-\tfrac{m_{g}}{\sqrt{2}}r}}{r}\cos\left(\frac{m_{g}}{\sqrt{2}}\,r\right)\
,$ (4.27)

whereas, by the optical theorem, its imaginary part relates to the total cross
section 666Note that the imaginary part of Eq. (4.26) remains finite at $r=0$
because $\sin(x)/x\to 1$.. At small distances $r\ll\sqrt{2}/m_{g}$, our result
recovers Newton’s potential:

$V(r)=-\frac{M}{16\pi
M_{p}\,r}\left(1-\frac{m_{g}}{\sqrt{2}}r+\frac{m_{g}^{3}}{6\sqrt{2}}r^{3}\right)+\cdots\
.$ (4.28)

The leading correction is constant, thus does not affect the dynamics, so the
first measurable new effect shows up only at next-to-leading order. We see
that Newton’s potential is modified at large distances $r\sim\sqrt{2}/m_{g}$.

The cosine function in Eq. (4.27) can turn the potential’s derivative
positive, thus creating islands of bounded motion of decreasing depth. At each
of these islands, gravity becomes repulsive. But because of the utterly small
graviton mass, such an effect is only felt at very large distances:

$r>\frac{\pi\sqrt{2}}{2}m_{g}^{-1}\ .$ (4.29)

At the first and highest barrier $r_{0}\sim 3.1\,m_{g}^{-1}$, the potential
height is given by (see Figure 1(a))

$V(r_{0})\sim 4.2\times 10^{-3}\,\frac{Mm_{g}}{M_{p}}\ .$ (4.30)

(a) First and highest barrier.

(b) First and deepest potential well.

Figure 1: Plot of the Newtonian potential corrected by the functional measure,
as given by Eq. (4.27).

If there is enough energy to overcome this barrier, the system alternates
between regions of attractive and repulsive gravity as the distance increases.
The system could get trapped between two consecutive barriers, thus forming a
gravitational bound state, should its energy be smaller than the potential
well (see Figure 1(b)). This effect, however, is rapidly weakened by the
exponential suppression in Eq. (4.27). Even at the deepest well depicted in
Figure 1(b), the existing energy from surrounding astrophysical events is
likely enough to keep such bound states from forming. On the other hand,
should the energy be smaller than Eq. (4.30), the system would not collapse as
the distance decreases, thus resolving the singularity at $r=0$. Note that the
ratio $m_{g}/M_{p}$ is utterly negligible in Eq. (4.30), thus only very
massive objects, such as black holes, could prevent such collapse from
happening.

##  5 Conclusions

Understanding the quantum nature of gravity requires, among other things,
grasping the graviton kinematics and dynamics. Quantum corrections trigger new
graviton self-interactions, which affect the graviton dynamics. As in any
gauge theory, it is generally believed that the diffeomorphism invariance
precludes the appearance of mass terms in the effective action of general
relativity. However, this follows by side-stepping the non-trivial Jacobian
that shows up in the path integral when diffeomorphisms (or, more generally,
field redefinitions) are performed. A proper regularization of the Jacobian
shows that it cannot be simply ignored.

A diffeomorphism-invariant functional measure thus requires the factor
$\sqrt{\text{Det}\,G_{ij}}$ to cancel out the functional Jacobian, in very
much the same way that $\sqrt{-g}$ is needed for integrations in curved
spacetimes. This additional factor contributes as a quantum effective
potential for the graviton, giving it a mass in the linear regime. Therefore,
general covariance is actually the reason for the existence, rather than for
the absence as commonly thought, of a mass for the graviton.

We stress that such a mass is not obtained by modifying general relativity.
One is simply quantizing general relativity by properly considering the
geometry of configuration space, which reflects in the definition of the
functional measure. One can view different configuration-space metrics as
different quantization approaches to the same classical theory. The naive
approach, where the functional Jacobian is disregarded, corresponds to a flat
configuration space. Unfortunately, there is no known physical principle other
than symmetry to help us determining the configuration-space metric. On
symmetry grounds, the lowest-order configuration-space metric is curved and
depends on the spacetime metric. This situation is dramatically different for
Yang-Mills fields, where the lowest-order configuration-space metric is
trivial and no gauge-invariant higher-order term exists. Therefore, gauge
invariance manifests itself quite differently in gravity than in the other
interactions, being the sole responsible for keeping Yang-Mills fields
massless (in the unbroken vacuum) and the graviton massive.

## Acknowledgments

IK is grateful to the National Council for Scientific and Technological
Development – CNPq (Grant No. 303283/2022-0) for partial financial support.

## References

  * [1] M. Fierz and W. Pauli, Proc. Roy. Soc. Lond. A 173 (1939), 211-232.
  * [2] G. R. Dvali and G. Gabadadze, Phys. Rev. D 63 (2001), 065007 [arXiv:hep-th/0008054 [hep-th]].
  * [3] G. R. Dvali, G. Gabadadze and M. Porrati, Phys. Lett. B 485 (2000), 208-214 [arXiv:hep-th/0005016 [hep-th]].
  * [4] G. R. Dvali, G. Gabadadze and M. Porrati, Phys. Lett. B 484 (2000), 112-118 [arXiv:hep-th/0002190 [hep-th]].
  * [5] E. A. Bergshoeff, O. Hohm and P. K. Townsend, Phys. Rev. Lett. 102 (2009), 201301 [arXiv:0901.1766 [hep-th]].
  * [6] C. de Rham, G. Gabadadze and A. J. Tolley, Phys. Rev. Lett. 106 (2011), 231101 [arXiv:1011.1232 [hep-th]].
  * [7] C. de Rham, Living Rev. Rel. 17 (2014), 7 [arXiv:1401.4173 [hep-th]].
  * [8] D. G. Boulware and S. Deser, Annals Phys. 89 (1975), 193.
  * [9] H. van Dam and M. J. G. Veltman, Nucl. Phys. B 22 (1970), 397-411.
  * [10] V. I. Zakharov, JETP Lett. 12 (1970), 312.
  * [11] A. I. Vainshtein, Phys. Lett. B 39 (1972), 393-394.
  * [12] E. Babichev and C. Deffayet, Class. Quant. Grav. 30 (2013), 184001 [arXiv:1304.7240 [gr-qc]].
  * [13] E. Babichev, C. Deffayet and R. Ziour, Phys. Rev. Lett. 103 (2009), 201102 [arXiv:0907.4103 [gr-qc]].
  * [14] C. Deffayet, G. R. Dvali, G. Gabadadze and A. I. Vainshtein, Phys. Rev. D 65 (2002), 044026 [arXiv:hep-th/0106001 [hep-th]].
  * [15] A. Nicolis, R. Rattazzi and E. Trincherini, Phys. Rev. D 79 (2009), 064036 [arXiv:0811.2197 [hep-th]].
  * [16] C. de Rham, A. J. Tolley and S. Y. Zhou, JHEP 04 (2016), 188 [arXiv:1602.03721 [hep-th]].
  * [17] M. A. Luty, M. Porrati and R. Rattazzi, JHEP 09 (2003), 029 [arXiv:hep-th/0303116 [hep-th]].
  * [18] A. Nicolis and R. Rattazzi, JHEP 06 (2004), 059 [arXiv:hep-th/0404159 [hep-th]].
  * [19] A. Adams, N. Arkani-Hamed, S. Dubovsky, A. Nicolis and R. Rattazzi, JHEP 10 (2006), 014 [arXiv:hep-th/0602178 [hep-th]].
  * [20] P. de Fromont, C. de Rham, L. Heisenberg and A. Matas, JHEP 07 (2013), 067 [arXiv:1303.0274 [hep-th]].
  * [21] E. Mottola, J. Math. Phys. 36 (1995) 2470 [arXiv:hep-th/9502109 [hep-th]].
  * [22] B. S. DeWitt, “ _The global approach to quantum field theory_ ”, Oxford University Press, Oxford, 2003.
  * [23] D. J. Toms, Phys. Rev. D 35 (1987) 3796.
  * [24] R. Casadio, A. Kamenshchik and I. Kuntz, Annals Phys. 449 (2023), 169203 [arXiv:2210.04368 [hep-th]].
  * [25] I. Kuntz and R. da Rocha, Nucl. Phys. B 993 (2023), 116258 [arXiv:2211.11913 [hep-th]].
  * [26] G. A. Vilkovisky, Nucl. Phys. B 234 (1984) 125.
  * [27] K. Meetz, J. Math. Phys. 10 (1969) 589.
  * [28] A. A. Slavnov and L. D. Faddeev, Teor. Mat. Fiz. 8 (1971) 297.
  * [29] E. S. Fradkin and G. A. Vilkovisky, Phys. Rev. D 8 (1973) 4241.
  * [30] E. S. Fradkin and G. A. Vilkovisky, Lett. Nuovo Cim. 19 (1977) 47.
  * [31] B. S. DeWitt, Phys. Rev. 160 (1967), 1113-1148.
  * [32] R. Casadio, A. Kamenshchik and I. Kuntz, Nucl. Phys. B 971 (2021), 115496 [arXiv:2102.10688 [hep-th]].
  * [33] I. Kuntz, R. Casadio and A. Kamenshchik, Mod. Phys. Lett. A 37 (2022) no.10, 2230007 [arXiv:2203.11259 [hep-th]].
  * [34] R. Casadio, A. Kamenshchik and I. Kuntz, Int. J. Mod. Phys. D 31 (2022) no.01, 2150130 [arXiv:2008.09387 [gr-qc]].
  * [35] P. Ellicott and D. J. Toms, Nucl. Phys. B 312 (1989), 700-714.
  * [36] B. S. DeWitt, “The effective action”, in Batalin, I.A. (Ed.) et al.: Quantum field theory and quantum statistics, Vol. 1, 191-222, CRC Press (1987).
  * [37] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116 (2016) no.22, 221101 [erratum: Phys. Rev. Lett. 121 (2018) no.12, 129902] [arXiv:1602.03841 [gr-qc]].


# Machine Learning guided high-throughput search of non-oxide garnets

Jonathan Schmidt  Institut für Physik, Martin-Luther-Universität Halle-
Wittenberg, D-06099 Halle, Germany  Hai-Chen Wang  Institut für Physik,
Martin-Luther-Universität Halle-Wittenberg, D-06099 Halle, Germany  Georg
Schmidt  Institut für Physik, Martin-Luther-Universität Halle-Wittenberg,
D-06099 Halle, Germany  Miguel A. L. Marques<EMAIL_ADDRESS>halle.de Institut für Physik, Martin-Luther-Universität Halle-Wittenberg,
D-06099 Halle, Germany

###### Abstract

Garnets, known since the early stages of human civilization, have found
important applications in modern technologies including magnetorestriction,
spintronics, lithium batteries, etc. The overwhelming majority of
experimentally known garnets are oxides, while explorations (experimental or
theoretical) for the rest of the chemical space have been limited in scope. A
key issue is that the garnet structure has a large primitive unit cell,
requiring an enormous amount of computational resources. To perform a
comprehensive search of the complete chemical space for new garnets, we
combine recent progress in graph neural networks with high-throughput
calculations. We apply the machine learning model to identify the potential
(meta-)stable garnet systems before systematic density-functional calculations
to validate the predictions. In this way, we discover more than 600 ternary
garnets with distances to the convex hull below 100 meV/atom with a variety of
physical and chemical properties. This includes sulfide, nitride and halide
garnets. For these, we analyze the electronic structure and discuss the
connection between the value of the electronic band gap and charge balance.

##  I Introduction

Garnets can be found throughout the world in diverse geological environments,
and have been known since prehistory mainly due to their use in jewelry as
gemstones. They are also relatively hard minerals, a property that makes them
useful for a series of industrial applications, such as in waterjet cutting or
as abrasives.

Generally, the garnets crystallize in a cubic structure (space group
$Ia\overline{3}d$) with chemical composition A3B2(B’C4)3, where the A atoms
are located in the $24c$ dodecahedral sites, the B atoms are in the $16a$
octahedral, and B’ atoms occupy the $24d$ tetrahedral sites. In ternary
garnets, B and B’ sites are occupied by the same chemical element. Around 1950
Yoder and Keith (1951); Bertaut and Forrat (1956); LeCraw  _et al._ (1958)
some rare-earth garnets, especially yttrium-based materials, started to
attract attention. Those garnets have a general formula of RE3B2(BO4)3 where
RE stands for rare-earth and B is a 3d magnetic transition metal (usually
iron) or a group IIIA element. Among these, one of the most used ones is
yttrium aluminum garnet (YAG), Y3Al2(AlO4)3, used as a synthetic simulant to
diamond due to its high refractive index ($>1.8$) Palik (1997). Doped YAGs
with other rare-earth elements have found numerous applications as lasing
media in modern medical laser devices Luke  _et al._ (2019) or in tunable
optical devices Basavad  _et al._ (2020); Sedmidubský  _et al._ (2019); Wako
_et al._ (2016); Liu and Vohra (1994).

Other important compounds, with interesting ferrimagnetic properties, are the
rare-earth iron garnets (RE3Fe2(FeO4)3, RIG). In the RIG structure, five Fe
atoms occupy two different sublattices, and the antiferromagnetic coupling
between sub-lattices and ferromagnetic coupling within the sublattice leads to
a ferrimagnetic configuration. RIGs can display a rather high Curie
temperature (around 560 K Harris (2012)), and some systems exhibit giant
magnetorestriction Sayetat (1986). Moreover, RIGs materials have a band gap
with values around 2.6 to 2.9 eV Gavriliuk  _et al._ (2005); Metselaar and
Larsen (1974). Among these materials yttrium iron garnet (YIG) stands out
because it has an exceptionally low Gilbert damping. YIG has first been used
as bulk material in optical insulators, circulators and Faraday rotators.
Since the last two decades, YIG is also more and more frequently used as thin
film material for spintronic applications Wu and Hoffmann (2013) because it
allows the transmission of spin currents although being an insulator by
itself. In recent years we have witnessed the attempt to replace yttrium by
lanthanides to increase the spin orbit coupling and introduce even
Dzyaloshinskii–Moriya interactions in hybrid systems.

Another interesting group of quaternary garnets is the lithium garnets
LN3M2(LiO4)3, where LN is a lanthanide and M is either Te, Ta, or Nb. With
partial Li-filling of the invasive positions of tetrahedral and distorted
octahedral sites, the stuffed lithium garnets (LN3M2Li2(LiO4)3) have a
promising lithium-ion conductivity and chemical stability, showing potential
as solid electrolytes in Li-batteries Subramanian  _et al._ (2021); Kim  _et
al._ (2020).

The many applications of garnets, and of YIG in particular, has increased the
need for garnets with different properties. Unfortunately the deposition of
garnet thin films with high quality is only possible on garnet substrates due
to the special crystal structure Schmidt  _et al._ (2020); Althammer (2018);
Serga  _et al._ (2010). Moreover for spintronics applications it would be a
huge step forward to be able to pair, for example, thin YIG films in hybrid
structures with other conducting films being either metallic or having such a
low band gap that reasonable electron conductivity is achievable at room
temperature Schmidt  _et al._ (2020); Althammer (2018); Serga  _et al._
(2010). Again to achieve a perfect interface in these structures one would
need to create an all-garnet hybrid which is currently prevented by the
obvious lack of room temperature conducting or metallic garnets. Such a
material class would dramatically extend the applicability of garnet thin
films.

However, despite these diverse applications and the technological relevance of
garnets, most of the research in garnets is confined to oxides Ye  _et al._
(2018), and only a few halides (also called cryolithionites) are known
experimentally Grew  _et al._ (2013). This is can be easily understood, as
oxides are usually simpler to work with under ambient experimental conditions.
Furthermore, a computational high-throughput search of new compositions for
the garnet prototype is challenging, as the garnet cubic primitive unit cell
contains 80 atoms, which is an order of magnitude larger than most structure
prototypes used in recent high-throughput searches such as
(double-)perovskites Schmidt  _et al._ (2021, 2017), (half)-Heuslers Oliynyk
_et al._ (2016), dichalcogenides Fukuda  _et al._ (2021), etc.

Fortunately, with the aid from state-of-the-art machine learning techniques,
the problem of searching through the entire combinatorial chemical space can
be significantly accelerated. With pre-trained machine models, we can filter
millions of compositions according to the predicted stability without
performing costly density functional theory (DFT) calculations for all the
compositions. Nevertheless, DFT validation of the stable compounds is still
necessary as a post-processing step.

In the present paper, we followed such a procedure to explore possible
(meta-)stable compounds beyond oxy-garnets. The rest of the paper is
structured as follows. In section II we explain the machine learning model and
the computational methods we applied. In section III we present the most
interesting crystal phases we uncovered in our work and discuss the potential
applications of the new proposed compounds. Finally we present our conclusions
and an outlook.

##  II Methods

###  II.1 Machine Learning Model

In this work we applied crystal graph attention networks, developed and
pretrained in Ref. Schmidt  _et al._ , 2021, to predict thermodynamically
stable materials. The networks use an attention-based message passing approach
based on the crystal graph representation of the crystal structure. Replacing
the normal distance information that is typically used as edge-representation
in crystal graph networks with solely the graph distance of the atoms to their
neighbors allows for precise predictions of unrelaxed structures. As garnets
crystallize in a cubic structure the list of neighbors and consequently the
graph distances are mostly constant throughout the geometry relaxation. This
removes the need to perform predictions with multiple cell constant ratios.

We study all garnet compositions A3B5C12 considering all elements up to Bi,
excluding only the rare-gases, and spanning a space of around 550k possible
compounds. Obviously, most of these hypothetical compounds are highly unstable
from the thermodynamic point of view. As discussed in Ref. Schmidt  _et al._
(2021), current machine learning approaches for the prediction of
thermodynamic stability suffer from a considerable error due to the strong
bias present in the existing datasets. To circumvent this problem, we use the
following workflow: (i) The machine-learning model is used to predict the
distance to the convex hull of stability for all 550k compounds. At the start
we use the pre-trained machine from Ref. Schmidt  _et al._ (2021); (ii) We
perform DFT geometry optimizations to validate all compounds predicted below
200 meV/atom from the hull; (iii) We add these calculations to a dataset
containing all DFT calculations for garnet systems. (iv) We use transfer
learning to train a new model using this dataset with a
training/validation/testing split of 80%/10%/10%. (v) The cycle is restarted.

In total we repeated the cycle three times. In the first, we performed DFT
calculations for 3320 compounds. The mean absolute error (MAE) of the initial
pre-trained machine was 0.497 eV/atom. This is a very high value, that was
expected as there were very few garnets in the dataset used in Ref. Schmidt
_et al._ (2021), and they spanned a very small chemical space. In the second
cycle we validated 7336 compounds, and the transfer-learning model performed
much better, with a MAE of 0.064 eV/atom. Finally, in the third cycle we
computed 3844 materials with DFT. The final model had a MAE of 0.058 eV/atom,
showing that the tranfer-learning workflow is meaningful and converges
quickly.

###  II.2 DFT Calculations

We perform DFT calculations using the package vasp with PAW Blöchl (1994)
datasets of version 5.2. and with the Perdew-Burke-Ernzerhof Perdew  _et al._
(1996) (PBE) exchange-correlation functional. Following the Materials Project
Jain  _et al._ (2013) recommendations, we use extra on-site corrections for
oxides, fluorides containing Co, Cr, Fe, Mn, Mo, Ni, V, and W. The on-site
corrections are repulsive and correct the d-states by respectively 3.32, 3.7,
5.3, 3.9, 4.38, 6.2, 3.25, and 6.2 eV. A cutoff of 520 eV is applied to the
planewaves, and $\Gamma$-centered $k$-point grids with a uniform density of
1000 k-points per reciprocal atom are used to sample the Brillouin zone. For
geometry optimizations all forces are converged to less than 0.005 eV/Å.

All calculations are performed with spin-polarization, starting from a
ferromagnetic ground-state as is customary in high-throughput searches.
Unfortunately, this means that, in most cases, antiferromagnetic or
ferrimagnetic systems will converge to an incorrect ferromagnetic ground-
state. This is important, in our context, particularly for ferrimagnetic
garnets having a 3d transition metal such as Fe, Ni, Co, Cr, Mn, or V in the B
position. We note that not only the spin-state of these garnets but also of
anti- and ferrimagnetic systems on the convex hull are treated incorrectly.
Consequently, the estimation of $E_{\text{hull}}$ for ferrimagnetic garnets is
far less accurate than for non-magnetic ones. For example, the experimentally
know Gd3Fe5O12 is predicted to be more than 1 eV above the hull according to
Materials Project database Jain  _et al._ (2013), a value that is certainly
grossly overestimated.

To properly estimate $E_{\text{hull}}$ for the ferrimagnetic garnets would
require obtaining the correct magnetic ordering for a portion of the convex
hull as well as for the garnets. This is a complex and computationally
expensive task, that is well beyond the scope of this work. Therefore, we made
the choice to restrict our discussion to systems not containing the 3d metals
mentioned above in the B-site.

It is well-known that the electronic band gaps calculated with the PBE
functional are severely underestimated Borlido  _et al._ (2019). Therefore, to
obtain a more reliable estimation of this important physical property we
performed calculations with the modified Becke-Johnson (mBJ) approximation
Tran and Blaha (2009), as this is by now recognized as one of the most
accurate functional for this task Borlido  _et al._ (2020).

To calculate the averaged carrier effective masses from the interpolated
eigenvalues we follow the approach of Ref. Hautier  _et al._ (2013).
Considering a temperature of 300 K, we deduce the chemical potential required
to reach a reference carrier concentration (1018 cm-3) by using BoltzTraP2. We
use a $k$-point mesh with a regular density of 2000 $k$-points per reciprocal
atom and interpolated the calculated eigenvalues using BoltzTraP2 Madsen and
Singh (2006); Madsen  _et al._ (2018). The calculated averaged carrier
effective masses can be seen as the intrinsic tendency for creating mobile
charge carriers in materials Hautier  _et al._ (2013).

##  III Results and Discussion

###  III.1 Stability and Chemistry

(a)

(b)

(c)

Figure 1: The distribution of the distance to the convex hull ($E_{\text{hull}}$) for (a) all calculated compositions, (b) halides, and (c) chalcogenides/nitrides. The cutoff to filter (meta-)stable systems (100 meV/atom) is shown as dashed vertical line.  Table 1: Number(meta-)stable systems ($N_{\text{stable}}$) below 100 meV/atom from the convex hull of thermodynamic stability and high-throughput success rate ($R$) for each category. Category | $N_{\text{stable}}$ |  $R$ (%)
---|---|---
Sulfides | 70 | 14
Selenides | 68 | 15
Tellurides | 28 | 7
Nitrides | 64 | 35
Fluorides | 62 | 17
Chlorides | 68 | 16
Bromides | 72 | 17
Iodides | 68 | 15
Hydrides | 69 | 8
Total | 569 | 14

In total the machine-learning model predicted around 12300 compositions below
200 meV/atom from the convex hull of thermodynamic stability that were not
present in the Materials Project database Jain  _et al._ (2013) nor in the
Inorganic Crystal Structure Database F. H. Allen and Sievers (1987) (ICSD).
All these calculations can be downloaded from the Materials Cloud repository
Schmidt  _et al._ (2022).

After the DFT validation calculations we re-evaluated the distance to the
convex hull ($E_{\text{hull}}$) of these candidates using the much more
complete convex hull of Ref. Schmidt  _et al._ (2021). The histogram of the
values $E_{\text{hull}}$ is shown in Fig. 1(a). We also separate the systems
into sulfides, selenides, tellurides, nitrides, chlorides, bromides, iodides,
and hydrides. These comprise the majority of all systems found. Most of the
candidates have an $E_{\text{hull}}$ larger than 100 meV/atom, but there are
still more than one thousand (about 9%) compositions below this threshold. The
high-throughput success rate, that we define by the number of compounds that
are within 100 meV/atom from the convex hull divided by the total number of
DFT calculations, stood at 14%, with a maximum of 35% for nitrides and a
minimum of 8% for hydrides. These numbers prove the efficiency of our machine-
learning assisted high-throughout search. The histogram of $E_{\text{hull}}$
for these categories as well as all calculated compositions are shown in Fig.
1.

The distribution of $E_{\text{hull}}$ for all systems follows the typical
skewed Gaussian with the peak located at around 200 meV/atom and a fat tail
that extends beyond 1 eV/atom, in agreement with the MAE errors for our
machine-learning models. Individual distributions for chalcogenides are also
skewed Gaussians peaking at around 200 meV/atom. The amount of (meta-)stable
compounds decreases from sulfides to tellurides, which is also expected as
this is the common trend of stability for the chalcogenides. Unlike the
situation for chalcogenides, the distribution curve for nitrides has multiple
peaks, and shows that there are plenty of potentially stable nitride garnets.
For halides and hydrides the histograms are again skewed Gaussian similar to
those of chalcogenides, but there is no clear trend in what concerns stability
across the group. For hydrides, the total number of systems is much larger,
but a lower percentage of them are (meta-)stable compared to the halides.

The total number of (quasi-)stable systems for each category is listed as
Table 1. A full list of the systems can be found in Table 1 of the
Supplementary Information (SI). We also selected a dozen of them to analyze
more closely in Table 3. Later discussions will mainly focus on these systems.

Table 2: The experimentally known oxy-garnets (not including a 3d metal in the B site), their ICSD ID, Materials Project ID, and predicted meta-stable counterparts with different C anions. The distance to the hull calculated with DFT is in parentheses (in meV/atom). Formula | ICSD ID | MP ID | Counterparts
---|---|---|---
Y3Al5O12 | 20090, 41144, 41145, 67102, 67103, | mp-3050 | Y3Al5S12(-4); Y3Al5Se12(36); Y3Al5Te12(84)
93634, 93635, 170157, 170158, 236589,
280104, 17687, 17688, 17689, 17690,
74607, 31496
La3Al5O12 |  | mp-780432 |  La3Al5S12(-4); La3Al5Se12(19); La3Al5Te12(49)
Eu3Al5O12 | 245326 | mp-21757 |  Eu3Al5S12(39); Eu3Al5Se12(93)
Tb3Al5O12 | 33602 | mp-14387 |  Tb3Al5S12(-4); Tb3Al5Se12(29); Tb3Al5Te12(79)
Er3Al5O12 | 170147, 280606, 170146, 62615 | mp-3384 |  Er3Al5S12(6); Er3Al5Se12(50); Er3Al5Te12(98)
Gd3Al5O12 | 192184 | mp-14133 |  Gd3Al5S12(-16); Gd3Al5Se12(21); Gd3Al5Te12(66)
Ho3Al5O12 | 409390, 33603 | mp-14388 |  Ho3Al5S12(-1); Ho3Al5Se12(44); Ho3Al5Te12(92)
Lu3Al5O12 | 259144, 17789, 182354 | mp-14132 |  Lu3Al5S12(28); Lu3Al5Se12(70)
Y3Ga5O12 | 80148, 14343, 185862, 23852 | mp-5444 |  Y3Ga5S12(47); Y3Ga5Se12(81)
La3Ga5O12 |  | mp-780561 |  La3Ga5S12(36); La3Ga5Se12(51); La3Ga5Te12(100)
Tb3Ga5O12 | 20831, 84875, 184934 | mp-5965 |  Tb3Ga5S12(37); Tb3Ga5Se12(75)
Sm3Ga5O12 | 9236, 84873, 291192 | mp-5800 |  Sm3Ga5S12(29); Sm3Ga5Se12(64); Sm3Ga5Te12(93)
Nd3Ga5O12 | 84872 | mp-15239 |  Nd3Ga5S12(26); Nd3Ga5Se12(54); Nd3Ga5Te12(96)
Gd3Ga5O12 | 9237, 37145, 192181, 84874, 184931 | mp-2921 |  Gd3Ga5S12(32); Gd3Ga5Se12(67); Gd3Ga5Te12(93)
Lu3Ga5O12 | 23850 | mp-14134 |  Lu3Ga5S12(74)
Dy3Ga5O12 | 409391 | mp-15576 |  Dy3Ga5S12(49); Dy3Ga5Se12(81)
Er3Ga5O12 | 9238 | mp-12236 |  Er3Ga5S12(55); Er3Ga5Se12(95)
Ho3Ga5O12 | 409390 | mp-15575 |  Ho3Ga5S12(52); Gd3Ga5Se12(89)

From our calculations we recover the majority of the oxide (and halide)
garnets that are already known experimentally, but we also obtain a wealth of
different compounds not present in available databases. Many of these systems
are oxides, that have been the subject of a recent high-throughput search Ye
_et al._ (2018) with results similar to ours. Interestingly, we also find a
wealth of other chalcogenides, nitrides, halides, and even hydrides as shown
in Figs. 1(b),1(c) and in Table1.

From the stable compounds, several are closely related to the oxide garnets by
the chemical substitution of oxygen by another chalcogen. These garnets along
with their counterparts are presented in Table 2. The existence of such
compounds is expected due to the chemical similarity among chalcogens. In
those chalcogenides, the dodecahedral sites (site A in A3B5S12) are mostly
occupied by rare earth elements, and according to the element occupying the
octahedral and tetrahedral sites (site B), the chalcogenides can be further
divided into several categories. The numbers of (meta-)stable systems for each
category decrease in the following order: occupying B with group IIIA elements
(Al, Ga, In, Tl), group IVA elements (Ge, Sn, Pb), group VA elements (As, Sb,
Bi), and transition metals (Ag, Cu, Sc, Ti). The preference of group IIIA
elements for site B can be understood by simple charge compensation arguments.
The most common oxidation state of the chalcogens is $-2$, while the rare
earth elements in sites A are $+3$: the composition reaches the “balanced” or
“compensated” state if B is in the oxidation state $+3$. Moreover, for
balanced A${}^{\text{III}}_{3}$B${}^{\text{III}}_{5}$C${}^{\text{-II}}_{12}$
chalcogenide compositions, one would expect the compounds to be
semiconductors. This is indeed what we find (see an example in section III.2).

Besides the chalcogenides, we discover 64 stable nitrides as seen in Table 1.
As discussed above, to reach a balanced oxidation state, elements with higher
oxidation state should be favored to combine with $-3$ oxydation state of
nitrogen. Indeed for nitride garnets A3B5N12, the position A is mostly
occupied by $+2$ or $+3$ chemical elements, and for position B the majority of
meta-stable systems have elements with oxidation state $+6$ (Mo and W).
Although a $+6$ element is required to achieve a balanced state, we also find
that relatively stable compounds are possible for $+5$ (Nb and Ta), and $+7$
(Re). We could argue that the balanced nitrides should be semiconducting.
However, due to the low electronegativity of nitrogen the gap may close (see
one such example in section III.2).

We can also identify several halides and hydrides from Figure 1(c). However,
halogens (and hydrogen in hydrides) have an oxidation state of $-1$, which
makes it more difficult to reach a balanced oxidation state with the A3B5C12
stoichiometry. One viable scheme is with an A element that is $+1$, while the
B elements in $16a$ and $24d$ Wycoff positions are respectively $+1$ and $+3$,
i.e. having the form of
A${}^{\text{I}}_{3}$B${}^{\text{I}}_{2}$(B${}^{\text{III}}$C${}^{\text{-I}}_{4}$)3.
Chemical elements exhibiting both $+1$ and $+3$ oxidation states are quite
rare. Nevertheless, we still discovered some meta-stable semiconducting halide
garnets, for example K3In5F12.

There are other compositions that do not belong to any of the discussed
classes, such as for example Sr12Zn3H5. Most of them are “inverted”-garnets,
i.e. with cations instead of anions occupying the C-sites, and have a
comparatively higher $E_{\text{hull}}$ than regular garnets. Furthermore, only
a few anti-garnets with the A and B sites both occupied by the chemical
elements of the nitrogen group could potentially become oxidation state
balanced.

We have to again emphasize that, in order to form (meta-)stable or
insulating/semiconducting compounds, charge compensation is neither a
necessary nor a sufficient condition, and we find many exceptions in Table S1
in the Supplementary Information. However, it gives us a simple, intuitive
argument to understand why a system is stabilized or has an electronic band
gap. Furthermore, we have to keep in mind that uncompensated systems might be
further stabilized through defects, such as vacancies. To simplify our
discussion, we leave such possibilities to future works, and focus in the
following on (meta-)stable regular garnets systems which could have balanced
charges, specifically chalcogenides (except oxides), halides, hydrides, and
nitrides.

Table 3: The calculated lattice constant ($a$, in Å), band gap calculated with PBE (Gap${}^{\text{PBE}}$ ) and MBJ (Gap${}^{\text{MBJ}}$) functional (in unit of eV)), distance to the convex hull ($E_{\text{hull}}$ in meV/atom), effective electron ($m_{\text{e}}^{*}$) and hole ($m_{\text{h}}^{*}$) masses (in unit of $m_{\text{e}}^{0}$), for some selected (meta-)stable sulfide and nitride garnets, data for Y3Al5O12 is also listed for comparison. Formula | $a$ | Gap${}^{\text{PBE}}$ | Gap${}^{\text{MBJ}}$ | $E_{\text{hull}}$ | $m_{\text{e}}^{*}$ | $m_{\text{h}}^{*}$
---|---|---|---|---|---|---
Y3Al5O12 | 12.125 | 4.53 | 6.12 | 0 | 1.3 | 6.8
Y3Al5S12 | 14.932 | 2.08 | 3.00 | 0 | 0.7 | 1.5
Y3Ga5S12 | 15.073 | 1.32 | 2.45 | 46 | 0.5 | 3.2
Y3In5S12 | 15.709 | 1.35 | 2.38 | 46 | 0.5 | 2.6
Y3Al5Se12 | 15.752 | 1.46 | 2.15 | 36 | 0.5 | 1.5
Y3Al5Te12 | 17.120 | 0.60 | 1.07 | 84 | 0.3 | 2.7
Y3Ge5S12 | 15.324 | 0.00 | 0.00 | 81 | – | –
Ca3W5N12 | 12.928 | 0.93 | 1.57 | 0 | 0.9 | 3.2
Ca3Re5N12 | 12.896 | 0.00 | 0.00 | 9 | – | –
La3Nb5N12 | 13.297 | 0.00 | 0.00 | 99 | – | –
K3In5F12 | 14.898 | 2.41 | 3.39 | 59 | 0.8 | 6.1
K3In5I12 | 19.977 | 0.00 | 0.00 | 61 | – | –
Mg3Rh5H12 | 11.008 | 0.00 | 0.00 | 51 | – | –
Y3Rh5H12 | 11.581 | 0.00 | 0.00 | 92 | – | –

###  III.2 Electronic Structure

We illustrate the possible electronic structures of our garnets through a few
selected examples depicted in Figs. 2–4.

Generally the electron states around the Fermi level in garnets can be
classified in three categories: (i) from s–p, p–p, or even d–p bonding
orbitals between B and C atoms, (ii) from corresponding anti-bonding between B
and C atoms, and, (iii) from d-electrons from A atoms if these are d-block or
f-block metals. From the simplest tight-binding model, and as expected from
$k\cdot p$ theory, we know that the larger the difference between the
electronegativity of B and C and the shorter B–C bond-length, the larger the
separations between states (i) and (ii). For the d (or f) block elements
occupying the A-site, the position of states (iii) can be between (i) and
(ii), overlapping with, or even above the latter. In charge compensated
situations, the bands (i) are completely filled, while (ii) and (iii) (if
applicable) bands are empty, resulting in insulating/semiconducting systems.
Otherwise, depletion of bands (i) or filling of bands (ii) or (iii) can
happen, leading to metallic systems. Although many deviations from this simple
tight-binding picture appear, we will see how these general patterns are
useful to understand the band-structures.

(a) Y3Al5O12

(b) Y3Al5S12

(c) Y3Ge5S12

(d) Y3Ga5S12

Figure 2: Calculated mBJ electronic band structures for selected chalgogenide
garnets. The Fermi level is set at zero.

(a) Ca3W5N12

(b) Ca3Re5N12

Figure 3: Calculated mBJ electronic band structures for selected nitride
garnets. The Fermi level is set at zero.

For chalcogenide garnets we chose as representative examples Y3Al5S12,
Y3Ge5S12, and Y3Ga5S12, while more examples can be found in the SI. The band-
structures of these compounds are shown in Fig 2, together with Y3Al5O12 for
comparison. The oxidation-state balanced Y3Al5CH12 (CH = S, Se, Te) garnets
are semiconductors as expected. For Y3Al5S12 the band structure is shown in
Fig. 2(b). Similar to its oxide counterpart in Fig. 2(a), the Y–$d$ states
dominate the conduction bands (CB) slightly hybridizing with Al(p) and
S(p)-states. These bands are from the type (iii) states as discussed above.
The valence bands (VB) around Fermi level are mainly composed of the localized
anionic $p$–states of type (i), also representing the typical situation
described above. From Y3Al5O12 to Y3Al5Te12 (see Fig. S1 in SI), following the
decreasing trend of electronegativity for chalcogens, the band gap shrinks and
the band edges become more dispersed.

The elements occupying the B position in charge-compensated chalcogenides also
have a crucial effect on the electronic structure. For example, Y3B5S12 (B =
Al, Ga, In, Tl, see Fig. S1 in the SI) all have a direct gap at $\Gamma$ with
the sole exception of Y3Al5S12 which has an indirect H–$\Gamma$ gap. Moreover,
going from Al to Tl the gap decreases, while the contibution of the s-states
from the B atoms to the bottom of the conduction bands increases, leading to
more extended bands and to lower effective electron masses. Furthermore, for
Y3Tl5S12 we can see that the Y(d) bands are above the Tl(s)-S(s) displaying
more disperse anti-bonding mixing, an resulting in very low electron effective
masses ($m^{*}_{\text{e}}=0.19$).

For chalcogenides with unbalanced oxidation states, such as Y3Ge5S12, similar
features can also be observed (see Fig. 2(c)). Between the empty Y(d)-S(p)
hybridization states of type (iii) and the full S(p) dominated type (i) bands
are the states from Ge(s) and S(p) anti-bonding of type (ii). As discussed
above they are partially occupied, so the system is metallic. It is possible
that the system might re-establish a balanced oxidation state by creating Ge
vacancies thus becoming semiconducting, though a detailed investigation of
such a possibility is beyond the scope of the present paper.

(a) K3In5F12

(b) K3In5I12

(c) Mg3Rh5H12

(d) Y3Rh5H12

Figure 4: Calculated mBJ electronic band structures for selected halide and
hydride garnets. The Fermi level is set at zero.

Some representatives of the band structures for nitrides garnets are shown in
Fig. 3. The oxidation state compensated Ca3W5N12 is semiconducting. The anti-
bonding $p-d$ bands forming the CB separate also into two manifolds: the lower
part is mainly constructed from the tungsten on the tetrahedral sites, and the
upper part mainly from the octahedral W. The VB are, as usual, mostly composed
of N(p) states. For charge unbalanced nitrides (such as Ca3Re5N12) that have
more valence electrons from Re compared to W, the Fermi level passes through
the partially occupied $p-d$ anti-bonding bands and the system becomes
metallic. However, the separation between $p-d$ bonding and anti-bonding bands
can still be seen as well as the double manifolds of $p-d$ anti-bonding
states. When a chemical element with less valence electrons replaces W, for
example, in La3Nb5N12 (see Fig. S1 in SI), the minority spin channel of $p-d$
bonding states is partially empty, and the system actually becomes a half-
metal.

For halides, we show the band structure of K3In5F12 and K3In5I12 as examples.
Apparently, the former compound has unbalance oxidation states, however it is
semiconducting. The Bader net charges of the In at octahedral and tetrahedral
sites are +1 and +2.4, respectively. The In in the B-site is therefore in the
+1 and +3 oxidation states, reaching charge compensation as discussed above.
The top of the valence bands are mainly composed of In(s)-F(p) anti-bonding
states, which separate from the lower valance bands, formed by localized
F(p)-In(p) bonding states, by around 4 eV. Unlike in most sulfides and
nitrides, where d-states from A atoms dominate the bottom of the conduction
bands, in K3In5F12 the bottom of the CB is mainly formed by In(p)-F(p) anti-
bonding hybridized states. This is because the K(s)-F(p) anti-bonding states
have a much higher energy than d-p anti-bonding states in those sulfides and
nitrides. When replacing F with I, the Bader net charges for the In atoms in
octahedral and tetrahedral sites are respectively +0.7 and +1.0, showing that
the +1/+3 oxidation states is not possible in K3In5I12, and the system becomes
metallic. The top valence band formed by In(p)-I(s) anti-bonding states
separates in two manifolds. The upper part comes from the In atoms in
tetrahedral sites and the lower belongs to the octahedral ones. This can be
partially explained by the fact that the In–I bond length is larger in the
latter.

In Fig 4 we show the band structure of selected examples of hydride garnets.
For both Mg3Rh5H12 and Y3Rh5H12 the separation between H(s)-Rh(d) mixed bands
and cation Rh(d) dominated bands are still present, similar to the situations
in other uncompensated garnets. These two systems are also both charge
uncompensated, the “extra” electrons partially filling the Rh(d) dominated
bands, leading to metallic systems. Another way to re-establish compensation
might be to force more H atoms to occupy the interstitial sites, but again, we
leave exploration of this possibility to future investigations.

##  IV Conclusions

We performed a machine-learning assisted high-throughput investigation of
ternary garnets. We concentrated in non-oxides (that have been studied
previously) and in ferromagnetic or paramagetic compounds. We find a wealth of
systems on the convex hull (i.e., thermodynamically stable) or close to it.
This includes chalcogenides (with the stability decreasing from S to Te),
nitrides, halides, hydrides, etc. We also found other possibilities, such as
“inverted” garnets, but these were slighlty less stable than the conventional
phase. The materials tend to be semiconducting/insulating when the composition
is charge compensated, otherwise we obtain metallic ground-states. The latter
ones could be especially relevant as, to out knowledge, no garnets conducting
at room temperature are know. Some of the metallic garnets even have lattice
constants that are suitable to create hetero-structures with YIG.

A few chalcogenide garnets, in particular the sulfides, are thermodynamic
stable, and are straightforward generalizations of common oxide garnets. Band-
gaps are as expected considerably smaller for the sulfides, and decrease
further across the periodic group. This opens up the possibility to engineer
the band gap of garnets by anionic alloying, from the extreme ultaviolet of
the oxide phase to UV-A regime or even into the visible. We predict several
nitride systems that have interesting electronic properties due to the
presence of transition metals in very high charge states. In view of the
recent synthesis of two exotic nitride perovskites Talley  _et al._ (2021);
Kloß  _et al._ (2021) that were predicted Sarmiento-Pérez  _et al._ (2015);
Flores-Livas  _et al._ (2019) with a method similar to the one used in this
paper, we are confident that also nitride garnets are accessible
experimentally. Finally, we find a few semiconducting halides where the
chemical element occupying the octahedral and the dodecahedral site is in two
different charge states.

Above all, we believe that our work proves that an exhaustive survey of the
ternary, and perhaps also of the quaternary, space of materials is now
accessible to high-throughput studies, even for large and complex unit cells.
This is made possible by machine learning methods, that already achieved an
outstanding maturity in the short time since their first appearance, and that
are reaching an unprecedented accuracy. We expect these methods to further
accelerate the discovery of new materials with exceptional properties.

##  V Acknowledgements

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V.
(www.gauss-centre.eu) for funding this project by providing computing time on
the GCS Supercomputer SuperMUC-NG at the Leibniz Supercomputing Centre under
the project pn25co.

## References

  * Yoder and Keith (1951) HS Yoder and ML Keith, “Complete substitution of aluminum for silicon: The system 3MnO·Al2O3· 3SiO2—3Y2O3·5Al2O3,” Am. Mineral. 36, 519–533 (1951).
  * Bertaut and Forrat (1956) F Bertaut and F Forrat, “Structure of ferrimagnetic rare-earth ferrites,” CR Acad. Sci 242, 382 (1956).
  * LeCraw  _et al._ (1958) R. C. LeCraw, E. G. Spencer, and C. S. Porter, “Ferromagnetic resonance line width in yttrium iron garnet single crystals,” Phys. Rev. 110, 1311–1313 (1958).
  * Palik (1997) ED Palik, in  _Handbook of Optical Constants of Solids_ (Elsevier, 1997).
  * Luke  _et al._ (2019) Alexander Maninagat Luke, Simy Mathew, Maram Majed Altawash, and Bayan Mohammed Madan, “Lasers: A review with their applications in oral medicine,” J Lasers Med Sci. 10, 324 (2019).
  * Basavad  _et al._ (2020) M. Basavad, H. Shokrollahi, H. Ahmadvand, and S.M. Arab, “Structural, magnetic and magneto-optical properties of the bulk and thin film synthesized cerium- and praseodymium-doped yttrium iron garnet,” Ceram. Int. 46, 12015–12022 (2020).
  * Sedmidubský  _et al._ (2019) David Sedmidubský, Vít Jakeš, Kateřina Rubešová, Pavla Nekvindová, Tomáš Hlásek, Roman Yatskiv, and Pavel Novák, “Magnetism and optical properties of Yb3Al5O12 hosted Er3+ – experiment and theory,” J. Alloys Compd. 810, 151903 (2019).
  * Wako  _et al._ (2016) A. H. Wako, F. B. Dejene, and H. C. Swart, “Effect of Ga3+ and Gd3+ ions substitution on the structural and optical properties of Ce3+-doped yttrium aluminium garnet phosphor nanopowders,” Luminescence 31, 1313–1320 (2016).
  * Liu and Vohra (1994) Jun Liu and Yogesh K. Vohra, “Sm:YAG optical pressure sensor to 180 GPa:calibration and structural disorder,” Appl. Phys. Lett. 64, 3386–3388 (1994).
  * Harris (2012) Vincent G. Harris, “Modern microwave ferrites,” IEEE Trans. Magn. 48, 1075–1104 (2012).
  * Sayetat (1986) F. Sayetat, “Huge magnetostriction in Tb3Fe5O12, Dy3Fe5O12, Ho3Fe5O12, Er3Fe5O12 garnets,” J. Magn. Magn. Mater. 58, 334–346 (1986).
  * Gavriliuk  _et al._ (2005) A. G. Gavriliuk, V. V. Struzhkin, I. S. Lyubutin, and I. A. Trojan, “Irreversible electronic transition with possible metallization in Y3Fe5O12 at high pressure,” JETP Lett. 82, 603–608 (2005).
  * Metselaar and Larsen (1974) R. Metselaar and P.K. Larsen, “High-temperature electrical properties of yttrium iron garnet under varying oxygen pressures,” Solid State Commun. 15, 291–294 (1974).
  * Wu and Hoffmann (2013) Mingzhong Wu and Axel Hoffmann, “Recent advances in magnetic insulators-from spintronics to microwave applications,” (2013).
  * Subramanian  _et al._ (2021) Kannan Subramanian, George V. Alexander, K. Karthik, Srabani Patra, M.S. Indu, O.V. Sreejith, Raja Viswanathan, Janani Narayanasamy, and Ramaswamy Murugan, “A brief review of recent advances in garnet structured solid electrolyte based lithium metal batteries,” J. Energy Storage 33, 102157 (2021).
  * Kim  _et al._ (2020) Abin Kim, Seungjun Woo, Minseok Kang, Heetaek Park, and Byoungwoo Kang, “Research progresses of garnet-type solid electrolytes for developing all-solid-state li batteries,” Front. Chem. 8 (2020), 10.3389/fchem.2020.00468.
  * Schmidt  _et al._ (2020) Georg Schmidt, Christoph Hauser, Philip Trempler, Maximilian Paleschke, and Evangelos Th. Papaioannou, “Ultra thin films of yttrium iron garnet with very low damping: A review,” Phys. Status Solidi B 257, 1900644 (2020).
  * Althammer (2018) Matthias Althammer, “Pure spin currents in magnetically ordered insulator/normal metal heterostructures,” J. Phys. D: Appl. Phys. 51, 313001 (2018).
  * Serga  _et al._ (2010) A A Serga, A V Chumak, and B Hillebrands, “YIG magnonics,” J. Phys. D: Appl. Phys. 43, 264002 (2010).
  * Ye  _et al._ (2018) Weike Ye, Chi Chen, Zhenbin Wang, Iek-Heng Chu, and Shyue Ping Ong, “Deep neural networks for accurate predictions of crystal stability,” Nat. Commun. 9 (2018), 10.1038/s41467-018-06322-x.
  * Grew  _et al._ (2013) E. S. Grew, A. J. Locock, S. J. Mills, I. O. Galuskina, E. V. Galuskin, and U. Halenius, “Nomenclature of the garnet supergroup,” Am. Mineral. 98, 785–811 (2013).
  * Schmidt  _et al._ (2021) Jonathan Schmidt, Love Pettersson, Claudio Verdozzi, Silvana Botti, and Miguel A. L. Marques, “Crystal graph attention networks for the prediction of stable materials,” Sci. Adv. 7, eabi7948 (2021).
  * Schmidt  _et al._ (2017) Jonathan Schmidt, Jingming Shi, Pedro Borlido, Liming Chen, Silvana Botti, and Miguel A. L. Marques, “Predicting the thermodynamic stability of solids combining density functional theory and machine learning,” Chem. Mater. 29, 5090–5103 (2017).
  * Oliynyk  _et al._ (2016) Anton O. Oliynyk, Erin Antono, Taylor D. Sparks, Leila Ghadbeigi, Michael W. Gaultois, Bryce Meredig, and Arthur Mar, “High-throughput machine-learning-driven synthesis of full-heusler compounds,” Chem. Mater. 28, 7324–7331 (2016).
  * Fukuda  _et al._ (2021) Masahiro Fukuda, Jingning Zhang, Yung-Ting Lee, and Taisuke Ozaki, “A structure map for AB2 type 2D materials using high-throughput DFT calculations,” Mater. Adv. 2, 4392–4413 (2021).
  * Blöchl (1994) P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B 50, 17953–17979 (1994).
  * Perdew  _et al._ (1996) John P Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865 (1996).
  * Jain  _et al._ (2013) Anubhav Jain, Shyue Ping Ong, Geoffroy Hautier, Wei Chen, William Davidson Richards, Stephen Dacek, Shreyas Cholia, Dan Gunter, David Skinner, Gerbrand Ceder, and Kristin a. Persson, “The materials project: A materials genome approach to accelerating materials innovation,” APL Mater. 1, 011002 (2013).
  * Borlido  _et al._ (2019) Pedro Borlido, Thorsten Aull, Ahmad W. Huran, Fabien Tran, Miguel A. L. Marques, and Silvana Botti, “Large-scale benchmark of exchange–correlation functionals for the determination of electronic band gaps of solids,” J. Chem. Theory Comput. 15, 5069–5079 (2019).
  * Tran and Blaha (2009) Fabien Tran and Peter Blaha, “Accurate band gaps of semiconductors and insulators with a semilocal exchange-correlation potential,” Phys. Rev. Lett. 102, 226401 (2009).
  * Borlido  _et al._ (2020) Pedro Borlido, Jonathan Schmidt, Ahmad W. Huran, Fabien Tran, Miguel A. L. Marques, and Silvana Botti, “Exchange-correlation functionals for band gaps of solids: benchmark, reparametrization and machine learning,” Npj Comput. Mater. 6, 96 (2020).
  * Hautier  _et al._ (2013) Geoffroy Hautier, Anna Miglio, Gerbrand Ceder, Gian-Marco Rignanese, and Xavier Gonze, “Identification and design principles of low hole effective mass p-type transparent conducting oxides,” Nat. Commun. 4 (2013), 10.1038/ncomms3292.
  * Madsen and Singh (2006) Georg K.H. Madsen and David J. Singh, “Boltztrap. a code for calculating band-structure dependent quantities,” Comput. Phys. Commun 175, 67 – 71 (2006).
  * Madsen  _et al._ (2018) Georg KH Madsen, Jesús Carrete, and Matthieu J Verstraete, “BoltzTraP2, a program for interpolating band structures and calculating semi-classical transport coefficients,” Comput. Phys. Commun. 231, 140–145 (2018).
  * F. H. Allen and Sievers (1987) G. Gergerhoff F. H. Allen and R. Sievers, eds.,  _Crystallographic databases_ (International Union of Crystallography, Chester, 1987).
  * Schmidt  _et al._ (2022) Jonathan Schmidt, Hai-Chen Wang, Georg Schmidt, and Miguel A. L. Marques, “Machine learning guided high-throughput search of non-oxide garnets,” Materials Cloud (2022), 10.24435/materialscloud:gm-n0.
  * Talley  _et al._ (2021) Kevin R. Talley, Craig L. Perkins, David R. Diercks, Geoff L. Brennecka, and Andriy Zakutayev, “Synthesis of LaWN3 nitride perovskite with polar symmetry,” Science 374, 1488–1491 (2021).
  * Kloß  _et al._ (2021) Simon D. Kloß, Martin L. Weidemann, and J. Paul Attfield, “Preparation of bulk-phase nitride perovskite LaReN3 and topotactic reduction to LaNiO2-type LaReN2,” Angew. Chem., Int. Ed. 60, 22260–22264 (2021).
  * Sarmiento-Pérez  _et al._ (2015) Rafael Sarmiento-Pérez, Tiago F. T. Cerqueira, Sabine Körbel, Silvana Botti, and Miguel A. L. Marques, “Prediction of stable nitride perovskites,” Chem. Mater. 27, 5957–5963 (2015).
  * Flores-Livas  _et al._ (2019) José A Flores-Livas, R Sarmiento-Pérez, Silvana Botti, Stefan Goedecker, and Miguel A L Marques, “Rare-earth magnetic nitride perovskites,” J. Phys.: Mater. 2, 025003 (2019).