Exploring the superwind mechanism for generating ultrahigh-energy cosmic rays using large-scale modeling of starbursts

The Pierre Auger Collaboration has provided a compelling indication for a possible correlation between the arrival directions of ultrahigh-energy cosmic rays and nearby starburst galaxies. Herein we show how the latest large-scale modeling of starburst galaxies is compatible with the cosmic rays producing the anisotropy signal being accelerated at the terminal shock of superwinds.


I. INTRODUCTION
Starburst-driven superwinds are complex, multiphase phenomena primarily powered by the momentum and energy injected by massive stars in the form of supernovae, stellar winds, and radiation [1]. According to the book, these superwinds are ubiquitous in galaxies where the starformation rate per unit area exceeds 10 −1 M ⊙ yr −1 kpc −2 . The deposition of mechanical energy by supernovae and stellar winds results in a bubble filled with hot (T ≲ 10 8 K) gas that is unbound by the gravitational potential because its temperature is greater than the local escape temperature. The overpressured bubble expands adiabatically, becomes supersonic at the edge of the starburst region, and eventually blows out of the disk into the halo forming a strong shock front on the contact surface with the cold gas in the halo. It was conjectured that charged particles can be accelerated by bouncing back and forth across this terminal shock up to extremely high energies [2,3]. However, criticisms on the choice of model parameters characterizing the speed of the shock and magnetic field strength were reported in [4][5][6][7].
Two recent results have further promoted a reanalysis, though. On the one hand, it was recently proven that the photodisintegration of 4 He on the cosmic microwave background (CMB) is less severe than earlier thought [8]. This implies that the physical survival probability of a 10 EeV helium nuclei coming from the nearest starbursts is 40% larger. On the other hand, full-blown simulations have been performed to accurately capture the hydrodynamic mixing and dynamical interactions between the hot and cold (T ∼ 10 4 K) phases in the outflow [9], as well as to provide a precise determination of the magnetic field in the halo [10]. Armed with the results of these Monte Carlo simulations, we reexamine here the acceleration mechanism in starburst superwinds and show that cosmic ray acceleration up to the highest observed energies is indeed feasible.

II. OBSERVATIONAL STATUS
The Pierre Auger collaboration reported a 4.5σ-significant correlation between the arrival direction of cosmic rays with energy E > 38 EeV and a model based on a catalog of bright starburst galaxies [11,12]. In the best-fit model, 11 þ5 −4 % of the ultrahigh-energy cosmic-ray (UHECR) flux originates from these objects and undergoes angular diffusion on a scale θ ∼ 15°AE 5°. The model signal is largely dominated by the four nearest starbursts: NGC 253 (which is 2.7 Mpc away), M82 (which is 3.6 Mpc away), and NGC 4945 and M83 (which are both 4 Mpc away and relatively close to one another); see Fig. 1. M82, however, is outside the field of view of Auger and therefore the anisotropy signal is dominated by the other three nearby starbursts. This is evident in the result of the test statistics which is almost independent of propagation effects [11,12]. The angular spread derives from a Fisher-Von Mises distribution, the equivalent of a Gaussian on the sphere, and corresponds to a top-hat scale ψ ∼ 1.59 × θ ¼ ð24 AE 8Þ°. Median deflections of particles in the Galactic magnetic field are estimated to be where Z is the charge of the UHECR in units of the proton charge [13]. Thus, the requirement θ G ≲ ψ implies that UHECRs contributing to the Auger anisotropy signal should have Z ≲ 10 and E=Z ∼ 10 EeV. Measurements of the primary energy, and shower properties like the atmospheric column depth at which the longitudinal development of a cosmic-ray shower reaches maximum, X max , and its fluctuations provide reliable information of the nuclear composition [15]. The first two moments of the X max distribution reported by the Pierre Auger collaboration, which are displayed in Fig. 2, seem to indicate that in the energy range of interest the nuclear composition is consistent with a mixture of CNO, F, and Ne, as required to accommodate the observed anisotropy signal, i.e., Z ≲ 10.
First generation UHECR experiments also pointed to a possible starburst origin for the highest energy events [17]. With current statistics, the Telescope Array (TA) collaboration cannot make a statistically significant corroboration or refutation of the UHECR⇋starburst connection [18]. However, TA data have revealed a pronounced directional hot spot [19] in arrival directions not far from the starburst galaxy M82 [20][21][22]. Indeed, the latest analysis carried out by the TA collaboration using ten years of data favors M82 FIG. 1. Model flux skymap of UHECR emission from the 23 brightest nearby starbursts, with a radio flux larger than 0.3 Jy (selected out the 63 objects within 250 Mpc search for γ-ray emission by the Fermi-LAT collaboration [14]) [11]. The continuum emission at 1.4 GHz has been used as a proxy to weight the UHECR emission and an attenuation factor accounting for energy losses incurred by UHECRs during propagation has been included. The map has been smoothed using a Fisher-Von Mises distribution with concentration parameter corresponding to a search radius of 15.0° [12]. The color scale indicates the probability density of the source skymap as a function of position on the sky. The hot spot near the Galactic South Pole originates in NGC 253, located at Galactic coordinates (l ¼ 97. as the origin of the TA hot spot [23]. Constraints based on the isotropic gamma-ray background at ≲TeV measured by the Fermi Large Area Telescope [24] seem to support the association of UHECRs and nearby starburst galaxies [25].
Whereas evidence is mounting, a definite answer as to whether a fraction of UHECRs indeed correlates with starburst galaxies will be given during this decade [26].

III. MODEL UPTAKE
By now, it is well established that charged particles can gain energy by diffusing in converging flows [27]. Diffusive shock acceleration has been widely accepted as the acceleration mechanism for UHECRs [28][29][30][31][32].
Consider a steady state collisionless shock propagating with velocity u s > 0 along a plasma at rest with a magnetic field B in the direction of the shock normal [33]. In the shock rest frame, the flow can be characterized by speeds u s and u s =ζ on different sides of the shock, where ζ is the density compression ratio at the shock. For a nonrelativistic shock with a high Mach number, ζ ¼ 4. Next, assume that particles of charge Ze and speed v ∼ c start off on the upstream (u s ) side of the shock, and diffuse around through collisions with plasma magnetic turbulence until they eventually cross the shock and get off downstream with lower velocity u s =ζ.
The diffusion coefficient D depends on the particle's energy E as well as on the strength and structure of the inhomogeneous magnetic field upstream. In the absence of a fully predictive theory of wave generation and particle diffusion, it is generally assumed that the diffusion coefficient is close to the so-called Bohm limit at all energies, defined by D ¼ r L v=3, where r L ¼ E=ðZeBÞ is the Larmor radius of the particle. This particular value of D is the lowest possible for isotropic turbulence and can only be achieved if the magnetic field is disordered on the scale of a particle's Larmor radius (i.e., if λ=r L ∼ 1, where λ is the mean-free path of the charged particles). Now, in the frame in which the upstream plasma is at rest the diffusion process isotropizes the angular distribution of particles. After a period of time Δt ¼ 4D=ðcu s Þ in which the particles have diffused in a magnetohydrodynamic plasma of speed u s , the particles experience the magnetic turbulence produced by the downstream plasma. Since the downstream plasma also gravitates to isotropize the particles elastically, the population of particles effectively sees a plasma moving towards it with speed u s ðζ − 1Þ=ζ upon arrival downstream (see e.g., [15] for details). This process of quasi-isotropization leads to a net increase in the average particle speed in the rest frame of the shock interface, yielding a mean fractional increase in energy of order ΔE=E ∼ u s =c.
Using the average time that particles spend upstream between shock crossing we find that the particles accelerate at a rate dE=dt ∼ ΔE=Δt ∼ Eu 2 s =ð4DÞ. The characteristic timescale for acceleration to energy E is then E=ðdE=dtÞ ¼ 4D=u 2 s . However, in our back-of-the-envelope estimate we have overlooked the time spent by the particles downstream. The downstream diffusion coefficient is expected to be much smaller than the upstream coefficient, because the downstream magnetic field is larger due to compression in the shock and probably highly turbulent. Therefore, inclusion of the downstream dwell time probably does not increase the acceleration timescale by more than a factor of 2 and following [34] hereafter we proceed by assuming that E=ðdE=dtÞ ¼ 8D=u 2 s . In the absence of energy losses it is straightforward to see that in this approximate model the maximum energy of the particles is whereB is the average magnetic field strength and τ the lifetime of the starburst. The maximum energy of the particles is constrained by the Hillas criterium, which states the Larmor radius of the particle should be no larger than the linear size of the accelerator, where R acc is the length scale of acceleration [35]. Largescale optical emission-line and x-ray nebulae oriented perpendicular to the stellar disks have been observed in several nearby starburst galaxies, suggesting that the superwind activity can be traced out to R acc ∼ 10 kpc. Indeed, the Hα and narrowband continuum optical images of M82 (taken with the KPNO 0.9 m telescope) are quite spectacular, with discernible Hα line emission extending to about 11 kpc above the plane of the galaxy; diffuse soft x-ray emission is seen over the same region by the ROSAT PSPC and HRI [36]. The Hα line emission has a very complex morphology that is highly suggestive of outgoing gas, with loops, tendrils, and filaments of line emission pointing out of the plane of the galaxy. In addition, XMMNewton EPIC observations of M82 show that the soft x-ray emission associated with the superwind extends out to a height of at least 14 kpc in the North and 7.5 kpc in the South [37]. In the case of NGC 253, data from ROSAT PSPC and HRI indicate that the superwind extends out to 10 kpc [38]. Near-infrared broadband and emission-line imaging have revealed the nucleus of NGC 4945 to be the site of a sizable starburst, the presence of which is illustrated by the conically shaped starburst superwind-blown cavity traced at many nearinfrared wavelengths [39,40]. However, optical filaments have only been traced out to about 2 kpc from the plane [41]. This perhaps could be justified by the fact that NGC 4945 lies at low Galactic latitude and consequently the foreground HI column is high. Indeed, as pointed out in [42], the current lack of evidence for the 10 kpc-scale emission that might be expected from a starburst-driven superwind is possibly due to the resulting high foreground optical depth for Hα and soft x-ray radiation. The x-ray data from Chandra is further biased against a detection of diffuse emission due to the significantly higher than normal background experienced during this observation [42]. Finally, M83 has a very complex morphology as inferred from Hα and soft x-ray emission, with an outflow of hot gas extending roughly 10 kpc into the halo above the disk of the galaxy [43,44]. Altogether, Oð10 kpcÞ scale superwinds have become the de facto standard while modeling nearby starburst galaxies [45][46][47][48].
We now turn to justify the fiducial parameters adopted in (2). The modeling of magnetic field is challenging because it is not directly visible to us and we have only a very limited idea about the structure of the fields that may exist in the acceleration region. Because no measurement technique is capable of capturing the entire complexity of the magnetic field, observations only reveal hints onB. Recently, however, a comprehensive study to understand the magnetic field structure in the core and halo of the starburst M82 has been carried out [10]. This study relies on the cosmic ray propagation code GALPROP that selfconsistently models the cosmic ray distribution and its diffuse emission. 1 The analysis focuses on two-dimensional axisymmetric models of the starburst core and superwind.
The starburst core is modeled as an oblate ellipsoid, with cylindrical-radius R core ¼ 0.2 kpc and thickness of z core ¼ 0.05 kpc. Throughout the core the magnetic field strength B and the gas number density n are constants B 0 and n 0 . Outside the core the field BðrÞ is defined as maxfB exp ; B pow g where and and where r ¼ ðR 2 þ z 2 Þ 1=2 is the spherical radius, with along the minor axis (i.e., R ¼ 0) the magnetic field strength becomes To a first approximation we can estimate the average magnetic field asB with z max ¼ 10 kpc. The normalization constants were determined from a simultaneous fit to the gamma-ray data from VERITAS [49] and Fermi [50], as well as the integrated radio data reported in [51][52][53][54]. Using these datasets, the cosmic ray propagation models were constrained in the parameter plane ðB 0 ; n 0 Þ [10]. There are three possible combinations (A, B, and B 0 ) of the two parameters, which provide excellent fits to the combined data and illustrate unique regions of the full parameter space. The numerical values for the parameters of models A, B, and B 0 are given in Table I; u h denotes the maximum assumed superwind velocity. We note in passing that the values of B 0 given in Table I are consistent with previous estimates and they are also similar to the magnetic field strength observed in the nearby starburst galaxy NGC 253 [55][56][57][58][59]. The power-law index β was inferred from the propagation of the cosmic rays out to several kpc, which must create the large radio halo seen at 22 and 92 cm [52]. As can be seen in Fig. 3 the condition B exp < B pow holds for the set of fiducial parameters of model B 0 . Substituting for B 0 and β into (7) we obtainB ¼ 180 μG. To visualize the behavior of B pow out of the minor axis, in Fig. 4 we show a comparison of the ratio B pow ðR; zÞ=B pow ðR ¼ 0; zÞ for various values of R. Since B ≥ B pow , we conclude that (7) provides a reasonable estimate ofB for model B 0 in the acceleration region. For models A and B, the values ofB given in Table I  the real value would not exceed this estimation by more than a factor of 2. Now, the magnetic fieldB carries with it an energy densityB 2 =ð8πÞ, and the flow carries with it an energy flux > u sB 2 =ð8πÞ. This sets a lower limit on the rate at which the energy is carried by the outflowing plasma, and which must be provided by the source. The flux carried by the outgoing plasma is a model-dependent parameter, which can be characterized within an order of magnitude. More concretely, 0.035 ≲ L B =L IR ≲ 0.35, where L IR ∼ 10 43.9 erg=s is the infrared luminosity [45]. The lower limit of L B corresponds to the estimate in [45] considering a supernova rate of 0.07 yr −1 , whereas the upper limit concurs with the estimate in [60], and could be obtained considering a supernova rate of 0.3 yr −1 [61] while pushing other model parameters to the most optimistic values. The relation (8) yields a magnetic field strength in the range 15 ≲B=μG ≲ 150, which reinforces theB interval given in Table I. For NGC 4945, the mechanical luminosity of the superwind is 10 43.8 erg=s and the supernova rate >0.3 yr [62]. The IR luminosity of M83 is measured to be 10 43.2 erg=s [63]. This suggests that both NGC 4945 and M83 can power a Oð100 μGÞ magnetic field, too. Very recently, a series of extremely high resolution global disk simulations of galaxy outflows have been carried out using the Cholla Galactic OutfLow Simulations (CGOLS) program [9]. A physically motivated prescription for clustered supernova served as a rough proxy for the full-blown simulation of a multiphase outflow from a disk galaxy. CGOLS is equipped with a two-phase analytic model capable of fitting the properties of the superwind as a function of radius. To accurately capture the hydrodynamic mixing and dynamical interactions between the hot and cold phases in the outflow the simulations were performed with high resolution (<5 pc) across a relatively large domain (20 kpc). CGOLS studies favor a volume average median velocity for the hot phase of the superwind of u h ∼ 1500 km=s. CGOLS studies also show that mixing between hot and cold gas in the wind is an effective way of transferring momentum from one phase to another and occurs at all radii. This process accelerates cold gas to u c ∼ 800 km=s on kpc scales.
The terminal velocities of the hot (h) and cold (c) phases of the superwind can be parametrized in terms of the thermalization efficiency ϵ (the fraction of energy of the central supernovae and stellar winds that goes into the outflow) and the mass loading factor α (a parameter commonly used to indicate the relative amount of mass loading; if α ¼ 1 the superwind is not mass loaded), where for this order of magnitude calculation, we have assumed that in total a 100 M ⊙ star injects Oð10 51 ergÞ into its surroundings during the wind phase [1]. Observational constraints derived from hard x-ray measurements of M82 pin down medium to high thermalization efficiencies (0.3 ≤ ϵ ≤ 1.0) and require the volume-filling superwind hot fluid that flows out of the starburst region to be only mildly centrally mass loaded (0.2 ≤ α ≤ 0.6) [48]. The allowed phase space then implies a terminal velocity of the M82 superwind in the range 1400 ≤ u h =ðkm=sÞ ≤ 2200, which is consistent with CGOLS simulations. The cold phase structure of the superwind must be examined using a tracer that can be directly related to the cold gas. α was estimated using observations of 12 CO ΔJ ¼ 1 − 0 transition lines in NGC 253, obtained with the Atacama Large Millimeter Array (ALMA) [64]. ALMA data favor a mass loading α ∼ 3 (with a lower limit of α ∼ 1). This yields 600 ≲ u c =ðkm=sÞ ≲ 1000, a range which is also consistent with CGOLS studies. The estimated superwind velocity of NGC 4945 is somewhat above 1000 km=s [65]. Altogether, these numbers justify the fiducial value, u s ∼ 1000 km=s, which agrees with the maximum superwind velocity assumed in the simulations with GALPROP; see Table I [ 10]. The age of the starbursts is also subject to large uncertainties. Observations seem to indicate that about 500 Myr ago M82 experienced a tidal encounter with its large spiral neighbor galaxy, M81, resulting in a large amount of gas being channeled into the core of the galaxy over the past 200 Myr, which produced a concentrated starburst and an associated pronounced peak in the cluster age distribution [66]. This starburst has continued for up to  about 50 Myr at a rate of 10 M ⊙ yr −1 [67]. A survey of bright evolved stars in the NGC 253 disk suggests that about 200 Myr ago this galaxy also interacted with a now defunct companion [68]. Some models predict that the age of the starburst in NGC 253 is similar to that of M82 [69]. Age estimates for the nuclear starburst of NGC 4945 are sparse and intrinsically uncertain. On the one hand, the excitation of the starburst core is very low, and this sets a lower limit on the starburst age of 5 Myr [70]. On the other hand, the mechanical luminosity associated with its high supernova rate is consistent with a starburst age of up to 30 Myr or even 50 Myr [40]. The age of the nuclear starburst in M83 was estimated to be between 10 and 30 Myr by comparing narrowband and broadband photometry of massive star clusters with theoretical population synthesis models [71]. However, the current data cannot exclude a longer duration of activity. Altogether this motivated our choice for the starburst age of τ ∼ 50 Myr. We now turn to discuss possible energy loss at the source. The accelerating nucleus will gain energy until it eventually suffers a photodisintegration. The photodisintegration rate depends on the energy density of the ambient radiation field. This is governed by the spatial distribution of photons, including both those from the CMB and stellar radiation fields. For compact regions near the galaxy core, starbursts exhibit an energy density in their stellar radiation fields which may exceed (or be comparable to) that of the CMB, but at the superwind scale starlight is expected to have a negligible energy density compared to that of the CMB [72,73].
We duplicate the statistical analysis given in [74] to determine the maximum energy at the source. The results for various nuclear species are shown in Fig. 5. We conclude that for the fiducial value of the acceleration time scale τ ∼ 50 Myr and magnetic field strengths in the range 15 ≲B=μG ≲ 150, the photodisintegration energy loss during the acceleration process can be safely neglected for nuclei with Z ≲ 10, and therefore (2) gives the maximum attainable energy.
Because the sources dominating the anisotropy signal are in our cosmic backyard, the energy losses during propagation are almost negligible. This is visible in the result of the test statistics reported by the Pierre Auger collaboration [11,12]. Note that the mean-free path of a 70 EeV nitrogen nucleus on the CMB is roughly 3 Mpc [15] and therefore CNO, F, and Ne can survive the trip from nearby sources in the energy range of interest. On the other hand, even if He nuclei could reach energies beyond 38 EeV in the acceleration process, because of their photodisintegration on the CMB [8] they would only contribute to the anisotropy signal near the minimum energy.

IV. CONCLUSIONS
A number of conclusions are in order connecting back with (2). If we assume the fiducial values for the velocity and age in (2), two of the three models examined (i.e., A and B, always in the notation of [10]) would produce too low averaged values of energies for the accelerated particles; E max ∼ Z EeV. Unless of course there is observational evidence for Z being sufficiently high, if models A or B are realized, starburst superwinds as a possible origin of the correlation commented in Sec. II will be challenging. Model B 0 , on the other hand, would naturally lead to particles whose averaged energy could be in excess of 38 EeV and beyond, E max ≳ Z × 10 EeV. In this case, relatively light nuclei (e.g., in the CNO, F, Ne region) would reach very high energies easily. This is consistent with observation; see Fig. 2. Even He nuclei may reach the needed UHECR regime. This would be helped by a less-dramatic photodisintegration along the way: new cross section parametrizations [8] have shown that, e.g., at E ∼ 10 1.8 EeV and a propagation distance of 3.5 Mpc, the 4 He intensity would be 35% larger than the output of the CRPropa 3 program and 42% larger than the output of the SimProp v2r4 program. Thus, 4 He could also contribute to the signal near threshold (E ∼ 38 EeV). The latter would of course be easier for larger velocities: We have assumed a FIG. 5. Energy loss during the acceleration process. The solid lines indicate the time independent cutoff energy for each particle species. This would be a brainwork situation in which the source produces an infinite number of identical shocks, with fresh injection at each shock and decompression between the shocks, such that the acceleration timescale ≫50 Myr [74][75][76][77]. It is very unlikely that particles above the cutoff energy would have not suffered at least one photodisintegration. The bands encompass 68% of the nuclei in the statistical sample. This means that only 16% of the nuclei would reach energies above the upper border of the shaded bands and 16% would only reach energies below the lower border. The dashed lines indicate the maximum acceleration energy as given by (2) in a time period of 50 Myr. This means that the maximum energy attained by the particles is the minimum between the solid and dashed lines: E max ¼ min (solid, dashed) lines, for each particle species. Note that for a timescale of 50 Myr, the particles cannot reach extreme energies unless the magnetic field is Oð100 μGÞ. For τ ∼ 50 Myr and 15 ≲B=μG ≲ 150, the photodisintegration energy loss for particle species with Z < 10 can be neglected, and the maximum energy is given by (2). fiducial value of 1000 km=s in (2), but in reality, u c < u s < u h and thus it is plausible to expect this could be greater than assumed. All in all, we have shown that UHECR acceleration in superwinds generated by nearby starbursts (M82, NGC 4945, NGC 253, and M83) is feasible.
Interestingly, we note that if with future data we indeed confirm the correlation with starbursts, in the context of superwinds, UHECR observations would directly feedback on models: The correlation would favor larger magnetic fields in the outflow. Observational guidance in assessing this conundrum is, at this point, essential. Tidal disruption events caused by black holes [22], low-luminosity gammaray bursts [78], and black hole unipolar induction [6] have been identified as potential UHECR accelerators inside starburst galaxies. However, all of these possible UHECR origins fail to explain which is the inherently unique feature(s) of starburst galaxies to account for their correlation [79]. A true smoking gun for all these proposals would be a correlation with the distribution of nearby matter or the more extreme active galactic nucleus as opposed to starburst galaxies. This is not what data seem to suggest. The data yielding the anisotropy signal favor a production mechanism of UHECRs above 38 EeV which is exclusive to starbursts. Hence, should future observations confirm the UHECR⇋starburst connection, Fermi-shock acceleration in starburst superwinds appears to be a feasible alternative for this distinguishing feature.