# Splitting of elliptic flow in a tilted fireball

Tribhuban Parida<EMAIL_ADDRESS>Sandeep Chatterjee
<EMAIL_ADDRESS>Department of Physical Sciences,
Indian Institute of Science Education and Research Berhampur,
Transit Campus (Govt ITI), Berhampur-760010, Odisha, India

###### Abstract

The splitting of elliptic flow measured in different regions of the momentum
space of produced hadrons has been recently studied in transport models and
proposed as a sensitive probe of the angular momentum carried by the fireball
produced in a relativistic heavy ion collision. The initial state angular
momentum also gives rise to rapidity odd directed flow which has been
measured. We consider a relativistic hydrodynamic framework with the initial
matter distribution suitably calibrated to describe the observed directed flow
and apply it to study the spilt in the elliptic flow. Our study suggests that
the split in the elliptic flow is mostly driven by directed and triangular
flows and may be used to constrain models of initial state rapidity
distribution of matter in the fireball.

##  I Introduction

The system of two non-central colliding relativistic heavy ion nuclei carries
large angular momentum. In the aftermath of the collision, a part of this
angular momentum is deposited in the locally thermalised fireball. The
hydrodynamic response results in fluid vorticity and possibly spin
polarisation which may be finally observed in the phase space occupation and
polarisation states of the emitted particles Liang and Wang (2005a, b);
Becattini  _et al._ (2008); Betz  _et al._ (2007); Ipp  _et al._ (2008);
Becattini  _et al._ (2013, 2015); Pang  _et al._ (2016a); Karpenko and
Becattini (2017).

There have been several attempts to model the initial longitudinal
distribution of various hydrodynamic fields like the energy density and fluid
velocity which after hydrodynamic evolution and particlization can leave their
imprint on different observables Snellings  _et al._ (2000); Betz  _et al._
(2007); Becattini  _et al._ (2013, 2015); Pang  _et al._ (2016a); Karpenko and
Becattini (2017); Bozek and Wyskiel (2010); Bozek  _et al._ (2011); Bożek  _et
al._ (2015); Bozek and Broniowski (2016); Bozek  _et al._ (2015); Broniowski
and Bożek (2016); Chatterjee and Bozek (2017); Bozek and Broniowski (2018);
Pang  _et al._ (2016b, 2015, 2018); Wu  _et al._ (2018); Shen and Alzhrani
(2020); Ryu  _et al._ (2021); Jiang  _et al._ (2021, 2022). While for quite
some time the rapidity dependence of directed flow has been used to
discriminate such models of the initial three dimensional matter distribution
Snellings  _et al._ (2000); Bozek and Wyskiel (2010); Shen and Alzhrani
(2020); Ryu  _et al._ (2021); Jiang  _et al._ (2021, 2022), recently it has
been demonstrated that even polarisation measurements of the final state
hadrons can constrain such ansatz of initial matter distribution Ryu  _et al._
(2021); Alzhrani  _et al._ (2022).

It has been pointed out that the non-zero angular momentum of the fireball
results in the splitting of elliptic flow in the momentum space of the final
state hadrons Chen  _et al._ (2021). This makes the splitting of elliptic flow
a sensitive observable to constrain models of the initial state three
dimensional matter distribution of the fireball. In a subsequent work, it was
argued that this splitting is mainly driven by the directed flow Zhang and Lin
(2021). These studies were conducted in a transport model framework and they
don’t describe the data on directed flow Abelev  _et al._ (2008). This raises
concern on their model prediction of the elliptic flow splitting.

The various models of the longitudinal profile of the fireball that has been
explored so far can be broadly divided into two categories: shifted initial
condition (SIC) Hirano and Tsuda (2002) and tilted initial condition (TIC)
Bozek and Wyskiel (2010). In SIC, the rapidity profile at each point on the
transverse plane is shifted according to the local centre of mass rapidity.
TIC is inspired from the ansatz that a participant nucleon deposits more
energy along its direction of motion Brodsky  _et al._ (1977); Back  _et al._
(2003); Bialas and Czyz (2005); Adil and Gyulassy (2005); Armesto  _et al._
(2007); Bzdak and Wozniak (2010). It has been shown that SIC fails to describe
the rapidity slope of directed flow at mid-rapidity in Au+Au at
$\sqrt{s_{\textrm{NN}}}=200$ GeV while TIC succeeds to describe the
experimental data on directed flow, albeit with a free parameter $\eta_{m}$
that parametrises the forward-backward asymmetry in the energy deposition of a
participant in the initial state Bozek and Wyskiel (2010). In this study, we
discuss the contrasting nature of the splitting of the elliptic flow in
different regions of the momentum space for both SIC as well as TIC,
underlying the significance of this observable in our efforts to comprehend
the longitudinal dynamics of the fireball.

##  II INITIAL RAPIDITY PROFILE

We study Au+Au collisions at $\sqrt{s_{\textrm{NN}}}=200$ GeV. The initial
condition for the hydrodynamic evolution of the fireball is obtained from the
Optical Glauber model. Here, the nucleus is modelled as a Woods-Saxon
distribution $\rho(x,y,z)=\frac{\rho_{0}}{1+\exp{(\frac{r-R}{a}})}$ where $\
r=\sqrt{x^{2}+y^{2}+z^{2}}$, $R=6.38$ fm and $a=0.535$ fm Shou  _et al._
(2015). The z-axis is taken along the beam axis while the x-axis is along the
impact parameter direction. The nuclear thickness function $T(x,y)$ is
obtained as

$T(x,y)=\int\rho(x,y,z)dz$ (1)

using which one can define the total forward ($N_{+}$) and backward ($N_{-}$)
going participants along the beam axis at a point $(x,y)$ on the transverse
plane

$N_{+}(x,y)=T(x-b/2,y)\left(1-(1-\sigma_{NN}T(x+b/2,y))\right)$ (2)
$N_{-}(x,y)=T(x+b/2,y)\left(1-(1-\sigma_{NN}T(x-b/2,y))\right)$ (3)

We have compared between two models of $\epsilon(x,y,\eta_{s})$, the initial
energy density deposited at a constant $\tau$ hypersurface at
$(x,y,\eta_{s})$: SIC Hirano and Tsuda (2002) and TIC Bozek and Wyskiel
(2010).

In case of SIC, the following ansatz is adopted for $\epsilon(x,y,\eta_{s})$

$\displaystyle\epsilon(x,y,\eta_{s})$ $\displaystyle=$
$\displaystyle\epsilon_{0}\left[\left(N_{+}(x,y)+N_{-}(x,y)\right)\frac{(1-\alpha)}{2}\right.$
(4)
$\displaystyle\left.+N_{coll}(x,y)\alpha\right]\times\epsilon_{\eta_{s}}(\eta_{s}-\eta_{sh}(x,y))$

where $\epsilon_{\eta_{s}}(\eta_{s}-\eta_{sh}(x,y))$ gives the $\eta_{s}$
distribution at $(x,y)$

$\epsilon_{\eta_{s}}(\eta_{s})=\exp\left(-\frac{\left(|\eta_{s}|-\eta_{0}\right)^{2}}{2\sigma_{\eta}^{2}}\theta(|\eta_{s}|-\eta_{0})\right)$
(5)

We have used $\epsilon_{0}=13.2~{}\text{GeV/fm}^{3}$, $\alpha=0.14$,
$\eta_{0}=1.3$ and $\sigma_{\eta}=1.5$ that provides a good description of the
($\eta-\frac{dN_{ch}}{d\eta}$) data. $\eta_{sh}(x,y)$ is given by

$\eta_{sh}=\frac{1}{2}\ln\frac{N_{+}(x,y)+N_{-}(x,y)+v_{N}(N_{+}(x,y)-N_{-}(x,y))}{N_{+}(x,y)+N_{-}(x,y)-v_{N}(N_{+}(x,y)-N_{-}(x,y))}$
(6)

Here, $v_{N}$ is the initial longitudinal velocity of each nucleon with mass
$m_{N}$.

$v_{N}=\sqrt{1-(4m_{N}^{2})/(s_{NN})}$ (7)

The second initial condition that we have studied is the TIC. In this case
$\epsilon(x,y,\eta_{s})$ is given by

$\displaystyle\epsilon(x,y,\eta_{s})$ $\displaystyle=$
$\displaystyle\epsilon_{0}\left[\left(N_{+}(x,y)f_{+}(\eta_{s})+N_{-}(x,y)f_{-}(\eta_{s})\right)\right.$
(8)
$\displaystyle\left.\times\left(1-\alpha\right)+N_{coll}(x,y)\epsilon_{\eta_{s}}\left(\eta_{s}\right)\alpha\right]$

where $\epsilon_{\eta_{s}}(\eta_{s})$ is the rapidity even profile as given in
Eq. 5 and $f_{+,-}(\eta_{s})$ introduce rapidity odd component in $\epsilon$

$f_{+,-}(\eta_{s})=\epsilon_{\eta_{s}}(\eta_{s})\epsilon_{F,B}(\eta_{s})$ (9)

where

$\epsilon_{F}(\eta_{s})=\begin{cases}0,&\text{if }\eta_{s}<-\eta_{m}\\\
\frac{\eta_{s}+\eta_{m}}{2\eta_{m}},&\text{if
}-\eta_{m}\leq\eta_{s}\leq\eta_{m}\\\ 1,&\text{if
}\eta_{m}<\eta_{s}\end{cases}$ (10)

and

$\epsilon_{B}(\eta_{s})=\epsilon_{F}(-\eta_{s})$ (11)

We have used $\eta_{m}=2.5$ to describe the directed flow data Abelev  _et
al._ (2008). For both the initial conditions, we have assumed the Bjorken flow
ansatz

$u^{\mu}(\tau_{0},x,y,\eta_{s})=(\cosh{\eta_{s}},0,0,\sinh{\eta_{s}})$

We evolve the above deposited initial energy distribution with the publicly
available MUSIC code Schenke  _et al._ (2010, 2012); Denicol  _et al._ (2018);
Paquet  _et al._ (2016) which implements evolution within the framework of 3+1
D relativistic hydrodynamics followed by Cooper Frye freezeout at $T=150$ MeV
and finally allowing all the resonances to decay to stable hadrons under
strong interaction. Thus, we obtain the momentum space probability
distribution of hadrons using which we compute various observables. We
consider the lattice QCD based equation of state, NEoS-B at zero baryon
density Monnai  _et al._ (2019); Bazavov  _et al._ (2014, 2012); Ding  _et
al._ (2015) and take the shear viscosity $\eta$ to entropy density $s$ ratio,
$\eta/s=0.08$. We have ignored the effects of bulk viscosity.

##  III SPLITTING OF THE ELLIPTIC FLOW

The azimuthal distribution of the hadrons in the plane transverse to the beam
axis can be expanded into Fourier components in the following way:

$\frac{dN}{d\phi}=\frac{1}{2\pi}\left(1+2\sum_{n}\left(v_{n}\cos(n(\phi-\psi_{\text{RP}}))+s_{n}\sin(n(\phi-\psi_{\text{RP}}))\right)\right)$
(12)

where $\psi_{\text{RP}}$ is the reaction plane angle in the laboratory frame.
$v_{n}$ and $s_{n}$ are the Fourier coefficients that characterise the
distribution.

There has been a recent proposal to measure the split $\Delta v_{2}$ in
$v_{2}=\langle\cos(2(\phi-\psi_{\text{RP}}))\rangle$, as measured in different
regions of the final hadron momentum space

$\Delta v_{2}={v_{2}}^{\text{R}}-{v_{2}}^{\text{L}}$ (13)

where
${v_{2}}^{\text{R}}=\langle\cos(2(\phi^{\text{R}}-\psi_{\text{RP}}))\rangle$
with $\phi^{\text{R}}\in((\psi_{\text{RP}}-\pi/2),(\psi_{\text{RP}}+\pi/2))$
and
${v_{2}}^{\text{L}}=\langle\cos(2(\phi^{\text{L}}-\psi_{\text{RP}}))\rangle$
with $\phi^{\text{L}}\in((\psi_{\text{RP}}+\pi/2),(\psi_{\text{RP}}+3\pi/2))$.
Here, $\langle...\rangle$ refers to averaging over the phase space of the
produced hadrons.

$\psi_{\text{RP}}$ is not directly measurable in experiments. The second order
event plane orientation $\psi_{2}$ and the first order spectator plane
$\psi_{\text{SP}}$ have been proposed as good proxies for $\psi_{\text{RP}}$
Adamczyk  _et al._ (2013); Abelev  _et al._ (2013). However, for the
determination of $\Delta v_{2}$, $\psi_{\text{SP}}$ alone is suitable as
$\psi_{2}=\pi$ is identified with $\psi_{2}=0$ and hence does not distinguish
between the phase spaces associated with $\phi^{\text{R}}$ and
$\phi^{\text{L}}$. Recently, the Event Plane Detector has been installed at
large rapidities which can also be used to estimate $\psi_{\text{RP}}$Adams
_et al._ (2020).

${v_{2}}^{\text{R}}$ and ${v_{2}}^{\text{L}}$ work out to be the following
Zhang and Lin (2021):

$\displaystyle v_{2}^{\text{R}}$ $\displaystyle=$
$\displaystyle\frac{\int_{\psi_{\text{RP}}-\frac{\pi}{2}}^{\psi_{\text{RP}}+\frac{\pi}{2}}\
\cos(2(\phi-\psi_{\text{RP}}))\frac{dN}{d\phi}\
d\phi}{\int_{\psi_{\text{RP}}-\frac{\pi}{2}}^{\psi_{\text{RP}}+\frac{\pi}{2}}\frac{dN}{d\phi}d\phi}$
(14) $\displaystyle\approx$
$\displaystyle\frac{v_{2}+\frac{4v_{1}}{3\pi}+\frac{12v_{3}}{5\pi}-\frac{20v_{5}}{21\pi}}{1+\frac{4v_{1}}{\pi}-\frac{4v_{3}}{3\pi}+\frac{4v_{5}}{5\pi}}$

$\displaystyle v_{2}^{\text{L}}$ $\displaystyle=$
$\displaystyle\frac{\int_{\psi_{\text{RP}}+\frac{\pi}{2}}^{\psi_{\text{RP}}+\frac{3\pi}{2}}\
\cos(2(\phi-\psi_{\text{RP}}))\frac{dN}{d\phi}\
d\phi}{\int_{\psi_{\text{RP}}+\frac{\pi}{2}}^{\psi_{\text{RP}}+\frac{3\pi}{2}}\frac{dN}{d\phi}d\phi}$
(15) $\displaystyle\approx$
$\displaystyle\frac{v_{2}-\frac{4v_{1}}{3\pi}-\frac{12v_{3}}{5\pi}+\frac{20v_{5}}{21\pi}}{1-\frac{4v_{1}}{\pi}+\frac{4v_{3}}{3\pi}-\frac{4v_{5}}{5\pi}}$

In Eqs. 14 and 15 we have omitted contributions from harmonics higher than the
fifth order. Further, from Eqs. 13, 14 and 15 we get

$\Delta
v_{2}\approx\frac{8v_{1}}{3\pi}+\frac{24v_{3}}{5\pi}-\frac{40v_{5}}{21\pi}$
(16)

In Eq. 16 only terms which are linear in the flow harmonics have been shown as
they are sufficient to estimate $\Delta v_{2}$ . Thus, $\Delta v_{2}$ is
sourced mainly by the odd flow harmonics. The collision geometry is such that
these odd flow harmonics have rapidity odd components with respect to the
reaction plane. Hence, $\Delta v_{2}$ is also rapidity odd following the odd
flow harmonics which is also evident from Eq. 13. There has been measurement
of directed flow with respect to $\psi_{\text{SP}}$ Abelev  _et al._ (2008).
Optical Glauber model with TIC followed by hydrodynamic expansion is able to
describe the $v_{1}$ measurement at $\sqrt{s_{\textrm{NN}}}=200$ GeV for Au+Au
collisionsBozek and Wyskiel (2010). Here, we use similar TIC within an optical
Glauber model followed by hydrodynamic expansion to compute the model
expectation for $\Delta v_{2}$. We expect similar results on including
fluctuations in the initial condition Bozek (2012).

##  IV RESULTS

Figure 1: (Color online) Phase space dependence of $v_{1}$ is computed both
for tilted initial condition (red solid line) and shifted initial condition
(blue dashed line) for 5-40 $\%$ centrality Au+Au collisions at
$\sqrt{s_{\textrm{NN}}}=200$ GeV. The model expectations are compared to
measurements from the STAR collaboration Abelev  _et al._ (2008). $v_{1}$ vs
$\eta$ is shown in panel (a) and $v_{1}$ vs $p_{T}$ is shown in panel (b).

We will now present the prediction for $\Delta v_{2}$ with respect to the
reaction plane as computed with TIC as well as SIC. As seen in Eq. 16, the
leading contributions to $\Delta v_{2}$ arise from the odd harmonics, out of
which there are STAR measurements on $v_{1}$ with respect to the spectator
plane Abelev  _et al._ (2008). Thus, we first compare the model results with
the STAR data for $v_{1}$ in Fig. 1. We have plotted the model expectations
for $v_{1}-\eta$ and $v_{1}-p_{T}$ in the panels (a) and (b) respectively and
compared them to the STAR measurements Abelev  _et al._ (2008). The model
results have been shown for both TIC (red solid line) as well as those from
SIC (dashed blue line). We find that for both $v_{1}-\eta$ and $v_{1}-p_{T}$,
TIC is able to describe the STAR data well while SIC fails. This is in
agreement with earlier studies Bozek and Wyskiel (2010); Jiang  _et al._
(2022). Thus, we expect TIC to provide correct prediction of $\Delta v_{2}$.

Figure 2: (Color online) The prediction for $\Delta v_{2}$ vs $\eta$ with tilt
initial condition has been plotted in red solid line. Further, the first three
leading flow harmonics that contribute to $\Delta v_{2}$ (see Eq. 16) are
shown as well.

The good description of $v_{1}$ by TIC motivates us further to compute $\Delta
v_{2}$ vs $\eta$ within the same scheme. In Fig. 2 we have shown the model
predictions for the $\eta$ dependence of $\Delta v_{2}$. Further, we have also
plotted the odd harmonics along with appropriate coefficients as suggested by
Eq. 16 that are the dominant contributors to $\Delta v_{2}$. Firstly, we note
that $\Delta v_{2}$ arises as a competition between $\frac{8v_{1}}{3\pi}$ and
$\frac{24v_{3}}{5\pi}$ as they are of opposite signs. The
$\frac{8v_{1}}{3\pi}$ term marginally wins and hence $\Delta v_{2}$ follows
its sign. For $|\eta|<1$, $\Delta v_{2}\sim 10^{-4}$ and it grows
substantially for larger $\eta$. The contribution from $v_{5}$ is negligible
in the entire $\eta$ range.

Figure 3: (Color online) The prediction for $\Delta v_{2}$ vs $p_{T}$ with
tilt initial condition has been plotted in red solid line. Further, the first
three leading flow harmonics that contribute to $\Delta v_{2}$ (see Eq. 16)
are shown as well.

The $p_{T}$ dependence is shown in Fig. 3. The predictions for $0<\eta<1.3$
and $1.3<\eta<2.5$ are shown separately in panels (a) and (b) respectively. As
expected from the $\eta$ dependence, the magnitude of the split is larger in
the $1.3<\eta<2.5$, otherwise, the trends are similar. $\Delta v_{2}$ is
negative for small $p_{T}$. There is a turning point around $p_{T}\sim 0.4$
GeV after which it crosses zero at $p_{T}\sim 1$ GeV and monotonically rises.
This $p_{T}$ dependence of $\Delta v_{2}$ is mainly borrowed from $v_{1}$
which also has a similar trend unlike $v_{3}$ which remains positive for all
$p_{T}$. It is interesting to note that for $p_{T}>1.5$ GeV, $v_{3}$ becomes
the dominant contributor.

Figure 4: (Color online) The predictions for phase space dependence of $\Delta
v_{2}$ has been compared between tilt initial condition and shift initial
condition. Further, the first three leading flow harmonics that contribute to
$\Delta v_{2}$ are also shown for each model to understand the origin in the
difference of their phase space dependence of $\Delta v_{2}$. $\Delta v_{2}$
vs $\eta$ is shown in panel (a) and $\Delta v_{2}$ vs $p_{T}$ is shown in
panel (b).

Further, we present a comparative study of $\Delta v_{2}$ computed with TIC
and SIC in Fig. 4. The results on the $\eta$ dependence of $\Delta v_{2}$ are
presented in panel (a) while the $p_{T}$ dependence is plotted in panel (b).
As seen in panel (a), for $|\eta|<1$, $\Delta v_{2}$ is around 10-20 times
larger in SIC as compared to TIC. Further they are of opposite signs- at
positive rapidities while the SIC gives a positive $\Delta v_{2}$, the TIC
yields a negative $\Delta v_{2}$. This may be traced to the fact that $v_{1}$
is of opposite sign in the two models. Further, from panel (b) we note that
the origin of this opposite sign is from the low $p_{T}$ region as for
$p_{T}>1.5$ GeV, both models yield positive values. While we noted earlier
that $\Delta v_{2}$ results out of tension between $v_{1}$ and $v_{3}$ in the
TIC, in the SIC we find that it is mostly controlled by $v_{1}$ as $v_{3}$
turns out to be small in this case. It is worth noting here that an earlier
study of $\Delta v_{2}$ within the transport framework of the AMPT model had
similar conclusions as the results here with SIC Zhang and Lin (2021). Thus,
the characteristics of $\Delta v_{2}$ can serve as a sensitive probe of
rapidity dependent initial condition of the fireball.

Figure 5: (Color online) The predictions for phase space dependence of $\Delta
v_{2}/v_{2}$ has been shown for tilt initial condition. Panel (a) shows the
$\eta$ dependence while panel (b) shows the $p_{T}$ dependence.

It is well known that the various model systematics of the transverse initial
condition can result in the variation of $v_{2}$ Shen  _et al._ (2011); Qiu
and Heinz (2011); Ruggieri  _et al._ (2013). Such model dependencies can creep
into our predictions of $\Delta v_{2}$ as well. We suggest to scale $\Delta
v_{2}$ with $v_{2}$ in order to cancel out such systematics of the transverse
initial condition which are not of interest here. We present the results of
$\Delta v_{2}$ scaled by $v_{2}$ in Fig. 5. The $\eta$ dependence is plotted
in panel (a). $|\Delta v_{2}/v_{2}|$ stays below 0.005 for $|\eta|<1$ beyond
which it has a rapid linear growth. The $p_{T}$ differential values are
plotted in panel (b). We obtain about $3\%$ and $5\%$ $\Delta v_{2}/v_{2}$ at
$p_{T}\sim 1.5$ GeV for $0<\eta<1.3$ and $1.3<\eta<2.5$ respectively. The
ratio of phase space integrated $\Delta v_{2}$ and $v_{2}$, $\Delta
v_{2}/v_{2}$ comes out to be -0.0019 and -0.0268 for $0<\eta<1.3$ and
$1.3<\eta<2.5$ respectively.

##  V SUMMARY

The collision geometry of a non-central relativistic heavy ion collision
introduces a large angular momentum in its initial state. Rapidity odd
directed flow is a natural consequence of this. Such observables probe the
longitudinal profile of the fireball. Recently, it has been observed that such
large angular momentum in the initial state also causes a split in the
magnitude of $v_{2}$, $\Delta v_{2}$ in different regions of the final state
hadron momentum space - parallel and anti-parallel to the impact parameter
direction Chen  _et al._ (2021); Zhang and Lin (2021). $\Delta v_{2}$ has been
proposed to be a sensitive probe of the initial rapidity profile of the
fireball. A transport model framework was adopted in these works that do not
describe the $v_{1}$ data.

In this work, we have revisited the estimation of $\Delta v_{2}$ within a 3+1
D relativistic hydrodynamic framework with tilted initial condition that is
known to describe the $v_{1}$ data Bozek and Wyskiel (2010). We find $v_{1}$
to be the leading contributor to $\Delta v_{2}$, as was reported in earlier
studies. However, unlike in those transport model based studies, we find that
the tilted initial condition gives rise to sizeable rapidity odd $v_{3}$ which
also contributes significantly to $\Delta v_{2}$, particularly for $p_{T}>1.5$
GeV, it becomes the dominant contributor to $\Delta v_{2}$. In order to
demonstrate the sensitivity of $\Delta v_{2}$ to the choice of the initial
condition, we have computed $\Delta v_{2}$ also with shifted initial
condition. In this case, similar to the earlier studies, $v_{3}$ comes out to
be negligible and hence $\Delta v_{2}$ gets contribution dominantly from
$v_{1}$. The $\eta$ dependence of $\Delta v_{2}$ is of opposite sign for
titled versus shifted initial conditions owing to the opposite signs of their
$v_{1}$. Finally, we have also presented the ratio of $\Delta v_{2}$ to
$v_{2}$ so as to get rid of the various model uncertainties that affect
$v_{2}$ and hence also $\Delta v_{2}$: for $p_{T}\sim 1.5$ GeV we obtain
$\Delta v_{2}/v_{2}\sim 3\%$ and $5\%$ for $0<\eta<1.3$ and $1.3<\eta<2.5$
respectively. Our study demonstrates that $\Delta v_{2}$ can play
complementary role to $v_{1}$ in constraining the rapidity profile of the
fireball.

##  VI ACKNOWLEDGEMENTS

SC acknowledges helpful discussions with Piotr Bozek and IISER Berhampur for
Seed Grant.

## References

  * Liang and Wang (2005a) Z.-T. Liang and X.-N. Wang, Phys. Rev. Lett. 94, 102301 (2005a), [Erratum: Phys.Rev.Lett. 96, 039901 (2006)], arXiv:nucl-th/0410079 .
  * Liang and Wang (2005b) Z.-T. Liang and X.-N. Wang, Phys. Lett. B 629, 20 (2005b), arXiv:nucl-th/0411101 .
  * Becattini  _et al._ (2008) F. Becattini, F. Piccinini, and J. Rizzo, Phys. Rev. C 77, 024906 (2008), arXiv:0711.1253 [nucl-th] .
  * Betz  _et al._ (2007) B. Betz, M. Gyulassy, and G. Torrieri, Phys. Rev. C 76, 044901 (2007), arXiv:0708.0035 [nucl-th] .
  * Ipp  _et al._ (2008) A. Ipp, A. Di Piazza, J. Evers, and C. H. Keitel, Phys. Lett. B 666, 315 (2008), arXiv:0710.5700 [hep-ph] .
  * Becattini  _et al._ (2013) F. Becattini, L. Csernai, and D. J. Wang, Phys. Rev. C 88, 034905 (2013), [Erratum: Phys.Rev.C 93, 069901 (2016)], arXiv:1304.4427 [nucl-th] .
  * Becattini  _et al._ (2015) F. Becattini, G. Inghirami, V. Rolando, A. Beraudo, L. Del Zanna, A. De Pace, M. Nardi, G. Pagliara, and V. Chandra, Eur. Phys. J. C 75, 406 (2015), [Erratum: Eur.Phys.J.C 78, 354 (2018)], arXiv:1501.04468 [nucl-th] .
  * Pang  _et al._ (2016a) L.-G. Pang, H. Petersen, Q. Wang, and X.-N. Wang, Phys. Rev. Lett. 117, 192301 (2016a), arXiv:1605.04024 [hep-ph] .
  * Karpenko and Becattini (2017) I. Karpenko and F. Becattini, Eur. Phys. J. C 77, 213 (2017), arXiv:1610.04717 [nucl-th] .
  * Snellings  _et al._ (2000) R. J. M. Snellings, H. Sorge, S. A. Voloshin, F. Q. Wang, and N. Xu, Phys. Rev. Lett. 84, 2803 (2000), arXiv:nucl-ex/9908001 .
  * Bozek and Wyskiel (2010) P. Bozek and I. Wyskiel, Phys. Rev. C 81, 054902 (2010), arXiv:1002.4999 [nucl-th] .
  * Bozek  _et al._ (2011) P. Bozek, W. Broniowski, and J. Moreira, Phys. Rev. C 83, 034911 (2011), arXiv:1011.3354 [nucl-th] .
  * Bożek  _et al._ (2015) P. Bożek, W. Broniowski, and A. Olszewski, Phys. Rev. C 91, 054912 (2015), arXiv:1503.07425 [nucl-th] .
  * Bozek and Broniowski (2016) P. Bozek and W. Broniowski, Phys. Lett. B 752, 206 (2016), arXiv:1506.02817 [nucl-th] .
  * Bozek  _et al._ (2015) P. Bozek, W. Broniowski, and A. Olszewski, Phys. Rev. C 92, 054913 (2015), arXiv:1509.04124 [nucl-th] .
  * Broniowski and Bożek (2016) W. Broniowski and P. Bożek, Phys. Rev. C 93, 064910 (2016), arXiv:1512.01945 [nucl-th] .
  * Chatterjee and Bozek (2017) S. Chatterjee and P. Bozek, Phys. Rev. C 96, 014906 (2017), arXiv:1704.02777 [nucl-th] .
  * Bozek and Broniowski (2018) P. Bozek and W. Broniowski, Phys. Rev. C 97, 034913 (2018), arXiv:1711.03325 [nucl-th] .
  * Pang  _et al._ (2016b) L.-G. Pang, H. Petersen, G.-Y. Qin, V. Roy, and X.-N. Wang, Eur. Phys. J. A 52, 97 (2016b), arXiv:1511.04131 [nucl-th] .
  * Pang  _et al._ (2015) L.-G. Pang, G.-Y. Qin, V. Roy, X.-N. Wang, and G.-L. Ma, Phys. Rev. C 91, 044904 (2015), arXiv:1410.8690 [nucl-th] .
  * Pang  _et al._ (2018) L.-G. Pang, H. Petersen, and X.-N. Wang, Phys. Rev. C 97, 064918 (2018), arXiv:1802.04449 [nucl-th] .
  * Wu  _et al._ (2018) X.-Y. Wu, L.-G. Pang, G.-Y. Qin, and X.-N. Wang, Phys. Rev. C 98, 024913 (2018), arXiv:1805.03762 [nucl-th] .
  * Shen and Alzhrani (2020) C. Shen and S. Alzhrani, Phys. Rev. C 102, 014909 (2020), arXiv:2003.05852 [nucl-th] .
  * Ryu  _et al._ (2021) S. Ryu, V. Jupic, and C. Shen, Phys. Rev. C 104, 054908 (2021), arXiv:2106.08125 [nucl-th] .
  * Jiang  _et al._ (2021) Z.-F. Jiang, C. B. Yang, and Q. Peng,  (2021), 10.1103/PhysRevC.104.064903, arXiv:2111.01994 [hep-ph] .
  * Jiang  _et al._ (2022) Z.-F. Jiang, S. Cao, X.-Y. Wu, C. B. Yang, and B.-W. Zhang, Phys. Rev. C 105, 034901 (2022), arXiv:2112.01916 [hep-ph] .
  * Alzhrani  _et al._ (2022) S. Alzhrani, S. Ryu, and C. Shen, (2022), arXiv:2203.15718 [nucl-th] .
  * Chen  _et al._ (2021) Z. Chen, Z. Wang, C. Greiner, and Z. Xu, (2021), arXiv:2108.12735 [hep-ph] .
  * Zhang and Lin (2021) C. Zhang and Z.-W. Lin, (2021), arXiv:2109.04987 [nucl-th] .
  * Abelev  _et al._ (2008) B. I. Abelev _et al._ (STAR), Phys. Rev. Lett. 101, 252301 (2008), arXiv:0807.1518 [nucl-ex] .
  * Hirano and Tsuda (2002) T. Hirano and K. Tsuda, Phys. Rev. C 66, 054905 (2002), arXiv:nucl-th/0205043 .
  * Brodsky  _et al._ (1977) S. J. Brodsky, J. F. Gunion, and J. H. Kuhn, Phys. Rev. Lett. 39, 1120 (1977).
  * Back  _et al._ (2003) B. B. Back _et al._ , Phys. Rev. Lett. 91, 052303 (2003), arXiv:nucl-ex/0210015 .
  * Bialas and Czyz (2005) A. Bialas and W. Czyz, Acta Phys. Polon. B 36, 905 (2005), arXiv:hep-ph/0410265 .
  * Adil and Gyulassy (2005) A. Adil and M. Gyulassy, Phys. Rev. C 72, 034907 (2005), arXiv:nucl-th/0505004 .
  * Armesto  _et al._ (2007) N. Armesto, L. McLerran, and C. Pajares, Nucl. Phys. A 781, 201 (2007), arXiv:hep-ph/0607345 .
  * Bzdak and Wozniak (2010) A. Bzdak and K. Wozniak, Phys. Rev. C 81, 034908 (2010), arXiv:0911.4696 [hep-ph] .
  * Shou  _et al._ (2015) Q. Y. Shou, Y. G. Ma, P. Sorensen, A. H. Tang, F. Videbæk, and H. Wang, Phys. Lett. B 749, 215 (2015), arXiv:1409.8375 [nucl-th] .
  * Schenke  _et al._ (2010) B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 82, 014903 (2010), arXiv:1004.1408 [hep-ph] .
  * Schenke  _et al._ (2012) B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 85, 024901 (2012), arXiv:1109.6289 [hep-ph] .
  * Denicol  _et al._ (2018) G. S. Denicol, C. Gale, S. Jeon, A. Monnai, B. Schenke, and C. Shen, Phys. Rev. C 98, 034916 (2018), arXiv:1804.10557 [nucl-th] .
  * Paquet  _et al._ (2016) J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 93, 044906 (2016), arXiv:1509.06738 [hep-ph] .
  * Monnai  _et al._ (2019) A. Monnai, B. Schenke, and C. Shen, Phys. Rev. C 100, 024907 (2019), arXiv:1902.05095 [nucl-th] .
  * Bazavov  _et al._ (2014) A. Bazavov _et al._ (HotQCD), Phys. Rev. D 90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Bazavov  _et al._ (2012) A. Bazavov _et al._ (HotQCD), Phys. Rev. D 86, 034509 (2012), arXiv:1203.0784 [hep-lat] .
  * Ding  _et al._ (2015) H. T. Ding, S. Mukherjee, H. Ohno, P. Petreczky, and H. P. Schadler, Phys. Rev. D 92, 074043 (2015), arXiv:1507.06637 [hep-lat] .
  * Adamczyk  _et al._ (2013) L. Adamczyk _et al._ (STAR), Phys. Rev. C 88, 064911 (2013), arXiv:1302.3802 [nucl-ex] .
  * Abelev  _et al._ (2013) B. Abelev _et al._ (ALICE), Phys. Rev. Lett. 110, 012301 (2013), arXiv:1207.0900 [nucl-ex] .
  * Adams  _et al._ (2020) J. Adams _et al._ , Nucl. Instrum. Meth. A 968, 163970 (2020), arXiv:1912.05243 [physics.ins-det] .
  * Bozek (2012) P. Bozek, Phys. Rev. C 85, 034901 (2012), arXiv:1110.6742 [nucl-th] .
  * Shen  _et al._ (2011) C. Shen, S. A. Bass, T. Hirano, P. Huovinen, Z. Qiu, H. Song, and U. Heinz, J. Phys. G 38, 124045 (2011), arXiv:1106.6350 [nucl-th] .
  * Qiu and Heinz (2011) Z. Qiu and U. W. Heinz, Phys. Rev. C 84, 024911 (2011), arXiv:1104.0650 [nucl-th] .
  * Ruggieri  _et al._ (2013) M. Ruggieri, F. Scardina, S. Plumari, and V. Greco, Phys. Lett. B 727, 177 (2013), arXiv:1303.3178 [nucl-th] .


# Identification of Sagittarius stream members in Angular Momentum space with
Gaussian mixture techniques

Jorge Peñarrubia1,2, Michael S. Petersen1
1Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
2Centre for Statistics, University of Edinburgh, School of Mathematics,
Edinburgh EH9 3FD, UK<EMAIL_ADDRESS>
###### Abstract

This paper uses Gaussian mixture techniques to dissect the Milky Way (MW)
stellar halo in angular momentum space. Application to a catalogue of 5389
stars near the plane of the Sagittarius (Sgr) stream with full 6D phase-space
coordinates supplied by Gaia EDR3 and SEGUE returns four independent dynamical
components. The broadest and most populated corresponds to the ‘smooth’ MW
halo. The narrowest and faintest contains 40 stars of the Orphan stream. We
find a component with little or no angular momentum likely associated with the
GSE substructure. We also identify 925 stars and 7 Globular Clusters with
probabilities $>90\%$ to be members of the Sgr stream. Comparison against
$N$-body models shows that some of these members trace the continuation of the
leading/trailing tails in the Southern/Northern hemispheres. The new
detections span $\sim 800^{\circ}$ on the sky, thus wrapping the Galaxy twice.

###### keywords:

Galaxy: halo–galaxies: kinematics and dynamics–galaxies: evolution

##  1 Introduction

Since its serendipitous discovery by Ibata et al. (1994), the Sagittarius
(Sgr) dwarf galaxy has become a poster child for the tidal disruption of
satellite galaxies and the hierarchical build up of galactic stellar haloes.
With early observations largely restricted to pencil-beam fields of view
(Ivezić et al., 2000; Vivas et al., 2001; Newberg et al., 2002; Martínez-
Delgado et al., 2002), the striking extent of the tidal tails was not revealed
until all-sky 2MASS survey became available (Majewski et al., 2003). It was
soon realized that dynamical modelling of the tidal tails provides powerful
constraints on the gravitational field of the Milky Way (MW). Yet, after more
than two decades of continuous research the Sgr stream still presents a number
of theoretical challenges. For example, dynamical models of the Sgr stream
that adopt a static MW potential have returned dark matter halo with prolate
(Helmi, 2004), oblate (Johnston et al., 2005) and triaxial (Law & Majewski,
2010) shapes. This mismatch may be due to the gravitational attraction induced
by a massive Large Magellanic Clound (LMC) on first infall, which has two
important effects: it perturbs the orbits of stream stars (Law & Majewski,
2010; Vera-Ciro & Helmi, 2013), and displaces the MW disc from the Galactic
barycentre (Gómez et al., 2015; Petersen & Peñarrubia, 2020), causing a reflex
motion in the kinematics of distant halo stars (Erkal et al., 2020; Petersen &
Peñarrubia, 2021, hereafter PP21) and Local Group galaxies (Peñarrubia et al.,
2016). Recently, Vasiliev et al. (2021, hereafter VBE21) show that the
disparate constraints on the dark matter halo shapes may be due to the
adoption of a static MW potential in earlier models. However, other issues,
like the origin of the bifurcated tails (Belokurov et al., 2006; Koposov et
al., 2012), remain poorly understood. E.g. it was suggested that the
bifurcations could be caused by rotation in the progenitor dwarf (Peñarrubia
et al., 2010). Yet, these models predict a significant amount of rotation in
the remnant core, which has been ruled out by different kinematic surveys
(Peñarrubia et al., 2011; Frinchaboy et al., 2012; Vasiliev & Belokurov, 2020;
del Pino et al., 2021). Furthermore, all published models predict the
existence of old stream wraps populated by stars lost during early pericentric
passages whose detection remains elusive, probably because they are more
phase-mixed than those stripped recently. Finding this material may constrain
the amount of dynamical friction experienced by the Sgr dwarf (Fardal et al.,
2019b), and its infall trajectory (Dierickx & Loeb, 2017).

This paper explores a statistical method to detect kinematic substructures in
the halo that are partially, or even fully phase mixed. In particular, we use
Gaussian mixture techniques to quantify clustering in angular momentum space.
Previous attempts to detect the Sgr stream using catalogues of halo stars with
full phase-space coordinates either apply cuts in energy and/or angular
momentum (Li et al., 2019; Johnson et al., 2020), or search for stars on
similar orbits with Friends-of-Friends algorithms (Yang et al., 2019). In
contrast, our method draws statistical associations between objects in the
halo without relying on assumptions on the form, shape or time-evolution of
the MW potential, or introducing subjective cuts in integrals of motion. Its
applicability is chiefly limited by the reduced number of halo stars with
measured 6D phase-space coordinates. For example, Ramos et al. (2020) provides
a catalogue of 182495 RR Lyrae in Gaia DR2 (Gaia Collaboration et al., 2018)
with 5D information, that is approximately 17 times larger than the 6D data
set introduced below.

Figure 1: Angular momentum components of stellar halo stars with available 6D
phase-space coordinates and located within $|B_{\rm sgr}|<20^{\circ}$ of the
Sgr dSph orbital plane (see text). Stars are colour-coded according to their
probability to belonging to the Sgr stream. Back stars mark the angular
momentum components of the Sgr dwarf galaxy. Black/white circles show the
exclusion limits applied by PP21 to remove Sgr stream stars from their
analysis. Stars on the left of the dotted line were labelled Sgr stream
members by Johnson et al. (2020). Orphan stream members ($p_{\rm mem}>0.9$)
are shown in green. Magenta crosses in the bottom-left of each panel show the
average uncertainty of angular momentum measurements, $\sim 600\,{\rm
kpc}\,{\rm km\,s^{-1}}$. Right panel shows the angular momentum distribution
of the stellar sample (black line). Coloured histograms weigh stars by
membership probabilities (see text).

##  2 Data

We use a parent catalogue of 6030 K Giants and 3952 Blue Horizontal Branch
(BHB) stars in the Milky Way stellar halo with available phase-space
information compiled by PP21 from public data sets (Xue et al., 2008; Yanny et
al., 2009; Xue et al., 2011) and updated with Gaia EDR3 astrometry (Gaia
Collaboration et al., 2021).

One of the main goal of this paper is to identify Sgr stream members. For a
better characterization of the angular momentum distribution of the Sgr
stream, it is convenient to minimize the fraction of MW interlopers in the
Bayesian fits of §3. To this aim, the reference frame is rotated such that the
average polar angle measured from the Galactic centre is at $B_{\rm sgr}=0$
(Majewski et al., 2003). Subsequently, stars within $|B_{\rm sgr}|<20^{\circ}$
at Galactocentric distances $r>20\,{\rm kpc}$ are chosen, which yields 3419 K
Giants and 1970 BHBs, with mean distance errors $\epsilon_{D}/D\simeq 0.05$
and 0.17, respectively.

In addition, we also use published positions, distances, radial velocities,
and proper motions of MW satellites (McConnachie & Venn, 2020) and globular
clusters (Vasiliev & Baumgardt, 2021; Baumgardt & Vasiliev, 2021) in order to
explore possible associations with the Sagittarius dwarf. To convert
heliocentric into Galactocentric quantities we place the sun at
$\vec{x}_{\odot}=(-8.17,0.0,0.02)$ kpc (Gravity Collaboration et al., 2019;
Bennett & Bovy, 2019), with $\vec{v}_{\odot}=(-12.9,245.6,7.78)$ km s-1
(Drimmel & Poggio, 2018).

##  3 Bayesian analysis

An efficient way to model the distribution of stars in angular momentum
$\vec{L}=\vec{r}\times\vec{v}$, with $\vec{r}$ and $\vec{v}$ measured in a
Galactocentric frame, is to adopt a Gaussian mixture likelihood (Kuhn &
Feigelson, 2017)

$\displaystyle\mathcal{L}(\vec{\theta}|\vec{S})=\sum_{i=1}^{N}\,\alpha_{i}\,\mathcal{N}_{i}(\vec{\theta}|\vec{S}),$
(1)

with individual weights normalized such that $\sum_{i=1}^{N}\alpha_{i}=1$.
Here, $\vec{\theta}=\\{D_{i},\ell_{i},b_{i},v_{\rm
los,i},\mu_{\ell,i},\mu_{b,i}\,\epsilon_{D,i},\epsilon_{v,i},\epsilon_{\ell,i},\epsilon_{b,i}\\}^{N_{\rm
sample}}_{i=1}$ is a sample stars with full phase-space information, and
$\vec{S}$ is the array of model parameters that define the multivariate
Gaussian probability functions

$\displaystyle\mathcal{N}(\vec{\theta}|\vec{S})=\frac{1}{[(2\pi)^{3}|\det(C)|]^{1/2}}\exp\big{[}-(\vec{L}-\langle\vec{L}\rangle)^{T}C^{-1}(\vec{L}-\langle\vec{L}\rangle)\big{]},$
(2)

where $C$ is the covariance matrix

$\displaystyle
C=\begin{bmatrix}\epsilon^{2}_{x}+\sigma_{x}^{2}&\epsilon_{x}\epsilon_{y}\rho_{xy}&\epsilon_{x}\epsilon_{z}\rho_{xz}\\\
\epsilon_{y}\epsilon_{x}\rho_{xy}&\epsilon^{2}_{y}+\sigma_{y}^{2}&\epsilon_{y}\epsilon_{z}\rho_{yz}\\\
\epsilon_{z}\epsilon_{x}\rho_{xz}&\epsilon_{z}\epsilon_{y}\rho_{yz}&\epsilon^{2}_{z}+\sigma_{z}^{2}\end{bmatrix}$
(3)

and $\rho_{ij}$ are correlation coefficients, while
$(\epsilon_{x},\epsilon_{y},\epsilon_{z})$ are uncertainties associated with
individual angular momentum values estimated by Monte-Carlo sampling
observational errors on heliocentric distance ($\epsilon_{D}$), line-of-sight
velocities ($\epsilon_{v}$), and proper motions
($\epsilon_{\ell},\epsilon_{b}$). The hyperparameters $\sigma^{2}_{i}$, with
$i=x,y,z$, account for the spread in angular momentum of stars in each
Gaussian mixture component beyond those introduced by statistical errors (e.g.
see Hobson et al., 2002). For a better exploration of the prior volume, it is
convenient to express the mean angular momentum in spherical coordinates
$\langle\vec{L}\rangle=(L_{x},L_{y},L_{z})=L(\cos\theta\cos\phi,\cos\theta\sin\phi,\sin\theta)$.
Each component is determined by 7 parameters $\vec{S_{i}}=\\{\log
L,\theta,\phi,\sigma_{x},\sigma_{y},\sigma_{z},\alpha\\}_{i}$, with
$i=1,...,N$, hence the total number of parameters in our fits is $7N-1$. Our
analysis uses $N=4$ components with 27 parameters.

We adopt flat priors on $\\{\alpha,\log L,\cos\theta,\phi\\}$, and Jeffreys
priors for the hyperparameters $\\{\sigma_{x},\sigma_{y}$, $\sigma_{z}\\}$,
with ranges that include reasonable values. Mixture models are fitted with the
code MultiNest (Feroz & Hobson, 2008; Feroz et al., 2009), which uses a
nested-sampling technique (Skilling, 2004) to calculate posterior
distributions and the evidence of the model.

From the posterior distributions, we estimate the probability that an object
(star, GC or dwarf galaxy) with observed quantities $\vec{\theta}$ belongs to
the Sgr stream as

$\displaystyle p_{\rm mem}(\vec{\theta})=\frac{\alpha_{\rm
sgr}\,\mathcal{N}_{\rm
sgr}(\vec{\theta}|\vec{S}_{p})}{\sum_{i=1}^{N}\,\alpha_{i}\,\mathcal{N}_{i}(\vec{\theta}|\vec{S}_{p})},$
(4)

where $\vec{S}_{p}$ is an array that contains the median of the posterior
distributions.

Table 1: Posteriors derived from a Gaussian mixture model with $N=4$ components. Only halo stars within $|B_{\rm sgr}|<20^{\circ}$ of the Sgr orbital plane and at Galactocentric distances $r>20\,{\rm kpc}$ are included in the fit. Angular momentum $L$ and the hyperparameters $\sigma_{i}$ are given in units of $\,{\rm kpc}\,{\rm km\,s^{-1}}$.  Param. | #1 (Sgr) | #2 (GSE) | #3 (smooth) | #4 (Orphan)
---|---|---|---|---
$\alpha$ | $0.322^{+0.002}_{-0.002}$ | $0.239^{+0.002}_{-0.002}$ | $0.438^{+0.004}_{-0.005}$ | $0.0077^{+0.0003}_{-0.0003}$
$\langle L_{x}\rangle$ | $+426^{+7}_{-6}$ | $-0.7^{+0.4}_{-0.4}$ | $+144^{18}_{-19}$ | $-3905^{+20}_{-17}$
$\langle L_{y}\rangle$ | $-4950^{+11}_{-11}$ | $+11^{+0.3}_{-0.4}$ | $-278^{+10}_{-9}$ | $-2322^{+11}_{-11}$
$\langle L_{z}\rangle$ | $-1436^{+7}_{-6}$ | $+3^{+0.2}_{-0.2}$ | $+122^{16}_{-16}$ | $-4664^{+19}_{-15}$
$\sigma_{x}$ | $655^{+5}_{-6}$ | $10.3^{+0.08}_{-0.08}$ | $1755^{+8}_{-8}$ | $13^{+1}_{-4}$
$\sigma_{y}$ | $1255^{+13}_{-14}$ | $11.3^{+0.3}_{-0.4}$ | $1926^{+10}_{-12}$ | $12^{+0.6}_{-0.9}$
$\sigma_{z}$ | $659^{+6}_{-5}$ | $10.2^{+0.06}_{-0.06}$ | $1733^{+16}_{-16}$ | $13^{+0.9}_{-1.2}$

##  4 Results

###  4.1 Halo stars in angular momentum space

Table 1 provides the median and 1-sigma uncertainties of the parameters of a
Gaussian mixture model with $N=4$ components. The posterior distributions on
individual parameters are well behaved, with little covariance between them.
Comparison against the angular momentum of the Sagittarius dwarf,
$\vec{L}_{\rm sgr}=\vec{r}_{\rm sgr}\times\vec{v}_{\rm
sgr}=(17.9,2.6,-6.6)\,{\rm kpc}\times(239.5,-29.6,213.5)\,{\rm
km\,s^{-1}}=(+360,-5402,-1152)\,{\rm kpc}\,{\rm km\,s^{-1}}$ (VBE21),
indicates that the first Gaussian component, with $\alpha\simeq 32\%$ of our
sample, corresponds to tidal debris from the Sgr dwarf galaxy. Notice,
however, that the median angular momentum of the stream is slightly different
from that of the remnant core. This is to be expected given that tidal
stripped stars escape from the Lagrange points L1/L2 with lower/higher angular
momentum than the disrupting progenitor (e.g. Peñarrubia et al., 2006). The
hyperparameters, $\sigma_{i}$, also contain important information. In
particular, we find that the spread in the Y-axis, which is roughly
perpendicular to the orbital plane of the dwarf, is approximately twice as
large as in the other two directions. We come back to this result in §5.

The remainder components of the Gaussian mixture model also show remarkable
properties. For example, Component #2, with $\sim 24\%$ of the fitted sample,
corresponds to stars with little or no angular momentum, thus moving on nearly
radial orbits. These properties are similar to those of the Gaia-Sausage-
Enceladus (GSE) substructure (Belokurov et al., 2018; Helmi et al., 2018). The
properties of this system are worth exploring in a separate contribution.
Component #3 contains the largest fraction of stars with $\simeq 44\%$ of the
sample. It also shows the broadest angular momentum distribution, suggesting
an association with the extended ‘smooth’ MW halo. Component #4 is the
faintest object in the mixture model with only $\alpha\simeq 0.8\%$ of the
sample (roughly 40 stars). The mean angular momentum and the sky coordinates
of these stars indicate membership to the Orphan stream (Belokurov et al.,
2006; Fardal et al., 2019a). The detection of such a faint substructure using
Gaussian mixture models opens up interesting follow-up applications (see also
§5). We note that adding an extra component ($N=5$) does not lead to the
identification of additional substructures.

Fig. 1 shows projections of the angular momentum vector of 5389 stars in the
Sagittarius stream plane ($|B_{\rm sgr}|<20^{\circ}$) colour-coded according
to their probability of belonging to the Sagittarius stream, Equation (4). For
reference, purple crosses on the bottom-left corner mark the average size of
individual error bars, $\simeq 600\,{\rm kpc}\,{\rm km\,s^{-1}}$. This Figure
reveals a number of interesting features. First, we find 171 BHBs and 754 K
Giants with a high probability, $p_{\rm mem}>0.9$, of being members of the Sgr
stream. As expected, these stars clump around the orbital angular momentum of
the Sgr dwarf, which is marked with black stars for reference. However, a
considerable amount of scatter is visible in the three panels, specially in
the $L_{y}$ component, which is roughly perpendicular to the orbital plane of
the Sgr dwarf. This is consistent with the hyperparameter $\sigma_{y}$ being
approximately twice as large as $\sigma_{x}$ and $\sigma_{z}$ (see Table 1).
Second, the distribution of non-members ($p_{\rm mem}<0.5$) is clearly non-
Gaussian. As a result, our mixture models need at least $N=4$ components in
order to find a good match. The angular momentum $L=|\vec{L}|$ distribution of
the sample is shown in the right panel with a black line. To illustrate the
location of the Gaussian mixture components we weigh each star by their
membership probability and re-compute the histrograms. As expected, we find
that stars with high probability of belonging to the GSE substructure have
little angular momentum. Stars labelled as members of the ‘smooth’ halo have a
broad angular momentum distribution, whereas Sgr and Orphan stars clump in a
relatively narrow region.

The angular momentum of the Sgr stream has been recently studied by Johnson et
al. (2020), who propose a simple criterion for membership. Namely, stars with
angular momentum coordinates $L_{z}<-(2.5+L_{y})/0.3$, with $L$ measured in
units of $10^{3}\,{\rm kpc}\,{\rm km\,s^{-1}}$ (dotted line in Fig. 1), are
labelled members. Application of this cut to our dataset yields 1500 stars.
All the 925 stream members detected in our analysis at $p_{\rm mem}>0.9$
satisfy this condition. We also find 317 likely interlopers ($p_{\rm
mem}<0.5$). Hence, these results suggest that Johnson et al. (2020) sample is
complete but has a low purity. Given the results from our sample, where
$317/1500\approx 21\%$ of the stars selected by the Johnson et al. (2020)
criterion are low-likelihood members, a similar fraction of the Johnson et al.
(2020) stream sample may be MW interlopers. However, we caution that the
Johnson et al. (2020) catalogue has different sky coverage and phase-space
densities, which complicates a direct comparison.

PP21 studied the angular momentum distribution of stars in the outskirts of
the MW with the opposite goal in mind, namely to remove Sgr stream members
from a catalogue of ‘smooth’ halo stars. To this aim, they exclude 1062 stars
with angular momenta similar to that of the Sgr dwarf, $|\vec{L}-\vec{L}_{\rm
sgr}|<3000\,{\rm kpc}\,{\rm km\,s^{-1}}$ (dotted-line black/white circles in
Fig. 1). One can see by eye that this cut is imperfect, as some Sgr members
are located outside the circled regions. Indeed, we find that 793 Sgr stream
members ($p_{\rm mem}>0.9$) and 72 MW interlopers ($p_{\rm mem}<0.5$) are
located within this volume. This implies that PP21 cut successfully removes a
large fraction ($793/925\simeq 86\%$) of stream members from their halo
sample, while the misidentification of smooth halo stars as stream members is
relatively low, $72/1062\simeq 7\%$.

Figure 2: Distribution and kinematics of Sgr stream members at a $p_{\rm
mem}>90\%$ confidence level (see text). Red and blue dots denote K Giant and
BHB stars. Dark/light grey dots show leading/trailing tail particles from
VBE21 $N$-body model. For reference, the disc plane is shown with a solid
black line, with dotted lines showing Galactic latitudes $b=\pm 20^{\circ}$.
The North Galactic Pole is marked with a black cross. Notice the presence of
(previously unknown) second wraps of the trailing/leading tails in the
Northern (R.A.$>100^{\circ}$)/Southern (R.A.$<100^{\circ}$) hemispheres,
respectively. We detect 7 Globular Clusters associated with the Sgr stream
(cyan squares). Green dots show members of the Orphan stream (component #4 in
Table 1).

###  4.2 The Sagittarius & Orphan streams

In Fig. 2 we plot the observational coordinates of Sgr stream stars with
membership probabilities $p_{\rm mem}>0.9$. Blue and red stars denote BHB and
K Giant stars, respectively. This plot reveals that while the Sgr stream can
be described by a simple multivariate normal function in angular momentum (see
Fig. 1), the distribution in phase space is extraordinarily more complex. To
guide the interpretation of the detections, we over-plot with dark/light grey
particles the leading/trailing tails of the Sgr stream model recently
published by VBE21, which accounts for both the gravitational attraction of
the LMC as well as the displacement of the MW disc from the Galactic
barycentre in response to the LMC infall (see detailed discussion in PP21, ).

Upper panel of Fig. 2 shows the sky projections of the stream members. The sky
coverage is incomplete, with large gaps at $R.A.\gtrsim 200^{\circ}$ and
$60^{\circ}\lesssim R.A.\lesssim 100^{\circ}$ owing to the limited SEGUE
footprint. Notice also that all BHBs are located within $r\lesssim 40\,{\rm
kpc}$ from the MW centre. This is an observational bias rather than a physical
effect, as BHBs are fainter than K Giants, and thus more difficult to detect
at large distances. They also have on average larger associated errors in
proper motions and radial velocities, thus complicating a statistical
association with the stream.

Comparison against the $N$-body model shows that most stars have distances,
line-of-sight velocities and proper motions consistent with those of the most
recent wraps of the leading ($120^{\circ}\lesssim R.A.\lesssim 220^{\circ}$)
and trailing ($-10^{\circ}\lesssim R.A.\lesssim 60^{\circ}$) tails,
respectively. In addition, we detect the continuation of the trailing tail in
the Northern Hemisphere ($R.A.\sim 170^{\circ}$ and $v_{\rm los}\sim+200\,{\rm
km\,s^{-1}}$) as well as the leading tail in the Southern Hemisphere
($R.A.\sim 30^{\circ}$ and $v_{\rm los}\sim+100\,{\rm km\,s^{-1}}$), thus
increasing the stream coverage by a factor $\sim 2$ with respect to current 5D
maps (e.g. Ramos et al., 2020). Given that the $N$-body model has not been
fitted to this particular data set the overall agreement with observations is
remarkable. Yet, a few discrepancies are noticeable. E.g. the continuation of
the leading tail at $R.A.\lesssim 60^{\circ}$ shows systematically lower
$v_{\rm los}$ values than predicted by the model. We will return to this issue
in §5.

Interestingly, material from the old wraps may have been detected by Johnson
et al. (2020) after applying a cut in $L_{y}$–$L_{z}$ (see dotted line in Fig.
1). These authors found metal-poor stars associated with the Sgr stream that
are both off-set and with a diffuse distribution of line-of-sight velocities
compared to the metal-rich component, which led them to speculate the possible
existence of a stripped stellar halo of the Sagittarius dwarf. However, §3
shows that the proposed cut in $L_{y}$–$L_{z}$ returns a non-negligible number
of MW contaminants ($\approx 21\%$ in our sample), suggesting that the broad
velocity distribution may be partially explained the presence of MW
interlopers. In contrast, selecting members from the Gaussian mixture of Fig.
2 shows that the older wraps of the trailing/leading tails are kinematically
cold, in agreement with VBE21 models.

The mixture decomposition presented in §3 can be used to draw statistical
associations between each of the Gaussian components and any object in our
Galaxy with available phase-space information. For example, of the 170 known
Globular Clusters in the Milky Way, we find 7 with high probability ($p_{\rm
mem}>0.9$) of being members of the Sgr system: M54, Arp 2, Pal 12, Whiting 1,
NGC 2419, Terzan 7 and 8 (cyan boxes in Fig. 2). The rest have membership
probabilities $p_{\rm mem}\lesssim 0.01$, suggesting that there are no
additional known GCs associated with the Sgr dwarf, in agreement with recent
analyses (e.g. Arakelyan et al., 2021; Johnson et al., 2020). Our results do
not support recent claims that NGC 5634 and 4147 trace ancient wraps of the
stream (Bellazzini et al., 2020). No known satellite galaxy appears to be
associated with the Sgr dwarf.

In addition, we also detect several members of the Orphan stream (component #4
in Table 1), shown in Fig. 2 with green dots. Comparison against the
positions, velocities and proper motions found by Fardal et al. (2019b) shows
excellent agreement, indicating that the Gaussian mixture decomposition of §3
can also uncover faint substructures in the stellar halo even when these
contribute to $\lesssim 1\%$ of the fitted sample.

##  5 Discussion & Summary

We show that modelling the distribution of stars in angular momentum space
with Gaussian mixture techniques provides a powerful method to detect accreted
substructures that are partially, or fully mixed in phase-space, without
making assumptions on the form, shape or time-evolution of the MW potential.
Application to a catalogue of 5389 stars in the plane of the Sgr stream with
available 6D phase-space coordinates reveals the presence of at least four
independent dynamical components. We associate the broadest and most numerous
with the ‘smooth’ stellar halo, and the faintest (with only $0.8\%$ of the
sample) with the Orphan stream. Interestingly, we also identify a third
component with little or no angular momentum that likely corresponds to the
GSE substructure.

Our statistical technique detects two older wraps of the Sagittarius stream
that correspond to the continuation of the leading/trailing tails in the
Southern/Northern hemispheres, showing that the Sgr stream circles the Galaxy
at least twice. The full extent of the tidal tails is shown in Fig. 3 in a
reference frame aligned with the orbital plane of the stream. Dots show the
orbital poles of individual stream members relative to that of the Sagittarius
dwarf galaxy as a function of the angular separation from the remnant core
(with notation $\Lambda>0$ for the leading, and $\Lambda<0$ trailing tails).
The discovery of older wraps double the known extent of the Sgr stream. While
published 5D data cover parts of the stream within
$-150^{\circ}\lesssim\Lambda\lesssim 200^{\circ}$ from the Sgr dwarf, 6D
detections increase this range out to $-300^{\circ}\lesssim\Lambda\lesssim
500^{\circ}$.

For comparison, the orbital poles of the $N$-body model are shown with grey
dots. To estimate the time at which different parts of the stream became
tidally unbound from the Sgr dwarf, we find the $N$-body particle closest to
each individual star and colour-code stream members according to their
stripping time. We distinguish between material lost at $t_{\rm
strip}<1.5\,{\rm Gyr}$ (cyan) and $t_{\rm strip}>1.5\,{\rm Gyr}$ (orange),
which approximately correspond to particles unbound during the last and
penultimate pericentric passages of the progenitor dwarf, respectively (for
details, see VBE21, ). As expected, the continuation of the leading & trailing
tails correspond to material stripped early. Interestingly, the orbital poles
of the old wraps remain roughly aligned with those of the Sgr dwarf. This
appears somewhat in tension with the $N$-body model, which predicts a strong
bending of the leading tail at $\Lambda\gtrsim+400^{\circ}$. The mismatch is
mainly due to the systematically higher line-of-sight velocities predicted for
the continuation of the leading tail at $R.A.\lesssim 60^{\circ}$ (see Fig.2).
Notice also that the orbital poles are considerably more scattered in the
longitudinal direction (which is roughly perpendicular to the orbital plane)
than predicted by the model. This may be connected with the notoriously large
spread of the component $L_{y}$ highlighted in Fig. 1 and quantified in Table
1 via the hyperparameter $\sigma_{y}$. The large scatter of the orbital poles
and $L_{y}$ values may be related to the bifurcation of the leading & trailing
tails (Yang et al., 2019). The detection of old stream wraps provide new
constraints to investiate the origin of the bifurcation in more detail.

The Gaussian mixture decomposition identifies $\sim 40$ members of the Orphan
stream, which contribute to $\approx 0.8\%$ of the fitted sample. It thus
appears that our statistical analysis is able to detect narrow clumps in
angular momentum with a small number of members –which are typically
associated with disrupting globular clusters or faint satellite galaxies– as
well as fully phase-mixed, hot substructures lurking in the stellar halo.
Follow-up contributions will explore whether adding metallicity and/or
chemical elements as independent dimensions in the Bayesian likelihood (eq. 1)
helps to uncover clustering in angular momentum at even finer detail.

Figure 3: Galactocentric coordinates of the orbital poles of individual stream
members ($p_{\rm mem}>0.9)$ relative to that of the Sgr dwarf as a function of
the angle along the stream, $\Lambda$. Stars are colour-coded according to the
time when they were tidally stripped from the progenitor (see text). The
location of the Sgr dwarf is marked with a white star.

## Acknowledgements

We thank Teresa Antoja, Pau Ramos, Eugene Vasiliev, Vasily Belokurov, Denis
Erkal, Benjamin Johnson and Charlie Conroy for helpful comments. MSP
acknowledges funding from a UK Science and Technology Facilities Council
(STFC) Consolidated Grant.

## Data Availability

The Sgr stream sample may be found on github: https://github.com/michael-
petersen/SgrL.

## References

  * Arakelyan et al. (2021) Arakelyan N. R., Pilipenko S. V., Sharina M. E., 2021, Astrophysical Bulletin, 75, 394
  * Baumgardt & Vasiliev (2021) Baumgardt H., Vasiliev E., 2021, arXiv e-prints, p. arXiv:2105.09526
  * Bellazzini et al. (2020) Bellazzini M., Ibata R., Malhan K., Martin N., Famaey B., Thomas G., 2020, A&A, 636, A107
  * Belokurov et al. (2006) Belokurov V., et al., 2006, ApJ, 642, L137
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Bennett & Bovy (2019) Bennett M., Bovy J., 2019, MNRAS, 482, 1417
  * Dierickx & Loeb (2017) Dierickx M. I. P., Loeb A., 2017, ApJ, 836, 92
  * Drimmel & Poggio (2018) Drimmel R., Poggio E., 2018, Research Notes of the American Astronomical Society, 2, 210
  * Erkal et al. (2020) Erkal D., et al., 2020, arXiv e-prints, p. arXiv:2010.13789
  * Fardal et al. (2019a) Fardal M. A., van der Marel R. P., Law D. R., Sohn S. T., Sesar B., Hernitschek N., Rix H.-W., 2019a, MNRAS, 483, 4724
  * Fardal et al. (2019b) Fardal M. A., van der Marel R. P., Sohn S. T., del Pino Molina A., 2019b, MNRAS, 486, 936
  * Feroz & Hobson (2008) Feroz F., Hobson M. P., 2008, MNRAS, 384, 449
  * Feroz et al. (2009) Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601
  * Frinchaboy et al. (2012) Frinchaboy P. M., Majewski S. R., Muñoz R. R., Law D. R., Łokas E. L., Kunkel W. E., Patterson R. J., Johnston K. V., 2012, ApJ, 756, 74
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Gómez et al. (2015) Gómez F. A., Besla G., Carpintero D. D., Villalobos Á., O’Shea B. W., Bell E. F., 2015, ApJ, 802, 128
  * Gravity Collaboration et al. (2019) Gravity Collaboration et al., 2019, A&A, 625, L10
  * Helmi (2004) Helmi A., 2004, ApJ, 610, L97
  * Helmi et al. (2018) Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85
  * Hobson et al. (2002) Hobson M. P., Bridle S. L., Lahav O., 2002, MNRAS, 335, 377
  * Ibata et al. (1994) Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature, 370, 194
  * Ivezić et al. (2000) Ivezić Ž., et al., 2000, AJ, 120, 963
  * Johnson et al. (2020) Johnson B. D., et al., 2020, ApJ, 900, 103
  * Johnston et al. (2005) Johnston K. V., Law D. R., Majewski S. R., 2005, ApJ, 619, 800
  * Koposov et al. (2012) Koposov S. E., et al., 2012, ApJ, 750, 80
  * Kuhn & Feigelson (2017) Kuhn M. A., Feigelson E. D., 2017, arXiv e-prints, p. arXiv:1711.11101
  * Law & Majewski (2010) Law D. R., Majewski S. R., 2010, ApJ, 714, 229
  * Li et al. (2019) Li J., et al., 2019, ApJ, 874, 138
  * Majewski et al. (2003) Majewski S. R., Skrutskie M. F., Weinberg M. D., Ostheimer J. C., 2003, ApJ, 599, 1082
  * Martínez-Delgado et al. (2002) Martínez-Delgado D., Zinn R., Carrera R., Gallart C., 2002, ApJ, 573, L19
  * McConnachie & Venn (2020) McConnachie A. W., Venn K. A., 2020, Research Notes of the American Astronomical Society, 4, 229
  * Newberg et al. (2002) Newberg H. J., et al., 2002, ApJ, 569, 245
  * Peñarrubia et al. (2006) Peñarrubia J., Benson A. J., Martínez-Delgado D., Rix H. W., 2006, ApJ, 645, 240
  * Peñarrubia et al. (2010) Peñarrubia J., Belokurov V., Evans N. W., Martínez-Delgado D., Gilmore G., Irwin M., Niederste-Ostholt M., Zucker D. B., 2010, MNRAS, 408, L26
  * Peñarrubia et al. (2011) Peñarrubia J., et al., 2011, ApJ, 727, L2
  * Peñarrubia et al. (2016) Peñarrubia J., Gómez F. A., Besla G., Erkal D., Ma Y.-Z., 2016, MNRAS, 456, L54
  * Petersen & Peñarrubia (2020) Petersen M. S., Peñarrubia J., 2020, MNRAS, 494, L11
  * Petersen & Peñarrubia (2021) Petersen M. S., Peñarrubia J., 2021, Nature Astronomy, 5, 251
  * Ramos et al. (2020) Ramos P., Mateu C., Antoja T., Helmi A., Castro-Ginard A., Balbinot E., Carrasco J. M., 2020, A&A, 638, A104
  * Skilling (2004) Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, American Institute of Physics Conference Series Vol. 735, American Institute of Physics Conference Series. pp 395–405, doi:10.1063/1.1835238
  * Vasiliev & Baumgardt (2021) Vasiliev E., Baumgardt H., 2021, arXiv e-prints, p. arXiv:2102.09568
  * Vasiliev & Belokurov (2020) Vasiliev E., Belokurov V., 2020, MNRAS, 497, 4162
  * Vasiliev et al. (2021) Vasiliev E., Belokurov V., Erkal D., 2021, MNRAS, 501, 2279
  * Vera-Ciro & Helmi (2013) Vera-Ciro C., Helmi A., 2013, ApJ, 773, L4
  * Vivas et al. (2001) Vivas A. K., et al., 2001, ApJ, 554, L33
  * Xue et al. (2008) Xue X. X., et al., 2008, ApJ, 684, 1143
  * Xue et al. (2011) Xue X.-X., et al., 2011, ApJ, 738, 79
  * Yang et al. (2019) Yang C., et al., 2019, ApJ, 886, 154
  * Yanny et al. (2009) Yanny B., et al., 2009, AJ, 137, 4377
  * del Pino et al. (2021) del Pino A., Fardal M. A., van der Marel R. P., Łokas E. L., Mateu C., Sohn S. T., 2021, ApJ, 908, 244