# Delays-induced Phase Transitions in Active Matter

Fatemeh Pakpour  Tamás Vicsek

###### Abstract

We consider the patterns of collective motion emerging when many aligning,
self-propelling units move in two dimensions while interacting through a
repulsive potential and are also subject to delays and random perturbations.
In this approach, delay plays the role analogous to reaction time so that a
given particle is influenced by the information about the velocity and the
position of the other particles in its vicinity with some time delay. To get
insight into the involved complex flows and the transitions between them we
use a simple model allowing – by fine-tuning of its few parameters – the
observation and analysis of behaviours that are less accessible by experiments
or analytic calculations and at the same time make the reproduction of
experimental results possible. We report for the first time about a transition
from an ordered, polarized collective motion to disorder as a function of the
increasing time delay. For a fixed intermediate value of the delay similar
transition (from order to disorder) is obtained as the repulsion radius is
increased. Our simulations show a transition from total polarization to two
kinds of states: fully disordered and a kind of state which is a mixture of
patches of fully disordered motion in the background of orderly moving other
particles. The transition occurs as the delay time is increased and is sharp,
indicating that the nature of this order-disorder transition is either of
first-order or is described by a sharply decreasing linear function. Our model
is a simplified version of a practical situation of quickly growing interest
because time delays are expected to play an increasingly important role when
the traffic of many, densely distributed autonomous drones will move around in
a quasi-two-dimensional air space.

††journal: Physica A

[inst1]organization=Science Department,addressline=Arak University of
Technology, city=Arak, postcode=3818146763, state=Markazi, country=Iran

[inst2]organization=Department of Biological Physics, Eötvös
University,addressline= Pázmány Péter sétány 1A, city=Budapest,
postcode=H-1117, country=Hungary

[inst3]organization=MTA-ELTE Statistical and Biological Physics Research Group
of HAS,addressline= Pázmány Péter sétány 1A, city=Budapest, postcode=H-1117,
country=Hungary

Almost all systems of collectively moving units are subject to some sort of
delays caused by finite time intervals needed for reacting to the positions
and velocities of the others. Still, a systematic study of the effect of
delays on the destruction of coherent flocking is lacking. Here we show that
if the members of a group of aligning particles travelling in two dimensions
have a finite size, i.e., repel each other, a large enough reaction time
results in either a complete disorder or in a new kind of motion state in
which patches of disordered and ordered states can coexist. Our findings are
expected to be relevant to situations such as the traffic of many units, e.g.,
drones in a confined area.

##  1 Introduction

Systems of many similar units exhibit collective behavioural patterns through
interactions among their units. In the case of collective motion, these
interactions can be looked at as a result of exchanging information about the
positions of the others and then accompanied by the corresponding reaction. In
nearly all systems the reaction to a particular set of information about the
positions of the others is not instantaneous due to several sources of non-
negligible reaction times resulting in delays. The mechanisms behind delayed
reactions are common in both living and technological systems. In this work,
we explore in which ways time delays can result in the destruction of coherent
behaviour during the motion of aligning self-propelled particles (SPP-s)
representing a typical manifestation of active matter.

The study of collective motion [1] of groups of biological and non-biological
agents involved also in the context of soft active matter [2] has been a
subject of quickly growing interest over the past two decades. Examples of
related studies include bacterial colonies [3], cell migration [4], swarms of
insects[5], schools of fish [6], flocks of birds [7], and crowds of people [8]
to nematic fluids [9], shaken metallic rods [10], and flying robots [11]. A
comprehensive understanding of the complex and nonlinear dynamics of a system
consisting of a large number of interacting entities in addition to
environmental impacts and important applications in biology [12], ecology
[13], ethology [6, 14], agriculture [5], economics [15], social science [16],
and cognitive science [17] can be functional in developing modern technologies
of artificial agents such as drones [18] and autonomous flying robots [19,
20]. The fascinating and extraordinary features of the collective motion in
diverse systems of interacting units over a wide range of system sizes such as
self-organization without any centralized control and leaderless coherent
motion [21, 22], hierarchical leader-follower relationships among the members
[23, 24], long-lasting giant number fluctuations [25, 26], turbulent flow and
vortex formation [27, 28, 29], and transitions to a variety of motions namely
marching, milling, rotating chains, and jamming [30, 16, 31, 32, 33, 34] have
been manifested through numerical and experimental investigations which can be
interpreted in the scope of equilibrium and non-equilibrium statistical
physics.

Various mathematical and computational models [21, 35, 36, 37, 38, 39, 40, 41,
42, 43] have been introduced to investigate the dynamics of self-propelled
particles. The simplest model accounting for the non-equilibrium phase
transition in flocking as a function of random perturbations was introduced by
Vicsek et.al. [36] which we shall refer to as SPP95. In the SPP95, the
collision avoidance effect is ignored and aligning with neighbouring particles
occurs instantly. However, in realistic systems, either natural or artificial,
in addition to the requirement of considering repulsive interactions to
prevent collisions, each agent needs some time to send, receive and respond to
the information signals from its neighbourhood. Thus, although the very simple
SPP95 model already displayed a non-trivial non-equilibrium phase transition
as a function of the level of random perturbations it was not designed to
account for the above-mentioned details of the interactions potentially
resulting in destroying the polarized motion which usually emerges if only
alignment is considered.

One of the very early studies considering the problem of delay time in the
collective behaviour of many-body systems appeared in the subject of vehicular
traffic. It was known that the observed density waves and instabilities in
vehicle traffic are not explained without considering the reaction time of
drivers [44, 45, 46]. Following that, the question of sensory temporal delay
provided a solution to the problem of gradient sensing mechanisms in bacterial
colonies. An experimental research showed that bacterial chemotaxis appears
using a time-sensing mechanism as if they have a kind of "memory device" to
compare the difference between environmental concentrations at successive
times [47]. The effect of time delay was also investigated in different many-
body physical systems such as physiological control systems [48], the ring
cavity system [49], cortical neurons [50], and lasers [51] which showed
oscillatory behaviours and chaos.

The emergence of instability in a two-dimensional system of self-propelled
agents interacting via a pair-wise linear attractive force and in the presence
of noise and a constant time delay was demonstrated in [52]. The study showed
that the delay induces a transition into a state which is an oscillating
aligned swarm that depends on the size of the interaction coupling
coefficient. In a work reporting on both a simulation of flocking birds with
two decentralized control algorithms and an experiment with real autonomous
flying robots, the realistic features of the flocking such as time delay,
inertial effects, sensor errors, and locality of communication besides a
friction-like interaction for alignment of the neighbouring particles, a
repulsive force, and noise considered [11]. The authors showed that although
time delay can reduce the stability of the flocking state the friction-like
term reduces the strength of this instability. They also demonstrated how time
delay leads to oscillations in the system and indicated that an additive
Gaussian noise term can decrease the instabilities caused by delay time.

In a one-dimensional numerical study on the directional switching of the self-
propelled agents with delayed interaction, it was shown that with increasing
delay time the average switching time of the particles increases [53]. The
authors examined three types of delay time: a particular time delay for all
particles, a random delay chosen from a normal distribution, and a state-
dependent form of delay time in which it increases with the level of disorder
in the vicinity of the particle. They showed the mean switching time for the
group with identical delay time is shorter than those with random delay time
and longer than the third group with heterogeneity. They also argued that the
time delay may facilitate the coherent motion of the latter group.

The effect of delay on the collective behaviour was also investigated
experimentally by using a small number of autonomous phototactic robots being
able to measure the intensity of a decaying light field emitted from the
others [54]. It was shown that sensorial delay - the time each particle needs
to change its velocity in response to the intensity of light received from
another particle - can be a control parameter for the type of collective
movement of particles from segregation to aggregation and clustering. To
confirm the scalability of the used mechanism, they simulated the behaviour of
a group of 100 robots and showed a positive delay results in forming clusters
while a negative delay can cause the robots to move away from each other.
Another numerical study on the role of time delay on the safety of an assembly
of up to 5 Unmanned Aerial Vehicle (UAVs) in a 3D space was conducted [55].
The research showed that communication delay has a major role during the
flock-building stage; longer delay time resulted in an extended time needed
for the flock to stabilize in the desired direction. The authors also showed a
very high delay time leads to an oscillatory behaviour which in turn increases
the possibility of collisions and decreases the safety of the group. They
argued that the number of agents building a flock is limited by the duration
of the communication delay.

The facilitation of coherent motion for a group of self-propelling agents at
short delay times was also demonstrated in a two-dimensional numerical
investigation based on the SPP95 and using a Voronoi tessellation tool for
more efficient measurement of clustering formation [56]. In contrast, the
study showed the order is destroyed when the delay time is around 25 times
longer. This phenomenon was also observed in another numerical work using the
SPP95-like model of active matter in three dimensions with periodic boundary
conditions [57]. The authors showed the delay time first facilitates the
ordered motion but with increasing time delay the characteristic oscillation
of delayed systems is observed. They observed an identical finite-size scaling
and correlation as natural swarms show. They also argued that the effect of
time delay in a system of SPPs is similar to inertia.

Although in practical situations the magnitude of the velocities of the SPPs
can vary over a broad range, in most of the related simulations the speed of
the particles was assumed to be a constant. A particularly detailed/complete
model [58] took into account numerous realistic factors playing an essential
role in the formation of flocks of birds. Their StarDisplay model and the
related simulations take into account the reaction times as well and are
capable of reproducing many realistic flocking scenarios. Perhaps the most
relevant example, in which the speed of the flocking units has to be
considered is a flock of drones moving in confined areas. Such a system was
investigated both numerically and experimentally in a study in which over 30
autonomous drones were used to explore the effects emerging due to the
criterion of avoidance (drones were not allowed to move too close to each
other) and the time delays resulting from the onboard calculations carried out
by the robots [18]. It was demonstrated that only a careful (evolutionary)
optimization of the parameters of the autopilot code could result in a
coherent and collisionless motion under such conditions.

In the present work, we focus on the statistical mechanics effects of taking
into account the two realistic conditions of delays and repulsion. Instead of
looking for reproducing observed patterns of flocking, we consider the general
question of a possible phase transition in the systems of SPPs as a result of
increasing delay times and repulsion radii.

##  2 Model, Parameters, and Definitions

In order to investigate the effect of time delays in systems in which they may
play a significant role, we consider $N$ particles moving continuously in an
$L\times L$ cell with periodic boundary conditions. At $t=0$, the particles
are randomly distributed within an area of a given size at positions
$\mathbf{x}_{i}(0)$. In addition, each particle is assigned a velocity
$\mathbf{v}_{i}(0)$ with the same magnitude, $v_{o}$, and randomly distributed
directions $\theta_{i}(0)$. The time interval between two updates of the
positions and directions is $\Delta t=1$. At each time step, the position of
the $i^{th}$ particle is updated according to:

Figure 1: Visualization of the model. $\mathit{(a)}$ is an explanatory panel,
demonstrating the way in which delay plus repulsion can result in a larger
angle between the departing particles as compared to the angle between their
trajectories before "collision". The interaction starts at $t=4$ and the delay
time for particle1 $\tau_{1}=3$, and for particle2 $\tau_{2}=4$. At $t=4$
particle1 remembers the information from particle2 at $t=1$ (the dotted
velocity vector), while particle2 remembers the information from particle1 at
$t=0$ (the dashed velocity vector). For each step, the green vector is the
alignment term, $\mathbf{v}^{align}$, the red one is the repulsion term,
$\mathbf{v}^{rep}$, the blue vector is the sum of two terms, and the black
vector is the final velocity of the particle after adding random noise to the
blue one. The repulsion term is in the direction of the distance vector
between two particles that was shown with the orange vector. $\mathit{(b)}$
shows the repulsion and alignment radii, $\mathit{(c)}$ illustrates the
repulsion force direction, $\mathit{(d)}$ displays the mathematical form of
the repulsion force. Supplementary movie 1 complements the above visualization
by illustrating the behaviour of the particles on a small scale compatible
with the range of the alignment radius.

$\mathbf{x}_{i}(t+1)=\mathbf{x}_{i}(t)+\mathbf{v}_{i}(t+1)\Delta t,$ (1)

where $\mathbf{x}_{i}(t)$ is the position at time $t$ and
$\mathbf{x}_{i}(t+1)$ and $\mathbf{v}_{i}(t+1)$ denote the position and the
velocity of the agent at time $t+1$. We use the term agent and particles
interchangeably and a collection of aligning agents represents active matter.
The velocity of each agent at time $t+1$ is calculated by considering three
terms: i) the effect of aligning with neighbouring particles located within a
certain alignment radius, $R$, ii) the repulsion from particles closer than a
certain repulsion radius, $2r_{rep}$, and a random noise $\eta_{i}(t)$. Each
particle has a delay time, $\tau_{i}$, it spends on processing information
about the positions and the velocities of the others, while it reacts without
delay regarding information about itself. The delay times in this simulation
were chosen from a Gaussian distribution [59] with mean $\mu_{\tau}$ and
standard deviation $\sigma_{\tau}$.

The alignment term for $i^{th}$ particle at time $t+1$ is calculated by taking
the average over the particle’s velocity at time $t$ and the velocity of
neighbouring particles located within radius $R$ at time $t-\tau_{i}$:

$\mathbf{v}^{align}_{i}(t+1)=\mathcal{N}^{v_{o}}\left(\frac{\sum_{j=1}^{N}a_{ij}(t-\tau_{i})\mathbf{v}_{j}(t-\tau_{i})+c_{align}\mathbf{v}_{i}(t)}{s_{i}(t-\tau_{i})+1}\right),$
(2)

where $\mathcal{N}^{v_{o}}()$ is a function which normalizes the magnitude of
its argument to $v_{o}$ and $s_{i}(t-\tau_{i})$ is the number of neighbours
inside radius $R$ at time $t-\tau_{i}$. $c_{align}$ is the coefficient of
alignment and $a_{ij}$ are the elements of the alignment matrix defined as
follows:

$\displaystyle a_{ij}=\begin{cases}1,\ \ i\neq j,\ \ i=1,...,N,\ \ j\in
S_{i}(t)\\\ 0,\ \ \text{otherwise,}\end{cases}$ (3)

where $S_{i}(t)$ is the set of particles inside the interaction radius $R$ and
for simplicity we assumed the simplest form of the contribution of the
velocity of a particle to the resulting direction (all $a_{ij}=1$ or $0$). The
alignment matrix is updated at each time step.

The repulsion term is a sum of pairwise repulsion forces depending on the
distance of particles according to:

$\begin{split}\mathbf{v}^{rep}_{i}(t+1)=c_{rep}\sum_{k}g(r_{i}+r_{k}-\lvert\mathbf{x}_{ik}(t-\tau_{i})\rvert)\frac{\mathbf{x}_{ik}(t-\tau_{i})}{\lvert\mathbf{x}_{ik}(t-\tau_{i})\rvert},\end{split}$
(4)

where $c_{rep}$ is the repulsion coefficient, $r_{i}$ is the repulsion radius
of particle $i$ and $r_{k}$ is the repulsion radius of the neighbouring
particle $k$ and $k$ runs over the particles located within a circle centred
at particle $i$ with the radius of $r_{i}+r_{k}$. We assume that both the
$r_{i}$ and $r_{k}$ are equal to $r_{rep}$. $\mathbf{x}_{ik}(t-\tau)$ is
$\mathbf{x}_{i}-\mathbf{x}_{k}$ at time $t-\tau$. The form of $g(x)$ can be
chosen quite freely, however, for the sake of simplicity we shall assume a
simple hard-core type dependence, meaning that repulsion acts only within a
given repulsion radius. Thus we choose the g(x) function as follows:

$\displaystyle g(x)=\begin{cases}0,\ \ \ |\mathbf{x}_{i}-\mathbf{x}_{k}|\geq
r_{i}+r_{k}\\\ x,\ \ \ \text{otherwise.}\end{cases}$ (5)

Finally, the velocity vector at time $(t+1)$ is obtained as follows:

$\mathbf{v}_{i}(t+1)=\mathcal{R}^{\eta_{i}(t+1)}\left[\mathcal{N}^{v_{max}}\left(\mathbf{v}^{align}_{i}(t+1)+\mathbf{v}^{rep}_{i}(t+1)\right)\right],$
(6)

where again $\mathcal{N}^{v_{max}}()$ normalizes its argument to a maximum
magnitude, $v_{max}$, and the function $\mathcal{R}^{\eta_{i}}()$ is a
rotation operator that changes the angle of its argument with $x$ axis by the
amount $\eta_{i}$ where $\eta_{i}$ is a random number chosen from a uniform
distribution from $[-\eta/2,\eta/2]$ interval. A visualization of the model is
shown in Figure 1 and a list of parameters used in our model is given in Table
1.

Table 1: The list of parameters together with their explored ranges used in
the simulations of the model of flocking with repulsion and delay. List of
Parameters
---
Parameter |  Definition |  Range
$N$ |  Number of particles |  $2000-16000$
$\rho$ |  Number density of the particles |  $1.0-3.0$
$v_{o}$ |  Initial velocity |  $0.01-0.03$
$R$ |  Alignment radius |  $1.0-1.5$
$c_{align}$ |  Alignment coefficient |  $1.0-3.0$
$r_{rep}$ |  Repulsion radius |  $0.2-0.45$
$c_{rep}$ |  Repulsion coefficient |  $0.1-1.0$
$v_{max}$ |  Maximum value of the velocity |  $3v_{o}$
$\eta$ |  Maximum noise angle |  $0-\pi$
$\mu_{\tau}$ |  Mean of time delay distribution |  $0-40$
$\sigma_{\tau}$ |  Standard deviation of the delay distribution |  $\mu_{\tau}/5,\mu_{\tau}/10$
$T$ |  The number of iterations in any run |  $6000-14000$

To describe the degree of global alignment in the system the normalized order
parameter is calculated using the following equation at each time step [1]:

$\noindent\Phi(t)=\frac{1}{N}\bigg{|}\sum_{j=1}^{N}\frac{\mathbf{v}_{j}(t)}{\lvert\mathbf{v}_{j}(t)\rvert}\bigg{|}.\\\
$ (7)

The value of $\Phi$ varies from nearly zero in a situation of disordered
motion and tends to $1.0$ for ordered/polarized motion.

##  3 Simulation and Results

The simulations were carried out using a GPU-accelerated code based on CUDA
Python making use of parallel processing. Each point in all plots was obtained
from averaging over 5 separate runs. The dynamics and the degree of
coordination in the system were tracked using equations $1-7$ through several
thousands of time steps depending on the system size. To account for the
effect of time delay we first allowed the particles to move freely from the
initial state for a short interval, normally a few time steps longer than the
maximum time delay, and then the new positions and velocities were calculated.
As mentioned in the previous section the time delays were drawn from a normal
distribution so that they were different for different particles. We explored
a wide range of parameters to find the most relevant patterns of the motion of
particles. The behaviour of the system and the phase transitions occurring
while changing the various parameters such as the mean time delay, noise, and
repulsion radius were studied.

Figure 2: A graph of normalized order parameter vs mean time delay for various
$N$. $\rho=1.0$, $R=1.0$, $c_{align}=c_{rep}=1.0$, $r^{rep}=0.42$,
$v_{o}=0.01$, and $\eta=0.31$. The phase transition curve is independent of
the system sizes.

Before deciding which value for N would be optimal to choose for the majority
of runs, we looked at the dependency of the order parameter versus delay
curves for a set of parameters and for a variety of N values. The results for
3 different numbers of particles, $N=2048$, $4096$, and $8192$ are displayed
in Figure 2 showing the transition from the ordered to the disordered states
as a function of the mean time delay. From this figure we conclude that there
is only a very weak dependence of the behaviour on the particle number $N$ for
$N>2048$ so - unless otherwise stated - the number of particles in our
simulations is 2048. Throughout this paper we report on behaviour in the
stationary state of the system which is obtained after the order parameter
values saturate (start to fluctuate around the given value)

Figure 3: Order parameter as a function of mean time delay for various
repulsion radii ranging from $0.30$ to $0.44$. The parameters used are:
$\rho=1.0$, $R=1.0$, $c_{align}=c_{rep}=1.0$, $v_{o}=0.01$, and $\eta=0.31$.
For $r_{rep}\geq 0.38$ there is a phase transition from an ordered phase to a
fully disordered phase while for $r_{rep}<0.38$ the system eventually reaches
a state which is a mixture of groups of fully disordered agents in the
background of orderly moving other particles.

To illustrate the effect of delays on the level of coherence of the motion of
particles the order parameters as a function of $\mu_{\tau}$ were plotted for
various repulsion radii in Figure 3 for moderate level of noise. This figure
demonstrates a behaviour that is very specific to our system. In spite of the
simplicity of our model it displays an atypical sequence of transitions
involving two distinct scenarios and a collective motion phase which has not
yet been observed in such simple models. Figure 3 indicates these kinds of
transitions: i) For relatively large repulsion radii (i.e., when two times
$r_{rep}$ is close to the alignment radius the system goes through a change
from a completely ordered to a fully disordered state as the average delay
time is increased. However, ii) for smaller repulsion radii the system -
although is also ordered for small delays - does not become fully disordered
even if the average delay is substantially increased. In this latter case, the
system is in a mixed state of nearly perfectly ordered and fully disordered
patches.

The transition i) to the fully disordered state is accompanied by a sharp
decrease in the order parameter indicating a phase change that is similar to a
first-order transition, but the decrease is less abrupt and leaves open the
possibility of a continuous transition with an exponent being $-1$. The
transition to the new, mixed state does not seem to have an analogy with
ordinary phase transitions, for example, the value of the order parameter in
the large delay limit is dependent on the repulsion radius and thus, does not
allow the definition of a single-parameter dependent order parameter
displaying values from $1$ to $0$. In Figure 4 four example screenshots of the
above-mentioned three states are shown. The videos of the corresponding
simulations can be found in Supplementary movies $(2-5)$. Two other videos for
$N=16384$ particles with similar parameters i.e. $\rho=1.0$, $R=1.0$,
$c_{align}=c_{rep}=1.0$, $v_{o}=0.01$, and $\eta=0.31$, and for
$(r_{rep}=0.38$, $\mu_{\tau}=5)$, and $(r_{rep}=0.32$, $\mu_{\tau}=20)$ is
found in Supplementary movies $6$ and $7$ respectively.

Figure 4: Four screenshots of different states of the system for $N=4096$. In
all shots $\rho=1.0$, $R=1.0$, $c_{align}=c_{rep}=1.0$, $v_{o}=0.01$, and
$\eta=0.31$. The instantaneous velocities are indicated by blue vectors of
length proportional to their magnitude. $\mathit{(a)}$ The fully ordered state
for small delay time, $\mu_{\tau}=5.0$ and $r_{rep}=0.4$ after $5000$ time
steps. $\mathit{(b)}$ Mixed state for $\mu_{\tau}=18$ and $r_{rep}=0.3$ after
$6000$ time steps. $\mathit{(c)}$ Mixed state for $\mu_{\tau}=35$ and
$r_{rep}=0.3$ after $6000$ time steps.$\mathit{(d)}$ The fully disordered
state for $\mu_{\tau}=15$ and $r_{rep}=0.4$ after $6000$ time steps. Figure 5:
Order parameter as a function of $r_{rep}$ for various mean delay time
$\mu_{\tau}$. The parameters used are: $\rho=1.0$, $R=1.0$,
$c_{align}=c_{rep}=1.0$, $v_{o}=0.01$, and $\eta=0.31$. Note that for large
delays and relatively small repulsion radii the system stays disordered while
it displays ordered collective motion for short average delay times.

The two main control parameters we consider in this work are the average delay
and the repulsion radius. (Although the model has several more parameters, see
Table 1, we mostly used for the other parameters either a simple trivial
value, e.g., 1, or some intermediate value. e.g., for the noise, which is far
from the given parameter resulting in complex behaviour.) In Figure 3 we
showed the behaviour of the order parameter as a function of delay, next we
also looked at how the order parameter changes if for a set of fixed delays
the radius of repulsion $r_{rep}$ is increased. Figure 5 demonstrates that, as
expected, increasing the repulsion radius results in a completely disordered
state when $r_{rep}$ approaches 0.5 being half of the alignment radius (2
times $r_{rep}$ getting close to $R=1$). Increased disorder follows from two
effects: first, the frequency of "collisions" becomes larger as the effective
density grows, and second, the area within which it is only the alignment term
which acts (within the annulus between the two circles of radii $r_{rep}$ and
$R=1$).

Figure 6: The phase diagram for fully ordered, fully disordered, and mixed
states with respect to the mean time delay $\mu_{\tau}$ and the repulsion
radius $r_{rep}$ for the two-dimensional model of flocking with delay and
repulsion. The other parameters used in the related simulations were; noise:
$\eta=0.31$, density: $\rho=1.0$, alignment radius $R=1.0$, repulsion
coefficient: $c_{rep}=c_{align}=1.0$, and default velocity $v_{o}=0.01$.

An important motivation for obtaining the curves in Figure 5 is that together
with those displayed in Figure 3 we can construct an approximate phase diagram
shown in Figure 6 depicting the regions of the two control parameters
$\mu_{\tau}$ and $r_{rep}$ for which a given collective motion state (Fully
ordered, Mixed and Fully-disordered) takes place. What Figure 6 demonstrates
is that the regions of the two main control parameters for which all three
phases can be observed are quite limited and the extension of the region
corresponding to the mixed phase dominates the behaviour.

##  4 Conclusions

Although macroscopic self-propelling (active) units of a system are typically
subject to interactions consisting of both a repulsion term and the effect of
time delays (reaction times) our simulation is the first study systematically
addressing the question of how these two factors affect the collective motion
of the units. In the limit of small delays and small repulsion, we find, as
expected, fully ordered, coherent motion in the system. Similarly, although
less trivially, we obtain a completely disordered motion state for large
delays and large repulsion radii. However, in addition to these, we have also
found a new kind of collective motion pattern: disordered patches embedded
into a globally ordered flow.

We argue that this complex behaviour can be understood by considering
competing effects on the behavior: As it is also illustrated by our 1 (a),
collisions in the presence of delays can result in trajectories of the
particles diverging by an angle being larger than the one along which they
approached each other. This mechanism results in an increased disorder and,
through a local increase of the magnitude of the "scattered" particles results
in an increased pushing of the neighbours deeper into their repulsion area
which, in turn, results in further sharp turns of the trajectories. This is an
instability which would result in completely disorderly motion without any
stabilizing effect. However, delay times - somewhat paradoxically, but in a
formerly shown manner (see [56, 57]) - may result in making the collective
motion more coherent. Perhaps the most straightforward interpretation of this
phenomenon is to look at very large delays as a limit in which the system
behaves quite like its mean-field version without delay (only already distant
particles interact). The interplay of the instability due to abrupt collisions
and the stabilizing effect of long delays results in a mixed phase containing
parts in which either the instability (resulting in disorder) or the
stabilizing long delay dominate the local behaviour. It should be noted that
an intermediate value of the order parameter in the present case is not
analogous to a partially ordered system of purely aligning particles because
in the latter case the level of disorder - due to noise - is distributed
homogeneously through the system.

##  5 Acknowledgements

The authors are grateful to Fahimeh Najafi and Máté Nagy for helpful
discussions. This research was partially supported by the following grants:
(Hungarian) National Research, Development and Innovation Office grants No:
K128780 and SNN 139598. It has also received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No 955576.

## References

  * [1] T. Vicsek, A. Zafeiris, Collective motion, Phys. Rep. 517 (2012) 71–140.
  * [2] M. Marchetti, J. Joanny, S. Ramaswamy, T. Liverpool, J. Prost, M. Rao, R. Simha, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85 (2013) 1143\.
  * [3] H. P. Zhang, A. Be’er, E. L. Florin, H. L. Swinney, Collective motion and density fluctuations in bacterial colonies, Proc. Natl. Acad. Sci. U. S. A. 107 (2010) 13626–13630.
  * [4] E. Méhes, T. Vicsek, Collective motion of cells: from experiments to models, Integr. Biol. 6 (2014) 831–854.
  * [5] J. Buhl, I. Sumpter, D. J.and Couzin, J. Hale, E. Despland, E. R. Miller, S. Simpson, From disorder to order in marching locusts, Science 312 (2006) 1402–1406.
  * [6] U. Lopez, J. Gautrais, I. D. Couzin, G. Theraulaz, From behaviuoral analyses to models of collective motion in fish schools, Interface Focus 2 (2012) 693–707.
  * [7] I. Bajec, F. Heppner, Organized flight in birds, Animal Behaviour 78 (2009) 777–789.
  * [8] M. Moussaid, N. Perozo, S. Garnier, D. Helbing, G. Theraulaz, The walking behavior of pedestrian social groups and its impact on crowd dynamics, PloS one 5 (2010) e10047.
  * [9] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, F. Sagués, Active nematics, Nat. Commun. 9 (2018) 3246.
  * [10] A. Kudrolli, Concentration dependent diffusion of self-propelled rods, Phys. Rev. Lett. 104 (2010) 088001.
  * [11] C. Virágh, G. Vásárhelyi, N. Tarcai, T. Szörényi, G. Somorjai, T. Nepusz, T. Vicsek, Flocking algorithm for autonomous flying robots, Bioinspiration Biomimetics 9 (2014) 025012.
  * [12] A. Deutsch, G. Theraulaz, T. Vicsek, Collective motion in biological systems, Interface Focus 2 (2012) 689–692.
  * [13] P. A. Westley, A. M. Berdahl, C. J. Torney, D. Biro, Collective movement in ecology: from emerging technologies to conservation and management, Philos. Trans. R. Soc. B 373 (2018) 20170004.
  * [14] D. J. T. Sumpter, Collective animal behavior, Princeton University Press (2010).
  * [15] R. Cont, J. Bouchaud, Herd behavior and aggregate fluctuations in financial markets, Macroecon. Dyn. 4 (2000) 170–196.
  * [16] D. Helbing, I. Farkas, T. Vicsek, Simulating dynamical features of escape panic, Nature 407 (2000) 487–490.
  * [17] I. D. Couzin, Collective cognition in animal groups, Trends in Cognitive Sciences 13 (2009) 36–43.
  * [18] G. Vásárhelyi, C. Virágh, G. Somorjai, T. Nepusz, A. E. Eiben, T.
Vicsek, Optimized flocking of autonomous drones in confined environments, Sci.
Robot. 3 (2018) 3536.

  * [19] M. Brambilla, E. Ferrante, M. Birattari, M. Dorigo, Swarm robotics: A review from the swarm engineering perspective, Swarm Intell. 7 (2013) 1–41.
  * [20] A. Turgut, H. Çelikkanat, F. Gökçe, E. Şahin, Self-organized flocking in mobile robot swarms, Swarm Intell. 2 (2008) 97–120.
  * [21] I. Aoki, A simulation study on the schooling mechanism in fish, Bull. Japan. Soc. Sci. Fish 48 (1982) 1081–1088.
  * [22] G. Grégoire, H. Chaté, Y. Tu, Moving and staying together without a leader, Phys. D: Nonlinear Phenom. 181 (2003) 157–170.
  * [23] M. Nagy, Z. Ákos, D. Biro, T. Vicsek, Hierarchical group dynamics in pigeon flocks, Nature 464 (2010) 890–893.
  * [24] I. D. Couzin, J. Krause, N. R. Franks, S. A. Levin, Effective leadership and decision-making in animal groups on the move, Nature 433 (2005) 513–516.
  * [25] H. Chaté, F. Ginelli, R. Montagne, Simple model for active nematics: Quasi-long-range order and giant fluctuations, Phys. Rev. Lett. 96 (2006) 180602\.
  * [26] V. Narayan, S. Ramaswamy, N. Menon, Long-lived giant number fluctuations in a swarming granular nematic, Science 317 (2007) 105–108.
  * [27] H. Cisneros, R. Cortez, C. Dombrowski, R. E. Goldstein, J. O. Kessler, Fluid dynamics of self-propelled microorganisms, from individuals to concentrated populations, Exp. Fluids 43 (2007) 737–753.
  * [28] A. Sokolov, R. Goldstein, F. Feldchtein, I. Aranson, Enhanced mixing and spatial instability in concentrated bacterial suspensions, Phys. Rev. E 80 (2009) 031903.
  * [29] C. W. Wolgemuth, Collective swimming and the dynamics of bacterial turbulence, Biophysical Journal 95 (2008) 1564–1574.
  * [30] A. Kudrolli, G. Lumay, D. Volfson, L. S. Tsimring, Swarming and swirling in self-propelled polar granular rods, Phys. Rev. Lett. 100 (2008) 058001.
  * [31] T. Ihle, Kinetic theory of flocking: Derivation of hydrodynamic equations, Phys. Rev. E 83 (2011) 030901.
  * [32] S. Mishra, A. Baskaran, M. C. Marchetti, Fluctuations and pattern formation in self-propelled particles, Phys. Rev. E 81 (2010) 061916.
  * [33] N. Tarcai, C. Virágh, D. Ábel, M. Nagy, P. Várkonyi, G. Vásárhelyi, T. Vicsek, Patterns, transitions and the role of leaders in the collective dynamics of a simple robotic flock, J. Stat. Mech.: Theory Exp. 4 (2011) P04010.
  * [34] D. Helbing, Traffic and related self-driven many-particle systems, Rev. Mod. Phys. 73 (2009) 1067.
  * [35] C. Reynolds, Flocks, herds, and schools: A distributed behavioral model, In Computer Graphics 21 (1987) 25–34.
  * [36] T. Vicsek, A. Czirok, E. Ben-Jacob, I. Cohen, O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75 (1955) 1226\.
  * [37] J. Toner, Y. Tu, Long range order in a two-dimensional dynamical xy model, Phys. Rev. Lett. 75 (1995) 4326.
  * [38] N. Shimoyama, K. Sugawara, T. Mizuguchi, Y. Hayakawa, M. Sano, Collective motion in a system of motile elements, Phys. Rev. Lett. 76 (1996) 3870.
  * [39] I. Couzin, J. Krause, R. James, G. Ruxton, N. Franks, Collective memory and spatial sorting in animal groups, J. Theor. Biol. 218 (2002) 1–11.
  * [40] J. Toner, Y. Tu, S. Ramaswamy, Hydrodynamics and phases of flocks, Ann. Phys. 318 (2005) 170–244.
  * [41] F. Cucker, S. Smale, Flocks, herds, and schools: A distributed behavioral model, IEEE Trans. Automat. 52 (2007) 852–862.
  * [42] A. Derzsi, G. Szölösi, T. Vicsek, Most minimal spp model, URL http://hal. elte. hu/ vicsek/SPP-minimal (2009).
  * [43] A. Cavagna, I. Giardina, T. Grigera, The physics of flocking: Correlation as a compass from experiments to theory, Phys. Rep. 728 (2018) 1–62.
  * [44] R. E. Chandler, R. Herman, E. W. Montroll, Traffic dynamics: studies in car following, Oper. Res. 6 (1958) 165–184.
  * [45] E. I. J. I. Kometani, T. S. U. N. A. Sasaki, On the stability of traffic flow (report-i), J. Oper. Res. Soc. Jpn. 2 (1958) 11–26.
  * [46] R. Herman, E. W. Montroll, R. B. Potts, R. W. Rothery, Traffic dynamics: analysis of stability in car following, Oper. Res. 7 (1959) 86–106.
  * [47] R. M. Macnab, D. E. Koshland Jr, The gradient-sensing mechanism in bacterial chemotaxis, Proc. Natl. Acad. Sci. U.S.A. 69 (1972) 2509–2512.
  * [48] M. C. Mackey, L. Glass, Oscillation and chaos in physiological control systems, Science 197 (1977) 287–289.
  * [49] K. Ikeda, Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system, Opt. commun. 30 (1979) 257–261.
  * [50] A. S. Landsman, I. B. Schwartz, Synchronized dynamics of cortical neurons with time-delay feedback, Nonlinear Biomed. Phys. 1 (2007) 1–9.
  * [51] T. W. Carr, I. B. Schwartz, M. Y. Kim, R. Roy, Delayed-mutual coupling dynamics of lasers: scaling laws and resonances, SIAM J. Appl. Math. 5 (2006) 699–725.
  * [52] E. Forgoston, I. B. Schwartz, Delay-induced instabilities in self-propelling swarms, Phys. Rev. E 77 (2008) 035203.
  * [53] Y. Sun, W. Lin, R. Erban, Time delay can facilitate coherence in self-driven interacting-particle systems, Phys. Rev. E 90 (2014) 062708.
  * [54] M. Mijalkov, A. McDaniel, J. Wehr, G. Volpe, Engineering sensorial delay to control phototaxis and emergent collective behaviors, Phys. Rev. X 6 (2016) 011008\.
  * [55] V. Casas, A. Mitschele-Thiel, On the impact of communication delays on uavs flocking behavior, IEEE Wireless Communications and Networking Conference Workshops (WCNCW) (2018) 67–72.
  * [56] R. Piwowarczyk, T. Selin, M.and Ihle, G. Volpe, Influence of sensorial delay on clustering and swarming, Phys. Rev. E 100 (2019) 012607.
  * [57] V. Holubec, D. Geiss, S. Loos, K. Kroy, F. Cichos, Finite-size scaling at the edge of disorder in a time-delay Vicsek model, Phys. Rev. Lett. 127 (2021) 258001\.
  * [58] C. Hemelrijk, L. van Zuidam, H. Hildenbrandt, What underlies waves of agitation in starling flocks, Behav. Ecol. Sociobiol. 69 (2015) 755–764.
  * [59] M. Islam, I. Faruque, Insect visuomotor delay adjustments in group flight support swarm cohesion, Sci. Rep. 13 (2023) 6407.


# Interstellar Bow Shocks around Fast Stars
Passing through the Local Interstellar Medium

J. Michael Shull  Department of Astrophysical and Planetary Sciences
University of Colorado, Boulder, CO 80309  S. R. Kulkarni  Owens Valley Radio
Observatory, 249-17
California Institute of Technology, Pasadena, CA 91125
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Bow-shocks are produced in the local interstellar medium by the passage of
fast stars from the Galactic thin-disk and thick-disk populations with
velocities $V_{*}=$ 40–80 km s-1. Stellar transits of local H I clouds occur
every 3500–7000 yr on average and last between $10^{4}$ and $10^{5}$ yr. There
could be 10–20 active bow shocks around low-mass stars inside clouds within 15
pc of the Sun. At local cloud distances of 3–10 pc, their turbulent wakes have
transverse radial extents $R_{\rm wake}\approx$ 100–300 AU , angular sizes
10–100 arcsec, and Ly$\alpha$ surface brightnesses of 2-8 Rayleighs in gas
with total hydrogen density $n_{\rm H}\approx 0.1~{}{\rm cm}^{-3}$ and
$V_{*}=$ 40-80 km s-1. These transit wakes may cover an area fraction
$f_{A}\approx(R_{\rm wake}/R_{\rm cl})\approx 10^{-3}$ of local H I clouds and
be detectable in IR (dust), UV (Ly$\alpha$, two-photon), or non-thermal radio
emission. Turbulent heating in these wakes could produce the observed elevated
rotational populations of H2 ($J\geq 2$) and influence the endothermic
formation of CH+ in diffuse interstellar gas at $T>10^{3}$ K.

##  1 Introduction

Recent observations have advanced our knowledge of gas clouds in the
interstellar medium (ISM) and the velocity dispersions of stars in the thin-
disk and thick-disk populations. This is an opportune moment for a new look at
the frequency and physical effects of bow shocks produced by stars passing
through low-density H I clouds near the Sun. Most of these transits will be by
low-mass stars, although more distant examples of bow shocks around runaway
O-type stars have been found in infrared emission (Van Buren & McCray 1988;
Peri et al. 2015; Kobulnicki et al. 2016). The Herschel far-infrared survey
(Cox et al. 2012) identified numerous bow shocks around AGB stars and red
supergiants, including Betelgeuse (Decin et al. 2012). Theoretical studies of
intermediate-velocity shocks (Kulkarni & Shull 2023) suggest that the far-UV
background radiation may include contributions from Ly$\alpha$ and two-photon
continuum emission from the metastable (2s) state of H I. To excite the 2s and
2p levels (10.2 eV) requires shock velocities $V_{s}>40$ km s-1 with post-
shock temperatures $T_{s}\approx(3\mu V_{s}^{2}/16k)\approx(36,300~{}{\rm
K})V_{40}^{2}$, where $V_{40}=V_{s}/40~{}{\rm km~{}s}^{-1}$. Magnetic fields
will reduce this post-shock temperature.

In addition to their astronomical importance, fast stars passing through local
clouds contribute to the thermal energy budget of interstellar gas in their
transit wakes and extended regions. Observations of interstellar H2 in
rotationally excited states with angular momentum quantum number $J\geq 2$
suggest additional excitation mechanisms besides UV-pumping (Gry et al. 2002;
Nehmé et al. 2008; Ingalls et al. 2011; Shull, Danforth, & Anderson 2021).
Elevated rotational populations could be produced by collisional excitation in
warm molecular gas ($T>500$ K) heated by MHD turbulent dissipation (Moseley et
al. 2021). The observed correlation of CH+ with rotationally-excited H2
(Jensen et al. 2010) has led to similar production schemes for CH+ involving
warm H2 (Falgarone & Puget 1995; Myers et al. 2015). Local clouds of diffuse
interstellar gas exhibit significant turbulent line widths of 2–3 km s-1
(Linsky et al. 2022) derived by comparing the widths of absorption lines from
atomic species of different masses. They found that the LISM is inhomogeneous
on length scales less than 5100 AU. The inferred temperature distribution
appears Gaussian, with some sight lines reaching $T\approx 12,000$ K. Table 1
lists a number of local “fast stars”, together with their velocities and
properties of the LISM (gas temperatures and turbulent velocities).

Strong shocks will be produced by stars moving at velocities (40–70 km s-1)
through the warm neutral medium of local clouds, in excess of the sound speed
$c_{s}$ (or magnetosonic wave speed $v_{m}$). For magnetized gas, magnetosonic
waves are characterized by a combination of adiabatic sound speed,
$c_{s}=(5kT/3\mu)^{1/2}=(9.30$ km s${}^{-1})T_{8000}^{1/2}$, and Alfvén
velocity, $v_{A}=(5.53~{}{\rm km~{}s}^{-1})(B/3\,{\mu G})n_{\rm H}^{-1/2}$.
Here, we scaled to gas temperature $T\approx(8000~{}{\rm K})T_{8000}$, total
hydrogen density $n_{\rm H}$ (in cm-3), mean particle mass $\mu=1.27m_{H}$ for
$y\equiv$ He/H = 0.1, and magnetic field strength $B$. Absorption-line studies
(Frisch et al. 2011; Linsky et al. 2022) of clouds in the local interstellar
medium (LISM) suggest $n_{\rm H}\approx 0.1$ cm-3 and $B\approx$ 3–4 $\mu$G
(McComas et al. 2012; Zank et al. 2013). For wave propagation perpendicular to
the magnetic field, the fast-mode magnetosonic wave speed,
$v_{m}=(v_{A}^{2}+c_{s}^{2})^{1/2}$. Considering spatial variations in LISM
parameters ($n_{H}$, $T$, $B$), magnetosonic wave speeds in the LISM could
range from 20 to 30 km s-1, with further uncertainties arising from
propagation angle and ionization fraction111Ionization fractions are important
at the shock interface, owing to charge exchange reactions with neutrals.
Scherer & Fichtner (2014) also pointed out the importance of including both
He+ and H+ in determining the speed of Alfvén waves. For wave propagation
(vector ${\bf k}$ at angle $\phi$ to the magnetic field, the fast-mode
(Kulsrud 2005) dispersion relation
$(\omega/k)^{2}=(v_{A}^{2}+c_{s}^{2})/2+\frac{1}{2}[(v_{A}^{2}-c_{s}^{2})^{2}+4v_{A}^{2}c_{s}^{2}\sin^{2}\phi]^{1/2}$
defines the phase and group velocities of the wave packets. Near the front of
the bow shock, we expect that $\phi\approx 90^{\circ}$ (${\bf k}\perp{\bf
B}$). Away from the forward region, the fast-mode wave speed could be somewhat
smaller than $(v_{A}^{2}+c_{s}^{2})^{1/2}$..

Observations of the kinematics of the Milky Way disk suggest that populations
of older stars in the Galactic disk have supersonic and super-Alfvénic
velocities relative to the LISM. The large stellar velocity dispersions
measured in the Galactic thin disk (Nordstrom et al. 2004; Anguiano et al.
2020) are likely the result of dynamical heating processes that increase the
stellar velocities. Massive Galactic objects (spiral arms, giant molecular
clouds, bars, infalling dwarf satellites) perturb the motions of the newborn
stars, causing their random velocities to increase with time. A recent study
(Viera et al. 2022) using Gaia-DR3 proper motions found two kinematically
distinct Gaussian components: a thin disk with velocity dispersions
($\sigma_{R}$, $\sigma_{\phi}$, $\sigma_{Z}$) = 30, 21, and 11 km s-1 and a
thick disk with dispersions of 49, 35, and 22 km s-1. These correspond to 2D
(in-plane) dispersions of $\sigma_{\rm 2D}=37$ km s-1 (thin disk) and 60 km
s-1 (thick disk). The vertical (out-of-plane) dispersions (11 and 22 km s-1)
and total 3D dispersions (38 to 64 km s-1) also differ between the components.
These velocities are considerably larger than those of stars forming inside
dense molecular gas, seen in stellar associations and star clusters with
radial velocity dispersions $\sigma_{\rm 1D}=$1–3 km s-1. Proper motion
measurements with Gaia (Kuhn et al. 2019; Shull, Darling, & Danforth 2021)
show that young star clusters are expanding at velocities $\sim 0.5$ km s-1
and will disperse in several Myr after the gas is expelled.

##  2 Stellar Bow Shocks in the Local ISM

###  2.1 Star-Cloud Encounters

Because of the large velocity dispersions for stars in Galactic disk
populations, H I clouds in the LISM will have frequent encounters with stars
at velocities 40–70 km s-1. Absorption-line studies along 157 stellar sight
lines (Redfield & Linsky 2008; Linsky et al. 2019) identified numerous warm H
I clouds within 15 pc of the Sun. The nearest structures within $d\sim 4$ pc
are labeled the Local Interstellar Cloud (LIC), the G Cloud, and the Blue
Cloud. Redfield & Linsky (2008) estimated that between 6 and 19% of the LISM
volume is filled with warm H I gas. An updated analysis (Linsky et al. 2022)
found that clouds within 4 pc are more tightly packed, with mean hydrogen
density $\langle n_{\rm HI}\rangle\approx 0.1$ cm-3 and typical volume factors
$f_{V}\approx 0.38-62$. For local clouds beyond 4 pc, they find $f_{V}<0.4$.

We now estimate the rate at which passing stars encounter warm H I gas clouds
in the LISM from cloud cross sections, stellar space densities, and velocity
dispersions. If the region within distance $R_{\rm LISM}$ of the Sun is filled
by $N_{\rm cl}$ spherical clouds of radius $R_{\rm cl}$, the volume filling
factor is $f_{V}=N_{\rm cl}(R_{\rm cl}/R_{\rm LISM})^{3}$, and the area
covering factor is $f_{A}=N_{\rm cl}(R_{\rm cl}/R_{\rm LISM})^{2}$. These
factors are related by $f_{A}=f_{V}(R_{\rm LISM}/R_{\rm cl})$ for spheres.
More generally, $f_{A}$ and $f_{V}$ depend on cloud geometry, although the
non-spherical or filamentary shapes of some local clouds do not change these
estimates substantially. For 15 clouds within 10 pc, the above approximation
suggests mean cloud radii, $R_{\rm cl}=(10~{}{\rm pc})(f_{V}/N_{\rm
cl})^{1/3}\approx$ 2–3 pc and area covering factors $f_{A}\approx$ 0.7–1.0 for
$f_{V}=$ 0.15–0.30.

The star-cloud encounter rate depends on the local density of stars, both
those currently seen within the local clouds and the stars in the disk
populations that transit the LISM at high velocity. The recent Gaia-DR3 survey
of the local volume within 10 pc (Reylé et al. 2022) found 541 stars in 336
systems, of which 286 were luminous stars (we do not include the 85 brown
dwarfs in the survey). This corresponds to a mean density of $0.065~{}{\rm
pc}^{-3}$, similar to that quoted in Table 4.5 of Mihalas & Binney (1981). The
normalized densities of stars in the thin-disk and thick-disk populations
should sum to similar values, although there is some disagreement as to their
relative local densities (Anguiano et al. 2020; Everall et al. 2022; Viera et
al. 2023). Not all of these stars will have velocities sufficient to produce
the type of shocks discussed here. Only a fraction will have the high
velocities needed to exceed magnetosonic speeds. The radial extent of the
astrospheres, bow shocks, and turbulent wakes depend on the presence of a
moderately strong stellar wind (Herbst et al. 2020).

To be conservative, we scale our encounter rate to a local density
$n_{*}=0.01~{}{\rm pc}^{-3}$ of stars of mean velocity
$V_{*}=(40~{}{\rm~{}km~{}s}^{-1})V_{40}$, passing through an extended LISM
with $R_{\rm LISM}=(15~{}{\rm pc})R_{15}$. The collision rate with local
clouds is

$k_{\rm coll}=n_{*}(\pi R_{\rm LISM}^{2})f_{A}V_{*}\approx(3\times
10^{-4}~{}{\rm yr}^{-1})\,V_{40}R_{15}^{2}f_{A}\;.$ (1)

For $f_{A}\approx$ 0.5–1.0, the encounter rate is one collision every
3500–7000 yr. A typical stellar crossing time through a cloud of mean length
$\ell\approx 4R_{\rm cl}/3\approx 3$ pc is

$t_{\rm cross}\approx\ell/V_{*}\approx(70,000~{}{\rm yr})(\ell/3~{}{\rm
pc})V_{40}^{-1}\;,$ (2)

and the number of active bow-shock encounters within 15 pc is

$N_{\rm bow}=k_{\rm coll}\,t_{\rm cross}\approx 20\,(\ell/3~{}{\rm
pc})f_{A}R_{15}^{2}\;.$ (3)

The chances of detection could be significantly higher for larger cloud
complexes beyond the LISM (Frisch et al. 2011). However, the local clouds lie
inside a low-density cavity, and the nearest star-forming clouds are more
distant, including gas in the Taurus molecular cloud (140 pc) and the
Scorpius-Centaurus association (130 pc).

###  2.2 Estimated Bow Shock Sizes

The two Voyager missions passed through the solar-wind termination shock at
distances of 94 and 84 AU (Stone et al. 2005; Richardson et al. 2008; Stone et
al. 2008). Low-mass stars like the Sun, with weak winds and low mass-loss
rates, are expected to produce astrospheres preceded by bow shocks with
transverse radial extents $R_{\rm wake}\approx$100–300 AU, similar to that of
the Sun’s heliosheath and termination shock222Observations from the
Interstellar Boundary Explorer (McComas et al. 2012; Zank et al. 2013; McComas
et al. 2019) suggest that the solar wind termination may be in a “bow wave”
rather than a strong shock, owing to the influence of inclined magnetic fields
in the heliosphere and LISM. This interpretation has been questioned by
Scherer & Fichtner (2014), who noted the importance of including both He+ and
H+ in the speeds of Alfvén and magnetosonic waves. of the solar wind. These
estimates are consistent with results from an HST survey (Wood et al. 2021) of
astrospheric Ly$\alpha$ absorption toward nearby M-dwarf stars. For stars with
super-magnetosonic motions through the ISM, the hydrogen wall and bow shock in
MHD models (e.g., Zank et al. 2013) precede the wind termination shock by a
substantial factor that depends on the stellar speed and magnetosonic Mach
number relative to the surrounding medium. This warm, extended gas is
influenced in density, temperature, and wave speeds. The trailing wakes
typically have transverse radial extent 2–3 times the size of the forward
radius of the termination shock. Most low-mass stars have mass-loss rates
comparable to (and sometimes less than) that of the Sun, $\dot{M}\approx
4\times 10^{-14}~{}M_{\odot}~{}{\rm yr}^{-1}$. At the typical 3–10 pc
distances to the local interstellar clouds, the bow-shock angular sizes could
range from 10–100 arcsec. More distant bow shocks would be smaller (0.1 arcsec
at 1 kpc) and probably undetectable.

A few, more distant fast-moving OB stars with much stronger winds have
observed arc-like features, attributed to bow shocks of parsec size (Van Buren
& McCray 1988; Mackey et al. 2016). An O-star moving at velocity $V_{*}\approx
40~{}{\rm km~{}s}^{-1}$ would produce a wind-termination shock of approximate
size,

$r_{b}\approx\left[\frac{\dot{M}_{w}V_{w}}{4\pi\mu n_{\rm
H}}\right]^{1/2}V_{*}^{-1}\approx(1.0~{}{\rm pc})\,\dot{M}_{-6}^{1/2}n_{\rm
H}^{-1/2}\left(\frac{V_{w}}{1000~{}{\rm
km~{}s}^{-1}}\right)^{1/2}\left(\frac{V_{*}}{40~{}{\rm
km~{}s}^{-1}}\right)^{-1}\;.$ (4)

This termination shock in O-star winds will be preceded a wind-driven bubble
of high-pressure gas (Weaver et al. 1977; Mackey et al. 2016), whose shape can
be distorted by supersonic stellar motion. Here, we express the stellar mass-
loss rate as $\dot{M}=(10^{-6}~{}M_{\odot}~{}{\rm yr}^{-1})\dot{M}_{-6}$ and
equate the wind ram pressure, $\dot{M}_{w}V_{w}/4\pi r^{2}$, with the
pressure, $\rho_{\rm ISM}V_{*}^{2}$, of the ISM of mass density $\rho_{\rm
ISM}=1.4m_{\rm H}n_{\rm H}$. At $d\approx 1$ kpc, the angular sizes would be
$\theta_{b}=r_{b}/d\approx 10^{-3}$ rad $(3.4^{\prime}~{}d_{\rm kpc}^{-1})$.
The above estimate of $r_{b}$ is uncertain by factors of 1.5–2.0 when one
includes the effects of magnetic fields.

###  2.3 Estimated Luminosity and Surface Brightness

Intermediate velocity shocks moving into H I gas will produce optical and UV
emission lines that could be signatures of fast-star passages through local
clouds. These include Ly$\alpha$, metal-ion emission lines, and two-photon
continuum from the metastable (2s) state of H I (Shull & McKee 1979). An
analysis of the bow-shock luminosity and surface brightness requires
calculations of the photon flux, $\Phi_{0}$, which is proportional to the
flux, $n_{\rm H}v_{s}$, of hydrogen atoms entering the shock front. The
radiative shock models (Shull & Silk 1979) of primeval material (H and He)
found photon fluxes that can be fitted to expressions rising steeply for
$V_{s}$ between 15 and 100 km s-1,

$\displaystyle\Phi_{0}(2\gamma)$ $\displaystyle=$ $\displaystyle(2.2\times
10^{4}~{}{\rm phot~{}cm}^{-2}~{}{\rm s}^{-1})(n_{\rm H}/0.1~{}{\rm
cm}^{-3})V_{40}^{2.18}\;,$ (5) $\displaystyle\Phi_{0}({\rm Ly}\alpha)$
$\displaystyle=$ $\displaystyle(2.7\times 10^{5}~{}{\rm phot~{}cm}^{-2}~{}{\rm
s}^{-1})(n_{\rm H}/0.1~{}{\rm cm}^{-3})V_{40}^{2.12}\;.$ (6)

We use these expressions to estimate the luminosities of Ly$\alpha$ and two-
photon emission produced in a bow shock of radius $r_{b}$ from a low-mass star
located at distance $d$. Because of the approximate nature of these estimates,
we consider plane-parallel shocks at small angles from the forward direction
of the star. More careful 3D models with MHD codes would clearly be warranted
in the future. Most of the emission comes from a polar cap with opening angle
$\theta_{\rm cap}\approx 30^{\circ}$ relative to the stellar velocity vector.
This cap occupies a fraction $(1-\cos\theta_{\rm cap})$ of the $2\pi
r_{b}^{2}$ area of the forward hemisphere, with a photon luminosity,

$S_{\rm phot}=2\pi r_{b}^{2}\Phi_{0}(1-\cos\theta_{\rm cap})\;,$ (7)

and its projected solid angle on the sky is approximately

$\Omega_{\rm arc}\approx\frac{2r_{b}\theta_{\rm cap}(\Delta
r)}{d^{2}}=\left(\frac{2r_{b}^{2}}{d^{2}}\right)\theta_{\rm
cap}(1-\cos\theta_{\rm cap})\;.$ (8)

Here, we adopt an angular length, $2\theta_{\rm arc}=2\theta_{\rm cap}$, for
the arc (angles in radians) and a projected arc width $\Delta
r=r_{b}(1-\cos\theta_{\rm cap})$. Using eq. (7), the arc surface brightness is

$B_{\rm arc}=\frac{S_{\rm phot}/4\pi d^{2}}{\Omega_{\rm
arc}}=\frac{\Phi_{0}}{4\theta_{\rm cap}}\;.$ (9)

As expected for optically thin emission, the brightness is independent of both
$d$ and $r_{b}$.

For solar-type stars passing through the LISM at $V_{*}\approx(40$ km
s${}^{-1})V_{40}$, we adopt $r_{b}=100$ AU, $n_{\rm H}\approx 0.1$ cm-3, and
$\theta_{\rm cap}\approx 30^{\circ}$ to find

$\displaystyle S_{\rm phot}(2\gamma)$ $\displaystyle=$
$\displaystyle(4.2\times 10^{34}~{}{\rm phot~{}s}^{-1})\left(\frac{n_{\rm
H}}{0.1~{}{\rm cm}^{-3}}\right)V_{40}^{2.18}\;,$ (10) $\displaystyle S_{\rm
phot}({\rm Ly}\alpha)$ $\displaystyle=$ $\displaystyle(5.1\times
10^{35}~{}{\rm phot~{}s}^{-1})\left(\frac{n_{\rm H}}{0.1~{}{\rm
cm}^{-3}}\right)V_{40}^{2.12}\;.$ (11)

With the 10.2 eV energies of the upper (2s, 2p) levels, these luminosity
coefficients sum to $9\times 10^{24}$ erg s-1, comparable to the kinetic
energy of interstellar gas entering the bow shock over the cap portion of the
front,

$\frac{\rho_{\rm ISM}V_{s}^{3}}{2}(2\pi r_{b}^{2})(1-\cos\theta_{\rm
cap})\approx(1.4\times 10^{25}~{}{\rm erg~{}s}^{-1})\left(\frac{n_{\rm
H}}{0.1~{}{\rm cm}^{-3}}\right)V_{40}^{3}\;.$ (12)

Additional heating of the dust in the gas assembled by the bow front will come
from stellar radiation, particularly in the case of OB stars. The Ly$\alpha$
photons will scatter in the low-density local H I, but owing to low H I column
densities only a small fraction will be absorbed by dust grains. However, some
of these bow shocks may also be seen in infrared emission, as observed for OB-
stars and young-star outflows with the Spitzer Space Telescope (Winston et al.
2012; Kobulnicki et al. 2016) and the Wide-field Infrared Survey Explorer
(Peri et al. 2015).

Most of the surface brightness of the bow shock comes from the material
directly in front of the star, visible on the sky as a small arc of emission
(Kobulnicky et al. 2016). Adopting $\theta_{\rm cap}=30^{\circ}$ ($\pi/6$
radians) and $\Delta r/r_{b}\approx 0.13$, and using equations (5) and (6) for
$\Phi_{0}$, we estimate photon surface brightnesses, $B_{\rm
arc}\approx(3\Phi_{0}/2\pi)$ of the arc,

$\displaystyle B_{\rm cap}(2\gamma)$ $\displaystyle\approx$
$\displaystyle(10,500~{}{\rm phot~{}cm}^{-2}~{}{\rm s}^{-1}~{}{\rm
sr}^{-1})V_{40}^{2.18}\left(\frac{n_{\rm H}}{0.1~{}{\rm cm}^{-3}}\right)\;,$
(13) $\displaystyle B_{\rm cap}({\rm Ly}\alpha)$ $\displaystyle\approx$
$\displaystyle(129,000~{}{\rm phot~{}cm}^{-2}~{}{\rm s}^{-1}~{}{\rm
sr}^{-1})V_{40}^{2.25}\left(\frac{n_{\rm H}}{0.1~{}{\rm cm}^{-3}}\right)\;.$
(14)

The above brightness coefficients (in units photons cm-2 s-1 sr-1) can also be
expressed as 0.13 and 1.6 Rayleighs. Thus, for stars in the thin-disk and
thick-disk populations, with rms velocities of 37 km s-1 and 60 km s-1
respectively, the Ly$\alpha$ surface brightnesses could range from 2–4
Rayleighs in gas with $n_{\rm H}\approx 0.1$ cm-3. Even greater brightnesses
might appear around several of the fast local stars listed in Table 1.

##  3 Turbulent Heating of the LISM gas

The LISM clouds penetrated by transiting low-mass stars will be affected by
turbulence generated by transiting fast stars at various shocks and in the
wake. In clouds of pathlength $\ell=$ 2–3 pc, stars moving at 40–80 km s-1
will spend 30,000–100,000 years, leaving wakes of dynamically disturbed gas.
Owing to downstream flows from the bow shock, the transverse sizes of the
wakes will be larger than the forward termination shock and astrosheath (e.g.,
Mac Low et al. 1991; Mackey et al. 2016). The fractional volume of the cloud
affected by the turbulent wake is

$f_{V}\approx\left(\frac{R_{\rm wake}}{R_{\rm cl}}\right)^{2}\approx(2\times
10^{-6})\left[\frac{R_{\rm wake}/300~{}{\rm AU}}{R_{\rm cl}/{\rm
pc}}\right]^{2}\;.$ (15)

Note that $R_{\rm wake}$ is the transverse size of the post-bowshock flow, not
its length. The total affected portion of local clouds includes both the
direct “bore holes” and extended regions influenced by heat generated by MHD
wave damping. As described by Spitzer (1982) and Ferriére et al. (1988), the
waves arise from reflected shocks produced by interactions with density
irregularities encountered by the bow shock. After a nonlinear cascade to
smaller eddies, much of the turbulence will decay into heat, on an eddy-
turnover time,

$t_{\rm eddy}\approx R_{\rm eddy}/V_{\rm eddy}\approx(10^{3}~{}{\rm
yr})\left(\frac{r_{b}}{100~{}{\rm AU}}\right)\left(\frac{3~{}{\rm
km/s}}{V_{\rm eddy}}\right)\;.$ (16)

Here, we assumed the initial eddy to have a characteristic size $R_{\rm
eddy}\approx r_{b}=100~{}{\rm AU}$ at the termination shock scale and an
initial turbulent velocity $V_{\rm eddy}\approx 3~{}{\rm km~{}s}^{-1}$. After
cascading to smaller scales through nonlinear effects, the turbulent energy is
dissipated. The turbulent heating rate can be estimated by dividing the
kinetic energy density ($\rho_{\rm gas}V_{\rm eddy}^{2}$) by the eddy turnover
time,

$\epsilon_{\rm turb}\approx\rho_{\rm gas}V_{\rm eddy}^{3}/R_{\rm
eddy}\approx(4\times 10^{-23}~{}{\rm erg~{}s}^{-1})n_{\rm
H}\left(\frac{r_{b}}{100~{}{\rm AU}}\right)\left(\frac{3~{}{\rm
km~{}s}^{-1}}{V_{\rm eddy}}\right)^{3}\;,$ (17)

where we assumed $\rho_{\rm gas}=1.4n_{\rm H}m_{\rm H}$. The turbulent heating
time would be

$t_{\rm heat}\approx\frac{3nkT/2}{\epsilon_{\rm turb}}\approx(10^{3}~{}{\rm
yr})\left(\frac{r_{b}}{100~{}{\rm AU}}\right)^{-1}\left(\frac{V_{\rm
eddy}}{3~{}{\rm km~{}s}^{-1}}\right)^{-3}\left(\frac{T}{6000~{}{\rm
K}}\right)\;.$ (18)

Because this time is much shorter than the stellar crossing time through the
cloud, the thermal state of gas in the bow-shock wake will be time dependent.
The situation will involve transient heating over multiple time scales, with
$t_{\rm heat}<t_{\rm cool}$ and $t_{\rm cool}\approx t_{\rm cross}$,

The post-shock wake will cool gradually, as the turbulent eddies quickly
cascade down to smaller scales. Gas temperatures will rise transiently to
$T\approx 20,000-30,000$ K, where H I can be collisionally ionized. Radiative
cooling by H I (Ly$\alpha$) and metal forbidden lines will reduce the
temperature to $T\leq 10^{4}$ K on an isobaric cooling timescale,

$t_{\rm cool}=\frac{5n_{\rm tot}kT/2}{n_{e}n_{H}{\cal
L}(T)}\approx(50,000~{}{\rm yr})\left(\frac{n_{\rm H}}{0.1~{}{\rm
cm}^{-3}}\right)^{-1}\left(\frac{T}{20,000~{}{\rm K}}\right){\cal
L}_{-22}^{-1}\;,$ (19)

comparable to the stellar crossing time. We assume a partially ionized medium
with $n_{e}\approx 0.5n_{\rm H}$ and $n_{\rm tot}\approx 1.6n_{\rm H}$, and a
cooling rate coefficient ${\cal L}\approx(10^{-22}~{}{\rm erg~{}cm}^{3}~{}{\rm
s}^{-1}){\cal L}_{-22}$ at 20,000 K. Analyzing the time-dependent evolution of
ionization and temperature behind bow shocks is a problem in hydrodynamics and
thermodynamics worth further study. The effects of magnetic fields in the
turbulent wakes (Baalmann et al. 2021 add to the complexities, including
reconnection and dissipation of magnetic field energy.

In our model we have ignored cooling through all elements other than hydrogen
and ignored external heat inputs. By $10^{4}$ K the cooling power from
hydrogen decreases with temperature and will compete with cooling by forbidden
and fine structure lines of metals. As cooling proceeds, the internal heat is
proportionally reduced and heating by external agents (e.g., ionization by
cosmic rays, photoelectric dust heating, ionizing soft X-ray radiation) will
become important. The recombination time scale is about $10^{5}n_{\rm H}^{-1}$
yr. The ionization time per atom from the diffuse EUV field which ionizes the
warm ionized medium is $10^{6}$ yr.

##  4 Conclusion and Future Prospects

This paper analyzes the interstellar bow-shocks in the LISM generated around
passing stars from the Galactic thin disk and thick disk with velocities of
40–70 km s-1. We estimate that bow shocks transit the local H I clouds every
3500–7000 yr. On average, there could be 10–20 active bow shocks inside clouds
within 15 pc of the Sun, with stellar transit times varying from
$10^{4}$–$10^{5}$ yr. At typical local cloud distances of 3–10 pc, the bow
shocks and astrospheres around low-mass stars likely have 100–300 AU
transverse radial extent and angular sizes of 10–100 arcsec. Their turbulent
transit wakes of passing stars might cover an area fraction
$f_{A}\approx(R_{\rm wake}/R_{\rm cl})^{2}\approx 10^{-3}$ of nearby clouds.
They could explain the 2–3 km s-1 turbulent widths and temperature variations
on scales less than 5100 AU in LISM clouds (Linsky et al. 2022), as well as
elevated rotational populations of H2 and endothermic production of CH+. Bow
shocks around rarer massive stars and AGB stars would produce shocks of parsec
scale and angular sizes of several arcmin (at $\sim 1$ kpc).

Detecting the low-mass bow shocks will require wide-field surveys, either in
the infrared (dust emission) or in the ultraviolet (for Ly$\alpha$ and two-
photon emission). The IR detections could be a task for the WISE data base,
and Ly$\alpha$ emission could be a project for the Alice spectrograph aboard
the New Horizons (NH) spacecraft during its extended mission. During six
great-circle scans (at 40-47 AU) the NH spectrograph observed interplanetary
Ly$\alpha$ (Gladstone et al. 2021) with a uniform background of $43\pm 3$
Rayleighs over the sky, interpreted as the Galactic Ly$\alpha$ background.
There is also resonantly scattered solar Ly$\alpha$, which at the current NH
distance (56 AU from the Sun) has fallen to about the same level. Thus, the
typical brightness in any direction (away from the Sun) would be 80–100 R. Our
calculations of far-UV radiation from the forward caps of 40–80 km s-1 bow
shocks suggest that arc-like structures around low-mass stars transiting local
clouds with $n_{\rm H}\approx 0.1~{}{\rm cm}^{-3}$ should have surface
brightnesses of 2–8 Rayleigh in Ly$\alpha$ and 0.1–0.4 Rayleigh in the UV two-
photon continuum. These could possibly be seen with NH and with slit
spectroscopy on the proposed UVEX mission (Kulkarni et al. 2021).

Examples of bow shocks around runaway O-type stars and AGB stars have already
been found in infrared emission (Kobulnicki et al. 2016; Peri et al. 2015; Cox
et al. 2012). Wide-angle searches in Ly$\alpha$, far-UV continua, and infrared
emission for lower luminosity examples could detect the more numerous cases of
low-mass stars transiting the LISM. Their signatures might also be detectable
in non-thermal radio emission from particles accelerated at the termination
shock and bow-shock vicinity. Although the predicted stellar transit wakes
have small covering fractions of LISM clouds, both by volume and on the sky,
the effects of down-stream turbulent dissipation, MHD wave heating, and
energetic particle acceleration could be more widespread throughout the local
clouds.

We thank Jeffrey Linsky and Seth Redfield for discussions on gas structures in
the local interstellar medium, Randy Gladstone for insights on Ly$\alpha$
scanning observations with the New Horizons spacecraft, and Brian Wood for
providing characteristics of the winds and astrospheric absorption around
nearby solar type stars. We also thank the anonymous referee for introducing
us to important papers on MHD models of astrospheres.

## References

  * (1)
  * (2) Anguiano, B., Majewski, S. R., Hayes, C. R., et al. 2020, AJ, 160, 43
  * (3)
  * (4) Baalmann, L. R., Scherer, K., Kleimann, J., et al. 2021, A&A, 650, A36
  * (5)
  * (6) Cox, N. L. J., Kerschbaum, F., van Marle, A.-J., et al. 2012, A&A, 537, A35
  * (7)
  * (8) Decin, L., Cox, N. L. J., Royer, P., et al. 2012, A&A, 548, A113
  * (9)
  * (10) Everall, A., Belokurov, V., Evans, N. W., et al. 2022, MNRAS, 511, 3863
  * (11)
  * (12) Falgarone, E., & Puget, J.-L. 1995, A&A, 293, 840
  * (13)
  * (14) Ferrière, K. M., Zweibel, E. G., & Shull, J. M. 1988, ApJ, 332, 984
  * (15)
  * (16) Frisch, P. C., Redfield, S., & Slavin, J. D. 2011, ARA&A, 49, 237
  * (17)
  * (18) Gladstone, G. R., Pryor, W. R., Hall, D. T., et al. 2021, AJ, 162, 241
  * (19)
  * (20) Gry, C., Boulanger, F., Nehmé, C., et al. 2002, A&A, 391, 675
  * (21)
  * (22) Herbst, K., Scherer, K., Ferreira, S. E. S., et al. 2020, ApJ, 897, L26
  * (23)
  * (24) Ingalls, J. G., Bania, T.M., Boulanger, F., et al. 2011, ApJ, 743, 174
  * (25)
  * (26) Jensen, A. G., Snow, T. P., Sonneborn, G., & Rachford, B. L. 2010, ApJ, 711, 1236
  * (27)
  * (28) Kobulnicky, H. A., Chick, W. T., Schurhammer, D. P., et al. 2016, ApJS, 227, 18
  * (29)
  * (30) Kuhn, M. A., Hillenbrand, L. A., Sills, A., et al. 2019, ApJ, 870, 32
  * (31)
  * (32) Kulkarni, S. R., Harrison, F. A., Grefenstette, B. W., et al. 2021, Science with the Ultraviolet Explorer (UVEX), arXiv:2111.15608
  * (33)
  * (34) Kulkarni, S. R., & Shull, J. M. 2023, PASP, submitted (arXiv: 2302.13963)
  * (35)
  * (36) Kulsrud, R. M. 2005, Plasma Physics for Astrophysics (Princeton: Princeton Univ. Press), p 108
  * (37)
  * (38) Linsky, J. L., Redfield, S., & Tilipman, D. 2019, ApJ, 886, 41
  * (39)
  * (40) Linsky, J. L., Redfield, S., Ryder, D., & Chasan-Taber, A. 2022, AJ, 164, 106
  * (41)
  * (42) Mackey, J., Haworth, T. J., Gvaramadze, V. V., et al. 2016, A&A, 586, A114
  * (43)
  * (44) Mac Low, M.-M., Van Buren, D., Wood, D. O. S., & Churchwell, E. 1991, ApJ, 369, 395
  * (45)
  * (46) McComas, D. J., Alexashov, D., Bzowski, M., et al. 2012, Science, 336, 1291
  * (47)
  * (48) McComas, D. J., Rankin, J. S., Schwadron, N. A., & Swaczyna, P. 2019, ApJ, 884, 145
  * (49)
  * (50) Mihalas, D., & Binney, J. 2005, Galactic Astronomy, Chapter 4 (San Francisco: W. H. Freeman)
  * (51)
  * (52) Moseley, E. R., Draine, B. T., Tomida, K., & Stone, J. M. 2021, MNRAS, 500, 3290
  * (53)
  * (54) Myers, A. F., McKee, C. F., & Li, P. S. 2015, MNRAS, 453, 2747
  * (55)
  * (56) Nehmé, C., Le Bourlot, J., Boulanger, F., et al. 2008, A&A, 483, 485
  * (57)
  * (58) Nordström, B., Major, M., Andersen, J., et al. 2004, A&A, 418, 989
  * (59)
  * (60) Peri, C. S., Benaglia, P., & Isequilla, N.L. 2015, A&A, 578, A45
  * (61)
  * (62) Redfield, S., & Linsky, J. L. 2008, ApJ, 673, 283
  * (63)
  * (64) Reylé, C., Jardine, K., Fouqué, P., et al. 2022, in $21^{\rm st}$ Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. A. S. Brun, J. Bouvier, P. Petit, arXiv:2302.02810
  * (65)
  * (66) Richardson, J. D., Kasper, J. C., Wang, C., et al. 2008, Nature, 454, 63
  * (67)
  * (68) Scherer, K., & Fichtner, H. 2014, ApJ, 782, 25
  * (69)
  * (70) Shull, J. M., & McKee, C. F. 1979, ApJ, 227, 131
  * (71)
  * (72) Shull, J. M., & Silk, J. 1979, ApJ, 234, 427
  * (73)
  * (74) Shull, J. M., Danforth, C. W., & Anderson, K. L. 2021, ApJ, 911, 55
  * (75)
  * (76) Shull, J. M., Darling, J., & Danforth, C. W. 2021, ApJ, 914, 18
  * (77)
  * (78) Slavin, J. D., & Frisch, P. C. 2008, A&A, 491, 53
  * (79)
  * (80) Spitzer, L. 1982, ApJ, 262, 315
  * (81)
  * (82) Stone, E. C., Cummings, A. C., McDonald, F. B., et al. 2005, Science, 309, 2017
  * (83)
  * (84) Stone, E. C., Cummings, A. C., McDonald, F. B., et al. 2008, Nature, 454, 71
  * (85)
  * (86) van Buren, D., & McCray, R. 1988, ApJ, 329, L93
  * (87)
  * (88) Vieira, K., Carraro, G., Korchagin, V., et al. 2022, ApJ, 932, 28
  * (89)
  * (90) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377
  * (91)
  * (92) Winston, E., Wolk, S.J., Bourke, T.L., et al. 2012, ApJ, 744, 126
  * (93)
  * (94) Wood, B. E., Müller, H.-R., Redfield, S., et al. 2021, ApJ, 915, 37
  * (95)
  * (96) Zank, G. P., Heerikhuisen, J., Wood, B. E., et al. 2013, ApJ, 763, 20
  * (97)

Table 1: Nearby Fast Starsaa List of local “fast stars” with spectral types and Gaia parallax distances. $V_{\rm ISM}$ is the ISM flow speed seen by the star in the stellar rest frame (Wood et al. 2021), and $V_{\rm rad}$ is the star’s heliocentric radial velocity. The last two columns list gas temperatures and turbulent velocities, inferred from an ensemble of absorption-line widths of hydrogen (both H I and D I) and various metal ions (Linsky et al. 2022).  Star | SpT | $V_{\rm ISM}$ | $V_{\rm rad}$ | $d$ | $T_{\rm gas}$ | $V_{\rm turb}$
---|---|---|---|---|---|---
|  | (km/s) | (km/s) | (pc) | (K) | (km/s)
61 Cyg A | K5 V | 86 | -65.2 | 3.50 | 6500 | 2.08
GJ 887 | M2 V | 85 | +8.8 | 3.29 | 6000 | 1.65
GJ 436 | M3 V | 79 | +9.6 | 9.78 | … | …
GJ 205 | M1.5 V | 70 | +8.6 | 5.70 | 4500 | 2.37
$\epsilon$ Ind | K5 V | 68 | -40.0 | 3.64 | 8340 | 1.97
61 Vir | G6.5 V | 51 | -7.8 | 8.53 | 6600 | 1.60
GJ 892 | K3 V | 49 | -18.5 | 6.54 | … | …
GJ 860A | M3 V | 47 | -24.0 | 4.01 | 4780 | 2.24
EV Lac | M3.5 V | 45 | +0.3 | 5.05 | 3030 | 2.21
$\pi^{1}$ UMa | G1.5 V | 43 | -12.8 | 14.4 | 2450 | 2.47
36 Oph A | K2 V | 40 | +0.4 | 5.95 | 5870 | 2.33
GJ 173 | M1 V | 38 | -6.8 | 11.2 | … | …
$\delta$ Eri | K0 IV | 37 | -6.2 | 9.09 | 3650 | 2.25
$\xi$ Boo A | G7 V | 32 | +1.6 | 6.71 | 5310 | 1.68
GJ 338 A | M0 V | 29 | +11.5 | 6.33 | 6650 | 2.29
$\delta$ Pav | G8 IV | 29 | -21.5 | 6.10 | 9100 | 2.44
GJ 15 A | M2 V | 28 | +11.7 | 3.56 | 7800 | 2.81
$\epsilon$ Eri | K2 V | 27 | +16.4 | 3.22 | 7410 | 2.03
$\alpha$ Cen B | K1 V | 25 | -22.6 | 1.35 | 5500 | 1.37
Sun | G2 V | 23 | … | … | 6400 | …
Means |  |  |  |  | 6638 | 2.54


# Measurements of $\boldsymbol{B\rightarrow K\pi}$ and
$\boldsymbol{B\rightarrow\pi\pi}$ Branching Fractions and
$\boldsymbol{\mathcal{A}_{CP}}$ Asymmetries at Belle II

Shu-Ping Lin on behalf of the Belle II Collaboration  University of Padova and
INFN, Via Marzolo 8, 35131, Padova, PD, Italy

###### Abstract

Analyses of $B$ meson decays to charmless hadronic final states are an
important part of the Belle II program. They are sensitive to effects from
non-standard model physics and provide experimentally precise constraints on
the weak interactions of quarks. We present recent Belle II results on
branching fractions and direct $CP$-violating asymmetries of the decays
$B^{0}\rightarrow K^{+}\pi^{-}$, $B^{+}\rightarrow K^{+}\pi^{0}$,
$B^{+}\rightarrow K^{0}\pi^{+}$, and $B^{0}\rightarrow K^{0}\pi^{0}$, and use
these to test the standard model through an isospin-based sum rule. In
addition, we measure the branching fraction and direct $CP$ asymmetry of the
decay $B^{+}\rightarrow\pi^{+}\pi^{0}$ and the branching fraction of the decay
$B^{0}\rightarrow\pi^{+}\pi^{-}$, which contribute towards the determination
of the Cabibbo–Kobayashi–Maskawa angle $\phi_{2}$. The data are collected with
the Belle II detector from the SuperKEKB asymmetric-energy $e^{+}e^{-}$
collider, consisting of $387\times 10^{6}$ $\Upsilon(4S)\rightarrow B\bar{B}$
events. We obtain $-0.03\pm 0.13\pm 0.04$ for the sum rule, in agreement with
the standard model expectation of zero and with a precision comparable to the
best existing determinations.

##  1 Motivation

###  1.1 $B\rightarrow K\pi$ Decays: Isospin Sum Rule

For $B\rightarrow K\pi$ decays, one may construct a sum rule exploiting
isospin symmetry using linear combinations of the branching fractions and $CP$
asymmetries [1]. This reduces the impact of theoretical uncertainties. The sum
rule parameter

$I_{K\pi}=\mathcal{A}_{K^{+}\pi^{-}}+\mathcal{A}_{K^{0}\pi^{+}}\cdot\frac{\mathcal{B}_{K^{0}\pi^{+}}}{\mathcal{B}_{K^{+}\pi^{-}}}\frac{\tau_{B^{0}}}{\tau_{B^{+}}}-2\mathcal{A}_{K^{+}\pi^{0}}\cdot\frac{\mathcal{B}_{K^{+}\pi^{0}}}{\mathcal{B}_{K^{+}\pi^{-}}}\frac{\tau_{B^{0}}}{\tau_{B^{+}}}-2\mathcal{A}_{K^{0}\pi^{0}}\cdot\frac{\mathcal{B}_{K^{0}\pi^{0}}}{\mathcal{B}_{K^{+}\pi^{-}}}$
(1)

is predicted to be zero in the Standard Model (SM) within $1\%$ [2, 3, 4],
therefore providing a null test for the SM. $\mathcal{A}_{K\pi}$ and
$\mathcal{B}_{K\pi}$ (neutral or charged $K$ and $\pi$) are the direct $CP$
asymmetry and the $CP$-averaged branching fraction of a $B\rightarrow K\pi$
decay, and $\tau_{B^{0}}$ and $\tau_{B^{+}}$ are the neutral and charged $B$
meson lifetimes, respectively. Charge conjugation is implied unless otherwise
indicated. The time-integrated $CP$ asymmetry parameter is defined as

$\mathcal{A}_{X}=\frac{\Gamma(\bar{B}\rightarrow\bar{X})-\Gamma(B\rightarrow
X)}{\Gamma(\bar{B}\rightarrow\bar{X})+\Gamma(B\rightarrow X)},$ (2)

where $\Gamma$ is the partial decay width to the final state $X$.
$\mathcal{A}_{X}$ corresponds to the direct $CP$ asymmetry of charged $B$
mesons and, in the absence of $CP$ violation in flavor oscillations, also
neutral $B$ mesons. This is a good approximation in the case of
$B^{0}$-$\bar{B}^{0}$ mixing.

The current world average of the sum rule parameter is $I_{K\pi}=-0.13\pm
0.11$ [5], and its precision is limited by the statistical uncertainty of
$\mathcal{A}_{K^{0}\pi^{0}}$ [6].

###  1.2 $B\rightarrow\pi\pi$ Decays: Towards the Angle $\phi_{2}$

The precision of the angle $\phi_{2}$ (or $\alpha$)
$=\arg(-V_{td}V^{*}_{tb}/V_{ud}V_{ub}^{*})$, where $V_{ij}$ are the elements
of the Cabibbo–Kobayashi–Maskawa (CKM) quark-mixing matrix [7, 8], is a
limiting factor of the global testing power of the CKM model. $\phi_{2}$ is
accessible through the combined information of
$B^{0}\rightarrow\pi^{+}\pi^{-}$, $B^{+}\rightarrow\pi^{+}\pi^{0}$, and
$B^{0}\rightarrow\pi^{0}\pi^{0}$ decays [9, 10]. Precise measurements of the
branching fraction and $CP$ asymmetries of the decays
$B^{0}\rightarrow\pi^{+}\pi^{-}$ and $B^{+}\rightarrow\pi^{+}\pi^{0}$ provide
constraints through isospin symmetry relations in the determination of
$\phi_{2}$.

##  2 Analysis Strategy

The analyses of all the decays follow a similar strategy. A common selection
is applied to all relevant final state particles. Multivariate algorithms are
developed for each mode to suppress the continuum background, which originates
from $e^{+}e^{-}\rightarrow q\bar{q}$ processes, where $q$ indicates a $u$,
$d$, $s$, or $c$ quark. Event shape variables, results of the $B$ meson decay
vertex fit, flavor-tagger outputs (except for $B^{0}\rightarrow
K_{S}^{0}\pi^{0}$ decay), and the beam-constrained mass
$M_{bc}=\sqrt{(E_{beam}^{*}/c^{2})^{2}-(p_{B}^{*}/c)^{2}}$ are used to train a
boosted decision tree (BDT). $E_{beam}^{*}$ and $p_{B}^{*}$ are the beam
energy and $B$ meson momentum, respectively, in the center-of-mass frame. A
loose requirement is applied to the BDT output, removing $90\%-99\%$ of the
continuum background. In order to improve the $M_{bc}$ resolution for decays
containing a neutral $\pi^{0}$, the measured $\pi^{0}$ momentum in the
calculation of $M_{bc}$ is replaced by

$\boldsymbol{p}_{\pi^{0}}^{*^{\prime}}=\sqrt{(E_{beam}^{*}-E_{h}^{*})^{2}/c^{2}-m_{\pi^{0}}^{2}c^{2}}\times\frac{\boldsymbol{p}_{\pi^{0}}^{*}}{|\boldsymbol{p}_{\pi^{0}}^{*}|},$
(3)

where $E_{h}^{*}$ is energy of the other hadron in the $B$ candidate decay,
$\boldsymbol{p}_{\pi^{0}}^{*}$ is the measured $\pi^{0}$ momentum, and
$m_{\pi^{0}}$ is the known $\pi^{0}$ mass. This substitution improves the
resolution of $M_{bc}$ and reduces its correlation to $\Delta E$.

The signal yields are measured with unbinned maximum likelihood fits to the
energy difference $\Delta E=E_{B}^{*}-E_{beam}^{*}$ and the transformed BDT
output $C^{\prime}$. $E_{B}^{*}$ is the energy of the $B$ meson in the center-
of-mass frame, and the BDT output is transformed using probability integral
transformation [12]. The two-dimensional probability density functions (PDF)
of $\Delta E$ and $C^{\prime}$ is assumed to be a product of two independent
one-dimensional PDFs, as the correlation between the two variables is small.
The $CP$ asymmetries are determined using the charge of the final state
particles for flavor-specific channels. For the $B^{0}\rightarrow
K_{S}^{0}\pi^{0}$ channel the $B$ mesons decay to the same $CP$ eigenstate,
and a category-based flavor tagger is used to determine the flavor of the
$B^{0}$, calibrated using $B\rightarrow D^{(*)}h^{+}$ decays [13].

The analysis procedures are first studied with simulated events before being
applied to the data sample. Several abundant control channels with similar
features to the signal channels are used to correct for discrepancies between
data and simulation, and to estimate systematic uncertainties.

##  3 Result

The fit results are listed in Table 1, and the fit projections to $\Delta E$
are shown in Figure 1. All results [14] agree with the world averages, and our
precisions are comparable to the best existing determinations from Belle [15,
16, 17] and BABAR [18, 19, 20, 21] despite using a smaller data sample.

The branching fraction and direct $CP$ asymmetry of $B^{0}\rightarrow
K^{0}\pi^{0}$ decays have also been measured in an analysis of the decay time
evolution [22] based on the same data sample, but featuring a different event
selection. The result presented in Table 1 is the combination of both
analyses, where the statistical and systematic correlations are taken into
account. The result for the $CP$ asymmetry of $B^{0}\rightarrow K^{0}\pi^{0}$
and the branching fraction of $B^{0}\rightarrow\pi^{+}\pi^{-}$ are the most
precise measurement made by a single experiment to date.

The sum rule of Equation 1 is calculated using our measurements of the
branching fractions and $CP$ asymmetries,

$I_{K\pi}=-0.03\pm 0.13\pm 0.04,$ (4)

where the measured ratio $\tau_{B^{0}}/\tau_{B^{+}}=0.9273\pm 0.0033$ [6] is
used, and the correlations between uncertainties are taken into account.
Common systematic uncertainties are cancelled in the ratios of the branching
fractions in Equation 1. The sum rule result is consistent with the
theoretical expectation. It is statistically limited and has a similar
precision to that from the average of measurements from Belle, BABAR and LHCb
collaborations [5].

Figure 1: $\Delta E$ fit projections and data distributions for $B^{0}\rightarrow K^{+}\pi^{-}$ (upper left), $B^{0}\rightarrow\pi^{+}\pi^{-}$ (upper right), $B^{+}\rightarrow K^{+}\pi^{0}$ (middle left), $B^{+}\rightarrow\pi^{+}\pi^{0}$ (middle right), $B^{+}\rightarrow K_{S}^{0}\pi^{+}$ (lower left) and $B^{0}\rightarrow K_{S}^{0}\pi^{0}$ (lower right) decays. Table 1: Results of the branching fraction and $CP$ asymmetries using $362\text{ fb}^{-1}$ Belle II data. The first uncertainties are statistical and the second are systematic. | $\mathcal{B}\ [10^{-6}]$ | $\mathcal{A}_{{\it CP}}$
---|---|---
$B^{0}\rightarrow K^{+}\pi^{-}$ | $20.67\pm 0.37\pm 0.62$ | $-0.072\pm 0.019\pm 0.007$
$B^{0}\rightarrow\pi^{+}\pi^{-}$ | $5.83\pm 0.22\pm 0.17$ | –
$B^{+}\rightarrow K^{+}\pi^{0}$ | $13.93\pm 0.38\pm 0.71$ | $0.013\pm 0.027\pm 0.005$
$B^{+}\rightarrow\pi^{+}\pi^{0}$ | $5.10\pm 0.29\pm 0.27$ | $-0.081\pm 0.054\pm 0.008$
$B^{+}\rightarrow K^{0}\pi^{+}$ | $24.37\pm 0.71\pm 0.86$ | $0.046\pm 0.029\pm 0.007$
$B^{0}\rightarrow K^{0}\pi^{0}$ | $10.73\pm 0.63\pm 0.62$ | $-0.01\pm 0.12\pm 0.04$

## References

## References

  * [1] M. Gronau, Phys. Lett. B 627, 82 (2005).
  * [2] M. Gronau and J.L. Rosner, Phys. Rev. D 59, 113002 (1999).
  * [3] G. Bell, M. Beneke, T. Huber, X.-Q. Li, Phys. Lett. B 750, 348 (2015).
  * [4] G. Bell, M. Beneke, T. Huber, X.-Q. Li, J. High Energy Phys. 04, 055 (2020).
  * [5] Y. Amhis et al. (HFLAV Collaboration), Phys. Rev. D 107, 052008 (2023).
  * [6] R. L. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
  * [7] N. Cabibbo, Phys. Rev. Lett. 10, 531 (1963).
  * [8] M. Kobayashi and T. Maskawa, Prog. Theor. Phys. 49, 652 (1973).
  * [9] M. Gronau and D. London, Phys. Rev. Lett. 65, 3381 (1990).
  * [10] J. Charles, et al., Eur. Phys. J. C 77, 574 (2017).
  * [11] F.Abudinén et al. (Belle II Collaboration), Phys. Rev. D 107, 112009 (2023).
  * [12] A. Hoecker et al., arXiv:physics/0703039.
  * [13] F.Abudinén et al. (Belle II Collaboration), Eur. Phys. J. C 82, 283 (2022).
  * [14] I. Adachi et al. (Belle II Collaboration), Phys. Rev. D 109, 012001 (2024)
  * [15] Y.-T. Duh et al. (Belle Collaboration), Phys. Rev. D 87, 031103 (2013).
  * [16] J. Dalseno et al. (Belle Collaboration), Phys. Rev. D 88, 092003 (2013).
  * [17] M. Fujikawa et al. (Belle Collaboration), Phys. Rev. D 81, 011101 (2010).
  * [18] J. P. Lees et al. (BABAR Collaboration), Phys. Rev. D 87, 052009 (2013).
  * [19] B. Aubert et al. (BABAR Collaboration), Phys. Rev. D 75, 012008 (2007).
  * [20] B. Aubert et al. (BABAR Collaboration), Phys. Rev. D 76, 091102 (2007).
  * [21] B. Aubert et al. (BABAR Collaboration), Phys. Rev. D 79, 052003 (2009).
  * [22] I. Adachi et al. (Belle II Collaboration), Phys. Rev. Lett. 131, 111803 (2023).