# Multiple Ultralight Axionic Wave Dark Matter and Astronomical Structures

Hoang Nhan Luu<EMAIL_ADDRESS>Institute for Advanced Study and
Department of Physics, Hong Kong University of Science and Technology, Hong
Kong  S.-H. Henry Tye<EMAIL_ADDRESS>Institute for Advanced Study and
Department of Physics, Hong Kong University of Science and Technology, Hong
Kong  Department of Physics, Cornell University, Ithaca, NY 14853, USA  Tom
Broadhurst<EMAIL_ADDRESS>Department of Theoretical Physics,
University of the Basque Country UPV/EHU, E-48080 Bilbao, Spain  Donostia
International Physics Center (DIPC), 20018 Donostia-San Sebastian (Gipuzkoa)
Spain.  Ikerbasque, Basque Foundation for Science, E-48011 Bilbao, Spain

###### Abstract

An ultralight scalar boson with mass $m_{1}\simeq 10^{-22}$ eV is gaining
credence as a Dark Matter (DM) candidate that explains the dark cores of dwarf
galaxies as soliton waves. Such a boson is naturally interpreted as an axion
generic in String Theory, with multiple light axions predicted in this
context. We examine the possibility of soliton structures over a wide range of
scales, accounting for galaxy core masses and the common presence of nuclear
star clusters. We present a diagnostic soliton core mass-radius plot that
provides a global view, indicating the existence of an additional axion with
mass $m_{2}\simeq 10^{-20}$ eV, with the possibility of a third axion with
mass $m_{3}\gtrsim 0.5\times 10^{-18}$ eV. We also argue that the relative
mass densities measured for these axions are consistent with their
cosmological production via the mis-alignment mechanism.

##  I Introduction

Dark Matter (DM) is well established from a wide set of astronomical evidence,
including dynamical, lensing and CMB data. However, its nature is far from
clear, requiring new physics beyond the known standard particle physics that
describes only $\simeq 17\%$ of the cosmological mass density Ade2015xua ;
Cyburt2015mya . It is understood that dark matter must be predominantly non-
relativistic, to the earliest limits of observation, otherwise large-scale
structure would be featureless on small scales. However, no evidence for
heavy, non-relativistic weakly interacting massive particles (WIMPs) has been
found, despite stringent laboratory searches. Alternatively, dark matter as a
Bose-Einstein condensate also satisfies the non-relativistic requirement, as
light bosons in the ground state behave as a condensate with high occupation
number. The uncertainty principle means bosons cannot be confined within the
de Broglie scale, thus avoiding the density cusp formation problem for WIMPs,
as raised in the fuzzy dark matter scenario Hu2000ke  with an ultralight
boson.

Ultralight bosonic dark matter has recently been successfully simulated for
the first time, revealing an unforeseen rich wavelike substructure that may be
termed “wave dark matter” ($\psi$DM) Schive2014a ; Schive2014b ; Schive2016 ;
Hui2016ltb ; Niemeyer2019aqm . A prominent solitonic wave at the base of every
virialised potential is predicted by $\psi$DM, representing the ground state,
where self gravity is matched by effective pressure from the uncertainty
principle. The solitons found in the simulations have flat cored profiles that
accurately match the known time independent solution of the Schrödinger-
Poisson equation Schive2014a ; Schive2014b ; Mocz2017wlg ; Veltmaat2018 , for
which the soliton mass scales inversely with its radius Guzman2006 .

Here we take seriously the generic String Theory prediction of multiple axion
fields on astronomical scales. An axion mass has been derived in this context
to be $m_{1}\simeq 10^{-22}$ eV Hu2000ke ; Schive2014a ; DeMartino2018zkx .
Here, we will advocate a solitonic origin for the puzzling dynamically
distinct, Nuclear Star Cluster (NSC) of $10^{7}\text{ M}_{\odot}$ in the Milky
Way that surrounds the central black hole on a scale of $\simeq 1$ pc. The
inner density profile of this NSC is fitted by a dense soliton of dark matter
corresponding to a heavier axion with mass $\simeq 10^{-20}$ eV. This inner
soliton amounts to a small dark matter contribution in addition to the
dominant soliton due to the $10^{-22}$ eV axion responsible for galaxy
formation in this context, forming this “soliton in soliton” structure within
the Milky Way on scales of 100 pc and 1 pc for axions of $10^{-22}$ eV and
$10^{-20}$ eV respectively.

Figure 1: The Soliton Core Mass-Radius Plot illustrates the map of soliton
core mass $M_{c}$ versus the core radius $r_{c}$ from cosmological structures
at hierarchically different scales of the observational data points plotted
together with theoretical prediction of solitons. The Milky Way (MW) and dSph
galaxies data suggest an axion mass $m_{1}\simeq 10^{-22}$ eV. The UFDs
require another axion with $m_{2}\simeq\times 10^{-20}$ eV and it is possible
that NSCs are consistent with this more massive axion as indicated here
tentatively for best studied examples, including the NSC in the Milky Way and
NGC205. The dark core of the well studied globular cluster Tuc-47 data has
also been suggested to harbour a soliton and the mass and scale of this
corresponds to and axion of approximately $m_{3}\gtrsim 0.5\times 10^{-18}$ eV
as shown lower left. Filled-colored bands are approximate mass ranges for the
three possible axions masses that may be inferred as above for the current
best data, spanning a wide range of astronomical objects that appear to
dominated or to contain significant dark matter in excess of the visible
stellar mass. The squares indicate objects for which axion mass constraints
have been obtained previously and circles those obtained here. The triangles
show predicted scale sizes of solitons in order that they lie in given axion
mass range. The Schwarzschild and the Jeans bounds are shown for reference.
Dashed error bars indicate preliminary contraints for which substantial
revision may be expected with a better determination of the scale of the
associated dark matter.

The available data suggests that the primordial density of the lighter axion
($m_{1}$) is substantially bigger than that of the heavier ($m_{2}$) axion,
even if the opposite should be true if they have equal number densities.
Fortunately, this is consistent with cosmological considerations where the
lighter scalar field is likely to have a higher density, if they are axions;
that is, they are massless (due to shift symmetries) until non-perturbative
dynamics generate the tiny masses.

In this paper, we first consider the solitonic properties of the multiple
axion scenario, with a solution for “nested” solitons corresponding to axions
of different masses. We find that, to a reasonable approximation, solitons
from different axions have negligible effect on each other. Based on this
fact, we then form a global view by displaying all claimed relevant data in a
plane of soliton core mass versus core radius, in Fig. 1; as shown, for a
given axion mass $m_{i}$, there is a scaling solitonic solution, where the
soliton core mass scales inverse to the core radius. The main result of this
paper, suggesting the existence of at least 2 distinct ultralight axions, with
masses $m_{1}\simeq 10^{-22}$ eV and $m_{2}\simeq 10^{-20}$ eV, and the
possibility of a third axion, $m_{3}\gtrsim 0.5\times 10^{-18}$ eV. In
addition to the data that we interpret as soliton-in-soliton phenomena, there
are also data indicating the presence of single $m_{2}$ solitons. Together,
they provide support for the existence of $m_{2}$ in addition to $m_{1}$.
Further focused testing of this unifying conclusion will clarify the extent to
which these widely different astronomical structures can be understood as a
manifestation of multi-ultralight axions. The data is also consistent with the
axionic cosmology where the lighter axion is expected to have a higher
primordial density. We conclude with predictions for new observations.

##  II Multiple Axions and Nested Solitons

It is straightforward to extend the $\psi$DM formalism to that for the
multiple axion case. At the first order perturbative theory and in the non-
relativistic limit, we have the Schr$\ddot{\text{o}}$dinger-Poisson equations
for N axion fields evolving on a Newtonian expanding background, which are
coupled via $\Phi$,

$\displaystyle\begin{split}i\hbar&\dfrac{\partial\psi_{1}}{\partial
t}=\left(-\dfrac{\hbar^{2}}{2m_{1}a^{2}}\nabla^{2}+m_{1}\Phi\right)\psi_{1},\\\
&\vdots\\\ i\hbar&\dfrac{\partial\psi_{N}}{\partial
t}=\left(-\dfrac{\hbar^{2}}{2m_{N}a^{2}}\nabla^{2}+m_{N}\Phi\right)\psi_{N},\\\
&\dfrac{\nabla^{2}\Phi}{4\pi
Ga^{2}}=\left|\psi_{1}\right|^{2}+...+\left|\psi_{N}\right|^{2}-\dfrac{3H^{2}}{8\pi
G}.\end{split}$ (1)

Let us focus on the $N=2$ case as this is relevant to this paper. For
simplicity, we take into account the fact that the characteristic time for
evolution of the system is short compared to the age of the universe, so $a$
becomes unity and $H$ vanishes. In addition, we also consider the system in
symmetrically spherical coordinate and find the stationary solution expressed
by $\psi_{i}(\textbf{x},t)=\psi_{i}(\textbf{x})e^{-iE_{i}t/\hbar}$.
Futhermore, to simplify physical constants, we rescale these quantities into
the dimensionless variables Ruffini1969

$\displaystyle r=\dfrac{\hbar^{2}}{2m^{2}GM}\tilde{r},$
$\displaystyle\psi_{i}=\left(\dfrac{2G^{3}m^{6}M^{4}}{\pi\hbar^{6}M_{i}}\right)^{1/2}\tilde{\psi_{i}},$
(2) $\displaystyle\Phi=\dfrac{2G^{2}M^{2}m^{2}}{\hbar^{2}}\tilde{\Phi},$
$\displaystyle E_{i}=\dfrac{2G^{2}M^{2}m^{2}}{\hbar^{2}}m_{i}\tilde{E_{i}}.$

where $M_{1},M_{2}$ are the total masses of the gravitational structure formed
by each axion and $m$, $M$ are the scale parameters which could be determined
in a specific system. $\Phi,E_{1},E_{2}$ are one-particle potential and one-
particle energy respectively, $\psi_{1},\psi_{2}$ are one-particle
wavefunctions which are normalized individually,
$\int|\psi_{i}|^{2}d^{3}\textbf{x}=1$. This normalization of wavefunctions
also implies $M_{tot}=\int\rho(\textbf{x})d^{3}\textbf{x}=M_{1}+M_{2}$, where
$\rho(\textbf{x})=M_{1}|\psi_{1}|^{2}+M_{2}|\psi_{2}|^{2}$. Finally, we obtain
a system of scale-invariant equations

$\displaystyle\begin{split}&\dfrac{\partial^{2}\tilde{\psi}_{1}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\psi}_{1}}{\partial\tilde{r}}+\left(\dfrac{m_{1}}{m}\right)^{2}(\tilde{\Phi}-\tilde{E}_{1})\tilde{\psi}_{1},\\\
&\dfrac{\partial^{2}\tilde{\psi}_{2}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\psi}_{2}}{\partial\tilde{r}}+\left(\dfrac{m_{2}}{m}\right)^{2}(\tilde{\Phi}-\tilde{E}_{2})\tilde{\psi}_{2},\\\
&\dfrac{\partial^{2}\tilde{\Phi}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\Phi}}{\partial\tilde{r}}+\left|\tilde{\psi}_{1}\right|^{2}+\left|\tilde{\psi}_{2}\right|^{2}.\end{split}$
(3)

These equations can be solved numerically under some necessary constraints

$\displaystyle\begin{split}&\tilde{\psi}^{\prime}_{1}(0)=\tilde{\psi}^{\prime}_{2}(0)=\tilde{\Phi}(0)=0\\\
&\tilde{\psi}_{1}(\infty)=\tilde{\psi}_{2}(\infty)=0,\\\
&\int^{\infty}_{0}|\tilde{\psi}_{i}|^{2}\tilde{r}^{2}d\tilde{r}=\dfrac{M_{i}}{M},\end{split}$
(4)

given the ratio $m_{1}/m_{2}$ and $M_{1}/M_{2}$, once we have solution in term
of $\tilde{\psi}_{i}(\tilde{r})$ we can find the corresponding density profile
at any scale by choosing appropriate value for $m$ and $M$. Because of the
existing shift symmetry in these equations, we can choose an arbitrary value
for the gravitational potential at the origin without changing the solution of
$\psi(r)$. As an illustration, we approach the problem by initially setting
the central values of wavefunctions $\tilde{\psi_{1}}(0)/\tilde{\psi_{2}}(0)$
and identify the normalization factors after the corresponding solution has
found, as a result, if
$\int^{\infty}_{0}|\tilde{\psi}_{i}|^{2}\tilde{r}^{2}d\tilde{r}=\alpha_{i}$,
the following density profile

$\displaystyle\rho(r)=\dfrac{2M^{4}m^{6}G^{3}}{\pi\hbar^{6}}\left(\left|\tilde{\psi_{1}}\right|^{2}+\left|\tilde{\psi_{2}}\right|^{2}\right),$
(5)

where $\tilde{\psi}_{1},\tilde{\psi}_{2}$ are solution of (LABEL:eq3),
describes a soliton with the total mass $M_{tot}=(\alpha_{1}+\alpha_{2})M$.

$\tilde{\psi}_{2}(0)$ / $\tilde{\psi}_{1}(0)$ | $\tilde{E}_{1}$ | $\tilde{E}_{2}$ | $M_{2}/M_{tot}$
---|---|---|---
0.5 | 1.68251566 | 0.69765202 | 2.7%
1.0 | 1.91670679 | 0.83895681 | 8.8%
2.0 | 2.59599377 | 1.24841522 | 23%
3.0 | 3.38136093 | 1.72370456 | 35%
Table 1: Stationary double-solution solutions in terms of $\tilde{E}_{1}$ and
$\tilde{E}_{2}$ considered for different central values of rescaled
wavefunctions with $m_{2}/m_{1}=3$. Notice that the accuracy of
$\tilde{E}_{1}$ and $\tilde{E}_{2}$ is extremely sensitive to the behavior of
the regular solution in this two-axion problem. Because the generic solution
blows up due to the existence of a movable-pole at a finite radius
$\tilde{r}_{sin}$, we need to impose the infinity boundary condition at some
$\tilde{r}_{max}<\tilde{r}_{sin}$ and gradually move this pole to a larger
radius. Hence, the solutions above are regular up to $r_{max}=9$.
$\tilde{E}_{1}$ and $\tilde{E}_{2}$ can also vary depending on the bins size
of $\tilde{r}$ in Runge-Kutte algorithm, here taken as $\Delta\tilde{r}=0.02$.
The detailed numerical method we used to solve for these values will be
presented in a future work (in preparation). Figure 2: Soliton mass density
profiles calculated in four cases of Table 1, with
$M_{tot}=M_{1}+M_{2}=10^{9}\text{ M}_{\odot}$ and $m_{2}=3m_{1}=2.4\times
10^{-22}$ eV, where $M_{i}$ is the total masses of individual soliton. We see
that as the fraction of light axion increases, the two-axion profile
approaches the single-axion one and vice versa.

As an illustration, assume $m_{1}/m_{2}=3$ and $m^{2}=m_{1}m_{2}$, a few
solutions are listed in Table 1. Fig. 2 compares these solutions with the
universal single-axion profile suggested in Schive2014a . Using the same setup
with $m_{2}=m_{1}/3=8\times 10^{-23}$ eV and
$\tilde{\psi}_{1}(0)=3\tilde{\psi}_{2}(0)$, we evaluate the magnitude of this
mutual influence on the nested soliton profile in terms of the residual
density deviation of the coupled case as a function of radius, as illustrated
in Fig. 3. The deviation from a simple sum of the two solitonic cores is
approximately 12%, as shown in the lower panel, and becomes much less than 1%
for our mass range of interest here, $m_{2}/m_{1}\sim 100$. We conclude that
soliton mass density profiles from different axions hardly influence each
other’s presence for the large axion hierarchies of relevance here (in Fig.
1), so we may accurately describe overlapping solitons from different axions
with the simple sum of single mass density profiles Schive2014a ,

$\displaystyle\rho_{s}(r)$
$\displaystyle=\rho_{s,1}(r)+\rho_{s,2}(r)\qquad\text{ where}$ (6)
$\displaystyle\rho_{s,i}(r)$ $\displaystyle\simeq\dfrac{1.9\times
10^{10}(m_{i}/10^{-22}\text{ eV})^{-2}(r_{c,i}/\text{pc})^{-4}}{[1+9.1\times
10^{-2}(r/r_{c,i})^{2}]^{8}}\text{ M}_{\odot}/\text{pc}^{3},$ (7)

where the core radius $r_{c}$ is the radius at which the density is half the
peak density, i.e., $\rho_{s}(r_{c})=\rho_{s}(0)/2$. This numerically-fitted
profile of $\rho_{s}(r)$ is well approximated up to $r\sim 3r_{c}$. Each
soliton follows its own individual core mass-core radius relation

$\displaystyle M_{c,i}$ $\displaystyle\simeq\dfrac{5.5\times
10^{10}}{(m_{i}/10^{-22}\text{ eV})^{2}(r_{c,i}/\text{pc})}\text{ M}_{\odot}.$
(8)

Here, the core mass $M_{c}$ is the mass enclosed within $r_{c}$ and the
(total) soliton mass is $M_{sol}\sim 4M_{c}$, with $r_{sol}\sim 5r_{c}$. Since
the independent solitons approximation in the soliton-in-soliton picture is
best when we restrict ourselves to radii smaller than the core radii, we shall
apply this result to within the core radii as much as possible.

Figure 3: Upper panel: Nested soliton profile calculated with three different methods with the same central density in each case, $\rho_{2}(0)=9\rho_{1}(0)\simeq 6.3\text{ M}_{\odot}/\text{pc}^{3}$, and $m_{2}=3m_{1}=8\times 10^{-23}$ eV. The red-solid curve shows the numerical solution of the coupled two-axion system (1), while the black-solid curve is obtained when simply adding two independent single-axion soliton profiles (7), neglecting the coupling between them. The solid-blue curve (that almost overlaps with the solid red curve) is obtained numerically assuming the heavier axion is solely influenced by the relatively flat gravitational potential of the lighter axion. Details can be found in Appendix A. Lower panel: Dashed-colored curves show the density difference $\Delta\rho$ between the solid black and blue curves above for different $m_{2}/m_{1}>1$ fractions. The corresponding total mass deviation $\Delta M/M$ characterizes the correction in each case, where $M$ is the mass scale corresponding to each case, which is computed from the normalization condition. Structure |  $m$ (eV) | $M_{c}(\text{M}_{\odot})$ |  $r_{c}$ (pc)
---|---|---|---
Carina | $1.9^{+1.9}_{-0.7}\times 10^{-22}$ | $2.5^{+7.7}_{-1.9}\times 10^{7}$ | $575^{+201}_{-195}$
Draco | $8.1^{+4.8}_{-3.0}\times 10^{-23}$ | $1.1^{+1.9}_{-0.7}\times 10^{8}$ | $724^{+231}_{-175}$
Fornax | $1.6^{+0.2}_{-0.2}\times 10^{-22}$ | $3.4^{+1.2}_{-0.9}\times 10^{7}$ | $617^{+75}_{-67}$
Leo I | $1.5^{+1.2}_{-0.5}\times 10^{-22}$ | $4.7^{+10.8}_{-3.3}\times 10^{7}$ | $513^{+179}_{-158}$
Leo II | $5.1^{+7.5}_{-3.2}\times 10^{-22}$ | $1.0^{+7.7}_{-0.9}\times 10^{7}$ | $200^{+268}_{-85}$
Sculptor | $2.0^{+0.5}_{-0.3}\times 10^{-22}$ | $3.0^{+1.08}_{-1.1}\times 10^{7}$ | $447^{+66}_{-66}$
Sextans | $2.5-9.5\times 10^{-23}$ | $2.6-17.8\times 10^{7}$ | $1.5-3\times 10^{3}$
Ursa Minor | $2.5^{+2.3}_{-1.0}\times 10^{-22}$ | $2.3^{+7.2}_{-1.8}\times 10^{7}$ | $389^{+200}_{-138}$
Segue 1 | $5.4^{+8.1}_{-4.0}\times 10^{-21}$ | $6.5^{+20.5}_{-4.2}\times 10^{5}$ | $31^{+114}_{-14}$
Reticulum II | $1.3^{+0.8}_{-1.0}\times 10^{-20}$ | $2.3^{+3.7}_{-1.4}\times 10^{5}$ | $36^{+1086}$
Carina II | $9.3^{+3.9}_{-5.7}\times 10^{-22}$ | $3.0^{+3.2}_{-1.4}\times 10^{7}$ | $\lesssim 676$
Hydrus 1 | $1.5^{+1.2}\times 10^{-20}$ | $1.3^{+4.8}_{-0.6}\times 10^{5}$ | $79^{+1469}$
Draco II* | $\sim 5.6\times 10^{-22}$ | $\sim 1.7\times 10^{7}$ | $\sim 105$
Triangulum II* | $\sim 3.8\times 10^{-22}$ | $\sim 2.4\times 10^{7}$ | $\sim 160$
Eridanus II* | $\gtrsim 6.0\times 10^{-20}$ | $\sim 1.5\times 10^{4}$ | $\lesssim 13$
Antlia 2 | $0.6-1.4\times 10^{-22}$ | $1.0-4.4\times 10^{7}$ | $1.8-3.4\times 10^{3}$
Milky Way | $0.6-1.4\times 10^{-22}$ | $1.5-6.6\times 10^{8}$ | $117-266$
NSC of MW | $2.2-10\times 10^{-20}$ | $1.7-14.5\times 10^{5}$ | $0.15-0.45$
NSC of NGC205 | $\gtrsim 2.5\times 10^{-20}$ | $<9.0\times 10^{4}$ | $0.9-20$
MGC | $\lesssim 10^{-22}$ | $\sim 5\times 10^{9}$ | $<1000$
Tucana 47 | $\gtrsim 5.4\times 10^{-19}$ | $\sim 6.9\times 10^{3}$ | $\lesssim 0.27$
Table 2: The axion mass $m$ (eV), the solitonic core mass $M_{c}$
($M_{\odot}$) and core radius $r_{c}$ (pc) fitting with different cosmological
structure including globular clusters, dwarf galaxies, elliptical galaxies,
clusters of galaxies and nuclear star clusters. *The core mass values for Tri-
II, Dra-II and Eri-II are inferred from different speculations on the total
halo masses as in Calabrese2016hmp ; Marsh2018zyw : $2\times 10^{10}\text{
M}_{\odot}$ and $1.2\times 10^{7}\text{ M}_{\odot}$ respectively.

##  III Solitonic Core Mass - Core Radius

Now we can discuss the observational data for the $\psi$DM model in Fig. 1 and
Table 2. This soliton core mass-core radius plot provides evidence that more
than one axion exists; it also allows us to make predictions (e.g., the size
of the massive galaxy clusters). Note that the masses involved spans 7 orders
of magnitude while the distance scale spans 5 orders of magnitude. The only
theoretical parameters are the axion masses.

$\bullet$ Dwarf Spheroidal Galaxies (dSphs): In the context of $\psi$DM, the
phase space distribution of stellar velocities and positions have been fitted
by Chen2017a ; Schive2014a  with a Jeans analysis to the classical dwarf
spheroidal galaxies, Carina, Draco, Fornax, Leo I, Sculptor, Ursa Minor,
providing strong evidence for an ultralight axion of $\simeq 10^{-22}$ eV,
that follows well the soliton mass - halo relation predicted by the $\psi$DM
simulationsSchive2016 . By assuming Gaussian distributions for the quoted
values of axion masses and cores radius in Table 1 of Chen2017a , the mass
cores of these galaxies are calculated using the relation (8). In addition we
include the Sextans dSph galaxy also analysed in the $\psi$DM
contextBroadhurst:2019fsl .

$\bullet$ Ultra-faint Dwarf Galaxies: These faint objects, now classified as
galaxies, are distinctly smaller and much less luminous than the classical
dwarf spheroidals with the largest mass-to-light among MW satellites,
dominated DM. Currently there is a paucity of accurate stellar velocity
measurements for this class of galaxy so that it has not been possible to
distinguish between a cored or cusped mass profile, making their origin quite
uncertain. Here we apply our Jeans analysis using individual stars as a first
step to obtain constraints for $\psi$DM for the four best studied examples
Simon2019 : Segue 1, Reticulum II, Carina II, Hydrus 1, which are chosen
because of their data availability and sufficient numbers of member stars
(more than 10), as described in Appendix B. Our values for scale radius and
mass of these four UFDs are shown in black in Figure 1, where it is clear that
although a relatively well defined values are obtained, they differ as a class
from the dwarf spheroidals in terms of their location in this plane, requiring
a significantly higher axion mass. Note, due to the large apparent size of
Carina II and its velocity dispersion, it can be classified as a dSph galaxy,
that we show falls nicely on the low axion mass track of Fig. 1.

We are aware of similar studies Calabrese2016hmp  for two other UFDs including
Triangulum II, Draco II and Marsh2018zyw  for Eridanus II with methods that
are tentative to infer axion masses. In particular, the authors of
Calabrese2016hmp  suggests a large maximum halo mass of Tri-II and Dra-II
about $2\times 10^{10}\text{ M}_{\odot}$ based on the dynamics of MW
satellites, while Marsh2018zyw  makes use of the stability of the old star
cluster residing in the center of Eri-II. With stellar extent of only $<20$
pc, Tri-II and Dra-II are more reasonably hosted by halos of $10^{8-9}\text{
M}_{\odot}$. The speculating disruption of the star cluster in Eri-II is also
likely to be compensated by the tidal disruption effect caused by MW
Schive2019rrw . Substantially increased spectroscopy will allow in the future
more accurately defined mass profiles for this class of small, dark matter
dominated galaxies.

$\bullet$ Ghostly Giant Galaxy: The new discovery by the Gaia satellite of a
nearly invisible, dark matter dominated galaxy, Antlia 2 (Ant-2), by tracking
its star motions Torrealba2018 , poses a strong conflict with cold dark
matter; but its very low mass and large size (more than double any known dwarf
galaxy) fit very nicely with a $m_{1}$ axionic soliton core, as shown in Fig 1
where it extends our predicted relation to lower mass, closer the limiting
Jeans scale. More detail can be found in Ref. Broadhurst:2019fsl .

$\bullet$ Massive Galaxy Cluster: Ref. Chen2018a  reports a compact dark mass
of $\sim 2\times 10^{10}\text{ M}_{\odot}$ near the center of a massive galaxy
cluster (MGC), from the radius of curvature of a small lensed structure in a
well resolved background galaxy lensed through the center of a massive lensing
cluster MACS1149 in recent deep Hubble Frontier Fields images. However, the
size is undetermined, with a radius bounded by about a few hundred pc. This
mass may be an offset black hole ejected from the central luminous galaxy or a
compact soliton with mass of $10^{10}\text{ M}_{\odot}$ that is expected at
the bottom of the potential of a massive cluster of $10^{15}\text{
M}_{\odot}$. Assuming that this is due to the $m_{1}$ axion, we predict an
extremely compact soliton with core radius $\simeq 10$ pc, corresponding to
the smaller de-Broglie wavelength of a massive cluster.

$\bullet$ Nested soliton in Milky Way and Central Nuclear Star Cluster: Using
a scaling relation derived from the simulations between the mass of soliton
core and its host virial mass, $M_{c}\propto M_{h}^{1/3}$ Schive2014a ;
Schive2014b ; Schive2016 ; Veltmaat2018 , Ref. DeMartino2018zkx  matches the
missing mass recently found within the central $\simeq 100$ pc of Milky Way
Portail2016 ; Zoccali2014  with a soliton of $\sim 10^{9}\text{ M}_{\odot}$
through Jeans analysis, corresponds to an axion with mass $\simeq 8\times
10^{-23}$ eV. This massive concentrated soliton explains well the projected
radial enhancement of bulge star velocity dispersion peaking at $\simeq 130$
km/s, that is $50$km/s above the general bulge level of $\simeq 80$ km/s.

The origin of the dynamically distinct nuclear star cluster present in the
Milky Way (and in most other galaxies) remains unclear, in particular the
presence of a core profile on a scale of $\simeq 1$ pc Buchholz2009  in the
old star population, for which a old stellar cusp is firmly predicted but not
evident Merritt2010 . The presence of excess unseen matter of $\simeq
10^{6}\text{ M}_{\odot}$ may be implied on a scale of $\simeq 0.4$ pc,
revealed by the high velocity orbit of the maser star IRS9, and excessive
proper motions of other stars at this radius Schodel2018 , that imply an
extended mass, additional to the central black hole on a parsec scale. Whether
this mass can be accounted for by stars of the surrounding NSC is unclear,
requiring a better understanding of the stellar mass function in this region.
An upper limit to the mass of the DM is set by the dynamics of the stellar
motion in the NSC which implies a total mass of $3\times 10^{7}\text{
M}_{\odot}$ on a scale of 3-5 pc radius of the NSC. Most of this dynamical
based mass is thought to be stellar, though the uncertain choice of initial
stellar mass function (IMF) means that as much as half this mass may not be
stellar in the case of a Chabrier IMF rather than the Salpeter form
Schodel2018 ; Alexander2005 . Under the conditions pointed above, we derive a
smaller soliton by matching the enclosed mass at the location of star IRS9
while also satisfying an upper bound closer to the center from the orbit of
star S2 Abuter2018drb . This “inner” soliton implies a heavier axion,
$m_{2}\simeq 10^{-20}$ eV, than the lighter axion $m_{1}$ responsible for the
“outer” soliton associated with bulge star dynamics of the Milky Way described
above and shown in Fig. 4. The solution of an intermediate black hole on top
of the outer soliton Brax:2019npi  for this structure is viable but less
likely providing that the two bounds are both satisfied.

$\bullet$ Early-type Galactic Nuclei of NGC 205: Using high resolution data
from Hubble Space Telescope, Ref. Nguyen2018  shows statistical evidence for
the lack of a central black hole inside the dwarf elliptical (dE) NGC 205,
where its nuclear star cluster resides. The decline of stellar mean velocities
towards the center might imply the existence of a dark soliton of $\sim\times
10^{5}\text{ M}_{\odot}$. However, uncertainties in the measured photometric
and dynamical masses of that nuclear star cluster does not rule out the
possibility of no soliton there, leaving only an upper bound of core mass. We
assume reasonably that the soliton radius to be of order of the star cluster’s
size and with the lower bound is not too small to avoid having a point-like
gravitational potential that is excluded by the declining central velocity
dispersion Nguyen2018 . Note that the extent of any solitonic DM profile is
very unclear in the data, but can be reasonably expected extend beyond the
stellar profile as it would have acted as a seed for gas attraction, so that
stars that subsequently form from the cooling of this gas should not be
expected to orbit beyond the soliton radius but to be generally of a smaller
scale. This means the total mass of any such dark solitonic component may be
substantially higher - nearly an order of magnitude may not be surprising if
core radius is twice as large as the scale length of the stars. Hence we treat
our constraint here cautiously, as indicated by the dashed error bars in Fig.
1 and we regard these values as underestimates of the radius and mass which we
may not be surprised to be revised towards a lighter axion mass than currently
indicated in Fig. 1.

$\bullet$ Globular Cluster 47 Tuc: Compact dark mass has been reported in one
of the best studied Globular Cluster within the Milky Way of $\simeq 3\times
10^{3}\text{ M}_{\odot}$ and naturally interpreted as long anticipated
intermediate mass black holes denBrok2013 ; Kiziltan2017 . However, support
for a point mass is shown to be lacking in the case of 47 Tuc in a new high
resolution stellar proper motion study Mann2018xkm  that prefers an extended
excess of binary stars and stellar mass black holes combined on a scale of 0.1
pc. A soliton explanation has been advocated for these dark excesses of 47 Tuc
Emami2018rxq  corresponding to an axion of mass $\simeq 0.5\times 10^{-18}$
eV, to account for the central dark core mass with an upper bound for the
soliton core radius $\sim 0.03$ pc. For a smaller $r_{c}$, the axion becomes
heavier. More data is needed to definitively test this tentative axion-soliton
explanation.

Figure 4: Fitting the soliton profile with the enclosed mass of the Nuclear
Star Cluster in the Milky Way, using the motion of the star IRS9 at $r\sim
0.365$ pc from the central black hole Reid2006cu . Here the enclosed mass
$M(r)$ inside radius $r$ is computed by integrating the soliton profile (7),
after subtracting a black hole mass of $\sim 4\times 10^{6}\text{ M}_{\odot}$
and stars at the center of NSC. As a reference for comparison, we also
highlight the mass range (the blue band) for the mass profile of solitons
corresponding to the lighter axion $m_{1}$ which have a much wider core radius
and relatively low central density DeMartino2018zkx . Note, the accurate orbit
of the closest orbiting star S2 provides a useful upper bound on any extended
matter additional to the black hole Abuter2018drb .

##  IV Cosmological evolution of axion relic abundance

Based on the astrophysical data, we expect the $m_{1}$ axion contributes
substantially more than the $m_{2}$ axion to the primordial dark matter
density, while the contribution from the $m_{3}$ axion, if present, is
negligible. This can be estimated roughly from the typical size needed to host
a nested local structure. For the one in Milky Way, the mean density over the
same volume is given by

$\displaystyle\dfrac{\bar{\rho}_{MW}}{\bar{\rho}_{NSC}}=\dfrac{M_{c,MW}}{M_{c,NSC}}\simeq
0.5\times 10^{3}.$ (9)

It is reasonable to expect that the relative primordial densities of the
$m_{1}$ axion and the $m_{2}$ axion is not far from this value,

$\displaystyle\frac{\bar{\rho}_{1}}{\bar{\rho}_{2}}\gg 1.$ (10)

If the primordial density $\bar{\rho}_{2}$ of the $m_{2}$ axion is much
bigger, one expects many more single $m_{2}$ solitons than observed.
Qualitatively, this is consistent with the axion physics in cosmology.

In the early universe, an axion potential is flat as the axion $\phi$ is
massless due to the presence of the shift symmetry, $m_{a}\sim 0$. So an axion
potential contributes an energy density like the vacuum energy whose density
stays more or less constant as the universe expands, until a mass is generated
dynamically. This happens when the Hubble parameter is close to its mass,
i.e., $3H\sim m_{a}$, leading to a potential for the axion,

$\displaystyle V(\phi)\simeq
m_{a}^{2}f_{a}^{2}\left[1-\cos(\phi/f_{a})\right]\sim m_{a}^{2}\phi^{2}/2$
(11)

where $f_{a}\gg m_{a}$ is the so called decay constant. Since there is no
reason that the expectation value of $\phi$ happens to sit precisely at the
bottom of the potential, it starts at $\phi_{ini}\sim f_{a}\theta_{ini}$ and
roll down towards the bottom and begin to oscillate coherently around the
bottom, which has the (simplified) equation of motion

$\displaystyle\ddot{\phi}+3H\phi+m^{2}_{a}\phi=0$ (12)

This is a damped simple harmonic oscillator. When coherent oscillation begins,
we can average over a period to obtain the average energy density $\rho$ (up
to a constant factor)

$\displaystyle\bar{\rho}\simeq[\dot{\phi}^{2}+m^{2}_{a}\phi^{2}]/2$ (13)

which leads to (recall $H=\dot{a}/a$ where $a$ is the cosmic scale factor)

$\displaystyle\dot{\bar{\rho}}=\left(-3H\right)\bar{\rho}\to{\bar{\rho}}\propto{1/a^{3}}$
(14)

so it behaves like matter density. This is the mis-alignment mechanism.

Now, dynamics clearly suggests that the conversion from vacuum energy density
to dark matter density at $a_{osc,j}$ is earlier for a heavier axion, so the
$m_{2}$ axion matter density is diluted earlier than the $m_{1}$ axion
density, ending with a lower $m_{2}$ axion energy density. In terms of the
initial densities $\rho_{ini,j}$, we have

$\displaystyle\dfrac{\bar{\rho}_{1}}{\bar{\rho}_{2}}=\dfrac{\rho_{ini,1}}{\rho_{ini,2}}\left(\dfrac{a_{osc,2}}{a_{osc,1}}\right)^{3},$
(15)

where $a_{osc,j}$ can be found by the condition $3H(a_{osc,j})=m_{j}$. Because
axion fields with the mass range of our interest started oscillating well
before matter-radiation equality, let us consider a radiation-dominated
universe. Assume that the initial densities of two axion fields are comparable
(see below), an exact solution can be solved in Fig. 5. On average, their
relic abundance ratio in late universe is roughly

$\displaystyle\dfrac{\bar{\rho}_{1}}{\bar{\rho}_{2}}\simeq\dfrac{\rho_{ini,1}}{\rho_{ini,2}}\left(\dfrac{m_{1}}{m_{2}}\right)^{-3/2}\simeq\dfrac{\rho_{ini,1}}{\rho_{ini,2}}\times
10^{3}.$ (16)

Let us consider the normalization for $\rho_{ini,j}$ for the $j$ axion,

$\rho_{ini,j}=f_{j}^{2}m_{j}^{2}\theta_{ini,j}$

where the initial angle $\theta_{ini,j}$ takes a random value of order
$\theta_{ini,j}\lesssim{\cal O}(1)$ radian. A priori, the $f_{j}$ are
unrelated. However, if the non-perturbative effect of the two axions are due
to the same dynamics, then it is reasonable to take
$f_{1}^{2}m_{1}^{2}=f_{2}^{2}m_{2}^{2}$ (recall the analogous situation for
the QCD axion, where it is related to the strong interaction pion via the
quark condensate, $f_{\pi}^{2}m_{\pi}^{2}=f_{a}^{2}m_{a}^{2}$.) If so, then

$\displaystyle{\rho_{ini,1}}/{\rho_{ini,2}}\sim
1\quad\rightarrow\quad{\bar{\rho}_{1}}/{\bar{\rho}_{2}}\sim 10^{3}$ (17)

This ratio is consistent with (10) despite being inferred from two different
epochs, showing that our two-axion picture is consistent throughout the
history of the Universe.

In fact, the ratio ${\rho_{ini,1}}/{\rho_{ini,2}}$ can be easily many orders
of magnitude away from unity if their masses have different dynamical origins.
That it is close to unity when we compare (9), (16) and (17) hints that these
two axion masses are generated by the same underlying non-perturbative
dynamics. Applying this to the $m_{3}$ axion implies another 2-3 orders of
magnitude lower in the relic density for the $m_{3}$ axion, which is
negligible compared to the dominant $m_{1}$ axion density, as the
observational data indicates.

Figure 5: The mis-alignment density evolution of the two axion fields with
mass $m_{1}$ and $m_{2}$. In the radiation-dominated universe, the equation
(12) is transformed to $d^{2}y/dx^{2}+(2/x)dy/dx+9x^{2}y=0$, where
$y=\phi/\phi_{ini}$ and $x=a/a_{osc}=a/\sqrt{3H_{0}/m_{a}}$; then it can be
solved numerically with the initial conditions $y(x=1)=1$,
$y^{\prime}(x=1)=0$. The normalized density is given by
$\rho/\rho_{ini}=\left[(1/3x)dy/dx\right]^{2}+y^{2}$. The fields behave as
dark matter with $\bar{\rho}\sim a^{-3}$ when averaging over its oscillating
period, which is roughly $m_{a}^{-1}$, as mentioned in the text.

##  V Discussions and Conclusion

We have focused on the soliton feature as a central prediction of the wave
dark matter picture and generalized to the case of multiple axions, motivated
by String Theory, deriving a stable time-independent, joint ground state
solution. We conclude that when the mass scales of two axions are spaced by
more than a factor of a few in mass, the combined solution asymptotes to a
nested pair of concentric solitons that are dynamically independent and this
allows us to make clear comparisons with the observed structure at the center
of our galaxy and other well resolved galaxies where a dense nuclear star
cluster is commonly present.

We have shown that the lightest axion of $m_{1}\simeq\times 10^{-22}$ eV
implied by both the Fornax Galaxy Schive2014a ; Schive2014b ; Schive2016  and
the central Milky Way DeMartino2018zkx , can also help interpret the massive
compact lens of $10^{10}\text{ M}_{\odot}$ uncovered recently near the center
of a massive galaxy cluster Chen2018a .

Our new Jeans-based dynamical analysis indicates strongly that second more
massive axion of $m_{2}\simeq 10^{-20}$ eV is required to explain dark matter
dominated class of newly discovered ”Ultra Faint Galaxies” (UFDs) found within
the Milky Way, because of the relatively small scale and high density of the
dark matter compared to all other galaxies, including the classical dwarf
Spheroidal class. This second axion, together with the very small orbits of
these UFDs within the Milky Way may imply these galaxy-like objects are more
like nuclear star clusters that may have survived extensive tidal destruction
of the host galaxy. This picture is tentatively supported by the possible
presence of an excess of dark matter in the two best studied nuclear star
clusters, that we have analysed here, at the center of the Milky way and at
the center of the nearby NGC205 dwarf galaxy, that we have shown indicate a
similar axion mass like the UFDs.

Finally, and more tentatively a compact inner mass of $\simeq 3\times
10^{3}\text{ M}_{\odot}$ found in the globular cluster 47 Tuc data may imply a
further axion mass of $m_{3}\gtrsim 0.5\times 10^{-18}$ eV. We have combined
this information in a single core mass-radius plot (Fig. 1) providing a
comprehensive “birds-eye” view of the axion role in structure formation,
concluding that these distinctive astronomical structures require more than
one axion in the wave dark matter context.

A definitive test for the presence of a nested soliton structure may be sought
using pulsar timing residuals imprinted on millisecond pulsars detected at the
Galactic Center Khmelnitsky2013lxt ; DeMartino2017qsa . Many thousands of
pulsars are expected Calore2015bsx , and can account for the GeV gamma-ray
excess Abazajian2011 ; Abazajian2014 , which are being searched with some
success Deneva2008rz . Detection of such pulsars is a major goal of the SKA,
including within the NSC region, given the high density star formation history
in this cluster. Such NSC pulsars located within the nested solitons may show
distinctive multi-frequency timing residuals on the respective Compton time
scales of these independently oscillating scalar fields, of a few hours and a
few months corresponding to $10^{-20}$ eV and $10^{-22}$ eV respectively,
providing a unique soliton signature as a definitive solution to the long
standing Dark Matter puzzle.

## Acknowledgments

We thank Tzihong Chiueh and George Smoot for valuable discussions. SHHT is
supported by the AOE grant AoE/P-404/18-6 issued by the Research Grants
Council (RGC) of the Government of the Hong Kong SAR China.

##  Appendix A Flat-density approximation

Although we have numerically solved the Schrodinger-Poisson equations of
nested system in case the masses of two axions are not quite different,
$m_{2}/m_{1}=3$, it is not practical to exactly solve these equations when
this mass difference is significant because we have to fine-tune two
eigenvalues $\tilde{E}_{1}$ and $\tilde{E}_{2}$ to maintain the stability of
the wavefunction of the heavier axion. This technical issue is generic in
finding numerical solution with multiple boundary conditions Moroz:1998dh . On
the other hand, we know the heavy-axion and light-axion wavefunction become
wider and narrower respectively as either mass or central density difference
increases. This fact implies that in the limit

$\displaystyle m_{2}\gg
m_{1}\quad\text{and}\quad\tilde{\psi_{2}}(0)\gg\tilde{\psi_{1}}(0),$ (18)

the heavy-axion wavefunction approximately evolves on a flat density
distribution background of the light axion. Therefore, we can separately write
the equations governing the evolution of two axions as follow

$\displaystyle\dfrac{\partial^{2}\tilde{\psi}_{2}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\psi}_{2}}{\partial\tilde{r}}+(\hat{\Phi}-\tilde{E}_{2})\tilde{\psi}_{2},$
$\displaystyle\dfrac{\partial^{2}\hat{\Phi}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\hat{\Phi}}{\partial\tilde{r}}+\left|\tilde{\psi}_{2}\right|^{2}+\left|\tilde{\psi}_{1}(0)\right|^{2},$

and

$\displaystyle\dfrac{\partial^{2}\tilde{\psi}_{1}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\psi}_{1}}{\partial\tilde{r}}+(\tilde{\Phi}-\tilde{E}_{1})\tilde{\psi}_{1},$
$\displaystyle\dfrac{\partial^{2}\tilde{\Phi}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\Phi}}{\partial\tilde{r}}+\left|\tilde{\psi}_{2\text{
(solved)}}\right|^{2}+\left|\tilde{\psi}_{1}\right|^{2}.$

They simply return to the familiar single-axion system of equations with an
external source of flat density. Here, we use $\hat{\Phi}$ and $\tilde{\Phi}$
to distinguish effective potentials in each system of equations.

##  Appendix B Jeans analysis for UFDs

The spherically symmetric Jeans equation when being projected along the line
of sight yields the following dispersion velocity Binney2008

$\displaystyle\sigma^{2}_{p}=\dfrac{2}{I(R)}\int_{R}^{\infty}\left(1-\beta\dfrac{R^{2}}{r^{2}}\right)\rho_{*}\sigma^{2}_{r}(r)\dfrac{rdr}{\sqrt{r^{2}-R^{2}}},$
(19)

where $I(R)$ is stellar surface density. From observations for UFDs, it is
standard to fit this density with the Plummer profile given by Walker2009b

$\displaystyle I(R)=\dfrac{L}{\pi
R_{h}^{2}}\left(1+\dfrac{R^{2}}{R_{h}^{2}}\right)^{-2},$ (20)

which deprojected in three-dimensional density of the form

$\displaystyle\rho_{*}(r)=\dfrac{3L}{4\pi
R_{h}^{3}}\left(1+\dfrac{r^{2}}{R_{h}^{2}}\right)^{-5/2}.$ (21)

Here, $R_{h}$ is called the half-light (effective) radius of the corresponding
galaxy and its projected value can be related to the deprojected one by
$r_{1/2}=(4/3)R_{h}$.

The high mass-to-light ratio of UFDs identifies them as extreme cases for
dark-matter dominated systems. In that case, stars are treated as ”tracers”
for an underlying gravitational potential generated by the corresponding mass
distribution and their radial velocity dispersion can be calculated via
Lokas:2003ks

$\displaystyle\rho_{*}\sigma^{2}_{r}(r)=Gr^{-2\beta}\int_{r}^{\infty}s^{2\beta-2}\rho_{*}(s)M(s)ds.$
(22)

The velocity anisotropy $\beta$ is ambiguous for most of UFDs. However, a
general radius-dependent profile can be taken into account as

$\displaystyle\beta(r)=(\beta_{\infty}-\beta_{0})\dfrac{r^{2}}{r^{2}+r^{2}_{\beta}}+\beta_{0},$
(23)

where $\beta_{0}$ and $\beta_{\infty}$ are asymptotic values at inner up to
outer edge of the dark matter halo, $r_{\beta}$ is the scale radius of the
stellar distribution. Because $\beta\equiv 1-\sigma^{2}_{r}/\sigma^{2}_{t}<1$
at all radius, we need to choose priors for $\beta_{0},\beta_{\infty}$
satisfied this condition later on.

In this work, we consider the dark matter profile with an inner soliton
extending to NFW profile
$\rho_{NFW}(r)=\rho_{0}(r/r_{s})^{-1}(1+r/r_{s})^{-2}$ at a transition radius
of $r_{trans}=3r_{c}$ Marsh:2015wka . The NFW characteristic density
$\rho_{0}$ is determined by setting
$\rho_{s}(r_{trans})=\rho_{NFW}(r_{trans})$; $r_{s}$ characterizes the scale
radius of NFW halos. The enclosed mass associated with this profile can be
computed explicitly Chen:2016unw ; Navarro:1995iw

$M(r)=M_{s}(r)=\dfrac{4.2\times 10^{7}\text{
M}_{\odot}}{m_{22}^{2}(r/\text{pc})}\dfrac{1}{(a^{2}+1)^{7}}\left(3465a^{13}\right.\\\
+23100a^{11}+65373a^{9}+101376a^{7}+92323a^{5}\\\
\left.+48580a^{3}-3465a+3465(a^{2}+1)^{7}\tan^{-1}(a)\right)\\\ \text{ with
}a=\left(2^{1/8}-1\right)^{1/2}(r/r_{c})\text{ when }r<r_{trans};$ (24)
$M(r)=M_{s}(r_{trans})+M_{NFW}(r)-M_{NFW}(r_{trans})\\\
=M_{s}(r_{trans})+4\pi\rho_{0}r^{3}_{s}\left[\ln\left(\dfrac{r_{s}+r}{r_{s}+r_{trans}}\right)\right.\\\
\left.-\dfrac{r}{r_{s}+r}+\dfrac{r_{trans}}{r_{s}+r_{trans}}\right]\text{when
}r>r_{trans}.$ (25)

Notice that the tail of NFW profile makes the enclosed mass divergent as
$\sim\ln r$ when $r\rightarrow\infty$.

Therefore, to calculate the mass we integrate this profile up to a cut-off
radius given by $r_{cut}=\text{Max}\left[r_{lim},r_{Roche}\right]$. Here,
$r_{lim}$ is about 10 times of the maximum stellar-extended radius;
$r_{Roche}$ is the standard Roche limit for a point mass of
$M_{halo}=10^{9}\text{ M}_{\odot}$ orbiting around Milky Way with a velocity
dispersion of $\sigma_{MW}=200\text{ km s}^{-1}$,

$\displaystyle
r_{Roche}\simeq\left(\dfrac{GM_{h}D^{2}}{2\sigma_{MW}^{2}}\right)^{1/3},$ (26)

and $D$ is the distance to the center of the dwarf satellite from the Earth.

Following the same method in Strigari2008i ; Wolf2010 , we have to determine a
viable range for 6 free parameters referred as

$\displaystyle\vec{\theta}=\\{m_{22},r_{c},r_{s},r_{\beta},\beta_{0},\beta_{\infty}\\}.$
(27)

The likelihood function is chosen as a two-component Gaussian distribution

$\displaystyle\mathcal{L}(\vec{\theta})=\prod_{i=1}^{N}\dfrac{1}{\sqrt{2\pi\sigma^{2}_{i}}}\text{exp}\left[-\dfrac{1}{2}\dfrac{(v_{i}-u)^{2}}{\sigma_{i}^{2}}\right].$
(28)

From observation, $v_{i}$ and $u$ are the heliocentric velocity for the
$i^{th}$ star and the systemic velocity of the satellite, respectively. The
velocity distribution is weighted by the total variance
$\sigma^{2}_{i}=\sigma^{2}_{i,th}+\sigma^{2}_{i,m}$, where $\sigma^{2}_{i,th}$
is computed theoretically from (19) at a certain projected radius $R_{i}$ for
an individual star while $\sigma^{2}_{i,m}$ refers to its velocity uncertainty
in measurement.

UFDs |  $R_{h}$ (pc) |  $D$ (kpc) |  $u$ (km s-1) | Ref
---|---|---|---|---
Segue 1 | 29 | 23 | 208.5 | Geha:2008zr ; Simon:2010zr
Reticulum II | 55 | 32 | 62.8 | Simon:2015fdw
Carina II | 91 | 37.4 | 477.2 | Li:2018fdw
Hydrus 1 | 53 | 27.6 | 80.4 | Koposov2018
Table 3: The structural parameters of UFDs and references to data sets used in
this work.

Uniform priors are assumed on the parameters over the following ranges

$\displaystyle-1<\log_{10}(m_{22})<5,$ $\displaystyle\quad
1<\log_{10}(r_{c}/\text{pc})<4,$ (29) $\displaystyle
1<\log_{10}(r_{s}/\text{pc})<4,$ $\displaystyle\quad
1<\log_{10}(r_{\beta}/\text{pc})<4,$ (30) $\displaystyle-10<\beta_{0}<1,$
$\displaystyle\quad-10<\beta_{\infty}<1,$ (31) $\displaystyle 0.2R_{h}<$
$\displaystyle r_{\beta}<r_{lim},$ (32) $\displaystyle
0.2R_{h}<r_{c}<2r_{cut},$ $\displaystyle\quad 0.2R_{h}<r_{c}<2r_{cut}.$ (33)

In Table 3, we show properties of UFD candidates and references to their
stellar data taken to analyze in this work. Most cases are confirmed as
unlikely to be under ongoing tidal disruption; hence they are taken to be in
dynamical equilibrium where Jeans analysis is applicable. Figures 6 illustrate
distributions for two relevant quantities: $M_{c},r_{c}$ which can show their
preferred values when being constrained from the data sets.

Figure 6: 1D and 2D posterior distributions of four UFDs taken from MCMC
chains using emcee Foremanmackey2012 , plotted using corner package corner2016
. Here, we simulate with 240 walkers and 10000 steps for each galaxy; only the
last $3000\times 240$ samples (thinning each 20 points) are chosen to plot.
The red dashed-vertical lines in each 1D diagram show the peak value of the
distributions and the black ones show errors at 60.6% of the peak height. 68%
and 95% confidence levels are shown by contours in the 2D correlation plots.

## References

  * (1) P. A. R. Ade et al. [Planck Collaboration], Astron. Astrophys. 594, A13 (2016) [arXiv:1502.01589 [astro-ph.CO]].
  * (2) R. H. Cyburt, B. D. Fields, K. A. Olive and T. H. Yeh, Rev. Mod. Phys. 88, 015004 (2016) [arXiv:1505.01076 [astro-ph.CO]].
  * (3) W. Hu, R. Barkana and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000) [astro-ph/0003365].
  * (4) H. Y. Schive, T. Chiueh and T. Broadhurst, Nature Phys. 10, 496 (2014) [arXiv:1406.6586 [astro-ph.GA]].
  * (5) H. Y. Schive, M. H. Liao, T. P. Woo, S. K. Wong, T. Chiueh, T. Broadhurst and W.-Y. P. Hwang, Phys. Rev. Lett. 113, no.26, 261302 (2014) [arXiv:1407.7762 [astro-ph.GA]].
  * (6) H. Y. Schive, T. Chiueh, T. Broadhurst and K. W. Huang, Astrophys. J. 818, no.1, 89 (2016) [arXiv:1508.04621 [astro-ph.GA]].
  * (7) L. Hui, J.P. Ostriker, S. Tremaine and E. Witten, Phys. Rev. D 95, no. 4, 043541 (2017) doi:10.1103/PhysRevD.95.043541 [arXiv:1610.08297 [astro-ph.CO]].
  * (8) J. C. Niemeyer, Prog. Part. Nucl. Phys., 103787 doi:10.1016/j.ppnp.2020.103787 [arXiv:1912.07064 [astro-ph.CO]].
  * (9) L. M. Widrow and N. Kaiser, Astrophys. J. 416, L71 (1993).
  * (10) P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov and L. Hernquist, Mon. Not. Roy. Astron. Soc. 471, no.4, 4559 (2017) [arXiv:1705.05845 [astro-ph.CO]].
  * (11) J. Veltmaat, J. C. Niemeyer and B. Schwabe, Phys. Rev. D 98, no.4, 043509 (2018) [arXiv:1804.09647 [astro-ph.CO]].
  * (12) F. S. Guzman and L. A. Urena-Lopez, Astrophys. J. 645, 814 (2006) [astro-ph/0603613].
  * (13) I. De Martino, T. Broadhurst, S. H. Tye, T. Chiueh and H. Schive, Phys. Dark Univ. 28, 100503 (2020) doi:10.1016/j.dark.2020.100503 [arXiv:1807.08153 [astro-ph.GA]].
  * (14) Joshua D. Simon, arXiv:1807.08153 [astro-ph.GA].
  * (15) D. J. E. Marsh and J. C. Niemeyer, arXiv:1810.08543 [astro-ph.CO].
  * (16) E. Calabrese and D. N. Spergel, Mon. Not. Roy. Astron. Soc. 460, no.4, 4397 (2016) [arXiv:1603.07321 [astro-ph.CO]].
  * (17) R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767–1783 (Nov, 1969).
  * (18) S. R. Chen, H. Y. Schive and T. Chiueh, Mon. Not. Roy. Astron. Soc. 468, no.2, 1338 (2017) [arXiv:1606.09030 [astro-ph.GA]].
  * (19) M. Portail, O. Gerhard, C. Wegg and M. Ness, Mon. Not. Roy. Astron. Soc. 465, no.2, 1621 (2017) [arXiv:1608.07954 [astro-ph.GA]].
  * (20) Zoccali, M. et al., Astron. Astrophys. 562 A66 (2014).
  * (21) M. C. Chen, T. Broadhurst, J. Lim, J. M. Diego, Y. Ohyama, H. Ford, and N. Benitez, Astrophys. J. 863 , no.2 135 (2018).
  * (22) D. Crnojevic, D. J. Sand, D. Zaritsky, K. Spekkens, B. Willman, and J. R. Hargis, Astrophys. J. 824 , no.1 L14 (2016).
  * (23) T. S. Li et al. [DES Collaboration], Astrophys. J. 838, no.1, 8 (2017) [arXiv:1611.05052 [astro-ph.GA]].
  * (24) H. Schive, T. Chiueh and T. Broadhurst, [arXiv:1912.09483 [astro-ph.GA]].
  * (25) R. M. Buchholz, R. Schodel and A. Eckart, Astron. Astrophys. 499, 483 (2009) [arXiv:0903.2135 [astro-ph.GA]].
  * (26) D. Merritt, Astrophys. J. 718 , no.2 739 (2010).
  * (27) R. Schodel et al. Astron. Astrophys. 609, A27 (2018) [arXiv:1701.03817 [astro-ph.GA]].
  * (28) T. Alexander, Phys. Rept. 419, 65 (2005) [astro-ph/0508106].
  * (29) R. Abuter et al. [GRAVITY Collaboration], Astron. Astrophys. 615, L15 (2018) [arXiv:1807.09409 [astro-ph.GA]].
  * (30) P. Brax, J. A. R. Cembranos and P. Valageas, Phys. Rev. D 101, no. 2, 023521 (2020) doi:10.1103/PhysRevD.101.023521 [arXiv:1909.02614 [astro-ph.CO]].
  * (31) C. Mann et al., arXiv:1807.03307 [astro-ph.GA].
  * (32) R. Emami, T. Broadhurst, G. Smoot, T. Chiueh and H. N. Luu, Phys. Rev. D 101, no.6, 063006 (2020) doi:10.1103/PhysRevD.101.063006 [arXiv:1806.04518 [astro-ph.CO]].
  * (33) G. Torrealba et al., arXiv:1811.04082 [astro-ph.CO].
  * (34) T. Broadhurst, I. de Martino, H. N. Luu, G. F. Smoot and S.-H. H. Tye, arXiv:1902.10488 [astro-ph.CO].
  * (35) Dieu. D. Nguyen et al. et al., Astrophys. J. 858, 118 (2018) [astro-ph/1711.04314].
  * (36) M. den Brok, G. van de Ven, R. van den Bosch, and L. Watkins, Mon. Not. Roy. Astron. Soc. 438 , no.1 487–493 (2014).
  * (37) B. Kiziltan, H. Baumgardt, and A. Loeb, Nature Phys. 542 , no.7640 203 (2017).
  * (38) M. J. Reid, K. M. Menten, S. Trippe, T. Ott and R. Genzel, Astrophys. J. 659, 378 (2007) [astro-ph/0612164].
  * (39) D. J. E. Marsh, Phys. Rept. 643, 1-79 (2016) doi:10.1016/j.physrep.2016.06.005 [arXiv:1510.07633 [astro-ph.CO]].
  * (40) A. Khmelnitsky and V. Rubakov, JCAP 1402, 019 (2014) [arXiv:1309.5888 [astro-ph.CO]].
  * (41) I. De Martino, T. Broadhurst, S. H. Henry Tye, T. Chiueh, H. Y. Schive and R. Lazkoz, Phys. Rev. Lett. 119, no.22, 221103 (2017) [arXiv:1705.04367 [astro-ph.CO]].
  * (42) F. Calore, M. Di Mauro, F. Donato, J. W. T. Hessels and C. Weniger, Astrophys. J. 827, no.2, 143 (2016) [arXiv:1512.06825 [astro-ph.HE]].
  * (43) K. N. Abazajian, S. Blanchet and J. P. Harding, Phys. Rev. D 84, 103007 (2011) [arXiv:1012.1247 [astro-ph.CO]].
  * (44) K. N. Abazajian, N. Canac, S. Horiuchi and M. Kaplinghat, Phys. Rev. D 90, no.2, 023526 (2014) [arXiv:1402.4090 [astro-ph.HE]].
  * (45) J. S. Deneva et al., Astrophys. J. 703, 2259 (2009) [arXiv:0811.2532 [astro-ph]].
  * (46) I. M. Moroz, R. Penrose and P. Tod, Class. Quant. Grav. 15, 2733 (1998). doi:10.1088/0264-9381/15/9/019
  * (47) Binney J., Tremaine S., 2008, Galactic Dynamics. Princeton Univ. Press, Princeton
  * (48) Walker M. G., Mateo M., Olszewski E. W., Penarrubia J., Wyn Evans N., Gilmore G., Astrophys. J. 704, 1274 (2009b) [arXiv:0906.0341 [astro-ph]]
  * (49) E. L. Lokas and G. A. Mamon, Mon. Not. Roy. Astron. Soc. 343, 401 (2003) doi:10.1046/j.1365-8711.2003.06684.x [astro-ph/0302461].
  * (50) D. J. E. Marsh and A. R. Pop, Mon. Not. Roy. Astron. Soc. 451, no. 3, 2479 (2015) doi:10.1093/mnras/stv1050 [arXiv:1502.03456 [astro-ph.CO]].
  * (51) J. F. Navarro, C. S. Frenk and S. D. M. White, Astrophys. J. 462, 563 (1996) doi:10.1086/177173 [astro-ph/9508025].
  * (52) S. R. Chen, H. Y. Schive and T. Chiueh, Mon. Not. Roy. Astron. Soc. 468, no. 2, 1338 (2017) doi:10.1093/mnras/stx449 [arXiv:1606.09030 [astro-ph.GA]].
  * (53) L. E. Strigari, J. S. Bullock, M. Kaplinghat, J. D. Simon, M. Geha, B. Willman and M. G. Walker, Nature 454, 1096 (2008) doi:10.1038/nature07222 [arXiv:0808.3772 [astro-ph]].
  * (54) Joe Wolf and Gregory D. Martinez and James S. Bullock and Manoj Kaplinghat and Marla Geha and Ricardo R. Munoz and Joshua D. Simon and Frank F. Avedo, Mon. Not. Roy. Astron. Soc. 406, 1220 (2010) [arXiv:0908.2995 [astro-ph.CO]]
  * (55) M. Geha, B. Willman, J. D. Simon, L. E. Strigari, E. N. Kirby, D. R. Law and J. Strader, Astrophys. J. 692, 1464 (2009) doi:10.1088/0004-637X/692/2/1464 [arXiv:0809.2781 [astro-ph]].
  * (56) Joshua D. Simon and Marla Geha et al., Astrophys. J. 733, no.1 (2011) doi:10.1088/0004-637X/733/1/46 [arXiv:1007.4198 [astro-ph]].
  * (57) J. D. Simon et al. [DES Collaboration], Astrophys. J. 808, no. 1, 95 (2015) doi:10.1088/0004-637X/808/1/95 [arXiv:1504.02889 [astro-ph.GA]].
  * (58) T. S. Li et al., Astrophys. J. 857, no. 2, (2018) doi:10.3847/1538-4357/aab666 [arXiv:1802.06810 [astro-ph.GA]].
  * (59) JSergey E. Koposov et al., Mon. Not. Roy. Astron. Soc. 479, 6343 (2018) [arXiv:1804.06430 [astro-ph.CO]]
  * (60) Daniel Foreman-Mackey and David W. Hogg and Dustin Lang and Jonathan Goodman, arXiv:1202.3665 [astro-ph.IM].
  * (61) Daniel Foreman-Mackey, The Journal of Open Source Software (2016) vol.1 no.2 doi:10.21105/joss.00024


# Similarity Based Label Smoothing For Dialogue Generation

Sougata Saha∗, Souvik Das, Rohini Srihari
Department of Computer Science and Engineering
University at Buffalo, New York
{sougatas, souvikda<EMAIL_ADDRESS>
Equal Contribution

###### Abstract

Generative neural conversational systems are generally trained with the
objective of minimizing the entropy loss between the training “hard” targets
and the predicted logits. Often, performance gains and improved generalization
can be achieved by using regularization techniques like label smoothing, which
converts the training “hard” targets to “soft” targets. However, label
smoothing enforces a data independent uniform distribution on the incorrect
training targets, which leads to an incorrect assumption of equi-probable
incorrect targets for each correct target. In this paper we propose and
experiment with incorporating data dependent word similarity based weighing
methods to transforms the uniform distribution of the incorrect target
probabilities in label smoothing, to a more natural distribution based on
semantics. We introduce hyperparameters to control the incorrect target
distribution, and report significant performance gains over networks trained
using standard label smoothing based loss, on two standard open domain
dialogue corpora.

##  1 Introduction

Generative conversational systems rely heavily on language modelling to be
able to generate an appropriate response to a user query. Given a context
which consist of multiple utterances in a conversation, a generative
conversational system can be formulated as a next utterance prediction
problem, where the task is to generate a response utterance conditioned on the
context utterances. With the advent of deep learning and availability of
sufficient training data, parametric models like recurrent neural networks and
transformers are generally implemented to achieve the language modelling task.
Trained by minimizing the expected cross entropy between the training hard
targets and the prediction logits, such models often overfit the training data
and does not generalize well on the test set. Label smoothing proposed by
Szegedy et al. (2015) to improve the performance of Inception net image
classifier on the ImageNet dataset, has widely gained acceptance in Natural
Language Processing tasks as a regularization technique to enhance the
generalization capability of deep neural networks. Vaswani et al. (2017) in
his work “Attention is all you need”, where he proposed the state-of-the-art
transformer architecture, had reported performance gains in machine
translation using label smoothing during training. Unlike other regularization
techniques which constrain the model parameters and hidden representations,
label smoothing augments the actual targets by reducing the actual target
probability, and assigning low probabilities to all classes, following a data
independent uniform distribution, thus preventing the model from predicting
the correct labels too confidently during training. However, as pointed out by
Pereyra et al. (2017) and Hinton et al. (2015), the probabilities assigned to
both the correct and incorrect classes constitute the knowledge of a network.
In language modelling, incorporating label smoothing and assigning a uniform
probability to all the incorrect classes can convey an incorrect knowledge to
the model. For example, in response to a user query “How are you doing ?”, if
we want to generate the sentence “I am doing good .” as the next utterance,
given that we have already generated the phrase “I am doing”, “great” and
“awesome” can also convey the same message as “good”. On the other hand, “bad”
would convey a different message, but would be logically correct. However a
random word like “aeroplane” would not make any sense. Hence, if we want to
use label smoothing, we would not want to follow a uniform distribution for
the incorrect classes, and rather weigh them using a weighing mechanism which
can present such knowledge to the model. In this paper we present ways of
imparting such information by modifying the data independent uniform
distribution in label smoothing with a more appropriate data dependent
distribution, which is proportional to the pre-trained word embedding
similarity between the actual target and the incorrect targets.

##  2 Related Work

Numerous techniques have been introduced to enhance the generalization
capability of neural networks. Although, as pointed out by Pereyra et al.
(2017), substantial advancements have been made in regularizing model
parameters, but not much work has been done in understanding external
regularization techniques like label smoothing or target data augmentation. We
can broadly categorize the most recent solutions to attain generalization in
conversations in the following two categories.
Loss function augmentation: Li et al. (2016) proposes using Maximum Mutual
Information along with the Cross Entropy loss, in order to penalize bland
responses like “I do not know”, which are frequent in conversational datasets.
Jiang et al. (2019) analysed that the Cross Entropy loss function prefers
frequent tokens, which leads to generating bland responses. Hence they
proposed augmenting the Cross Entropy Loss with a frequency based corpus
dependent weighing mechanism, in order to yield diverse responses. Wang et al.
(2020) experiments with using optimal transport to match sequences generated
in the teacher and student modes, and increasing performance of student forced
networks on the test dataset by reducing the gap between the two modes.
Data augmentation: Cai et al. (2020) demonstrated that conversational datasets
generally don’t exhibit coherence in query response pairs, which affect the
Cross Entropy loss. They propose a training data augmentation module, which
can not only replace words in the actual target response with similar words
using BERT Devlin et al. (2019), but also augment the style of the response,
preserving the meaning. They further introduced a neural weighting mechanism,
which can assign weights or importance to the augmented and golden training
data, and report significant performance gains. Kang and Hashimoto (2020)
demonstrated that the log loss is not robust to noise, and hence proposed
truncating the distribution of the training targets to achieve an easy to
optimize and more robust loss function. He and Glass (2020) introduced a
network which can provide negative generated samples, and train the generation
model to maximize the log likelihood of training data while minimizing the
likelihood of negative samples. Our proposed method falls in the first
category, as we do not augment the training data, and instead augment the
probability of incorrect labels for each correct label.

##  3 Methods and Experiments

We experiment with ways to augment the data independent uniform distribution
enforced by Label Smoothing. Let $U_{i}$ be an utterance consisting of words
$\\{w_{j}\\}_{j=1}^{N}$, where $N$ is the number of words in the utterance.
For each word $w_{j}$, in label smoothing a probability $1-s$ is assigned to
the true label $w_{j}$, and a probability of $s$ (smoothing factor) is
distributed uniformly among the rest of the $k$ words in the vocabulary. We
augment the distribution of the incorrect class by weighting the smoothing
factor $s$ according to the cosine similarity between the Glove Pennington et
al. (2014) word embedding of the correct word in the training data and all the
words in the vocabulary. Thus, if the correct word to be predicted is “good”,
then the words “great” and “awesome” in the vocabulary would get a higher
proportion of the smoothing factor $s$, compared to an unrelated word like
“aeroplane”, thus presenting a more correct knowledge to the model.
Mathematically, let $\vec{w_{j}}$ be the Glove word embedding of word $w_{j}$,
$\mathbf{W_{k}}$ be a matrix containing the Glove word embedding for all the
words in the vocabulary (including $w_{j}$), $\vec{w}_{j\\_sim}$ be the vector
of cosine similarity between the word $w_{j}$ and all the words in the
vocabulary. Since Glove word embeddings are learned representations, they can
be noisy. Hence, we introduce a threshold $t$, below which we set the cosine
similarity value in $\vec{w}_{j\\_sim}$ as 0. We achieve this by introducing a
mask $mask_{j}$, and multiply the similarity vector with the mask. The
resulting vector is normalised to lie between 0 and 1, and finally multiplied
by $s$. We treat $t$ as a model hyperparameter, and is tuned using grid
search. We further reason that although Glove embeddings are learned from text
corpora, there are possibilities that dissimilar words can lie in close
proximity in the embedding space, resulting in a high cosine similarity score,
and presenting an incorrect knowledge to the model. To circumvent this
problem, we further experiment with filtering out the cosine similarities of
dissimilar words based on WordNet sysnets pri (2010), which we achieve by
implementing another mask $mask_{j}^{{}^{\prime}}$.

$\displaystyle\vec{w}_{j\\_sim}$
$\displaystyle=\frac{\vec{w_{j}}\cdot\mathbf{W_{k}}}{||\vec{w_{j}}||\cdot||\mathbf{W_{k}}||}$
$\displaystyle\vec{w}_{j\\_dist}$
$\displaystyle=\frac{\vec{w}_{j\\_sim}*mask_{j}*mask_{j}^{{}^{\prime}}}{\sum(\vec{w}_{j\\_sim}*mask_{j}*mask_{j}^{{}^{\prime}})}*s$
$\displaystyle\vec{w}_{j\\_dist}[j^{*}]$ $\displaystyle=1-s$
$\displaystyle\text{where, }mask_{j}$
$\displaystyle=\left\\{\begin{array}[]{@{}ll@{}}0,&\text{if}\
\vec{w}_{j\\_sim}<=t\\\ 0,&\text{if}\ \vec{w}_{j\\_sim}=1\\\
1,&\text{otherwise}\end{array}\right.$ $\displaystyle\text{and,
}mask_{j}^{{}^{\prime}}$
$\displaystyle=\left\\{\begin{array}[]{@{}ll@{}}1,&\text{if $w_{k}$ is a
synonym of $w_{j}$}\\\ 0,&\text{otherwise}\end{array}\right.$

###  3.1 Dataset

We perform experiments on (i) The DailyDialog dataset Li et al. (2017) : A
multi-turn open domain dialogue dataset which has 13,118 pertaining to diverse
day-to-day topics , and (ii) The Empathetic Dialogues dataset Rashkin et al.
(2019) : An open domain multi-turn dataset consisting of 25,000 conversations
grounded in emotional situations. We use the same training, validation and
testing splits as mentioned in the datasets. We concatenate all the turns in
the query in one long text, and use two special tokens: “[speaker1]” and
“[speaker2]” to distinguish the speakers. In order to speed up computation, we
restrict the context to the most recent 50 tokens, which is determined
analytically from the corpora. Please refer to the supplementary material for
the code and the dataset.

###  3.2 Model

Since the primary scope of this paper is to experiment with different loss
functions, we used a standard transformer encoder-decoder architecture as
proposed by Vaswani et al. (2017), where the encoder encodes the most recent
utterance in the conversation, along with context from the previous turns. The
encoder-decoder comprises of 3 layers each, with 300 dimensional hidden
representation, with 6 attention heads in each multi-headed attention layer.
The embedding layer is populated with 300 dimensional Glove embeddings, which
are trained along with the entire network. Finally, a fully connected linear
layer predicts the next word.

###  3.3 Experiments

We treat the vanilla Cross Entropy (CE) loss, CE loss with label smoothing,
Kullback–Leibler (KL) divergence loss and KL loss with label smoothing as the
baselines. We experiment with different smoothing values
$s\in\\{NA,0.1,0.2\\}$ , cosine similarity thresholds
$t\in\\{NA,0.0,0.5,0.8\\}$, and also perform ablation study to analyze the
usefulness of the WordNet similarity mask $w\in\\{NA,0,1\\}$. In total, we
experiment with 30 different settings for each dataset.

Figure 1: Illustration of probability of the incorrect words. The x-axis =
vocabulary, y-axis = probability. Setting t = 0.8, or using WordNet mask
filters out most words, making the target distribution equivalent to vanilla
CE loss targets. Using t = 0.5 or 0.0 yields a less dramatic effect, and helps
preserving the knowledge of the incorrect labels.

| s = NA | s = 0.1 | s = 0.2 | s = 0.1 | s = 0.2
---|---|---|---|---|---
t = NA | t = NA | t = NA | t = 0 | t = 0.5 | t = 0 | t = 0.5
Dataset | Metric | Loss | w = NA | w = NA | w = NA | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1
DD | SacreBLEU | CE | 1.6625 | 1.8523 | 1.7251 | 1.9896 | 1.7627 |  | 2.1936
---
(+ 12.67 %)
1.8575 | 2.1158 | 1.8020 | 2.0536 | 1.9302
KL | 1.9469 | 1.7536 | 1.7931 | 1.8459 | 1.8181 | 1.9128 | 1.8858 | 1.9387 | 1.8092 | 1.7292 | 1.8856
ROUGE L | CE | 0.1209 | 0.1243 | 0.1243 | 0.1238 | 0.1209 |  | 0.1270
---
(+ 0.57 %)
0.1260 | 0.1217 | 0.1204 | 0.1238 | 0.1244
KL | 0.1263 | 0.1223 | 0.1233 | 0.1227 | 0.1264 | 0.1243 | 0.1242 | 0.1223 | 0.1253 | 0.1232 | 0.1234
METEOR | CE | 0.1244 | 0.1324 | 0.1286 | 0.1342 | 0.1287 |  | 0.1379
---
(+ 4.16 %)
0.1314 | 0.1344 | 0.1279 | 0.1346 | 0.1313
KL | 0.1324 | 0.1303 | 0.1303 | 0.1346 | 0.1324 | 0.1327 | 0.1311 | 0.1319 | 0.1310 | 0.1298 | 0.1296
ED | SacreBLEU | CE | 2.2794 | 2.4084 | 2.1903 |  | 2.4427
---
(+ 1.42 %)
2.2082 | 2.1922 | 2.2164 | 2.3187 | 2.2622 | 2.3125 | 2.2569
KL | 2.2715 | 2.1682 | 2.2797 | 2.2774 | 2.2781 | 2.3370 | 2.3615 | 2.4319 | 2.2749 | 2.4393 | 2.1431
ROUGE L | CE | 0.1382 | 0.1437 | 0.1373 | 0.1443 | 0.1409 | 0.1425 | 0.1385 | 0.1416 | 0.1398 | 0.1411 | 0.1381
KL | 0.1406 | 0.1395 | 0.1426 | 0.1441 | 0.1454 | 0.1435 | 0.1387 |  | 0.1465
---
(+ 1.95 %)
0.1401 | 0.1430 | 0.1394
METEOR | CE | 0.1254 | 0.1287 | 0.1257 |  | 0.1324
---
(+ 2.08 %)
0.1266 | 0.1248 | 0.1245 | 0.1291 | 0.1266 | 0.1278 | 0.1243
KL | 0.1250 | 0.1233 | 0.1297 | 0.1272 | 0.1302 | 0.1290 | 0.1246 | 0.1323 | 0.1253 | 0.1283 | 0.1234

Table 1: Comparison of sacreBLEU, ROUGE L and METEOR scores using variants of
Cross Entropy (CE) loss and Kullback–Leibler (KL) divergence loss on
DailyDialog (DD) and EmpatheticDialogues (ED) datasets.

##  4 Results and Analysis

We compare the (i) sacreBLEU score Post (2018): a standardised version of the
BLEU score Papineni et al. (2002), (ii) ROUGE L score Lin (2004): which
compares Longest Common Subsequence (LCS), and automatically takes into
account sentence level structure similarity and identifies longest co-
occurring in sequence n-grams, (iii) METEOR score Banerjee and Lavie (2005):
an improvement over BLEU score, which incorporates stemming and synonymy
matching along with exact word matching. Table 1 summarizes the results we
obtained in each of the experiments. The supplementary material contains
results for all the 60 experiments along with additional evaluation metrics
like BERTscore Zhang* et al. (2020) and ROUGE 1 & 2\. In Table 1, the columns
containing “NA” are the baseline results, against which improvements are
measured.
Observations From the experiments we observe that, (i) Using a data dependent
cosine similarity based distribution for label smoothing significantly
outperforms the baseline (vanilla entropy based loss with or without label
smoothing). We observe 12.67 % increase in BLEU score, 0.57 % increase in
ROUGE L score, and 4.16 % increase in METEOR score for the DailyDialog
dataset, and 1.42 % increase in BLEU score, 1.95 % increase in ROUGE L score,
and 2.08 % increase in METEOR score for the EmpatheticDialogues dataset. (ii)
Using additional WordNet synonym based filtering ($w$) does not help
performance. To understand why this is happening, we plotted the distribution
of the smoothing factor $s$ for the randomly selected word “fun”, and observed
that the word had only one overlapping WordNet synonym in our vocabulary:
“play”. This caused the word “play” to be assigned a probability of 0.1, while
all the other words are assigned a probability of 0, except for “fun”, which
was assigned a probability of 0.9. We reason that the sparsity in synonyms
does not help in reducing the overconfidence of the model, as the final
distribution is very similar to non-smoothing targets. Figure 1 illustrates
the probabilities assigned to the incorrect labels of the word “fun”, by each
of the methods discussed in this paper. (iii) Using CE loss instead of KL
generally improves performance while using label smoothing. We reason that
this happens because in case of label smoothing, the constant entropy
coefficient in KL loss reduces the overall loss, thus reducing the gradients
during back propagation, which results in slower learning. (iv) Generally,
using high smoothing value ($s$) does not help in learning. (v) The cosine
similarity threshold $t$ should be treated as a hyperparameter, and will
require tuning depending on the vocabulary of the dataset used. (vi) We also
noticed that a cosine similarity threshold $t$ as high as 0.8 does not help in
learning. We reason that using a high threshold creates a scenario similar to
using WordNet synonyms, where the smoothing probability is distributed among
very few (or no) words. Note that in order to enhance readability, the results
with 0.8 threshold are omitted from Table 1, and are presented in the
additional supplementary materials.

##  5 Conclusion

Label smoothing has an undesirable property of assigning uniform probability
to incorrect labels, which present an incorrect knowledge to learn from. In
this paper we propose ways to convert the uniform distribution to a data
dependent distribution by weighing the smoothing probability using cosine
similarity of word embeddings between the correct and incorrect labels. We
further experiment with WordNet synonyms as an additional filtering criteria,
and report our findings. Using our proposed methodology, we attain significant
improvements over the baseline metrics across all datasets. However, one
drawback that we notice in the proposed system is the inability to factor in
context, while weighing the distribution of the incorrect labels. As future
research, we intend to address this drawback using more contextualised
representations instead of static embeddings.

## References

  * pri (2010) 2010\.  About wordnet.
  * Banerjee and Lavie (2005) Satanjeev Banerjee and Alon Lavie. 2005.  METEOR: An automatic metric for MT evaluation with improved correlation with human judgments.  In _Proceedings of the ACL Workshop on Intrinsic and Extrinsic Evaluation Measures for Machine Translation and/or Summarization_ , pages 65–72, Ann Arbor, Michigan. Association for Computational Linguistics.
  * Cai et al. (2020) Hengyi Cai, Hongshen Chen, Yonghao Song, Cheng Zhang, Xiaofang Zhao, and Dawei Yin. 2020.  Data manipulation: Towards effective instance learning for neural dialogue generation via learning to augment and reweight.  _arXiv preprint arXiv:2004.02594_.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * He and Glass (2020) Tianxing He and James Glass. 2020.  Negative training for neural dialogue response generation.
  * Hinton et al. (2015) Geoffrey Hinton, Oriol Vinyals, and Jeff Dean. 2015.  Distilling the knowledge in a neural network.
  * Jiang et al. (2019) Shaojie Jiang, Pengjie Ren, Christof Monz, and Maarten de Rijke. 2019.  Improving neural response diversity with frequency-aware cross-entropy loss.  _The World Wide Web Conference on - WWW ’19_.
  * Kang and Hashimoto (2020) Daniel Kang and Tatsunori Hashimoto. 2020.  Improved natural language generation via loss truncation.
  * Li et al. (2016) Jiwei Li, Michel Galley, Chris Brockett, Jianfeng Gao, and Bill Dolan. 2016.  A diversity-promoting objective function for neural conversation models.  In _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 110–119, San Diego, California. Association for Computational Linguistics.
  * Li et al. (2017) Yanran Li, Hui Su, Xiaoyu Shen, Wenjie Li, Ziqiang Cao, and Shuzi Niu. 2017.  DailyDialog: A manually labelled multi-turn dialogue dataset.  In _Proceedings of the Eighth International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 986–995, Taipei, Taiwan. Asian Federation of Natural Language Processing.
  * Lin (2004) Chin-Yew Lin. 2004.  ROUGE: A package for automatic evaluation of summaries.  In _Text Summarization Branches Out_ , pages 74–81, Barcelona, Spain. Association for Computational Linguistics.
  * Papineni et al. (2002) Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu. 2002.  Bleu: a method for automatic evaluation of machine translation.  In _Proceedings of the 40th Annual Meeting of the Association for Computational Linguistics_ , pages 311–318, Philadelphia, Pennsylvania, USA. Association for Computational Linguistics.
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher Manning. 2014.  GloVe: Global vectors for word representation.  In _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1532–1543, Doha, Qatar. Association for Computational Linguistics.
  * Pereyra et al. (2017) Gabriel Pereyra, George Tucker, Jan Chorowski, Łukasz Kaiser, and Geoffrey Hinton. 2017.  Regularizing neural networks by penalizing confident output distributions.
  * Post (2018) Matt Post. 2018.  A call for clarity in reporting BLEU scores.  In _Proceedings of the Third Conference on Machine Translation: Research Papers_ , pages 186–191, Belgium, Brussels. Association for Computational Linguistics.
  * Rashkin et al. (2019) Hannah Rashkin, Eric Michael Smith, Margaret Li, and Y-Lan Boureau. 2019.  Towards empathetic open-domain conversation models: A new benchmark and dataset.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 5370–5381, Florence, Italy. Association for Computational Linguistics.
  * Szegedy et al. (2015) Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jonathon Shlens, and Zbigniew Wojna. 2015.  Rethinking the inception architecture for computer vision.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.
  * Wang et al. (2020) Guoyin Wang, Chunyuan Li, Jianqiao Li, Hao Fu, Yuh-Chen Lin, Liqun Chen, Yizhe Zhang, Chenyang Tao, Ruiyi Zhang, Wenlin Wang, Dinghan Shen, Qian Yang, and Lawrence Carin. 2020.  Improving text generation with student-forcing optimal transport.
  * Zhang* et al. (2020) Tianyi Zhang*, Varsha Kishore*, Felix Wu*, Kilian Q. Weinberger, and Yoav Artzi. 2020.  Bertscore: Evaluating text generation with bert.  In _International Conference on Learning Representations_.

##  6 Supplementary Material

###  6.1 All Experiment Results

Table 2 shows the different variants of the baselines that were computed for
both the DailyDialog and EmpatheticDialogues datasets. All performance
improvements are compared against these baselines. For a metric, the best
baseline score among all the hyperparameter settings is chosen to report
improvements. Table 3 shows the results of using different hyperparameter
settings and loss function in the DailyDialog dataset, and Table 4 shows the
results obtained on the EmpatheticDialogues dataset. The best results with
detailed comparison against baselines are already discussed in the main paper.

|  | DailyDialog Dataset | EmpatheticDialogue Dataset
---|---|---|---
|  | s = NA | s = 0.1 | s = 0.2 | s = NA | s = 0.1 | s = 0.2
|  | t = NA | t = NA | t = NA | t = NA | t = NA | t = NA
Metric | Loss | w = NA | w = NA | w = NA | w = NA | w = NA | w = NA
sacreBLEU | CE | 1.6625 | 1.8523 | 1.7251 | 2.2794 | 2.4084 | 2.1903
KL | 1.9469 | 1.7536 | 1.7931 | 2.2715 | 2.1682 | 2.2797
BERTScore | CE | 0.8522 | 0.8520 | 0.8520 | 0.8539 | 0.8544 | 0.8527
KL | 0.8529 | 0.8520 | 0.8510 | 0.8540 | 0.8531 | 0.8541
ROUGE 1 | CE | 0.1272 | 0.1312 | 0.1319 | 0.1536 | 0.1592 | 0.1527
KL | 0.1336 | 0.1298 | 0.1300 | 0.1560 | 0.1545 | 0.1587
ROUGE 2 | CE | 0.0282 | 0.0303 | 0.0299 | 0.0251 | 0.0292 | 0.0251
KL | 0.0305 | 0.0283 | 0.0282 | 0.0267 | 0.0259 | 0.0271
ROUGE L | CE | 0.1209 | 0.1243 | 0.1243 | 0.1382 | 0.1437 | 0.1373
KL | 0.1263 | 0.1223 | 0.1233 | 0.1406 | 0.1395 | 0.1426
METEOR | CE | 0.1244 | 0.1324 | 0.1286 | 0.1254 | 0.1287 | 0.1257
KL | 0.1324 | 0.1303 | 0.1303 | 0.1250 | 0.1233 | 0.1297

Table 2: Baseline results of diverse automatic text generation metrics on the
DailyDialog and EmpatheticDialogues datasets. The hyperparameters s, t and w
control the usage of Label Smoothing, Cosine similarity threshold and WordNet
filtering respectively. For the baseline, t and w were not used, which is
indicated by NA. s = NA signifies vanilla entropy based loss without Label
Smoothing.

###  6.2 Model Training and Parameters

All the models were trained on a single Nvidia V-100 GPU, for 15 epochs each
with a learning rate of 2e-4, batch size of 64, and using AdamW optimizer. The
gradients of the model were clipped with a value of 1, and dropout with
probability 0.1 was applied during training. The average run-time of each
experiment is 60 minutes, with each of the trained models having 17.7 M
parameters. The code, dataset and best performing models are publicly
available through this link: download link.

|  | s = 0.1 | s = 0.2
---|---|---|---
|  | t = 0 | t = 0.5 | t = 0.8 | t = 0 | t = 0.5 | t = 0.8
Metric | Loss | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1
sacreBLEU | CE | 1.9896 | 1.7627 | 2.1936 | 1.8575 | 1.6676 | 1.8859 | 2.1158 | 1.8020 | 2.0536 | 1.9302 | 1.5674 | 1.8502
KL | 1.8459 | 1.8181 | 1.9128 | 1.8858 | 1.7957 | 1.7453 | 1.9387 | 1.8092 | 1.7292 | 1.8856 | 1.5874 | 1.9707
BERTScore | CE | 0.8518 | 0.8529 | 0.8515 | 0.8527 | 0.8507 | 0.8509 | 0.8513 | 0.8525 | 0.8525 | 0.8519 | 0.8507 | 0.8512
KL | 0.8520 | 0.8527 | 0.8517 | 0.8515 | 0.8520 | 0.8516 | 0.8509 | 0.8525 | 0.8522 | 0.8518 | 0.8518 | 0.8515
ROUGE 1 | CE | 0.1309 | 0.1279 | 0.1353 | 0.1326 | 0.1260 | 0.1280 | 0.1298 | 0.1271 | 0.1315 | 0.1317 | 0.1250 | 0.1290
KL | 0.1301 | 0.1332 | 0.1318 | 0.1311 | 0.1281 | 0.1276 | 0.1301 | 0.1328 | 0.1310 | 0.1312 | 0.1263 | 0.1325
ROUGE 2 | CE | 0.0282 | 0.0276 | 0.0309 | 0.0300 | 0.0287 | 0.0280 | 0.0286 | 0.0286 | 0.0308 | 0.0310 | 0.0276 | 0.0305
KL | 0.0283 | 0.0312 | 0.0300 | 0.0294 | 0.0297 | 0.0291 | 0.0285 | 0.0312 | 0.0292 | 0.0299 | 0.0277 | 0.0299
ROUGE L | CE | 0.1238 | 0.1209 | 0.1270 | 0.1260 | 0.1200 | 0.1203 | 0.1217 | 0.1204 | 0.1238 | 0.1244 | 0.1183 | 0.1222
KL | 0.1227 | 0.1264 | 0.1243 | 0.1242 | 0.1213 | 0.1207 | 0.1223 | 0.1253 | 0.1232 | 0.1234 | 0.1185 | 0.1252
METEOR | CE | 0.1342 | 0.1287 | 0.1379 | 0.1314 | 0.1270 | 0.1319 | 0.1344 | 0.1279 | 0.1346 | 0.1313 | 0.1223 | 0.1280
KL | 0.1346 | 0.1324 | 0.1327 | 0.1311 | 0.1262 | 0.1275 | 0.1319 | 0.1310 | 0.1298 | 0.1296 | 0.1247 | 0.1330

Table 3: Results of diverse automatic text generation metrics on the
DailyDialog dataset, trained with variants of Entropy based loss with
different hyperparameter settings: cosine similarity threshold (t), Label
Smoothing (s) and WordNet filtering (w).

|  | s = 0.1 | s = 0.2
---|---|---|---
|  | t = 0 | t = 0.5 | t = 0.8 | t = 0 | t = 0.5 | t = 0.8
Metric | Loss | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1 | w = 0 | w = 1
sacreBLEU | CE | 2.4427 | 2.2082 | 2.1922 | 2.2164 | 2.3467 | 2.2596 | 2.3187 | 2.2622 | 2.3125 | 2.2569 | 2.3944 | 2.2767
KL | 2.2774 | 2.2781 | 2.3370 | 2.3615 | 2.2347 | 2.2769 | 2.4319 | 2.2749 | 2.4393 | 2.1431 | 2.2566 | 2.2652
BERTScore | CE | 0.8543 | 0.8539 | 0.8547 | 0.8528 | 0.8536 | 0.8544 | 0.8544 | 0.8536 | 0.8539 | 0.8532 | 0.8531 | 0.8544
KL | 0.8541 | 0.8543 | 0.8544 | 0.8528 | 0.8536 | 0.8528 | 0.8544 | 0.8528 | 0.8544 | 0.8526 | 0.8535 | 0.8543
ROUGE 1 | CE | 0.1612 | 0.1564 | 0.1577 | 0.1531 | 0.1558 | 0.1551 | 0.1589 | 0.1550 | 0.1575 | 0.1531 | 0.1553 | 0.1590
KL | 0.1594 | 0.1613 | 0.1596 | 0.1540 | 0.1549 | 0.1552 | 0.1619 | 0.1554 | 0.1588 | 0.1545 | 0.1564 | 0.1569
ROUGE 2 | CE | 0.0287 | 0.0270 | 0.0271 | 0.0250 | 0.0274 | 0.0267 | 0.0269 | 0.0265 | 0.0267 | 0.0264 | 0.0261 | 0.0262
KL | 0.0270 | 0.0290 | 0.0273 | 0.0266 | 0.0256 | 0.0253 | 0.0288 | 0.0269 | 0.0274 | 0.0257 | 0.0251 | 0.0268
ROUGE L | CE | 0.1443 | 0.1409 | 0.1425 | 0.1385 | 0.1402 | 0.1388 | 0.1416 | 0.1398 | 0.1411 | 0.1381 | 0.1396 | 0.1423
KL | 0.1441 | 0.1454 | 0.1435 | 0.1387 | 0.1397 | 0.1393 | 0.1465 | 0.1401 | 0.1430 | 0.1394 | 0.1404 | 0.1416
METEOR | CE | 0.1324 | 0.1266 | 0.1248 | 0.1245 | 0.1267 | 0.1264 | 0.1291 | 0.1266 | 0.1278 | 0.1243 | 0.1247 | 0.1292
KL | 0.1272 | 0.1302 | 0.1290 | 0.1246 | 0.1235 | 0.1254 | 0.1323 | 0.1253 | 0.1283 | 0.1234 | 0.1257 | 0.1269

Table 4: Results of diverse automatic text generation metrics on the
EmpatheticDialogues dataset, trained with variants of Entropy based loss with
different hyperparameter settings: cosine similarity threshold (t), Label
Smoothing (s) and WordNet filtering (w).