# Simulating Pedestrian Avoidance: The Humans vs Zombies Scenario

Juan P. Oriana<EMAIL_ADDRESS>Instituto Tecnológico de Buenos Aires
(ITBA), Lavardén 315, 1437 C. A. de Buenos Aires, Argentina  German A.
Patterson<EMAIL_ADDRESS>Daniel R. Parisi<EMAIL_ADDRESS>Instituto Tecnológico de Buenos Aires (ITBA), CONICET, Lavardén 315, 1437 C.
A. de Buenos Aires, Argentina

###### Abstract

This study introduces a unique active matter system as an application of the
pedestrian collision avoidance paradigm, that proposes dynamically adjusting
the desired velocity. We present a fictitious human-zombie scenario set within
a closed geometry, combining prey-predator behavior with a one-way contagion
process that can transform prey into predators. The system demonstrates varied
responses: in cases where agents have the same maximum speeds, a single zombie
always catches a human, whereas two zombies never catch a single human. As the
number of human agents increases, observables, such as the final fraction of
zombie agents and total conversion times, exhibit a significant change in the
system’s behavior at intermediate density values. Most notably, there is
evidence of a first-order phase transition when the mean population speed is
analyzed as an order parameter.

pedestrian dynamics; collision avoidance; pedestrian navigation; predator-
pray; active matter; microscopic simulation

††preprint: APS/123-QED

##  I Introduction

A fundamental aspect in pedestrian dynamics modeling is describing how
pedestrians navigate and avoid collisions in various environments. Initially,
force-based models like the Social Force Model [1, 2] and the Predictive
Collision Avoidance Model [3] simulate collision avoidance maneuvers by
applying elusive forces directly on the agents.

An alternative framework, as proposed by Wang et al. [4] and expanded in [5,
6, 7, 8], involves dynamically adjusting the agents’ desired velocity for
navigation and avoidance maneuvers. This method closely resembles the natural
walking behavior, wherein individuals adapt their speed and direction
according to their surroundings, especially in scenarios where physical
contact is absent. Moreover, this approach can be implemented in any
microscopic agent model that has a specified target velocity, irrespective of
whether the movement is defined by forces or rules.

As an alternative experimentation platform for this framework, we suggest
exploring a fictional and basic system that imitates a predator-prey scenario.
This model involves agents resembling zombies ($Z$) attacking humans ($H$),
drawing inspiration from prevalent depictions in popular culture, including
movies, series, comics, and novels.

Other fictional systems involving zombies have previously been explored as an
analogy for understanding the transmission of infectious diseases [9, 10, 11].
These studies often consider macroscopic population variables governed by
differential equations to model their dynamics. On the other hand, adopting a
microscopic approach, Libál et al. [12] investigate an $SCZR$ (Susceptible-
Cleric-Zombie-Recovered) model by employing a population of self-propelled
particles performing random motion.

In our proposed zombie-human system, particles exhibit movement directions
determined by a specific heuristic. Zombie agents possess a desired velocity
pointing toward the nearest human agent, while human agents aim to avoid
zombies and collisions with walls and other humans. Additionally, our system
is simplistic in terms of infection, following a $SI$ (Susceptible-Infected)
model. The agents have only two states, and the only possible transition is
unidirectional: $H\rightarrow Z$.

Beyond offering a distinctive and imaginative framework for the pedestrian
navigation paradigm mentioned earlier, our proposed setting intriguingly
amalgamates two concepts within a unified system: a) predator-prey dynamics
and b) infectious disease contagion. This scenario stands out for its
distinctive aspect: the contagion not only affects the prey but also
transforms it into a predator.

The paper is structured as follows: Sec. II introduces the microscopic model,
while Sec. III outlines the simulation scenario. Sec. IV presents the outcomes
of numerous simulations conducted to characterize the system. Finally, Sec. V
contains our concluding remarks.

##  II The Human-Zombie Microscopic Model

The model considers two types of agents: zombies and humans. Common to both,
the rule-based operational model used is the CPM [13], plus a new avoidance
layer that modifies the agent’s desired velocity to prevent collisions. For
human agents, the goal is to evade zombies, walls and other humans, while for
zombies, the objective is to target the nearest human.

The original CPM assigns each agent a radius that can vary within the range
from ${r_{min}}$ to $r_{max}$ following specific rules. These limiting values
are shared by all agents in the simulation, regardless of whether they are
zombies or humans.

The agent’s radius is directly tied to the agent’s speed $v_{i}$ as

$|\mathbf{v}_{i}|=v_{i}=v_{i}^{max}\left(\frac{r-r_{min}}{r_{max}-r_{min}}\right)^{\beta}\
,$ (1)

where $\beta$ is a constant. The agent’s maximum velocity $v^{max}$ is the
free velocity and it is distinctive for each agent category:

$v_{i}^{max}=\begin{cases}v_{z}^{max}&\text{if agent $i$ is a zombie}\ ,\\\
v_{h}^{max}&\text{if agent $i$ is a human}\ .\end{cases}$ (2)

The velocity of each agent is adjusted in each discrete temporal increments
$\Delta t$ by Eq. (1). Simultaneously, the radius and position of each agent
are updated according to:

$\displaystyle r_{i}(t+\Delta{t})$
$\displaystyle=\min{\left[r_{max},r_{i}(t)+\frac{r_{max}}{\tau/\Delta{t}}\right]}\
,$ (3) $\displaystyle\mathbf{x}_{i}(t+\Delta{t})$
$\displaystyle=\mathbf{x}_{i}(t)+\mathbf{v}_{i}\Delta{t}\ ,$ (4)

where $\tau$ is a characteristic time constant describing the duration it
takes for a particle to reach its maximum speed after starting from rest.

On the other hand, when a collision occurs, the state of the involved
particles is updated as follows. First, they contract their radii to the
minimal value, $r_{min}$. Then, at the next time step, the agents undergo
repulsion using an escape velocity, $\textbf{v}^{e}$, in a direction
determined by the center of two particles, each moving away from each other.
The magnitude of this escape velocity corresponds to the maximum velocity,
$v^{e}=v^{max}$.

Moreover, in our specific scenario, the interaction between humans and zombies
results in the immobilization of both agents for a period $t_{c}$,
representing the contagion process. Upon completion, the human agent
transforms into a zombie and emerges with a radius of $r=r_{min}$.

Additionally, to counteract artifacts arising from symmetric configurations,
we add a small angular noise to the direction of the desired velocity. Let
$\theta_{i}$ represent the angle formed by $\textbf{v}_{i}$ with respect to
the positive x-semiaxis, and let $\theta^{t}_{i}$ denote the angle of the
desired direction vector of agent $i$, thus,
$\theta_{i}=\theta_{i}^{t}+\Delta\theta$, with $\Delta\theta$ defined as a
uniform random value within the range of $[-\mu/2,\mu/2]$.

In what follows, we describe the heuristic that allows the avoidance layer to
dynamically adjust the direction of the desired velocity of agents:
$\mathbf{v}_{i}=v_{i}\mathbf{e}_{i}^{t}$. This direction is represented by the
unit vector
$\mathbf{e}^{t}_{i}=\frac{\mathbf{x}^{t}_{i}-\mathbf{x}_{i}}{{|\mathbf{x}^{t}_{i}-\mathbf{x}_{i}|}}$,
which is determined by a temporal target denoted as $\mathbf{x}_{i}^{t}$. This
target is computed differently depending on whether agent $i$ is a human or a
zombie.

If agent $i$ is a human, $\mathbf{{x}}_{i}^{t}$ is obtain from the sum of the
collision vectors $\mathbf{n}_{c}^{ij}$. Each $\mathbf{n}_{c}^{ij}$ quantifies
the extent to which a human $i$ should adjust its direction to move away from
an agent $j$. Therefore, this desired target is defined by Eq. 5.

$\mathbf{x}_{i}^{t}=\sum_{z\in Z^{i}}\mathbf{n}_{c}^{iz}+\sum_{h\in
H^{i}}\mathbf{n}_{c}^{ih}+\mathbf{n}_{c}^{iw}\ ,$ (5)

where $Z^{i}$ represents the set of the $N^{e}_{z}$ nearest zombies to agent
$i$, $H^{i}$ the set of the $N^{e}_{h}$ nearest humans to agent $i$.

The vectors $\mathbf{n}_{c}^{ij}$ are computed as:

$\mathbf{n}_{c}^{ij}=\mathbf{e}^{ij}Ae^{-\frac{d_{ij}}{B}}\ ,$ (6)

where $d_{ij}$ is the distance between agents,
$\mathbf{e}^{ij}=\frac{\mathbf{x}^{i}-\mathbf{x}^{j}}{|\mathbf{x}^{i}-\mathbf{x}^{j}|}$
is the direction from agent $i$ to agent $j$, and $A$ and $B$ are constants
denoting the intensity and a characteristic length, respectively. The latter
two parameters vary depending on whether agent $j$ is a human, a zombie, or
the wall itself.

The last term, $\mathbf{n}_{c}^{iw}$, accounts for the interaction of agent
$i$ with walls ($j=w$). In this specific case, we will consider
$\mathbf{x}^{w}$ as the nearest point on the boundary to agent $i$.

If agent $i$ is a zombie, its target is the nearest human agent $h$, defined
by $\mathbf{x}_{i}^{t}=\mathbf{e}^{ih}$.

All the previously defined vectors can be visualized in Fig. 1

$\mathbf{a}_{k}$$\mathbf{a}_{i}$$\mathbf{a}_{j}$$\mathbf{n}_{c}^{ij}$$\mathbf{n}_{c}^{ik}$$\mathbf{n}_{c}^{iw}$$\mathbf{x}_{i}^{t}$$\mathbf{x}_{k}^{t}$
Figure 1: Selection of the dynamic target $\textbf{x}_{i}^{t}$ for a human
agent $a_{i}$ depends on the positions of the wall, another human $a_{j}$
(blue) and a zombie agent $a_{k}$ (green). Also depicted, selection of dynamic
target $\textbf{x}_{k}^{t}$ for a zombie agent $a_{k}$, solely dependent on
human agent $a_{i}$.

##  III Simulations

The agents move within a circular arena with a radius $R_{a}=11$ m, as shown
in Fig. 2. The simulations start with a solitary zombie agent positioned at
the center of the domain, and $N_{h}(t=0)\equiv N_{h}$ human agents randomly
distributed, ensuring that their distance from both the center and the wall
exceeds 1 m (1 m $\leq|\mathbf{x}_{i}|\leq R_{a}-$ 1 m), as illustrated by the
dotted circles in Fig. 2.

$R_{a}$$\mathbf{x}_{i}$ Figure 2: Illustration of the simulation domain and
the initial configuration featuring a zombie at the center (green agent) and
$N_{h}=5$ human agents (blue agents).

Considering the model described in Sec. II, the fixed parameters for the
agents are the following: $\Delta t=0.05$ s; $r_{min}=0.15$ m; $r_{max}=0.35$
m; $v_{h}^{max}=4.0$ m/s; $\beta=0.9$; $\tau=0.5$ s; $\mu=0.052$ rad (3°);
$A_{h}=4$; $B_{h}=1$ m; $N^{e}_{h}=1$; $A_{z}=8$; $B_{z}=4$ m; $N^{e}_{z}=2$;
$A_{w}=8$; $B_{w}=1$ m; $t_{c}=3$ s.

The initial number of humans ($N_{h}$) will vary between 10 and 100, in steps
of 5, and is considered the control parameter. It is important to note that
the total number of agents is $N=N_{h}+1$ ($N_{h}$ humans plus one zombie),
and this number remains constant in each simulation as human agents turn into
zombies. Moreover, the free speed of the zombie, $v_{z}^{max}$, is swept from
3.8 m/s to 4.2 m/s in steps of 0.1 m/s, being the reference case
$v_{z}^{max}=4$ m/s . Overall, there are a total of 95 system configurations
to explore. The statistical analysis of each was conducted on 2500
realizations of random initial conditions. Each realization continued until
either all agents had become zombies or the simulated time reached
$t_{max}=2000$ s. We define the final simulation time as $T$, and the number
of simulation steps as $S$. Hence, $S=\frac{T}{\Delta{t}}$. In the case that
$T=t_{max}$, then $S=40,000$.

We will now define the observables under investigation. First, the fraction of
zombies at any given time $t$ in a realization $k$ is given by:

$\phi^{k}_{z}(t)=\frac{N^{k}_{z}(t)}{N^{k}_{h}+1}=\frac{N^{k}_{z}(t)}{N}$ (7)

and compute $\left\langle\phi_{z}(t)\right\rangle$, which represents the
average of $\phi_{z}^{k}$ across the realizations. We focus our study on the
endpoint of the average final zombie fraction (hereafter referred to simply as
the final zombie fraction), denoted as
$\left\langle\phi_{z}^{final}\right\rangle$, which is defined when
$\left\langle\phi_{z}(t)\right\rangle$ reaches a stationary value.

The second observable we calculate is the time required to achieve total
conversion, denoted as the ’total conversion time’, $T_{c}$, and we average it
across all realizations as the ’average total conversion time’, $\left\langle
T_{c}\right\rangle$. Any realization that did not achieve total conversion
within $t_{max}$ was excluded from this particular analysis. If fewer than 10%
of realizations achieve total conversion, that specific configuration will not
be considered for analysis.

The last observable we define is the mean velocity, used as a type of order
parameter. This is based on the understanding that zombie-human interactions
typically occur at zero velocity, consequently reducing the mean velocity of
the agents. This reduction is even more significant than the usual decrease
caused by density in pedestrian dynamics and the CPM model [13].

Considering a specific realization $k$ at time $t$, the system’s velocity,
accounting for both agent types, is expressed as follows:

$v^{k}(t)=\frac{1}{N}\sum_{i=1}^{N}|\mathbf{v}_{i}^{k}(t)|\ ,$ (8)

where $\mathbf{v}_{i}$ is the vector velocity of particle $i$. Then, we
compute the average across the realizations as:

$\bar{v}=\bar{v}(t)=\frac{1}{K}\sum_{k=1}^{K}v^{k}(t)\ ,$ (9)

with $K=2500$, the number of realizations per configuration. Lastly, in the
case of stationary states, the time average of mean population velocity is
computed as $\left\langle\bar{v}\right\rangle_{t}$.

##  IV Results

###  IV.1 Dynamics involving 2 and 3 agents

The behaviors of the system when $N_{h}=1$ and $N_{h}=2$ are interesting and
notably different from one another. In the simplest scenario, the system
consists of a single human agent ($N_{h}=1$) in addition to the initial zombie
($N_{z}(t=0)=1$) within the simulation domain. We will show that, under
conditions of equal maximum speed ($v_{z}^{max}=v_{h}^{max}=v^{max}$), the
zombie will always catch the human agent.

Figure 3: Dynamics with 2 agents. A: Human agent evading a zombie while moving
along the curved wall. B: Distance of agents from the center of the enclosure.

In the long-term pursuit between a human agent and a zombie, the human’s
trajectory becomes circular, running parallel to the enclosure’s wall. This
behavior is illustrated in Fig. 3A, showing the trajectories of both agents.
Fig. 3B depicts the distance of both agents from the center of the enclosure
over time. The zombie agent pursues the prey with a smaller turning radius,
resulting in a higher angular velocity. Consequently, the zombie always
catches up with its prey. See Supplemental Material [14] for a video showing
this dynamics.

This outcome indicates that in our simulations, when the speed of zombies is
similar to that of humans, the initial zombie will inevitably succeed in
converting at least one human. This claim is strengthened by the fact that
human agents must evade not just walls, but also other agents, which, on top
of that, can reduce their speed.

An important consideration is the potential for a human to escape if a model
incorporating inertia and the capability for zig-zag evasive maneuvers is
employed. Under these conditions, the aforementioned conclusion regarding an
isolated human escaping from a zombie would no longer hold.

Let us now examine the immediately more intricate case ($N_{h}=2$), showcasing
behavior markedly distinct from the $N_{h}=1$ scenario. Under this condition,
the zombie chases the nearest human until it catches him, which will
definitely happen, as demonstrated above. After converting the human into a
zombie, the system’s state changes to $N_{z}(t)=2$ and $N_{h}(t)=1$. Given
these conditions, with both zombies having only one remaining human target,
they both move towards that individual. With the two zombies being in close
proximity following the conversion, they collide while moving towards the same
target. This collision decelerate the zombies, granting the human an
opportunity to accelerate away from both of them at maximum speed. As they
resume their approach towards the target, after a few time steps, they collide
once more, again diminishing their speed and allowing the human to distance
themselves again. This pattern repeats indefinitely, ensuring that two
zombies, with maximum desired speeds equal to that of the human, can never
catch him, within the context of this operational model. An animation of this
dynamics can be seen on Supplemental Material [15].

Summing up, under the given conditions, when $N_{h}=1$ the conversion will
always be complete, i.e.: $\phi_{z}=1$. In contrast, this will never occur
when $N_{h}=2$, resulting in a stable zombie fraction of $\phi_{z}=2/3$.

###  IV.2 Dynamics involving more than 3 agents

Having previously analyzed the fundamental cases of one and two humans, we now
turn our attention to more densely populated initial conditions.

To illustrate the system’s dynamics, we present several representative
snapshots in Fig. 4 for different initial values of $N_{h}$, taken at specific
moments corresponding to distinct zombie fractions.

It is crucial to highlight that after the first conversion (when reaching
$N_{z}=2$), the two zombies may either target the same human agent or
different ones, depending on the positions of other human agents. In the
former scenario, a loop could arise where the remaining humans can escape from
the two zombies that continuously collide with each other (similar to the
$N_{h}=2$ scenario discussed in the previous Sec. IV.1). This scenario is
termed ’convergent pursuit,’ and a video example is available in the
Supplemental Material [16]. In the latter scenario, where the two zombies
split to pursue different human agents, they will succeed, mirroring the
$N_{h}=1$ situation described in Sec. IV.1. This is referred to as ’divergent
pursuit,’ and a video is available in the Supplemental Material [17].
Additionally, snapshots of both types of pursuit can be seen in Fig. 4B for
’divergent’ and Fig. 4E and H for ’convergent’.

A

B

C

D

E

F

G

H

I

Figure 4: Temporal evolution of three exemplary configurations. A: $N_{h}=20$
at $t=0$. B: $N_{h}=20$ at first conversion. C: $N_{h}=20$ at $\phi_{z}=0.5$.
D: $N_{h}=50$ at $t=0$. E: $N_{h}=50$ at first conversion. F: $N_{h}=100$ at
$\phi_{z}=0.5$. G: $N_{h}=100$ at $t=0$. H: $N_{h}=100$ at first conversion.
I: $N_{h}=100$ at $\phi_{z}=0.5$.

After establishing these definitions, simulations are categorized based on the
following criteria: if, post the initial conversion, both zombies remain
within a distance of less than 2 meters for 3 seconds, the scenario is
classified as convergent pursuit. Otherwise, it falls under divergent pursuit.

Figure 5 shows that the distribution of total simulations between each
category is fairly balanced. However, when considering only simulations that
result in total conversion, the only mechanism leading to total conversion for
$N_{h}\leq 25$ is the divergent pursuit type. As $N_{h}$ increases, a
transition is observed, reaching a balance after $N_{h}\geq 45$. In this
region of elevated densities, both types of initial pursuits exhibit
comparable probabilities of achieving total conversion due to the decelerating
effect caused by the increased density of human agents.

Figure 5: Fraction of divergent pursuit cases as a function of the number of
initial human agents in the system with a maximum zombie speed of
$v_{z}^{max}=v_{h}^{max}=4$ m/s.

####  IV.2.1 Final Zombie Fraction

In this section, we examine the dynamics of the population, focusing on the
average final zombie fraction, $\left\langle\phi_{z}^{final}\right\rangle$
[Eq. 7]. Since the number of agents remains constant in each simulation, the
fraction of human agents is complementary to the zombie fraction:
$\left\langle\phi_{z}^{final}\right\rangle=1-\left\langle\phi_{h}^{final}\right\rangle$.
Consequently, we focus our discussion solely on the zombie fraction.

First, we examine the time evolution of $\left\langle\phi_{z}(t)\right\rangle$
in Fig. 6A for the scenario where $v_{z}^{max}=v_{h}^{max}=4$ m/s. The mean
zombie fraction is observed to reach stationary values within the simulation
time, clustering around low and high values of $N_{h}$. A notable transition
occurs around intermediate values of $N_{h}\sim 40$, requiring more time to
reach the stationary regime.

A

B

Figure 6: Zombie Fraction. A: Time evolution of $\phi_{z}$ for $v_{z}^{max}=4$
m/s and the specified $N_{h}$ values. B: Final zombie fraction
$\phi_{z}^{final}$ as a function of the initial number of human agents,
categorized by the type of pursuit.

This transition becomes clearer in Fig .6A, where
$\left\langle\phi_{z}^{final}\right\rangle$ is plotted against the initial
number of human agents. The curve labeled ’All cases’ represents the
stationary data points from Fig. 6A. For values up to $N_{h}=30$, the final
zombie fraction is $\left\langle\phi_{z}^{final}\right\rangle\approx 0.5$.
However, beyond $N_{h}=45$, the system invariably results in total conversion.

The intermediate values of $\left\langle\phi_{z}^{final}\right\rangle$ at
lower $N_{h}$ arise from two extreme possible outcomes of the simulations: one
with a very high zombie fraction and the other with a very low one. In the
latter scenario, humans have ample space to outpace the zombies, who collide
with each other, impeding their pursuit (convergent pursuit). Figure 6B
illustrates that these two extremes correspond to divergent and convergent
pursuit scenarios. This confirms our earlier observations from the analysis of
Fig. 5: total conversion at lower $N_{h}$ values is solely observed in cases
of divergent pursuit.

On the other hand, when the initial human population exceeds $N_{h}\sim 45$,
the constrained space inhibits human escape without mutual interference. This
heightened density reduces their speed, invariably leading to zombies catching
them, regardless of pursuit type.

Now, we explore the same observable for varying maximum zombie speeds,
maintaining a constant human maximum speed of $v_{h}^{max}=4$ m/s. Figure 7A
displays five curves corresponding to $v_{z}^{max}\in[3.8,4.2]$ m/s.
Observations reveal that the curves qualitatively exhibit a similar pattern,
showing a transition from partial to total conversion at intermediate $N_{h}$
values. However, as $v_{z}^{max}$ increases, the final zombie fraction at
lower $N_{h}$ values decreases, and the $N_{h}$ value at which the system
saturates becomes higher. Both observations support the concept that higher-
speed zombies can capture all human agents at lower critical densities.
Additionally, the likelihood of converting the entire population becomes
higher at these lower densities.

AB Figure 7: Final Conversion for $v_{z}^{max}$ Lower and Higher than
$v_{h}^{max}$. A: Final zombie fraction as a function of $N_{h}$ for various
$v_{z}^{max}$ values. B: Ratio of realizations resulting in total conversion
to the total number of simulations conducted for each $N_{h}$ and
$v_{z}^{max}$ value.

In Fig.7B the fraction of simulations achieving total conversion is shown. For
the cases of $v_{z}^{max}\geq v_{h}^{max}=4$ m/s, the trend mirrors that of
the final zombie fraction, albeit with noticeable differences observed for
lower values of $v_{z}^{max}$. This occurs because a small portion of the
realizations fail to achieve total conversion. This is primarily caused by
multiple zombies colliding among themselves while pursuing a single human
agent, preventing them from catching it. Furthermore, in instances where total
conversion does not occur, the final fraction tends to be significantly high,
usually around $\phi_{z}\sim\frac{N_{h}}{N_{h}+1}$. This explains the
disparities between panels A and B of Fig. 7 when the maximum zombie speed is
lower than the human speed. The video of a specific realization for this
scenario is available in the Supplemental Material [18], with
$v_{z}^{max}=3.8$ m/s being less than $v_{h}^{max}$.

####  IV.2.2 Total Conversion Times

Now, we turn our attention to the study of the time needed to achieve total
conversion. As Fig. 7B clearly shows, total conversion is not always attained,
being more probable at higher values of $N_{h}$ and $v_{z}^{max}$. Hence, the
pertinent simulations within this context constitute a subset of those
discussed in the preceding section. With this in mind, Fig. 8A illustrates the
average total conversion time as a function of the initial number of human
agents and for different maximum zombie speeds. Intriguingly, all curves
exhibit a distinct maximum for all zombie speeds. In scenarios where zombie
speeds exceed human speed, the peak conversion time occurs near $N_{h}\sim
30$, while for zombie speeds lower than or equal to human speed, the peak is
around $N_{h}\sim 40$.

This suggests that achieving complete contagion occurs relatively quickly,
irrespective of whether the initial number of humans is small or large. With
fewer humans, they are less easily captured, yet the total conversion time
remains low due to the smaller count of agents requiring conversion.
Conversely, scenarios with a high number of humans result in easier capture
(as density reduces their speed). At intermediate human numbers, a balance
between these effects emerges: they are not easily captured, yet their
adequate numbers contribute to slowing down the process.

Additional insight is presented in Fig. 8B, where scenarios of convergent and
divergent pursuit are delineated for the case of zombie speed equal to human
speed ($v_{z}^{max}=v_{h}^{max}$). Convergent pursuit results in total
conversion only for $N_{h}\geq 30$, and, as expected, it represents a slow
process that results in prolonged conversion times for intermediate values of
$N_{h}$. In contrast, divergent pursuit serves as the sole mechanism for
achieving total conversion at lower $N_{h}$ values, occurring more rapidly.
The peak observed in total conversion times distinctly arises from averaging
the swift conversions in divergent pursuit with the slower, diminishing
conversion times linked to convergent pursuit.

AB Figure 8: Average Total Conversion Time. A: Across different configurations
of $v_{z}^{max}$. B: Categorized by the type of pursuit for the scenario where
$v_{z}^{max}=v_{h}^{max}=4$ m/s.

####  IV.2.3 Mean Population Velocities

Previous studies have suggested a shift in the system’s behavior at
intermediate $N_{h}$ values. To delve deeper into this phenomenon, we will now
investigate the system’s mean speed.

The determination of mean velocity becomes less straightforward when humans
and zombies possess different maximum speeds. In cases where
$v_{z}^{max}<v_{h}^{max}$, the mean speed at low $N_{h}$ takes an indefinitely
long time to stabilize. Conversely, in scenarios where
$v_{z}^{max}>v_{h}^{max}$, simulations with high $N_{h}$ values typically
result in total conversions before the mean speed can stabilize.

Stationary mean population speeds are evident solely when humans and zombies
share equal maximum speeds ($v_{z}^{max}=v_{h}^{max}$), as depicted in Fig.
9A. This figure also demonstrates that the stationary speed for low $N_{h}$ is
close to the maximum ($\bar{v}\approx v_{z}^{max}\approx v_{h}^{max}\approx 4$
m/s), while for high $N_{h}$, the mean speed falls below $2$ m/s. At
intermediate $N_{h}$ values, there is a noticeable transition in mean speed
between these two groups.

This transition becomes clearer when examining the temporal mean in the
stationary regime of mean speed (considered after $t=250$ s), as shown in Fig.
9B. The standard deviations of the mean speed peak at $N_{h}=50$, suggesting a
potential phase transition at the critical value $N_{hc}\approx 50$. In Fig.
9C, the fourth-order Binder cumulant
$U_{4}=1-\frac{\left\langle\bar{v}^{4}\right\rangle_{t}}{3\left\langle\bar{v}^{2}\right\rangle_{t}^{2}}$
is computed for the same data, exhibiting a distinct minimum at $N_{hc}$,
indicative of a first-order transition [19].

ABC Figure 9: Population velocities for the scenario
$v_{z}^{max}=v_{h}^{max}=4$ m/s with various $N_{h}$ values. A: Evolution of
population velocity averaged across all realizations. B: Mean population
velocities after $t=250$ s, accompanied by the corresponding standard
deviation. C: Binder parameter computed using the same data as in panel B.

##  V Conclusions

Employing an operational model of pedestrian dynamics [13], known for
accurately replicating the decrease in average speeds based on density
(fundamental diagram of pedestrian traffic), we enhanced it by integrating
evasion and pursuit mechanisms contingent upon the agent’s type. The prey
agents (humans) transition to the predator state (zombies) upon being captured
by the predator.

This novel system not only acts as a test scenario for implementing the
collision avoidance paradigm, where alterations are made solely to the agents’
desired velocity, but it also exhibits a non-trivial emergent dynamic.

In basic scenarios involving two or three agents, it was observed that within
a confined space, a predator can consistently catch its prey, even when they
possess identical maximum speeds. However, when two predators pursue the same
prey, their interference diminishes their speed, preventing them from ever
capturing the sole prey.

The mutual collision mechanism can be disrupted when multiple preys exist, and
each predator pursues a different prey. This scenario, termed ’divergent
pursuit’, becomes a crucial factor in achieving total conversion in
environments characterized by lower population density.

In scenarios involving numerous agents, we regard the global density ($N_{h}$,
the initial number of humans) as the controlling parameter. Both the final
zombie fraction and the average times for total conversion manifest a distinct
behavioral change at intermediate $N_{h}$ values. In the case of the final
zombie fraction, it abruptly transitions from $\phi_{z}^{final}=0.5$ to
$\phi_{z}^{final}=1$, and for the total conversion times, a peak is observed
at values of $N_{h}$ that depend on the maximum zombie speed.

Remarkably, when considering the average population velocity as an order
parameter in scenarios where the maximum speeds of humans and zombies are
equal, the results present evidence of a first-order phase transition
occurring at the critical value of $N_{hc}\approx 50$.

The proposed human-zombie scenario could be expanded into more complex
versions by including several probabilistic elements: the chance that a human
can eliminate a zombie, the likelihood of a human being immune, and the
possibility of a zombie being cured and reverting to a human state. These
additions would introduce new dynamics and complexities to the system.

##  VI Acknowledgments

This work was funded by project PICT 2019-00511 (Agencia Nacional de Promoción
Científica y Tecnológica, Argentina).

## References

  * Helbing and Molnár [1995] D. Helbing and P. Molnár, Physical Review E 51, 4282 (1995).
  * Helbing et al. [2000] D. Helbing, I. Farkas, and T. Vicsek, Nature 407, 487 (2000).
  * Karamouzas et al. [2009] I. Karamouzas, P. Heil, P. Van Beek, and M. H. Overmars, in _Motion in Games: Second International Workshop, MIG 2009, Zeist, The Netherlands, November 21-24, 2009. Proceedings 2_ (Springer, 2009), pp. 41–52.
  * Wang et al. [2015] Q.-L. Wang, Y. Chen, H.-R. Dong, M. Zhou, and B. Ning, Chinese Physics B 24, 038901 (2015).
  * Martin and Parisi [2020a] R. Martin and D. Parisi, Collective Dynamics 5, 324 (2020a).
  * Martin and Parisi [2020b] R. Martin and D. Parisi, Collective Dynamics 5 (2020b).
  * Martin and Parisi [2020c] R. F. Martin and D. R. Parisi, in _Traffic and Granular Flow 2019_ (Springer, 2020c), pp. 205–210.
  * Martin and Parisi [2023] R. F. Martin and D. R. Parisi, Physica A: Statistical Mechanics and its Applications p. 129414 (2023).
  * Munz et al. [2009] P. Munz, I. Hudea, J. Imad, and R. J. Smith, Infectious disease modelling research progress 4, 133 (2009).
  * Alemi et al. [2015] A. A. Alemi, M. Bierbaum, C. R. Myers, and J. P. Sethna, Physical Review E 92, 052801 (2015).
  * Bauer [2019] H. Bauer, Undergraduate thesis, University of Lynchburg (2019), URL https://digitalshowcase.lynchburg.edu/cgi/viewcontent.cgi?article=1122&context=utcp.
  * Libál et al. [2023] A. Libál, P. Forgács, Á. Néda, C. Reichhardt, N. Hengartner, and C. Reichhardt, Physical Review E 107, 024604 (2023).
  * Baglietto and Parisi [2011] G. Baglietto and D. R. Parisi, Physical Review E 83 (2011).
  * SMV [a] https://youtu.be/ZCp7hJBvf14.
  * SMV [b] https://youtu.be/84liICIE0Ek.
  * SMV [c] https://youtu.be/YDJkbw_ayiQ.
  * SMV [d] https://youtu.be/g8WlDL1h65A.
  * SMV [e] https://youtu.be/yIMXLRIe1mY.
  * Puzzo et al. [2022] M. L. R. Puzzo, E. S. Loscar, A. De Virgiliis, and T. S. Grigera, Journal of Physics: Condensed Matter 34, 314001 (2022).


# Close Binary Fractions in _accreted_ and _in-situ_ Halo Stars

Dolev Bashi,1 Vasily Belokurov,2 Simon Hodgkin2
1Astrophysics Group, Cavendish Laboratory, University of Cambridge, JJ Thomson
Avenue, Cambridge CB3 0HE, UK
2Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The study of binary stars in the Galactic halo provides crucial insights into
the dynamical history and formation processes of the Milky Way. In this work,
we aim to investigate the binary fraction in a sample of _accreted_ and _in-
situ_ halo stars, focusing on short-period binaries. Utilising data from Gaia
DR3, we analysed the radial velocity (RV) uncertainty $\sigma_{\mathrm{RV}}$
distribution of a sample of main-sequence stars. We used a novel Bayesian
framework to model the dependence in $\sigma_{\mathrm{RV}}$ of single and
binary systems allowing us to estimate binary fractions $F$ in a sample of
bright ($G_{\mathrm{RVS}}$< 12) Gaia sources. We selected the samples of _in-
situ_ and _accreted_ halo stars based on estimating the 6D phase space
information and affiliating the stars to the different samples on an action-
angle vs energy ($L_{\mathrm{z}}-E$) diagram. Our results indicate a higher,
though not significant, binary fraction in _accreted_ stars compared to the
_in-situ_ halo sample. We further explore binary fractions using cuts in $E$
and $L_{z}$, and find a higher binary fractions in both high-energy and
prograde orbits that might be explained by difference in metallicity. By
cross-matching our Gaia sample with APOGEE DR17 catalogue, we confirm the
results of previous studies on higher binary fractions in metal-poor stars and
find the fractions of _accreted_ and _in-situ_ halo stars consistent with this
trend. Our finding provides new insights into binary stars’ formation
processes and dynamical evolution in the primordial Milky Way Galaxy and its
_accreted_ dwarf Galaxies.

###### keywords:

binaries: close – Galaxy: halo – Galaxy: kinematics and dynamics – methods:
statistical – techniques: radial velocities

††pubyear: 2024††pagerange: Close Binary Fractions in _accreted_ and _in-situ_
Halo Stars–Close Binary Fractions in _accreted_ and _in-situ_ Halo Stars

##  1 Introduction

The Galactic halo, a vast and sparsely populated region enveloping the Milky
Way, presents a unique environment for the study of stellar populations
(Helmi, 2008, 2020; Deason & Belokurov, 2024). It comprises predominantly old
stars, remnants of both the primordial galactic formation era and the
subsequent accretion events. This dichotomy in the halo’s star formation
history is epitomised in the two distinct groups of stars: ’ _accreted_ ’ halo
stars, originating from smaller, consumed galaxies and the ’ _in-situ_ ’ halo
stars, formed within the nascent Milky Way (Deason & Belokurov, 2024).

The _accreted_ component of the halo in the solar neighbourhood, is
exemplified by the Gaia-Sausage-Enceladus (GS/E; Belokurov et al., 2018; Helmi
et al., 2018; Haywood et al., 2018), a prominent substructure discerned
through analyses of Gaia data (Gaia Collaboration et al., 2016, 2018)
indicating the remnant of the Milky Way’s most consequential recent merger, an
event estimated to have concluded around $8-10$ billion years ago. The GS/E
appears to dominate the halo around the Sun, but a wide variety of smaller
halo substructures have also been detected (Helmi et al., 1999; Myeong et al.,
2018a, b, 2019; Matsuno et al., 2019; Koppelman et al., 2019; Naidu et al.,
2020; O’Hare et al., 2020; Malhan, 2022; Dodd et al., 2023).

Conversely, the _in-situ_ halo population includes the ’splashed’
high-$\alpha$ disc (Bonaca et al., 2017; Gallart et al., 2019; Di Matteo et
al., 2019; Belokurov et al., 2020), and the Aurora stars (Belokurov &
Kravtsov, 2022; Conroy et al., 2022; Chandra et al., 2023; Zhang et al.,
2023), both formed in the Milky Way proper before the arrival of the GS/E
progenitor. _Splash_ stars, with their relative metal-rich composition ([Fe/H]
$\gtrsim-1$) and highly eccentric orbits, are thought to have been displaced
from the Galactic disc during the cataclysmic merger event that formed the
GS/E structure. Having been formed before the merger, _Splash_ stars are on
average as old as the stars in the GS/E (see e.g. Belokurov et al., 2020).
_Aurora_ is characterized by its metal-poor nature ([Fe/H] $\leq$ -1) and a
wide range of azimuthal velocities Belokurov & Kravtsov (2022); Chandra et al.
(2023); Zhang et al. (2023). Identified through chemical tagging, _Aurora_
represents the oldest in-situ stellar population of the Milky Way, formed
before the Galaxy had developed a coherently rotating disc (see discussion in
Semenov et al., 2024; Dillamore et al., 2024). While first spotted in the
vicinity of the Sun, _Aurora_ likely follows a steep radial density, with most
of its mass near the Galactic centre (see Belokurov & Kravtsov, 2023; Rix et
al., 2022).

Alongside the _accreted_ stars, these components of the _in-situ_ halo,
_Splash_ and _Aurora_ , underscore the Milky Way’s dynamic history, from its
earliest star formation episodes to the significant reshaping triggered by
galactic mergers. This layered approach facilitated by the Gaia space
observatory data, provides a nuanced understanding of how different stellar
populations within the Galactic halo offer insight into the chronological
sequence of events that shaped the Milky Way.

In the context of the Galactic halo, the study of binary stars acquires
additional significance. Binary stars are a fundamental component of our
Galaxy’s architecture. They serve as critical probes for understanding stellar
formation and evolution (Carney & Latham, 1987; Tokovinin et al., 2006; Latham
et al., 2002; Raghavan et al., 2010; Moe et al., 2019; Price-Whelan et al.,
2020; Mazzola et al., 2020; Bashi et al., 2023; Gaia Collaboration et al.,
2023b; Chen et al., 2024). The frequency and characteristics of binary systems
in this region provide insights into the conditions of star formation and the
dynamical processes in the early Galaxy and the smaller, assimilated dwarf
galaxies. For a recent review of results in the study of binary star research,
particularly those facilitated by Gaia, we refer the reader to El-Badry
(2024).

Previous studies have predominantly focused on binary systems within the disc
and bulge of the Milky Way or have treated the halo as a homogeneous entity.
One exception might be the recent work of Hwang et al. (2022), which focused
on wide-binary stars fractions.

In recent years, the field of astrophysics has seen an influx of data from
wide-field spectroscopic surveys such as LAMOST (Cui et al., 2012), GALAH (De
Silva et al., 2015) and APOGEE (Majewski et al., 2017). These surveys have
been instrumental in measuring the Radial-Velocities (RVs) and chemical
abundances of hundreds of thousands of stars. A key aspect of some of these
surveys is the collection of multi-epoch data for individual stars. This
approach has enabled a statistical analysis of RV variability as a function of
various stellar properties, including metallicity (Tian et al., 2018; Moe et
al., 2019; Price-Whelan et al., 2020) and $\alpha$ elements (Badenes et al.,
2018; Mazzola et al., 2020; Niu et al., 2021).

For this study, we meticulously selected our sample from the Gaia Data Release
3 (DR3 Gaia Collaboration et al., 2023a), focusing on specific criteria to
ensure the relevance and accuracy of our analysis. While Gaia DR3 does not
include each source’s individual epoch RV measurements (Katz et al., 2023), we
can exploit the available RV variability information such as the standard
deviation of the epoch radial velocity measurements of the individual source
$\sigma_{\mathrm{RV}}$, or the $\Delta\mathrm{RV_{max}}$, the maximum detected
shift in the radial velocities to assess binary frequency in a statistical
framework (Maoz et al., 2012; Badenes et al., 2018; Moe et al., 2019; Mazzola
et al., 2020; Andrew et al., 2022).

While it is generally hypothesized that _accreted_ and _in-situ_ halo stars
might exhibit different binary fractions due to their distinct formation
histories, empirical evidence supporting this notion is sparse. This research
aims to fill this gap by employing a novel methodological approach to estimate
the close binary fractions among these stars thereby shedding light on the
diverse evolutionary pathways that have shaped the Galactic halo.

The paper is structured as follows. In Section 2, we present our sample
selection of bright Gaia Main-Sequence (MS) sources. Section 3 discusses the
dependence between $\sigma_{\mathrm{RV}}$ and the stellar magnitude seen in
the Gaia sample. In Section 4, we demonstrate our selection procedure of
_accreted_ and _in-situ_ halo stars. Section 5 presents our methodological
approach using Bayesian modelling of the binary fractions. Sections 6 and 7
give our main results on the binary fractions in the _accreted_ and _in-situ_
, as well as their dependence on metallicity. Section 8 summarizes our results
and discusses their possible implications.

##  2 selecting a sample of bright main-sequence sources

We targeted stars within $3$ kpc of the solar vicinity with measurable radial
velocities (Katz et al., 2023). The query used to obtain data from the Gaia
DR3, was designed to filter and isolate a specific subset of MS and sub-Giant
objects, focusing on those pertinent to our study of binary fractions in the
Galactic halo. Consequently, we limited our sample to stars with magnitudes
estimated from the integration of the Radial Velocity Spectrometer
$G_{\mathrm{RVS}}$ brighter than $12$ mag using rv_method_used=1.
$G_{\mathrm{RVS}}$ is estimated using the flux in the RVS spectra (Sartoretti
et al., 2023). This magnitude constraint was applied to ensure the reliability
and quality of the RV data, as the epoch radial velocities of stars beyond
this threshold were not considered reliable enough to derive the combined
radial velocities as the median of the epoch radial velocities (Katz et al.,
2023).

We further selected stars with well-constrained parallax measurements
(parallax_over_error > 5 and parallax > 0) , ensuring high-quality distance
estimates. To minimize contamination from disc stars, we selected stars with a
galactic latitude ($b$) satisfying the condition $|b|>10^{\circ}$. This
criterion was essential to focus our study on halo stars and reduce the
inclusion of stars from the galactic disc with different dynamical properties.

We used Gaia extinction values to accurately plot the Color-Magnitude Diagram
(CMD) of our sample. The inclusion of stars with non-null ag_gspphot and
ebpminrp_gspphot values ensure that we have reliable photometric extinction-
corrected data, crucial for accurate luminosity and temperature determination.

Within the CMD, we employed a polygonal selection method to isolate mainly
dwarf stars. Our decision to focus on MS stars stemmed from our objective to
study stellar variability caused by the motion in a binary system. Using this
selection, we avoided the added complexity in cases of evolved stars with
their variable and evolutionary transient stages as well as the impact of
stellar evolution on binary fraction.

To further constrain our sample to select stars that were observed with good
coverage in epoch phase-space, we selected sources with
rv_visibility_periods_used larger than 8 where a visibility period is a group
of transits separated from other groups by a gap of at least $4$ days (Gaia
Collaboration et al., 2023a; Katz et al., 2023).

Figure 1: Extinction corrected Colour-Magnitude Diagram (CMD) of available
Gaia sources with combined radial velocities. Left panel displays points
density with colour scale represents the number of sources falling into each
bin. The green polygon-shaped dashed line (see text) delineates our border of
the selected MS sample. Right panel shows a zoomed-in view with MIST Isochrone
for single (q=0) and twin binary (q=1) systems in different range of
metallicities and ages.

In the left panel of Fig. 1, we show a CMD of our sample of $N=1,153,051$ MS
and Sub-Giants sources (out of a total stellar sample of $N=2,174,836$
sources) delineated by green polygonal borders.111The used polygon bounds are
the following (BP_RP, Abs. G): (0.4, 4.5), (0.4, 2.5), (0.75, 2.5), (1, 4),
(1.75, 7), (1, 7), (0.4, 4.5). The colour scale represents the number of
sources falling into each bin. The right panel of Fig. 1 displays a zoomed
selection of the polygonal area, where we have superimposed the MESA (Paxton
et al., 2011) Isochrones and Stellar Tracks models (MIST ; Dotter, 2016; Choi
et al., 2016) Isochrones interpolation for single stars (i.e. binary mass
ratio, $q=0$) and twin binaries ($q=1$). For single stars, the tracks span a
range of metallicities, accurately tracing the progression of metal-poor to
metal-rich stars within the MS. For the [Fe/H]=-1 and [Fe/H]=-2 a $12$ Gyrs
Isochrone was selected while for the [Fe/H]=0 we selected a $4.5$ Gyrs
Isochrone to better characterise the current population of stars in these
metallicity bins. Additionally, for the case of twin binaries, the
corresponding track is plotted $0.75$ magnitude above the single star line
(Hurley & Tout, 1998) representing the photometric appearance of an equally
bright companion in the CMD, effectively mimicking the luminosity of a normal
metal-rich dwarf.

As can be seen by Fig. 1, the polygon borders are meticulously chosen to
include mainly dwarf stars and capture the sequence of metal-poor stars,
typically found at the lower extremity of the MS. Notably, the binary sequence
of metal-rich stars (at BP_RP > 1) is clearly discernible within our selection
and corresponds nicely to the twin stars [Fe/H] = 0 Isochrone. While most
stars in this region are not halo stars, their inclusion was deemed essential
for completeness. As a final remark, we note that in the case of a metal-poor
star with an equally bright companion (twin) the photometrically position on
the CMD, appear as a normal metal-rich dwarf.

##  3 Stellar variability: radial velocity uncertainty

To get the standard deviation of the epoch radial velocity measurements of the
individual bright source $\sigma_{\mathrm{RV}}$, we used the following
expression:

$\sigma_{\mathrm{RV}}=\sqrt{\frac{2N_{\mathrm{RV}}}{\pi}\left(\zeta_{\mathrm{RV}}^{2}-0.11^{2}\right)}$
(1)

where $\zeta_{\mathrm{RV}}$ is the uncertainty on the median of the epoch
radial velocities (radial_velocity_error) to which a constant shift of $0.11$
km/s was added to take into account a calibration floor contribution (Andrew
et al., 2022; Katz et al., 2023), and $N_{\mathrm{RV}}$ is the number of
transits used to derive the median radial velocity (rv_nb_transits).

Fig. 2 shows the distribution of $N=1,153,051$ MS sources in
$\sigma_{\mathrm{RV}}$ as function of $G_{\mathrm{RVS}}$. As might be
expected, the typical standard deviation of the epoch radial velocity
measurements increases towards fainter sources. Above this observable curve,
we find a further wider spread of sources of suspected binary systems with
high $\sigma_{\mathrm{RV}}$ values.

Figure 2: RV uncertainty $\sigma_{\mathrm{RV}}$ vs $G_{\mathrm{RVS}}$ of our
MS sample (sources inside green polygon of Fig. 1). The colour scale
represents the number of sources falling into each bin. A clear monotonic
trend between $\sigma_{\mathrm{RV}}$ and $G_{\mathrm{RVS}}$ is evident, with
points well above this trend being suspected binary systems.

To demonstrate this expected trend, we show in Fig. 3 the
$\sigma_{\mathrm{RV}}$-$G_{\mathrm{RVS}}$ scatter of a sample of three types
of systems: (i) fraction of planet-host stars listed in the NASA Exoplanet
Archive222https://exoplanetarchive.ipac.caltech.edu/index.html with a flag
suggesting detection by radial-velocity and can be approximated as a well-
constrained sample of single stars; (ii) a sample of single-lined
spectroscopic binaries (SB1) listed in the Gaia DR3 NSS catalogue (Gaia
Collaboration et al., 2023b) with a clean-score above 0.578 based on Bashi et
al. (2022) metric and (iii) eclipsing binaries (EB) listed in the first Gaia
catalogue of EB candidates (Mowlavi et al., 2023). As expected, the single
stars sample sits nicely on the exponential curve and have the lowest
$\sigma_{\mathrm{RV}}$ for a given $G_{\mathrm{RVS}}$, while the SB1s are
found above this curve, and the EBs with typically shorter periods and larger
semi-amplitudes have the highest $\sigma_{\mathrm{RV}}$ values.

In section 5, we will further elaborate on this dependence in developing our
model to estimate binary fractions and the typical (mean)
$\sigma_{\mathrm{RV}}$ value as a function of $G_{\mathrm{RVS}}$.

Figure 3: RV uncertainty $\sigma_{\mathrm{RV}}$ vs $G_{\mathrm{RVS}}$ of
single stars (cyan), single-line spectroscopic binaries listed in the Gaia DR3
NSS catalogue (magenta) and eclipsing binaries listed in the first Gaia
catalogue of eclipsing-binary candidates (green) after cross match with our MS
sample. Gray background colour scale represents the number of sources falling
into each bin of the general sample similar to Fig. 2.

##  4 _accreted_ and _in situ_ samples

In addition to the previously mentioned criteria, we implemented further cuts
to refine our selection of halo stars. Specifically, we applied a velocity
threshold, selecting stars with a total velocity
($V_{\mathrm{tot}}=\sqrt{U^{2}+V^{2}+W^{2}}$) greater than $220$ km/s. This
threshold is based on the understanding that the halo shows close to net zero
rotation and therefore in the Solar neighborhood, the velocity of its stars is
on average offset from that of the disc’s stars by the local value of the
circular velocity (Nissen & Schuster, 2010; Helmi et al., 2018; Belokurov et
al., 2018; Bashi & Zucker, 2022; Lane et al., 2022).

To accurately estimate the dynamical properties of our selected stars, we used
the Galpy package (Bovy, 2015). Using Galpy, we calculated the orbital energy
($E$) and the z-component of angular momentum ($L_{z}$) for each star in our
sample. These calculations were based on the assumption of the Milky Way’s
potential as described in the MWPotential2014 model which is a well-
established model that approximates the gravitational potential of the Milky
Way, considering the contributions from the bulge, disc, and halo. We plot all
energies offset by the potential of MWPotential2014 at infinity so that $E<0$
are bound orbits (e.g. Lane et al., 2022).

In the $E-L_{z}$ plane, the samples are confined by a bounding parabola
defined by the escape velocity as a function of changing energy at a fixed
Galactocentric radius.

In Fig. 4 we show $N=1,928$ stars on the $E-L_{z}$ plane that were left after
applying the above cuts and a background grey contours marking the position of
the overall MS sample which mainly comprised of disc stars. We mark the Sun’s
position at $E_{\odot}=-1.04\times 10^{5}~{}\mathrm{km^{2}/s^{2}}$ and use it
to gauge the boundary energy we select to separate between the _accreted_ and
_in-situ_ samples. Follow Belokurov & Kravtsov (2023), and to reduce
contamination, we selected conservative cuts on $E$ and affiliated _accreted_
stars as sources with $E>E_{\odot}+0.2$ and $L_{z}<2$ and _in-situ_ stars as
$E<E_{\odot}-0.2$ and $L_{z}<0$.

By incorporating these velocity and dynamic cuts, we were able to identify
stars that are characteristic of the Milky Way’s halo more precisely, and were
left with $N=181$ _accreted_ and $N=291$ _in-situ_ halo stars. This refined
selection ensures that our analysis of binary fractions is focused on the
intended stellar population, thereby enhancing the reliability of our binary
fractions estimates.

Figure 4: $E-L_{z}$ values of our MS Gaia sources, with contour background
grey points mark the main location of disc stars and blue points mark stars
with $V_{\mathrm{tot}}>220$ km/s (see text for full selection cuts). Blue
points inside the dashed green box mark our selected sample of ’ _in-situ_ ’
sources, while blue points inside the dashed orange box mark the sample of
_accreted_ stars. The black fine-dashed curve corresponds to the maximal
angular momentum at fixed energy, i.e. orbits with circular velocity. $\odot$
marks the location of the Sun in the chosen potential.

##  5 Methodology: Modeling Binary Fractions using the RV uncertainty
distribution

Single stars typically exhibit minimal RV variation, primarily due to factors
like stellar activity and instrumental noise, resulting in low
$\sigma_{\mathrm{RV}}$ values that follow a trend similar to the cyan points
of Fig. 3 with fainter stars having larger variation. On the other hand, in
binary systems, the orbital motion of each star causes periodic changes in
their radial velocities. For a given binary semi-amplitude, $K$, in the most
simplified scenario of a circular orbit, assuming the measurements are
uniformly distributed across the orbital phase, the added RV uncertainty from
the binary nature of the system would be
$\sigma_{\mathrm{RV}}^{b}=K/\sqrt{2}$. Thereby, the standard deviation of the
epoch radial velocity measurements value is a direct indicator of the binary
nature of the system (Maoz et al., 2012; Badenes et al., 2018; Moe et al.,
2019; Mazzola et al., 2020; Andrew et al., 2022), where a larger
$\sigma_{\mathrm{RV}}$ typically suggests a closer or more massive companion
as $K\propto P^{-1/3}q(1+q)^{-2/3}$ where $P$ is the binary orbital period and
$q$ is the binary mass ratio.

By analyzing the $\sigma_{\mathrm{RV}}$-$G_{\mathrm{RVS}}$ distribution for a
given sample of Gaia stars, we can statistically estimate the binary fraction.
Given a typical $\sigma_{\mathrm{RV}}$ of $\sim 5$ km/s at
$G_{\mathrm{RVS}}$=12, our method is sensitive to close-binaries with orbital
period $P\lesssim 1000$ days for a solar mass twin binary system.

A dataset consisting predominantly of single stars would exhibit a narrow
distribution centred around, and following the trend of dense sources shown in
Fig. 2, thereby reflecting minimal radial velocity variation. In contrast, a
dataset with a high binary fraction would show a broader distribution,
extending towards higher $\sigma_{\mathrm{RV}}$ values for a given
$G_{\mathrm{RVS}}$, indicative of the significant radial velocity changes
caused by binary companions. A mixed population of single and binary stars
would likely display a bimodal distribution, with two distinct peaks
corresponding to the two populations.

In most cases, our method treats higher-order multiple star systems, such as
triples, as binary stars. In this scenario, we might be sensitive to the close
binary substructure where the ratio of the outer period $P_{\mathrm{out}}$ to
the inner period $P_{\mathrm{in}}$, adjusted for eccentricity
$e_{\mathrm{out}}$, is expected to have
$P_{\mathrm{out}}\left(1-e_{\mathrm{out}}\right)^{3}/P_{\mathrm{in}}>5$
(Tokovinin, 2004) while we may miss the RV variation of the outer companion.

To model the population of sources in log $\sigma_{\mathrm{RV}}$, we define a
density function, which is a sum of two Gaussian distributions to model the
populations of single and binary stars. To properly model the fractions of
single and binary stars, we multiply each Gaussian by the expected fractions,
where $F$ marks the binary fraction. The final density function will then be

$\displaystyle f(x|G_{\mathrm{RVS}};\theta)$
$\displaystyle=(1-F)\cdot\mathcal{N}_{s}\left(x|\mu_{s}(G_{\mathrm{RVS}}),\sigma_{s}^{2}\right)$
(2)
$\displaystyle+F\cdot\mathcal{N}_{b}\left(x|\mu_{b}(G_{\mathrm{RVS}}),\sigma_{b}^{2}\right).$

where $\theta={\left(F,a,b,G_{\mathrm{min}},d,\sigma_{s},\sigma_{b}\right)}$,
defines our parameter space and
$\mathcal{N}_{s}\left(x|\mu_{s}(G_{\mathrm{RVS}}),\sigma_{s}^{2}\right)$,
$\mathcal{N}_{b}\left(x|\mu_{b}(G_{\mathrm{RVS}}),\sigma_{b}^{2}\right)$ are
the Gaussian density functions of the single and binary populations with means
$\mu_{s}(G_{\mathrm{RVS}})$, $\mu_{b}(G_{\mathrm{RVS}})$ and standard-
deviations $\sigma_{s}$, $\sigma_{b}$ respectively.

To model the dependence between $\sigma_{\mathrm{RV}}$ and $G_{\mathrm{RVS}}$,
we used an exponent function and defined $\mu_{s}(G_{\mathrm{RVS}})$ as

$\mu_{s}(G_{\mathrm{RVS}};a,,b,G_{\mathrm{min}})=a+e^{b\cdot(G_{\mathrm{RVS}}-G_{\mathrm{min}})}$
(3)

where $a,~{}b,~{}G$ are free parameters which control the shape of the
expotential function. Similarly, we defined the binary mean,
$\mu_{b}(G_{\mathrm{RVS}})$ in quadrature as

$\mu_{b}(G_{\mathrm{RVS}};a,,b,G_{\mathrm{min}},d)=\log_{10}{\left(\sqrt{\left(10^{\mu_{s}(G_{\mathrm{RVS}};a,,b,G_{\mathrm{min}})}\right)^{2}+d^{2}}\right)}\\\
$ (4)

where $d$ is a free parameter that accounts for the extra variability in
$\sigma_{\mathrm{RV}}$ as a result of the system’s binary nature.

In our analysis, we sought to account for potential variations in
$\sigma_{\mathrm{RV}}$ attributed to line-broadening where early type stars at
a fixed $G_{\mathrm{RVS}}$ having larger errors (Andrew et al., 2022). To
explore this dependence, we extend our model by introducing extra factor to
account for the relation between $\sigma_{\mathrm{RV}}$ and BP_RP. Despite
this, the outcomes from our ’simplified’ model were found to align with those
from a model that explicitly incorporates spectral type. This agreement
suggests that the parameters $\sigma_{s}$ and $\sigma_{b}$ within our model
effectively encapsulate the influence of spectral type, thereby moderating the
dispersion of data points around the fitted single and binary star models.

Table 1: Prior distributions of model parameters of binary fraction ($F)$, and single and binary Gaussian’s means ($a,~{}b,~{}G_{\mathrm{min}},~{}d$) and stds ($\sigma_{s},~{}\sigma_{b}$).  Parameter | Prior
---|---
$F$ | $\log\mathcal{U}(0.01,1)$
$a$ | $\mathcal{U}(-5,5)$
$b$ | $\log\mathcal{U}(10^{-3},1)$
$G_{\mathrm{min}}$ | $\mathcal{U}(0,12)$
$d$ | $\log\mathcal{U}(10^{-4},100)$
$\sigma_{s}$ | $\log\mathcal{U}(10^{-3},10)$
$\sigma_{b}$ | $\log\mathcal{U}(10^{-2},10)$

Finally, as the stars whose combined RV formal uncertainty was larger than
$40$ km/s were flagged as invalid and removed from the final Gaia DR3
catalogue (Katz et al., 2023), we used a truncated normal distribution to cope
with this abrupt cut in the observed distributions.

The final likelihood expression would then be

$\mathcal{L}=\prod_{i=1}^{N}f(x^{i}|G_{\mathrm{RVS}}^{i};\theta).$ (5)

Figure 5: Posterior histograms of binary fractions in _accreted_ (orange) and
_in-situ_ (green) samples based on the selection presented in section 4.

To test our algorithm, we applied it on different sets of simulated datasets
and used an MCMC routine with $40$ walkers and $10^{4}$ steps, using
uninformative priors on the seven free parameters of our model as listed in
Table 1. We used the Python emcee package (Foreman-Mackey et al., 2013) to
find the parameter values (and their uncertainties) that maximize the sample
likelihood.

In the first test set, we attempted to find the best-fit model parameters
${\left(F,a,b,G_{\mathrm{min}},d,\sigma_{s},\sigma_{b}\right)}$ from a
simulated sample based on the defined density function of Equation 2.

In the second test set, we simulated a population of single vs short-period
binary systems using a random selection of binary separations, mass ratios,
eccentricity distribution inclination and phase (e.g. Mazzola et al., 2020)
and attempted to get the best-fit parameters of
${\left(F,a,b,G_{\mathrm{min}}\right)}$.

In the third test set, we used a mixed sample of known single stars and Gaia
SB1s (Bashi et al., 2022; Gaia Collaboration et al., 2023b) based on the
subsamples used in Fig. 3 and attempted to get the best-fit parameters of
${\left(F\right)}$.

We found good agreement between the posterior distributions and the actual
simulated parameters in all test cases, thereby supporting our Bayesian model
and encouraging us to proceed and estimate binary fractions on real datasets.

##  6 Binary fractions in halo stars

###  6.1 The _accreted_ and _in-situ_ subsamples

Given our selection of $N=181$ and $N=291$ _accreted_ and _in-situ_ sources
respectively we list in Table 2 the best fit posterior values of our Bayesian
model (median and $\%16,\%84$ quantiles uncertainties). We present in Fig. 5
the posterior distributions of the _accreted_ (orange) and _in-situ_ (green)
binary fractions. We find a $\sim 2\sigma$ difference in favor of a higher
binary fraction in _accreted_ stars, $F_{\mathrm{accreted}}=0.238\pm 0.038$,
compared to _in-situ_ stars, $F_{\mathrm{in-situ}}=0.182\pm 0.032$.

In Fig. 6, we show a scatter plot of $\sigma_{\mathrm{RV}}$ vs.
$G_{\mathrm{RVS}}$ for the _accreted_ and _in-situ_ samples. The red dotted
lines mark our best-fit posterior distributions of the mean values using
Equation 3 which describes the relation between $\sigma_{\mathrm{RV}}$ and
$G_{\mathrm{RVS}}$ in the case of a single source solution. Overall, the best
fit posteriors ($a$, $b$, $G_{\mathrm{min}}$), which determine the curves
shape of the _accreted_ and _in-situ_ samples are consistent within $1\sigma$
as can be seen in Table. 2.

To explore the distributions of $\sigma_{\mathrm{RV}}$ after removing the
effect of stellar $G_{\mathrm{RVS}}$, we show in Fig. 7 the ’flattened’
normalised distributions after substracting the observed
$\sigma_{\mathrm{RV}}$ values with the best fit model
$\mu_{s}(G_{\mathrm{RVS}};a,,b,G_{\mathrm{min}})$. Comparing these
distributions might also shed some light on the nature of these binary
systems. As expected both histograms have clear peaks around zero, marking the
distributions of single sources with standard-deviation $\sigma_{s}$. The
extended tail of binaries is also evident where we find the _accreted_ sample
having more sources at high
$\sigma_{\mathrm{RV}}-\mu_{s}(G_{\mathrm{RVS}};a,,b,G_{\mathrm{min}})$ values
compared to the _in-situ_ sub-sample. This observed trend is also supported by
our posterior model results where the _accreted_ best-fit model on $d$ and
$\sigma_{b}$, which are responsible to mark the position and spread of the
binary sequence, is larger than in the _in-situ_ case.

Figure 6: $\sigma_{\mathrm{RV}}$ vs $G_{\mathrm{RVS}}$ of the _accreted_ (left
panel) and _in-situ_ (right panel) samples. Dotted red lines mark the best-fit
Bayesian model of our posterior distribution based on Equation 3.

A significant limitation encountered in utilizing the Gaia sample pertains to
the obscured insights into the metallicity content of the two halo samples.
Metallicity is known to anti-correlate with binary fraction (Badenes et al.,
2018; Moe et al., 2019; Price-Whelan et al., 2020; Mazzola et al., 2020).
Consequently, any disparity in the metallicity content between the _accreted_
and _in-situ_ samples may lead to variations in binary fractions that are not
inherently related to their dynamical properties. To further explore this, we
examine in section 7 a subset of our samples with reported metallicity
estimates using the APOGEE DR17 catalogue (Abdurro’uf et al., 2022).

Table 2: Posterior values of _accreted_ and _in-situ_ samples.  Parameter | _accreted_ | _in-situ_
---|---|---
$N$ | $181$ | $291$
$F$ | $0.2384\pm 0.038$ | $0.182\pm 0.032$
$a$ | $1.84_{-0.55}^{+0.44}$ | $1.59_{-0.39}^{+0.45}$
$b$ | $0.163_{-0.027}^{+0.059}$ | $0.174_{-0.034}^{+0.042}$
$G_{\mathrm{min}}$ | $6.0_{-2.3}^{+2.5}$ | $7.0_{-2.4}^{+1.8}$
$d$ | $20.0_{-5.0}^{+8.2}$ | $14.6_{-3.2}^{+4.0}$
$\sigma_{s}$ | $0.1150_{-0.0091}^{+0.0094}$ | $0.1252_{-0.0082}^{+0.0092}$
$\sigma_{b}$ | $0.601_{-0.075}^{+0.097}$ | $0.519_{-0.055}^{+0.066}$
Figure 7: ’Flattened’ normalised distributions of _accreted_ (orange) and _(_
in-situ) (green) sources after substracting the observed
$\sigma_{\mathrm{RV}}$ values with the best fit model
$\mu_{s}(G_{\mathrm{RVS}};a,,b,G_{\mathrm{min}})$.

###  6.2 Cuts in $E$ and $L_{z}$

We next attempted to estimate binary fractions using different cuts in $E$ and
$L_{z}$ that might be affiliated with different Galactic components (Helmi &
White, 1999; Belokurov et al., 2018; Myeong et al., 2019; Matsuno et al.,
2019; Naidu et al., 2020; Belokurov & Kravtsov, 2022, 2023). While _accreted_
stars and especially the GS/E, are more common at higher $E$ values in the
solar neighbourhood (Belokurov et al., 2018; Helmi et al., 2018; Haywood et
al., 2018; Naidu et al., 2020), other substructures also have their expected
places on the $E-L_{z}$ diagram. For example, at high $L_{z}$, we are evident
to the extended high-$\alpha$ sequence of disc stars as well as _Splash_ stars
(Bonaca et al., 2017; Gallart et al., 2019; Di Matteo et al., 2019; Belokurov
et al., 2020). This sample might be contaminated by a higher fraction of
metal-rich stars. Similarly, at low energy values of retrograde orbits, we
should expect to have a higher population of _Aurora_ , _in-situ_ stars
(Belokurov & Kravtsov, 2022; Conroy et al., 2022; Chandra et al., 2023; Zhang
et al., 2023).

In the upper panel of Fig. 8 we show the binary fractions using different cuts
in energy, $E^{\mathrm{cut}}$ that separates our sample to high- and low-
energy sources (see red and violet points respectively). Our analysis
demonstrates that below a critical energy threshold
($E<E^{\mathrm{cut}}\approx E_{\odot}$), binary fractions remain relatively
similar. This consistency can be explained by the high contamination rate of
_in-situ_ stars which dominates the sample of high-energy stars, making the
two subsamples of low- and high-energy similar in stellar properties.
Conversely, above this energy limit ($E>E^{\mathrm{cut}}\approx E_{\odot}$),
we observe a noticeable increase in binary fractions. This trend intensifies
at higher energy levels, where the sample likely comprises a purer selection
of _accreted_ stars. These results align with findings from our previous
section where we estimated binary fractions in a conservative subsample of
_accreted_ and _in-situ_ stars. The increase in binary fraction at higher
energy bins suggests a more significant presence of dynamically heated,
_accreted_ structures especially of the GS/E, with lower metallicities.

The variation of binary fractions with angular momentum also presents an
interesting trend. In the bottom panel of Fig. 8 we show the binary fractions
as function of the angular momentum cut-off ($L_{z}^{\mathrm{cut}}$) that
separates our sample to retrograde-prograde and slow-fast rotation sources
(see red and violet points respectively). We find a consistent higher binary
fractions in sources with $L_{z}<L_{z}^{\mathrm{cut}}$ compared to the
complement subgroup. In stark contrast, the most prograde sources,
predominantly associated with the high-$\alpha$ disc and _Splash_ (Bonaca et
al., 2017; Gallart et al., 2019; Di Matteo et al., 2019; Belokurov et al.,
2020), exhibit significantly lower binary fractions than the overall halo
sample. Possible explanations to this difference might be the higher metal
content associated with this subsample compared to the slow and retrograde
sources.

Figure 8: Binary fractions in halo stars and their corresponding uncertainties
using different cuts in energy ($E^{\mathrm{cut}}$, upper panel) and angular
momentum ($L_{z}^{\mathrm{cut}}$, bottom panel). Red points display sources
with $E$ ($L_{z}$) higher (lower) than the threshold, while violet points
display the complement subsets.

##  7 Metallicity dependence

In this section, we delve into the relationship between metallicity and binary
star fractions while considering the samples of _accreted_ and _in-situ_ halo
stars. We used the APOGEE DR 17 catalogue (Abdurro’uf et al., 2022) to get
estimates of stellar metallicity and cross-matched our Gaia MS sample with it.
To minimise effects by stellar mass on binary fraction and to select a
representative sample similar to our halo sample, we limited our selection to
sources with APOGEE effective temperature (TEFF) between $4600-6200$ K and
stellar mass, approximated using the Torres et al. (2010) relation, within the
$0.5-1.2M_{\odot}$ range.

Figure 9: The Gaia -APOGEE cross-match sample. Left panel shows the $E-L_{z}$
values while the right panel the [Fe/H]-[Al/Fe] scatter. Green points, mark
_in-situ_ stars based on previous section kinematic selection, orange points
mark the _accreted_ sample. Grey contour background points mark the main
location of the disc and other _in-situ_ substructures.

On the one hand, while there are larger catalogues with reported metallicity
estimates (e.g. Andrae et al., 2023), it is extremely important when dealing
with metal-poor binary stars to use a high spectroscopic derivation of the
metal content to avoid over-estimation of the derived metallicities. On the
other hand, performing our binary search analysis solely using APOGEE stars
would not be applicable as most targets had only a single exposure with a
small fraction of sources having 2-4 RV epoch observations (Abdurro’uf et al.,
2022; Mazzola et al., 2020). These limited number of epochs, is insufficient
to get good estimates of $\sigma_{\mathrm{RV}}$ compared to Gaia’s DR3 typical
$20-30$ epochs (Katz et al., 2023).

As part of our initial filtering, we removed stars with measurements that are
believed to be problematic, as indicated by the following general flags:
VERY_BRIGHT_NEIGHBOR, LOW_SNR, as well as these element pipeline analysis
flags: FE_H_FLAG, MG_FE_FLAG, AL_FE_FLAG. However, we found that even while
using these cuts, there were a few cases where the extracted metallicity
estimates listed in APOGEE where underestimated. In most of these cases, the
outliers correlated strongly with high values of $\sigma_{\mathrm{RV}}$,
suggesting that the binary nature of theses sources, and especially in cases
of high mass-ratios, have biased the elemental estimation. These cases are
common particularly in blended sources and multiple-star systems with more
than one star that contribute significantly to the overall flux. A similar
issue was also highlighted in Price-Whelan et al. (2020), where, in their
analysis aimed at characterizing close binaries, it was assumed that their
sample consisted exclusively of single-line sources.

Consequently, to remove these outliers, we found the value on RV_CCFWHM, which
marks the Full-Width Half-Maximum (FWHM) of the cross-correlation peak from
the combined vs. best-match synthetic spectrum, as a good indicator of ’well-
behaved’ solutions in cases with value below $150$ km/s. The cross-correlation
technique is employed to match observed stellar spectra against a library of
synthetic spectra. The FWHM of the cross-correlation peak, denoted by
RV_CCFWHM, quantifies the width of this peak at half of its maximum height,
reflecting the degree of match between the observed and synthetic spectra.
Thereby a peak that is too broad might indicate an imprecise match.

The cross-match yielded $N=48,484$ sources (disc and halo) in common, with
only $N=9$ and $N=27$ _accreted_ and _in-situ_ sources respectively having
using section 4 subsamples selection. On the left panel of Fig. 9, we show a
scatter plot of the sources on $E-L_{z}$ plane similar to Fig. 4 while on the
right panel we show the distribution on the [Fe/H]-[Al/Fe] diagram. As
expected, our kinematic selected halo sub-samples dominate the metal-poor
fraction of stars, with _accreted_ sources having lower [Al/Fe] values
compared to _in-situ_ stars (Hasselquist et al., 2021; Belokurov & Kravtsov,
2022, 2023; Deason & Belokurov, 2024). In addition, we note on a few other
_in-situ_ metal-poor prograde stars seen in this plot (grey points), probably
affiliated with the _Aurora_ population which while assumed to be
kinematically hot, it has an approximately isotropic velocity ellipsoid and a
modest net rotation (Belokurov & Kravtsov, 2022; Zhang et al., 2023).

Using the APOGEE sample, we can attempt to estimate binary fractions in bins
of metallicity. The black crosses in Fig 10 present a detailed analysis of the
binary fraction across various metallicity bins. We focus on the median values
and the $16\%$, $84\%$ quantile values of the binary fractions posteriors. In
addition we plot by red and violet crosses the fraction of retrograde
($L_{z}<0$; $N=51$) and prograde ($L_{z}\geq 0$; $N=83$) halo sources
respectively.

The low number of _accreted_ ($N=9$) and _in-situ_ ($N=27$) halo stars in our
Gaia-APOGEE cross-matched subsample is not sufficient to produce valuable
estimates of binary fractions using our Bayesian model. Instead, we decided to
estimate the median and $16\%$, $84\%$ quantile values on the metallicty
values in each halo subsample and place the binary fractions found in section
6 as the expected values. Consequently we mark by orange and green boxes the
_accreted_ and _in-situ_ binary fractions to stress that these are not a
direct calculation of the binary fractions of corresponding samples shown in
Fig. 9.

Examining the trend in filed stars’ binary fraction as a function of [Fe/H],
Fig 10 corroborates previous research findings that have consistently shown
higher binary fractions in metal-poor stars (Badenes et al., 2018; Moe et al.,
2019; Mazzola et al., 2020; Price-Whelan et al., 2020). This trend is
prominently visible across the different metallicity bins, reinforcing the
significant role of metallicity in the formation of binary stars.

In our context, one further notable observation from Fig 10 is the consistency
of the _accreted_ and _in-situ_ with the field stars trend, suggesting that
the anti-correlation in close-binary fractions with metallicity continues into
the domain of halo stars. Secondly, we find an interesting discrepancy in
field binary fractions compared to the halo prograde stars. Similar
discrepancy is also evident when comparing the trend of field binary fractions
of Fig. 10 to the the highest prograde bin ($L_{z}\geq
L_{z}^{\mathrm{cut}}=0.5$) in the bottom panel of Fig. 8. Assuming the
prograde bin includes high-$\alpha$ disc and _Splash_ stars, we might expect
it to be populated mainly by sources with [$\alpha$/Fe] $\geq 0.2$. However,
our bin estimate of a binary fraction of $F\approx 0.10$, is lower than the
expected value observed in the trend of Fig. 10 for similar metallicity.
Nonetheless, given the Gaia-APOGEE sample is dominated by thin disc stars,
most of the source in this diagram have lower-$\alpha$ values and have a
higher stellar mass compared to thick disc stars. Both of these parameters
inspire a higher binary fraction (Raghavan et al., 2010; Moe et al., 2019;
Mazzola et al., 2020) that can mitigate this effect.

##  8 Conclusions

This study provides valuable insights into close binary star fractions within
the Galactic _accreted_ and _in-situ_ halo populations, which stands as a
pivotal component in understanding the dynamical history and formation of the
Milky Way Galaxy.

By leveraging the comprehensive data from Gaia DR3 (Gaia Collaboration et al.,
2023a), our methodological approach, utilized a novel Bayesian framework to
analyse the radial velocity uncertainty distribution of Gaia RVS sources (Katz
et al., 2023; Sartoretti et al., 2023). Using this approach, we were able to
model the populations of single and binary stars as a sum of two Gaussian
distributions allowing a nuanced understanding of close-binary star fractions
in different Galactic environments using only Gaia data.

Our findings reveal a higher binary fraction in _accreted_ stars, compared to
_in-situ_ halo stars. On the one hand, this observation suggests that the
processes influencing binary star formation and evolution might differ across
these distinct populations. On the other hand, given the slight lower metal
content in _accreted_ stars compared to _in-situ_ halo stars, the difference
in binary fractions might be explained by the difference in metallicity. A
similar conclusion was also made for wide-binary fractions in Hwang et al.
(2022), although their selection of _accreted_ and _in-situ_ sources was
different than the one we used. Furthermore, their selection of stars on the
$E-L_{z}$ plane was more contaminated causing the the binary fractions to be
much more consistent in any case.

Figure 10: Posterior values of binary fractions using APOGEE data in bins of
[Fe/H]. Black crosses mark the whole Gaia-APOGEE sample. Red and violet
crosses mark the fractions in retrograde ($L_{z}<0$) and prograde ($L_{z}\geq
0$) halo sources respectively. We also mark with orange and green boxes the
_accreted_ and _in-situ_ binary fractions based on the values estimated in
section 6.

In our study, we also explored binary fractions across different galactic
components, employing varied energy and angular momentum cuts to differentiate
these groups. We found that, below a specific energy threshold
($E<E^{\mathrm{cut}}\approx E_{\odot}$), binary fractions remained relatively
uniform due to a high contamination rate of _in-situ_ stars that dominate the
high-energy sample, resulting in similar stellar properties between low- and
high-energy stars. Above this energy threshold, binary fractions increased,
especially in higher energy levels where the sample likely consists more of
_accreted_ stars, of the GS/E, indicating a prevalence of dynamically heated,
metal-poor structures. Additionally, when examining the variation of binary
fractions with angular momentum, we observed consistently higher binary
fractions in retrograde sources compared to their prograde counterparts. This
disparity is stark, as the most prograde sources, which are typically metal-
rich, exhibited significantly lower binary fractions than those in the overall
halo sample. These findings suggest that the observed trends in binary
fractions can provide insights into the dynamic and compositional differences
within the Galaxy.

Overall, our research corroborates previous studies, indicating higher binary
fractions in metal-poor stars (Moe et al., 2019; Price-Whelan et al., 2020;
Mazzola et al., 2020). Our analysis reveals that the binary fractions of
_accreted_ and _in-situ_ samples are consistent with the field stars trend, as
evidenced in Fig. 10. This observation suggests that the anti-correlation
between close-binary fractions and metallicity continues into the domain of
halo stars. Particularly, it underscores the influence of metallicity during
the early stages of the Milky Way’s formation and in the _accreted_ dwarf
galaxies. The higher prevalence of binary systems in metal-poor environments
suggests a fundamental aspect of binary formation in the early universe, where
lower metallicity was more common. Such subtleties in the data not only
enhance our understanding of the binary star dynamics but also pose intriguing
questions about the historical processes that have shaped these populations.

Building on this, the impact of primary mass on binary fraction adds another
layer of complexity to our analysis. Since binary fraction depends on the
primary mass of the system (Raghavan et al., 2010; Moe et al., 2019),
exploring this relationship could bias our results. To address this, we
conducted a comparative analysis of the Absolute Magnitude (Abs. G)
distributions between the two halo samples, hypothesizing that it could serve
as a proxy for primary mass. Our findings revealed that the distributions in
Abs. G between the _accreted_ and _in-situ_ samples were comparable,
suggesting that despite the challenges posed by metallicity, the primary mass,
as inferred from Abs. G., does not significantly differentiate the binary
fractions of the two halo samples. Similarly, in the analysis of the Gaia-
APOGEE metallicity bins we intentionally restricted our sample to narrow
region of effective temperature $4600-6200$ K and stellar mass within the
$0.5-1.2M_{\odot}$ range to decrease possible biases with primary mass.
However, this in any case, would not change dramatically, nor reverse, the
observed trend of higher binary fraction towards metal-poor stars. In fact, as
Moe et al. (2019) pointed out, as metal-poor stars tend to be of later type,
the trend seen of higher binary fraction in metal-poor stars should become
more significant.

Latest discoveries of dormant compact objects, such as Gaia BH3 (Gaia
Collaboration et al., 2024) and a list of Neutron star candidates (El-Badry et
al., 2024), have highlighted the critical importance of exploring binary
fractions in metal-poor regions, particularly among halo stars. These findings
suggest that metal-poor environments, often associated with older galactic
structures, may harbor a significant number of these elusive compact remnants.
The presence of such objects in binary systems can provide valuable insights
into the evolutionary paths that lead to their formation.

Future work could aim to expand the halo sample utilized in this study to
obtain more accurate estimates of the mean binary fraction values and their
uncertainties. For instance, we limited our analysis to stars with
$G_{\mathrm{RVS}}$<12 to ensure the reliability and quality of the RV data.
Nonetheless, the Gaia RVs are available to stars up-to $G_{\mathrm{RVS}}$=14.
Another potential strategy involves leveraging the RV scatter, incorporating
data from various ground-based spectroscopic surveys such as LAMOST (Cui et
al., 2012), GALAH (De Silva et al., 2015) and APOGEE (Majewski et al., 2017).
This approach would allow us to extend our sample, including a broader range
of stars suspected to be in binary systems.

As for extending our methodical approach, future works might also exploit our
Bayesian model to estimate the likelihood of individual Gaia sources to be
affiliated as single stars or as binary stars.This refined focus could
significantly enhance our understanding of the distinct characteristics and
behaviors of specific systems allowing further study of their specific nature.

## Acknowledgements

D.B. acknowledges the support of the Blavatnik family and the British council
as part of the Blavatnik Cambridge Fellowship. D.B. also extends his gratitude
to Didier Queloz for his hospitality and encouragement. We wish to thank Shay
Zucker for a helpful discussion on the Bayesian methodology and Jo Bovy for
assisting with the definition of the Galpy’s galactic potential. This work has
also made use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC; ht tps://www.cosmos.esa.int/web/gaia/
dpac/consortium). Funding for DPAC has been provided by national institutions,
in particular the institutions participating in the Gaia Multilateral
Agreement. This research also made use of TOPCAT (Taylor, 2005), an
interactive graphical viewer and editor for tabular data. This research has
made use of the NASA Exoplanet Archive, which is operated by the California
Institute of Technology, under contract with the National Aeronautics and
Space Administration under the Exoplanet Exploration Program.

## Data Availability

The research presented in this article predominantly relies on data that is
publicly available and accessible online through the Gaia  DR3 and APOGEE
archives.

## References

  * Abdurro’uf et al. (2022) Abdurro’uf et al., 2022, ApJS, 259, 35
  * Andrae et al. (2023) Andrae R., Rix H.-W., Chandra V., 2023, ApJS, 267, 8
  * Andrew et al. (2022) Andrew S., Penoyre Z., Belokurov V., Evans N. W., Oh S., 2022, MNRAS, 516, 3661
  * Badenes et al. (2018) Badenes C., et al., 2018, ApJ, 854, 147
  * Bashi & Zucker (2022) Bashi D., Zucker S., 2022, MNRAS, 510, 3449
  * Bashi et al. (2022) Bashi D., Shahaf S., Mazeh T., Faigler S., Dong S., El-Badry K., Rix H. W., Jorissen A., 2022, MNRAS, 517, 3888
  * Bashi et al. (2023) Bashi D., Mazeh T., Faigler S., 2023, MNRAS, 522, 1184
  * Belokurov & Kravtsov (2022) Belokurov V., Kravtsov A., 2022, MNRAS, 514, 689
  * Belokurov & Kravtsov (2023) Belokurov V., Kravtsov A., 2023, MNRAS, 525, 4456
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Belokurov et al. (2020) Belokurov V., Sanders J. L., Fattahi A., Smith M. C., Deason A. J., Evans N. W., Grand R. J. J., 2020, MNRAS, 494, 3880
  * Bonaca et al. (2017) Bonaca A., Conroy C., Wetzel A., Hopkins P. F., Kereš D., 2017, ApJ, 845, 101
  * Bovy (2015) Bovy J., 2015, ApJS, 216, 29
  * Carney & Latham (1987) Carney B. W., Latham D. W., 1987, AJ, 93, 116
  * Chandra et al. (2023) Chandra V., et al., 2023, arXiv e-prints, p. arXiv:2310.13050
  * Chen et al. (2024) Chen X., Liu Z., Han Z., 2024, Progress in Particle and Nuclear Physics, 134, 104083
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Conroy et al. (2022) Conroy C., et al., 2022, arXiv e-prints, p. arXiv:2204.02989
  * Cui et al. (2012) Cui X.-Q., et al., 2012, RAA, 12, 1197
  * De Silva et al. (2015) De Silva G. M., et al., 2015, MNRAS, 449, 2604
  * Deason & Belokurov (2024) Deason A. J., Belokurov V., 2024, arXiv e-prints, p. arXiv:2402.12443
  * Di Matteo et al. (2019) Di Matteo P., Haywood M., Lehnert M. D., Katz D., Khoperskov S., Snaith O. N., Gómez A., Robichon N., 2019, A&A, 632, A4
  * Dillamore et al. (2024) Dillamore A. M., Belokurov V., Kravtsov A., Font A. S., 2024, MNRAS, 527, 7070
  * Dodd et al. (2023) Dodd E., Callingham T. M., Helmi A., Matsuno T., Ruiz-Lara T., Balbinot E., Lövdal S., 2023, A&A, 670, L2
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * El-Badry (2024) El-Badry K., 2024, arXiv e-prints, p. arXiv:2403.12146
  * El-Badry et al. (2024) El-Badry K., et al., 2024, arXiv e-prints, p. arXiv:2405.00089
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., et al., 2013, emcee: The MCMC Hammer, Astrophysics Source Code Library, record ascl:1303.002 (ascl:1303.002)
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A2
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2023a) Gaia Collaboration et al., 2023a, A&A, 674, A1
  * Gaia Collaboration et al. (2023b) Gaia Collaboration et al., 2023b, A&A, 674, A34
  * Gaia Collaboration et al. (2024) Gaia Collaboration et al., 2024, arXiv e-prints, p. arXiv:2404.10486
  * Gallart et al. (2019) Gallart C., Bernard E. J., Brook C. B., Ruiz-Lara T., Cassisi S., Hill V., Monelli M., 2019, Nature Astronomy, 3, 932
  * Hasselquist et al. (2021) Hasselquist S., et al., 2021, ApJ, 923, 172
  * Haywood et al. (2018) Haywood M., Di Matteo P., Lehnert M. D., Snaith O., Khoperskov S., Gómez A., 2018, ApJ, 863, 113
  * Helmi (2008) Helmi A., 2008, A&ARv, 15, 145
  * Helmi (2020) Helmi A., 2020, ARA&A, 58, 205
  * Helmi & White (1999) Helmi A., White S. D. M., 1999, MNRAS, 307, 495
  * Helmi et al. (1999) Helmi A., White S. D. M., de Zeeuw P. T., Zhao H., 1999, Nature, 402, 53
  * Helmi et al. (2018) Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85
  * Hurley & Tout (1998) Hurley J., Tout C. A., 1998, MNRAS, 300, 977
  * Hwang et al. (2022) Hwang H.-C., et al., 2022, MNRAS, 513, 754
  * Katz et al. (2023) Katz D., et al., 2023, A&A, 674, A5
  * Koppelman et al. (2019) Koppelman H. H., Helmi A., Massari D., Price-Whelan A. M., Starkenburg T. K., 2019, A&A, 631, L9
  * Lane et al. (2022) Lane J. M. M., Bovy J., Mackereth J. T., 2022, MNRAS, 510, 5119
  * Latham et al. (2002) Latham D. W., Stefanik R. P., Torres G., Davis R. J., Mazeh T., Carney B. W., Laird J. B., Morse J. A., 2002, AJ, 124, 1144
  * Majewski et al. (2017) Majewski S. R., et al., 2017, AJ, 154, 94
  * Malhan (2022) Malhan K., 2022, ApJ, 930, L9
  * Maoz et al. (2012) Maoz D., Badenes C., Bickerton S. J., 2012, ApJ, 751, 143
  * Matsuno et al. (2019) Matsuno T., Aoki W., Suda T., 2019, ApJ, 874, L35
  * Mazzola et al. (2020) Mazzola C. N., et al., 2020, MNRAS, 499, 1607
  * Moe et al. (2019) Moe M., Kratter K. M., Badenes C., 2019, ApJ, 875, 61
  * Mowlavi et al. (2023) Mowlavi N., et al., 2023, A&A, 674, A16
  * Myeong et al. (2018a) Myeong G. C., Evans N. W., Belokurov V., Amorisco N. C., Koposov S. E., 2018a, MNRAS, 475, 1537
  * Myeong et al. (2018b) Myeong G. C., Evans N. W., Belokurov V., Sanders J. L., Koposov S. E., 2018b, MNRAS, 478, 5449
  * Myeong et al. (2019) Myeong G. C., Vasiliev E., Iorio G., Evans N. W., Belokurov V., 2019, MNRAS, 488, 1235
  * Naidu et al. (2020) Naidu R. P., Conroy C., Bonaca A., Johnson B. D., Ting Y.-S., Caldwell N., Zaritsky D., Cargile P. A., 2020, ApJ, 901, 48
  * Nissen & Schuster (2010) Nissen P. E., Schuster W. J., 2010, A&A, 511, L10
  * Niu et al. (2021) Niu Z., Yuan H., Wang S., Liu J., 2021, ApJ, 922, 211
  * O’Hare et al. (2020) O’Hare C. A. J., Evans N. W., McCabe C., Myeong G., Belokurov V., 2020, Phys. Rev. D, 101, 023006
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Price-Whelan et al. (2020) Price-Whelan A. M., et al., 2020, ApJ, 895, 2
  * Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
  * Rix et al. (2022) Rix H.-W., et al., 2022, ApJ, 941, 45
  * Sartoretti et al. (2023) Sartoretti P., et al., 2023, A&A, 674, A6
  * Semenov et al. (2024) Semenov V. A., Conroy C., Chandra V., Hernquist L., Nelson D., 2024, ApJ, 962, 84
  * Taylor (2005) Taylor M. B., 2005, in Shopbell P., Britton M., Ebert R., eds, Astronomical Society of the Pacific Conference Series Vol. 347, Astronomical Data Analysis Software and Systems XIV. p. 29
  * Tian et al. (2018) Tian Z.-J., et al., 2018, Research in Astronomy and Astrophysics, 18, 052
  * Tokovinin (2004) Tokovinin A., 2004, in Allen C., Scarfe C., eds, Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 21, Revista Mexicana de Astronomia y Astrofisica Conference Series. pp 7–14
  * Tokovinin et al. (2006) Tokovinin A., Thomas S., Sterzik M., Udry S., 2006, A&A, 450, 681
  * Torres et al. (2010) Torres G., Andersen J., Giménez A., 2010, A&ARv, 18, 67
  * Zhang et al. (2023) Zhang H., Ardern-Arentsen A., Belokurov V., 2023, arXiv e-prints, p. arXiv:2311.09294