# Yang–Mills Monopoles in
Extremal Reissner–Nordström Black Hole Metric

Xiangqin Zhang
School of Mathematics and Statistics
Henan University
Kaifeng, Henan 475004, PR China

Yisong Yang111Email address for correspondence<EMAIL_ADDRESS>
Courant Institute of Mathematical Sciences
New York University
New York, NY 10012, USA

###### Abstract

We show that the closed-form monotone solution linking two different vacuum
states, for the exterior Yang–Mills wave equation over the extremal
Reissner–Nordström spacetime found in a recent work of Bizoń and Kahl, is the
unique solution to an associated domain-wall type energy minimization problem.
We also derive the accompanying interior Yang–Mills wave equation within the
same formalism. We obtain a few more closed-form solutions, regular and
singular, among other solutions of oscillatory behavior, and discuss several
interesting features of the solutions based on some energy consideration.

Keywords: Reissner–Nordström black hole metric, Yang–Mills fields, monopoles,
domain-wall solitons, minimization problem, existence and uniqueness of
solution.

PACS numbers: 02.30.Hq, 11.15.q, 11.27.+d, 12.39.Ba

MSC numbers: 34B40, 35J50, 81T13

##  1 Introduction

The concept of magnetic monopoles was first conceived by P. Curie [9] and
later theoretically formulated by Dirac [10] based on electromagnetic duality
of the Maxwell theory. A remarkable feature of this concept is that its
existence could lead to an explanation why electric charges are multiples of a
unit charge [10, 13, 25, 27, 39]. Such a charge quantization phenomenon was
then extended by Schwinger [28] to dyons, namely, particles of both electric
and magnetic charges. In modern development, thought experiments with
monopoles immersed in a type-II superconductor led Mandelstam [19, 20], Nambu
[23], and t’ Hooft [33, 33] to come up with a linear confinement mechanism for
monopoles and their colored, non-Abelian or the Yang–Mills, field theory
versions, which has propelled enormous activities and progress towards an
understanding [29, 30] of the unsettled puzzle of quark confinement [14]. Like
that in the classical Maxwell theory, a monopole in the pure Yang–Mills
theory, for example, the celebrated Wu–Yang monopole [1, 13], is of infinite
energy [16]. There are several well-recognized paths which may be taken to
overcome this infinity or divergence problem. These include adopting the model
of the Born–Infeld nonlinear electrodynamics [4, 5], coupling the Yang–Mills
fields with matter fields [1, 13, 17, 24, 25, 32, 35, 38, 40], namely,
considering the Yang–Mills–Higgs theory, and hosting the Yang–Mills fields
over curved spacetime backgrounds due to gravity. In the context of this last
approach, if a gravitational black hole is present and fields are considered
outside the event horizon so that possible singularities leading to energy
divergence are concealed inside the horizon, as in the scenario of the cosmic
censorship hypothesis [22, 37], the total energy of the fields outside the
horizon could remain finite, as a result. The present study belongs to this
category of investigation. Specifically our work is motivated by the recent
paper of Bizoń and Kahl [3] who considered spherically symmetric $SU(2)$
Yang–Mills fields in the exterior of the event horizon of the extremal
Reissner–Nordström black hole and obtained a countable family of static
regular solutions to the Yang–Mills equations, given in terms of a scalar
field sequence $\\{W_{n}\\}$, say, among which the first two, $W_{1}$ and
$W_{2}$, are in closed forms. These solutions resemble the Bartnick–McKinnon
solutions [2, 31] featured as soliton-like configurations connecting vacuum
states at infinity and oscillating about the unit-charge monopole state in
local domains. The main result of our present work is to show that, viewed as
a connecting orbit between different vacuum states which passes through the
unit-charge monopole state midway, $W_{1}$ is the unique global energy
minimizer of the problem. The solution $W_{2}$, however, enjoys a different
signature – it exhibits itself as a local energy maximizer in a well-
formulated sense. Besides, we shall also obtain a new closed-form solution
which has two singular points (in fact, two singular shells), near which
energy blow-up takes place. Furthermore, we adapt the formalism of [3] to
consider interior monopoles (inside the event horizon) and derive a governing
wave equation which may be viewed as a companion of the wave equation of
Bizoń–Kahl [3] for exterior monopoles. For this interior monopole equation, we
obtain three closed-form solutions, similar to $W_{1}$, $W_{2}$, and the third
singular solution to the exterior monopole equation described above, among
which two are regular everywhere away from the curvature singularity of the
Reissner–Nordström spacetime and the third one is singular at one point (in
fact, one shell). It is noted that this interior monopole equation possesses a
family of oscillatory solutions as the exterior equation, as well.

The content of the rest of the paper is outlined as follows. In Section 2, we
formulate the problems to be studied with a review of the Yang–Mills fields
over an extremal Reissner–Nordström spacetime and the exterior monopole
equation derived in [3], followed with an extension of such a formulation to
get the interior, ‘accompanying’, monopole equation. Subsequently, we focus on
static equations. In Section 3, we review the two closed-form regular
solutions, $W_{1}$ and $W_{2}$, obtained in [3], and present a third,
singular, one, for the exterior equation. We obtain the exact values of the
energies of these two regular solutions to be further explored in later
sections. We then present two closed-form regular solutions and a third,
singular, one, for the interior equation, similar to those for the exterior
equation. In Section 4, we study several issues of the exterior solutions
based on some energy consideration. In particular, we formulate a domain-wall
type minimization problem that associates a distinguishing significance to
$W_{1}$ as the unique global minimizer of the problem. We also present a
virial identity which displays an energy partition relation for finite-energy
solutions. In Section 5, we first prove the existence of a solution to the
domain-wall minimization problem formulated in Section 4. Note that, due to
the behavior of the weight function in the energy functional, it is not a
trivial matter to preserve the asymptotic limit of the weak limit of a
minimizing sequence. Fortunately, such a difficulty may be overcome by a local
convergence method leading to realizing the weak limit as a monotone solution
to the exterior equation and then by an energy comparison. We next prove the
uniqueness of a monotone solution to the exterior equation. Since the energy
functional is not convex, such a uniqueness result is usually not ensured [11]
and seems surprising in the present context. Consequently, we recognize that
$W_{1}$ is actually this energy minimizer whose energy calculated in Section 3
is thus the energy minimum of the formulated problem. We also illustrate that
$W_{2}$ is a local energy maximizer among a class of testing configurations,
which seems to suggest that, energetically, the closed-form solution $W_{2}$
of the exterior monopole equation has a certain ‘significance’ associated as
well. In Section 6, we briefly comment on how to obtain oscillatory solutions
for both exterior and interior monopole equations using their polar-variable
representations.

##  2 Yang–Mills monopole equations in Reissner–Nordström metric

Consider the Reissner–Nordström metric [22, 37]

$\mathrm{d}s^{2}=-\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)\,\mathrm{d}t^{2}+\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)^{-1}\,\mathrm{d}r^{2}+r^{2}(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2}),$
(2.1)

describing a charged black hole of mass $M$ and charge $Q$ in the
Schwarzschild spherical coordinates. At the critical situation

$M=Q,$ (2.2)

the metric (2.1) assumes the form

$\mathrm{d}s^{2}=-\frac{(r-M)^{2}}{r^{2}}\mathrm{d}t^{2}+\frac{r^{2}}{(r-M)^{2}}\,\mathrm{d}r^{2}+r^{2}(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2}).$
(2.3)

known as the extremal Reissner–Nordström black hole metric with its event
horizon at $r=M$. We will be interested in both the exterior region, $r>M$,
and interior one, $r<M$, possibly extended to its curvature singularity at
$r=0$ as exhibited by the associated Kretschmann invariant given by

$K=R_{\alpha\beta\mu\nu}R^{\alpha\beta\mu\nu}=\frac{8M^{2}}{r^{8}}(7M^{2}-12Mr+6r^{2}),$
(2.4)

composed from the Riemann tensor induced by the gravitational metric form
(2.1).

###  2.1 Exterior region

We first consider the space region exterior to the event horizon, $r>M$. As in
[3], we use the variables

$\tau=\frac{t}{4M},\quad x=\ln\left(\frac{r}{M}-1\right),$ (2.5)

to recast (2.3) into

$\displaystyle\mathrm{d}s^{2}$ $\displaystyle=$
$\displaystyle\Omega(x)\left(-\mathrm{d}\tau^{2}+C^{4}(x)(\mathrm{d}x^{2}+\mathrm{d}\theta^{2}+\sin^{2}\theta\,\mathrm{d}\phi^{2})\right),$
(2.6) $\displaystyle\Omega(x)$ $\displaystyle=$
$\displaystyle\frac{16M^{2}}{(1+\mathrm{e}^{-x})^{2}}\quad
C(x)=\cosh\left(\frac{x}{2}\right).$ (2.7)

Thus the exterior region in terms of the radial variable, $r>M$, is converted
in terms of the variable $x$ to the full line $-\infty<x<\infty$, such that
the ‘midway’ spot, $x=0$, where the conformal factor $C(x)$ is minimized,
corresponds to the classical Schwarzschild radius

$r_{\mbox{\small s}}=2M,$ (2.8)

which will have a specific meaningfulness in our study.

We consider the $SU(2)$ Yang–Mills theory over such a spacetime. Use
$\sigma_{a}$ ($a=1,2,3$) to denote the Pauli matrices and set

$t_{a}=\frac{\sigma_{a}}{2\mathrm{i}},\quad a=1,2,3.$ (2.9)

Then $\\{t_{a}\\}$ generates $SU(2)$ and satisfies the commutator relation
$[t_{a},t_{b}]=\epsilon_{abc}t_{c}$. As in [3], we choose the gauge field to
be the 1-form

$A=W(\tau,x)\,\omega+t_{3}\cos\theta\,\mathrm{d}\phi,\quad\omega=t_{1}\,\mathrm{d}\theta+t_{2}\sin\theta\,\mathrm{d}\phi,$
(2.10)

where $W$ is a scalar field. Then we compute to get the Yang–Mills curvature
field

$\displaystyle F$ $\displaystyle=$ $\displaystyle\mathrm{d}A+A\wedge A$ (2.11)
$\displaystyle=$ $\displaystyle
t_{1}(W_{\tau}\,\mathrm{d}\tau+W_{x}\,\mathrm{d}x)\wedge\mathrm{d}\theta+t_{2}\sin\theta(W_{\tau}\,\mathrm{d}\tau+W_{x}\,\mathrm{d}x)\wedge\mathrm{d}\phi+t_{3}\sin\theta(W^{2}-1)\,\mathrm{d}\theta\wedge\mathrm{d}\phi.\quad$

Therefore, with $x^{0}=\tau,x^{1}=x,x^{2}=\theta,x^{3}=\phi$ and
$F=F^{a}_{\mu\nu}t_{a}\,\mathrm{d}x^{\mu}\wedge\mathrm{d}x^{\nu}$, we can read
(2.11) to obtain the nonzero and independent components of the 2-form $F$ to
be

$F^{1}_{02}=\frac{1}{2}W_{\tau},\quad F^{1}_{12}=\frac{1}{2}W_{x},\quad
F^{2}_{03}=\frac{1}{2}\sin\theta W_{\tau},\quad
F^{2}_{13}=\frac{1}{2}\sin\theta W_{x},\quad
F^{3}_{23}=\frac{1}{2}\sin\theta(W^{2}-1).$ (2.12)

On the other hand, if we rewrite line element (2.6)–(2.7) as
$\mathrm{d}s^{2}=g_{\mu\nu}\,\mathrm{d}x^{\mu}\mathrm{d}x^{\nu}$, then the
Yang–Mills action density reads

$\displaystyle{\cal L}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\,F^{a}_{\mu\nu}g^{\mu\alpha}g^{\nu\beta}F^{a}_{\alpha\beta}\sqrt{|\det(g_{\alpha\beta})|}$
(2.13) $\displaystyle=$
$\displaystyle\left(-\frac{C^{2}}{4}W_{\tau}^{2}+\frac{1}{4C^{2}}W_{x}^{2}+\frac{1}{8C^{2}}(1-W^{2})^{2}\right)\sin\theta,$

in view of (2.12). Hence we can rewrite the Yang–Mills action as

$L=\int{\cal L}\,\mathrm{d}^{4}x=2\pi\int{\cal
L}_{W}\,\mathrm{d}x\mathrm{d}\tau,$ (2.14)

where

${\cal
L}_{W}=-\frac{C^{2}}{2}W_{\tau}^{2}+\frac{1}{2C^{2}}W_{x}^{2}+\frac{1}{4C^{2}}(1-W^{2})^{2},$
(2.15)

as obtained in [3], so that the Yang–Mills equation is the Euler–Lagrange
equation of (2.15):

$W_{\tau\tau}=C^{-2}(C^{-2}W_{x})_{x}+C^{-4}W(1-W^{2}),$ (2.16)

whose static limit reads

$W^{\prime\prime}-\tanh\left(\frac{x}{2}\right)W^{\prime}+W(1-W^{2})=0,$
(2.17)

where and in the sequel, we use the notation $W^{\prime}=W_{x}$, etc.,
interchangeably.

###  2.2 Interior region

We next consider the region that is interior to the event horizon, $0<r<M$.
Thus, using the variables

$\tau=\frac{t}{4M},\quad x=\ln\left(1-\frac{r}{M}\right),\quad 0<r<M,$ (2.18)

we have $-\infty<x<0$ and we see that the line element (2.3) becomes

$\displaystyle\mathrm{d}s^{2}$ $\displaystyle=$
$\displaystyle\Omega(x)\left(-\mathrm{d}\tau^{2}+S^{4}(x)(\mathrm{d}x^{2}+\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2})\right),$
(2.19) $\displaystyle\Omega(x)$ $\displaystyle=$
$\displaystyle\frac{16M^{2}}{(1-\mathrm{e}^{-x})^{2}},\quad
S(x)=\sinh\left(\frac{x}{2}\right).$ (2.20)

Therefore, from (2.10) and (2.11), we analogously obtain the effective
Yang–Mills action density

${\cal
L}_{W}=-\frac{S^{2}}{2}W_{\tau}^{2}+\frac{1}{2S^{2}}W_{x}^{2}+\frac{1}{4S^{2}}(1-W^{2})^{2},$
(2.21)

over the spatial interval $-\infty<x<0$, with the associated Euler–Lagrange
equation

$W_{\tau\tau}={S^{-2}}(S^{-2}W_{x})_{x}+{S^{-4}}(1-W^{2})W,$ (2.22)

whose static limit reads

$W^{\prime\prime}-\coth\left(\frac{x}{2}\right)W^{\prime}+W(1-W^{2})=0.$
(2.23)

The equations (2.22) and (2.23) are a new pair of equations, accompanying
(2.16) and (2.17).

##  3 Exact solutions

In [3], the following two nontrivial exact solutions to the exterior
Yang–Mills monopole equation (2.17) are obtained:

$W_{1}(x)=\tanh\left(\frac{x}{2}\right),\quad W_{2}(x)=\frac{2\cosh
x-2-\sqrt{6}}{2\cosh x+4+3\sqrt{6}}.$ (3.1)

Here we find that there is a third exact solution,

$W_{3}(x)=\frac{2\cosh x-2+\sqrt{6}}{2\cosh x+4-3\sqrt{6}},$ (3.2)

which is apparently singular at

$x_{\pm}=\pm\ln\left(\frac{3\sqrt{6}}{2}+2\sqrt{3}-2-\frac{3\sqrt{2}}{2}\right)\approx\pm
1.104268225,$ (3.3)

giving rise to the corresponding singular radii

$r_{\pm}=\left(1+\mathrm{e}^{x_{\pm}}\right)M,$ (3.4)

in the original radial variable, which are distributed about the Schwarzschild
radius $r_{\mbox{\small s}}$ following

$M<r_{-}<r_{\mbox{\small
s}}=2M<r_{+},\quad\quad\frac{r_{-}+r_{+}}{2}=\left(1+\cosh(x_{+})\right)M>r_{\mbox{\small
s}}.$ (3.5)

For the interior Yang–Mills monopole equation (2.23), we likewise find the
following three nontrivial exact solutions:

$\displaystyle w_{1}(x)$ $\displaystyle=$
$\displaystyle\coth\left(\frac{x}{2}\right),$ (3.6) $\displaystyle w_{2}(x)$
$\displaystyle=$ $\displaystyle\frac{2\cosh x+2-\sqrt{6}}{2\cosh
x-4+3\sqrt{6}},$ (3.7) $\displaystyle w_{3}(x)$ $\displaystyle=$
$\displaystyle\frac{2\cosh x+2+\sqrt{6}}{2\cosh x-4-3\sqrt{6}},$ (3.8)

for $-\infty<x<0$. The solutions $w_{1}$ and $w_{2}$ given in (3.6) and (3.7),
respectively, are regular but $w_{3}$ in (3.6) is singular at

$x_{0}=-\ln\left(2[1+\sqrt{3}]+\frac{3}{\sqrt{2}}[1+\sqrt{3}]\right),$ (3.9)

whereby the irrelevant positive root is discarded.

We now compute the energies of these solutions.

Recall that the energy-momentum tensor induced from (2.13) is

$T_{\mu\nu}=-F^{a}_{\mu\alpha}g^{\alpha\beta}F^{a}_{\nu\beta}+\frac{1}{4}g_{\mu\nu}(F^{a}_{\gamma\delta}g^{\gamma\alpha}g^{\delta\beta}F^{a}_{\alpha\beta}),$
(3.10)

which via ${\cal H}=T^{0}_{0}=g^{00}T_{00}$ gives rise to the Hamiltonian
energy density associated with (2.6)–(2.7) and (2.12) to be

${\cal
H}_{W}=\frac{C^{2}}{2}W_{\tau}^{2}+\frac{1}{2C^{2}}W_{x}^{2}+\frac{1}{4C^{2}}(1-W^{2})^{2}.$
(3.11)

Thus, evaluating the static total energy

$E(W)=2\pi\int_{-\infty}^{\infty}\left(\frac{1}{2C^{2}}W_{x}^{2}+\frac{1}{4C^{2}}(1-W^{2})^{2}\right)\,\mathrm{d}x,$
(3.12)

we obtain the exact value

$E(W_{1})=\frac{8\pi}{5},$ (3.13)

for the solution $W_{1}$ given in (3.1). For the solution $W_{2}$ given in
(3.1), utilizing MAPLE, we obtain the exact value

$E(W_{2})=2\pi\frac{(12198704\sqrt{2}\sqrt{3}+30388736)\tanh^{-1}\left(\frac{5}{3\sqrt{2}+\sqrt{3}}\right)+15707196\sqrt{2}+13222288\sqrt{3}}{35368572\sqrt{3}+43079049\sqrt{2}},$
(3.14)

which is approximately $2\pi(0.9664005272)$. These two exact results are
consistent with those found in [3] based on numerical evaluation. Not
surprisingly, the singular solution $W_{3}(x)$ is of infinite energy. In other
words, this solution behaves singularly in both point-wise and energy-wise
manners at $x_{\pm}$.

We now turn our attention to the interior solutions listed in (3.6)–(3.8). In
view of (3.10), (2.12), and (2.19)–(2.20), we get the Hamiltonian energy
density

${\cal
H}_{W}=\frac{S^{2}}{2}W_{\tau}^{2}+\frac{1}{2S^{2}}W_{x}^{2}+\frac{1}{4S^{2}}(1-W^{2})^{2},$
(3.15)

which renders the total energy for a static interior solution in the full
region $0<r<M$ to be

$E=2\pi\int_{-\infty}^{0}\left(\frac{1}{2S^{2}}W_{x}^{2}+\frac{1}{4S^{2}}(1-W^{2})^{2}\right)\,\mathrm{d}x.$
(3.16)

All the solutions given in (3.6)–(3.8) are of infinity energy due to the
divergence of the energy density (3.15) as a consequence of the curvature
singularity of the Reissner–Nordström black hole spacetime at $r=0$ displayed
by (2.4). Nevertheless, we may assume that we are in a situation that gravity
is generated from a bulk of matter which occupies a finite region within the
event horizon,

$r\leq r_{\mbox{\tiny B}},\quad 0<r_{\mbox{\tiny B}}<M.$ (3.17)

Then we can focus on the ‘regularized’ region

$-\infty<x<x_{\mbox{\tiny B}}<0,\quad x_{\mbox{\tiny
B}}=\ln\left(1-\frac{r_{\mbox{\tiny B}}}{M}\right),\quad r_{\mbox{\tiny
B}}\leq r<M,$ (3.18)

instead. It is seen that the solutions (3.6) and (3.7) are of finite energies
over (3.18) for any $0<r_{\mbox{\tiny B}}<M$ and the solution (3.8) is of
finite energy only when $x_{\mbox{\tiny B}}$ is below its singular point
$x_{0}$ spelled out in (3.9). Thus, we are led to the finite-energy condition

$x_{\mbox{\tiny
B}}<-\ln\left(2[1+\sqrt{3}]+\frac{3}{\sqrt{2}}[1+\sqrt{3}]\right)\approx-2.421226122,$
(3.19)

or equivalently,

$r_{\mbox{\tiny
B}}>\left(\frac{(3+3\sqrt{3})\sqrt{2}+4\sqrt{3}+2}{(3+3\sqrt{3})\sqrt{2}+4\sqrt{3}+4}\right)M\approx(0.9111873445)M.$
(3.20)

Over a given interval (3.18), the energy of the solution (3.6) is generically
much greater than that of the solution (3.7), since the former is singular at
its right boundary point, $x=0$, and the latter is everywhere regular. As a
comparison, we list our calculation of the energy

$E(W)=2\pi\int_{-\infty}^{x_{\mbox{\tiny
B}}}\left(\frac{1}{2S^{2}}W_{x}^{2}+\frac{1}{4S^{2}}(1-W^{2})^{2}\right)\,\mathrm{d}x,$
(3.21)

for $x_{\mbox{\tiny B}}=-1$ for the solutions $w_{1}$ and $w_{2}$, given in
(3.6) and (3.7), respectively, as follows:

$E(w_{1})=2\pi(3.273953196),\quad E(w_{2})=2\pi(0.3220989496).$ (3.22)

##  4 Energy consideration of exterior monopoles

We now investigate the problem whether there is a nontrivial energetically
stable exterior monopole solution linking two vacuum states, $W=\pm 1$,
characterized by the boundary condition $W(-\infty)=\pm 1$ and $W(\infty)=\pm
1$. It is clear that the problem of linking the same vacuum state with
$W(-\infty)=W(\infty)=\pm 1$ is not well posed without some further
elaboration since the vacuum states themselves minimize the energy and may be
approached arbitrarily by nontrivial field configurations. Thus it remains to
consider the problem of linking two different vacuum states, say,

$\eta=\inf\\{E(W)\,|\,W(-\infty)=-1,W(\infty)=1\\},$ (4.1)

and ask whether this minimization problem has a solution, where $E(W)$ is as
given in (3.12). This is typically a domain wall problem [6, 7, 8, 26, 36],
interpolating two ground state domains, $W=\pm 1$. However, it is not hard to
see that this problem is not well posed either.

In fact, let $u(x)$ be a smooth function satisfying

$u(x)=-1,\quad x\leq-1;\quad u(x)=1,\quad x\geq 1.$ (4.2)

Then, using $u_{(a)}=u(x-a)$ as a testing function, we have

$\displaystyle E(u_{(a)})$ $\displaystyle=$
$\displaystyle\pi\int_{a-1}^{a+1}\frac{1}{C^{2}(x)}\left((u^{\prime}_{(a)})^{2}+\frac{1}{2}(u^{2}_{(a)}-1)^{2}\right)\,\mathrm{d}x$
(4.3) $\displaystyle\leq$ $\displaystyle
K_{0}\int_{a-1}^{a+1}\frac{\mathrm{d}x}{C^{2}(x)}\to 0\quad\mbox{as
}|a|\to\infty.$

Here $K_{0}>0$ is a constant depending on the properties of $u$ over $[-1,1]$.
As a consequence, this shows the quantity $\eta$ in (4.1) is zero which is
therefore not attainable.

On the other hand, since the boundary condition in (4.1) dictates that the
unit magnetic charge monopole configuration, $W=0$, is to occur somewhere as
an intermediate state, we may conveniently impose the following additional
‘intermediate state condition’

$W(0)=0.$ (4.4)

It is interesting that this condition implies that the potential part of the
energy density

${\cal H}_{W}=\frac{1}{2C^{2}(x)}W_{x}^{2}+\frac{1}{4C^{2}(x)}(1-W^{2})^{2}$
(4.5)

is maximized at the spot where the intermediate state occurs, $x=0$. In Figure
1, we plot the energy density (4.5) for the solutions $W_{1}$ and $W_{2}$
given in (3.1) as an illustration. It is interesting to notice that the former
peaks sharply at $x=0$ where the solution passes the intermediate state and
the latter also peaks at the same spot which is not where the solution passes
the intermediate state but is the midpoint between the two spots,

$x^{\pm}_{0}=\ln\left(1+\frac{\sqrt{6}}{2}\pm\frac{\sqrt{6+4\sqrt{6}}}{2}\right),$
(4.6)

where the solution passes the intermediate state. Such a property of energy
concentration is typical for solitons in field theories.

Figure 1: Plots of the energy distributions of the solutions linking two
different vacuum states and the same vacuum state and the unit-charge monopole
state viewed as an intermediate state. The former is depicted by a solid-line
curve showing the concentration of energy at the spot where the solution
passes through the intermediate state. The latter is drawn by a dash-line
curve exhibiting the concentration of energy in the middle of the two spots
where the solution passes through the intermediate state. Energetically both
solutions demonstrate themselves as domain-wall solitons.

These pictures coincide with what we know about a domain wall soliton
configuration in general.

Thus, we are led to modifying (4.1) into the minimization problem

$\eta_{0}\equiv\inf\\{E(W)\,|\,W(-\infty)=-1,W(\infty)=1,W(0)=0\\}.$ (4.7)

We shall show that $W_{1}$ given in (3.1) is the unique solution to this
problem. Consequently, as a by-product, we have

$\eta_{0}=\frac{8\pi}{5},$ (4.8)

in view of (3.13).

Let $W$ be a finite-energy critical point of the energy functional (3.12). Set
$W_{\lambda}(x)=W(\lambda x)$. Then, from
$\frac{\mathrm{d}}{\mathrm{d}\lambda}E(W_{\lambda})|_{\lambda=1}=0$, we get
the virial identity

$\int_{-\infty}^{\infty}\frac{1}{C^{2}(x)}\left(\left[1+x\tanh\left(\frac{x}{2}\right)\right]W_{x}^{2}+\frac{x}{2}\,\tanh\left(\frac{x}{2}\right)(W^{2}-1)^{2}\right)\,\mathrm{d}x=\int_{-\infty}^{\infty}\frac{(W^{2}-1)^{2}}{2C^{2}(x)}\,\mathrm{d}x,$
(4.9)

which is an energy partition relation. This identity demonstrates that the
‘potential energy’ is much greater than the ‘elastic energy’ for an exterior
monopole. For example, for the exact solution $W=W_{1}$ in (3.1), we have

$\int_{-\infty}^{\infty}\frac{(W^{2}-1)^{2}}{4C^{2}(x)}\,\mathrm{d}x=\frac{8}{15},\quad\int_{-\infty}^{\infty}\frac{W_{x}^{2}}{2C^{2}(x)}\,\mathrm{d}x=\frac{4}{15},$
(4.10)

and, for $W=W_{2}$ given in (3.1), these quantities are approximately
$0.8765154575$ and $0.08988507022$, respectively.

##  5 Existence and uniqueness of exterior energy minimizer

The symmetry of the functional (3.12) indicates that the problem (4.7) amounts
to considering the minimization problem

$I_{0}\equiv\inf\left\\{I(W)\equiv\int_{0}^{\infty}\frac{1}{C^{2}(x)}\left((W^{\prime})^{2}+\frac{1}{2}(W^{2}-1)^{2}\right)\,\mathrm{d}x\,\bigg{|}\,W(0)=0,W(\infty)=1\right\\},$
(5.1)

over the set of admissible functions which are absolutely continuous on all
compact subintervals of the half line $[0,\infty)$.

Let $\\{{\cal W}_{n}\\}$ be a minimizing sequence of (5.1) satisfying

$\lim_{n\to\infty}I({\cal W}_{n})=I_{0};\quad I({\cal W}_{n})\leq
I_{0}+\frac{1}{n},\quad n=1,2,\dots.$ (5.2)

By modifying the sequence $\\{{\cal W}_{n}\\}$ if necessary, we may assume
that $\\{{\cal W}_{n}\\}$ enjoys the property that each function ${\cal
W}_{n}$ is nondecreasing and $0\leq{\cal W}_{n}(x)\leq 1$. Furthermore, for
each ${\cal W}_{n}$, we may modify it such that it minimizes the partial
functional

$I_{n}(W)\equiv\int_{0}^{n}\frac{1}{C^{2}(x)}\left((W^{\prime})^{2}+\frac{1}{2}(W^{2}-1)^{2}\right)\,\mathrm{d}x,$
(5.3)

among functions fulfilling the boundary condition $W(0)=0,W(n)={\cal
W}_{n}(n)$. Therefore ${\cal W}_{n}$ when restricted to the interval $[0,n]$
is a critical point of $I_{n}$. As a consequence, we have

$\int_{0}^{n}\frac{1}{C^{2}(x)}\left({\cal W}^{\prime}_{n}w^{\prime}+({\cal
W}_{n}^{2}-1){\cal W}_{n}w\right)\,\mathrm{d}x=0,\quad\forall w\in
C_{0}^{1}(0,n).$ (5.4)

Extracting a suitable subsequence if necessary, we may assume that $\\{{\cal
W}_{n}\\}$ is weakly convergent in the Sobolev space $W^{1,2}(0,a)$, with the
weighted measure $C^{-2}(x)\,\mathrm{d}x$, for any $a>0$. Let ${{\cal W}}$
denote such a weak limit. Fix $a>0$. Neglecting the first few terms of the
sequence if necessary, we may also assume $n>a$. Thus (5.4) gives us

$\int_{0}^{a}\frac{1}{C^{2}(x)}\left({\cal W}^{\prime}_{n}w^{\prime}+({\cal
W}_{n}^{2}-1){\cal W}_{n}w\right)\,\mathrm{d}x=0,\quad\forall w\in
C_{0}^{1}(0,a).$ (5.5)

Using the compact embedding from the Sobolev space into the space $C[0,a]$ and
letting $n\to\infty$ in (5.5), we arrive at

$\int_{0}^{a}\frac{1}{C^{2}(x)}\left({{\cal W}}^{\prime}w^{\prime}+({{\cal
W}}^{2}-1){{\cal W}}w\right)\,\mathrm{d}x=0,\quad\forall w\in C_{0}^{1}(0,a),$
(5.6)

which leads us to conclude by the standard elliptic regularity theory [12, 18]
that ${{\cal W}}$ is a classical solution to the Euler–Lagrange equation of
the functional $I$. Hence ${\cal W}$ satisfies

$\left(\frac{{\cal
W}^{\prime}}{C^{2}(x)}\right)^{\prime}+\frac{1}{C^{2}(x)}(1-{\cal W}^{2}){\cal
W}=0,\quad{\cal W}(0)=0.$ (5.7)

Since $\\{{\cal W}_{n}\\}$ is a nondecreasing sequence with values confined in
$[0,1]$, so is its limit ${\cal W}$. However, since ${\cal W}$ satisfies
(5.7), we see that ${\cal W}\equiv 0$ or ${\cal W}$ is strictly increasing. We
now show that the former does not occur.

In fact, by using weak convergence, we have

$I_{a}({\cal W})\leq\liminf_{n\to\infty}I_{a}({\cal
W}_{n})\leq\lim_{n\to\infty}I({\cal W}_{n})=I_{0},$ (5.8)

for any $a>0$. Thus, letting $a\to\infty$ in (5.8), we get

$I({\cal W})\leq I_{0}.$ (5.9)

Recall the result

$I_{0}\leq I\left(\tanh\,\frac{x}{2}\right)=\frac{4}{5}.$ (5.10)

Besides, we also have $I(0)=1$. Hence ${\cal W}\not\equiv 0$. Consequently,
${\cal W}$ strictly increases. In particular,

$\lim_{x\to\infty}{\cal W}(x)=L,$ (5.11)

for some $L\in(0,1]$. From (5.11), we deduce that there is a sequence
$\\{x_{n}\\}$, $x_{n}\to\infty$ as $n\to\infty$, such that ${\cal
W}^{\prime}(x_{n})\to 0$ as $n\to\infty$. Using this result as the boundary
condition and integrating the differential equation in (5.7), we have

$\frac{{\cal
W}^{\prime}(x)}{C^{2}(x)}=\int^{\infty}_{x}\frac{1}{C^{2}(y)}(1-{\cal
W}^{2}(y)){\cal W}(y)\,\mathrm{d}y.$ (5.12)

Thus, if $L<1$ in (5.11), we infer from (5.12) that ${\cal
W}^{\prime}(x)\to\infty$ as $x\to\infty$, which is false. Therefore $L=1$ in
(5.11). This establishes that ${\cal W}$ lies in the admissible space of the
problem (5.1). Thus $I({\cal W})\geq I_{0}$. Combining this result with (5.9),
we arrive at

$I({\cal W})=I_{0}.$ (5.13)

In other words, the existence of a solution to the minimization problem (5.1)
follows.

We next prove that, actually,

${\cal W}(x)=\tanh\left(\frac{x}{2}\right).$ (5.14)

Thus, in particular, equality in (5.10) holds. We achieve this goal by showing
that a positive solution to the initial value problem (5.7), with a limiting
value at infinity, is unique.

In fact, let ${\cal W}_{1}$ and ${\cal W}_{2}$ be two such solutions. Then
$w={\cal W}_{1}-{\cal W}_{2}$ satisfies

$\left(\frac{w^{\prime}}{C^{2}(x)}\right)^{\prime}+\frac{w}{C^{2}(x)}(1-{\cal
W}_{1}^{2}-{\cal W}_{2}^{2}-{\cal W}_{1}{\cal W}_{2})=0.$ (5.15)

On the other hand, for any differentiable function $h$, we have

$\displaystyle C^{2}(x)\left(\frac{{\cal W}_{1}{\cal
W}^{\prime}_{1}h^{2}}{C^{2}(x)}\right)^{\prime}+{{\cal
W}_{1}^{2}(h^{\prime})^{2}}+{({\cal W}_{1}{\cal W}_{2}+{\cal W}_{2}^{2}){\cal
W}_{1}^{2}h^{2}}$ $\displaystyle=([{\cal W}_{1}h]^{\prime})^{2}-(1-{\cal
W}_{1}^{2}-{\cal W}_{2}^{2}-{\cal W}_{1}{\cal W}_{2})({\cal W}_{1}h)^{2}.$
(5.16)

Now substituting $h=\frac{w}{{\cal W}_{1}}$ into (5), we have

$\displaystyle\frac{(w^{\prime})^{2}}{C^{2}}-\frac{1}{C^{2}}(1-{\cal
W}_{1}^{2}-{\cal W}_{2}^{2}-{\cal W}_{1}{\cal W}_{2})w^{2}$
$\displaystyle=\left(\frac{{\cal W}_{1}^{\prime}w^{2}}{C^{2}{\cal
W}_{1}}\right)^{\prime}+\frac{(w^{\prime}{\cal W}_{1}-w{\cal
W}_{1}^{\prime})^{2}}{C^{2}}+\frac{({\cal W}_{1}{\cal W}_{2}+{\cal
W}_{2}^{2})w^{2}}{C^{2}}.$ (5.17)

Moreover, using the Cauchy–Kovalevskaya theorem, we know that the solution of
(5.7) has the asymptotic form

${\cal W}(x)=bx+cx^{3}+\mbox{O}(x^{5}),\quad x\sim 0.$ (5.18)

As a consequence of (5.18), we may integrate (5) and drop the resulting
boundary terms to arrive at

$\displaystyle\int_{0}^{\infty}\left(\frac{(w^{\prime})^{2}}{C^{2}}-\frac{1}{C^{2}}(1-{\cal
W}_{1}^{2}-{\cal W}_{2}^{2}-{\cal W}_{1}{\cal
W}_{2})w^{2}\right)\,\mathrm{d}x$
$\displaystyle=\int_{0}^{\infty}\frac{1}{C^{2}(x)}\left({(w^{\prime}{\cal
W}_{1}-w{\cal W}_{1}^{\prime})^{2}}+{({\cal W}_{1}{\cal W}_{2}+{\cal
W}_{2}^{2})w^{2}}\right)\,\mathrm{d}x.$ (5.19)

However, by virtue of (5.15), we may rewrite the left-hand side of (5) as

$\displaystyle\int_{0}^{\infty}\left(\frac{(w^{\prime})^{2}}{C^{2}}-\frac{1}{C^{2}}(1-{\cal
W}_{1}^{2}-{\cal W}_{2}^{2}-{\cal W}_{1}{\cal
W}_{2})w^{2}\right)\,\mathrm{d}x$
$\displaystyle=\left(\frac{ww^{\prime}}{C^{2}}\right)_{0}^{\infty}-\int_{0}^{\infty}w\left(\left[\frac{w^{\prime}}{C^{2}}\right]^{\prime}+\frac{w}{C^{2}}(1-{\cal
W}_{1}^{2}-{\cal W}_{2}^{2}-{\cal W}_{1}{\cal W}_{2})\right)\,\mathrm{d}x=0.$
(5.20)

Combining (5) and (5), we conclude with $w\equiv 0$.

Hence the existence and uniqueness of an exterior energy-minimizing monopole
solution linking different vacuum states at $x\pm\infty$ and passes through
the unit-charge monopole state at $x=0$ is established which is the exact
solution $W_{1}$ stated in (3.1) and found in [3] which determines the minimum
energy, $E_{\min}=\frac{8\pi}{5}$, as stated in (3.13).

It is interesting to note that a rich family of even and odd solutions of
(2.17) linking the same and different vacuum states is obtained in [3] which
are characterized by the number $n$ of zeros they possess. When $n=1$, their
calculated energy coincides with our exact value; when $n\geq 2$, their
calculated energies are all significantly greater than that with $n=1$.

Motivated by the construction in this section and the even solutions with
multiple zeros linking to the same vacuum state, $W(-\infty)=W(\infty)=\pm 1$,
obtained in [3], we may impose the following minimization problem

$\sigma_{0}\equiv\inf\\{E(W)\,|\,W(-\infty)=W(\infty)=1,W(0)\leq 0\\},$ (5.21)

over the admissible set of functions which are absolutely continuous on all
compact subintervals of the real line. Thus, with the function (5.14) as a
testing function, i.e., choosing $W(x)=|\tanh(\frac{x}{2})|$ for all $x$, we
have the estimate

$\sigma_{0}\leq\frac{8\pi}{5}.$ (5.22)

With the problem (5.21) in mind, we consider the exact solution $W_{2}$ given
in (3.1) whose value at $x=0$ is

$W_{2}(0)=-h_{0},\quad h_{0}=\frac{1}{3+\sqrt{6}},$ (5.23)

and we ask what this ‘height’ $h_{0}$ means with regard to the problem (5.21).
So we take the trial function

$W_{h}(x)=\frac{2\cosh x-d}{2\cosh x+4+3\sqrt{6}},\quad d\equiv
2+3(2+\sqrt{6})h,\quad 0<h<1,$ (5.24)

which satisfies $W_{h}(0)=-h$ and is suggested by the form of $W_{2}(x)$ in
(3.1). Hence we can compute the energy $E(W_{h})$. Surprisingly, we find that
$E(W_{h})$ is maximized at $h=h_{0}$ as illustrated in Figure 2.

Figure 2: A plot of the normalized energy $\frac{E(W_{h})}{2\pi}$ against
$d=2+3(2+\sqrt{6})h$ for $\frac{1}{6}<h<\frac{1}{5}$ such that the critical
height $h_{0}=\frac{1}{3+\sqrt{6}}$ corresponds to
$d_{0}=\frac{12+5\sqrt{6}}{3+\sqrt{6}}\approx 4.45$ where the energy peaks
with the signature value, $0.9664$, first obtained in [3].

Therefore, the quantity $h_{0}$ may be regarded as a ‘canonical height’ at
which the energy confined over the profile family of functions of the form
(5.24) is maximized and brings forth a critical point of the energy
functional.

##  6 Remarks on oscillatory solutions

In [3], a family of solutions of (2.17) with prescribed numbers of zeros are
obtained. In this section, we revisit such a construction of oscillatory
solutions from the viewpoint of recasting the equation into a form of its
‘polar variable’ representation such that the appearance of multiple zeros of
the solutions becomes somewhat more transparent as we now show.

For our purpose, we use the following polar variable anzatz

$W=\rho\cos{\omega},\quad W^{\prime}=\rho\sin{\omega},$ (6.1)

to represent a solution of (2.17). Thus, we see that the pair
$(\rho,{\omega})$ satisfies the coupled system of equations

$\displaystyle\rho^{\prime}$ $\displaystyle=$
$\displaystyle\tanh\left(\frac{x}{2}\right)\,\rho\sin^{2}{\omega}+\rho^{3}\cos^{3}{\omega}\sin{\omega},$
(6.2) $\displaystyle{\omega}^{\prime}$ $\displaystyle=$
$\displaystyle-1+\tanh\left(\frac{x}{2}\right)\,\sin{\omega}\cos{\omega}+\rho^{2}\cos^{4}{\omega}.$
(6.3)

To justify this representation, we observe in view of the uniqueness theorem
for solutions of the initial value problems of ordinary differential equations
that if $\rho$ satisfying (6.2) is positive at a point then it will stay
positive where it exists. So, we may impose the initial condition

$\rho(0)=a>0,\quad{\omega}(0)=b.$ (6.4)

Thus, in particular, to produce an even solution of (2.17), we take $b=0$, and
an odd solution, we impose $b=\frac{\pi}{2}$. In the subsequent discussion, we
fix $b$ and adjust $a$ to achieve multiple zeros for the field $W$.
Accordingly, in view of (6.1), such a goal may be achieved by making the
angular variable ${\omega}$ circulate through the full circle multiple times,
in a sense to be made precise below.

Indeed, we fix $b,K_{0}>0$. Then the theorem of continuous dependence of
solutions on the initial data allows us to get some number $a_{0}>0$ such that
the solution to (6.2)–(6.4) satisfies

$0<\rho(x)\leq\frac{1}{2},\quad 0\leq x\leq K_{0},$ (6.5)

when $a\in(0,a_{0}]$. Applying (6.5) to (6.3) and using
$\sin{\omega}\cos{\omega}\leq\frac{1}{2}$, we get

${\omega}^{\prime}(x)\leq-\frac{1}{4},\quad 0\leq x\leq K_{0}.$ (6.6)

Therefore, for any integer $n\geq 1$, we may set $K_{0}=4n\pi$ (say). Then
(6.6) indicates that there are $n$ points $x_{1},\dots,x_{n}$ in $(0,K_{0}]$
such that

${\omega}(x_{i})=\mbox{an odd multiple of }\,\frac{\pi}{2},\quad i=1,\dots,n.$
(6.7)

As a consequence of this observation and (6.1), we see that if $b$ is not an
odd multiple of $\frac{\pi}{2}$, then $x_{1},\dots,x_{n}$ are the zeros of
$W$, and if $b$ is an odd multiple of $\frac{\pi}{2}$, then
$b,x_{1},\dots,x_{n}$ are the zeros of $W$. Hence a construction of a local
solution with a prescribed number of zeros is obtained. Note that such a
solution may be even, odd, or neither.

For the interior monopole equation (2.23), if we use the polar representatiion
(6.1) with the flippled variable $x=-y$, then we get the transformed governing
system

$\displaystyle\rho^{\prime}$ $\displaystyle=$
$\displaystyle\coth\left(\frac{y}{2}\right)\,\rho\sin^{2}{\omega}-\rho^{3}\cos^{3}{\omega}\sin{\omega},$
(6.8) $\displaystyle{\omega}^{\prime}$ $\displaystyle=$ $\displaystyle
1+\coth\left(\frac{y}{2}\right)\,\sin{\omega}\cos{\omega}-\rho^{2}\cos^{4}{\omega},$
(6.9)

where $0<y<\infty$ and $\rho^{\prime}=\frac{\mathrm{d}\rho}{\mathrm{d}y}$,
etc. Although $\coth(\frac{y}{2})$ is unbounded as $y\to 0$, corresponding to
the curvature singularity, we may again restrict our attention to a region
away from the singularity in order to get solutions with multiple zeros in
$W$. For example, it suffices to request $\coth(\frac{y}{2})<2$ for our method
here to work. Such a condition translates into

$y>\ln 3.$ (6.10)

In view of (2.18), this restriction gives us the interior region

$\frac{2M}{3}<r<M,$ (6.11)

to accommodate oscillatory monopole solutions similar to those in the exterior
region.

In order to obtain a global solution with a prescribed number of zeros with
desired asymptotic limits at infinity, $W(\pm\infty)=\pm 1$, additional work
is needed for the choice of the initial data. Here we omit the details since
much of the idea along the line is presented in [3] mathematically and
numerically. See also [15, 21, 31].

## References

  * [1] A. Actor, Classical solutions of $SU(2)$ Yang–Mills theories, Rev. Mod. Phys. 51 (1979) 461–525.
  * [2] R. Bartnik and J. McKinnon, Particlelike solutions of the Einstein–Yang–Mills equations Phys. Rev. Lett. 61 (1988) 141–144.
  * [3] P. Bizoń and M. Kahl, A Yang–Mills field on the extremal Reissner–Nordström black hole, Class. Quantum Grav.  33 (2016) 175013.
  * [4] M. Born and L. Infeld, Foundation of the new field theory, Nature 132 (1933) 1004.
  * [5] M. Born and L. Infeld, Foundation of the new field theory, Proc. Roy. Soc. (London) A 144 (1934) 425–451.
  * [6] L. Cao, S. Chen, and Y. Yang, Domain wall solitons arising in classical gauge field theories, Commun. Math. Phys. 369 (2019) 317–349.
  * [7] S. Chen and Y. Yang, Phase transition solutions in geometrically constrained magnetic domain wall models, J. Math. Phys. 51 (2010) 023504.
  * [8] S. Chen and Y. Yang, Domain wall equations, Hessian of superpotential, and Bogomol’nyi bounds, Nucl. Phys. B 904 (2016) 470–493.
  * [9] P. Curie, Sur la possibilité d’existence de la conductibilité magnétique et du magnétisme libre, Séances de la Société Française de Physique (Paris) (1894) 76–77.
  * [10] P. Dirac, Quantised singularities in the electromagnetic field, Proc. Roy. Soc. (London) A 133 (1931) 60–72.
  * [11] L. C. Evans, Partial Differential Equations, Graduate Studies in Mathematics 19, Amer. Math. Soc., Providence, Rhode Island, 2002.
  * [12] D. Gilbarg and N. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer, Berlin and New York, 1977.
  * [13] P. Goddard and D. I. Olive, Magnetic monopoles in gauge field theories, Rep. Prog. Phys. 41 (1978) 1357–1437.
  * [14] J. Greensite, An Introduction to the Confinement Problem, Lecture Notes in Physics 821, Springer-Verlag, Berlin and New York, 2011.
  * [15] S. P. Hastings, J. B. McLeod, and W. C. Troy, Static spherically symmetric solutions of a Yang–Mills field coupled to a dilaton, Proc. Roy. Soc. (London) A 449 (1995) 479–491.
  * [16] A. Jaffe and C. H. Taubes, Vortices and Monopoles, Birkhäuser, Boston, 1980.
  * [17] B. Julia and A. Zee, Poles with both magnetic and electric charges in non-Abelian gauge theory, Phys. Rev. D 11 (1975) 2227–2232.
  * [18] O. A. Ladyzhenskaya and N. N. Uraltseva, Linear and Quasilinear Elliptic Equations, Academic Press, New York and London, 1968.
  * [19] S. Mandelstam, Vortices and quark confinement in non-Abelian gauge theories, Phys. Lett. B 53 (1975) 476–478.
  * [20] S. Mandelstam, General introduction to confinement, Phys. Rep. C 67 (1980) 109–121.
  * [21] K. McLeod, W. C. Troy, and F. B. Weissler, Radial solutions of $\Delta u+f(u)=0$ with prescribed numbers of zeros, J. Diff. Eqs. 83 (1990) 368–378.
  * [22] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation, Freeman, New York, 1973.
  * [23] Y. Nambu, Strings, monopoles, and gauge fields, Phys. Rev. D 10 (1974) 4262–4268.
  * [24] A. M. Polyakov, Particle spectrum in the quantum field theory, JETP Lett. 20 (1974) 194–195.
  * [25] J. Preskill, Magnetic monopoles, Annu. Rev. Nucl. Part. Sci. 34 (1984) 461–530.
  * [26] R. Rajaraman, Solitons and Instantons, North Holland, Amsterdam, 1982.
  * [27] L. H. Ryder, Quantum Field Theory, 2nd ed., Cambridge U. Press, Cambridge, U. K., 1996.
  * [28] J. Schwinger, A magnetic model of matter, Science 165 (1969) 757–761.
  * [29] M. Shifman and A. Yung, Supersymmetric solitons and how they help us understand non-Abelian gauge theories, Rev. Mod. Phys. 79 (2007) 1139–1196.
  * [30] M. Shifman and A. Yung, Supersymmetric Solitons, Cambridge U. Press, Cambridge, U. K., 2009.
  * [31] J. Smoller, A. Wasserman, S. T. Yau, and J. B. McLeod, Smooth static solutions of the Einstein–Yang/Mills equations, Bull. Amer. Math. Soc. (N.S.) 27 (1992) 239–242; Commun. Math. Phys. 143 (1991) 115–147.
  * [32] G. ’t Hooft, Magnetic monopoles in unified gauge theories, Nucl. Phys. B 79 (1974) 276–284.
  * [33] G. ’t Hooft, On the phase transition towards permanent quark confinement, Nucl. Phys. B 138 (1978) 1–25.
  * [34] G. ’t Hooft, Topology of the gauge condition and new confinement phases in non-Abelian gauge theories, Nucl. Phys. B 190 (1981) 455–478.
  * [35] Yu. S. Tyupkin, V. A. Fateev, and A. S. Shvarts, Particle-like solutions of the equations of gauge theories, Theoret. Math. Phys. 26 (1976) 270–273.
  * [36] T. Vachaspati, Kinks and Domain Walls: An Introduction to Classical and Quantum Solitons, Cambridge U. Press, 2006, Cambridge, U. K.
  * [37] R. M. Wald, General Relativity, U. Chicago Press, Chicago and London, 1984.
  * [38] E. J. Weinberg, Classical Solutions in Quantum Field Theory: Solitons and Instantons in High Energy Physics, Cambridge U. Press, Cambridge, U. K., 2012.
  * [39] Y. Yang, Solitons in Field Theory and Nonlinear Analysis, Springer-Verlag, Berlin and New York, 2001.
  * [40] X. Zhang and Y. Yang, Solutions to the minimization problem arising in a dark monopole model in gauge field theory, Nucl. Phys. B 951 (2020) 114851.


# Automated all-functionals infrared and Raman spectra

Lorenzo Bastonero<EMAIL_ADDRESS>U Bremen Excellence Chair, Bremen
Center for Computational Materials Science, and MAPEX Center for Materials and
Processes, University of Bremen, D-28359 Bremen, Germany  Nicola Marzari  U
Bremen Excellence Chair, Bremen Center for Computational Materials Science,
and MAPEX Center for Materials and Processes, University of Bremen, D-28359
Bremen, Germany  Theory and Simulation of Materials (THEOS), and National
Centre for Computational Design and Discovery of Novel Materials (MARVEL),
École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Laboratory for Materials Simulations, Paul Scherrer Institut, 5232 Villigen
PSI, Switzerland

###### Abstract

Infrared and Raman spectroscopies are ubiquitous techniques employed in many
experimental laboratories, thanks to their fast and non-destructive nature
able to capture materials’ features as spectroscopic fingerprints.
Nevertheless, these measurements frequently need theoretical support in order
to unambiguously decipher and assign complex spectra. Linear-response theory
provides an effective way to obtain the higher-order derivatives needed, but
its applicability to modern exchange-correlation functionals remains limited.
Here, we devise an automated, open-source, user-friendly approach based on
ground-state density-functional theory and the electric enthalpy functional to
allow seamless calculations of first-principles infrared and Raman spectra. By
employing a finite-displacement and finite-field approach, we allow for the
use of any functional, as well as an efficient treatment of large low-symmetry
structures. Additionally, we propose a simple scheme for efficiently sampling
the Brillouin zone with different electric fields. To demonstrate the
capabilities of our approach, we provide illustrations using the ferroelectric
LiNbO3 crystal as a paradigmatic example. We predict infrared and Raman
spectra using various (semi)local, Hubbard corrected, and hybrid functionals.
Our results also show how PBE0 and extended Hubbard functionals yield in this
case the best match in term of peak positions and intensities, respectively.

## Introduction

Vibrational spectroscopies are one of the most powerful, fast and reliable
methods for materials’ structure identification and characterization. Among
the different techniques probing atomic motions, infrared (IR) absorption and
Raman experiments are widespread in most experimental set-ups, thanks to their
low cost and precision in identifying microscopic features. Despite these
advantages, data interpretation for many cases calls for theoretical support.
Density-functional theory (DFT) is most often employed as a reliable framework
for accurate predictions of such spectra, providing in principle exact
estimates of the second- and third-order derivatives of the energy functional
that are needed to predict vibrational spectra. Over the past decade, many
techniques have been developed to efficiently obtain non-resonant spectra of
insulating crystals, such as linear response theory [1, 2] for infrared
absorption, second-order response of the density matrix [3] and the $2n+1$
theorem [4] for Raman couplings. While being elegant and attractive, these
methods come at the cost of rather complex code implementations, limiting in
practice applications to selected semi-local functionals and often only to
norm-conserving pseudopotentials. This is especially the case for the third-
order derivatives of the total energy, required for the Raman calculations.
One common expedient for the latter is to resort to finite differences of the
dielectric tensor upon atomic displacements, computed using density-functional
perturbation theory (DFPT) [5, 6]. However, this approach scales in general as
$6N_{at}$ with the number of atoms $N_{at}$, a burden in low-symmetry systems,
and becomes out of reach for the non-linear optical susceptibility tensor,
which has been shown to be fundamental to achieve good agreement both in
single crystals [7] and powder spectra [8]. An alternative framework for
carrying out such derivatives is the finite numerical differentiation of
forces and polarization upon the application of a homogeneous electric
field[9, 10, 11, 4, 12]. By this technique, along with finite displacements
[13], one can compute straightforwardly all the quantities needed for IR and
Raman spectra, as well as allowing to use any functional (Hubbard corrected,
hybrids, …), and pseudopotential formalism (norm-conserving, ultrasoft, PAW)
with little effort. Furthermore, the use of finite fields overcomes the
expensive computation and poor scaling of the variational dielectric tensor
approach, allowing for the study of complex materials such as amorphous and
defective structures [12, 14]. In fact, within this scheme, one would need to
perform typically only 12 ground-state calculations in the presence of a
homogeneous electric field (see Eq. 13 and related discussion). While
appealing, this approach has seen scarce applicability, possibly due to the
many steps and pitfalls related to the underlying algorithm, which can be
highly non-trivial even for expert users.

A key solution is to streamline all of these operations through a modern
workflow manager, adhering in passing to the FAIR [15] principles for data
sharing, while being able of automatically submitting, parsing, and processing
the outputs generated by multiple calculations. In this study, we present a
comprehensive formulation for a unified finite-displacement and finite-field
approach, along with efficient Brillouin zone sampling, and encode this
formalism in a self-contained Python package that operates within the AiiDA
framework [16, 17], offering a scalable computational infrastructure for
automated and reproducible workflows and data provenance.

To showcase the capabilities of the approach, we demonstrate its effectiveness
in predicting the infrared and Raman spectra of LiNbO3 using seven different
functionals: LDA, PBE, PBEsol, PBEsol+U, PBEsol+U+V, PBE0, and HSE06.

## Results

### The finite-displacement and finite-field method

Infrared and Raman spectroscopy are investigation methods probing the atomic
vibrations of materials. In both techniques, a source of light is shone on the
sample, and the incident photons will be either absorbed or scattered by the
material. In IR experiments, peaks of the absorption spectra correspond to
collective atomic motions in resonance with that frequency, corresponding to
long-wavelength polar optical phonons in crystals. Raman spectroscopy records
instead a scattered frequency $\omega_{S}$; in this case, the excited lattice
vibration corresponds to the shift from the incoming laser frequency
$\omega_{L}$. Far from resonance the two phenomena are described as the change
in polarization $\mathbfcal{P}$ and the change in polarizability $\chi$, for
IR and Raman respectively, leading to different selection rules for the
activated phonon modes which ultimately provide a comprehensive fingerprint of
the material. The former mechanism is coupled to the Born effective charge
$\mathbf{Z}_{I}^{*}$ [18, 1] of an atom $I$ in the unit cell, and the latter
instead to the Raman tensor defined as $d\mathbf{\chi}/d\mathbf{\tau}_{I}$ in
the Placzek approximation [19, 20, 21], where $\mathbf{\chi}$ is the
electronic susceptibility tensor and $\tau_{I}$ the atomic displacement.

In terms of these quantities, the scattering intensities of a normal mode
$\nu$ are

$I_{\text{IR}}^{\nu}\propto\frac{|\bar{\mathbf{Z}}^{\nu}|^{2}}{\omega_{\nu}}$
(1)

for infrared spectra, and

$I_{\text{Raman}}^{\nu}\propto\frac{\omega_{S}^{4}}{\omega_{\nu}}(n_{\nu}(T)+1)|\mathbf{e}_{S}\cdot\mathbf{\alpha}^{\nu}\cdot\mathbf{e}_{L}|^{2}$
(2)

for Raman Stokes processes (i.e. when the scattered light frequency is lower
than the incident one). The tensors defining these amplitudes are defined via
the polarization vector

$\bar{Z}^{\nu}_{i}=\sum_{I,k}Z^{*}_{I,ik}\frac{e^{\nu}_{I,k}}{\sqrt{m_{I}}}~{},$
(3)

and the Raman susceptibility tensor

$\alpha^{\nu}_{ij}=\sqrt{\Omega}\sum_{I,k}\frac{\partial\chi^{(1)}_{ij}}{\partial\tau_{I,k}}\frac{e^{\nu}_{I,k}}{\sqrt{m_{I}}}~{},$
(4)

where $\omega_{\nu}$ and $\mathbf{e}^{\nu}$ are the frequency and the
eigenvector of the phonon $\nu$, $n_{\nu}(T)$ is the Bose-Einstein occupation
function, $m_{I}$ the mass of the atom $I$, and $\Omega$ the volume of the
unit cell In the harmonic approximation, the phonon normal modes at zone
center are obtained from the interatomic force constants (IFCs)

$\Phi_{IJ,ij}=\frac{\partial^{2}E}{\partial\tau_{I,i}\partial\tau_{J,j}}~{},$
(5)

where $E$ is the potential energy surface. The IFCs can be written in terms of
first-order derivative of atomic forces $\mathbf{F}_{I}$ with respect to
atomic displacements, and thus it can be obtained using a set of small
displacement configurations, as discussed in the literature [13, 1]. The
diagonalization of the mass-scaled IFCs, describing the harmonic ion dynamics,
gives the phonon frequencies and eigenvectors:

$\frac{1}{\sqrt{m_{I}m_{J}}}\Phi_{IJ,ij}e_{J,j}=\omega_{\nu}^{2}e_{I,i}~{}.$
(6)

While the IFCs extend over an arbitrary (but decaying) number of unit cells,
these photon spectroscopies utilize photons which carry a negligible momentum
$\mathbf{q}\approx 0$, thus only probing phonons in the vicinity of the zone
center, and justifying the evaluation via finite displacements in the
primitive cell only. The Born effective charges and Raman tensors can also be
defined as derivatives of forces with respect to a macroscopic electric field
$\mathbfcal{E}$ [12, 4]. For the former we have:

$Z^{*}_{I,ij}=\frac{\partial F_{I,j}}{\partial\mathcal{E}_{i}}~{}$ (7)

while for the latter, the Raman tensor reads:

$\frac{\partial\chi_{ij}}{\partial\tau_{I,k}}=\frac{1}{\Omega}\frac{\partial^{2}F_{I,k}}{\partial\mathcal{E}_{i}\partial\mathcal{E}_{j}}~{}.$
(8)

In a finite-difference scheme, this means that a small homogeneous electric
field must be applied in order to carry out the numerical differentiation. In
periodic-boundary conditions, this can be accomplished by extending the energy
functional $E[\psi]$ to include a term [22, 9, 10] accounting for the coupling
of the electric field and the total polarization of the system (often referred
as electric-enthalpy functional):

$\mathcal{F}[\psi]=E[\psi]-\Omega~{}\mathbfcal{E}\cdot[\mathbfcal{P}^{ion}+\mathbfcal{P}^{el}[\psi]]~{},$
(9)

where $\mathbfcal{P}^{ion}$ and $\mathbfcal{P}^{el}$ are the ionic and
electronic polarizations, respectively; the latter being described by the
modern theory of polarization [23, 24]. Such functional, while not bounded
from below, admits local minima representing self-consistent stationary
solutions [25, 26], if $\mathbfcal{E}$ is properly chosen; this guarantees the
simultaneous evaluation of forces and electronic polarization at convergence.
As a consequence, the electronic polarization (that comes for free) can be
exploited to extract its first and second-order derivatives with respect to an
electric field [4, 10]. These correspond respectively to the high-frequency
dielectric tensor $\mathbf{\epsilon}^{\infty}$, used to account for the non-
analytical contribution to the IFCs [2], and the non-linear optical
susceptibility $\mathbf{\chi}^{(2)}$, key in non-centrosymmetric crystals
where the Frölich contribution to the Raman tensor could be sizable [7, 8].
Their relationship with the electronic polarization derivatives is as follows:

$\displaystyle\epsilon^{\infty}_{ij}$
$\displaystyle=\delta_{ij}+\frac{4\pi}{\Omega}\frac{\partial\mathcal{P}^{el}_{i}}{\partial\mathcal{E}_{j}}~{},$
(10) $\displaystyle\chi^{(2)}_{ijk}$
$\displaystyle=\frac{4\pi}{\Omega}\frac{\partial^{2}\mathcal{P}^{el}_{i}}{\partial\mathcal{E}_{j}\partial\mathcal{E}_{k}}~{}.$
(11)

### Numerical differentiation

The possibility of consistently describing homogenoues electric fields allows
to access any derivative in atomic displacements and electric field, while
being easily extensible to advanced functionals and pseudopotential
formulations [12]. We calculate here the tensors for infrared and Raman cross-
sections using a central derivative formulation. These derivatives are
performed around $\mathbfcal{E}=\mathbf{0}$, which resembles the common
experimental setup where an external electromagnetic source is negligible,
either natural or artificial. Indicating (by dropping the indices) a computed
quantity as $A$, either a force or the electronic polarization, we calculate
its $m$-th derivatives with a numerical accuracy of $n$-th order via the
discretization formula:

$\frac{\partial^{m}A}{\partial\mathcal{E}_{i}^{m}}=\frac{1}{\Delta\mathcal{E}^{m}_{i}}\sum_{l=-n/2}^{n/2}~{}c_{l}^{m}~{}A\left(l\cdot\Delta\mathcal{E}_{i}\right)+\mathcal{O}(\Delta\mathcal{E}_{i}^{n}),$
(12)

where $\Delta\mathcal{E}_{i}>0$ is the small electric field, which we refer to
as the finite difference step, along the $i$-th Cartesian direction, and
$c_{l}^{m}$ the central derivative coefficients that can be automatically
found using the Fornberg’s algorithm[27]. The discretization is performed
using uniform steps such that $l=-\frac{n}{2},\dots,0,\dots,\frac{n}{2}$ is
always an integer number and the applied fields are evenly spaced (see Fig.
1(a)). For mixed second-order derivatives of $A$, i.e. for the Raman and non-
linear susceptibility tensors, we exploit a formula devised in [12] to reduce
the number of calculations. We assign
$\mathcal{E}_{i}=\mathcal{E}_{j}=\lambda$ and express the mixed derivative as:

$\frac{\partial
A}{\partial\mathcal{E}_{i}\partial\mathcal{E}_{j}}=\frac{1}{2}\left[\frac{\partial^{2}A}{\partial\lambda^{2}}-\frac{\partial^{2}A}{\partial\mathcal{E}_{i}^{2}}-\frac{\partial^{2}A}{\partial\mathcal{E}_{j}^{2}}\right]~{},$
(13)

where each second-order derivative is evaluated through Eq. 12; this can be
visualized schematically in Fig. 1(a). For example using a 2nd order numerical
accuracy, we would need only 12 self-consistent field (SCF) calculations with
a non-zero electric ($\Delta\mathcal{E}$-SCF) field to carry out all the
needed tensors both for IR and Raman spectra, even in low-symmetry systems
such as amorphous materials [12, 14]; this number can be reduced if one
exploits the symmetries of the material. In fact, there could exist symmetries
belonging to the point group of the crystal that transform an electric field
direction into an other. For example, if the crystal satisfy inversion
symmetry, this is sufficient to reduce by a half the number of
$\Delta\mathcal{E}$-SCF calculations needed, since forces and electronic
polarization produced by $\Delta\mathcal{E}_{i}$ are connected by a similarity
transformation to the one produced by $-\Delta\mathcal{E}_{i}$; keeping the
previous numerical example, we would need then just 6
$\Delta\mathcal{E}$-SCFs. For highly-symmetric crystals, such as cubic Si,
only two electric field directions must be evaluated; in this case, the number
of calculations reduces further to only 2. In general, the number of
$\Delta\mathcal{E}$-SCF calculations needed scales with the numerical accuracy
$n$ of the central formula as $\mathcal{N}_{n}=2\cdot\mathcal{N}_{n-2}$, where
$n$ is a positive multiple of 2.

While this approach is advantageous, we also highlight that the choice of the
field magnitude is restricted. As already stated before the electric enthalpy
functional 9 has no global minima, and long-lived meta-stable states exist
only below a certain critical field $\mathbfcal{E}^{c}$, where the Zener
tunneling is suppressed [22, 9]. An estimate of such critical value can be
performed using [9]:

$e|\mathbfcal{E}^{c}\cdot\mathbf{a}_{i}|\simeq E_{\text{gap}}/N_{i}$ (14)

where $e$ is the electric charge, $\mathbf{a}_{i}$ the lattice vectors
defining the unitcell, $E_{\text{gap}}$ the electronic band gap, and $N_{i}$
the number of k-points sampling the Brillouin zone in the $i$-th crystal
direction. For very small band gap semiconductors this critical field can be
very small, limiting in practice the convergence due to numerical noise. In
these circumstances, a balance must be found; nevertheless, for the majority
of cases this does not represent an actual limitation.

### Directional sampling

A subtle issue related to the Berry phase formalism used in the electric-
enthalpy formulation is the dense sampling of the Brillouin zone required to
achieve well converged values of the polarization, needed to carry out the
numerical derivatives. In fact, one can write the electronic polarization
using the string-averaged discretized Berry phase formulation [23] along the
direction of a primitive reciprocal lattice vector $\mathbf{b}_{i}$ as:

$\mathbfcal{P}_{el}\cdot\mathbf{b}_{i}=\frac{2e}{N^{(i)}_{\perp}\Omega}\sum_{l=1}^{N^{(i)}_{\perp}}\text{Im}\ln\prod_{j=0}^{N_{\parallel}^{(i)}-1}\det[S(\mathbf{k}_{j}^{(i)},\mathbf{k^{\prime}}_{j+1}^{(i)})]~{}.$
(15)

Here
$S_{nm}(\mathbf{k},\mathbf{k}^{\prime})=\braket{u_{n\mathbf{k}}}{u_{m\mathbf{k}^{\prime}}}$
represents the overlap matrix between Bloch’s states $u_{l\mathbf{k}}$
belonging to occupied bands $l$, and $N_{\perp}^{(i)}$ is the number of
strings (i.e. the number of k-points orthogonal to $\mathbf{b}_{i}$), each
having $N_{\parallel}^{(i)}$ k-points. Thus, when we apply a static electric
field along lattice vectors, an electronic polarization is induced according
to the field direction, properly sampled by increasing the number of k-points
within each string $N_{\parallel}^{(i)}$. While the derivatives can be
performed in the crystal reference system [28], we describe in the following
some possible implementations using a Cartesian reference system (the one
employed in the code). In this coordinates system one can follow two
approaches, namely either performing an orthogonal transformation of the
lattice to have the chosen field direction along one of the new basis vectors,
or generalise Eq. 15 to a generic reciprocal vector. It is clear it is more
convenient to choose the latter, as it requires only to modify Eq. 15 to
account for different $N_{\parallel}^{(i)}$ for each $N_{\perp}^{(i)}$, while
the former would require in general large supercells to accomodate the rotated
crystal. Nevertheless, many codes use automatically generated Monkhorst-Pack
(MP) meshes, and it is not intuitive how to properly choose a regular grid
that samples correctly an arbitrary direction. To automate the choice of such
grids, we propose here the concept of directional sampling targeting the
increase of the number of k-points in any desired direction and crystal
lattice, while still retaining a regular MP mesh. The main idea is to define
the BZ mesh by two quantities: an orthogonal and a parallel k-point distance,
$d_{\perp}$ and $d_{\parallel}$ respectively. The orthogonal distance
represents a minimum distance between k-points in the BZ that are commonly
used to achieve convergence of total energies and forces, thus usually
producing rather sparse meshes \- this quantity is then used as a reference
minimum for the generation of the grid. The parallel distance instead is used
to define a denser number of k-points in the desired direction, which in our
case is referring to the applied electric field. To do so, we transform the
electric field direction in crystal coordinates. This step allows us to assign
weights for the sampling in reciprocal space. Let
$\tilde{\mathbfcal{E}}=(\tilde{\mathcal{E}}_{1},\tilde{\mathcal{E}}_{2},\tilde{\mathcal{E}}_{3})$
be the electric field in the lattice reference system, we then define the
weights as:

$w_{i}=\frac{|\tilde{\mathcal{E}}_{i}|}{\max\\{|\tilde{\mathbfcal{E}}|\\}}~{},$
(16)

aiming at maximizing the sampling for lattice directions where the electric
dipole is expected to be more relevant. Finally, we can assign the number of
k-points $N_{i}$ of an MP mesh as:

$N_{i}=\max\left\\{w_{i}\frac{|\mathbf{b}_{i}|}{d_{\parallel}},\frac{|\mathbf{b}_{i}|}{d_{\perp}}\right\\}~{}.$
(17)

This methodology is easy to implement and generates meshes that better sample
the direction where the induced polarization is expected, significantly
reducing the number of k-points needed instead when using a constant distance
to define the entire MP grid.

We point out that the parallel distance defined may not resemble the true
distance among k-points within strings along a defined direction. This is
easily understood with a trivial example reported in Fig. 1, which represents
a simple two-dimensional BZ lattice sampling with an applied field along a
mixed direction. It is clear that the distance among adjacent k-points within
strings will have a greater distance $d^{true}_{\parallel}$ than the
$d_{\parallel}$ introduced due to geometrical reasons.

Figure 1:  (a) Schematic illustration of the central discretization formula 12
and 13 for derivatives with respect to the electric field in the directions
$i=\\{x,y,\lambda=(x,y)\\}$. The black dots represent the applied electric
fields $l\Delta\mathcal{E}_{i}$ used in the SCF calculations to evaluate the
electric-field-dependent forces and polarization to employ in 12 and 13, while
the different coloured circles are associated to the numerical accuracy $n$
achieved in the numerical differentiation (here shown up to $n=6$). It can be
noticed that the evaluation of the numerical derivative at $n>2$ provides
automatically derivatives at lower accuracy $n$; this is helpful to check the
convergence a posteriori in respect to the finite step size
$\Delta\mathcal{E}$. (b) Square reciprocal lattice showing a shifted regular
k-points sampling (black dots) of the BZ generated via Eqs. 16-17 for the
shown electric field $\tilde{\mathbfcal{E}}$ in reciprocal coordinates. The
strings of k-points are indicated by straight black lines, and are only shown
on part of the BZ for clarity. In this particular case
$\tilde{\mathcal{E}}_{1}=\tilde{\mathcal{E}}_{2}$, hence $w_{1}=w_{2}=1$ from
Eq. 16, meaning the generated mesh via Eq. 17 have k-points along $b_{1}$ and
$b_{2}$ equally spaced by $d_{\parallel}$, indicated on the figure by an
horizontal arrow. It can be seen that this introduced distance differs from
the actual distance $d^{true}_{\parallel}$ (diagonal arrow) among string-
connected k-points.

It is understood that our simplistic solution might be in the future
substituted with generalised MP grids [29] that can be ad hoc fine tuned for
effective and inexpensive samplings. However, this is far from being a trivial
task and thus out of scope for the present work.

### Computational workflow

In this section, we finally present the computational workflows that allow to
simulate the desired spectra, which can be infrared, Raman, or both, in an
automated fashion. We designed two main independent workflows as AiiDA
WorkChains [16]: one carrying out solely the finite displacement calculations
for the phonon modes to evaluate the IFCs, which we called PhononWorkChain,
and the other for the mixed total-energy derivatives via the homogeneous
electric fields, the DielectricWorkChain. When joined together, the above two
provide the IRamanSpectraWorkChain, as schematically outlined in Fig. 2.

Figure 2:  Computational workflow scheme of the IRamanSpectraWorckChain and of
its sub-workflows, the DielectricWorkChain and the PhononWorkChain. The inputs
can be easily, fully defined via the get_builder_from_protocol while
specifying only a StructureData and a Code. The black circles on top of boxes
mean that the calculation starts using charge density and wavefunctions from a
previous calculation.

The IRamanSpectraWorkChain takes as inputs a structure, a code and other
details on how to run the workflows (e.g. SCF parameters, wall-time,
parallelization options, and so on). These information are then split and/or
shared between the two sub-workflows. Both start running an initial SCF
ground-state calculation, meant to produce the charge density $\rho_{0}$ and
wavefunctions $\\{\psi_{i}\\}$ that will be used as a starting point for the
next runs. We generally represent this restart throughout the flowchart of
Fig. 2 by a curved white arrow in a black circle, and we use the dashed arrows
to indicate that the files have been copied over to a fresh new folder. When
these are finished, the PhononWorkChain generates the structures with
irreducible atomic displacements and computes their forces. The latter are
then used to produce the matrix of IFCs, which will be employed in a post-
processing step to obtain the phonon frequencies and eigenvectors. In
parallel, the DielectricWorkChain evaluates $\mathcal{E}^{c}$ through a non-
SCF calculation and extracts the numerical accuracy $n$ and the electric field
step $\Delta\mathcal{E}$, as follows: when the critical value is larger than
$10^{-4}~{}\text{(Ry)a.u.}$ ($Rydberg$ atomic units; 1 (Ry)a.u. $\approx
36.3609$ V/Å), a 4th order numerical expansion is chosen, otherwise a 2nd
order is used. The 4th order allows to achieve a greater numerical accuracy
for the finite differences, thus it is utilized whenever it is possible to
apply an electric field above threshold. In cases the electric field that can
be applied is smaller than $10^{-4}~{}\text{(Ry)a.u.}$, a 2nd order expansion
is sufficient, thus avoiding calculations that would be affected by too much
numerical noise [21]; the amplitude of the electric field step in fact may
result too small to produce significant digit changes in forces and
polarization, hence possibly worsening the numerical derivative even if adding
extra points to the formula. The step $\Delta\mathcal{E}$ is instead selected
according to:

$\Delta\mathcal{E}=\begin{cases}\frac{10^{-3}}{n}&\text{if
$\mathcal{E}^{c}>10^{-3}$},\\\ \frac{\text{round($\mathcal{E}^{c}$,
5)}}{n}&\text{if $10^{-4}<\mathcal{E}^{c}\leq 10^{-3}$},\\\
\frac{\text{round($\mathcal{E}^{c}$, 6)}}{n}&\text{if $\mathcal{E}^{c}\leq
10^{-4}$},\end{cases}$ (18)

where the electric field units are again expressed in (Ry)a.u. and round(x,y)
rounds the number $x$ to the $y$ digit after the decimal point. These choices
should guarantee larger values of the finite step to reduce the numerical
noise on forces and polarization, and at the same time higher accuracy to
remove the step size dependence of the finite difference formula, which can be
verified a posteriori when $n>2$ (see also Fig. 1(a)). We point out here that
eventually the user can still enforce the accuracy $n$ and/or the
$\Delta\mathcal{E}$ value that will be used in Eqs. 12 and 18. Once an
electric field is chosen and the ground-state SCF calculation is finished, a
symmetry analysis is performed to find the inequivalent directions
$\hat{\mathbfcal{E}}$ of the electric fields to apply, where
$\hat{\mathbfcal{E}}$ is a versor. For each such inequivalent direction an MP
mesh is determined through Eq. 17 and a series of SCF calculations are then
launched at $\mathbfcal{E}=\mathbf{0}$ with the corresponding mesh restarting
from the ground-state density of the previous SCF. The reason for computing
the SCF at $\mathbfcal{E}=\mathbf{0}$ with these meshes is twofold: first, it
is crucial to evaluate the polarization operator, set up in following finite
electric-field calculations at the corresponding mesh reference, and, second,
it guarantees the same level of convergence on forces and polarization for the
finite differentiation. In fact, $\hat{\mathbfcal{E}}$ and
$-\hat{\mathbfcal{E}}$ produce the same k-point meshes (as can be verified via
Eqs. 16-17), hence all the calculated values, at ${\mathbfcal{E}}=\mathbf{0}$
included, employed in the finite formula 12 have the same level of convergence
with respect to the particular k-point grid. Once finished, electric-enthalpy
calculations start for each inequivalent direction $\hat{\mathbfcal{E}}$ at
the smallest step value
$\Delta\mathbfcal{E}=\Delta\mathcal{E}\cdot\hat{\mathbfcal{E}}$. If the
numerical accuracy $n$ is greater than 2, further SCF calculations with
homogeneous electric fields are run, restarting from the previous steps till
the maximum value $(n/2)\Delta\mathbfcal{E}$; this loop corresponds to
increasing/decresaing the value of $l$ in Eq. 12. At this point, the computed
forces and polarizations are gathered together to perform the numerical
differentiation according to Eqs. 12 and 13 for achieving second and third
rank tensors 7-8-10. As a final step, the IRamanSpectraWorkChain joins all the
information in a unique AiiDA data type, named VibrationalData, able to
efficiently store large arrays and information in the AiiDA repository and
database, and with provenance graph for full reproducibility [16]. This data
type is a self-contained Python class, and it is provided with a number of
post-processing methods and features that allow the user to compute powder and
single crystal polarized vibrational intensities and frequencies, along with
the symmetry labels of the respective phonon modes. We also implemented an
IntensitiesAverageWorkChain which performs the spherical average spectra,
important in Raman for non-centrosymmetric crystals [8].

As mentioned earlier, the inputs of these workflows can be fully customized by
the user. We do not report here the full list of these inputs, which can be
explored through the usual command line interface of AiiDA [16] and through
the on-line documentation made available. However, we would like to emphasize
the possibility of fully defining the entire set of inputs via a limited set
of them, as already been done in other AiiDA packages [17]. In similar fashion
to Ref. [17], we identified the fundamental inputs to be the code, i.e. the
Quantum ESPRESSO [30, 31, 32] pw.x binary, associated with the computer, that
runs the calculations, the relaxed structure and the protocol, a string
defining in a general and user-friendly fashion the accuracy of the
simulation. This minimal set, used through the get_builder_from_protocol
method [17] designed for all the workchains, provides the user with a pre-
filled builder (the main unit in AiiDA used to run workflows) ready to be
submitted. This is adequate for non experts in the field that might be
interested in simulating their structure’s spectra with minimal knowledge, but
also in the context of high-throughput searches or for teaching or testing.

### All-functionals LiNbO3 spectra

We demonstrate the power of this approach studying the ferroelectric phase of
LiNbO3 (space group R3c, no. 161) for which the vibrational spectra have been
long debated in the literature [33, 34, 35, 36] and exhibit large LO-TO phonon
splitting at $\Gamma$ [33] that are ideal for exploring different functional
flavours. We selected three commonly used local and semi-local funcationals,
LDA, PBE and PBEsol, then Hubbard-corrected PBEsol in the standard and
extended version, namely DFT+U [37, 38] and DFT+U+V [39, 40, 41], and two of
the most used hybrids functionals for inorganic crystalline materials, HSE06
[42, 43] and PBE0 [44].

Through the use of the IRamanSpectraWorkChain, we carried out the simulation
of the vibrational spectra for all the functionals above and we compare our
results to experimental measurements from the most recent literature [8, 35,
36]. In particular, we compare our simulations in Fig. 3 to (a) different
Raman single-crystal polarization setups [35] and (b) two polarized infrared
absorption spectra [36], and in Fig. 4 to polycristalline Raman powder spectra
[8] for which we use the spherical average procedure accounting for the
$\chi^{(2)}$ contribution, as outlined in Ref. [8]. The single-crystal data
allowed us to compare directly the frequencies and, in the case of Raman, the
intensities. In (a) and (b), we ranked the spectra in the plots in function of
the statistical estimator
$\overline{|\Delta\omega|}=\sum_{\nu=1}^{N}|{\omega_{\nu}^{\text{exp}}-\omega_{\nu}^{\text{theo}}}|/N$,
where $N$ is the number of active modes detected.

Figure 3: Raman and IR spectra of ferroelectric LiNbO3 as obtained with seven
different functional, and compared to experimental measurements. (a) Raman
single-crystals spectra for different polarization setups (reported in
paranthesis) for pure TO modes corresponding to E and A1 symmetries. (b)
Infrared absorption using polarized light for E (upper panel) and A1 (lower
panel) symmetry modes. The theoretical spectra are ordered according to the
error $\overline{|\Delta\omega|}$ (shown in Tab. 1), from smallest (pink,
PBEsol+Usc+Vsc) to largest (blue, LDA). Hubbard-corrected PBEsol$+U_{sc}$ and
PBEsol$+U_{sc}+V_{sc}$ are labeled just $U_{sc}$ and $U_{sc}+V_{sc}$, the
subscript sc refers to the self-consistent calculation of $U$ and $V$ [45,
46]. Theoretical intensities are smeared with a constant 8 cm-1 broadened
Lorentzian.  Figure 4:  Raman powder spectra using the spherical average
formula [8] and accounting for the Frölich term to reproduce the
characteristic asymmetric line shapes. Theoretical intensities are smeared
with a constant 14 cm-1 broadened Lorentzian. The theoretical spectra are
reported in order of decreasing error $\overline{|\Delta\omega|}$, and
Hubbard-corrected PBEsol are abbreviated only with their corrective
parameters, as in Fig. 3.

All the functionals show good agreement with experiments and with the results
of the literature [33, 35]. PBEsol+Usc+Vsc, PBE, and the two hybrids show
competitive results among each other, while PBEsol, PBEsol+Usc and LDA have
quantitative larger discrepancies. In particular in Tab. 1 we report
statistical estimators for the comparison of the theoretical and experimental
phonon frequencies and intensities for Fig. 3(a).

| LDA | PBEsol+Usc | PBEsol | PBE | PBEsol+Usc+Vsc | HSE06 | PBE0
---|---|---|---|---|---|---|---
$\overline{|\Delta\omega|}$ | 21.3 | 16.5 | 14.4 | 10.8 | 9.7 | 9.3 | 8.9
$\overline{\Delta\omega}$ | 20.8 | 13.2 | 13.0 | 7.7 | 7.5 | 5.8 | 5.3
$|\Delta\omega|_{max}$ | 37.2 | 33.8 | 27.3 | 22.0 | 21.6 | 24.0 | 23.0
$\overline{|\Delta\omega|^{*}}$ | 13.8 | 13.0 | 8.8 | 9.2 | 8.7 | 9.8 | 9.6
$\overline{|\Delta I|}$ | 37 | 22 | 19 | 21 | 15 | 27 | 25
Table 1:  Statistical error estimators between functional predictions and
experiments, given in cm-1 for frequencies $\omega$ and in percentage for the
intensities $I$. The estimators are:
$\overline{|\Delta\omega|}=\sum_{\nu=1}^{N}|{\Delta\omega_{\nu}}|/N$ is the
mean absolute error;
$\overline{\Delta\omega}=\sum_{\nu=1}^{N}{\Delta\omega_{\nu}}/N$ is the mean
error; $|\Delta\omega_{max}|$ is the absolute max difference;
$\overline{|\Delta\omega|^{*}}=\sum_{\nu=1}^{N}|{\Delta\omega_{\nu}-\overline{\Delta\omega}}|/N$
is the absolute error for frequencies shifted by their mean difference;
$\overline{|\Delta I|}=100\cdot\min_{f}\\{\Delta I(f)\\}$ represents the
percentage error of the optimal scaled intensities; $N$ is the number of
considered modes and $\Delta\omega_{\nu}$ is the difference between
experimental and theoretical functional mode frequency.

These spectra setups allowed to accurately estimate the error of the
intensities evaluated through the minimization of the relative error:

$\Delta
I(f)=\sqrt{\frac{\sum_{\nu=1}^{N}({I_{\nu}^{\text{theo}}-fI_{\nu}^{\text{exp}}})^{2}}{\sum_{\nu=1}^{N}(fI_{\nu}^{\text{exp}})^{2}}}$
(19)

where $f$ represents the unknown scaling factor between theoretical $I^{theo}$
and experimental $I^{exp}$ intensities, due to the use of arbitrary units in
the experimental spectra. The smallest absolute error on frequencies
$\overline{|\Delta\omega|}$ is provided by the PBE0 functional with 8.9 cm-1,
and on intensities $\overline{|\Delta I|}$ by the PBEsol+Usc+Vsc functional
with 15%. It is to note that the latter also shows a remarkable
$\overline{|\Delta\omega|}$ of 9.7 cm-1, less than 1 cm-1 apart from the best
performing functional for this material, whereas the former has 25% error for
$\overline{|\Delta I|}$, 10% more with respect to PBEsol+Usc+Vsc. We also
reported the renormalized absolute error $\overline{|\Delta\omega|}^{*}$, for
which the frequencies are shifted by the bare mean error
$\overline{\Delta\omega}$ before the modulus average; this gives some insights
whether the discrepancy arises from a systematic error, e.g. due to the
geometry used. Interestingly, after renormalization, most of the functionals
get rather close error values. While this might be due to the renormalization
of the phonon frequencies due to the thermal expansion which tends to soften
the vibrational modes, all the theoretical spectra were performed at
constrained experimental lattice (see the Method section), meaning the
systematic error is unrelated to this effect. While other thermal
contributions should in principle be accounted for, their inclusion would
require more sophisticated level of theory [47, 48, 49] which goes beyond the
scope of the present work. Moreover, at the experimental conditions of our
references, such contributions are expected to produce only minor phonon
shifts. Given these points we might conclude in this case that such errors are
intrinsic to the functionals.

## Discussion

In this paper, we have presented a comprehensive and efficient approach that
utilizes finite displacements and finite fields to calculate infrared and
Raman spectra with any modern exchange-correlation functional and
pseudopotential formulation (norm-conserving, ultrasoft, PAW). Our approach
overcomes the current expensive alternative of the variational dielectric
tensor technique [6] for Raman coupling calculations, which remains
impractical in large or low symmetry systems. This is also very relevant when
generating coupling tensors in low-symmetric supercells to incorporate
temperature and quantum effects [48, 49, 50], or to train a tensorial machine-
learning potential useful for molecular dynamics simulations.

The formalism has been implemented as a highly optimized, automatic workflow
package for the prediction of vibrational spectra. The technical and practical
challenges have been successfully addressed with the help of the AiiDA
infrastructure [16], which empowered the design of fully reproducible,
reusable, user-friendly workflows for one of the most important class of
spectroscopic techniques for materials characterisation. As such, we believe
it will be broadly useful to both the computational community and the
experimental community at large.

Furthermore, we have shown that the workflows can effectively operate with
state-of-the-art functionals and can be employed for a comprehensive analysis
of both single crystals and powder spectra. In conclusion, we believe that the
present approach will pave the way for accelerated materials characterization
of materials displaying complex, challenging chemistry.

## Methods

The infrared and Raman calculations were carried out using the workflows
developed, relying on aiida-phonopy [51] for pre- and post-processing and
aiida-quantumespresso [16, 17], having Quantum ESPRESSO [30, 31, 32] as DFT
quantum engine. Before performing the vibrational spectra calculations, the
positions of the atoms in the LiNbO3 primitive cell are optimized at fixed
experimental lattice geometry [52] for all functionals until the total energy
and forces acting on atoms are less than $10^{-4}~{}\text{Ry/atom}$ and
$5\cdot 10^{-5}~{}\text{Ry/Bohr}$, respectively. Norm-conserving
pseudopotentials from the Pesudo-Dojo library [53] are used for the LDA and
hybrid functionals, employing PBE pseudopotentials for the latter, and
pseudopotentials from the precision SSSP (version 1.1) [54, 55, 56, 57, 58,
59, 53] library for the others. A cutoff of 90 Ry is used for the
wavefunctions; 4 and 12 times greater cutoffs are used for the charge density
(respectively for Pseudo-Dojo and SSSP based pseudopotentials). The Brillouin
zone is sampled uniformly with a 0.3 Å-1 k-points distance. The Hubbard
parameters are computed following the self-consistent procedure exploiting
density-functional perturbation theory, as explained in [45, 46, 60], using a
q-point distance of 0.4 Å-1, orthogonalized atomic orbitals [61, 46] for the
occupation matrices, and with onsite (intersite) parameters converged within
0.01 (0.005) eV (a very conservative threshold). The hybrid calculations were
performed using a 90 Ry cutoff and a single q-point mesh for the exact
exchange operator, and the adaptively compressed exchange (ACE) technique [62]
was used to speed up the convergence of the calculations. Finally, finite
differences were carried out using a 0.01 Å displacement distance and a denser
k-point distance of 0.15 Å-1 for the phonon calculations (in the
PhononWorkChain), and an electric field step of $5\cdot 10^{-4}$ (Ry)a.u., a
2nd-order numerical accuracy, and a k-point distance of 0.15 Å-1 in the
parallel direction of the applied electric field for the electric-field
derivatives (as carried out in the DielectricWorkChain).

## References

  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. D. Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Reviews of Modern Physics 73, 515 (2001).
  * 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, Physical Review B 55, 10355 (1997).
  * Lazzeri and Mauri [2003] M. Lazzeri and F. Mauri, First-principles calculation of vibrational Raman spectra in large systems: Signature of small rings in crystalline, Physical Review Letters 90, 036401 (2003).
  * Veithen  _et al._ [2005] M. Veithen, X. Gonze, and P. Ghosez, Nonlinear optical susceptibilities, Raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory, Physical Review B 71, 125107 (2005).
  * Baroni and Resta [1986] S. Baroni and R. Resta, Ab initio calculation of the low-frequency Raman cross section in silicon, Phys. Rev. B 33, 5969 (1986).
  * Liang  _et al._ [2019] Q. Liang, S. Dwaraknath, and K. A. Persson, High-throughput computation and evaluation of Raman spectra, Scientific Data 6, 10.1038/s41597-019-0138-y (2019).
  * Calzolari and Nardelli [2013] A. Calzolari and M. B. Nardelli, Dielectric properties and Raman spectra of ZnO from a first principles finite-differences/finite-fields approach, Scientific Reports 3, 10.1038/srep02999 (2013).
  * Popov  _et al._ [2020] M. N. Popov, J. Spitaler, V. K. Veerapandiyan, E. Bousquet, J. Hlinka, and M. Deluca, Raman spectra of fine-grained materials from first principles, npj Computational Materials 6, 10.1038/s41524-020-00395-3 (2020).
  * Souza  _et al._ [2002] I. Souza, J. Íñiguez, and D. Vanderbilt, First-principles approach to insulators in finite electric fields, Physical Review Letters 89, 117602 (2002).
  * Umari and Pasquarello [2002] P. Umari and A. Pasquarello, Ab initio molecular dynamics in a finite homogeneous electric field, Physical Review Letters 89, 157602 (2002).
  * Umari and Pasquarello [2003] P. Umari and A. Pasquarello, Polarizability and dielectric constant in density-functional supercell calculations with discretek-point samplings, Physical Review B 68, 085114 (2003).
  * Umari and Pasquarello [2005] P. Umari and A. Pasquarello, Density functional theory with finite electric field, International Journal of Quantum Chemistry 101, 666 (2005).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scripta Materialia 108, 1 (2015).
  * Niu  _et al._ [2008] L. Niu, J. Zhu, W. Gao, A. Liu, X. Han, and S. Du, First-principles calculation of vibrational Raman spectra of tetrahedral amorphous carbon, Physica B: Condensed Matter 403, 3559 (2008).
  * Wilkinson  _et al._ [2016] M. D. Wilkinson, M. Dumontier, I. J. Aalbersberg, G. Appleton, M. Axton, A. Baak, N. Blomberg, J.-W. Boiten, L. B. da Silva Santos, P. E. Bourne, J. Bouwman, A. J. Brookes, T. Clark, M. Crosas, I. Dillo, O. Dumon, S. Edmunds, C. T. Evelo, R. Finkers, A. Gonzalez-Beltran, A. J. Gray, P. Groth, C. Goble, J. S. Grethe, J. Heringa, P. A. ’t Hoen, R. Hooft, T. Kuhn, R. Kok, J. Kok, S. J. Lusher, M. E. Martone, A. Mons, A. L. Packer, B. Persson, P. Rocca-Serra, M. Roos, R. van Schaik, S.-A. Sansone, E. Schultes, T. Sengstag, T. Slater, G. Strawn, M. A. Swertz, M. Thompson, J. van der Lei, E. van Mulligen, J. Velterop, A. Waagmeester, P. Wittenburg, K. Wolstencroft, J. Zhao, and B. Mons, The FAIR guiding principles for scientific data management and stewardship, Scientific Data 3, 10.1038/sdata.2016.18 (2016).
  * Huber  _et al._ [2020] S. P. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. V. Yakutovich, C. W. Andersen, F. F. Ramirez, C. S. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari, B. Kozinsky, and G. Pizzi, AiiDA 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance, Scientific Data 7, 10.1038/s41597-020-00638-4 (2020).
  * Uhrin  _et al._ [2021] M. Uhrin, S. P. Huber, J. Yu, N. Marzari, and G. Pizzi, Workflows in AiiDA: Engineering a high-throughput, event-based engine for robust and modular computational workflows, Computational Materials Science 187, 110086 (2021).
  * Max Born [1998] K. H. Max Born,  _Dynamical Theory of Crystal Lattices_ (OUP Oxford, 1998).
  * Placzek [1934] G. Placzek, Rayleigh-streuung und Raman-effekt, Akademische Verlagsgesellschaft (1934).
  * Brüesch [1986] P. Brüesch,  _Phonons: Theory and Experiments II_ (Springer Berlin Heidelberg, 1986).
  * Porezag and Pederson [1996] D. Porezag and M. R. Pederson, Infrared intensities and Raman-scattering activities within density-functional theory, Physical Review B 54, 7830 (1996).
  * Nunes and Gonze [2001] R. W. Nunes and X. Gonze, Berry-phase treatment of the homogeneous electric field perturbation in insulators, Physical Review B 63, 155107 (2001).
  * Vanderbilt and King-Smith [1993] D. Vanderbilt and R. D. King-Smith, Electric polarization as a bulk quantity and its relation to surface charge, Physical Review B 48, 4442 (1993).
  * Resta [1994] R. Resta, Macroscopic polarization in crystalline dielectrics: the geometric phase approach, Reviews of Modern Physics 66, 899 (1994).
  * Souza  _et al._ [2004] I. Souza, J. Íñiguez, and D. Vanderbilt, Dynamics of Berry-phase polarization in time-dependent electric fields, Physical Review B 69, 085106 (2004).
  * Umari and Pasquarello [2004] P. Umari and A. Pasquarello, Finite electric field in density functional calculations with periodic boundary conditions, Computational Materials Science 30, 116 (2004).
  * Fornberg [1988] B. Fornberg, Generation of finite difference formulas on arbitrarily spaced grids, Mathematics of Computation 51, 699 (1988).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Maximally localized generalized wannier functions for composite energy bands, Physical Review B 56, 12847 (1997).
  * Wang  _et al._ [2021] Y. Wang, P. Wisesa, A. Balasubramanian, S. Dwaraknath, and T. Mueller, Rapid generation of optimal generalized monkhorst-pack grids, Computational Materials Science 187, 110100 (2021).
  * 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, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with quantum ESPRESSO, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Giannozzi  _et al._ [2020] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. F. Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, Quantum ESPRESSO toward the exascale, The Journal of Chemical Physics 152, 154105 (2020).
  * Veithen and Ghosez [2002] M. Veithen and P. Ghosez, First-principles study of the dielectric and dynamical properties of lithium niobate, Physical Review B 65, 214302 (2002).
  * Margueron  _et al._ [2012] S. Margueron, A. Bartasyte, A. M. Glazer, E. Simon, J. Hlinka, I. Gregora, and J. Gleize, Resolved E-symmetry zone-centre phonons in LiTaO3 and LiNbO3, Journal of Applied Physics 111, 104105 (2012).
  * Sanna  _et al._ [2015] S. Sanna, S. Neufeld, M. Rüsing, G. Berth, A. Zrenner, and W. G. Schmidt, Raman scattering efficiency in LiTaO3 and LiNbO3 crystals, Physical Review B 91, 224302 (2015).
  * Kojima  _et al._ [2016] S. Kojima, K. Kanehara, T. Hoshina, and T. Tsurumi, Optical phonons and polariton dispersions of congruent LiNbO3 studied by far-infrared spectroscopic ellipsometry and Raman scattering, Japanese Journal of Applied Physics 55, 10TC02 (2016).
  * Himmetoglu  _et al._ [2011] B. Himmetoglu, R. M. Wentzcovitch, and M. Cococcioni, First-principles study of electronic and structural properties of CuO, Physical Review B 84, 115108 (2011).
  * Himmetoglu  _et al._ [2013] B. Himmetoglu, A. Floris, S. de Gironcoli, and M. Cococcioni, Hubbard-corrected DFT energy functionals: The LDA+U description of correlated systems, International Journal of Quantum Chemistry 114, 14 (2013).
  * Jr and Cococcioni [2010] V. L. C. Jr and M. Cococcioni, Extended DFT+U+V method with on-site and inter-site electronic interactions, Journal of Physics: Condensed Matter 22, 055602 (2010).
  * Cococcioni and Marzari [2019] M. Cococcioni and N. Marzari, Energetics and cathode voltages of LiMPO4 olivines ( M=Fe , Mn) from extended Hubbard functionals, Physical Review Materials 3, 033801 (2019).
  * Timrov  _et al._ [2020] I. Timrov, F. Aquilante, L. Binci, M. Cococcioni, and N. Marzari, Pulay forces in density-functional theory with extended Hubbard functionals: From nonorthogonalized to orthogonalized manifolds, Physical Review B 102, 235159 (2020).
  * Heyd  _et al._ [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened coulomb potential, The Journal of Chemical Physics 118, 8207 (2003).
  * Heyd  _et al._ [2006] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Erratum: “Hybrid functionals based on a screened coulomb potential” [J. Chem. Phys. 118, 8207 (2003)], The Journal of Chemical Physics 124, 219906 (2006).
  * Perdew  _et al._ [1996] J. P. Perdew, M. Ernzerhof, and K. Burke, Rationale for mixing exact exchange with density functional approximations, The Journal of Chemical Physics 105, 9982 (1996).
  * Timrov  _et al._ [2018] I. Timrov, N. Marzari, and M. Cococcioni, Hubbard parameters from density-functional perturbation theory, Physical Review B 98, 085127 (2018).
  * Timrov  _et al._ [2021] I. Timrov, N. Marzari, and M. Cococcioni, Self-consistent Hubbard parameters from density-functional perturbation theory in the ultrasoft and projector-augmented wave formulations, Physical Review B 103, 045141 (2021).
  * Monacelli  _et al._ [2018] L. Monacelli, I. Errea, M. Calandra, and F. Mauri, Pressure and stress tensor of complex anharmonic crystals within the stochastic self-consistent harmonic approximation, Physical Review B 98, 024106 (2018).
  * Monacelli and Mauri [2021] L. Monacelli and F. Mauri, Time-dependent self-consistent harmonic approximation: anharmonic nuclear quantum dynamics and time correlation functions, Physical Review B 103, 104305 (2021).
  * Monacelli  _et al._ [2021] L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea, and F. Mauri, The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects, Journal of Physics: Condensed Matter 33, 363001 (2021).
  * Siciliano  _et al._ [2023] A. Siciliano, L. Monacelli, G. Caldarelli, and F. Mauri, Wigner Gaussian dynamics: simulating the anharmonic and quantum ionic motion, Physical Review B 107, 174307 (2023).
  * Bastonero [2022] L. Bastonero, AiiDA-Phonopy code v1.0.0, https://github.com/aiida-phonopy/aiida-phonopy (2022).
  * Weigel  _et al._ [2020] T. Weigel, C. Funke, M. Zschornak, T. Behm, H. Stöcker, T. Leisegang, and D. C. Meyer, X-ray diffraction using focused-ion-beam-prepared single crystals, Journal of Applied Crystallography 53, 614 (2020).
  * 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, Computer Physics Communications 226, 39 (2018).
  * Prandini  _et al._ [2018] G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, Precision and efficiency in solid-state pseudopotential calculations, npj Computational Materials 4, 10.1038/s41524-018-0127-2 (2018).
  * Garrity  _et al._ [2014] K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Pseudopotentials for high-throughput DFT calculations, Computational Materials Science 81, 446 (2014).
  * Schlipf and Gygi [2015] M. Schlipf and F. Gygi, Optimization algorithm for the generation of ONCV pseudopotentials, Computer Physics Communications 196, 36 (2015).
  * Willand  _et al._ [2013] A. Willand, Y. O. Kvashnin, L. Genovese, Á. Vázquez-Mayagoitia, A. K. Deb, A. Sadeghi, T. Deutsch, and S. Goedecker, Norm-conserving pseudopotentials with chemical accuracy compared to all-electron calculations, The Journal of Chemical Physics 138, 104109 (2013).
  * Corso [2014] A. D. Corso, Pseudopotentials periodic table: From H to Pu, Computational Materials Science 95, 337 (2014).
  * Topsakal and Wentzcovitch [2014] M. Topsakal and R. Wentzcovitch, Accurate projected augmented wave (PAW) datasets for rare-earth elements (RE=La-Lu), Computational Materials Science 95, 263 (2014).
  * Timrov  _et al._ [2022] I. Timrov, N. Marzari, and M. Cococcioni, HP - A code for the calculation of Hubbard parameters using density-functional perturbation theory, Computer Physics Communications 279, 108455 (2022).
  * Floris  _et al._ [2020] A. Floris, I. Timrov, B. Himmetoglu, N. Marzari, S. de Gironcoli, and M. Cococcioni, Hubbard-corrected density functional perturbation theory with ultrasoft pseudopotentials, Physical Review B 101, 064305 (2020).
  * Lin [2016] L. Lin, Adaptively compressed exchange operator, Journal of Chemical Theory and Computation 12, 2242 (2016).
  * Bastonero and Marzari [2023] L. Bastonero and N. Marzari, Automated all-functionals infrared and raman spectra, Materials Cloud Archive 2023.118, 10.24435/materialscloud:90-36 (2023).

## Data availability

The data produced in this work can be found on the Materials Cloud
Archive[63].

## Code availability

The source code is made available open-source on GitHub
(https://github.com/bastonero/aiida-vibroscopy). It is also distributed as an
installable package through the Python Package Index
(https://pypi.org/project/aiida-vibroscopy/).

## Acknowledgements

L.B. would like to thank (in alphabetic order) Marnik Bercx, Chiara
Cignarella, Virginie de Mestral, Michele Kotiuga, Eric Macke, Giovanni Pizzi,
Norma Rivano and Iurii Timrov for fruitful discussions. The authors gratefully
acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) under
Germany’s Excellence Strategy (EXC 2077, No. 390741603, University Allowance,
University of Bremen) and Lucio Colombi Ciacchi, the host of the “U Bremen
Excellence Chair Program,” as well as computing time granted by the Resource
Allocation Board and provided on the supercomputer Lise and Emmy at NHR@ZIB
and NHR@Göttingen as part of the NHR infrastructure. The calculations for this
research were conducted with computing resources under the project hbc00053
and hbi00059.

## Author contributions

L.B. designed and developed the code, wrote the first draft of the manuscript,
ran and analysed the calculations. N.M. supervised the project, discussed the
results, and contributed to the final version of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Correspondence should be addressed to L.B.