pdftexdestination with the same hyperrefOption hyperrefToken pdftex(dest)

# Effect of the Solar dark matter wake on planets

Indranil Banik1 and Pavel Kroupa1,2
1Helmholtz-Institut für Strahlen und Kernphysik (HISKP), University of Bonn,
Nussallee 14$-$16, D-53115 Bonn, Germany
2Charles University, Faculty of Mathematics and Physics, Astronomical
Institute, V Holešovičkách 2, CZ-18000 Praha 8, Czech Republic Email:
<EMAIL_ADDRESS>(Indranil Banik)
<EMAIL_ADDRESS>(Pavel Kroupa)

###### Abstract

The Galaxy is conventionally thought to be surrounded by a massive dark matter
(DM) halo. As the Sun goes through this halo, it excites a DM wake behind it.
This local asymmetry in the DM distribution would gravitationally affect the
motions of Solar System planets, potentially allowing the DM wake to be
detected or ruled out. Hernandez (2019) recently calculated that the DM-
induced perturbation to Saturn’s position is 252 metres net of the effect on
the Sun. No such anomaly is seen in Saturn’s motion despite very accurate
tracking of the Cassini spacecraft, which orbited Saturn for ¿13 years.

Here, we revisit the calculation of how much Saturn would deviate from
Keplerian motion if we fix its position and velocity at some particular time.
The DM wake induces a nearly resonant perturbation whose amplitude grows
almost linearly with time. We show that the Hernandez (2019) result applies
only for an observing duration comparable to the ${\approx 250}$ million year
period of the Sun’s orbit around the Galaxy. Over a 100 year period, the
perturbation to Saturn’s orbit amounts to ¡1 cm, which is quite consistent
with existing observations. Even smaller perturbations are expected for the
terrestrial planets.

###### keywords:

ephemerides – celestial mechanics – space vehicles – dark matter – gravitation
– planets and satellites: dynamical evolution and stability

††pubyear: 2019††pagerange: Effect of the Solar dark matter wake on
planets–References

##  1 Introduction

Over the last century, one of the big ‘elephants in the room’ for astronomers
is the fact that very large dynamical discrepancies often exist between the
observed rotation curves of galaxies and the predictions of Newtonian gravity
applied to their luminous matter distributions (e.g. Babcock, 1939; Rubin &
Ford, 1970; Rogstad & Shostak, 1972). These acceleration discrepancies are
usually attributed to halos of dark matter (DM) surrounding each galaxy
(Ostriker & Peebles, 1973). However, the discrepancies follow some remarkable
regularities (Famaey & McGaugh, 2012) that can be summarised as a unique
relation between the acceleration inferred from the rotation curve and that
expected from the baryonic distribution (McGaugh et al., 2016). Such a radial
acceleration relation was predicted several decades earlier using Milgromian
dynamics (MOND, Milgrom, 1983). In this model, the dynamical effects usually
attributed to DM are instead provided by an acceleration dependence of the
gravity law.

It is important to test both MOND and DM in regimes different to those which
gave rise to the paradigms in the first place. For MOND, the dynamics of wide
binary stars is a promising ‘experimentum crucis’ in the near future (e.g.
Banik & Zhao, 2018). If the observed dynamical discrepancies are due to DM,
then massive objects moving through the Galactic DM halo should experience
dynamical friction (Chandrasekhar, 1943). This is also true for the Sun, which
would be expected to create a trailing DM wake (Hernandez, 2019). Because the
DM would be overdense behind the Sun, its distribution would not be
spherically symmetric, contrary to the assumption of several previous works
(e.g. Pitjeva & Pitjev, 2013). This asymmetry might allow for much stronger
constraints on the DM density. Following this idea, Hernandez (2019) found
that the Sun’s DM wake would cause the position of Saturn to deviate by 252
metres. Such an effect is not seen in Cassini radio tracking data despite it
being accurate to 32 metres (Viswanathan et al., 2017). This led Hernandez
(2019) to rule out the DM hypothesis.

In this work, we conduct more detailed calculations of how Saturn would
deviate from Keplerian expectations if the DM wake is present (Section 2). Our
results are shown in Figure 2 for a century of observations, assuming Saturn’s
initial position and velocity are known exactly. We derive a perturbation
amplitude many orders of magnitude smaller than suggested by Hernandez (2019).
In Section 3, we explain why this is. Our conclusions are given in Section 4.

##  2 Methods and results

Figure 1: The Solar System geometry assumed in our work. Saturn orbits the Sun
on a nearly circular orbit with instantaneous heliocentric position $\bm{R}$,
which defines an orbital phase angle $\phi$ relative to some initial
direction. We define this to be our $x$-axis and take it to be the direction
within Saturn’s orbital plane most closely aligned with $\widehat{\bm{w}}$,
the velocity of the Sun with respect to the local DM. A non-zero velocity
would cause a density enhancement towards $-\widehat{\bm{w}}$, the so-called
‘DM wake’. Its density is axisymmetric with respect to $\widehat{\bm{w}}$,
making the DM-induced tidal acceleration of Saturn very nearly parallel to
this direction (see text).

In this contribution, we use a co-ordinate system in which the $xy$ plane
corresponds to the orbital plane of Saturn, whose orbital pole is along $z$
and instantaneous position is $\bm{R}$. Assuming near-circular motion, we
consider the orbit of Saturn in the epicyclic approximation. If its motion
were purely Keplerian, its position and velocity would imply that its guiding
centre radius (semi-major axis) is $R_{0}$. At any time $t$, its orbital
radius $R\equiv\left|\bm{R}\right|$ is slightly larger by an amount $r$, whose
evolution is governed by

$\displaystyle\ddot{r}\,+\,\Omega^{2}r~{}=~{}0\,,$ (1)

where $\Omega^{2}=\frac{GM_{\odot}}{{R_{0}}^{3}}$ is the Keplerian orbital
frequency of Saturn and $\dot{q}\equiv\frac{dq}{dt}$ for any quantity $q$.

In Figure 1, we show the relation between Saturn’s orbital plane and the DM
wake, which we take to lie in the direction
$\widehat{\bm{w}}\equiv\left(\cos\theta_{0},0,\sin\theta_{0}\right)\propto\bm{v}_{\odot}$.
In our notation, $\theta_{0}$ is the minimum angle between any vector within
the orbital plane of Saturn and the velocity of the Sun, $\bm{v}_{\odot}$,
with respect to the local DM. We assume this DM has no ordered motion in the
Galactocentric frame, implying that $\bm{v}_{\odot}$ can be measured in this
frame.

The Sun moves through the Galactic DM halo with a circular velocity of
$v_{c,\odot}$ in addition to some non-circular velocity
$\left(U_{\odot},V_{\odot},W_{\odot}\right)$. This is defined in the usual
Galactic Cartesian frame in which $x$ points towards the Galactic Centre, $z$
towards the North Galactic Pole and $y$ in the direction necessary to make the
co-ordinate system right-handed. Fortunately, $y$ points along the local
direction of the large scale ordered rotation of the Galactic disk. Thus, the
Sun has a Galactocentric velocity

$\displaystyle\bm{v}_{\odot}~{}\equiv~{}\begin{bmatrix}U_{\odot}\\\
V_{\odot}+v_{c,\odot}\\\ W_{\odot}\end{bmatrix}\,.$ (2)

In the rest of this contribution, we use our previously mentioned co-ordinate
system aligned with the orbital plane of Saturn. As the direction
$\widehat{\bm{R}}\equiv\bm{R}/R$ towards Saturn changes over the course of its
orbit, the angle $\theta$ between $\widehat{\bm{w}}$ and $\widehat{\bm{R}}$ is

$\displaystyle\cos\theta~{}\equiv~{}\widehat{\bm{w}}\cdot\widehat{\bm{R}}~{}=~{}\cos\theta_{0}\cos\phi\,,$
(3)

where $\phi=\Omega t$ is the orbital phase angle of Saturn such that
$\bm{R}\equiv R\left(\cos\phi,\sin\phi,0\right)$. This is valid for a nearly
circular orbit on which $\cos\theta$ has a maximum at $t=0$. In this case, the
component of $\bm{R}$ parallel to $\widehat{\bm{w}}$ is

$\displaystyle
w~{}\equiv~{}\widehat{\bm{w}}\cdot\bm{R}~{}=~{}R\cos\theta~{}=~{}R\cos\theta_{0}\cos\phi\,.$
(4)

The gravity due to the DM wake is approximately parallel to $\bm{v}_{\odot}$
(Mulder, 1983). This is due to the large-scale asymmetry caused by the motion
of the Sun. Some force is also expected in the direction orthogonal to
$\bm{v}_{\odot}$, but due to axisymmetry such a force must vanish very close
to the Sun. This is not true for a force parallel to $\bm{v}_{\odot}$, which
we therefore take to be the dominant effect of the DM wake.

After subtracting the gravity exerted by the DM wake on the Sun, the residual
acceleration of a test particle close to it can be approximated as (equation
AIV.13 in Mulder, 1983)

$\displaystyle\bm{F}_{w}$ $\displaystyle=$
$\displaystyle-\left(F_{1}+F_{0}\frac{\cos\theta}{\left|\cos\theta\right|}\right)\widehat{\bm{w}}\,,\,\text{where}$
(5) $\displaystyle\frac{F_{1}}{F_{typ}}$ $\displaystyle=$ $\displaystyle
0.21\ln\left(\frac{R}{2r_{min}}\right)\,,$ (6)
$\displaystyle\frac{F_{0}}{F_{typ}}$ $\displaystyle=$ $\displaystyle
0.44\,,\,\text{with}$ (7) $\displaystyle F_{typ}$ $\displaystyle=$
$\displaystyle\frac{4\mathrm{\pi}G^{2}\rho_{{}_{DM}}M_{\odot}}{\sigma^{2}}\quad\text{and}$
(8) $\displaystyle\sqrt{\frac{GM_{\odot}}{r_{min}}}$ $\displaystyle=$
$\displaystyle v_{\odot}\,,$ (9)

where $G$ is Newton’s gravitational constant, $\rho_{{}_{DM}}$ is the
unperturbed local density of DM, $\sigma$ is its one-dimensional velocity
dispersion, $v_{\odot}\equiv\left|\bm{v}_{\odot}\right|$ and $r_{min}$ is the
heliocentric distance at which $v_{\odot}$ is equal to the circular velocity
about the Sun (equation 11 in Mulder, 1983). Note that these results only hold
if the unperturbed DM has a Maxwell-Boltzmann velocity distribution with
$\sigma=v_{\odot}/\sqrt{2}$ (Hernandez, 2019). In this case, the most likely
DM velocity is $v_{\odot}\approx 250$ km/s, very close to the values reported
for Milky Way analogues in cosmological hydrodynamical simulations (Bozorgnia
& Bertone, 2017, figure 3). The numerical coefficients in Equation 5 need to
be adjusted for a different $v_{\odot}/\sigma$ (appendix IV in Mulder, 1983).

When estimating the DM wake-induced displacements, we adopt the parameter
values listed in Table 1. At the orbital radius of Saturn, this leads to a
typical tidal acceleration of

$\displaystyle F_{typ}~{}=~{}4.5\times 10^{-21}\,m/s^{2}\,.$ (10)

This is subject to a ${\approx 10\%}$ uncertainty due to imprecise knowledge
of the DM density (Hagen & Helmi, 2018). Its uncertain velocity distribution
also influences our results, with simulations indicating deviations from our
assumed Maxwell-Boltzmann form (e.g. Bozorgnia & Bertone, 2017). There is also
a ${\approx 1\%}$ uncertainty in $v_{\odot}$, corresponding to a ${\approx
2\%}$ uncertainty in $F_{typ}$.

Parameter | Meaning | Value and units
---|---|---
R | Orbital radius of Saturn | 9.58 AU
$v_{c,\odot}$ | See Equation 2 | 233.3 km/s
$U_{\odot}$ | See Equation 2 | 11.1 km/s
$V_{\odot}$ | See Equation 2 | 12.24 km/s
$W_{\odot}$ | See Equation 2 | 7.25 km/s
$\rho_{{}_{DM}}$ | Local dark matter density | 0.018 $M_{\odot}/pc^{3}$
Table 1: Our adopted values of the parameters relevant to this contribution.
We obtain $v_{c,\odot}$ from McGaugh (2018) and the non-circular velocity of
the Sun from Schönrich et al. (2010). The resulting Solar velocity of 245.9
km/s (Equation 2) implies that $\theta_{0}=60.6^{\circ}$ for a planet orbiting
within the Ecliptic (Figure 1). The local DM density is constrained by
kinematic observations of the Solar neighbourhood (Xia et al., 2016; Hagen &
Helmi, 2018). We assume that the unperturbed DM has a Maxwell-Boltzmann
velocity distribution with local one-dimensional velocity dispersion
$\sigma=v_{\odot}/\sqrt{2}=173.9$ km/s.

Over a period of $T=10$ years, the maximum DM-induced deviation is $\approx
F_{typ}\,T^{2}/2=0.22$ mm. Given that the ephemerides of Saturn are accurate
to ${\approx 30}$ metres (Viswanathan et al., 2017), it is unlikely that the
motion of Saturn would be noticeably affected by the Solar DM wake.
Nonetheless, we determine the expected perturbation more accurately in order
to reveal its ‘resonant’ nature (Hernandez, 2019). The precise meaning of this
will become clear in the following sections.

###  2.1 Motion within the orbital plane

The force from the DM wake varies only a little for a planet on a near-
circular orbit with $R_{0}\gg r_{min}$, as is the case for all Solar System
planets. Thus, Equation 5 approximately leads to a wake potential within the
orbital ($xy$) plane of

$\displaystyle\Phi_{w}$ $\displaystyle=$
$\displaystyle\overbrace{F_{1}w}^{\Phi_{w,1}}+F_{0}\frac{\cos\theta}{\left|\cos\theta\right|}w$
(11) $\displaystyle=$ $\displaystyle
R\cos\theta_{0}\left(F_{1}\cos\phi+F_{0}\left|\cos\phi\right|\right)$ (12)

On the last line, we restrict attention to the behaviour of $\Phi_{w}$ within
the $xy$ plane.

We now use Chapter 3.3.3 of Binney & Tremaine (2008) to consider the evolution
of Saturn’s orbit in a weakly non-axisymmetric potential. The derivation there
requires decomposition of $\Phi_{w}$ into Fourier modes $\propto\cos m\phi$.
The $F_{1}$ term contributes only to the $m=1$ mode while the $F_{0}$ term
contributes to all modes with even $m$, including the case $m=0$. For reasons
that we clarify in Section 3.1, we consider only the $m=0$ and $m=1$ modes,
which we denote $\Phi_{w,0}$ and $\Phi_{w,1}$, respectively. Using equation
3.146 in Binney & Tremaine (2008), we get that

$\displaystyle\ddot{r}+\Omega^{2}r$ $\displaystyle=$
$\displaystyle-\left(\frac{\partial\Phi_{w,1}}{\partial
R}+\overbrace{\frac{2\Omega\Phi_{w,1}}{R\widetilde{\Omega}}}\right)\cos\widetilde{\Omega}t$
(13)
$\displaystyle-\cos\theta_{0}F_{0}\langle\left|\cos\phi\right|\rangle\,,\text{
where}$ $\displaystyle\widetilde{\Omega}$ $\displaystyle\equiv$
$\displaystyle\Omega-\epsilon\,.$ (14)

The term marked with an overbrace arises because the angular momentum
oscillates with time, but this effect is only caused by the non-axisymmetric
part of $\Phi_{w}$. To handle the axisymmetric ($m=0$) part, we use the
superposition principle, exploiting the fact that the governing equations are
linear for small perturbations. The mean value of $\left|\cos\phi\right|$ is
$\langle\left|\cos\phi\right|\rangle=2/\mathrm{\pi}$. The wake potential
oscillates with a very slow frequency $\epsilon\ll\Omega$ due to the Galactic
orbit of the Sun. Thus, a very accurate approximation to Equation 13 is

$\displaystyle\ddot{r}+\Omega^{2}r$ $\displaystyle=$
$\displaystyle-\cos\theta_{0}\left(3F_{1}\cos\widetilde{\Omega}t+\frac{2}{\mathrm{\pi}}F_{0}\right)\,.$
(15)

We now define scaled force constants

$\displaystyle\widetilde{F}_{0}$ $\displaystyle=$
$\displaystyle-\frac{2\cos\theta_{0}}{\mathrm{\pi}}F_{0}\,,$ (16)
$\displaystyle\widetilde{F}_{1}$ $\displaystyle=$
$\displaystyle-3\cos\theta_{0}F_{1}\,.$ (17)

Thus, we get that

$\displaystyle\ddot{r}+\Omega^{2}r~{}=~{}\widetilde{F}_{1}\cos\widetilde{\Omega}t+\widetilde{F}_{0}\,.$
(18)

The particular integral of this can be guessed as

$\displaystyle
r~{}=~{}C\cos\widetilde{\Omega}t+\frac{\widetilde{F}_{0}}{\Omega^{2}}\,.$ (19)

This is a valid solution if the oscillation amplitude is

$\displaystyle
C~{}=~{}\frac{\widetilde{F}_{1}}{\Omega^{2}-{\widetilde{\Omega}}^{2}}\,.$ (20)

Now, suppose we start observing the orbit of Saturn at some time ${t=0}$. Its
position and velocity may have been affected by $\Phi_{w}$ at earlier times,
but we have no way of knowing this due to the lack of prior observations.
Thus, all we see is that Saturn behaves as if it has some $r$ and $\dot{r}$ at
$t=0$, which defines a particular (osculating) Keplerian orbit. To detect the
DM wake, it is necessary that the behaviour of $r$ deviate from simple
harmonic motion satisfying $\ddot{r}+\Omega^{2}r=0$ with $r$ and $\dot{r}$
having their observed values at ${t=0}$. Thus, we define a further radial
perturbation ${\delta r}$ which must satisfy $\delta r=\dot{\delta r}=0$ when
${t=0}$. In standard Keplerian motion, ${\delta r=0}\forall t$.

We have seen that Equation 19 solves the governing Equation 18, but it does
not satisfy this initial condition. To ensure that it does so, we add
appropriate multiples of the complementary functions $\cos\Omega t$ and
$\sin\Omega t$.

$\displaystyle\delta
r~{}=~{}\frac{\widetilde{F}_{0}}{\Omega^{2}}\left(1-\cos\Omega
t\right)\,+\,\frac{\widetilde{F}_{1}}{\Omega^{2}-{\widetilde{\Omega}}^{2}}\left(\cos\widetilde{\Omega}t-\cos\Omega
t\right)\,.$ (21)

We now take the limit that ${\epsilon\to 0}$, which implies that
${\widetilde{\Omega}\to\Omega}$. This is valid because $\epsilon$ corresponds
to the Galactic orbit of the Sun. Our observations span only a tiny fraction
of this ${\approx 250}$ Myr period, making it safe to assume that ${\epsilon
t\ll 1}$. In this limit, Equation 21 becomes

$\displaystyle\delta
r~{}=~{}\frac{\widetilde{F}_{0}}{\Omega^{2}}\left(1-\cos\Omega
t\right)\,+\,\frac{\widetilde{F}_{1}}{2\Omega}t\sin\Omega t\,.$ (22)

Our derivation could alternatively have been done by assuming ${\epsilon=0}$
from the outset. In this case, we would need to solve an equation of the form
$\ddot{x}+\Omega^{2}x=\cos\Omega t$. Given the initial conditions
$x=\dot{x}=0$, the solution to this is $x=t\sin\Omega t/\left(2\Omega\right)$,
which would yield exactly the same result. Thus, the actual value of
$\epsilon$ is completely irrelevant to our analysis as long as $\epsilon$ is
sufficiently small that ${\epsilon T\ll 1}$ over the period $T$ covered by
accurate observations.

To simplify the notation, we define oscillation amplitudes

$\displaystyle
x_{i}~{}\equiv~{}\frac{\widetilde{F}_{i}}{\Omega^{2}}\,,~{}i=0,1\,.$ (23)

Combining our previous results, it is easy to see that

$\displaystyle x_{0}$ $\displaystyle=$
$\displaystyle-\frac{2\cos\theta_{0}F_{0}}{\mathrm{\pi}\Omega^{2}}=-0.014\,\text{mm,}$
(24) $\displaystyle x_{1}$ $\displaystyle=$
$\displaystyle-\frac{3\cos\theta_{0}F_{1}}{\Omega^{2}}=-0.18\,\text{mm}.$ (25)

In terms of these $x_{i}$, Equation 22 can be written as

$\displaystyle\delta
r~{}=~{}x_{0}\left(1-\cos\phi\right)\,+\,\frac{x_{1}}{2}\phi\sin\phi\,.$ (26)

The $x_{i}$ are more than just convenient shorthand $-$ they capture the very
essence of our whole derivation. If the observing duration is comparable to
the orbital period $2\mathrm{\pi}/\Omega$, we expect to see a perturbation of
order ${F_{i}/\Omega^{2}}$. Thus, the $x_{i}$ tell us roughly how large the
deviation from Keplerian motion would be after a significant fraction of the
planetary orbital period. This also follows from assuming $\delta r\approx
F_{1}t^{2}/2$ over a duration $t=1/\Omega$, a typical orbital timescale. Over
longer periods, $\delta r\propto t$ due to the effects of orbital mechanics
(Equation 26).

####  2.1.1 Tangential motion

Similarly to ${\delta r}$, we define a perturbation ${\delta p}$ parallel to
the circular velocity of Saturn. ${\delta p}$ is governed by equation 3.145 of
Binney & Tremaine (2008).

$\displaystyle\dot{\delta p}~{}=~{}-2\Omega\delta
r-\frac{\widetilde{F}_{1}}{\Omega}\left(\cos\Omega t-1\right)\,.$ (27)

As with our solution for $\delta r$, the arbitrary additive constant is chosen
to ensure that $\dot{\delta p}=0$ at $t=0$. Thus, we get that

$\displaystyle\delta p$ $\displaystyle=$ $\displaystyle
2x_{0}\left(\sin\phi-\phi\right)+x_{1}\left(\left(\phi\cos\phi-\sin\phi\right)-\left(\sin\phi-\phi\right)\right)$
(28) $\displaystyle=$ $\displaystyle
2x_{0}\left(\sin\phi-\phi\right)+x_{1}\left(\phi\cos\phi-2\sin\phi+\phi\right)\,.$

###  2.2 Vertical motion

Applying Equation 5 in the $z$-direction orthogonal to the orbital plane of
Saturn, we get that

$\displaystyle
F_{z}~{}=~{}-\sin\theta_{0}\left(F_{1}+F_{0}\operatorname{sign}\left(\cos\phi\right)\right)$
(29)

The oscillatory term involving $F_{0}$ has a mean value of zero but clearly
contributes to the ${m=1}$ mode, which leads to a resonant perturbation
(Section 2.1). Its amplitude can be found by Fourier transforming
$\operatorname{sign}\left(\cos\phi\right)$, which yields a first non-zero term
of $\left(4\cos\phi\right)/\mathrm{\pi}$. Neglecting ${m>1}$ terms as before
and noting that the vertical epicyclic frequency is equal to the circular
orbit frequency $\Omega$, we get that

$\displaystyle\ddot{z}+\Omega^{2}z$ $\displaystyle=$ $\displaystyle
F_{z,0}+F_{z,1}\cos\Omega t\,,\,\text{where}$ (30) $\displaystyle F_{z,0}$
$\displaystyle=$
$\displaystyle-\sin\theta_{0}F_{1}=\frac{\tan\theta_{0}\widetilde{F}_{1}}{3}\,~{}\text{and}$
(31) $\displaystyle F_{z,1}$ $\displaystyle=$
$\displaystyle-\frac{4\sin\theta_{0}F_{0}}{\mathrm{\pi}}=2\tan\theta_{0}\widetilde{F}_{0}\,.$
(32)

This is directly analogous to the resonant limit of Equation 18, whose
solution we have just found (Equation 22). Therefore, the vertical
displacement is

$\displaystyle\delta z$ $\displaystyle=$
$\displaystyle\frac{F_{z,0}}{\Omega^{2}}\left(1-\cos\Omega
t\right)\,+\,\frac{F_{z,1}}{2\Omega}t\sin\Omega t$ (33) $\displaystyle=$
$\displaystyle\tan\theta_{0}\left(\frac{x_{1}}{3}\left(1-\cos\phi\right)\,+\,x_{0}\phi\sin\phi\right)\,.$
(34)

###  2.3 The total displacement

In Figure 2, we show the expected deviation of Saturn from Keplerian motion
over a period of 100 years. Presently, our observations span approximately
half its orbital period of 30 years thanks to the Cassini mission (Matson,
1992). Given an observing accuracy of ${\approx 30}$ metres (Viswanathan et
al., 2017), it is clear that the DM wake cannot be detected using current
Solar System observations. Indeed, it would be extremely challenging to detect
even if we had similarly accurate observations over a full century.

Observations over such a long period would cover much more than a single
planetary orbit, allowing us to approximate that ${\Omega t\gg 1}$. Making
this approximation in Equation 22, we see that the expected displacement
${\propto t/\Omega}$, which is smaller for a less distant planet like Mars.
This justifies our focus on Saturn since it is the most distant planet to
which we have sent an orbiting spacecraft.

Figure 2: The deviation of Saturn from Keplerian motion due to the Sun’s DM
wake, assuming a known initial position and velocity. We show the deviation in
the radial (solid red), azimuthal (dotted blue) and vertical (dot-dashed pink)
directions, with the total shown in black. Observations are assumed to start
when the Sun-Saturn line is maximally aligned with the DM wake (${\phi=0}$ in
Figure 1). The perturbations oscillate with an amplitude that grows almost
linearly (Section 2). Even so, the total deviation is ¡1 cm after a century of
observations.

##  3 Discussion

###  3.1 Higher harmonics of the wake potential

In Section 2, we showed that the ${m=1}$ mode in $\Phi_{w}$ leads to a
resonant perturbation in the limit that $\epsilon\to 0$. The ${m=0}$ mode
leads to a constant radial gravity which causes Saturn to have a smaller time-
averaged $R$ than would otherwise be the case (Equation 26). For a fixed
angular momentum, this increases Saturn’s average angular velocity, causing
${\delta p}$ to grow linearly with time.

Higher frequency modes $\left(m>1\right)$ are also present in Equation 12. We
now investigate the response to these modes in more detail. For this purpose,
Equation 20 is particularly useful since it has a much simpler interpretation
in the non-resonant case.

####  3.1.1 Radial motion

Equation 12 leads to ${m>1}$ modes only because it contains a term of the form
$\left|\cos\phi\right|$. We begin by giving the Fourier representation of
$\Phi_{w}$ for modes $m\geq 2$, thus neglecting modes we have already
considered in Section 2.

$\displaystyle\Phi_{w}~{}=~{}-\frac{4RF_{0}\cos\theta_{0}}{\mathrm{\pi}}\sum_{j=1}^{\infty}\frac{\left(-1\right)^{j}\cos\left(m\phi\right)}{m^{2}-1}\,,\,m\equiv
2j\,.$ (35)

Note that only even $m$ modes are non-zero. Neglecting the $m=0$ mode,
Equation 13 now becomes

$\displaystyle\ddot{\delta
r}+\Omega^{2}r~{}=~{}\frac{12F_{0}\cos\theta_{0}}{\Omega^{2}\mathrm{\pi}}\sum_{j=1}^{\infty}\frac{\left(-1\right)^{j}\cos
m\phi}{m^{2}-1}\,,\,m\equiv 2j\,.$ (36)

Applying the standard solution for a forced harmonic oscillator (Equation 20)
and adding appropriate multiples of the complementary function to satisfy the
boundary condition $\delta{r}=\dot{\delta{r}}=0$ at ${\phi=0}$, the solution
is

$\displaystyle\delta
r~{}=~{}6x_{0}\sum_{j=1}^{\infty}\frac{\left(-1\right)^{j}\left(\cos
m\phi-\cos\phi\right)}{\left(m^{2}-1\right)^{2}}\,,\,m\equiv 2j\,.$ (37)

Comparing this with the dominant $x_{1}$ term in Equation 26, we can now
estimate the fractional error in that equation due to neglecting $m>1$ modes.

$\displaystyle\text{Fractional error in }\delta
r~{}\approx~{}\frac{12\,x_{0}}{\left(m^{2}-1\right)^{2}x_{1}}\approx
0.10\left(m=2\right)\,.$

####  3.1.2 Tangential motion

The generalisation of Equation 27 to ${m>1}$ is

$\displaystyle\dot{\delta p}~{}=~{}-2\Omega
r-\frac{\Phi_{w}}{R\Omega}\left(\cos m\phi-1\right)\,.$ (38)

Substituting in our solution for ${\delta r}$ due to modes with ${m>1}$
(Equation 37) and applying our usual boundary conditions, this becomes

$\displaystyle\dot{\delta p}~{}=~{}-\sum_{j=1}^{\infty}\left[\frac{2\Omega
x_{0}\left(-1\right)^{j}}{\left(m^{2}-1\right)^{2}}\right]$ (39)
$\displaystyle\left[6\left(\cos
m\phi-\cos\phi\right)\right.+\left(m^{2}-1\right)\left(\cos
m\phi-1\right)\left.\right]\,,\,m\equiv 2j\,.$

Requiring $\delta p=0$ fixes the solution to

$\displaystyle\delta
p~{}=~{}-\sum_{j=1}^{\infty}\frac{2\,x_{0}\left(-1\right)^{j}}{\left(m^{2}-1\right)^{2}}$
$\displaystyle\left[\left(m^{2}-5\right)\frac{\sin
m\phi}{m}-6\sin\phi-\left(m^{2}-1\right)\phi\right]\,.$ (40)

Similarly to Section 3.1.1, we can estimate the fractional effect of ${m>1}$
modes on ${\delta p}$ by comparing Equations 28 and 40. In the limit that
$\phi\gg 1$, we get that

$\displaystyle\text{Fractional error in }\delta
p~{}\approx~{}\frac{2x_{0}}{\left(m^{2}-1\right)x_{1}}\approx
0.05\,\left(m=2\right)\,.$

####  3.1.3 Vertical motion

Equation 29 induces ${m>1}$ modes due to the
$\operatorname{sign}\left(\cos\phi\right)$ term, which yields only odd $m$
modes. By finding the strengths of these Fourier modes and neglecting the
previously considered case ${m=1}$, we get that

$\displaystyle\ddot{\delta
z}+\Omega^{2}z~{}=~{}-\frac{4F_{0}\sin\theta_{0}}{\mathrm{\pi}}\sum_{j=1}^{\infty}\frac{\left(-1\right)^{j+1}}{2j-1}\cos\overbrace{\left(2j-1\right)}^{m}\phi\,.$
(41)

Following our derivation in Section 3.1.2, the solution is

$\displaystyle\delta
z~{}=~{}\sum_{j=1}^{\infty}\frac{2\,x_{0}\tan\theta_{0}\left(-1\right)^{j+1}}{m\left(m^{2}-1\right)}\left(\cos
m\phi-\cos\phi\right)\,.$ (42)

Comparing this with Equation 34 shows that it has an error from $m>1$ modes of

$\displaystyle\text{Fractional error in }\delta
z~{}\approx~{}\frac{6\,x_{0}}{m\left(m^{2}-1\right)x_{1}}\approx
0.02\,\left(m=3\right)\,.$

####  3.1.4 Combined effect

Figure 2 indicates that ${\delta p}$ is responsible for the largest part of
the total DM-induced displacement. Our results in Section 3.1.2 indicate that
${m>1}$ modes only affect ${\delta p}$ by ${\approx 5\%}$. Similarly small
effects can be expected for the sub-dominant contributions from ${\delta r}$
(Section 3.1.1) and ${\delta z}$ (Section 3.1.3). Thus, the total DM-induced
displacement of Saturn is only affected a few percent by $m>1$ harmonics in
$\Phi_{w}$. This is much smaller than uncertainties from other parameters like
the ${\approx 10\%}$ error in the local DM density (Hagen & Helmi, 2018).
Therefore, it is currently not very important to consider these higher
harmonics.

###  3.2 Comparison with Hernandez (2019)

Our results in Figure 2 are much smaller than those obtained by Hernandez
(2019). This is because we properly take the resonant limit of Equation 21
while Hernandez (2019) found the maximum amplitude of the induced
perturbations. This maximum is reached after order the Galactic orbital period
of the Sun. For a much shorter observing duration, the maximum possible
amplitude of $\delta r$ is not relevant. The fact that this is only a few
hundred metres (Hernandez, 2019) shows just how little the DM wake would
really influence the Solar System.

The difference between our approach and that of Hernandez (2019) is perhaps
most evident in his equation 6, which states that the perturbation is order
$\widetilde{F}_{1}/\left(\epsilon\Omega\right)$ in our notation. In reality,
our observations span a duration $T\approx 1/\Omega$, so the maximum deviation
that we could see due to the $\widetilde{F}_{1}$ term is
$\widetilde{F}_{1}T^{2}/2\approx\widetilde{F}_{1}/\left(2\Omega^{2}\right)$.
Hernandez (2019) obtained a much larger result because he implicitly assumes
that our observations are long enough to cover all relevant oscillation
periods, including the Galactic orbit of the Sun. This overestimates the
observing duration by $\approx 7$ orders of magnitude, explaining why our
calculated perturbation amplitudes are smaller than his by roughly this
factor.

Nevertheless, we emphasize the great importance of using the Solar System as a
laboratory for testing fundamental physics. Radio tracking data from Cassini
has placed strong constraints on the MOND interpolating function, though some
choices are still consistent with observations (Hees et al., 2014, 2016).
Nearby wide binary stellar systems could be even more important as the MOND
circular velocity exceeds the Newtonian value by ${\approx 20\%}$ (Banik &
Zhao, 2018). Such systems do appear to show a departure from Newtonian
expectations (Hernandez et al., 2019), though it was later shown that the two
can be reconciled with a more careful rejection of outliers (Banik, 2019).

The Chandrasekhar dynamical friction expected from DM remains a promising way
to show its reality, independently of direct and indirect detection
experiments. Although the effects are small in the Solar System, galactic-
scale tests are much more promising since the amount of DM in a system roughly
scales with the cube of its size, more than compensating for the inverse
square law of Newtonian gravity. Results from several galactic tests argue
against the reality of the putative DM halos (Angus et al., 2011; Kroupa,
2015; Oehm et al., 2017). More accurate observations of the Pisces Overdensity
could shed light on this issue if it is indeed a feature in the Galactic
stellar halo caused by the orbit of the Large Magellanic Cloud (Belokurov et
al., 2019).

##  4 Conclusions

The hypothetical DM wake behind the Sun would have effects on planets within
the Solar System (Hernandez, 2019). However, that work overestimates the
effect due to implicitly assuming that our observations cover a duration
comparable to the Galactic orbital period of the Sun. This is most evident in
his equation 6, which shows a larger deviation than the product of force and
observing duration squared. Over the finite ($\la 30$ year) period covered by
accurate space-age observations, the DM wake of the Sun has an imperceptible
effect on the motion of Solar System planets. Thus, the DM hypothesis is not
in tension with presently available Solar System ephemerides. Galactic-scale
tests appear more promising, with the nearby M81 group providing evidence
against the expected dynamical friction (Oehm et al., 2017). In addition, wide
binary systems should soon clarify the true cause of the observed dynamical
discrepancies in galaxies.

## Acknowledgements

IB is supported by an Alexander von Humboldt postdoctoral fellowship. The
graph was prepared using matlab${}^{\text{\textregistered}}$.

## References

  * Angus et al. (2011) Angus G. W., Diaferio A., Kroupa P., 2011, MNRAS, 416, 1401
  * Babcock (1939) Babcock H. W., 1939, Lick Observatory Bulletin, 19, 41
  * Banik (2019) Banik I., 2019, MNRAS, 487, 5291
  * Banik & Zhao (2018) Banik I., Zhao H., 2018, MNRAS, 480, 2660
  * Belokurov et al. (2019) Belokurov V., Deason A. J., Erkal D., Koposov S. E., Carballo-Bello J. A., Smith M. C., Jethwa P., Navarrete C., 2019, MNRAS, 488, L47
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Bozorgnia & Bertone (2017) Bozorgnia N., Bertone G., 2017, International Journal of Modern Physics A, 32, 1730016
  * Chandrasekhar (1943) Chandrasekhar S., 1943, ApJ, 97, 255
  * Famaey & McGaugh (2012) Famaey B., McGaugh S. S., 2012, Living Reviews in Relativity, 15, 10
  * Hagen & Helmi (2018) Hagen J. H. J., Helmi A., 2018, A&A, 615, A99
  * Hees et al. (2014) Hees A., Folkner W. M., Jacobson R. A., Park R. S., 2014, Physical Review D, 89, 102002
  * Hees et al. (2016) Hees A., Famaey B., Angus G. W., Gentile G., 2016, MNRAS, 455, 449
  * Hernandez (2019) Hernandez X., 2019, preprint, Arxiv (arXiv:1901.10605v2)
  * Hernandez et al. (2019) Hernandez X., Cortés R. A. M., Allen C., Scarpa R., 2019, International Journal of Modern Physics D, 28, 1950101
  * Kroupa (2015) Kroupa P., 2015, Canadian Journal of Physics, 93, 169
  * Matson (1992) Matson D. L., 1992, in Kaldeich B., ed., ESA Special Publication Vol. 338, Symposium on Titan.
  * McGaugh (2018) McGaugh S. S., 2018, Research Notes of the American Astronomical Society, 2, 156
  * McGaugh et al. (2016) McGaugh S., Lelli F., Schombert J., 2016, Phys. Rev. Lett., 117, 201101
  * Milgrom (1983) Milgrom M., 1983, ApJ, 270, 365
  * Mulder (1983) Mulder W. A., 1983, A&A, 117, 9
  * Oehm et al. (2017) Oehm W., Thies I., Kroupa P., 2017, MNRAS, 467, 273
  * Ostriker & Peebles (1973) Ostriker J. P., Peebles P. J. E., 1973, ApJ, 186, 467
  * Pitjeva & Pitjev (2013) Pitjeva E. V., Pitjev N. P., 2013, MNRAS, 432, 3431
  * Rogstad & Shostak (1972) Rogstad D. H., Shostak G. S., 1972, ApJ, 176, 315
  * Rubin & Ford (1970) Rubin V. C., Ford Jr. W. K., 1970, ApJ, 159, 379
  * Schönrich et al. (2010) Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829
  * Viswanathan et al. (2017) Viswanathan V., Fienga A., Gastineau M., Laskar J., 2017, NSTIM, 108
  * Xia et al. (2016) Xia Q., et al., 2016, MNRAS, 458, 3839


Now, we assemble all the results of this paragraph and formulate the theorem
about the eigenvalues that lie near zero.

######  Theorem 3.24.

Consider Assumption 3.22. Then for any $\lambda_{0}\in\mathscr{O}$, there
exists a complex neighborhood $U$ of $\lambda_{0}$ such that $\lambda\in U$ is
an eigenvalue of $\mathfrak{D}_{\epsilon}$ if and only if

$n(\lambda,\epsilon)e^{2\xi_{c}(\hat{c})/\epsilon}=1$

where as $\epsilon\downarrow 0$, the function $n(\lambda,\epsilon)$ satisfies

$n(\lambda,\epsilon)=-1+\mathcal{O}(\epsilon).$

###### Proof.

Figure 22 shows the geometric configuration in the complex $x$-plane for our
case, along with the labeling of the corresponding principal regions. In
particular, there are two regions, namely
$\mathscr{D}_{0}=\hat{\mathscr{D}}_{1}$ and
$\mathscr{D}_{1}=\hat{\mathscr{D}}_{0}$, that have both $c$ and $\hat{c}$ on
their boundaries, each being a principal region with respect to either turning
point.

We would like to connect the WKB form of a solution on the segment
$(-\infty,b_{2}]_{\mathscr{P}}$ of a progressive path
$\mathscr{P}\in\mathscr{D}_{2}$, to the the WKB form of the same solution on
the segment $[\hat{b}_{2},+\infty)_{\hat{\mathscr{P}}}$ of a progressive path
$\hat{\mathscr{P}}\in\hat{\mathscr{D}}_{2}$. The essential observation for one
to notice is the relation between $\xi_{c}$ and $\xi_{\hat{c}}$ in the
intermediate area. We have

$\xi_{c}(x)-\xi_{\hat{c}}(x)=\xi_{c}(\hat{c})=-\xi_{\hat{c}}(c),\quad
x\in\mathscr{D}_{0}\cup\mathscr{D}_{1}.$

We take an arbitrary point $a_{0}$ on the common Stokes curve linking $c$ and
$\hat{c}$. Then, $a_{0}$ can be joined to $b_{2}$ by an extension of
$\mathscr{P}$ that passes through $c$ and is progressive; now
$[c,a_{0}]_{\mathscr{P}}$ coincides with the common Stokes curve. By applying
Theorem 3.20 (with $k=0$ and $j=2$) we find that the WKB form of
$w(x,\lambda,\epsilon)$ at $a_{0}$ is given by

$w(x,\lambda,\epsilon)=w_{\rm I,1}^{(c)}(x,\lambda,\epsilon)+w_{\rm
L}^{(c)}(x,\lambda,\epsilon)$

where as $\epsilon\downarrow 0$ we have the asymptotics

$\displaystyle V_{0}^{1/4}(x,\lambda)w_{\rm I,1}^{(c)}(x,\lambda,\epsilon)$
$\displaystyle\doteq
i^{-3}[\lambda_{2,0}+\mathcal{O}(\hat{\chi})]e^{\xi_{c}(x,\lambda)/\epsilon}$
$\displaystyle V_{0}^{1/4}(x,\lambda)w_{\rm L}^{(c)}(x,\lambda,\epsilon)$
$\displaystyle\doteq
i^{-2}[\lambda_{2,1}+\mathcal{O}(\hat{\chi})]e^{-\xi_{c}(x,\lambda)/\epsilon}.$

where $\hat{\chi}(\epsilon)$ is given by (3.18) and (3.16).

To prepare for passage through $\hat{c}$, we observe that
$\mathscr{D}_{0}=\hat{\mathscr{D}}_{1}$. Consequently, we relabel $a_{0}$ now
naming it $\hat{a}_{1}$ and extend $\hat{\mathscr{P}}$, still progressive, to
pass through $\hat{c}$ and continue along the common Stokes curve until
$\hat{a}_{1}$ is reached. In the WKB form of $w(x,\lambda,\epsilon)$ at
$\hat{a}_{1}$ found by passage through $c$, there are two terms: $w_{\rm
I,1}^{(c)}(x,\lambda,\epsilon)$ and $w_{\rm L}^{(c)}(x,\lambda,\epsilon)$. Let
us find how these two behave when they pass through $\hat{c}$.

The contribution from $w_{\rm L}^{(c)}(x,\lambda,\epsilon)$ to the WKB form on
$[\hat{b}_{2},+\infty)_{\hat{\mathscr{P}}}$ is obtained by replacing
$e^{-\xi_{c}(x,\lambda)/\epsilon}$ by
$e^{-\xi_{c}(\hat{c},\lambda)/\epsilon}e^{-\xi_{\hat{c}}(x,\lambda)/\epsilon}$
and applying Theorem 3.20 (with $c$ replaced by $\hat{c}$ and $k=2,j=1$). We
find

$w_{\rm L}^{(c)}(x,\lambda,\epsilon)=w_{\rm
I}^{(\hat{c})}(x,\lambda,\epsilon)$

where as $\epsilon\downarrow 0$ we have

$V_{0}^{1/4}(x,\lambda)w_{\rm I}^{(\hat{c})}(x,\lambda,\epsilon)\doteq
i^{0}[\lambda_{1,2}+\mathcal{O}(\hat{\chi})]e^{-\xi_{c}(\hat{c},\lambda)/\epsilon}e^{\xi_{\hat{c}}(x,\lambda)/\epsilon}.$

To handle the contribution from $w_{\rm I,1}^{(c)}(x,\lambda,\epsilon)$, the
key step is to regard $a_{0}$ as a member of $\hat{\mathscr{D}}_{0}$ and
relabel it as $\hat{a}_{0}$. Since this entails crossing a Stokes curve,
$\xi_{\hat{c}}(x,\lambda)$ is replaced by $-\xi_{\hat{c}}(x,\lambda)$. Thus,
$e^{\xi_{c}(x,\lambda)/\epsilon}$ becomes
$e^{\xi_{c}(\hat{c},\lambda)/\epsilon}e^{-\xi_{\hat{c}}(x,\lambda)/\epsilon}$,
where $\hat{c}$ is regarded as a member of $\mathscr{D}_{0}$ in calculating
$\xi_{c}(\hat{c},\lambda)\in i{\mathbb{R}}^{+}$ and
$x\in\hat{\mathscr{D}}_{0}$. Now, by an application of Theorem 3.20 (with $c$
replaced by $\hat{c}$ and $k=2,j=0$), we find

$w_{\rm I,1}^{(c)}(x,\lambda,\epsilon)=w_{\rm
I}^{(\hat{c})}(x,\lambda,\epsilon)$

where as $\epsilon\downarrow 0$, the following asymptotic form is satisfied

$V_{0}^{1/4}(x,\lambda)w_{\rm
I}^{(\hat{c})}(x,\lambda,\epsilon)\doteq-i^{2}[\lambda_{0,2}+\mathcal{O}(\hat{\chi})]e^{\xi_{c}(\hat{c},\lambda)/\epsilon}e^{\xi_{\hat{c}}(x,\lambda)/\epsilon}.$

Now, combining the last two results, we are capable of obtaining the behavior
of $w(x\lambda,\epsilon)$ on $[\hat{b}_{2},+\infty)_{\hat{\mathscr{P}}}$. As
$\epsilon\downarrow 0$ we find that for
$x\in[\hat{b}_{2},+\infty)_{\hat{\mathscr{P}}}$ we have

$V_{0}^{1/4}(x,\lambda)w(x,\lambda,\epsilon)\doteq\Big{\\{}[1+\mathcal{O}(\hat{\chi})]e^{-\xi_{c}(\hat{c},\lambda)/\epsilon}+[1+\mathcal{O}(\hat{\chi})]e^{\xi_{c}(\hat{c},\lambda)/\epsilon}\Big{\\}}e^{\xi_{\hat{c}}(x,\lambda)/\epsilon}.$

But $\lambda$ is an eigenvalue of $\mathfrak{D}_{\epsilon}$ if and only if
$w(x,\lambda,\epsilon)$ is decaying as $x\rightarrow+\infty$ which is
equivalent with the fact that the quantity inside the braces in the last
result is zero. This easily leads to

$n(\lambda,\epsilon)e^{2\xi_{c}(\hat{c})/\epsilon}=1$

for some function $n(\lambda,\epsilon)$ so that we have

$n(\lambda,\epsilon)=-1+\mathcal{O}(\hat{\chi})\quad\text{as}\quad\epsilon\downarrow
0.$

This completes our proof. ∎

##  4\. Norming Constants

Here, we present the results for the semiclassical behavior of the norming
constants that correspond to the eigenvalues (see sections 2 and 3) of our
Dirac operator. These norming constants, are the proportionality constants
that relate the two eigenfunctions corresponding to a particular eigenvalue of
our operator. We have the following result.

######  Corollary 4.1.

Consider a $\lambda_{0}$ on the asymptotic spectral arcs satisfying the
condition (H1) of §2.4 and for some complex neighborhood $U$ of $\lambda_{0}$
take $\lambda\in U$ to be an eigenvalue of $\mathfrak{D}_{\epsilon}$. Then the
norming constant corresponding to $\lambda$ satisfies the following
asymptotics as $\epsilon\downarrow 0$

$\pm 1+\mathcal{O}(\epsilon).$

###### Proof.

Suppose that ${\bf u}_{0}^{\alpha}$ and ${\bf u}_{0}^{\beta}$ are the
eigenfunctions corresponding to the eigenvalue $\lambda$. These were defined
explicitly in §2.4 as solutions of the problem (2.2). From Lemma 2.17, and the
formulae (2.35), (2.36) (or equivalently Theorem 2.18) we find that

$\frac{{\bf u}_{0}^{\beta}}{{\bf u}_{0}^{\alpha}}=\frac{{\mathcal{W}}[{\bf
u}_{2}^{\beta},{\bf u}_{0}^{\beta}]}{{\mathcal{W}}[{\bf u}_{2}^{\alpha},{\bf
u}_{0}^{\alpha}]}e^{z(\beta,\lambda,\alpha)/\epsilon}=\frac{{\mathcal{W}}[{\bf
u}_{0}^{\beta},{\bf u}_{1}^{\beta}]}{{\mathcal{W}}[{\bf u}_{0}^{\alpha},{\bf
u}_{1}^{\alpha}]}e^{-z(\beta,\lambda,\alpha)/\epsilon}.$

This gives rise to the fact

$\bigg{(}\frac{{\bf u}_{0}^{\beta}}{{\bf
u}_{0}^{\alpha}}\bigg{)}^{2}=\frac{{\mathcal{W}}[{\bf u}_{0}^{\beta},{\bf
u}_{1}^{\beta}]{\mathcal{W}}[{\bf u}_{2}^{\beta},{\bf
u}_{0}^{\beta}]}{{\mathcal{W}}[{\bf u}_{0}^{\alpha},{\bf
u}_{1}^{\alpha}]{\mathcal{W}}[{\bf u}_{2}^{\alpha},{\bf u}_{0}^{\alpha}]}.$

Now using the asymptotic formulas (2.40) and (2.41), we find that $({\bf
u}_{0}^{\beta}/{\bf u}_{0}^{\alpha})^{2}$ behaves like
$-\delta(\alpha_{0},\beta_{0})+\mathcal{O}(\epsilon)$ as $\epsilon\downarrow
0$. Hence as $\epsilon\downarrow 0$, the norming constant satisfies

$\frac{{\bf u}_{0}^{\beta}}{{\bf u}_{0}^{\alpha}}=\pm 1+\mathcal{O}(\epsilon)$

since in all of our cases the turning points in our analysis are zeros of
$g_{-}$ which consequently force $\delta$ to be $-1$ [cf. (2.39)]. ∎

######  Remark 4.2.

The signs of the norming constants change consecutively from one eigenvalue to
the next one as we move along the asymptotic spectral arcs on the spectral
plane.

######  Remark 4.3.

For the norming constants of the eigenvalues that lie near zero, Theorem 3.24
shows that we obtain the same result.

##  5\. Reflection Coefficient

In this section we will consider the reflection coefficient for our Dirac
operator

$\mathfrak{D}_{\epsilon}=\begin{bmatrix}-\frac{\epsilon}{i}\frac{d}{dx}&-iA(x)\exp\\{iS(x)/\epsilon\\}\\\
-iA(x)\exp\\{-iS(x)/\epsilon\\}&\frac{\epsilon}{i}\frac{d}{dx}\end{bmatrix}.$

where again $A(x)=\sech(2x)$ and also the phase function
$S(x)=\sech(2x),x\in{\mathbb{R}}$. Recall that the continuous spectrum of our
Dirac operator is the whole real line. So in this section we are considering
$\lambda\in{\mathbb{R}}$.

###  5.1. Reflection away from zero.

We begin with the case where this $\lambda\in{\mathbb{R}}$ is independent of
$\epsilon$ and consider a $\delta>0$ so that $|\lambda|\geq\delta$. Under the
(different) change of variables

$W_{\pm}=\frac{u_{2}\exp\\{iA/(2\epsilon)\\}\pm
u_{1}\exp\\{-iA/(2\epsilon)\\}}{\sqrt{A\pm i(\lambda+A^{\prime}/2)}}$

equation (2.2) is transformed to the following equation (actually, by applying
the transformation above, we get two independent equations; we only consider
the case for the lower index and set $W=W_{-}$).

$\frac{d^{2}W}{dx^{2}}=[-\epsilon^{-2}f(x,\lambda)+g(x,\lambda)]y$ (5.1)

where $f$ and $g$ are given by the following formulae

$f(x,\lambda)=A^{2}(x)+[\lambda+\tfrac{A^{\prime}(x)}{2}]^{2}$ (5.2)

and

$\displaystyle g(x,\lambda)$
$\displaystyle=\frac{1}{2}\Big{[}\big{(}\log\\{A(x)-i[\lambda+\tfrac{A^{\prime}(x)}{2}]\\}\big{)}^{\prime}\Big{]}^{2}-\frac{\Big{(}\\{A(x)-i[\lambda+\tfrac{A^{\prime}(x)}{2}]\\}^{1/2}\Big{)}^{\prime\prime}}{\\{A(x)-i[\lambda+\tfrac{A^{\prime}(x)}{2}]\\}^{1/2}}$
$\displaystyle=\frac{3}{4}\Bigg{\\{}\frac{A^{\prime}(x)-i\frac{A^{\prime\prime}(x)}{2}}{A(x)-i[\lambda+\tfrac{A^{\prime}(x)}{2}]}\Bigg{\\}}^{2}-\frac{1}{2}\frac{A^{\prime\prime}(x)-i\frac{A^{\prime\prime\prime}(x)}{2}}{A(x)-i[\lambda+\tfrac{A^{\prime}(x)}{2}]}$
(5.3)

We have the following definitions.

######  Definition 5.1.

We define an error-control function $H(x,\lambda)$ for equation (5.1), to be a
primitive of

$\frac{1}{f^{1/4}(x,\lambda)}\frac{d^{2}}{dx^{2}}\Big{[}\frac{1}{f^{1/4}(x,\lambda)}\Big{]}-\frac{g(x,\lambda)}{f^{1/2}(x,\lambda)}.$
(5.4)

######  Definition 5.2.

We define the variation of $H$ in the interval
$(x_{1},x_{2})\subseteq(0,+\infty)$ to be given by

$\mathcal{V}_{x_{1},x_{2}}[H](\lambda)=\int_{x_{1}}^{x_{2}}\bigg{|}\frac{1}{f^{1/4}(t,\lambda)}\frac{d^{2}}{dx^{2}}\Big{[}\frac{1}{f^{1/4}(t,\lambda)}\Big{]}-\frac{g(t,\lambda)}{f^{1/2}(t,\lambda)}\bigg{|}dt.$
(5.5)

Observe that (5.2) implies $f>0$ in ${\mathbb{R}}$. Consequently, equation
(5.1) has no turning points. Furthermore, notice that

  * •

$g$ is complex-valued

  * •

$f$ is twice continuously differentiable with respect to $x$ and

  * •

$g$ is continuous.

These properties allow one (cf. Theorem $2.2$ of §$2.4$ from chapter 6 of [25]
along with the remarks from §5.1 of the same chapter) to arrive at the
following theorem.

######  Theorem 5.3.

Take an arbitrary interval $(x_{1},x_{2})\subseteq{\mathbb{R}}$. Then,
equation (5.1) has in the above interval, two twice continuously
differentiable solutions $w_{\pm}$ with

$w_{\pm}(x,\epsilon)=f^{-1/4}(x,\lambda)\exp\Big{\\{}\pm\frac{i}{\epsilon}\int
f^{1/2}(t,\lambda)dt\Big{\\}}[1+r_{\pm}(x,\epsilon)].$

The remainders $r_{\pm}$ satisfy

$|r_{\pm}(x,\epsilon)|\hskip 5.0pt,\hskip 5.0pt\epsilon
f^{-1/2}(x,\epsilon)\Big{|}\frac{\partial r_{\pm}}{\partial
x}(x,\epsilon)\Big{|}\leq\exp\\{\epsilon\mathcal{V}_{\kappa,x}[H](\lambda)\\}-1$
(5.6)

where $\kappa$ is an arbitrary (finite or infinite) point in the closure of
$(x_{1},x_{2})$, provided that $\mathcal{V}_{\kappa,x}[H](\lambda)<+\infty$.

######  Remark 5.4.

Since $g$ is not real, we cannot expect the solutions $w_{\pm}$ to be complex
conjugates.

######  Remark 5.5.

It follows that

  * •

$r_{\pm}(x,\epsilon)\rightarrow 0$ as $x\to\kappa$ and

  * •

$\epsilon f^{-1/2}(x,\epsilon)\frac{\partial r_{\pm}}{\partial
x}(x,\epsilon)\rightarrow 0$ as $x\to\kappa$.

Using (5.4), notice that $H$ is independent of $\epsilon$ whence the right-
hand side of (5.6) is $\mathcal{O}(\epsilon)$ as $\epsilon\downarrow 0$ and
fixed $x$. But $\mathcal{V}_{x_{1},x_{2}}[H](\lambda)<+\infty$ which implies
that this $\mathcal{O}$-term is uniform with respect to $x$ since
$\mathcal{V}_{\kappa,x}[H](\lambda)\leq\mathcal{V}_{x_{1},x_{2}}[H](\lambda)$.
Hence

$w_{\pm}(x,\epsilon)\sim
f^{-1/4}(x,\lambda)\exp\Big{\\{}\pm\frac{i}{\epsilon}\int
f^{1/2}(t,\lambda)dt\Big{\\}}\quad\text{as}\quad\epsilon\downarrow 0$

uniformly in $(x_{1},x_{2})$.

Next we define the Jost solutions. The Jost solutions are defined as the
components of the bases $\\{J_{-}^{l},J_{+}^{l}\\}$ and
$\\{J_{-}^{r},J_{+}^{r}\\}$ of the two-dimensional linear space of solutions
of equation (5.1), which satisfy the asymptotic conditions

$\displaystyle J_{\pm}^{l}(x,\lambda)$ $\displaystyle\sim\exp\big{\\{}\pm
i\lambda x/\epsilon\big{\\}}\quad\text{as}\quad x\to-\infty$ $\displaystyle
J_{\pm}^{r}(x,\lambda)$ $\displaystyle\sim\exp\big{\\{}\pm i\lambda
x/\epsilon\big{\\}}\quad\text{as}\quad x\to+\infty.$

From scattering theory, we know that the reflection coefficient
$\mathcal{R}(\lambda,\epsilon)$ for the waves incident on the potential from
the right, can be expressed in terms of wronskians of the Jost solutions. More
presicely, we have

$\mathcal{R}(\lambda,\epsilon)=\frac{W[J_{-}^{l},J_{-}^{r}]}{W[J_{+}^{r},J_{-}^{l}]}.$
(5.7)

The next step is to construct the Jost solutions as WKB solutions. For this,
we define the following four WKB solutions

$\displaystyle\bar{w}_{\pm}^{l}(x,\epsilon)=f^{-1/4}(x,\lambda)\exp\Big{\\{}\pm\frac{i}{\epsilon}\Big{(}\lambda
x+\int_{-\infty}^{x}[f^{1/2}(t,\lambda)-\lambda]dt\Big{)}\Big{\\}}[1+\bar{r}_{\pm}^{l}(x,\epsilon)]$
$\displaystyle\bar{w}_{\pm}^{r}(x,\epsilon)=f^{-1/4}(x,\lambda)\exp\Big{\\{}\pm\frac{i}{\epsilon}\Big{(}\lambda
x+\int_{+\infty}^{x}[f^{1/2}(t,\lambda)-\lambda]dt\Big{)}\Big{\\}}[1+\bar{r}_{\pm}^{r}(x,\epsilon)]$

which we are going to modify slightly in a while. If we take the limits as
$x\to\pm\infty$ of the above, we instantly notice the following relations
between $\bar{w}_{\pm}^{l},\bar{w}_{\pm}^{l}$ and the Jost solutions
$J_{\pm}^{l},J_{\pm}^{r}$; we have

$\displaystyle J_{\pm}^{l}$ $\displaystyle=\lambda^{1/2}\bar{w}_{\pm}^{l}$
$\displaystyle J_{\pm}^{r}$ $\displaystyle=\lambda^{1/2}\bar{w}_{\pm}^{r}.$

Let now $w_{\pm}^{l},w_{\pm}^{r}$ be four WKB solutions satisfying

$\displaystyle
w_{\pm}^{l}(x,\epsilon)=f^{-1/4}(x,\lambda)\exp\Big{\\{}\pm\frac{i}{\epsilon}\int_{0}^{x}f^{1/2}(t,\lambda)dt\Big{\\}}[1+r_{\pm}^{l}(x,\epsilon)]$
(5.8) $\displaystyle
w_{\pm}^{r}(x,\epsilon)=f^{-1/4}(x,\lambda)\exp\Big{\\{}\pm\frac{i}{\epsilon}\int_{0}^{x}f^{1/2}(t,\lambda)dt\Big{\\}}[1+r_{\pm}^{r}(x,\epsilon)].$
(5.9)

Once again, the connnection between $w_{\pm}^{l},w_{\pm}^{r}$ and
$\bar{w}_{\pm}^{l},\bar{w}_{\pm}^{r}$ is evident. Indeed, as
$\epsilon\downarrow 0$, we have

$\displaystyle\bar{w}_{\pm}^{l}$
$\displaystyle=\exp\Big{\\{}\pm\frac{i}{\epsilon}\int_{-\infty}^{0}[f^{1/2}(t,\lambda)-\lambda]dt\Big{\\}}w_{\pm}^{l}\big{(}1+\text{o}(1)\big{)}$
$\displaystyle\bar{w}_{\pm}^{r}$
$\displaystyle=\exp\Big{\\{}\mp\frac{i}{\epsilon}\int_{0}^{+\infty}[f^{1/2}(t,\lambda)-\lambda]dt\Big{\\}}w_{\pm}^{r}\big{(}1+\text{o}(1)\big{)}.$

Subsequently, as $\epsilon\downarrow 0$, for the Jost solutions we have

$\displaystyle J_{\pm}^{l}$
$\displaystyle=\lambda^{1/2}\exp\Big{\\{}\pm\frac{i}{\epsilon}\int_{-\infty}^{0}[f^{1/2}(t,\lambda)-\lambda]dt\Big{\\}}w_{\pm}^{l}\big{(}1+\text{o}(1)\big{)}$
(5.10) $\displaystyle J_{\pm}^{r}$
$\displaystyle=\lambda^{1/2}\exp\Big{\\{}\mp\frac{i}{\epsilon}\int_{0}^{+\infty}[f^{1/2}(t,\lambda)-\lambda]dt\Big{\\}}w_{\pm}^{r}\big{(}1+\text{o}(1)\big{)}.$
(5.11)

Recall that the properties of $A$ show that the function $t\mapsto
f^{1/2}(t,\lambda)-\lambda$ is in $L^{1}({\mathbb{R}})$. Furthermore, we have

$\int_{-\infty}^{0}[f^{1/2}(t,\lambda)-\lambda]dt=\int_{0}^{+\infty}[f^{1/2}(t,\lambda)-\lambda]dt=\frac{1}{2}\|f^{1/2}(\cdot,\lambda)-\lambda\|_{L^{1}({\mathbb{R}})}$
(5.12)

and we define

$\sigma(\lambda)=\|f^{1/2}(\cdot,\lambda)-\lambda\|_{L^{1}({\mathbb{R}})}.$
(5.13)

Substituting (5.10), (5.11), (5.12) and (5.13) in (5.7) we have

$\mathcal{R}(\lambda,\epsilon)=e^{i\sigma(\lambda)/\epsilon}\frac{W[w_{-}^{l},w_{-}^{r}]}{W[w_{+}^{r},w_{-}^{l}]}.$
(5.14)

Finally, using (5.8), (5.9) and (5.6) we find that

  * •

$W[w_{-}^{l},w_{-}^{r}]=\mathcal{O}(1)$ as $\epsilon\downarrow 0$ and

  * •

$W[w_{+}^{r},w_{-}^{l}]=-\frac{2i}{\epsilon}+\mathcal{O}(1)$ as
$\epsilon\downarrow 0$.

Substituting these last results in (5.14) we get that

$\mathcal{R}(\lambda,\epsilon)=e^{i\sigma(\lambda)/\epsilon}\mathcal{O}(\epsilon)\quad\text{as}\quad\epsilon\downarrow
0.$

Hence we have just proved the following theorem.

######  Theorem 5.6.

Consider an arbitrary $\delta>0$. Then the reflection coefficient of equation
(5.1) as defined by (5.7), satisfies

$\mathcal{R}(\lambda,\epsilon)=\mathcal{O}(\epsilon)\quad\text{as}\quad\epsilon\downarrow
0$ (5.15)

uniformly for $|\lambda|\geq\delta$.

###  5.2. Reflection near zero.

Now we turn to the case where $\lambda$ depends on $\epsilon$
[$\lambda=\lambda(\epsilon)$] and particularly we let $\lambda$ approach $0$
like $\epsilon^{b}$ for an $\epsilon$-independent positive constant $b$. Using
(5.2) and (5.3), we see that (5.4) can be written as

$\displaystyle\frac{1}{f^{1/4}(x,\lambda)}$
$\displaystyle\frac{d^{2}}{dx^{2}}\Big{[}\frac{1}{f^{1/4}(x,\lambda)}\Big{]}-\frac{g(x,\lambda)}{f^{1/2}(x,\lambda)}=$
$\displaystyle\hskip
10.0pt\frac{5}{16}\\{A(x)^{2}+[\lambda+\tfrac{A^{\prime}(x)}{2}]^{2}\\}^{-5/2}\\{2A(x)A^{\prime}(x)+[\lambda+\tfrac{A^{\prime}(x)}{2}]A^{\prime\prime}(x)\\}^{2}$
$\displaystyle-\frac{1}{4}\\{A(x)^{2}+[\lambda+\tfrac{A^{\prime}(x)}{2}]^{2}\\}^{-3/2}$
$\displaystyle\hskip
61.17325pt\cdot\\{2[A(x)A^{\prime\prime}(x)+A^{\prime}(x)^{2}]+[\lambda+\tfrac{A^{\prime}(x)}{2}]A^{\prime\prime\prime}(x)+\tfrac{A^{\prime\prime}(x)^{2}}{2}\\}$
$\displaystyle-\frac{3}{4}\\{A(x)^{2}+[\lambda+\tfrac{A^{\prime}(x)}{2}]^{2}\\}^{-1/2}\Bigg{\\{}\frac{A^{\prime}(x)-i\frac{A^{\prime\prime}(x)}{2}}{A(x)-i[\lambda+\tfrac{A^{\prime}(x)}{2}]}\Bigg{\\}}^{2}$
$\displaystyle+\frac{1}{2}\\{A(x)^{2}+[\lambda+\tfrac{A^{\prime}(x)}{2}]^{2}\\}^{-1/2}\frac{A^{\prime\prime}(x)-i\frac{A^{\prime\prime\prime}(x)}{2}}{A(x)-i[\lambda+\tfrac{A^{\prime}(x)}{2}]}.$
(5.16)

We can easily see that each of the terms in the sum (5.2) is less than

$C\lambda^{-5/2}e^{-\gamma|x|}$

where $C,\gamma>0$, since $A(x)=\sech(2x)$ and
$A(x),A^{\prime}(x),A^{\prime\prime}(x),A^{\prime\prime\prime}(x)\sim\pm
e^{-2|x|}$ as $x\rightarrow\pm\infty$. Recalling that
$\lambda(\epsilon)\in[\epsilon^{b},+\infty)$ where $b>0$ is independent of
$\epsilon$, and (5.6) we see that the variation in (5.5) behaves like

$\mathcal{V}_{0,+\infty}[H](\lambda(\epsilon))=\mathcal{O}\Big{(}\epsilon^{-5b/2}\Big{)}\quad\text{as}\quad\epsilon\downarrow
0.$

Hence using (5.8), (5.9) we get

  * •

$W[w_{-}^{l},w_{-}^{r}]=\mathcal{O}(\epsilon^{-5b/2})$ as $\epsilon\downarrow
0$ and

  * •

$W[w_{+}^{r},w_{-}^{l}]=-\frac{2i}{\epsilon}+\mathcal{O}(\epsilon^{-1-5b/2})$
as $\epsilon\downarrow 0$.

Substituting these last results in (5.14), we finally obtain that

$\mathcal{R}(\lambda(\epsilon),\epsilon)=e^{i\sigma(\lambda(\epsilon))/\epsilon}\mathcal{O}(\epsilon^{1-5b/2})\quad\text{as}\quad\epsilon\downarrow
0.$ (5.17)

So, we have showed the following.

######  Theorem 5.7.

Consider $0<b<\tfrac{2}{5}$ (independent of $\epsilon$). Then the reflection
coefficient of equation (5.1) as defined by (5.7), satisfies

$\mathcal{R}(\lambda(\epsilon),\epsilon)=\mathcal{O}(\epsilon^{1-5b/2})\quad\text{as}\quad\epsilon\downarrow
0$ (5.18)

uniformly for $\lambda(\epsilon)$ in any closed interval of
$[\epsilon^{b},+\infty)$.

######  Remark 5.8.

A similar (and better) estimate can be achieved by using the exact WKB method.
See Corollary 2.8 in [14].

##  6\. Conclusion: Application to the semiclassical focusing NLS

In a sense what we have been able to prove in the previous sections is that
the behavior of the solution of (1.1) in the regime $\epsilon\downarrow 0$ can
be approximated by a “soliton ensemble”, i.e a purely soliton solution. This
solution is uniquely defined by the eigenvalues of the associated Dirac
operator and their corresponding norming constants: the eigenvalues are
exactly the WKB approximations of the previous sections while the norming
constants alternate between -1 and 1. The role of the continuous part of the
spectrum is negligible when epsilon is small, since the reflection coefficient
is small away from zero 888We refer to [16] for the actual details justifying
the approximation.. The remaining question is the semiclassical behavior of
the soliton ensemble itself.

This is a generalization of the problem considered in [19] and [20], the
difference being that here the relevant WKB-eigenvalues lie in the Y-shaped
set of Figure 10 (together with its reflection with respect to the real axis),
while in [19] the spectrum is purely imaginary (the four non-imaginary
“branches” are not there). Even though the proofs of [19] and [20] are long
and complicated they can be readily transferred from one case to another.

We simply recall that the starting point of the strategy in [19] consists of
considering a fairly arbitrary contour around the imaginary discrete spectrum
and transforming the meromorphic Riemann-Hilbert problem (with no jump
contour) given by the inverse scattering theory to a holomorphic Riemann-
Hilbert problem with jump given on that fairly arbitrary contour. $Eventually$
this contour has to be deformed to an optimal contour on which the asymptotic
analysis of the holomorphic Riemann-Hilbert problem will be feasible.

The limiting density of eigenvalues is a crucial quantity in the analysis. In
the special case $A(x)=S(x)=\sech(2x)$, for example,

$\rho(\lambda)=\frac{\lambda}{\pi}\int_{x_{1}(\lambda)}^{x_{2}(\lambda)}\frac{dx}{(A(x)^{2}+(\lambda+\frac{1}{2}S^{\prime}(x))^{2})^{1/2}}$
(6.1)

for $\lambda\in\tilde{\Lambda}_{12}$, while

$\rho(\lambda)=\frac{\lambda}{\pi}\int_{x_{1}(\lambda)}^{x_{6}(\lambda)}\frac{dx}{(A(x)^{2}+(\lambda+\frac{1}{2}S^{\prime}(x))^{2})^{1/2}}$
(6.2)

for $\lambda\in\tilde{\Lambda}_{16}$ and

$\rho(\lambda)=\frac{\lambda}{\pi}\int_{x_{2}(\lambda)}^{x_{6}(\lambda)}\frac{dx}{(A(x)^{2}+(\lambda+\frac{1}{2}S^{\prime}(x))^{2})^{1/2}}$
(6.3)

for $\lambda\in\tilde{\Lambda}_{26}$. By symmetry, the last two integrals are
equal and it is not hard to see that all three integrals are analytic in the
upper half-plane, the top integral one being equal to the sum of the two
others999Just write every integral as a contour integral on a closed curve
wrapping around each $\tilde{\Lambda}_{ij}$ and then deform to a curve
independent of $\lambda$ to get local analyticity. Then, as $\lambda$ varies,
start deforming those closed curves along the moving cuts (actually Stokes
curves moving with $\lambda$) to achieve global analyticity. The geometry of
turning points is known numerically and this is crucial. One can deform such
that the number of turning points along each deformation remains two and cover
the whole upper half-plane. The only possibly problematic point is the
bifurcation point $\lambda_{\otimes}$ where a third turning point appears. But
the linear relation between the three integrals above shows that analyticity
holds at that point as well.. So $\rho$ is not continuous at the bifurcation
point $\lambda_{\otimes}$ if defined as a function on the Y-shaped set of
Figure 10 but it is a sum of two functions $\rho_{1}$ and $\rho_{2}$ such that
$\rho_{1}$ is continuous in $\tilde{\Lambda}_{12}\cup\tilde{\Lambda}_{16}$ and
$\rho_{2}$ is continuous in $\tilde{\Lambda}_{12}\cup\tilde{\Lambda}_{26}$. In
fact, each $\rho_{j}$ admits (the same) analytic extension, which is of course
$\frac{1}{2}\rho$. We also note that because of the nice behavior of $A$ and
$S$ at infinity, $\rho$ is regular at $0$.

It follows eventually that the strategy of [19] and [20] goes through by
performing the first step of [19] for a sum of two measures, one supported on
$\tilde{\Lambda}_{12}\cup\tilde{\Lambda}_{16}$ and the other on
$\tilde{\Lambda}_{12}\cup\tilde{\Lambda}_{26}$ each with an analytic density
with an entire extension. Even though the Y-shaped set of Figure 10 is more
complicated than the imaginary interval appearing as an asymptotic spectrum in
the zero phase case, the symmetries involved (with respect to conjugation and
taking negatives) ensure the properties of the density $\rho$ needed in the
analysis of [20]. We leave the details to the interested reader 101010Note
that the imaginary part of $\rho$ is positive on the real line. This is the
main condition in [20] for the admissibility of the contour maximizing the
equilibrium energy. .

##  Appendix A Modified Bessel Functions

###  A.1. Primary solutions

Let us in this section of the appendix describe some properties of Modified
Bessel Functions (or MBFs) which are used in certain parts of this survey. For
the proofs and a more detailed analysis on the material presented in this
section, we refer the interested reader to [26] and §2 of [24].

We start by investigating solutions of the following equation

$\frac{d^{2}w(t)}{dt^{2}}=\frac{9}{4}tw(t).$ (A.1)

When $t$ is real we shall adopt the function

$U(t)=\Big{(}\frac{2t}{\pi}\Big{)}^{1/2}K_{1/3}(t^{3/2})$ (A.2)

as a standard solution, where $K_{1/3}$ denotes the MBF of order $1/3$ in the
usual notation. We continue to employ this solution in the more general case
of $t$ being a complex number. The functions on the RHS of (A.2) are
understood to assume their principal values when $\mathrm{ph}t=0$ and be
defined by continuity for other values of $\mathrm{ph}t$.

###  A.2. Secondary solutions

Let us fix some notation and define some sectors on the complex plane. So, for
$j\in\\{0,1,2\\}$ we set

$S_{j}=\\{\,x\in{\mathbb{C}}\mid(2j-1)\tfrac{\pi}{3}\leq\mathrm{ph}x\leq(2j+1)\tfrac{\pi}{3}\,\\}.$

If there exists a number $0<\delta<1$ we say that the set

${\rm
int}_{\delta}(S_{j})=\\{\,x\in{\mathbb{C}}\mid(2j-\delta)\tfrac{\pi}{3}\leq\mathrm{ph}x\leq(2j+\delta)\tfrac{\pi}{3}\,\\}$

is an internal part of the sector $S_{j}$.

We now introduce the secondary solutions for equation (A.1). Here, we treat
$t$ as a complex number. These are defined by

$U_{j}(t)=U(te^{-i2\pi j/3}),\quad j\in\\{0,1,2\\}.$

It can be easily observed that $U_{j}$ is a solution of equation (A.1) that is
decaying in the sector $S_{j}$, as $t\rightarrow\infty$.

For $j,k\in\\{0,1,2\\}$, the wronskian of the pair $\\{U_{j},U_{k}\\}$ is
given by

$\displaystyle W[U_{j},U_{k}]$
$\displaystyle=3ie^{-i(j+k)\pi/3}\lambda_{j,k}\quad\text{where}$
$\displaystyle\lambda_{j,k}$
$\displaystyle=\frac{\sin[(k-j)\pi/3]}{\sin\pi/3}$

From this, we conclude that $U_{j}$, $U_{k}$ are linearly independent when
$j,k\in\\{0,1,2\\}$ with $j\neq k$.

###  A.3. Auxiliary functions

Note: throughout this section, we consider $j,k\in\\{0,1,2\\}$ with $j\neq k$
(so that $\lambda_{j,k}\neq 0$).

First of all, it is well known that as $t\rightarrow\infty$, the solution
$U_{j}$ is decaying in $S_{j}$ and dominant in $S_{k}$ and by symmetry,
$U_{k}$ is dominant in $S_{j}$ and decaying in $S_{k}$. This in turn, implies
that $\\{U_{j},U_{k}\\}$ comprises a numerically satisfactory pair of
solutions in the closed region $S_{j}\cup S_{k}$ except possibly in a
neighborhood of $0$. Hence, the pair above constitutes an appropriate solution
basis in $S_{j}\cup S_{k}$ and in order to majorize these solutions, we shall
introduce auxiliary weight, modulus and phase functions.

We start by defining the function $E(t)$ by the formula

$E(t)=|\exp\\{(-1)^{j}t^{3/2}\\}|\quad\text{for}\quad t\in S_{j}$

where the branch of $t^{3/2}$ is $|t|^{3/2}\exp\\{i3(\mathrm{ph}t)/2\\}$.On
the boundaries of $S_{j}$, $E(t)=1$ and so $E$ is continuous everywhere with
$E(t)\geq 1$. As weight functions in $S_{j}\cup S_{k}$ we assign the functions
$\mathsf{E}_{j,k}(t)$ defined by

$\mathsf{E}_{j,k}(t)=\begin{cases}1/E(t),\quad t\in S_{j}\\\ E(t),\quad t\in
S_{k}\end{cases}$

It is a continuous function with $\mathsf{E}_{j,k}(0)=1$. Also it satisfies
the following symmetry

$\mathsf{E}_{k,j}(t)=1/\mathsf{E}_{j,k}(t)\quad\text{for}\quad t\in S_{j}\cup
S_{k}.$

For the modulus and phase functions in $S_{j}\cup S_{k}$ we define

$\displaystyle|U_{j}(t)|$
$\displaystyle=\mathsf{E}_{j,k}(t)\mathsf{M}_{j,k}(t)\cos\mathsf{\theta}_{j,k}(t)$
$\displaystyle|U_{j}^{\prime}(t)|$
$\displaystyle=\mathsf{E}_{j,k}(t)\mathsf{N}_{j,k}(t)\cos\mathsf{\omega}_{j,k}(t)$

and

$\displaystyle|U_{k}(t)|$
$\displaystyle=\frac{1}{\mathsf{E}_{j,k}(t)}\mathsf{M}_{j,k}(t)\sin\mathsf{\theta}_{j,k}(t)$
$\displaystyle|U_{k}^{\prime}(t)|$
$\displaystyle=\frac{1}{\mathsf{E}_{j,k}(t)}\mathsf{N}_{j,k}(t)\sin\mathsf{\omega}_{j,k}(t)$

From these, we arrive at

$\displaystyle\mathsf{M}_{j,k}(t)$
$\displaystyle=\bigg{[}\frac{|U_{j}(t)|^{2}}{\mathsf{E}_{j,k}(t)^{2}}+\mathsf{E}_{j,k}(t)^{2}|U_{k}(t)|^{2}\bigg{]}^{1/2}$
$\displaystyle\mathsf{N}_{j,k}(t)$
$\displaystyle=\bigg{[}\frac{|U_{j}^{\prime}(t)|^{2}}{\mathsf{E}_{j,k}(t)^{2}}+\mathsf{E}_{j,k}(t)^{2}|U_{k}^{\prime}(t)|^{2}\bigg{]}^{1/2}$
$\displaystyle\mathsf{\theta}_{j,k}(t)$
$\displaystyle=\arctan\bigg{[}\mathsf{E}_{j,k}(t)^{2}\bigg{|}\frac{U_{k}(t)}{U_{j}(t)}\bigg{|}\bigg{]}$
$\displaystyle\mathsf{\omega}_{j,k}(t)$
$\displaystyle=\arctan\bigg{[}\mathsf{E}_{j,k}(t)^{2}\bigg{|}\frac{U_{k}^{\prime}(t)}{U_{j}^{\prime}(t)}\bigg{|}\bigg{]}$

The above formulae show that each one of
$\mathsf{M}_{j,k}(t),\mathsf{N}_{j,k}(t),\mathsf{\theta}_{j,k}(t)$ and
$\mathsf{\omega}_{j,k}(t)$ is continuous in $S_{j}\cup S_{k}$. Also we have
the following symmetry relations

$\displaystyle\mathsf{M}_{j,k}(t)$ $\displaystyle=\mathsf{M}_{k,j}(t)$
$\displaystyle\mathsf{N}_{j,k}(t)$ $\displaystyle=\mathsf{N}_{k,j}(t)$
$\displaystyle\mathsf{\theta}_{j,k}(t)+$
$\displaystyle\mathsf{\theta}_{k,j}(t)=\frac{\pi}{2}$
$\displaystyle\mathsf{\omega}_{j,k}(t)+$
$\displaystyle\mathsf{\omega}_{k,j}(t)=\frac{\pi}{2}$

In internal parts of $S_{j}$, $S_{k}$, for large $|t|$ the asymptotic behavior
of the modulus functions is given by

$\displaystyle\mathsf{M}_{j,k}(t)$
$\displaystyle\sim(1+\lambda_{j,k}^{2})^{1/2}|t|^{-1/4}$
$\displaystyle\mathsf{N}_{j,k}(t)$
$\displaystyle\sim\frac{3}{2}(1+\lambda_{j,k}^{2})^{1/2}|t|^{1/4}.$

But in the full domain $S_{j}\cup S_{k}$, the modulus and phase functions
fluctuate as $t\rightarrow\infty$. For $\mathsf{M}_{j,k}(t)$ we have the
following bound

$\mathsf{M}_{j,k}(t)\leq
C_{j,k}|t|^{-1/4}[1+\Theta(t^{3/2})]\quad\text{for}\quad t\in S_{j}\cup S_{k}$

where $\Theta(t)$ and $C_{j,k}$ are given by

$\displaystyle\Theta(t)$ $\displaystyle=\exp\\{5\pi/(72|t|)\\}-1$ (A.3)
$\displaystyle C_{j,k}$
$\displaystyle=\max\\{[1+(|\lambda_{j,k}|+|\lambda_{j,k+1}|)^{2}]^{1/2},[1+(|\lambda_{j,k}|+|\lambda_{j,k-1}|)^{2}]^{1/2}\\}.$
(A.4)

Unfortunately, no such simple bound is available for $\mathsf{N}_{j,k}(t)$ and
for this we shall introduce another pair of modulus and phase functions that
will be convenient to work with. This is defined by

$\displaystyle\bigg{|}\frac{d[t^{1/4}U_{j}(t)]}{dt}\bigg{|}$
$\displaystyle=|t|^{1/4}\mathsf{E}_{j,k}(t)\hat{\mathsf{N}}_{j,k}(t)\cos\hat{\mathsf{\omega}}_{j,k}(t)$
$\displaystyle\bigg{|}\frac{d[t^{1/4}U_{k}(t)]}{dt}\bigg{|}$
$\displaystyle=|t|^{1/4}\frac{1}{\mathsf{E}_{j,k}(t)}\hat{\mathsf{N}}_{j,k}(t)\sin\hat{\mathsf{\omega}}_{j,k}(t)$

Immediately, we find that

$\displaystyle\hat{\mathsf{N}}_{j,k}(t)$
$\displaystyle=|t|^{-1/4}\bigg{\\{}\frac{1}{\mathsf{E}_{j,k}(t)^{2}}\bigg{|}\frac{d[t^{1/4}U_{j}(t)]}{dt}\bigg{|}^{2}+\mathsf{E}_{j,k}(t)^{2}\bigg{|}\frac{d[t^{1/4}U_{k}(t)]}{dt}\bigg{|}^{2}\bigg{\\}}^{1/2}$
$\displaystyle\hat{\mathsf{\omega}}_{j,k}(t)$
$\displaystyle=\arctan\bigg{\\{}\mathsf{E}_{j,k}(t)^{2}\bigg{|}\frac{d[t^{1/4}U_{k}(t)]}{dt}\bigg{|}\Big{/}\bigg{|}\frac{d[t^{1/4}U_{j}(t)]}{dt}\bigg{|}\bigg{\\}}.$

Symmetry now tell us that

$\displaystyle\hat{\mathsf{N}}_{j,k}(t)$
$\displaystyle=\hat{\mathsf{N}}_{k,j}(t)$
$\displaystyle\hat{\mathsf{\omega}}_{j,k}(t)+$
$\displaystyle\hat{\mathsf{\omega}}_{k,j}(t)=\frac{\pi}{2}$

while as $t\rightarrow\infty$ in internal parts of $S_{j}$ or $S_{k}$ we have

$\hat{\mathsf{N}}_{j,k}(t)\sim\frac{3}{2}(1+\lambda_{j,k}^{2})^{1/2}|t|^{1/4}.$

Lastly, $\hat{\mathsf{N}}_{j,k}(t)$ satisfies the following bound

$\hat{\mathsf{N}}_{j,k}(t)\leq\frac{3}{2}C_{j,k}|t|^{1/4}[1+\Theta(t^{3/2})]\quad\text{for}\quad
t\in S_{j}\cup S_{k}.$

where $\Theta(t)$ and $C_{j,k}$ are given by (A.3) and (A.4) respectively.

## Data Availability

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## Acknowledgements

We would like to thank Nikos Efremidis for Figures 1 and 2 in the
Introduction. Figures 9 through 21 have been taken from [22]. Also

  * •

The research of the first author was supported by a JSPS grant-in-aid No.
21K03303.

  * •

The second author acknowledges the support of the Institute of Applied and
Computational Mathematics of the Foundation of Research and Technology -
Hellas (FORTH), via grant MIS 5002358 and the support of the University of
Crete via grant 10753. He also expresses his sincere gratitude to the
Independent Power Transmission Operator(IPTO) for a scholarship through the
School of Sciences and Engineering of the University of Crete. Finally he
acknowledges the support of the Fulbright Institute in Greece through a
Visiting Researcher Fellowship at the Department of Mathematics, Princeton
University.

  * •

The third author acknowledges the generous support of Ritsumeikan University
during a visit in $2018$ and the support of the Greek Secretariat of Research
and Technology under ARISTEIA II Grant No. 3964.

## References

  * [1] J. C. Bronski: Semiclassical eigenvalue distribution of the Zakharov–Shabat eigenvalue problem, Physica D 97, pp. 376–397 (1996)
  * [2] J. C. Bronski: Spectral instability of the semiclassical Zakharov–Shabat eigenvalue problem, Physica D 152-153, pp. 163–170 (2001)
  * [3] T. Brooke Benjamin, J.E. Feir: The disintegration of wave trains on deep water. Part 1. Theory, Journal of Fluid Mechanics, 27 (3), pp. 417–430 (1967)
  * [4] G. Biondini, D. Mantzavinos: Universal Nature of the Nonlinear Stage of Modulational Instability, Phys. Rev. Lett., 116, 043902 (2016)
  * [5] M. Bertola, A. Tovbis: Universality for the focusing nonlinear Schrödinger equation at the gradient catastrophe point: Rational breathers and poles of the tritronquée solution to Painlevé I, Communications in Pure and Applied Mathematics, 66, n.5, pp. 678–752 (2013)
  * [6] P. Deift, X. Zhou: A Steepest Descent Method for Oscillatory Riemann-Hilbert Problems; Asymptotics for the MKdV Equation, Annals of Mathematics, Second Series, 137, n.2, pp. 295-368 (1993)
  * [7] P. Deift, S. Venakides, X. Zhou: New Results in Small Dispersion KdV by an Extension of the Steepest Descent Method for Riemann-Hilbert Problems, Inter.Math.Res.Notices, 1997, pp.286-299.
  * [8] N. Dencker, The pseudospectrum of systems of semiclassical operators, Analysis & PDE, 1(3) (2008), pp.323-373.
  * [9] J. Ecalle: Cinq applications des fonctions résurgentes, prépublication Orsay 1981.
  * [10] S. Fujiié, C. Lasser, L. Nédélec: Semiclassical resonances for a two-level Schrödinger operator with a conical intersection, Asymptotic Analysis, 65, n.1-2 (2009), pp. 17-58.
  * [11] S. Fujiié, T. Ramond: Matrice de scattering et résonances asociées à une orbite hétérocline, Ann. Inst. H. Poincaré Phys. Théor., 69, n.1 (1998), pp. 31-82.
  * [12] L. Faddeev, L. Takhtajan: Hamiltonian Methods in the Theory of Solitons, Springer 1987.
  * [13] C. Gérard, A. Grigis: Precise Estimates of Tunneling and Eigenvalues near a Potential Barrier, J.Differential Equations, 72 (1988), pp.149-177.
  * [14] S. Fujiié and S. Kamvissis, Semiclassical WKB Problem for the Non-self-adjoint Dirac Operator with Analytic Potential, Journal of Mathematical Physics 61, no. 1, 011510 (2020).
  * [15] N. Hatzizisis and S. Kamvissis, Semiclassical WKB Problem for the Non-Self-Adjoint Dirac Operator with a Decaying Potential, Journal of Mathematical Physics 62, no. 3, 033510 (2021).
  * [16] N. Hatzizisis and S. Kamvissis, Semiclassical WKB Problem for the Non-Self-Adjoint Dirac Operator with a Multi-Humped Decaying Potential, arXiv 2106.07253.
  * [17] K. Hirota and J. Wittsten, Complex eigenvalue splitting for the Dirac operator, Communications in Mathematical Physics, v. 383, pp. 1527-1558 (2021).
  * [18] N. Korneev, J.A. Catana Castellanos, V.A. Vysloukh: Multisoliton pulse breakup in WKB approximation, Optik 207 (2020) 164359.
  * [19] S. Kamvissis, Kenneth D. T.-R. McLaughlin, P. D. Miller : Semiclassical Soliton Ensembles for the Focusing Nonlinear Schrödinger Equation, Annals of Mathematics, 154 (2003), Princeton University Press, Princeton, NJ.
  * [20] S. Kamvissis, E. A. Rakhmanov: Existence and Regularity for an Energy Maximization Problem in Two Dimensions, Journal of Mathematical Physics, 46, n.8 (2005); addendum in Journal of Mathematical Physics, 50, n.9 (2009).
  * [21] E. Madelung : Quantentheorie in Hydrodynamischer Form,  Z. Phys., 40 (3-4), pp. 322 - 326 (1927)
  * [22] P. D. Miller, Some Remarks on a WKB Method for the Non-selfadjoint Zakharov-Shabat Eigenvalue Problem with Analytic Potentials and Fast Phase, Physica D, Nonlinear Phenomena 152 (2001), pp. 145-162.
  * [23] F. W. J. Olver, Second-Order Linear Differential Equations with Two Turning Points, Philosophical Transactions of the Royal Society of London, Series A, Mathematical and Physical Sciences 278, no. 1279 (1975), pp. 137-174.
  * [24] F. W. J. Olver, General Connection Formulae for Liouville-Green Approximations in the Complex Plane, Philosophical Transactions of the Royal Society of London, Series A, Mathematical and Physical Sciences, Vol. 289, No. 1364, pp. 501-548, (1978).
  * [25] F. W. J. Olver, Asymptotics and Special Functions, AK Peters/CRC Press (1997).
  * [26] F. W. J Olver, D. W. Lozier, R. F. Boisvert and C. W. Clark, eds. NIST Handbook of Mathematical Functions (Hardback and CD-ROM), Cambridge University Press (2010).
  * [27] A. Voros : The return of the quartic oscillator. The complex W.K.B. method, Ann. Inst. H. Poincaré, 29, pp. 211-338 (1983)
  * [28] V. E. Zakharov and A. B. Shabat : Exact Theory of Two-dimensional Self-focusing and One-dimensional Self-modulation of Wave in Nonlinear Media, Journal of Experimental and Theoretical Physics 34 n.1, pp. 62-69 (1972).