11institutetext: Facultad de Ciencias Exactas, UNLP, Calle 47 y 115, CP(1900),
La Plata, Buenos Aires, Argentina.
11email<EMAIL_ADDRESS>22institutetext: Instituto Argentino de
Radioastronomía (CCT La Plata, CONICET), C.C.5, (1894) Villa Elisa, Buenos
Aires, Argentina.
22email<EMAIL_ADDRESS>33institutetext: Departament de Física
Quántica i Astrofísica, Institut de Ciéncies del Cosmos (ICCUB), Universitat
de Barcelona (IEEC-UB), Martí i Franquès 1, E-08028 Barcelona, Spain.
33email<EMAIL_ADDRESS>
# Non-thermal emission in hyper-velocity and semi-relativistic stars

J. R. Martinez 1122 S. del Palacio 22 V. Bosch-Ramon & G. E. Romero 223322

###### Abstract

Context. There is a population of runaway stars that move at extremely high
speeds with respect to their surroundings. The fast motion and the stellar
wind of these stars, plus the wind-medium interaction, can lead to particle
acceleration and non-thermal radiation.

Aims. We characterise the interaction between the winds of fast runaway stars
and their environment, in particular to establish their potential as cosmic-
ray accelerators and non-thermal emitters.

Methods. We model the hydrodynamics of the interaction between the stellar
wind and the surrounding material. We self-consistently calculate the
injection and transport of relativistic particles in the bow shock using a
multi-zone code, and compute their broadband emission from radio to
$\gamma$-rays.

Results. Both the forward and reverse shocks are favourable sites for particle
acceleration, although the radiative efficiency of particles is low and
therefore the expected fluxes are in general rather faint.

Conclusions. We show that high-sensitivity observations in the radio band can
be used to detect the non-thermal radiation associated with bow shocks from
hypervelocity and semi-relativistic stars. Hypervelocity stars are expected to
be modest sources of sub-TeV cosmic rays, accounting perhaps for a $\sim 0.1$%
of that of galactic cosmic rays.

###### Key Words.:

Radiation mechanisms: non-thermal – Stars: winds, outflows – Acceleration of
particles – Shock waves

##  1 Introduction

Massive stars have intense ultraviolet (UV) radiation fields that accelerate
the surface material, launching powerful supersonic winds. These stellar winds
interact with the interstellar medium (ISM) generating two shock fronts: a
forward shock (FS) that propagates through the ISM and a reverse shock (RS)
that propagates through the stellar wind (Weaver et al., 1977). These shocks
are potential sites for non-thermal phenomena; they have been detected in a
few occasions (Prajapati et al., 2019; Sánchez-Ayaso et al., 2018), and have
been suggested to produce galactic cosmic rays up to PeV energies (e.g.
Aharonian et al., 2019; Morlino et al., 2021).

Stars that have a supersonic velocity with respect to the local ISM are known
as runaway stars. In this scenario, the interaction region becomes bow-shaped
(van Buren & McCray, 1988), and thus the whole interaction structure is
usually called bow shock (BS). The FS compresses and heats up dust and gas
that emit mostly infrared (IR) and optical radiation (Peri et al., 2012;
Kobulnicky et al., 2016), but BSs can also accelerate particles up to
relativistic energies via diffusive shock acceleration (DSA). These particles,
in turn, can interact with matter and electromagnetic fields, producing
broadband non-thermal (NT) radiation (as predicted by, e.g., del Valle &
Romero, 2012; del Palacio et al., 2018; del Valle & Pohl, 2018). Nonetheless,
despite more than 700 BSs have been cataloged, NT emission has been clearly
detected only in one of them, BD+43°3654 (Benaglia et al., 2010, 2021). We
note that it has also been detected radio emission from the BS produced by the
high-mass X-ray binary Vela X-1, but in this case the nature (thermal or NT)
of the emission is still uncertain (van den Eijnden et al., 2021).

Hyper-velocity stars are the subclass of runaway stars with velocities of
hundreds to a few thousands of km s-1 (Brown, 2015). Recent observational
campaigns have catalogued hundreds of HVSs with the data provided by Gaia, and
models predict an ejection rate of HVSs of 10-4–10-5 yr-1 (Zhang et al.,
2013). In the case of massive HVSs, they are expected to generate strong BSs,
being promissory NT sources. The Hills mechanism (Hills, 1988) explains the
HVSs origin via a 3-body exchange between a stellar binary and a supermassive
black hole. The black hole disrupts the binary, ejecting one of its components
at great velocities. The velocity of ejection depends on the supermassive
black hole mass, and the total mass and semimajor axis of the binary.
Recently, this mechanism gained great support by the discovery of a $\sim
1700\leavevmode\nobreak\ {\rm km\;s^{-1}}$ A-type star ejected from Sgr A∗
(Koposov et al., 2020).

Tutukov & Fedorova (2009) predicted a putative subclass of HVSs with semi-
relativistic velocities, called semi-relativistic stars (SRSs). Numerical
simulations support this prediction (Loeb & Guillochon, 2016), and according
to Dremova et al. (2017) a modified Hills mechanism can explain their origin.
This mechanism consists of the gravitational interaction of two supermassive
black holes that eject stars located in their central clusters. This mechanism
predicts velocities of tens of thousands of km s-1, the maximum speed of
ejection of the SRSs being determined by the mass of the secondary black hole
and the mass of the ejected star (e.g. Guillochon & Loeb, 2015).

In this work, we aim to characterise the NT particle production and associated
emission spectra of BSs produced by massive stars that propagate with extreme
velocities. In particular, we focus on massive HVSs and a putative SRS.

The paper is organised as follows. In Sec. 2, we present a multi-zone emission
model that is suitable for extreme velocity stars for which the FS is also
relevant. We present and discuss our results in Sec. 3, and finally, we
conclude with a summary of the main findings of our work in Sec. 4.

##  2 Model

Figure 1: Sketch of the model considered. The position of the CD is
represented by a black solid line, while the orange and blue regions represent
the FS and the RS, respectively. The solid lines with arrows represent
different streamlines in each shock, injected in different positions separated
by an angle $\Delta\theta$. We also show the orientation of perpendicular and
tangential vectors to the shocks in different positions, alongside
$\@vec{v}_{\rm w}$ and $\@vec{V}_{\rm ISM}=-\@vec{V}_{\star}$, and the angle
$\alpha$ between them. Adapted from del Palacio et al. (2018).

###  2.1 Scenarios studied

We define the fiducial cases to study keeping a compromise between the
potential detectability of the sources and the feasibility of finding such
objects. Regarding the luminosity of a BS produced by a massive star, the most
important parameters are those associated with the properties of the stellar
wind, and how they relate to the properties of the medium (del Palacio et al.,
2018). We focus here particularly on massive stars in the main sequence, as
this is the evolutionary stage in which they spend most of their life. These
stars produce more powerful winds (and therefore more luminous BSs) for
younger spectral types, although younger and more massive stars are less
numerous (Salpeter, 1955). In addition to stellar mass dependence, the
energetics and radiation efficiency of the BS can also increase with the
stellar spatial velocity (Martínez et al., 2021).

Observations with the Gaia satellite have detected several B-type stars with
speeds over 100 km s-1, but not a substantial quantity of O-type stars with
those velocities (Kreuzer et al., 2020). According to Marchetti et al. (2019),
the highest velocities found in HVSs are of the order of $V_{\star}\sim
10\,000$ km s-1, although estimations above 1 000 km s-1 are unreliable.
Putting everything together, we decide to study putative B0 and B1 type HVSs
with velocities between 500 and 1 000 km s-1.

Within the scenarios of interest, we set the most promissory and less common
stellar spectral type (B0) for the most conservative case, that is, a star
with the lowest spatial velocity (500 km s-1) propagating through the Galactic
disk (d). For a B0 star, plausible wind parameters are $\dot{M}_{\rm
w}=10^{-8}\leavevmode\nobreak\ $M☉ yr-1 and $v_{\rm w}=$1 500 km s-1 (Krtička,
2014; Kobulnicky et al., 2019). In addition, we consider HVSs with
$V_{\star}=1\,000$ km s-1, $\dot{M}_{\rm w}=10^{-9}\leavevmode\nobreak\ $M☉
yr-1 and $v_{\rm w}=$1 200 km s-1 propagating through three different media:
the Galactic disk, the Galactic halo (h) and a molecular cloud (mc). This is
motivated by the different densities in each of these media, easily reached by
the HVSs, and the expectation that the FS is more luminous in a denser medium
(Martínez et al., 2021). Lastly, we consider a B2 SRS with a spatial velocity
of $V_{\star}=60\,000$ km s-1, $\dot{M}_{\rm w}=10^{-8}\leavevmode\nobreak\
$M☉ yr-1 and $v_{\rm w}=$1 000 km s-1, a viable scenario according to Dremova
et al. (2017).

We summarise the characteristics of the selected scenarios in Table 1.
Henceforth, we will refer to the systems studied as <spectral type>-<velocity
of the star [in units of $10^{3}$ km s-1]>-<medium of propagation>. For
instance, B1–1–d represents a B1 star propagating at 1 000 km s-1 through the
Galactic disk.

Table 1: Parameters of the systems modelled. The values of $v_{\mathrm{w}}$ and $\dot{M}_{\rm w}$ are taken from Krtička (2014) and Kobulnicky et al. (2019), and $R_{\star}$ and $T_{\star}$ from Harmanec (1988). We assume solar abundances. |  | Scenario
---|---|---
Parameter | Symbol | B0–0.5–d | B1–1–h | B1–1–d | B1–1–mc | B2–60–d
Peculiar velocity |  $V_{\star}$ [km s-1] | 500 | 1000 | 1000 | 1000 | 60 000
Ambient density |  $n_{\mathrm{ISM}}$ [cm-3] | 10 | 0.1 | 10 | 100 | 10
Ambient mean molecular weight | $\mu_{\rm ISM}$ | 1.28 | 1.28 | 1.28 | 2.35 | 1.28
Ambient temperature [K] | TISM | 100 | 100 | 100 | 10 | 100
Spectral type |  | B0 | B1 | B1 | B1 | B2
Wind velocity |  $v_{\mathrm{w}}$ [km s-1] | 1500 | 1200 | 1200 | 1200 | 1000
Wind mass-loss rate |  $\dot{M}_{\rm w}$ [M☉ yr-1] | $10^{-8}$ | $10^{-9}$ | $10^{-9}$ | $10^{-9}$ | $10^{-10}$
Stellar radius |  $R_{\star}$ [R☉] | 5.5 | 4.8 | 4.8 | 4.8 | 4
Stellar temperature |  $T_{\star}$ [kK] | 29 | 25 | 25 | 25 | 20
Distance to BS apex |  $R_{0}$ [AU] | 792 | 1120 | 112 | 26.2 | 0.54

###  2.2 Geometry

The BS forms as the result of the collision between the stellar wind and the
ISM material. The former can be modelled as a spherical wind and the latter as
a planar wind in the reference frame of the star. The FS propagates through
the ISM and the RS propagates through the stellar wind, both separated by a
contact discontinuity (CD). The shape of the CD is defined by the condition
that the total momentum flows tangential to the shocked region, so that the
flux of mass through the CD surface is null (Wilkin, 1996; Christie et al.,
2016). Each system is then divided in four regions: free-flowing stellar wind,
unshocked ISM, shocked ISM, and shocked stellar wind, as represented in Fig.
1.

The DSA mechanism operates in strong and adiabatic shocks. For massive stars,
the RS always fulfils these conditions given that radiative cooling is not
efficient there (del Valle & Romero, 2012). The FS, on the other hand, is
radiative for typical runaway stars, but for the HVSs and SRSs considered in
this work the FS is strong and adiabatic instead (see forthcoming Sec. 2.3).
Thus, both the RS and the FS are promissory accelerators of cosmic rays and
are included in the model.

The closest position of the CD to the star, known as the stagnation point, is
located along the axis of symmetry where the total pressures of the ISM and
the stellar wind cancel each other out. The thermal pressure in both the ISM
and the wind is negligible. Thus, the ISM total pressure ($P_{\rm ISM}$) is:

$P_{\rm ISM}=\rho_{\rm ISM}\,V_{\star}^{2},$ (1)

where $\rho_{\rm ISM}$ is the density of the ISM. On the other hand, the total
pressure of the stellar wind is given in terms of the mass-loss rate, the
distance to the star, and the velocity of the wind:

$P_{\rm w}=\frac{\dot{M}_{\rm w}v_{\rm w}}{4\pi R(\theta)^{2}}.$ (2)

Matching Eqs. (1) and (2), the stagnation point is located at

$R_{0}=\left[\frac{\dot{M}_{\rm w}v_{\rm w}}{4\pi\rho_{\rm
ISM}V_{\star}^{2}}\right]^{1/2}.$ (3)

As it will be shown in the forthcoming Sec. 3, in the case of BSs from HVSs
and SRSs, NT processes are relevant even at distant regions from the apex.
Thus, a one-zone model approximation in which the emitter is considered
homogeneous is not appropriate if one aims at more detailed quantitative
predictions. Therefore, we adopt a multi-zone emission model based on the one
developed by del Palacio et al. (2018), with the major difference being the
incorporation of the FS in the model. For this, we follow an analogous
approach as the one used for the RS in del Palacio et al. (2018) with some
small modifications in the hydrodynamics.

Each shock in the BS is treated as a two-dimensional (2D) structure.
Neglecting the width of the shocked gas layers, these are co-spatial with the
CD111Fig. 1 shows a shocked gas layer of non-zero width for illustrative
purposes only.. The shape of the CD is calculated using the formulae given by
Wilkin (1996). We assume that the shocked gas flows downstream at a fixed
angle $\phi$ around the symmetry axis, dragging away NT particles. We model
the 2D structure as a sum of 1D linear emitters embedded in the 3D physical
space. For a given angle $\phi$ there are several 1D emitters, each starting
at a different angle $\theta$ with respect to the symmetry axis where
relativistic particles are injected, and consisting of multiple cells located
along the path of the 1D emitter on the CD. Particles enter in the BS and are
accelerated in the first cell of each 1D emitter, and then they move to the
following cells up to an angle $\theta_{\rm max}$ (Fig. 1). All the 1D
emitters at a certain angle $\phi$ are summed up, thereby obtaining a 1D
structure that contains all the relativistic particles along a shock at that
particular angle. Finally, we rotate this 1D structure made of all the 1D
emitters with the same $\phi$ value around the symmetry axis to get the full
2D structure of the BS. At last, we note that in all the scenarios we study,
the stellar velocity is significantly below the speed of light. Therefore,
relativistic effects in the hydrodynamics, as well as in the calculation of
the NT particle distribution at each cell and the radiation they produce
(e.g., Doppler boosting), can be neglected.

###  2.3 Hydrodynamics

We introduce here the semi-analytical hydrodynamical model used to
characterise the properties of the shocked gas in the BS. In del Palacio et
al. (2018), Rankine-Hugoniot jump conditions were assumed, which near the apex
yields correct values up to first order. Here, we adopt a more consistent
approach that is suitable to model both the RS and the FS up to distant
regions from the BS apex. Firstly, we assume that the total pressure at
$R_{0}$ is the ambient total pressure. Secondly, we consider Bernoulli’s
equation across the shock:

$\frac{1}{2}v_{\rm u}^{2}+\left(\frac{\gamma_{\rm ad}}{\gamma_{\rm
ad}-1}\right)\frac{P_{\rm th,u}(\theta)}{\rho_{\rm
u}(\theta)}=\frac{1}{2}v(\theta)^{2}+\left(\frac{\gamma_{\rm ad}}{\gamma_{\rm
ad}-1}\right)\frac{P_{\rm th}(\theta)}{\rho(\theta)},$ (4)

where we labelled the upstream region with the subscript u. Given that
$R_{0}\gg R_{\star}$ for the scenarios considered here (Table 1), we adopt
$v_{\rm w}=v_{\infty}$. Then, $v_{\rm u}$ is equal to $v_{\infty}$ in the RS,
and equal to $V_{\star}$ in the FS. Moreover, we neglect the thermal pressure
in the upstream region when compared to the ram (kinetic) pressure, since
$P_{\rm th,u}\ll P_{\rm ram,u}\propto v_{\rm u}^{2}$. The incoming stellar
wind and ISM impact perpendicularly to the BS apex, where the fluid halts.
Hereafter, we represent quantities at the stagnation point ($\theta=0$) with a
subscript 0. Using Eq. (4), we determine the density at the stagnation point
as:

$\rho_{0}=\left(\frac{\gamma_{\rm ad}}{\gamma_{\rm ad}-1}\right)\frac{2P_{{\rm
th}0}}{v_{\rm u}^{2}}=5\,\frac{P_{{\rm th}0}}{v_{\rm u}^{2}}=5\rho_{\rm u},$
(5)

assuming that the fluid behaves like an ideal gas with adiabatic coefficient
$\gamma_{\rm ad}=5/3$.

We assume that the shocked gas moves parallel to the CD. In the regions where
the shocked fluid is subsonic, the flux of momentum that crosses the shock
perpendicularly heats the downstream material, converting the total pressure
in the upstream into thermal pressure in the downstream. We can thus calculate
the total pressure at each point as:

$P(\theta)=P_{\rm th}(\theta)=\rho_{\rm u}(\theta)v_{\rm u}v_{{\rm
u}\perp}(\theta).$ (6)

On the other hand, if the fluid becomes supersonic at an angle $\theta_{\rm
c}$, we consider that only the momentum density component perpendicular to the
BS that crosses the BS also perpendicularly is converted to thermal energy.
This two-region approach requires adopting a soft transition of the
thermodynamic quantities between both regimes. We thus adopt a prescription
for the total pressure for $\theta>\theta_{\rm c}$ of the form:

$P(\theta)=P_{\rm th}(\theta)=\rho_{\rm u}(\theta)\left(\frac{v_{\rm
u}}{v_{{\rm u}\perp}(\theta_{\rm c})}\right)v_{{\rm u}\perp}(\theta)^{2}.$ (7)

We can then calculate the density using the politropic relation:

$\rho(\theta)=\rho_{0}\left(\frac{P(\theta)}{P_{0}}\right)^{1/\gamma_{\rm
ad}}.$ (8)

However, at large values of $\theta$ this underestimates the fluid density, so
we impose the condition $\rho(\theta)\geq\rho_{\rm min}(\theta)$, where
$\rho_{\rm min}$ is defined as

$\rho_{\rm min}(\theta)=\begin{dcases}\rho_{\rm u}\quad{\rm FS}\\\
\left(\frac{4\pi}{\Omega(\pi-\theta)}-1\right)\rho_{\rm u}\quad{\rm
RS},\end{dcases}$ (9)

with $\Omega(x)=2\pi\left(1-\cos{(x)}\right)$, which takes into account the
fraction of stellar wind accumulated in the RS up to an angle $\theta$. To
derive these expressions we have assumed that the shocked flow behaves as only
one stream line which at each location adapts to the impact of the incoming
upstream material, having homogeneous conditions in the direction
perpendicular to the CD. This approximation is valid in the subsonic region,
although it weakens in the supersonic region as the shocked flow is not
causally connected along the CD.

The assumption of a laminar flow requires the suppression of dynamical
instabilities, such as Rayleigh-Taylor (triggered by density differences
across the CD) and Kelvin-Helmholtz instabilities (triggered by tangential
velocity differences across the CD). Such instabilities can arise in stellar
BSs, especially when the ambient medium is dense (e.g. Comeron & Kaper, 1998;
Meyer et al., 2016), although a high stellar velocity inhibits their
appearance up to $\theta\gtrsim 135\degree$ (Comeron & Kaper, 1998; Christie
et al., 2016). Neglecting instabilities is further justified in the presence
of adiabatic shocks (e.g. Falceta-Gonçalves & Abraham, 2012), as it is the
case in the context of HVSs and SRSs.

A shock is considered adiabatic when the gas escapes from the shock region
before it cools significantly. This condition can be expressed as $t_{\rm
th}/t_{\rm conv}\gtrsim 1$, with

$t_{\rm th}(\theta)=\frac{k_{\rm B}\,\mu_{\rm u}\,m_{\rm
p}\,T}{\zeta(\theta)\,\rho_{\rm u}\,\lambda(T)},$ (10)

where the function $\lambda(T)$ depends on the temperature of the shock
(Myasnikov et al., 1998), and $m_{\rm p}$ is the proton mass. For the
convection timescale we use the expression given in del Palacio et al. (2018),
$t_{\mathrm{conv}}=R(\theta)/v_{\parallel}(\theta)$. In Fig. 2 we show the
logarithm of the ratio $t_{\rm th}/t_{\rm conv}$ as a function of the angle
$\theta$ for the FSs. This ratio increases with $\theta$ up to $\theta\sim
45\degree$ in the HVSs and up to $\theta\sim 70\degree$ in the SRS as the
fluid accelerates. After this, the ratio slowly decreases as $t_{\rm
conv}\propto R(\theta)$ increases. As a result, the FS is adiabatic up to
$\theta\sim 160\degree$ in the systems B1–1–d and B1–1–mc, and up to
$\theta\sim 135\degree$ in the system B0–05–d (see Fig. 2). Moreover, the FS
of the systems B1–1–h and B2–60–d, as well as the RSs in all cases studied,
always fulfil the adiabaticity condition.

Figure 2: Logarithm of the ratio $t_{\rm th}/t_{\rm conv}$ for the different
FSs studied in this work. For reference, we plot a black dotted horizontal
line at zero (where $t_{\mathrm{th}}=t_{\mathrm{conv}}$); curves above this
line corresponds to adiabatic shocks.

Lastly, the magnetic field in the subsonic regime
($\theta<\theta_{\mathrm{c}}$) is obtained imposing that, at each position,
its pressure is a fraction $\eta_{B}$ of the thermal pressure of the shocked
region:

$B(\theta)=\left(\eta_{B}\,8\pi P_{\rm th}(\theta)\right)^{1/2}.$ (11)

In the supersonic regime ($\theta>\theta_{\mathrm{c}}$), we assume that the
magnetic field remains frozen to the plasma, and that it is tangent to the
shock surface. The later assumption is motivated by the fact that the magnetic
field component perpendicular to the shock normal is amplified by adiabatic
compression, and thus in general it should be the dominant component in both
shocks (RS and FS).

We can then obtain the magnetic field as:

$B(\theta)=B(\theta_{\mathrm{c}})\,\left(\frac{\rho(\theta)}{\rho(\theta_{\mathrm{c}})}\frac{v(\theta_{\mathrm{c}})}{v(\theta)}\right)^{1/2}.$
(12)

In Fig. 3 we show the dependence of the thermodynamic quantities with $\theta$
for both shocks. In the apex the conversion of kinetic energy to internal
energy is maximised as the shock is perpendicular. The tangential velocity
increases monotonically with $\theta$, going from zero in the stagnation point
to $v_{\rm t}\sim v_{\infty}$ in the RS and $v_{\rm t}\sim V_{\star}$ in the
FS. The pressure slowly decays with $\theta$ in the FS, while it drops more
abruptly in the RS because $P_{\rm RS}\propto\rho_{\rm w}\propto R^{-2}$; as a
consequence, the other quantities that depend on $P$ also decay gradually. We
highlight that the magnetic field decreases slowly in the FS, which favours
synchrotron emission up to large values of $\theta$. Lastly, we note that this
hydrodynamic model yields densities along the shocks that are slightly higher
than the ones obtained assuming Rankine-Hugoniot jump conditions. The
discrepancy is a factor $\sim 1.5$ for angles $\theta\lesssim 60\degree$, and
a factor $\sim 2$ in the distant regions with $\theta>60\degree$ (see Fig. 8).
The reason for this is that, as explained above, different regions of the
shocked layer affect each other making the local hydrodynamical conditions
depart from the Rankine-Hugoniot ones.

Figure 3: _Left panel:_ Thermodynamic variables in the RS as a function of the
position angle along the shock. We also give as a reference the linear
distance along the shock, $\Delta s=\int_{0}^{\theta}\mathrm{d}l$. The sub-
index 0 refers to values in the apex. _Right panel:_ The same but for the FS.

###  2.4 Non-thermal particles

Relativistic particles can be accelerated via DSA in hypersonic and adiabatic
shocks, such as the ones present in the BS (Sect. 2.3). Additionally, both
electrons and protons could be accelerated via shock drift acceleration (SDA)
in the RS, as the stellar magnetic field lines are expected to be parallel to
the shock surface (Marcowith et al., 2016). Nevertheless, acquiring
relativistic energies through SDA requires multiple interactions with the
shock front, similarly to DSA (Matthews et al., 2020). Given that both
mechanisms lead to the injection of a power-law distribution of particles and
we treat the acceleration details phenomenologically, we refer in what follows
only to acceleration via DSA although SDA might be also involved.

The energy distribution of the injected particles at the i-th cell is assumed
to be

$Q(E,\theta_{i})=Q_{0}E^{-p}\exp(-E/E_{\rm max}),$ (13)

where $Q_{0}$ is a normalisation factor, $p$ is the spectral index, and
$E_{\rm max}$ is the cut-off energy, all dependent on $\theta_{i}$.

The normalisation constant $Q_{0}$ is set by the condition $\int
EQ(E,\theta_{i})\,{\rm d}E=\Delta L_{\rm NT}(\theta_{i})$, being $\Delta
L_{\rm NT}(\theta_{i})$ the power available to accelerate NT particles at each
position. That is,

$\Delta L_{\rm NT}(\theta_{i})=f_{\rm NT}\,\Delta L_{\perp}(\theta_{i})=f_{\rm
NT}\,S_{\rm E}(\theta_{i})\,A_{\perp}(\theta_{i}),$ (14)

where $S_{\rm E}(\theta_{i})$ is the energy flux per unit volume of the
corresponding fluid, and $A_{\perp}(\theta_{i})$ is the area of the cell
surface projected perpendicular to $\@vec{v}_{\rm u}$. This area is calculated
as $A_{\perp}(\theta_{i})=R(\theta_{i})\,\sin(\theta_{i})\,\Delta
l(\theta_{i})\,\sin{(\alpha_{i})}\,\Delta\phi$, with $\Delta l(\theta_{i})$
being the length of the cell and $\alpha_{i}$ the angle between $\@vec{v_{\rm
u}}$ and the tangent to the shock, $\@vec{T}$ (Fig. 1). The parameter $f_{\rm
NT}$ is defined as the fraction of the power injected to the BS that goes to
NT particles. We adopt $f_{\rm NT}=0.1$ and assume that 95% of this power goes
to protons, while the remaining 5% goes to electrons. Finally, the energy flux
per unit volume in the subsonic regime is $S_{\rm E}=0.5\,\rho_{\rm u}\,v_{\rm
u}^{3}$, whereas in the supersonic regime only the perpendicular velocity
component is converted into thermal energy and therefore $S_{\rm
E}=0.5\,\rho_{\rm u}\,v_{\rm u}^{2}\,v_{\rm u\perp}$. Contrary to what happens
in the RS, the injection of energy in the FS is relevant even at large values
of $\theta$, where the flow is more susceptible to develop instabilities. That
could potentially increase the area of the shock, thus enhancing the injected
power and, consequently, the emitted luminosity (e.g. de la Cita et al.,
2017). Nevertheless, the variability induced by this effect is not expected to
dominate the average luminosities predicted in our model.

We can determine $p$ in terms of the compression factor $\zeta$ as (e.g.
Caprioli & Spitkovsky, 2014):

$p(\theta_{i})=\frac{\zeta(\theta_{i})+2}{\zeta(\theta_{i})-1},\quad\zeta(\theta_{i})=\left(\frac{\gamma_{\rm
ad}+1}{\gamma_{\rm ad}}\right)\frac{M(\theta_{i})^{2}}{M(\theta_{i})^{2}+2},$
(15)

where $M$ is the Mach number. Finally, $E_{\rm max}$ is obtained by equating
the acceleration and energy loss timescales.

The steady-state particle distribution at the injection cell is:

$N_{0}(E,\theta_{i})\approx Q(E,\theta_{i})\times{\rm min}(t_{\rm cell},t_{\rm
cool}),$ (16)

being $t_{\rm cell}$ the cell convection time and $t_{\rm cool}$ the cooling
timescale. In the BS, electrons cool mainly by synchrotron and inverse Compton
(IC) interactions, the latter with both the stellar UV (IC⋆) and dust IR
(ICIR) radiation fields (Martínez et al., 2021). Adiabatic losses can also be
relevant, whereas relativistic Bremsstrahlung losses are negligible. Protons
cool by proton-proton inelastic collisions, a rather minor effect, and
adiabatic losses.

The NT particles are confined within the shock and so they are dragged by the
fluid. This occurs because the particles gyro-radii, $r_{\rm
g}(E,\theta)\propto E/B(\theta)$, even for $E\sim E_{\mathrm{max}}$ are much
smaller than the shocked layer width (the shock typical scale) $H(\theta)$.
The latter is calculated considering mass conservation across the shock

$H(\theta)=\frac{\int_{0}^{\theta}\rho_{\rm u}(\theta^{\prime})\,v_{\rm
u}\,A_{\perp}(\theta^{\prime})\,{\rm
d}\theta^{\prime}}{2\pi\,R(\theta)\,\sin{(\theta)}\,\rho(\theta)\,v_{\parallel}(\theta)},$
(17)

yielding a typical value of $H(\theta)\sim 0.3R(\theta)$ for $\theta<\pi/2$
(Christie et al., 2016). Under these conditions, let us consider that there is
a certain number of NT particles with energy $E$ in the $i$-th cell. By the
time they reach the ($i$+1)-th cell, their energy will be $E^{\prime}\leq E$,
and the size and convection velocity of the cell will also be different.
Nonetheless, in the steady-state, the flux of particles in the position and
energy space must be conserved. Considering this, we obtain the evolution of
the particle energy distribution along a linear emitter222We note that the
inclusion of $t_{\mathrm{cell}}$ in this expression is a small correction to
the one used in del Palacio et al. (2018).:

$N(E^{\prime},i+1)=N(E,i)\,\frac{|\dot{E}(E,i+1)|}{|\dot{E}(E^{\prime},i+1)|}\,\frac{t_{\mathrm{cell}}(i+1)}{t_{\mathrm{cell}}(i)},$
(18)

with $\left|\dot{E}(E,i)\right|=E/t_{\rm cool}(E,i)$ the cooling rate for
particles of energy $E$ at the position $\theta_{i}$ and $t_{\rm cell}(i)$ the
convection time of the i-th cell. Finally, the energy $E^{\prime}$ is given by
the condition $t_{\rm cell}=\int_{E}^{E^{\prime}}\dot{E}(\tilde{E},i)\,{\rm
d}\tilde{E}$.

We use the formulae given by Khangulyan et al. (2014) to calculate the
isotropic IC cooling timescales, given that the electron distribution is
isotropic at each position and the fluid is non-relativistic. These
expressions take into account the Klein-Nishina (K-N) cross section for the
interaction at high energies. For the case of the stellar radiation field, we
consider the star as a black body emitter with temperature $T_{\star}$ and a
dilution factor of the photon field
$\kappa_{\star}=\left[R_{\star}/(2R(\theta))\right]^{2}$. For the case of the
IR photon field produced by the dust, we assume isotropy within the NT
emitter, and we approximate its spectrum with a Planck law of temperature
$T_{\rm IR}=100$ K. The corresponding dilution factor is $\kappa_{\rm
IR}(\theta)=U_{\rm IR}(\theta)/U_{\rm BB}$, where $U_{\rm BB}$ is the energy
density of the radiation of a black body with temperature $T_{\mathrm{IR}}$.
Since $U_{\rm IR}\approx L_{\rm IR}/\left(4\pi R(\theta)^{2}c\right)$
(considering that the dust is concentrated in a thin shell surrounding the BS)
and $U_{\rm BB}=4\sigma T_{\rm IR}^{4}/c$, we obtain:

$\kappa_{\rm IR}=\frac{L_{\rm IR}}{16\pi\sigma T_{\rm IR}^{4}R(\theta)^{2}}.$
(19)

##  3 Results

First, we compute the timescales and particle energy distributions for the
scenarios studied. Then, we use a one-zone model to estimate the luminosity
scaling with the relevant parameters of the systems. Finally, we present the
spectral energy distributions (SEDs), and discuss the detectability of each
system.

The magnetic field in the shocked region can be generated by adiabatic
compression of the ISM (star) magnetic field in the FS (RS) and/or be
amplified by the action of cosmic rays. Adopting $\eta_{\rm B}=0.1$ in Eq.
(11) yields values of $B$ consistent with a ratio between NT energy density
and magnetic energy density of $U_{\rm NT}/U_{\rm B}\gtrsim 1$. Thus, in both
shocks the magnetic field could be the result of amplification by cosmic rays
(Bell, 2004). Alternatively, in the RS the magnetic field could come from
adiabatic compression of the stellar magnetic field. Under these conditions,
and adopting an Alfvén radius $r_{\rm A}\sim R_{\star}$, the stellar magnetic
field in the stellar surface is $B_{\star}\sim
0.25\,B(\theta)\,\left(R(\theta)/R_{\star}\right)\,\left(v_{\infty}/v_{\rm
rot}\right)$, with the star rotation speed being $v_{\rm rot}\sim
0.1\,v_{\infty}$ (del Palacio et al., 2018). We obtain $B_{\star}\sim 25$ G
for the B0 star, $B_{\star}\sim 10$ G for the B1 stars, and $B_{\star}\sim 3$
G for the B2 star. Given that these are plausible values (Parkin et al.,
2014), adopting $\eta_{\rm B}=0.1$ is a reasonable assumption. Additionally,
we consider an equipartition scenario ($\eta_{B}=1$) to set an upper limit to
the predicted radio fluxes. At last, given that the majority of power is
injected within $\theta<170\degree$, and that instabilities can be significant
for large values of $\theta$, we fix $\theta_{\rm max}=170\degree$.

Table 2: Power injected in NT particles (electrons and protons) in each system, assuming $f_{\mathrm{NT}}=0.1$ and $\eta_{\rm B}=1$ (see text in Sec. 2.4 for details). |  $L_{\rm NT}$ [erg s-1]
---|---
Scenario | RS | FS | Total
B0–05–d | $5.3\times 10^{32}$ | $6.9\times 10^{32}$ | $1.2\times 10^{33}$
B1–1–h | $3.4\times 10^{31}$ | $1.1\times 10^{32}$ | $1.4\times 10^{32}$
B1–1–d | $3.4\times 10^{31}$ | $1.1\times 10^{32}$ | $1.4\times 10^{32}$
B1–1–mc | $3.4\times 10^{31}$ | $1.1\times 10^{32}$ | $1.4\times 10^{32}$
B2–60–d | $2.3\times 10^{30}$ | $5.5\times 10^{32}$ | $5.5\times 10^{32}$

###  3.1 Relativistic particle population

In Table 2 we show the power injected in NT particles (both electrons and
protons) for each scenario. This quantity increases with younger spectral
types and the velocity of the star, as expected from Eq. (14). Considering
that there are at most tens of thousands of HVSs in the Milky Way (Marchetti
et al., 2019), the total luminosity injected in NT particles by these sources
in the Galaxy is likely below $10^{38}$ erg s-1. Then, HVSs do not
significantly contribute to the Galactic cosmic ray population, which has a
much larger (by a factor $\sim 10^{3}$) contribution from supernova remnants.

####  3.1.1 Forward shock

An example of the timescales considered is shown in Figs. 4 and 5 for the FS
of the system B1–1–d. For electrons near the apex, convection (escape)
dominates up to $E_{\rm e}\leq 4$ GeV. In the range $4\leavevmode\nobreak\
{\rm GeV}\leq E_{\rm e}\lesssim 160$ GeV, IC losses with the stellar photon
field are dominant, also versus IR IC losses; these interactions occur in the
Thomson regime for $E_{\rm e}\lesssim 30$ GeV, and in the K-N regime for
$E_{\rm e}\gtrsim 30$ GeV. At last, diffusion (escape) losses dominate above
160 GeV, and electrons reach energies of $E_{\rm e,max}\sim 500$ GeV. In a
scenario with a lower ISM density (system B1–1–h), the stagnation point is
further from the star, and so the density of stellar photons in the BS is
smaller. As a consequence, IC timescales are larger, and convection losses are
dominant up to $E_{\rm e}\lesssim 100$ GeV, and diffusion losses are dominant
in the range $100\leavevmode\nobreak\ {\rm GeV}\lesssim E_{\rm e}\lesssim
300\leavevmode\nobreak\ {\rm GeV}\sim E_{\rm e,max}$. Similarly, $R_{0}$ also
increases for main-sequence stars with younger spectral types, as
$\dot{M}_{\rm w}$ and $v_{\rm w}$ increase. Finally, for a SRS (B2–60–d),
$R_{0}$ is significantly closer to the star ($R_{0}\approx
30\leavevmode\nobreak\ R_{\star}$). Then, IC⋆ interactions are dominant for
$3\leavevmode\nobreak\ {\rm GeV}\leq E_{\rm e}\lesssim 600\leavevmode\nobreak\
{\rm GeV}$, and ICIR interactions are dominant in the range
$600\leavevmode\nobreak\ {\rm GeV}\lesssim E_{\rm e}\lesssim
2\leavevmode\nobreak\ {\rm TeV}$. Above $E_{\rm e}\geq 2\leavevmode\nobreak\
{\rm TeV}$ and up to $E_{\rm max}$ diffusion dominates. Since
$B(\theta)\propto P(\theta)\propto V_{\star}^{2}\gg 1\,000$ km s-1, the
acceleration of NT particles is very efficient, yielding maximum energies of
$E_{\rm e,max}\sim 5$ TeV for electrons and $E_{\rm p,max}\sim 10$ TeV for
protons.

In distant regions from the apex, the stellar photon field is more diluted and
therefore IC losses are less important. Additionally, the convection timescale
shortens as $v_{\parallel}$ increases. As a consequence, convection losses are
completely dominant in the FS of all systems for angles $\theta\gtrsim
100\degree$. Moreover, particle acceleration is less efficient in distant
regions since both $v_{\perp}$ and $B$ decrease with $\theta$ (Fig. 3). Then,
$E_{\rm max}$ diminishes with $\theta$ for both electrons and protons.

In contrast, protons are convected away from the FS without radiating a
significant fraction of their energy. For the HVSs, close to the apex
convection losses are the dominant process for protons with $E_{\rm p}<160$
GeV. In the case of the SRS (scenario B2–60–d), convection dominates up to
$E_{\rm p}\gtrsim 1$ TeV333 We note that after convecting away from the BS,
protons diffuse in the surrounding medium and, in dense environments, can
produce significant $\gamma$-ray radiation via proton-proton collisions (del
Valle et al., 2015).. Above the energies mentioned, and up to $E_{\rm p,max}$,
diffusion losses are dominant.

Protons reach their maximum energy at the apex ($\theta=0$). We estimate the
scaling of $E_{\rm p,max,0}$ with the system parameters by matching $t_{\rm
ac,0}(E_{\rm p})=t_{\rm diff,0}(E_{\rm p})$. Considering diffusion in the Bohm
regime, these timescales are:

$\displaystyle t_{\rm ac,0}(E_{\rm p})$ $\displaystyle=\frac{2\pi
c}{q}\frac{E_{\rm p}}{B_{0}\,v_{{\perp}0}}=\frac{2\pi c}{q}\frac{E_{\rm
p}}{B_{0}\,v_{\rm u}},$ (20) $\displaystyle t_{\rm diff,0}(E_{\rm p})$
$\displaystyle=\frac{R_{0}^{2}}{2D_{\rm
B,0}}=\frac{3q}{2c}\frac{R_{0}^{2}\,B_{0}}{E_{\rm p}},$ (21)

where $c$, $q$ and $D_{\rm B}$ are the speed of light, proton charge and
diffusion coefficient, respectively. Finally, using Eqs. (6) and (11), we get
$E_{\rm max,p,0}\propto R_{0}\,B_{0}\,V_{\star}^{0.5}\propto\dot{M}_{\rm
w}^{0.5}\,v_{\rm w}^{0.5}\,V_{\star}^{0.5}$. This maximum energy is $E_{\rm
max,p}\sim 1$ TeV for the HVSs and $E_{\rm max,p}\sim 10$ TeV for the SRS.

When convection dominates, particles move along the BS with an energy
distribution that keeps the same spectral index as the injected distribution.
On the other hand, when IC in the Thomson regime or diffusion dominates, the
particle energy distribution is softened. In Fig. 6 we show the particle
energy distribution of electrons and protons in different regions of the RS
and the FS for the system B1–1–d.

####  3.1.2 Reverse shock

Characteristic timescales near the apex for the RS are similar to the ones of
the FS of the HVSs. Nonetheless, the behaviour is different for the SRS.
Despite the acceleration timescale decreases as the magnetic field is
stronger, the IC cooling timescale decreases more drastically, yielding
maximum electron energies of $E_{\rm e,max}\sim 10$ GeV in the RS (the RS
shock velocity is much lower than in the FS). Moreover, we highlight that IC
cooling is still relevant up to $\theta\sim 160\degree$ for the RS of this
system.

Finally, matching Eq. (20) and (21) for the RS we find for protons that
$E_{\rm max,p}(0)\propto R_{0}\,B_{0}\,V_{\rm w}^{0.5}\propto\dot{M}_{\rm
w}^{0.5}\,v_{\rm w}$. In this case, the maximum energy is $E_{\rm
max,p,0}\gtrsim 1$ TeV for the HVSs and $E_{\rm max,p,0}\sim 200$ GeV for the
SRS.

Figure 4: Cooling times for three different positions of electrons for the FS
of the system B1-1-d. The IC cooling timescale of electrons takes into account
both the stellar and dust photon fields. Figure 5: Cooling times for protons
of the system B1–1–d near the apex.

Figure 6: Particle energy distribution for electrons (top panels) and protons
(bottom panels). The left panels are for the RS and the right panels are for
the FS. The colour scale represents five different regions of the emitter that
correspond to intervals of length $\Delta\theta=\theta_{\rm max}/5=34\degree$.
The black dashed line corresponds to the total particle energy distribution.

###  3.2 Analytical estimates on emissivity scaling

A one-zone approximation is suitable to obtain order-of-magnitude estimates of
the radiative outputs, with the advantage that the dependencies of the
emission with respect to different system parameters become explicit (e.g. del
Palacio et al., 2018). We therefore apply this formalism to derive the scaling
of the BS luminosity with the relevant system parameters.

We consider that the BS has an effective surface of size
$R_{\mathrm{eff}}=a\,R_{0}$. This is mostly relevant for the FS, for which
$a\sim 5$ (as explained below), whereas for the RS it is $a\sim 1$. If $t_{\rm
rad}$ is the cooling timescale of the dominant NT mechanism in a certain
energy range, we estimate the NT power emitted in that range by each shock as

$L_{\rm rad}\sim L_{\rm NT}\,\left(\frac{t_{\rm conv}}{t_{\rm
rad}}\right)\propto\begin{dcases}f(a)\,\dot{M}_{\rm w}\,v_{\rm
w}^{2}\,\left(\frac{t_{\rm conv}}{t_{\rm rad}}\right)&{\rm RS}\\\ \rho_{\rm
ISM}\,R_{\rm eff}^{2}\,V_{\star}^{3}\left(\frac{t_{\rm conv}}{t_{\rm
rad}}\right)&{\rm FS},\end{dcases}$ (22)

where $f(a)$ is an order unity function that tends to $f(a)\approx 0.6$ as $a$
increases. Qualitatively, the ratio $t_{\rm conv}/t_{\rm syn}$ determines the
NT luminosity emitted in the radio band. On the other hand, $\gamma$-ray
emission depends on the ratio $t_{\rm conv}/t_{\rm IC\star}$, as IC process
with the stellar radiation field is the dominant NT process for those
energies. We mote that these timescale-ratio dependencies of the radiation
luminosity are strictly valid for a dominant $t_{\rm conv}$.

The convection (escape) timescale is roughly determined by the effective
radius and the sound speed in the shocked gas:

$t_{\rm conv}\sim\frac{R_{\rm eff}}{c_{\rm
s}}\propto\begin{dcases}a\,\dot{M}_{\rm w}^{0.5}\,v_{\rm w}^{-0.5}\,n_{\rm
ISM}^{-0.5}\,\mu_{\rm ISM}^{-0.5}\,V_{\star}^{-1}&{\rm RS}\\\ \\\
a\,\dot{M}_{\rm w}^{0.5}\,v_{\rm w}^{0.5}\,n_{\rm ISM}^{-0.5}\,\mu_{\rm
ISM}^{-0.5}\,V_{\star}^{-2}&{\rm FS}.\end{dcases}$ (23)

On the other hand, the synchrotron cooling timescale is (e.g. Blumenthal &
Gould, 1970)

$t_{\rm syn}\propto B^{-2}\propto\begin{dcases}\left(\dot{M}_{\rm w}\,v_{\rm
w}\,R_{\rm eff}^{-2}\right)^{-1}\propto a^{2}\,n_{\rm ISM}^{-1}\,\mu_{\rm
ISM}^{-1}\,V_{\star}^{-2}&{\rm RS}\\\ \left(\rho_{\rm
ISM}\,V_{\star}^{2}\right)^{-1}\propto n_{\rm ISM}^{-1}\,\mu_{\rm
ISM}^{-1}\,V_{\star}^{-2}&{\rm FS}.\end{dcases}$ (24)

Thus, from Eqs. (22)–(24) we obtain

$L_{\rm syn}\sim L_{\rm NT}\,\left(\frac{t_{\rm conv}}{t_{\rm
syn}}\right)\propto\begin{dcases}f(a)\,a^{-1}\,\dot{M}_{\rm w}^{1.5}\,v_{\rm
w}^{1.5}\,n_{\rm ISM}^{0.5}\,\mu_{\rm ISM}^{0.5}\,V_{\star}&{\rm RS}\\\ \\\
a^{3}\,\dot{M}_{\rm w}^{1.5}\,v_{\rm w}^{1.5}\,n_{\rm ISM}^{0.5}\,\mu_{\rm
ISM}^{0.5}\,V_{\star}&{\rm FS}.\end{dcases}$ (25)

Therefore, synchrotron emission depends mostly on the stellar parameters:
massive stars are promissory radio emitters, as they have fast, powerful
winds. To a lesser extent, synchrotron luminosity increases with the speed of
the star and the medium density. Finally, we note that the FS luminosity has a
stronger dependency on the effective radius of the BS than the RS luminosity.
This supports the need to incorporate the factor $a$ when using one-zone
models to estimate the luminosity of BSs in systems for which the FS
contribution can be significant or even dominant. A value of $a\sim 5$ yields
luminosities that match within a factor two or three with those obtained using
a more precise multi-zone model.

The IC⋆ cooling timescale depends on the energy density of the stellar photon
field: $t_{{\rm IC}\star}^{-1}\propto U_{\star}\propto L_{\star}R_{\rm
eff}^{-2}$. Assuming $L_{\star}\propto\dot{M}_{\rm w}^{0,5}v_{\rm w}^{0.5}$
(Kobulnicky et al., 2017), and using Eqs. (22) and (23), we estimate

$L_{\rm IC_{\star}}\sim L_{\rm NT}\,\left(\frac{t_{\rm conv}}{t_{\rm
IC_{\star}}}\right)\propto\begin{dcases}f(a)\,a^{-1}\,\dot{M}_{\rm w}\,v_{\rm
w}\,n_{\rm ISM}^{0.5}\,\mu_{\rm ISM}^{0.5}\,V_{\star}&{\rm RS}\\\ \\\
a\,\dot{M}_{\rm w}\,v_{\rm w}\,n_{\rm ISM}^{0.5}\,\mu_{\rm
ISM}^{0.5}\,V_{\star}&{\rm FS}.\end{dcases}$ (26)

As a consequence, high energy emission depends mostly on the mass-loss rate
and the velocities of the star and the wind. Massive stars moving with high
velocities in a dense medium are the most promising $\gamma$-ray sources.
However, as shown in the following section, even for these sources the
expected fluxes are significantly below the sensitivity threshold of current
$\gamma$-ray observatories.

###  3.3 Spectral energy distribution

In Fig. 7 we show the SEDs obtained for all the systems studied. As expected,
the NT spectrum is dominated by synchrotron emission in the radio band and IC⋆
emission in $\gamma$-rays. In the X-ray band, the IC⋆ component is usually
dominant, although a non-negligible contribution from synchrotron radiation
can also be expected for high magnetic fields, which can even be dominant in
SRSs. On the other hand, relativistic Bremsstrahlung and hadronic NT emission
are not relevant.

In all scenarios, the NT radiation of the FS is brighter than that of the RS.
In addition, the system B0–05–d is the most luminous, in agreement with the
discussion in Sec. 3.2, where we showed that synchrotron emission depends
mostly on the stellar wind parameters. On the other hand, the system B1–1–mc
has the most luminous BS among systems with
$V_{\star}=1\,000\leavevmode\nobreak\ $km s-1. Thus, as expected from Sec.
3.2, a denser medium favours NT emission.

In Table 3 we show the fluxes predicted at $\nu=1.4$ GHz, assuming a distance
to the star of 1 kpc. We conclude that a system with the characteristics of
B0–05–d could be detected by the new generation of radio interferometers, such
as the SKA (Cassano et al., 2018) and the ngVLA (McKinnon et al., 2019), if it
is at a distance $<3$ kpc. On the other hand, assuming $\eta_{B}=1$ increases
the fluxes by a factor of five, giving room for detection of slightly less
favourable sources.

The SRS is the most luminous X-ray source among the systems considered. Given
that electrons are accelerated up to energies $E_{\rm max}\sim 10$ TeV, the
synchrotron spectrum reaches energies of $\epsilon\lesssim 1$ MeV in the FS.
Nevertheless, we predict luminosities of $L_{\rm X}\sim 10^{28}$ erg s-1
between 1 keV and 10 keV, undetectable by current instruments unless at a
highly unlikely short distance. However, we highlight that by considering
$\eta_{\rm B}=1$ and that if a 50% of $L_{\rm NT}$ went to electrons, this
luminosity could increase by a factor of $\sim 200$, making it detectable by
Chandra or future instruments like Lynx, if the system is at a distance
$\lesssim 1$ kpc.

Finally, we note that IC emission reaches energies of $\epsilon\gtrsim 10$ TeV
in the FS of the system B2–60–d, for which IC with both the stellar and the IR
field occur in the K-N regime. Nonetheless, the predicted $\gamma$-ray
radiation is undetectable with present or forthcoming instrumentation in all
systems. The fluxes predicted are at least two orders of magnitude below the
detection threshold of CTA and Fermi-LAT, even assuming distances of 1 kpc to
the source (Bruel et al., 2018; Maier, 2019).

Figure 7: SEDs for the systems considering $f_{\mathrm{NT}}=0.1$ and $\eta_{B}=0.1$. Solid lines correspond to the RS, while dotted lines correspond to the FS. Table 3: Energy fluxes predicted at $\nu=1.4$ GHz assuming a distance of 1 kpc to the stars. | S${}_{\mathrm{1.4\,GHz}}$ [$\mu$Jy]
---|---
| $\eta_{B}=0.1$ | $\eta_{B}=1$
Scenario | RS | FS | Total | RS | FS | Total
B0–05–d | 24.3 | 317.4 | 341.7 | 132.5 | 1784 | 1916.5
B1–1–h | 0.3 | 3.0 | 3.3 | 1.4 | 16.6 | 18.0
B1–1–d | 0.9 | 10.9 | 11.8 | 4.6 | 56.8 | 61.4
B1–1–mc | 1.9 | 24.0 | 25.9 | 10.2 | 124.6 | 134.8
B2–60–d | 0.3 | 1.6 | 1.9 | 1.4 | 8.6 | 10.0

##  4 Conclusions

We investigated the shocks produced by massive stars that move with hypersonic
or semi-relativistic velocities with respect to their surrounding medium. We
introduced a refined model for calculating the emission from stellar BSs,
especially relevant for systems moving with very high velocities
($V_{\star}>300$ km s-1). Our results show that in these systems both the RS
and the FS are adiabatic and hypersonic, making them promising NT particle
accelerators, potentially contributing with a $\sim 0.1$% component to the
galactic cosmic rays at energies $\lesssim 1$ TeV. Protons and electrons are
accelerated in these systems up to energies $\gtrsim 1$ TeV. Nonetheless, the
detection of NT emission associated with these relativistic particles is
likely to remain elusive, as the predicted fluxes are too faint for current
observatories. We suggest that, in the near future, the most promising
observational breakthrough could be achieved by the next generation of
interferometers operating at low radio frequencies, which can potentially
detect the NT emission produced by BSs of early-type HVSs. We also note that
under rather optimistic conditions for the production of leptonic radiation
(i.e. $\eta_{\rm B}=1$ and with a $\sim 5$% of the available energy injected
into non-thermal electrons) sources within 1 kpc from Earth might be
detectable in X-rays with future instruments.

###### Acknowledgements.

V.B-R. & G.E.R. acknowledge financial support from the _State Agency for
Research_ of the Spanish Ministry of Science and Innovation under grant
PID2019-105510GB-C31 and through the _Unit of Excellence María de Maeztu
2020-2023_ award to the Institute of Cosmos Sciences (CEX2019-000918-M).
V.B-R. is also supported by the Catalan DEC grant 2017 SGR 643, and is
Correspondent Researcher of CONICET, Argentina, at the IAR.

## References

  * Aharonian et al. (2019) Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nature Astronomy, 3, 561
  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550
  * Benaglia et al. (2021) Benaglia, P., del Palacio, S., Hales, C., & Colazo, M. E. 2021, MNRAS, 503, 2514
  * Benaglia et al. (2010) Benaglia, P., Romero, G. E., Martí, J., Peri, C. S., & Araudo, A. T. 2010, A&A, 517, L10
  * Blumenthal & Gould (1970) Blumenthal, G. R. & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237
  * Brown (2015) Brown, W. R. 2015, ARA&A, 53, 15
  * Bruel et al. (2018) Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, arXiv e-prints, arXiv:1810.11394
  * Caprioli & Spitkovsky (2014) Caprioli, D. & Spitkovsky, A. 2014, ApJ, 783, 91
  * Cassano et al. (2018) Cassano, R., Fender, R., Ferrari, C., et al. 2018, arXiv e-prints, arXiv:1807.09080
  * Christie et al. (2016) Christie, I. M., Petropoulou, M., Mimica, P., & Giannios, D. 2016, MNRAS, 459, 2420
  * Comeron & Kaper (1998) Comeron, F. & Kaper, L. 1998, A&A, 338, 273
  * de la Cita et al. (2017) de la Cita, V. M., Bosch-Ramon, V., Paredes-Fortuny, X., Khangulyan, D., & Perucho, M. 2017, A&A, 598, A13
  * del Palacio et al. (2018) del Palacio, S., Bosch-Ramon, V., Müller, A. L., & Romero, G. E. 2018, A&A, 617, A13
  * del Valle & Pohl (2018) del Valle, M. V. & Pohl, M. 2018, ApJ, 864, 19
  * del Valle & Romero (2012) del Valle, M. V. & Romero, G. E. 2012, A&A, 543, A56
  * del Valle et al. (2015) del Valle, M. V., Romero, G. E., & Santos-Lima, R. 2015, MNRAS, 448, 207
  * Dremova et al. (2017) Dremova, G. N., Dremov, V. V., & Tutukov, A. V. 2017, Astronomy Reports, 61, 573
  * Falceta-Gonçalves & Abraham (2012) Falceta-Gonçalves, D. & Abraham, Z. 2012, MNRAS, 423, 1562
  * Guillochon & Loeb (2015) Guillochon, J. & Loeb, A. 2015, ApJ, 806, 124
  * Harmanec (1988) Harmanec, P. 1988, Bulletin of the Astronomical Institutes of Czechoslovakia, 39, 329
  * Hills (1988) Hills, J. G. 1988, Nature, 331, 687
  * Khangulyan et al. (2014) Khangulyan, D., Aharonian, F. A., & Kelner, S. R. 2014, ApJ, 783, 100
  * Kobulnicky et al. (2019) Kobulnicky, H. A., Chick, W. T., & Povich, M. S. 2019, AJ, 158, 73
  * Kobulnicky et al. (2016) Kobulnicky, H. A., Chick, W. T., Schurhammer, D. P., et al. 2016, ApJS, 227, 18
  * Kobulnicky et al. (2017) Kobulnicky, H. A., Schurhammer, D. P., Baldwin, D. J., et al. 2017, AJ, 154, 201
  * Koposov et al. (2020) Koposov, S. E., Boubert, D., Li, T. S., et al. 2020, MNRAS, 491, 2465
  * Kreuzer et al. (2020) Kreuzer, S., Irrgang, A., & Heber, U. 2020, A&A, 637, A53
  * Krtička (2014) Krtička, J. 2014, A&A, 564, A70
  * Loeb & Guillochon (2016) Loeb, A. & Guillochon, J. 2016, Annals of Mathematical Sciences and Applications, 1, 183
  * Maier (2019) Maier, G. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 733
  * Marchetti et al. (2019) Marchetti, T., Rossi, E. M., & Brown, A. G. A. 2019, MNRAS, 490, 157
  * Marcowith et al. (2016) Marcowith, A., Bret, A., Bykov, A., et al. 2016, Reports on Progress in Physics, 79, 046901
  * Martínez et al. (2021) Martínez, J. R., del Palacio, S., & Romero, G. E. 2021, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 62, 274
  * Matthews et al. (2020) Matthews, J. H., Bell, A. R., & Blundell, K. M. 2020, New A Rev., 89, 101543
  * McKinnon et al. (2019) McKinnon, M., Beasley, A., Murphy, E., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 81
  * Meyer et al. (2016) Meyer, D. M. A., van Marle, A. J., Kuiper, R., & Kley, W. 2016, MNRAS, 459, 1146
  * Morlino et al. (2021) Morlino, G., Blasi, P., Peretti, E., & Cristofari, P. 2021, MNRAS, 504, 6096
  * Myasnikov et al. (1998) Myasnikov, A. V., Zhekov, S. A., & Belov, N. A. 1998, MNRAS, 298, 1021
  * Parkin et al. (2014) Parkin, E. R., Pittard, J. M., Nazé, Y., & Blomme, R. 2014, A&A, 570, A10
  * Peri et al. (2012) Peri, C. S., Benaglia, P., Brookes, D. P., Stevens, I. R., & Isequilla, N. L. 2012, A&A, 538, A108
  * Prajapati et al. (2019) Prajapati, P., Tej, A., del Palacio, S., et al. 2019, ApJ, 884, L49
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Sánchez-Ayaso et al. (2018) Sánchez-Ayaso, E., del Valle, M. V., Martí, J., Romero, G. E., & Luque-Escamilla, P. L. 2018, ApJ, 861, 32
  * Tutukov & Fedorova (2009) Tutukov, A. V. & Fedorova, A. V. 2009, Astronomy Reports, 53, 839
  * van Buren & McCray (1988) van Buren, D. & McCray, R. 1988, ApJ, 329, L93
  * van den Eijnden et al. (2021) van den Eijnden, J., Heywood, I., Fender, R., et al. 2021, arXiv e-prints, arXiv:2111.10159
  * Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377
  * Wilkin (1996) Wilkin, F. P. 1996, ApJ, 459, L31
  * Zhang et al. (2013) Zhang, F., Lu, Y., & Yu, Q. 2013, ApJ, 768, 153

##  Appendix A Comparison of the hydrodynamical model

We compare the values of the pressure, density, and speed of the shocked fluid
obtained with the prescriptions used in this work (Sec. 2.3) with those
obtained by considering Rankine-Hugoniot jump conditions. We summarise the
results in Fig. 8, which shows the ratio between these quantities calculated
with each formalism for different angles $\theta$. Both prescriptions yield
very similar speeds for the shocked gas, although the prescription used in
this work yield greater values for the remaining thermodynamic quantities. The
discrepancy on the pressure increases with the angle $\theta$ until
$\theta\sim 60\degree$, while the discrepancy on the density keeps increasing
with $\theta$ and it reaches a difference of a factor two at $\theta\sim
140\degree$.

Figure 8: Ratio between the values given by the hydrodynamic prescription used
in this work and the values using Rankine-Hugoniot jump conditions for strong
shocks for the RS. The green, blue, and orange lines correspond to the
pressure, density, and speed of the shocked fluid, respectively. The black
dotted line at unity is marked as a reference.