EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-EP-2020-161 LHCb-PAPER-2020-023 14 September 2020

Observation of multiplicity-dependent prompt $\chi_{c1}(3872)$ and
${\psi{(2S)}}$ production

in $pp$ collisions

LHCb collaboration†††Authors are listed at the end of this Letter.

The production of $\chi_{c1}(3872)$ and ${\psi{(2S)}}$ hadrons is studied as a
function of charged particle multiplicity in $pp$ collisions at a center-of-
mass energy of $8\text{\,Te\kern-1.00006ptV}$, corresponding to an integrated
luminosity of 2$\text{\,fb}^{-1}$. For both states, the fraction that is
produced promptly at the collision vertex is found to decrease as charged
particle multiplicity increases. The ratio of $\chi_{c1}(3872)$ to
${\psi{(2S)}}$ cross-sections for promptly produced particles is also found to
decrease with multiplicity, while no significant dependence on multiplicity is
observed for the equivalent ratio of particles produced away from the
collision vertex in $b$-hadron decays. This behavior is consistent with a
calculation that models the $\chi_{c1}(3872)$ structure as a compact
tetraquark. Implications for the binding energy of the $\chi_{c1}(3872)$ state
are discussed.

Published in Phys. Rev. Lett. 126 (2021) 092001

© 2024 CERN for the benefit of the LHCb collaboration. CC BY 4.0 licence.

In recent years, multiple new resonances containing heavy quarks have been
observed that do not fit into the framework of conventional hadrons, see Ref.
[1] for a recent review. The most studied of these exotic hadrons is the
${\chi_{c1}(3872)}$ state, also known as $X(3872)$. It was first discovered in
the mass spectrum of
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}$ in $B$-meson decays
by the Belle collaboration [2], and has since been confirmed by multiple other
experiments [3, 4, 5, 6]. Despite intense scrutiny, the exact nature of the
${\chi_{c1}(3872)}$ state is still unclear.

Multiple explanations of the $\chi_{c1}(3872)$ structure have been proposed.
Shortly after its discovery, it was considered as one of several possible
charmonium states [7]. However, LHCb has since measured the quantum numbers to
be $J^{PC}=1^{++}$ [8], which disfavors its assignment as conventional
charmonium because no compatible charmonium states with these quantum numbers
are expected near the measured mass [9]. Other models consider the
$\chi_{c1}(3872)$ state to be a tetraquark, which may have further
substructure, composed of a diquark-antidiquark bound state [10, 11, 12] or a
hadrocharmonium state where two light quarks orbit a charmonium core [13].
Mixtures of various exotic and conventional states have also been studied [14,
15, 16, 17]. The remarkable proximity of the $\chi_{c1}(3872)$ mass to the sum
of the ${D}^{0}$ and ${\kern 1.79993pt\overline{\kern-1.79993ptD}}{}^{*0}$
meson masses have led to the consideration of its structure as a hadronic
molecule, a state comprising these two mesons bound via pion exchange [18,
19]. In this case the binding energy of the $\chi_{c1}(3872)$ hadron would be
small, as the mass difference $(M_{{{D}^{0}}}+M_{{{\kern
1.25995pt\overline{\kern-1.25995ptD}}{}^{*0}}})-M_{\chi_{c1}(3872)}=0.07\pm
0.12$$\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ is consistent with zero [20].
Consequently, these models assign the $\chi_{c1}(3872)$ state a large radius
of $O$(10 fm) [17, 21]. Results from recent LHCb studies of the
$\chi_{c1}(3872)$ line shape are compatible with the molecular interpretation
but do not exclude other possibilities [20, 22].

Techniques developed to study quarkonium production in proton-nucleus ($p$A)
collisions can be used to probe the binding energy of hadrons. Measurements of
charmonium production in $p$A collisions at fixed target experiments [23, 24]
and colliders [25, 26, 27, 28, 29, 30] showed that ${\psi{(2S)}}$ production
is suppressed more than ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}$ production in
rapidity regions where a relatively large number of charged particles are
produced. Similarly, measurements of ${\mathchar 28935\relax}$ production at
the Large Hadron Collider (LHC) revealed that the ${\mathchar
28935\relax{(2S)}}$ and ${\mathchar 28935\relax{(3S)}}$ states are suppressed
more than the ${\mathchar 28935\relax{(1S)}}$ [31, 32]. As the effects
governing heavy quark production and transport through the nucleus are assumed
to be similar for states with the same quark content, the mechanism for the
suppression of excited states is expected to occur in the late stages of the
collision, after the heavy quark pair has hadronized into a final state.
Models incorporating final-state effects, such as heavy quark pair breakup via
interactions with comoving hadrons, describe the relative suppression of
excited quarkonium states in $p$A collisions[33, 34, 35, 36, 37]. Similar
final-state effects can also disrupt formation of the $\chi_{c1}(3872)$ state
via interactions with pions produced in the underlying event [38], and would
be especially significant if the $\chi_{c1}(3872)$ structure is a large,
weakly bound hadronic molecule.

High-multiplicity $pp$ collisions provide a hadronic environment that
approaches heavy ion collisions in many respects. Recently, phenomena
typically thought only to occur in collisions of large nuclei have been
observed in high-multiplicity $pp$ collisions, including a near-side ridge in
two-particle angular correlations [39], strangeness enhancement [40], and
collective flow [41]. Multiplicity-dependent modification of ${\mathchar
28935\relax}$ production has also been observed [42]. Therefore, high-
multiplicity $pp$ collisions provide a testing ground for examining final-
state effects observed on quarkonium in $p$A and AA collisions. Measurements
of such effects can provide new constraints on the structure of the
$\chi_{c1}(3872)$ [43].

In this Letter, measurements of the fractions of $\chi_{c1}(3872)$ and
${\psi{(2S)}}$ states, $f_{\mathrm{prompt}}$, that are produced promptly at
the $pp$ collision vertex as a function of charged particle multiplicity are
presented. The $\chi_{c1}(3872)$ and ${\psi{(2S)}}$ states are compared by
measuring the ratio of the $\chi_{c1}(3872)$ to ${\psi{(2S)}}$ cross-sections,
as a function of multiplicity. The $\chi_{c1}(3872)$ and ${\psi{(2S)}}$
candidates are reconstructed through their decays to the
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}$ final state, where
the ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}$ meson subsequently decays to
${\mu^{+}}\mu^{-}$ pairs. This study uses data collected with the LHCb
detector at a center-of-mass energy $\sqrt{s}=8$ TeV, corresponding to an
integrated luminosity of 2$\text{\,fb}^{-1}$.

The LHCb detector [44, 45] is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$, designed for the study of particles
containing $b$ or $c$ quarks. The detector elements comprise a silicon-strip
vertex detector (VELO) surrounding the $pp$ interaction region that allows $b$
hadrons to be identified from their characteristically long flight distance; a
tracking system that provides a measurement of the momentum, $p$, of charged
particles; two ring-imaging Cherenkov (RICH) detectors that discriminate
between different species of charged hadrons, and a series of tracking
detectors interleaved with hadron absorbers for identifying muons. In this
analysis, multiplicity is represented by the number of charged particle tracks
reconstructed in the VELO, $\mathrm{N_{tracks}^{VELO}}$. The VELO track-
reconstruction efficiency has been measured to be about $99\%$ [46].

Simulation is required to model the effects of the detector acceptance and the
imposed selection requirements. In the simulation, $pp$ collisions are
generated using Pythia [47, 48] with a specific LHCb configuration [49].
Decays of unstable particles are described by EvtGen [50]. The interaction of
the generated particles with the detector, and its response, are implemented
using the Geant4 toolkit [51, *Agostinelli:2002hh] as described in Ref. [53].

Events considered in this analysis are selected by a set of triggers designed
to record events containing the decay
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\rightarrow{\mu^{+}}\mu^{-}$. Tracks from
triggered events that are identified as good muon candidates are retained. The
muons are required to have momentum $p>10\text{\,Ge\kern-1.00006ptV\\!/}c$ and
a momentum component transverse to the beam direction
$p_{\mathrm{T}}>650\text{\,Me\kern-1.00006ptV\\!/}c$. Candidate
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}$ mesons are formed from a pair of
oppositely charged muons with an invariant mass within $\pm
39\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ (corresponding to 3 times the
resolution on the mass) of the known ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}$ mass
and combined $p_{\mathrm{T}}>3\text{\,Ge\kern-1.00006ptV\\!/}c$. Charged pion
candidates are selected using particle identification information from the
RICH detectors. They are required to have
$p>3\text{\,Ge\kern-1.00006ptV\\!/}c$ to ensure that the pions are above
threshold in one of the RICH detectors, and have
$p_{\mathrm{T}}>500\text{\,Me\kern-1.00006ptV\\!/}c$ to reduce combinatorial
background.

Selected ${\mu^{+}}\mu^{-}{{\pi}^{+}}{{\pi}^{-}}$ combinations that form a
good-quality common vertex are fitted with kinematic constraints that require
all tracks to originate from a common vertex and constrain the dimuon mass to
the known ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}$ mass [54]. The decay kinematics
are required to satisfy
$M_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}}-M_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}}-M_{{{\pi}^{+}}{{\pi}^{-}}}<300\text{\,Me\kern-1.00006ptV\\!/}c^{2}$
and the candidates must have
$p_{\mathrm{T}}>5\text{\,Ge\kern-1.00006ptV\\!/}c$ and be within the
pseudorapidity range $2<\eta<4.5$. The resulting
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}$ invariant-mass
spectrum is shown in Fig. 1.

To avoid multiplicity biases arising from tracks produced in multiple
collisions that occur in the same beam crossing, only events with a single
reconstructed collision vertex are considered. The position of collision
vertices is restricted to a range along the beam direction
$-60<z<120\text{\,mm}$, to avoid biases from missing tracks that fall outside
the VELO acceptance.

Figure 1: The ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\pi^{+}\pi^{-}$ invariant-
mass spectrum. The inset shows the region of the $\chi_{c1}(3872)$ resonance.

Both $\chi_{c1}(3872)$ and ${\psi{(2S)}}$ hadrons can be produced promptly at
the $pp$ collision vertex, either directly or in strong decays of higher
charmonia states, or in the decays of $b$ hadrons, which travel several
millimeters before decaying. The prompt component of the signal is separated
from the component originating from $b$ decays by performing a simultaneous
fit to the ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}$
invariant-mass spectrum and the pseudo decay-time spectrum. The pseudo decay-
time $t_{z}$ is defined as

$t_{z}\equiv\frac{(z_{\rm decay}-z_{\rm PV})\times M}{p_{z}},$ (1)

where $z_{\rm decay}-z_{\rm PV}$ is the difference between the positions of
the reconstructed vertex of the
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}$ and the collision
vertex along the beam axis, $M$ is the known mass[54] of the reconstructed
${\psi{(2S)}}$ or $\chi_{c1}(3872)$ candidate, and $p_{z}$ is the candidate’s
momentum along the beam axis. The signal in the $t_{z}$ spectrum is fit with a
delta function representing the prompt component and an exponential decay
function representing the component from $b$ decays, which are convolved with
a double Gaussian resolution function. Two different parameterizations of the
$t_{z}$ background components using mass sidebands above and below the mass
peak of interest are employed. The first is an empirically determined
analytical function as was done in Ref. [55] and the second directly uses the
$t_{z}$ shape templates taken from the mass sidebands in the data.

In the fit to the invariant-mass spectrum, the ${\psi{(2S)}}$ peak is
represented by a sum of two Crystal Ball functions, as in a previous LHCb
analysis at $7\text{\,Te\kern-1.00006ptV}$ [56]. The measured
$\chi_{c1}(3872)$ peak is well described by a Gaussian function. The
background contribution is studied by examining the invariant-mass spectrum
constructed by like-sign pion pairs, and is found to be well described by
third-order Chebyshev polynomials; this shape is used to represent the
background when fitting the
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}$ mass spectra. The
invariant mass and $t_{z}$ spectra are divided into bins of
$\mathrm{N_{tracks}^{VELO}}$, and the fit is performed in each bin, separately
for the ${\psi{(2S)}}$ meson and $\chi_{c1}(3872)$ state. An example is shown
in Fig. 2.

Figure 2: The ${\psi{(2S)}}$ (left) invariant-mass and (right) $t_{z}$
spectrum in the $p_{\mathrm{T}}$ and multiplicity ranges
$p_{\mathrm{T}}>5\text{\,Ge\kern-1.00006ptV\\!/}c$ and
$60<\mathrm{N_{tracks}^{VELO}}<80$, with the simultaneous fit superimposed.

The total yield and the measured fraction of the inclusive signal that is
produced at the collision vertex for each state is determined by a fit. Due to
the different production mechanisms, the observed $\chi_{c1}(3872)$ and
${\psi{(2S)}}$ hadrons from these sources have different momentum
distributions, which may lead to differences in acceptance and reconstruction
efficiencies. To account for this effect, the $p_{\mathrm{T}}$ distributions
of prompt and displaced signal candidates are extracted from the data using
the $sPlot$ technique [57], and the $p_{\mathrm{T}}$ distributions of the
simulated particles are weighted to match those of the data. Since these
measurements are binned in multiplicity and effects of multiplicity-dependent
breakup may depend on $p_{\mathrm{T}}$, the simulation is reweighted to match
$p_{\mathrm{T}}$ distributions extracted from low- and high-multiplicity
samples. The difference in the acceptance and reconstruction efficiencies
found using these different parameterizations of the $p_{\mathrm{T}}$
distributions is taken as a systematic uncertainty. Corrections are applied to
account for the relative acceptance of the LHCb spectrometer between particles
produced at the primary vertex and in ${b}$ decays
$\varepsilon_{\mathrm{prompt}}^{\mathrm{acc}}/\varepsilon_{\mathrm{{b}}}^{\mathrm{acc}}$,
which is determined via simulation to be 1.00$\pm$0.01 for the ${\psi{(2S)}}$
state and 1.02$\pm$0.01 for the $\chi_{c1}(3872)$ state, and for the relative
reconstruction and selection efficiencies
$\varepsilon_{\mathrm{prompt}}^{\mathrm{reco}}/\varepsilon_{\mathrm{{b}}}^{\mathrm{reco}}$,
which are 0.99$\pm$0.03 for the ${\psi{(2S)}}$ state and 1.11$\pm$0.04 for the
$\chi_{c1}(3872)$ state. The central value of $f_{\mathrm{prompt}}$ is taken
as the average of the values obtained from the two fitting methods, while the
difference is taken as a systematic uncertainty, which ranges from 1$\%$ to
2$\%$ for the ${\psi{(2S)}}$ fits and from 2$\%$ to 6$\%$ for the
$\chi_{c1}(3872)$ fits. This uncertainty is uncorrelated between bins and is
comparable to the statistical uncertainty on the $\psi(2S)$ data, while the
statistical uncertainty dominates on the $\chi_{c1}(3872)$ data. The
uncertainty on the relative efficiency is taken as a systematic uncertainty on
$f_{\mathrm{prompt}}$, which is correlated between the data points for each
species. The resulting values of $f_{\mathrm{prompt}}$ as a function of
multiplicity are shown in Fig. 3, up to $\mathrm{N_{tracks}^{VELO}}$ = 200.
The fraction of events with $\mathrm{N_{tracks}^{VELO}}>$200 is negligible and
is not included in the analysis. The horizontal position of each point is the
average value of $\mathrm{N_{tracks}^{VELO}}$ for signal events within that
bin.

Figure 3: The fraction, $f_{\mathrm{prompt}}$, of promptly produced
$\chi_{c1}(3872)$ and ${\psi{(2S)}}$ hadrons, as a function of the number of
tracks reconstructed in the VELO. The vertical error bars (boxes) represent
the uncorrelated (correlated) uncertainties, while the horizontal error bars
indicate bin widths.

A clear decrease of $f_{\mathrm{prompt}}$ is seen as the multiplicity
increases, for both the ${\psi{(2S)}}$ and $\chi_{c1}(3872)$ hadrons. This
could be due to a combination of several effects: the average multiplicity is
higher in events containing a ${b}\bar{{b}}$ pair due to their fragmentation
into hadrons and subsequent decays [58, 59]; or the suppression of prompt
${\psi{(2S)}}$ and $\chi_{c1}(3872)$ production via interactions with other
particles produced at the vertex, which decreases the prompt production in
high-multiplicity events, but does not affect production in $b$ decays.

The prompt and $b$-decay components are examined directly by calculating the
ratio of the $\chi_{c1}(3872)$ and ${\psi{(2S)}}$ cross-sections,
$\sigma_{\chi}/\sigma_{\psi}$, times their respective branching fractions to
the ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}$ final state,
${\mathcal{B}}_{\chi}$ and ${\mathcal{B}}_{\psi}$. This ratio is given by

$\displaystyle\frac{\sigma_{\chi}}{\sigma_{\psi}}\frac{{\mathcal{B}}_{\chi}}{{\mathcal{B}}_{\psi}}=\frac{N_{\chi}\,f_{\mathrm{prompt}}^{\chi}}{N_{\psi}\,f_{\mathrm{prompt}}^{\psi}}\frac{\varepsilon_{\psi}^{\mathrm{acc}}}{\varepsilon_{\chi}^{\mathrm{acc}}}\frac{\varepsilon_{\psi}^{\mathrm{reco}}}{\varepsilon_{\chi}^{\mathrm{reco}}}\frac{\varepsilon_{\psi}^{\mathrm{PID}}}{\varepsilon_{\chi}^{\mathrm{PID}}}.$
(2)

Here, $N$ is the signal yield, $f_{\mathrm{prompt}}$ is the prompt fraction
and the $\varepsilon$ terms represent various efficiency corrections of the
corresponding state. The ratio of cross-sections from ${b}$ decays is found by
replacing $f_{\mathrm{prompt}}$ with (1$-f_{\mathrm{prompt}}$) in Eq. 2.
Correlated systematic uncertainties largely cancel in the ratio, and the
result is dominated by uncorrelated uncertainties. The ratio of efficiencies
for four charged decay products to fall within the LHCb acceptance,
$\varepsilon_{\psi}^{\mathrm{acc}}/\varepsilon_{\chi}^{\mathrm{acc}}$, is
found via simulation to be consistent with one with an uncertainty of
approximately 1$\%$ that is determined by varying the $p_{\mathrm{T}}$
distributions of the simulated ${\psi{(2S)}}$ and $\chi_{c1}(3872)$ hadrons.
Control samples of identified muons and pions obtained from data are used to
measure the ratio of muon and pion identification efficiencies,
$\varepsilon_{\psi}^{\mathrm{PID}}/\varepsilon_{\chi}^{\mathrm{PID}}$, which
is near one with an uncertainty of about 1$\%$ due to the finite size of the
control sample. The only relative efficiency that has a significant deviation
from unity is the ratio of reconstruction efficiencies,
$\varepsilon_{\psi}^{\mathrm{reco}}/\varepsilon_{\chi}^{\mathrm{reco}}$, which
is found via simulation to be 0.58$\pm$0.02 (0.65$\pm$0.04) for particles that
are produced promptly (in ${b}$-decays). This is due to the different
kinematic properties of the pion pair produced in the decays: pions from
$\chi_{c1}(3872)$ hadron decays proceed through an intermediate
$\rho^{0}(770)$ resonance [60] and have a higher reconstruction efficiency
than pions from the ${\psi{(2S)}}$ decay due to their higher $p_{T}$. The
uncertainty on the ratio of reconstruction efficiencies is taken from the
variations observed when weighting the $p_{\mathrm{T}}$ distributions of the
simulated ${\psi{(2S)}}$ and $\chi_{c1}(3872)$ hadrons to match those in the
data in different multiplicity bins, as previously discussed.

Figure 4: The ratio of the $\chi_{c1}(3872)$ and ${\psi{(2S)}}$ cross-sections
measured in the ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}$
channel as a function of the number of tracks reconstructed in the VELO. The
point-to-point uncorrelated (correlated) uncertainties are shown as vertical
error bars (boxes) and the bin widths are shown as horizontal error bars. See
text for details on calculations from Ref. [43]

The ratio of cross-sections is shown in Fig. 4. A decrease in the prompt
production of $\chi_{c1}(3872)$ hadrons relative to prompt $\psi(2S)$ mesons
is observed as the charged particle multiplicity increases. To illustrate this
effect, a linear fit to this data, which considers only the uncorrelated
uncertainties, is performed and returns a negative slope that differs from
zero by 5 standard deviations.

After preliminary LHCb results on multiplicity-dependent $\chi_{c1}(3872)$
production were presented [61], calculations of these observables based on the
Comover Interaction Model [34, 35] were performed [43]. In this model,
promptly produced $\chi_{c1}(3872)$ and ${\psi{(2S)}}$ hadrons interact with
other produced particles, with a breakup cross-section $\sigma_{\rm br}$ that
is determined by their radius and binding energy. The model assumes no
interactions at low multiplicity, and the calculations are normalized to the
data in the lowest multiplicity bin. A purely molecular $\chi_{c1}(3872)$ has
a large radius and correspondingly high $\sigma_{\rm br}$, and is quickly
dissociated as multiplicity increases. If coalescence provides an additional
formation mechanism for molecular $\chi_{c1}(3872)$, the ratio
$\sigma_{\chi_{c1}(3872)}/\sigma_{{\psi{(2S)}}}$ rises with multiplicity.
Neither of these calculations are consistent with the data. A compact
tetraquark $\chi_{c1}(3872)$ has a slightly larger radius and $\sigma_{\rm
br}$ than the ${\psi{(2S)}}$, and in this scenario
$\sigma_{\chi_{c1}(3872)}/\sigma_{{\psi{(2S)}}}$ gradually decreases with
multiplicity, matching the measured trend.

In contrast to the prompt data, the ratio of cross-sections for production in
${b}$ decays shows a slight increase, which is not statistically significant.
A linear fit to these data points, again without considering the correlated
systematic uncertainty, gives a positive slope that is consistent with zero
within 1.6 standard deviations. Since these hadrons originate from displaced
decay vertices of ${b}$ hadrons, they are not subject to suppression via
interactions with other particles produced at the primary vertex.
Consequently, this ratio is set only by the branching fractions of ${b}$
decays to $\chi_{c1}(3872)$ and ${\psi{(2S)}}$ hadrons. The multiplicity
dependence of $b$ hadron production has not been studied in detail and
modification of the ${b}$ hadron admixture could affect $\chi_{c1}(3872)$
production, as different ${b}$ hadron species may have different decay
probabilities to $\chi_{c1}(3872)$ states [62, 63]. However, the uncertainties
preclude drawing any firm conclusions on multiplicity-dependent modifications
of ${b}$ hadronization from this data.

In conclusion, the prompt $\chi_{c1}(3872)$ and prompt ${\psi{(2S)}}$
production cross-sections decrease relative to their production via $b$ decays
as the charged particle multiplicity increases in $pp$ collisions at
$8\text{\,Te\kern-1.00006ptV}$. A comparison between the $\chi_{c1}(3872)$ and
${\psi{(2S)}}$ states shows that, in contrast to production from $b$-decays,
which display no significant dependence on multiplicity, prompt production of
$\chi_{c1}(3872)$ is suppressed relative to prompt ${\psi{(2S)}}$ production
as multiplicity increases. This observation is an important ingredient to
obtain a full understanding of the nature of the $\chi_{c1}(3872)$ state.

## Acknowledgements

We express our gratitude to our colleagues in the CERN accelerator departments
for the excellent performance of the LHC. We thank the technical and
administrative staff at the LHCb institutes. We acknowledge support from CERN
and from the national agencies: CAPES, CNPq, FAPERJ and FINEP (Brazil); MOST
and NSFC (China); CNRS/IN2P3 (France); BMBF, DFG and MPG (Germany); INFN
(Italy); NWO (Netherlands); MNiSW and NCN (Poland); MEN/IFA (Romania); MSHE
(Russia); MICINN (Spain); SNSF and SER (Switzerland); NASU (Ukraine); STFC
(United Kingdom); DOE NP and NSF (USA). We acknowledge the computing resources
that are provided by CERN, IN2P3 (France), KIT and DESY (Germany), INFN
(Italy), SURF (Netherlands), PIC (Spain), GridPP (United Kingdom), RRCKI and
Yandex LLC (Russia), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil), PL-
GRID (Poland) and OSC (USA). We are indebted to the communities behind the
multiple open-source software packages on which we depend. Individual groups
or members have received support from AvH Foundation (Germany); EPLANET, Marie
Skłodowska-Curie Actions and ERC (European Union); A*MIDEX, ANR, Labex P2IO
and OCEVU, and Région Auvergne-Rhône-Alpes (France); Key Research Program of
Frontier Sciences of CAS, CAS PIFI, Thousand Talents Program, and Sci. & Tech.
Program of Guangzhou (China); RFBR, RSF and Yandex LLC (Russia); GVA, XuntaGal
and GENCAT (Spain); the Royal Society and the Leverhulme Trust (United
Kingdom).

## References

  * [1] S. L. Olsen, T. Skwarnicki, and D. Zieminska, _Nonstandard heavy mesons and baryons: Experimental evidence_ , Rev. Mod. Phys. 90 (2018) 015003, arXiv:1708.04012
  * [2] Belle collaboration, S. K. Choi et al., _Observation of a narrow charmonium - like state in exclusive ${{B}^{\pm}}\\!\rightarrow{{K}^{\pm}}{{\pi}^{+}}{{\pi}^{-}}{{J\mskip-3.0mu/\mskip-2.0mu\psi}}$ decays_, Phys. Rev. Lett. 91 (2003) 262001, arXiv:hep-ex/0309032
  * [3] CDF collaboration, D. Acosta et al., _Observation of the narrow state $X(3872)\rightarrow J/\psi\pi^{+}\pi^{-}$ in $\bar{p}p$ collisions at $\sqrt{s}=1.96\text{\,Te\kern-1.00006ptV}$_, Phys. Rev. Lett. 93 (2004) 072001, arXiv:hep-ex/0312021
  * [4] D0 collaboration, V. M. Abazov et al., _Observation and properties of the $X(3872)$ decaying to $J/\psi\pi^{+}\pi^{-}$ in $p\bar{p}$ collisions at $\sqrt{s}=1.96\text{\,Te\kern-1.00006ptV}$_, Phys. Rev. Lett. 93 (2004) 162002, arXiv:hep-ex/0405004
  * [5] BaBar collaboration, B. Aubert et al., _Study of the $B\rightarrow J/\psi K^{-}\pi^{+}\pi^{-}$ decay and measurement of the $B\rightarrow X(3872)K^{-}$ branching fraction_, Phys. Rev. D71 (2005) 071103, arXiv:hep-ex/0406022
  * [6] LHCb collaboration, R. Aaij et al., _Observation of $X(3872)$ production in $p$ $p$ collisions at $\sqrt{s}=$7 TeV_, Eur. Phys. J. C72 (2012) 1972, arXiv:1112.5310
  * [7] E. J. Eichten, K. Lane, and C. Quigg, _Charmonium levels near threshold and the narrow state $X(3872)\rightarrow\pi^{+}\pi^{-}J/\psi$_, Phys. Rev. D69 (2004) 094019, arXiv:hep-ph/0401210
  * [8] LHCb collaboration, R. Aaij et al., _Determination of the $X(3872)$ meson quantum numbers_, Phys. Rev. Lett. 110 (2013) 222001, arXiv:1302.6269
  * [9] T. Barnes, S. Godfrey, and E. S. Swanson, _Higher charmonia_ , Phys. Rev. D72 (2005) 054026, arXiv:hep-ph/0505002
  * [10] L. Maiani, F. Piccinini, A. D. Polosa, and V. Riquer, _Diquark-antidiquarks with hidden or open charm and the nature of X(3872)_ , Phys. Rev. D71 (2005) 014028, arXiv:hep-ph/0412098
  * [11] R. D. Matheus, S. Narison, M. Nielsen, and J. M. Richard, _Can the X(3872) be a $1^{++}$ four-quark state?_, Phys. Rev. D75 (2007) 014005, arXiv:hep-ph/0608297
  * [12] S. Dubnicka, A. Z. Dubnickova, M. A. Ivanov, and J. G. Korner, _Quark model description of the tetraquark state X(3872) in a relativistic constituent quark model with infrared confinement_ , Phys. Rev. D81 (2010) 114007, arXiv:1004.1291
  * [13] P. C. Wallbott, G. Eichmann, and C. S. Fischer, _$X(3872)$ as a four-quark state in a Dyson-Schwinger/Bethe-Salpeter approach_, Phys. Rev. D100 (2019) 014033, arXiv:1905.02615
  * [14] P. Colangelo, F. De Fazio, and S. Nicotri, _$X(3872)\rightarrow D\bar{D}\gamma$ decays and the structure of $X(3872)$_, Phys. Lett. B 650 (2007) 166, arXiv:hep-ph/0701052
  * [15] R. D. Matheus, F. S. Navarra, M. Nielsen, and C. M. Zanetti, _QCD Sum Rules for the X(3872) as a mixed molecule-charmoniun state_ , Phys. Rev. D80 (2009) 056002, arXiv:0907.2683
  * [16] M. Butenschoen, Z.-G. He, and B. A. Kniehl, _Deciphering the $X(3872)$ via its polarization in prompt production at the CERN LHC_, Phys. Rev. Lett. 123 (2019) 032001, arXiv:1906.08553
  * [17] S. Coito, G. Rupp, and E. van Beveren, _X(3872) is not a true molecule_ , Eur. Phys. J. C73 (2013) 2351, arXiv:1212.0648
  * [18] N. A. Tornqvist, _Isospin breaking of the narrow charmonium state of Belle at 3872-MeV as a deuson_ , Phys. Lett. B590 (2004) 209, arXiv:hep-ph/0402237
  * [19] E. Braaten and M. Kusunoki, _Exclusive production of the X(3872) in B meson decay_ , Phys. Rev. D71 (2005) 074005, arXiv:hep-ph/0412268
  * [20] LHCb collaboration, R. Aaij et al., _Study of the $\psi_{2}(3823)$ and ${\chi_{c}}(3872)$ states in $B^{+}\rightarrow({{J\mskip-3.0mu/\mskip-2.0mu\psi}}\pi^{+}\pi^{-})K^{+}$ decays_, JHEP 08 (2020) 123, arXiv:2005.13422
  * [21] M. Cardoso, G. Rupp, and E. van Beveren, _Unquenched quark-model calculation of $X(3872)$ electromagnetic decays_, Eur. Phys. J. C75 (2015) 26, arXiv:1411.1654
  * [22] LHCb collaboration, R. Aaij et al., _Study of the line shape of the ${\chi_{c}}(3872)$ meson_, arXiv:2005.13419, submitted to Phys. Rev. D.
  * [23] NuSea collaboration, M. J. Leitch et al., _Measurement of differences between $J/\psi$ and $\psi^{\prime}$ suppression in p-A collisions at 800-GeV/c_, Phys. Rev. Lett. 84 (2000) 3256, arXiv:nucl-ex/9909007
  * [24] NA50 collaboration, B. Alessandro et al., _J/psi and psi-prime production and their normal nuclear absorption in proton-nucleus collisions at 400-GeV_ , Eur. Phys. J. C48 (2006) 329, arXiv:nucl-ex/0612012
  * [25] PHENIX collaboration, A. Adare et al., _Nuclear Modification of $\psi^{\prime},\chi_{c}$, and J/$\psi$ Production in d+Au Collisions at $\sqrt{s_{\scriptscriptstyle\text{NN}}}$ = 200 GeV_, Phys. Rev. Lett. 111 (2013) 202301, arXiv:1305.5516
  * [26] ALICE collaboration, B. B. Abelev et al., _Suppression of $\psi$(2S) production in p-Pb collisions at $\sqrt{s_{\scriptscriptstyle\text{NN}}}$ = 5.02 TeV_, JHEP 12 (2014) 073, arXiv:1405.3796
  * [27] LHCb collaboration, R. Aaij et al., _Study of $\psi{(2S)}$ production cross-sections and cold nuclear matter effects in $p$Pb collisions at $\sqrt{s_{\scriptscriptstyle\text{NN}}}=$5 TeV_, JHEP 03 (2016) 133, arXiv:1601.07878
  * [28] PHENIX collaboration, A. Adare et al., _Measurement of the relative yields of $\psi(2S)$ to $\psi(1S)$ mesons produced at forward and backward rapidity in $p$$+$$p$, $p$$+$Al, $p$$+$Au, and 3He$+$Au collisions at $\sqrt{s_{\scriptscriptstyle\text{NN}}}$ = 200 GeV_, Phys. Rev. C95 (2017) 034904, arXiv:1609.06550
  * [29] CMS collaboration, A. M. Sirunyan et al., _Measurement of prompt $\psi$(2S) production cross sections in proton-lead and proton-proton collisions at $\sqrt{s_{\scriptscriptstyle\text{NN}}}=$ 5.02 TeV_, Phys. Lett. B790 (2019) 509, arXiv:1805.02248
  * [30] ALICE collaboration, S. Acharya et al., _Measurement of nuclear effects on $\psi\rm{(2S)}$ production in p-Pb collisions at $\sqrt{\textit{s}_{\rm NN}}=8.16$ TeV_, JHEP 07 (2020) 237, arXiv:2003.06053
  * [31] CMS collaboration, S. Chatrchyan et al., _Event Activity Dependence of Y(nS) Production in $\sqrt{s_{\scriptscriptstyle\text{NN}}}$=5.02 TeV pPb and $\sqrt{s}$=2.76 TeV pp Collisions_, JHEP 04 (2014) 103, arXiv:1312.6300
  * [32] LHCb collaboration, R. Aaij et al., _Study of $\mathchar 28935\relax$ production in $p$Pb collisions at $\sqrt{s_{\scriptscriptstyle\text{NN}}}=$8.16 TeV_, JHEP 11 (2018) 194, arXiv:1810.07655
  * [33] A. Capella, A. Kaidalov, A. Kouider Akil, and C. Gerschel, _$J/\psi$ and $\psi^{\prime}$ suppression in heavy ion collisions_, Phys. Lett. B393 (1997) 431, arXiv:hep-ph/9607265
  * [34] A. Capella, E. G. Ferreiro, and A. B. Kaidalov, _Nonsaturation of the ${J\mskip-3.0mu/\mskip-2.0mu\psi}$ suppression at large transverse energy in the comovers approach_, Phys. Rev. Lett. 85 (2000) 2080, arXiv:hep-ph/0002300
  * [35] E. G. Ferreiro, _Excited charmonium suppression in proton–nucleus collisions as a consequence of comovers_ , Phys. Lett. B749 (2015) 98, arXiv:1411.0549
  * [36] Y.-Q. Ma, R. Venugopalan, K. Watanabe, and H.-F. Zhang, _$\psi(2S)$ versus $J/\psi$ suppression in proton-nucleus collisions from factorization violating soft color exchanges_, Phys. Rev. C97 (2018) 014909, arXiv:1707.07266
  * [37] E. G. Ferreiro and J.-P. Lansberg, _Is bottomonium suppression in proton-nucleus and nucleus-nucleus collisions at LHC energies due to the same effects?_ , JHEP 10 (2018) 094, Erratum ibid. 03 (2019) 063, arXiv:1804.04474
  * [38] E. Braaten, H.-W. Hammer, and T. Mehen, _Scattering of an Ultrasoft Pion and the X(3872)_ , Phys. Rev. D82 (2010) 034018, arXiv:1005.1688
  * [39] CMS collaboration, V. Khachatryan et al., _Observation of Long-Range Near-Side Angular Correlations in Proton-Proton Collisions at the LHC_ , JHEP 09 (2010) 091, arXiv:1009.4122
  * [40] ALICE collaboration, J. Adam et al., _Enhanced production of multi-strange hadrons in high-multiplicity proton-proton collisions_ , Nature Phys. 13 (2017) 535, arXiv:1606.07424
  * [41] CMS collaboration, V. Khachatryan et al., _Evidence for collectivity in pp collisions at the LHC_ , Phys. Lett. B765 (2017) 193, arXiv:1606.06198
  * [42] CMS collaboration, A. M. Sirunyan et al., _Investigation into the event-activity dependence of $\Upsilon$(nS) relative production in proton-proton collisions at $\sqrt{s}=$ 7 TeV_, arXiv:2007.04277
  * [43] A. Esposito et al., _The nature of $X(3872)$ from high-multiplicity $pp$ collisions_, arXiv:2006.15044
  * [44] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [45] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [46] R. Aaij et al., _Performance of the LHCb Vertex Locator_ , JINST 9 (2014) P09007, arXiv:1405.7808
  * [47] T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * [48] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [49] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , J. Phys. Conf. Ser. 331 (2011) 032047
  * [50] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [51] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [52] Geant4 collaboration, S. Agostinelli et al., _Geant4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [53] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [54] Particle Data Group, M. Tanabashi et al., _Review of particle physics_ , Phys. Rev. D98 (2018) 030001
  * [55] LHCb collaboration, R. Aaij et al., _Measurement of $\psi{(2S)}$ production cross-sections in proton-proton collisions at $\sqrt{s}=$7 and 13 TeV_, Eur. Phys. J. C80 (2020) 185, arXiv:1908.03099
  * [56] LHCb collaboration, R. Aaij et al., _Measurement of $\psi{(2S)}$ meson production in $p$ $p$ collisions at $\sqrt{s}=$7 TeV_, Eur. Phys. J. C72 (2012) 2100, Erratum ibid. C80 (2020) 49, arXiv:1204.1258
  * [57] M. Pivk and F. R. Le Diberder, _sPlot: A statistical tool to unfold data distributions_ , Nucl. Instrum. Meth. A555 (2005) 356, arXiv:physics/0402083
  * [58] P. C. Rowson et al., _Charged Multiplicity of Hadronic Events Containing Heavy Quark Jets_ , Phys. Rev. Lett. 54 (1985) 2580
  * [59] OPAL collaboration, G. Abbiendi et al., _Charged particle multiplicities in heavy and light quark initiated events above the Z0 peak_ , Phys. Lett. B550 (2002) 33, arXiv:hep-ex/0211007
  * [60] CDF collaboration, A. Abulencia et al., _Measurement of the dipion mass spectrum in $X(3872)\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{\pi}^{+}}{{\pi}^{-}}$ decays._, Phys. Rev. Lett. 96 (2006) 102002, arXiv:hep-ex/0512074
  * [61] J. Matthew Durham for the LHCb collaboration, _LHCb measurements of the exotic tetraquark candidate $\chi_{c1}(3872)$ in high-multiplicity $pp$ and $p$Pb collisions_, in 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions, 2020, arXiv:2002.01551
  * [62] ATLAS collaboration, M. Aaboud et al., _Measurements of $\psi(2S)$ and $X(3872)\rightarrow J/\psi\pi^{+}\pi^{-}$ production in $pp$ collisions at $\sqrt{s}=8$ TeV with the ATLAS detector_, JHEP 01 (2017) 117, arXiv:1610.09303
  * [63] CMS collaboration, A. M. Sirunyan et al., _Observation of the B ${}^{0}_{\mathrm{s}}\rightarrow$X(3872)$\phi$ decay_, arXiv:2005.04764

LHCb collaboration

R. Aaij31, C. Abellán Beteta49, T. Ackernley59, B. Adeva45, M. Adinolfi53, H.
Afsharnia9, C.A. Aidala84, S. Aiola25, Z. Ajaltouni9, S. Akar64, J.
Albrecht14, F. Alessio47, M. Alexander58, A. Alfonso Albero44, Z. Aliouche61,
G. Alkhazov37, P. Alvarez Cartelle47, S. Amato2, Y. Amhis11, L. An21, L.
Anderlini21, A. Andreianov37, M. Andreotti20, F. Archilli16, A. Artamonov43,
M. Artuso67, K. Arzymatov41, E. Aslanides10, M. Atzeni49, B. Audurier11, S.
Bachmann16, M. Bachmayer48, J.J. Back55, S. Baker60, P. Baladron Rodriguez45,
V. Balagura11, W. Baldini20, J. Baptista Leite1, R.J. Barlow61, S. Barsuk11,
W. Barter60, M. Bartolini23,h, F. Baryshnikov80, J.M. Basels13, G. Bassi28, B.
Batsukh67, A. Battig14, A. Bay48, M. Becker14, F. Bedeschi28, I. Bediaga1, A.
Beiter67, V. Belavin41, S. Belin26, V. Bellee48, K. Belous43, I. Belov39, I.
Belyaev38, G. Bencivenni22, E. Ben-Haim12, A. Berezhnoy39, R. Bernet49, D.
Berninghoff16, H.C. Bernstein67, C. Bertella47, E. Bertholet12, A. Bertolin27,
C. Betancourt49, F. Betti19,e, M.O. Bettler54, Ia. Bezshyiko49, S. Bhasin53,
J. Bhom33, L. Bian72, M.S. Bieker14, S. Bifani52, P. Billoir12, M. Birch60,
F.C.R. Bishop54, A. Bizzeti21,r, M. Bjørn62, M.P. Blago47, T. Blake55, F.
Blanc48, S. Blusk67, D. Bobulska58, J.A. Boelhauve14, O. Boente Garcia45, T.
Boettcher63, A. Boldyrev81, A. Bondar42,u, N. Bondar37, S. Borghi61, M.
Borisyak41, M. Borsato16, J.T. Borsuk33, S.A. Bouchiba48, T.J.V. Bowcock59, A.
Boyer47, C. Bozzi20, M.J. Bradley60, S. Braun65, A. Brea Rodriguez45, M.
Brodski47, J. Brodzicka33, A. Brossa Gonzalo55, D. Brundu26, A. Buonaura49, C.
Burr47, A. Bursche26, A. Butkevich40, J.S. Butter31, J. Buytaert47, W.
Byczynski47, S. Cadeddu26, H. Cai72, R. Calabrese20,g, L. Calefice14, L.
Calero Diaz22, S. Cali22, R. Calladine52, M. Calvi24,i, M. Calvo Gomez83, P.
Camargo Magalhaes53, A. Camboni44, P. Campana22, D.H. Campora Perez47, A.F.
Campoverde Quezada5, S. Capelli24,i, L. Capriotti19,e, A. Carbone19,e, G.
Carboni29, R. Cardinale23,h, A. Cardini26, I. Carli6, P. Carniti24,i, K.
Carvalho Akiba31, A. Casais Vidal45, G. Casse59, M. Cattaneo47, G.
Cavallero47, S. Celani48, J. Cerasoli10, A.J. Chadwick59, M.G. Chapman53, M.
Charles12, Ph. Charpentier47, G. Chatzikonstantinidis52, C.A. Chavez
Barajas59, M. Chefdeville8, C. Chen3, S. Chen26, A. Chernov33, S.-G. Chitic47,
V. Chobanova45, S. Cholak48, M. Chrzaszcz33, A. Chubykin37, V. Chulikov37, P.
Ciambrone22, M.F. Cicala55, X. Cid Vidal45, G. Ciezarek47, P.E.L. Clarke57, M.
Clemencic47, H.V. Cliff54, J. Closier47, J.L. Cobbledick61, V. Coco47, J.A.B.
Coelho11, J. Cogan10, E. Cogneras9, L. Cojocariu36, P. Collins47, T.
Colombo47, L. Congedo18, A. Contu26, N. Cooke52, G. Coombs58, G. Corti47, C.M.
Costa Sobral55, B. Couturier47, D.C. Craik63, J. Crkovská66, M. Cruz Torres1,
R. Currie57, C.L. Da Silva66, E. Dall’Occo14, J. Dalseno45, C. D’Ambrosio47,
A. Danilina38, P. d’Argent47, A. Davis61, O. De Aguiar Francisco61, K. De
Bruyn77, S. De Capua61, M. De Cian48, J.M. De Miranda1, L. De Paula2, M. De
Serio18,d, D. De Simone49, P. De Simone22, J.A. de Vries78, C.T. Dean66, W.
Dean84, D. Decamp8, L. Del Buono12, B. Delaney54, H.-P. Dembinski14, A.
Dendek34, V. Denysenko49, D. Derkach81, O. Deschamps9, F. Desse11, F.
Dettori26,f, B. Dey72, P. Di Nezza22, S. Didenko80, L. Dieste Maronas45, H.
Dijkstra47, V. Dobishuk51, A.M. Donohoe17, F. Dordei26, A.C. dos Reis1, L.
Douglas58, A. Dovbnya50, A.G. Downes8, K. Dreimanis59, M.W. Dudek33, L.
Dufour47, V. Duk76, P. Durante47, J.M. Durham66, D. Dutta61, M. Dziewiecki16,
A. Dziurda33, A. Dzyuba37, S. Easo56, U. Egede68, V. Egorychev38, S.
Eidelman42,u, S. Eisenhardt57, S. Ek-In48, L. Eklund58, S. Ely67, A. Ene36, E.
Epple66, S. Escher13, J. Eschle49, S. Esen31, T. Evans47, A. Falabella19, J.
Fan3, Y. Fan5, B. Fang72, N. Farley52, S. Farry59, D. Fazzini24,i, P. Fedin38,
M. Féo47, P. Fernandez Declara47, A. Fernandez Prieto45, F. Ferrari19,e, L.
Ferreira Lopes48, F. Ferreira Rodrigues2, S. Ferreres Sole31, M. Ferrillo49,
M. Ferro-Luzzi47, S. Filippov40, R.A. Fini18, M. Fiorini20,g, M. Firlej34,
K.M. Fischer62, C. Fitzpatrick61, T. Fiutowski34, F. Fleuret11,b, M.
Fontana47, F. Fontanelli23,h, R. Forty47, V. Franco Lima59, M. Franco
Sevilla65, M. Frank47, E. Franzoso20, G. Frau16, C. Frei47, D.A. Friday58, J.
Fu25, Q. Fuehring14, W. Funk47, E. Gabriel31, T. Gaintseva41, A. Gallas
Torreira45, D. Galli19,e, S. Gambetta57, Y. Gan3, M. Gandelman2, P. Gandini25,
Y. Gao4, M. Garau26, L.M. Garcia Martin55, P. Garcia Moreno44, J. García
Pardiñas49, B. Garcia Plana45, F.A. Garcia Rosales11, L. Garrido44, D.
Gascon44, C. Gaspar47, R.E. Geertsema31, D. Gerick16, L.L. Gerken14, E.
Gersabeck61, M. Gersabeck61, T. Gershon55, D. Gerstel10, Ph. Ghez8, V.
Gibson54, M. Giovannetti22,j, A. Gioventù45, P. Gironella Gironell44, L.
Giubega36, C. Giugliano20,g, K. Gizdov57, E.L. Gkougkousis47, V.V. Gligorov12,
C. Göbel69, E. Golobardes83, D. Golubkov38, A. Golutvin60,80, A. Gomes1,a, S.
Gomez Fernandez44, M. Goncerz33, G. Gong3, P. Gorbounov38, I.V. Gorelov39, C.
Gotti24, E. Govorkova31, J.P. Grabowski16, R. Graciani Diaz44, T.
Grammatico12, L.A. Granado Cardoso47, E. Graugés44, E. Graverini48, G.
Graziani21, A. Grecu36, L.M. Greeven31, P. Griffith20, L. Grillo61, S.
Gromov80, L. Gruber47, B.R. Gruberg Cazon62, C. Gu3, M. Guarise20, P. A.
Günther16, E. Gushchin40, A. Guth13, Y. Guz43,47, T. Gys47, T. Hadavizadeh68,
G. Haefeli48, C. Haen47, J. Haimberger47, S.C. Haines54, T. Halewood-leagas59,
P.M. Hamilton65, Q. Han7, X. Han16, T.H. Hancock62, S. Hansmann-Menzemer16, N.
Harnew62, T. Harrison59, C. Hasse47, M. Hatch47, J. He5, M. Hecker60, K.
Heijhoff31, K. Heinicke14, A.M. Hennequin47, K. Hennessy59, L. Henry25,46, J.
Heuel13, A. Hicheur2, D. Hill62, M. Hilton61, S.E. Hollitt14, P.H. Hopchev48,
J. Hu16, J. Hu71, W. Hu7, W. Huang5, X. Huang72, W. Hulsbergen31, R.J.
Hunter55, M. Hushchyn81, D. Hutchcroft59, D. Hynds31, P. Ibis14, M. Idzik34,
D. Ilin37, P. Ilten52, A. Inglessi37, A. Ishteev80, K. Ivshin37, R.
Jacobsson47, S. Jakobsen47, E. Jans31, B.K. Jashal46, A. Jawahery65, V.
Jevtic14, M. Jezabek33, F. Jiang3, M. John62, D. Johnson47, C.R. Jones54, T.P.
Jones55, B. Jost47, N. Jurik47, S. Kandybei50, Y. Kang3, M. Karacson47, J.M.
Kariuki53, N. Kazeev81, M. Kecke16, F. Keizer54,47, M. Kenzie55, T. Ketel32,
B. Khanji47, A. Kharisova82, S. Kholodenko43, K.E. Kim67, T. Kirn13, V.S.
Kirsebom48, O. Kitouni63, S. Klaver31, K. Klimaszewski35, S. Koliiev51, A.
Kondybayeva80, A. Konoplyannikov38, P. Kopciewicz34, R. Kopecna16, P.
Koppenburg31, M. Korolev39, I. Kostiuk31,51, O. Kot51, S. Kotriakhova37,30, P.
Kravchenko37, L. Kravchuk40, R.D. Krawczyk47, M. Kreps55, F. Kress60, S.
Kretzschmar13, P. Krokovny42,u, W. Krupa34, W. Krzemien35, W. Kucewicz85,33,k,
M. Kucharczyk33, V. Kudryavtsev42,u, H.S. Kuindersma31, G.J. Kunde66, T.
Kvaratskheliya38, D. Lacarrere47, G. Lafferty61, A. Lai26, A. Lampis26, D.
Lancierini49, J.J. Lane61, R. Lane53, G. Lanfranchi22, C. Langenbruch13, J.
Langer14, O. Lantwin49,80, T. Latham55, F. Lazzari28,s, R. Le Gac10, S.H.
Lee84, R. Lefèvre9, A. Leflat39, S. Legotin80, O. Leroy10, T. Lesiak33, B.
Leverington16, H. Li71, L. Li62, P. Li16, X. Li66, Y. Li6, Y. Li6, Z. Li67, X.
Liang67, T. Lin60, R. Lindner47, V. Lisovskyi14, R. Litvinov26, G. Liu71, H.
Liu5, S. Liu6, X. Liu3, A. Loi26, J. Lomba Castro45, I. Longstaff58, J.H.
Lopes2, G. Loustau49, G.H. Lovell54, Y. Lu6, D. Lucchesi27,l, S. Luchuk40, M.
Lucio Martinez31, V. Lukashenko31, Y. Luo3, A. Lupato61, E. Luppi20,g, O.
Lupton55, A. Lusiani28,q, X. Lyu5, L. Ma6, S. Maccolini19,e, F. Machefert11,
F. Maciuc36, V. Macko48, P. Mackowiak14, S. Maddrell-Mander53, O. Madejczyk34,
L.R. Madhan Mohan53, O. Maev37, A. Maevskiy81, D. Maisuzenko37, M.W.
Majewski34, S. Malde62, B. Malecki47, A. Malinin79, T. Maltsev42,u, H.
Malygina16, G. Manca26,f, G. Mancinelli10, R. Manera Escalero44, D.
Manuzzi19,e, D. Marangotto25,n, J. Maratas9,t, J.F. Marchand8, U. Marconi19,
S. Mariani21,47,21, C. Marin Benito11, M. Marinangeli48, P. Marino48, J.
Marks16, P.J. Marshall59, G. Martellotti30, L. Martinazzoli47, M.
Martinelli24,i, D. Martinez Santos45, F. Martinez Vidal46, A. Massafferri1, M.
Materok13, R. Matev47, A. Mathad49, Z. Mathe47, V. Matiunin38, C. Matteuzzi24,
K.R. Mattioli84, A. Mauri31, E. Maurice11,b, J. Mauricio44, M. Mazurek35, M.
McCann60, L. Mcconnell17, T.H. Mcgrath61, A. McNab61, R. McNulty17, J.V.
Mead59, B. Meadows64, C. Meaux10, G. Meier14, N. Meinert75, D. Melnychuk35, S.
Meloni24,i, M. Merk31,78, A. Merli25, L. Meyer Garcia2, M. Mikhasenko47, D.A.
Milanes73, E. Millard55, M.-N. Minard8, L. Minzoni20,g, S.E. Mitchell57, B.
Mitreska61, D.S. Mitzel47, A. Mödden14, R.A. Mohammed62, R.D. Moise60, T.
Mombächer14, I.A. Monroy73, S. Monteil9, M. Morandin27, G. Morello22, M.J.
Morello28,q, J. Moron34, A.B. Morris74, A.G. Morris55, R. Mountain67, H. Mu3,
F. Muheim57, M. Mukherjee7, M. Mulder47, D. Müller47, K. Müller49, C.H.
Murphy62, D. Murray61, P. Muzzetto26, P. Naik53, T. Nakada48, R. Nandakumar56,
T. Nanut48, I. Nasteva2, M. Needham57, I. Neri20,g, N. Neri25,n, S. Neubert74,
N. Neufeld47, R. Newcombe60, T.D. Nguyen48, C. Nguyen-Mau48, E.M. Niel11, S.
Nieswand13, N. Nikitin39, N.S. Nolte47, C. Nunez84, A. Oblakowska-Mucha34, V.
Obraztsov43, D.P. O’Hanlon53, R. Oldeman26,f, C.J.G. Onderwater77, A.
Ossowska33, J.M. Otalora Goicochea2, T. Ovsiannikova38, P. Owen49, A.
Oyanguren46, B. Pagare55, P.R. Pais47, T. Pajero28,47,q, A. Palano18, M.
Palutan22, Y. Pan61, G. Panshin82, A. Papanestis56, M. Pappagallo18,d, L.L.
Pappalardo20,g, C. Pappenheimer64, W. Parker65, C. Parkes61, C.J. Parkinson45,
B. Passalacqua20, G. Passaleva21, A. Pastore18, M. Patel60, C. Patrignani19,e,
C.J. Pawley78, A. Pearce47, A. Pellegrino31, M. Pepe Altarelli47, S.
Perazzini19, D. Pereima38, P. Perret9, K. Petridis53, A. Petrolini23,h, A.
Petrov79, S. Petrucci57, M. Petruzzo25, A. Philippov41, L. Pica28, M.
Piccini76, B. Pietrzyk8, G. Pietrzyk48, M. Pili62, D. Pinci30, J. Pinzino47,
F. Pisani47, A. Piucci16, Resmi P.K10, V. Placinta36, S. Playfer57, J.
Plews52, M. Plo Casasus45, F. Polci12, M. Poli Lener22, M. Poliakova67, A.
Poluektov10, N. Polukhina80,c, I. Polyakov67, E. Polycarpo2, G.J. Pomery53, S.
Ponce47, A. Popov43, D. Popov5,47, S. Popov41, S. Poslavskii43, K. Prasanth33,
L. Promberger47, C. Prouve45, V. Pugatch51, A. Puig Navarro49, H. Pullen62, G.
Punzi28,m, W. Qian5, J. Qin5, R. Quagliani12, B. Quintana8, N.V. Raab17, R.I.
Rabadan Trejo10, B. Rachwal34, J.H. Rademacker53, M. Rama28, M. Ramos
Pernas55, M.S. Rangel2, F. Ratnikov41,81, G. Raven32, M. Reboud8, F. Redi48,
F. Reiss12, C. Remon Alepuz46, Z. Ren3, V. Renaudin62, R. Ribatti28, S.
Ricciardi56, D.S. Richards56, K. Rinnert59, P. Robbe11, A. Robert12, G.
Robertson57, A.B. Rodrigues48, E. Rodrigues59, J.A. Rodriguez Lopez73, A.
Rollings62, P. Roloff47, V. Romanovskiy43, M. Romero Lamas45, A. Romero
Vidal45, J.D. Roth84, M. Rotondo22, M.S. Rudolph67, T. Ruf47, J. Ruiz Vidal46,
A. Ryzhikov81, J. Ryzka34, J.J. Saborido Silva45, N. Sagidova37, N. Sahoo55,
B. Saitta26,f, D. Sanchez Gonzalo44, C. Sanchez Gras31, C. Sanchez
Mayordomo46, R. Santacesaria30, C. Santamarina Rios45, M. Santimaria22, E.
Santovetti29,j, D. Saranin80, G. Sarpis61, M. Sarpis74, A. Sarti30, C.
Satriano30,p, A. Satta29, M. Saur5, D. Savrina38,39, H. Sazak9, L.G.
Scantlebury Smead62, S. Schael13, M. Schellenberg14, M. Schiller58, H.
Schindler47, M. Schmelling15, T. Schmelzer14, B. Schmidt47, O. Schneider48, A.
Schopper47, M. Schubiger31, S. Schulte48, M.H. Schune11, R. Schwemmer47, B.
Sciascia22, A. Sciubba30, S. Sellam45, A. Semennikov38, M. Senghi Soares32, A.
Sergi52,47, N. Serra49, J. Serrano10, L. Sestini27, A. Seuthe14, P. Seyfert47,
D.M. Shangase84, M. Shapkin43, I. Shchemerov80, L. Shchutska48, T. Shears59,
L. Shekhtman42,u, V. Shevchenko79, E.B. Shields24,i, E. Shmanin80, J.D.
Shupperd67, B.G. Siddi20, R. Silva Coutinho49, G. Simi27, S. Simone18,d, I.
Skiba20,g, N. Skidmore74, T. Skwarnicki67, M.W. Slater52, J.C. Smallwood62,
J.G. Smeaton54, A. Smetkina38, E. Smith13, M. Smith60, A. Snoch31, M.
Soares19, L. Soares Lavra9, M.D. Sokoloff64, F.J.P. Soler58, A. Solovev37, I.
Solovyev37, F.L. Souza De Almeida2, B. Souza De Paula2, B. Spaan14, E. Spadaro
Norella25,n, P. Spradlin58, F. Stagni47, M. Stahl64, S. Stahl47, P. Stefko48,
O. Steinkamp49,80, S. Stemmle16, O. Stenyakin43, H. Stevens14, S. Stone67,
M.E. Stramaglia48, M. Straticiuc36, D. Strekalina80, S. Strokov82, F.
Suljik62, J. Sun26, L. Sun72, Y. Sun65, P. Svihra61, P.N. Swallow52, K.
Swientek34, A. Szabelski35, T. Szumlak34, M. Szymanski47, S. Taneja61, Z.
Tang3, T. Tekampe14, F. Teubert47, E. Thomas47, K.A. Thomson59, M.J. Tilley60,
V. Tisserand9, S. T’Jampens8, M. Tobin6, S. Tolk47, L. Tomassetti20,g, D.
Torres Machado1, D.Y. Tou12, M. Traill58, M.T. Tran48, E. Trifonova80, C.
Trippl48, A. Tsaregorodtsev10, G. Tuci28,m, A. Tully48, N. Tuning31, A.
Ukleja35, D.J. Unverzagt16, A. Usachov31, A. Ustyuzhanin41,81, U. Uwer16, A.
Vagner82, V. Vagnoni19, A. Valassi47, G. Valenti19, N. Valls Canudas44, M. van
Beuzekom31, H. Van Hecke66, E. van Herwijnen80, C.B. Van Hulse17, M. van
Veghel77, R. Vazquez Gomez45, P. Vazquez Regueiro45, C. Vázquez Sierra31, S.
Vecchi20, J.J. Velthuis53, M. Veltri21,o, A. Venkateswaran67, M. Veronesi31,
M. Vesterinen55, D. Vieira64, M. Vieites Diaz48, H. Viemann75, X. Vilasis-
Cardona83, E. Vilella Figueras59, P. Vincent12, G. Vitali28, A. Vollhardt49,
D. Vom Bruch12, A. Vorobyev37, V. Vorobyev42,u, N. Voropaev37, R. Waldi75, J.
Walsh28, C. Wang16, J. Wang3, J. Wang72, J. Wang4, J. Wang6, M. Wang3, R.
Wang53, Y. Wang7, Z. Wang49, D.R. Ward54, H.M. Wark59, N.K. Watson52, S.G.
Weber12, D. Websdale60, C. Weisser63, B.D.C. Westhenry53, D.J. White61, M.
Whitehead53, D. Wiedner14, G. Wilkinson62, M. Wilkinson67, I. Williams54, M.
Williams63,68, M.R.J. Williams57, F.F. Wilson56, W. Wislicki35, M. Witek33, L.
Witola16, G. Wormser11, S.A. Wotton54, H. Wu67, K. Wyllie47, Z. Xiang5, D.
Xiao7, Y. Xie7, H. Xing71, A. Xu4, J. Xu5, L. Xu3, M. Xu7, Q. Xu5, Z. Xu5, Z.
Xu4, D. Yang3, Y. Yang5, Z. Yang3, Z. Yang65, Y. Yao67, L.E. Yeomans59, H.
Yin7, J. Yu70, X. Yuan67, O. Yushchenko43, K.A. Zarebski52, M. Zavertyaev15,c,
M. Zdybal33, O. Zenaiev47, M. Zeng3, D. Zhang7, L. Zhang3, S. Zhang4, Y.
Zhang47, A. Zhelezov16, Y. Zheng5, X. Zhou5, Y. Zhou5, X. Zhu3, V.
Zhukov13,39, J.B. Zonneveld57, S. Zucchelli19,e, D. Zuliani27, G. Zunica61.

1Centro Brasileiro de Pesquisas Físicas (CBPF), Rio de Janeiro, Brazil
2Universidade Federal do Rio de Janeiro (UFRJ), Rio de Janeiro, Brazil
3Center for High Energy Physics, Tsinghua University, Beijing, China
4School of Physics State Key Laboratory of Nuclear Physics and Technology,
Peking University, Beijing, China
5University of Chinese Academy of Sciences, Beijing, China
6Institute Of High Energy Physics (IHEP), Beijing, China
7Institute of Particle Physics, Central China Normal University, Wuhan, Hubei,
China
8Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IN2P3-LAPP, Annecy,
France
9Université Clermont Auvergne, CNRS/IN2P3, LPC, Clermont-Ferrand, France
10Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France
11Université Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France
12LPNHE, Sorbonne Université, Paris Diderot Sorbonne Paris Cité, CNRS/IN2P3,
Paris, France
13I. Physikalisches Institut, RWTH Aachen University, Aachen, Germany
14Fakultät Physik, Technische Universität Dortmund, Dortmund, Germany
15Max-Planck-Institut für Kernphysik (MPIK), Heidelberg, Germany
16Physikalisches Institut, Ruprecht-Karls-Universität Heidelberg, Heidelberg,
Germany
17School of Physics, University College Dublin, Dublin, Ireland
18INFN Sezione di Bari, Bari, Italy
19INFN Sezione di Bologna, Bologna, Italy
20INFN Sezione di Ferrara, Ferrara, Italy
21INFN Sezione di Firenze, Firenze, Italy
22INFN Laboratori Nazionali di Frascati, Frascati, Italy
23INFN Sezione di Genova, Genova, Italy
24INFN Sezione di Milano-Bicocca, Milano, Italy
25INFN Sezione di Milano, Milano, Italy
26INFN Sezione di Cagliari, Monserrato, Italy
27Universita degli Studi di Padova, Universita e INFN, Padova, Padova, Italy
28INFN Sezione di Pisa, Pisa, Italy
29INFN Sezione di Roma Tor Vergata, Roma, Italy
30INFN Sezione di Roma La Sapienza, Roma, Italy
31Nikhef National Institute for Subatomic Physics, Amsterdam, Netherlands
32Nikhef National Institute for Subatomic Physics and VU University Amsterdam,
Amsterdam, Netherlands
33Henryk Niewodniczanski Institute of Nuclear Physics Polish Academy of
Sciences, Kraków, Poland
34AGH - University of Science and Technology, Faculty of Physics and Applied
Computer Science, Kraków, Poland
35National Center for Nuclear Research (NCBJ), Warsaw, Poland
36Horia Hulubei National Institute of Physics and Nuclear Engineering,
Bucharest-Magurele, Romania
37Petersburg Nuclear Physics Institute NRC Kurchatov Institute (PNPI NRC KI),
Gatchina, Russia
38Institute of Theoretical and Experimental Physics NRC Kurchatov Institute
(ITEP NRC KI), Moscow, Russia, Moscow, Russia
39Institute of Nuclear Physics, Moscow State University (SINP MSU), Moscow,
Russia
40Institute for Nuclear Research of the Russian Academy of Sciences (INR RAS),
Moscow, Russia
41Yandex School of Data Analysis, Moscow, Russia
42Budker Institute of Nuclear Physics (SB RAS), Novosibirsk, Russia
43Institute for High Energy Physics NRC Kurchatov Institute (IHEP NRC KI),
Protvino, Russia, Protvino, Russia
44ICCUB, Universitat de Barcelona, Barcelona, Spain
45Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de
Santiago de Compostela, Santiago de Compostela, Spain
46Instituto de Fisica Corpuscular, Centro Mixto Universidad de Valencia -
CSIC, Valencia, Spain
47European Organization for Nuclear Research (CERN), Geneva, Switzerland
48Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL),
Lausanne, Switzerland
49Physik-Institut, Universität Zürich, Zürich, Switzerland
50NSC Kharkiv Institute of Physics and Technology (NSC KIPT), Kharkiv, Ukraine
51Institute for Nuclear Research of the National Academy of Sciences (KINR),
Kyiv, Ukraine
52University of Birmingham, Birmingham, United Kingdom
53H.H. Wills Physics Laboratory, University of Bristol, Bristol, United
Kingdom
54Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom
55Department of Physics, University of Warwick, Coventry, United Kingdom
56STFC Rutherford Appleton Laboratory, Didcot, United Kingdom
57School of Physics and Astronomy, University of Edinburgh, Edinburgh, United
Kingdom
58School of Physics and Astronomy, University of Glasgow, Glasgow, United
Kingdom
59Oliver Lodge Laboratory, University of Liverpool, Liverpool, United Kingdom
60Imperial College London, London, United Kingdom
61Department of Physics and Astronomy, University of Manchester, Manchester,
United Kingdom
62Department of Physics, University of Oxford, Oxford, United Kingdom
63Massachusetts Institute of Technology, Cambridge, MA, United States
64University of Cincinnati, Cincinnati, OH, United States
65University of Maryland, College Park, MD, United States
66Los Alamos National Laboratory (LANL), Los Alamos, United States
67Syracuse University, Syracuse, NY, United States
68School of Physics and Astronomy, Monash University, Melbourne, Australia,
associated to 55
69Pontifícia Universidade Católica do Rio de Janeiro (PUC-Rio), Rio de
Janeiro, Brazil, associated to 2
70Physics and Micro Electronic College, Hunan University, Changsha City,
China, associated to 7
71Guangdong Provencial Key Laboratory of Nuclear Science, Institute of Quantum
Matter, South China Normal University, Guangzhou, China, associated to 3
72School of Physics and Technology, Wuhan University, Wuhan, China, associated
to 3
73Departamento de Fisica , Universidad Nacional de Colombia, Bogota, Colombia,
associated to 12
74Universität Bonn - Helmholtz-Institut für Strahlen und Kernphysik, Bonn,
Germany, associated to 16
75Institut für Physik, Universität Rostock, Rostock, Germany, associated to 16
76INFN Sezione di Perugia, Perugia, Italy, associated to 20
77Van Swinderen Institute, University of Groningen, Groningen, Netherlands,
associated to 31
78Universiteit Maastricht, Maastricht, Netherlands, associated to 31
79National Research Centre Kurchatov Institute, Moscow, Russia, associated to
38
80National University of Science and Technology “MISIS”, Moscow, Russia,
associated to 38
81National Research University Higher School of Economics, Moscow, Russia,
associated to 41
82National Research Tomsk Polytechnic University, Tomsk, Russia, associated to
38
83DS4DS, La Salle, Universitat Ramon Llull, Barcelona, Spain, associated to 44
84University of Michigan, Ann Arbor, United States, associated to 67
85AGH - University of Science and Technology, Faculty of Computer Science,
Electronics and Telecommunications, Kraków, Poland

aUniversidade Federal do Triângulo Mineiro (UFTM), Uberaba-MG, Brazil
bLaboratoire Leprince-Ringuet, Palaiseau, France
cP.N. Lebedev Physical Institute, Russian Academy of Science (LPI RAS),
Moscow, Russia
dUniversità di Bari, Bari, Italy
eUniversità di Bologna, Bologna, Italy
fUniversità di Cagliari, Cagliari, Italy
gUniversità di Ferrara, Ferrara, Italy
hUniversità di Genova, Genova, Italy
iUniversità di Milano Bicocca, Milano, Italy
jUniversità di Roma Tor Vergata, Roma, Italy
kAGH - University of Science and Technology, Faculty of Computer Science,
Electronics and Telecommunications, Kraków, Poland
lUniversità di Padova, Padova, Italy
mUniversità di Pisa, Pisa, Italy
nUniversità degli Studi di Milano, Milano, Italy
oUniversità di Urbino, Urbino, Italy
pUniversità della Basilicata, Potenza, Italy
qScuola Normale Superiore, Pisa, Italy
rUniversità di Modena e Reggio Emilia, Modena, Italy
sUniversità di Siena, Siena, Italy
tMSU - Iligan Institute of Technology (MSU-IIT), Iligan, Philippines
uNovosibirsk State University, Novosibirsk, Russia


# Optimal State Estimation in the Presence of Non-Gaussian Uncertainty via
Wasserstein Distance Minimization

Himanshu Prabhat and Raktim Bhattacharya Himanshu Prabhat and Raktim
Bhattacharya are with the Department of Aerospace Engineering, Texas A&M
University, TX, USA. Their email addresses are<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>respectively.

###### Abstract

This paper presents a novel distribution-agnostic Wasserstein distance-based
estimation framework. The goal is to determine an optimal map combining prior
estimate with measurement likelihood such that posterior estimation error
optimally reaches the Dirac delta distribution with minimal effort. The
Wasserstein metric is used to quantify the effort of transporting from one
distribution to another. We hypothesize that minimizing the Wasserstein
distance between the posterior error and the Dirac delta distribution results
in optimal information fusion and posterior state uncertainty. Framework
validation is demonstrated by the successful recovery of the classical Kalman
filter for linear systems with Gaussian uncertainties. Notably, the proposed
Wasserstein filter does not rely on particle representation of uncertainty.
Furthermore, the classical result for the Gaussian Sum Filter (GSF) is
retrieved from the Wasserstein framework. This approach analytically exhibits
the suboptimality of GSF and enables the use of nonlinear optimization
techniques to enhance the accuracy of the Gaussian sum estimator.

##  I INTRODUCTION

The necessity for nonlinear estimation algorithms arises from various
challenges inherent to complex systems and their associated dynamics.
Nonlinear dynamics often evolve in nonlinear manifolds, necessitating the
development of estimation techniques that can capture these complexities.
Sensor models are frequently nonlinear and exhibit time-varying, state-
dependent noise, complicating the estimation process further. While linearized
dynamics with frequent measurements have proven effective in practice, they
tend to fail when measurements are sparse spatio-temporally. In such cases,
nonlinearities in the dynamics dominate system trajectories, resulting in non-
Gaussian state uncertainties that cannot be resolved using linear
approximations. Particle-based methods have been proposed to address these
issues; however, they are often plagued with concerns related to unbiasedness,
lack of convergence guarantees, and computational inefficiencies.

As a result, nonlinear estimation remains an active area of research, with
numerous novel approaches emerging in the literature. Recently, the use of
optimal transport has been gaining popularity in the field of nonlinear
estimation. Developments in this area show promise, offering potential
solutions to the challenges posed by nonlinear dynamics and complex sensor
models. By incorporating these techniques, researchers are advancing the state
of the art in nonlinear estimation, enhancing its accuracy and reliability in
a wide range of applications.

The first challenge in nonlinear filtering is the propagation of uncertainty
to get the prior state uncertainty, which is governed by advection-diffusion
equations. Wasserstein gradient flows provide a powerful means of
understanding and solving many diffusion equations. Specifically, Fokker-
Planck equations, which model the diffusion of probability measures, can be
understood as gradient descent over entropy functionals in Wasserstein space.
This simplifies the computation of the prior uncertainty in nonlinear
filtering [1, 2, 3].

The second challenge is determining the optimal posterior state uncertainty
with available measurements. This step is difficult as the prior and the
measurement uncertainties can potentially be in different manifolds,
represented with a mix of Cartesian and circular variables. Information fusion
in the optimal transport framework is primarily performed using particle
representation of distributions. Danu et al. [4] treat the problem of
evaluating the cost of particle filter clouds association based on the
Wasserstein distance of different orders, analyzing the implications of clouds
cardinality (for weighted particles), and of various resampling methods (for
unweighted particles). More recently, Reich et al. [5, 6] leveraged the
optimal transport property of Wasserstein distance to resample weighted point
clouds from Bayesian updates and generate unweighted particles. The
formulation works quite well in nonlinear manifolds, as shown by Das et al.
[7, 8]. Taghvaei et al. [9] apply an optimal transport framework to develop
algorithms that feedback measurements in the evolution of the state
probability density function, resulting in an elegant solution to the ensemble
Kalman filtering problem. The resulting system is a controlled interaction
between a finite number of particles that minimizes the Wasserstein distance
between distributions over consecutive time steps. The work of Reich et al.
and Taghvaei et al. overcame the problem ensemble Kalman filtering has with
biases with finite particles and particle degeneracy. These early
breakthroughs have led to considerable research in this field [10, 11, 12, 13,
14].

### Main Contributions

In this paper, we introduce an alternative formulation for data assimilation,
which seeks to determine an optimal map that combines both the prior
distribution and the likelihood. The objective is to minimize the Wasserstein
distance between the posterior error distribution and the Dirac distribution
centered at the origin. Utilizing the Dirac distribution for error enables the
representation of a perfect state estimate characterized by zero error and
dispersion. We hypothesize that minimizing the Wasserstein distance between
the posterior error and the Dirac distribution results in optimal information
fusion and posterior state uncertainty. The Wasserstein distance, which
defines a metric on the manifold of probability density functions (PDFs), is
well-suited for formulating general estimation algorithms applicable to
arbitrary uncertainty representations. It is important to note that the work
presented in this paper does not rely on particle representations of
uncertainty. Instead, it focuses on Gaussian and the mixture of Gaussian
representations, which offers a more tractable and computationally efficient
approach to the problem.

Specifically, the following are the key contributions of the paper:
i) In this paper, we introduce a novel general state-estimation formulation
that aims to minimize the Wasserstein distance between the posterior error
distribution and the Dirac distribution centered at the origin. The framework
is distribution agnostic and can be applied to information fusion over general
manifolds with arbitrary state uncertainty.
ii) To validate the formulation, we demonstrate that the new framework
successfully recovers the classical Kalman filter for linear systems
characterized by Gaussian uncertainty.
iii) The proposed Wasserstein filter recovers the classical result for the
Gaussian Sum Filter. Notably, an analytical proof for the suboptimality
condition of Gaussian Sum Estimators in minimum mean square error (MMSE) sense
is presented using the new formulation.
iv) Lastly, we showcase that by employing a nonlinear optimization technique,
which is warm-started with the suboptimal solution, the accuracy of Gaussian
Sum Estimators can be enhanced.

Key results are presented in Theorems 1 and 2, along with supporting
corollaries.

##  II MATHEMATICAL PRELIMINARIES

######  Definition 1.

(Wasserstein distance) Given $z_{1},z_{2}\in\mathbb{R}^{n}$, with
$\mathcal{P}_{2}(p_{1},p_{2})$ containing the collection of all probability
measures $p$ supported on the product space $\mathbb{R}^{2n}$, having a finite
second moment, with first marginal $p_{1}$ and second marginal $p_{2}$.
2-Wasserstein distance, denoted as $W_{2}$, between two probability measures
$p_{1},p_{2}$, is defined as

$\displaystyle W_{2}(p_{1},p_{2})\triangleq$
$\displaystyle\left(\displaystyle\inf_{p\in\mathcal{P}_{2}(p_{1},p_{2})}\displaystyle\int_{\mathbb{R}^{2n}}\parallel
z_{1}-z_{2}\parallel_{\ell_{2}\left(\mathbb{R}^{n}\right)}^{2}\>dp(z_{1},z_{2})\right)^{\frac{{1}}{2}}.$
(1)

######  Remark 1.

The Wasserstein distance can be interpreted as the cost for the Monge-
Kantorovich optimal transportation plan[15]. This adds an intuitive meaning to
the Wasserstein distance in terms of the least effort needed to morph from one
distribution to another. The particular choice of $\ell_{2}$ norm with order 2
is motivated by [16]. Further, one can prove (p. 208, [15]) that $W_{2}$
defines a metric on the manifold of PDFs, which makes it suitable for
formulating distribution agnostic estimation algorithms.

######  Proposition 1.

The Wasserstein distance $W_{2}$ between two multivariate Gaussians
$\mathcal{N}(\mu_{1},\Sigma_{1})$ and $\mathcal{N}(\mu_{2},\Sigma_{2})$ in
$\mathcal{R}^{n}$ is given by

$\displaystyle
W_{2}^{2}\left(\mathcal{N}(\mu_{1},\Sigma_{1}),\mathcal{N}(\mu_{2},\Sigma_{2})\right)=$
$\displaystyle\|\mu_{1}-\mu_{2}\|^{2}+\textbf{tr}\left(\Sigma_{1}+\Sigma_{2}-2\left(\sqrt{\Sigma_{1}}\Sigma_{2}\sqrt{\Sigma_{1}}\right)^{\frac{1}{2}}\right).$
(2)

######  Corollary 1.

Square of 2-Wasserstein distance between Gaussian distribution
$\mathcal{N}(\mu,\Sigma)$ and the Dirac delta function $\delta(x-\mu_{c})$ is
given by,

$\displaystyle
W_{2}^{2}(\mathcal{N}(\mu,\Sigma),\delta(x-\mu_{c}))=\|\mu-\mu_{c}\|^{2}+\textbf{tr}\left(\Sigma\right).$
(3)

###### Proof.

Defining the Dirac delta function as (see e.g., p. 160-161, [17])

$\displaystyle\delta(x-\mu_{c})=\lim_{\mu\to\mu_{c},\Sigma\to
0}\mathcal{N}(\mu,\Sigma),$

and substituting in (2), we get the result. ∎

The Wasserstein distance in (3) can also be written as,

$\displaystyle
W_{2}^{2}(\mathcal{N}(\mu,\Sigma),\delta(x-\mu_{c}))=\|\mu-\mu_{c}\|^{2}+\textbf{tr}\left(\Sigma\right),$
$\displaystyle={\textbf{tr}\left((\mu-\mu_{c})(\mu-\mu_{c})^{T}+\Sigma\right)}.$
(4)

##  III WASSERSTEIN FILTER

Here we consider the problem of updating the prior state estimate with
available measurement data to arrive at the posterior state estimate. We
assume $x\in\mathcal{R}^{n}$ is the true state. The prior estimate of $x$ is
denoted by $x^{-}$, which is a random variable with associate probability
density function $p_{x^{-}}(x^{-})$. Similarly, the posterior estimate is
denoted by $x^{+}$, also a random variable with associated probability density
function $p_{x^{+}}(x^{+})$. We next assume that measurement
$y\in\mathcal{R}^{m}$ is a function of the true state $x$, but corrupted by an
additive noise $n$. It is modeled as

$\displaystyle y:=g(x)+n.$ (5)

The noise is assumed to be a random variable with associated probability
density function $p_{n}(n)$.

The errors associated with the prior and posterior estimates are defined as

$\displaystyle e^{-}$ $\displaystyle:=x^{-}-x,$ (6) $\displaystyle e^{+}$
$\displaystyle:=x^{+}-x,$ (7)

which are also random variables with probability density functions
$p_{e^{-}}(e^{-})$ and $p_{e^{+}}(e^{+})$ respectively.

The objective of optimal estimation is to determine $x^{+}$ from $x^{-}$ and
$y$ by minimizing some statistical quantity associated with $e^{+}$. More
generally, we can define the optimal estimation problem as

$\displaystyle\min_{T}d(e^{+})$ (8)

where $T:(x^{-},y)\mapsto x^{+}$, $d(e^{+})$ is some cost function to be
minimized.

In optimal Wasserstein filtering, we define
$d(e^{+}):=W^{2}_{2}(p_{e^{+}}(e^{+}),\delta(e^{+}))$, i.e. minimize the
Wasserstein distance of the posterior-error’s probability density function
from the Dirac at the origin $\delta(e^{+})$. The distribution $\delta(e^{+})$
is the degenerate probability density function that represents the posterior
error with zero mean and variance. The optimization therefore determines the
map $T(\cdot,\cdot)$ that minimizes the Wasserstein distance between
$p_{e^{+}}(e^{+})$ and $\delta(e^{+})$, resulting in the best posterior
estimate of $x$ in this sense.

In the next subsections, we present a few specific cases of (8), which are of
engineering significance.

###  III-A Linear Measurement with Gaussian Uncertainty

Given a linear sensor model $y:=Cx+n,n\sim\mathcal{N}(0,R)$ with Gaussian
prior estimate $x^{-}$, we define a posterior linear map as
$x^{+}=T(x^{-},y):=Gx^{-}+Hy$. It can be shown that $e^{+}$ is Gaussian
distribution $\mathcal{N}(\mu_{e}^{+},\Sigma_{e}^{+})$, with

$\displaystyle e^{+}$ $\displaystyle=x^{+}-x,$
$\displaystyle=Gx^{-}+(HC-I)(x^{-}-e^{-})+Hn,$
$\displaystyle=\begin{bmatrix}(G+HC-I)&-(HC-I)&H\end{bmatrix}\begin{pmatrix}x^{-}\\\
e^{-}\\\ n\end{pmatrix}.$ (9)

If $\mathbb{E}\left[n\right]=0$ and $\mathbb{E}\left[e^{-}\right]=0$, then

$\displaystyle\mu_{e}^{+}$
$\displaystyle:=\mathbb{E}\left[e^{+}\right]=(G+HC-I)\mu^{-},$ (10a)
$\displaystyle\Sigma_{e}^{+}$
$\displaystyle:=\mathbb{E}\left[e^{+}{e^{+}}^{T}\right]-\mu_{e}^{+}{\mu_{e}^{+}}^{T}$
$\displaystyle=\begin{bmatrix}(G+HC-I)&-(HC-I)&H\end{bmatrix}.$
$\displaystyle\mathbb{E}\left[\begin{pmatrix}x^{-}\\\ e^{-}\\\
n\end{pmatrix}\begin{pmatrix}x^{-}\\\ e^{-}\\\
n\end{pmatrix}^{T}\right]\begin{bmatrix}(G+HC-I)^{T}\\\ -(HC-I)^{T}\\\
H^{T}\end{bmatrix}$ $\displaystyle-\mu_{e}^{+}{\mu_{e}^{+}}^{T}.$ (10b)

In the case when, the joint probability density function associated with
$e^{+}$ is Gaussian, denoted by $\mathcal{N}(\mu_{e}^{+},\Sigma_{e}^{+})$, the
Wasserstein distance between $\mathcal{N}(\mu_{e}^{+},\Sigma_{e}^{+})$ and the
Dirac delta function $\delta(e^{+})$ is given by,

$\displaystyle
W_{2}^{2}(\mathcal{N}(\mu_{e}^{+},\Sigma_{e}^{+}),\delta(e^{+}))$
$\displaystyle=\textbf{tr}\left(\mu_{e}^{+}{\mu_{e}^{+}}^{T}+\Sigma_{e}^{+}\right).$
(11)

Therefore, minimizing
$W_{2}^{2}(\mathcal{N}(\mu_{e}^{+},\Sigma_{e}^{+}),\delta(e^{+}))$ is the same
as minimizing the posterior error variance, as it is done in Kalman filtering.
The following theorem validates our argument that the proposed Wasserstein
filter for a linear map between Gaussian distributed posterior and prior is
the same as the Kalman filter.

######  Proposition 2.

(Orthogonality Condition) Consider a linear estimator with generic linear
mapping between posterior and prior at the $k^{\text{th}}$ time step as

$\displaystyle{x_{k}}^{+}:=G{x_{k}}^{-}+Hy_{k}.$ (12)

Let the posterior error be defined as ${e_{k}}^{+}:={x_{k}}^{+}-x_{k}$. The
necessary conditions for a linear estimator to be optimal in the MMSE sense is
given by

$\displaystyle\mathbb{E}\left[{e_{k}^{+}}{x_{k}^{-}}^{T}\right]$
$\displaystyle=0$ (13)
$\displaystyle\mathbb{E}\left[{e_{k}^{+}}{y_{k}}^{T}\right]$
$\displaystyle=0.$ (14)

This implies that at any time step $k$, the posterior error must be orthogonal
to the prior state and current measurement.

###### Proof.

Given the linear estimator (12), the cost function for MMSE optimization can
be written as

$\displaystyle J$
$\displaystyle=\textbf{tr}\left\\{\mathbb{E}\left[{e_{k}^{+}}{e_{k}^{+}}^{T}\right]\right\\}$
$\displaystyle=\textbf{tr}\left\\{\mathbb{E}\left[{e_{k}^{+}}{x_{k}^{-}}^{T}\right]G^{T}+\right.$
$\displaystyle\left.{}\mathbb{E}\left[{e_{k}^{+}}{y_{k}}^{T}\right]H^{T}-\mathbb{E}\left[{e_{k}^{+}}{x_{k}}^{T}\right]\right\\}.$
(15)

Applying necessary first-order optimal conditions,$\frac{\partial J}{\partial
G}=0;\frac{\partial J}{\partial H}=0$, (13) and (14) are obtained. ∎

######  Corollary 2.

Given an optimal linear estimator in the MMSE sense, consider a linear time-
invariant state propagation model $x_{k+1}=A{x}_{k}+w_{k}$, with propagation
noise $w_{k}\sim\mathcal{N}(0,Q)$. It can be shown that prior state and error
estimate at any time step $k$, $x_{k}^{-}$ and $e_{k}^{-}:=x_{k}^{-}-x$
respectively, are orthogonal.

###### Proof.

Following the given definitions, we simplify
$\mathbb{E}\left[{e_{k}^{-}}{x_{k}^{-}}^{T}\right]$ as:

$\displaystyle\mathbb{E}\left[{e_{k}^{-}}{x_{k}^{-}}^{T}\right]$
$\displaystyle\quad=\mathbb{E}\left[{({x_{k}^{-}}-x_{k})}{x_{k}^{-}}^{T}\right]$
$\displaystyle\quad=\mathbb{E}\left[{({A{x}_{k-1}^{+}}-{A{x_{k-1}}}-w_{k-1})}{A{x}_{k-1}^{+}}^{T}\right]$
$\displaystyle\quad=\mathbb{E}\left[{({A{e}_{k-1}^{+}}-w_{k-1})}{(A{x}_{k-1}^{+})}^{T}\right]$
$\displaystyle\quad=\mathbb{E}\left[{({A{e}_{k-1}^{+}})}{(A{x}_{k-1}^{+})}^{T}\right]-\mathbb{E}\left[({w_{k-1})}{(A{x}_{k-1}^{+})}^{T}\right]$
$\displaystyle\quad=A{\mathbb{E}\left[{{e}_{k-1}^{+}}{x_{k-1}^{+}}^{T}\right]}A^{T}-{\mathbb{E}\left[{w_{k-1}}{{x}_{k-1}^{+}}^{T}\right]}A^{T}$
(16)

Using the uncorrelated property of white noise, (16) reduced to

$\displaystyle\mathbb{E}\left[{{e}_{k}^{-}}{{x}_{k}^{-}}^{T}\right]$
$\displaystyle\quad=A{\mathbb{E}\left[{{e}_{k-1}^{+}}{x_{k-1}^{+}}^{T}\right]}A^{T}-{\mathbb{E}\left[{w_{k-1}}\right]\mathbb{E}\left[{{x}_{k-1}^{+}}^{T}\right]}A^{T}$
$\displaystyle\quad=A{\mathbb{E}\left[{{e}_{k-1}^{+}}{x_{k-1}^{+}}^{T}\right]}A^{T}.$
(17)

Following the results of proposition 2, the remaining term in (17) drops off.
Thus,

$\displaystyle\mathbb{E}\left[{{e}_{k}^{-}}{{x}_{k}^{-}}^{T}\right]=0.$ (18)

Eq.(18) establishes the orthogonal relation between prior state estimate and
prior error at any time step $k$. ∎

######  Theorem 1.

Given linear sensor model $y:=Cx+n$, with measurement noise
$n\sim\mathcal{N}(0,R)$, and defining a linear map between posterior and prior
state estimates as $x^{+}=T(x^{-},y):=Gx^{-}+Hy$, minimization of Wasserstein
distance $d(e^{+}):=W^{2}_{2}(p_{e^{+}}(e^{+}),\delta(e^{+}))$ results in the
Kalman filter with $G^{\ast}:=I-H^{\ast}C$ and
$H^{\ast}:=\Sigma_{e}^{-}C^{T}(C\Sigma_{e}^{-}C^{T}+R)^{-1}$.

###### Proof.

Using (4), the optimization problem can be formally stated as:

$\displaystyle\min_{T}\textbf{tr}\left(\mu_{e}^{+}{\mu_{e}^{+}}^{T}+\Sigma_{e}^{+}\right).$
(19)

Substituting (10a) and (10b) in (19), we get the cost function $J$ as:

$\displaystyle\footnotesize J$
$\displaystyle=\textbf{tr}\Bigg{\\{}\begin{bmatrix}(G+HC-I)&-(HC-I)&H\end{bmatrix}\times$
$\displaystyle\mathbb{E}\left[\begin{pmatrix}x^{-}\\\ e^{-}\\\
n\end{pmatrix}\begin{pmatrix}x^{-}\\\ e^{-}\\\
n\end{pmatrix}^{T}\right]\begin{bmatrix}(G+HC-I)^{T}\\\ -(HC-I)^{T}\\\
H^{T}\end{bmatrix}\Bigg{\\}}$ (20)

With $\Sigma_{e}^{-}=\mathbb{E}\left[e^{-}{e^{-}}^{T}\right]$, we expand the
cost function and further simplify the expression using properties of
measurement white noise as:

$\displaystyle J$
$\displaystyle=\textbf{tr}\\{(G+HC-I){\mathbb{E}\left[x^{-}{x^{-}}^{T}\right]}{(G+HC-I)^{T}}-$
$\displaystyle(G+HC-I){\mathbb{E}\left[x^{-}{e^{-}}^{T}\right]}{(HC-I)^{T}}-$
$\displaystyle(HC-I){\mathbb{E}\left[x^{-}{e^{-}}^{T}\right]}{(G+HC-I)^{T}}+$
$\displaystyle(HC-I){\Sigma_{e}^{-}}{(HC-I)^{T}}+HR{H}^{T}\\}$ (21)

Enforcing the first order optimal conditions $\frac{\partial J}{\partial
G}=0;\frac{\partial J}{\partial H}=0$, we get:

$\displaystyle(G+HC-I){\mathbb{E}\left[x^{-}{x^{-}}^{T}\right]}$
$\displaystyle-(HC-I){\mathbb{E}\left[x^{-}{e^{-}}^{T}\right]}=0$ (22a)
$\displaystyle(HC-I){\Sigma_{e}^{-}}{C^{T}}-(G+HC-I){\mathbb{E}\left[x^{-}{e^{-}}^{T}\right]}{C^{T}}$
$\displaystyle+HR=0$ (22b)

∎

An analytical solution for (22a), (22b) can be obtained by setting $G+HC-I=0$.
Following the assumption, (22a) reduced to
${\mathbb{E}\left[x^{-}{e^{-}}^{T}\right]}=0$. Indeed, this result holds for
linear estimator optimal in MMSE sense as shown in corollary 2. Solving (22b),
we obtain corresponding stationary point,

$\displaystyle H^{\ast}$
$\displaystyle={\Sigma_{e}^{-}}{C^{T}}(C{\Sigma_{e}^{-}}C^{T}+R)^{-1}$ (23a)
$\displaystyle G^{\ast}$
$\displaystyle=I-{\Sigma_{e}^{-}}{C^{T}}(C{\Sigma_{e}^{-}}C^{T}+R)^{-1}C$
(23b)

The positive definiteness of the Hessian matrix
$\textbf{H}=(C{\Sigma_{e}^{-}}C^{T}+R)$ affirms the optimality of the obtained
stationary point. Further, posterior error mean $(\mu_{e}^{+})$, covariance
$(\Sigma_{e}^{+})$, and measurement update are obtained as:

$\displaystyle\mu_{e}^{+}$ $\displaystyle=0,$ (24a)
$\displaystyle\Sigma_{e}^{+}$
$\displaystyle=(HC-I){\Sigma_{e}^{-}}(HC-I)^{T}+HRH^{T},$
$\displaystyle=H(C{\Sigma_{e}^{-}}C^{T}+R)H^{T}+{\Sigma_{e}^{-}}-$
$\displaystyle\qquad HC{\Sigma_{e}^{-}}-{\Sigma_{e}^{-}}{C^{T}}{H^{T}},$
$\displaystyle={\Sigma_{e}^{-}}-{\Sigma_{e}^{-}}{C^{T}}(C{\Sigma_{e}^{-}}C^{T}+R)^{-1}C{\Sigma_{e}^{-}}$
(24b) $\displaystyle x^{+}$
$\displaystyle=x^{-}+{\Sigma_{e}^{-}}{C^{T}}(C{\Sigma_{e}^{-}}C^{T}+R)^{-1}(y-Cx^{-}).$
(24c)

Therefore, for linear measurements and Gaussian uncertainty models, the
Wasserstein filter is equivalent to the Kalman filter.

###  III-B Mixture of Gaussian (MoG) Prior

In the previous section, we considered Gaussian prior estimate with a linear
measurement model with Gaussian noises(5). In this section, we consider the
case of fusing multimodal Gaussian mixture prior i.e.,

$\displaystyle
x^{-}\sim\sum_{i}^{M}\lambda_{i}^{-}\mathcal{N}(\mu_{i}^{-},\Sigma_{i}^{-}),$
(25)

where $0\leq\lambda_{i}\leq 1$ and

$\sum_{i}^{M}\lambda_{i}=1,$

with a linear measurement model containing Gaussian measurement noise.
Alternatively, the $M$-order Gaussian mixture prior can be expressed using the
generative model approach with a random vector
$\boldsymbol{\beta}^{-}:=[\beta_{1}^{-},..,\beta_{M}^{-}],$ indicating on each
generating Gaussian, with pdf given by [18],

$\displaystyle\rho(\boldsymbol{\beta}^{-})$
$\displaystyle=\sum_{j=1}^{M}{\lambda_{j}^{-}}\delta({\beta_{j}^{-}}-1)$ (26)
where, $\displaystyle\beta_{j}^{-}$ $\displaystyle=\begin{cases}1,&\text{if
}x^{-}\text{generated by }j^{th}\text{ Gaussian}\\\
0,&\text{otherwise}\end{cases}$ (27) and, $\displaystyle P$
$\displaystyle(\beta_{j}^{-}=1):=P(\beta_{j}^{-})=\lambda_{j}^{-}.$ (28)

We introduce a latent random variable
${z_{j}^{-}}:=(x^{-}|\beta_{j}^{-})\sim\mathcal{N}(\mu_{j}^{-},\Sigma_{j}^{-})$
characterizing the $j^{th}$ Gaussian node. Using the generative model
approach, the multimodal Gaussian prior can thus be written as

$\displaystyle x^{-}\sim\sum_{i}^{M}P(\beta_{i}^{-})\rho(z_{i}^{-}).$ (29)

With linear Gaussian measurement, we can assume the posterior distribution to
be of the form

$\displaystyle x^{+}$
$\displaystyle\sim\sum_{i}^{M}\lambda_{i}^{+}N(\mu_{i}^{+},\Sigma_{i}^{+})$
or, $\displaystyle x^{+}$
$\displaystyle\sim\sum_{i}^{M}P(\beta_{i}^{+})\rho(z_{i}^{+})$ (30) where,
$\displaystyle z_{i}^{+}$
$\displaystyle:=(x^{+}|\beta_{i}^{+})\sim\mathcal{N}(\mu_{i}^{+},\Sigma_{i}^{+}).$

This interpretation of a mixture of Gaussian(MoG) suggests that there are $M$
latent variables representing each Gaussian. Fusion of a MoG prior with
Gaussian likelihood results in posterior MoG with the same number of nodes.
Each prior Gaussian node is one-to-one mapped with another posterior node and
we assume that this mapping is linear, but the map’s parameters are dependent
upon the mode and the measurement values. The linear mapping between each
corresponding $j^{th}$ Gaussian node of posterior and prior is given by

$\displaystyle z_{i}^{+}=G_{i}z_{i}^{-}+H_{i}y$ (31)

As in (7), we define the posterior error variable as,

$\displaystyle e^{+}=x^{+}-x,$ (32)

where $x$ is the deterministic true state variable. $e^{+}$ is a shifted
version of the distribution of $x^{+}$ with pdf given by

$\displaystyle e^{+}$
$\displaystyle\sim\sum_{i}^{M}\lambda_{i}^{+}\mathcal{N}(\mu_{i}^{+}-x,\Sigma_{i}^{+})$
(33) $\displaystyle\text{ or},$ $\displaystyle e^{+}$
$\displaystyle\sim\sum_{i}^{M}P(\beta_{i}^{+})\rho(e_{i}^{+})$ (34)
$\displaystyle\text{where},$ $\displaystyle e_{i}^{+}$
$\displaystyle:=z_{i}^{+}-x.$ (35)

######  Proposition 3.

($W_{2}$ between MoG and Dirac reference PDF) The Wasserstein distance between
a multimodal Gaussian Mixture Distribution and the Dirac delta at origin is
given by,

$\displaystyle W_{2}^{2}=\sum_{i}^{M}\lambda_{i}W_{i}^{2}.$ (36)

###### Proof.

refer Thm. 2 [19] ∎

######  Remark 2.

Given the Gaussian mixture prior alongside linear Gaussian measurement, the
Wasserstein distance between posterior error $e^{+}$ and Dirac delta
distribution centered at the origin is given by,

$\displaystyle W_{2}^{2}$
$\displaystyle=\sum_{i}^{M}\lambda_{i}^{+}W_{i}^{2}(N(\mu_{i}^{+}-x,\Sigma_{i}^{+}),\delta(y))$
$\displaystyle=\sum_{i}^{M}tr(\lambda_{i}^{+}\\{(\mu_{i}^{+}-x)({\mu_{i}^{+}-x})^{T}+\Sigma_{i}^{+}\\})$
$\displaystyle=\sum_{i}^{M}tr(\lambda_{i}^{+}\mathbb{E}\left[{e_{i}^{+}}{e_{i}^{+}}^{T}\right])$
(37)

Using (31), (35) in (37) we get,

$\displaystyle W_{i}^{2}$
$\displaystyle=\textbf{tr}\\{\begin{bmatrix}(G_{i}+H_{i}C-I)&-(H_{i}C-I)&H_{i}\end{bmatrix}.$
$\displaystyle\mathbb{E}\left[\begin{pmatrix}z_{i}^{-}\\\ e_{i}^{-}\\\
n\end{pmatrix}\begin{pmatrix}z_{i}^{-}\\\ e_{i}^{-}\\\
n\end{pmatrix}^{T}\right]\begin{bmatrix}(G_{i}+H_{i}C-I)^{T}\\\
-(H_{i}C-I)^{T}\\\ H_{i}^{T}\end{bmatrix}\\}$ (38)

Wasserstein distance minimization for MoG prior results in the following
optimization problem:

$\displaystyle\min_{\\{\lambda_{i},G_{i},H_{i}\\}_{\forall
i\in[1,M]}}\sum_{i}^{M}\textbf{tr}\\{\lambda_{i}^{+}\mathbb{E}\left[{e_{i}^{+}}{e_{i}^{+}}^{T}\right]\\}$
(39)

subject to

$\sum_{i}^{M}\lambda_{i}^{+}=1,\lambda_{i}\geq 0.$

Although the objective function is non-convex, there exists a convex upper
bound

$\displaystyle\sum_{i}^{M}\textbf{tr}\\{\lambda_{i}^{+}\mathbb{E}\left[{e_{i}^{+}}{e_{i}^{+}}^{T}\right]\\}\leq\sum_{i}^{M}\textbf{tr}\\{\mathbb{E}\left[{e_{i}^{+}}{e_{i}^{+}}^{T}\right]\\},$
(40)

which can be conveniently optimized. It can be shown that bound minimization
leads to Gaussian sum filter (GSF)[18].

######  Theorem 2.

Given the prior distribution represented by MoG alongside linear sensor model
$y:=Cx+n;n\sim\mathcal{N}(0,R)$, we assume linear maps between corresponding
posterior and prior latent variables $z_{i}s$ as
$z_{i}^{+}=T_{i}(z_{i}^{-},y):=G_{i}z_{i}^{-}+H_{i}y$. The minimization of the
upper bound of Wasserstein distance
$d(e^{+})\leq\sum_{i}^{M}W^{2}_{2}(p_{e_{i}^{+}}(e_{i}^{+}),\delta(e_{i}^{+}))$
given in (40) results in the Gaussian sum filter (GSF) with
$G_{i}^{\ast}:=I-H_{i}^{\ast}C$ and
$H_{i}^{\ast}:=\Sigma_{i}^{-}C^{T}(C\Sigma_{i}^{-}C^{T}+R)^{-1}$.

###### Proof.

Relaxed minimization problem can be derived from (39) as,

$\displaystyle\sum_{i}^{M}\textbf{tr}\\{\lambda_{i}^{+}\mathbb{E}\left[{e_{i}^{+}}{e_{i}^{+}}^{T}\right]\\}\leq\sum_{i}^{M}\textbf{tr}\\{\mathbb{E}\left[{e_{i}^{+}}{e_{i}^{+}}^{T}\right]\\}.$
(41)

Since the cost function is a sum of squares with decoupled associated
parameters, each term can be independently optimized. Following the expansion
of the $i^{th}$ term using (38), we have

$\displaystyle J_{i}$
$\displaystyle=\textbf{tr}\\{(G_{i}+H_{i}C-I){\mathbb{E}\left[x_{i}^{-}{x_{i}^{-}}^{T}\right]}{(G_{i}+H_{i}C-I)^{T}}$
$\displaystyle-(G_{i}+H_{i}C-I){\mathbb{E}\left[x_{i}^{-}{e_{i}^{-}}^{T}\right]}{(H_{i}C-I)^{T}}-$
$\displaystyle(H_{i}C-I){\mathbb{E}\left[x_{i}^{-}{e_{i}^{-}}^{T}\right]}{(G_{i}+H_{i}C-I)^{T}}+$
$\displaystyle(H_{i}C-I){\Sigma_{i}^{-}}{(H_{i}C-I)^{T}}+H_{i}R{H_{i}}^{T}\\}.$
(42)

Considering the linear map between Gaussian marginalized prior and posterior
with minimum mean squared marginalized error as the objective function, the
orthogonality conditions mentioned in proposition 2 holds well, i.e.

$\displaystyle\mathbb{E}\left[{{e_{i}}_{k}^{-}}{{x_{i}}_{k}^{-}}^{T}\right]=0.$
(43)

The objective function minimization follows the same approach detailed in
theorem 1, thus leading to the following stationary point:

$\displaystyle H_{i}^{\ast}$
$\displaystyle={\Sigma_{i}^{-}}{C^{T}}(C{\Sigma_{i}^{-}}C^{T}+R)^{-1}$ (44a)
$\displaystyle G_{i}^{\ast}$
$\displaystyle=I-{\Sigma_{i}^{-}}{C^{T}}(C{\Sigma_{i}^{-}}C^{T}+R)^{-1}C$
(44b)

A reasonable estimate of mixing coefficients $\lambda_{i}s$ can be obtained
using the Bayesian approach, [20]

$\displaystyle\lambda_{i}^{+}=\frac{{\lambda_{i}^{-}}\mathcal{N}(y;C\mu_{i}^{-},C^{T}\Sigma_{i}^{-}C+R)}{\sum_{j}^{M}\lambda_{j}^{-}\mathcal{N}(y;C\mu_{j}^{-},C^{T}\Sigma_{j}^{-}C+R)}$
(45)

∎

Algorithm 1 GSF Measurement Update

$\textbf{Given:}\ \mu_{i}^{-},\Sigma_{i}^{-}\ \forall i=1,..M$

$\textbf{Posterior Parameters:}\ \\{H_{i},\lambda_{i}^{+}\\}$

$\ H_{i}={\Sigma_{i}^{-}}{C^{T}}(C{\Sigma_{i}^{-}}C^{T}+R)^{-1}$

$\
\lambda_{i}^{+}=\frac{{\lambda_{i}^{-}}\mathcal{N}(y;C\mu_{i}^{-},C^{T}\Sigma_{i}^{-}C+R)}{\sum_{j}^{M}\lambda_{j}^{-}\mathcal{N}(y;C\mu_{j}^{-},C^{T}\Sigma_{j}^{-}C+R)}$

Posterior Update:

$\ \mu_{i}^{+}=\mu_{i}^{-}+H_{i}(y-C\mu_{i}^{-})$

$\ {\Sigma_{i}^{+}}={\Sigma_{i}^{-}}-H_{i}C{\Sigma_{i}^{-}}$

$\textbf{Posterior GMM:}\
x^{+}\sim\sum_{i}^{M}\lambda_{i}^{+}N(\mu_{i}^{+},\Sigma_{i}^{+})$

It should be noted that obtained $H_{i}^{\ast},\ \lambda_{i}^{+}$ for GSF are
optimal only for the relaxed minimization problem in (41) and sub-optimal for
exact optimization problem in (39). Considering the sub-optimal nature of the
solution, an intuitive idea is to perform a local search in the vicinity of
$(H_{i}^{\ast},\ \lambda_{i}^{+}),$ with the goal of improvement over the GSF.
The cost function in (42) can be simplified using (43) and trivial first order
optimal condition $G_{i}+H_{i}C-I=0$,

$\displaystyle J_{i}$ $\displaystyle=$
$\displaystyle\qquad\textbf{tr}\\{(H_{i}C-I){\Sigma_{i}^{-}}{(H_{i}C-I)^{T}}+H_{i}R{H_{i}}^{T}\\}.$
(46)

Using (46), we formulate the nonlinear Gaussian sum filter (nGSF) which solves
the following nonlinear constrained optimization with sub-optimal solutions in
(44a), (44b), and (45) as initial guess,

$\displaystyle\hat{J}$ $\displaystyle=min_{\\{\lambda_{i}^{+},H_{i}\\}}$
$\displaystyle\sum_{i}^{M}\lambda_{i}^{+}\textbf{tr}\\{(H_{i}C-I){\Sigma_{i}^{-}}{(H_{i}C-I)^{T}}+H_{i}R{H_{i}}^{T}\\},$
(47)

subject to

$\sum_{i}^{M}\lambda_{i}^{+}=1;\lambda_{i}^{+}\geq 0.$

Algorithm 2 nGSF Measurement Update

$\textbf{Given:}\ \mu_{i}^{-},\Sigma_{i}^{-}\ \forall i=1,..M$

$\textbf{Initial Guess:}\ \\{H_{i0},\lambda_{i0}^{+}\\}$

$\ H_{i0}={\Sigma_{i}^{-}}{C^{T}}(C{\Sigma_{i}^{-}}C^{T}+R)^{-1}$

$\
\lambda_{i0}^{+}=\frac{{\lambda_{i}^{-}}\mathcal{N}(y;C\mu_{i}^{-},C^{T}\Sigma_{i}^{-}C+R)}{\sum_{j}^{M}\lambda_{j}^{-}\mathcal{N}(y;C\mu_{j}^{-},C^{T}\Sigma_{j}^{-}C+R)}$

Nonlinear Optimization:

$\ \textit{Cost:}\ \min_{\\{\lambda_{i}^{+},H_{i}\\}_{\forall
i\in[1,M]}}\sum_{i}^{M}\textbf{tr}\\{\lambda_{i}^{+}\mathbb{E}\left[{e_{i}^{+}}{e_{i}^{+}}^{T}\right]\\}$

$\ \textit{Constraints:}\ \sum_{i}^{M}\lambda_{i}^{+}=1;\lambda_{i}^{+}\geq 0$

Posterior Update:

$\ \mu_{i}^{+}=\mu_{i}^{-}+H_{i}^{\ast}(y-C\mu_{i}^{-})$

$\ {\Sigma_{i}^{+}}={\Sigma_{i}^{-}}-H_{i}^{\ast}C{\Sigma_{i}^{-}}$

$\textbf{Posterior GMM:}\
x^{+}\sim\sum_{i}^{M}\lambda_{i}^{+}N(\mu_{i}^{+},\Sigma_{i}^{+})$

##  IV EXAMPLE

We consider the classical duffing oscillator as an example to compare the
performance of the nonlinear Gaussian sum filter (nGSF) with baseline GSF. The
nonlinear dynamics in terms of two states $x=[x_{1},x_{2}]^{\prime}$ is given
by:

$\displaystyle\dot{x}_{1}=x_{2},\ \dot{x}_{2}=-x_{1}-0.25x_{2}-x_{1}^{3}.$
(48)

We propagate Gaussian distributed uncertainty in the initial condition, i.e.
$\Delta_{x_{0}}\sim\mathcal{N}(0,I)$, using particle ensemble. Next, we obtain
the Gaussian mixture prior by fitting the MoG over the propagated point cloud.
It is evident from Fig.1 that despite the Gaussian uncertainty in the initial
condition, the propagated distribution is non-Gaussian. In the simulated
cases, the prior state distribution at each time step is approximated by a
Gaussian mixture with 10 components.

Considering a linear scalar measurement corrupted by Gaussian noise, i.e.

$\displaystyle y=x+v_{k},\ v_{k}\sim\mathcal{N}(0,0.1),$ (49)

we obtain the posterior MoG using the GSF and nGSF algorithms. Random sampling
from the obtained posterior gives us a prior point cloud for the next time
step. The estimation performance for baseline GSF and nGSF is presented in
Fig.2 and Fig. 3.

Figure 1: Initial and propagated state uncertainty after 1 time step (0.5
seconds). Figure 2: Position ($x_{1}$) estimates with filter sample time:
0.5s. Figure 3: Velocity ($x_{2}$) estimates with filter sample time: 0.5s.

We observe the performance enhancement in terms of reduced error variance in
the case of nGSF compared to GSF.

##  V CONCLUSIONS

A novel distribution agnostic estimation framework based on Wasserstein
distance is formulated in this work. It is demonstrated that for linear
Gaussian uncertainty, the optimal filter in the Wasserstein metric is
essentially the Kalman filter. Moreover, the Gaussian Sum Filter is derived in
the Wasserstein framework with explicit proof of sub-optimality in the MMSE
sense. Notably, other approaches do not provide such explicit guarantees.
Furthermore, the proposed nGSF is shown to enhance the accuracy of the
Gaussian sum estimator using a nonlinear optimization technique, at the cost
of additional computation. In our future work, we will relax the linear
Gaussian measurement restriction and incorporate nonlinear measurement models
with non Gaussian noise.

##  VI ACKNOWLEDGMENTS

This work is supported by AFOSR grant FA9550-22-1-0539 with Dr. Erik Blasch as
the program director.

## References

  * [1] Petr Mokrov, Alexander Korotin, Lingxiao Li, Aude Genevay, Justin M Solomon, and Evgeny Burnaev.  Large-scale wasserstein gradient flows.  Advances in Neural Information Processing Systems, 34:15243–15256, 2021.
  * [2] Charlie Frogner and Tomaso Poggio.  Approximate inference with wasserstein gradient flows.  In International Conference on Artificial Intelligence and Statistics, pages 2581–2590. PMLR, 2020.
  * [3] Soroosh Shafieezadeh Abadeh, Viet Anh Nguyen, Daniel Kuhn, and Peyman M Mohajerin Esfahani.  Wasserstein distributionally robust kalman filtering.  Advances in Neural Information Processing Systems, 31, 2018.
  * [4] Daniel Danu, Thia Kirubarajan, and Thomas Lang.  Wasserstein distance for the fusion of multisensor multitarget particle filter clouds.  In 2009 12th International Conference on Information Fusion, pages 25–32. IEEE, 2009.
  * [5] Sebastian Reich and Colin J Cotter.  Ensemble filter techniques for intermittent data assimilation.  Large Scale Inverse Problems: Computational Methods and Applications in the Earth Sciences, 13:91–134, 2013.
  * [6] Sebastian Reich.  A nonparametric ensemble transform method for bayesian inference.  SIAM Journal on Scientific Computing, 35(4):A2013–A2024, 2013.
  * [7] Niladri Das, Vedang Deshpande, and Raktim Bhattacharya.  Optimal-transport-based tracking of space objects using range data from a single ranging station.  Journal of Guidance, Control, and Dynamics, 42(6):1237–1249, 2019.
  * [8] Niladri Das, Riddhi Pratim Ghosh, Nilabja Guha, Raktim Bhattacharya, and Bani Mallick.  Optimal transport based tracking of space objects in cylindrical manifolds.  The Journal of the Astronautical Sciences, 66(4):582–606, 2019.
  * [9] Amirhossein Taghvaei and Prashant G Mehta.  An optimal transport formulation of the ensemble kalman filter.  IEEE Transactions on Automatic Control, 66(7):3052–3067, 2020.
  * [10] Rahul Singh, Isabel Haasler, Qinsheng Zhang, Johan Karlsson, and Yongxin Chen.  Inference with aggregate data in probabilistic graphical models: An optimal transport approach.  IEEE Transactions on Automatic Control, 67(9):4483–4497, 2022.
  * [11] Yongxin Chen, Tryphon T Georgiou, and Michele Pavon.  Optimal transport in systems and control.  Annual Review of Control, Robotics, and Autonomous Systems, 4:89–113, 2021.
  * [12] Amirhossein Taghvaei and Bamdad Hosseini.  An optimal transport formulation of bayes’ law for nonlinear filtering algorithms.  In 2022 IEEE 61st Conference on Decision and Control (CDC), pages 6608–6613. IEEE, 2022.
  * [13] Adrien Corenflos, James Thornton, George Deligiannidis, and Arnaud Doucet.  Differentiable particle filtering via entropy-regularized optimal transport.  In International Conference on Machine Learning, pages 2100–2111. PMLR, 2021.
  * [14] Yongxin Chen, Johan Karlsson, and Axel Ringh.  Optimal transport for applications in control and estimation.  IEEE Control Systems Magazine, 41(4):28–33, 2021.
  * [15] C. Villani.  Topics in optimal transportation, volume 58.  Amer Mathematical Society, 2003.
  * [16] A. Halder and R. Bhattacharya.  Further results on probabilistic model validation in wasserstein metric.  In 51st IEEE Conference on Decision and Control, Maui, 2012.
  * [17] Sadri Hassani.  Mathematical physics, a modern introduction to its foundations, 1999.
  * [18] I. Bilik and J. Tabrikian.  Optimal recursive filtering using gaussian mixture model.  In IEEE/SP 13th Workshop on Statistical Signal Processing, 2005, pages 399–404, 2005.
  * [19] Kooktae Lee, Abhishek Halder, and Raktim Bhattacharya.  Performance and robustness analysis of stochastic jump linear systems using wasserstein metric.  Automatica, 51:341–347, 2015.
  * [20] Anderson B. and Moore J.  Optimal filtering.  Prentice Hall, 1979.