# Remarks on the linear wave equation

John M. Ball  Heriot-Watt University and Maxwell Institute for Mathematical
Sciences, Edinburgh  Hong Kong Institute for Advanced Study

###### Abstract.

We make some remarks on the linear wave equation concerning the existence and
uniqueness of weak solutions, satisfaction of the energy equation, growth
properties of solutions, the passage from bounded to unbounded domains, and
reconciliation of different representations of solutions.

Dedicated to Bob Pego in admiration.

##  1\. Introduction

Let $\Omega\subset\mathbb{R}^{n}$ be open with boundary $\partial\Omega$. In
this paper we consider the linear wave equation

(1.1) $\displaystyle u_{tt}=\Delta u,$

for $u=u(x,t)$, $x\in\Omega$, $t\geqslant 0$, with boundary condition

(1.2) $\displaystyle\left.u\right|_{\partial\Omega}=0,$

(interpreted appropriately if $\Omega$ is unbounded) and initial conditions

(1.3) $\displaystyle u(x,0)=u_{0}(x),\,u_{t}(x,0)=v_{0}(x),$

where $u_{0}\in H_{0}^{1}(\Omega),\,v_{0}\in L^{2}(\Omega)$. Our aim is to
make some remarks which are hard to find in the literature, although largely
implicit in it, concerning the existence and uniqueness of weak solutions to
(1.1)-(1.3), satisfaction of the energy equation

(1.4) $\displaystyle\frac{1}{2}\int_{\Omega}(|\nabla
u(x,t)|^{2}+u_{t}(x,t)^{2})\,dx=\frac{1}{2}\int_{\Omega}(|\nabla
u_{0}(x)|^{2}+v_{0}(x)^{2})\,dx,\;t\geqslant 0,$

growth properties of solutions, the passage from bounded to unbounded domains,
and the reconciliation of different representations of solutions. Although
these remarks are perhaps original only in various details, it is hoped that
readers may find their combination useful.

If $\Omega$ is bounded, perhaps the easiest method for proving existence and
uniqueness is to represent $u$ as a Fourier expansion in the eigenfunctions
$\omega_{j}\in H_{0}^{1}(\Omega)$ of $-\Delta$, with corresponding real
eigenvalues $\lambda_{j}>0$, in terms of which

(1.5) $\displaystyle
u(x,t)=\sum_{j=1}^{\infty}\left(u_{0j}\cos(\sqrt{\lambda_{j}}t)+v_{0j}\frac{\sin(\sqrt{\lambda_{j}}t)}{\sqrt{\lambda_{j}}}\right)\omega_{j},$

where $u_{0j}=(u_{0},\omega_{j}),v_{0j}=(v_{0},\omega_{j})$ and
$(\cdot,\cdot)$ denotes the inner product in $L^{2}(\Omega)$.

In the case $\Omega=\mathbb{R}^{n}$ the solution can be given explicitly in
terms of the Fourier transforms $\hat{u}_{0},\,\hat{v}_{0}$ of the initial
data as

(1.6) $\displaystyle
u(x,t)=\frac{1}{(2\pi)^{n/2}}\int_{\mathbb{R}^{n}}\left(\hat{u}_{0}(\xi)\cos(|\xi|t)+\hat{v}_{0}(\xi)\frac{\sin(|\xi|t)}{|\xi|}\right)e^{ix\cdot\xi}\,d\xi.$

Alternatively one can use Poisson’s method of spherical means (see, for
example, [19]), which in the case $n=3$ leads to Kirchhoff’s solution

(1.7) $\displaystyle u(x,t)=\frac{1}{4\pi
t}\int_{S(x,t)}v_{0}(y)\,dS_{y}+\frac{\partial}{\partial t}\left(\frac{1}{4\pi
t}\int_{S(x,t)}u_{0}(y)\,dS_{y}\right),$

where $S(x,t)=\\{y\in\mathbb{R}^{3}:|y-x|=t\\}$ and $S_{y}$ is the usual
$(n-1)$-dimensional surface measure. Yet another method is to represent the
solution as a superposition of plane waves via the Radon transform [12, 21].

For $\Omega$ a general (possibly unbounded) open set these methods do not
apply, and a natural approach that we review in Section 2 is via the Hille-
Yosida theorem. Among various remarks we comment on the growth in time of the
$L^{2}$ norm of solutions when $\Omega$ is unbounded. In order to make the
connection with weak solutions in the sense of distributions, and thus to
establish uniqueness of weak solutions, it is convenient to calculate the
adjoint of the wave operator, which we do in Section 5. We show that the
adjoint is injective and that there are no nontrivial linear constants of
motion. The use of the phase space $H=H_{0}^{1}(\Omega)\times L^{2}(\Omega)$
for unbounded domains $\Omega$ is motivated by taking a limit of the semiflows
corresponding to an increasing sequence of bounded open sets $\Omega_{j}$ with
union $\Omega$ (Theorem 9).

It is not completely obvious how to reconcile the different representations of
solutions described above. We explore this for Kirchhoff’s solution in Section
5, showing that the derivation of Kirchhoff’s solution implies smoothing
properties under the taking of averages over spheres and balls that are
familiar in harmonic analysis. We illustrate the related harmonic analysis
methods by showing (Theorem 10) that if $B=B(0,1)$ denotes the unit ball in
$\mathbb{R}^{n}$ then for $f\in L^{2}(\mathbb{R}^{n})$ and $t>0$ the average

(1.8) $\displaystyle{\mathcal{N}}_{t}(f)(x)=\int_{B}f(x+tz)\,dz$

belongs to $H^{\frac{n+1}{2}}(\mathbb{R}^{n})$.

##  2\. Existence and uniqueness via the Hille-Yosida theorem

Let $H=H_{0}^{1}(\Omega)\times L^{2}(\Omega)$. $H$ is a Hilbert space with
inner product

$\langle\left(\begin{array}[]{ll}u\\\
v\end{array}\right),\left(\begin{array}[]{ll}p\\\
q\end{array}\right)\rangle=\int_{\Omega}(up+\nabla u\cdot\nabla p+vq)\,dx.$

We briefly review how, for a general (possibly unbounded) domain
$\Omega\subset\mathbb{R}^{n}$, it can be proved using the Hille-Yosida theorem
that (1.1)-(1.3) generates a semiflow on $H$, and that the energy equation
(1.4) holds. For equivalent treatments see [1, Chapter 7], [8, Chapter IV],
both using cosine families, and for the case of bounded domains [7, p444]. We
write (1.1) in the form

(2.1) $\displaystyle\dot{w}=Aw,$

where $v=u_{t}$, $w=\left(\begin{array}[]{ll}u\\\ v\end{array}\right)$ and

(2.4) $\displaystyle A=\left(\begin{array}[]{ll}0&{\mathbf{1}}\\\
\Delta&0\end{array}\right).$

We regard $A$ as an unbounded linear operator on $H$ with domain

(2.7) $\displaystyle D(A)=\\{\left(\begin{array}[]{l}u\\\
v\end{array}\right)\in H:\Delta u\in L^{2}(\Omega),v\in H_{0}^{1}(\Omega)\\}.$

Then $D(A)$ is dense in $H$ and it is easily checked that $A:D(A)\to H$ is
closed. We apply the following special case of the Hille-Yosida theorem (see,
for example, [23, Corollary 3.8], [7, p441]):

######  Theorem 1.

A closed densely defined linear operator $A$ on a real Banach space $X$ is the
generator of a $C^{0}$-semigroup $\\{T(t)\\}_{t\geqslant 0}$ satisfying for
some $\omega\in\mathbb{R}$

$\|T(t)\|\leqslant e^{\omega t},\;t\geqslant 0,$

if and only if $(\omega,\infty)\subset\rho(A)$ and
$\|R_{\lambda}\|\leqslant\frac{1}{\lambda-\omega}\text{ for }\lambda>\omega.$

Here a $C^{0}$-semigroup $\\{T(t)\\}_{t\geqslant 0}$ is a family of bounded
linear operators $T(t):X\to X$ satisfying (i) $T(0)=$identity, (ii)
$T(s+t)=T(s)T(t)$ for all $s,t\geqslant 0$ and (iii) $t\mapsto T(t)p$ is
continuous from $[0,\infty)\to X$ for all $p\in X$. The resolvent set
$\rho(A)$ is the set of $\lambda$ such that $\lambda{\mathbf{1}}-A:D(A)\to X$
is one-to-one and onto, and for $\lambda\in\rho(A)$ the resolvent operator
$R_{\lambda}:X\to X$ is defined by
$R_{\lambda}w=(\lambda{\mathbf{1}}-A)^{-1}w.$

######  Theorem 2.

Let $A$ be given by (2.4),(2.7). Then $A$ is the generator of a
$C^{0}$-semigroup $\\{T(t)\\}_{t\geqslant 0}$ on $H$, and the energy equation

(2.8) $\displaystyle E(T(t)p)=E(p),\;\;t\geqslant 0$

is satisfied for all $p=\left(\begin{array}[]{l}u_{0}\\\
v_{0}\end{array}\right)\in H$, where

(2.9) $\displaystyle E(w):=\frac{1}{2}\int_{\Omega}(|\nabla
u(x)|^{2}+v(x)^{2})\,dx.$

###### Proof.

We first show that $(0,\infty)\subset\rho(A)$. Thus we need to prove that for
any $\lambda>0$ and $f\in H_{0}^{1}(\Omega)$, $g\in L^{2}(\Omega)$ there
exists a unique solution $\left(\begin{array}[]{l}u\\\ v\end{array}\right)\in
H$ with $\Delta u\in L^{2}(\Omega)$, $v\in H_{0}^{1}(\Omega)$ to

(2.10) $\displaystyle\lambda u-v$ $\displaystyle=f,$ (2.11)
$\displaystyle\lambda v-\Delta u$ $\displaystyle=g.$

Since $v=\lambda u-f$ we just need to show that there is a unique solution
$u\in H_{0}^{1}(\Omega)$ to

(2.12) $\displaystyle-\Delta u+\lambda^{2}u-(\lambda f+g)=0.$

The existence follows by minimization of the functional

(2.13) $\displaystyle I(u)=\int_{\Omega}\left(\frac{1}{2}|\nabla
u|^{2}+\frac{\lambda^{2}}{2}u^{2}-(\lambda f+g)u\right)\,dx$

over $H_{0}^{1}(\Omega)$ via the direct method of the calculus of variations
and showing that the minimizer satisfies the Euler-Lagrange equation, and the
uniqueness follows since the difference $z$ of two solutions to (2.12)
satisfies $\int_{\Omega}(z^{2}+|\nabla z|^{2})\,dx=0$. To prove the resolvent
estimate, note that $R_{\lambda}\left(\begin{array}[]{l}f\\\
g\end{array}\right)=\left(\begin{array}[]{l}u\\\ v\end{array}\right)$, and
thus taking the inner product in $H$ of (2.10), (2.11) with
$\left(\begin{array}[]{c}u\\\ v\end{array}\right)$ we obtain

(2.14) $\displaystyle\lambda\int_{\Omega}(u^{2}+|\nabla
u|^{2}+v^{2})\,dx=(v,u)+\int_{\Omega}(fu+\nabla f\cdot\nabla u+gv)\,dx,$

where $(\cdot,\cdot)$ denotes the inner product in $L^{2}(\Omega)$. But
$\displaystyle(v,u)\leqslant\frac{1}{2}\int_{\Omega}(u^{2}+|\nabla
u|^{2}+v^{2})\,dx$, and hence

(2.15) $\displaystyle(\lambda-\frac{1}{2})\int_{\Omega}(u^{2}+|\nabla
u|^{2}+v^{2})\,dx\leqslant\int_{\Omega}(f,\nabla f,g)\cdot(u,\nabla u,v)\,dx,$

from which the estimate

(2.16)
$\displaystyle\|R_{\lambda}\|\leqslant\frac{1}{\lambda-\frac{1}{2}}\text{ for
}\lambda>\frac{1}{2}$

follows.

It remains to prove that the energy equation (2.8) holds, this not being
immediately obvious since the formal derivation of it via multiplication of
(1.1) by $u_{t}$ and integration is not valid for $u_{0}\in
H_{0}^{1}(\Omega),v_{0}\in L^{2}(\Omega)$. To this end we note that

(2.17) $\displaystyle E(w):=\frac{1}{2}\int_{\Omega}(|\nabla
u(x)|^{2}+v(x)^{2})\,dx$

is a $C^{1}$ function of $w=\left(\begin{array}[]{c}u\\\
v\end{array}\right)\in H$, and that if $w\in D(A)$ then

(2.18) $\displaystyle E^{\prime}(w)(Aw)=\int_{\Omega}(\nabla u\cdot\nabla
v+v\cdot\Delta u)\,dx=0.$

But by a well-known result for linear semigroups (see, for example, [23,
Theorem 2.4]), the map $t\mapsto T(t)p$ is $C^{1}$ for $p\in D(A)$ with
derivative $AT(t)p$. Thus if $p\in D(A)$ then $t\mapsto E(T(t)p)$ is $C^{1}$
with derivative $E^{\prime}(T(t)p)(AT(t)p)=0$. Hence $E(T(t)p)=E(p)$ for all
$t\geqslant 0$, $p\in D(A)$, and thus, since $D(A)$ is dense in $H$ and $E$
and $T(t)$ are continuous, also for $p\in X$ . ∎

######  Remark 1.

An alternative approach to proving existence, involving much the same
calculations, is to apply the Lumer-Phillips theorem (see e.g. [23, Theorem
4.3]) by showing that the operator $-A+\lambda{\mathbf{1}}$ is maximal
monotone for $\lambda>\frac{1}{2}$.

######  Remark 2.

The fact that the time reversibility of (1.1) implies that $A$ generates a
group $\\{T(t)\\}_{t\in\mathbb{R}}$ of bounded linear operators can be proved
either by checking that $-A$ generates a $C^{0}$-semigroup or by verifying
that

(2.23) $\displaystyle T(-t):=\left(\begin{array}[]{cc}{\mathbf{1}}&0\\\
0&-{\mathbf{1}}\end{array}\right)T(t)\left(\begin{array}[]{cc}{\mathbf{1}}&0\\\
0&-{\mathbf{1}}\end{array}\right)$

satisfies $T(-t)T(t)={\mathbf{1}}$ by calculating $\frac{d}{dt}(T(-t)T(t)p)=0$
for $p\in D(A)$.

######  Remark 3.

Using [23, Theorem 2.4] we also have the regularity result (see [1, Theorem
7.2]) that if $u_{0}\in H_{0}^{1}(\Omega)$, $\Delta u_{0}\in L^{2}(\Omega)$
and $u_{1}\in H_{0}^{1}(\Omega)$ then $u_{tt}=\Delta u\in
C([0,\infty);L^{2}(\Omega))$, $u_{t}\in C([0,\infty);H_{0}^{1}(\Omega))$.

######  Remark 4.

Note that if $\left(\begin{array}[]{c}u_{0}\\\ v_{0}\end{array}\right)\in H$
then, since $\left(\begin{array}[]{c}u_{0}\\\
0\end{array}\right)=A\left(\begin{array}[]{c}0\\\ u_{0}\end{array}\right)$,

(2.30) $\displaystyle T(t)\left(\begin{array}[]{c}u_{0}\\\
v_{0}\end{array}\right)=T(t)\left(\begin{array}[]{c}0\\\
v_{0}\end{array}\right)+AT(t)\left(\begin{array}[]{c}0\\\
u_{0}\end{array}\right),$

so that the semigroup is determined by its action on initial data with first
component zero (for the same observation for more general hyperbolic equations
see [18, p15]). Said differently, the solution $u$ with initial data
$u(x,0)=u_{0}(x),u_{t}(x,0)=0$ is given by $u=h_{t}$, where $h$ is the
(strong) solution of the wave equation with initial data
$h(x,0)=0,h_{t}(x,0)=u_{0}(x)$.

##  3\. Growth of $L^{2}$ norm as $t\to\infty$.

Theorem 2 gives the extra information that $\|u(\cdot,t)\|_{2}\leqslant
e^{t/2}\|p\|_{H}$, where $\|\cdot\|_{2}$ denotes the norm in $L^{2}(\Omega)$.
However we have the better estimate

(3.1)
$\displaystyle\|u(\cdot,t)\|_{2}^{2}\leqslant\|u_{0}\|_{2}^{2}+2t(u_{0},v_{0})+2E(u_{0},v_{0})t^{2}.$

This follows from the energy equation (2.8) by integrating the identity
$\frac{d}{dt}(u_{t},u)=\|u_{t}\|_{2}^{2}-\|\nabla u\|_{2}^{2}$ to deduce that

(3.2)
$\displaystyle\|u(\cdot,t)\|_{2}^{2}=\|u_{0}\|_{2}^{2}+2t(u_{0},v_{0})+\int_{0}^{t}\int_{0}^{s}(\|u_{t}(\cdot,\tau)\|_{2}^{2}-\|\nabla
u(\cdot,\tau)\|_{2}^{2})\,d\tau\,ds.$

The identity (3.2) follows straightforwardly for $p\in D(A)$, and then for an
arbitrary $p\in H$ via approximation of $p$ by a sequence $p^{(j)}\in D(A)$.

If $\Omega$ is bounded then the energy equation (2.8) and Poincaré inequality
imply that $\|u(\cdot,t)\|_{2}$ is bounded, but for unbounded $\Omega$ it is
possible that $\|u(\cdot,t)\|_{L^{2}(\Omega)}\to\infty$ as $t\to\infty$. For
example, in the case $\Omega=\mathbb{R}^{n}$ with $u_{0}=0$ we have that

$\displaystyle\hat{u}(\xi,t)=\hat{v}_{0}(\xi)\frac{\sin(|\xi|t)}{|\xi|}.$

For $\varepsilon>0$ and $r=|\xi|$ let

$\displaystyle\hat{v}_{0}(\xi)=\left\\{\begin{array}[]{ll}r^{-\frac{n}{2}+\varepsilon},&r\in[0,1),\\\
0,&r\geqslant 1.\end{array}\right.$

Then

$\int_{\mathbb{R}^{n}}|\hat{v}_{0}|^{2}d\xi=\omega_{n}\int_{0}^{1}r^{-1+2\varepsilon}dr<\infty,$

where $\omega_{n}={\mathcal{H}}^{n-1}(S^{n-1})$, so that $v_{0}\in
L^{2}(\mathbb{R}^{n})$ and is real and radially symmetric (because the Fourier
transform of a radially symmetric function is radially symmetric – see, for
example, [27, Corollary 1.2]). But

$\displaystyle\|\hat{u}\|^{2}_{L^{2}(\mathbb{R}^{n})}$
$\displaystyle=\omega_{n}\int_{0}^{1}r^{2\varepsilon-3}\sin^{2}(rt)\,dr$
$\displaystyle=t^{2(1-\varepsilon)}\omega_{n}\int_{0}^{t}s^{2\varepsilon-1}\left(\frac{\sin
s}{s}\right)^{2}\,ds,$

so that by Plancherel’s theorem
$\lim_{t\to\infty}\|u(\cdot,t)\|_{2}t^{\varepsilon-1}=C_{\varepsilon}>0$.

Note, however, that for $\Omega=\mathbb{R}^{n}$ we always have

(3.4) $\displaystyle\lim_{t\to\infty}t^{-1}\|u(\cdot,t)\|_{2}=0.$

This follows from Plancherel’s theorem since

(3.5) $\displaystyle t^{-2}\|\hat{u}(\cdot,t)\|_{2}^{2}\leqslant
2\int_{\mathbb{R}^{n}}\left(|\hat{u}_{0}|^{2}\frac{\cos^{2}(|\xi|t)}{t^{2}}+|\hat{v}_{0}|^{2}\frac{\sin^{2}(|\xi|t)}{|\xi|^{2}t^{2}}\right)\,d\xi,$

which tends to zero as $t\to\infty$ by dominated convergence as
$\frac{\sin\tau}{\tau}$ is bounded.

In fact, for $\Omega=\mathbb{R}^{n},n=1,2,$ estimates in [17] show that blow-
up of the $L^{2}$ norm as $t\to\infty$ occurs when $v_{0}\in
L^{1}(\mathbb{R}^{n})$ and $\int_{\mathbb{R}^{n}}v_{0}\,dx=0$. In the case
$n=1$ the blow-up follows easily from d’Alembert’s formula

(3.6) $\displaystyle
u(x,t)=\frac{1}{2}\left(u_{0}(x+t)+u_{0}(x-t)+\int_{x-t}^{x+t}v_{0}(z)\,dz\right),$

since if $v_{0}\in L^{2}(\mathbb{R})$ with
$\int_{\mathbb{R}}v_{0}^{+}\,dx\neq\int_{\mathbb{R}}v_{0}^{-}\,dx$ then by
Fatou’s Lemma

$\displaystyle\infty$
$\displaystyle=\int_{\mathbb{R}}\lim_{t\to\infty}\left(\int_{x-t}^{x+t}v_{0}(z)\,dz\right)^{2}dx$
$\displaystyle\leqslant\liminf_{t\to\infty}\|2u(\cdot,t)-u_{0}(\cdot+t)-u_{0}(\cdot-t)\|^{2}_{2}$
$\displaystyle\leqslant
8\left(\|u_{0}\|_{2}^{2}+\lim_{t\to\infty}\|u(\cdot,t)\|_{2}^{2}\right).$

If, on the other hand, $v_{0}\in L^{1}(\mathbb{R})$ and
$\int_{\mathbb{R}}v_{0}\,dx=0$ then it can happen either that
$\|u(\cdot,t)\|_{2}$ remains bounded as $t\to\infty$, or that
$\lim_{t\to\infty}\|u(\cdot,t)\|_{2}=\infty$. For the first case one can take
$v_{0}=0$, and for the second take $v_{0}$ odd with

$\displaystyle v_{0}(x)=\left\\{\begin{array}[]{cc}0,&x\in[0,1],\\\
x^{-\alpha},&x\in(1,\infty),\end{array}\right.$

with $1<\alpha<\frac{3}{2}$. Then $v_{0}\in L^{1}(\mathbb{R})\cap
L^{2}(\mathbb{R})$ with $\int_{\mathbb{R}}v_{0}\,dx=0$, and (without loss of
generality taking $u_{0}=0$)

$\displaystyle 4\|u(\cdot,t)\|^{2}_{2}$
$\displaystyle=\int_{\mathbb{R}}\left(\int_{x-t}^{x+t}v_{0}(s)\,ds\right)^{2}dx$
$\displaystyle\geqslant\int_{t+1}^{\infty}\left(\int_{x-t}^{x+t}s^{-\alpha}ds\right)^{2}dx=\int_{t+1}^{\infty}\left(\frac{(x+t)^{1-\alpha}-(x-t)^{1-\alpha}}{1-\alpha}\right)^{2}dx$
$\displaystyle=t^{3-2\alpha}\int_{1+\frac{1}{t}}^{\infty}\left(\frac{(1+{y})^{1-\alpha}-(y-1)^{1-\alpha}}{1-\alpha}\right)^{2}dy\geqslant
Ct^{3-2\alpha}$

for $t\geqslant 1$ and some constant $C>0$.

For any $n$ the set of $v_{0}\in L^{2}(\mathbb{R}^{n})$ such that
$\|u(\cdot,t)\|_{2}$ is bounded as $t\to\infty$ is dense in
$L^{2}(\mathbb{R}^{n})$. Indeed if $\hat{v}_{0}=0$ in some neighbourhood
$B(0,\varepsilon)$ of $0$, then

(3.8) $\displaystyle\|\hat{u}(\cdot,t)\|^{2}_{2}\leqslant
2\left(\|\hat{u}_{0}\|^{2}_{2}+\varepsilon^{-2}\int_{|\xi|>\varepsilon}|\hat{v}_{0}|^{2}d\xi\right)\leqslant
C<\infty$

for all $t\geqslant 0$. The linear subspace of such functions is dense in
$L^{2}(\mathbb{R}^{n})$, since otherwise there would be some nonzero
$\theta\in L^{2}(\mathbb{R}^{n})$ with $(\theta,v_{0})=0$ whenever
$\hat{v}_{0}\in
C_{0}^{\infty}(\mathbb{R}^{n}\setminus\overline{B(0,\varepsilon)})$. But by
Plancherel’s theorem this implies that $(\hat{\theta},\hat{v}_{0})=0$ for all
such $\hat{v}_{0}$, so that $\hat{\theta}(\xi)=0$ for $|\xi|>\varepsilon$ and
any $\varepsilon>0$, thus $\hat{\theta}=0$ and hence $\theta=0$, a
contradiction.

It is not clear to the author whether there can be a solution with
$\|u(\cdot,t)\|_{2}$ unbounded but not tending to infinity as $t\to\infty$.
This seems to depend delicately on the behaviour of $\hat{v}_{0}(\xi)$ as
$|\xi|\to 0$.

##  4\. The adjoint of the wave operator and weak solutions

In order to show that $T(t)p$ in Theorem 2 is the unique weak solution,
appropriately defined, of (1.1) we first recall the definition and properties
(see e.g. [9]) of the adjoint $A^{*}$ of a closed densely defined linear
operator $A$ on a real Banach space $X$ with dual space $X^{*}$. Let
$D(A^{*})$ be the set of those $v\in X^{*}$ for which there exists $v^{*}\in
X^{*}$ such that

$\langle w,v^{*}\rangle=\langle Aw,v\rangle\text{ for all }w\in D(A).$

Since $D(A)$ is dense, $v^{*}$ is unique. Then $A^{*}:D(A^{*})\to X^{*}$ is
the linear operator defined on $X^{*}$ by $A^{*}v=v^{*}$, so that

$\displaystyle\langle w,A^{*}v\rangle=\langle Aw,v\rangle\text{ for all }w\in
D(A).$

$A^{*}$ is closed, and if $X$ is reflexive then $D(A^{*})$ is dense in
$X^{*}$.

######  Definition 4.1.

Let $\tau>0$. A function $w:[0,\infty)\to X$ is a weak solution of
$\dot{w}=Aw$ on $[0,\infty)$ if
(i) $w:[0,\infty)\to X$ is weakly continuous,
(ii) for every $z\in D(A^{*})$ the function $t\mapsto\langle w(t),z\rangle$ is
continuously differentiable on $[0,\infty)$ and

$\frac{d}{dt}\langle w(t),z\rangle=\langle w(t),A^{*}z\rangle,\;t\geqslant 0.$

######  Remark 5.

This is weaker than the definition in [3] of a weak solution in that we do not
assume that $w:[0,\infty)\to X$ is (strongly) continuous.

We will use the following uniqueness result.

######  Theorem 3.

Let $A$ generate the $C^{0}$-semigroup $\\{T(t)\\}_{t\geqslant 0}$ of bounded
linear operators on $X$. Then for any $p\in X$, $w(t):=T(t)p$ is the unique
weak solution of $\dot{w}=Aw$ on $[0,\infty)$ satisfying $w(0)=p$.

###### Proof.

We first show that $w(t)=T(t)p$ is a weak solution. Let $p_{j}\in D(A)$ with
$p_{j}\to p$ in $X$. Then for $t\geqslant 0$ and $z\in D(A^{*})$ we have that

$\displaystyle\langle T(t)p_{j},z\rangle$ $\displaystyle=\langle
p_{j},z\rangle+\int_{0}^{t}\langle AT(s)p_{j},z\rangle ds$ (4.1)
$\displaystyle=\langle p_{j},z\rangle+\int_{0}^{t}\langle
T(s)p_{j},A^{*}z\rangle ds,$

and passing to the limit $j\to\infty$ using the continuity of $T(t)$ and the
boundedness of $\|T(s)p_{j}\|$ on $[0,t]$ we get

(4.2) $\displaystyle\langle T(t)p,z\rangle=\langle
p,z\rangle+\int_{0}^{t}\langle T(s)p,A^{*}z\rangle ds,$

from which (ii) follows.

To prove the uniqueness, suppose that there are two weak solutions
$w,\tilde{w}$ with initial data $p$, and let $W=w-\tilde{w}$. Then
$W:[0,\infty)\to X$ is weakly continuous, hence strongly measurable and
bounded in norm on compact subsets of $[0,\infty)$ [14, pp 59,75,84]. In
particular $W$ is (Bochner) integrable on $[0,t]$ for any $t>0$. Hence, for
any $z\in D(A^{*})$,

(4.3) $\displaystyle\langle W(t),z\rangle=\int_{0}^{t}\langle
W(s),A^{*}z\rangle ds=\langle\int_{0}^{t}W(s)\,ds,A^{*}z\rangle.$

Integrating (4.3) with respect to $t$ we have that

(4.4)
$\displaystyle\langle\int_{0}^{t}W(s)\,ds,z\rangle=\langle\int_{0}^{t}\int_{0}^{s}W(\tau)\,d\tau\,ds,A^{*}z\rangle$

for all $z\in D(A^{*})$, so that by a lemma in [3] (see also [9, p127])
$\int_{0}^{t}\int_{0}^{s}W(\tau)\,d\tau\,ds\in D(A)$ and

(4.5)
$\displaystyle\int_{0}^{t}W(s)\,ds=A\int_{0}^{t}\int_{0}^{s}W(\tau)\,d\tau\,ds,\;t\geqslant
0.$

Hence $r(t):=\int_{0}^{t}\int_{0}^{s}W(\tau)\,d\tau$ is differentiable in $t$
and solves $\dot{r}(t)=Ar(t)$ with $r(0)=0$, so that by a well-known result
(see e.g. [20, p483], [23, Chapter 4], [24, Theorem 35.2]) $r(t)=0$. Hence
also $\int_{0}^{t}W(s)\,ds=0$, thus $\int_{0}^{t}\langle W(s)\,ds,z\rangle=0$
for any $z\in X^{*}$. Differentiating with respect to $t$ and using the
continuity of $\langle W(s),z\rangle$ we have that $\langle W(t),z\rangle=0$
for all $t\geqslant 0$ and $z$. Hence $W=0$ and $w=\tilde{w}$. ∎

######  Remark 6.

This result is given for $X$ a Hilbert space in [2, Corollary 4.8.1]. We note
that the theorem [2, Theorem 4.8.3] of which it is a corollary appears to omit
the hypothesis that $w$ is weakly continuous, and its proof requires some
consequent adjustment.

######  Remark 7.

It is natural to conjecture that for a closed densely defined operator $A$ the
existence of a unique weak solution $w(t)=S(t)p$ satisfying $w(0)=p$ for each
initial data $p\in X$ implies that $A$ generates a $C^{0}$-semigroup
$\\{T(t)\\}_{t\geqslant 0}$ and that $S(t)=T(t)$. This would be a
generalization of the result in [3] to the case when weak solutions are only
required to be weakly continuous in $t$. A crucial step would be to show that
each linear map $S(t)$ is continuous, which was proved in [3] using the closed
graph theorem. However, to generalize this step would seem to require a closed
graph theorem for a linear map from $X$ to the space $C([0,T];X_{w})$, the
space of weakly continuous maps from $[0,T]$ to $X$ with the compact open
topology. Although there are many generalizations of the closed graph theorem
to maps between topological vector spaces (see e.g. [16]), the author was
unable to find one which applies to this case.

Thus we need to calculate the adjoint of the wave operator (2.4) on the
Hilbert space $H=H_{0}^{1}(\Omega)\times L^{2}(\Omega)$.

######  Lemma 4.

The Laplace operator $\Delta$ with $D(\Delta)=\\{u\in H_{0}^{1}(\Omega):\Delta
u\in L^{2}(\Omega)\\}$ is self-adjoint on $L^{2}(\Omega)$.

###### Proof.

This is proved in [1, Example 7.2.1]. Alternatively, one can first note that
if $v\in H_{0}^{1}(\Omega)$ there exists a sequence $\varphi^{(j)}\in
C_{0}^{\infty}(\Omega)$ with $\varphi^{(j)}\to v$ in $H^{1}(\Omega)$, so that
if $u\in D(\Delta)$ we have that

$(-\Delta u,v)=\lim_{j\to\infty}(-\Delta
u,\varphi^{(j)})=\lim_{j\to\infty}(\nabla u,\nabla\varphi^{(j)})=(\nabla
u,\nabla v),$

and hence $(-\Delta u,v)=(u,-\Delta v)$ for $u,v\in D(\Delta)$. Hence
$-\Delta$ is symmetric. As in the proof of Theorem 2, for any $z\in
L^{2}(\Omega)$ there exists a unique $u\in D(\Delta)$ satisfying $-\Delta
u+u=z$. Hence $-\Delta$ is also maximal monotone, so that the self-adjointness
follows from [6, Proposition 7.6]. ∎

For $\theta\in L^{2}(\Omega)$ denote by $({\mathbf{1}}-\Delta)^{-1}\theta$ the
unique solution $u\in D(\Delta)$ of $-\Delta u+u=\theta$.

######  Theorem 5.

The adjoint of the wave operator $A$ is given by

$A^{*}=\left(\begin{array}[]{cc}0&({\mathbf{1}}-\Delta)^{-1}-{\mathbf{1}}\\\
{\mathbf{1}}-\Delta&0\end{array}\right),$

with

$D(A^{*})=\\{\left(\begin{array}[]{l}\chi\\\ \psi\end{array}\right)\in
H:\Delta\chi\in L^{2}(\Omega),\psi\in H_{0}^{1}(\Omega)\\}.$

###### Proof.

By the definition of $A^{*}$ we have that $\left(\begin{array}[]{ll}\chi\\\
\psi\end{array}\right)\in D(A^{*})\text{ and
}A^{*}\left(\begin{array}[]{ll}\chi\\\
\psi\end{array}\right)=\left(\begin{array}[]{ll}p\\\ q\end{array}\right)$ if
and only if $\left(\begin{array}[]{ll}\chi\\\ \psi\end{array}\right)\in H$ and

(4.16) $\displaystyle\langle\left(\begin{array}[]{ll}u\\\
v\end{array}\right),\left(\begin{array}[]{ll}p\\\
q\end{array}\right)\rangle=\langle\left(\begin{array}[]{ll}v\\\ \Delta
u\end{array}\right),\left(\begin{array}[]{ll}\chi\\\
\psi\end{array}\right)\rangle\text{ for all }\left(\begin{array}[]{ll}u\\\
v\end{array}\right)\in D(A),$

that is

(4.17) $\displaystyle\int_{\Omega}(pu+\nabla p\cdot\nabla
u+qv)\,dx=\int_{\Omega}(\chi v+\nabla\chi\cdot\nabla v+\psi\Delta u)\,dx$

for all $u,v\in H_{0}^{1}(\Omega)$ with $\Delta u\in L^{2}(\Omega)$. (4.17) is
equivalent to the two equations

(4.18) $\displaystyle\int_{\Omega}(\chi v+\nabla\chi\cdot\nabla v)\,dx$
$\displaystyle=\int_{\Omega}qv\,dx\;\text{ for all }v\in H_{0}^{1}(\Omega),$
(4.19) $\displaystyle\int_{\Omega}(pu+\nabla p\cdot\nabla u)\,dx$
$\displaystyle=\int_{\Omega}\psi\Delta u\,dx\;\text{ for all }u\in D(\Delta).$

But (4.18) says that $q=({\mathbf{1}}-\Delta)\chi$, while (4.19) can be
written as

(4.20) $\displaystyle\int_{\Omega}pu\,dx=\int_{\Omega}(p+\psi)\Delta
u\,dx\;\text{ for all }u\in D(\Delta),$

so that, by Lemma 4, $p+\psi\in D(\Delta)$ and $\Delta(p+\psi)=p$, from which
it follows that $p=[({\mathbf{1}}-\Delta)^{-1}-{\mathbf{1}}]\psi$. ∎

######  Remark 8.

For a bounded domain $\Omega$ one can use the equivalent inner product
$((u,v))=\int_{\Omega}\nabla u\cdot\nabla v\,dx$ on $H_{0}^{1}(\Omega)$, when
the adjoint takes the simpler form (see [5, Lemma 3.1])
$A^{*}=-\left(\begin{array}[]{cc}0&{\mathbf{1}}\\\
\Delta&0\end{array}\right)$.

By Definition 4.1, $w=\left(\begin{array}[]{l}u\\\ v\end{array}\right)$ is a
weak solution of the wave equation $\dot{w}=Aw$ on $[0,\infty)$ if and only if
$w:[0,\infty)\to H$ is weakly continuous, and, for any
$z=\left(\begin{array}[]{l}\chi\\\ \psi\end{array}\right)\in D(A^{*})$, the
function $t\mapsto\langle w(t),z\rangle$ is continuously differentiable with
derivative

(4.21) $\displaystyle\frac{d}{dt}\langle w(t),z\rangle=\langle
w(t),A^{*}z\rangle,\;t\geqslant 0.$

Equivalently, for $t\geqslant 0$,

(4.22) $\displaystyle\frac{d}{dt}\int_{\Omega}(u\chi+\nabla
u\cdot\nabla\chi+v\psi)\,dx=\int_{\Omega}(up+\nabla u\cdot\nabla p+vq)\,dx,$

where $\Delta(p+\psi)=p$ and $q=\chi-\Delta\chi$, or

(4.23) $\displaystyle\frac{d}{dt}\int_{\Omega}uq\,dx$
$\displaystyle=\int_{\Omega}vq\,dx\;\text{ for all }q\in L^{2}(\Omega),$
(4.24) $\displaystyle\frac{d}{dt}\int_{\Omega}v\psi\,dx$
$\displaystyle=-\int_{\Omega}\nabla u\cdot\nabla\psi\,dx\;\text{ for all
}\psi\in H_{0}^{1}(\Omega),$

since for any $q\in L^{2}(\Omega)$ there exists a unique solution $\chi\in
H_{0}^{1}(\Omega)$ to $\chi-\Delta\chi=q$. But it is easily checked that
(4.23) holds if and only if $u$ is weakly differentiable with respect to $t$
with $u_{t}(\cdot,t)=v(\cdot,t)$ for all $t\geqslant 0$, that is

(4.25)
$\displaystyle\int_{0}^{\infty}\varphi^{\prime}(t)u(\cdot,t)\,dt=-\int_{0}^{\infty}\varphi(t)v(\cdot,t)\,dt\text{
for all }\varphi\in C_{0}^{\infty}(0,\infty).$

Then (4.24) becomes

(4.26)
$\displaystyle\frac{d}{dt}\int_{\Omega}u_{t}\psi\,dx=-\int_{\Omega}\nabla
u\cdot\nabla\psi\,dx\text{ for all }\psi\in H_{0}^{1}(\Omega),t\geqslant 0.$

Hence $w$ is a weak solution if and only if $u:[0,\infty)\to
H_{0}^{1}(\Omega)$, $v:[0,\infty)\to L^{2}(\Omega)$ are weakly continuous,
$v=u_{t}$, and (4.26) holds. In particular weak solutions in this sense are
unique.

######  Remark 9.

For the case $\Omega=\mathbb{R}^{n}$ it is possible to prove uniqueness of
weak solutions for $u_{0},v_{0}$ merely distributions using properties of the
fundamental solution of the wave equation (see [28, Theorem 13.1]).

Next we give some further properties of $A,A^{*}$ from which we deduce the
absence of nontrivial linear conserved quantities.

######  Lemma 6.

$(i)$ $\Delta D(\Delta)$ is dense in $L^{2}(\Omega)$.
$(ii)$ $R(A)=\\{Az:z\in D(A)\\}$ is dense in $H$.
$(iii)$ $A^{*}:D(A^{*})\to H=H^{*}$ is one-to-one.

###### Proof.

(i). Suppose that $\int_{\Omega}z\Delta u\,dx=0$ for some $z\in L^{2}(\Omega)$
and all $u\in D(\Delta)$. By Lemma 4, $z\in D(\Delta)$ and $\Delta z=0$.
Choosing $u=z$ we thus have $\int_{\Omega}|\nabla z|^{2}dx=0$ and so $\nabla
z=0$ in $\Omega$. Hence $z=0$ (for example because the extension $\tilde{z}$
of $z$ by zero belongs to $H_{0}^{1}(\mathbb{R}^{n})$ and $\nabla\tilde{z}=0$,
so that $\tilde{z}$ is constant and thus zero).
(ii) If not there would exist a nonzero $\left(\begin{array}[]{l}\chi\\\
\psi\end{array}\right)\in H$ with $\langle\left(\begin{array}[]{l}\chi\\\
\psi\end{array}\right),A\left(\begin{array}[]{l}u\\\
v\end{array}\right)\rangle=0$ for all $\left(\begin{array}[]{l}u\\\
v\end{array}\right)\in D(A)$, that is

(4.27) $\displaystyle\int_{\Omega}(\chi v+\nabla\chi\cdot\nabla v+\psi\Delta
u)\,dx=0$

for all $u,v\in H_{0}^{1}(\Omega)$ with $\Delta u\in L^{2}(\Omega)$. Taking
first $u=0$ and $v=\chi$ we get that $\chi=0$. Then $\psi=0$ by (i).
(iii) Now suppose $A^{*}\left(\begin{array}[]{l}\chi\\\
\psi\end{array}\right)=0$, for some $\left(\begin{array}[]{l}\chi\\\
\psi\end{array}\right)\in D(A^{*})$. Then

$\langle\left(\begin{array}[]{l}u\\\
v\end{array}\right),A^{*}\left(\begin{array}[]{l}\chi\\\
\psi\end{array}\right)\rangle=0$

for all $\left(\begin{array}[]{l}u\\\ v\end{array}\right)\in D(A)$, so that by
(ii) $\chi=\psi=0$. Hence $A^{*}$ is one-to-one. ∎

######  Theorem 7.

There is no nontrivial linear constant of motion for (2.1), that is there is
no nonzero $z\in H$ such that $\langle T(t)p,z\rangle=\langle p,z\rangle$ for
all $t\geqslant 0$ and $p\in D(A)$.

###### Proof.

If $\langle T(t)p,z\rangle$ were constant in $t$ for all $p\in D(A)$ then

$0=\left.\frac{d}{dt}\langle T(t)p,z\rangle\right|_{t=0}=\langle Ap,z\rangle$

for all $p\in D(A)$, so that $z=0$ by Lemma 6 (ii). ∎

Finally in this section we motivate the use of the phase space
$H=H_{0}^{1}(\Omega)\times L^{2}(\Omega)$ for an unbounded domain
$\Omega\subset\mathbb{R}^{n}$. To this end we assume that
$\Omega=\bigcup_{j=1}^{\infty}\Omega_{j}$ is the union of an increasing
($\Omega_{j}\subset\Omega_{j+1}$) sequence of bounded open sets
$\Omega_{j}\subset\mathbb{R}^{n}$. In Theorem 9 below we show that the
semiflow $\\{T(t)\\}_{t\geqslant 0}$ on $H$ for the wave equation is the limit
of the corresponding semiflows $\\{T_{j}(t)\\}_{t\geqslant 0}$ for the wave
equation on $\Omega_{j}$ with phase space $H_{j}=H_{0}^{1}(\Omega_{j})\times
L^{2}(\Omega_{j})$. We use the following lemma, which is a slight
generalization of [4, Lemma 5.12].

######  Lemma 8.

Let $X$ be a reflexive Banach space, $T>0$, and $w^{(j)}:[0,T]\to X$ satisfy
$\rm(i)$ $\|w^{(j)}(t)\|_{X}\leqslant M<\infty$ for all $j=1,2,\ldots$ and
$t\in[0,T]$,
$\rm(ii)$ there is a dense subset $E$ of $X^{*}$ such that for any $v\in E$
the functions $t\mapsto\langle w^{(j)},v\rangle$ are equicontinuous on $[0,T]$
for $j$ sufficiently large, that is given $\varepsilon>0$, there exists
$\delta(v)>0$ and $J(v)$ such that

(4.28) $\displaystyle|\langle w^{(j)}(t),v\rangle-\langle
w^{(j)}(s),v\rangle|<\varepsilon\text{ for
}|s-t|\leqslant\delta(v),\;j\geqslant J(v).$

Then there exists a weakly continuous $w:[0,T]\to X$ and a subsequence
$w^{(j_{k})}$ of $w^{(j)}$ converging uniformly to $w$ in the weak topology,
i.e. for any sequence $s_{k}\to s$ in $[0,T]$ we have
$w^{(j_{k})}(s_{k})\rightharpoonup w(s)$ in $X$.

###### Proof.

By (i) and a diagonal argument we can extract a subsequence $w^{(j_{k})}$ such
that $w^{(j_{k})}(\tau)$ converges weakly to a limit for any rational
$\tau\in[0,T]$. We claim that $w^{(j_{k})}(t)$ converges weakly to a limit
$w(t)$ for any $t\in[0,T]$. This follows provided $\langle
w^{(j_{k})}(t),v\rangle$ is a Cauchy sequence for any $v\in X^{*}$, since then
by (i) the limit is a bounded linear functional of $v$, so that since $X$ is
reflexive it defines an element of $X^{**}=X$. To prove this, let
$\varepsilon>0$ and choose $\tilde{v}\in E$ with
$\|\tilde{v}-v\|_{X^{*}}<\frac{\varepsilon}{2M}$, and then a rational
$\tau\in[0,T]$ with $|\tau-t|\leqslant\delta(\tilde{v})$. Then for $k,l$
sufficiently large we have by (i) and (4.28) that

$\displaystyle|\langle w^{(j_{k})}(t)-w^{(j_{l})}(t),v\rangle|$
$\displaystyle\leqslant|\langle
w^{(j_{k})}(t)-w^{(j_{l})}(t),\tilde{v}\rangle|+\varepsilon$
$\displaystyle\leqslant|\langle
w^{(j_{k})}(\tau)-w^{(j_{l})}(\tau),\tilde{v}\rangle|+3\varepsilon\leqslant
4\varepsilon,$

as required. Similar arguments then show that $w$ is weakly continuous and
that $w^{(j_{k})}(s_{k})\rightharpoonup w(s)$ if $s_{k}\to s$. ∎

For a function $f\in H_{j}$ set

(4.31)
$\displaystyle\bar{f}(x)=\left\\{\begin{array}[]{cl}f(x),&x\in\Omega_{j}\\\
0,&x\in\Omega\setminus\Omega_{j}.\end{array}\right.$

Note that if $f\in H_{j}$ then $\bar{f}\in H$.

######  Theorem 9.

Let $p_{j}\in H_{j}$ and $\bar{p}_{j}\to p$ in $H$. Then
$\overline{T_{j}(t)p_{j}}\to T(t)p$ in $H$ uniformly on compact subsets of
$[0,\infty)$.

###### Proof.

We denote by $A_{j}$ the infinitesimal generator of $T_{j}(t)$, that is

(4.34) $\displaystyle A_{j}=\left(\begin{array}[]{ll}0&{\mathbf{1}}\\\
\Delta_{j}&0\end{array}\right),$

with

$D(A_{j})=\\{\left(\begin{array}[]{l}u\\\ v\end{array}\right)\in
H:\Delta_{j}u\in L^{2}(\Omega_{j}),v\in H_{0}^{1}(\Omega_{j})\\},$

where $\Delta_{j}=\Delta$ with domain $D(\Delta_{j})=\\{u\in
H_{0}^{1}(\Omega_{j}):\Delta u\in L^{2}(\Omega_{j})\\}$. Then by Theorem 5 we
have that

$A_{j}^{*}=\left(\begin{array}[]{cc}0&({\mathbf{1}}-\Delta_{j})^{-1}-{\mathbf{1}}\\\
{\mathbf{1}}-\Delta_{j}&0\end{array}\right),$

with

$D(A_{j}^{*})=\\{\left(\begin{array}[]{l}\chi\\\ \psi\end{array}\right)\in
X_{j}:\Delta\chi\in L^{2}(\Omega_{j}),\psi\in H_{0}^{1}(\Omega_{j})\\}.$

By Theorem 3 we have that

(4.35) $\displaystyle\langle T_{j}(t)p_{j},v\rangle=\langle
p_{j},v\rangle+\int_{0}^{t}\langle T_{j}(s)p_{j},A^{*}_{j}v\rangle ds$

for all $t\geqslant 0$ and $v\in D(A_{j}^{*})$. Furthermore, if $T>0$ there
exists $M>0$ such that

(4.36)
$\displaystyle\|\overline{T_{j}(t)p_{j}}\|_{H}=\|T_{j}(t)p_{j}\|_{H_{j}}\leqslant
M\text{ for all }t\in[0,T],$

this following from the energy equation (2.8) and the estimate (3.1). Let
$E=\\{\left(\begin{array}[]{cc}\chi\\\ \psi\end{array}\right):\chi,\psi\in
C_{0}^{\infty}(\Omega)\\}$, which is a dense subset of $H^{*}=H$. Given any
$v\in E$ we have that $v\in D(A_{j}^{*})$ for large enough $j$. Thus from
(4.35), (4.36) we have that for $v\in E$ and large enough $j$

(4.37)
$\displaystyle|\langle{T_{j}(t)p_{j}},v\rangle-\langle{T_{j}(s)p_{j}},v\rangle|\leqslant
C|t-s|,\text{ for }s,t\in[0,T],$

where $C=C(v)$ is a constant, and
$\langle\overline{T_{j}(t)p_{j}},v\rangle=\langle{T_{j}(t)p_{j}},v\rangle$ for
each $t$. Hence $w^{(j)}(t):=\overline{T_{j}(t)p_{j}}$ satisfies the
hypotheses of Lemma 8 for any $T>0$ and so there is a subsequence
$w^{(j_{k})}$ and a weakly continuous $w:[0,\infty)\to H$ such that
$w^{(j_{k})}$ converges uniformly to $w$ on compact subsets of $[0,\infty)$ in
the weak topology. Writing $p_{j}=\left(\begin{array}[]{c}u_{0j}\\\
v_{0j}\end{array}\right)$,
$T_{j}(t)p_{j}=\left(\begin{array}[]{c}u_{j}(\cdot,t)\\\
v_{j}(\cdot,t)\end{array}\right)$, we have from (4.23), (4.24) that for any
$q,\psi\in C_{0}^{\infty}(\Omega)$ and $t\geqslant 0$, and for $k$
sufficiently large to ensure ${\rm supp}\,q\subset\Omega_{j_{k}},\,{\rm
supp}\,\psi\subset\Omega_{j_{k}}$,

(4.38) $\displaystyle\int_{\Omega}u_{j_{k}}(x,t)q(x)\,dx$
$\displaystyle=\int_{\Omega}u_{0j_{k}}(x)q(x)\,dx+\int_{0}^{t}\int_{\Omega}v_{j_{k}}(x,s)q(x)\,dx\,ds,$
(4.39) $\displaystyle\int_{\Omega}v_{j_{k}}(x,t)\psi(x)\,dx$
$\displaystyle=\int_{\Omega}v_{0j_{k}}(x)\psi(x)\,dx-\int_{0}^{t}\int_{\Omega}\nabla
u_{j_{k}}(x,s)\cdot\nabla\psi(x)\,dx\,ds.$

Passing to the limit $k\to\infty$ and setting
$w(t)=\left(\begin{array}[]{c}u(\cdot,t)\\\ v(\cdot,t)\end{array}\right)$,
$p=\left(\begin{array}[]{c}u_{0}\\\ v_{0}\end{array}\right)$, we get that

(4.40) $\displaystyle\int_{\Omega}u(x,t)q(x)\,dx$
$\displaystyle=\int_{\Omega}u_{0}(x)q(x)\,dx+\int_{0}^{t}\int_{\Omega}v(x,s)q(x)\,dx\,ds,$
(4.41) $\displaystyle\int_{\Omega}v(x,t)\psi(x)\,dx$
$\displaystyle=\int_{\Omega}v_{0}(x)\psi(x)\,dx-\int_{0}^{t}\int_{\Omega}\nabla
u(x,s)\cdot\nabla\psi(x)\,dx\,ds.$

By approximation (4.40), (4.41) hold for all $q\in L^{2}(\Omega),\psi\in
H_{0}^{1}(\Omega)$, so that $w$ is a weak solution on $[0,\infty)$ with
$w(0)=p$. Hence, by Theorem 3, $w(t)=T(t)p$. The uniqueness also implies by a
standard argument that the whole sequence $w^{(j)}$ converges, so that if
$s_{j}\to s$ in $[0,\infty)$ then $\overline{T_{j}(s_{j})p_{j}}\rightharpoonup
T(s)p$ in $H$. Since by the energy equation, for all $t\geqslant 0$

$\displaystyle\int_{\Omega}\left(|\nabla
u_{j}(x,t)|^{2}+|v_{j}(x,t)|^{2}\right)\,dx$
$\displaystyle=\int_{\Omega}\left(|\nabla u_{0j}|^{2}+|v_{0j}|^{2}\right)\,dx$
(4.42) $\displaystyle\to\int_{\Omega}\left(|\nabla
u_{0}|^{2}+|v_{0}|^{2}\right)\,dx$ $\displaystyle=\int_{\Omega}\left(|\nabla
u(x,t)|^{2}+|v(x,t)|^{2}\right)\,dx,$

we have that $\nabla u_{j}(\cdot,s_{j})\to\nabla u(\cdot,s)$ strongly in
$L^{2}(\Omega)^{n}$ and $v_{j}(\cdot,s_{j})\to v(\cdot,s)$ strongly in
$L^{2}(\Omega)$. From (3.2) we deduce that also $u_{j}(\cdot,s_{j})\to
u(\cdot,s)$ strongly in $L^{2}(\Omega)$ and hence
$\overline{T_{j}(s_{j})p_{j}}\to T(s)p$ strongly in $X$ as claimed. ∎

######  Remark 10.

See [10] for a discussion on solving the wave equation in $\mathbb{R}^{n}$ in
a different energy space.

##  5\. Kirchhoff’s formula and smoothing via averaging

We return to Kirchhoff’s formula for a $C^{2}$ solution of the wave equation
for $n=3$ and initial data $u(\cdot,0)=u_{0}$, $u_{t}(\cdot,0)=v_{0}$, namely

(5.1) $\displaystyle u(x,t)=\frac{1}{4\pi
t}\int_{S(x,t)}v_{0}(y)\,dS_{y}+\frac{\partial}{\partial t}\left(\frac{1}{4\pi
t}\int_{S(x,t)}u_{0}(y)\,dS_{y}\right).$

For simplicity, and bearing in mind Remark 4, we suppose that $u_{0}=0$, so
that

$\displaystyle u(x,t)$ $\displaystyle=\frac{1}{4\pi
t}\int_{S(x,t)}v_{0}(y)\,dS_{y}$ (5.2)
$\displaystyle=\frac{t}{4\pi}\int_{S^{2}}v_{0}(x+tz)\,dS_{z}$

Hence, formally we also have

$\displaystyle u_{t}(x,t)$
$\displaystyle=\frac{1}{t}u(x,t)+\frac{t}{4\pi}\int_{S^{2}}\nabla
v_{0}(x+tz)\cdot z\,dS_{z}$ (5.3)
$\displaystyle=\frac{1}{t}u(x,t)+\frac{t^{2}}{4\pi}\Delta_{x}\int_{B}v_{0}(x+tz)\,dz,$

where $B=B(0,1)$.

The representations (5.2), (5.3) are at first sight puzzling, because in view
of Theorem 2 we expect them to be meaningful when we just have $v_{0}\in
L^{2}(\Omega)$. However they can be understood because of the smoothing
properties of averaging over spheres and balls that are well known in harmonic
analysis. Indeed Stein [25] proved in particular that if we define the
spherical mean

(5.4) $\displaystyle{\mathcal{M}}_{t}(f)(x)=\int_{S^{2}}f(x+tz)\,dS_{z},$

then there is a constant $C>0$, independent of $t$, such that

(5.5) $\displaystyle\|{\mathcal{M}}_{t}(f)\|_{2}\leqslant C\|f\|_{2}$

for all $f\in C_{0}^{\infty}(\mathbb{R}^{3})$ and $t>0$. Approximating $v_{0}$
in $L^{2}(\mathbb{R}^{3})$ by functions $f^{(j)}\in
C_{0}^{\infty}(\mathbb{R}^{3})$ and applying (5.5) with $f=f^{(j)}-f^{(k)}$ we
find that ${\mathcal{M}}_{t}(f^{(j)})$ is a Cauchy sequence in
$L^{2}(\mathbb{R}^{3})$, so that the integral in (5.2) can be defined for each
$t$ as an element of $L^{2}(\mathbb{R}^{3})$. Furthermore we can estimate
derivatives of $M_{t}(f)$, as is done in [29]. Thus we find that

(5.6) $\displaystyle\|\frac{\partial}{\partial
x_{i}}{\mathcal{M}}_{t}(f)\|_{2}\leqslant C_{1}t^{-1}\|f\|_{2}$

for some constant $C_{1}$ and $i=1,2,3$. In [29] it is also shown that for
$B=\\{x\in\mathbb{R}^{3}:|x|<1\\}$ the ball average

(5.7) $\displaystyle{\mathcal{N}}_{t}(f)(x):=\int_{B}f(x+tz)\,dz,$

satisfies the estimates

(5.8) $\displaystyle\|{\mathcal{N}}_{t}(f)\|_{2}\leqslant
C_{2}\|f\|_{2},\;\|\frac{\partial}{\partial
x_{i}}{\mathcal{N}}_{t}(f)\|_{2}\leqslant
C_{2}t^{-1}\|f\|_{2},\;\|\Delta{\mathcal{N}}_{t}(f)\|_{2}\leqslant
C_{2}t^{-2}\|f\|_{2},$

for some constant $C_{2}$.

In fact the derivation of Kirchhoff’s solution combined with Theorem 2 gives
an alternative proof of the estimates (5.5),(5.6),(5.8). It suffices to take
$\varphi^{(j)}\in C_{0}^{\infty}(\mathbb{R}^{n})$ with $\varphi^{(j)}\to
v_{0}$ in $L^{2}(\mathbb{R}^{3})$. Then the solution $u^{(j)}$ with
$u^{(j)}(\cdot,0)=0,\,u_{t}^{(j)}(\cdot,0)=\varphi^{(j)}$ is smooth and given
by

(5.9) $\displaystyle u^{(j)}(x,t)=\frac{1}{4\pi
t}\int_{S(x,t)}\varphi^{(j)}(y)\,dS_{y},$

and satisfies

(5.10) $\displaystyle
u^{(j)}_{t}(x,t)=\frac{1}{t}u^{(j)}(x,t)+\frac{t^{2}}{4\pi}\Delta_{x}\int_{B}\varphi^{(j)}(x+tx)\,dz.$

By Theorem 2 we have that $u^{(j)}\to u$ in
$C([0,\tau];H_{0}^{1}(\mathbb{R}^{3}))$ and $u^{(j)}_{t}\to u_{t}$ in
$C([0,\tau];L^{2}(\mathbb{R}^{3}))$ for any $\tau>0$, where $u$ is the unique
weak solution with initial data $u(\cdot,0)=0,u_{t}(\cdot,0)=v_{0}$ , and thus
is given by (5.2). Hence, setting $f=v_{0}$, we obtain (5.5) from (3.1), and
(5.6) from (2.8). The first estimate in (5.8) is immediate since

$\|{\mathcal{N}}_{t}(f)\|_{2}^{2}\leqslant\int_{B}\left(\int_{B}1^{2}dz\int_{\mathbb{R}^{3}}f(x+tz)^{2}dz\right)dx.$

By (5.10), (3.1), (2.8) we have that
$t^{2}\|\Delta{\mathcal{N}}_{t}(\varphi^{(j)})\|_{2}\leqslant M<\infty$, from
which we get the third estimate. The middle estimate then follows using the
relation $(g,\Delta g)=-\|\nabla g\|_{2}^{2}$.

As an illustration of the harmonic analysis methods used to derive estimates
such as (5.5),(5.6),(5.8) we prove the following result.

######  Theorem 10.

Let $f\in L^{2}(\mathbb{R}^{n})$. Then for $B=\\{x\in\mathbb{R}^{n}:|x|<1\\}$
the average ${\mathcal{N}}_{t}(f)$ defined for $t>0$ by

(5.11) $\displaystyle{\mathcal{N}}_{t}(f)(x)=\int_{B}f(x+tz)\,dz$

belongs to $H^{\frac{n+1}{2}}:=H^{\frac{n+1}{2}}(\mathbb{R}^{n})$ with

(5.12) $\displaystyle\|{\mathcal{N}}_{t}(f)\|_{H^{\frac{n+1}{2}}}\leqslant
C\max(1,t^{-\frac{n+1}{2}})\|f\|_{2},$

for some constant $C>0$ independent of $t$.

###### Proof.

We use the fact (see [11, p175], [22, pp 605-606], [26, p 338]) that the
Fourier transform $\hat{\chi}_{B}$ of the characteristic function $\chi_{B}$
of the unit ball $B$ satisfies

(5.13) $\displaystyle\hat{\chi}_{B}(\xi)\leqslant
C_{n}(1+|\xi|)^{-\frac{n+1}{2}}$

for some positive constant $C_{n}$. We note that $h:={\mathcal{N}}_{1}(f)$
satisfies $h(x)=(\chi_{B}*f)(x)$, so that
$\hat{h}(\xi)=(2\pi)^{\frac{n}{2}}\hat{\chi}_{B}(\xi)\hat{f}(\xi),$ and thus
by (5.13)

(5.14) $\displaystyle(1+|\xi|^{\frac{n+1}{2}})|\hat{h}(\xi)|\leqslant
C|\hat{f}(\xi)|.$

Hence
$\|h\|_{H^{\frac{n+1}{2}}}=\|(1+|\xi|^{\frac{n+1}{2}})\hat{h}\|_{2}\leqslant
C\|\hat{f}\|_{2}=C\|f\|_{2}$, and hence (5.12) holds for $t=1$.

Set $f_{t}(x)=f(tx)$. Since
${\mathcal{N}}_{t}(f)(tx)={\mathcal{N}}_{1}(f_{t})(x)$ we have
$\widehat{{\mathcal{N}}_{t}(f)}(\xi)=t^{n}\widehat{{\mathcal{N}}_{1}(f_{t})}(t\xi)$,
so that by (5.14)

$\displaystyle(1+|\xi|^{\frac{n+1}{2}})|\widehat{{\mathcal{N}}_{t}(f)}|(\xi)$
$\displaystyle\leqslant
C\left(\frac{1+|\xi|^{\frac{n+1}{2}}}{1+|t\xi|^{\frac{n+1}{2}}}\right)t^{n}|\widehat{f_{t}}(t\xi)|$
(5.15) $\displaystyle\leqslant C\max(1,t^{-\frac{n+1}{2}})|\hat{f}(\xi)|,$

giving (5.12) for any $t$. ∎

######  Remark 11.

Because of results of Hlawka [15], Herz [13] the same result holds if $B$ is
replaced by a bounded convex set $C\subset\mathbb{R}^{n}$ with sufficiently
smooth boundary having everywhere positive Gaussian curvature. However if $B$
is replaced by the cube $Q=(-1,1)^{n}$ then
${\mathcal{N}}^{Q}_{1}(f)(x):=\int_{Q}f(x+z)\,dz$ has less regularity than
${\mathcal{N}}_{1}(f)$. In fact

$\displaystyle\hat{\chi}_{Q}(\xi)$
$\displaystyle=\frac{1}{(2\pi)^{\frac{n}{2}}}\int_{Q}e^{-i\xi\cdot x}\,dx$
$\displaystyle=\frac{2^{n}}{(2\pi)^{\frac{n}{2}}}\prod_{j=1}^{n}\frac{\sin\xi_{j}}{\xi_{j}}.$

Hence if $\alpha\geqslant 0$ then $(1+|\xi|^{\alpha})\hat{\chi}_{Q}(\xi)\in
L^{\infty}(\mathbb{R}^{n})$ iff $\alpha\leqslant 1$. Hence
${\mathcal{N}}^{Q}_{1}(f)\in H^{1}(\mathbb{R}^{n})$ but in general
${\mathcal{N}}^{Q}_{1}(f)\not\in H^{\alpha}(\mathbb{R}^{n})$ for $\alpha>1$.

Acknowledgement. This paper was completed while visiting the Hong Kong
Institute for Advanced Study as a Senior Fellow.

## References

  * [1] W. Arendt, C. J. K. Batty, M. Hieber, and F. Neubrander.  Vector-valued Laplace transforms and Cauchy problems, volume 96 of Monographs in Mathematics.  Birkhäuser/Springer Basel AG, Basel, second edition, 2011.
  * [2] A. V. Balakrishnan.  Applied functional analysis, volume 3 of Applications of Mathematics.  Springer-Verlag, New York, second edition, 1981.
  * [3] J. M. Ball.  Strongly continuous semigroups, weak solutions, and the variation of constants formula.  Proc. Amer. Math. Soc., 64:370–373, 1977.
  * [4] J. M. Ball.  On the asymptotic behaviour of generalized processes, with applications to nonlinear evolution equations.  J. Differential Eqns, 27:224–265, 1978.
  * [5] J. M. Ball.  Global attractors for damped semilinear wave equations.  Discrete and Continuous Dynamical Systems, 10:31–52, 2004.
  * [6] H. Brezis.  Functional analysis, Sobolev spaces and partial differential equations.  Universitext. Springer, New York, 2011.
  * [7] L. C. Evans.  Partial differential equations, volume 19 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, second edition, 2010.
  * [8] H. O. Fattorini.  Second order linear differential equations in Banach spaces, volume 108 of North-Holland Mathematics Studies.  North-Holland Publishing Co., Amsterdam, 1985.  Notas de Matemática, 99. [Mathematical Notes].
  * [9] S. Goldberg.  Unbounded linear operators: Theory and applications.  McGraw-Hill Book Co., New York-Toronto-London, 1966.
  * [10] J. A. Goldstein and M. Wacker.  The energy space and norm growth for abstract wave equations.  Appl. Math. Lett., 16(5):767–772, 2003.
  * [11] L. Grafakos.  Classical Fourier analysis, volume 249 of Graduate Texts in Mathematics.  Springer, New York, second edition, 2008.
  * [12] S. Helgason.  The Radon transform, volume 5 of Progress in Mathematics.  Birkhäuser Boston, Inc., Boston, MA, second edition, 1999.
  * [13] C. S. Herz.  Fourier transforms related to convex sets.  Ann. of Math. (2), 75:81–92, 1962.
  * [14] E. Hille and R. S. Phillips.  Functional analysis and semi-groups, volume 31 of Colloq. Publ. Col.. Amer. Math. Soc., Providence, 1957.
  * [15] E. Hlawka.  Über Integrale auf konvexen Körpern. I.  Monatsh. Math., 54:1–36, 1950.
  * [16] T. Husain.  The open mapping and closed graph theorems in topological vector spaces. Clarendon Press, Oxford,, 1965.
  * [17] R. Ikehata.  $L^{2}$-blowup estimates of the wave equation and its application to local energy decay.  J. Hyperbolic Differ. Equ., 20(1):259–275, 2023.
  * [18] F. John.  Plane waves and spherical means applied to partial differential equations.  Interscience Publishers, New York-London, 1955.
  * [19] F. John.  Partial differential equations, volume 1 of Applied Mathematical Sciences.  Springer-Verlag, New York, fourth edition, 1982.
  * [20] T. Kato.  Perturbation theory for linear operators. Springer-Verlag, Berlin-New York,,, second edition, 1976.
  * [21] P. D. Lax.  Hyperbolic partial differential equations, volume 14 of Courant Lecture Notes in Mathematics.  New York University, Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, 2006.  With an appendix by Cathleen S. Morawetz.
  * [22] B. Makarov and A. Podkorytov.  Real analysis: measures, integrals and applications.  Universitext. Springer, London, 2013.  Translated from the 2011 Russian original.
  * [23] A. Pazy.  Semigroups of linear operators and applications to partial differential equations, volume 44 of Applied Mathematical Sciences.  Springer-Verlag, New York, 1983.
  * [24] G. R. Sell and Y. You.  Dynamics of evolutionary equations.  Springer Verlag.  Book in preparation.
  * [25] E. M. Stein.  Maximal functions: Poisson integrals on symmetric spaces.  Proc. Nat. Acad. Sci. U.S.A., 73(8):2547–2549, 1976.
  * [26] E. M. Stein.  Harmonic analysis: real-variable methods, orthogonality, and oscillatory integrals, volume 43 of Princeton Mathematical Series.  Princeton University Press, Princeton, NJ, 1993.  With the assistance of Timothy S. Murphy, Monographs in Harmonic Analysis, III.
  * [27] E. M. Stein and G. Weiss.  Introduction to Fourier analysis on Euclidean spaces.  Princeton Mathematical Series, No. 32. Princeton University Press, Princeton, N.J., 1971.
  * [28] F. Trèves.  Basic linear partial differential equations.  Pure and Applied Mathematics, Vol. 62. Academic Press [Harcourt Brace Jovanovich, Publishers], New York-London, 1975.
  * [29] J. Zhao and D. Fan.  Derivative estimates of averaging operators and extension.  Front. Math. China, 14(2):475–491, 2019.


# Confronting the neutron star population with inverse cascades

Nikhil Sarin  Nordita, Stockholm University and KTH Royal Institute of
Technology
Hannes Alfvéns väg 12, SE-106 91 Stockholm, Sweden  Oskar Klein Centre for
Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova,
Stockholm SE-106 91, Sweden<EMAIL_ADDRESS>Axel Brandenburg  Nordita,
Stockholm University and KTH Royal Institute of Technology
Hannes Alfvéns väg 12, SE-106 91 Stockholm, Sweden  The Oskar Klein Centre,
Department of Astronomy, Stockholm University, AlbaNova, Stockholm SE-106 91,
Sweden  Brynmor Haskell  Nicolaus Copernicus Astronomical Centre, Polish
Academy of Sciences, Bartycka 18, 00-716 Warsaw, Poland

###### Abstract

The origin and evolution of magnetic fields of neutron stars from birth has
long been a source of debate. Here, motivated by recent simulations of the
Hall cascade with magnetic helicity, we invoke a model where the large-scale
magnetic field of neutron stars grows as a product of small-scale turbulence
through an inverse cascade. We apply this model to a simulated population of
neutron stars at birth and show how this model can account for the evolution
of such objects across the $P\dot{P}$ diagram, explaining both pulsar and
magnetar observations. Under the assumption that small-scale turbulence is
responsible for large-scale magnetic fields, we place a lower limit on the
spherical harmonic degree of the energy-carrying magnetic eddies of $\approx
40$. Our results favor the presence of a highly resistive pasta layer at the
base of the neutron star crust. We further discuss the implications of this
paradigm on direct observables, such as the nominal age and braking index of
pulsars.

stars:magnetars, pulsars: general, stars: magnetic field

##  1 Introduction

Neutron stars harbor the strongest known magnetic fields in the universe, with
the large-scale poloidal field strength $B_{\rm p}$ in the so-called magnetars
exceeding values of $B_{\rm p}\approx 10^{15}{\rm G}$. The magnetic field
strength plays a crucial role in determining the observational properties of a
neutron star and whether, for example, it will be seen as a standard radio
pulsar or a magnetar (e.g., Borghese, 2023). Furthermore, the large-scale
dipolar field is thought to be primarily responsible for the spindown of
neutron stars (Ostriker & Gunn, 1969), leading to a spindown law of the form
$\dot{\Omega}\approx\Omega^{n}$, where $\Omega=2\pi/P$ is the spin frequency
of the star with $P$ being the rotation period, and $n$ is the so-called
braking index that can be observationally constrained from the combination
$n=\ddot{\Omega}\Omega/\dot{\Omega}^{2}$. For dipole magnetic spindown, one
has $n=3$, with deviations from this value expected if other mechanisms are
driving the spindown (e.g., $n=1$ for winds, and $n=5$ for gravitational wave
emission from quadrupolar ‘mountains’).

In this picture, neutron stars are assumed to be born rapidly rotating, in the
upper left corner of the $P\dot{P}$ diagram, and evolve towards the bottom
right corner and the so-called ‘death’ line as they spin down with a constant
magnetic field. Observed braking indices are, however, often different from
$n=3$, and there is significant observational and theoretical evidence for
magnetic field evolution during the lifetime of a neutron star; see Igoshev et
al. (2021b) for a recent review. Therefore, understanding neutron star
formation and evolution mechanisms is essential to understanding the observed
neutron star population.

The origin of neutron star magnetic fields is, however, still debated. It is
generally agreed that the fossil field inherited from the progenitor star is
insufficient to explain the strongest neutron star fields; therefore, some
field amplification is necessary. Most models invoke a large-scale dynamo at
birth, powered by convection or differential rotation (Thompson & Duncan,
1993). The source of the turbulence driving the dynamo comes from neutrino-
driven convection. An additional source of turbulence is the magnetorotational
instability (Guilet et al., 2022; Reboul-Salze et al., 2022), which draws
energy from the differential rotation, i.e., ultimately from potential energy.
This instability requires the presence of a magnetic field, which is
accomplished through a positive feedback loop whereby magnetic energy can be
amplified further by dynamo action (Brandenburg et al., 1995). Other possible
dynamo mechanisms tap energy through fallback accretion (Barrère et al., 2022)
or through precession (Lander, 2021). In all those cases, small-scale and
large-scale dynamo action appear together. In this Letter, we define a large-
scale field as the part contained in the dipole component, while the small-
scale field is the remainder.

Magnetic field amplification can also occur later in the life of a neutron
star due to the re-emergence of a buried magnetic field due to Ohmic
dissipation (Muslimov & Page, 1996; Ho, 2011) or Hall drift (Gourgouliatos &
Cumming, 2015). Most of these models of magnetic field evolution in neutron
stars assume initially low-order multipoles (Gourgouliatos & Hollerbach, 2018;
Gourgouliatos et al., 2020), although the true magnetic field structure is a
matter of considerable debate (Igoshev et al., 2021b). Alternatively, the
magnetic field could be predominantly of small scales at a very young age of
the neutron star, but the field could then undergo what is known as inverse
cascading (Wareing & Hollerbach, 2009, 2010; Igoshev et al., 2021a). The
spectral magnetic energy at small wavenumbers or low multiples would then
increase with time rather than decrease. This idea was advanced in a recent
paper (Brandenburg, 2020, hereafter B20), where it was shown that the
resulting large-scale field, $B_{\rm LS}$, increases approximately linearly
with time, typically by three orders of magnitude in the models of B20, while
the thermal emission continues to decrease.

---
Figure 1: Evolution of neutron star properties for varying initial conditions
(red curves). From left to right, the top panels show the evolution across the
$P\dot{P}$ diagram, and then the evolution with time of the period, and the
period derivative. The bottom panels show the evolution with time of the
observed braking index, the large-scale magnetic field, and the characteristic
age.

In this Letter, we investigate the implications of the scenario proposed in
B20 on a simulated population of neutron stars and observables such as the
braking index, characteristic age $\tau_{\rm char}$, and the viability of this
idea to explain the neutron star population. We begin in Section 2 discussing
the motivation for a model where the magnetic field grows from a small-scale
turbulent field at birth to a large-scale field. In Section 3, we describe the
spin evolution of a neutron star under this hypothesis and perform simulations
showcasing the evolution of a population of neutron stars for various
observables in Section 4. We discuss implications of our results and conclude
in Section 5.

##  2 Growth of large-scale field

A proto-neutron star is typically fully convective at birth and remains so for
tens of thousands of turnover times (Epstein, 1979), corresponding to a time
of $\mathcal{O}(1\,\mathrm{min})$ (Burrows & Lattimer, 1986). It is commonly
assumed that, within this time, the progenitor magnetic field is destroyed and
then regenerated by a dynamo (Thompson & Duncan, 1993). However, given the
limited time between the collapse of the progenitor star and the end of the
convective period, it is conceivable that only a small-scale magnetic field
with scales comparable to the size of the turbulent eddies is generated. At
the neutron star surface, the scale of such a field can be estimated by the
density scale height $H$, which we expect to coincide roughly with the
thickness of the crust of $\approx 5\%$ of the neutron star radius (Cumming et
al., 2004). This corresponds to a typical wavenumber $k=1/H$, and thus to a
typical spherical harmonics degree of $\ell=(H/R)^{-1}\approx 20$, where $R$
is the radius of the neutron star. Note, however, that the thickness of the
crust will depend on the star’s mass, the details of the equation of state,
and the scale height will be some fraction of this thickness. Therefore, we
will explore a range $10\leq\ell_{\rm max}\leq 40$.

Eventually, when the neutron star solidifies, the magnetic field in the crust
is still randomly distributed in space. It then decays in a way analogous to a
turbulent decay, but here the dynamics are governed by the Hall cascade
(Goldreich & Reisenegger, 1992). In this case, similar to the case of ordinary
magnetohydrodynamic (MHD) turbulence, the magnetic field undergoes an inverse
cascade, where not only the peak of the magnetic energy spectrum moves toward
smaller $k$, but also the spectral energy at the smallest wavenumber increases
in time. Numerical simulations in B20 found the large-scale field $B_{\rm LS}$
to increase with decreasing rms magnetic field strength, $B_{\rm rms}$, like
$B_{\rm LS}\propto B_{\rm rms}^{-5}$; see his Figure 12(b).

In the Hall cascade, $B_{\rm rms}$ decays with time like $B_{\rm rms}\propto
t^{-p/2}$, where $p=2/5$ is the decay exponent of the magnetic energy for a
helical field (B20). Thus,

$B_{\rm LS}\propto B_{\rm rms}^{-5}\propto t,$ (1)

so the large-scale magnetic field increases linearly with time; see Sarin et
al. (2023) for supplemental material showing that the growth of $B_{\rm LS}$
is actually closer to $t^{1.2}$, although this minor departure from the
expected linear growth will be ignored here. This linear growth applies to the
case when the magnetic field is fully helical (see Runs E and F in Figure 12
of B20). However, it also applies to the case of nearly nonhelical fields (see
Runs C and D in Figure 12 of B20) because the magnetic helicity is nearly
unchanged, but the magnetic energy decays, so the ratio increases, making the
magnetic field in the end nearly fully helical; see Tevzadze et al. (2012) for
similar behavior in the MHD case.

The increase of $B_{\rm LS}$ continues until the largest scale in the system
has been reached. Mathematically, the scale increases with $t$ like $\xi_{\rm
M}\propto t^{q}$, where $q=2/5=0.4$ in the helical case, and $4/13\approx
0.31$ in the nonhelical case (Brandenburg, 2023). When $\xi_{\rm M}$ becomes
comparable to $R$, we must adopt global spherical geometry.

So far, inverse cascading has not yet been seen in global simulations. There
are multiple reasons for this. First, it helps to initialize the simulations
with a magnetic field having a broken power law energy spectrum with a proper
inertial range and a sufficiently steep subinertial range. Second, large
numerical resolution is required to obtain inverse cascading, at least in the
nonhelical case (B20 used $1024^{3}$ mesh points). The global simulations of
Gourgouliatos et al. (2020), for example, did not have power laws, extended
only to spherical harmonic degrees of eighty, and did not include magnetic
helicity, which is the case considered here. Dehman et al. (2023) also
presented global simulations with initially complex magnetic field
configurations, but again not with power-law initial fields or with helicity.

In the absence of suitable global simulations, we argue here that we can
substitute the spherical harmonics representation of the Laplacian,
$\ell(\ell+1)/R^{2}$, by just $k^{2}$, as in Cartesian geometry. Simplifying
this further to $k\approx\ell/R$, and using $k\approx\xi_{\rm M}^{-1}\propto
t^{-q}$ with $q=2/5$ in the helical case, we have

$\ell\approx kR\propto Rt^{-q}.$ (2)

As argued above, the initial value of $\ell$ is expected to be somewhere
around $\ell_{\rm max}\approx 20$, and will then decrease to $\ell=1$ like

$\ell=\ell_{\rm max}\,(t/\tau_{1})^{-q},$ (3)

where $\tau_{1}$ marks the start of the decay. Thus, we have $\ell=1$ at
$t\equiv\tau_{2}=\tau_{1}\ell_{\rm max}^{1/q}$.

After that time, the inverse cascade stops, and then even the large-scale
field can only decay. This decay is expected to be either exponential or in a
power-law fashion. We assume here the latter and thus propose

$B_{\rm LS}^{2}=B_{\rm p,0}^{2}\frac{(1+t/\tau_{1})^{2}}{(1+t/\tau_{2})^{m}}$
(4)

for the time evolution of the large-scale magnetic field, where $B_{\rm p,0}$
is the initial large-scale poloidal (dipolar $\ell=1$) component of the
neutron star magnetic field. Here, $m$ dictates the slope of the large-scale
field decay, which is not well constrained. Finally, let us note that the
value of $\tau_{1}$ is linked, in the simulations of B20, to the diffusive
timescale $t_{\rm d}\approx\xi_{\rm M}^{2}/\eta$, where $\eta$ is the magnetic
resistivity. Low values of $\tau_{1}$ generally correspond to a high
resistivity, as we will discuss in the following sections.

##  3 Spin evolution

With the behavior of the magnetic field described above, we can continue to
describe the evolution of the neutron star spin. For the sake of simplicity,
we start with the standard vacuum spin-down formula (Ostriker & Gunn, 1969)
ignoring other spin-down mechanisms such as gravitational-wave emission. Using
Equation (4), our neutron star spin-evolution is then governed by

$\dot{\Omega}=-\frac{B_{\rm
p,0}^{2}R^{6}\Omega(t)^{3}}{6c^{3}I}\frac{(1+t/\tau_{1})^{2}}{(1+t/\tau_{2})^{m}}.$
(5)

Here, $c$ is the speed of light, and $I$ is the moment of inertia of the
neutron star. We emphasize that the above model is effectively the common
vacuum-dipole spin-down equation (Ostriker & Gunn, 1969) except for accounting
of the evolution of the large-scale magnetic field following Eq. (4). Under
the assumption that the magnetic field evolves from small scales to large
scales in a fully helical manner, the timescales $\tau_{1}$ and $\tau_{2}$ can
be related to each other via

$\tau_{2}=\tau_{1}\ell_{\rm max}^{5/2}.$ (6)

Here, $\ell_{\rm max}$ can be estimated as the inverse ratio of the electron
density scale height to the radius of the neutron star, i.e., $(H_{\rm
e}/R)^{-1}$. As discussed above, the magnetic field grows until $\ell\to 1$,
i.e., until the wavenumber approaches a dipole from the initial small-scale
wavenumber $\ell_{\rm max}$. For later reference, values of $\ell_{\rm
max}=40$, 20, and 10 correspond to $\tau_{2}/\tau_{1}\approx 10^{4}$, 1800,
and 300. Owing to the linear scaling of $B_{\rm LS}$ in Eq. (1), the ratios
$\tau_{2}/\tau_{1}$ also correspond to the amplification factor of the large-
scale field.

---
Figure 2: Real $P\dot{P}$ diagram of pulsars obtained from the ATNF catalog
while the middle panel shows a simulated $P\dot{P}$ diagram with our model for
a $1000$ different neutron star initial conditions using
$\tau_{1}=10^{2}\,{\rm yr}$ and $\tau_{2}=10^{6}\,{\rm yr}$, corresponding to
$\ell_{\rm max}=40$. Both panels include lines for constant magnetic fields
and characteristic ages and the pulsar death line, calculated using the psrqpy
software package. The color of the dots in the middle panel indicates the age
of the selected neutron star, with neutron stars older than
$10^{9}\,\mathrm{yr}$ in black, red for neutron stars between
$10^{5}$–$10^{9}\,\mathrm{yr}$, blue for neutron stars between
$10^{2}$–$10^{5}\,\mathrm{yr}$, and green for neutron stars younger than
$10^{2}\,\mathrm{yr}$. The last panel shows the histogram of the large-scale
magnetic field of the $1000$ neutron stars at their selected ages.

##  4 Simulation

With the model for the spin and magnetic field evolution of neutron stars
described above, we now consider the effect of such a model on several neutron
star properties and observables through a series of simulations. For each of
our simulations, we draw the initial spin-period $P_{0}$ and initial magnetic
field $B_{\rm p,0}$ from astrophysically motivated distributions $p$ (Igoshev
et al., 2022; Pagliaro et al., 2023), in particular,

$\displaystyle p(\ln P_{0}/\mathrm{s})$
$\displaystyle=\mathcal{N}(-1.25,0.99),$ $\displaystyle p(\log_{10}B_{\rm
p,0}/\mathrm{G})$ $\displaystyle=\mathcal{U}(8,11),$

where $\mathcal{N}$ and $\mathcal{U}$ denote normal and uniform distributions
in given ranges.

We also fix $m=2.3$ throughout this work, consistent with the decay of the
magnetic field under ohmic dissipation (e.g., Pons & Geppert, 2007). We note
here that a precise value of $m$ is not important to this work as we focus on
the growth of the magnetic field rather than the decay. We have verified that
the results are robust to our choices of $m$, the initial spin, and large-
scale magnetic field. Our motivation for choosing these ranges were to show 1)
that one does not need large scale magnetic fields at birth to make magnetars
and 2) that a regular neutron star spin population can make magnetars.

We start by investigating the evolution with time; we set
$\tau_{1}=10^{2}\,\mathrm{yr}$ to $\tau_{2}=10^{6}\,\mathrm{yr}$, i.e.,
$\ell_{\rm max}\approx 40$, and evolve $100$ neutron stars drawn with initial
periods and magnetic fields from the distributions above forward in time.
Figure 1 shows the evolution for different observables. In particular, we show
the evolution of neutron stars across the $P\dot{P}$ diagram, the period and
period derivative as functions of time, the observable braking index, dipole
component of the magnetic field and the characteristic age. Figure 1 helps
illustrate the model’s main features. For example, large-scale initial fields
as small as $10^{11}\,\mathrm{G}$ can be amplified to magnetar strengths
($10^{14}\,\mathrm{G}$, alleviating the concern of needing a large-scale
poloidal field at birth. In particular, the age at such field strengths
depends on the choice of $\tau_{1}$ and $\ell_{\rm max}$. Similarly, the model
can naturally produce a range of braking indices despite being built only on
the assumption of vacuum-dipole emission, and naturally returns values of the
braking index $n\lesssim 2$ for young pulsars, as measured in the observed
population (Espinoza et al., 2017). We also show the evolution of the
characteristic age with time. We find that at young ages, there is a huge
discrepancy between the characteristic age and the actual age of the neutron
star, while after some time these ages tend to agree. This is expected as our
spin evolution is only different from $n=3$ at early times. The evolution
across the $P\dot{P}$ diagram also demonstrates that the model can accommodate
neutron stars in various parts of the $P\dot{P}$ diagram despite a narrow
starting position and explain the different locations we observe for real
neutron stars.

---
Figure 3: Simulated $P\dot{P}$ diagrams for the same initial conditions in all
panels and ages but different $\tau_{1}$ and $\ell_{\rm max}$. The colors
correspond to the ages of the neutron stars, with the same ranges as described
in Fig. 2.

To illustrate this last point, we now simulate a hypothetical $P\dot{P}$
diagram with $\tau_{1}=1\,\mathrm{yr}$ and $\ell_{\rm max}=40$ corresponding
to $\tau_{2}=10^{4}\,\mathrm{yr}$, ignoring selection effects and drawing ages
randomly throughout the lifespan of the neutron star. In Figure 2, the left
panel shows a real $P\dot{P}$ diagram of pulsars obtained from the ATNF
catalog (Manchester et al., 2005) through the psrqpy software package (Pitkin,
2018). The middle panel shows a $1000$ simulated neutron stars with our model
drawn randomly from the above distributions and selected at different ages.
Both panels include lines for constant magnetic fields and characteristic ages
and the pulsar death line (Zhang et al., 2000), calculated through psrqpy with
the region below shaded in yellow. The last panel shows the histogram of the
large-scale magnetic field of the $1000$ simulated neutron stars at their
selected ages. We emphasize that we do not model processes such as mass
transfer that is likely responsible for neutron stars observed in the lower
left of the $P\dot{P}$ diagram, i.e., the millisecond pulsars. We caution the
reader against a direct comparison of the observed population and our
simulated population. We do not model the radio luminosity of any of the
neutron stars and therefore cannot account for any selection effects that may
limit our ability to see neutron stars in certain parts of the $P\dot{P}$
diagram (Faucher-Giguère & Kaspi, 2006; Szary et al., 2014).

Figure 2 also demonstrates that this evolutionary model for the magnetic field
can produce neutron stars in the “magnetar” region of the $P\dot{P}$ diagram
while also accounting for the distribution of neutron stars across the
$P\dot{P}$ diagram for a small range of initial conditions. Furthermore, it
helps illustrate the significant differences in neutron stars’ actual ages and
magnetic fields compared to the characteristic spin-down age and constant
magnetic-field models. This may be important to reconcile the often
significant discrepancies in inferred ages of neutron stars in supernova
remnants with their location in the $P\dot{P}$ diagram.

The above $P\dot{P}$ diagram represented one single case of $\tau_{1}$ and
$\ell_{\rm max}$. However, as we mentioned in Section 2, these values are
unknown. In Figure 3, we show simulated $P\dot{P}$ diagrams for a range of
these values. We note that each panel has the same initial conditions (i.e.,
$P_{0}$ and $B_{\rm p,0}$) and all neutron stars are selected at the same
ages, and the difference in their “observed” locations is only due to
different $\tau_{1}$ and $\ell_{\rm max}$. Figure 3 can be interpreted in two
primary ways. First, for almost any choice of $\tau_{1}$ and $\ell_{\rm max}$,
this model can accommodate a large portion of the $P\dot{P}$ diagram for a
small range of initial conditions and provide a natural explanation for the
often discrepant age estimates from spindown and a supernova remnant. Second,
to explain magnetars with this model, we need large $\ell_{\rm max}$ and small
$\tau_{1}$, i.e., we need the inverse cascades to begin early and need the
initial field to be confined to small scales.

##  5 Conclusions

Our results demonstrate that it is indeed possible to explain the production
of large-scale magnetic fields through an inverse cascade. This means that an
initial small-scale magnetic field gets gradually converted into a large-scale
one. The efficiency of this process increases with increasing dynamic range in
space and time, i.e., it becomes more efficient, the smaller the scale of the
initial field (larger $\ell_{\rm max}$), and the earlier the inverse cascade
begins (smaller $\tau_{1}$).

Our models for the growth of the magnetic field early on in the life of the
star can explain the low ($n\lesssim 2$) braking indices measured in young
pulsars, and their observed evolution in the $P\dot{P}$ diagram (Espinoza et
al., 2017). This is consistent with numerical results presented in
Gourgouliatos & Cumming (2015), where the low braking indices of Vela and
other young pulsars were interpreted as a byproduct of the increase in large-
scale magnetic fields in Hall MHD. However, note that Gourgouliatos & Cumming
(2015), started with initially strong large-scale magnetic fields, while we
can reproduce magnetars with much smaller initial large-scale fields.

The most extreme case studied in this Letter is $\ell_{\rm max}=40$ and
$\tau_{1}=1\,{\rm yr}$. We estimated $\ell_{\rm max}=20$ based on the scale
height, but this provided a rough estimate, because the convective eddies
during the first minute of the neutron star’s life can well be smaller than
the local scale height as discussed in Section 2. Somewhat larger values of
$\ell_{\rm max}=40$ appear therefore feasible. However, the situation with the
value of $\tau_{1}$ is less obvious. Once the crust has been established, the
relevant timescale is the magnetic diffusion time, which was also used as the
natural time unit in B20 and Brandenburg (2023). In their Cartesian models, it
was defined as $\tau_{\rm d}=\xi_{\rm M}^{2}/\eta$. Using $\eta=4\times
10^{-8}\,{\rm m}^{2}\,{\rm s}^{-1}$ and $\xi_{\rm M}=R/\ell_{\rm max}\approx
700\,{\rm m}$ with $\ell_{\rm max}=20$, we have $\tau_{\rm
d}=0.4\,\mathrm{Myr}$. Somewhat smaller values are expected as we increase
$\ell$, for example to 100, in which case we can have $\xi_{\rm M}=150\,{\rm
m}$, and therefore $\tau_{\rm d}=20\,\mathrm{kyr}$. In addition, the time
$\tau_{1}$ when the inverse cascade begins is a certain fraction of $\tau_{\rm
d}$, for example $\tau_{1}/\tau_{\rm d}=0.01$, corresponds to the time when a
certain characteristic quantity (the Hosking integral) reached a maximum; see
Figure 4 of Brandenburg (2023). Thus, the smallest conceivable value of
$\tau_{1}$ is 200 yr for a standard crustal composition. It should be noted,
however, that recent calculations of the conductivity in the presence of a
nuclear pasta phase (Pelicer et al., 2023) have yielded significantly smaller
conductivities, which could accommodate values of $\tau_{1}\approx 1$ yr.

Our results show that low values of $\tau_{1}$ (together with high values of
$\ell_{\rm max}$) will create stars in the magnetar range in most of our
models, and therefore point towards the presence of a thin, highly resistive
pasta layer at the base of the crust, as suggested also by Pons et al. (2013),
to explain the absence of long period isolated X-ray pulsars. We note that all
models produce neutron stars below the death line, but those would not be
observable as pulsars owing to their low luminosity.

Another implication of this model is on the ages of magnetars. If the magnetic
field is confined to small-scales at birth and the initial dipole component is
small, then there needs to be a minimum time such that the field can grow to
magnetar strengths. For the choice of $\tau_{1}=10\,\mathrm{yr}$ and
$\ell_{\rm max}=40$, this is $\geq 500\,\mathrm{yr}$ for our simulation, i.e.,
there should not exist any neutron stars with $B_{\rm p}\gtrsim
10^{14}\,\mathrm{G}$, younger than $\approx 500\,\mathrm{yr}$. A comparison
between characteristic age and true age is presented in Figure 1, where we see
that at early times most of our systems are much younger than they appear (we
caution, however, that the systems we produce with large characteristic ages
have very small values of $\dot{P}$ at birth, which would make them
challenging to observe), while in the range between $10^{3}$ and $10^{5}$
years, a small number of systems, which correspond to the magnetars, are
slightly older than they appear. However, we note that we do not capture the
true diversity of initial conditions such as spin or considerations of the
equation of state which could produce magnetar-like field strength earlier in
time. We will explore the full effect of different initial conditions,
equations of state, and compare directly to observations in future work.

##  6 Acknowledgments

We are grateful to Clara Dehman, Cristobal Espinoza, and Kostas Gourgouliatos
for helpful comments. We also thank the referee for their helpful feedback
that has improved the presentation of this work. N.S. is supported by a
Nordita fellowship. A.B. acknowledges support from the Swedish Research
Council (Vetenskapsrådet, grant number 2019-04234). Nordita is supported in
part by NordForsk. B.H. acknowledges support from the National Science Centre
Poland (NCN) via grant OPUS 2018/29/B/ST9/02013.

## References

  * Barrère et al. (2022) Barrère, P., Guilet, J., Reboul-Salze, A., Raynaud, R., & Janka, H. T. 2022, A&A, 668, A79, doi: 10.1051/0004-6361/202244172
  * Borghese (2023) Borghese, A. 2023, IAUS, 363, 51, doi: 10.1017/S1743921322000357
  * Brandenburg (2020) Brandenburg, A. 2020, ApJ, 901, 18, doi: 10.3847/1538-4357/abad92
  * Brandenburg (2023) —. 2023, JPlPh, 89, 175890101, doi: 10.1017/S0022377823000028
  * Brandenburg et al. (1995) Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741, doi: 10.1086/175831
  * Burrows & Lattimer (1986) Burrows, A., & Lattimer, J. M. 1986, ApJ, 307, 178, doi: 10.1086/164405
  * Cumming et al. (2004) Cumming, A., Arras, P., & Zweibel, E. 2004, ApJ, 609, 999, doi: 10.1086/421324
  * Dehman et al. (2023) Dehman, C., Viganò, D., Ascenzi, S., Pons, J. A., & Rea, N. 2023, MNRAS, 523, 5198, doi: 10.1093/mnras/stad1773
  * Epstein (1979) Epstein, R. I. 1979, MNRAS, 188, 305, doi: 10.1093/mnras/188.2.305
  * Espinoza et al. (2017) Espinoza, C. M., Lyne, A. G., & Stappers, B. W. 2017, MNRAS, 466, 147, doi: 10.1093/mnras/stw3081
  * Faucher-Giguère & Kaspi (2006) Faucher-Giguère, C.-A., & Kaspi, V. M. 2006, ApJ, 643, 332, doi: 10.1086/501516
  * Goldreich & Reisenegger (1992) Goldreich, P., & Reisenegger, A. 1992, ApJ, 395, 250, doi: 10.1086/171646
  * Gourgouliatos & Cumming (2015) Gourgouliatos, K. N., & Cumming, A. 2015, MNRAS, 446, 1121, doi: 10.1093/mnras/stu2140
  * Gourgouliatos & Hollerbach (2018) Gourgouliatos, K. N., & Hollerbach, R. 2018, ApJ, 852, 21, doi: 10.3847/1538-4357/aa9d93
  * Gourgouliatos et al. (2020) Gourgouliatos, K. N., Hollerbach, R., & Igoshev, A. P. 2020, MNRAS, 495, 1692, doi: 10.1093/mnras/staa1295
  * Guilet et al. (2022) Guilet, J., Reboul-Salze, A., Raynaud, R., Bugli, M., & Gallet, B. 2022, MNRAS, 516, 4346, doi: 10.1093/mnras/stac2499
  * Ho (2011) Ho, W. C. G. 2011, MNRAS, 414, 2567, doi: 10.1111/j.1365-2966.2011.18576.x
  * Igoshev et al. (2022) Igoshev, A. P., Frantsuzova, A., Gourgouliatos, K. N., Tsichli, S., & et al. 2022, MNRAS, 514, 4606, doi: 10.1093/mnras/stac1648
  * Igoshev et al. (2021a) Igoshev, A. P., Gourgouliatos, K. N., Hollerbach, R., & Wood, T. S. 2021a, ApJ, 909, 101, doi: 10.3847/1538-4357/abde3e
  * Igoshev et al. (2021b) Igoshev, A. P., Popov, S. B., & Hollerbach, R. 2021b, Univ, 7, 351, doi: 10.3390/universe7090351
  * Lander (2021) Lander, S. K. 2021, MNRAS, 507, L36, doi: 10.1093/mnrasl/slab086
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993, doi: 10.1086/428488
  * Muslimov & Page (1996) Muslimov, A., & Page, D. 1996, ApJ, 458, 347, doi: 10.1086/176817
  * Ostriker & Gunn (1969) Ostriker, J. P., & Gunn, J. E. 1969, ApJ, 157, 1395, doi: 10.1086/150160
  * Pagliaro et al. (2023) Pagliaro, G., Alessandra Papa, M., Ming, J., Lian, J., & et al. 2023, arXiv e-prints, arXiv:2303.04714, doi: 10.48550/arXiv.2303.04714
  * Pelicer et al. (2023) Pelicer, M. R., Antonelli, M., Menezes, D. P., & Gulminelli, F. 2023, MNRAS, 521, 743, doi: 10.1093/mnras/stad562
  * Pitkin (2018) Pitkin, M. 2018, JOSS, 3, 538, doi: 10.21105/joss.00538
  * Pons & Geppert (2007) Pons, J. A., & Geppert, U. 2007, A&A, 470, 303, doi: 10.1051/0004-6361:20077456
  * Pons et al. (2013) Pons, J. A., Viganò, D., & Rea, N. 2013, NatPh, 9, 431, doi: 10.1038/nphys2640
  * Reboul-Salze et al. (2022) Reboul-Salze, A., Guilet, J., Raynaud, R., & Bugli, M. 2022, A&A, 667, A94, doi: 10.1051/0004-6361/202142368
  * Sarin et al. (2023) Sarin, N., Brandenburg, A., & Haskell, B. 2023, Zenodo, doi: 10.5281/zenodo.8088084
  * Szary et al. (2014) Szary, A., Zhang, B., Melikidze, G. I., Gil, J., & Xu, R.-X. 2014, ApJ, 784, 59, doi: 10.1088/0004-637X/784/1/59
  * Tevzadze et al. (2012) Tevzadze, A. G., Kisslinger, L., Brandenburg, A., & Kahniashvili, T. 2012, ApJ, 759, 54, doi: 10.1088/0004-637X/759/1/54
  * Thompson & Duncan (1993) Thompson, C., & Duncan, R. C. 1993, ApJ, 408, 194, doi: 10.1086/172580
  * Wareing & Hollerbach (2009) Wareing, C. J., & Hollerbach, R. 2009, PhPl, 16, 042307, doi: 10.1063/1.3111033
  * Wareing & Hollerbach (2010) —. 2010, JPlPh, 76, 117, doi: 10.1017/S0022377809990158
  * Zhang et al. (2000) Zhang, B., Harding, A. K., & Muslimov, A. G. 2000, ApJ, 531, L135, doi: 10.1086/312542