Diffusive shock acceleration at EeV and associated multimessenger flux from ultra-fast outflows driven by active galactic nuclei

ABSTRACT

Active galactic nuclei (AGN) can launch and sustain powerful winds featuring mildly relativistic velocity and wide opening angle. Such winds, known as ultra-fast outflows (UFOs), can develop a bubble structure characterized by a forward shock expanding in the host galaxy and a wind termination shock separating the fast, cool wind from the hot shocked wind. In this work, we explore whether diffusive shock acceleration can take place efficiently at the wind termination shock of UFOs. We calculate the spectrum of accelerated particles and find that protons can be energized up to the EeV range promoting UFOs to promising candidates for accelerating ultra-high energy cosmic rays (UHECRs). We also compute the associated gamma-ray and neutrino fluxes and compare them with available data in the literature. We observe that high-energy (HE) neutrinos are efficiently produced up to hundreds of PeV while the associated gamma rays could be efficiently absorbed beyond a few tens of GeV by the optical-ultraviolet AGN photon field. By assuming a typical source density of non-jetted AGN, we expect that UFOs could play a dominant role as diffuse sources of UHECRs and HE neutrinos. We finally apply our model to the recently observed NGC1068 and we find out that under specific parametric conditions an obscured UFO could provide a sizeable contribution to the observed gamma-ray flux while only contributing up to ∼10 per cent to the associated neutrino flux.

acceleration of particles, neutrinos, cosmic rays, galaxies: active, gamma-rays: galaxies

1 INTRODUCTION

Fast wide angle winds are one of the most intriguing feedback mechanisms of active galactic nuclei (AGN; Silk & Rees 1998). Their impact on the host galaxies has long been considered to affect the dynamical evolution of the interstellar medium (ISM) and act as a regulator of star formation (Crenshaw, Kraemer & George 2003). The discovery of blue-shifted Fe K absorption lines in X-ray spectra of AGN (see e.g. Chartas et al. 2002) brought compelling evidence of mildly relativistic velocities typically ranging from |$0.1 \, c$| to |$0.3 \, c$| (where c is the speed of light) in such winds, which thereafter have been often referred to as ultra-fast outflows (UFOs). Mildly relativistic flows were already known to be present in the vicinity of AGN from broad absorption lines observed in their spectra (Weymann et al. 1991). The discovery of UFOs allowed us to understand that high kinetic luminosity (⁠|$\dot{E} \simeq 10^{41}\!-\!10^{45}\,{\rm erg\,s}^{-1}$| Tombesi et al. 2015; Fiore et al. 2017) was also possible in fast AGN-driven winds with wide opening angles. Recently UFOs have been systematically detected in both radio-quiet and radio-loud AGN (see e.g. Markowitz, Reeves & Braito 2006; Braito et al. 2007; Cappi et al. 2009; Tombesi et al. 2010a, b, 2015). The number of observations keeps increasing with time, despite of the observational challenges, also thanks to high resolution grating spectra in the soft X-ray (see e.g. Pounds et al. 2003) and the discovery of ultraviolet (UV) lines in addition to those already known in the X-ray (Kriss et al. 2018; Mehdipour et al. 2022). The search for a launching mechanism of the UFOs has not found a definitive answer yet, although there is a general agreement to ascribe this phenomenon to the accretion activity of the super massive black hole (SMBH) (see Laha et al. 2021, for an updated review).

The dynamics of the wind and the associated feedback on the host galaxy depend on whether the outflow conserves energy or momentum (King & Pounds 2015). In particular, if the wind plasma does not cool efficiently when it shocks with the surrounding medium, the system is energy-conserving, and the momentum flux is boosted while the wind sweeps up the external matter. In the opposite scenario, namely, if the wind radiates most of its thermal energy, it evolves conserving the momentum. In spite of the strong radiation fields that could strongly affect the cooling of electrons, Faucher-Giguère & Quataert (2012) showed that two-temperature plasma effects are likely to slow down radiative losses for protons thereby favouring an energy-conserving dynamics.

A recent Fermi-LAT analysis (Ajello et al. 2021) showed that UFOs are a new class of gamma-ray emitters. In this analysis, the average gamma-ray emission from a sample of 11 nearby (⁠|$z < 0.1$|⁠) radio-quiet AGN with an UFO is derived by adopting a stacking analysis. The average best-fit gamma-ray spectral slope is measured to be 2.1 ± 0.3, and the gamma-ray luminosity is found to scale with the AGN bolometric luminosity.

AGN-driven outflows, similar to stellar winds (see e.g. Weaver et al. 1977; Koo & McKee 1992a, b), are expected to develop a structure characterized by an inner wind termination shock (hereafter wind shock), a contact discontinuity, and an outer forward shock. The forward shock has been proposed as a plausible site for particle acceleration (see e.g. Lamastra et al. 2016; Wang & Loeb 2016a; Lamastra et al. 2019; Ajello et al. 2021), where ideal conditions for efficient production of gamma rays and high-energy (HE) neutrinos are expected (see also McDaniel, Ajello & Karwin 2023, for a recent study on molecular outflows). Wang & Loeb (2017) highlighted also the possibility that, in somewhat extreme conditions, a fast AGN-driven wind could have the energy budget to accelerate CRs up to the ultra-high-energy (UHE) range. The associated cumulative contribution of AGN-driven winds to the diffuse gamma-ray and neutrino flux has been also explored (see e.g. Wang & Loeb 2016a, b; Lamastra et al. 2017; Liu et al. 2018). In addition, the amplitude of the recently observed spectrum of the diffuse neutrino flux (Abbasi et al. 2021), in light of the constraints imposed by the diffuse gamma-ray flux observed by Fermi-LAT (Ackermann et al. 2015), suggests that there could be a class of sources at least partially opaque to gamma rays (see e.g. Murase, Guetta & Ahlers 2016).

Indeed, a search for time-integrated point-like neutrino sources (Aartsen et al. 2020) highlighted an excess in the direction of the Seyfert galaxy NGC1068. The most intriguing aspect of the emission of such a galaxy is the lack of gamma rays in the TeV band (Acciari et al. 2019), where the neutrino flux is observed. The natural implication of a highly opaque cosmic particle accelerator triggered several studies on the multimessenger implications of HE particles populating the innermost region of AGN, such as discs and accretion flows (see e.g. Kimura, Murase & Mészáros 2019; Gutiérrez, Vieyro & Romero 2021), combined emission from successful and failed AGN winds (Inoue et al. 2022), and AGN corona (see e.g. Inoue, Khangulyan & Doi 2020; Murase, Kimura & Mészáros 2020; Kheirandish, Murase & Kimura 2021; Eichmann et al. 2022; Murase 2022). Interestingly, we recently witnessed a growth in the statistical significance of the signal from NGC1068 (Abbasi et al. 2022a).

Even though there is an increasing evidence for HE particles populating the innermost regions of active galaxies, our understanding of the acceleration mechanisms at play in such environments is still incomplete. The mechanisms capable of energizing HE particles in the vicinity of SMBH is indeed a partially unexplored field we aim to assess in this work together with its multimessenger consequences.

Therefore, we develop a model for particle acceleration by exploring the diffusive shock acceleration (DSA) mechanism at the wind shock of UFOs. At this shock, unique conditions for acceleration of protons at energies as high as ∼10¹⁸ eV can be found. In addition, the medium property could make such sources bright in HE neutrinos while being partially opaque to gamma rays beyond 10–10² GeV. In particular, the optical thickness to gamma-rays depends on the specific parametric conditions. We discuss the role of UFOs as UHE cosmic ray (UHECR) sources in light of the spectral behaviour of the particle flux escaping the AGN wind bubble. In particular, we show that at the highest energies spectral features harder than E⁻² could appear in the spectrum of escaping particles due to the interplay of diffusion-advection and energy losses. This result is of great interest in light of phenomenological studies aiming at modelling the spectrum and mass composition observed by the Pierre Auger Observatory, which suggest UHECRs to be characterized by hard spectra (see e.g. Unger, Farrar & Anchordoqui 2015, and references therein). Moreover, if UFOs were common in galaxies, this would produce important implications for their diffuse multimessenger emission. We provide here an order of magnitude estimate of the potential role of UFOs for a diffuse flux of HE neutrinos and CR at the ankle. Finally, we explore whether an UFO could play a role in the neutrino emission observed in NGC1068 discussing possible realizations of such a system, which better agree with observed data.

The manuscript is organized as follows: in Section 2, we describe the model for the wind bubble and the associated particle acceleration and transport formalism; in Section 3, we discuss our results in terms of spectra of accelerated and escaping particles, we perform a parameter space scan focusing on the maximum energy. In Section 4, we discuss the multimessenger implications in terms of HE photons, and neutrinos and in Section 5, we specialize our model by applying it to the case of NGC1068 and discuss possible model improvements and alternative scenarios. We draw our conclusions in Section 6.

2 MODEL FOR PARTICLE ACCELERATION AND MULTIMESSENGER EMISSION IN UFO

The fast wind launched and sustained by AGN expands with large opening angle and mildly relativistic velocity. At such speed, the outflow is supersonic therefore it drives a forward shock expanding in the external medium, while a contact discontinuity separates the shocked wind (SW) material from the shocked ambient medium (SAM). The impact of the wind on the external medium creates a shock inside the wind material, the wind shock, which is trailing behind the contact discontinuity and it is oriented towards the central engine. The outflow, characterized by these three discontinuities, features a bubble structure as sketched in Fig. 1.

Figure 1.

Structure of the wind bubble. The SMBH (BH) responsible for the wind launching is located on the left of the sketch. The blue (red) arrows correspond to the cool (shocked) wind of the upstream (downstream) region. The wind shock (R_sh) separates these two regions. The SAM is located between the contact discontinuity (R_cd) and the forward shock (R_fs), which bounds the system (credit: I. Peretti).

The environment surrounding active SMBHs, where UFOs expand is extremely complex: from the innermost parts of AGN, one can find the accretion disc, the broad-line region (BLR) with highly dense clouds of density up to |$n_{\rm c} \lesssim 10^{10} \, \rm cm^{-3}$| (see e.g. Ricci & Trakhtenbrot 2022), a dusty torus, the narrow-line region with clouds of typical density of about |$10^4 \, \rm cm^{-3}$| and often a larger circum-nuclear disc (see e.g. Urry & Padovani 1995). For the purposes of modelling the dynamics of the UFO, we assume a uniform circum-nuclear medium (CNM) of effective density n₀. A plausible range of values for n₀ can be derived from the column density N_H, based on the work by Ricci et al. (2017): considering |$N_{\rm H} \lesssim 10^{25} \,{\rm cm}^{-2}$| as a typical range for the column density of AGN and assuming |$10^6\!-\!10^9 \,{\rm M}_{\odot }$| as a plausible mass range for SMBHs, we obtain 1.5–63 pc as a typical radius for the sphere of influence of the black hole (⁠|$GM_{\rm SMBH}/\sigma _*^2$|⁠) and an upper limit to the external medium density of |$n_0 \lesssim 5 \cdot 10^4\!-\!10^6\,{\rm cm}^{-3}$|⁠. Here, we adopted the relation between black hole mass and stellar velocity dispersion, |$M_{\rm SMBH}/10^9\,{\rm M}_{\odot}\simeq 0.309 \times (\sigma _{*}/ 200\,{\rm km\,s^{-1}})^{4.38}$|(Kormendy & Ho 2013).

In a uniform medium of density n₀, the wind expands with constant velocity up to the radius at which the swept-up mass roughly balances the whole mass of the outflow. After the swept-up mass becomes dynamically relevant, the outflow starts decelerating. During the deceleration phase, the forward shock R_fs and the wind shock R_sh evolve self-similarly according to different scaling laws: R_fs ∼ t^3/5 and R_sh ∼ t^2/5 (see also Weaver et al. 1977; Koo & McKee 1992a, b, for detailed discussions and Appendix  A for additional details). Since the wind shock decelerates faster than the forward shock, the hot bubble, namely the spherical shell between the wind shock and the contact discontinuity, grows with time while remaining approximately adiabatic. In this context, the wind bubble evolution can be considered as energy-conserving (see e.g. Faucher-Giguère & Quataert 2012). On the other hand, radiative losses in the SAM can be efficient (see e.g. Nims, Quataert & Faucher-Giguère 2015), so that the whole swept-up mass is eventually compressed into a relatively thin layer between the contact discontinuity, R_cd, and R_fs. While the wind bubble slows down its expansion, the relative velocity between the plasma and the wind shock remains high, namely the shock stays strong. Therefore, we refer to the innermost region of free expanding wind and to the SW, respectively, as upstream and downstream. Different from the wind shock, the Mach number of the forward shock strongly depends on the temperature and conditions of the surrounding medium. Therefore, it is not guaranteed that the forward shock is strong for a sufficient amount of time.

We assume a spherically symmetric geometry for the outflow and, since the wind launching region has a negligible size compared to the whole bubble, we also assume a constant upstream velocity u₁. The SW is adiabatic therefore the velocity profile reads: u₂(r) = u₂(R_sh/r)², where u₂ = u₁/4. In agreement with observations of UFOs, we limit our investigation to a plasma velocity lower than ∼0.3c. We thus neglect relativistic effects due to the mildly relativistic motion of the plasma. The gas density in the upstream region scales as |$n_1(r) = \dot{M}/\big[4 \pi r^2 u_1 m_p \big]$|⁠, while in the downstream region it is constant and equal to n₂ = 4n₁(R_sh). The gas density between R_cd and R_fs depends on the amount of matter accumulated during the outflow evolution, |$n_{\rm SAM} = n_0 R_{\rm fs}^3/(R_{\rm fs}^3-R_{\rm cd}^3)$|⁠, where we assume |$R_{\rm cd} \simeq 0.9 \, R_{\rm fs}$| (Sharma et al. 2014). Under this assumption, one obtains n_SAM ≈ 4n₀ as a typical gas density in the SAM layer. We postulate a turbulent nature for the magnetic field, and we estimate its amplitude in the upstream region under the assumption that a fraction |$\epsilon_{\rm B} \lesssim 10\ \hbox{per cent}$| of the ram pressure is converted into magnetic field energy density, namely |$U_{\rm B}(r) = \epsilon_{\rm B} m_p n_1(r) u_1^2$|⁠. At the wind shock, we assume that the magnetic field gets compressed by a factor |$\sqrt{11}$|⁠, typical of strong shocks, and remains constant throughout the whole downstream region. We adopt the quasi-linear theory of diffusion, |$D(r,p)= v(p) r_{\rm L}^{2-\delta }(r,p) l_{\rm c}^{\delta -1}/3$|⁠, where v is the particle velocity, r_L is the Larmor radius, δ is the slope of the turbulence power spectrum, and l_c is the coherence length of the magnetic field that we assume to be comparable in size with the launching radius of the wind. In addition, we account for the small-angle scattering regime of diffusion, namely |$D\propto r_{\rm L}^2$|⁠, which takes place when r_L > l_c (see e.g. Subedi et al. 2017; Dundovic et al. 2020).

The average lifetime of AGN is inferred to be ≲10⁷ yr (Yu & Tremaine 2002). During this time, the AGN is expected to show multiple episodes of activity with duty cycles of ≲10⁵ yr duration (Schawinski et al. 2015). This suggests that, even though their typical age is not known at present, UFOs could have a lifetime ranging from hundred years up to several thousands of years. In this work, we explore UFOs under the assumption that they can be powered for a sufficient amount of time so as to allow them to reach the deceleration phase. Therefore, |$\rm O(10^3 \, \rm yr)$| is a conservative assumption, while we comment in Section 3.1, the impact of assuming a different lifetime up to values comparable with the AGN duty cycle. In this context, the dynamical evolution of the system becomes soon slower than all relevant time-scales involving HE particles. Hence, the process of particle acceleration and transport can be treated as stationary (see Fig. 2 in Section 3, where the typical time-scales for HE particles in a prototype UFO are discussed). We assume a spherically symmetric transport, where particles are injected via DSA at the wind shock, whereas, once they reach the forward shock location, they freely escape the wind bubble. The transport equation reads:

$$\begin{eqnarray} r^2 u \partial _r f = \partial _r[r^2 D \partial _r f] + \frac{p}{3} \partial _p f \, \partial _r[r^2 u] + r^2 [Q-\lambda f]\, , \end{eqnarray}$$

(1)

where u = u(r) is the wind velocity profile, D = D(r, p) is the diffusion coefficient, Q = Q(r, p) is the injection term, and λ = λ(r, p) is the loss rate accounting for pp and pγ interactions (see Appendix  B for additional details). As boundary conditions, consistently with the spherical symmetry, we assume a null net flux at the centre of the system, |$u f-D \partial _r f|_{r=0}=0$|⁠, while, as mentioned earlier, we regard the forward shock as a free escape boundary, f(R_fs, p) = 0. The injection term reads:

$$\begin{eqnarray} Q(r,p) = Q_0(p) \delta [r-R_{\rm sh}] = \frac{\eta _{\rm CR} u_1 n_1}{4 \pi p^2} \delta [p-p_{\rm inj}] \delta [r-R_{\rm sh}], \end{eqnarray}$$

(2)

where |$p_{\rm inj} = 1 \, \rm GeV/c$| is the injection momentum of particles (picked up from the plasma) that enter the DSA process and η_CR is the efficiency factor, normalized such that the CR pressure at the shock is a small fraction (⁠|${\lesssim}10\ \hbox{per cent}$|⁠) of the plasma ram pressure.

Figure 2.

Typical time-scales regulating the transport of HE particles compared with the age of the system (thick black line). The blue dashed line represents the acceleration time-scale, while energy losses via photomeson and Bethe–Heitler pair-production are represented, respectively, by red and magenta dot–dashed lines. Advective and diffusive escape are represented by orange and green dotted lines.

We solve equation (1) following the same procedure developed in Morlino et al. (2021) and Peretti et al. (2022) with the modification that, in the present work, energy losses in the downstream region are also accounted for due to their possible observational relevance. In fact, the relative small distance from an active SMBH makes the environment potentially hostile for HE particles, where energy losses could affect the acceleration and limit the escape.

The details of the calculation are reported in Appendix  B, while the general form of the solution at the wind shock reads:

$$\begin{eqnarray} f_{\rm sh}(p) = C p^{-s} {\rm e}^{-\Gamma _{\rm cut}(p)}, \end{eqnarray}$$

(3)

where C is a constant (see Appendix  B), s is the spectral index (s = 4 in strong shocks), and Γ_cut is a HE cut-off function (Γ_cut = Γ_l + Γ_e as detailed in Appendix  B), which increases with momentum.

As the background photon field, we use the spectral energy distribution (SED) model shown in the top panel of Fig. A1 provided in Appendix  A. Such a SED is characterized by the big blue bump and an X-ray power-law component as described in Marconi et al. (2004). Such a field is assumed to decrease with the second power of distance from the central engine. Consistently with the AGN SED, we also account for the far infrared (FIR) component produced by a dusty torus (Mullaney et al. 2011). While the r⁻² profile is adequate for the photon field produced by the accretion disc and AGN corona (optical and X-rays), we point out that the FIR radiation field produced by the dusty torus is probably characterized by a more complex spatial structure (Błażejowski et al. 2000), possibly with a uniform behaviour up to a radius comparable with the typical size of the torus (≈1 pc). In Section 3.1, we discuss the implications and possible limits of such approximation.

Gamma rays and neutrinos from pp and pγ interactions are computed following Kelner, Aharonian & Bugayov (2006) and Kelner & Aharonian (2008), respectively. The energy losses due to Bethe–Heitler pair-production are taken into account as energy loss mechanism taking place at a rate |$t^{-1}_{\rm BH}$| (see Gao, Asano & Mészáros 2012). Gamma–gamma absorption on the AGN SED, including the torus, is also accounted for by adopting the cross section appropriate for the case of an isotropic photon field (see Aharonian 2004). Finally, the associated flux-density of HE protons escaping the system is computed self-consistently from the solution to the transport equation as |$j_{\rm esc}(p) = - D_2 \partial _r f|_{r=R_{\rm fs}}$|⁠.

The calculations illustrated here have been carried out in the context of the thin shell approximation, in which the SW is perfectly separated from the SAM. This implies that the cold gas acting as a target for pp interactions (see below) is all located close to the R_fs. It is worth pointing out that the spatial distribution of the gas in the SW might be affected by instabilities and mixing that may lead to a more pervasive, though clumpy, structure of the gas, which in turn might affect the spatial and spectral properties of the secondary emission. These effects will be investigated in a future dedicated work.

3 RESULTS

In order to present our model and discuss its physical implications, we assume a set of typical parameters, hereafter referred to as our benchmark scenario. The parameters defining our benchmark scenario, summarized in Table 1, have been chosen according to the following criteria: we assume |$u_1=0.2 \, c$| as the average value for the terminal wind speed of UFOs; |$\dot{M}= 10^{-1} \, \rm {\rm M}_{\odot } \, yr^{-1}$| has been chosen such that the total kinetic power |$\dot{E} = \dot{M} u_1^2/2$| matches about ∼3 per cent of the total AGN bolometric luminosity (Fiore et al. 2017) characterized by an X-ray luminosity |$L_X=10^{44} \, \rm erg \, s^{-1}$|⁠; |$l_{\rm c} = 10^{-2} \, \rm pc$| is compatible with the launching radius of the wind as predicted by accretion disc wind models (see e.g. Murray et al. 1995); the age of the system |$t_{\rm age} = 10^3 \, \rm yr$| has been chosen in order to assure stationary conditions, which are not guaranteed for much younger systems (the resulting shock radii at such an age are |$R_{\rm sh}\approx 0.8 \, \rm pc$| and |$R_{\rm fs}\approx 3 \, \rm pc$|⁠); ϵ_B = 0.05 guarantees a minor dynamical impact of the turbulent magnetic field; δ = 3/2 is motivated by an MHD-like (Kraichnan) turbulence cascade.

Table 1.

Parameters of the benchmark UFO and of the three alternative scenario considered for the multimessenger emission.

Parameter	benchmark	A	B	C
u1/c	0.2	–	–	–
|$\dot{M} [\rm {\rm M}_{\odot } \, yr^{-1}]$|	10−1	–	–	–
ξCR	0.05	0.087	0.1	0.12
ϵB	0.05	–	–	–
|$l_{\rm c} [\rm pc]$|	10−2	–	–	–
δ	3/2	–	–	–
|$L_X [\rm erg \, s^{-1}]$|	1044	–	–	–
|$n_{0} [\rm cm^{-3}]$|	104	2 · 103	5 · 102	2 · 102
|$t_{\rm age} [\rm yr]$|	103	3 · 103	104	2 · 104

In agreement with the upper limits presented in Section 2, we assume |$n_{0} = 10^4 \, \rm cm^{-3}$| as a typical value for the external ambient medium that can also be found in the core of luminous infrared galaxies (see e.g. Downes & Solomon 1998; Faucher-Giguère & Quataert 2012), and we discuss in Section 3.1 the impact of adopting densities up to |$n_0=10^6 \, \rm cm^{-3}$|⁠. Higher densities could be reached if the bulk of the column density were concentrated within a smaller volume consistent with the size of the BLR, |$R_{\rm BLR}\lesssim 0.1 \,{\rm pc}$| (Bentz et al. 2009). However, the investigation of such a scenario goes beyond the limits of applicability of our model because, in such dense and compact environment, the system would be in calorimetric conditions. Nevertheless, we discuss possible implications of ultra-dense environments in Section 5.2.

Even though we study particle acceleration and transport by solving the stationary transport equation, equation (1), a more direct understanding of the physics property of the solution can be obtained by analysing the typical time-scales of the different competing processes. Fig. 2 illustrates the typical time-scales for HE particles as computed at R_sh for the benchmark scenario. Here, the age of the system (τ_age, thick black line) is compared with the following time-scales: acceleration (⁠|$\tau _{\rm acc} \approx s D_1/u_1^2$|⁠, blue dashed line), diffusion (τ_diff = (R_esc − R_sh)²/D₂, green dotted line the diffusion), advection (τ_adv = (R_esc − R_sh)/ <u₂>, orange dotted), the pγ photomeson (τ_pγ, red dot–dashed line), and Bethe–Heitler pair-production (τ_BH, magenta dot–dashed). Inelastic pp collisions were also taken into account, however at the wind shock, the target density is of the order of |$20 \, \rm cm^{-3}$|⁠, so that the associated time-scale exceed |$10^5 \, \rm yr$|⁠, the upper limit of the plot. Therefore, pp interactions are irrelevant for the acceleration. On the other hand, this does not exclude them as relevant loss mechanisms in the SAM, where the density is orders of magnitude higher. From the interplay between the different time-scales, it is possible to observe: 1) τ_acc ≪ τ_age and the minimum between losses and escape is also smaller than the age supporting the stationary assumption; 2) τ_acc as well as τ_diff feature a break between 10²PeV and 1 EeV due to the transition in the diffusion coefficient from the QLT behaviour (r_L < l_c) to the small-angle scattering regime (r_L > l_c); 3) energy losses via pγ photomeson production play a dominant role at the highest energies and are expected to set the maximum energy; 4) τ_pγ and τ_BH increase with the second power of the distance moving outward from R_sh to R_fs therefore the transport in the downstream region is characterized by a close competition between energy losses and escape; 5) the Bethe–Heitler pair production does not play a dominant role since at low energy the transport is advection-dominated while at the highest energies is regulated by the photomeson production.

The spatial transport of particles is regulated by advection at low energy and by diffusion at the highest energies while energy losses, both adiabatic (in the upstream region) and inelastic collisions, can affect the normalization and/or introduce spectral features. The top panel of Fig. 3 shows the spatial distribution for three different CR energies. The upstream region (R/R_sh < 1) is characterized by the competition between diffusion, which tries to homogenize spatially the particles, and advection which prevents low energy particles to diffuse upwind. In particular, one can see that the higher the energy the stronger the impact of diffusion. The red dotted line illustrates the spatial distribution at low energies while the blue and black lines show results for an intermediate energy and near the exponential cut-off, respectively. In the downstream region (R/R_sh > 1) one can see that, different from the upstream one, advection-dominated transport leads to a spatially homogenized solution whereas diffusion-dominated transport leads to a number suppression while approaching the free escape boundary. This behaviour is a natural result of the spherical geometry of the system.

Figure 3.

Top panel: spatial distribution of the CR phase space density. Low-energy particles behave in the system as illustrated by the red dotted line, high-energy particle behaviour is represented by the blue dot–dashed curve, while the black curve shows the behaviour of particles at the maximum energy. Bottom panel: spectrum of particles at the shock (thick green line) compared to the spectral shape of the escaping flux (dashed magenta line). The dotted curves represent the particle spectra in the downstream region. From red to blue, the dotted lines are computed at r/R_fs = 0.29, 0.33, 0.44, 0. 57, 0.75, 0.9.

The bottom panel of Fig. 3 illustrates the spectrum of accelerated particles at the shock (thick green line), at different radii in the downstream region (dotted curves where the red one is the closest to R_sh while the blue one approaches R_fs) and the associated spectrum of the escaping flux (purple dashed line). The spectrum of accelerated particles at the wind termination shock, as suggested by equation (3) and as naturally predicted by DSA in a finite system, is a power-law of index s with maximum energy |$E_{\rm max} \simeq 1 \, \rm EeV$| and does not show any relevant additional spectral feature. On the other hand, the particle spectrum gradually steepens in the downstream region moving from the wind shock to the forward shock as a result of escape and pγ energy losses. In the downstream region, energy losses play a crucial role in shaping the spectrum of particles escaping the system. In particular, the photomeson interactions on the big blue bump occur faster than escape at |${\sim}10^{17}\,{\rm eV}$| as one can also deduce from Fig. 2. This results in a dip in the spectrum at such energy, whereas at higher energies the escape is more efficient so that the spectrum hardens at the highest energies.

A comment on the maximum energy is in order: the exponential function regulating the cut-off, Γ_cut, accounts for the geometry of the system and loss mechanisms, so that it cannot be simplified as a ratio E/E_max. Therefore, here we define E_max as the energy where p^sf_sh is suppressed by one e-fold. In what follows we describe in detail the impact of different realizations of the system to the maximum energy.

3.1 Impact of parameters on the maximum energy

A qualitative estimate of the maximum energy set by the geometry of the system can be obtained by comparing the upstream diffusion length, D₁/u₁, with the size of such region, R_sh (see also Morlino et al. 2021; Peretti et al. 2022, for additional discussion). Since at the highest energies r_L is already larger than l_c one can write the maximum energy as follows:

$$\begin{eqnarray} E_{\max} &=& q_{\rm B}\sqrt{\frac{6}{c}} \left[\frac{\epsilon_{\rm B} \dot{M} l_{\rm c}}{R_{\rm sh}} \right]^{1/2} u_1 \\ & \simeq& 1.4 \, {\rm EeV} \, \bigg(\frac{\epsilon_{\rm B}}{0.05} \frac{\dot{M}}{10^{-1}\,{\rm M}_{\odot}\,{\rm yr}^{-1}} \frac{l_{\rm c}}{10^{-2} \, {\rm pc}} \frac{1 \, \rm pc}{R_{\rm sh}} \bigg)^{1/2} \frac{u_1}{0.2 \, c}. \end{eqnarray}$$

(4)

As one can see from equation (4), the maximum energy for DSA at the wind shock of UFOs turns out to be of the order of EeV for standard values of parameters.

Table 2 highlights the impact of different parametric assumptions on the maximum energy. In particular, we see that, according to equation (4), E_max scales roughly linearly with u₁ and with the square root of |$\dot{M}$| and ϵ_B. The impact of l_c on E_max can be understood as follows: when |$l_{\rm c} \gg 10^{-2} \, \rm pc$|⁠, the diffusion coefficient is much larger than the benchmark scenario so that the diffusion length reaches the size of the system at lower energies; when |$l_{\rm c} \ll 10^{-2} \, \rm pc$| the energy at which the diffusion coefficient changes regime (from the standard quasi-linear theory ∼E^{2 − δ} to the small pitch-angle scattering regime ∼E²) shifts to lower energies thereby resulting in a larger value of D at the highest energies. Therefore, since at the highest energies, diffusion dominates, a local maximum in E_max appears for |$l_{\rm c} \simeq 10^{-2} \,{\rm pc}$|⁠. Different assumptions on the slope of the turbulence cascade (Kolmogorov-like) have a negligible impact on E_max because at the highest energies, where r_L > l_c, diffusion proceeds in a different regime.

Table 2.

Impact on the maximum energy of a parameter variations. All parameters are set to the benchmark UFO values shown in Table 1 except for those indicated in the first two columns. The last row shows the result for benchmark values for comparison.

Parameter(s)	Variation(s)	Emax [EeV]
u1/c	0.03/0.1/0.3	|$0.03/0.31/1.86$|
|$\dot{M} [{\rm M}_{\odot } \, {\rm yr}^{-1}]$|	10−2/1	|$0.29/2.82$|
ϵB	0.01/0.1	|$0.53/1.41$|
lc[pc]	|$3 \cdot 10^{-3}/10^{-1}$|	|$0.81/0.24$|
δ	5/3 (Kolmogorov)	1.02
tage[yr]	|$10^2/10^4/10^5/10^6$|	|$0.58/1.12/0.88/0.63$|
|$n_{0} [\rm cm^{-3}]$|	103/105/106	1.11/0.75/0.26
Urad	none/double	|$2.04/0.77$|
(⁠|$\dot{M},u_1$|⁠)	pessimistic/optimistic	|$0.01/4.53$|
no variations (benchmark)		1.06

We explored a wide range of possible ages of the system: from 10² yr, where it overcomes all relevant time-scales being still consistent with the stationary assumption, up to 1 Myr, where it becomes comparable with the AGN duty cycle. We observe that the age of the system does not have a strong impact on E_max, which is affected by less than a factor 2 for the wide range of alternatives considered. Interestingly, we notice that for an age |${\gtrsim}10^4 \,{\rm yr}$|⁠, the matter accumulated in the SAM layer becomes calorimetric for low energy particles. This could result in a further hardening of the escaping flux.

The density of the external medium has a direct impact on the size of the system, while the effect on the maximum energy is non-trivial. For densities n₀ ≲ 10⁵, the effect on the maximum energy is moderate and follows approximately equation(4) (see also equation (A1)). On the other hand, higher densities (n₀ ≳ 10⁵) result in a smaller size of the system (⁠|$R_{\rm sh} \lesssim 0.4 \, \rm pc$|⁠) and greater amount of matter accumulated over the age of the system. For such a range of density, we observe that the SAM layer becomes calorimetric for low energy particles similar to old UFOs. Finally, we observe a deviation from equation (4) when the external density is as large as |$n_0= 10^6 \, \rm cm^{-3}$| and the wind shock radius reduces to |$R_{\rm sh}\approx 0.2 \, \rm pc$|⁠. In fact, for such a reduced size, the pγ interactions start to be efficient on the big blue bump photon field of the AGN.

Interestingly, as also highlighted in Fig. 2, different assumptions in the photon field highlight a trend which suggests that the pγ interactions on the infrared field of the torus regulate the maximum energy. In particular, E_max increases by a factor 2 when the photon field is removed, while it decreases when a stronger photon field is considered. This suggests that the infrared field of the torus could play a crucial role in regulating the maximum energy achievable in UFOs. We finally explore the combined effect of maximum (minimum) values of u₁ and |$\dot{M}$| corresponding to a plausible optimistic (pessimistic) scenario. In this context, one can see that UFOs can be responsible for proton acceleration with E_max ranging from 10 PeV up to 5 EeV. In particular, the objects at the high luminosity end of a hypothetical luminosity function of UFOs are candidate acceleration sites of UHECRs where protons could reach a few EeV. Heavier nuclei could be accelerated to higher total energies, provided they survive photodisintegration. The latter possibility depends on the photon background present at the acceleration site and on the relative distance between R_sh and R_fs.

Our benchmark scenario as well as the scenarios characterized by larger size (resulting from smaller n₀ or larger t_age) are not affected by our approximate treatment of the spatial behaviour of the torus photon field because |$R_{\rm sh} \gtrsim 1\,{\rm pc}$|⁠. On the other hand, scenarios characterized by high density and/or younger age could, in principle, be partially affected since R_sh would be located at distances smaller than the typical size of the dusty torus. Nevertheless, we do not expect a strong impact of different radial dependence of the torus photon field for the following reasons: 1) HE particles cannot access the innermost upstream region due to the dominant effect of advection; 2) since it also scales as r⁻², the optical big blue bump becomes the dominant thermal photon field for energy losses and limits the maximum energy to be lower than 300 PeV for |$R_{\rm sh}\lesssim 0.2 \, \rm pc$|⁠. Finally, the gamma-ray absorption is moderately affected by different assumptions on the torus photon field because sub-TeV photons are absorbed on the optical-UV field while beyond 10 TeV they are absorbed by the EBL (Franceschini & Rodighiero 2017) en route to Earth.

4 GAMMA-RAYS AND HE NEUTRINOS FROM UFOS AND CONSTRAINTS TO THEIR LOCAL DENSITY

The gas swept-up from the dense environment of the SMBH as well as the strong radiation field of the AGN can make hadronic interactions observationally relevant in UFOs. Since interactions are copiously taking place in the UFO wind bubble, a high level of gamma-ray and HE neutrino emission can be expected. Fig. 4 illustrates the gamma-ray (thick blue line) and the single-flavour neutrino flux (dotted red for pp and dot–dashed orange for pγ) flux expected from the benchmark UFO placed at a redshift |$z = 0.013$|⁠. In particular, the gamma-ray flux is compared with the typical UFO spectral energy distribution (SED) as found in Ajello et al. (2021).

Figure 4.

Gamma-ray (thick blue line) and HE neutrino flux produced in the benchmark UFO (⁠|$R_{\rm sh}\simeq 0.8 \, \rm pc$|⁠) via pp (dotted red line) and pγ (orange dot–dashed line) interactions. The acronym SAM (SW) refers to the SAM SW. The thin blue-to-cyan lines represents the gamma-ray flux computed, respectively, for scenario A, B, and C (see Table 1) in order to illustrate the dependence of the gamma-ray absorption on the bubble expansion (where |$R_{\rm sh} \simeq 2, 5, {\rm and}\ 8 \, {\rm pc}$|⁠, respectively). The UFO is assumed to be located at |$z = 0.013$| in order to be directly compared with the best-fit UFO SED provided in Ajello et al. (2021), where the grey band represent the 1σ uncertainty band of such a best-fit UFO SED.

Despite the fact that the benchmark scenario represents an average UFO in terms of power and maximum energy, the gamma-ray flux of the benchmark UFO (⁠|$R_{\rm sh} \simeq 0.8 \, \rm pc$| and |$R_{\rm fs} \simeq 3 \, \rm pc$|⁠) cannot be representative of the whole class due to the strong radial-dependence of the γγ absorption. Therefore, in Fig. 4, we also compare it with the expected gamma-ray fluxes as predicted from scenario A (⁠|$R_{\rm sh} \simeq 2 \, \rm pc$| and |$R_{\rm fs} \simeq 8 \, \rm pc$|⁠), B (⁠|$R_{\rm sh} \simeq 5 \, \rm pc$| and |$R_{\rm fs} \simeq 22 \, \rm pc$|⁠), and C (⁠|$R_{\rm sh} \simeq 8 \, \rm pc$| and |$R_{\rm fs} \simeq 40 \, \rm pc$|⁠) as described in Table 1. Scenarios A, B, and C do not differ from the benchmark scenario in terms of total power but illustrate alternative realizations of it, having a larger size resulting from a longer evolution in a less dense environment. One can observe that these scenarios enhance the gamma-ray emission above ∼10 GeV due to a weaker gamma-gamma absorption and allow a better agreement with the UFO sample observed by Ajello et al. (2021).

Regardless of the age of the system, gamma-rays of energy greater than a few TeV are completely absorbed by the infrared radiation field of the torus and on the EBL. Therefore, pp neutrinos in the 10 TeV–1 PeV band as well as pγ neutrinos in the energy band 10² TeV–10² PeV would be produced without their gamma-ray counterpart. UFOs are thus expected to be bright neutrino sources featuring spectra as hard as ∼E⁻², while being highly opaque to TeV (and possibly 10–10² GeV) gamma-rays (Inoue et al. 2022).

If the AGN activity during the duty cycle were intermittent, particles leaving the active UFO could end up in a larger scale expanding slower wind that might result from a previous UFO. In this scenario, depending on the local diffusion coefficient experienced by HE particles, adiabatic losses might play a role. However, if the local diffusion were comparable to the one in the undisturbed ISM or somehow in between the undisturbed ISM in our Galaxy (⁠|$D \approx 10^{28} E_{\rm GeV}^{1/3} \, \rm cm^2 \, s^{-1}$|⁠) and the active UFO, the effect on the highest energy particles would be moderate. On the other hand, if the UFO lifetimes were longer than our assumption for the benchmark scenario, the SAM layer could increase luminosity with time (see also Peretti et al. 2022) and eventually turning calorimetric for low-energy particles.

UFOs could be common in nearby luminous infrared galaxies (LIRGs), such as active starburst galaxies and Seyfert galaxies. However, the abundance and distribution of these objects throughout the Universe as well as their luminosity function are poorly known. Therefore, in what follows, we estimate the order of magnitude of their diffuse multimessenger emission in terms of EeV cosmic rays and associated HE neutrinos and gamma rays. Since the horizon for the Bethe–Heitler pair-production suffered by UHECRs on the cosmic microwave background (CMB) is placed beyond z > 2, we neglect such a loss mechanism in our calculations.

We first assume as a prototype UFO the EeV-atron described by the benchmark scenario presented in Table 1. As discussed in Appendix  B, the escaping flux is regulated by the interplay between diffusion, advection, and energy losses. However, despite its complex analytic expression, assuming an ∼E⁻² spectrum the power contained by the escaping particles can be approximated as follows:

$$\begin{eqnarray} L_{\rm CR} & =& \int {\rm d}p \, 4 \pi p^2 [{\rm pc}] 4 \pi R_{\rm esc}^2 j_{\rm esc} \simeq \, \frac{3}{4} \, \xi _{\rm CR} \, \eta _{\rm loss} \, \dot{M} \, u_1^2 \\ & \simeq& 2 \cdot 10^{43} \, {\eta _{\rm loss}} \frac{\xi _{\rm CR}}{0.05} \frac{\dot{M}} {0.1\,{\rm M}_{\odot}\,{\rm yr}^{-1}} \left(\frac{u_{1}}{0.2 \, c} \right)^2 \frac{\rm erg}{\rm s}, \end{eqnarray}$$

(5)

where j_esc is the escaping flux of protons as defined in Appendix  B, η_loss ≤ 1 is an age-dependent parameter, which accounts for the relative reduction in the escaping flux due to energy losses, while the other parameters are normalized to the values shown in Table 1. The CR luminosity is related to the CR spectral injection rate |$\mathcal {Q}_{\rm CR}$| (units of GeV⁻¹s⁻¹) as |$L_{\rm CR}=\int {\rm d}E \, E \mathcal {Q}_{\rm CR}(E)$|⁠. For simplicity, we assume in the following that the CR emission follows E⁻² from GeV to EeV, such that |$E_p^2 \mathcal {Q}_{\rm CR} \simeq \chi L_{\rm CR}$| with χ ≡ 1/ln (EeV/GeV) ≃ 0.05.

In general, the locally observed CR spectrum ϕ_CR (units of GeV⁻¹cm⁻²s⁻¹sr⁻¹) is related to the spectral emission rate of extragalactic sources via a set of transport equations. For CR protons in the EeV energy range, we can assume that the transport is dominated by continuous energy loss due to the expansion of the Universe, while we neglect the effect of intergalactic magnetic fields. Following the notation of Ahlers & Halzen (2018), we can estimate the local contribution of UFO EeV-atrons as:

$$\begin{eqnarray} \big[E^2_p \phi _{\rm CR}\big]_{{\rm EeV}}\simeq \frac{\xi_z}{4 \pi}\frac{c}{H_0} \rho_0 \big[E^2_{p}\mathcal {Q}_{\rm CR}\big]_{{\rm EeV}}. \end{eqnarray}$$

(6)

The factor ξ_z is of order unity accounting for the integral in redshift of the source distribution. In particular, ξ_z ≃ 0.5 (2.6, 7.8) for a flat (star-formation rate, AGN) distribution, where the factor for the AGN distribution has been computed for Log L_X = 44–45 adopting data from Ueda et al. (2014). The parameter ρ₀ represents the local comoving density of sources for which we assume |$\rho _0 = 10^{-5} \, \rm Mpc^{-3}$| as a reference. Such a value is often quoted as a typical density of AGN with X-ray luminosity of the order of |$L_X \simeq 10^{44} \, \rm erg \, s^{-1}$| (see e.g. Ueda et al. 2014; Murase & Waxman 2016; Fiore et al. 2017). It is worth mentioning that a source density |${\sim}10^{-4}\!-\!10^{-5}\,{\rm Mpc}^{-3}$| matches also the number density inferred for powerful starbursts such as luminous and ultra-luminous infrared galaxies (see Gruppioni et al. 2013; Peretti et al. 2020; Condorelli et al. 2023), which are currently also considered to be plausible hosts for UHECR accelerators (Aab et al. 2018). Expressing the contribution of EeV CR protons to the spectral emission as |$\big[E_p^2 \mathcal {Q}_{\rm CR}\big]_{{\rm EeV}} = \chi L_{\rm CR}$| with χ ≃ 0.05 we arrive at:

$$\begin{eqnarray} \big[E^2_p\phi _{\rm CR}\big]_{{\rm EeV}} \simeq 3\cdot 10^{-7} \frac{\rho _{0}}{10^{-5}\,{\rm Mpc^{-3}}} \, \frac{\xi _{z}}{2.6} \, \frac{L_{\rm CR}}{10^{43}\,{\rm erg\,s}^{-1}} \, \rm \frac{GeV}{cm^2\, s \, sr}\, . \end{eqnarray}$$

(7)

In comparison, the observed CR spectrum has value of [E²ϕ_CR]_EeV ≃ 2 · 10⁻⁷ GeV cm⁻²s⁻¹ sr⁻¹, very close to our estimate.

The associated all-flavour neutrino spectral injection rate |$\mathcal {Q}_{\nu }$| resulting from pp interactions can be related to |$\mathcal {Q}_{\rm CR}$| as:

$$\begin{eqnarray} E^2_\nu \mathcal {Q}_{\rm all\, \nu }&\simeq& \frac{1}{2} \kappa \, \sigma _{pp} \, c \, t_{\rm age}\, n_{\rm ISM} \, \eta _{\rm loss}^{-1} \, E_p^2 \mathcal {Q}_{\rm CR} \\&\simeq& 0.07 \frac{n_{\rm ISM,b}}{10^{4}\,{\rm cm}^{-3} } \frac{t_{\rm age}}{10^{3}\,\rm yr} \, \eta _{\rm loss}^{-1} \, E_p^2 \mathcal {Q}_{\rm CR}, \end{eqnarray}$$

(8)

where κ ≃ 0.5 is the inelasticity of pp interactions with cross section σ_pp ≃ 3 · 10⁻²⁶ cm², and the CR proton energy is related to neutrino energy as E_ν ≃ 0.05E_p. Notice that in this order of magnitude estimate, we consider only neutrinos resulting from pp interactions since they dominate up to at least |$10^2 \, \rm TeV$|⁠. Neutrinos from pγ as previously shown, do not feature a different order of magnitude in the HE SED therefore the pp neutrino flux can be assumed as a good approximation also for the order of magnitude flux reached by the pγ component.

Again, following Ahlers & Halzen (2018), we can now estimate the isotropic neutrino flux as:

$$\begin{eqnarray} E^2_\nu \phi _{\rm all \nu }= \frac{\xi _z}{4 \pi } \frac{c}{H_0}\rho _0 E^2_\nu \mathcal {Q}_{\rm all \nu}. \end{eqnarray}$$

(9)

Adopting in equation (9), the expressions for the CR and neutrino luminosity as described in equations (5) and (8) together with |$E_p^2 \mathcal {Q}_{\rm CR} \simeq \chi L_{\rm CR}$|⁠, one obtains an estimate of the total single flavour neutrino flux at the Earth:

$$\begin{eqnarray} E^2_\nu \phi _{\rm all \nu} &\simeq& 2\cdot 10^{-8} \frac{\rm GeV}{\rm cm^2\,s \, sr}\\ &&\times \, \frac{\rho _{0}}{10^{-5}\,{\rm Mpc}^{-3}} \frac{\xi _{z}}{2.6} \frac{n_{\rm ISM,b}} {10^{4}\,{\rm cm}^{-3}} \frac{t_{\rm age}}{10^{3}\,\rm yr} \frac{L_{\rm CR}}{10^{43}\,{\rm erg\,s}^{-1}}. \end{eqnarray}$$

(10)

This flux is comparable to the level of the diffuse neutrino flux observed by IceCube (see Abbasi et al. 2022b) for energies larger than 100 TeV.

These estimates indicate that if UFOs were as abundant as typical non-jetted AGN and powerful starbursts, they could potentially be the dominant sources of cosmic rays at EeV, while also strongly contributing to the diffuse neutrino flux observed by IceCube above 100 TeV. In this context, the associated gamma-ray flux would not be expected to contribute substantially to the diffuse gamma-ray flux observed by Fermi-LAT (Ackermann et al. 2015) due to the flat spectral shape and the strong absorption taking place inside the source environment above ∼10² GeV that would limit the energy budget of the electromagnetic cascade in the propagation of the gamma rays to the Earth. In a scenario of minimal absorption (scenario C), one can expect at most a gamma-ray flux at ≲1 TeV at the same level of the neutrino one, namely |$E_{\gamma }^2 \phi _{\gamma } (E_{\gamma }) \sim E_{\nu _{\mu }}^2 \phi _{\nu _{\mu }} (E_{\nu _{\mu }})$|⁠.

Directly observed UFOs as well as X-ray AGN and luminous infrared galaxies, where a UFO can be obscured, represent a promising source class where UHECRs could be accelerated. The detection of a flat (∼E⁻²) neutrino spectrum extending in the PeV range from the innermost core of these galaxies could serve as a hint to discriminate whether DSA is taking place as we propose in this work.

5 APPLICATION TO NGC1068

In what follows, we specialize our calculations to the nearby Seyfert galaxy NGC1068 located at a distance |$D_{\rm L} = 14 \,{\rm Mpc}$|⁠. Notice that NGC1068 is Compton thick AGN (Matt et al. 1997). Therefore, a clear detection of an UFO in its nuclear region is extremely challenging, even though Mizumoto, Izumi & Kohno (2019) reported some indications in that sense. Despite one cannot have a compelling evidence of an UFO in NGC1068 in the following, we explore the possibility that this galaxy hosts an obscured UFO in its nuclear region. We compute the gamma-ray and associated neutrino emission through pp and pγ interactions by adopting |$\dot{M} = 2 \cdot 10^{-1} \, \rm {\rm M}_{\odot } \, yr^{-1}$|⁠, |$u_1 = 0.1 \, c$|⁠, and |$l_{\rm c} =10 \, \rm pc$|⁠, while assuming all other parameters as in the benchmark scenario. In particular, the coherence length assumed here to be as large as the system size (⁠|$2 \, R_{\rm fs} \simeq l_{\rm c} \simeq 10 \, \rm pc$|⁠) results in a maximum energy of approximately |$E_{\rm max} = 5 \, \rm PeV$|⁠. This is found in agreement with the trend presented in the previous section (Section 3.1). Notice that the value of l_c required for this calculation is considerably larger than in our benchmark scenario discussed in Section 3. In fact, assuming smaller values of l_c leads to a neutrino spectrum extending up to several PeV with an approximately ∼E⁻² spectrum, in contradiction with the IceCube sensitivity.

Fig. 5 illustrates the multimessenger flux produced by the accelerated particles in the system and under the assumption that the pressure of accelerated particles is ∼5 per cent of the ram pressure at the wind termination shock. One can see that the gamma-ray emission (thick lines) is dominated by the pp contribution in the SAM (blue line), whereas the emission from the SW (cyan line) is more than two order of magnitude dimmer. As discussed in Section 4, the strong photon field associated with the accretion disc and torus makes the source opaque to gamma rays above a few tens GeV, so that all TeV photons are completely absorbed. The neutrino flux is dominated by the pp channel in the range from GeV up to ∼10² TeV, where it features a reduction due to the maximum energy. Photomeson interactions take also place in the UFO environment. However, they produce a negligible impact on the spectrum. Overall, the neutrino flux shows a remarkable flat spectrum over more than five orders of magnitude, where the associated gamma-ray counterpart gets absorbed beyond ∼10 GeV. It is possible to notice that 1) a UFO could dominate the gamma-ray flux observed by Fermi-LAT and 2) in the TeV range the UFO could contribute from a few up to ∼10 per cent of the flux measured by IceCube (Abbasi et al. 2022a) leaving room for other possible sources such as the AGN corona, the molecular outflow or the starburst ring. A standard UFO seems therefore disfavoured for explaining the level of neutrino flux observed by IceCube in light of the stringent upper limits imposed by MAGIC as well as the detected Fermi-LAT flux at lower energies.

Figure 5.

Multimessenger emission for the case of NGC1068. The thick blue line represents the gamma-ray emission dominated by pp interactions in the SAM, while the red dotted line is the associated neutrino flux. The emission from the SW (gamma-ray in cyan dashed and neutrino in magenta dot–dot–dashed) as well as the photomeson (orange dot–dashed) emission are subdominant. The model prediction is compared with Fermi-LAT (Abdollahi et al. 2020), MAGIC (Acciari et al. 2019), and IceCube (Abbasi et al. 2022a) data.

5.1 Forward shock scenario for NGC1068

As an alternative scenario to the acceleration at the wind termination shock one could explore the same UFO at an early stage during which the forward shock is expected to play the most relevant role in terms of particle acceleration.

The forward shock expanding in the unperturbed external medium of density n₀ ≃ 10⁴ cm⁻³ can indeed foster particle acceleration through DSA. Let us consider a young UFO that has not entered the deceleration phase yet, namely, it would be expanding with constant velocity ∼u₁. In this context, the forward shock is expected to be strong (⁠|${\cal M}\gg 1$|⁠) thereby efficient in accelerating particles. The spectrum of accelerated particles at such shock can be written as f(p) = A(p/p₀)^−α exp (− p/p_max), where the normalization A is estimated assuming that a fraction ξ_fs ∼ 0.1 of the ram pressure |$m_{p}n_{0} u_1^2$| is converted into energized particles. In the test-particle limit, α = 4, while an upper limit on the maximum momentum can be estimated with the Hillas criterion, leading to:

$$\begin{eqnarray} p_{\rm max, fs} \simeq \frac{\xi _{\rm fs}}{0.1}\, \frac{R_{\rm fs}}{\text{pc}}\, \frac{u_1}{10^8\,{\rm cm\,s^{-1}}}\, \frac{B}{\mu \text{G}} 10^3\,{\rm GeV\,c}^{-1}. \end{eqnarray}$$

(11)

An acceleration site with |$R_{\rm fs} \sim 10^{-2} \,{\rm pc}$| and |$B \sim 10^2\!-\!10^3 \,\mu$|G allows for p_max ≳ 10⁴–10⁵ GeV c⁻¹, in principle sufficient for the pp and pγ processes to produce gamma rays and neutrinos in an energy range accessible to current instruments.

Assuming that the accelerated protons fill at most a volume V ∼ 4 · 10⁻⁶(R_fs/10⁻²pc)³ pc³, in which the interactions with the gas and the AGN photon field are taking place, the UFO is expected to produce a gamma-ray luminosity (before any absorption effects are taken into account) of the order of L_γ(10³–10⁴GeV) ≃ 10⁻¹⁵–10⁻¹⁴ erg cm⁻¹² s⁻¹, and a neutrino luminosity of the order of L_ν(10³–10⁴ GeV) ≃ 10⁻¹⁶–10⁻¹⁵ erg cm⁻² s⁻¹, indicating that the emission from the forward shock is expected to be subdominant compared to the one from the wind termination shock happening at later times.

Late stage UFOs, as discussed in Ajello et al. (2021), could be luminous enough to be detected in the GeV range by space-based telescopes such as Fermi-LAT. However, this would require the central engine to be active for a time |$t \gg 10^3 \, \rm yr$|⁠. On the other hand, high acceleration efficiency at the wind termination shock can be expected to be found soon after the deceleration phase has begun.

5.2 Discussion on NGC1068

UFOs are characterized by a prominent gamma-ray (up to ∼10–10² GeV) and neutrino emission (typically up to ∼10² TeV) resulting from pp interactions. Therefore, in the standard test particle regime considered in this work, the resulting gamma-ray and neutrino spectra will feature roughly the spectral slope (∼E⁻²) of their parent protons. We observe that from the energetic point of view, the level of neutrino flux observed is compatible with the amount of power that UFOs can supply on average suggesting the AGN as a plausible origin for such an emission. However, the strongest limitation comes from the level of gamma rays measured by Fermi-LAT. In fact, even though the AGN radiation field can absorb efficiently TeV photons, GeV gamma-rays come basically unabsorbed. Indeed, as pointed out by several authors (see e.g. Inoue et al. 2020; Murase et al. 2020; Kheirandish et al. 2021; Eichmann et al. 2022; Inoue et al. 2022; Murase 2022), such level of neutrinos is not compatible with a pp scenario with a ∼E⁻² spectrum unless the emission comes from a region optically thick to GeV and sub-GeV gamma rays, namely the nearest neighbourhood of the SMBH or AGN-corona having a size of |${\lesssim}10^2 \, R_{\rm s}$| (where R_s is the Schwarzschild radius).

Even though it is not clear whether the launching radius of a UFO can be localized at such a small distance from the SMBH, we explored under which conditions one could expect a UFO to be confined close to the AGN corona for a sufficient amount of time during its deceleration phase. We found that, in order to reproduce the level of neutrino flux inferred by IceCube (Abbasi et al. 2022a) without exceeding the gamma-ray flux, the energy budget would not exceed standard values typical of UFOs (⁠|$\dot{E}\lesssim 10^{45} \, \rm erg \, s^{-1}$|⁠) but the requirement in terms of average gas density of the external medium would be |$n_{\rm ISM} \gtrsim 10^{10} \, \rm cm^{-3}$|⁠. Such a density could be compatible with the high column density inferred for this source (see e.g. Matt et al. 1997) therefore this could be a plausible scenario to account for the observed neutrino flux. A detailed modelling of a confined wind expanding in the ultra-dense medium goes beyond the scope of this work therefore it is left for a follow-up investigation.

The recent anisotropy study carried out by the Pierre Auger Observatory suggests that Seyfert and Starburst galaxies (or objects with a spatial distribution related to these sources) may be related to this anisotropy (Aab et al. 2018). It is intriguing that NGC1068 is the fourth most relevant object in such catalogs, providing a contribution of about 5–10 per cent. However, our calculations suggests that if the gamma-ray flux detected by Fermi-LAT is associated with the UFO activity in this source, then the maximum energy is bound to be too low to be relevant for UHECRs. Viceversa, if UHECRs are to be produced in the UFO, a different source for the gamma rays has to be found.

6 CONCLUSIONS

In this work, we investigated the potential of DSA at the wind termination shocks of ultra-fast outflows (UFOs) in the core of active galaxies. We developed a model of acceleration and transport of particles in the wind bubbles excavated by UFOs and we studied the multimessenger implications in terms of escaping particles and HE photons and neutrinos produced through hadronic interactions.

We found that protons can be accelerated up to the EeV range and that the far infrared photon field of the torus could play a dominant role in setting the maximum energy in such sources. In addition, the transport condition in the UFO environment could result in a spectral hardening of the escaping flux of the cosmic rays at the highest energies. Such energetic and spectral properties are crucial in light of the recent results obtained by the Pierre Auger Observatory, which suggest that the sources of UHECRs could be characterized by hard spectra (see e.g. Aab et al. 2017, and references therein).

UFOs are extremely interesting in terms of multimessenger emission since, in standard conditions, they are expected to shine in gamma rays up to ∼10 GeV, while being opaque beyond 10–10² GeV depending on the parametric configuration. While gamma rays are efficiently absorbed, HE neutrinos from pp and pγ interactions will be copiously produced with a spectrum extending up to ∼10² PeV featuring a spectral slope as hard as the one of their parent protons. Such a property is particularly interesting since UFOs have the potential and possibly the number density in the Universe to simultaneously dominate the CR spectrum at the ankle and the diffuse neutrino flux observed beyond ∼10² TeV.

We finally applied our model to the UFO that could be present in the nuclear region of NGC1068, and we found out that, if confirmed, in favourable parametric conditions, a well-developed spherical UFO wind bubble could dominate the gamma-ray flux observed by Fermi-LAT, while it is unlikely to contribute more than ∼10 per cent of the total neutrino flux observed by IceCube at TeV. On the other hand, the base of the outflow could be a plausible region, where protons can be injected and produce HE neutrinos while the associated GeV gamma-ray counterpart can be efficiently reprocessed to lower energies (see also Inoue et al. 2022).

The present model focused on the acceleration and transport of protons in UFOs, while the study of potential implications for primary and secondary electrons as well as electron–positron pairs is left for future work. UFOs are in fact expected to be perfect electron calorimeters given the dominant synchrotron and inverse Compton time-scales compared to the escape time-scale, and the case of leptons would require a specific treatment. We computed the typical time-scales for electrons, and we found that primary electrons cool very rapidly without reaching TeV energies. Secondaries as well as electron–positron pairs are also expected to cool mostly via synchrotron losses since beyond ∼1 TeV, the interaction with the big blue bump takes place in the Klein–Nishina regime. In such energy range, the electromagnetic cascade is expected to be synchrotron-dominated for standard parametric assumptions (ϵ_B ≳ 10⁻²). Therefore, one could expect the leptonic emission to produce some possible signatures in the hard X-ray to soft gamma-ray energy band and in radio.

Heavy nuclei should also be implemented in a future work since their transport, similar to electrons, requires the inclusion of calorimetric conditions, continuous energy losses as well as fragmentation and re-injection at lower masses. While our focus has been on the proton population, it is worth mentioning that nuclei can be expected to be co-accelerated in UFOs with a rigidity dependence in the maximum energy and possibly a higher efficiency in the injection of energetic heavier elements (see e.g. Caprioli, Yi & Spitkovsky 2017). In particular, while protons only partially suffer energy losses, heavy nuclei of electric charge Z are expected to efficiently fragment on the AGN photon field. Therefore, one can expect that at early time, the majority of heavy nuclei would be reprocessed, while at later time they could start escaping efficiently with an energy as large as |${\sim}Z \cdot \rm EeV$|⁠.

With the present work, we propose UFOs as candidate sources of UHECRs and efficient HE neutrino emitters possibly opaque to gamma rays beyond a few tens GeV. Such properties make UFOs a remarkable source class not only for the UHECRs detected by the Pierre Auger Observatory but also for the diffuse HE gamma-ray and neutrino fluxes observed, respectively, by Fermi-LAT and IceCube.

ACKNOWLEDGEMENTS

The research activity of EP and MA was supported by Villum Fonden (project No. 18994). EP was also supported by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 847523 ‘INTERACTIONS’. FGS acknowledges financial support from the PRIN MIUR project ‘ASTRI/CTA Data Challenge’ (PI: P. Caraveo), contract 298/2017. GM is partially supported by the INAF Theory grant 2022 ‘Star Clusters as cosmic ray factories’. EP is grateful to Antonio Condorelli for insightful discussions.

DATA AVAILABILITY

No data has been analysed or produced in this work.

APPENDIX A: PROPERTIES OF THE AGN WIND BUBBLE

The dynamics of the wind shock and forward shock during the deceleration phase of a wind bubble have been studied by Koo & McKee (1992a, b). The location of the wind shock is described by the following relation:

$$\begin{eqnarray} R_{\rm sh} &=& 23.0 \, { \rm pc}\bigg(\frac{t_{\rm age}}{1\,{\rm Myr}}\bigg)^{2/5}\bigg(\frac{\dot{E}}{10^{38}\,{\rm erg}\,{\rm s}^{-1}}\bigg)^{3/10} \\ &&\times \, \bigg(\frac{n_{0}}{1\,{\rm cm}^{-3}}\bigg)^{-3/10} \bigg(\frac{u_{1}}{10^8\,{\rm cm}\,{\rm s^{-1}}}\bigg)^{-1/2}, \end{eqnarray}$$

(A1)

where t_age is the age of the system, |$\dot{E}$| the kinetic power (⁠|$\dot{E}= \dot{M}u_1^2 /2$|⁠), n₀ the external medium density, and u₁ the upstream wind speed. The forward shock radius is found as follows:

$$\begin{eqnarray} R_{\rm fs} = 76.0 \, {\rm pc}\bigg(\frac{t_{\rm age}}{1\,{\rm Myr}}\bigg)^{3/5} \bigg(\frac{\dot{E}}{10^{38}\,{\rm erg}\,{\rm s}^{-1}}\bigg)^{1/5}\left(\frac{n_{0}}{1\,{\rm cm}^{-3}}\right)^{-1/5}. \end{eqnarray}$$

(A2)

The radiation field (top panel) and the corresponding opacity (bottom panel) of the wind bubble are shown in Fig. A1. In particular, the typical AGN radiation field features a combination of thermal and non-thermal components. At around ∼10 meV, one can identify the infrared thermal field of the torus, while the second thermal component peaking at ∼10 eV is the thermal radiation from the accretion disc known as the big blue bump. Finally, starting from the keV range, the non-thermal tail produced by electrons in the AGN corona is present. Such a non-thermal tail can additionally feature a somehow different shape due to the mechanism known as Compton reflection. In this work, we decided to ignore such an effect since it has a strong dependence on the local parameters of the AGN, while having a negligible impact on the high-energy particles and their radiation. The gamma-ray opacity shown in the bottom panel clearly highlights the prominent impact of the big blue bump and the torus absorbing, respectively, in the energy ranges 10 GeV–TeV and beyond 10 TeV.

Figure A1.

Top panel: the AGN photon field for three different X-ray luminosity. Dotted red, dashed green, and dot–dashed blue represent the typical spectral energy distribution of an AGN having an X-ray luminosity of 10⁴², 10⁴⁴, and 10⁴⁶ |$\rm erg \, s^{-1}$|⁠, respectively. Bottom panel: gamma–gamma opacity of the wind bubble, where the different curves highlight the scenarios discussed in Section 3. In particular: the thick black line represents the benchmark scenario, while the thin lines describe the absorption in the Scenario A, B, and C.

APPENDIX B: SOLUTION TO THE TRANSPORT EQUATION

The transport equations (1) are solved separately in the upstream and the downstream regions and joined at the wind shock location as described in Morlino et al. (2021); Peretti et al. (2022). In what follows, we highlight the main analytical steps required to obtain the formal solution and the key aspects of the iterative algorithm we adopt to find the numerical solution.

B1 Upstream region

By integrating equation (1) from r = 0 to r < R_fs, one obtains:

$$\begin{eqnarray} r^2u_1f_1(r,p) = r^2D_1(r,p) \partial _rf_1(r,p){\,-\,}G_1(r,p){\,-\,}H_1 (r,p), \end{eqnarray}$$

(B1)

where the subscript ‘1’ refers to the upstream region, while the functions G₁ and H₁ have the following expressions:

$$\begin{eqnarray} G_1(r,p) = \frac{1}{3} \int _0^r\,{\rm d}r^{\prime} \left[- \frac{\partial \ln [p^3f_1(r^{\prime },p)]}{\partial \ln p} \right] f_1(r^{\prime },p) \, \partial _{r^{\prime}}\big[r^{\prime 2}u_1\big] , \\ \end{eqnarray}$$

(B2)

$$\begin{eqnarray} H_1(r,p) = \int _0^r\,{\rm d}r^{\prime} \, \lambda _1(r^{\prime },p) \, r^{\prime 2} \, f_1 (r^{\prime },p), \end{eqnarray}$$

(B3)

where |$\lambda _1= \tau _{pp}^{-1} + \tau _{{\rm p}\gamma }^{-1}+ \tau _{\rm BH}^{-1}$| is the CR energy loss accounting for pion production in pγ and pp interactions (see Kelner & Aharonian 2008) as well as Bethe–Heitler (BH) pair-production (see e.g. Gao et al. 2012). Defining the effective upstream velocity as:

$$\begin{eqnarray} V_{\rm eff,1} (r,p) = u_1 \left[1 + \frac{G_1(r,p)+H_1(r,p)}{u_1 r^2 f_1(r,p)} \right], \end{eqnarray}$$

(B4)

the solution of equation (B1) is straightforwardly obtained:

$$\begin{eqnarray} f_1(r,p) = f_{\rm sh}(p) \exp \left[- \int _r^{R_{\rm sh}}\,{\rm d}r^{\prime } \frac{V_{\rm eff,1}(r^{\prime },p)}{D_1(r^{\prime },p)} \right]. \end{eqnarray}$$

(B5)

B2 Downstream region

Similarly to the upstream case, also in the downstream, equation (1) is first approached through a spatial integral exploiting the boundary condition. Integrating from r > R_fs to r = R_fs, one obtains:

$$\begin{eqnarray} r^2u_2(r)f_2(r,p) = r^2D_2(r,p) \partial _rf_2(r,p) + R_{\rm fs}^2 j_{\rm esc}(p) + H_2 (r,p),\\ \end{eqnarray}$$

(B6)

where the subscript ‘2’ refers to the downstream region, while the function H₂ reads:

$$\begin{eqnarray} H_2(r,p) = \int _r^{R_{\rm fs}}{\rm d}r^{\prime } \, \lambda _2(r^{\prime },p) \, r^{\prime 2} f_2 (r^{\prime },p)\, . \end{eqnarray}$$

(B7)

As in the upstream region, it is convenient to define the effective velocity in the downstream region as:

$$\begin{eqnarray} V_{\rm eff,2}(r,p) = u_2(r) \bigg[1-\frac{H_2(r,p)}{r^2u_2(r)f_2(r,p)} \bigg]. \end{eqnarray}$$

(B8)

It is also useful to define the integral function I₂ as:

$$\begin{eqnarray} I_2(r,p) = \int _{R_{\rm sh}}^r \frac{{\rm d}r^{\prime }}{r^{\prime 2}} {\rm e}^{-\varphi _2(r^{\prime },p)}, \end{eqnarray}$$

(B9)

where φ₂(r, p) has the following expression:

$$\begin{eqnarray} \varphi _2(r,p) = \int _{R_{\rm sh}}^{r} {\rm d}r^{\prime } \frac{V_{\rm eff,2}(r^{\prime },p)}{D_2(p)}. \end{eqnarray}$$

(B10)

Integrating equation (B6), one obtains the expression for the escaping flux and the downstream solution as:

$$\begin{eqnarray} f_2(r,p) = f_{\rm sh}(p) \,{\rm e}^{\varphi _2(r,p)} \bigg[1-\frac{I_2(r,p)}{I_2(R_{\rm fs},p)} \bigg] \end{eqnarray}$$

(B11)

$$\begin{eqnarray} j_{\rm esc}(p) = f_{\rm sh}(p) \frac{D_2(p)}{R_{\rm fs}^2 I_2(R_{\rm fs},p)}. \end{eqnarray}$$

(B12)

We finally notice that the escaping flux j_esc can be rewritten as

$$\begin{eqnarray} j_{\rm esc}(p) = \eta _{\rm loss} \frac{u_2 f_{\rm sh}(p)}{1-\exp [-R_{\rm sh}u_2(1-R_{\rm sh}/R_{\rm fs})/D_2(p)]} \frac{R_{\rm sh}^2}{R_{\rm fs}^2}, \end{eqnarray}$$

(B13)

where η_loss is a parameter ≲1 accounting for energy losses.

B3 Solution at the shock

The shock solution is obtained by integrating equation (1) across an infinitely small layer embedding the wind shock. The result is the following:

$$\begin{eqnarray} [D_2 \partial _r f_2-D_1 \partial _r f_1]_{r=R_{\rm sh}}-\frac{u_1-u_2}{3} p \partial _p f_{\rm sh}(p) + Q_0 (p) = 0.\\ \end{eqnarray}$$

(B14)

Substituting equations (B1) and (B6) in the first term on the left-hand side, equation (B14) can be rewritten as:

$$\begin{eqnarray} \frac{s Q_0(p)}{u_1} = p \partial _p f_{\rm sh}(p) + sf_{\rm sh}(p) + s \Psi _{ l}(p) f_{\rm sh}(p) + s \Psi _{ e}(p) f_{\rm sh}(p),\\ \end{eqnarray}$$

(B15)

where the functions Ψ_k (k = l, e) are defined as:

$$\begin{eqnarray} \Psi _{l}(p) = \frac{{G}_1(R_{\rm sh},p)+{H}_1(R_{\rm sh},p)+{H}_2(R_{\rm sh},p)} {u_1 R_{\rm sh}^2 f_{\rm sh}(p)}, \end{eqnarray}$$

(B16)

$$\begin{eqnarray} \Psi _{e}(p) = \frac{\big[D_2 {I}_2^{-1}(R_{\rm sh},p)-R_{\rm sh}^2 u_2\big]}{u_1 R_{\rm sh}^2}. \end{eqnarray}$$

(B17)

Here, the subscripts l and e stands for loss and escape, respectively. Finally, by recognizing a total derivative on the right-hand side of equation (B15), the solution at the shock can be obtained as:

$$\begin{eqnarray} f_{\rm sh}(p) = \frac{s \eta n_1}{4 \pi p_{\rm inj}^3} \bigg(\frac{p_{\rm inj}}{p}\bigg)^s {\rm e}^{-\Gamma _l(p)}\,{\rm e}^{-\Gamma _e(p)}, \end{eqnarray}$$

(B18)

where:

$$\begin{eqnarray} \Gamma _{l(e)}(p) = s \int _{p_{\rm inj}}^p \frac{{\rm d}p^{\prime }}{p^{\prime }} \Psi _{l(e)}(p^{\prime }). \end{eqnarray}$$

(B19)

Notice that equation (B18) can be re-written in the compact form (see equation (3)) |$f_{\rm sh}(p) = C \, p^{-s} {\rm exp}[-\Gamma _{\rm cut}(p)]$|⁠, where |$C= s\, \eta \, n_1\, p_{\rm inj}^{s-3}/(4 \pi)$| and Γ_cut = Γ_l + Γ_e.

B4 Iteration algorithm

The solution to the transport equation on the two sides of the shock, f₁ (equation (B5)) and f₂ (equation (B11)), and at the shock, f_sh (equation (B18)), do not have a simple analytic form since they depend on each other through the functions V_eff, 1, V_eff, 2, and Ψ_l(e), respectively. A solution can be found found via an iterative algorithm.

We initialize the solutions for the set of functions (⁠|$f_{\rm sh}^{(0)},f_1^{(0)},f_2^{(0)}$|⁠) by the solutions resulting from thefollowing no-loss conditions: |$G_1^{(0)} = H_1^{(0)} = H_2^{(0)} = 0$|⁠; |$V_{\rm eff,1}^{(0)}(r,p)= u_1$|⁠; |$V_{\rm eff,2}^{(0)}(r,p)= u_2(r)$|⁠. This results in |$\Psi _l^{(0)} = 0$|⁠, while |$\Psi _e^{(0)}$| reduces to the following analytic form:

$$\begin{eqnarray} \Psi _e^{(0)}(p) = \frac{u_2/u_1}{ {\rm exp}\left[\frac{R_{\rm sh}u_2}{D_2(p)}\left(1-\frac{R_{\rm sh}}{R_{\rm fs}}\right)\right] - 1 }\, . \end{eqnarray}$$

(B20)

We start from this initial approximation and find iterative solutions by re-computing all functions with the set of solution of the previous iteration, namely:

$$\begin{eqnarray} && \left(f_{\rm sh}^{(i)},f_1^{(i)},f_2^{(i)}\right) \rightarrow \left(G_1^{(i+1)}, H_1^{(i+1)}, H_2^{(i+1)}\right) \\ &&\quad \rightarrow \left(V_{\rm eff,1}^{(i+1)}, V_{\rm eff,2}^{(i+1)}, \Psi _l^{(i+1)},\Psi _e^{(i+1)}\right) \\ &&\quad \rightarrow f_{\rm sh}^{(i+1)} \\ &&\quad \rightarrow \left(f_1^{(i+1)},f_2^{(i+1)} \right), \end{eqnarray}$$

where (i) and (i + 1) indicate the i-th and (i + 1)-th iteration. This algorithm is repeated until the phase space density at the n-th iteration f⁽ⁿ⁾ is indistinguishable from the solution found at the iteration (n − 1)-th, f^{(n − 1)}, namely when a convergence condition has been obtained.