# Longtime Dynamics of Irrotational Spherical Water Drops: Initial Notes

Chengyang Shao

In this note, we propose several unsolved problems concerning the irrotational
oscillation of a water droplet under zero gravity. We will derive the
governing equation of this physical model, and convert it to a quasilinear
dispersive partial differential equation defined on the sphere, which formally
resembles the capillary water waves equation but describes oscillation defined
on curved manifold instead. Three types of unsolved mathematical problems
related to this model will be discussed in observation of hydrodynamical
experiments under zero gravity111There are numbers of visual materials on such
experiments conducted by astronauts. See for example
https://www.youtube.com/watch?v=H_qPWZbxFl8&t or
https://www.youtube.com/watch?v=e6Faq1AmISI&t.: (1) Strichartz type
inequalities for the linearized problem (2) existence of periodic solutons (3)
normal form reduction and generic lifespan estimate. It is pointed out that
all of these problems are closely related to certain Diophantine equations,
especially the third one.

##  1\. Capillary Spherical Water Waves Equation: Derivation

###  1.1. Water Waves Equation for a Bounded Water Drop

Comparing to gravity water waves problems, the governing equation for a
spherical droplet of water under _zero gravity_ takes a very different form.
At a first glance it looks similar to the water waves systems as mentioned
above, but some crucial differences do arise after careful analysis. To the
author’s knowledge, besides those dealing with generic free-boundary Euler
equation ([CS07], [SZ08]), the only reference on this problem is Beyer-Günther
[BG98], in which the local well-posedness of the equation is proved using a
Nash-Moser type implicit function theorem. We will briefly discribe known
results for gravity water waves problems in the next subsection.

To start with, let us pose the following assumptions on the fluid motion that
we try to describe:

  * •

(A1) The perfect, irrotational fluid of constant density $\rho_{0}$ occupies a
smooth, compact region in $\mathbb{R}^{3}$.

  * •

(A2) There is no gravity or any other external force in presence.

  * •

(A3) The air-fluid interface is governed by the Young-Laplace law, and the
effect of air flow is neglected.

We assume that the boundary of the fluid region has the topological type of a
smooth compact orientable surface $M$, and is described by a time-dependent
embedding $\iota(t,\cdot):M\to\mathbb{R}^{3}$. We will denote a point on $M$
by $x$, the image of $M$ under $\iota(t,\cdot)$ by $M_{t}$, and the region
enclosed by $M_{t}$ by $\Omega_{t}$. The outer normal will be denoted by
$N(\iota)$. We also write $\bar{\nabla}$ for the flat connection on
$\mathbb{R}^{3}$.

Adopting assumption (A3), we have the Young-Laplace equation:

$\sigma_{0}H(\iota)=p_{i}-p_{e},$

where $H(\iota)$ is the (scalar) mean curvature of the embedding, $\sigma_{0}$
is the surface tension coefficient (which is assumed to be a constant), and
$p_{i},p_{e}$ are respectively the inner and exterior air pressure at the
boundary; they are scalar functions on the boundary and we assume that $p_{e}$
is a constant. Under assumptions (A1) and (A2), we obtain Bernoulli’s
equation, sometimes referred as the pressure balance condition, on the
evolving surface:

(1.1) $\left.\frac{\partial\Phi}{\partial
t}\right|_{M_{t}}+\frac{1}{2}|\bar{\nabla}\Phi|_{M_{t}}|^{2}-p_{e}=-\frac{\sigma_{0}}{\rho_{0}}H(\iota),$

where $\Phi$ is the velocity potential of the velocity field of the air. Note
that $\Phi$ is determined up to a function in $t$, so we shall leave the
constant $p_{e}$ around for convenience reason that will be explained shortly.
According to assumption (A1), the function $\Phi$ is a harmonic function
within the region $\Omega_{t}$, so it is uniquely determined by its boundary
value, and the velocity field within $\Omega_{t}$ is $\bar{\nabla}\Phi$. The
kinematic equation on the free boundary $M_{t}$ is naturally obtained as

(1.2) $\frac{\partial\iota}{\partial t}\cdot
N(\iota)=\bar{\nabla}\Phi|_{M_{t}}\cdot N(\iota).$

Finally, we would like to discuss the conservation laws for (1.1)-(1.2). The
preservation of volume $\text{Vol}(\Omega_{t})=\text{Vol}(\Omega_{0})$ is a
consequence of incompressibility. The system describes an Eulerian flow
without any external force, so the center of mass moves at a uniform speed
along a fixed direction, i.e.

(1.3)
$\frac{1}{\mathrm{Vol}(\Omega_{0})}\int_{\Omega_{t}}Pd\mathrm{Vol}(P)=V_{0}t+C_{0},$

with Vol being the Lebesgue measure, $P$ marking points in $\mathbb{R}^{3}$,
$V_{0}$ and $C_{0}$ being the velocity and starting position of center of mass
respectively. Furthermore, the total momentum is conserved, and since the flow
is a potential incompressible one, the conservation of total momentum is
expressed as

(1.4) $\int_{M_{t}}\rho_{0}\Phi
N(\iota)d\mathrm{Area}(M_{t})\equiv\rho_{0}\mathrm{Vol}(\Omega_{0})V_{0}.$

Most importantly, it is not surprising that (1.1)-(1.2) is a Hamilton system
(the Zakharov formulation for water waves; see Zakharov [Zak68]), with
Hamiltonian

(1.5)
$\sigma_{0}\mathrm{Area}(\iota)+\frac{1}{2}\int_{\Omega_{t}}\rho_{0}|\bar{\nabla}\Phi|^{2}d\mathrm{Vol}=\sigma_{0}\mathrm{Area}(M_{t})+\frac{1}{2}\int_{M_{t}}\rho_{0}\Phi|_{M_{t}}\left(\bar{\nabla}\Phi|_{M_{t}}\cdot
N(\iota)\right)d\mathrm{Area},$

i.e. potential proportional to surface area plus kinetic energy of the fluid.

###  1.2. Converting to a Differential System

It is not hard to verify that the system (1.1)-(1.2) is invariant if $\iota$
is composed with a diffeomorphism of $M$; we may thus regard it as a
_geometric flow_. If we are only interested in perturbation near a given
configuration, we may reduce system (1.1)-(1.2) to a non-degenerate dispersive
differential system concerning two scalar functions defined on $M$, just as
Beyer and Günther did in [BG98]. In fact, during a short time of evolution,
the interface can be represented as the graph of a function defined on the
initial surface: if $\iota_{0}:M\to\mathbb{R}^{3}$ is a fixed embedding close
to the initial embedding $\iota(0,x)$, we may assume that
$\iota(t,x)=\iota_{0}(x)+\zeta(t,x)N_{0}(x)$, where $\zeta$ is a scalar
“height” function defined on $M_{0}$ and $N_{0}$ is the outer normal vector
field of $M_{0}$.

Figure 1. The shape of the surface

With this observation, we shall transform the system (1.1)&(1.2) into a non-
local system of two real scalar functions $(\zeta,\phi)$ defined on $M$, where
$\zeta$ is the “height” function described as above, and
$\phi(t,x)=\Phi(t,\iota(t,x))$ is the boundary value of the velocity
potential, pulled back to the underlying manifold $M$.

The operator

$B_{\zeta}:\phi\to\bar{\nabla}\Phi|_{M_{t}}$

maps the pulled-back Dirichlet boundary value $\phi$ to the boundary value of
the gradient of $\Phi$. We shall write

$(D[\zeta]\phi)N(\iota)$

for its normal part, where $D[\zeta]$ is the Dirichlet-Neumann operator
corresponding to the region enclosed by the image of $\iota_{0}+\zeta N_{0}$.
Thus

$\frac{\partial\zeta}{\partial t}N_{0}\cdot N(\iota)=D[\zeta]\phi.$

We also need to calculate the restriction of $\partial_{t}\Phi$ on $M_{t}$ in
terms of $\phi$ and $\iota$. By the chain rule,

$\displaystyle\left.\frac{\partial\Phi}{\partial t}\right|_{M_{t}}$
$\displaystyle=\frac{\partial\phi}{\partial
t}-\bar{\nabla}\Phi|_{M_{t}}\cdot\frac{\partial\iota}{\partial t}$
$\displaystyle=\frac{\partial\phi}{\partial
t}-\left(\bar{\nabla}\Phi|_{M_{t}}\cdot
N_{0}\right)\frac{\partial\zeta}{\partial t}$
$\displaystyle=\frac{\partial\phi}{\partial t}-\frac{1}{N_{0}\cdot
N(\iota)}\left(\bar{\nabla}\Phi|_{M_{t}}\cdot N_{0}\right)\cdot D[\zeta]\phi.$

We thus arrive at the following nonlinear system:

(EQ(M)) $\left\\{\begin{aligned} \frac{\partial\zeta}{\partial
t}&=\frac{1}{N_{0}\cdot N(\iota)}D[\zeta]\phi,\\\ \frac{\partial\phi}{\partial
t}&=\frac{1}{N_{0}\cdot N(\iota)}\left(B_{\zeta}\phi\cdot N_{0}\right)\cdot
D[\zeta]\phi-\frac{1}{2}|B_{\zeta}\phi|^{2}-\frac{\sigma_{0}}{\rho_{0}}H(\iota)+p_{e},\end{aligned}\right.$

where $\iota=\iota_{0}+\zeta N_{0}$.

######  Remark 1.

_We may obtain an explicit expression of
$B_{\zeta}\phi=\bar{\nabla}\Phi|_{M_{t}}$ in terms of $\phi$ (together with
the connection $\nabla_{0}$ on the fixed embedding $\iota_{0}(M)$), just as
standard references did for Euclidean or periodic water waves, but that is not
necessary for our discussion at the moment. It is important to keep in mind
that the preservation of volume and conservation of total momentum (1.3)-(1.4)
convert to integral equalities of $(\zeta,\phi)$. These additional
restrictions are not obvious from the differential equations (EQ(M)), though
they can be deduced from (EQ(M)) since they are just rephrase of the original
physical laws (1.1)-(1.2). _

For $M=S^{2}$, the case that we shall discuss in detail, we refer to the
system as the _capillary spherical water waves equation_. To simplify our
discussion, we shall be working under the center of mass frame, and require
the eigenmode $\Pi^{(0)}\phi$ to vanish for all $t$. This could be easily
accomplished by absorbing the eigenmode into $\phi$ since the equation is
invariant by a shift of $\phi$. In a word, from now on, we will be focusing on
the non-dimensional capillary spherical water waves equation

(EQ) $\left\\{\begin{aligned} \frac{\partial\zeta}{\partial
t}&=\frac{1}{N_{0}\cdot N(\iota)}D[\zeta]\phi,\\\ \frac{\partial\phi}{\partial
t}&=\frac{1}{N_{0}\cdot N(\iota)}\left(B_{\zeta}\phi\cdot N_{0}\right)\cdot
D[\zeta]\phi-\frac{1}{2}|B_{\zeta}\phi|^{2}-H(\iota),\end{aligned}\right.$

where $\iota=(1+\zeta)\iota_{0}$, and $\Pi^{(0)}\phi\equiv 0$. We assume that
the total volume of the fluid is $4\pi/3$, so that the preservation of volume
is expressed as

(1.6) $\frac{1}{3}\int_{S^{2}}(1+\zeta)^{3}d\mu_{0}\equiv\frac{4\pi}{3},$

where $\mu_{0}$ is the standard measure on $S^{2}$. The inertial movement of
center of mass (1.3) and conservation of total momentum (1.4) under our center
of mass frame are expressed respectively as

(1.7) $\int_{S^{2}}(1+\zeta)^{4}N_{0}d\mu_{0}=0,\quad\int_{S^{2}}\phi
N(\iota)d\mu(\iota)=0,$

where $\mu(\iota)$ is the induced surface measure. Further, the Hamiltonian of
the system is

(1.8)
$\mathbf{H}[\zeta,\phi]=\mathrm{Area}(\iota)+\frac{1}{2}\int_{S^{2}}\phi\cdot
D[\zeta]\phi\cdot d\mu(\iota),$

and for a solution $(\zeta,\phi)$ there holds $\mathbf{H}[\zeta,\phi]\equiv
4\pi$.

Up to this point, we are still working within the realm of well-established
frameworks. We already know that the general free-boundary Euler equation is
locally well-posed due to the work of [CS07] or [SZ08], and due to the curl
equation the curl free condition persists during the evolution. On the other
hand, the Cauchy problem of system (1.1) and (1.2) is known to be locally
well-posed, due to Beyer and Günther in [BG98]. They used an iteration
argument very similar to a Nash-Moser type argument in the sense that it
involves multiple scales of Banach spaces and “tame” maps in a certain sense.
Finally, it is not hard to transplant the potential-theoretic argument of Wu
[Wu99] to prove the local well-posedness.

To sum up, we already know that the system (EQ(M)) for a compact orientable
surface $M$ (hence (EQ) specifically) is locally well-posed. But this is all
we can assert for the motion of a water droplet under zero gravity. In the
following part of this note, we will propose several questions and conjectures
concerning the long-time behaviour of water droplets under zero gravity.

###  1.3. Previous Works on Water Waves

There has already been several different approaches to describe the motion of
perfect fluid with free boundary. One is to consider the motion of perfect
fluid occupying an arbitrary domain in $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$,
either with or without surface tension. The motion is described by a free
boundary value problem of Euler equation. This generic approach was employed
by Countand-Shkoller [CS07] and Shatah-Zeng [SZ08]. Both groups proved the
local-wellposedness of the problem. This approach has the advantage of being
very general, applicable to all geometric shapes of the fluid.

On the other hand, when coming to potential flows of a perfect fluid, the curl
free property results in _dispersive_ nature of the problem. The motion of a
curl free perfect fluid under gravity and a free boundary value condition is
usually referred to as _the gravity water waves problem_. The first
breakthroughs in understanding local well-posedness were works of Wu [Wu97]
[Wu99], who proved local well-posedness of the gravity water waves equation
without any smallness assumption. Lannes [Lan05] extended this to more generic
bottom shapes. Taking surface tension into account, the problem becomes
_gravity-capillary water waves_. Schweizer [Sch05] proved local well-posedness
with small Cauchy data of the gravity-capillary water waves problem, and Ming-
Zhang [MZ09] proved local well-posedness without smallness assumption.
Alazard-Metevier [AM09] and Alazard-Burq-Zuily [ABZ11] used para-differential
calculus to obtain the optimal regularity for local well-posedness of the
water waves equation, either with or without surface tension.

For discussion of long time behavior, it is important to take into account the
dispersive nature of the problem. For linear dispersive properties, there has
been work of Christianson-Hur-Staffilani [CHS10]. For gravity water waves
living in $\mathbb{R}^{2}$, works on lifespan estimate include Wu [Wu09] and
Hunter-Ifrim-Tataru [HIT16] (almost global result), Ionescu-Pusateri [IP15]
and Alazard-Delort [AD13] and Ifrim-Tataru [IT14] (global result). For gravity
water waves living in $\mathbb{R}^{3}$, there are works of Germain-Masmoudi-
Shatah [GMS12] and Wu [Wu11] (no surface tension), Germain-Masmoudi-Shatah
[GMS15] (no gravity), Deng-Ionescu-Pausader-Pusateri [DIPP17] (gravity-
capillary water waves) and Wang [Wan20] (gravity-capillary water waves with
finite depth). These results all employed different forms of decay estimates
derived from dispersive properties.

As for long time behavior of periodic water waves, Berti-Delort [BD+18]
considered gravity-capillary water waves defined on $\mathbb{T}^{1}$, Berti-
Feola-Pusateri [BFP22] considered gravity water waves defined on
$\mathbb{T}^{1}$, Ionescu-Pusateri [IP19] considered gravity-capillary water
waves defined on $\mathbb{T}^{2}$, and obtained an estimate on the lifespan
beyond standard energy method. All of the three groups used para-differential
calculus and suitable normal form reduction; the results of Berti-Delort and
Ionescu-Pusateri were proved for physical data of full Lebesgue measure.

To sum up, all results on the gravity water waves problem listed above are
concerned with an equation for two scalar functions defined on a fixed flat
manifold, being one of the following: $\mathbb{R}^{1}$, $\mathbb{R}^{2}$,
$\mathbb{T}^{1}$, $\mathbb{T}^{2}$, sometimes called the “bottom” of the
fluid. These two functions represents the geometry of the liquid-gas interface
and the boundary value of the velocity potential, respectively. The manifold
itself is considered as the bottom of the container in which all dynamics are
performed. We observe that the differential equation (EQ(M)) is,
mathematically, fundamentally different from the water waves equations that
have been well-studied.

##  2\. Initial Notes on Unsolved Problems

In this section, we propose unsolved problems related to the spherical
capillary water waves system (EQ(M)) with $M=S^{2}$. Not surprisingly, these
problems all have deep backgrounds in number theory.

###  2.1. Linearization Around the Static Solution

We are mostly interested in the stability of the static solution of (EQ(M)). A
static solution should be a fluid region whose shape stays still, with motion
being a mere shift within the space. In this case, we have $p_{e}=0$ since the
reference is relatively static with respect to the air. Moreover, the velocity
field $\bar{\nabla}\Phi$ and “potential of acceleration”
$\partial\Phi/\partial t$ must both be spatially uniform, so that the left-
hand-side of the pressure balance condition (1.1) is a function of $t$ alone.
It follows that $\iota(M)$ is always a compact embedded surface of constant
mean curvature, hence in fact always an Euclidean sphere by the Alexandrov
sphere theorem (see [MP19]), and we may just take $M=S^{2}$. Moreover, since
$\iota(S^{2})$ should enclose a constant volume by incompressibility, the
radius of that sphere does not change. After suitable scaling, we may assume
that the radius is always 1, and $\rho_{0}=1$, $\sigma_{0}=1$, to make
(1.1)-(1.2) non-dimensional. Finally, by choosing the center of mass frame, we
may simply assume that the spatial shift is always zero, so that the velocity
potential $\Phi\equiv\Phi_{0}$, a real constant. It is harmless to fix it to
be zero.

Thus, under our convention, a static solution of (1.1)-(1.2) takes the form

(2.1) $\left(\begin{array}[]{c}\iota(t,x)\\\
\Phi(t,x)\end{array}\right)=\left(\begin{array}[]{c}\iota_{0}(x)\\\
\displaystyle{0}\end{array}\right),$

where $a\in\mathbb{R}^{3}$ is a constant vector, and $\iota_{0}$ is the
standard embedding of $S^{2}$ as $\partial B(0,1)\subset\mathbb{R}^{3}$.
Equivalently, this means that a static solution of (EQ(M)) under our
convention must be $(\zeta,\phi)=(0,0)$. Note here that the Gauss map of
$\iota_{0}$ coincides with itself.

We can now start our perturbation analysis around a static solution at the
linear level. Let $\mathcal{E}^{(n)}$ be the space of spherical harmonics of
order $n$, normalized according to the standard surface measure on $S^{2}$. In
particular, $\mathcal{E}^{(1)}$ is spanned by three components of $N_{0}$. Let
$\Pi^{(n)}$ be the orthogonal projection on $L^{2}(S^{2})$ onto
$\mathcal{E}^{(n)}$, $\Pi_{\leq n}$ be the orthogonal projection on
$L^{2}(S^{2})$ onto $\bigoplus_{k\leq n}\mathcal{E}^{(k)}$, $\Pi_{\geq n}$ be
the orthogonal projection on $L^{2}(S^{2})$ onto $\bigoplus_{k\geq
n}\mathcal{E}^{(k)}$. For $\iota=(1+\zeta)\iota_{0}$, the linearization of
$-H(\iota)$ around the sphere $\zeta\equiv 0$ is $\Delta\zeta+2\zeta$, where
$\Delta$ is the Laplacian on the sphere $S^{2}$; cf. the standard formula for
the second variation of area in [BC12]. Then $H^{\prime}(\iota_{0})$ acts on
$\mathcal{E}^{(n)}$ as the multiplier $-(n-1)(n+2)$. Note that even if we
consider the dimensional form of (EQ), there will only be an additional
scaling factor $\sigma_{0}/(\rho_{0}R^{2})$, where $R$ is the radius of the
sphere. On the other hand, the following solution formula for the Dirichlet
problem on $B(0,1)$ is well-known: if $f\in L^{2}(S^{2})$, then the harmonic
function in $B(0,1)$ with Dirichlet boundary value $f$ is determined by

$u(r,\omega)=\sum_{n\geq 0}r^{n}(\Pi^{(n)}f)(\omega),$

where $(r,\omega)$ is the spherical coordinate in $\mathbb{R}^{3}$. Thus the
Dirichlet-Neumann operator $D[\iota_{0}]$ acts on $\mathcal{E}^{(n)}$ as the
multiplier $n$. Note again that even if we consider the dimensional form (EQ),
there will only be an additional scaling factor $R^{-1}$.

Thus, setting

$u=\Pi^{(0)}\zeta+\Pi^{(1)}\zeta+\sum_{n\geq
2}\sqrt{(n-1)(n+2)}\cdot\Pi^{(n)}\zeta+i\sum_{n\geq
1}\sqrt{n}\cdot\Pi^{(n)}\phi,$

we find that the linearization of (EQ) around the static solution (2.1) is a
linear dispersive equation

(2.2) $\frac{\partial u}{\partial t}+i\Lambda u=0,$

where the $3/2$-order elliptic operator $\Lambda$ is given by a multiplier

$\Lambda=\sum_{n\geq 2}\sqrt{n(n-1)(n+2)}\cdot\Pi^{(n)}=:\sum_{n\geq
0}\Lambda(n)\Pi^{(n)}.$

Note that $(\zeta,\phi)$ is completely determined by $u$. At the linear level,
there must hold $\Pi^{(0)}u\equiv 0$ because the first variation of volume
must be zero; and $\Pi^{(1)}u\equiv 0$ because of the conservation laws (1.7).

Let us also re-write the original nonlinear system (EQ) into a form that
better illustrates its perturbative nature. For simplicity, we use
$O(u^{\otimes k})$ to abbreviate a quantity that can be controlled by
$k$-linear expressions in $u$, and disregard its continuity properties for the
moment. For example, $\|u\|_{H^{1}}^{2}+\|u\|_{H^{2}}^{4}$ is an expression of
order $O(u^{\otimes 2})$ when $u\to 0$.

Since the operator $\Lambda$ acts degenerately on
$\mathcal{E}^{(0)}\oplus\mathcal{E}^{(1)}$, we should be more careful about
the eigenmodes $\Pi^{(0)}u$ and $\Pi^{(1)}u$. The volume preservation equation
(1.6) implies $\partial_{t}\Pi^{(0)}\zeta=O(u^{\otimes 2})$. Projecting (EQ)
to $\mathcal{E}^{(1)}$, which is spanned by the components of $N_{0}$, we
obtain $\partial_{t}\Pi^{(1)}\zeta=\Pi^{(1)}\phi+O(u^{\otimes 2})=O(u^{\otimes
2})$ since the conservation law (1.7) implies $\Pi^{(1)}\phi=O(u^{\otimes
2})$; and $\partial_{t}\Pi^{(1)}\phi=O(u^{\otimes 2})$ since
$H^{\prime}(\iota_{0})=-\Delta-2$ annihilates $\mathcal{E}^{(1)}$. We can thus
formally re-write the nonlinear system (EQ) as the following:

(2.3) $\frac{\partial u}{\partial t}+i\Lambda u=\mathfrak{N}(u),$

with $\mathfrak{N}(u)=O(u^{\otimes 2})$ vanishing quadratically as $u\to 0$.
Note that we are disregarding all regularity problems at the moment.

###  2.2. Question at Linear Level

At the linear level, our first unanswered question is

######  Question 1.

Does the solution of the linear capillary spherical water waves equation (2.2)
satisfy a Strichartz type estimate of the form

$\|e^{it\Lambda}f\|_{L^{p}_{T}L^{q}_{x}}\lesssim_{T}\|f\|_{H^{s}},$

where $L^{p}_{T}L^{q}_{x}=L^{p}([0,T];L^{q}(S^{2}))$, and the admissible
indices $(p,q)$ and $s$ should be determined?

Answer to Question 1 should be important in understanding the _dispersive_
nature of linear capillary spherical water waves. For Schrödinger equation on
a compact manifold, a widly cited result was obtained by Burq-Gérard-Tzvetkov
[BGT04]:

######  Theorem 2.1.

On a general compact Riemannian manifold $(M^{d},g)$, there holds

$\|e^{it\Delta_{g}}f\|_{L^{p}_{t}L^{q}_{x}([0,T]\times
M)}\lesssim_{T}\|f\|_{H^{1/p}(M)},$

where

$\frac{2}{p}+\frac{d}{q}=\frac{d}{2},\quad p>2.$

The authors used a time-localization argument for the parametrix of
$\partial_{t}-i\Delta_{g}$ to prove this result. For the sphere $S^{d}$, this
inequality is not optimal. The authors further used a Bourgain space argument
to obtain the optimal Strichartz inequality:

######  Theorem 2.2.

Let $(S^{d},g)$ be the standard n-dimensional sphere. For a function $f\in
C^{\infty}(S^{d})$, there holds the Strichartz inequality

$\|e^{it\Delta_{g}}f\|_{L^{p}_{t}L^{q}_{x}([0,T]\times
M)}\lesssim_{T}\|f\|_{H^{s}(M)},\quad s>s_{0}(d)$

where

$s_{0}(2)=\frac{1}{8},\quad s_{0}(d)=\frac{d}{4}-\frac{1}{2},\,d\geq 3.$

Furthermore these inequalities are optimal in the sense that the Sobolev index
$s$ cannot be less than or equal to $s_{0}(d)$.

The proof is the consequence of two propositions. The first one is the
“decoupling inequality on compact manifolds”, in particular the following
result proved by Sogge [Sog88]:

######  Proposition 2.1.

Let $\Pi_{k}$ be the spectral projection to eigenspaces with eigenvalues in
$[k^{2},(k+1)^{2}]$ on $S^{d}$. Then there holds

$\|\Pi_{k}\|_{L^{2}\to L^{q}}\leq C_{q}n^{s(q)},$

where

$s(q)=\left\\{\begin{matrix}\frac{d-1}{2}\left(\frac{1}{2}-\frac{1}{2q}\right),\quad&2\leq
q\leq\frac{2(d+1)}{d-1}\\\
\frac{d-1}{2}-\frac{d}{q},\quad&\frac{2(d+1)}{d-1}\leq
q\leq\infty.\end{matrix}\right.$

These estimates are sharp in the following sense: if $h_{k}$ is a zonal
spherical harmonic function of degree $k$ on $S^{d}$, then as $k\to\infty$,

$\|h_{k}\|_{L^{q}}\simeq C_{q}k^{s(q)}\|h_{k}\|_{L^{2}}.$

The second one is a Bourgain space embedding result:

######  Proposition 2.2.

For a function $f\in C^{\infty}_{0}(\mathbb{R}\times S^{d})$, define the
Bourgain space norm

$\|f(t,x)\|_{X^{s,b}}:=\left\|\langle\partial_{t}+i\Delta_{g}\rangle^{b}f(t,x)\right\|_{L^{2}_{t}H^{s}_{x}}.$

Then for $b>1/2$ and $s>s_{0}(d)$, there holds

$\|f\|_{L^{4}(\mathbb{R}\times S^{d})}\leq C_{s,b}\|f\|_{X^{s,b}}.$

The key ingredient for proving this proposition is the following number-
theoretic result:

$\\#\\{(p,q)\in\mathbb{N}^{2}:p^{2}+q^{2}=A\\}=O(A^{\varepsilon}).$

As for optimality of the Strichartz inequality, the authors of [BGT04]
implemented standard results of Gauss sums.

The parametrix and Bourgain space argument can be repeated without essential
change for the linear capillary spherical water waves equation (2.2), but this
time the Bourgain space argument would be more complicated: the Bourgain space
norm is now

$\|f(t,x)\|_{X^{s,b}}:=\left\|\langle\partial_{t}+i\Lambda\rangle^{b}f(t,x)\right\|_{L^{2}_{t}H^{s}_{x}},$

and the embedding result becomes

$\|f\|_{L^{4}(\mathbb{R}\times S^{2})}\leq C_{s,b}\|f\|_{X^{s,b}}$

for $f\in C^{\infty}_{0}(\mathbb{R}\times S^{2})$ and all $s>3\rho/8+1/8$,
$b>1/2$, where $\rho$ is the infimum of all exponents $\rho^{\prime}$ such
that when $A\to\infty$, the number

$\\#\left\\{(n_{1},n_{2})\in\mathbb{N}^{2}:\,\frac{1}{2}\leq\frac{n_{2}}{n_{1}}\leq
2,\,|\Lambda(n_{1})+\Lambda(n_{2})-A|\leq\frac{1}{2}\right\\}\leq
C_{\rho^{\prime}}A^{\rho^{\prime}}.$

Some basic analytic number theory implies $\rho=1/3$, and thus the range of
$s$ is $s>1/4$. Surprisingly this index is not better than that predicted by
the parametrix method. It remains unknown whether this index could be further
optimized.

To close this subsection, we note that the capillary spherical water wave
lives on a compact region, so the dispersion does not take away energy from a
locality to infinity. This is a crucial difference between waves on compact
regions and waves in Euclidean spaces. In particular, we do not expect decay
estimate for $e^{it\Lambda}f$. For the nonlinear problem (EQ), techniques like
vector field method (Klainerman-Sobolev type inequalities) do not apply.

###  2.3. Rotationally Symmetric Solutions: Bifurcation Analysis

Illuminated by observations in hydrodynamical experiments under zero gravity,
and suggested by the existence of standing gravity capillary water waves due
to Alazard-Baldi [AB15], we propose the following conjecture:

######  Conjecture 2.1.

There is a Cantor family of small amplitude periodic solutions to the
spherical capillary water waves system (EQ).

Let us conduct the bifurcation analysis that suggests why this conjecture
should be true. By time rescaling, we aim to find solution
$(\zeta,\phi,\omega_{0})$ of the following system that is $2\pi$-peiodic in
$t$:

(2.4) $\left\\{\begin{aligned} \omega_{0}\frac{\partial\zeta}{\partial
t}&=\frac{1}{N_{0}\cdot N(\iota)}D[\zeta]\phi,\\\
\omega_{0}\frac{\partial\phi}{\partial t}&=\frac{1}{N_{0}\cdot
N(\iota)}\left(B_{\zeta}\phi\cdot N_{0}\right)\cdot
D[\zeta]\phi-\frac{1}{2}|B_{\zeta}\phi|^{2}-H(\iota),\end{aligned}\right.$

together with the conservation laws (1.6)-(1.8). Here we refer $\omega_{0}>0$
as the _fundamental frequency_. The linearization of this system at the
equilibrium $(\zeta,\phi)=(0,0)$ is

(2.5) $L_{\omega_{0}}\left(\begin{matrix}\zeta\\\
\phi\end{matrix}\right):=\left(\begin{matrix}\omega_{0}\partial_{t}&-D[0]\\\
-\Delta-2&\omega_{0}\partial_{t}\end{matrix}\right)\left(\begin{matrix}\zeta\\\
\phi\end{matrix}\right)=0,\quad\Pi^{(0)}\zeta=\Pi^{(1)}\zeta=\Pi^{(1)}\phi=0.$

We restrict to _rotationally symmetric_ solutions of the system: that is,
water droplets which are always rotationally symmetric with a fixed axis. In
addition, we require $\zeta$ to be even in $t$ and $\phi$ to be odd in $t$.
The solution thus should take the form

$\zeta(t,x)=\sum_{j,n\geq
0}\zeta_{jn}\cos(jt)Y_{n}(x),\quad\phi(t,x)=\sum_{j\geq 1,n\geq
0}\phi_{jn}\sin(jt)Y_{n}(x),$

where $Y_{n}$ is the $n$’th zonal spherical harmonic, i.e. the (unique)
normalized spherical harmonic of degree $n$ that is axially symmetric. In
spherical coordinates this means that
$Y_{n}(\theta,\varphi)=P_{n}(\cos\theta)$, where $P_{n}$ is the $n$’th
Legendre polynomial. Since $\phi_{0n}$ are irrelevant we fix them to be 0.
Then

$L_{\omega_{0}}\left(\begin{matrix}\zeta\\\
\phi\end{matrix}\right)=\sum_{j,n\geq
0}\left(\begin{matrix}(-\omega_{0}j\zeta_{jn}-n\phi_{jn})\sin(jt)Y_{n}(x)\\\
\left((n-1)(n+2)\zeta_{jn}+\omega_{0}j\phi_{jn}\right)\cos(jt)Y_{n}(x)\end{matrix}\right).$

In order that $(\zeta,\phi)^{\mathrm{T}}\in\mathrm{Ker}L_{\omega_{0}}$, at the
level $n=0$, we must have $\zeta_{j0}=\phi_{j0}=0$ for all $j\geq 0$. At the
level $n=1$, we have $\zeta_{01}=\phi_{01}=0$, and for $j\geq 1$ there holds
$\omega_{0}j\zeta_{j1}-\phi_{j1}=0$ and $\omega_{0}j\phi_{j1}=0$, so
$\zeta_{j1}=\phi_{j1}=0$ for all $j\geq 0$. Hence $\zeta_{jn},\phi_{jn}$ can
be nonzero only for $j\geq 1$ and $n\geq 2$.

Consequently, $L_{\omega_{0}}$ has a one-dimensional kernel if and only if the
Diophantine equation

(2.6) $\omega_{0}^{2}j^{2}=n(n-1)(n+2),\quad j\geq 1,\,n\geq 2$

has exactly one solution $(j_{0},n_{0})$. If $\omega_{0}$ has this property,
then at the linear level, the lowest frequency of oscillation is

$\omega_{0}j_{0}=\sqrt{n_{0}(n_{0}-1)(n_{0}+2)}.$

We look into this Diophantine equation. Obviously $\omega_{0}^{2}$ has to be a
rational number. The equation is closely related to a family of elliptic
curves over $\mathbb{Q}$:

$E_{c}:y^{2}=x(x-c)(x+2c)=x^{3}+c^{2}x^{2}-2c^{2}x,\quad c\in\mathbb{N}.$

If we set $a/b=\omega_{0}^{2}$ (irreducible fraction), then integral solutions
of (2.6) are in 1-1 correspondence with integral points with natural number
coordinates on the elliptic curve $E_{ab}$, under the following map:

$(j_{0},n_{0})\to(abn_{0},a^{2}bj_{0})\in E_{ab}.$

Thus we just need to find natural numbers $a,b$ such that there is a unique
(up to negation) integral point $(x,y)\in E_{ab}$, where $x>0$ is divided by
$ab$, and $y$ is divided by $a^{2}b$. We seek for $n_{0}$ as small as possible
with such property, which gives lowest frequency of oscillation as small as
possible. For $ab=1,\cdots,50$, we find that if $ab=15$, then the integral
points on elliptic curve $E_{15}$ (up to negation of Mordel-Weil group) are

$(-30,0),\,(-5,50),\,(0,0),\,(15,0),\,(24,108),\,(90,900).$

The only point $(x,y)$ with $ab|x$ and $ab^{2}|y$ is (90,900), which gives
$n_{0}=6$, and the lowest frequency
$\Lambda(n_{0})=\sqrt{n_{0}(n_{0}-1)(n_{0}+2)}$ of oscillation is
$4\sqrt{15}\simeq 15.49\cdots$.

There are other choices of $a,b$. We list down the value of $ab$ below 50, the
corresponding $n_{0}$ and the lowest frequency $\Lambda(n_{0})$:

$ab$ | $n_{0}$ | $\Lambda(n_{0})$
---|---|---
15 | 6 | $4\sqrt{15}$
17 | 49 | $4\sqrt{323}$
22 | 9 | $6\sqrt{22}$
26 | 50 | $15\sqrt{78}$
42 | 7 | $3\sqrt{42}$
46 | 576 | $2040\sqrt{46}$
50 | 25 | $90\sqrt{2}$

See Appendix A for the MAGMA code used to find these values. This list
suggests that $n_{0}=6$ might be the smallest order that meets the
requirement, but this remains unproved. We summarize these into the following
number-theoretic question:

######  Question 2.

For the family of elliptic curves

$E_{ab}:\,y^{2}=x(x-ab)(x+2ab),\quad a,b\in\mathbb{N},$

how many choices of $a,b\in\mathbb{N}$ are there such that, there is exactly
one integral point $(x,y)\in E_{ab}$ with $x,y>0$ and $ab|x$, $a^{2}b|y$? For
such $a,b$ and integral point $(x,y)$, is the minimal value of $x/(ab)$
exactly 6, or is it smaller?

Of course, a complete answer of Question 2 should imply very clear
understanding of periodic solutions of the spherical capillary water waves
equation constructed using bifurcation analysis. But at this moment we are
satisfied with existence, so we may pick any $\omega_{0}=a/b$ and
$(j_{0},n_{0})$ that meets the requirement, for example the simplest case
$n_{0}=6$, and any of the following choices of $\omega_{0}$ and $j_{0}$:

$\omega_{0}=\sqrt{15},\,j_{0}=4;\quad\omega_{0}=\sqrt{\frac{1}{15}},\,j_{0}=60;\quad\omega_{0}=\sqrt{\frac{3}{5}},\,j_{0}=20;\quad\omega_{0}=\sqrt{\frac{5}{3}},\,j_{0}=12.$

We thus refine our conjecture as follows:

######  Conjecture 2.2.

Let $(n_{0},j_{0})$ be a pair of natural numbers with $n_{0}\geq 2$,
$j_{0}\geq 1$, and set $\omega_{0}=\sqrt{\Lambda(n_{0})}/j_{0}$. Suppose that
the only natural number solution of the Diophantine equation

$\omega_{0}^{2}j_{0}=\Lambda(n_{0})^{2}=n_{0}(n_{0}-1)(n_{0}+1)$

is $(j_{0},n_{0})$. Then there is a Cantor set with positive measure of
parameters $\omega$, clustered near $\omega_{0}$, such that the spherical
capillary water waves equation (2.4) admits small amplitude periodic solution
with frequency $\omega$.

The counterpart of Conjecture 2.2 for gravity capillary standing water waves
was proved by Alazard-Baldi [AB15] using a Nash-Moser type theorem. The key
technique in their proof was to find a conjugation of the linearized operator
of the gravity capillary water waves system on $\mathbb{T}^{1}$ to an operator
of the form

$\omega\partial_{t}+iT+i\lambda_{1}|D_{x}|^{1/2}+i\lambda_{-1}|D_{x}|^{-1/2}+\text{Operator
of order}\leq-\frac{3}{2},$

where $T$ is an elliptic Fourier multiplier of order $3/2$, and
$\lambda_{1},\lambda_{-1}$ are real constants. The frequency $\omega$ lives in
a Cantor type set that clusters around a given frequency so that the kernel of
the linearized operator is 1-dimensional. With this conjugation, they were
able to find periodic solutions of linearized problems required by Nash-Moser
iteration.

It is expected that this technique could be transplanted to the equation
(2.4), since our analysis for (2.6) suggests that the 1-dimensional kernel
requirement for bifurcation analysis is met. It seems that the greatest
technical issue is to find a suitable conjugation that takes the linearized
operator of (2.4) to an operator of the form

$\omega\partial_{t}+i(T_{3/2}+T_{1/2}+T_{-1/2})+\text{Operator of
order}\leq-\frac{3}{2},$

where each $T_{k}$ is a real Fourier multiplier acting on spherical harmonics.
The difficulty is that, since we are working with pseudo-differential
operators on $S^{2}$, the formulas of symbolic calculus are not as neat as
those on flat spaces. It seems necessary to implement some global harmonic
analysis for compact homogeneous spaces, e.g. extension of results collected
in Ruzhansky-Turunen [RT09]. Unfortunately, those results do not include
pseudo-differential operators with “rough coefficients” and para-differential
operators, so it seems necessary to re-write the whole theory.

###  2.4. Number-Theoretic Obstruction with Normal Form Reduction

As pointed out in Section 1, the system (2.3) is locally well-posed due to a
result in [BG98]. General well-posedness results [CS07], [SZ08] for free
boundary value problem of Euler equation also apply. As for lifespan estimate
for initial data $\varepsilon$-close to the static solution (2.1), it should
not be hard to conclude that the lifespan should be bounded below by
$1/\varepsilon$. The result relates to the fact that the sphere is a _stable_
critical point of the area functional, cf. [BC12]. This is nothing new: a
suitable energy inequality should imply it. However, although the clue is
clear, the implementation is far from standard since we are working on a
compact manifold. A rigorous proof still calls for hard technicalities.

Now we will be looking into the nonlinear equation (2.3) for its longer time
behavior. Although appearing similar to the well-studied water waves equation
in e.g. [IP19], [Lan05], [Wu97], [Wu99], there is a crucial difference between
the dispersive relation in (2.3) and the well-studied water waves equations:
the dispersive relation exhibits a strong _rigidity property_ , i.e. the
arbitrary physical constants enter into the dispersive relation $\Lambda$ only
as _scaling factors_. For the gravity-capillary water waves, the linear
dispersive relation reads

$\sqrt{g|\nabla|+\sigma|\nabla|^{3}},$

where $g$ is the gravitational constant and $\sigma$ is the surface tension
coefficient. For the Klein-Gordon equation on a Riemannian manifold, the
linear dispersive relation reads

$\sqrt{-\Delta+m^{2}},$

where $m$ is the mass. In [IP19], Ionescu and Pusateri referred such
dispersive relations as having _non-degenerate dependence on physical
parameter_ , while in our context it is appropriate to refer to the dependence
as _degenerate_. We will see that this crucial difference brings about severe
obstructions for the long-time well-posedness of the system.

Following the idea of Delort and Szeftel [DS04], we look for a normal form
reduction of (2.3) and explain why the rigidity property could cause
obstructions. Not surprisingly, the obstruction is due to resonances, and
strongly relates to the solvablity of a Diophantine equation. Delort and
Szeftel cast a normal form reduction to the small-initial-data problem of
quasilinear Klein-Gordon equation on the sphere and obtained an estimate on
the lifespan longer than the one provided by standard energy method. After
their work, normal form reduction has been used by mathematicians to
understand water waves on flat tori, for example [BD+18], [BFP22] and [IP19].
The idea was inspired by the normal form reduction method introduced by Shatah
[Sha85]: for a quadratic perturbation of a linear dispersive equation

$\partial_{t}u+iLu=N(u)=O(u^{\otimes 2}),$

using a new variable $u+B(u,\bar{u})$ with a suitably chosen quadratic
addendum $B(u,\bar{u})$ can possibly eliminate the quadratic part of $N$, thus
extending the lifespan estimate beyond the standard $1/\varepsilon$.

So we shall write the quadraticr part of $\mathfrak{N}(u)$ as

$\sum_{n_{3}\geq 0}\Pi^{(n_{3})}\left[\sum_{n_{1}\geq 0}\sum_{n_{2}\geq
0}\mathcal{M}_{1}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}u\right)+\mathcal{M}_{2}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}\bar{u}\right)+\mathcal{M}_{3}\left(\Pi^{(n_{1})}\bar{u},\Pi^{(n_{2})}\bar{u}\right)\right],$

where $\mathcal{M}_{1},\mathcal{M}_{2},\mathcal{M}_{3}$ are complex bi-linear
operators, following the argument of Section 4 in [DS04]. They are independent
of $t$ since the right-hand-side of the equation does not depend on $t$
explicitly. Let’s look for a diffeomorphism

$u\to v:=u+\mathbf{B}[u,u]$

in the function space $C^{\infty}(S^{2})$, where $\mathbf{B}$ is a bilinear
operator, so that the equation (2.3) with quadratic nonlinearity reduces to an
equation with cubic nonlinearity. The $\mathbf{B}[u,u]$ is supposed to take
the form
$\mathbf{B}[u,u]=\mathbf{B}_{1}[u,u]+\mathbf{B}_{2}[u,u]+\mathbf{B}_{3}[u,u]$,
with

$\mathbf{B}_{1}[u,u]=\sum_{n_{3}\geq 0}\sum_{n_{1},n_{2}\geq
0}b_{1}(n_{1},n_{2},n_{3})\Pi^{(n_{3})}\mathcal{M}_{1}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}u\right),$
$\mathbf{B}_{2}[u,u]=\sum_{n_{3}\geq 0}\sum_{n_{1},n_{2}\geq
0}b_{2}(n_{1},n_{2},n_{3})\Pi^{(n_{3})}\mathcal{M}_{2}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}\bar{u}\right),$
$\mathbf{B}_{3}[u,u]=\sum_{n_{3}\geq 0}\sum_{n_{1},n_{2}\geq
0}b_{3}(n_{1},n_{2},n_{3})\Pi^{(n_{3})}\mathcal{M}_{3}\left(\Pi^{(n_{1})}\bar{u},\Pi^{(n_{2})}\bar{u}\right),$

where the $b_{j}(n_{1},n_{2},n_{3})$’s are complex numbers to be determined.
Implementing (2.3), we find

$\displaystyle(\partial_{t}$ $\displaystyle+i\Lambda)(u+\mathbf{B}[u,u])$
$\displaystyle=\mathfrak{N}(u)+\sum_{n_{3}\geq 0}\sum_{n_{1},n_{2}\geq
0}b_{1}(n_{1},n_{2},n_{3})\Pi^{(n_{3})}\mathcal{M}_{1}\left(\Pi^{(n_{1})}\partial_{t}u,\Pi^{(n_{2})}u\right)$
$\displaystyle\quad+\sum_{n_{3}\geq 0}\sum_{n_{1},n_{2}\geq
0}b_{1}(n_{1},n_{2},n_{3})\Pi^{(n_{3})}\mathcal{M}_{1}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}\partial_{t}u\right)+(\text{similar
terms})$ $\displaystyle\quad+\sum_{n_{3}\geq 0}\sum_{n_{1},n_{2}\geq
0}i\Lambda(n_{3})b_{1}(n_{1},n_{2},n_{3})\Pi^{(n_{3})}\mathcal{M}_{1}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}u\right)$
$\displaystyle=\mathfrak{N}(u)+\sum_{n_{3}\geq 0}\sum_{\min(n_{1},n_{2})\leq
1}i\left[\Lambda(n_{3})-\Lambda(n_{1})-\Lambda(n_{2})\right]b_{1}(n_{1},n_{2},n_{3})\Pi^{(n_{3})}\mathcal{M}_{1}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}u\right)$
$\displaystyle\quad+\sum_{n_{3}\geq 0}\sum_{n_{1},n_{2}\geq
2}i\left[\Lambda(n_{3})-\Lambda(n_{1})-\Lambda(n_{2})\right]b_{1}(n_{1},n_{2},n_{3})\Pi^{(n_{3})}\mathcal{M}_{1}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}u\right)$
$\displaystyle\quad+(\text{similar terms})+O(u^{\otimes 3}).$

We aim to eliminate most of the second order portions of $\mathfrak{N}(u)$.
The coefficients $b_{j}(n_{1},n_{2},n_{3})$ are fixed as follows:

(2.7) $\displaystyle b_{1}(n_{1},n_{2},n_{3})$
$\displaystyle=i\left[\Lambda(n_{3})-\Lambda(n_{1})-\Lambda(n_{2})\right]^{-1},\quad
n_{1},n_{2},n_{3}\geq 2$ $\displaystyle b_{2}(n_{1},n_{2},n_{3})$
$\displaystyle=i\left[\Lambda(n_{3})-\Lambda(n_{1})+\Lambda(n_{2})\right]^{-1},\quad
n_{1},n_{2},n_{3}\geq 2$ $\displaystyle b_{3}(n_{1},n_{2},n_{3})$
$\displaystyle=i\left[\Lambda(n_{3})+\Lambda(n_{1})+\Lambda(n_{2})\right]^{-1},\quad
n_{1},n_{2},n_{3}\geq 2$ $\displaystyle b_{1,2,3}(n_{1},n_{2},n_{3})$
$\displaystyle=0,\quad\text{if
}\Lambda(n_{3})\pm\Lambda(n_{1})\pm\Lambda(n_{2})=0\text{ or
}\min(n_{1},n_{2},n_{3})\leq 1,$

then a large portion of the second order part of $\Pi_{\geq
2}\mathfrak{N}\left(u\right)$ will be eliminated. In fact, for
$n_{1},n_{2},n_{3}\geq 2$, if
$\Lambda(n_{3})\pm\Lambda(n_{1})\pm\Lambda(n_{2})\neq 0$, then the term

$\Pi^{(n_{3})}\left[\sum_{n_{1},n_{2}\geq
2}\mathcal{M}_{1}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}u\right)+\mathcal{M}_{2}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}\bar{u}\right)+\mathcal{M}_{3}\left(\Pi^{(n_{1})}\bar{u},\Pi^{(n_{2})}\bar{u}\right)\right]$

is cancelled out. On the other hand, by the volume preservation equality (1.6)
and conservation law (1.7), there holds $\Pi^{(0)}u=O(u^{\otimes 2})$,
$\Pi^{(1)}u=O(u^{\otimes 2})$, so the low-low interaction

$\Pi_{\geq 2}\left[\sum_{\min(n_{1},n_{2})\leq
1}\mathcal{M}_{1}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}u\right)+\mathcal{M}_{2}\left(\Pi^{(n_{1})}u,\Pi^{(n_{2})}\bar{u}\right)+\mathcal{M}_{3}\left(\Pi^{(n_{1})}\bar{u},\Pi^{(n_{2})}\bar{u}\right)\right]$

is automatically $O(u^{\otimes 3})$.

Thus the existence and continuity of the normal form $\mathbf{B}[u,u]$ depends
on the property of the 3-way resonance equation

(2.8) $\Lambda(n_{3})-\Lambda(n_{1})-\Lambda(n_{2})=0,\quad
n_{1},n_{2},n_{3}\geq 2.\quad n_{3}\leq n_{1}+n_{2}$

which is equivalent to the Diophantine equation

(2.9) $[F(n_{1})+F(n_{2})-F(n_{3})]^{2}-4F(n_{1})F(n_{2})=0,\quad
n_{1},n_{2},n_{3}\geq 2,$

where $F(X)=X(X-1)(X+2)$.

If the tuple $(n_{1},n_{2},n_{3})$ is non-resonant, i.e. it is such that
$b_{1,2,3}(n_{1},n_{2},n_{3})\neq 0$, then some elementary number theoretic
argument will give a lower bound on $|b_{1,2,3}(n_{1},n_{2},n_{3})|$ in terms
of a negative power (can be fixed as $-9/2$) of $n_{1},n_{2},n_{3}$. This is
usually referred as _small divisor estimate_.

To study the distribution of resonant frequencies, we propose the following
unsolved question:

######  Question 3.

Does the Diophantine equation (2.9) have finitely many solutions?

However, the Diophantine equation (2.9) does admit non-trivial solutions
$(5,5,8)$ and $(10,10,16)$. In other words, the second order terms e.g.

$\Pi^{(8)}\mathcal{M}_{1}(\Pi^{(5)}u\cdot\Pi^{(5)}u),\quad\Pi^{(5)}\mathcal{M}_{1}(\Pi^{(8)}u\cdot\Pi^{(5)}\bar{u}),$

in the quadratic part of $\mathfrak{N}(u)$ cannot be eliminated by normal form
reduction. On the other hand, it seems to be very hard to determine whether
(2.9) still admits any other solution. We have the following proposition (the
author would like to thank Professor Bjorn Poonen for the proof):

######  Proposition 2.3.

The Diophantine equation (2.9) has no solution with $n_{1}\leq 10^{4}$ other
than $(5,5,8)$ and $(10,10,16)$.

The proof of this proposition is computer-aided. The key point is to use the
so-called Runge’s method to show that if $(n_{1},n_{2},n_{3})$ is a solution,
then there must hold $n_{2}=O(n_{1}^{2})$. For a given $n_{1}$, this reduces
the proof to numerical verification for _finitely many possibilities_. The
algorithm can of course be further optimized, but due to some algebraic
geometric considerations, it is reasonable to conjecture that the solution of
(2.9) should be very rare. In fact, there are two ways of viewing the problem.
We observe that if $(n_{1},n_{2},n_{3})$ is a solution, then
$F(n_{1})F(n_{2})$ must be a square, and the square free part of
$F(n_{1}),F(n_{2}),F(n_{3})$ must be the same. Further reduction turns the
problem into finding integral points on a family of elliptic curves

$Y^{2}=cF(X),\quad c\text{ is square-free},$

which is of course difficult, but since Siegel’s theorem asserts that there
are only finitely many integer points on an elliptic curve over $\mathbb{Q}$,
it is reasonable to conjecture that there are not “too many” solutions to
(2.9). We may also view the problem as finding integral (rational) points on a
given algebraic surface. The complex projective surface corresponding to (2.9)
is given by

$\displaystyle\mathfrak{V}:[X(X-W)(X+2W)$
$\displaystyle+Y(Y-W)(Y+2W)-Z(Z-W)(Z+2W)]^{2}$
$\displaystyle\quad=4X(X-W)(X+2W)Y(Y-W)(Y+2W),$

where $[X,Y,Z,W]$ is the homogeneous coordinate on $\mathbb{CP}^{3}$. With the
aid of computer, we obtain

######  Proposition 2.4.

The complex projective surface $\mathfrak{V}\subset\mathbb{CP}^{3}$ has
Kodaira dimension 2 (i.e. it is of general type under Kodaira-Enriquez
classification), and its first Betti number is 0.

See Appendix A for the code.

The first part of the proposition suggests that the rational points of
$\mathfrak{V}$ should be localized on finitely many algebraic curves laying on
$\mathfrak{V}$; nevertheless, this seemingly simple suggestion is indeed a
special case of the _Bomberi-Lang conjecture_ , a hard problem in number
theory (its planar case is known as the celebrated Faltings’s theorem). The
second part suggests that the rational points of $\mathfrak{V}$ should be rare
since its Albanese variety is a single point. But these are just heuristics
that solutions to (2.9) should be rare. In general, determining the
solvability of a given Diophantine equation is very difficult 222For example,
the seemingly simple Diophantine equation $x^{3}+y^{3}+z^{3}=42$ is in fact a
puzzle of more than 60 years, and its first solution was found recently by
Booker-Sutherland [BS21]. It is of extremely large magnitude: $42=(-80\ 538\
738\ 812\ 075\ 974)^{3}+80\ 435\ 758\ 145\ 817\ 515^{3}+12\ 602\ 123\ 297\
335\ 631^{3}$. Another example is the equation of same type
$x^{3}+y^{3}+z^{3}=3$. Beyond the easily found solutions (1, 1, 1), (4, 4,
-5), (4, -5, 4), (-5, 4, 4), the next solution reads $(569\ 936\ 821\ 221\
962\ 380\ 720,-569\ 936\ 821\ 113\ 563\ 493\ 509,-472\ 715\ 493\ 453\ 327\
032)$., as number theorists and arithmetic geometers generally believe.

The reason that such issues do not occur for water waves in the flat setting
or nonlinear Klein-Gordon equations is twofold. First of all, the resonance
equation is easily understood even in the degenerate case in the flat setting.
For example, the capillary water waves without gravity on $\mathbb{T}^{2}$ has
dispersive relation $|\nabla|^{1/2}$, and the 3-way resonance equation is

$\sqrt[4]{k_{1}^{2}+k_{2}^{2}}+\sqrt[4]{l_{1}^{2}+l_{2}^{2}}=\sqrt[4]{m_{1}^{2}+m_{2}^{2}},$

with the additional requirement $m=k+l$. We already know that the resonance
equation has no non-trivial solution at all, cf. [BFP22]. But even without
using $m=k+l$, we would be able to conclude that there are at most finitely
many non-trivial solutions from the celebrated Faltings’s theorem on rational
points on high-genus algebraic projective curves (although this is like using
a sledge hammer to crack a nut). Secondly, for the non-degenrate case, for
example the gravity-capillary waves, the dispersive relation reads
$\sqrt{g|\nabla|+\sigma|\nabla|^{3}}$, so if the ratio $\sigma/g$ is a
transcedental number then the 3-way resonance equation has no solution.
Furthermore, using some elementary calculus and a measure-theoretic argument,
it can be shown, not without technicalities, that the resonances admit certain
small-divisor estimates for almost all parameters. This is exactly the
argument employed by Delort-Szeftle [DS04], Berti-Delort [BD+18] and Ionescu-
Pusateri [IP19], so that their results were stated for _almost all_
parameters. These parameters are, roughly speaking, badly approximated by
algebraic numbers.

However, the resonance equation (2.8) is inhomogeneous and allows no arbitrary
physical parameter at all. Furthermore, since product of spherical harmonics
are no longer spherical harmonics in general, Fourier series techniques
employed by [BFP22] [BD+18] [IP19] that works for the torus are never valid
for $S^{2}$; for example, we cannot simply assume $n_{3}=n_{1}+n_{2}$ in
(2.8), as already illustrated by the solutions (5,5,8) (10,10,16). These are
the crucial differences between the capillary spherical water waves and all
known results for water waves in the flat setting.

###  2.5. Heuristics for Lifespan Estimate

To summarize, almost global lifespan estimate of (2.3) depends on the
difficult number theoretic question 3. Before it is fully resolved, we can
only expect partial results regarding the normal form transformation.

If there are only finitely many solutions to the Diophantine equation (2.9),
then under the normal form reduction $u\to v=u+\mathbf{B}[u,u]$ with
coefficients given by (2.7), the equation (2.3) is transformed into the
following system:

$\displaystyle\frac{\partial}{\partial t}\Pi_{c}v$ $\displaystyle=O(v^{\otimes
2}),$ $\displaystyle\frac{\partial}{\partial t}(1-\Pi_{c})v$
$\displaystyle=O(v^{\otimes 3}),$

where $\Pi_{c}$ is the orthogonal projection to
$\bigoplus_{n_{3}}\mathcal{E}^{(n_{3})}\subset L^{2}(S^{2})$, with $n_{3}$
being either 0 or 1, or exhausting the third component of all nontrivial
solutions of (2.9).

There is no reasonable assertion to be made if the conjecture fails. However,
if the conjecture does hold true, then we can expect that the lifespan
estimate for $\varepsilon$-Cauchy data of (EQ) goes beyond $\varepsilon^{-1}$,
as what we expect for gravity water waves in the periodic setting, e.g. in
[IP19]:

######  Conjecture 2.3.

If the Diophantine equation (2.9) has only finitely many solutions, then there
is some $\alpha>0$ such that for $\varepsilon$-Cauchy data of (EQ), the
lifespan goes beyond $\varepsilon^{-(1+\alpha)}$ as $\varepsilon\to 0$.

Let’s explain the heuristic as follows. The argument we aim to implement is
the standard “continuous induction method”, i.e. for some suitably large $s$
and $K$ and suitable $\alpha>0$, assuming $T=\varepsilon^{-(1+\alpha)}$ and
$\sup_{t\in[0,T]}\|v\|_{H^{s}(g_{0})}\leq K\varepsilon$, we try to prove a
better bound $\sup_{t\in[0,T]}\|v\|_{H^{s}(g_{0})}\leq K\varepsilon/2$. Here
$g_{0}$ is the standard metric on $S^{2}$. It is intuitive to expect such a
result for the cubic equation ${\partial_{t}}(1-\Pi_{c})v=O(v^{\otimes 3})$.
As for the quadratic equation ${\partial_{t}}\Pi_{c}v=O(v^{\otimes 2})$, it is
crucial to implement the conservation of energy $\mathbf{H}[\zeta,\phi]\equiv
4\pi$ for a solution. We summarize it as

######  Proposition 2.5.

Fix $T>0$. Let $u$ be a smooth solution of (2.3) and $v=u+\mathbf{B}[u,u]$ be
as above. Suppose for some suitably large $s$ and $K$, there holds

$\sup_{t\in[0,T]}\|v\|_{H^{s}(g_{0})}\leq K\varepsilon$

with $\varepsilon$ sufficiently small. Then there is in fact a better bound
for the low frequency part $\Pi_{c}v$:

$\sup_{t\in[0,T]}\|\Pi_{c}v\|_{H^{s}(g_{0})}\leq K\varepsilon/4.$

###### Proof.

We consider the “approximate” energy functional

$\mathbf{H}_{0}[\zeta,\phi]=4\pi+\int_{S^{2}}2\zeta\cdot
d\mu_{0}+\frac{1}{2}\int_{S^{2}}(|\nabla_{0}\zeta|^{2}+2|\zeta|^{2})d\mu_{0}+\frac{1}{2}\int_{S^{2}}\left||\nabla_{0}^{1/2}|\phi\right|^{2}d\mu_{0},$

where $\mu_{0}$ is the standard area measure on $S^{2}$ and $\nabla_{0}$ is
the standard connection on $S^{2}$. This is nothing but the quadratic
approximation to $\mathbf{H}[\zeta,\phi]$ in (1.8) at $(0,0)$, so there holds
$\mathbf{H}_{0}[\zeta,\phi]=\mathbf{H}[\zeta,\phi]+O(u^{\otimes 3})$. Using
volume preservation (1.6), we obtain $\int_{S^{2}}\zeta
d\mu_{0}=-\|\zeta\|_{L^{2}(g_{0})}^{2}+O(\zeta^{\otimes 3})$, so summarizing
we have

(2.10)
$\frac{1}{2}\int_{S^{2}}(|\nabla_{0}\zeta|^{2}-2|\zeta|^{2})d\mu_{0}+\frac{1}{2}\int_{S^{2}}\left||\nabla_{0}^{1/2}|\phi\right|^{2}d\mu_{0}=O(u^{\otimes
3}).$

Note that we used the conservation law $\mathbf{H}[\zeta,\phi]\equiv 4\pi$. By
spectral calculus on $S^{2}$, we have

$\|\zeta\|_{H^{1}(g_{0})}^{2}\simeq\|\Pi^{(0)}\zeta\|_{L^{2}(g_{0})}^{2}+\|\Pi^{(1)}\zeta\|_{L^{2}(g_{0})}^{2}+\int_{S^{2}}(|\nabla_{0}\zeta|^{2}-2|\zeta|^{2})d\mu_{0},$

and by volume preservation and conservation of momentum (1.7) we find

(2.11)
$\|\zeta\|_{H^{1}(g_{0})}^{2}\simeq\int_{S^{2}}(|\nabla_{0}\zeta|^{2}-2|\zeta|^{2})d\mu_{0}+O(\zeta^{\otimes
4}).$

Now, for some $N_{0}>0$ relating to the loss of regularity caused by
$\mathbf{B}$, we may choose $s>>2N_{0}$. Then if $\|v\|_{H^{s}(g_{0})}\leq
K\varepsilon$, it follows that $\|u\|_{H^{s-N_{0}}(g_{0})}\leq
K^{\prime}\varepsilon$, so by (2.10) and (2.11) we have

$\|u\|_{L^{2}(g_{0})}^{2}\simeq\|\zeta\|_{H^{1}(g_{0})}^{2}+\|\phi\|_{H^{1/2}(g_{0})}^{2}\leq
K^{\prime}\varepsilon^{3}.$

Thus

$\displaystyle\|v\|_{L^{2}(g_{0})}$ $\displaystyle\leq
C\|u\|_{L^{2}(g_{0})}+C\|u\|_{H^{N_{0}}(g_{0})}^{2}$ $\displaystyle\leq
C\varepsilon^{3/2}+K^{\prime}\varepsilon^{2}.$

Since the spectrum of $\Pi_{c}v$ is bounded, by Bernstein type inequality we
have

$\|\Pi_{c}v\|_{H^{s}(g_{0})}\leq C\|v\|_{L^{2}(g_{0})}\leq
K^{\prime}\varepsilon^{3/2}(1+\varepsilon^{1/2}).$

If $\varepsilon$ is sufficiently small then this implies
$\|\Pi_{c}v\|_{H^{s}(g_{0})}\leq K\varepsilon/4$. ∎

We point out that the above proof is independent from the magnitude of the
lifespan $T$, so it is always applicable as long as the cubic equation
${\partial_{t}}(1-\Pi_{c})v=O(v^{\otimes 3})$ is well-understood. There are
two crucial points in the proof of Proposition 2.5: the conservation of
energy, and that the projection $\Pi_{c}$ is of finite rank, so that a
Bernstein type inequality holds. The last fact holds only if there are
finitely many 3-way resonances, i.e. there are only finitely many solutions to
the Diophantine equation (2.8).

Finally, we propose an even more ambitious conjecture concerning global
dynamical properties of spherical water droplets, which is again illuminated
by observation in hydrodynamical experiments under zero gravity, and also
suggested by the results of Berti-Montalto [BM20]:

######  Conjecture 2.4.

If the Diophantine equation (2.9) has only finitely many solutions, then a KAM
type result holds for (EQ): there is a family of infinitely many quasi-
periodic solutions of (EQ), depending on a parameter which takes value in a
Cantor-type set.

##  Appendix A MAGMA Code

MAGMA is a large, well-supported software package designed for computations in
algebra, number theory, algebraic geometry, and algebraic combinatorics. In
this appendix, we give the MAGMA code used to conduct computations on
Diophantine equations related to the spherical capillary water waves system.

###  A.1. Integral Points on Elliptic Curve

We can find all integral points on a given elliptic curve over $\mathbb{Q}$
using MAGMA. For a monic cubic polynomial $f(x)$, the function
`EllipticCurve(f)` creates the elliptic curve

$E:y^{2}=f(x),$

and the function `IntegralPoints(E)` returns a sequence containing all the
integral points on $E$ under the homogeneous coordinate of $\mathbb{QP}^{2}$,
modulo negation. We use this to find out all integral points on the elliptic
curve

$E_{c}:\,y^{2}=x^{3}+cx^{2}-2c^{2}x=x(x-c)(x+2c).$

for natural number $c\leq 50$. The MAGMA code is listed below, which excludes
all the $c$’s such that there are only trivial integral points
$\\{(-c,0),(0,0),(c,0)\\}$ on $E_{c}$.

¿ Qx¡x¿ := PolynomialRing(Rationals()); ¿ for c in [1..50] do ¿ E :=
EllipticCurve(x^3+c*x^2-2*c^2*x); ¿ S, reps := IntegralPoints(E); ¿ if # S gt
3 then ¿ print c, E; ¿ print S; ¿ end if; ¿ end for; 2 Elliptic Curve defined
by y^2 = x^3 + 2*x^2 - 8*x over Rational Field [ (-4 : 0 : 1), (-2 : 4 : 1),
(-1 : -3 : 1), (0 : 0 : 1), (2 : 0 : 1), (4 : 8 : 1), (8 : -24 : 1), (50 : 360
: 1) ] 8 Elliptic Curve defined by y^2 = x^3 + 8*x^2 - 128*x over Rational
Field [ (-16 : 0 : 1), (-8 : 32 : 1), (-4 : -24 : 1), (0 : 0 : 1), (8 : 0 :
1), (9 : -15 : 1), (16 : 64 : 1), (32 : -192 : 1), (200 : 2880 : 1) ] 13
Elliptic Curve defined by y^2 = x^3 + 13*x^2 - 338*x over Rational Field [
(-26 : 0 : 1), (0 : 0 : 1), (13 : 0 : 1), (121 : 1386 : 1) ] 15 Elliptic Curve
defined by y^2 = x^3 + 15*x^2 - 450*x over Rational Field [ (-30 : 0 : 1), (-5
: -50 : 1), (0 : 0 : 1), (15 : 0 : 1), (24 : 108 : 1), (90 : -900 : 1) ] 17
Elliptic Curve defined by y^2 = x^3 + 17*x^2 - 578*x over Rational Field [
(-34 : 0 : 1), (-32 : -56 : 1), (0 : 0 : 1), (17 : 0 : 1), (833 : 24276 : 1) ]
18 Elliptic Curve defined by y^2 = x^3 + 18*x^2 - 648*x over Rational Field [
(-36 : 0 : 1), (-32 : -80 : 1), (-18 : 108 : 1), (-9 : -81 : 1), (0 : 0 : 1),
(18 : 0 : 1), (36 : 216 : 1), (72 : -648 : 1), (450 : 9720 : 1) ] 22 Elliptic
Curve defined by y^2 = x^3 + 22*x^2 - 968*x over Rational Field [ (-44 : 0 :
1), (-32 : -144 : 1), (0 : 0 : 1), (22 : 0 : 1), (198 : 2904 : 1) ] 23
Elliptic Curve defined by y^2 = x^3 + 23*x^2 - 1058*x over Rational Field [
(-46 : 0 : 1), (0 : 0 : 1), (23 : 0 : 1), (50 : -360 : 1) ] 26 Elliptic Curve
defined by y^2 = x^3 + 26*x^2 - 1352*x over Rational Field [ (-52 : 0 : 1),
(-49 : -105 : 1), (0 : 0 : 1), (26 : 0 : 1), (1300 : 47320 : 1) ] 30 Elliptic
Curve defined by y^2 = x^3 + 30*x^2 - 1800*x over Rational Field [ (-60 : 0 :
1), (-50 : 200 : 1), (-45 : -225 : 1), (-24 : -216 : 1), (-20 : 200 : 1), (-6
: 108 : 1), (0 : 0 : 1), (30 : 0 : 1), (36 : 144 : 1), (40 : -200 : 1), (75 :
-675 : 1), (90 : 900 : 1), (300 : 5400 : 1), (324 : -6048 : 1), (480 : -10800
: 1), (7290 : 623700 : 1), (10830 : -1128600 : 1), (226875 : 108070875 : 1) ]
32 Elliptic Curve defined by y^2 = x^3 + 32*x^2 - 2048*x over Rational Field [
(-64 : 0 : 1), (-32 : 256 : 1), (-16 : -192 : 1), (0 : 0 : 1), (32 : 0 : 1),
(36 : -120 : 1), (64 : 512 : 1), (128 : -1536 : 1), (800 : 23040 : 1) ] 33
Elliptic Curve defined by y^2 = x^3 + 33*x^2 - 2178*x over Rational Field [
(-66 : 0 : 1), (0 : 0 : 1), (33 : 0 : 1), (81 : -756 : 1) ] 35 Elliptic Curve
defined by y^2 = x^3 + 35*x^2 - 2450*x over Rational Field [ (-70 : 0 : 1),
(-49 : 294 : 1), (-45 : -300 : 1), (-40 : 300 : 1), (-14 : -196 : 1), (0 : 0 :
1), (35 : 0 : 1), (50 : 300 : 1), (175 : -2450 : 1), (224 : 3528 : 1), (280 :
-4900 : 1), (4410 : -294000 : 1), (14450 : -1739100 : 1) ] 39 Elliptic Curve
defined by y^2 = x^3 + 39*x^2 - 3042*x over Rational Field [ (-78 : 0 : 1), (0
: 0 : 1), (39 : 0 : 1), (147 : -1890 : 1) ] 42 Elliptic Curve defined by y^2 =
x^3 + 42*x^2 - 3528*x over Rational Field [ (-84 : 0 : 1), (-56 : -392 : 1),
(-12 : 216 : 1), (0 : 0 : 1), (42 : 0 : 1), (63 : -441 : 1), (294 : 5292 : 1)
] 43 Elliptic Curve defined by y^2 = x^3 + 43*x^2 - 3698*x over Rational Field
[ (-86 : 0 : 1), (-32 : -360 : 1), (0 : 0 : 1), (43 : 0 : 1) ] 46 Elliptic
Curve defined by y^2 = x^3 + 46*x^2 - 4232*x over Rational Field [ (-92 : 0 :
1), (0 : 0 : 1), (46 : 0 : 1), (26496 : 4316640 : 1) ] 50 Elliptic Curve
defined by y^2 = x^3 + 50*x^2 - 5000*x over Rational Field [ (-100 : 0 : 1),
(-50 : 500 : 1), (-25 : -375 : 1), (-4 : 144 : 1), (0 : 0 : 1), (50 : 0 : 1),
(100 : 1000 : 1), (200 : -3000 : 1), (1250 : 45000 : 1) ] Running Magma
V2.27-7. Seed: 821911319; Total time: 2.430 seconds; Total memory usage:
85.16MB.

###  A.2. Classification of Projective Surface

Some basic geometric parameters of the complex projective surface

$\displaystyle\mathfrak{V}:[X(X-W)(X+2W)$
$\displaystyle+Y(Y-W)(Y+2W)-Z(Z-W)(Z+2W)]^{2}$
$\displaystyle\quad=4X(X-W)(X+2W)Y(Y-W)(Y+2W)$

in $\mathbb{CP}^{3}$ can be computed using MAGMA. The author would like to
thank Professor Bjorn Poonen for introducing MAGMA and providing the code
listed below.

    > Q:=Rationals();
    > P<x,y,z,w>:=ProjectiveSpace(Q,3);
    > fx:=x*(x-w)*(x+2*w);
    > fy:=y*(y-w)*(y+2*w);
    > fz:=z*(z-w)*(z+2*w);
    > V:=Surface(P,(fz-fx-fy)^2-4*fx*fy);
    > KodairaEnriquesType(V);
    2 0 General type
    Running Magma V2.27-7.
    Seed: 1338492807; Total time: 0.670 seconds; Total memory usage: 32.09MB.

Note that the Kodaira dimension is invariant regardless of the choice of base
field, so it is legitimate to choose the base field to be $\mathbb{Q}$ in the
above code. The variable $w$ is used to homogenize the equation. The function
`KodairaEnriquesType(V)` returns three values for the given projective surface
$V$: the first is the Kodaira dimension, the second is irrelevant when the
Kodaira dimension is not $-\infty$, 1 or 0, and the third is the Kodaira-
Enriquez classification of the surface $X$.

## References

  * [AB15] Thomas Alazard and Pietro Baldi.  Gravity capillary standing water waves.  Archive for Rational Mechanics and Analysis, 217:741–830, 2015\.
  * [ABZ11] Thomas Alazard, Nicolas Burq, and Claude Zuily.  On the water-wave equations with surface tension.  2011\.
  * [AD13] Thomas Alazard and Jean-Marc Delort.  Global solutions and asymptotic behavior for two dimensional gravity water waves.  arXiv preprint arXiv:1305.4090, 2013.
  * [AM09] Thomas Alazard and Guy Métivier.  Paralinearization of the dirichlet to neumann operator, and regularity of three-dimensional water waves.  Communications in Partial Differential Equations, 34(12):1632–1704, 2009.
  * [BC12] Joao Lucas Barbosa and Manfredo do Carmo.  Stability of hypersurfaces with constant mean curvature.  In Manfredo P. do Carmo–Selected Papers, pages 221–235. Springer, 2012.
  * [BD+18] Massimiliano Berti, Jean-Marc Delort, et al.  Almost global solutions of capillary-gravity water waves equations on the circle.  Springer, 2018.
  * [BFP22] Massimiliano Berti, Roberto Feola, and Fabio Pusateri.  Birkhoff normal form and long time existence for periodic gravity water waves.  Communications on Pure and Applied Mathematics, 2022.
  * [BG98] Klaus Beyer and Matthias Günther.  On the cauchy problem for a capillary drop. part i: irrotational motion.  Mathematical methods in the applied sciences, 21(12):1149–1183, 1998.
  * [BGT04] Nicolas Burq, Pierre Gérard, and Nikolay Tzvetkov.  Strichartz inequalities and the nonlinear schrödinger equation on compact manifolds.  American Journal of Mathematics, 126(3):569–605, 2004.
  * [BM20] Massimiliano Berti and Riccardo Montalto.  Quasi-periodic standing wave solutions of gravity-capillary water waves, volume 263.  American mathematical society, 2020.
  * [BS21] Andrew R Booker and Andrew V Sutherland.  On a question of mordell.  Proceedings of the National Academy of Sciences, 118(11):e2022377118, 2021.
  * [CHS10] Hans Christianson, Vera Mikyoung Hur, and Gigliola Staffilani.  Strichartz estimates for the water-wave problem with surface tension.  Communications in Partial Differential Equations, 35(12):2195–2252, 2010.
  * [CS07] Daniel Coutand and Steve Shkoller.  Well-posedness of the free-surface incompressible euler equations with or without surface tension.  Journal of the American Mathematical Society, 20(3):829–930, 2007\.
  * [DIPP17] Y Deng, A. D. Ionescu, B. Pausader, and F. Pusateri.  Global solutions of the gravity-capillary water-wave system in three dimensions.  Acta Mathematica, 219(2):213–402, 2017.
  * [DS04] J-M Delort and Jeremie Szeftel.  Long-time existence for small data nonlinear klein-gordon equations on tori and spheres.  International Mathematics Research Notices, 2004(37):1897–1966, 2004.
  * [GMS12] Pierre Germain, Nader Masmoudi, and Jalal Shatah.  Global solutions for the gravity water waves equation in dimension 3.  Annals of Mathematics, pages 691–754, 2012.
  * [GMS15] Pierre Germain, Nader Masmoudi, and Jalal Shatah.  Global solutions for capillary waves equation in dimension 3.  Comm. Pure Appl. Math, 68(4):625–687, 2015.
  * [HIT16] John K Hunter, Mihaela Ifrim, and Daniel Tataru.  Two dimensional water waves in holomorphic coordinates.  Communications in Mathematical Physics, 346(2):483–552, 2016.
  * [IP15] Alexandru D Ionescu and Fabio Pusateri.  Global solutions for the gravity water waves system in 2d.  Inventiones mathematicae, 199:653–804, 2015.
  * [IP19] AD Ionescu and F Pusateri.  Long-time existence for multi-dimensional periodic water waves.  Geometric and Functional Analysis, 29(3):811–870, 2019.
  * [IT14] Mihaela Ifrim and Daniel Tataru.  Two dimensional water waves in holomorphic coordinates ii: global solutions.  arXiv preprint arXiv:1404.7583, 2014.
  * [Lan05] David Lannes.  Well-posedness of the water-waves equations.  Journal of the American Mathematical Society, 18(3):605–654, 2005\.
  * [MP19] Rolando Magnanini and Giorgio Poggesi.  On the stability for alexandrov’s soap bubble theorem.  Journal d’Analyse Mathématique, 139(1):179–205, 2019.
  * [MZ09] Mei Ming and Zhifei Zhang.  Well-posedness of the water-wave problem with surface tension.  Journal de mathématiques pures et appliquées, 92(5):429–455, 2009.
  * [RT09] Michael Ruzhansky and Ville Turunen.  Pseudo-differential operators and symmetries: background analysis and advanced topics, volume 2.  Springer Science & Business Media, 2009.
  * [Sch05] Ben Schweizer.  On the three-dimensional euler equations with a free boundary subject to surface tension.  In Annales de l’IHP Analyse non linéaire, volume 22, pages 753–781, 2005.
  * [Sha85] Jalal Shatah.  Normal forms and quadratic nonlinear klein-gordon equations.  Communications on Pure and Applied Mathematics, 38(5):685–696, 1985\.
  * [Sog88] Christopher D Sogge.  Concerning the lp norm of spectral clusters for second-order elliptic operators on compact manifolds.  Journal of functional analysis, 77(1):123–138, 1988.
  * [SZ08] Jalal Shatah and Chongchun Zeng.  Geometry and a priori estimates for free boundary problems of the euler’s equation.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 61(5):698–744, 2008.
  * [Wan20] X Wang.  Global regularity for the 3d finite depth capillary water waves.  Annales Scientifiques de l’École Normale Supérieure, 53(4):847–943, 2020.
  * [Wu97] Sijue Wu.  Well-posedness in sobolev spaces of the full water wave problem in 2-d.  Inventiones mathematicae, 130(1):39–72, 1997.
  * [Wu99] Sijue Wu.  Well-posedness in sobolev spaces of the full water wave problem in 3-d.  Journal of the American Mathematical Society, 12(2):445–495, 1999\.
  * [Wu09] Sijue Wu.  Almost global wellposedness of the 2-d full water wave problem.  Inventiones mathematicae, 177(1):45–135, 2009.
  * [Wu11] Sijue Wu.  Global wellposedness of the 3-d full water wave problem.  Inventiones mathematicae, 184(1):125–220, 2011.
  * [Zak68] Vladimir E Zakharov.  Stability of periodic waves of finite amplitude on the surface of a deep fluid.  Journal of Applied Mechanics and Technical Physics, 9(2):190–194, 1968.


#  Ab initio prediction of an order-disorder transition in Mg2GeO4:
implication for the nature of super-Earth’s mantles

Koichiro Umemoto1,2, and Renata M. Wentzcovitch3,4 1Earth-Life Science
Institute, Tokyo Institute of Technology, Tokyo, Japan
2Theoretical Quantum Physics Laboratory, Cluster for Pioneering Research,
RIKEN, Wako-shi, Saitama 351-0198, Japan
3Department of Applied Physics and Applied Mathematics, Columbia University,
New York, NY, USA.
4Department of Earth and Environmental Sciences, Lamont Doherty Earth
Observatory, Columbia University, New York, NY, USA.

###### Abstract

Here we present an ab initio prediction of an order-disorder transition (ODT)
from $I\bar{4}2d$-type to a Th3P4-type phase in the cation sublattices of
Mg2GeO4, a post-post-perovskite (post-PPv) phase. This uncommon type of
prediction is achieved by carrying out a high-throughput sampling of atomic
configurations in a 56-atom supercell followed by a Boltzmann ensemble
statistics calculation. Mg2GeO4 is a low-pressure analog of $I\bar{4}2d$-type
Mg2SiO4, a predicted major planet-forming phase of super-Earths’ mantles.
Therefore, a similar ODT is anticipated in $I\bar{4}2d$-type Mg2SiO4 as well,
which should impact the internal structure and dynamics of these planets. The
prediction of this Th3P4-type phase in Mg2GeO4 enhances further the
relationship between the crystal structures of Earth/planet-forming silicates
and oxides at extreme pressures and those of rare-earth sesquisulfides at low
pressures.

## Introduction

Ab initio quasiharmonic (QHA) calculations of polymorphic phase transitions at
extreme pressure and temperature conditions have proven to be highly
predictive for nearly two decades RIMG2010 . Combined with materials discovery
methods (e.g., Refs. Glass2006 ; Lonie2011 ; Pickard2011 ; Wang2012 ; Wu2014 ;
Curtis2018 ), these simulation tools offer a powerful approach to
investigating phase transition phenomena at challenging experimental
conditions typical of planetary interiors. Mineral physics and geophysics have
benefitted immensely from these developments in materials simulations in the
past two decades. For example, in 2004, ab initio predictions played a crucial
role in discovering and elucidating the major post-perovskite (PPv) transition
in MgSiO3 bridgmanite at deep Earth interior conditions, e.g., 2,500K at 125
GPa Murakami2004 ; Oganov2004 ; Tsuchiya2004 . More recently, ab initio QHA
calculations have explored pressure and temperature conditions expected in
interior of super-Earths, terrestrial-type exoplanets more massive than Earth.
The interest in these planets stems from their similarities and differences
with Earth and their potential habitability. Besides, ab initio methods are
highly predictive when addressing planet-forming silicates and oxides,
motivating experiments. These phases are oxides involving Mg, Si, ferrous and
ferric Fe, Al, and Ca, primarily. In deep interiors of large super-Earths with
up to $\sim$13 Earth mass ($M_{\oplus}$), the range of pressures and
temperatures can reach tens of tera-Pascals (TPa) (hundreds of Mbar) and
$10^{4}$–$10^{5}$ K Hakim2018 ; vandenBerg2019 . Despite remarkable
developments in experimental techniques Smith2014 ; Dubrovinsky2012 ;
Dubrovinskaia2016 ; Dewaele2018 ; Sakai2018 , these pressure-temperature
conditions are still very challenging to experiments making ab initio
predictions of phase transition phenomena in planet-forming phases critical to
advancing planetary modeling. This progress has been registered in a series of
ab initio discoveries concerning the nature of super-Earth’s mantle-forming
phases in the important MgO-SiO2 system Wu2014 ; Umemoto2006a ;
UmemotoWentzcovitch2011 ; Niu2015 ; Umemoto2017 . They have culminated in
partial experimental confirmations in a low-pressure analog system, NaF-MgF2
Dutta2018 , and detailed modeling of these planets’ internal structure and
dynamics Hakim2018 ; vandenBerg2019 . In particular, a sequence of “post-post-
perovksite” (post-PPv) transitions in the MgO-SiO2 system Wu2014 ;
Umemoto2006a ; UmemotoWentzcovitch2011 ; Niu2015 ; Umemoto2017  was predicted
to occur up to $\sim$3 TPa and 10,000 K, starting from MgSiO3 PPv and ending
in its dissociation into the elementary oxides MgO and SiO2. This dissociation
process was predicted to occur in three stages: 1) a dissociation reaction,
MgSiO3 PPv $\to$ $I\bar{4}2d$-type Mg2SiO4 \+ $P2_{1}/c$-type MgSi2O5; 2)
further dissociation of $P2_{1}/c$-type MgSi2O5 $\to$ $I\bar{4}2d$-type
Mg2SiO4 \+ Fe2P-type SiO2; 3) final dissociation of $I\bar{4}2d$-type Mg2SiO4
$\to$ CsCl-type MgO + Fe2P-type SiO2. If MgSiO3 coexists with MgO or SiO2,
other intermediate recombination reactions producing Mg2SiO4 or MgSi2O5 can
intervene in the three-stage dissociation process. These post-PPv transitions
were shown to have profound effects on super-Earths’ mantle dynamics Hakim2018
; vandenBerg2019 . However, the exceedingly high predicted transition
pressures ($>\sim 500$ GPa) make their experimental confirmation quite
challenging.

Two potential low-pressure analog systems, i.e., MgO-GeO2 and NaF-MgF2, were
proposed as viable experimental alternates UmemotoWentzcovitch2015 ;
UmemotoWentzcovitch2019 . Both displayed some of these novel high-pressure
phases found in the MgO-SiO2 system. NaMgF3 PPv was predicted to exhibit the
novel $P2_{1}/c$-type structure of MgSi2O5 under pressure before its full
dissociation into NaF and MgF2. The predicted phases and transformations were
experimentally confirmed Dutta2018 . In contrast, MgGeO3 PPv under pressure
was expected to produce the novel $I\bar{4}2d$-type Mg2SiO4 structure. The
predicted reactions in the MgO-GeO2 system have not been experimentally
confirmed yet since they happen at higher pressures. Also, both systems were
expected to display the respective analog recombination reactions. This state
of affairs brings us to our present study.

Using ab initio techniques, we predict another type of phase transition in the
MgO-GeO2 system, a temperature-induced change from $I\bar{4}2d$-type to
Th3P4-type structure in Mg2GeO4. This is not a regular polymorphic phase
transition but an order-disorder transition (ODT) in the cation sublattices.
The crystal structure of $I\bar{4}2d$-type Mg2GeO4 is shown in Fig. 1. In this
phase, Mg and Ge cations are regularly ordered. Disorder in the cation
sublattices makes the cation sites indistinguishable and results in the
Th3P4-type structure. Details of both crystal structures are given in
Supplementary Information. The result predicted here is unexpected since
disorder occurs between two sublattices containing cations with nominally
different valences. Specifically, Mg and Ge are known to form stoichiometric
end-member phases (MgO and GeO2), preserving their nominal valence at
comparable pressures. This type of prediction is also uncommon since it
requires reliable methods to compute free energy in disordered solid-
solutions, which is computationally much more costly than regular polymorphic
transition. The predicted pressure and temperature transition conditions for
this ODT are more easily achievable in laboratory experiments than the analog
one in Mg2SiO4. Still, a similar ODT is also expected to occur in
$I\bar{4}2d$-type Mg2SiO4 and should be crucial for modeling interiors of
super-Earths.

## Statistical treatment

The ODT critical temperature, $T_{c}$, was calculated using the same approach
previously used to compute the ice-VII to -VIII ODT boundary Umemoto2010 . The
method consists in sampling, if not all, the most significant symmetrically
distinct (irreducible) atomic configurations of a chosen supercell. To
represent the disordered Th3P4-type phase we chose a 56-atom supercell (8
Mg2GeO4, $\sqrt{2}\times\sqrt{2}\times 1$ supercell of the conventional unit
cell) and generated an ensemble of 125 irreducible configurations using the
scheme described below. We then computed the static partition function for
this ensemble of 125 configurations:

$Z_{static}(V,T)=\sum^{125}_{i=1}w_{i}\exp\left(-\frac{E_{i}(V)}{k_{B}T}\right),$
(1)

where $E_{i}(V)$ and $w_{i}$ are the total energy and multiplicity of the
$i$th irreducible configuration ($\sum_{i=1}^{125}w_{i}=3489$), and $k_{B}$ is
the Boltzmann’s constant. The static partition function is then extended to
include zero-point motion (ZPM) and phonon thermal excitation energies within
the QHA Umemoto2010 ; Qin2019 . From $Z(V,T)$, all thermodynamic potentials
and functions can be calculated: Helmholtz free energy $F(V,T)=-k_{B}\ln Z$,
pressure $P(V,T)=-(\partial F/\partial V)_{T}$, Gibbs free energy
$G(V,T)=F+PV$ (converted to $G(P(V,T),T)$), entropy $S(V,T)=-(\partial
F/\partial T)_{V}$, constant-volume heat capacity
$C_{V}(V,T)=-T(\partial^{2}S/\partial T^{2})_{V}$, constant-pressure heat
capacity $C_{P}(P,T)=-T(\partial^{2}S/\partial T^{2})_{P}$, and so forth.
Finally the ODT is obtained by locating a peak in $C_{P}{(T)}$. The procedure
used to generate the 125 cation configurations in Eq. 1 is schematically
depicted in Fig. S1. We started with the ordered structure containing the 24
cations (16 Mg and 8 Si ions) in their respective Wyckoff sites of the
$I\bar{4}2d$-type phase. This lowest enthalpy configuration is shown in Fig.
1(a) and corresponds to the leftmost configuration in Fig. S1. We refer to it
as the “zero-interchange structure”. Then we sequentially interchanged Mg/Ge
pairs once. These single interchanges produce 128 (16 $\times$ 8)
configurations where only 4 are irreducible, each with its distinct
multiplicity. We refer to this first generation of structure as “one-
interchange” configurations. Starting from these 128 configurations, we repeat
this process and produce 3360 configurations among which 120 are irreducible
“two-interchange” configurations. With zero, one, and two cation interchanges,
a total of $1+4+120=125$ configurations were generated and used to compute the
partition function. Discussion of convergence issues related to supercell size
and number of configurations is offered in the Supplementary Information
section.

## Results and discussion

$C_{P}$ profiles at several pressures obtained using static free energy
calculations ($C_{P}^{st}$) are shown in Fig. 2(a). The peak in $C_{P}^{st}$
can be easily identified. Including the vibrational free energy contribution
to the total free energy, a Debye-like contribution is added to $C_{P}^{st}$
producing $C_{P}^{st+vib}$. The peak in $C_{P}^{st+vib}$ appears as a hump in
the Dulong-Petit regime of $C_{P}^{st+vib}$ (Fig. 2(b)). The peak would have
appeared very sharp if the calculation could have been carried out with an
infinitely large supercell and number of configurations. By adding only the
ZPM energy, $E_{ZPM}$, the heat capacity ($C_{P}^{st+zpm}$) still resembles
$C_{P}^{st}$, except that $T_{c}$ lowers by $\sim$100 K at 200 GPa (Fig.
2(c)). This temperature shift is essentially a volume effect caused by the
expansion of the equilibrium volume upon inclusion of $E_{ZPM}$.

With increasing pressure, $T_{c}$, i.e., the peak temperature in $C_{P}$,
increases producing the phase boundary shown in Fig. 2(d). This happens
because the enthalpy difference between all configurations and the ground
state one, the ordered $I\bar{4}2d$-type phase, increases with pressure as
shown in Fig. S3. Among the 125 configurations, the $I\bar{4}2d$-type phase
has the lowest enthalpy at all pressures investigated here. The four
irreducible configurations generated by one-interchange of cations are more
similar to the ordered ground state structure and tend to have lower
enthalpies than the 120 irreducible configurations produced by two-
interchanges, with some exceptions. This behavior of $T_{c}$, i.e., the
positive Clapeyron slope is opposite to that observed in the ice-VII to -VIII
ODT. In the case of ice, all configuration enthalpies converge to a single
value under pressure, that of ice X. This is because the ice ODT precedes and
turns into a hydrogen-bond symmetrization transition under pressure
Umemoto2010 . Besides, in the present ODT $T_{c}$ is only slightly altered by
quantum effects, e.g, ZPM or thermal excitation effects. It takes place above
the Debye temperature, $\theta_{\rm Debye}$, in the Dulong-Petit regime of
$C_{V}$. Therefore, it has a classical origin and can be reasonably well
addressed using the static partition function in Eq. 1.

The following dissociation and recombination post-PPv transitions were
predicted in the MgO-GeO2 system UmemotoWentzcovitch2019 :

Dissociation -

MgGeO3 (PPv) $\to$ Mg2GeO4 ($I\bar{4}2d$-type) + GeO2 (pyrite-type) at $\sim$
175 GPa followed by the transition of GeO2 from pyrite- to cotunnite-type,
which is not affected by the order-disorder transition;

Recombination -

MgGeO3 (PPv) + MgO (B1-type) $\to$ Mg2GeO4 ($I\bar{4}2d$-type) at $\sim$ 173
GPa.

These transitions remain valid at low temperatures, but above the ODT’s,
$T_{c}$ one should replace the $I\bar{4}2d$-type phase by the Th3P4-type one.
The newly computed phase boundaries for these transitions are shown in Fig. 3.
Both reactions have negative Clapeyron slopes, a common behavior in pressure
induced structural transitions involving an increase in cation coordination
Navrotsky1980 . The ODT widens the stability fields of the dissociation
products of MgGeO3 PPv (GeO2 and Mg2GeO4) and of the recombination product
(Mg2GeO4) because the configuration entropy lowers the Gibbs free energy of
Mg2GeO4 in the disordered Th3P4-type phase at higher temperatures. The
configuration entropy also decreases in magnitude the negative Clapeyron
slopes of dissociation and recombination transitions $(dT/dP=\Delta V/\Delta
S)$ at high temperatures. As pointed out earlier UmemotoWentzcovitch2019 ,
hysteresis might prevent experimental observation of these dissociation and
recombination transitions from MgGeO3 PPv, in which case a polymorphic from
PPv to a Gd2S3-type phase might intervene. A similar phenomenon was observed
experimentally in the NaF-MgF2 system in which NaMgF3 PPv transformed to a
U2S3-type post-PPv phase at low temperatures Dutta2018 .

## Conclusion

We have predicted an order-disorder transition (ODT) in the cation sublattices
of the $I\bar{4}2d$-type Mg2GeO4, a post-PPv phase in the Mg-Ge-O system. This
type of prediction is uncommon and is not accomplished simply using modern
materials discovery techniques (e.g., Refs. Glass2006 ; Lonie2011 ;
Pickard2011 ; Wang2012 ; Wu2014 ; Curtis2018 ). In addition to structural
prediction, the ODT prediction requires effective statistical sampling of
atomic configurations Umemoto2010 . Besides, $T_{c}$ cannot be calculated by
direct comparison of Gibbs free energy as in a regular first order transition.
Instead, $T_{c}$ is obtained by calculating the position of a peak in
$C_{P}(T)$ throughout this second order transition. In this study, this was
accomplished using ab initio quasiharmonic (QHA) calculations on a 56-atom
supercell. Although anharmonic effects may play a role in first order
transitions at the high temperatures investigated here, harmonic or anharmonic
vibrational effects play a secondary role in the present ODT.

The predicted $I\bar{4}2d$-type to Th3P4-type phase change expands toward
lower pressures and temperatures the stability fields of the post-PPv
dissociation/recombination products containing Mg2GeO4. The MgO-GeO2 system is
a partly low-pressure analog of the Earth/planet-forming MgO-SiO2 system. Both
$I\bar{4}2d$-type Mg2GeO4 and Mg2SiO4 were predicted to occur as post-PPv
dissociation/recombination transition products. Therefore, the present study
strongly suggests that a similar ODT should also occur in the Mg and Si cation
sublattices of $I\bar{4}2d$-type Mg2SiO4 at the high temperatures typical of
the deep interiors of super-Earths Hakim2018 ; vandenBerg2019 .

Finally, it should be noted that the Th3P4-type structure is a high
temperature form of several rare-earth sesquisulfides, $RR^{\prime}$S3
($R$,$R$’=lanthanoid or actinoid), with vacancies at cation sites
Zachariasen1949 ; Flahaut1979 . For example, Gd2S3 is stabilized in the
orthorhombic $\alpha$ phase at low temperature and transforms to the
Th3P4-type $\gamma$ phase with vacancies at high temperature. As pointed out
earlier Umemoto2008 , the $RR^{\prime}$S3 family of structures form an analog
system to high-pressure phases of MgSiO3 and Al2O3 (bridgmanite, PPv, and
U2S3-type). Hence, the prediction of the Th3P4-type phase in this study
strengthens further the structural relationship between Earth/planet-forming
phases at ultrahigh pressures and rare-earth sesquisulfides.

## Computational details

Calculations for the 125 cation configurations were performed using the local-
density approximation Perdew1981  to density-functional theory. For all atomic
species, Vanderbilt-type pseudopotentials Vanderbilt1990  were generated. The
valence electron configurations and cutoff radii for the pseudopotentials were
$2s^{2}2p^{6}3s^{2}$ and 1.6 a.u. for Mg, $4s^{2}4p^{1}3d^{10}$ and 1.6 a.u.
for Ge, and $2s^{2}2p^{4}$ and 1.4 a.u. for O, respectively. Cutoff energies
for the plane-wave expansion are 70 Ry. The $2\times 2\times 2$ k-point mesh
was used for the 56-atom supercell. For structural optimization under
arbitrary pressures between 100 and 800 GPa, we used the variable-cell-shape
damped molecular dynamics Wentzcovitch1991 ; Wentzcovitch1993 . Dynamical
matrices were calculated on the $2\times 2\times 2$ q mesh using density-
functional perturbation theory Giannozzi1991 ; Baroni2001 . The vibrational
contribution to the partition function was taken into account within the
quasi-harmonic approximation (QHA) Wallace1972 ; Umemoto2010 ; Qin2019 . The
$8\times 8\times 8$ q-point mesh was used for the QHA summation. All
calculations were performed using the Quantum-ESPRESSO Giannozzi2009  and the
qha software Qin2019 .

## Acknowledgments

K.U. acknowledges support of a JSPS Kakenhi Grant # 17K05627. R.M.W.
acknowledges support of a US Department of Energy Grant DE-SC0019759. All
calculations were performed at Global Scientific Information and Computing
Center and in the ELSI supercomputing system at the Tokyo Institute of
Technology, the HOKUSAI system of RIKEN, and the Supercomputer Center at the
Institute for Solid State Physics, the University of Tokyo.

## References

  * (1) Wentzcovitch RM, Stixrude S (eds.) (2010) Theoretical and computational methods in mineral physics: Geophysical Applications, Reviews in Mineralogy and Geochemistry, 71.
  * (2) Glass CW, Oganov AR, Hansen N (2006) USPEX – Evolutionary crystal structure prediction, Comp. Phys. Commun. 175, 713-720.
  * (3) Lonie DC, Zurek E (2011) XTALOPT: An open-source evolutionary algorithm for crystal structure prediction, Comp. Phys. Commun. 175, 713-720.
  * (4) Pickard CJ, Needs RJ (2011) Ab initio random structure searching, J. Phys. Condens. Matter 23, 053201.
  * (5) Wang Y, Lv J, Zhu L, Ma Y (2012) CALYPSO: A method for crystal structure prediction, Comp. Phys. Commun. 183, 2063-2070.
  * (6) Wu SQ, et al. (2014) An adaptive genetic algorithm for crystal structure prediction,J. Phys.: Condens. Matter 26, 035402.
  * (7) Curtis F, et al. (2018) GAtor: A First-Principles Genetic Algorithm for Molecular Crystal Structure Predction, J. Chem. Theory Comput. 14, 2246-2264.
  * (8) Murakami M, Hirose K, Kawamura K, Sata N, Ohishi Y (2004) Post-perovskite phae transition in MgSiO3. Science 304, 855-858.
  * (9) Oganov AR, Ono S (2004) Theoretical and experimental evidence for a post-perovsite phase of MgSiO3 in Earth’s D” layer, Nature 430, 445-448.
  * (10) Tsuchiya T, Tsuchiya J, Umemoto K, Wentzcovitch RM (2004) Phase transition i MgSiO3 perovskite in the earth’s lower mantle, Earth Planet. Sci. Lett. 224, 241-248.
  * (11) Hakim K, et al. (2018) A new ab initio equation of state of hcp-Fe and its implication on the interior tructure and mass-radius relations of rocky super-Earths, Icarus 313, 61-78.
  * (12) van den Berg AP, Yuen DA, Umemoto K, Jacobs MHG, Wentzcovich RM (2019) Mass-dependent dynamics of terrestrial exoplanets using ab initio mineral prperties, Icarus 317, 412-426.
  * (13) Smith RF, et al. (2014) Ramp compression of diamond to five terapascals, Nature 511, 330-333.
  * (14) Dubrovinsky L, Dubrovinskaia N, Prakapenka VB, Abakumov AM (2012) Implementation of micro-ball nanodiamond anvils for high-pressure studies above Mbar, Nat. Commun. 3, 1163.
  * (15) Dubrovinskaia N, et al. (2016) Terapascal static pressure generation with ultrahigh yield strength nanodiamond, Sci Adv. 2, e1600341.
  * (16) Dewaele A, Loubeyre P, Occelli F, Marie O, Mezouar M (2018) Toroidal diamond anvil cell for detailed measurements under extreme static pressres, Nat. Commun. 9, 2913.
  * (17) Sakai T, et al. (2018) High pressure generation using double-stage diamond anvil technique: problems an equations of state of rhenium, High Pressure Res. 38, 107-119.
  * (18) Umemoto K, Wentzcovitch RM, Allen PB (2006) Dissociation of MgSiO3 in th Cores of Gas Giants and Terrestrial Exoplanets, Science 311, 983-986.
  * (19) Umemoto K, Wentzcovitch RM (2011) Two-stage dissociation in MgSiO3 post-peroskite, Earth Plenet. Sci. Lett. 311, 225-229.
  * (20) Niu H, Oganov AR, Chen XQ, Li D (2015) Prediction of novel stable compound in the Mg-Si-O system under exoplanet pressures, Sci. Rep. 5, 18347.
  * (21) Umemoto K, et al. (2017) Phase tansitions in MgSiO3 post-perovskite in super-Earth mantles, Earth Planet. Sci Lett. 478, 40-45.
  * (22) Dutta R, Greenberg E, Prakapenka VB, Duffy TS (2019) Phase Transitions beynd post-perovskite in NaMgF3 to 160 GPa, Proc. Nat. Acad. Sci. 116, 1934-19329.
  * (23) Umemoto K, Wentzcovitch RM (2015) Two-stages Dissociation of NaMgF3 Post-Pervskite: A Potential Low-Pressure Analog of MgSiO3 at Multi-Mbar Pressures, JP Conf. Proc. 4, 011002.
  * (24) Umemoto K, Wentzvocitch RM (2019) Ab initio exploration of post-PPv transitionsin low-pressure analogs of MgSiO3, Phys. Rev. Mat. 3, 123601.
  * (25) Umemoto K, Wentzcovitch RM, de Gironcoli S, Baroni S (2010) Order–disorder phase boundary between ice VII and VIII obtained by first princiles, Chem. Phys. Lett. 499, 236-240.
  * (26) Qin T, Zhang Q, Wentzcovitch RM, Umemoto K (2019) qha: A Python package for quasiharmonic free energy calculation for multi-configration systems, Comp. Phys. Commun. 237, 199-207.
  * (27) Navrotsky A (1980) Lower mantle phase transitions may generally have negative pressur-temperature slople, Geophys. Res. Lett. 7, 709-711.
  * (28) Zachariasen WH (1949) Crystal Chemical Studies of the $5f$-Series of Elements. VI. he Ce2S3-Ce3S4 Type of Structure, Acta Cryst. 2, 57-60.
  * (29) Flahaut J (1979) Sulfides, selenides, and tellurides. Handbook on the Physics and Cheistry of Rare-Earths, eds. Gschneidner K. A. Jr. and Eyring L. R. (North-Holland Amsterdam), 4, 1.
  * (30) Umemoto K, Wentzcovitch RM (2008) Prediction of an U2S3-type polymorph of Al2O3 at 3.7 Mbar, Proc. Nat. Acad. Sci. 105, 6526-6530.
  * (31) Perdew JP, Zunger A (1981) Self-interaction correction to density-functional aproximations for many-electron systems, Phys. Rev. B, 23, 5048.
  * (32) Vanderbilt D (1990) Soft self-consistent pseudopotentials in a generalized eigenvalueformalism, Phys. Rev. B, 41, R7892.
  * (33) Wentzcovitch RM (1991) Invariant molecular-dynamics approach to structural phase trnsitions, Phys. Rev. B 44, 2358.
  * (34) Wentzcovitch RM, Martins JL, Price GD (1993) Ab Initio Molecular Dynamicswith Variable Cell Shape: Application to MgSiO3, Phys. Rev. Lett. 70, 347.
  * (35) Baroni S, de Gironcoli S, Dal Corso A, Giannozzi P (2001) Phonons and related crystal properties from density-functional perturbation theoy, Rev. Mod. Phys. 73, 515.
  * (36) Giannozzi P, de Gironcoli S, Pavone P, Baroni S (1991) Ab initio calculation o phonon dispersions in semiconductors, Phys. Rev. B 43, 7231.
  * (37) Wallace D (1972) Thermodynamics of Crystals, John Wiley, Hoboken, N. J..
  * (38) Giannozzi P, et al. (2009) Quantum ESPRESSO: a modular and open-source software projet for quantum simulations of materials, J. Phys.: Condens. Matter, 21, 39502.

Figures

Figure 1:  (a) Crystal structure of $I\bar{4}2d$-type Mg2GeO4. Yellow, grey,
and red spheres denote Mg, Ge, and O ions, respectively. (b) Coodination
polyhedra around Mg and Ge ions. Numbers next to the O atoms represent Mg-O
and Ge-O bond-lengths in Angstroms at 400 GPa.

Figure 2:  Constant-pressure heat capacities, $C_{P}$, calculated using 125
configurations in the 56-atom supercell and (a) static free energies only (no
vibrational effects included) ($C_{P}^{st}$), (b) free energies calculated
using the QHA ($C_{P}^{st+vib}$), and (c) static free energy plus $E_{ZPM}$
($C_{P}^{st+zpm}$). Vertical dashed lines indicate the peak temperatures of
$C_{P}^{st+zpm}$. (d) Phase boundary of the ODT between $I\bar{4}2d$-type and
Th3P4-type Mg2GeO4 obtained from the peak temperatures in (c).

Figure 3:  Phase transition boundaries in the MgO-GeO2 system. The blue solid
lines show the ODT in Mg2GeO4 given by the peak temperature in $C_{P}$ (Fig.
2(c)). The red lines denote the post-PPv phase boundaries involving Th3P4-type
Mg2GeO4, while the black lines denote those involving $I\bar{4}2d$-type
Mg2GeO4 calculated in Ref. UmemotoWentzcovitch2019 . The dashed lines
represent the metastable continuation of these phase boundaries.