G. N. Koutsokostas et al

*D. J. Frantzeskakis, Department of Physics, National and Kapodistrian University of Athens, Panepistimiopolis, Zografos, Athens 15784, Greece.

# Nonlinear instability and solitons in a self-gravitating fluid

G. N. Koutsokostas  S. Sypsas  O. Evnin  T. P. Horikis  D. J. Frantzeskakis
Department of Physics, National and Kapodistrian University of Athens, Athens,
Greece  High Energy Physics Research Unit, Faculty of Science, Chulalongkorn
University, Bangkok, Thailand  NARIT, Don Kaeo, Mae Rim, Chiang Mai, Thailand
Theoretische Natuurkunde, Vrije Universiteit Brussel and International Solvay
Institutes, Brussels, Belgium  Department of Mathematics, University of
Ioannina, Ioannina, Greece<EMAIL_ADDRESS>
(<day> <Month>, <year>; <day> <Month>, <year>; <day> <Month>, <year>)

###### Abstract

[Abstract]

We study a spherical, self-gravitating fluid model, which finds applications
in cosmic structure formation. We argue that since the system features
nonlinearity and gravity-induced dispersion, the emergence of solitons becomes
possible. We thus employ a multiscale expansion method to study, in the weakly
nonlinear regime, the evolution of small-amplitude perturbations around the
equilibrium state. This way, we derive a spherical nonlinear Schrödinger (NLS)
equation that governs the envelope of the perturbations. The effective NLS
description allows us to predict a “nonlinear instability” (occurring in the
nonlinear regime of the system), namely, the modulational instability which,
in turn, may give rise to spherical soliton states. The latter feature a very
slow (polynomial) curvature-induced decay in time. The soliton profiles may be
used to describe the shape of dark matter halos at the rims of the galaxies.

††articletype: Article Type

##  1 Introduction

A set of strongly-interacting particles on scales much larger than their mean
free-path, i.e., where the motion of each individual particle relative to its
neighbors is negligible, can be described via conservation laws. Within this
description, the collective behavior of the system will be governed by the
equations of fluid mechanics.

One of the simplest possible cases in this context is an isentropic one-
dimensional (1D) flow, described by the following system of partial
differential equations (PDEs) —see, e.g., Refs.1, 2: the continuity equation
(conservation of mass),

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x}(\rho
u)=0,$ (1)

where $\rho(x,t)$ is the mass density and $u(x,t)$ is the fluid velocity, as
well as the Euler equation (conservation of momentum),

$\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}=-\frac{1}{\rho}\frac{\partial p}{\partial x},$ (2)

where $p(x,t)$ is the pressure. The latter is connected to the density via an
equation of state,

$p=p(\rho),$ (3)

where it is assumed that $dp/d\rho=c_{s}^{2}(\rho)>0$, with $c_{s}$ being the
(generally density-dependent) speed of sound. The nonlinear system (1)-(3) can
be expressed in the following conservation form:

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x}(\rho
u)=0,$ (4) $\displaystyle\frac{\partial}{\partial t}(\rho
u)+\frac{\partial}{\partial x}\left[\rho u^{2}+p(\rho)\right]=0,$ (5)

which is in fact identical to the compressible isentropic gas dynamical
equations, and can be solved by the ingenious method introduced by Riemann 3.
Indeed, the hyperbolic Eqs. (4)-(5) can be diagonalized by the Riemann
invariants:

$\displaystyle
r_{1}=u-\int^{\rho}\frac{c_{s}(\rho^{\prime})}{\rho^{\prime}}d\rho^{\prime},\quad
r_{2}=u+\int^{\rho}\frac{c_{s}(\rho^{\prime})}{\rho^{\prime}}d\rho^{\prime},$
(6)

with characteristic velocities

$v_{1}=u-c_{s}(\rho),\quad v_{2}=u+c_{s}(\rho),$ (7)

so that

$\frac{\partial r_{j}}{\partial t}+v_{j}\frac{\partial r_{j}}{\partial
x}=0,\quad j=1,2.$ (8)

Due to the monotonicity 111Generally, for a polytropic gas,
$p(\rho)=A\rho^{\gamma}$ (where $\gamma$ is the polytropic index and $A={\rm
const.}$), one finds $g(\rho)=[2/(\gamma-1)]c_{s}(\rho)$, with $c_{s}(\rho)$
being a monotone function. of
$g(\rho)=\int^{\rho}\left(c_{s}(\rho^{\prime})/\rho^{\prime}\right)d\rho^{\prime}$,
Eqs. (6) can be inverted leading to:

$u=\frac{1}{2}(r_{1}+r_{2}),\quad\rho=g^{-1}\left(\frac{1}{2}(r_{2}-r_{1})\right).$
(9)

The above purely nonlinear hyperbolic system is a prototypical setting where
smooth initial data may develop singularities in finite time, i.e., shock
waves. In this regard, the existence of Riemann invariants allows for
estimates of the breaking times at which these may occur4.

Nevertheless, apart from the intrinsic nonlinearity which may give rise to
purely nonlinear waves in the form of shocks, the fluid model under
consideration may also display the important phenomenon of dispersion if
gravity is included. Indeed, let us follow Ref.5, 6, and introduce a 1D self-
gravitating gas, such that the Euler equation is modified as

$\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}=-\frac{1}{\rho}\frac{\partial p}{\partial x}-\frac{\partial\Phi}{\partial
x},$ (10)

where $\Phi$ is the gravitational potential obeying the Poisson equation:

$\frac{\partial^{2}\Phi}{\partial x^{2}}=4\pi G\left(\rho-\rho_{0}\right).$
(11)

Here, $G$ is Newton’s gravitational constant and $\rho_{0}$ is an equilibrium
density. Since we are interested in fluctuations propagating on top of the
equilibrium background density, let us define the (linear) overdensity as
$\rho_{1}\equiv(\rho-\rho_{0})/\rho_{0}$. Then, the pertinent PDEs describing
this self-gravitating gas, namely Eqs. (1), (10), (11), together with the
state-equation (3), when linearized around the equilibrium state, admit plane
wave solutions $\rho_{1}\propto\exp[i(kx-\omega t)]$, with the frequency
$\omega$ and the wavenumber $k$ satisfying the dispersion relation

$\omega^{2}(k)=c_{s}^{2}(\rho_{0})\left(k^{2}-k_{J}^{2}\right),$ (12)

where $k_{J}=(4\pi G\rho_{0})^{1/2}/c_{s}(\rho_{0})$ is the so-called Jeans
wavenumber. Since $\partial^{2}\omega(k)/\partial k^{2}\neq 0$, the nonlinear
system is also dispersive2, a fact that suggests the possibility of soliton
formation.

Apart from highlighting the presence of gravity-induced dispersion, the
relation (12) also suggests the occurrence of the (linear) Jeans instability
7: perturbations of wavenumbers $|k|<k_{J}$ are unstable and grow
exponentially in time. Thus, at scales larger than the Jeans length, where
gravity dominates the internal pressure, overdensities collapse under their
gravitational pull leading to the emergence of structures.

Recently, there has been a substantial increase of interest in self-
gravitating fluid models, as they generically appear in studies in
astrophysics and plasmas 8, as well as in cosmology, and particularly in the
context of cold dark matter —see, e.g., Refs. 9, 10, 11, 12. In such settings,
obviously, pertinent self-gravitating fluid systems are considered in higher-
dimensions (and in many cases in a spherical geometry). However, once the
dynamics is considered in the fully nonlinear regime, due to their complicated
character, multidimensional self-gravitating fluid models are usually studied
by means of numerical simulations —see, e.g., Refs. 12, 13.

Here, motivated by these works, our scope is to study analytically a radial
three-dimensional (3D) fluid model [see Eqs. (14)-(16) and (31) below]. Upon
employing a multiscale expansion method we are able to approximate the system,
under certain conditions, with a much more tractable equation. In particular,
we show that, in the weakly nonlinear regime, small-amplitude excitations
around the equilibrium state of the fluid are described, at slow length/time
scales, by an effective spherical nonlinear Schrödinger (NLS) equation of the
form:

$i\frac{\partial
A}{\partial\tilde{t}}+i\frac{A}{\tilde{t}}+\frac{1}{2}\frac{\partial^{2}\omega(k)}{\partial
k^{2}}\frac{\partial^{2}A}{\partial\tilde{r}^{2}}+\gamma(k)|A|^{2}A=0,$ (13)

where $A$ is the envelope of the perturbations, $\tilde{t}$ and $\tilde{r}$
denote proper slow time and slow radial variable, $\omega(k)$ is given by Eq.
(12), and $\gamma(k)$ is the nonlinearity coefficient, which we specify below.
Notice that the spherical NLS (13), which has been used to describe radial
solitons in plasmas 14, 15 and fluids 16, is actually a traditional NLS
equation, which also incorporates the curvature term $iA/\tilde{t}$ (in fact,
as we will see, $r^{-1}\propto\tilde{t}^{-1}$) stemming from the consideration
of the radial $3$D geometry. This term suggests that, although the initial
system is conservative, the underlying wave structures feature a curvature-
induced decay (see details in Sec. 5 below). This effect is reminiscent to the
situation occurring, e.g., in the usual radially symmetric 2nd-order wave
equation, whose spherical (D’Alembert) wave solutions feature a curvature-
induced decay $\propto 1/r$.

With the effective spherical NLS model at hand, we are able to deduce the
following. First, we investigate a certain “nonlinear instability”, i.e., one
occurring in the nonlinear regime of the system, the so-called modulational
instability (MI) 17 that can manifest itself in addition to the Jeans
instability. This way, we find that the weak, spatially uniform and time-
varying solution of the NLS is —under certain conditions— subject to MI, for
wavenumbers outside the Jeans instability band, namely for $k_{J}<|k|<\alpha
k_{J}$, with $\alpha\approx 1.34$. As mentioned, this is an example of
nonlinearity (i.e., the nonlinearity-induced MI) competing with the gravity-
induced dispersion, giving rise to the formation of solitons. The latter may
appear as regular sech-shaped solitons18 or rogue waves19, 20 and they are
characterized by an exponential or an algebraic decay in space, respectively,
while both experience a slow curvature- induced polynomial decay in time.

Physically, our results could be relevant for cosmic structure formation. The
system (1), (10), (11), extended to three dimensions is a fair description of
the universe during its (dark) matter-dominated era. Small density
fluctuations along the initially homogeneous background, provide the seeds
that grow to form local dark-matter overdensities, objects which then act as
gravitational wells that attract (luminous) baryonic matter, eventually
forming the galaxies we observe. The radial distribution of dark matter (DM)
in these galactic halos is inferred by the rotational motion of their baryonic
component21 and seems to hint upon a layered structure: a central, smooth,
high-density core followed by a decaying tail22. In pure dark matter models,
the presence of a core can be attributed to the Jeans instability supported by
pressure23. In this context, the main claim of the present work is the
appearance of MI-induced solitonic structures at scales larger than the Jeans
length; for a cored DM halo this might appear as a slightly overdense outer
shell.

The paper is organized as follows. First, in Sec. 2 we introduce the fluid
model, and study the linear regime, where the Jeans instability occurs. In
Sec. 3, using a multiscale expansion method, we derive the effective spherical
NLS equation. In Sec. 4, we analyze the nonlinear instability (MI) that is
predicted in the framework of the NLS, while in Sec. 5 we present spherical
soliton solutions that can be formed due to MI. Finally, in Sec. 6 we
summarize our findings and discuss possibilities for future work.

##  2 The model and Jeans instability

We consider the dynamical equations for an ideal, multi-dimensional fluid, in
the context of non-relativistic Einstein gravity. These include the continuity
equation:

$\frac{\partial\rho}{\partial t}+\mathbf{\nabla}\cdot(\rho\mathbf{u})=0,$ (14)

where $\rho(\mathbf{r},t)$ and $\mathbf{u}(\mathbf{r},t)$ denote the mass
density and velocity, respectively, as well as the Euler equation:

$\frac{\partial\mathbf{u}}{\partial
t}+\left(\mathbf{u}\cdot\mathbf{\nabla}\right)\mathbf{u}=-\frac{1}{\rho}\mathbf{\nabla}p-\mathbf{\nabla}\Phi.$
(15)

Here, the gravitational potential $\Phi(\mathbf{r},t)$ is coupled to the mass
density via the Poisson equation,

$\Delta\Phi=4\pi G\left(\rho-\rho_{0}\right),$ (16)

where, as above, $G$ is Newton’s gravitational constant, $\rho_{0}$ is the
equilibrium density and $\Delta=\nabla^{2}$. Finally, the pressure $p$ in Eq.
(19) is related to the mass density $\rho$ by an equation of state,
$p=p(\rho,T)$, with $T$ being the temperature. We consider the case of
constant temperature and adopt a phenomenological model, where the pressure
and the density obey the following equation:

$p(\rho)=c_{s}^{2}\rho,$ (17)

with $c_{s}$ being the speed of isothermal sound, which is assumed to be
constant in what follows. This assumption corresponds to Boyle’s law for ideal
gases, and is relevant to the considered dark matter problem for the following
reasons. First, we note that the state equation (17) corresponds to the case
of a spherical polytropic gas, and suggests that $c_{s}^{2}=kT/m$, where $k$
is the Boltzmann’s constant, $T$ is the temperature of the dark matter
(assumed to be constant and very low), and $m$ is the mass of the dark matter
particles. In the standard cosmological model, known as Lambda-cold-dark-
matter ($\Lambda$CDM), dark matter is assumed to consist of cold, non-
interacting particles, the mass of which is currently unknown. In a fluid
language, “cold” translates to a vanishing sound speed and, consequently, the
Jeans length vanishes too. However, a speed of sound can be attributed24 due
to (isotropic) velocity dispersion: $c_{s}=\sqrt{\langle v_{i}^{2}\rangle}$,
which suggests that $c_{s}=$ const. On the other hand, another possibility
with a non-vanishing Jeans length is the so-called fuzzy dark matter25, which
considers the dynamics of a nonrelativistic, ultra-light scalar field coupled
to gravity. In fact, our qualitative results depend only on the existence of
the Jeans length and are realization independent.

It is convenient for our analysis to express the model (14)-(16) and (17) in
dimensionless form. This can be done upon measuring the density $\rho$ and
pressure $p$ in units of their equilibrium values $\rho_{0}$ and
$p_{0}=p(\rho_{0})$, and normalize the velocity $\mathbf{u}$, the
gravitational potential $\Phi$, the space variable $\mathbf{r}$ and time $t$
with respect to $c_{s}$, $c_{s}^{2}$, $(p_{0}/4\pi G\rho_{0}^{2})^{1/2}$ and
$(4\pi G\rho_{0})^{-1/2}$, respectively. In what follows, we set the sound
speed to $c_{s}=10^{-6}c_{0}=3\times 10^{2}\leavevmode\nobreak\ {\rm m/s}$,
with $c_{0}$ the speed of light; such a value is compatible with
observations26. Then, the above units take the following values: the
equilibrium density is $\rho_{0}\approx 2.5\times 10^{-27}\leavevmode\nobreak\
{\rm kg/m^{3}}$, while the characteristic length and time scales are
$7.16\times 10^{22}\leavevmode\nobreak\ {\rm m}\approx 7\leavevmode\nobreak\
{\rm kpc}$ and $6.9\times 10^{17}\leavevmode\nobreak\ {\rm s}\approx 2.2\times
10^{10}\leavevmode\nobreak\ {\rm Julian\leavevmode\nobreak\ years}$ [the
gravitational constant is $G=6.67\times 10^{-11}\leavevmode\nobreak\ {\rm
m^{3}/(kg\cdot s^{2})}]$. As we will show later, the structures we will be
focusing on have characteristic sizes a few times larger than the Jeans
radius, placing them at the edge of typical galactic halos. Finally, as per
common practice, since the gravitationally bound objects we will be discussing
decouple from the expanding background, in what follows we will neglect the
cosmological expansion and treat spacetime as static.

Under this assumption and in these units, the considered fluid model is
written in the following dimensionless form:

$\displaystyle\frac{\partial\rho}{\partial
t}+\mathbf{\nabla}\cdot(\rho\mathbf{u})=0,$ (18)
$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\mathbf{\nabla})\mathbf{u}=-\frac{1}{\rho}\mathbf{\nabla}p-\mathbf{\nabla}\Phi,$
(19) $\displaystyle\Delta\Phi=\rho-1,$ (20) $\displaystyle p(\rho)=\rho.$ (21)

A simple nontrivial solution of the above system is

$\rho=1,\quad p=1,\quad\mathbf{u}=0,\quad\Phi=0,$ (22)

which corresponds to a stationary background solution.

To study the stability of the background, we introduce the perturbation ansatz

$\rho=1+\epsilon\rho_{1},\quad p=1+\epsilon
p_{1},\quad\mathbf{u}=\epsilon\mathbf{u}_{1},\quad\Phi=\epsilon\Phi_{1},$ (23)

where $0<\epsilon\ll 1$ is a formal small parameter denoting the amplitude of
fluctuations with respect to the background, while $\rho_{1}$, $p_{1}$,
$\mathbf{u}_{1}$, and $\Phi_{1}$ are arbitrary perturbations of order ${\cal
O}(1)$. Substituting Eqs. (23) into Eqs. (18)-(21) we obtain at ${\cal
O}(\epsilon)$ the following linear system:

$\displaystyle\frac{\partial\rho_{1}}{\partial
t}+\mathbf{\nabla}\cdot\mathbf{u}_{1}=0,$ (24)
$\displaystyle\frac{\partial\mathbf{u}_{1}}{\partial
t}=-\mathbf{\nabla}p_{1}-\mathbf{\nabla}\Phi_{1},$ (25)
$\displaystyle\Delta\Phi_{1}=\rho_{1},$ (26) $\displaystyle p_{1}=\rho_{1},$
(27)

which can be reduced to a single equation for $\rho_{1}$, namely:

$\frac{\partial^{2}\rho_{1}}{\partial t^{2}}-\Delta\rho_{1}-\rho_{1}=0.$ (28)

Equation (28) admits plane wave solutions, namely
$\rho_{1}\propto\exp[i(\mathbf{k}\cdot\mathbf{r}-\omega t)]+{\rm c.c.}$ (where
c.c. denotes complex conjugate), with the perturbation’s frequency $\omega$
and the amplitude $k=|\mathbf{k}|$ of the wavevector $\mathbf{k}$ obeying the
dispersion relation:

$\omega^{2}(k)=k^{2}-k_{J}^{2},$ (29)

where the parameter $k_{J}$ is the Jeans wavenumber and takes the value
$k_{J}=1$ in our units (recall that, in physical units, the Jeans wavenumber
reads $k_{J}=(4\pi G\rho_{0})^{1/2}/c_{s}(\rho_{0})$, as mentioned in the
Introduction). Importantly, the dispersion relation (29) implies that for
perturbation wavenumbers less than $k_{J}$, i.e., for $|k|<k_{J}$, the
frequency $\omega(k)$ becomes purely imaginary and, hence, the perturbation
behaves as $\exp[{\rm Im}\,\omega(k)t]\exp(i\mathbf{k}\cdot\mathbf{r})$ (the
other mode decays exponentially in time), featuring a growth rate:

${\rm Im}\,\omega(k)=\sqrt{1-k^{2}}.$ (30)

In such a case, the system is unstable, and the relevant instability is the
so-called Jeans instability. Notice that the growth rate becomes maximum for
$k=0$ (i.e., for infinite wavelengths) and its value is ${\rm
Im}\,\omega(0)=1$.

Figure 1: The Jeans instability growth rate ${\rm Im}\,\omega(k)$ as a
function of the wavenumber $k>0$. According to the linear theory, the
instability band is defined by the interval $k<k_{\rm J}=1$, highlighted with
the blue segment. The red segment depicts the extension of the instability
band once nonlinearity is taken into account —see discussion in Sec. 4.

The growth rate, ${\rm Im}\,\omega(k)$, as a function of $k$ is shown in Fig.
1. Importantly, the growth rate drops to zero for $k=k_{J}=1$, implying that
the Jeans instability ceases to exist for $|k|\geq 1$ (indeed, as seen from
Eq. (29), in this case $\omega(k)\in\mathbb{R}\leavevmode\nobreak\ \forall
k$). Nevertheless, as we will show in the next section, an analysis of the
(weakly) nonlinear regime of the problem reveals that a nonlinear instability
may occur for $k\gtrsim 1$. The relevant extension of the instability band
beyond $k_{J}$ is depicted as a red horizontal line in Fig. 1 (see details
below).

##  3 Derivation of the effective spherical NLS equation

The above discussion on the Jeans instability actually stems from the
consideration of the linearized version of the self-gravitating fluid model.
Our scope in this section is to show that it is also possible to study
analytically the full, nonlinear problem. The latter, apart from being
nonlinear, is also dispersive, since the dispersion relation (29) implies that
$\partial^{2}\omega(k)/\partial k^{2}\neq 0$. It is then our intention to
employ the method of multiple scales27 to asymptotically reduce the original
self-gravitating fluid model to a simpler nonlinear evolution equation,
namely, a spherical NLS equation, that admits soliton solutions.

We start by noticing that, in the case of a spherical system under
consideration, it is natural to seek for solutions of Eqs. (18)-(21) with
angular symmetry; in such a case, the Laplacian in Eq. (20) may be expressed
as

$\Delta=\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial r},$ (31)

while the fluid velocity reduces to $\mathbf{u}=(u,0,0)$. Notice that in this
work we focus on the simplest nontrivial setup: a spherically symmetric
overdensity. Here, $r=0$ represents the center of the object that will finally
be formed. In a realistic scenario, there will be a collection of such objects
and, hence, there, the spherical symmetry assumption breaks down. Strictly
speaking, we are going to solve the problem on a ball of radius $r_{0}$ which
is much larger than the radius of the object in question and much smaller than
the distance between two objects.

In the wavenumber regime $k>k_{J}=1$, where the system is still not prone to
the Jeans instability, we seek solutions of the system (18)-(20) in the form
of the following asymptotic expansions:

$\displaystyle\rho=1+\sum_{j=1}^{\infty}\epsilon^{j}\rho_{j}(r,t,R,T_{2}),\quad
u=\sum_{j=1}^{\infty}\epsilon^{j}u_{j}(r,t,R,T_{2}),\quad\Phi=\sum_{j=1}^{\infty}\epsilon^{j}\Phi_{j}(r,t,R,T_{2}).$
(32)

Here, the unknown functions, $\rho_{j}$, $u_{j}$ and $\Phi_{j}$, depend on the
fast variables $r$ and $t$, as well as the slow variables $R_{1}$, $T_{1}$,
and $T_{2}$, which are defined as follows:

$R=R_{1}-VT_{1},\quad R_{1}=\epsilon(r-r_{0}),\quad T_{1}=\epsilon t,\quad
T_{2}=\epsilon^{2}t,$ (33)

where $r_{0}$ is the radius of the initial configuration (e.g., the initial
radius of a spherical soliton —see below), while the unknown velocity $V$ in
the definition of $R$ will be determined in a self- consistent manner. Notice
that $R=R_{1}-VT_{1}=\epsilon(r-Vt-r_{0})$ is actually a traveling wave
coordinate associated with the expansion ($V>0$) or contraction ($V<0$) of the
sought spherical structure. As such, it describes the increase or decrease of
the spherical structure’s radius with time. The domain of the variable $R$ is
$-\epsilon(r_{0}+Vt)<R<+\infty$ (since $r>0$), and hence it can take both
positive and negative values, for radii outside or inside the sphere of radius
$r_{0}+Vt$. Furthermore, the above unknown fields are assumed to satisfy
vanishing boundary conditions as $r\rightarrow\infty$.

It is also relevant to make the following remarks. First, since the system
features quadratic nonlinearities, the above slow variables are the only ones
that appear at the third order of approximation [i.e., at ${\cal
O}(\epsilon^{3})$], where the first nonlinear equation emerges —as we will see
this will be the NLS equation. In addition, the choice of the slow traveling
wave coordinate $R$ is made to facilitate the calculations at ${\cal
O}(\epsilon^{2})$ and ${\cal O}(\epsilon^{3})$ —see below. Furthermore, it is
noticed that Eqs. (33) suggest that $1/r\sim\epsilon^{2}/(VT_{2})$, which
means that curvature effects will become important at ${\cal
O}\\!\left(\epsilon^{3}\right)$. Finally, it should be mentioned that our
analysis is similar to the one in the $(1+1)$-dimensional setting 6;
nevertheless, in that work, the reductive perturbation method —instead of the
multiple scales technique [see Eqs. (33)]— is used, while the important, to
our spherical setting, above mentioned curvature term was not taken into
regard. As we will see, the inclusion of this term will lead us to a spherical
NLS equation, rather the traditional one appearing in the $1$D setting.

Substituting the expansions (32) into Eqs. (18)-(20) we obtain, at each order
in $\epsilon$, the following system:

$\displaystyle\frac{\partial\rho_{j}}{\partial t}+\frac{\partial
u_{j}}{\partial r}=f_{j},$ (34) $\displaystyle\frac{\partial u_{j}}{\partial
t}+\frac{\partial\rho_{j}}{\partial r}+\frac{\partial\Phi_{j}}{\partial
r}=g_{j},$ (35) $\displaystyle\frac{\partial^{2}\Phi_{j}}{\partial
r^{2}}-\rho_{j}=h_{j},$ (36)

where $j=1,2,\dots,$ denotes the order of the perturbation scheme, while the
inhomogeneous parts $f_{j}$, $g_{j}$ and $h_{j}$ will be defined below.
Similarly to the analysis of the linearized version of the model, it is
possible to derive from Eqs. (34)-(36) a single equation for the fields
$\rho_{j}$, namely:

$\frac{\partial^{2}\rho_{j}}{\partial
t^{2}}-\frac{\partial^{2}\rho_{j}}{\partial r^{2}}-\rho_{j}=\frac{\partial
f_{j}}{\partial t}-\frac{\partial g_{j}}{\partial r}+h_{j}.$ (37)

Below we focus on the first three orders of approximation, for which the above
system leads to the following results.

First, at $\mathcal{O}(\epsilon)$ (i.e., for $j=1$), the system (34)-(36) is
linear and homogeneous, namely, $f_{1}=g_{1}=h_{1}=0$; accordingly, Eq. (37)
(for $j=1$) becomes identical to the 1D version of Eq. (28), with
$\Delta=\partial_{r}^{2}$. Thus, a solution of Eq. (37) for $j=1$ is of the
form

$\rho_{1}\sim A\exp(i\theta)\,,\quad\theta=kr-\omega t,$ (38)

with $\omega$ and $k$ obeying the Jeans dispersion relation (29) and $A$ being
an unknown complex envelope function, which depends on the slow variables $R$
and $T_{2}$. Notice that for the above solution, we have considered the case
of right-going waves, with $\omega=\sqrt{k^{2}-1}$. In addition, it can
readily be seen that the remaining fields $u_{1}$ and $\Phi_{1}$ take a
similar form, and thus a solution of the system (34)-(36) reads

$\displaystyle\rho_{1}$ $\displaystyle=$
$\displaystyle-k^{2}A\exp(i\theta)+{\rm c.c.},$ (39) $\displaystyle u_{1}$
$\displaystyle=$ $\displaystyle-k\omega A\exp(i\theta)+{\rm c.c.},$ (40)
$\displaystyle\Phi_{1}$ $\displaystyle=$ $\displaystyle A\exp(i\theta)+{\rm
c.c.}.$ (41)

Here, the unknown complex function $A$ represents, at this order of
approximation, the envelope of the perturbation of the background, which
modulates the carrier wave $\exp(i\theta)$; as we will show below, the
function $A$ satisfies a spherical NLS equation.

Proceeding to the next order, $\mathcal{O}(\epsilon^{2})$ (i.e., $j=2$), the
system (34)-(36) becomes inhomogeneous, with the right-hand side terms given
by:

$\displaystyle f_{2}$ $\displaystyle=$
$\displaystyle-\frac{\partial\rho_{1}}{\partial T_{1}}-\frac{\partial
u_{1}}{\partial R_{1}}-\frac{\partial}{\partial r}\left(\rho_{1}u_{1}\right),$
(42) $\displaystyle g_{2}$ $\displaystyle=$ $\displaystyle-\frac{\partial
u_{1}}{\partial T_{1}}-\frac{\partial\rho_{1}}{\partial
R_{1}}-\frac{\partial\Phi_{1}}{\partial R_{1}}-\rho_{1}\frac{\partial
u_{1}}{\partial t}-u_{1}\frac{\partial u_{1}}{\partial
r}-\rho_{1}\frac{\partial\Phi_{1}}{\partial r},$ (43) $\displaystyle h_{2}$
$\displaystyle=$ $\displaystyle-2\frac{\partial^{2}\Phi_{1}}{\partial
r\partial R_{1}}.$ (44)

The solvability condition for the system (34)-(36) at this order (or
equivalently for Eq. (37) for $j=2$) is obtained by demanding that there be no
terms proportional to $\exp(i\theta)$ on the r.h.s. of (37), which in turn
guarantees the absence of secular terms in the solution for $\rho_{2}$. In
particular, the secular term is of the form 222Notice that if we had only used
the variables $R_{1}$ and $T_{1}$, the relevant term would have been of the
form
$\exp(i\theta)\left(\partial_{T_{1}}-\omega^{\prime}\partial_{R_{1}}\right)A$.
This secular term vanishes as long as $A$ depends on
$R=R_{1}-\omega^{\prime}T_{1}=R_{1}-VT_{1}$.
$(V-\omega^{\prime}(k))A_{R}\exp(i\theta)$, and its removal yields the
following expression for the unknown velocity $V$:

$V=\frac{k}{\omega}=\omega^{\prime}(k).$ (45)

This means that $V$ is actually equal to the group velocity
$\omega^{\prime}(k)\equiv\partial\omega/\partial k$, as derived from the Jeans
dispersion relation (29). In addition, at the same order, we obtain the
following expressions for the unknown fields $\rho_{2}$, $u_{2}$, $\Phi_{2}$:

$\displaystyle\rho_{2}$ $\displaystyle=$
$\displaystyle-\frac{2}{3}k^{4}(2k^{2}-3)A^{2}\exp(2i\theta)+{\rm c.c.},$ (46)
$\displaystyle u_{2}$ $\displaystyle=$
$\displaystyle-\frac{k^{3}(4k^{4}-7k^{2}+3)}{3\omega}A^{2}\exp(2i\theta)+\frac{i}{\omega}\frac{\partial
A}{\partial R}\exp(i\theta)+{\rm c.c.},$ (47) $\displaystyle\Phi_{2}$
$\displaystyle=$
$\displaystyle\frac{1}{6}k^{4}(2k^{2}-3)A^{2}\exp(2i\theta)+\frac{2i}{k}\frac{\partial
A}{\partial R}\exp(i\theta)+{\rm c.c.}.$ (48)

Finally, at order $\mathcal{O}(\epsilon^{3})$ (i.e., for $j=3$), the right-
hand side terms of the system (34)-(36) are given by:

$\displaystyle f_{3}=$ $\displaystyle-$
$\displaystyle\frac{2}{VT_{2}}u_{1}-\frac{\partial\rho_{1}}{\partial
T_{2}}-\frac{\partial\rho_{2}}{\partial T_{1}}-\frac{\partial}{\partial
R_{1}}(\rho_{1}u_{1})-\frac{\partial}{\partial
r}(\rho_{1}u_{2}+\rho_{2}u_{1})-\frac{\partial u_{2}}{\partial R_{1}},$ (49)
$\displaystyle g_{3}=$ $\displaystyle-$ $\displaystyle\frac{\partial
u_{1}}{\partial T_{2}}-\rho_{1}\frac{\partial u_{1}}{\partial
T_{1}}-\frac{\partial u_{2}}{\partial T_{1}}-u_{1}\frac{\partial
u_{1}}{\partial R_{1}}-\frac{\partial\rho_{2}}{\partial
R_{1}}-\rho_{1}\frac{\partial\Phi_{1}}{\partial
R_{1}}-\frac{\partial\Phi_{2}}{\partial R_{1}}$ (50) $\displaystyle-$
$\displaystyle\rho_{1}\frac{\partial u_{2}}{\partial
t}-\rho_{1}u_{1}\frac{\partial u_{1}}{\partial r}-\frac{\partial}{\partial
r}(u_{1}u_{2})-\rho_{1}\frac{\partial\Phi_{2}}{\partial
r}-\rho_{2}\frac{\partial\Phi_{1}}{\partial r},$ $\displaystyle h_{3}=$
$\displaystyle-$
$\displaystyle\frac{2}{VT_{2}}\frac{\partial\Phi_{1}}{\partial
r}-\frac{\partial^{2}\Phi_{1}}{\partial
R_{1}^{2}}-2\frac{\partial^{2}\Phi_{2}}{\partial r\partial R_{1}}.$ (51)

At this order too, the solvability condition for Eqs. (34)-(36) (or
equivalently for Eq. (37)), is also obtained by removing the secular terms,
$\propto\exp(i\theta$). This way, we derive the following NLS equation for the
envelope function $A$:

$i\frac{\partial A}{\partial
T_{2}}+\frac{i}{T_{2}}A+\frac{1}{2}\omega^{\prime\prime}(k)\frac{\partial^{2}A}{\partial
R^{2}}+\gamma(k)|A|^{2}A=0,$ (52)

where the dispersion and nonlinearity coefficients, $\omega^{\prime\prime}(k)$
and $\gamma(k)$, are given by:

$\displaystyle\omega^{\prime\prime}(k)$ $\displaystyle\equiv$
$\displaystyle\frac{\partial^{2}\omega}{\partial
k^{2}}=-\frac{1}{\omega^{3}},$ (53) $\displaystyle\gamma(k)$ $\displaystyle=$
$\displaystyle\frac{k^{4}(8k^{4}-21k^{2}+12)}{6\omega}.$ (54)

The functional forms of $\omega^{\prime\prime}(k)$ and $\gamma(k)$ are
depicted in Fig. 2, in the wavenumber regime $k>k_{J}$ (i.e., in the regime of
interest, free of the Jeans instability). It is observed that the dispersion
coefficient is $\omega^{\prime\prime}(k)<0\leavevmode\nobreak\ \forall k$,
while the nonlinearity coefficient becomes zero for $k=\alpha k_{J}$, where
the value of parameter $\alpha$ is given by $\alpha=(1/4)\sqrt{21+\sqrt{57}}$.
For $k_{J}<k<\alpha k_{J}$ we have $\gamma(k)<0$, while for $k>\alpha k_{J}$
we have $\gamma(k)>0$.

Figure 2: The dependence of the dispersion [solid (blue) line] and
nonlinearity [dashed (red) line] coefficients, $\omega^{\prime\prime}$ and
$\gamma$, on the wavenumber $k$, in the regime $k>k_{J}$. The dispersion
coefficient is always negative, while the nonlinearity coefficient becomes
zero at $k=\alpha k_{J}$, with $\alpha=(1/4)\sqrt{21+\sqrt{57}}$, and it is
negative (positive) for $k<\alpha k_{J}$ ($k>\alpha k_{J}$).

##  4 The nonlinear instability

We proceed with our analysis upon identifying an instability which stems from
the NLS Eq. (52). We will show that the latter possesses a spatially
homogeneous and time-dependent solution, which is prone to modulational
instability. This instability that occurs in the nonlinear regime of the
system also depends on the curvature and, in some cases, may be arrested. This
picture complements the one obtained in the linear regime of the problem,
where solely the Jeans instability is present.

We start by introducing the transformation:

$A(R,T_{2})=\frac{T_{0}}{T_{2}}\leavevmode\nobreak\ \Psi(R,T_{2}),$ (55)

where $T_{0}=\epsilon^{2}t_{0}$ is an initial time (see below). Then, by
substituting in the NLS Eq. (52), it is straightforward to find that the
function $\Psi$ satisfies the following equation:

$i\frac{\partial\Psi}{\partial
T_{2}}+\frac{1}{2}\omega^{\prime\prime}(k)\frac{\partial^{2}\Psi}{\partial
R^{2}}+\gamma(k)\frac{T_{0}^{2}}{T_{2}^{2}}\leavevmode\nobreak\
|\Psi|^{2}\Psi=0.$ (56)

Observe that Eq. (56) is actually a traditional $1$D NLS equation, but with
the important difference that the nonlinearity coefficient $\gamma(k)$ is now
replaced by $\gamma(k)(T_{0}/T_{2})^{2}$. Hence, in other words, the spherical
NLS equation (52) can be viewed as a traditional NLS, which gets linearized,
becoming a linear Schrödinger equation, in the asymptotic limit of
$T_{2}\rightarrow+\infty$.

Next, we proceed by identifying an elementary, spatially homogeneous and time-
dependent solution of the NLS Eq. (56), which is of the form:

$\Psi_{b}(T_{2})=\Psi_{0}\exp[i\Theta(T_{2})],$ (57)

where $\Psi_{0}$ is the constant amplitude of the solution, assumed to be
real, and $\Theta(T_{2})$ is a phase, with the associated nonlinear frequency
shift, $d\Theta/dT_{2}$, satisfying the equation:

$\frac{d\Theta}{dT_{2}}=\gamma(k)\Psi_{0}^{2}\frac{T_{0}^{2}}{T_{2}^{2}}.$
(58)

Next, we wish to study the linear stability of the solution $\Psi_{b}(T_{2})$.
For this purpose, we introduce the linearization ansatz:

$\Psi(R,T_{2})=[\Psi_{0}+(\delta\Psi)(R,T_{2})]\exp[i\Theta(T_{2})],$ (59)

where $|(\delta\Psi)|\ll\Psi_{0}$ is a small, generally complex, perturbation.
Substituting into Eq. (56), linearizing with respect to $\delta\Psi$, and also
using Eq. (58), we obtain:

$i\frac{\partial(\delta\Psi)}{\partial
T_{2}}+\frac{1}{2}\omega^{\prime\prime}(k)\frac{\partial^{2}(\delta\Psi)}{\partial
R^{2}}+\gamma(k)\Psi_{0}^{2}\frac{T_{0}^{2}}{T_{2}^{2}}[(\delta\Psi)+(\bar{\delta\Psi})]=0,$
(60)

where $\bar{(\delta\Psi)}$ is the complex conjugate of $(\delta\Psi)$.
Separating real and imaginary parts in $(\delta\Psi)$, namely
$(\delta\Psi)=(\delta\Psi)_{r}+i(\delta\Psi)_{i}$, Eq. (60) leads to the
following system of linear PDEs:

$\displaystyle-\frac{\partial(\delta\Psi)_{i}}{\partial
T_{2}}+\frac{1}{2}\omega^{\prime\prime}(k)\frac{\partial^{2}(\delta\Psi)_{r}}{\partial
R^{2}}+2\gamma(k)\Psi_{0}^{2}\frac{T_{0}^{2}}{T_{2}^{2}}(\delta\Psi)_{r}=0,$
(61) $\displaystyle\frac{\partial(\delta\Psi)_{r}}{\partial
T_{2}}+\frac{1}{2}\omega^{\prime\prime}(k)\frac{\partial^{2}(\delta\Psi)_{i}}{\partial
R^{2}}=0.$ (62)

Then, the compatibility condition
$\partial^{3}(\delta\Psi)_{i})/(\partial^{2}R\partial
T_{2})=\partial^{3}(\delta\Psi)_{i})/(\partial T_{2}\partial^{2}R)$ leads to
the following equation for $(\delta\Psi)_{r}$:

$\frac{\partial^{2}(\delta\Psi)_{r}}{\partial
T_{2}^{2}}+\frac{1}{4}(\omega^{\prime\prime}(k))^{2}\frac{\partial^{4}(\delta\Psi)_{r}}{\partial
R^{4}}+\gamma(k)\omega^{\prime\prime}(k)\Psi_{0}^{2}\frac{T_{0}^{2}}{T_{2}^{2}}\frac{\partial^{2}(\delta\Psi)_{r}}{\partial
R^{2}}=0.$ (63)

To this end, considering solutions of the form,

$(\delta\Psi)_{r}(R,T_{2})=\psi_{0}(T_{2})\exp(-iKR)+{\rm c.c.},$ (64)

we cast Eq. (63) into the form:

$\frac{\partial^{2}\psi_{0}}{\partial T_{2}^{2}}+\Omega^{2}(T_{2})\psi_{0}=0,$
(65)

where

$\Omega^{2}(T_{2})=\frac{1}{4}(\omega^{\prime\prime}(k))^{2}K^{2}\left(K^{2}-K_{c}^{2}\right),\quad
K_{c}^{2}(T_{2})=\frac{4\gamma(k)A_{0}^{2}}{\omega^{\prime\prime}(k)}\leavevmode\nobreak\
\frac{T_{0}^{2}}{T_{2}^{2}},$ (66)

with $A_{0}=\Psi_{0}$ (see Eq. (55) for $T_{2}=T_{0}$). Equation (65) is a
2nd-order ordinary differential equation (ODE) with a time-dependent
coefficient, $\Omega^{2}(T_{2})$. Obviously, if $\Omega^{2}(T_{2})<0$ or, in
other words, if

$\gamma(k)\omega^{\prime\prime}(k)>0,\quad{\rm and}\quad|K|<K_{c},$ (67)

then the elementary solution (57) is prone to modulational instability (MI).
Since $\omega^{\prime\prime}<0$ for every $k$ [see Eq. (53) and Fig. 2], the
first of the above MI conditions is satisfied when $\gamma(k)<0$. As discussed
above, the latter condition is fulfilled in the interval [see Eq. (54)]:

$k_{J}<|k|<\alpha k_{J},\quad\alpha=\frac{1}{4}\sqrt{21+\sqrt{57}}\approx
1.34,$ (68)

which defines the nonlinear instability band; this is highlighted by a red
line-segment in Fig. 1. Hence, our analysis reveals that the Jeans instability
band, $|k|<k_{J}$, is extended by a factor $\alpha$ (i.e., up to $\approx
33.6\%$), when the dynamics of the full nonlinear problem is taken into
account.

The dependence of the critical wavenumber $K_{c}$ on time renders the
instability growth rate, Im$(\Omega)$, also a function of time (this is a
typical situation occurring in the context of the spherical and cylindrical
NLS equations, appearing in studies on plasma physics —see, e.g., Refs. 14,
15, 28). In particular, the growth rate Im$(\Omega)$, as can be found from Eq.
(66), is of the form:

${\rm
Im}(\Omega)(K,T_{2})=\frac{1}{2}|\omega^{\prime\prime}(k)|K\sqrt{K_{c}^{2}(T_{2})-K^{2}},$
(69)

and is depicted in Fig. 3 as a function of $|K|<K_{c}$ for three different
values of the slow time $T_{2}$, namely for $T_{2}=0.005$, $0.008$ and $0.01$,
corresponding to the later stages of the evolution, e.g., for $t=0.5$, $0.8$,
$1$ and $\epsilon=0.1$. (The carrier wavenumber has been fixed to
$k=(1+\alpha)/2\approx 1.17$, i.e., in the middle of the nonlinear instability
band.) It is important to point out that, as shown in Fig. 3, the maximum
growth rate of the nonlinear instability is of the same order of magnitude as
the linear Jeans instability (see Fig. 1). Furthermore, it is observed that,
due to the dependence of $K_{c}$ on the slow time $T_{2}$, both the maximum
growth rate and the range of the nonlinear instability band decrease with
$T_{2}$.

Figure 3: The instability growth rate, ${\rm Im}(\Omega)$ as a function of $K$
[see Eq. (69)], for $k=(1+\alpha)k_{J}/2\approx 1.17$ [i.e., in the middle of
the nonlinearity-induced extension of the instability band —see Eq. (68)], and
for three different times: $T_{2}=5\times 10^{-3}$ [upper, dashed-dotted
(blue) curve], $T_{2}=8\times 10^{-3}$ [middle, dashed (red) curve], and
$T_{2}=10^{-2}$ [lower, solid (orange) curve]. Notice that both the maximum
value of ${\rm Im}(\Omega)$ and the range of the nonlinear instability band
decrease with $T_{2}$.

An important remark at this point is the following. Contrary to the 1D case
[when the curvature term $iA/T_{2}$ in Eq. (52) is absent], the nonlinear
instability ceases to exist for $K_{c}^{2}-K^{2}=0$, i.e., above some critical
time $T_{2}^{\rm(c)}$ (a similar result was obtained in relevant studies in
the context of plasma physics —see, e.g., Refs. 14, 15, 28). In particular,
the instability growth will cease for $T_{2}\geq T_{2}^{\rm(c)}$, where, for a
given wavenumber $K$, the critical time is given by:

$\displaystyle T_{2}^{\rm(c)}$ $\displaystyle=$
$\displaystyle\frac{T_{0}A_{0}}{K}\sqrt{\frac{4|\gamma(k)|}{|\omega^{\prime\prime}(k)|}}.$
(70)

Hence, the growth $\Gamma$ of the perturbation that takes place as long as the
nonlinear instability occurs, i.e., from the initial time $T_{0}$ up to the
time $T_{2}^{\rm(c)}$, is given by:

$\Gamma(K)=\exp\left(\int_{T_{0}}^{T_{2}^{\rm(c)}}{\rm
Im}(\Omega)dT_{2}^{\prime}\right)=\exp\left(A_{0}^{2}T_{0}|\gamma(k)|F(w)\right),$
(71)

where the function $F(w)$ is defined as follows:

$\displaystyle
F(w)=\frac{2}{w}\Bigg{[}\sqrt{w}\ln\left(\sqrt{w}+\sqrt{w-1}\right)-\sqrt{w-1}\Bigg{]},\quad\text{with}\quad
w\equiv\dfrac{K_{c}^{2}(T_{0})}{K^{2}}=\dfrac{4\gamma(k)A_{0}^{2}}{\omega^{\prime\prime}(k)K^{2}}.$
(72)

The function $F(w)$ is shown in Fig. 4. As seen therein, $F(w)$ initially
increases, then reaches a maximum, namely

${\rm max}F(w_{c})=\frac{2\sqrt{w_{c}-1}}{w_{c}},$ (73)

where $w_{c}$ is determined by the transcendental equation

$2\sqrt{w_{c}-1}=\sqrt{w_{c}}\ln\left(\sqrt{w_{c}}+\sqrt{w_{c}-1}\right),$
(74)

while $F(w)$ decreases slowly thereafter, with
$\lim_{w\rightarrow\infty}F(w)=0$. Note that ${\rm max}F(w_{c})\approx 0.55$
and $w_{c}\approx 12$ (see Fig. 4). Thus, according to the above discussion,
while the growth of the MI is known to be always an increasing function of
time in the 1D problem 17, this is not the case in the spherical geometry
under consideration; hence, in the spherical geometry, wave patterns are
expected to be, in principle, structurally more stable than the 1D ones. This
can also be explained by the formal reduction of the spherical NLS [via the
transformation of Eq. (55)] to the 1D NLS (56), with a time-dependent
nonlinearity coefficient, $\propto 1/T_{2}^{2}$. This suggests that, at
sufficiently large times, the NLS becomes a linear Schrödinger equation as
mentioned above, which in turn explains the fact that (a) the MI instability
threshold is time- dependent [see Eq. (69)], and (b) wave structures of the
spherical model are more stable than ones of the $1$D model.

Let us now estimate the typical length scale of the instability as the
wavenumber of the maximally growing mode. Since the only $K$-dependence of the
growth factor comes via $F(w)$, which peaks at $w_{c}=12$, we have that
$K_{c}^{2}(T_{0})/K_{\rm MI}^{2}=12$ and substituting $K_{c}(T_{0})$ from
(66), we obtain (recall that $k=\epsilon K$)

$k_{\rm
MI}(k)=\epsilon\sqrt{\frac{\gamma(k)}{3\omega^{\prime\prime}(k)}}A_{0}.$ (75)

For a carrier wave with $k=1.17$ the dispersion/nonlinearity coefficients read
$|\omega^{\prime\prime}|=4.4$ and $|\gamma|=0.9$. Suppose that the background
(57) becomes order one at some time $T_{2}=T_{0}A_{0}$. Then for $T_{2}$ close
to the end of the dynamics, say $T_{2}=10^{-2}$, and $T_{0}=10^{-3}$, with
$A_{0}\sim 50$, we obtain a total growth of an order of magnitude, while the
MI scale reads $k_{\rm MI}\approx 0.1$. Upon defining the characteristic
length as $\lambda_{\rm MI}=2\pi/k_{\rm MI}$, we see that structures may
emerge at scales larger than the Jeans length. This is compatible with the
small-amplitude approximation that we have used, which holds true at the tails
of galactic halos where the density is close to its equilibrium value. In
relation to this let us also note that, while the function $F(w)$ clearly
attains a maximum at $w_{c}\approx 12$, it decreases quite slowly thereafter.
This means that the instability persists for scales much larger than the Jeans
length.

Figure 4: The function $F(w)$ given in Eq. (72). It is observed that this
function is initially increasing, then features a maximum ${\rm
max}F(w_{c})\approx 0.55$ at $w=w_{c}\approx 12$, and then it asymptotically
decreases to zero.

##  5 Solitons and rogue waves

The derivation of the effective NLS Eq. (52) that was presented above allows
us to predict the existence of coherent wave structures, namely soliton
solutions, which may describe large scale structures in the universe. Exact
soliton solutions exist only in the 1D setting, i.e., in the absence of the
curvature term $iA/T_{2}$. Nevertheless, as we will discuss below, approximate
soliton solutions, analogous to the 1D case, are possible in the spherical
geometry as well.

Generally, soliton structures are supported by the NLS in the case where
dispersion and nonlinearity are counterbalanced, i.e., when
$|\omega^{\prime\prime}(k)/\gamma(k)|={\cal O}(1)$. Given the form of
$\omega^{\prime\prime}(k)$ and $\gamma(k)$, this requirement is fulfilled
around the middle of the interval $k_{J}<k<\alpha k_{J}$ where the nonlinear
instability occurs, as seen in Fig. 5. Furthermore, the emergence of solitons
is naturally connected with the initial conditions associated with the NLS
equation. In that regard, it is necessary to make the following relevant
comments.

Figure 5: The ratio $|\omega^{\prime\prime}(k)/\gamma(k)|$ as a function of
the wavenumber $k$ in the interval $k_{J}<k<\alpha k_{J}$, i.e., in the
nonlinear instability band. As seen, this ratio is ${\cal O}(1)$ for
wavenumbers around the middle of this interval.

First, it is mentioned that coherent patterns can result from the onset of the
nonlinear instability, i.e., the modulational instability (MI) mechanism
discussed in the previous Section. Although MI was shown to occur for small-
amplitude perturbations on top of the spatially uniform and time-dependent
elementary solution (57), any initial condition —also in the form of a random
perturbation (which is the physically relevant condition in our problem)— can
be viewed as a Fourier superposition of plane waves. Such a superposition will
naturally contain wavenumbers that undergo MI 333Here it should be stressed
that we are referring to initial conditions characterized by wavenumbers which
are Jeans stable; such a situation may occur in physically relevant parameter
regimes, as explained in the end of Section 4. which, in turn, during its
evolution, may give rise to localized coherent structures. Indeed, it can be
shown (see, e.g., Ref.29) that at the early stage of its evolution, the MI
results in the excitation of integer multiples of the perturbation wavenumbers
of the unstable wave solutions, i.e., causes a widening of the wavenumber
spectrum, which results in localization in space. This localization mechanism
suggests the possibility of soliton formation.

The most celebrated solutions of the NLS equation, which were found by means
of the Inverse Scattering Transform method18, are the so-called “bright
solitons” (this term is coined from the field of nonlinear optics30). Such
bell-shaped envelope solitons play a key role in the long-time dynamics of the
NLS model, since they act as statistical attractors to which the system
relaxes31, 32 (see also numerical evidence in Ref.33). Nevertheless, another
important waveform that may be generated due to MI, when seeded from random
perturbations of a background, is a single or a series of high-contrast peaks
of random intensity; these peaks can be identified as “extreme events” or
“rogue waves” resulting from the unstable growth of weak wave modulations,
which evolve into short groups of steep waves (see, e.g., the reviews19, 20
and references therein). These waves can be mathematically described by the
so-called Peregrine soliton34. In contrast to the “regular” sech-shaped NLS
soliton solution discussed above, the rogue wave is represented by a rational
function, with the property of being localized in both time and space. Both of
these important solutions will be discussed below.

We now proceed by presenting soliton solutions of the effective NLS Eq. (52).
First, in the purely 1D setting, the NLS equation is a completely integrable
system, and possesses exact soliton solutions18. The NLS bright soliton is
characterized by a pulse-shaped envelope wave structure. In our problem, it is
relevant to consider the stationary form of the NLS soliton, namely:

$\displaystyle A_{\rm sol}(R,T_{2})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{|\omega^{\prime\prime}|}{|\gamma|}}\leavevmode\nobreak\
\eta\leavevmode\nobreak\ {\rm
sech}\left[\eta\left(R-R_{0}\right)\right]\exp\left(i\eta^{2}|\omega^{\prime\prime}|T_{2}\right),$
(76)

where $\eta$ is a free ${\cal O}(1)$ parameter, characterizing the amplitude,
the inverse width and the frequency of the soliton, while $R_{0}$ represents
the location of the soliton center.

In addition, in the 1D setting, there exists another exact soliton solution,
which is localized both in space and time and decays algebraically rather than
exponentially, contrary to the soliton of Eq. (76). The pertinent solution is
the so-called Peregrine soliton34 mentioned above, which is a rational
solution of Eq. (52), of the following form:

$\displaystyle A_{\rm RW}(R,T_{2})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{|\omega^{\prime\prime}|}{|\gamma|}}\eta\left[1-\frac{4(1+2i\eta^{2}|\omega^{\prime\prime}|T_{2})}{1+4\eta^{2}(R-R_{0})^{2}+4\eta^{4}|\omega^{\prime\prime}|^{2}T_{2}^{2}}\right]\exp(i\eta^{2}|\omega^{\prime\prime}|T_{2}),$
(77)

where the free parameter $\eta$ is connected, as before, with the amplitude,
inverse width and frequency of the rogue wave, while $R_{0}$ represents again
the location of the center. Observe that the soliton waveform is on top of a
finite background, decays to it asymptotically for either large $R$ or $T_{2}$
but exhibits a nontrivial behavior over a small region of $(R,T_{2})$. That is
why these localized solutions, that “appear from nowhere and disappear without
a trace”35 are considered as possible prototypes of rogue waves.

Here, a couple of comments are in order. First, although the Peregrine soliton
does not decay to zero as $R\rightarrow\infty$ (as was assumed in our
multiscale perturbation scheme), the soliton core may well describe extreme
wave events that may occur on top of the equilibrium density of the system
(see relevant results and pertinent discussions in Refs.36, 37, 38). Second,
in our problem, the time needed for the Peregrine soliton to disappear is very
large (recall that $T_{2}=\epsilon^{2}t$). Lastly, the solution (77) suggests
that fluctuations around the equilibrium density decay algebraically [rather
than exponentially, as is the case corresponding to the soliton (76)].

Next, we proceed with the spherical geometry where, as mentioned above, exact
analytical soliton solutions of Eq. (52) are not available. Nevertheless,
approximate soliton solutions can be found in specific limiting cases,
depending on the order of magnitude of the terms involved in the NLS equation
(52). In particular, substituting either the soliton (76) or the rogue wave
(77) into the NLS and considering a specific time interval $\tau>0$ (i.e.,
$T_{0}\leq T_{2}\leq T_{0}+\tau$), we obtain:

$\displaystyle\frac{1}{T_{2}}A\propto\frac{\eta}{\tau},\quad\frac{1}{2}|\omega^{\prime\prime}(k)|A_{RR}\approx|\gamma(k)||A|^{2}A\approx\eta^{3},$
(78)

in the aforementioned regime $|\omega^{\prime\prime}(k)/\gamma(k)|={\cal
O}(1)$, where 1D solitons do exist. Then, two limiting cases can readily be
identified:

  * (a)

The curvature effect is much stronger than the effects of dispersion and
nonlinearity, namely,

$\frac{1}{T_{2}}A\gg\frac{1}{2}|\omega^{\prime\prime}(k)|A_{RR}\approx|\gamma(k)||A|^{2}A,$

for which case we have $\eta^{2}\tau\ll 1$. In other words, in a situation
where the initial amplitude $\eta$ is sufficiently small and the time interval
$\tau$ is sufficiently short, the dispersive and nonlinear terms of Eq. (52)
can be neglected. Then, in this limit, the NLS is reduced to the linear
equation:

$A_{T_{2}}+\frac{1}{T_{2}}A=0,$ (79)

which, upon integration, leads to the following expression for $A(R,T_{2})$:

$A(R,T_{2})=\frac{T_{0}}{T_{2}}F(R),$ (80)

where $F(R)$ is an arbitrary function of $R$. Then, we may construct an
approximate soliton solution upon assuming that $F(R)$ has initially the
profile of a planar soliton or a rogue wave [i.e., $F(R)$ has the functional
form given by Eq. (76) or Eq. (77) for $T_{2}=T_{0}$]. In such a case, Eq.
(80) dictates that the relevant waveform decays in time, according to the law
$T_{0}/T_{2}$. For instance, the approximate functional form of the planar
soliton reads:

$\displaystyle A(R,T_{2})$ $\displaystyle\approx$
$\displaystyle\eta\frac{T_{0}}{T_{2}}\sqrt{\frac{|\omega^{\prime\prime}|}{|\gamma|}}\leavevmode\nobreak\
{\rm
sech}\left[\eta\left(R-R_{0}\right)\right]\exp\left(i\eta^{2}|\omega^{\prime\prime}|T_{0}\right),$
(81)

while the form of the rogue wave can be found in a similar manner.

  * (b)

The curvature effect is much weaker than the effects of dispersion and
nonlinearity, namely,

$\frac{1}{T_{2}}A\ll\frac{1}{2}|\omega^{\prime\prime}(k)|A_{RR}\approx|\gamma(k)||A|^{2}A,$

for which case we have $\eta^{2}\tau\gg 1$. Assuming now that $\eta={\cal
O}(1)$, the latter inequality suggests that for sufficiently large times
($\tau\gg 1$), the effects of dispersion and nonlinearity dominate over the
curvature effect. This means that the latter can be considered as a small
perturbation, whose effect can be studied in the framework of the so-called
adiabatic approximation of the perturbation theory for solitons39. According
to this approach, which can be applied in the case of the planar soliton (76),
the functional form of the soliton remains unchanged but its amplitude becomes
time-dependent and evolves according to the law
$\eta=\eta_{0}(T_{0}/T_{2})^{2}$; here, $\eta_{0}$ is a constant setting the
solitary wave amplitude at the initial time $T_{2}=T_{0}$. Thus, an
approximate spherical soliton solution, valid in the regime of long times
($\tau\gg 1$), takes the form:

$\displaystyle
A(R,T_{2})\approx\eta_{0}\left(\frac{T_{0}}{T_{2}}\right)^{2}\sqrt{\frac{|\omega^{\prime\prime}|}{|\gamma|}}\leavevmode\nobreak\
{\rm
sech}\left[\eta_{0}\left(\frac{T_{0}}{T_{2}}\right)^{2}\left(R-R_{0}\right)\right]\exp\left\\{i|\omega^{\prime\prime}|\int_{T_{0}}^{T_{2}}\eta_{0}^{2}\left(\frac{T_{0}}{T_{2}}\right)^{4}dT_{2}^{\prime}\right\\}.$
(82)

To conclude this section, it is found that, in both aforementioned asymptotic
regimes (a) and (b), the solitons undergo a slow (algebraic) decay when the
curvature effect is taken into account. It is also pointed out that the
analytical estimates for the bright soliton decay presented above have also
been discussed in the context of plasma physics; there, direct numerical
simulations performed in the framework of Eq. (52), were found to be in good
agreement with the analytical results14, 15. Similar findings were also
reported in Ref.40 for the case of rogue waves. Both solutions describe
spherical shells of radius $R_{0}$ which decay exponentially/algebraically as
$R\rightarrow\infty$; in the cosmological context, such structures may appear
at the rims of galactic halos, where the density $\rho$ is close to the
background density $\rho_{0}$, as per our considerations in the multiscale
expansion method.

##  6 Conclusions

In this work we have studied the nonlinear dynamics of a self-gravitating
fluid model in spherical geometry, which may find applications in cosmic
structure formation.

The 1D version of the model without gravity, being a purely nonlinear
hyperbolic system, gives rise to shock waves, i.e., the development of
discontinuities in finite time. When gravity is included, the model features
—in addition to nonlinearity— dispersion. In particular, small-amplitude
(linear) excitations $\propto\exp[i(kr-\omega t)]$ (where $\omega$ and $k$
represent the frequency and the wavenmber) on top of the steady-state
background are characterized by a dispersion relation of the form
$\omega^{2}(k)=k^{2}-k_{J}^{2}$, where $k_{J}$ is the so-called Jeans
wavenumber. Obviously, since $\partial^{2}\omega(k)/\partial k^{2}\neq 0$, the
system features dispersion and nonlinearity simultaneously, a fact suggesting
the possibility of soliton formation.

In the linear regime described above, the particular form of the dispersion
relation indicates the presence of an important instability, namely, the Jeans
instability, which manifests itself for wavenumbers $|k|<k_{J}$. The Jeans
instability occurs when gravity dominates the internal pressure, resulting in
gravitational collapse and the formation of large-scale structures in the
universe.

Based on the above, we have tried to address two important questions: (a)
apart from the Jeans instability that is predicted in the linear regime of the
model, is there any other instability that may occur in the full nonlinear
version of the model, and (b) since we are dealing with a nonlinear dispersive
system, can we predict and derive soliton solutions? Our analysis suggests a
positive answer to both questions.

In particular, we have employed a multiscale expansion method and derived from
the original spherically symmetric self-gravitating fluid model an effective
spherical nonlinear Schrödinger (NLS) equation for the envelope of small-
amplitude perturbations around the equilibrium state of the system. This NLS
model allowed us to predict an instability occurring in the nonlinear regime
of the system, namely the modulational instability (MI) of the weak, spatially
uniform, time-varying solution of the NLS (this solution is the simplest
nontrivial perturbation satisfying the NLS on top of the equilibrium state of
the fluid model). We have found that, while the Jeans instability occurs for
$|k|<k_{J}$, the MI sets in for $k_{J}<|k|<\alpha k_{J}$, with $\alpha\approx
1.34$. This means that Jeans stable wavenumbers may become unstable in the
aforementioned MI band. Since at the early stage of the evolution the MI gives
rise to widening of the wavenuber spectrum —which, in turn, causes
localization in real space— soliton formation is possible even outside the
Jeans instability domain. In other words, in the physical problem under
consideration, nonlinearity (and, in particular, the “nonlinear instability”
that manifests itself in the nonlinear regime), when competing with the
gravity-induced dispersion, may give rise to large scale structures in the
universe.

Indeed, the derived effective NLS equation allowed us to predict such
structures, which may be of the form of solitons or rogue waves. These
structures feature, respectively, an exponential or an algebraic decay in
space, while both of them experience a slow curvature-induced polynomial decay
in time. The functional form of the predicted soliton states could be relevant
for the description of the edges of dark-matter halos.

Our analysis and results suggest a number of interesting research directions
for future studies. First of all, it would be relevant to perform direct
numerical simulations in the spherical self-gravitating fluid model to better
clarify the various possibilities investigated in this work. In addition, it
would also be interesting to consider the effect of viscosity in the self-
gravitating fluid. Furthermore, a study of a self-gravitating two-fluid model,
with each fluid representing the luminous and the dark matter, is particularly
challenging. The above themes are currently under investigation and relevant
results will be presented elsewhere.

## Acknowledgments

This work is dedicated to Professor Ioannis G. Stratis, a close friend and
collaborator, with whom we had the chance to have numerous discussions about
science, research and academic life during his great career. Constructive
discussions with T. A. Apostolatos are kindly acknowledged. The research work
of G.N.K. was supported by the Hellenic Foundation for Research and Innovation
(HFRI) under the HFRI PhD Fellowship grant (Fellowship Number: $5860$). O.E.
and S.S. have been supported by Thailand NSRF via PMU-B (grant numbers
B01F650006 and B37G660013).

## References

  * 1 Courant R, Friedrichs KO. Supersonic Flow and Shock Waves.  Springer, New York .  1948\.
  * 2 Whitham GB. Linear and Nonlinear Waves.  Wiley, New York .  1974\.
  * 3 Riemann B. Ueber die Fortpflanzung ebener Luftwellen von endlicher Schwingungsweite. Abhandlungen der Königlichen Gesellschaft der Wissenschaften zu Göttingen 1860; Retrieved 2012-08-08: 1-23.
  * 4 Lax PD. Development of Singularities of Solutions of Nonlinear Hyperbolic Partial Differential Equations. J. Math. Phys. 1964; 5: 611-613.
  * 5 Kates RE. Self-focusing as a mechanism for the formation of dense structures in Newtonian cosmology. Astronomy and Astrophysics 1986; 168: 1-5.
  * 6 Ono H, Nakata I. Soliton Formation in a Self-Gravitating Gas. Progr. Theor. Phys. 1994; 92: 9-15.
  * 7 Jeans JH. The stability of a spherical nebula. Phil. Trans. Roy. Soc. (London) 1902; 199A: 1-53.
  * 8 Chavanis P. Dynamical stability of infinite homogeneous self-gravitating systems and plasmas: application of the Nyquist method. Eur. Phys. J. B 2012; 85: 229-259.
  * 9 Hui L, Ostriker JP, Tremaine S, Witten E. Ultralight scalars as cosmological dark matter. Phys. Rev. D 2017; 95: 043541.
  * 10 Böhmer CG, Harko T. Can dark matter be a Bose-Einstein condensate?. J. Cosmology and Astroparticle Phys. 2007; 06: 025.
  * 11 Chavanis P. Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions. I. Analytical results. Phys. Rev. D 2011; 84: 043531.
  * 12 Mocz P, Vogelsberger M, Robles VH, et al. Galaxy formation with BECDM - I. Turbulence and relaxation of idealized haloes. Monthly Notices Roy. Astr. Soc. 2017; 471: 4559-4570.
  * 13 Schive HY, Chiueh T, Broadhurst T. Cosmic structure as the quantum interference of a coherent dark wave. Nat. Phys. 2014; 10: 496-499.
  * 14 Sabry R, Moslem WM, Shukla PK, Saleem H. Cylindrical and sperical ion-acoustic envelope solitons in multicomponent plasmas with positrons. Phys. Rev. E 2009; 79: 056402.
  * 15 Xue JK. Propagation of nonplanar dust-acoustic envelope solitary waves in a two-ion-temperature dusty plasma. Phys. Plasmas 2004; 11: 1860-1865.
  * 16 Huang G, Lou S, Dai X. Cylindrical envelope solitary waves. Chin. Phys. Lett. 1990; 7: 398-401.
  * 17 Zakharov V, Ostrovsky L. Modulation instability: The beginning. Physica D: Nonlinear Phenomena 2009; 238: 540-548.
  * 18 Zakharov VE, Shabat AB. Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media. Sov. Phys. JETP 1971; 34: 62-69.
  * 19 Onorato M, Residoric S, Bortolozzoc U, Montinad A, Arecchie FT. Rogue waves and their generating mechanisms in different physical contexts. Phys. Reports 2013; 528: 47-89.
  * 20 Tikan A, Randoux S, El G, Tovbis A, Copie F, Suret P. Local emergence of Peregrine solitons: experiments and theory. Frontiers in Phys. 2021; 8: 599435.
  * 21 Rubin VC. A brief history of dark matter: 1–13; Space Telescope Science Institute Symposium Series. Cambridge University Press .  2004
  * 22 Blok dWJG. The Core-Cusp Problem. Advances in Astronomy 2010; 2010: 1–14.  doi: 10.1155/2010/789293
  * 23 Schive HY, Chiueh T, Broadhurst T. Cosmic Structure as the Quantum Interference of a Coherent Dark Wave. Nature Phys. 2014; 10: 496–499.  doi: 10.1038/nphys2996
  * 24 Armendariz-Picon C, Neelakanta JT. How Cold is Cold Dark Matter?. JCAP 2014; 03: 049.  doi: 10.1088/1475-7516/2014/03/049
  * 25 Hu W, Barkana R, Gruzinov A. Cold and fuzzy dark matter. Phys. Rev. Lett. 2000; 85: 1158-1161.  doi: 10.1103/PhysRevLett.85.1158
  * 26 Muller CM. Cosmological bounds on the equation of state of dark matter. Phys. Rev. D 2005; 71: 047302.  doi: 10.1103/PhysRevD.71.047302
  * 27 Jeffrey A, Kawahara T. Asymptotic Methods in Nonlinear Wave Theory.  Pitman, Boston .  1982\.
  * 28 Demiray H. Modulation of cylindrical (spherical) waves in a plasma with vortex electron distribution. Phys. Plasmas 2018; 25: 072113.
  * 29 Hasegawa A, Kodama Y. Solitons in Optical Communications.  Clarendon Press, Oxford .  1995\.
  * 30 Kivshar YS, Agrawal GP. Optical solitons: from fibers to photonic crystals.  Academic Press, San Diego .  2003\.
  * 31 Zakharov VE, Pushkarev AN, Shvets VF, V.Yankov V. Soliton turbulence. Pisma Zh. Eksp. Teor. Fiz. 1988; 48: 79-82.
  * 32 Dyachenko S, Zakharov VE, Pushkarev AN, Shvets VF, Yankov VV. Soliton turbulence in nonintegrable wave systems. Zh. Eksp. Teor. Fiz. 1989; 96: 2026-2031.
  * 33 Jordan R, Josserand C. Self-organization in nonlinear wave turbulence. Phys. Rev. E 2000; 61: 1527-1539.
  * 34 Peregrine DH. Water waves, nonlinear Schrödinger equations and their solutions. J. Austral. Math. Soc. Ser. B 1983; 25: 16-43.
  * 35 Akhmediev N, Ankiewicz A, Taki M. Waves that appear from nowhere and disappear without a trace. Phys. Lett. A 2009; 373: 675-678.
  * 36 Anastassi Z, Fotopoulos G, Frantzeskakis DJ, et al. Spatiotemporal algebraically localized waveforms for a nonlinear Schrödinger model with gain and loss. Physica D 2017; 355: 24-33.
  * 37 Horikis TP, Ablowitz MJ. Rogue waves in nonlocal media. Phys. Rev. E 2017; 95: 042211.
  * 38 Ablowitz MJ, Horikis TP. Rogue waves in birefringent optical fibers: elliptical and isotropic fibers. J. Opt. 2017; 19: 065501.
  * 39 Kivshar YS, Malomed BA. Dynamics of solitons in nearly integrable systems. Rev. Mod. Phys. 1989; 61: 763-915.
  * 40 Sabry R. Freak waves in Saturn’s magnetosphere. Astrophys. Space Sci. 2015; 355: 33-41.