Stellar wind bubbles of OB stars as Galactic cosmic-ray re-accelerators

Cosmic rays are highly energetic messengers propagating in magnetized plasma, which are, possibly but not exclusively, accelerated at astrophysical shocks. Amongst the variety of astrophysical objects presenting shocks, the huge circumstellar stellar wind bubbles forming around very massive stars, are potential non thermal emitters. We present the 1D magnetohydrodynamical simulation of the evolving magnetized surroundings of a single, OB type main sequence 60 Mo star, which is post processed to calculate the re-acceleration of preexisting non-thermal particles of the Galactic cosmic ray background. It is found that the forward shock of such circumstellar bubble can, during the early phase (1 Myr) of its expansion, act as a substantial reaccelerator of pre existing interstellar cosmic rays. This results in an increasing excess emission flux by a factor of 5, the hadronic component producing gamma-rays by pion0 decay being more important than those by synchrotron and inverse Compton radiation mechanisms. We propose that this effect is at work in the circumstellar environments of massive stars in general and we conjecture that other nebulae such as the stellar wind bow shocks of runaway massive stars also act as Galactic cosmic-ray re-accelerators. Particularly, this study supports the interpretation of the enhanced hadronic emission flux measured from the surroundings of kappa Ori as originating from the acceleration of pre-existing particles at the forward shock of its wind bubble.


INTRODUCTION
The detailed components of the non-thermal spectrum of the interstellar medium (ISM) of the Milky Way is far from being understood.A source of Galactic cosmic rays is the so-called remnants nebulae, produced by the fast expansion of a supernova blastwave, released by the death of a star or a stellar system, into its local ambient medium (Shklovskii 1954;Weiler & Sramek 1988;Blasi 2011;Vink 2012Vink , 2020)).Such cosmic rays, principally emitting in radio, X-rays, and -rays via inverse Compton, non-thermal synchrotron mechanisms and hadronic processes, are well documented thanks to the plethora of data collected by high-energy observatories such as Fermi and H.E.S.S., which permitted to confront theoretical models like the diffusive shock acceleration (Schatzman 1963;Axford et al. 1977;Bell 1978a,b;Blandford & Ostriker 1978;Krymsky et al. 1979) with observations (Abdo et al. 2010;H. E. S. S. Collaboration 2018).The forthcoming Cherenkov Telescope Array (CTA ground-based observatory will continue this ongoing effort (Acharyya 2023;Acero 2023).However, if shocks and turbulence in supernova remnants are identified accelerators of charged particles, at least in the sub-PeV band of the cosmic ray spectrum, additional unidentified sources able to bring electrons and protons of the ISM to very high kinetic energies must exist to fully explain the Galactic cosmic ray spectrum before and beyond the knee (Cristofari et al. 2021).Amongst all potential cosmic-ray accelerators, wind collision in binary systems (Reimer et al. 2007;Hamaguchi et al. 2018;Pittard et al. 2021), pulsars nebulae (Bednarek & Bartosik 2004) but also superbubbles around stellar clusters with massive stars (Bykov 2001;Butt & Bykov 2008;Morlino et al. 2021;Morlino 2021) have been theoretically demonstrated to present the ability to substantially accelerate particles within the diffusive shock acceleration, which turned out to be consistent with high-energy observations (Joubaud et al. 2020).Particularly, there is a growing suspicion for the circumstellar shocks found in nebulae generated by the interaction between the strong winds of high-mass stars and the ISM, to be a major source of Galactic cosmic rays (Cardillo 2019;Meyer et al. 2020).
The basic picture of the circumstellar medium of massive stars is the wind-ISM interaction model of Weaver et al. (1977), in which the supersonic wind of a main-sequence massive star collides with its constant local ambient medium.This results in a structured nebula, delimited by an inner termination shock (the so-called reverse shock) and an outer so-called forward shock, and centered onto the region of free-steaming stellar wind.The gas between the reverse and forward shocks is made of two layers of shocked diluted hot stellar wind and shocked dense cold ISM material, separated by a contact discontinuity.Those simple objects motivated numerous studies, from the beginning of numerical fluid dynamics to recent computing-intensive calculations, revealing their complex functioning, moderated by radiative cooling (van Marle & Keppens 2010), thermal conduction (Zhekov & Myasnikov 1998), stellar evolution (Avedisova 1972;Garcia-Segura & Mac Low 1995a,b;Garcia-Segura et al. 1996) and their sensitivity to the large-scale magnetization of the ISM (van Marle et al. 2015).If the star moves fast through the ISM, the wind bubble is turned into a bow shock (Gull & Sofia 1979;Kobulnicky et al. 2010;Cox et al. 2012;Gvaramadze et al. 2014;Kobulnicky et al. 2017), characterized by an arc-like shape, an internal organization similar to circumstellar wind bubbles, and by a radiative forward shock that further stimulated the analytic and numerical exploration of the surroundings of massive runaway stars (Wilkin 1996;Comerón & Kaper 1998;van Marle et al. 2011;Meyer et al. 2016;Acreman et al. 2016).Stellar wind bubbles and bow shocks have been largely studied as the host site where the most massive stars die as a core-collapse supernova event (Dwarkadas 2005(Dwarkadas , 2007)), giving birth to complex supernova remnants (van Veelen et al. 2009;Meyer et al. 2015Meyer et al. , 2021;;Velázquez et al. 2023;Villagran et al. 2024) or -ray bursts (van Marle et al. 2005(van Marle et al. , 2006(van Marle et al. , 2007)).Nevertheless, their role as cosmic-rays accelerators still deserve investigations.
The Earth's closest circumstellar structure is the termination shock of the solar wind interacting with the ISM.It has been theoretically suggested (Jokipii 1968;Jokipii et al. 2004) and interpreted in the context of measures of the spacecraft Voyager to produce a local component of cosmic rays (Ferreira et al. 2008;Jokipii & Kóta 2014;Pogorelov et al. 2017).This idea has been extended to the termination (reverse) shock of circumstellar wind bubbles around main-sequence, OB-type massive stars (Casse & Paul 1980), although this results have been questioned (Voelk & Forman 1982) but also confirmed (Webb et al. 1985).Equivalently, the raising number of detected stellar wind bow shocks of supersonically-moving high-mass stars (Peri et al. 2012(Peri et al. , 2015;;Kobulnicky et al. 2016Kobulnicky et al. , 2018) ) spurred the interest of the heliospheric community to massive stellar objects such as the OB star  Cephei, predicting local cosmic-ray anisotropies to be generated inside of such stellar wind cavities (Scherer et al. 2015;Herbst et al. 2022).
The observational and experimental quest for cosmic rays originating from the surroundings of massive stars is a vivid field of research.On the one hand, the astrospheres of runaway massive stars are predicted to be the site of particle acceleration (del Valle & Romero 2012, 2014;del Valle et al. 2015;del Valle & Pohl 2018).Such emissions have been detected from the bow shock EB27 around BD+43 • 3654 (Benaglia et al. 2010;del Palacio et al. 2018;Benaglia et al. 2021) and possible association with  -rays data have been reported in Sánchez-Ayaso et al. (2018).Recently, a second stellar wind bow shock from a massive runaway star revealed in its turn radio non-thermal signature (Moutzouri et al. 2022).On the other hand, more systematic attempts to measure, e.g.high-energy X-rays (Schulz et al. 2014;Toalá et al. 2016bToalá et al. , 2017;;De Becker et al. 2017;Binder et al. 2019) and very high-energy emission (Schulz et al. 2014) from the surroundings of stellar wind bow shocks turned to be unfruitful.At the same time, the surroundings of evolved Wolf-Rayet stars seem to be non-thermal emitters (Prajapati et al. 2019), as predicted by Zirakashvili & Ptuskin (2018) and phenomenological studies pointing out the importance of cosmic-rays re-acceleration in the circumstellar shocks (Cardillo 2019).This puzzling situation leaves wide open the question of the existence of cosmic rays from circumstellar shock.This paper is a numerical investigation of the non-thermal phenomenon in the context of the circumstellar medium of a very massive (⩾ 60 M ⊙ ) star.We explore, by means of 1D magnetohydrodynamical (MHD) simulations completed by particle acceleration calculations, the possibility of expanding stellar wind bubbles around a very massive star to generate non-thermal emission by the re-acceleration of pre-existing Galactic background cosmic rays.We established time-dependent predictions for the non-thermal synchrotron radiation, inverse Compton emission, and hadronic -rays feedback by   decay of the surroundings of a very massive star, which we further consider in the context of observational data.
Our study is organized as follows.We first introduce the reader to the various numerical methods that we have used to simulate the 1D MHD stellar wind bubble of a 60 M ⊙ star and to calculate its non-thermal emission in Section 2. In Section 3 we present results for the MHD structure of the surroundings of very massive stars, the particle acceleration at the shocks generated by wind-ISM interaction, and the resulting non-thermal emission.The caveats of our results are further detailed in Section 4, discussed in the context of observations and we conclude in Section 5.

METHODS
In this section, we present the numerical methods used to model the magnetized circumstellar medium of a very massive star and to calculate the acceleration of particles at work therein.

Governing equations
The problem of the plasma in MHD circumstellar bubbles of massive stars is described as follows by the non-ideal MHD equations, and,

𝜕 𝑩 𝜕𝑡
with  the gas density,  the plasma velocity vector,  the thermal pressure of the gas, Î the identity matrix, the linear momentum vector,  the magnetic field vector and the total energy of the system.This system of equations is closed by the definition of the adiabatic sound speed in the medium, with  = 5/3 the adiabatic index for ideal gas.
Losses by optically-thin radiative processes are represented by the term, with, the gas temperature where  B is the Boltzmann constant,  H the mass of the proton,  is mean molecular weight, set to 0.61 for gas ionized by the radiation of hot massive stars.The function Λ() stands for the cooling rate, while Γ() represents the heating rate of the gas, respectively.We refer the reader interested to the detailed description of the Γ() and Λ() to Meyer et al. (2014).
The MHD simulations of circumstellar bubbles of massive stars are carried out with the pluto code (Mignone et al. 2007(Mignone et al. , 2012;;Vaidya et al. 2018) using the Harten-Lax-van Leer Riemann solver (Harten et al. 1983) with a 2 nd -order numerical scheme together with linear reconstruction of the primitive variables.We make use of the eight-wave MHD formulation by Powell (1997) which ensures that the divergence-free condition of the magnetic field vector, is satisfied everywhere in the computational domain, throughout the whole calculation.Last, the timestep of the numerical simulations is controlled by the Courant-Friedrich-Levy criterion, which we set initially to  cfl = 0.1.

Initial conditions
We perform our MHD wind bubble models using a 1D sphericallysymmetric coordinate system and a [0,  rmax ] computational domain with  max = 150 pc, that is mapped with a uniform grid of 2000 zones.To avoid spurious MHD effects at the stellar wind boundary, we first model the early 0.01 Myr of the wind-ISM interaction in a pure hydrodynamical manner and using a smaller domain of  max = 15 pc uniformly discretised with 1000 zones, before mapping the solution onto the larger domain.As in Dwarkadas (2005Dwarkadas ( , 2007)), we perform the simulations in the reference frame of the star.The stellar wind is launched within a spherical zone of 20 grid zones centered onto the origin of the computational domain.
The wind density is as follows, with  the distance to the star,  () the mass-loss rate and  w () the time-dependent stellar wind terminal velocity.The various surface properties of the star (  (),  w ()) are taken from the tabulated stellar evolutionary track for the 60 M ⊙ star of Groh et al. (2014), which has already been used to simulate wind bubbles a two-dimensional fashion (Meyer et al. 2020).The ambient medium has properties corresponding to the warm phase of the ISM, with constant number density  ISM = 0.79 cm −3 and temperature 8000 K, respectively.Particle acceleration being of intrinsic MHD nature, a magnetic field is added into the outflowing stellar wind.The magnetic field is incorporated into the model at the time 0.01 Myr, at the moment of the mapping of the hydrodynamical wind-ISM interaction to the MHD grid.The stellar boundary conditions are made of a radial component, and of a toroidal component, where  ★ is the stellar surface magnetic field,  ★ the stellar radius and   =  rot is the toroidal stellar rotation velocity at the equator (Meyer et al. 2021;Meyer 2021a).Since we perform 1D MHD simulations, any latitude-dependence of the stellar wind has disappeared and we calculate the models in the equatorial plane of the spherical coordinate system (, = /2,).In this study, we adopt the stellar surface properties of  ★ = 50 G, and  rot = 50 km −1 , which are characteristic of OB stars (Hubrig et al. 2013) Inside of the computational domain, we impose a reconstructed magnetic field, taking the form of a piece-wise function.In the remaining part of the paper, we will note the termination shock  R , the contact discontinuity  CD and the forward shock  F .The freely-expanding wind region ( ⩽  R ) and the shocked stellar wind ( R ⩽  ⩽  CD ), are filled with, while the region of shocked ISM ( CD ⩽  ⩽  F ) is filled as,  r =   =   =  ISM / √ 3, where  is the compression ratio of the shock, taken to  = 4. Hence, the value of the compressed ISM magnetic field in the outer layer of the wind bubble is, with  r =   =   =  ISM / √ 3 and  ISM = 7 G. Similarly, we consider the magnetization of the ISM, also with  r =   =   =  ISM / √ 3 in each grid zone of the computational domain  F < .After the mapping, the 1D MHD simulations are restarted at a time  = 0.01 and further calculated over the part of the main-sequence phase for which the compression ratio of the forward shock propagating into the ISM is larger than 4.This corresponds to the early 1 Myr of the star's evolution (Groh et al. 2014) and we temporally discretize our simulation so that about 100 MHD simulation outputs are produced.

Particle acceleration calculations
The 1D MHD stellar wind bubble model simulated with the pluto code is post-processed with the particle acceleration code ratpac.The Radiation Acceleration Transport Parallel Code (ratpac) is a numerical framework solving the acceleration and subsequent transport of electrons and protons in a magnetized plasma presenting shocks.It is a modular, parallelized code written in python and making use of the fipy library of finite volume solvers for partial differential equations.Originally developed in the context of supernova remnants, ratpac post-processes arbitrary pre-calculated or analytic plasma flows including a pattern of MHD shocks (Telezhinsky et al. 2012a(Telezhinsky et al. ,b, 2013;;Wilhelm et al. 2020), although it can self-consistently be used on the fly, coupled together with an MHD solver such as the pluto code (Brose et al. 2016;Sushch et al. 2018;Brose et al. 2021).
The MHD profiles (, ,  ) are loaded into ratpac and, at each timestep of the calculation, the MHD flow profile is interpolated between the temporal neighboring snapshots.The coordinate transformation, with  = / F , permits us to greatly reduce numerical costs and to reach spatial resolution near the shock region much finer than that of the pre-calculated MHD profiles.Once the interpolated MHD profile is transformed to the new coordinate system, a shock finder algorithm looks for the position of the forward shock of the stellar wind bubble  F .The position of the shock is found by looping through the radial direction, from the ISM to the star,  F being caught via the compression ratio  ⩾ 2 of the density.To compensate for the poor resolution in the vicinity of the shock, the velocity profile is fitted with a logistic function using 5 grid zones upstream and downstream of the forward shock.The velocity flow is finally resharpened to ensure a clear sharp jump around  F .The time-dependent numerical method is based on kinetic calculations within the test-particle approximation, which neglects any feedback of the cosmic rays to the gas dynamics.The main output of ratpac is the non-thermal particle number density  (, , ) with  the radial coordinate,  the momentum (or equivalently the energy ) of the particle and  the time, respectively.The evolution of  is determined in spherical coordinates, with a grid that is co-moving with the considered shock, where the following diffusion-advection equation is solved, where  is the spatial diffusion coefficient,  is the vector velocity,  / represents the energy losses.

Diffusion model
A Böhm-like diffusion model is adopted in this study, controlled by a free parameter  .We adopt a coefficient that is a piece-wise function of the distance to the central star (Telezhinsky et al. 2012a(Telezhinsky et al. ,b, 2013;;Wilhelm et al. 2020), which includes a free parameter  accounting for sub-grid microphysical diffusion processes at work in the interior of the circumstellar bubble, such as magnetic turbulence cascade, not included in our simulations.It reads, where is the diffusion coefficient of the Galactic cosmic ray background acting as a bath of non-thermal particles in which the wind bubble evolves, with  = 1/3 and  0 = 10 29 cm 2 s −1 , and where, represents the exponential transition between the Böhm-based diffusion coefficient in the bubble and that in the ISM.Our work explores the effects of the diffusion coefficient by varying it in the range 10 ⩽  ⩽ 400.Our diffusion coefficient strategy is further illustrated in Fig. 1.

Energy threshold and particle re-acceleration
We include the re-acceleration of Galactic cosmic rays into the simulations by switching-off the injection of particles (Wilhelm et al. 2020) and by imposing a reservoir of non-thermal particle distribution taken the Galactic cosmic-ray background spectrum.The non-thermal particles are initially distributed using a uniform spatial distribution from the outer boundary of the computational domain up to the termination shock of the stellar wind bubble.The initial time of the particle acceleration simulation is performed when the stellar wind bubble is 0.5 Myr old, time instance from which the termination shock, contact discontinuity, and forward shock of the circumstellar bubble, is clearly distinguishable and capturable by the shock/discontinuity-finder algorithm.Neglecting the early phase of the bubble's expansion underestimate ISM particle re-acceleration, potentially important as the compression ratio of the forward shock is large.The goal of this work is not to estimate the total, precise non-thermal budget of the wind bubble, but rather to qualitatively corroborate previous observations, see Section 4. The cosmic-ray spectrum that we adopt is that of Brose et al. (2016) and Wilhelm et al. (2020), modified for the purpose of the present study: since the Galactic cosmic rays do not reach the termination shock of wind bubbles around massive stars (Voelk & Forman 1982), the pre-existing non-thermal particles are imposed from the ISM to the contact discontinuity of the bubble  CD .It reads, where  Gal () is taken from Wilhelm et al. (2020) for the electrons, and from Moskalenko et al. (2002) and Jaffe et al. (2011) for the protons, respectively.To be re-accelerated at the shock by diffusive shock acceleration, the particles need to have energy sufficiently high to cross the shock the first time and begin the Fermi cycles.We call  th the threshold energy necessary to trigger the re-acceleration process of a particle.This lower-energy cut-off of the distribution of available cosmic ray is the governing parameter of diffusive re-acceleration (Thornbury & Drury 2014;Drury & Strong 2015).We introduce a corresponding cut-off in the Galactic cosmic-ray background in order not to artificially take into account the numerous low-energy particles which would bias the solution by being considered as accelerable by the code while they can not cross the forward shock.
The quantity  th is evaluated by forcing the pre-existing nonthermal particles at the forward shock to have properties such that their acceleration characteristic lengthscale (the precursor length  acc ) is larger than the shock's radiative characteristic lengthscale (the cooling length  cool ).They read, and respectively, with  F the speed of the forward shock expanding through the ISM and  the gas velocity in the post-shock region at the forward shock.In the above relation,  cool is the cooling timescale, as defined in Meyer et al. (2017), with  u =  ISM the number density upstream the shock,  d =  u the number density downstream of the magnetized forward shock of the bubble measured from the simulation, and with, the compression ratio at the forward shock, and, where  and  A are the Mach and Alfvenic Mach numbers upstream the shock, respectively (Shu 1992).Similarly, the tempera-ture downstream of the shock, is the gas temperature downstream the shock, while  ISM = 8000 K is the upstream shock temperature, and Λ() is the rate for opticallythin cooling for ionized gas that is interpolated from the cooling curve presented in Meyer et al. (2014).The Böhm diffusion coefficient reads, with  the speed of light,  is the elementary charge, and 3 the compressed magnetic field downstream of the forward shock.
The value of the cut-off momentum  th is obtained by solving  cool =  acc for , which straightforwardly provides the cut-off energy  th , for both electrons and protons.The evolution of  th for several values of  is plotted in Fig. 2. It shows that the cut-off energy is of the same order of magnitude during an early 1 Myr of the star's life, and decreases with time as the wind bubble expands and the compression ratio of its forward shock gets smaller.Our results are consequently a lower limit estimate of the non-thermal feedback for the outer circumstellar shocks around early-type, mainsequence massive stars.Without including this initial non-thermal population in the ISM, the process of accelerating particles at shocks during the wind bubble initial expansion phase is inefficient and its associated emission is consequently negligible.This capability to accelerate particles decreases as the bubble further grows and its compression ratio decreases, as accelerated protons escape and as the electrons lose energy by radiation mechanisms.

RESULTS
We hereby present the 1D MHD simulations of the wind-ISM interaction of a 60 M ⊙ massive star, together with the cosmic-ray distribution from the particle acceleration calculation and the bubble's emission spectra.

MHD wind bubble of magnetized very massive star
In Fig. 3 we show the number density field (in cm −3 ) and the magnetic field structure (in  G) in the stellar wind bubble produced by wind-ISM interaction around a 60 M ⊙ massive star.The initial conditions at the moment of the mapping of the initial hydrodynamical wind-ISM calculation as a MHD simulation on a larger grid.The density field already exhibits the typical structure of a wind bubble as described in Weaver et al. (1977), with a reverse shock, a contact discontinuity and the forward shock solid blue line).The freely-expanding stellar wind characterised by a number density decreasing ∝ 1/ 2 from the origin of the domain to the termination shock.A region of low-density shocked wind establishes itself between the termination shock and the contact discontinuity, while a dense region of shocked ISM material forms between the contact discontinuity and the forward shock, respectively (Fig. 3a).The reconstructed magnetic field consists in a Parker wind in the section between the origin of the domain and the contact discontinuity.The compressed magnetic field is imposed in the dense region of ISM gas, and the ISM is filled with the constant field  ISM .As the central star evolves through the main-sequence phase of its life and blows its strong wind, corresponding to a mass-loss rate  ≈ 10 −6 M ⊙ yr −1 and wind velocity  w ≈ 2000 km s −1 .Under the effect of the momentum injected into the circumstellar medium, the ram pressure of the wind pushes the structured bubbles to larger radial distances in the ISM (Fig. 3b).The toroidal component of the magnetic field is compressed at the termination shock, and, as it governs the total magnetic field since the radial field decreases ∝ 1/ 2 according to the stellar wind density.It then increases as a function of the radius from the star  and is it proportional to the density jump in the region of shocked ISM.This magnetic field structure is conserved throughout the whole simulation (Fig. 3c).Fig. 4 plots the properties of the expanding forward shock of the stellar wind bubble over the early 1 Myr of the evolution of the mainsequence 60 M ⊙ .At time 0.4 Myr, the position of the forward shock is ≈ 26 pc from the star, and this distance quasi-linearly increases as a function of time up to times > 1 Myr, as the volume of the wind bubble augments (Fig. 4a).The shock speed, i.e. the velocity at which the forward shock goes through the local ambient medium monotonically decreases from  F ≈ 38 km s −1 at time 0.4 Myr to  F ≈ 27 km s −1 at time 1.0 Myr, respectively (Fig. 4b).Similarly, the value compression ratio of the forward shock decreases from  ≈ 5.5 to  ≈ 3.8 over the same time interval (Fig. 4c).At later times, the expansion of the wind bubble slows down, and the forward shock becomes weaker as the bubble reaches its final size (Meyer et al. 2020).We start the post-processing of the MHD wind bubble with the particle acceleration code at 0.4 Myr, because at earlier times, the number of grid zones resolving the layer of shocked ISM between the contact discontinuity and the forward shock, is not sufficient to permit the shock finder algorithm to properly track the forward shock.The cosmic ray calculation is continued up to the moment the compression ratio at the forward shock is  ⩽ 4, at times ≈ 1 Myr.

Non-thermal particle distribution
In Fig. 5 we display the distribution of non-thermal particles for electrons (left) and protons (right) as a function of the normalized momentum, for the different diffusion coefficient models, differing by the parameter  , see Eq. 19.The number density of non-thermal particles (in cm −3 ) is shown for several time instances, 0.5 Myr (dotted lines), 0.7 Myr (dashed lines) and 1.0 Myr (solid line).The blue line is the initial Galactic cosmic-ray background that we assume, and which is cut-off at lower momentum for a given threshold energy.The energy threshold (orange vertical lines) is the energy under which no acceleration is possible, and it depends on the adopted diffusion model, i.e. on the parameter  .
As the system evolves, a modification of the non-thermal particles distribution with respect to the initial Galactic cosmic-ray spectrum happens.An excess of non-thermal electrons above the initial population is produced for a diffusion model  = 10, together with a shift of its maximum momentum to lower values, because of the efficient synchrotron cooling at work therein (Fig. 5a).A similar excess of particles is visible in the proton distribution at momenta larger than the initial cut-off of the Galactic background cosmic-ray spectrum, which does not have its maximum momentum shifted, since no efficient cooling is at work regarding this particle specie (Fig. 5d).The same trend can be seen in the evolution of the nonthermal electron and proton distribution calculated with different diffusion laws, i.e. a larger Böhm-like diffusion factor and a smaller momentum cut-off of the initial Galactic spectrum, providing more injected pre-existing particle to re-accelerate, because of the smaller threshold energy (Fig. 5b,e,c,f).
Fig. 6 compares the number density distribution of non-thermal particles at the end of the simulations (1 Myr) for different diffusion models (10 ⩽  ⩽ 400).The accelerated particles of the model with  = 10 reach a higher maximum energy than the models with  = 200 and  = 400, respectively, because their energy threshold is higher, and their acceleration timescale is shorter as a result of the smaller diffusion coefficient.Hence, a more important accumulation of non-thermal particles with  = 10 than with  = 200 or with  = 400 happens.The maximum energy  max ∝  2  reached by the non-thermal electrons is similar in all models, regardless of the parameter  because of the losses by synchrotron emission that are independent of the diffusion coefficient, and cool them to the same energy (Fig. 6a).On the other hand, the distribution of non-thermal protons do not exhibit a deficit of particles with respect to the initial value of the Galactic cosmic-ray background, as they do not experience cooling mechanisms, see Fig. 6b.
To better understand the acceleration process, one can have a closer look at the evolution of the distribution of non-thermal protons in the vicinity of the forward shock the stellar wind bubble, for given selected energies, as a function of time.Fig. 7 plots the radial particle distribution in the calculation for the diffusion model with  = 200, for energies  th and  max throughout the simulation, for several time instances.The spatial coordinate is shown as normalised to the radius of the forward shock.The black line corresponds to the values of the pre-existing Galactic non-thermal bath.In panel Fig. 7b, there is a peak in the particle distribution at the forward shock because of the confinement of particles due to the re-acceleration of Galactic particles, followed by an exponential decrease behind the forward shock.The slope of the radial distribution behind the shock increases since the magnetic field in the post-shock region decreases at the forward shock of the expanding wind bubble.The Böhm-like diffusion coefficient goes as / (see solid and dashed lines of Fig. 7b).In panel Fig. 7a, the acceleration of particles at the threshold energy  th induces a decrease of protons at the forward shock.There are no more particles with energy lower than  th to fill the void: the Böhm-like mechanism is energy-dependent, therefore the lower the energy and the longer the timsecale to diffuse such protons, hence the re-filling of this void is only possible by replenishing it with lower-energy cosmic rays brought by spatial advection towards the forward shock.

Emission spectra
Fig. 8 plots the evolution of synchrotron emission of the forward shock of the circumstellar stellar wind bubble forming around the massive star, at selected time instances during the early 1 Myr of its main-sequence evolution.The colors distinguish between several adopted Boehm-based diffusion models, and the several lines of panels distinguish between the emission mechanisms, i.e., synchrotron (top panels), inverse Compton scattering with cosmic microwave background photons (second line of panels),   decays (third line of panels), and all processes together (bottom panels).The emission flux  is calculated as in Das et al. (2022) at the forward shock of the wind bubble, where the particles are re-accelerated.It is converted into luminosities by integrating the flux over a volume.The luminosity represents a shell with an inner radius at the contact discontinuity and an outer radius extending 0.5 pc beyond the forward shock.
The fluxes of the various components of the emission spectrum increase as a function of time, which translates the continuous reacceleration of pre-existing non-thermal particles at the forward shock of the bubble, see differences between the left-hand and right-hand columns of panels.The fluxes are more important for the synchrotron and the hadronic emission mechanism.The effects of cooling diminish the maximum energy of the synchrotron flux distribution (Fig. 8,1a-1c), while the maximum energy reached by the photon generated via inverse Compton or by hadronic processes is constant over time and only differs according to the used diffusion model (Fig. 8,2a-2c"3a-3c).The differences caused by the adopted diffusion models show that, in large part, the Boehm-like coefficient corresponds to a higher emission flux.The emission flux by synchrotron radiation is in the 10 −6 -10 2 eV with a peak at about 10 −4 eV and which maximum increases, for the several efficient diffusion models, by a factor of 5 during the 1 Myr of the bubble's life.Similarly, the emission flux by the inverse Compton mechanism lies within the 10 −4 -10 2 GeV energy band.It peaks at about 10 −3 GeV but displays a plateau up to about 1 GeV, and its maximum flux increases by a factor of 6.6 over the early main-sequence of the star.The emission flux by   decays is in the 10 −5 -10 1 TeV energy band, peaks at about 10 −3 TeV and its maximum increases by a factor of 5 during the simulation (Fig. 8).The forward shock of the expanding wind bubble is therefore able to increase the non-thermal emission by a factor of 5 compare to the local cosmic ray background.The hadronic flux is more important than that by inverse Compton, which is itself larger than the synchrotron emission flux.The important hadronic flux comes from the cut-off energy of the Galactic cosmic-ray spectrum that is smaller for protons than for the electrons responsible for the emission by synchrotron and inverse Compton mechanisms.

DISCUSSION
This section presents the limitations of our model.Finally, we compare our findings with observational data and discuss our results in the context of massive stars and other observations.

Caveats of the model
First, although we time-dependently interpolate the stellar surface properties of its central massive star, the simulation that we perform for its circumstellar medium are one-dimensional, see also Weaver et al. (1977); Garcia-Segura et al. (1996); Zhekov & Myasnikov (1998); Dwarkadas (2005Dwarkadas ( , 2007)).This study is the first work presenting a 1D MHD model of stellar wind bubble including all three components of the velocity and magnetic fields, which is a step forward compare to the simulations of Weaver et al. (1977); Dwarkadas (2005Dwarkadas ( , 2007)).Nevertheless, our adopted 1D computing strategy still suffers from an evident loss of realism compare to multi-dimensional simulations (Freyer et al. 2003(Freyer et al. , 2006)).Especially, the development of instabilities at the contact discontinu- ity (Meyer et al. 2020), but also the turbulence of the gas in the region of shocked wind (Dwarkadas 2007) which change the properties in the layer of shocked ISM, are not included in our 1D MHD models.
A number of physical processes at work are also neglected.The most prominent of all is heat conduction, affecting both the thermal transfer and the properties at the termination shock and at the forward shock of the wind bubble.This could potentially change the compression ratio of the shocks, and, consequently, the efficiency of the acceleration of electrons and protons at this emplacement.The ambient medium hosting the central massive star is setup up in its most simple description: a uniform region with properties corresponding to the warm phase of the ISM ( ISM = 0.79 cm −3 and  ≈ 8000 K).A more realistic ambient medium, e.g.accounting for the native turbulence of the gas (Rogers & Pittard 2013;Mackey et al. 2015), its intrinsic clumpiness (Wareing et al. 2017;Pittard 2019), or the large-scale magnetic field of the ISM, see study of van Marle et al. (2015).
A few more points can be mentioned.In particular, the adiabatic compression of cosmic rays in the post-shock region of radiative shocks, resulting from density waves travelling with highdensity regions of the ISM, is a mechanism which could matter in the context of the circumstellar wind bubble (  1982).The pre-existing population of Galactic comic rays would then be increased by a factor proportional to the compression ratio of the forward shock (Uchiyama et al. 2010), increasing the local population of non-thermal particles available for re-acceleration.This is particularly important during the early expansion phase of the bubble, and could magnify the hadronic emission of the stellar surroundings of massive stars, as it is the case in some supernova remnants (Tang & Chevalier 2014;Cardillo et al. 2016;Tang 2019;Sushch & Brose 2023).This strengthens the results presented in this paper.
Secondly, the target photon field for the inverse Compton emission in our calculation does not include thermal infrared photons produced by dust trapped in the circumstellar medium and heated by the starlight.This has been show to matter in the context of stellar wind bow shocks around OB stars (del Valle & Pohl 2018).We produce our emission spectra at times about 0.5 Myr after the onset of the main-sequence phase, when the bubble is extended and its dusty ISM gas has already reached distances ⩾ 20 pc (Fig. 3b).Since the starlight infrared flux by unit solid angle diminishes by the 1/ 2 , with  the distance to the star, it is reasonably acceptable to neglect it in the present work than in del Valle & Pohl (2018).A more detailed study of that problem would be highly desirable in the future.
Third, our model does not account for non-thermal Bremsstrahlung, an emission mechanism of relativistic electrons  emitting when their path is deviated by the local ISM electrical field of charged particles.This produces from hard X-rays to -rays emission.Since it affects supernova shocks propagating the ISM at densities similar to that of the warm phase of the ISM, this mechanism could affect the non-thermal X-rays to -rays emission spectrum of stellar wind bubbles of OB star.A back-of-the-envelope estimate of the non-thermal relativistic Bremsstrahlung and inverse Compton emission timescales can further educate us on the respective importance of these emission processes.We adopt the prescriptions of Berezinskii et al. (1990) for the Bremsstrahlung timescale and of Longair (2011) for the inverse Compton electron lifetime.Results are plotted in Fig. 9, indicating that the Bremsstrahlung losses are quicker than the inverse Compton at energies < 10 4 GeV, while the inverse Compton electron lifetimes are shorter at higher energies > 10 4 GeV.Consequently, one should, in the future, consider in more detail the relative importance of this emission mechanism.Its luminosity might nevertheless be negligible compared to the inverse Compton integrated flux, as suggested by the study on bow shocks from runaway massive stars in dense environments by del Valle & Romero (2014).Our treatment of the particle acceleration mechanism is, in its turn, subject to limitations.We make use of the robust, but simple code ratpac developed for the non-thermal spectra of blastwave from supernova remnants expanding in a pre-shaped medium (Telezhinsky et al. 2012a;Wilhelm et al. 2020).However, our approach neglects some features implemented into ratpac, such as the amplification of the magnetic field ahead of the forward shock (Brose et al. 2016(Brose et al. , 2019(Brose et al. , 2021)).Last, one should underline that spatial resolution plays a preponderant tole in particle acceleration simulations.Hence, both the 1D MHD stellar wind bubble (Pittard et al. 2021) and the shock reconstruction during the particle acceleration calculation would benefit a high number of cells discretizing the respective regions of interest, in the vicinity of the forward shock.Potential solutions could be either a new coordinate transformation as that adopted in this study (see Eq. 17), or the use of so-called Particle-In-Cell simulations to simulate the propagation of the forward shock (Bohdan et al. 2021).

Generic comparison to other works
The first investigation of the non-thermal phenomenon at the reverse shock of circumstellar wind bubbles is that of Casse & Paul (1980), which concluded on its possible but overall inefficiency, see also (Voelk & Forman 1982;Webb et al. 1985).Because the stellar magnetic field decreases abruptly as a function of distance to the star because of a Parker spiral, its value at the termination shock is rather small, and the efficiency of the corresponding particle acceleration smaller.However, these series of studies assumed a Weaver-like bubble, which stellar wind parameters are today known to be far too high in terms of mass-loss rate (Brott et al. 2011).The acceleration of non-thermal particles in supernova remnants has been shown to add-up when multiple core-collapse explosions happen in OB associations, suggesting that, in the stellar clusters where the wind and supernova blastwave of several generations of high-mass stars collide, these effects should be very efficient Bykov & Toptygin (2001); Bykov (2001); Butt & Bykov (2008).Similarly, the study of Morlino et al. (2021) demonstrates the efficient production of cosmic rays at the termination shock of superbubbles around clusters of massive stars, using semi-analytic estimates.These studies concern the cosmic-ray acceleration at the termination shock of circumstellar nebulae surrounding single massive stars and/or superbubbles forming around a stellar cluster of high-mass stars, which is different from our approach since we show how native cosmic rays in the ISM are re-accelerated by the forward shocks of stellar wind bubbles around massive stars, which acts as a magnifier of the local non-thermal particle distribution.

Comparison to the case of 𝜅 Ori
This work would not be complete without a comparison between our results and the observation of non-thermal emission from the surroundings of the massive main-sequence B-type star  Ori (Cardillo et al. 2019).The historical star  Ori is located in the apparent constellation of Orion, also hosting the red supergiant star Betelgeuse.It is also a massive star, classified as a supergiant star of spectral type B0.5 Ia that is believed to be a variable, with moderately fluctuating magnitude (Marchili et al. 2018).It is embedded into a large stellar wind bubble, produced by its own wind-ISM interaction happening during its past main-sequence phase and which generated a large circumstellar structure according to the theory of Weaver et al. (1977).The fate of  Ori is that of a core-collapse supernova, which will take place inside of its wind-blown bubble and eventually produce a supernova remnant following the scenario that has been investigated in the studies of Dwarkadas (2005Dwarkadas ( , 2007) ) and Das et al. (2022).Interestingly, the region of  Ori exhibits an hadronic emission excess which has been constrained to originate from the re-acceleration of cosmic rays pre-existing as an ISM non-thermal bath in which the circumstellar medium grows, see data acquired by the Astro rivelatore Gamma a Immagini Leggero (AGILE) high-energy facility operating in the -ray energy band (Marchili et al. 2018;Cardillo et al. 2019).
In Fig. 10 we compare the hadronic emission spectra in our model with the AGILE data of the bubble around  Ori.The best fit obtained between the observational data and the present numerical results is that calculated with a modified Böhm-like diffusion model  = 400, assuming a distance to the source of 260 pc, which is of the order of that measured to  Ori Cardillo et al. (2019).This good agreement between our simulations and the observations from the vicinity of  Ori should not erase the fact the scenario that is investigated in our study is slightly different from that in the data reporting cosmic-rays re-acceleration.Indeed,  Ori is into an ambient medium that is denser than the assumed warm phase of the Galactic plane, as its main-sequence stellar wind interacted with a dense cloud (see paragraph below).Having such consistency between numerical model and observational data within the context of acceleration of pre-existing cosmic rays, fully supports the fact that the standard Böhm diffusion is not appropriate to depict the nonthermal evolution of the bubble of  Ori, as demonstrated by Cardillo (2019); Cardillo et al. (2019).
There is a degeneracy in terms of stellar surface properties (stellar mass-loss rate, wind velocity) and ISM density, respectively, in the shaping of wind-blown bubbles.A similar scaling relation exists in the shaping of stellar wind bow shocks around massive runaway stars, see discussion in section 3.1.3of Meyer et al. (2017).These elements authorise direct comparison between our model and the AGILE data.The theory of Weaver et al. (1977) gives the timedependent evolution of the forward shock of a wind bubble around massive stars, with  0 the ISM background number density,  36 = /10 36 erg s −1 the stellar mechanical luminosity and  6 = /10 6 yr the age of the star since the onset of the stellar wind bubble.With, where  w and  w are the wind mass-loss rate and terminal wind velocity, respectively, one obtains, which we can express for both the bubbles of  Ori and of our 60 M ⊙ star, and, respectively.By equalling them, Using the values of  ,60 M ⊙ = 2 × 10 −6 M ⊙ yr −1 and  ,60 M ⊙ = 2020 km s −1 that we use for the main-sequence phase of the massive star (Groh et al. 2014;Meyer et al. 2020) and similarly taking  w, = 2 × 10 −7 M ⊙ yr −1 and  w, = 2000 km s −1 as in Lennon et al. (1991); Cardillo et al. (2019), one finds, Hence, our adopted prescription for the wind power and the modelled time window of the star's main-sequence phase are consistent with the current knowledge that we have regarding  Ori, in terms of dimensions of its circumstellar wind bubble.Since many uncertainties persist regarding the values of the mass-loss rate and the wind velocity of  Ori, it could be blowing a stronger wind and still have the same size, as long as the compression ratio at the forward shock is more important.In the context of a radiative shock, the compression ratio can reach high values > 10 as suggested in Cardillo et al. (2019); however, the magnitude of the ISM magnetic field is also unconstrained, which could lead to diminishing effects at the forward shock (Meyer et al. 2017(Meyer et al. , 2021)), thereby altering the re-acceleration mechanism of pre-existing particles.Further, more detailed simulations are necessary to fine-tune these particular situations and better understand the excess -ray emission from massive star wind bubbles.

The X-rays surroundings of 𝜅 Ori
The overall surroundings of  Ori is a complex environment that is accessible for study through multi-wavelength observations with high-angular resolution data acquisition techniques.In particular, the study by Pillitteri et al. (2016) presents data obtained with the X-ray Multi-Mirror (XMM)-Newton space-based observatory, operating in the soft X-ray waveband, of the immediate vicinity of  Ori.This work proposes that a large and dense shell of swept-up material has formed around the central massive star.The X-ray fluxes are estimated to be in the range of 0.8 − 7.01 × 10 −14 erg s −1 cm −2 ≈ 6.24 × 10 −12 −4.43 × 10 −11 GeV s −1 cm −2 , and are interpreted as originating from the nearby hundreds of low-mass protostars that have formed in a dusty ring around  Ori.These protostars were found to obscure their reflected infrared starlight emission on the dust particles in past near-infrared observational surveys.Since the X-ray fluxes of the  Ori region are much larger than those that we predict, we conclude that our modeling is consistent with the observations, suggesting that the real data may consist of two components: one from the young stellar objects of the dusty ring and another from the re-accelerated cosmic rays, the latter being several orders of magnitude fainter than the former.This additional element supports the conclusions of Cardillo et al. (2019), stating that the circumstellar medium of  Ori is a cosmic-ray re-accelerator.

The radio synchrotron emission of the surroundings of 𝜅 Ori
The numerical simulations presented in this study predict nonthermal emission produced by the synchrotron mechanism (see panel 1a of Fig. 8).In the energy range 10 −12 -10 −7 GeV, which falls within the radio waveband (mm and longer wavelengths), these luminosities reach 10 31 GeV s −1 , corresponding to fluxes on the order of 1.0-1.1 × 10 −10 GeV s −1 cm −2 , or, in other cgs units, 1.60 × 10 −13 erg s −1 cm −2 .Since synchrotron radio emission has been detected from the Orion region, this allows for a qualitative comparison between such measurements and our predictions.These multi-wavelength observations have, among others, been performed with the Very Large Array (VLA) and the Atacama Large Millimeter/submillimeter Array (ALMA) in the radio waveband, reporting variability in the light curves of young stellar objects (Massi et al. 2006;Forbrich et al. 2008;Rivilla et al. 2015;Forbrich et al. 2017;Vargas-González et al. 2021), providing a high-energy non-thermal counterpart to the pre-main-sequence flares affecting protostars, whose origin has been proposed within the burst mode of accretion in star formation scenarios (Vorobyov & Basu 2006;Vorobyov 2009;Vorobyov & Basu 2010;Vorobyov et al. 2018;Meyer et al. 2019Meyer et al. , 2021Meyer et al. , 2022)).Additionally, several other radio sources in Orion, observed with the Arecibo telescope and the VLA, have been constrained as the outcome of synchrotron emission mechanisms occurring at shocks generated by outflows from higher-mass stars, such as a Herbig-Haro object (Yusef-Zadeh et al. 1990), the edge of the Orion-Eridanus superbubble (Bracco et al. 2023), or a more complex structure, e.g.consisting of relic supernova remnants and bow shock nebulae around runaway pulsars (West et al. 2022).The interpretation of the synchrotron radio emission in the Orion region as originating from particles accelerated at the edges of protostellar jets or other circumstellar shocks, like the expanding forward shock of an old supernova remnant or a stellar wind bubble, is not contradictory to the alternative scenario in which pre-existing cosmic rays would be re-accelerated, or even to a third hypothesis in which both mechanisms occur simultaneously.Hence, our model is consistent with the presence of synchrotron radio emission in the Orion region and might contribute, at least in part, to their explanation as originating from the shocked environments of massive stars like  Ori.

Comparison to stellar wind bow shocks
This work is consistent with the interpretation of observational results on stellar wind bubbles and stellar wind bow shocks who either reported or claimed the importance of cosmic-ray acceleration in massive stellar surroundings, because of the acceleration of particles at the termination shock of circumstellar structures such as stellar wind bow shocks (del Valle & Romero 2012, 2014;del Valle et al. 2015;del Valle & Pohl 2018).Particularly, the bow shock EB27 around the massive OB star BD+43 • 3654 has been predicted and shown to be a non-thermal emitter (Benaglia et al. 2010;del Valle et al. 2013;del Palacio et al. 2018;Sánchez-Ayaso et al. 2018;Benaglia et al. 2021), although other searches for non-thermal emission in stellar surroundings have been unfruitful (Schulz et al. 2014;Toalá et al. 2016bToalá et al. , 2017;;De Becker et al. 2017;Binder et al. 2019), suggesting that if the mechanisms of particle re-/acceleration is at work therein, it must be of lower efficiency than in other astrophysical sources such as AGN, pulsars winds, or supernova remnants.BD+60 • 2522 is an O-type massive star driving the nebula NGC 7635, well-known for its multi-wavelength emission, such as optical and infrared (Moore et al. 2002).Its stellar wind power and its space velocity are weaker than that of BD+43 • 3654, however, its ambient medium is denser than that of BD+43 • 3654, so that the numerical work of Green et al. (2019) suggested first the Bubble Nebula could be a potential non-thermal emitter.Previous obser-vational campaigns of the XMM-Newton space-born telescope of the Bubble Nebula revealed neither thermal nor non-thermal X-ray emission, see Toalá et al. (2016aToalá et al. ( , 2020)).The study of Moutzouri et al. (2022) deepens this effort using original interferometric radio continuum data at frequencies 4-12 GHz with the Very Large Array, rectified with single dish polarimetric observation at 4-8 GHz using the Effelsberg radio.The observational studies of Moutzouri et al. (2022) thus revealed that a second stellar wind bow shock around the surroundings of the runaway massive star NGC 7635 also emits detectable non-thermal radiation.
The nature of the non-thermal emission is determined by mapping the spectral index of the bow shocks and by comparing the flux of their spectral energy distribution with one-zone models for the acceleration/compression of wind/ISM non-thermal particles.Two different mechanisms have been invoked when considering non-thermal emission of stellar wind bow shocks of massive runaway stars: on the one hand, the emission of particles accelerated at the termination of the circumstellar medium and emitting in the shocked magnetic field (del Valle & Pohl 2018), on the other hand, the re-acceleration of pre-existing particles of the ISM at the forward shock of the star's bow shocks (Cardillo et al. 2019;Moutzouri et al. 2022).The second process is consistent with the surroundings of BD+43 • 3654, while unrealistic parameters at the forward shock are necessary to explain those of BD+60 • 2522.The first process could not fit the data for both bow shocks.Since bow shocks are compact nebulae with strong shocks, the ambient medium of BD+43 • 3654 has a density of  ≈ 15 cm −3 , which is less than that measured for  Ori, and both wind terminal velocities are similar (≈ 2020 km s −1 ), but its mass-loss rate is larger (  = 9 × 10 −6 M ⊙ yr −1 ), one can suggest that the radiation of accelerated particles might not be the appropriate mechanism to explain the -ray excess of the wind bubble of  Ori.This supports the conclusions of  Ori being a cosmic rays re-accelerator Cardillo et al. (2019).

CONCLUSION
In this study, we explore the possibility of shocks in a stellar wind bubble (Weaver et al. 1977) forming around very massive stars to be the site of particle acceleration and to generate non-thermal emission (Casse & Paul 1980;Voelk & Forman 1982;Webb et al. 1985).The wind-ISM interaction shaping the circumstellar medium of a 60 M ⊙ (Groh et al. 2014) is modeled by means of 1D MHD numerical simulations using the pluto code (Mignone et al. 2007(Mignone et al. , 2012;;Vaidya et al. 2018), and further post-processed with the particle acceleration code ratpac (Telezhinsky et al. 2012a(Telezhinsky et al. ,b, 2013)).The 1D MHD simulation account for the self-consistent time-dependent evolution of the stellar surface magnetic field, based on the Parker solution (Parker 1958) and on the stellar properties of Groh et al. (2014).The MHD structure of the wind bubble is further analyzed by a particle acceleration code, within the test particle approximation, which runs a shock-finder to catch its forward shock and calculate the acceleration of particles there.We obtain the time-dependent non-thermal particle density of accelerated electrons and protons at the shock, together with their radiation properties by synchrotron, inverse Compton and hadronic emission, respectively.
We show that the forward shock of a young stellar wind bubble around an early-type main-sequence very massive star, that is expanding into its local ambient medium, is able to re-accelerate pre-existing Galactic cosmic rays.The shock at the outer young bubble is sufficiently compressed to have a compression ratio ⩾ 4 and to accelerate the non-thermal electrons and protons, which have diffused from the ISM to the vicinity of the forward shock.The stellar wind bubble radiates via synchrotron, inverse Compton mechanisms and emits -ray emission by   decay.The emission flux spans from the sub-eV to the TeV energy bands, with a dominant component of the hadronic component at the GeV energy band and are multiplied by a factor of 5 during the early expansion of the bubble.This process is at work in the initial expansion phase of the bubble, since the ISM magnetic field damps the shock at later times and the shock ceases to be radiative at later times.
Our results are qualitatively in accordance with the enhanced -ray emission at about 1 GeV measured from the surroundings of the massive star  Ori, This emission has been shown to be inconsistent with diffusive shock acceleration at the forward shock of its circumstellar bubble, but compatible with the re-acceleration of ISM cosmic rays (Cardillo et al. 2019).This phenomenon likely applies to all bubbles of static massive stars, as well as the bow shocks of runaway stars (Moutzouri et al. 2022).Therefore, such a process should not be omitted when seeking to better understand the Galactic cosmic-ray spectrum (Cardillo 2019).This study needs to be extended to a broader scope, exploring both the stellar properties and the local conditions of the ISM.

Figure 1 .
Figure 1.Number density of a 1 Myr old stellar wind bubble, formed by a 60 M ⊙ star (black lines in cm −3 ).The figure highlights its principal discontinuities (termination shock, contact discontinuity, forward shock, dashed red lines), but also how the Böhm-like diffusion coefficient  ( ) of Eq. 19 is treated as a function of the distance to the star, into the different regions of the circumstellar medium.

Figure 2 .
Figure 2. Time evolution of the energy threshold  th used to cut off the Galactic cosmic-ray background, for electrons (a) and for protons (b).The quantity are plotted as a function of time (in Myr) for different multiplying factor  of the Böhm diffusion coefficient  B .

Figure 3 .
Figure 3. Number density (solid blue line, in cm −3 ) and magnetic field strength (dashed red line, in B) in our simulation of the stellar wind bubble of a 60 M ⊙ star, for different times of its evolution, from 0.01 Myr (a) to 1.0 Myr (c).

Figure 4 .
Figure 4. Properties of the forward shock of the stellar wind bubble of a 60 M ⊙ star.The panels show the shock position  F (in pc), the shock speed  F and the compression ratio of the gas at the forward shock  in the time interval between 0.4 Myr and 1.0 Myr.

Figure 5 .
Figure5.Distribution of the non-thermal particles as a function of momentum (electrons, left panels), and (protons, right panels) as a function of the adopted diffusion coefficient in the ambient medium, controlled by the parameter  .The number density of non-thermal particles (in cm −3 ) is shown as a function of the normalized momentum, for several time instances, 0.5 Myr (dotted lines), 0.7 Myr (dashed lines) and 1.0 Myr (solid line).The blue lines are the initial Galactic cosmic-ray background and the orange vertical line is the momentum threshold.

Figure 6 .
Figure 6.Distribution of the non-thermal particles as a function of momentum (electrons, top panel; protons, bottom), as a function of the adopted diffusion coefficient in the ambient medium, for  = 10 (thick solid red line),  = 200 (dotted black line) and  = 400 (thin solid dark green line).The number density of non-thermal particles (in cm −3 ) is shown as a function of the normalized momentum, at a time instance corresponding to the end of the simulation at time 1.0 Myr.

Figure 7 .
Figure 7. Evolution of the distribution of non-thermal protons (in cm −3 ) for the threshold (top) and maximum (bottom) energies, at selected times of the wind bubble evolution, in the diffusion model with  = 200.The black line is the Galactic cosmic ray background.

Figure 8 .
Figure 8. Time-sequence evolution of the non-thermal emission spectra of the forward shock of the stellar wind bubble created by a 60 M ⊙ star during the early 1 Myr of its main-sequence phase.The luminosities are plotted as a function of the energy at times 0.5 Myr (left), 0.7 Myr (middle) and 1.0 Myr (right).The emission is displayed for the synchrotron (SY, top panels), inverse Compton (IC, second line of panels), -rays by   decay mechanisms (PD, third line of panels) and for all three processes (bottom panels).

Figure 9 .
Figure 9.Comparison between the relativistic Bremsstrahlung timescale (solid red line) and inverse Compton scattering of CMB (cosmological microwave background) photons (dashed blue line) timescale (in Myr), for the non-thermal emission at the forward shock of the stellar wind bubble around a massive star, 1 Myr after the onset of its central star main-sequence phase.

Figure 10 .
Figure10.Comparison between the -rays spectrum of the stellar wind bubble created by a 60 M ⊙ star after 1 Myr of main-sequence evolution, for the different adopted diffusion models.The plot units, frame box size and the orange dots are data from the AGILE -rays excess observed from the vicinity of the B-type star  Ori and reported in fig.3-4 ofCardillo et al. (2019).