DPUR/TH/72

October, 2021

Scale Invariance and Dilaton Mass

Ichiro Oda111 E-mail address<EMAIL_ADDRESS>
Department of Physics, Faculty of Science, University of the Ryukyus,

Nishihara, Okinawa 903-0213, Japan

###### Abstract

We consider a generic scale invariant scalar quantum field theory and its
symmetry breakdown. Based on the dimension counting identity, we give a
concise proof that dilaton is exactly massless at the classical level if scale
invariance is broken spontaneously. On the other hand, on the basis of the
generalized dimension counting identity, we prove that the dilaton becomes
massive at the quantum level if scale invariance is explicitly broken by
quantum anomaly.

It is pointed out that a subtlety occurs when scale invariance is
spontaneously broken through a scale invariant regularization method where the
renormalization scale is replaced with the dilaton field. In this case, the
dilaton remains massless even at the quantum level after spontaneous symmetry
breakdown of scale symmetry, but when the massless dilaton couples non-
minimally to the Einstein-Hilbert term and is applied for cosmology, it is
phenomenologically ruled out by solar system tests unless its coupling to
matters is much suppressed compared to the gravitational interaction.

##  1 Introduction

A dilaton is defined as a Nambu-Goldstone boson of scale (or dilatation)
invariance and appears ubiquitously in various areas of modern theoretical
physics such as string theory [1] and scalar-tensor gravity [2]. Since scale
symmetry prohibits the appearance of dimensionful parameters in an action, it
is not a symmetry in our world. Thus, in any case scale invariance must be
broken spontaneously or explicitly by quantum anomaly.

One of the important problems associated with the dilaton is the magnitude of
its mass. If the dilaton is exactly massless, it mediates a long-range force
between massive objects as in the Newtonian force. This fact imposes a
phenomenological constraint on parameters in a theory including the dilaton
such as the Brans-Dicke model [3] since the long-range force stemming from the
massless dilaton could affect the perihelion advance of Mercury, for instance.
On the other hand, if the dilaton has a small mass, there might be an
observational possibility of a non-Newtonian force whose force range depends
on the magnitude of the mass and the detail of a model [4].

The aims of this article are two-fold. One aim is related to a recent paper
discussing a no-go theorem that a massless dilaton cannot exist in
renormalizable, unitary scalar quantum field theory (QFT) in four dimensions
[5]. In this article, we wish to present an alternative and very concise proof
of this no-go theorem on the basis of the dimension counting identity and its
anomalous generalization. One of attractive points in our proof is that we can
deal with non-renormalizable QFT as well.

The other aim is to point out a subtlety of the massless dilaton when it is
applied for cosmology. In the manifestly scale invariant regularization method
[6]-[10], scale invariance is spontaneously broken while maintaining quantum
scale invariance, that is, no quantum anomaly. In this case, the dilaton is
precisely massless even at the quantum level, but in order to carry out the
scale invariant regularization method, there appear some problems to be solved
in future, those are, non-renormalizability [11], fine-tuning of coupling
constants [12] and treatment of the unbroken phase [9]. In this article, we
wish to point out that there is another problem when the massless dilaton is
applied for cosmology since the scalar force generated by the massless dilaton
could in general affect the solar system tests and the late-time cosmology, so
it must be screened by the still-unknown mechanism in order to be consistent
with observations and predictions obtained from general relativity.

We close this section with an overview of this article. In Section 2, we
discuss classical scale invariance and show that the dilaton mass as well as
the cosmological constant are identically vanishing at the stationary point of
a potential. In Section 3, we work with a quantum theory with scale anomaly.
It is then easy to show that the dilaton, which is massless at the classical
level, becomes massive owing to the anomaly. In Section 4, we consider the
case where scale invariance is spontaneously broken and quantum scale
invariance is maintained via the scale invariant regularization method. The
final section is devoted to the conclusion.

##  2 Classical scale invariance

In this section, we will present a concise proof that the dilaton is precisely
massless when a theory has strict scale invariance at the classical level.
This proof is based on the dimension counting identity of a potential and
holds for both classical and quantum scale invariances. However, as will be
argued in a later section, one has to pay additional attention to the case of
the quantum scale invariance.

Let us start with a classically scale invariant theory with $n$ real scalar
fields $\phi_{i}(i=1,2,\cdots,n)$222We use a Minkowskian metric of signature
$(-+++)$.:

$\displaystyle{\cal{L}}=-\frac{1}{2}\sum_{i=1}^{n}\partial_{\mu}\phi_{i}\partial^{\mu}\phi_{i}-V(\phi),$
(1)

where $V(\phi)$ is a generic scale invariant potential of scalar fields
$\phi_{i}$. The classical scale invariance then leads to the dimension
counting identity:

$\displaystyle\sum_{i=1}^{n}\phi_{i}\frac{\partial
V(\phi)}{\partial\phi_{i}}=4V(\phi).$ (2)

At the stationary point of the potential $\phi_{i}=v_{i}$, we have

$\displaystyle\left.\frac{\partial
V(\phi)}{\partial\phi_{i}}\right|_{\phi=v}=0.$ (3)

Together with Eq. (3), Eq. (2) implies that the vacuum energy, or equivalently
the cosmological constant when the theory couples to gravity, is vanishing at
the classical level:

$\displaystyle V(v)=0.$ (4)

Taking the partial derivative of Eq. (2) with respect to the scalar fields
$\phi_{j}$, we obtain

$\displaystyle\sum_{i=1}^{n}\phi_{i}\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}=3\frac{\partial
V(\phi)}{\partial\phi_{j}}.$ (5)

At the stationary point of the potential, this equation reads

$\displaystyle\sum_{j=1}^{n}\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}v_{j}=0.$
(6)

In a similar manner, one can obtain a series of equations holding at the
stationary point

$\displaystyle\sum_{j,k}\left.\frac{\partial^{3}V(\phi)}{\partial\phi_{i}\partial\phi_{j}\partial\phi_{k}}\right|_{\phi=v}v_{j}v_{k}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\sum_{j,k,l}\left.\frac{\partial^{4}V(\phi)}{\partial\phi_{i}\partial\phi_{j}\partial\phi_{k}\partial\phi_{l}}\right|_{\phi=v}v_{j}v_{k}v_{l}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\sum_{m}\left.\frac{\partial^{5}V(\phi)}{\partial\phi_{i}\partial\phi_{j}\partial\phi_{k}\partial\phi_{l}\partial\phi_{m}}\right|_{\phi=v}v_{m}$
$\displaystyle=$ $\displaystyle 0,$ (7)

and so on.

At this point, let us focus our attention on Eq. (6). The existence of a
nonzero solution $v_{i}\neq 0$ leads us to the equation

$\displaystyle\det\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}=0.$
(8)

For simplicity, let us take $v_{i}$ to

$\displaystyle v_{i}=(0,0,\cdots,0,v).$ (9)

Then, Eq. (6) states that

$\displaystyle\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{n}}\right|_{\phi=v}=0.$
(10)

Since we assume that the potential $V(\phi)$ is a generic one in the sense
that

$\displaystyle\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{I}\partial\phi_{J}}\right|_{\phi=v}\neq
0,$ (11)

with being $I,J=1,\cdots,n-1$, we can obtain the result

$\displaystyle{\rm
rank}\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}=n-1.$
(12)

We can now make clear the meaning of Eq. (12) by expanding the potential
around $v_{i}$ as follows:

$\displaystyle V(\phi)=V(v)+\sum_{i}\left.\frac{\partial
V(\phi)}{\partial\phi_{i}}\right|_{\phi=v}\varphi_{i}+\frac{1}{2}\sum_{i,j}\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}\varphi_{i}\varphi_{j}+\cdots,$
(13)

where $\dots$ describes higher-order terms in $\varphi_{i}$ and we have
assumed that

$\displaystyle\varphi_{i}\equiv\phi_{i}-v_{i},\qquad|\varphi_{i}|\ll 1.$ (14)

With the help of Eqs. (3) and (4), Eq. (13) can be cast to the form

$\displaystyle V(\phi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i,j=1}^{n}\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}\varphi_{i}\varphi_{j}+\cdots$
(15) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{I=1}^{n-1}\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{I}^{2}}\right|_{\phi=v}\varphi_{I}^{2}+\cdots,$

where in the last equality the symmetric matrix
$\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}$
with the rank (12) is transformed to an $(n-1)\times(n-1)$ diagonal matrix by
a suitable orthogonal matrix.333The sufficient condition such that the
stationary point $\phi_{i}=v_{i}$ is the local minimum of a potential is that
the Hessian matrix
$\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{I}\partial\phi_{J}}\right|_{\phi=v}(I,J=1,2,\cdots,n-1)$,
is the positive definite matrix. The final expression in Eq. (15) clearly
shows that spontaneous symmetry breakdown (SSB) of scale invariance occurs,
thereby yielding $n-1$ massive scalar fields and one massless scalar field,
which is nothing but the dilaton predicted by the Nambu-Goldstone (NG)
theorem. It is of interest that all massless scalar fields except for the
dilaton become massive as a result of the SSB of scale invariance.

Now we will present two simple examples, one of which has a property (11) and
we therefore obtain $n-1$ massive scalar fields and one massless dilaton after
the SSB of scale invariance whereas the other example has a potential such
that there is an internal $SO(n)$ global symmetry, which is spontaneously
broken to an $SO(n-1)$ symmetry, and consequently we have $n-1$ and $1$ Nambu-
Goldstone bosons corresponding to the SSB of the internal symmetry and scale
invariance, respectively.

The first example has a potential:

$\displaystyle
V(\phi)=\frac{1}{2}\sum_{I,J=1}^{n-1}\lambda_{IJ}\phi_{I}\phi_{J}\phi_{n}^{2},$
(16)

where the coupling constant is symmetric, $\lambda_{IJ}=\lambda_{JI}$. Then,
the first derivatives of the potential become

$\displaystyle\frac{\partial
V(\phi)}{\partial\phi_{I}}=\sum_{J=1}^{n-1}\lambda_{IJ}\phi_{J}\phi_{n}^{2},\qquad\frac{\partial
V(\phi)}{\partial\phi_{n}}=\sum_{I,J=1}^{n-1}\lambda_{IJ}\phi_{I}\phi_{J}\phi_{n}.$
(17)

Note that we can see that Eq. (2) is certainly satisfied in this case.

As a nontrivial vacuum which triggers the SSB of scale invariance, let us take
a stationary point

$\displaystyle(\phi_{I},\phi_{n})=(0,0,\cdots,0,v).$ (18)

At the vacuum (18), the second derivatives of the potential read

$\displaystyle\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{I}\partial\phi_{J}}\right|_{\phi=v}$
$\displaystyle=$
$\displaystyle\left.\lambda_{IJ}\phi_{n}^{2}\right|_{\phi=v}=\lambda_{IJ}v^{2},\qquad\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{I}\partial\phi_{n}}\right|_{\phi=v}=\left.2\sum_{J=1}^{n-1}\lambda_{IJ}\phi_{J}\phi_{n}\right|_{\phi=v}=0,$
$\displaystyle\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{n}^{2}}\right|_{\phi=v}$
$\displaystyle=$
$\displaystyle\left.\sum_{I,J=1}^{n-1}\lambda_{IJ}\phi_{I}\phi_{J}\right|_{\phi=v}=0.$
(19)

Using this result, we find that Eqs. (11) and (12) are valid so we can obtain
$n-1$ massive scalar fields and one massless dilaton after the SSB of scale
invariance.

Of course, this result can be checked directly by expanding the fields about
the vacuum (18) as

$\displaystyle\phi_{I}=\varphi_{I},\qquad\phi_{n}=v+\varphi_{n},$ (20)

by which the potential (16) can be rewritten as

$\displaystyle V(\phi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}v^{2}\sum_{I,J=1}^{n-1}\lambda_{IJ}\varphi_{I}\varphi_{J}+v\sum_{I,J=1}^{n-1}\lambda_{IJ}\varphi_{I}\varphi_{J}\varphi_{n}+\frac{1}{2}\sum_{I,J=1}^{n-1}\lambda_{IJ}\varphi_{I}\varphi_{J}\varphi_{n}^{2}.$
(21)

This expression means that we have $n-1$ massive scalar fields $\varphi_{I}$
with the mass matrix squared $M_{IJ}^{2}=v^{2}\lambda_{IJ}$ and a massless
dilaton $\varphi_{n}$.

Next let us consider the second example whose potential is of form:

$\displaystyle
V(\phi)=\frac{\lambda}{4!}\left(\sum_{i=1}^{n}\phi_{i}^{2}-\phi_{n+1}^{2}\right)^{2},$
(22)

where $n+1$ massless scalar fields are introduced, one of which has negative
sign in front of the $\phi_{n+1}^{2}$ term whereas the other $n$ fields have
positive sign and give us an $SO(n)$ internal symmetry.

The stationary condition is defined as in Eq. (3) and the cosmological
constant is vanishing owing to scale invariance as in Eq. (4). Using Eq. (22),
we can also check explicitly the validity of the identity (2). Furthermore,
Eq. (5) yields

$\displaystyle\sum_{i=1}^{n+1}\phi_{i}\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}=3\frac{\partial
V(\phi)}{\partial\phi_{j}}.$ (23)

At the stationary point of the potential, this equation takes the form

$\displaystyle\sum_{j=1}^{n+1}\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}v_{j}=0.$
(24)

As before, the existence of a nonzero solution $v_{i}\neq 0$ provides us with
Eq. (8). As a stationary point, let us take $v_{i}$ to

$\displaystyle
v_{i}\equiv(v_{1},v_{2},\cdots,v_{n},v_{n+1})=(0,0,\cdots,0,v,v).$ (25)

Then, on account of Eqs. (22) and (25), it turns out that the rank of the
matrix in Eq. (24) is given by

$\displaystyle{\rm
rank}\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}=1.$
(26)

This result tells us that we have only one massive particle and the remaining
$n$ particles become massless after SSB. In other words, the vacuum (25)
triggers the SSB of an $SO(n)$ symmetry to an $SO(n-1)$ symmetry in addition
to the scale symmetry, and consequently we have $n$ massless NB bosons, one of
which is the dilaton generated by the SSB of scale invariance and the other
$n-1$ NB bosons come from that of the internal symmetry.

The above result can be also investigated by expanding the fields around the
nontrivial vacuum (25) obeying the stationary condition as follows. Taking the
variation of the potential (22) with respect to $\phi_{i}$ with
$i=1,2,\cdots,n$ and $\phi_{n+1}$, respectively, yields

$\displaystyle\frac{\partial V(\phi)}{\partial\phi_{i}}$ $\displaystyle=$
$\displaystyle\frac{\lambda}{3!}\phi_{i}\left(\sum_{i=1}^{n}\phi_{i}^{2}-\phi_{n+1}^{2}\right),$
$\displaystyle\frac{\partial V(\phi)}{\partial\phi_{n+1}}$ $\displaystyle=$
$\displaystyle-\frac{\lambda}{3!}\phi_{n+1}\left(\sum_{i=1}^{n}\phi_{i}^{2}-\phi_{n+1}^{2}\right),$
(27)

from which we can see the existence of a flat direction defined as

$\displaystyle\sum_{i=1}^{n}\phi_{i}^{2}-\phi_{n+1}^{2}=0.$ (28)

Since the nontrivial vacuum configuration (25) is a solution to Eq. (28), we
can expand the fields around this configuration as follows:

$\displaystyle\phi_{I}=\varphi_{I},\qquad\phi_{n}=v+\varphi_{n},\qquad\phi_{n+1}=v+\varphi_{n+1},$
(29)

where $I=1,2,\cdots,n-1$. Then, the potential (22) can be cast to the form

$\displaystyle V(\phi)$ $\displaystyle=$ $\displaystyle\frac{\lambda
v^{2}}{6}(\varphi_{n}-\varphi_{n+1})^{2}+\frac{\lambda
v}{6}(\varphi_{n}-\varphi_{n+1})\left(\sum_{i=1}^{n}\varphi_{i}^{2}-\varphi_{n+1}^{2}\right)$
(30) $\displaystyle+$
$\displaystyle\frac{\lambda}{4!}\left(\sum_{i=1}^{n}\varphi_{i}^{2}-\varphi_{n+1}^{2}\right)^{2}.$

Introducing new fields defined as

$\displaystyle\chi_{n}=\frac{\varphi_{n}+\varphi_{n+1}}{\sqrt{2}},\qquad\chi_{n+1}=\frac{\varphi_{n}-\varphi_{n+1}}{\sqrt{2}},$
(31)

the potential (30) can be rewritten as

$\displaystyle V(\phi)=\frac{\lambda
v^{2}}{3}\chi_{n+1}^{2}+\frac{\sqrt{2}\lambda
v}{6}\chi_{n+1}\left(\sum_{I=1}^{n-1}\varphi_{I}^{2}+2\chi_{n}\chi_{n+1}\right)+\frac{\lambda}{4!}\left(\sum_{I=1}^{n-1}\varphi_{I}^{2}+2\chi_{n}\chi_{n+1}\right)^{2}.$
(32)

Note that at the same time the kinetic term can be expressed in terms of the
new fields in a canonical form:

$\displaystyle-\frac{1}{2}\sum_{i=1}^{n+1}(\partial_{\mu}\phi_{i})^{2}=-\frac{1}{2}\sum_{I=1}^{n-1}(\partial_{\mu}\varphi_{I})^{2}-\frac{1}{2}(\partial_{\mu}\chi_{n})^{2}-\frac{1}{2}(\partial_{\mu}\chi_{n+1})^{2}.$
(33)

As a result, we have one massive scalar field $\chi_{n+1}$ with mass squared
$M_{\chi_{n+1}}^{2}=\frac{2}{3}\lambda v^{2}$, a massless dilaton $\chi_{n}$
and $n-1$ massless NG bosons $\varphi_{I}$ corresponding to the SSB from an
$SO(n)$ to an $SO(n-1)$, as expected.

To close this section, it is worthwhile to point out the following facts
obtained above: In classically scale invariant theories, before the SSB of
scale invariance, all particles are massless. Once scale invariance is
spontaneously broken, the massless particles become massive except for the
dilaton which is a NB boson corresponding to the SSB of scale invariance.
Provided that there are internal symmetries in addition to scale invariance in
a classical action, we have additional NB bosons associated with the SSB of
the internal symmetries. In this case, all particles become massive but the NB
bosons and the magnitude of mass is proportional to the scale $v$ of the SSB.

##  3 Quantum anomaly

In this section, we work with quantum field theory (QFT) where we encounter a
quantum anomaly called scale anomaly. It is well known that the presence of
the ultra-violet divergences in QFT makes it necessary to introduce a
dimensionful renormalization scale $\mu$ in a theory, which could violate
classical scale invariance. Provided that the renormalization scale $\mu$ is
introduced, the dimension counting identity must be generalized to

$\displaystyle\left(\mu\frac{\partial}{\partial\mu}+\sum_{i=1}^{n}\phi_{i}\frac{\partial}{\partial\phi_{i}}\right)V(\phi)=4V(\phi),$
(34)

where the dimension counting identity picks up an anomalous piece
corresponding to a derivative with respect to the renormalization scale $\mu$.

Now, at the stationary point $\phi_{i}=v_{i}$ obeying Eq. (3)444Precisely
speaking, in QFT the stationary point should be understood as the vacuum
expectation value of the field operator, i.e., $\langle\phi_{i}\rangle=v_{i}$,
and a classical potential must be replaced by an effective potential., the
cosmological constant is not zero and given by

$\displaystyle V(v)=\frac{\mu}{4}\left.\frac{\partial
V(\phi)}{\partial\mu}\right|_{\phi=v}\neq 0.$ (35)

Taking the partial derivative of Eq. (34) with respect to the scalar fields
$\phi_{j}$ produces555Here, for simplicity, we assume that the renormalized
field $\phi_{i}$ is independent of $\mu$, i.e., the anomalous dimension is
zero, which holds at one loop level.

$\displaystyle\sum_{i=1}^{n}\phi_{i}\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}=3\frac{\partial
V(\phi)}{\partial\phi_{j}}-\mu\frac{\partial^{2}V(\phi)}{\partial\mu\partial\phi_{j}}.$
(36)

At the stationary point, this equation becomes

$\displaystyle\sum_{j=1}^{n}\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}v_{j}=-\left.\mu\frac{\partial^{2}V(\phi)}{\partial\mu\partial\phi_{i}}\right|_{\phi=v}.$
(37)

Since the RHS is not in general vanishing because of scale anomaly, for
nonvanishing $v_{i}$ we must have the equation:

$\displaystyle\det\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}\neq
0.$ (38)

Then, since the matrix
$\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}$
is a symmetric matrix and can be therefore diagonalized by an orthogonal
matrix, Eq. (38) implies that

$\displaystyle{\rm
rank}\left.\frac{\partial^{2}V(\phi)}{\partial\phi_{i}\partial\phi_{j}}\right|_{\phi=v}=n.$
(39)

Eq. (39) shows that all the massless scalar fields become massive after
symmetry breaking of scale invariance owing to the scale anomaly. To put it
differently, we have no more a massless dilaton under the presence of the
quantum anomaly at the quantum level although the theory was scale invariant
at the classical level.

In order to verify this fact explicitly, let us consider an example with a
single massless scalar field whose Lagrangian density is of form:

$\displaystyle{\cal{L}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-V(\phi)$ (40)
$\displaystyle=$
$\displaystyle-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-\frac{\lambda}{4!}\phi^{4},$

which is manifestly invariant under scale transformation.

By following Coleman and Weinberg [13] and using the modified minimal-
subtraction (or $\overline{MS}$ scheme), we calculate the effective potential
to one loop level:

$\displaystyle
V_{eff}(\phi)=\frac{\lambda}{4!}\phi^{4}+\frac{\lambda^{2}}{256\pi^{2}}\phi^{4}\left(\ln\frac{\lambda\phi^{2}}{2\mu^{2}}-\frac{3}{2}\right),$
(41)

where $\mu$ is the renormalization scale. It is easy to see that this
effective potential satisfies the generalized dimension counting identity (34)
in the sense that

$\displaystyle\left(\mu\frac{\partial}{\partial\mu}+\phi\frac{\partial}{\partial\phi}\right)V_{eff}(\phi)=4V_{eff}(\phi).$
(42)

Taking the derivatives of the effective potential over $\phi$ leads to

$\displaystyle\frac{\partial V_{eff}(\phi)}{\partial\phi}$ $\displaystyle=$
$\displaystyle\lambda\phi^{3}\left[\frac{1}{3!}+\frac{\lambda}{64\pi^{2}}\left(\ln\frac{\lambda\phi^{2}}{2\mu^{2}}-1\right)\right],$
$\displaystyle\frac{\partial^{2}V_{eff}(\phi)}{\partial\phi^{2}}$
$\displaystyle=$
$\displaystyle\frac{\lambda}{2}\phi^{2}\left[1+\frac{\lambda}{32\pi^{2}}\left(3\ln\frac{\lambda\phi^{2}}{2\mu^{2}}-1\right)\right].$
(43)

From these equations, we can find $\langle\phi\rangle=v$ which minimizes the
effective potential, and the dilaton mass squared:

$\displaystyle v$ $\displaystyle=$
$\displaystyle\mu\sqrt{\frac{2}{\lambda}}\exp\left({\frac{3\lambda-32\pi^{2}}{6\lambda}}\right),$
$\displaystyle m_{D}^{2}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial^{2}V_{eff}(\phi)}{\partial\phi^{2}}\right|_{\phi=v}=\frac{\lambda^{2}v^{2}}{32\pi^{2}}.$
(44)

Note that since $v\propto\mu$, the dilaton mass is proportional to the
renormalization scale $\mu$, which describes the magnitude of a violation of
scale invariance.

##  4 Quantum scale invariance

It is a familiar fact that some fields have built-in mechanisms to keep them
massless; vector fields have gauge invariance and fermions have chiral
invariance. Perhaps, it is scale invariance that would make scalar fields be
massless, but the existence of the ultra-violet divergences in QFT forces us
to introduce a dimensionful renormalization scale $\mu$ in a theory, which
violates classical scale invariance inevitably at the quantum level. Without
scale invariance, scalar fields would have no immunity to protect them from
acquiring a nonzero mass.

However, if there were no scale anomaly at the quantum level, we would also
have the dimension counting identity (2). Then, by following our argument in
Section 2, we find that the dilaton remains massless even at the quantum
regime as well as the classical one. In fact, Englert et al. [6] have already
proposed an interesting approach such that scale invariance can be maintained
even at the quantum level if we introduce a scalar field, which was called
“dilaton”, instead of the renormalization scale $\mu$. Their reasoning is so
simple that let us follow it briefly: By using the dilaton, the Einstein-
Hilbert action in general relativity, which has no scale anomaly, can be
transformed to a locally scale invariant scalar-tensor gravity [14, 15] 666A
relation between the locally scale invariant scalar-tensor gravity and Weyl
transverse gravity was examined in Ref. [16] and the Higgs mechanism in scale
invariant gravitational theories has been investigated in Ref. [17]., which
could have a Weyl anomaly corresponding to local scale invariance.777Here the
difference between global and local scale invariances is not important since
any globally scale invariant theories can be promoted to local ones by
coupling gravity in an appropriate manner. With the order reversed, if we
start with the locally scale invariant scalar-tensor gravity and take the
unitary gauge, i.e., a constant dilaton, we reach the Einstein-Hilbert action.
This quantum equivalence between general relativity and the locally scale
invariant scalar-tensor gravity then suggests that the anomaly associated with
local scale invariance, which exists in the scale invariant scalar-tensor
gravity, might be fake since there is no such an anomaly in general
relativity.

It is true that the approach advocated by Englert et al. [6] is indeed free
from the quantum anomaly and we can explicitly check that the dilaton is
exactly massless even at the quantum level [10]. However, in this approach,
which we will call the scale invariant regularization method, there are
several problems we have to pay the price of success. First, such scale
invariant theories are non-renormalizable as in general relativity and make
sense only as effective field theories [11]. Second, there is a fine-tuning
problem for keeping a flat direction stable against radiative corrections,
which can be understood as the incarnation of the cosmological constant
problem [12]. Finally, it is not clear how to treat with the unbroken phase
($\langle\phi\rangle=0$ where $\phi$ is a dilaton) since the formalism makes
sense only in the broken phase ($\langle\phi\rangle\neq 0$) [9]. In this
section, we wish to point out that there is another phenomenological problem
when we attempt to apply the scale invariant regularization method for
cosmology in addition to the above three problems.

It is natural that the massless dilaton $\phi$ couples to the Einstein-Hilbert
term with the non-minimal coupling term888We follow the conventions and
notation of the MTW textbook [18].

$\displaystyle\frac{1}{\sqrt{-g}}{\cal{L}}=\frac{1}{2}\xi\phi^{2}R-\frac{1}{2}\epsilon
g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi+L_{m},$ (45)

where $\epsilon=\pm 1=\frac{|\omega|}{\omega}$ with $\omega$ being the Brans-
Dicke parameter which will appear shortly. Here $\epsilon=1$ corresponds to a
normal field with a positive energy while $\epsilon=-1$ does to a ghost field
with a negative energy. Moreover, $\xi$ is a constant and in particular the
case of $\xi=\frac{1}{6}$ and $\epsilon=-1$ produces a locally scale invariant
scalar-tensor gravity, and $L_{m}$ denotes the Lagrangian for a generic matter
field without mass terms. Finally, it is worthwhile to stress that the
Langrangian density (45) is invariant under scale transformation.

By defining a new scalar field $\varphi$, which is sometimes called the Brans-
Dicke scalar field, and the Brans-Dicke parameter $\omega$ as

$\displaystyle\varphi=\frac{1}{2}\xi\phi^{2},\qquad\omega=\frac{\epsilon}{4\xi},$
(46)

we can rewrite (45) as a well-known Brans-Dicke Lagrangian:

$\displaystyle\frac{1}{\sqrt{-g}}{\cal{L}}=\varphi
R-\omega\frac{1}{\varphi}g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi+L_{m}.$
(47)

The field equation for the Brans-Dicke scalar field then turns out to be

$\displaystyle\Box\varphi=\frac{1}{4\omega+6}T,$ (48)

where $T$ denotes the trace part of the energy-momentum tensor $T_{\mu\nu}$
for matter field defined as
$T_{\mu\nu}\equiv-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}L_{m})}{\delta
g^{\mu\nu}}$.

In the manifestly scale invariant regularization method, scale invariance is
expected to be broken by Planck-scale physics [10], so the dilaton will
acquire a large vacuum expectation value (VEV), $\langle\phi\rangle\sim
M_{Pl}$ where $M_{Pl}$ is the reduced Planck mass defined in terms of the
Newton constant $G$ as $M_{Pl}\equiv\frac{1}{\sqrt{8\pi G}}$. Indeed, this
expectation can be realized by the following argument: First, let us expand
the dilaton $\phi$ about its VEV as

$\displaystyle\phi=\langle\phi\rangle+\tilde{\phi},$ (49)

and substitute this expansion into the Lagrangian (45), from which we can
obtain

$\displaystyle\frac{1}{\sqrt{-g}}{\cal{L}}=\frac{1}{2}\xi\langle\phi\rangle^{2}R+\xi\left(\langle\phi\rangle\tilde{\phi}+\frac{1}{2}\tilde{\phi}^{2}\right)R-\frac{1}{2}\epsilon
g^{\mu\nu}\partial_{\mu}\tilde{\phi}\partial_{\nu}\tilde{\phi}+L_{m},$ (50)

Next, to make the first term on the RHS of Eq. (50) coincide with the
Einstein-Hilbert term, we have to choose
$\xi\langle\phi\rangle^{2}=M_{Pl}^{2}$. With $\xi={\cal{O}}(1)$, we have
$\langle\phi\rangle\sim M_{Pl}$. This completes our confirmation that scale
invariance is broken around the Planck scale.999Incidentally, as will be seen
in Eq. (51), $\xi$ is a tiny quantity. If we take $\xi=10^{-6}$ for instance,
we have $\langle\phi\rangle=10^{3}M_{Pl}$, which is beyond the Planck mass so
it is not physically interesting.

Eq. (48) implies that the Brans-Dicke scalar field mediates a long-range force
between massive objects as in the Newtonian force and its coupling strength is
given by $\frac{1}{4\omega+6}$, which is roughly of the same order of
magnitude as that of the Newtonian force since
$\omega=\frac{\epsilon}{4\xi}={\cal{O}}(1)$, that is, $\omega$ is of the order
of unity. Note that as in the Brans-Dicke theory, most of models related to
the dilaton have couplings with matters which are of gravitational strength.

On the other hand, the currently best time delay measurements by the Cassini
probe in the solar system gives us the bound [19]

$\displaystyle\omega>4.3\times 10^{3}\ (\mbox{or}\quad\xi<5.8\times
10^{-5}),\qquad\epsilon=1.$ (51)

Unfortunately, the magnitude of the Brans-Dicke parameter $\omega$ is too
large to be no room that it is close to anywhere around the unity. This
problem is known as a long-standing curse in the Brans-Dicke theory and its
generalizations [2]. To solve this problem, some people have advocated that
the dilaton might have small mass in such a way that it mediates a finite-
range force [4]. If its Compton wave-length is smaller than the size of the
solar system, we can neglect the finite-range force from the massive dilaton,
thereby making the constraint (51) be irrelevant phenomenologically.

However, with quantum scale invariance, the dilaton is exactly massless at the
both classical and quantum levels, so it would be not impossible but very
difficult to construct a realistic cosmological model which is consistent with
the Cassini bound (51).

##  5 Conclusions

In this article, on the basis of the dimension counting identity and its
generalization, we have given a concise proof that the dilaton is exactly
massless at the classical level if scale invariance is broken spontaneously
whereas it becomes massive at the quantum level if scale invariance is
explicitly broken by quantum anomaly. Moreover, we have shown that when
quantum scale invariance is maintained through the scale invariant
regularization method, the dilaton remains massless at the quantum level.
Compared to the proof in Ref. [5], our proof is simpler and has an advantage
of being able to apply for non-renormalizable quantum field theory.

It is of interest to notice that before the spontaneous symmetry breakdown of
scale invariance, all the particles are massless because of scale invariance,
and after the SSB they become massive except for the dilaton unless there are
internal symmetries such as $SO(n)$. Here all masses arise from a nonzero
vacuum expectation value of the dilaton field.

Furthermore, we have pointed out that in addition to three problems of the
manifestly scale invariant regularization scheme where the renormalization
scale is replaced with a dilaton field, there is the fourth problem which
emerges when we couple the massless dilaton to the Einstein-Hilbert term in
the non-minimal way; a difficulty of constructing a realistic cosmological
model by using the massless dilaton. It seems that a crucial point is to
develop a new mechanism for nullifying the scalar force stemming from the
dilaton. We wish to find such a mechansim in future.101010While preparing for
this paper, an approach for screening dilaton interaction was proposed in
[20].

Acknowledgements

This work is partly supported by the JSPS Kakenhi Grant No. 21K03539.

## References

  * [1] J. Polchinski, ”String Theory”, Cambridge University Press, 1998.
  * [2] Y. Fujii and K. Maeda, ”The Scalar-Tensor Theory of Gravitation”, Cambridge University Press, 2003.
  * [3] C. Brans and R. H Dicke, “Mach’s Principle and a Relativistic Theory of Gravitation”, Phys. Rev. 124 (1961) 925.
  * [4] Y. Fujii, “Dilaton and Possible Non-Newtonian Gravity”, Nature Phys. Sci. 234 (1971) 5.
  * [5] D. Nogradi and B. Ozsvath, ”Dilaton in Scalar QFT: A No-go Theorem”, arXiv: 2109.09822 [hep-th].
  * [6] F. Englert, C. Truffin and R. Gastmans, ”Conformal Invariance in Quantum Gravity”, Nucl. Phys. B 117 (1976) 407.
  * [7] M. Shaposhnikov and D. Zenhausern, “Quantum Scale Invariance, Cosmological Constant and Hierarchy Problem”, Phys. Lett. B 671 (2009) 162.
  * [8] A. Codello, G. D’Odorico, C. Pagani and R. Percacci, “The Renormalization Group and Weyl-invariance”, Class. Quant. Grav. 30 (2013) 115015.
  * [9] C. Tamarit, “Running Couplings with A Vanishing Scale Anomaly”, JHEP12 (2013) 098.
  * [10] D. M. Ghilencea, “Manifestly Scale-invariant Regularization and Quantum Effective Operators”, Phys. Rev. D 93 (2016) 105006.
  * [11] M. Shaposhnikov and F. Tkachov, “Quantum Scale-invariant Models as Effective Field Theories”, arXiv:0905.4857 [hep-th].
  * [12] T. Kugo, “Necessity and Insufficiency of Scale Invariance for Solving Cosmological Constant Problem”, PoS CORFU2019 (2020) 071.
  * [13] S. R. Coleman and E. J. Weinberg, ”Radiative Corrections as the Origin of Spontaneous Symmetry Breaking”, Phys. Rev. D 7 (1973) 1888.
  * [14] P. A. M. Dirac, ”Long Range Forces and Broken Symmetries”. Proc. R. Soc. A 333 (1973) 403.
  * [15] S. Deser, ”Scale Invariance and Gravitational Coupling”, Ann. Phys. Rev. (N. Y) 59 (1970) 248.
  * [16] I. Oda, “Classical Weyl Transverse Gravity”, Eur. Phys. J. C 77 (2017) 284.
  * [17] I. Oda, “Higgs Mechanism in Scale-Invariant Gravity”, Adv. Stud. Theor. Phys. 8 (2014) 215\.
  * [18] C. W. Misner, K. S. Thorne and J. A. Wheeler, “Gravitation”, W H Freeman and Co (Sd), 1973.
  * [19] B. Bertotti, L. Iess and P. Tortora, “A Test of General Relativity Using Radio Links with the Cassini Spacecraft”, Nature 425 (2003) 374.
  * [20] C. P. Burgess and F. Quevedo, ”Axion Homeopathy: Screening Dilaton Interactions”, arXiv:2110.10352 [hep-th].


# Unified Differentiable Learning of the Electric Enthalpy and Dielectric
Properties with Exact Physical Constraints

Stefano Falletta<EMAIL_ADDRESS>John A. Paulson School of
Engineering and Applied Sciences, Harvard University, Cambridge MA, USA
Andrea Cepellotti  John A. Paulson School of Engineering and Applied Sciences,
Harvard University, Cambridge MA, USA  Chuin Wei Tan  John A. Paulson School
of Engineering and Applied Sciences, Harvard University, Cambridge MA, USA
Anders Johansson  John A. Paulson School of Engineering and Applied Sciences,
Harvard University, Cambridge MA, USA  Albert Musaelian  John A. Paulson
School of Engineering and Applied Sciences, Harvard University, Cambridge MA,
USA  Cameron J. Owen  Department of Chemistry and Chemical Biology, Harvard
University, Cambridge MA, USA  Boris Kozinsky  John A. Paulson School of
Engineering and Applied Sciences, Harvard University, Cambridge MA, USA
Robert Bosch LLC Research and Technology Center, Watertown, MA, USA

###### Abstract

The modern theory of polarization and electric enthalpy functionals have paved
the way to first-principles studies of polarization in crystalline, disordered
and liquid materials. However, calculating the polarization in large-scale ab-
initio molecular dynamics remains prohibitive due to the high computational
cost involved. In this work, we introduce a machine-learning formulation that
simultaneously learns to predict electric enthalpy, atomic forces,
polarization, Born charges, and polarizability for a given atomic
configuration. These quantities are derived through differentiation of the
electric enthalpy with respect to atomic positions and electric field within a
local equivariant representation of the atomic environment. Our new method
exactly constrains the polarization to be a conservative vector field, obeys
the acoustic sum rule for Born charges, and ensures the conservation of
electric enthalpy in machine-learning molecular dynamics. Importantly, the
formulation is general and does not necessarily require equivariance. The
locality of the model enables the efficient calculation of vibrational and
dielectric properties of materials over large length-scale and long time-scale
molecular dynamics at quantum-mechanical accuracy. We apply our approach to
determine the infrared spectrum and frequency-dependent dielectric constant of
$\alpha$-SiO2 using molecular dynamics, finding excellent agreement with
results obtained from density functional perturbation theory. Additionally, we
show that our model accurately captures the electric-field contributions to
the electric enthalpy, forces, and polarization in reference to density
functional theory calculations. This demonstrates the capability of our
formulation to perform large-scale molecular dynamics under electric fields,
thus paving the way to efficient and accurate studies of dielectric response
of complex materials.

††preprint: APS/123-QED

##  I Introduction

Polarization is a crucial physical quantity for understanding the properties
of dielectrics and ferroelectrics [1, 2]. The modern theory of polarization
[3, 1, 4, 5, 6, 7, 8, 2, 9] and the development of electric enthalpy
functionals [10, 11] have enabled first-principles studies of polarization
with density functional theory (DFT). Predicting the polarization over a
molecular dynamics (MD) trajectory allows for the calculation of vibrational
and dielectric properties of crystalline, disordered and liquid materials [12,
13, 14, 15, 16, 17, 18, 19], such as infrared spectra and dielectric
constants. However, performing large length-scale and long time-scale MD
simulations using first-principles methods remains challenging due to the
associated high computational cost.

Machine-learning (ML) force fields bridge the gap between accuracy and
efficiency of MD simulations. The accuracy of ML force fields has been
significantly enhanced in recent years by incorporating exact $E(3)$ symmetry
equivariance (translations, rotations and inversion), starting with the NequIP
model and subsequent approaches [20, 21, 22, 23, 24, 25]. Recently, ML
approaches have been used to predict properties beyond energy and forces [26],
including polarization [27, 28, 29, 30, 31, 32, 33, 34, 35, 36], Born charges
[37, 38, 35], polarizability [39, 40, 41, 42, 35], dielectric constants [37,
43, 29, 31, 35, 36], infrared spectra [37, 33, 34, 44, 35, 36], and Raman
spectra [40, 41, 44, 42, 45], with applications to molecules [39, 30, 41, 34,
28, 27, 42], liquid water [37, 40, 29, 33, 44, 35, 36], and solids [31, 38,
41, 44]. These methods are based on a variety of techniques, including
Gaussian process regression [39, 28, 31, 42, 45], kernel ridge regression [30,
34, 41], random forest regression [43], other kernel-based regression models
[27, 30, 37], deep neural networks [40, 29, 38, 35], and message-passing
equivariant neural networks [33, 44]. To study polarization properties during
a MD, current ML force fields are utilized in conjunction with specialized ML
models. Specifically, the ML force field is used to predict energies and
forces for the time evolution, while separate dedicated models are used to
predict polarization properties for each structure. However, this does not
guarantee the enforcement of physical constraints on electric enthalpy,
polarization, and Born charges. This includes ensuring that polarization is a
conservative vector field, Born charges adhere to the acoustic sum rule, and
the electric enthalpy being conserved during machine-learning molecular
dynamics (MLMD). These constraints are relevant as they relate to the quantum
nature of polarization, to the translational invariance of polarization and
charge neutrality of the system, and to electric enthalpy conservation.

In this work, we introduce a self-consistent neural network formulation for
learning polarization, Born charges, and polarizability in addition to the
electric enthalpy, forces, and stress. The polarization is obtained by
differentiating the electric enthalpy with respect to the electric field, the
Born charges by differentiating the polarization with respect to the atomic
positions, the polarizability by differentiating the polarization with respect
to the electric field. These derivative relations ensure the enforcement of
physical constraints, namely the polarization being a conservative vector
field, the acoustic sum rule for Born charges, and the conservation of the
electric enthalpy in machine-learning molecular dynamics. Our approach is
based on a strictly local equivariant description of the atomic environment,
thereby overcoming the limitations of message-passing neural networks. The
locality of the model enables efficient calculations of electric enthalpy,
forces, stress, polarization, Born charges, and polarizability over large
length-scale and long time-scale molecular dynamics simulations. We apply the
method to $\alpha$-SiO2, and obtain the infrared spectrum and frequency-
dependent dielectric constant in excellent agreement with results computed
from density functional perturbation theory. Furthermore, we demonstrate that
our formulation can be used to perform molecular dynamics under finite
electric fields, since it captures the electric-field contributions to the
electric enthalpy, forces, and polarization. This paves the way to large-scale
studies of vibrational and dielectric properties of crystalline, disordered
and liquid materials, and to large-scale molecular dynamics of materials under
electric fields at the first-principles level of accuracy.

##  II Methodology

Figure 1: a, Schematics outlining our neural network formulation that
simultaneously learns to predict electric enthalpy $E$, forces
$\mathbf{F}_{i}$, polarization $\mathbf{P}$, Born charges
$\mathbf{Z}^{*}_{i}$, and polarizability $\bm{\alpha}$ given in input the
atomic positions $\mathbf{r}_{i}$ and the electric field $\bm{\mathcal{E}}$.
Here, $i$ denotes the atom index. b, Workflow for generating DFT data through
calculations with finite electric fields. The DFT labels used for training and
validation are calculated in the limit of zero electric field, and are circled
for ease of reference.

For a system subject to a uniform electric field, one can define the electric
enthalpy functional as follows [10, 46, 11]

$E=E^{0}-\bm{\mathcal{E}}\cdot\mathbf{P},$ (1)

where $E^{0}$ is the energy functional in the absence of an electric field,
$\bm{\mathcal{E}}$ the electric field, and $\mathbf{P}$ the polarization as
calculated within the modern theory of polarization [4, 5, 2]. The
polarization is related to the electric enthalpy through the following
differential relation,

$P_{\mu}=-\frac{dE}{d\mathcal{E}_{\mu}},$ (2)

where $\mu$ denotes the Cartesian component. In the modern theory of
polarization, $\mathbf{P}$ is a multi-valued quantity defined modulo the
quantum of polarization, namely $\Delta\mathbf{P}=e\mathbf{R}$, where
$\mathbf{R}$ is a lattice vector. In the limit of a vanishing electric field,
the derivative of the polarization with respect to the atomic displacement
defines the Born charges [47]:

$Z^{*}_{i\mu\nu}=\left.\frac{1}{e}\frac{dP_{\mu}}{dr_{i\nu}}\right|_{\mathcal{E}_{\mu}=0},$
(3)

where $r_{i\nu}$ is the Cartesian coordinate $\nu$ of the atom $i$. Born
charges obey the acoustic sum rule, namely $\sum_{i}Z^{*}_{i\mu\nu}=0$ [48].
This reflects the translation invariance of polarization and the charge
neutrality of the system. The derivative of the polarization with respect to
the electric field defines the polarizability, namely

$\alpha_{\mu\nu}=\frac{dP_{\mu}}{d\mathcal{E}_{\nu}}.$ (4)

We use the differential relations in Eqs. (2), (3) and (4) to formulate the
learning approach for polarization, Born charges and polarizability in
addition to electric enthalpy, forces, and stress. This is based on performing
derivatives of the electric enthalpy, a scalar quantity, with respect to
atomic positions and electric field. The framework is general and can be
implemented in a variety of ML methods, including both invariant and
equivariant neural networks [20, 21] and kernel-based methods [49]. In this
work, we consider the equivariant neural network Allegro [21, 22], which
offers the accuracy advantage of an equivariant neural network that does not
involve message passing. This overcomes the efficiency and scalability
limitations of message-passing neural networks while maintaining the same
level of accuracy [21, 22]. In addition, Allegro minimizes the amount of data
in input for training compared to other existing methods [50, 51, 52].

Our unified framework for force field and polarization properties is
implemented as follows. First, we include the electric field as an input of
the network along with the atomic positions. In particular, the spherical
harmonics describing the electric field and interatomic distances are
concatenated and treated on the same footing within the neural network
architecture [21]. Thanks to the modularity of the Allegro code, this
incorporation is achieved without significant modifications of its core
architecture. Then, we determine the polarization by differentiating the
output of the model, i.e. the electric enthalpy, with respect to the electric
field. This procedure is analogous to that used to derive atomic forces as
derivatives of the electric enthalpy with respect to the atomic positions (cf.
Fig. 1a). Then, following Eqs. (3) and (4), we derive Born charges and
polarizability by differentiating the polarization with respect to atomic
positions and electric field, respectively. Since Born charges are per-atom
quantities, they provide a large amount of data for learning the polarization,
thus increasing the data efficiency. Polarization, Born charges and
polarizability can be learned along with the electric enthalpy, atomic forces,
and stress by adding the following extra contribution to the conventional
force field loss function:

$\displaystyle\Delta\mathcal{L}$
$\displaystyle=\frac{\lambda_{P}}{3N}\sum_{\mu=1}^{3}\left|\left(-\frac{d\hat{E}}{d\mathcal{E}_{\mu}}-P_{\mu}\right)\text{mod}\
\Delta P_{\mu}\right|$
$\displaystyle\qquad+\frac{\lambda_{Z}}{9N}\sum_{i=1}^{N}\sum_{\mu=1}^{3}\sum_{\nu=1}^{3}\left|-\frac{d^{2}\hat{E}}{d\mathcal{E}_{\mu}dr_{i\nu}}-Z^{*}_{i\mu\nu}\right|$
$\displaystyle\qquad+\frac{\lambda_{\alpha}}{9N}\sum_{\mu=1}^{3}\sum_{\nu=1}^{3}\left|-\frac{d^{2}\hat{E}}{d\mathcal{E}_{\mu}d\mathcal{E}_{\nu}}-\alpha_{\mu\nu}\right|,$
(5)

where $\lambda_{P}$, $\lambda_{Z}$, and $\lambda_{\alpha}$ are the loss
weights, $N$ is the number of atoms, and $\hat{E}$ is the predicted electric
enthalpy. In Eq. (5), the loss contributions of extensive quantities, namely
polarization and polarizability, are normalized by the number of atoms.
Moreover, the loss contribution related to the polarization is constructed to
give the minimal image difference between the predicted and the DFT values of
polarization, taking into account the multi-value nature of the polarization.

Large-scale MLMD for determining electric enthalpy, forces and polarization
properties can be performed with LAMMPS [53], through a dedicated interface
that we developed for our ML model. Our interface generally allows to perform
MLMD in the presence of a time-dependent and space-dependent electric field.
This requires accounting for the corresponding variation in electric enthalpy,
forces and polarization due to the electric field. For simplicity, we here
consider the case of a constant electric field. In this case, the electric
enthalpy is obtained as in Eq. (1). The forces are calculated as

$\mathbf{F}_{i}=\mathbf{F}_{i}^{0}+e\bm{\mathcal{E}}\cdot\mathbf{Z}_{i},$ (6)

where $\mathbf{F}_{i}^{0}$ is the force on atom $i$ in the absence of the
electric field. The polarization is given by

$\mathbf{P}=\mathbf{P}^{0}+\bm{\mathcal{E}}\cdot\bm{\alpha},$ (7)

where $\mathbf{P}^{0}$ is the polarization in the absence of the electric
field. In Eqs. (6) and (7), $\mathbf{F}_{i}$ and $\mathbf{P}$ obtained in the
presence of an electric field are calculated using quantities determined in
the limit of zero electric field, namely $\mathbf{F}^{0}$, $\mathbf{P}^{0}$,
$\mathbf{Z}_{i}$ and $\bm{\alpha}$. This stems from the linearity of the
electric enthalpy functional with respect to small electric fields in the
modern theory of polarization. Moreover, learning the Born charges is
essential for the correct evaluation of forces under an electric field.
Indeed, the change in forces due to small electric fields can become
comparable to the accuracy of the ML model for forces, and thus may not be
captured by training on only energy and forces in the presence of an electric
field.

Our method offers several physical and computational advantages. First, it is
physically exact as it predicts electric enthalpy and related derivatives
within a single unified model, in contrast to previous approaches that used
multiple independent models. In this regard, the polarization is guaranteed to
be a conservative vector field, since it is calculated as the gradient of the
electric enthalpy with respect to the electric field [cf. Eq. (2)]. This
implies that the integral of polarization over a cyclic adiabatic evolution of
the Hamiltonian is either zero or the quantum of polarization [2]. Moreover,
this ensures that differences of polarization between two states do not depend
on the adiabatic path connecting the states. Second, the acoustic sum rule for
Born charges is satisfied, as these are calculated as derivatives of
polarization with respect to atomic positions [cf. Eq. (3)]. This reflects the
translation invariance of polarization and the exact charge neutrality of the
system. Third, the electric enthalpy obtained in MLMD is conserved, as atomic
forces in the presence of a constant electric field are obtained as a gradient
of the electric enthalpy. Furthermore, our method can be applied to disordered
systems, as the polarization is calculated as derivative of the electric
enthalpy. Specifically, evaluating polarization via integration of Eq. (3)
solely captures polarization differences and requires a reference
configuration, which is not definable for disordered systems. Lastly, our
scheme can be easily generalized to determine additional quantities, such as
piezoelectric tensors ($e_{ijk}=-d\sigma_{ij}/d\mathcal{E}_{k}$, where
$\sigma_{ij}$ is the stress tensor), and to other forms of perturbations, such
as magnetic fields.

Our method also presents several practical advantages. First, it is based on
an equivariant local description of the atomic environment, eliminating the
need for message passing. This is crucial for performing large-scale MPI-
parallelized MD simulations, as the receptive field of message-passing neural
networks can grow excessively, constraining parallel computation and the scale
of MD simulations [21]. Second, the model predicts a scalar, i.e. the electric
enthalpy, and derives therefrom all quantities through automatic
differentiation. This ensures the involvement of a minimal number of tensorial
paths and favors the efficiency of the model, as opposed to predicting
directly vectorial or tensorial properties, as done in earlier approaches that
do not satisfy the exact physical constraints discussed above. Also, it is
possible to even use invariant scalar representations in a Gaussian process or
scalar neural network regression frameworks. Third, the equivariance of our
model reduces the amount of input data required for training [50, 51, 52],
conferring robustness to the model [20, 21, 22]. This is particularly valuable
when learning polarization and Born charges, which require computationally
expensive DFT calculations in the presence of electric fields. Furthermore,
MLMD can be conducted with or without electric fields using the model trained
under zero electric field conditions, eliminating the necessity for training
separate models at different electric field magnitudes. This offers advantages
in terms of computational cost and memory. Lastly, the learning of
polarization properties within the force field does not necessitate
equivariance and can be implemented in other existing approaches, such as
those based on Gaussian process regression [49].

##  III Results

We use our unified model for electric enthalpy, forces, and polarization
properties to predict the vibrational and dielectric properties of
$\alpha$-SiO2. First, we generate DFT data for training. We collect
uncorrelated frames thorugh classical MD simulations using the Vashishta
potential [54] in the NVT ensemble. These simulations are performed at both
300 K and 600 K, each for 100 ps. We take 100 uncorrelated snapshots from each
MD at intervals of 1 ps, yielding a total of 200 frames. Next, for each frame,
we perform DFT calculations to determine energy, forces, and polarization in
the absence of electric field. To calculate Born charges and polarizability,
we perform DFT calculations in the presence of small uniform electric fields,
and use finite-differences expressions involving forces [cf. Eq. (8)] and
polarization [cf. Eq. (4)], respectively. This procedure is illustrated in
Fig. 1b. Then, we train our ML models and exhaustively test them over various
hyperparameters, including cutoff radius, maximum order of spherical
harmonics, loss coefficients, number of tensor features, and network
architecture. We find that generally a maximum order of spherical harmonics of
3 or greater improves the model, especially for simpler network architectures.
Further computational details are provided in the Supplementary Information
(SI).

Figure 2: a, Autocorrelation function of polarization as a function of time
$t$ and b, infrared spectrum of $\alpha$-SiO2 at 300 K as a function of
frequency $\omega$ calculated through MLMD employing our unified model for
electric enthalpy, forces, and polarization properties. In b, the infrared
spectrum obtained with DFPT is reported for comparison. The infrared spectra
are smeared using a Gaussian broadening of 20 cm-1. The amplitudes of the main
infrared peak obtained in the two cases are matched as a guide to the eye.

We then construct a 24696-atom supercell and perform MLMD for 200 ps after 40
ps of equilibration in the absence of an electric field within the NVT
ensemble. This large size scale ensures excellent convergence of results and
gives a demonstration of the scalability of our method. By determining the
autocorrelation function of polarization during the MLMD (cf. Fig. 2a), we can
then calculate the infrared spectrum of $\alpha$-SiO2 [cf. Eq. (9)]. As
illustrated in Fig. 2b, we find the infrared spectrum peaks at frequencies of
1041, 420, and 765 cm-1. For comparison, we obtain the infrared spectrum of
$\alpha$-SiO2 through density functional perturbation theory (DFPT)
calculations [55, 56, 57, 58, 59], and find excellent agreement with the MLMD
result, as shown in Fig. 2b. The frequencies of the main infrared peaks
obtained in the two cases are given in Table 1. Next, the polarizability
during the MLMD is used to calculate the high-frequency dielectric constant
[cf. Eq. (10)]. Then, the static dielectric constant can be determined through
the fluctuation-dissipation theorem, which relies on the polarization values
obtained during the MLMD [cf. Eq. (11)]. By employing the autocorrelation
function of polarization from the MLMD [cf. Eq. (12)], we calculate the
frequency-dependent dielectric constant, which we illustrate in Fig. 3. For
comparison, we determine the frequency-dependent dielectric constant using the
DFPT approach in Refs. [60, 61], and find a very good agreement with the MLMD
results. This agreement in infrared spectra and frequency-dependent dielectric
constant is expected, due to the harmonicity of $\alpha$-SiO2.

Figure 3: a, Real part and b, imaginary part of the frequency-dependent dielectric constant of $\alpha$-SiO2 at 300 K as a function of frequency $\omega$ calculated through MLMD employing our unified model for electric enthalpy, forces, and polarization properties. Here, $\varepsilon=(\varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33})/3$. The results obtained with DFPT are reported for comparison. A Gaussian broadening of 20 cm-1 is used. Table 1: Main frequencies (in cm-1) of the infrared activity in $\alpha$-SiO2 calculated through MLMD and DFPT. | MLMD | DFPT
---|---|---
$\omega_{1}$ | 1041 | 1044
$\omega_{2}$ | 420 | 423
$\omega_{3}$ | 765 | 763

Furthermore, we verify that our formulation can be used to perform MLMD under
electric fields. Specifically, we study the extent to which our ML model
captures the electric-field contributions to electric enthalpy, forces and
polarization, with reference to DFT calculations in $\alpha$-SiO2. These
contributions can be directly quantified by determining the change in
dielectric constant upon structural relaxation due to the presence of an
electric field [cf. Eq. (13)]. Indeed, accurate electric-field contributions
on forces and polarization allow to correctly determine the structural
displacements due to the electric field and the associated change in
polarization. This requires accurate Born charges [cf. Eq. (6)] and
polarizability [cf. Eq. (7)], and yields accurate electric enthalpy [cf. Eq.
(1)]. To this aim, we first perform DFT calculations under a finite electric
field [11] applied along the $z$ direction starting from a pristine bulk
structure. As illustrated in Fig. 4a, we find high-frequency and static
dielectric constants of $\varepsilon_{zz}^{\infty}=2.41$ and
$\varepsilon_{zz}^{0}=4.76$, respectively. Next, starting from the same
initial structure, we perform a MLMD simulation at fixed volume using our ML
model under the same electric field, with temperature ramping down from 10 K
to 0.01 K. As illustrated in Fig. 4b, our model achieves the same values of
dielectric constants as obtained with DFT. This validates the accuracy of the
electric-field contributions to electric enthalpy, forces, and polarization,
and confirms the suitability of our formulation for performing MLMD under
finite electric fields.

Figure 4: Dielectric constant $\varepsilon_{zz}$ of $\alpha$-SiO2 as obtained
through DFT structural relaxations (in red) or through MLMD at fixed volume
(in blue) as a function of the relaxation steps in the presence of an electric
field along $\mathbf{z}$. The temperature of the MLMD is initialized to 10 K
and ramps down to 0.01 K. In both cases, the electric field is set to
$10^{-2}$ Ry a.u., and the initial structure is obtained from a DFT relaxation
in the absence of the electric field. The horizontal turquoise line represents
the average dielectric constant in the last steps of the MLMD. The high-
frequency dielectric constant $\varepsilon_{zz}^{\infty}$ and the static
dielectric constant $\varepsilon_{zz}^{0}$ obtained with DFT and MLMD
essentially coincide and are reported.

##  IV Conclusions

In conclusion, we introduced a unified equivariant neural network framework to
learn and predict polarization, Born charges, and polarizability in addition
to electric enthalpy, forces, and stress. This is achieved through
differentiation of the electric enthalpy with respect to the input variables
of the network, namely the atomic positions and the electric field. Our
approach enforces exact physical constraints, namely the polarization being a
conservative vector field, the acoustic sum rule for Born charges, and the
conservation of the electric enthalpy in machine-learning molecular dynamics.
Our method is based on an equivariant local description of the atomic
environment, which allows to efficiently determine polarization properties in
large size-scale and time-scale molecular dynamics simulations at first-
principles accuracy. We applied the model to study the vibrational and
dielectric properties of $\alpha$-SiO2, and found that the infrared spectrum
and the frequency-dependent dielectric constant are in excellent agreement
with the results obtained with density functional perturbation theory. In
addition, we demonstrated that our formulation can be used to perform large-
scale molecular dynamics under electric fields as it captures electric-field
contributions to electric enthalpy, forces and polarization in close agreement
with reference DFT calculations. Overall, this work offers a promising new
direction in using machine learning techniques to accelerate the
investigations of dielectric and vibrational spectroscopic properties of
complex materials, including crystalline, disordered, and liquid materials,
defects dynamics under electric field, and ferroelectric response under
electric fields. Finally, ideas introduced in this work generalize readily to
differentiable learning strategies for a wide variety of thermodynamic
potential state functions, such as the free energy, grand potential, enthalpy
and others.

Acknowledgements
We thank L. Giacomazzi for insight on DFPT calculations. This work was
supported primarily by the US Department of Energy, Office of Basic Energy
Sciences Award No. DE-SC0022199 as well as by the Camille and Henry Dreyfus
Foundation Grant No. ML-22-075, the Department of Navy award N00014-20-1-2418
issued by the Office of Naval Research and Robert Bosch LLC. Computational
resources were provided by the Harvard University FAS Division of Science
Research Computing Group. C.J.O. was supported by the National Science
Foundation Graduate Research Fellowship Program under Grant No. (DGE1745303).
A.J. was supported by Aker Scholarship. Computational resources were provided
by the FAS Division of Science Research Computing Group at Harvard University.
Additional resources include the National Energy Research Scientific Computing
Center (NERSC), a DOE Office of Science User Facility supported by the Office
of Science of the U.S. Department of Energy under Contract No. DE-
AC02-05CH11231 using NERSC award BES-ERCAP0024206.

## References

  * Resta [1994] R. Resta, Macroscopic polarization in crystalline dielectrics: the geometric phase approach, Rev. Mod. Phys. 66, 899 (1994).
  * Resta and Vanderbilt [2007] R. Resta and D. Vanderbilt, Theory of polarization: a modern approach, in  _Physics of ferroelectrics: a modern perspective_ (Springer, 2007).
  * Resta [1992] R. Resta, Theory of the electric polarization in crystals, Ferroelectrics 136, 51 (1992).
  * King-Smith and Vanderbilt [1993] R. D. King-Smith and D. Vanderbilt, Theory of polarization of crystalline solids, Phys. Rev. B 47, 1651 (1993).
  * Vanderbilt and King-Smith [1993] D. Vanderbilt and R. D. King-Smith, Electric polarization as a bulk quantity and its relation to surface charge, Phys. Rev. B 48, 4442 (1993).
  * Tagantsev [1991] A. K. Tagantsev, Electric polarization in crystals and its response to thermal and elastic perturbations, Phase Transit. 35, 119 (1991).
  * Posternak  _et al._ [1990] M. Posternak, A. Baldereschi, A. Catellani, and R. Resta, Ab initio study of the spontaneous polarization of pyroelectric BeO, Phys. Rev. Lett. 64, 1777 (1990).
  * Vanderbilt [2000] D. Vanderbilt, Berry-phase theory of proper piezoelectric response, J. Phys. Chem. Solids 61, 147 (2000).
  * Spaldin [2012] N. A. Spaldin, A beginner’s guide to the modern theory of polarization, J. Solid State Chem. 195, 2 (2012).
  * Nunes and Gonze [2001] R. W. Nunes and X. Gonze, Berry-phase treatment of the homogeneous electric field perturbation in insulators, Phys. Rev. B 63, 155107 (2001).
  * Umari and Pasquarello [2002] P. Umari and A. Pasquarello, Ab initio molecular dynamics in a finite homogeneous electric field, Phys. Rev. Lett. 89, 157602 (2002).
  * Kirkwood [1939] J. G. Kirkwood, The dielectric polarization of polar liquids, J. Chem. Phys. 7, 911 (1939).
  * Fröhlich [1958] H. Fröhlich, Theory of dielectrics: dielectric constant and dielectric loss, Clarendon Press: Oxford (1958).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Rep. Progr. Phys. 29, 255 (1966).
  * Neumann [1983] M. Neumann, Dipole moment fluctuation formulas in computer simulations of polar systems, Molecular Physics 50, 841 (1983).
  * Neumann and Steinhauser [1983] M. Neumann and O. Steinhauser, On the calculation of the frequency-dependent dielectric constant in computer simulations, Chem. Phys. Lett. 102, 508 (1983).
  * Neumann  _et al._ [1984] M. Neumann, O. Steinhauser, and G. S. Pawley, Consistent calculation of the static and frequency-dependent dielectric constant in computer simulations, Mol. Phys. 52, 97 (1984).
  * Evans and Powles [1982] W. Evans and J. Powles, The computer simulation of the dielectric properties of polar liquids. the dielectric constant and relaxation of liquid hydrogen chloride, Mol. Phys. 45, 695 (1982).
  * McQuarrie [2000] D. McQuarrie, Statistical Mechanics; University Science books, Sausalito, CA , 222 (2000).
  * Batzner  _et al._ [2022] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nat. Commun. 13, 2453 (2022).
  * Musaelian  _et al._ [2023] A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth, and B. Kozinsky, Learning local equivariant representations for large-scale atomistic dynamics, Nat. Commun. 14, 579 (2023).
  * Kozinsky  _et al._ [2023] B. Kozinsky, A. Musaelian, A. Johansson, and S. Batzner, Scaling the leading accuracy of deep equivariant models to biomolecular simulations of realistic size, in  _Int. Conf. High Performa. Comput._ (Association for Computing Machinery, New York, NY, USA, 2023).
  * Batatia  _et al._ [2022] I. Batatia, D. P. Kovacs, G. Simm, C. Ortner, and G. Csanyi, MACE: higher order equivariant message passing neural networks for fast and accurate force fields, Adv. Neural Inf. Process. Syst. 35, 11423 (2022).
  * Nigam  _et al._ [2022] J. Nigam, M. J. Willatt, and M. Ceriotti, Equivariant representations for molecular Hamiltonians and N-center atomic-scale properties, J. Chem. Phys. 156, 014115 (2022).
  * Geiger and Smidt [2022] M. Geiger and T. Smidt, e3nn: Euclidean neural networks, (2022), arXiv:2207.09453 [cs.LG] .
  * Fedik  _et al._ [2022] N. Fedik, R. Zubatyuk, M. Kulichenko, N. Lubbers, J. S. Smith, B. Nebgen, R. Messerly, Y. W. Li, A. I. Boldyrev, K. Barros, _et al._ , Extending machine learning beyond interatomic potentials for predicting molecular properties, Nat. Rev. Chem. 6, 653 (2022).
  * Christensen  _et al._ [2019] A. S. Christensen, F. A. Faber, and O. A. von Lilienfeld, Operators in quantum machine learning: Response properties in chemical space, J. Chem. Phys. 150, 064105 (2019).
  * Veit  _et al._ [2020] M. Veit, D. M. Wilkins, Y. Yang, J. DiStasio, Robert A., and M. Ceriotti, Predicting molecular dipole moments by combining atomic partial charges and atomic dipoles, J. Chem. Phys. 153, 024113 (2020).
  * Krishnamoorthy  _et al._ [2021] A. Krishnamoorthy, K.-i. Nomura, N. Baradwaj, K. Shimamura, P. Rajak, A. Mishra, S. Fukushima, F. Shimojo, R. Kalia, A. Nakano, and P. Vashishta, Dielectric constant of liquid water determined with neural network quantum molecular dynamics, Phys. Rev. Lett. 126, 216403 (2021).
  * Staacke  _et al._ [2022] C. G. Staacke, S. Wengert, C. Kunkel, G. Csányi, K. Reuter, and J. T. Margraf, Kernel charge equilibration: efficient and accurate prediction of molecular dipole moments with a machine-learning enhanced electron density model, Mach. Learn.: Sci. Technol. 3, 015032 (2022).
  * Gigli  _et al._ [2022] L. Gigli, M. Veit, M. Kotiuga, G. Pizzi, N. Marzari, and M. Ceriotti, Thermodynamics and dielectric response of BaTiO3 by data-driven modeling, npj Comput. Mater. 8, 209 (2022).
  * Schütt  _et al._ [2021] K. Schütt, O. Unke, and M. Gastegger, Equivariant message passing for the prediction of tensorial properties and molecular spectra, in  _Int. Conf. Mach. Learn._, Vol. 139, edited by M. Meila and T. Zhang (PMLR, 2021) pp. 9377–9388.
  * Schienbein [2023] P. Schienbein, Spectroscopy from machine learning by accurately representing the atomic polar tensor, J. Chem. Theory Comput. 19, 705 (2023).
  * Shao  _et al._ [2023] X. Shao, L. Paetow, M. E. Tuckerman, and M. Pavanello, Machine learning electronic structure methods based on the one-electron reduced density matrix, Nat. Commun. 14, 10.1038/s41467-023-41953-9 (2023).
  * Zhang  _et al._ [2022] L. Zhang, H. Wang, M. C. Muniz, A. Z. Panagiotopoulos, R. Car, and W. E, A deep potential model with long-range electrostatic interactions, J. Chem. Phys. 156, 124107 (2022), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0083669/16725582/124107_1_online.pdf .
  * Joll  _et al._ [2024] K. Joll, P. Schienbein, K. M. Rosso, and J. Blumberger, Molecular dynamics simulation with finite electric fields using perturbed neural network potentials (2024), arXiv:2403.12319 [physics.chem-ph] .
  * Choudhary  _et al._ [2020] K. Choudhary, K. F. Garrity, V. Sharma, A. J. Biacchi, A. R. Hight Walker, and F. Tavazza, High-throughput density functional perturbation theory and machine learning predictions of infrared, piezoelectric, and dielectric responses, npj Comput. Mater. 6, 64 (2020).
  * Shimizu  _et al._ [2023] K. Shimizu, R. Otsuka, M. Hara, E. Minamitani, and S. Watanabe, Prediction of born effective charges using neural network to study ion migration under electric fields: applications to crystalline and amorphous Li3PO4, (2023), arXiv:2305.19546 [cond-mat.mtrl-sci] .
  * Wilkins  _et al._ [2019] D. M. Wilkins, A. Grisafi, Y. Yang, K. U. Lao, R. A. DiStasio, and M. Ceriotti, Accurate molecular polarizabilities with coupled cluster theory and machine learning, Proc. Natl. Acad. Sci. 116, 3401 (2019).
  * Sommers  _et al._ [2020] G. M. Sommers, M. F. Calegari Andrade, L. Zhang, H. Wang, and R. Car, Raman spectrum and polarizability of liquid water from deep neural networks, Phys. Chem. Chem. Phys. 22, 10592 (2020).
  * Grumet  _et al._ [2023] M. Grumet, C. von Scarpatetti, T. Bučko, and D. A. Egger, Delta machine learning for predicting dielectric properties and Raman spectra (2023), arXiv:2307.10578 [cond-mat.mtrl-sci] .
  * Fang  _et al._ [2023] M. Fang, S. Tang, Z. Fan, Y. Shi, N. Xu, and Y. He, A study of simulating Raman spectra for alkanes with a machine learning-based polarizability model (2023), arXiv:2301.13490 [physics.chem-ph] .
  * Takahashi  _et al._ [2020] A. Takahashi, Y. Kumagai, J. Miyamoto, Y. Mochizuki, and F. Oba, Machine learning models for predicting the dielectric constants of oxides based on high-throughput first-principles calculations, Phys. Rev. Mater. 4, 103801 (2020).
  * Kapil  _et al._ [2024] V. Kapil, D. P. Kovács, G. Csányi, and A. Michaelides, First-principles spectroscopy of aqueous interfaces using machine-learned electronic and quantum nuclear effects, Faraday Discuss. 249, 50 (2024).
  * Berger  _et al._ [2024] E. Berger, J. Niemelä, O. Lampela, A. H. Juffer, and H.-P. Komsa, Raman spectra of amino acids and peptides from machine learning polarizabilities (2024), arXiv:2401.14808 [physics.comp-ph] .
  * Souza  _et al._ [2002] I. Souza, J. Íñiguez, and D. Vanderbilt, First-principles approach to insulators in finite electric fields, Phys. Rev. Lett. 89, 117602 (2002).
  * Neto and Salinas [2007] A. M. F. Neto and S. R. Salinas,  _The physics of ferroelectrics-a modern perspective_ (Springer, 2007).
  * Pick  _et al._ [1970] R. M. Pick, M. H. Cohen, and R. M. Martin, Microscopic theory of force constants in the adiabatic approximation, Phys. Rev. B 1, 910 (1970).
  * Vandermause  _et al._ [2020] J. Vandermause, S. B. Torrisi, S. Batzner, Y. Xie, L. Sun, A. M. Kolpak, and B. Kozinsky, On-the-fly active learning of interpretable bayesian force fields for atomistic rare events, npj Comput. Mater. 6, 20 (2020).
  * Loose  _et al._ [2023] T. D. Loose, P. G. Sahrmann, T. S. Qu, and G. A. Voth, Coarse-graining with equivariant neural networks: A path toward accurate and data-efficient models, J. Phys. Chem. B 127, 10564 (2023).
  * Fu  _et al._ [2023] X. Fu, Z. Wu, W. Wang, T. Xie, S. Keten, R. Gomez-Bombarelli, and T. Jaakkola, Forces are not enough: Benchmark and critical evaluation for machine learning force fields with molecular simulations (2023), arXiv:2210.07237 [physics.comp-ph] .
  * Maxson and Szilvasi [2024] T. Maxson and T. Szilvasi, Transferable water potentials using equivariant neural networks (2024), arXiv:2402.16204 [physics.chem-ph] .
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun. 271, 108171 (2022).
  * Broughton  _et al._ [1997] J. Q. Broughton, C. A. Meli, P. Vashishta, and R. K. Kalia, Direct atomistic simulation of quartz crystal oscillators: Bulk properties and nanoscale devices, Phys. Rev. B 56, 611 (1997).
  * de Gironcoli  _et al._ [1989] S. de Gironcoli, S. Baroni, and R. Resta, Piezoelectric properties of III-V semiconductors from first-principles linear-response theory, Phys. Rev. Lett. 62, 2853 (1989).
  * Giannozzi  _et al._ [1991] P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni, Ab initio calculation of phonon dispersions in semiconductors, Phys. Rev. B 43, 7231 (1991).
  * Gonze [1997] X. Gonze, First-principles responses of solids to atomic displacements and homogeneous electric fields: Implementation of a conjugate-gradient algorithm, Phys. Rev. B 55, 10337 (1997).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B 55, 10355 (1997).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73, 515 (2001).
  * Pasquarello and Car [1997] A. Pasquarello and R. Car, Dynamical charge tensors and infrared spectrum of amorphous SiO2, Phys. Rev. Lett. 79, 1766 (1997).
  * Giacomazzi  _et al._ [2009] L. Giacomazzi, P. Umari, and A. Pasquarello, Medium-range structure of vitreous SiO2 obtained through first-principles investigation of vibrational spectra, Phys. Rev. B 79, 064202 (2009).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117, 1 (1995).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter 21, 395502 (2009).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * van Setten  _et al._ [2018] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, Comput. Phys. Commun. 226, 39 (2018).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, PyTorch: An imperative style, high-performance deep learning library, in  _Adv. Neural Inf. Process. Syst._, Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).
  * Kingma and Ba [2017] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2017), arXiv:1412.6980 [cs.LG] .
  * Marsalek and Markland [2017] O. Marsalek and T. E. Markland, Quantum dynamics and spectroscopy of ab initio liquid water: The interplay of nuclear and electronic quantum effects, J. Phys. Chem. Lett. 8, 1545 (2017).
  * Ramírez  _et al._ [2004] R. Ramírez, T. López-Ciudad, P. Kumar P, and D. Marx, Quantum corrections to classical time-correlation functions: Hydrogen bonding and anharmonic floppy modes, J. Chem. Phys. 121, 3973 (2004).

## Supplementary Information

###  IV.1 Computational details

Figure 5: Parity plots of energy, forces, polarization, Born charges, and
polarizability comparing the results obtained with DFT and with our unified ML
model. Each $\mu$ component of polarization is folded in the range
$[-\tfrac{1}{2}\Delta P_{\mu},+\tfrac{1}{2}\Delta P_{\mu}]$, where
$\Delta\mathbf{P}$ is the quantum of polarization. Mean absolute errors (MAE)
are provided in the inset of the plots. Energy, forces and polarization are
determined in the absence of the electric field.

We model $\alpha$-SiO2 with a 72-atom orthorhombic supercell ($a=9.97$ Å,
$b=8.63$ Å, $c=10.96$ Å). The classical molecular dynamics for generating the
data frames are performed with lammps [62]. We use the classical Vashishta
potential [54], and set the time step to 2 fs. The DFT calculations are
performed using a plane-wave density functional approach as implemented in the
Quantum Espresso suite [63]. We use the semilocal Perdew-Burke-Ernzerhof (PBE)
functional [64], and describe the core-valence interactions using
normconserving pseudopotentials [65]. We sample the Brillouin zone with a
2$\times$2$\times$2 k-grid and set the energy cutoff to 100 Ry. Born charges
and polarizability can be calculated through DFT calculations with finite
electric fields. In particular, the finite electric field is set to 10-4 Ry
a.u. along each Cartesian direction, which ensures a linear regime of
polarization with respect to electric field. The Born charges are calculated
through the following expression

$Z^{*}_{i\mu\nu}=\left.\frac{1}{e}\frac{dF_{i\nu}}{d\mathcal{E}_{\mu}}\right|_{\mathcal{E}_{\mu}=0},$
(8)

which derives from combining Eq. (2) and the definition of atomic forces
$F_{i\nu}=-dE/dr_{i\nu}$. Equation (8) can be evaluated by taking finite
differences of forces in the presence and absence of electric field,
respectively. Similarly, the polarizabilities are determined through Eq. (4)
by taking finite differences of polarization in the presence and absence of
the electric field, respectively. The DFPT calculations are performed using a
9-atom supercell and a 4$\times$4$\times$4 kgrid. The infrared spectrum and
the frequency-dependent dielectric constant from DFPT are computed at the
$\Gamma$-point of the Brillouin zone.

The Allegro code is compiled using PyTorch 1.11 [66]. We use 2-layers models
with full $O(3)$ symmetry, radial cutoff $r_{\text{max}}=6$ Å, maximum order
of spherical harmonics $\ell_{\text{max}}=3$, 64 tensor features, and 8 Bessel
functions to represent interatomic distances. In particular,
$r_{\text{max}}=6$ Å is found to capture most of the relevant peaks in the
radial distribution function of $\alpha$-SiO2. The network is constituted of
two-body latent multi-layer perceptron (MLP) with hidden dimensions
$[64,64,64]$, a latent MLP with dimensions $[64,64,64]$, a linear projection
as embedding MLP, and a single layer of dimension 128 as final MLP. SiLU
nonlinearities are applied to the two-body latent MLP and to the latent MLP.
The loss function weights are $\lambda_{E}=1$, $\lambda_{F}=1$,
$\lambda_{P}=1$, $\lambda_{Z}=10$, and $\lambda_{\alpha}=1000$. We fix the
batch size to 1, and the learning rate to 0.001. The scheduler patience is set
to 250. We use the Adam optimizer [67] in PyTorch, with default parameters of
$\nu_{1}=0.9$, $\nu_{2}=0.999$, and $\epsilon=10^{-8}$ without weight decay.
The dataset is split in 80% for training, and 20% for validation. The
validation of our model $\alpha$-SiO2 on the entire dataset yields mean
absolute errors of 0.21 eV on the energy, 27 meV$\cdot$Å-1 on forces, 10
me$\cdot$Å on polarization, 5$\cdot$10-3 on Born charges, and 12
me$\cdot$A${}^{2}\cdot$V-1 on polarizabilities. Parity plots comparing the DFT
labels with the ML predictions are illustrated in Fig. 5.

###  IV.2 Vibrational and dielectric properties from molecular dynamics

We show how vibrational and dielectric properties can be determined from
polarization and polarizability throughout a MD trajectory in the absence of
the electric field at a given temperature [12, 13, 14, 15, 16, 17, 18, 19].
Specifically, the infrared absorption spectrum can be calculated as [68, 69]

$I(\omega)\propto\omega^{2}\text{Re}\left[\int_{0}^{+\infty}dt\,e^{-i\omega
t}\braket{\mathbf{P}^{0}(t)\cdot\mathbf{P}^{0}(0)}\right],$ (9)

where $\omega$ is the frequency, $t$ the time, $\mathbf{P}^{0}$ the
polarization in the absence of an electric field, and
$\braket{\mathbf{P}^{0}(t)\cdot\mathbf{P}^{0}(0)}$ the average of the
autocorrelation function of polarization. The polarizability can be used to
determine the high-frequency dielectric constant through the following
expression

$\varepsilon^{\infty}_{\mu\nu}=1+\frac{4\pi}{\Omega}\braket{\alpha_{\mu\nu}},$
(10)

where $\Omega$ is the volume, and $\braket{\alpha_{\mu\nu}}$ the average
polarizability. The static dielectric constant can then be determined by
adding an ionic contribution related to the polarization during the MD though
the fluctuation-dissipation theorem, namely

$\varepsilon^{0}_{\mu\nu}=\varepsilon^{\infty}_{\mu\nu}+\frac{4\pi}{\Omega}\frac{\text{cov}(P^{0}_{\mu},P^{0}_{\nu})}{k_{\text{B}}T},$
(11)

where $k_{\text{B}}$ is the Boltzmann constant, $T$ the temperature, and
$\text{cov}(P^{0}_{\mu},P^{0}_{\nu})$ the covariance of $P^{0}_{\mu}$ and
$P^{0}_{\nu}$. Then, one can determine the frequency-dependent dielectric
constant using the autocorrelation function of polarization as follows:

$\displaystyle\varepsilon_{\mu\nu}(\omega)=1+(\varepsilon^{0}_{\mu\nu}-1)\left[1-i\omega\int_{0}^{+\infty}dt\,e^{-i\omega
t}\frac{\braket{P^{0}_{\mu}(t)P^{0}_{\nu}(0)}}{\text{cov}(P^{0}_{\mu},P^{0}_{\nu})}\right].$
(12)

In Fig. 4, the dielectric constant is calculated as

$\varepsilon_{\mu\nu}=1+\frac{4\pi}{\Omega}\frac{P_{\mu}(\mathbf{R},\bm{\mathcal{E}})-P_{\mu}(\mathbf{R}^{0},\mathbf{0})}{\mathcal{E}_{\nu}},$
(13)

where $P_{\mu}(\mathbf{R},\bm{\mathcal{E}})$ is the polarization obtained in
the presence of an electric field $\bm{\mathcal{E}}$ for a given structure
$\mathbf{R}$, and $P_{\mu}(\mathbf{R}^{0},\bm{0})$ is the polarization
obtained in the absence of the field for the initial structure
$\mathbf{R}^{0}$. The structure $\mathbf{R}^{0}$ is found by performing DFT
relaxations in the absence of the electric field. Equations (10), (11) and
(13) are written in atomic units. In eV units, $4\pi$ is replaced with
$1/\epsilon_{0}$, where $\epsilon_{0}$ is the vacuum permittivity.