Can superbubbles accelerate ultra-high energy protons?

We critically assess limits on the maximum energy of protons accelerated within superbubbles around massive stellar clusters, considering a number of different scenarios. In particular, we derive under which circumstances acceleration of protons above peta-electronvolt (PeV) energies can be expected. While the external forward shock of the superbubble may account for acceleration of particles up to 100 TeV, internal primary shocks such as supernova remnants expanding in the low density medium or the collective wind termination shock which forms around a young compact cluster provide more favourable channels to accelerate protons up to 1 PeV, and possibly beyond. Under reasonable conditions, clustered supernovae launching powerful shocks into the magnetised wind of a young and compact massive star cluster are found to be the most promising systems to accelerate protons above 10 PeV. On the other hand, stochastic re-acceleration in the strongly turbulent plasma is found to be much less effective than claimed in previous works, with a maximum proton energy of at most a few hundred TeV.


INTRODUCTION
The search for Galactic cosmic-ray (CR) sources able to accelerate above PeV energies has been ongoing for decades. The importance of this energy range was established early on, as it presented an obstacle to CR origin theories via diffusive shock acceleration (DSA) at supernova remnant (SNR) shocks (Lagage & Cesarsky 1983b;Heavens 1984;Axford 1994). For standard ISM magnetic field values, and assuming the most optimistic conditions for transport of CRs near the shock (i.e. the Bohm diffusion limit) the maximum energy for typical SNRs was shown to fall close to 100 TeV (e.g. Cesarsky & Lagage 1981;Lagage & Cesarsky 1983a). Such low maximum energies could not account for the knee feature in the CR spectrum at a few PeV, let alone the steeper second Galactic population extending above the knee to ∼ 10 18 eV energies. Alternative scenarios, in particular the acceleration around massive star wind-termination shocks, were shown to barely accelerate up to PeV (Casse & Paul 1980;Cesarsky & Montmerle 1983a), and to suffer several issues such as the problem of the injection and confinement at a perpendicular shock (the magnetic field being azimuthal far away from the massive star).
Significant progress followed the demonstration by Lucek & Bell (2000), and in particular Bell (2004), that CRs could drive substantial non-linear amplification of magnetic fields around young SNR shocks. This scenario of particle acceleration however requires rather special conditions to reach PeV bands (e.g. Zirakashvili & Ptuskin 2008;Bell et al. 2013;Gabici et al. 2019;Cristofari 2021) and lacks observational confirmation (though see HESS Collaboration et al. 2022).
To date, a convincing scenario for the production of multi-PeV ★ E-mail: thibault@mpi-hd.mpg.de energy protons in our Galaxy is lacking. Such extremely energetic CRs are not only needed to reproduce the CR Galactic-extragalactic transition as measured by ground-based detectors (Parizot 2014), but also to account for recent UHE gamma-ray observations (Cao et al. 2021a;Cao et al. 2021b). Besides, a number of sources are detected up to hundreds of TeV with no significant cut-off within the detection range of the instrument (e.g. Abramowski et al. 2012a;HESS Collaboration et al. 2016;Abeysekara et al. 2020;Abeysekara et al. 2021). Interestingly, a number of these observations are correlated with massive star clusters or OB associations born within star-forming regions. It was already pointed out in the early 80s that massive stars are often found in clusters and that collective effects between stellar winds and possibly SNRs (referred to as "SNOBs" at that time) could play a key role in the acceleration of high energy CRs (e.g. Montmerle 1979;Cesarsky & Montmerle 1983b). Although the SNR origin of Galactic CRs has been extensively explored, the peculiarities of the environment within or around a cluster or association have been scarcely considered.
Here we discuss the importance of massive stellar clusters (MSCs) which form within dense molecular clouds. The stellar feedback from the cluster heats the surrounding medium which inflates a cavity around the cluster. The cavity keeps expanding during the whole cluster lifetime. After several Myr, it reaches a size of about 100 pc and is commonly referred to as a superbubble. A number of arguments support the superbubble (SB) origin of CRs (see e.g. Lingenfelter 2018; Bykov et al. 2020, for recent reviews). In particular these environments are expected to excite a very turbulent plasma which effectively confines the particles while a number of shocks sweep up the medium and accelerate the particles, such that a substantial fraction of the stellar power can be efficiently transferred into CRs (Montmerle 1979;Bykov 2001). This was confirmed in the last decade by an increasing number of gamma-ray observations (e.g. Ackermann et al. 2011;Abramowski et al. 2012b;H. E. S. S. Collaboration et al. 2015;Yang et al. 2018;Aharonian et al. 2019;Abeysekara et al. 2020;Mestre et al. 2021). This class of sources have generically a spectrum compatible with the local CR flux when propagation effects are taken into account (d /d ∝ −2.2...−2.5 ). Besides, the composition of the plasma around massive star clusters or associations is expected to differ from that of the standard interstellar medium as it is enriched by the winds of Wolf-Rayet (WR) stars. This may explain a number of composition anomalies measured in the local CR flux (e.g. Gupta et al. 2020;Tatischeff et al. 2021).
The LHAASO observatory has reported the detection of a PeV photon coincident with the Cygnus region (Cao et al. 2021a), implying CRs of several PeV are produced therein. The maximum energy of CRs produced within MSCs, OB associations or SBs is expected to increase compared to the case of an isolated SNR, because the stellar outflows should in principle be able to excite a strong turbulent magnetic field, which is needed to confine the particles. This follows from the Hillas criterion (Hillas 1984), which is commonly invoked to set an upper bound on the maximum energy of the particles accelerated in a given source as max < / , where , and are respectively the typical velocity, magnetic field and size of the source. Assuming ≈ 10 µG the magnetic field inside the SB, ≈ 100 pc the size of the SB and ≈ 3000 km/s the typical velocity of the SNR shocks and stellar winds (Bykov et al. 1995;Bykov & Toptygin 2001), one obtains max ∼ 10 PeV. This simple estimate has long motivated the candidacy of SBs as sources of super PeV protons.
This estimate however relies on the implicit assumption that the particle accelerator is operating on the full scale of 100 pc. This has some support from gamma-ray observations (e.g. in the Cygnus region, Abeysekara et al. 2021) which display an extended emission over several tens of pc. However, it is also very possible that CRs are accelerated in a subsystem, e.g. a SNR expanding within the SB, and only produce gamma-rays after they have diffused within the cavity and reached the dense and strongly magnetised shell which surrounds the SB. This demonstrates the importance of distinguishing SBs as gamma-ray sources and SBs as particle accelerators. In particular, it is not clear which acceleration mechanism would take place over the whole SB scale, or if there is such a mechanism at all.
The present paper aims at clarifying this question, considering a number of detailed scenarios for the acceleration of particles within SB environments (including compact MSCs or looser stellar associations). The Hillas limitation is justified by theoretical arguments which motivate credible values for the relevant source size, velocity, and magnetic field. In Section 2 we anticipate the results using simple estimates in order to get benchmark values for the maximum proton energies. In Section 3 we detail the properties of the SB, in particular how its size, pressure, density, mean magnetic field and mean turbulence velocity evolve. Using physically motivated acceleration mechanisms, we then compute directly the maximum energy considering three candidate accelerators: the SB forward shock (Section 4), the supernovae and wind termination shocks (Section 5), and the turbulence (Section 6). In each case, we discuss the possibility of accelerating UHE (about 10 PeV) protons. We conclude in Section 7.

BENCHMARK ESTIMATES
The Hillas criterion is sometimes formulated as a confinement criterion: max ≈ , i.e. particles with Larmor radius larger than the size of the accelerator escape. This is the least stringent constraint one can put on the maximum energy and it does not give any clue about the acceleration mechanism. As particles are accelerated in electric fields (or equivalently the motion of scattering centres), considering a specific acceleration mechanism introduces a factor / in the above estimate, where is some characteristic velocity specific to the acceleration mechanism under investigation. The maximum energy is then that achieved in the electric potential Φ = / , such that max ≈ / . It should be noted that this should still be understood as an upper bound, as other processes, such as losses, may decrease the maximum achievable energy.
While the confinement criterion is rather universal, introducing the "characteristic velocity of the source" is more challenging as it is not always a well-defined quantity. In order to estimate its value, one needs to specify the mechanism driving the acceleration of particles. In this introductory section, we consider a number of possible mechanisms which may work within SBs, in order to derive preliminary limitations on the maximum energy. Each mechanism will be considered in more detail in the following sections.
SBs are complex environments. The hot rarefied plasma is spanned by primary shocks, which decay into turbulent motions and MHD waves. Collective effects such as re-acceleration processes contribute to the acceleration of particles. SBs are also delimited by a forward shock which expands into the ISM. While the size of this shock can be very large ( 100 pc, Weaver et al. 1977), it is too slow ( ∼ 30 km/s) to accelerate PeV protons, even assuming large magnetic fields in the ISM. More generally, in the case of strong shocks, the velocity is identified as the velocity of the shock, as it is precisely the velocity jump at the shock discontinuity which drives the acceleration of the particles via the DSA mechanism (e.g. Drury 1983). Inside SBs, there are two types of strong primary shocks: the time-dependent SNR shocks which expand after a supernova (SN) explosion, and the wind termination shocks (WTS) which surround the individual massive stars, or the entire stellar cluster if it is compact enough. SNRs and WTSs have a typical velocity of several 1000 km/s. In the low-density SB interior, SNRs expand to a typical radius of a few tens of pc before reaching the Sedov-Taylor phase, which is generally larger than the radius of the WTS, even in the case of a WTS powered by a very massive compact cluster. The size of the latter shock depends on the mechanical power of the stellar cluster and is typically of the order of 10 pc (Weaver et al. 1977). It is believed that, due to in situ CR acceleration, efficient magnetic field amplification takes place upstream of SNR shocks, leading to a magnetic field of up to several tens of µG. On the other hand, it is less clear if the magnetic field can be as efficiently amplified in stellar winds. In both cases, the maximum energy is generally inferred to be only ∼ 1 PeV (e.g. Gupta et al. 2020;Morlino et al. 2021;Vieu et al. 2022b). While this is slightly larger than the maximum energy achieved at isolated SNR, atypical conditions would be required for protons to be accelerated well beyond PeV by primary shocks embedded in SBs. A promising situation combining the advantages of both WTS and SNRs might be that of a SNR expanding in a wind profile close to a compact cluster. A powerful cluster may indeed convert a substantial amount of its mechanical energy into turbulence, amplifying magnetic fields up to hundreds of µG in its vicinity, such that particles could be accelerated up to 10 PeV by a powerful SNR propagating in the wind, even in the absence of additional magnetic-field amplification. Considering a similar scenario of SNR shock -cluster wind interaction with efficient turbulence generation, Bykov et al. (2015) found that proton energies up to 40 PeV could be reached in the case of fast shocks ( = 10 4 km/s).
The winds and SNRs collectively sustain turbulence in the SB. Supersonic hydrodynamic turbulence, composed of random largescale motions of the plasma, provides compression and rarefaction waves which may also be able to accelerate particles. Indeed, the random waves are scattering centres and the original stochastic Fermi acceleration mechanism takes place (Skilling 1975;Ptuskin 1988;Bykov & Toptygin 1990). In this case, the velocity should be identified as the mean velocity of the flows, that is, the square root of the energy stored in the turbulence. This provides typically ∼ 100 km/s, roughly an order of magnitude smaller than the velocity of primary shocks. Even if, in this case, the turbulence is distributed inside the whole SB, such that the relevant size of the accelerator is of the order of 100 pc, the expected maximum energy is estimated as a few hundreds of TeV. If a substantial fraction of the mechanical power of the cluster goes into magnetic field amplification, 1 PeV could be reached, although this is already an optimistic value. Finally, an ensemble of individual winds, gathered in a loose cluster, that is, a cluster around which no collective WTS forms, may be viewed as an ensemble of stochastic shocks from the point of view of high energy particles and act as scattering centres, leading to a stochastic acceleration process over strong primary shocks (Bykov 2001). Although the velocity of the scattering centres is of the order of the wind velocity, the scattering cross-section is rather small, i.e. particles diffusing inside the SB have a low probability of encountering an individual WTS multiple times. Indeed, the size of an individual WTS around a massive star is about 1 pc, that is at least several times smaller than the radius of a loose stellar cluster. The average plasma velocity in the stellar cluster is correspondingly about (1 pc/ ) 3 , where is the radius of the cluster and ∼ 2000 km/s the average wind velocity. In the case of an extended cluster, 5 pc and the average velocity is of the order of a few tens km/s. The maximum energy achieved by protons via this acceleration mechanism is not expected to exceed 500 TeV. In a linear treatment of the problem, and assuming optimistic parameters ( = 3000 km/s, = 30 µG, over a region of radius = 5 pc), Klepach et al. (2000) showed that PeV energies could indeed be hardly reached via this mechanism of stochastic acceleration. It is furthermore difficult to justify how an average velocity of 3000 km/s could be maintained within an ensemble of stars surrounded by individual wind termination shocks, each of them having a radius much smaller than 5 pc. Table 1 summarises the general expectations on the maximum energy which have been discussed in the present Section. The remainder of the paper is devoted to a more detailed analysis of the possible acceleration mechanisms in the peculiar SB environment.

SUPERBUBBLE PROPERTIES
To motivate the calculations that follow, we introduce the key assumptions of SB environments. We consider a massive stellar cluster which has formed inside a dense molecular cloud. Massive stars (which are defined in the following as the stars with mass > 8 ), emit powerful winds during their lives. If the cluster is compact enough, a collective WTS will form around the MSC. Otherwise, each individual star is surrounded by its own WTS. Massive stars explode as SNe at the end of their lives and turn into expanding SNRs. SNRs and WTS heat the surrounding medium, which inflates a low-density cavity called a superbubble (SB). Thermodynamical considerations provide the expansion law of the SB as (Castor et al. 1975;Weaver et al. 1977): where P is the average mechanical power of the stellar cluster, 0 is the hydrogen number density of the surrounding medium, which typically ranges from 1 cm −3 (the average ISM) to 100 cm −3 (giant molecular clouds), is the age of the MSC and ∼ 0.1 a power conversion factor which accounts for the fact that not all the stellar power is converted into superbubble shell expansion (e.g. Yadav et al. 2017;Vieu et al. 2022b). In deriving the expression above, it was assumed that is much larger than the radius of the star cluster and that SNRs become subsonic before reaching the edge of the SB. The average mechanical power of the MSC is roughly constant in time and can be written as P = 10 36 * erg/s, where * is the number of massive stars ( > 8 ) in the MSC (Vieu et al. 2022b). This includes the main sequence (MS), Wolf-Rayet (WR) and supernova remnant (SNR) phases of the stellar evolution (the slow winds of the red supergiant phase have a negligible contribution Seo et al. 2018).
The MSC not only blows the SB, it also drives the generation of large-scale MHD turbulence via the decay of primary shocks (e.g. Bykov & Toptygin 1987;Padoan et al. 2016) or various instabilities operating at a macroscopic scale (e.g. Inoue et al. 2009), along with CR-driven instabilities (e.g. Marcowith et al. 2016, and references therein). Without delving into details of the turbulence injection mechanism, one can write a simple phenomenological equation for the generation of turbulence as (Zhou & Matthaeus 1990): where ( ) is the turbulent spectral energy density at wavenumber , is the spectral energy flux from the large scales to the small scales and the source of the turbulence at the largest scale 2 / 0 . Assuming that the magnetic and hydrodynamic energies are in equipartition, the turbulence spectrum can be normalised such that ∫ d ( ) = 2 /(4 ) =¯2, where and¯are respectively the turbulent magnetic field and turbulent velocity averaged over the perturbations (or, equivalently, over the volume in the case of homogeneous turbulence). Assuming that the stellar feedback is homogeneous over a volume , one can write the source term in this volume as (Vieu et al. 2022b): where < 1 is the efficiency of turbulence generation. Numerical simulations suggest 30% (Gallegos-Garcia et al. 2020). Further numerical work is needed to better constrain this parameter. In the following numerical estimates, we will assume a reasonable value = 10%. Solving Eq. 2 provides typically: where is the average density in the volume . Depending on the nature of the turbulence (i.e. the expression of the flux in Eq. 2), this estimate varies by a small numerical factor of order unity. The average density and temperature in the SB interior are mainly driven by the evaporation and mass loading at the shell (Weaver et al. 1977). Thermodynamical calculations provide (Mac Low & McCray 1988): The scalings derived in this section describe the thermodynamic, kinetic and magnetic properties of the SB, which we will now use to estimate the maximum energy of the protons accelerated in SB cavities, considering a number of possible scenarios.

PARTICLE ACCELERATION BY THE SUPERBUBBLE FORWARD SHOCK
SBs are delimited by a forward shock, whose size and velocity are readily given by Eq. 1. Even assuming optimistic values for the mechanical power (P = 10 39 erg/s) and ISM density ( 0 = 1 cm −3 ), the forward shock rapidly becomes weak, e.g. ≈ 65 km/s at = 1 Myr. In the case of a weak shock, the maximum energy achieved by the particles via DSA is generally limited by the available acceleration time. Assuming that the diffusion coefficient is homogeneous around the SB shock, the acceleration rate reads (Lagage & Cesarsky 1983b): where is the shock compression ratio, with M the shock sonic Mach number and the adiabatic index of the medium. In order to find the maximum momentum, we assume Bohm diffusion, = /( ), with ≤ 3, and integrate the acceleration rate. The Larmor radius of particles at the maximum energy (or rigidity) is plotted in Figure 1 (solid curves) as function of the age of the SB for optimistic values of the parameters (see caption).
For completeness we also plot on Figure 1 the size limitation on the maximum energy. Indeed, particles whose diffusion length is larger than the shock radius are expected to escape upstream of the shock. Assuming Bohm diffusion, this provides another limitation as ,max = / . This limitation is generically less stringent than the time limitation discussed above.
Even in the case of a very massive cluster (P = 10 39 erg/s -about 1000 massive stars), the maximum Larmor radius never exceeds 0.02 pc. With more reasonable parameters, e.g. 0 10 cm −3 , it falls below 0.01 pc.
In order to estimate the maximum energy from the maximum = 10 38 erg/s = 10 39 erg/s Figure 1. Evolution of the maximum Larmor radius achieved in the SB forward shock around a canonical cluster (P = 10 38 erg/s, blue curve) and around a very massive cluster (P = 10 39 erg/s, red curve). The solid curves show the limitation due to the finite age of the SB while the dashed curves show the limitation due to the finite size of the SB. We assumed optimistic parameters: = 1, 0 = 1 cm −3 , and an ISM temperature of 10 4 K.
Larmor radius, one has to specify the magnetic field. If the magnetic field experienced by the particles diffusing around the forward shock is self-generated via non-resonant streaming instability, one expects ∝ 1/2 3/2 and therefore decreases with time (Bell 2004). In fact the self-generated magnetic field is already negligible at early times as the shock is too slow ( ≈ 0.05 µG for * = 1000 and 0 = 1 cm −3 ). The fact that after a few Myr the radius of the forward shock becomes of the order of 100 pc is irrelevant in this scenario.
This shows that particles diffusing upstream of the forward shock likely experience the ISM background magnetic field rather than the self-generated magnetic field, and thus the Bohm assumption is already optimistic. The ISM magnetic field is typically a few µG such that even in the most optimistic case displayed in Figure 1, the particles are barely accelerated up to 100 TeV. The ISM magnetic field may be higher in peculiar regions of our Galaxy, such as the Galactic centre, but the ambient density is then also higher than 1 cm −3 , which slows down the SB expansion. The SB forward shock is therefore generically unable to accelerate PeV protons above 100 TeV, because it rapidly becomes very weak.

PARTICLE ACCELERATION BY INTERNAL PRIMARY SHOCKS
As the large-scale collective shock surrounding the SB is unable to accelerate protons to ultra high energies, we now consider the acceleration at internal shocks, which are smaller but faster. Two cases must be distinguished. We first discuss the case of a "loose" cluster, where the separation between the clustered stars is too large for a collective outflow to form. In this case, the individual WTS surrounding each massive stars are too small to accelerate very high energy CRs and the only relevant primary shocks are those associated to SNRs. We then discuss the case of a "compact" cluster, where the stars are confined in a small region, allowing the formation of a large collective WTS. In this case we consider the acceleration of particles around the WTS itself, and then at SNR expanding in the "free wind" upstream of the WTS.

Supernova remnants in loose clusters
As stated above, the average stellar separation in a loose cluster is too large for a WTS to form. In this case, the SNRs expand in the low-density cavity for times longer than typical in standard ISM conditions (Parizot et al. 2004). In the resulting turbulent medium, it is reasonable to assume a Bohm diffusion regime with a uniform diffusion coefficient. Magnetic fields can be efficiently amplified via two mechanisms: the generation of a MHD turbulent cascade from the largest scales where the massive stars inject a fraction of their kinetic power (Eq. 4), or CR driven instability due to strong CR currents (Bell 2004). Neglecting the latter for the moment, the maximum energy, achieved at the end of the free expansion phase, can be estimated (Vieu et al. 2022b where ej is the mass of the ejecta. Using the magnetic field and density calculated in Section 3, we plot in the left panel of Figure 2 the evolution of the maximum momentum as a function of the cluster age. We adopt an ISM density of 1 cm −3 and turbulence generation efficiency = 10%, as stated above. The maximum momentum decreases with time as the SB expands and the turbulence is diluted over a larger volume. For these generally optimistic parameters, PeV energies are approached, though more typical values will provide max ∼ 0.3 − 1 PeV, similar to what is obtained in isolated SNR shocks. Let us then account for the self-generation of the magnetic field upstream of the SNR shock. Under favourable conditions, this may enhance the maximum energy in the case of isolated SNR. It is indeed expected that a fraction of the shock pressure will be converted in magnetic turbulence due to CR-driven instabilities (e.g. Dorfi & Drury 1985;Bell 2004;Beresnyak et al. 2009;Downes & Drury 2014;Xu & Lazarian 2017;Zweibel 2020, etc.). We write: where ∼ 10 % is the acceleration efficiency. Plugging this into Eq. 8 provides: This is plotted in the right panel of Figure 2 for 0 = 100 cm −3 and = 10%. Interestingly, in this case, max weakly depends on the density inside the SB, and thus is almost constant throughout the SB lifetime.
This shows that even though the CR-driven instabilities are less efficient in a low-density medium, PeV protons could be generically expected from SNR shocks expanding inside SBs. However the maximum proton energy is not found to reach 10 PeV, even around very fast SNR shocks. Besides, as the cavity expands, the turbulence level decreases: if CR-driven instabilities are not effective, PeV protons are only expected around young clusters ( 10 Myr). Note that from the scalings provided in Section 3, one can show that the magnetic field is more efficiently generated from the stellar outflows if the ISM density is low. On the other hand, if the magnetic field is generated via CR-driven instabilities, a high ISM density is more favourable. Therefore, the two mechanisms of magnetic field amplification do not generally apply simultaneously.

Wind termination shock (WTS)
We now consider the case of a compact cluster. Before investigating the properties of the SNRs exploding within compact clusters, we discuss the WTS which surrounds the cluster. Indeed, the outflows expelled by massive stars gathered in a compact cluster pile up onto a collective shock around the cluster. The size of the WTS is determined by the balance between the wind pressure upstream of the WTS and the pressure inside the SB downstream of the WTS. The latter is derived from the average temperature and pressure inside the SB as given by Eqs. 5 and 6. Using the fits of stellar evolution models provided by e.g. Seo et al. (2018), we compute the wind mechanical power, and thus the wind pressure, during the entire lifetime of 1000 random MSC samples generated using a Salpeter initial mass function (Salpeter 1955) for 100 massive stars. The evolution of the WTS radius averaged over the samples is shown in Figure 3.
The radius of the WTS reaches a maximum at 3 Myr, then decreases as the most massive stars start to explode, and eventually drops dramatically after 10 Myr, since no WR stars remain in the MSC after this time. At 3 Myr, the wind mechanical power is comparable to the energy deposited by supernova remnants (about 10 51 erg/Myr for a cluster of 100 massive stars). An upper limit on the size of the wind termination shock can be computed from a pressure balance as (Vieu et al.  pc . Note that the collective WTS will disappear if its hypothetical radius becomes smaller than the extent of the MSC. This typically happens once all WR stars have exploded, after 10 Myr. The acceleration mechanism discussed in this section is therefore only applicable in general to young and compact MSCs. The dynamical timescale of the WTS is much larger than the DSA acceleration time. Therefore the maximum energy achievable by particles accelerated at the WTS is limited by the shock size. It was shown by Morlino et al. (2021) that the standard Hillas estimate gives an accurate zeroth-order estimate of the maximum energy. The maximum energy is achieved when the diffusion length becomes of the order of the size of the shock, for beyond this limit the particles escape in the upstream region. In the case of Bohm diffusion, = 2 /( ), the maximum momentum is then: where = 3 g / mfp , with mfp the mean free path due to magnetic deflections. The most optimistic case, the Bohm limit, corresponds to = 3. Even in the scenario where the ISM density is small ( 0 ∼ 1 cm −3 ) and the star cluster very massive ( * ≈ 1000), PeV protons can only be accelerated if the magnetic field is amplified to typically 10 µG. Given that the efficiency of CR-driven instabilities in the wind density profile is still an open question, we conclude that in the absence of other field amplification mechanisms, the proton energies are not expected to reach more than a few PeV, in the most optimistic estimate. This conclusion is in agreement with the in-depth computation of Morlino et al. (2021) and the simulations of Gupta et al. (2020). The magnetic field is also amplified by the generation of largescale MHD turbulence from the decay of primary outflows emitted by the stars (winds and SNRs within the cluster). As discussed in Section 3, it is expected that a fraction of the mechanical power of the MSC is converted into MHD turbulence within the compact cluster. A major difference with the case of a loose cluster is that the magnetic field is confined. Indeed, the magnetic field advected in the free wind region has a 1/ radial profile and the associated energy density is negligible. Let us assume that the stars are evenly distributed within a sphere of radius , which is the "extent" of the cluster. We then expect the turbulence to be homogeneous in the cluster. The turbulence scale is taken as the distance between the stars, which is = (2/3) 1/3 −1/3 * . The magnetic field in the cluster is estimated from Eq. 4 as: Because the magnetic field decreases as 1/ outside of the cluster core, the intensity of the magnetic field at the position of the WTS is inversely proportional to the radius of the WTS and the maximum momentum does not depend anymore on the size of the WTS, but rather on the cluster size: max ≈ 3 ( / ) which leads to the following estimate: This suggests that WTS around compact clusters can generically accelerate protons up to several PeV. Remarkably, the maximum energy weakly depends on the properties of the cluster, for instance its mass and extension, provided that it is compact and powerful enough to allow the formation of a collective WTS. The numerical values input in Eq. 14 correspond to those typically assumed for Westerlund 1, where recent HESS observations show a ring-shaped gamma-ray emission enclosing the cluster, which could theoretically correlate with the position of the WTS (Mohrmann et al. 2021). To accelerate protons beyond 10 PeV at a MSC WTS requires a very massive ( * ∼ 1000, i.e. a total cluster mass of about 5 × 10 4 ) extended ( ∼ 3 pc) cluster, with highly efficient turbulence generation. The only known clusters which displayed such properties are RSGC01, RSGC02, RSGC03 (Krumholz et al. 2019) but with an age above 10 Myr they are too old to display a WTS: the Wolf-Rayet stars have already exploded several Myrs ago.

SNRs in the cluster
SNRs exploding within young compact clusters are not directly launched in the low-density SB cavity. It is in fact yet unclear what happens to SNR shocks expanding deep inside the cluster. It is very possible that the shock, whose pressure and size are comparable with that of the stellar outflows of the most massive stars, will be disrupted before leaving the cluster. This is expected in particular to enhance the turbulence level and the magnetic fields in the core of the cluster (Eq. 13). Keeping these reservations in mind, let us nevertheless assume that a SNR shock is launched in the free wind region with a relative velocity at = (e.g. a massive star explodes close to the edge of the cluster). In the 1/ 2 wind profile, the ejecta-dominated SNR shock follows a self-similar evolution, ∝ , with ≈ 7/8 (e.g. Finke & Dermer 2012;Gaggero et al. 2018). Particles will be accelerated while the shock expands upstream of the WTS (Gupta et al. 2020). The maximum energy is limited either by the size or the age of the shock. Because the magnetic field decreases as 1/ in the free wind region, the size limitation gives: which decreases with time as the shock slows down. On the other hand, the time limitation on the maximum momentum is inferred by considering the acceleration rate (Eq. 7) in the case of a strong shock ( = 4). Assuming again Bohm diffusion and integrating on time leads to (for < 1): The maximum energy is achieved when both criteria 15 and 16 are equal (earlier the SNR is too young, later the SNR is too slow), i.e. when the SNR shock reaches a radius of: However, for = 7/8, * = 44 which is expected to exceed the radius of the WTS. Therefore the maximum energy will be reached when the SNR collides with the WTS and can be computed by setting = in Eq. 16. We get, assuming = 7/8: max ≈ 10 5000 km/s 3 1 − PeV (18) For reasonable cluster sizes, / ∼ 5 − 30, we get 1 − ( / ) 1/7 ∼ 0.2 − 0.4 and therefore a maximum energy of the order of several PeV for the parameters used in Eq. 18. Therefore, if a supernova launches a fast shock (with velocity relative to the wind 10 000 km/s) in the outer parts of a compact cluster which efficiently excites MHD turbulence to generate magnetic fields of about 100 µG in its vicinity, protons could be accelerated around the SNR beyond 10 PeV. Note that the collision between the SNR and the WTS is not expected to enhance the maximum energy (Vieu et al. 2020), unless non-trivial amplification effects arise, which is beyond the scope of the present work.
A recent SNR shock expanding in a collective wind around a MSC might be responsible for the acceleration of particles in e.g. Westerlund 1, the Quintuplet, [ 2003] 179. Other very young MSCs such as Westerlund 2, NGC 3603, or the Arches are a priori excluded as no supernova is expected to have exploded within the cluster yet.

Reacceleration by SNR shocks
Around massive star clusters, SNR shocks span the same region successively. The question of particle acceleration by successive timedependent shocks was investigated in Vieu et al. (2022a), where the maximum achievable energy was left as a free parameter. In general, the reacceleration of the particles around an isolated shock is not expected to dramatically enhance the maximum energy compared to the first injection step. Even in the case of interacting winds, efficient reacceleration requires non trivial confinement enhancement (Vieu et al. 2020). Eventually, when a SNR launched at the edge of a compact cluster converges toward the collective WTS surrounding the cluster, it will reaccelerate the particles pre-accelerated by the WTS, which may leave detectable spectral features such as bumps or hardenings in the TeV bands (Bykov et al. 2019;Vieu et al. 2020). After the collision between the SNR and the WTS, further particle acceleration may proceed around the transmitted and reflected waves. Although these secondary shocks are expected to be weak (Landau & Lifshitz 1987), they provide additional reacceleration time which may slightly enhance the maximum energy, by a factor ∼ 2.

Acceleration by compression and rarefaction waves
Inside the SB, the particles probe large-scale compression and rarefaction waves. This gives rise to a stochastic acceleration process (second order Fermi mechanism). In order to infer the maximum momentum achievable via this acceleration mechanism, we must investigate the transport of particles in a strongly turbulent medium.
In a strongly turbulent medium, particles diffuse due to their interactions with MHD waves (referred to as "small-scale diffusion"), but also due to their advection in the random large-scale flows. The characteristic diffusion, acceleration and escape times read: where is the largest (energy-containing) scale of the turbulence ( = 2 / 0 ), is the spatial diffusion coefficient, is the momentum diffusion coefficient and is the extent of the turbulent region. A critical parameter driving the acceleration of particles is the ratio of acceleration over escape time, denoted ≡ / . We have: In the following we seek an expression for . In order to do so, we compute the transport parameters and associated timescales from first principles following Bykov & Toptygin (1990, 1993. The starting point is the transport equation in the diffusion approximation: where is the total velocity field including all perturbations and is the small-scale diffusion coefficient which only includes the interaction between the particles and the magnetic waves, not the effect of the large-scale flows (in contrast to the total diffusion coefficient defined above). Eq. 21 describes the small-scale diffusion over MHD perturbations, the advection in the large-scale plasma flows and the energy change due to the interactions with compression and rarefaction waves. It is not possible to solve directly Eq. 21 within a semi-analytical framework. Instead, Eq. 21 must be first averaged over the HD fluctuations, in order to obtain a transport equation in terms of the mean velocity¯and not the fluctuating velocity . This average is not a trivial task. Indeed, when the small-scale diffusion length / is smaller than the largest eddy turnover length , the particle distribution function is expected to strongly fluctuate over small-scales and perturbation theory breaks down (within one eddy turnover length, one cannot write =¯+ with ¯) . One must then renormalise the transport coefficients in order to account for strong long-wavelength fluctuations. Bykov & Toptygin (1990) provides the result as (for isotropic turbulence): The interaction with compression and rarefaction waves leads, on average, to a diffusion in momentum space. Besides, the advection in the random large-scale flows leads, on average, to a spatial diffusion, which is included, along with the small-scale coefficient , in the transport parameter . The transport coefficients are calculated as (Bykov & Toptygin 1990): where and are respectively the longitudinal and transverse wave spectra to be specified. We chose the Fourier convention ( , ) = ∫ d d˜( , ) ( · − ) . In the following we assume a simple isotropic power-law turbulence spectrum: normalised such that 4 ∫ d d 2 ( , ) =¯2. is the spectral index of the turbulence energy spectrum (e.g. = 5/3 for a Kolmogorov spectrum, = 3/2 for a Iroshnikov-Kraichnan spectrum , = 2 in the case of supersonic turbulence), and [ ] is the Heaviside step-function.
Using the spectrum 25, one can compute numerically the spatial and momentum diffusion coefficients as given by Eqs. 23 and 24. The evolution of the parameter as function of the small-scale diffusion coefficient is shown in Figure 4  = 3/2 + 3/2 (1 + 4 /9) 1/2 .
When is momentum-dependent, the bulk of the particle spectrum accelerated via stochastic acceleration in turbulence is not a powerlaw and the maximum energy is ill-defined. In the following we determine the maximum momentum as that producing ∼ 10. Above this value the differential energy spectrum becomes steeper than −3 , which is at the limit of what has been observed in terms of steep gamma-ray spectra (e.g. Cao et al. 2021b). The plot on the right panel of Figure 4 is the same plot as that on the left, but where = 10 was fixed, and the axis reversed. This readily shows the value of above which the acceleration stops, as function of the geometry of the system. If increases monotonically with momentum, the maximum momentum will always be achieved in the linear regime, as we expect < in superbubbles. Similar conclusions hold when one considers more realistic turbulence spectra than the ansatz (25), e.g. with an exponential decay of the time correlations and a smooth cut-off at 0 .

Perturbative approach
The above shows that it is in fact not necessary to implement a nonperturbative procedure to infer a maximum momentum. It is straightforward to show that in the linear regime ( ≈ ) Eq. 24 becomes, for any turbulence spectrum: This is a generic result of stochastic particle acceleration in perturbation theory (e.g. Ptuskin 1988). A similar result is obtained when considering the stochastic acceleration on MHD waves characterised by a velocity ≈¯as we assumed equipartition between kinetic and magnetic energy in the turbulence (e.g. Thornbury & Drury 2014). In this case, the total momentum diffusion, including both the interaction with magnetic and hydrodynamic waves, is twice that given by Eq. 27. We get: and the acceleration stops when ≈¯, that is, the expected Hillas limit.

Diluted turbulence (loose association)
We now investigate the criterion max ∼¯derived above. The hydrodynamic turbulence is generated in the SB by the massive stars, with an average flow velocity given by Eq. 4. The volume of the SB is computed according to Eq. 1. This procedure implicitely assumes that the stellar energy input is diluted in the whole SB, which is only expected is the stellar cluster is not compact, but rather loose, i.e. no collective wind termination shock forms. The estimate of the maximum momentum in the case of smallscale Bohm diffusion ( ∝ ) is shown in the right panel of Figure 5, as function of the age of the stellar cluster. We took 0 = 1 cm −3 , = 10%, = 10 pc, which are optimistic parameters. Even in this case, only the most massive stellar clusters would produce PeV protons in the few first Myr of their lives. Besides, this result has been obtained assuming a Bohm diffusion regime over the whole SB volume. Harder turbulence spectra such as Kolmogorov or Kraichnan laws are more likely to be generated on this scale (e.g. Gallegos-Garcia et al. 2020). This would lead to a much larger diffusion coefficient and dramatically decrease the maximum energy of the acceleration mechanism. For instance, it was shown by Ferrand & Marcowith (2010) and Vieu et al. (2022a) that in the case of a Kolmogorov turbulence spectrum, the stochastic acceleration mechanism is irrelevant in SBs beyond GeV bands. Therefore, to accelerate protons to PeV energies and above in a SB filled with a realistic turbulence spectrum would require either an extremely small ambient density, or an unrealistically large turbulent energy density, well above the stellar energy input.

Confined turbulence (compact cluster)
In the scenario discussed above the turbulence was assumed to be homogeneous over the whole SB, which is only relevant if the stars are evenly spread over a large region. In the case of a compact cluster, the turbulence is efficiently generated within the MSC and rapidly decreases in the cold wind surrounding the cluster. Thus the velocity and magnetic field fluctuations are confined within the MSC. The acceleration region is then the extent of the cluster and not the whole SB. The turbulence scale is taken as the distance between the stars, = (2/3) 1/3 −1/3 * . Using Eq 4 we can estimate the properties of the turbulence within the cluster as: ≈ 11 0.1 cm −3 1 6 2/3 * 1 3 10 pc where is the density within the boundaries of the cluster. Then the maximum proton energy reads, assuming Bohm diffusion: PeV .
Even in the case of a very massive ( * = 1000) and compact ( = 1 pc) cluster, the maximum energy barely reaches 0.5 PeV in the most optimistic case. The maximum energy achievable via stochastic acceleration in the turbulent SB is found to be much smaller than previous estimates (e.g. Bykov & Toptygin 2001). In the latter, the region of acceleration was assumed to extend over a size = 150 pc and to be filled with strong outflows of mean velocity¯= 3000 km/s. Applying these numbers would imply a maximum proton energy up to 100 PeV. However, the velocity¯which is relevant in this acceleration scenario is not that of the primary SNR and wind termination shocks, but that experienced on average by the particles at each scattering. Most of the time, the particles scatter off secondary waves and not primary winds or SNR shocks. In fact, the primary wind outflows contribute to the average outflow velocity in a loose association of * massive stars as = * ( / ) 3 , where is the average wind velocity and is the average extent of a single star WTS. In other words, the scattering cross section between the particles and the strong primary shocks is very small. On the other hand, to generate homogeneous turbulence with characteristic velocity of the order of 1000 km/s via the decay of the primary shocks would require an unrealistic stellar power or an extremely low density, as discussed above.
Overall, PeV protons might only be reaccelerated via this mechanism in the Galactic centre region:¯∼ 30 km/s and ∼ 100 µG over ∼ 100 pc (Morris & Serabyn 1996) may indeed provide max ∼ 1 PeV assuming a Bohm diffusion regime. Such high turbulence level in the Galactic centre is in particular powered by massive star outflows from the Arches, Quintuplet and Nuclear clusters, whose collective feedback might excite strong turbulence and magnetic fields to stochastically accelerate PeV protons.

CONCLUSIONS
Young MSCs and SBs have been proposed as promising candidates sources for protons at 1 PeV and above (Montmerle 1979;Aharonian et al. 2019;Cao et al. 2021a). In this work, we have shown that this is in fact not straightforward to accomplish. We considered a number of scenarios with typical environmental parameters, consistently ensuring that the global energy balance is conserved, in particular between the mechanical energy input of the star cluster and the turbulence excited around it. Using physically motivated arguments we have justified that the Hillas limit ¯can be applied, provided the characteristic velocity of the system¯is accurately identified, as well as the magnetic field.
Even using optimistic values for the key parameters, the SB forward shock surrounding the whole cavity and the MHD turbulence diluted inside the cavity are generically too weak to accelerate PeV protons. The same conclusion holds when accounting for the contribution of a distribution of strong primary stellar winds, for in this case the scattering cross-section between the particles and the winds is low, even within very massive stellar clusters.
Although primary SNR shocks expanding in a low-density cavity blown by an extended cluster may be able to accelerate protons up to a few PeV if the cluster is young or if CR-driven instabilities operate at their limit, more promising scenarios are found in the case of a compact cluster which powers a collective wind. In this case, the massive stars evolve in a peculiar medium. Indeed, the turbulence is expected to be confined within the cluster as the turbulent energy density released in the collective wind around the cluster is negligible. This implies the generation of large disordered magnetic fields and strong outflows in the cluster core. If the cluster is powerful and compact enough, a large WTS will form during about 10 Myr and accelerate particles. We found that protons up to several PeV are to be expected. Furthermore, SNR shocks may be launched from the edge of the cluster within the free wind region. These will also efficiently accelerate CRs before over-taking the WTS. In this scenario, the most massive compact stellar clusters in our Galaxy such as Westerlund 1 provide channels to efficiently accelerate protons up to about 10 PeV, without the need for additional self-generated magnetic field amplification. Early computations (e.g. Voelk & Biermann 1988) already suggested that SNR shocks expanded in the massive star winds of their progenitors could accelerate PeV protons. The mechanism that we discussed is closely related, though it alleviates the issue of the acceleration efficiency at a perpendicular shock.
Remarkably, our estimate for the maximum energy achieved in this last scenario weakly depends on the mass of the cluster, as long as it is powerful enough to form a collective wind. As most massive stars are expected to be born within compact MSCs, this class of CR sources is expected to be rather common in the Milky Way. Our knowledge of the open cluster population in the Galaxy is still very sparse (Cantat-Gaudin 2022). Considerable progress is currently made in the course of the ongoing analysis of the Gaia Data Release 2 and 3, with thousands of open clusters identified and characterised, among which about 100 are younger than 10 Myr (Cantat-Gaudin et al. 2020;Tarricq et al. 2021), and therefore expected to power a strongly magnetised collective wind into which SNR shocks will expand. With an average rate of one powerful clustered supernova per 100 kyr per young MSC, we should expect at least one of such event per millennium in the Galaxy, which may be more than enough to account for the CR population in PeV bands (Cristofari et al. 2020).
Interestingly, the acceleration of PeV particles around primary shocks embedded in superbubbles always proceeds in a low density environment. Subsequent hadronic interactions are only expected to take place in the dense shell delimiting the superbubble cavity. The incident flux of CR reaching the shell after diffusion in the cavity may generically be too low to produce gamma-ray signatures detectable by the current 100 TeV gamma-ray sensitive instruments such as LHAASO.

DATA AVAILABILITY
No new data were generated or analysed in support of this research.