Cosmic ray acceleration by multiple shocks in the jets of Active Galactic Nuclei

. Active galactic nuclei are one of the most promising sources for accelerating particles up to the highest energies. In this contribution, we present a scenario in which cosmic rays are accelerated in multiple shocks created by the interaction of relativistic AGN jets with the winds of embedded massive stars. We solve the Fokker-Planck equation considering escape and radiative losses as well as the collective e ﬀ ect of the shocks and the reacceleration of the particles. Finally, we calculate the maximum energies that the particles can achieve and discuss the possibility of producing ultra-high energy cosmic rays in this astrophysical situation.


Introduction
Starburst galaxies and active galactic nuclei (AGNs) have been in the center of the discussion about the origin of ultra-high-energy cosmic rays (UHECRs) since the last years. Arguments in favour and against each of these astrophysical objects have been exposed by several authors based on theoretical and observational studies [1][2][3].
When proposing cosmic-ray sources, it is not only necessary to explain the high energies, but also the mass composition inferred from observations, which agrees with intermediate mass nuclei dominating at the higher energies [4,5]. In the case of starburst galaxies, there is a large amount of massive stars and supernova explosions injecting heavy nuclei in the galactic disk, however, their strong photon fields enhance the photo-disintegration of these particles, and therefore, the superwind region in the halo was proposed as the most suitable place in these galaxies to accelerate particles up to the highest energies [1,3]. Nevertheless, from the theoretical point of view to achieve energies of 10 19 − 10 20 eV has been demonstrated to be very unlikely there [see e.g., [6][7][8][9]. Besides that, starbursts do not satisfied the minimum energetic requirement to be sources of UHECRs [10,11]. Jets from powerful radiogalaxies, on the other hand, fulfil the energetic requirements, but justifying the presence of atomic nuclei inside them is not trivial [see e.g., [12][13][14][15][16][17][18]. Moreover, the mechanism and efficiency of relativistic jets to accelerate particles are continuously a matter of debate.
The situation of a jet colliding with stellar winds was proposed a long time ago as a way to mass-loading jets with baryons [19]. Star-jet interactions were also suggested as a possible scenario to explain the bright Xray knots observed in Centaurus A jet and as sources of * e-mail: analaura.muller@eli-beams.eu * * e-mail: anabella.araudo@eli-beams.eu gamma-rays in AGNs [20][21][22][23][24]. Additionally, the wind of massive stars is very rich in intermediate-mass nuclei. In particular, the presence of oxygen, nitrogen, and carbon makes the atmosphere of Wolf-Rayet (WR) stars to have a composition similar to the one observed in cosmic rays [25,26]. Motivated by these ideas, we propose in this work a scenario where cosmic rays are accelerated by multiple shocks. Those shocks are produced by the interaction of the relativistic jet with massive stars, as it is sketched in Fig. 1. In Section 2 we describe our scenario. Section 3 contains the description of the procedure that we follow to calculate the particle distribution of the reaccelerated particles. In Section 4 we present our results and finally, Section 5 is dedicated to our conclusions.

Model
When the jet collides with the wind of a massive star a pair of shock waves is created (see Fig. 1). A non-relativistic shock propagates in the wind with velocity ∼ v wind and a relativistic shock propagates backwards into the jet with speed ∼ v jet . A fraction of the kinetic energy of the jet is dissipated in each interaction and, therefore, the jet slowly decelerates. We assume that its power evolves as L jet (z) = L 0 jet e −τ , where L 0 jet is the initial jet kinetic luminosity and where n ⋆ (z) is the density of stars, and σ sp and σ j are the cross sections of the bow shock and the jet, respectively [22]. The lower limit of the integral is given by the minimum distance at which the first star can be located, we adopt z = 1 pc for this value. The size of the bow shock around the star is determined by the position of the stagna- tion point R sp . By equating the jet and wind kinetic pressures we obtain [20] R sp (z) = whereṀ ⋆ is the mass-loss rate of the star, R jet (z) = z tan θ is the jet radius, and θ is the jet opening angle. Taking into account the effect of the deceleration, R sp (z) can be written in first approximation as with R 0 sp the the location of the stagnation point when the jet does not decelerate and N ⋆ (z) = ∫ z 1 pc n ⋆ (z ′ ) π R jet (z ′ ) 2 dz ′ is the cumulative number of stars into the jet. In Figure 2 we plot N ⋆ using the stellar distribution from [22] n ⋆ pc −3 ∼ 4.6 × 10 −4 (497) y η 0.89 with M BH = 10 8 M ⊙ the mass of the central black hole, η accr = 0.1 the AGN accretion rate, and y = 1 and 2. The opening angle of the jet is assumed to be θ = 5 • . We can see in Fig 2 that for y = 1 one star is expected in the first 50 pc of the jet, ten stars at 100 pc, and fifty stars at 300 pc. The stellar distributions in Eq. (4) correspond to the total distribution of massive stars. In the present contribution we are interested only in WR stars given that they have more powerful winds and a chemical composition rich in intermediate-mass nuclei, although the number of WR stars is expected to be lower than the total number of massive stars. We assume that the number of WR stars is n WR = 0.1n ⋆ . We adopt standard parameters for a WR star:Ṁ ⋆ = 10 −4 M ⊙ yr −1 , luminosity L ⋆ = 10 38 erg s −1 , temperature T ⋆ = 3 × 10 4 K, and v wind = 3000 km s −1 . For the jet we adopt L 0 jet = 10 42 erg s −1 , Lorentz factor Γ 0 = 5, and θ = 5 • . For calculating the deceleration of the jet we consider all the other stars to be OB stars with typical parameters v wind = 2000 km s −1 andṀ ⋆ = 10 −6 M ⊙ yr −1 . In Fig. 3 we plot R sp (z) (and R jet for comparison). Note that R sp increases significantly when the number of WR stars into the jet is much larger than 10 (z > 180 pc for y = 1 and z > 63 pc for y = 2). Additionally, the deceleration of the jet makes also the Lorentz factor to be lower than Γ 0 . Its change can be found by using the definition , with ρ 0 jet the jet density, and solving the transcendental equation In Fig. 4 we plot Γ for the two stellar distributions considered in the present study. We can see that Γ ∼ Γ 0 until z ∼ 60 pc for y = 2 and z ∼ 200 pc for y = 1, where N WR = 10 is large enough to affect the dynamics of the jet and also R sp (z) (See Fig. 3). The magnetic field in the jet is assumed to be B jet (z) = , with z 0 = 5 × 10 −5 (M BH /10 7 M ⊙ ) pc the position where the jet is launched [22,27,28]. The magnetic field of the stellar wind is described as in [29]. Using these expressions, log(z/pc) where Z is the atomic number [30]. Note that R sp = 5 pc and B jet = 125 µG at z = 338 pc for y = 2. In Fig. 5 we show that E Hillas increases as a result of the faster increment of R sp (z) with respect to the decrease of B jet and Γ. At z ∼ 1 kpc, E Hillas = 10 18 eV and 2.6 × 10 19 eV, for protons and iron nuclei, respectively. Nevertheless, E Hillas is an upper limit. The actual maximum energies expected at each shock will be determined by the radiative and escape losses, as we discuss in the next section.

Multiple shock acceleration
Particles can be accelerated in the bow shocks created by the interaction of the jet with the winds of WR stars. Acceleration at single bow shocks has been considered before in order to study the gamma-ray emission from misaligned AGN jets [20][21][22]. Because a large number of stars can be simultaneously into the jet, particles accelerated in a shock can be reaccelerated at another shock at larger z. We follow the next steps to obtain the particle distributions after multiple-shock interactions 1. As we are interested in accelerating nuclei, first we consider the injection of particles by the shock into the wind, because heavy particles are present there. This shock is non-relativistic and strong, with a velocity v sh ∼ v wind . We assume that particles are accelerated by the Fermi I mechanism whose typical injection is S (p) = S 0 p −α e −p/p max , with α = 4 for strong shocks, and the maximum momentum p max determined by equalling the acceleration and losses timescales. The only relevant radiative process for the particles are proton-proton inelastic collisions, and we consider also diffusion and convection escape losses. The normalisation is with L nt the luminosity injected into non-thermal particles. We adopt L nt = 0.1Ṁ ⋆ v 2 wind /2, i.e., 10% of the kinetic energy of the stellar wind. The particle distribution f (p) is obtained by solving the Fokker-Planck equation as in [31].

2.
The particles distribution f (p) produced at the first shock propagates to the next one. During propagation, we only consider that particles can lose energy because of the radiative processes or escape the system. To obtain the particle distribution reaching the next shock, we solve the Fokker-Planck equation with temporal and energy dependency and assume that the particles move along the jet convected by the outflow material, thus time and position inside the jet relates as t = ∫ z 3. The amount of accelerated particles encountering the next bow shock is f prop ∝ (R sp (z)/R jet ) 2 . These cosmic rays will be reaccelerated at the relativistic shocks in the jet. We estimate the distribution after reacceleration in one shock using the following expression [32][33][34] f reacc (p) = α rel where α rel = 3/(1 − ξ −1 ) = 4.5 when the compression factor is ξ = 3. This expression is formally valid for the case of reacceleration in non-relativistic shocks. We consider it as good approximation in our scenario since Γ 0 = 5 and the shocks are mildly relativistic [35]. We calculate p max at each shock and introduce a cut-off e −p/p max to the distribution.
4. All these steps are repeated until there are not more bow shocks inside the jet.

Preliminary results
As a preliminary calculation, we consider the shocks produced by three WR stars arbitrarily located at z 1 = 100, z 2 = 150, and z 3 = 200 pc. At z 1 we consider only acceleration in the wind bow shock, because we are interested in accelerating heavy nuclei. At z 2 and z 3 we consider reacceleration at the shocks in the jet. For the preliminary results presented in this contribution, we do not consider acceleration of particles in the wind shocks at z 1 nor at z 3 . Fig. 6 shows the timescales related to the parameters mentioned before for a wind shock at z 1 (top panel), and jet shocks at z 2 (middle panel) and z 3 (bottom panel), respectively. Adopting the Bohm diffusion regime and the escape distance 0.5R sp (z), the maximum energy of the particles (E max ) is determined by the particle escape out of the acceleration region. In the case of the wind shock, it can be seen that E max is much smaller than in the case of the shocks in the jet, mostly because B jet > B wind . As mentioned in Section 2, R sp (z) increases faster with z than the decay of the magnetic field of the jet, and therefore E max slowly increases with z, obtaining E max (150 pc) = 5 × 10 16 eV and E max (200 pc) = 6 × 10 16 eV for protons.

Particle distributions
We calculate the particle distributions following the procedure described above. In Fig. 7 we exhibit the proton distribution after each step. The distribution injected in one shock arrives to the next one after propagation almost without being modified. This is a consequence of the low jet matter density, which prevents the cooling through proton-proton collisions.After reacceleration, the distribution becomes slightly harder (α = −3.8). Figure 8 displays the injected distribution at z 1 and the result after the interaction with two successive shocks at  z 2 , and z 3 . The distribution extends from shock to shock up to higher energies because of the higher E max at each shock. After the interaction with two shocks, the spectral index of the particle distribution becomes α = −3.6.

Conclusions
We investigate the acceleration of particles at multiple shocks inside AGN jets. These shocks are formed by the interaction of the jet with the winds of WR stars. The jet is charged with protons and heavier nuclei which are removed from the stellar winds. In particular, WR stars con-  tribute with oxygen, nitrogen, and carbon, alike the mass composition inferred from the observation of UHECRs. Particles are accelerated and reaccelerated in the bow shocks by the Fermi I mechanism. We find that the spectrum of the cosmic rays is just slightly modified by propagation between the shocks because of the low jet matter and photon densities. Accounting for the jet deceleration after each interaction, we find that the Hillas energy at the shocks increases with z allowing to reach higher energies. The energy dissipated in the collisions also reduces significantly the Lorentz factor of the jet after the interaction with tens of stars, causing the shocks into the jet to become mildly relativistic. Reacceleration by multiple shocks, on the other hand, causes the spectrum to become flatter than the typical injection slope through first order Fermi acceleration. Overall, we show the potential of this scenario to produce UHECR in AGN jets. Nevertheless, better estimations require to improve the characterisation of the stellar populations and their spatial distribution in the surroundings of AGN jets.