#  Chains of Boson Stars

C. A. R. Herdeiro†, J. Kunz‡, I. Perapechka⋆, E. Radu† and Y. Shnir⋄

†Department of Mathematics, University of Aveiro and CIDMA, Campus de
Santiago, 3810-183 Aveiro, Portugal
‡ Institute of Physics, University of Oldenburg, Germany Oldenburg D-26111,
Germany
⋆ Department of Theoretical Physics and Astrophysics, Belarusian State
University, Minsk 220004, Belarus
⋄ BLTP, JINR, Dubna 141980, Moscow Region, Russia

(January 2021)

###### Abstract

We study axially symmetric multi-soliton solutions of a complex scalar field
theory with a sextic potential, minimally coupled to Einstein’s gravity. These
solutions carry no angular momentum and can be classified by the number of
nodes of the scalar field, $k_{z}$, along the symmetry axis; they are
interpreted as chains with $k_{z}+1$ boson stars, bound by gravity, but kept
apart by repulsive scalar interactions. Chains with an odd number of
constituents show a spiraling behavior for their ADM mass (and Noether charge)
in terms of their angular frequency, similarly to a single fundamental boson
star, as long as the gravitational coupling is small; for larger coupling,
however, the inner part of the spiral is replaced by a merging with the
fundamental branch of radially excited spherical boson stars. Chains with an
even number of constituents exhibit a truncated spiral pattern, with only two
or three branches, ending at a limiting solution with finite values of ADM
mass and Noether charge.

##  1 Introduction

Many non-linear physical systems support non-topological solitons, which
represent spatially localized field configurations. One of the simplest
examples in flat space is given by $Q$-balls, which are particle-like
configurations in a model with a complex scalar field possessing a harmonic
time dependence and a suitable self-interaction potential [1, 2, 3]. When
$Q$-balls are coupled to gravity, the so-called Boson Stars (BSs) emerge,
which represent solitonic solutions with a topologically trivial and globally
regular geometry. The simplest such configurations are static and spherically
symmetric, the scalar field possessing a mass term only, without self-
interaction [4, 5]. These solutions are usually dubbed mini-Boson Stars (mBS),
being regarded as macroscopic quantum states, which are prevented from
gravitationally collapsing by Heisenberg’s uncertainty principle; also, they
do not have a regular flat spacetime limit.

Both $Q$-balls and BSs carry a Noether charge associated with an unbroken
continuous global $U(1)$ symmetry. The charge $Q$ is proportional to the
angular frequency of the complex boson field and represents the boson particle
number of the configurations [2, 3].

In flat spacetime $Q$-balls exist only within a certain frequency range for
the scalar field: between a maximal value $\omega_{\rm max}$, which
corresponds to the mass of the scalar excitations, and some lower non-zero
critical value $\omega_{\rm min}$, that depends on the form of the potential.
On the one hand, as the frequency $\omega$ is approaching its extremal values,
both the mass $M$ and the Noether charge $Q$ of the configurations diverge. On
the other hand, $M$ and $Q$ attain a minimum at some critical value
$\omega_{\rm cr}\in[\omega_{\rm min},\omega_{\rm max}]$. Away from
$\omega_{\rm cr}$ the mass and the charge of the configurations increase
towards the divergent value at the boundary of the domain. Within
$[\omega_{\rm min},\omega_{\rm cr}]$ the configurations become unstable
against decay.

The situation is different for BSs: gravity modifies the critical behavior
pattern of the configurations. The fundamental branch of the BS solutions
starts off from the perturbative excitations at $\omega\sim\omega_{\rm max}$,
at which $M,Q$ trivialise (rather than diverge). Then, the BSs exhibit a
spiraling (or oscillating) pattern of the frequency dependence of both charge
and mass, where both tend to some finite limiting values at the centers of the
spirals. Qualitatively, the appearance of the frequency-mass spiral may be
related to oscillations in the force-balance between the repulsive scalar
interaction and the gravitational attraction in equilibria. Indeed, radially
excited rotating BSs do not exhibit a spiraling behavior; instead the second
branch extends all the way back to the upper critical value of the frequency
$\omega_{\rm max}$, forming a loop [6].

The main purpose of this paper is to report on the existence of a new type of
solutions, which correspond to chains of BSs.111Similar chains, but for a
scalar field without self-interactions, were recently studied in [7], in the
context of multipolar BSs. Here, we emphasise the interpretation of multiple
BSs, rather than a single multipolar BS. These are static, axially symmetric
equilibrium configurations interpreted as a number of BSs located
symmetrically with respect to the origin, along the symmetry axis. We
construct these solutions and investigate their physical properties for a
choice of the scalar field potential with quartic and sextic self-interaction
terms, which was employed in most of the $Q$-balls literature. We note that
similar configurations of chains of constituents are known to exist both for
gravitating and flat space non-Abelian monopoles and dyons [8, 9, 10, 11, 12,
14, 15, 13, 16, 17, 18, 19, 20], Skyrmions [21, 22, 23], electroweak
sphalerons [24, 25, 26, 27, 28], $SU(2)$ non-self dual configurations [29, 30]
and Yang-Mills solitons in ADS4 spacetime [31, 32].

In these multi-component BSs configurations the constituents form a chain
along the symmetry axis and, consequently, the scalar field is required to
possess ’nodes’ (zeros of the scalar field amplitude). Configurations with
$k_{z}$ nodes on the symmetry axis possess $k_{z}+1$ constituents.
Configurations with even and odd numbers of constituents can show a different
pattern, when their domain of existence is mapped out. In particular, we find
that the pattern exhibited by a mass-frequency diagram of the chains of BSs
can differ both from the typical spiraling picture and from the closed loop
scenario. For chains with an even number of constituents the pattern always
terminates at critical solutions. For chains with an odd number of
constituents, the pattern depends on the strength of the gravitational
interaction. The configurations then either merge with the corresponding
radial excitation of the fundamental solution, or the central constituent of
the configurations exhibits oscillations while retaining smaller satellite
constituents.

This paper is organized as follows. In Section II the theoretical setting is
specified. This includes the action, the equations of motion, the Ansatz and
the boundary conditions for the BS chains. The numerical results for these new
equilibrium configurations are shown in Section III . We give our conclusions
in Section IV.

##  2 The model

###  2.1 Action, field equations and global charges

We consider a self-interacting complex scalar field $\Phi$, which is minimally
coupled to Einstein’s gravity in a $(3+1)$-dimensional space-time. The
corresponding action of the system is

$\mathcal{S}=\int d^{4}x\sqrt{-g}\left[\frac{R}{16\pi
G}-\frac{1}{2}g^{\mu\nu}\left(\Phi_{,\,\mu}^{*}\Phi_{,\,\nu}+\Phi_{,\,\nu}^{*}\Phi_{,\,\mu}\right)-U(|\Phi|^{2})\right],$
(2.1)

where $R$ is the Ricci scalar curvature, $G$ is Newton’s constant, the
asterisk denotes complex conjugation, and $U$ denotes the scalar field
potential.

Variation of the action (2.1) with respect to the metric leads to the Einstein
equations

$E_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R-8\pi G~{}T_{\mu\nu}=0\ ,$
(2.2)

where

$T_{\mu\nu}\equiv\Phi_{,\mu}^{*}\Phi_{,\nu}+\Phi_{,\nu}^{*}\Phi_{,\mu}-g_{\mu\nu}\left[\frac{1}{2}g^{\sigma\tau}(\Phi_{,\sigma}^{*}\Phi_{,\tau}+\Phi_{,\tau}^{*}\Phi_{,\sigma})+U(|\Phi|^{2})\right]\,,$
(2.3)

is the stress-energy tensor of the scalar field.

The corresponding equation of motion of the scalar field is the non-linear
Klein-Gordon equation

$\left(\Box-\frac{dU}{d|\Phi|^{2}}\right)\Phi=0\,,$ (2.4)

where $\Box$ represents the covariant d’Alembert operator.

The solutions considered in this work have a static line-element (with a
timelike Killing vector field $\xi$), being topologically trivial and globally
regular, $i.e.$ without an event horizon or conical singularities, while the
scalar field is finite and smooth everywhere. Also, they approach
asymptotically the Minkowski spacetime background. Their mass $M$ can be
obtained from the respective Komar expressions [33],

${M}=\frac{1}{{4\pi G}}\int_{\Sigma}R_{\mu\nu}n^{\mu}\xi^{\nu}dV~{}.$ (2.5)

Here $\Sigma$ denotes a spacelike hypersurface (with the volume element $dV$),
while $n^{\mu}$ is normal to $\Sigma$, with $n_{\mu}n^{\mu}=-1$ [33]. After
replacing in (2.5) the Ricci tensor by the stress-energy tensor, via the
Einstein equations (2.2), one finds the equivalent expression

$\displaystyle
M=\,2\int_{\Sigma}\left(T_{\mu\nu}-\frac{1}{2}\,g_{\mu\nu}\,T_{\gamma}^{\
\gamma}\right)n^{\mu}\xi^{\nu}dV\ .$ (2.6)

The action (2.1) is invariant with respect to the global $\mathrm{U}(1)$
transformations of the complex scalar field, $\Phi\to\Phi e^{i\alpha}$, where
$\alpha$ is a constant. The following Noether 4-current is associated with
this symmetry

$j_{\mu}=-i(\Phi\partial_{\mu}\Phi^{\ast}-\Phi^{\ast}\partial_{\mu}\Phi)\,.$
(2.7)

It follows that integrating the timelike component of this 4-current in a
spacelike slice $\Sigma$ yields a second conserved quantity – the Noether
charge:

$\displaystyle Q=\int_{\Sigma}j^{\mu}n_{\mu}dV.$ (2.8)

###  2.2 The ansatz and equations

Apart from being static, the configurations in this work are also axially
symmetric222The scalar field possesses a time dependence (with
$\partial_{t}\Phi=-i\omega\Phi$), which disappears at the level of the energy-
momentum tensor. However, the scalar field is axially symmetric,
$\partial_{\varphi}\Phi=0$. . Thus, in a system of adapted coordinates, the
space-time possesses two commuting Killing vector fields

$\displaystyle\xi=\partial/\partial t~{}~{}{\rm
and}~{}~{}\eta=\partial/\partial\varphi,$ (2.9)

with $t$ and $\varphi$ the time and azimuthal coordinates, respectively. Line
elements with these symmetries are usually studied by employing a metric
ansatz [34]

$\displaystyle
ds^{2}=-e^{-2U(\rho,z)}dt^{2}+e^{2U(\rho,z)}\left[e^{2k(\rho,z)}(d\rho^{2}+dz^{2})+P(\rho,z)^{2}d\varphi^{2}\right],$
(2.10)

where $(\rho,z)$ correspond, asymptotically, to the usual cylindrical
coordinates333 In the Einstein-Maxwell theory, it is always possible to set
$P\equiv\rho$ , such that only two independent metric functions appear in the
equations, with $(\rho,z)$ the canonical Weyl coordinates [34]. For a
(complex) scalar field matter content, however, the generic metric ansatz
(2.10) with three independent functions is needed. . However, it the numerical
treatment of the Einstein-Klein-Gordon equations, it is useful to employ
‘quasi-isotropic’ spherical coordinates $(r,\theta)$, defined by the
coordinate transformation in (2.10)

$\displaystyle\rho=r\sin\theta,~{}~{}z=r\cos\theta~{},$ (2.11)

with the usual range $0\leqslant r<\infty$, $0\leqslant\theta\leq\pi$. Then
the metric can then be written in a Lewis-Papapetrou form, with

$ds^{2}=-fdt^{2}+\frac{m}{f}\left(dr^{2}+r^{2}d\theta^{2}\right)+\frac{l}{f}r^{2}\sin^{2}\theta
d\varphi^{2}~{}.$ (2.12)

The three metric functions $f$, $l$, and $m$ are functions of the variables
$r$ and $\theta$ only, chosen such that the trivial angular and radial
dependence of the (asymptotically flat) line element is already factorized.
The relation between the metric functions in the above line-element and those
in (2.10) is $f=e^{-2U}$, $lr^{2}\sin^{2}\theta=P^{2}$, $m=e^{2k}$. The
symmetry axis of the spacetime is located at the set of points with vanishing
norm of $\eta$, $||\eta||=0$; it corresponds to the $z$-axis in (2.10) or the
set with $\theta=0,\pi$ in (2.12). The Minkowski spacetime background is
approached for $r\to\infty$, with $f=l=m=1$.

For the scalar field we adopt an Ansatz with a harmonic time dependence, while
the amplitude depends on $(r,\theta)$,

$\Phi=\phi(r,\theta)e^{-i\omega t},$ (2.13)

where $\omega\geq 0$ is the angular frequency. The corresponding stress-energy
tensor is static, with the nonvanishing components

$\displaystyle
T_{r}^{r}=\frac{f}{m}\bigg{(}\phi_{,r}^{2}-\frac{\phi_{,\theta}^{2}}{r^{2}}\bigg{)}+\frac{\omega^{2}\phi^{2}}{f}-U(\phi^{2}),~{}~{}T_{\theta}^{\theta}=\frac{f}{m}\bigg{(}-\phi_{,r}^{2}+\frac{\phi_{,\theta}^{2}}{r^{2}}\bigg{)}+\frac{\omega^{2}\phi^{2}}{f}-U(\phi^{2}),~{}~{}T_{r}^{\theta}=\frac{2f}{r^{2}m}\phi_{,r}\phi_{,\theta}$
$\displaystyle
T_{\varphi}^{\varphi}=-\frac{f}{m}\bigg{(}\phi_{,r}^{2}+\frac{\phi_{,\theta}^{2}}{r^{2}}\bigg{)}+\frac{\omega^{2}\phi^{2}}{f}-U(\phi^{2}),~{}~{}T_{t}^{t}=-\frac{f}{m}\bigg{(}\phi_{,r}^{2}+\frac{\phi_{,\theta}^{2}}{r^{2}}\bigg{)}-\frac{\omega^{2}\phi^{2}}{f}-U(\phi^{2}).$
(2.14)

After inserting the ansatz (2.12), (2.13) into the field equations (2.2),
(2.4) we find a system of six coupled partial differential equations that
needs to be solved. There are three equations for the metric functions
$f,l,m$, found by taking the following suitable combinations of the Einstein
equations, $E_{r}^{r}+E_{\theta}^{\theta}=0$, $E_{\varphi}^{\varphi}=0$ and
$E_{t}^{t}=0$, yielding

$\displaystyle
f_{,rr}+\frac{f_{,\theta\theta}}{r^{2}}+\frac{2f_{,r}}{r}+\frac{\cot\theta
f_{,\theta}}{r^{2}}-\frac{1}{f}\bigg{(}f_{,r}^{2}+\frac{f_{,\theta}^{2}}{r^{2}}\bigg{)}+\frac{1}{2l}\bigg{(}f_{,r}l_{,r}+\frac{f_{,\theta}l_{,\theta}}{r^{2}}\bigg{)}+16\pi
G\left(U(\phi^{2})-\frac{2\omega^{2}\phi^{2}}{f}\right)m=0~{},$ $\displaystyle
l_{,rr}+\frac{l_{,\theta\theta}}{r^{2}}+\frac{3l_{,r}}{r}+\frac{2\cot\theta
l_{,\theta}}{r^{2}}-\frac{1}{2l}\bigg{(}l_{,r}^{2}+\frac{l_{,\theta}^{2}}{r^{2}}\bigg{)}+32\pi
G\left(U(\phi^{2})-\frac{\omega^{2}\phi^{2}}{f}\right)\frac{lm}{f}=0~{},$
(2.15) $\displaystyle
m_{,rr}+\frac{m_{,\theta\theta}}{r^{2}}+\frac{m_{,r}}{r}+\frac{m}{2f^{2}}\bigg{(}f_{,r}^{2}+\frac{f_{,\theta}^{2}}{r^{2}}\bigg{)}-\frac{1}{m}\bigg{(}m_{,r}^{2}+\frac{m_{,\theta}^{2}}{r^{2}}\bigg{)}$
$\displaystyle\ \ \ \ \ \ \ +16\pi
G\left[\frac{f}{m}\bigg{(}\phi_{,r}^{2}+\frac{\phi_{,\theta}^{2}}{r^{2}}\bigg{)}+U(\phi^{2})-\frac{\omega^{2}\phi^{2}}{f}\right]\frac{m^{2}}{f}=0~{},$

and one equation for the scalar field amplitude

$\displaystyle\phi_{,rr}+\frac{\phi_{,\theta\theta}}{r^{2}}+\left(\frac{2}{r}+\frac{l_{,r}}{2l}\right)\phi_{,r}+\left(\cot\theta+\frac{l_{,\theta}}{2l}\right)\frac{\phi_{,\theta}}{r^{2}}+(\frac{\omega^{2}}{f}-\frac{\partial
U}{\partial\phi^{2}})\frac{m}{f}\phi=0\ .$ (2.16)

Apart from these, there are two more Einstein equations,
$E_{\theta}^{r}=0,~{}E_{r}^{r}-E_{\theta}^{\theta}=0$, which are not solved in
practice, being treated as constraints and used to check the numerical
accuracy of the solutions.

The mass $M$ is computed from the Komar expression (2.5), where we insert the
metric ansatz (2.12), with unit vector $n=-\sqrt{f}dt$, the volume element
$dV=1/\sqrt{f}\,\sqrt{-g}\,dr\,d\theta\,d\varphi$, and
$\sqrt{-g}=r^{2}\sin\theta\frac{\sqrt{l}m}{f}$. In evaluating (2.5), we use
the fact that $R_{t}^{t}$ is a total derivative:

$\displaystyle\sqrt{-g}R_{t}^{t}$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial
r}\left(\frac{r^{2}\sin\theta\sqrt{l}f_{,r}}{2f}\right)-\frac{\partial}{\partial\theta}\left(\frac{\sin\theta\sqrt{l}f_{,\theta}}{2f}\right).$

Then it follows that $M$ can be read off from the asymptotic expansion of the
metric function $f$

$\displaystyle f=1-\frac{2MG}{r}+\mathcal{O}\left(\frac{1}{r^{2}}\right)\ .\ \
\ $ (2.17)

Alternatively, the mass $M$ can be obtained by direct integration of (2.6) ,

$\displaystyle
M=\int\left(T_{t}^{t}-T_{r}^{r}-T_{\theta}^{\theta}-T_{\varphi}^{\varphi}\right)\,\sqrt{-g}\,dr\,d\theta\,d\varphi=4\pi\int_{0}^{\infty}dr\int_{0}^{\pi}d\theta~{}r^{2}\sin\theta\frac{\sqrt{l}m}{f}\left(U(\phi^{2})-\frac{2\omega^{2}\phi^{2}}{f}\right)~{}.$
(2.18)

In terms of the above Ansatz the Noether charge $Q$ is given by

$Q=4\pi\omega\int\limits_{0}^{\infty}dr\int\limits_{0}^{\pi}d\theta~{}\frac{\sqrt{l}m}{f^{2}}r^{2}\sin\theta~{}\phi^{2}~{}.$
(2.19)

Also, let us note that the solutions in this work have no horizon. Therefore
they are zero entropy objects, without an intrinsic temperature. Still, in
analogy to black holes, one may write a “first law of thermodynamics” [35],
albeit without the entropy term, which reads

$dM=\omega dQ~{}.$ (2.20)

This gives a relation between the mass and Noether charge of neighbouring BS
solutions which can be used to check the numerical accuracy of the solutions.

###  2.3 The potential and scaling properties

The solutions in this work were found for a potential originally proposed in
[36, 37] ,

$U(|\Phi|^{2})=\nu|\Phi|^{6}-\lambda|\Phi|^{4}+\mu^{2}|\Phi|^{2}$ (2.21)

which is chosen such that nontopological soliton solutions - $Q$-balls - exist
in the absence of the Einstein term in the action (2.1), $i.e.$ on a Minkowski
spacetime background [38, 39, 40, 41]. Also, at least in the spherically
symmetric case, this choice of the potential follows for the existence of very
massive and highly compact objects approaching the black hole limit [38].

The parameter $\mu$ in (2.21) corresponds to the mass of the scalar
excitations around the $|\Phi|=0$ vacuum. Apart from that, the potential
possesses two more parameters $(\nu,\lambda)>0$, determining the self-
interactions, which are chosen in such a manner that it possesses a local
minimum at some finite non-zero value of the field $|\Phi|$, besides the
global minimum at $|\Phi|=0$.

Two of the constants in (2.21) can actually be absorbed into a redefinition of
the radial coordinate together with a rescaling of the scalar field,

$\displaystyle r\to\frac{u_{0}}{\mu}r\
,~{}~{}~{}\phi\to\frac{\sqrt{\mu}}{\nu^{1/4}\sqrt{u_{0}}}\phi\ ,$ (2.22)

with $u_{0}>0$ an arbitrary constant. Note that the scalar field frequency
changes accordingly, $\omega\to\frac{u_{0}}{\mu}\omega.$

Then, the potential (2.21) becomes (up to an overall factor),

$\displaystyle U=\phi^{6}-\bar{\lambda}\phi^{4}+u_{0}^{2}\phi^{2},\qquad{\rm
with}~{}~{}\bar{\lambda}=\frac{\lambda u_{0}}{\mu\sqrt{\nu}}\ .$ (2.23)

The choice employed in most of the $Q$-ball literature is

$\displaystyle u_{0}^{2}=1.1\ ,~{}~{}{\rm and}~{}~{}\bar{\lambda}=2\ ,$ (2.24)

which are the values used also in this work.

For completeness, let us mention that for the specific ansatz (2.12), (2.13)
with the above scalings, the equations for gravitating $Q$-balls (which are
effectively solved) can also be derived by extremizing the following reduced
action444In (2.26) we use the dimensionless radial variable $r$ and the
dimensionless frequency $\omega$ given in units set by $\mu$.

$\displaystyle S_{red}=\int drd\theta\left({\cal L}_{g}-4\alpha^{2}{\cal
L}_{s}\right),$ (2.25)

with [for ease of notation, we denote $(\nabla S)\cdot(\nabla T)\equiv
S_{,r}T_{,r}+\frac{1}{r^{2}}S_{,\theta}T_{,\theta}$]:

$\displaystyle{\cal L}_{g}$ $\displaystyle=$ $\displaystyle
r^{2}\sin\theta\sqrt{l}\bigg{[}\frac{1}{2lm}(\nabla l)\cdot(\nabla
m)-\frac{1}{2f^{2}}(\nabla f)^{2}$ (2.26)
$\displaystyle{~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}}-\frac{1}{rl}\left(l_{,r}+\frac{\cot\theta}{r}l_{,\theta}\right)+\frac{1}{rm}\left(m_{,r}+\frac{\cot\theta}{r}m_{,\theta}\right)\bigg{]},$
$\displaystyle{\cal L}_{s}$ $\displaystyle=$ $\displaystyle
r^{2}\sin\theta\frac{m\sqrt{l}}{f}\bigg{[}\frac{f}{m}(\nabla\phi)^{2}-\frac{\omega^{2}\phi^{2}}{f}+\phi^{6}-2\phi^{4}+1.1\phi^{2}\bigg{]}.$

The (dimensionless) coupling constant reads

$\displaystyle\alpha^{2}=\frac{4\pi G\mu}{\sqrt{\nu}u_{0}},$ (2.27)

which determines the strength of the gravitational coupling of the solutions.

As with the known spherically symmetric configurations, solutions exist for
$0\leqslant\alpha<\infty$. The limit $\alpha\to 0$ corresponds to the non-
backreacting case, $i.e.$ $Q$-balls on a fixed Minkowski background. To
understand the large $\alpha$ limit, we define $\hat{\phi}=\phi/\alpha$. Then
for large values of the effective gravitational coupling $\alpha$, the non-
linearity of the potential (2.21) becomes suppressed and the system approaches
the usual Einstein-Klein-Gordon model [4, 5], with its corresponding mBS
solutions. That is, the resulting equations are identical to those found for a
model with a non-self-interacting, massive complex scalar field $\hat{\phi}$.
However, in this work we shall restrict our study to the case of a finite,
nonzero $\alpha$. The basic properties of the mBS chains were studied (in a
more general context) in the recent work [7], while the issue of flat
spacetime $Q$-ball chains will be discussed elsewhere.

###  2.4 Boundary conditions

The solutions studied in this work are globally regular and asymptotically
flat, possessing finite mass and Noether charge. Appropriate boundary
conditions guarantee that these conditions are satisfied.

Starting with the boundary conditions for the metric functions, regularity of
the solutions at the origin requires

$\partial_{r}f\bigl{.}\bigr{|}_{r=0}=\partial_{r}m\bigl{.}\bigr{|}_{r=0}=\partial_{r}l\bigl{.}\bigr{|}_{r=0}=0\,.$
(2.28)

Demanding asymptotic flatness at spatial infinity yields

$f\bigl{.}\bigr{|}_{r\to\infty}=m\bigl{.}\bigr{|}_{r\to\infty}=l\bigl{.}\bigr{|}_{r\to\infty}=1.$
(2.29)

Axial symmetry and regularity impose on the symmetry axis at $\theta=0,\pi$
the conditions

$\partial_{\theta}\phi\bigl{.}\bigr{|}_{\theta=0,\pi}=\partial_{\theta}f\bigl{.}\bigr{|}_{\theta=0,\pi}=\partial_{\theta}m\bigl{.}\bigr{|}_{\theta=0,\pi}=\partial_{\theta}l\bigl{.}\bigr{|}_{\theta=0,\pi}=0\,.$
(2.30)

Further, the condition of the absence of a conical singularity requires that
the solutions should satisfy the constraint
$m\bigl{.}\bigr{|}_{\theta=0,\pi}=l\bigl{.}\bigr{|}_{\theta=0,\pi}$. In our
numerical scheme we explicitly verified (within the numerical accuracy) this
condition on the symmetry axis.

Turning now to the boundary conditions for the scalar field amplitude, we
mention first that $\phi$ approaches asymptotically the global minimum,

$\phi\bigl{.}\bigr{|}_{r\to\infty}=0~{},$ (2.31)

while on the symmetry axis we impose

$\partial_{\theta}\phi\bigl{.}\bigr{|}_{\theta=0,\pi}=0\ .$ (2.32)

The behaviour of the scalar field at the origin is more complicated, depending
on the considered parity ${\cal P}$. As mentioned in the Introduction, the
solutions split into two classes, distinguished by their behaviour $w.r.t.$ a
reflection along the equatorial plane $\theta=\pi/2$. The geometry is left
invariant under this transformation,

$\displaystyle
f(r,\pi-\theta)=f(r,\theta),~{}~{}l(r,\pi-\theta)=l(r,\theta),~{}~{}m(r,\pi-\theta)=m(r,\theta)~{};$
(2.33)

for the scalar field, however, there are two possibilities:

$\displaystyle{\cal
P}=1~{}~{}~{}{\rm(even~{}parity)}:~{}~{}\phi(r,\pi-\theta)={\phantom{-}}\phi(r,\theta)\
,$ (2.34) $\displaystyle{\cal
P}=-1~{}{\rm(odd~{}parity)}:~{}~{}\phi(r,\pi-\theta)=-\phi(r,\theta)~{}.$
(2.35)

We use this symmetry to reduce the domain of integration to
$[0,\infty)\times[0,\pi/2]$, with

$\partial_{\theta}f\bigl{.}\bigr{|}_{\theta=\pi/2}=\partial_{\theta}m\bigl{.}\bigr{|}_{\theta=\pi/2}=\partial_{\theta}l\bigl{.}\bigr{|}_{\theta=\pi/2}=0\,,$
(2.36)

and

${\cal
P}=1:~{}\partial_{\theta}\phi\bigl{.}\bigr{|}_{\theta=\pi/2}=0,~{}~{}~{}{\cal
P}=-1:~{}\phi\bigl{.}\bigr{|}_{\theta=\pi/2}=0,~{}~{}$ (2.37)

together with

$\displaystyle{\cal P}=1:~{}\partial_{r}\phi\bigl{.}\bigr{|}_{r=0}=0~{}\
,~{}~{}~{}{\cal P}=-1:~{}\phi\bigl{.}\bigr{|}_{r=0}=0\,.$ (2.38)

Finally, let us mention that, for both ${\cal P}=\pm 1$, one can formally
construct an approximate expression of the solutions compatible with the above
boundary conditions, $e.g.$ by assuming the existence of a power series form
close to $r=0$. The resulting expressions are of little help in understanding
the properties of the solutions, and a numerical approach is necessary555We
recall that, more than fifty years after their discovery [4, 5], the static
mBS are still not known in closed form. . However, one important result of
this study is the bound-state condition

$\displaystyle\omega\leqslant\mu,$ (2.39)

which emerges from the finite mass requirement. No similar result is found for
the minimal value of frequency.

###  2.5 Numerical method

The set of four coupled non-linear elliptic partial differential equations for
the functions $(f,l,m;\phi)$ has been solved numerically subject to the
boundary conditions defined above. In a first stage, a new compactified radial
variable $x$ is introduced, replacing $r$, with

$\displaystyle x\equiv\frac{r}{c+r},$ (2.40)

with $c$ an arbitrary parameter, typically of order one. With this choice, the
semi-infinite region $[0,\infty)$ is mapped to the finite interval $[0,1]$.
This avoids the use of a cutoff radius.

The numerical calculations have been performed by using a using a sixth-order
finite difference scheme. The system of equations is discretized on a grid
with a typical number of points $129\times 89$. The underlying linear system
is solved with the Intel MKL PARDISO sparse direct solver [50] using the
Newton-Raphson method. Calculations are performed with the packages
FIDISOL/CADSOL [51] and CESDSOL666Complex Equations – Simple Domain partial
differential equations SOLver is a C++ package being developed by one of us
(I.P.).  library. In all cases, the typical errors are of order of $10^{-4}$.

For the choice (2.24) of the potential’s parameters, the only input parameters
are the gravitational coupling constant $\alpha$ together with the scalar
field’s frequency $\omega$. The parity of the solutions is imposed via the
boundary conditions (2.34), (2.38). The number of individual constituents
results from the numerical output. Also, in all plots below, quantities are
given in natural units set by $\mu,G$.

##  3 The solutions

###  3.1 Nodal structure and energy distribution

The choice of the parity ${\cal P}$ is related to the number of distinct
constituents of the solutions (as resulting $e.g.$ from the spatial
distribution of the energy density). Moreover, denoting the number of nodes on
the symmetry axis by $k_{z}$, the solutions can be classified by the number of
nodes $k_{z}$. The number of constituents of the chains is given by $k_{z}+1$.
The even-parity configurations (${\cal P}=1$) have an even $k_{z}$, while the
solutions with an odd $k_{z}$ are found for ${\cal P}=-1$. For example, the
spherically symmetric solutions have $k_{z}=0$ and one single constituent
localized at the origin, $r=0$. The simplest non-spherical configuration has
$k_{z}=1$ and represents a pair of static BSs777 In principle, the $k_{z}=1$
solutions can be thought of as the static limit of negative parity spinning
configurations considered in [41].  with opposite phases, the inversion of the
sign of the scalar field function $\phi$ under reflections
$\theta\to\pi-\theta$ corresponds to the shift of the phase $\omega t\to\omega
t+\pi$.

It was pointed out that the character of the interaction between Q-balls in
Minkowski spacetime depends on their relative phase [42, 43]. If the Q-balls
are in phase, the interaction is attractive, if they are out of phase, there
is a repulsive force between them. Thus, an axially symmetric $k_{z}=1$
solution can be balanced by the gravitational attraction.

Solutions with $k_{z}>1$ exist as well; the maximal value we have reached so
far is $k_{z}=5$; however, they are likely to exist for an arbitrarily large
$k_{z}$.

To illustrate these aspects, we display in Fig. 1 several functions of
interest for the five types of representative configurations. Both odd and
even parity configurations are shown there, with the node numbers $k_{z}=0-4$;
also the solutions have the same values of the input parameters, $\alpha=0.25$
and $\omega/\mu=0.8$, being located on the 1st branch in the
$(\omega,M)$-diagram (as described below). These chains possess one to five
constituents (from top to bottom), as seen by the number of peaks of the
charge density, shown in the left panels. The middle panels represent the
scalar field amplitude $\phi$, and the right panels show the metric function
$f=-g_{tt}$. For the sake of clarity we have chosen to exhibit these figures
in polar coordinates $(\rho,z)$, as given by eq. (2.11).

Figure 1: Chains of BSs with one to five constituents (from top to bottom) on
the first (a.k.a. fundamental) branch for $\alpha=0.25$ at $\omega/\mu=0.80$:
$3d$ plots of the $U(1)$ scalar charge distributions (left plots), the scalar
field functions $\phi$ (middle plots) and the metric functions $f$ (right
plots) versus the coordinates $\rho=r\sin\theta$ and $z=r\cos\theta$.

The first row shows a single spherically symmetric BS for comparison. The
second row exhibits the pair of BSs. The charge density has two symmetric
peaks, the metric function has two symmetric troughs, while the scalar field
function is anti-symmetric, featuring a peak and a trough. The triplet,
quartet and quintet in the next few rows feature $k_{z}+1$ very similar (in
size and shape) peaks for the charge density and troughs for the metric
function, while the scalar field shows alternating peaks and troughs, all
located symmetrically along the $z$-axis. Thus on the fundamental branch we
basically encounter a chain consisting of $k_{z}+1$ BSs, all possessing
similar size, shape and distance from their next neighbors.

This picture partially changes as we move along the domain of solutions, for a
given coupling $\alpha$. This is seen in Fig. 2, where these chains
($k_{z}=0-4$) are now shown for illustrative solutions sitting on the second
branch of the $M$ $vs.$ $\omega$ domain of existence (see the discussion in
the next subsection), with $\alpha=0.25$ for $k_{z}=0-3$, $\alpha=0.5$ for
$k_{z}=4$, and $\omega/\mu=\\{0.43,0.47,0.57,0.7,0.7\\}$ for $k_{z}=0-4$,
respectively. As we look at the charge density of the configurations, we see a
dominant peak at the center for odd chains and a dominant inner pair for even
chains, while the other peaks of the triplet, quartet and quintet have turned
into slight elevations, hardly visible in the figures. In the metric functions
the troughs at the center of the odd chains dominate, and likewise the (almost
merged) inner pair of troughs of the even chains. All other troughs are
weakened substantially. The scalar field itself, however, retains the outer
peaks and troughs to a somewhat greater extent, still reflecting clearly the
number of constituents of the chains.

Figure 2: Chains of BSs with one to five constituents (from top to bottom) on
the second branch for differing values of $\alpha$ and $\omega$: $3d$ plots of
the $U(1)$ scalar charge distributions (left plots), the scalar field
functions $\phi$ (middle plots) and the metric functions $f$ (right plots)
versus the coordinates $\rho=r\sin\theta$ and $z=r\cos\theta$.

###  3.2 The $\omega$-dependence and the branch structure

Recall the frequency dependence for the single BSs and for fixed coupling
$\alpha$ [38]. The set of BSs emerges from the vacuum at the maximal
frequency, given by the boson mass $\omega_{\text{max}}=\mu$. Thus, unlike the
case of $Q$-balls in flat space, where mass and charge diverge, these
quantitites vanish in this limit. Decreasing the frequency $\omega$ spans the
first or fundamental branch, which terminates at the first backbending of the
curve, at which point it moves towards larger frequencies. The curve then
follows a spiraling/oscillating pattern, with successive backbendings.

The mass and charge form a spiral, as $\omega$ is varied, while the minimum of
the metric function $f$ and the maximum of the scalar function $\phi$ shows
damped oscillations. The set of solutions tends to a limiting solution at the
center of the spiral which has finite values of the mass and charge. However,
the values of the scalar field function $\phi$ and the metric function $f$ at
the center of the star, which represent the maximal and minimal values of
these functions, $\phi_{\text{max}}$ and $f_{\text{min}}$, respectively, do
not seem to be finite, with $\phi_{\text{max}}$ diverging and $f_{\text{min}}$
vanishing in the limit.

Let us now consider the frequency dependence of the BS chains, when the
coupling $\alpha$ is kept fixed. Like the single BSs, all chains emerge
similarly from the vacuum at the maximal frequency, given by the boson mass
$\omega_{\text{max}}=\mu$, where their mass and charge vanish in the limit.

When $\omega$ is decreased, mass and charge rise and the chains follow along
their fundamental branch. As for the single BS, for all chains this
fundamental branch ends at a minimal value of the frequency, from where a
second branch arises. But then even and odd chains will in general exhibit
different patterns, which will also depend on the coupling strength $\alpha$.

####  3.2.1 Even chains

Figure 3: Comparison of the $k_{z}=0$ single BSs (dashed curves), the
$k_{z}=1$ pair of BSs (solid curves), and the $k_{z}=3$ quartet of BSs (dash-
dotted curves): scaled ADM mass $M$ (upper left panel), scaled charge $Q$
(upper right panel), minimal value of the metric function $f_{\text{min}}$
(middle left panel), maximal value of the scalar field $\phi_{\text{max}}$ of
the $k_{z}=1$ pair only (middle right panel), and separation $z_{\text{d}}$
between the two components of the $k_{z}=1$ pair only (lower panel) $vs.$
frequency $\omega$ for several values of the coupling $\alpha$. Note the
quadratic scale for $f_{\text{min}}$ and $z_{\text{d}}$.

Let us consider the BS pair and the higher even chains in more detail. We
illustrate the $\omega$-dependence for even chains in Fig. 3, selecting the BS
pair ($k_{z}=1$) and the BS quartet ($k_{z}=3$), and comparing with the single
BS. In the upper panels we show the $k_{z}=1$ pair (solid curves), the
$k_{z}=3$ quartet (dash-dotted curves), as well as the $k_{z}=0$ single BS
(dashed curves). For the latter we only show the first few branches. In the
two upper panels we exhibit the scaled ADM mass $M$ (left) and the scaled
charge $Q$ (right). The different colors refer to different values of the
coupling $\alpha$. The middle left panel shows in an analogous manner the
minimal value $f_{\text{min}}$ of the metric function $f$. Restricting to the
$k_{z}=1$ pair, we then exhibit the maximal value $\phi_{\rm max}$ of the
scalar field function $\phi$ on right middle panel, and the separation
$z_{\text{d}}$ between the components of the $k_{z}=1$ pair in the lower
panel. All quantities are shown versus the frequency $\omega$.

Whereas the single BSs constitute an infinite set of branches, that form a
spiral or a damped oscillation, depending on the quantity of interest [38],
the even chains seem to end in a limiting configuration quite abruptly
somewhere in the middle of a branch. The number of branches before this
limiting configuration is encountered depends on the strength of the
gravitational coupling $\alpha$. For small $\alpha$ there are more branches,
for larger $\alpha$, the limiting configuration is already encountered on the
second branch, as illustrated in Fig. 3 for the BS pair.

The mass and the charge exhibit only the onset of a spiraling behavior with
two branches for the larger $\alpha$ values ($\alpha=0.5,1$), and three
branches for the smaller ones ($\alpha=0.15,0.2,0.25$). While the minimum
$f_{\text{min}}$ the metric function (middle left) and the maximum
$\phi_{\text{max}}$ of the scalar function (middle right) exhibit only two or
three oscillations. The coordinate distance $z_{\text{d}}$ between the two
components of the pair, as given by twice the value of the $z$-coordinate of
the maximum $\phi_{\text{max}}$, exhibits both types of behavior (lower). For
small $\alpha$ ($\alpha=0.15$) $z_{\text{d}}$ shows three oscillations, while
it decreases continuously. For larger $\alpha$ ($\alpha=0.25$) it exhibits the
onset of a spiral, with again larger values on the third branch.

Let us now address the limiting behavior of the pairs in more detail, and its
dependence on the coupling $\alpha$. To that end, we illustrate in Fig. 4 the
profiles of the metric function $f$ (left panels) and the scalar field
function $\phi$ (right panels) on the $z$-axis, choosing a large (upper
panels) and a small (lower panels) value of the coupling $\alpha$. Since for
the large $\alpha$, the coordinate distance $z_{\text{d}}$ does not decrease
monotonically, but increases again towards the limiting solution, we retain
two well-separated components in the limit. The minimum of the metric function
$f_{\text{min}}$ tends to zero in the limit, while the maximum (minimum) of
the scalar field function $\phi_{\text{max}}$ ($\phi_{\text{min}}$) becomes
extremely sharp.

Figure 4: Profiles on the symmetry axis of the metric function $f$ (left) and
the scalar field function $\phi$ of the $k_{z}=1$ (almost) critical solution
on the second branch at $\omega=0.85$ and $\alpha=1$ (upper plots), and on the
third branch at $\omega=0.52$ and $\alpha=0.5$ (lower plots). Also,
$z=r\cos\theta$, with $\theta=0,\pi$.

In fact, the scalar field amplitude $\phi$ acquires the shape of two
antisymmetric peak(on)s associated with the locations of the minima of the
metric function $f$, where the name peakons refers to field configurations
with extremely large absolute values of the second derivative at the maxima
[52, 53]. Further numerical investigation of such singular solutions is,
unfortunately, plagued by severe numerical problems. We remark, that in the
case of odd chains the choice of spherical coordinates and the presence of the
major peak(on) at the origin alleviate the numerical problems considerably.

For smaller $\alpha$, the coordinate distance $z_{\text{d}}$ between the two
components of the pair decreases monotonically towards the limiting solution.
Fig. 4 shows that the closeness of the extrema of the scalar field function
coincides with a very steep rise of the scalar field function at the origin.
The metric function, however, retains only a small peak at the origin. Here
the numerical grid allows for better resolution of this extremal behavior, and
thus a closer approach to the limiting solution. Still, the approach towards
the limiting solution is restricted by numerical accuracy.

The $k_{z}=3$ quartet represents a bound state of four BSs, located
symmetrically along the symmetry axis. It may also be viewed as a bound state
of two bound pairs of BSs, since this configuration is not fully symmetric.
The two inner BSs are slightly larger than the two outer BSs, though there is
not too much distinction between the stars as long as the configuration
resides on the fundamental branch. As seen in Fig. 3, these quartet
configurations share most of the properties with $k_{z}=1$ pairs. For the
chosen values of the coupling $\alpha$ we find two branches of solutions, the
fundamental branch connected to the perturbative excitations at
$\omega\to\omega_{\text{max}}$, and the second branch leading to a limiting
solution. While approaching this limiting solution the outer extrema of the
metric and of the scalar field function become less pronouced, leaving
basically the two inner extrema, which evolve completely analogously to the
extrema of the pair towards a limiting solution. We expect this scenario to
represent a generic pattern seen for all chains with odd $k_{z}$ (although so
far we have checked it for $k_{z}=1,3,5$ only).

####  3.2.2 Odd chains

As noticed above, the odd chains always possess a BS centered at the origin,
with the remaining BSs are located symmetrically with respect to the origin,
on the symmetry axis. The presence of the central BS constituent suggests that
the $(\omega,M)$-pattern of the odd chains could resemble that found for a
single BS. To scrutinize this conjecture, let us consider the branch structure
for the $k_{z}=2$ triplet of BSs, exhibited in Fig. 5 for four values of the
coupling $\alpha$ (with different colors). Analogously to the even chains, we
show in the two upper panels the scaled ADM mass $M$ (left) and the scaled
charge $Q$ (right). The middle panels show the minimal value $f_{\text{min}}$
of the metric function (left) and the maximal value $\phi_{\rm max}$ of the
scalar field function (right), while the lower panel shows the separation
$z_{\text{d}}$ between the neighboring components of the $k_{z}=2$ triplet.
Again, all quantities are shown as a function of the frequency $\omega$.

Figure 5: $k_{z}=2$ triplet of BSs: scaled ADM mass $M$ (upper left panel),
scaled charge $Q$ (upper right panel), minimal value of the metric function
$f_{\text{min}}$ (middle left panel), maximal value of the scalar field
function $\phi_{\text{max}}$ (middle right panel), and separation
$z_{\text{d}}$ between neighboring components of the $k_{z}=2$ triplet (lower
panel) $vs.$ frequency $\omega$ for several values of the coupling $\alpha$.
Note the quadratic scale for $Q$, $f_{\text{min}}$ and $z_{\text{d}}$.

Let us first consider small values of the coupling $\alpha$. Then our
expectation is borne out, and we observe the triplet forming spirals for the
mass and charge, while damped oscillations for the extremal values
$f_{\text{min}}$ and $\phi_{\rm max}$. The extremal values always reside at
the center of the configurations, $i.e.$, at the origin. Also the separation
$z_{\text{d}}$ between the neighboring components of the $k_{z}=2$ triplet
forms a spiral. Here, along this spiral, the contribution of the central BS to
the full configurations becomes dominant, while the outer BSs contributions
diminish.

When increasing the coupling constant beyond $\alpha\gtrsim 0.45$, however,
the scenario changes completely. While the configurations still follow a
similar pattern along the fundamental branch, the further part of the
$(\omega,M)$-diagram does not involve a spiraling/oscillating behavior.
Instead there is a single second branch, that leads all the way back to the
vacuum configuration. Thus all physical quantities exhibit loops, beginning
and ending at $\omega_{\text{max}}$. This may be interpreted as follows. Along
the second branch the configurations again feature a dominant central BS, but
the ’satellite’ BSs dissolve into a (sort of) boson shell. Moreover, the
system tends more and more towards spherical symmetry. Now we recall that a
central BS surrounded by a boson shell is precisely what constitutes a
radially excited spherically symmetric BS with a single radial node,
$n_{r}=1$. This suggests that the $k_{z}=2$ triplet might merge with a
$n_{r}=1$ single BS.

Let us therefore compare the $M-\omega$ diagram of the $k_{z}=2$ triplets at
large coupling $\alpha$ with that of the radially excited $n_{r}=1$ single
BSs. Such a comparison is shown in Fig. 5, where the $k_{z}=2$ BS triplets are
marked by solid curves and the radially excited $n_{r}=1$ single BSs by dashed
curves for two values of $\alpha$. The upper panels show the scaled ADM mass
$M$ (left) and the scaled charge $Q$ (right), while the lower panels show the
minimal value $f_{\text{min}}$ of the metric function $f$ (left) and the
maximal value $\phi_{\rm max}$ of the scalar field function $\phi$ (right).

Figure 6: Comparison of the radially excited $n_{r}=1$ single BSs (dashed
curves), and $k_{z}=2$ triplet of BSs (solid curves): scaled ADM mass $M$
(upper left panel), scaled charge $Q$ (upper right panel), minimal value of
the metric function $f_{\text{min}}$ (lower left panel), maximal value of the
scalar field function $\phi_{\text{max}}$ (lower right panel) $vs.$ frequency
$\omega$ for two values of the coupling $\alpha$. Note the quadratic scale for
$f_{\text{min}}$.

These figures are indeed very telling. The radially excited $n_{r}=1$ single
BSs exhibit the typical curve pattern of spherically symmetric BSs. They
emerge from the vacuum, form the fundamental branch and end in a
spiraling/oscillating pattern. The $k_{z}=2$ triplet likewise emerges from the
vacuum, forms a fundamental branch, and a second branch, but this second
branch of the triplet nicely overlaps with the fundamental branch of the
$n_{r}=1$ single BSs at some critical value of the frequency
$\omega_{\text{cr}}$. The overlap happens when mass and charge of the
$n_{r}=1$ single BSs have already passed their maximal values and the radially
excited stars are descending into the spiral.

It is now clear how the domain of existence of odd chains with more
constituents should be, for sufficiently large values of the coupling
$\alpha$. Let us consider a chain with $k_{z}$ nodes and thus $k_{z}+1$
constituents. Then this chain will feature on its fundamental branch $k_{z}+1$
more or less equal BSs, located on the symmetry axis. Subsequently the central
BS will start to dominate while the satellite BSs will dissolve into boson
shells. As the system tends towards spherical symmetry, it will overlap with a
$n_{r}=k_{z}/2$ single BS at some critical value of the frequency. We have
checked this behavior for the $k_{z}=4$ BS quintet, which indeed overlaps with
the radially excited $n_{r}=2$ single BS.

##  4 Conclusions

The main purpose of this paper was to report the existence of a new type of
solitonic configurations for a model with a gravitating self-interacting
complex scalar field. These configurations represent chains of BSs, with
$k_{z}+1$ constituents located symmetrically along the symmetry axis. The
number $k_{z}\geqslant 0$ represents the number of nodes on the symmetry axis
of the scalar field amplitude $\phi$.

The chains emerge from the vacuum $\phi=0$ at a maximal value of the boson
field frequency $\omega$, which is given by the field’s mass. For any $k_{z}$,
a fundamental branch of solutions is found emerging from this vacuum, in a
($\omega,M$)-diagram (with $M$ the ADM mass), ending at a minimal value of the
frequency $\omega_{\rm min}$, whose value is determined by the self-
interaction potential and the gravitational coupling strength $\alpha$. The
subsequent solution curve depends on the number of constituents and the
coupling $\alpha$. A single spherical BS has been argued to form an infinite
number of branches, leading to spirals or damped oscillations (depending on
the quantities considered) as its limiting configuration is approached. For
even chains we do not see such an endless spiraling/oscillating pattern.
Instead we observe only two to three branches, depending on the coupling
$\alpha$.

As the limiting configuration is approached the even chains retain basically
two of their $k_{z}+1$ constituents, whose metric function $g_{00}$ exhibits
two sharp peaks, reaching a very small value, while the scalar field features
two sharp opposite extrema located right at the location of these peaks. The
resulting configurations then feature huge second derivatives of the
functions, which impede further numerical analysis towards the limiting
solution.

The odd chains show a similar pattern as the single BS, when the coupling
$\alpha$ is small. This may be interpreted as the central BS dominating the
configurations on the higher branches. For larger $\alpha$, however, the
pattern changes totally, and the chains overlap on their second branch with a
radially excited spherical single BS with $n_{r}=k_{z}/2$ (radial) nodes. In
this case the central dominant BS will be surrounded by $n_{r}$ ‘boson
shells’.888 Boson shells without a central BS are also known [54, 55].

These solutions can be generalized in various directions. The most obvious
generalization is to include rotation. For the scalar field this means to
include an explicit harmonic dependence on the azimuthal angle $\varphi$. The
rotating single-BSs were obtained long ago [56, 57, 58, 40]. The rotating BSs
with odd parity and $k_{z}=1$, representing the rotating generalizations of
the pair of BSs, were also discussed in the literature [41]. We predict the
existence of rotating generalizations for the triplet and the higher chains
discussed in this work.

Single non-rotating BSs cannot be endowed with a black hole at the center; the
no-hair theorems forbid their existence [59]. This result is, however,
circumvented in the presence of spin, hairy generalizations of the Kerr black
hole (BH) with a complex scalar field being reported in literature [44, 60,
61, 62, 63, 64, 45, 65]. These hairy BHs obey a synchronization condition
relating the angular velocity of the event horizon and the boson field
frequency. Most studies so far considered only an even parity scalar field,
see $e.g.$ [66, 67, 68, 69, 70, 71, 72, 73, 74, 47, 75, 46, 76]. Synchronized
hairy BHs with an odd parity scalar field were obtained in [47, 46, 76]. While
these solutions represent only the simplest type of generalizations,
containing a single black hole at the center of configurations with one
(parity even) constituent or two (parity odd) constituents one can easily
image to put a black hole either at the center of rotating configurations with
more than two components, or to put a black hole at the center of each of the
components along the symmetry axis. Such configurations should correspond to
hairy double Kerr solutions or hairy multi-Kerr solutions. It would be
interesting to see, whether the presence of the scalar hair can regularize
such solutions, so that no conical singularity would be needed to hold them in
equilibrium.

Along similar lines, but replacing rotation by an electric charge, one could
investigate chains of electromagnetically charged BSs, generalizing the
results for a single charged BSs in Einstein-Klein-Gordon-Maxwell model [78,
54, 79]. Some results in this direction were reported in the recent work [82],
where (flat space) $Q$-chains were constructed in a model with a $U(1)$ gauged
scalar field, for a particular choice of the constants in the potential
(2.21). Gravitating generalizations of these solutions should also exist,
sharing some of the properties of the (ungauged) BS chains in this work. In
this context, we mention the recent results in [79, 80], showing that the no-
hair theorem in [81] does not apply to a single charged static BS in a model
with a $Q$-ball type potential (2.21), with the existence of BH
generalizations. For chains of charged BSs, it would be again of particular
interest to see whether they would support regular static multi-BH solutions.

## Acknowledgements

We gratefully acknowledge support by the DFG funded Research Training Group
1620 “Models of Gravity” and by the DAAD. We would also like to acknowledge
networking support by the COST Actions CA16104 and CA15117. Ya.S. gratefully
acknowledges the support by the Ministry of Sceince and High Education of
Russian Federation, project FEWF-2020-0003. and by the BLTP JINR Heisenberg-
Landau program 2020. This work is also supported by the Center for Research
and Development in Mathematics and Applications (CIDMA) through the Portuguese
Foundation for Science and Technology (FCT - Fundacao para a Ciência e a
Tecnologia), references UIDB/04106/2020 and UIDP/04106/2020 and by national
funds (OE), through FCT, I.P., in the scope of the framework contract foreseen
in the numbers 4, 5 and 6 of the article 23, of the Decree-Law 57/2016, of
August 29, changed by Law 57/2017, of July 19. We acknowledge support from the
projects PTDC/FIS-OUT/28407/2017, CERN/FIS-PAR/0027/2019 and PTDC/FIS-
AST/3041/2020. This work has further been supported by the European Union’s
Horizon 2020 research and innovation (RISE) programme H2020-MSCA-RISE-2017
Grant No. FunFiCO-777740. The authors would like to acknowledge networking
support by the COST Action CA16104.

## References

  * [1] G. Rosen, J. Math. Phys. 9 (1968) 996, 999
  * [2] R. Friedberg, T. D. Lee and A. Sirlin, Phys. Rev. D 13 (1976) 2739
  * [3] S. R. Coleman, Nucl. Phys. B 262 (1985) 263 Erratum: [Nucl. Phys. B 269 (1986) 744].
  * [4] D. J. Kaup, Phys. Rev. 172 (1968) 1331.
  * [5] R. Ruffini and S. Bonazzola, Phys. Rev. 187 (1969) 1767.
  * [6] L. G. Collodel, B. Kleihaus and J. Kunz, Phys. Rev. D 96 (2017) no.8, 084066
  * [7] C. A. R. Herdeiro, J. Kunz, I. Perapechka, E. Radu and Y. Shnir, Phys. Lett. B 812 (2021), 136027 [arXiv:2008.10608 [gr-qc]].
  * [8] B. Kleihaus and J. Kunz, Phys. Rev. D 61 (2000) 025003
  * [9] B. Kleihaus and J. Kunz, Phys. Rev. Lett. 85 (2000), 2430-2433
  * [10] B. Kleihaus, J. Kunz and Y. Shnir, Phys. Lett. B 570 (2003) 237
  * [11] B. Kleihaus, J. Kunz and Y. Shnir, Phys. Rev. D 68 (2003) 101701
  * [12] B. Kleihaus, J. Kunz and Y. Shnir, Phys. Rev. D 70 (2004) 065010
  * [13] B. Kleihaus, J. Kunz and Y. Shnir, Phys. Rev. D 71 (2005) 024013
  * [14] R. Teh and K. Wong, J. Math. Phys. 46 (2005), 082301
  * [15] V. Paturyan, E. Radu and D. Tchrakian, Phys. Lett. B 609 (2005), 360-366
  * [16] B. Kleihaus, J. Kunz and U. Neemann, Phys. Lett. B 623 (2005) 171
  * [17] J. Kunz, U. Neemann and Y. Shnir, Phys. Lett. B 640 (2006) 57
  * [18] J. Kunz, U. Neemann and Y. Shnir, Phys. Rev. D 75 (2007) 125008
  * [19] K. Lim, R. Teh and K. Wong, J. Phys. G 39 (2012), 025002
  * [20] R. Teh, A. Soltanian and K. Wong, Phys. Rev. D 89 (2014) no.4, 045018
  * [21] S. Krusch and P. Sutcliffe, J. Phys. A 37 (2004) 9037
  * [22] Y. Shnir and D. H. Tchrakian, J. Phys. A 43 (2010) 025401
  * [23] Y. Shnir, Phys. Rev. D 92 (2015) no.8, 085039
  * [24] B. Kleihaus, J. Kunz and M. Leissner, Phys. Lett. B 663 (2008), 438-444
  * [25] R. Ibadov, B. Kleihaus, J. Kunz and M. Leissner, Phys. Lett. B 663 (2008), 136-140
  * [26] R. Ibadov, B. Kleihaus, J. Kunz and M. Leissner, Phys. Lett. B 686 (2010), 298-306
  * [27] R. Ibadov, B. Kleihaus, J. Kunz and M. Leissner, Phys. Rev. D 82 (2010), 125037
  * [28] R. Teh, B. Ng and K. Wong, Annals Phys. 362 (2015), 170-195
  * [29] E. Radu and D. H. Tchrakian, Phys. Lett. B 636 (2006) 201
  * [30] Y. M. Shnir, EPL 77 (2007) no.2, 21001.
  * [31] O. Kichakova, J. Kunz, E. Radu and Y. Shnir, Phys. Rev. D 86 (2012) 104065.
  * [32] O. Kichakova, J. Kunz, E. Radu and Y. Shnir, Phys. Rev. D 90 (2014) no.12, 124012.
  * [33] R. M. Wald, General Relativity, (University of Chicago Press, Chicago, 1984).
  * [34] D. Kramer, H. Stephani, E. Herlt, and M. MacCallum, Exact Solutions of Einstein’s Field Equations, Cambridge University Press, Cambridge, (1980).
  * [35] T. D. Lee and Y. Pang, Phys. Rept. 221 (1992) 251.
  * [36] W. Deppert and E. W. Mielke, Phys. Rev. D 20 (1979) 1303.
  * [37] E. W. Mielke and R. Scherzer, Phys. Rev. D 24, 2111 (1981).
  * [38] R. Friedberg, T. Lee and Y. Pang, Phys. Rev. D 35 (1987), 3658
  * [39] M.S. Volkov and E. Wohnert, Phys. Rev. D 66 (2002) 085003.
  * [40] B. Kleihaus, J. Kunz and M. List, Phys. Rev. D 72 (2005) 064002 .
  * [41] B. Kleihaus, J. Kunz, M. List and I. Schaffer, Phys. Rev. D 77 (2008) 064025
  * [42] R. Battye and P. Sutcliffe, Nucl. Phys. B 590, (2000) 329
  * [43] P. Bowcock, D. Foster and P. Sutcliffe, J. Phys. A 42 (2009), 085403
  * [44] S. Hod, Phys. Rev. D 86 (2012) 104026 Erratum: [Phys. Rev. D 86 (2012) 129902]
  * [45] C. Herdeiro, E. Radu and H. Runarsson, Phys. Lett. B 739 (2014) 302
  * [46] J. Kunz, I. Perapechka and Y. Shnir, Phys. Rev. D 100 (2019) no.6, 064032
  * [47] Y. Q. Wang, Y. X. Liu and S. W. Wei, Phys. Rev. D 99 (2019) no.6, 064036
  * [48] Y. Brihaye and B. Hartmann, Phys. Rev. D 79 (2009) 064013
  * [49] A. Bernal, J. Barranco, D. Alic and C. Palenzuela, Phys. Rev. D 81 (2010) 044031
  * [50] N.I.M. Gould, J.A. Scott and Y. Hu, ACM Transactions on Mathematical Software 33 (2007) 10;
O. Schenk and K. Gärtner Future Generation Computer Systems 20 (3) (2004) 475.

  * [51] W. Schönauer and R. Weiß, J. Comput. Appl. Math. 27, 279 (1989) 279;
M. Schauder, R. Weiß and W. Schönauer, Universität Karlsruhe, Interner Bericht
Nr. 46/92 (1992).

  * [52] R. Camassa and D. D. Holm, Phys. Rev. Lett. 71 (1993), 1661
  * [53] J. Lis, Phys. Rev. D 84 (2011), 085030
  * [54] B. Kleihaus, J. Kunz, C. Lammerzahl and M. List, Phys. Lett. B 675 (2009), 102-115
  * [55] B. Kleihaus, J. Kunz, C. Lammerzahl and M. List, Phys. Rev. D 82 (2010), 104050
  * [56] F. E. Schunck and E. W. Mielke, Phys. Lett. A 249, 389 (1998).
  * [57] S. Yoshida and Y. Eriguchi, Phys. Rev. D 56, 762 (1997).
  * [58] F. D. Ryan, Phys. Rev. D 55, 6081 (1997).
  * [59] C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D 24 (2015) no.09, 1542014
  * [60] C. A. R. Herdeiro and E. Radu, Phys. Rev. Lett. 112 (2014) 221101
  * [61] C. Herdeiro and E. Radu, Phys. Rev. D 89 (2014) no.12, 124018
  * [62] C. Herdeiro and E. Radu, Class. Quant. Grav. 32 (2015) no.14, 144001
  * [63] S. Hod, Phys. Rev. D 90 (2014) no.2, 024051
  * [64] C. L. Benone, L. C. B. Crispino, C. Herdeiro and E. Radu, Phys. Rev. D 90, no. 10, 104024 (2014)
  * [65] C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D 23 (2014) no.12, 1442014
  * [66] B. Kleihaus, J. Kunz and S. Yazadjiev, Phys. Lett. B 744 (2015) 406
  * [67] C. A. R. Herdeiro, E. Radu and H. Runarsson, Phys. Rev. D 92 (2015) no.8, 084059
  * [68] C. Herdeiro, J. Kunz, E. Radu and B. Subagyo, Phys. Lett. B 748, 30 (2015)
  * [69] C. Herdeiro, E. Radu and H. Runarsson, Class. Quant. Grav. 33 (2016) no.15, 154001
  * [70] Y. Brihaye, C. Herdeiro and E. Radu, Phys. Lett. B 760 (2016) 279
  * [71] S. Hod, Phys. Lett. B 751 (2015) 177
  * [72] C. Herdeiro, J. Kunz, E. Radu and B. Subagyo, Phys. Lett. B 779, 151 (2018)
  * [73] C. Herdeiro, I. Perapechka, E. Radu and Y. Shnir, JHEP 1810 (2018) 119
  * [74] C. Herdeiro, I. Perapechka, E. Radu and Y. Shnir, JHEP 02 (2019), 111
  * [75] J. F. M. Delgado, C. A. R. Herdeiro and E. Radu, Phys. Lett. B 792 (2019), 436-444
  * [76] J. Kunz, I. Perapechka and Y. Shnir, JHEP 07 (2019), 109
  * [77] B. Kleihaus and J. Kunz, Phys. Lett. B 494 (2000), 130-134
  * [78] P. Jetzer and J. van der Bij, Phys. Lett. B 227 (1989), 341-346
  * [79] C. A. R. Herdeiro and E. Radu, Eur. Phys. J. C 80 (2020) no.5, 390 [arXiv:2004.00336 [gr-qc]].
  * [80] J. P. Hong, M. Suzuki and M. Yamada, Phys. Rev. Lett. 125 (2020) no.11, 111104 [arXiv:2004.03148 [gr-qc]].
  * [81] A. E. Mayo and J. D. Bekenstein, Phys. Rev. D 54 (1996), 5059-5069 [arXiv:gr-qc/9602057 [gr-qc]].
  * [82] V. Loiko, I. Perapechka and Y. Shnir, [arXiv:2012.01052 [hep-th]].


# Assimilation of SWOT Altimetry and Sentinel-1 Flood Extent Observations for
Flood Reanalysis - A Proof-of-Concept

###### Abstract

In spite of astonishing advances and developments in remote sensing
technologies, meeting the spatio-temporal requirements for flood hydrodynamic
modeling remains a great challenge for Earth Observation. The assimilation of
multi-source remote sensing data in 2D hydrodynamic models participates to
overcome such a challenge. The recently launched Surface Water and Ocean
Topography (SWOT) wide-swath altimetry satellite provides a global coverage of
water surface elevation at a high resolution. SWOT provides complementary
observation to radar and optical images, increasing the opportunity to observe
and monitor flood events. This research work focuses on the assimilation of 2D
flood extent maps derived from Sentinel-1 C-SAR imagery data, and water
surface elevation from SWOT as well as in-situ water level measurements. An
Ensemble Kalman Filter (EnKF) with a joint state-parameter analysis is
implemented on top of a 2D hydrodynamic TELEMAC-2D model to account for errors
in roughness, input forcing and water depth in floodplain subdomains. The
proposed strategy is carried out in an Observing System Simulation Experiment
based on the 2021 flood event over the Garonne Marmandaise catchment. This
work makes the most of the large volume of heterogeneous data from space for
flood prediction in hindcast mode paves the way for nowcasting.

Index Terms—  Fluvial floods, Data assimilation, EnKF, TELEMAC-2D, Garonne,
Sentinel-1, SWOT, OSSE

##  1 Introduction

The inclusion of Earth Observations (EO) from space into flood risk prediction
strategies presents a great opportunity to improve the ability to monitor,
anticipate flooding, mitigate its impacts, as well as protecting assets
worldwide. The properties of inland water bodies are monitored by altimetry
missions that provide along-track water surface elevation (WSE) from nadir
(e.g., TOPEX/Poseidon, Jason, SARAL/AltiKa, Sentinel-3, Sentinel-6) or large-
swath altimeters (SWOT), as well as from radar/optical satellites
(Sentinel-1/Sentinel-2) that provide high-resolution water extent maps.
Nevertheless, remote-sensing (RS) data can provide spatially dense information
but with a sparse revisit frequency. On the other hand, in-situ gauge
measurements offer a high temporal sampling but they are limited to low
spatial coverage. Neither suffice to fully capture both spatial and temporal
dynamics of flooding.

Data assimilation (DA) has emerged as an efficient tool to deal with
uncertainties in hydrology and hydrodynamics, such as in roughness, inflow
discharge, channel and floodplain geometry and/or hydraulic state. It combines
physics-based equations that govern the hydrodynamic models with ground
measurements and/or multi-mission satellite data that offer different
characteristics in space and time. Therefore, it allows to correct initial,
boundary conditions and model parameters to issue improved nowcasts and
forecasts. Leveraging multi-source flood observations allows densifying the
observing network, both spatially and temporally, as well as diversifying
their characteristics. Previous works [1, 2, 3] have shown that this allows
for a better performance of the EnKF—that relies on the stochastic computation
of forecast error covariance matrix among a limited number of perturbed
simulations—to represent the flow dynamics in the riverbed and floodplain.
Indeed, flood extent maps derived from Sentinel-1 (S1) C-band SAR imagery data
were assimilated jointly with in-situ water depth measurements from observing
stations of the VigiCrue network in France, over the Garonne Marmandaise
catchment; this improves the representation of the major flood events in 2019
and 2021 [2, 3]. Launched in December 2022, the SWOT wide-swath altimetry
satellite, with the Ka-band Radar Interferometer onboard, shall provide a
global coverage of WSE at a high resolution presented as rasters and pixel
clouds from which river products are delivered along the predefined river
centerlines.

The present work investigates the addition, to this workflow, of river
products derived from SWOT satellite at nodes located over the centerline of
the river. Given the recent launch of SWOT and the limited availability of
real SWOT river product, the experiment is carried out in an Observing System
Simulation Experiment (OSSE) framework using synthetical S1, SWOT and in-situ
observations generated from a reference simulation.

##  2 Method

Fig. 1: Observations for OSSE based on the 2021 event.

The study is carried out with the TELEMAC-2D111https://opentelemac.org/ (T2D)
hydrodynamic solver, over a 50-km reach of the Garonne River between Tonneins
and La Réole, for a major flood event in 2021. The upstream forcing is
provided by the observed stage height at Tonneins translated into discharge
using a local rating curve. The OSSE observing setting was based on a real
scenario for in-situ and S1 data, whereas the SWOT pass plan was enhanced,
compared to the 21-day repeat cycle of the nominal science orbit, by tripling
the number of passes during the flood period. This allows a temporal sampling
of 1-2 days for SWOT over the studied catchment, which is more ideal for flood
studies. Fig. 1 depicts the synthetical in-situ water-level for the 2021 flood
event at VigiCrue observing stations, namely Tonneins (blue), Marmande
(orange) and La Réole (green). The overpass times of S1 and SWOT are indicated
as vertical black dashed lines and yellow dash-dotted lines, respectively.

Fig. 2: Workflow of EnKF DA strategy in OSSE. Fig. 3: Generation of node-level
WSE SWOT river products on 2021-01-16. Color-coded aggregated pixels for each
node (left panel) and pixel class (right panel).

The data assimilation workflow in OSSE is presented in Fig. 2. An EnKF with a
joint state-parameter analysis is implemented on top of T2D [2]. The EnKF
control vector is composed of spatially-distributed friction coefficients
$K_{s}$, a corrective parameter of the inflow discharge $Q$, typically
necessary when dealing with highly uncertain inflow [4]. In addition, it is
augmented with a uniform correction of hydraulic states within the floodplain
subdomains $\delta H$. The ensemble of T2D simulations for the EnKF algorithm
are generated from a priori values for the friction and upstream forcing
perturbed with errors with gaussian statistical properties [1, 2]. The
settings of the OSSE’s reference simulation are based on a previous
assimilation of real in-situ data. This reference simulation are used to
generate synthetic observations, thanks to the dedicated observation
operators. Such operators are used to compute the model equivalent of the
observations for each member of the ensemble. True WSE values are extracted
from the reference simulation at all observation times and locations to
generate synthetical in-situ data at stream-gauge stations. On the other hand,
a threshold of 5 cm is applied on reference water depth maps to generate
synthetical flood extents at S1 overpass times, which are used to compute wet
surface ratios (WSR) in selected subdomains of the floodplain [2]. Such non-
Gaussian ratio is further treated with a Gaussian anamorphosis [3, 5].

The synthetical SWOT river product are generated from the reference water
depth maps at SWOT overpass times, first processed with the Tools4SWOTsims
toolbox [6]. This toolbox consists of a set of Python scripts to map 1D/2D
hydrodynamic model outputs into 2D WSE rasters that are compatible with SWOT-
HR simulator222https://github.com/CNES/swot-hydrology-toolbox, dedicated to
hydrology science [7]. The outputs of the SWOT-HR simulator are SWOT
observations in pixel cloud format. Such pixel cloud data is then processed by
the RiverObs333https://github.com/SWOTAlgorithms/RiverObs which aggregates the
WSEs from a selection of pixels to issue WSEs at nodes every 200 m along river
centerline and at river reach every 10 km (which are defined by the SWOT River
Database [8]). The selection of pixels is based on a class assignment related
to the nature of the water pixels (i.e. open water, water near land, dark
water, etc), and weighted by their respective uncertainties to issue node-
level WSE with an error below 10 cm as prescribed in the SWOT requirements.
The details for SWOT WSE, width, and slope computation are given in SWOT River
Single Pass Product Description Document [9] and example data
products444https://podaac.jpl.nasa.gov/swot?tab=datasets.

Fig. 3 depicts a typical SWOT node-level river product over a meander of a
reach of the Garonne River observed on 2021-01-16, each node is an aggregation
of pixels. The pixels are color-coded either by their associated node on the
left panel, or by their class on the right panel. For each member of the
ensemble, the model equivalent of the SWOT node-level WSE is computed using
the SWOT-dedicated observation operator. It is the average of the WSEs of
certain associated pixels selected according to their quality flag; each pixel
WSE is interpolated from the T2D mesh of simulated water depth.

##  3 Experimental Results

Four EnKF DA experiments are carried out considering different combinations of
in-situ, S1 and SWOT synthetical observations, as presented in Table 1. In the
following, the deterministic simulation using the control vector true values
is represented as Truth, and the simulation without DA using the calibrated
parameters is called Open Loop (OL). IDA assimilates only in-situ water-level
measurements, IGDA assimilates in-situ WSE and S1-derived WSR observations,
FDA assimilates all available flood observations whereas RSDA assimilates only
remote-sensing data.

Quantitative performance assessments are carried out in the observational
spaces by comparing the DA analysis WSE time-series at observing stations with
the true water-level observations, and by the means of the resulting
contingency maps and the overall Critical Success Index (CSI) metric computed
for the DA analysis flood extent maps, with respect to the truth flood extent
maps.

Table 1: Summary of experimental settings. | Assimilated observations |
---|---|---
Exp. | In-situ | Sentinel-1 | SWOT |
name | WSE | WSR | WSE | Control vector
OL | $\square$ | $\square$ | $\square$ | -
IDA |  $\square$ $\checkmark$ | $\square$ | $\square$ | $K_{s},Q$
IGDA |  $\square$ $\checkmark$ |  $\square$ $\checkmark$ | $\square$ | $K_{s},Q,\delta H$
RSDA | $\square$ |  $\square$ $\checkmark$ |  $\square$ $\checkmark$ | $K_{s},Q,\delta H$
FDA |  $\square$ $\checkmark$ |  $\square$ $\checkmark$ |  $\square$ $\checkmark$ | $K_{s},Q,\delta H$
Fig. 4: WSEs at Tonneins for all experiments (top panel), and errors computed
with respect to the Truth WSE (bottom panel).

Fig. 4 depicts the WSE simulated at Tonneins from the five experiments on the
top panel, and the anomalies to the truth on the bottom panel. The SWOT-
derived WSE at the closest nodes to Tonneins are represented by black crosses.
The RMSEs of the WSE at all three observing stations computed over the entire
flood event are provided in Table 2. The assimilation of in-situ data with
high temporal sampling suffices to represent the WSE (in the riverbed) at the
gauge stations. It reduces the RMSE from 26 cm in the OL to 5-7 cm in IDA,
IGDA and FDA.

The improvement brought by FDA demonstrates the feasibility of multi-sources
DA. It should be noted that for RSDA, in-situ observations are independent
observations. The discrepancies between in-situ and SWOT WSE are most likely
due to the space aggregation of SWOT pixels, thus the assimilation of SWOT
does not guarantee an improvement with respect to the WSEs at gauge stations.
In RSDA, the correction at Tonneins is implicitly influenced by misfits to
SWOT at various locations, including those precede Tonneins. Here, the
assimilation of SWOT degrades the results with respect to in-situ observations
at Tonneins. These results from RSDA should be further investigated, with a
focus on the consistency between the SWOT node-level WSE and the truth WSE.
The contingency maps and CSI computed with respect to the truth flood extents
are shown in Fig. 5 for all experiments. OL tends to heavily underestimate the
flooding at the flood peak and overestimate the flooding during recess. The
assimilation of synthetical S1 WSR and/or SWOT node-level WSE (in IGDA, RSDA
and FDA) improves the representation of the floodplain dynamics and thus have
better CSI score with respect to OL and IDA.

Fig. 5: Contingency maps for all experiments with respect to truth flood extents at flood peak and recess. Table 2: Root-mean-square errors of simulated WSE at observing stations. RMSE [m] | Tonneins | Marmande | La Réole
---|---|---|---
OL | 0.260 | 0.398 | 0.578
IDA | 0.052 | 0.042 | 0.053
IGDA | 0.054 | 0.044 | 0.049
RSDA | 0.410 | 0.435 | 0.413
FDA | 0.075 | 0.056 | 0.053

##  4 Conclusions and Perspectives

This paper demonstrates the feasibility of multi-sources DA and presents the
merits of assimilating different types of flood observations, in the framework
of an EnKF built on top of a T2D hydrodynamic model. The FDA experiment that
assimilates all observations presents moderate performance, suggesting
potential inconsistency between the in-situ and SWOT-nodes observations. It
also shows that the assimilation of remote-sensing only observations (in RSDA)
fails to improve the OL results in the riverbed, in exchange for improving the
floodplain dynamics. This trade-off will be investigated in further studies
with a spatial analysis along the river SWOT nodes. Lastly, the good
performance of the joint assimilation (in IGDA) of in-situ WSE data and
S1-derived flood extent maps demonstrates the benefits of these two data
sources, which are the ones with the lowest temporal sampling in real-world
applications. This work heralds toward a reliable methodology for flood
prediction for poorly-gauged catchments, making the most of innovative EO data
such as real SWOT data.

## References

  * [1] T. H. Nguyen, S. Ricci, C. Fatras, A. Piacentini, A. Delmotte, E. Lavergne, and P. Kettig, “Improvement of flood extent representation with remote sensing data and data assimilation,” IEEE Transactions on Geoscience and Remote Sensing, vol. 60, pp. 1–22, 2022.
  * [2] T. H. Nguyen, S. Ricci, A. Piacentini, C. Fatras, P. Kettig, G. Blanchet, S. Peña Luque, and S. Baillarin, “Dual state-parameter assimilation of sar-derived wet surface ratio for improving fluvial flood reanalysis,” Water Resources Research, vol. 58, no. 11, p. e2022WR033155.  e2022WR033155 2022WR033155.
  * [3] T. H. Nguyen, S. Ricci, A. Piacentini, E. Simon, R. Rodriguez Suquet, and S. Peña Luque, “Gaussian Anamorphosis for Ensemble Kalman Filter Analysis of SAR-Derived Wet Surface Ratio Observations,” IEEE Transactions on Geoscience and Remote Sensing, vol. 62, pp. 1–21, 2024.
  * [4] T. H. Nguyen, S. Ricci, A. Piacentini, Q. Bonassies, R. Rodriquez Suquet, S. Peña Luque, K. Marlis, and C. H. David, “Reducing uncertainties of a chained hydrologic-hydraulic models to improve flood forecasting using multi-source earth observation data,” in IGARSS 2023 - 2023 IEEE International Geoscience and Remote Sensing Symposium, pp. 1525–1528, 2023.
  * [5] T. H. Nguyen, S. Ricci, A. Piacentini, E. Simon, R. R. Suquet, and S. Peña Luque, “Dealing with non-gaussianity of sar-derived wet surface ratio for flood extent representation improvement,” in IGARSS 2023 - 2023 IEEE International Geoscience and Remote Sensing Symposium, pp. 1595–1598, 2023.
  * [6] C. Emery, S. Ricci, and A. Piacentini, “Tools4SWOTsims and SMURF–Python Libraries Supporting the Generation of Synthetical SWOT-like data in Anticipation of Real Products For River Hydrodynamics Studies with the Mascaret Hydraulic Model,” in AGU Fall Meeting Abstracts, vol. 2022, pp. OS22A–21, 2022.
  * [7] R. P. d. M. Frasson, R. Wei, M. Durand, J. T. Minear, A. Domeneghetti, G. Schumann, B. A. Williams, E. Rodriguez, C. Picamilh, C. Lion, T. Pavelsky, and P.-A. Garambois, “Automated river reach definition strategies: Applications for the surface water and ocean topography mission,” Water Resources Research, vol. 53, no. 10, pp. 8164–8186, 2017.
  * [8] G. H. Allen and T. M. Pavelsky, “Global extent of rivers and streams,” Science, vol. 361, no. 6402, pp. 585–588, 2018.
  * [9] C. Stuurman and C. Pottier, “SWOT product description: Level-2 KaRIn high rate river single pass vector product,” Technical Report, California Institute of Technology, 2020.