Dynamical complexity in micro-scale disk-wind systems

Powerful winds at accretion-disk scales have been observed in the past 20 years in many AGN. These are the so-called ultrafast outflows (UFOs). Outflows are intimately related to mass accretion through the conservation of angular momentum, and they are therefore a key ingredient of most accretion disk models around black holes (BHs). Nuclear winds and outflows can also provide the feedback that regulates the joint BH and galaxy growth. We reconsidered UFO observations in the framework of disk-wind scenarios, both magnetohydrodynamic disk winds and radiatively driven winds. We studied the statistical properties of observed UFOs from the literature and derived the distribution functions of the ratio $\bar \omega$ of the mass-outflow and -inflow rates and the ratio $\lambda_w$ of the mass-outflow and the Eddington accretion rates. We studied the links between $\bar \omega$ and $\lambda_w$ and $\lambda={L_{bol}}/{L_{Edd}}$. We derived the typical wind-activity history in our sources by assuming that it can be statistically described by population functions. We find that the distribution functions of $\bar \omega$ and $\lambda_w$ can be described as power laws above some thresholds, suggesting that there may be many wind subevents for each major wind event in each AGN activity cycle, which is a fractal behavior. We then introduced a cellular automaton to investigate how the dynamical properties of an idealized disk-wind system change following the introduction of feedback rules. We find that without feedback, the system is overcritical. Conversely, when feedback is present, regardless of whether it is magnetic or radiation driven, the system can be driven toward a self-organized critical state.Our results corroborate the hypothesis that AGN feedback is a necessary key ingredient in disk-wind systems, and following this, in shaping the coevolution of galaxies and supermassive BHs.


Introduction
The coevolution of black holes (BH) and galaxies is one of the most discussed topics in modern extragalactic astrophysics.BHgalaxy coevolution aims to understand the interplay between active galactic nuclei (AGN), BH accretion, and the evolution of their host galaxies.At the very root of the problem lies the crucial question whether the supermassive BHs (SMBH) systematically found in massive galaxies nuclei are a coincidence or a necessity.If the growth of a nuclear SMBH is the result of independent processes that produce both BH and galaxy growth, then BHs are a coincidence.Alternatively, if the two systems influence their respective growth in a nonlinear feedback loop significantly, BHs may be a necessity.
Today, BHs are an observational evidence in quite a broad mass range, from a few to billions of solar masses (The LIGO/ VIRGO collaboration 2019; The Event Horizon Telescope Collaboration 2019; GRAVITY Collaboration 2022).We know that SMBHs in the nuclei of bright central galaxies (BCGs) are the drivers of radio bubbles that heat the intracluster matter (ICM) of cool-core clusters and groups (see McNamara & Nulsen 2007;Cattaneo et al. 2009;Fabian 2012;Gaspari et al. 2020).This is the main observational evidence so far in favor of the BH necessity scenario.However, BCGs are rare and peculiar galaxies, and their properties cannot be extrapolated straightforwardly to an enormously larger population of galaxies, for which BHs can still be considered freaks of nature.Indirect hints that the BH necessity scenario can be applied to the population of more common galaxies do exist and include a) the similarity of the star formation and BH accretion histories (Madau & Dikinson 2014;Aird et al. 2015;Fiore et al. 2017, hereafter F17), b) BH feedback is often invoked as the main driver of the observed correlations between BH and galaxy bulge properties (e.g.Gaspari et al. 2019;Marsden et al. 2020 and references therein), and c) BH feedback is also invoked as the cause of the low growth efficiency of high-mass galaxies (e.g., F17 and references therein).However, similar star formation and BH accretion histories can be the result of the BH coincidence scenario as well, and BH feedback is by no means the only or established main driver of the BH-bulge correlation and star formation quenching.In conclusion, despite very significant observational and theoretical efforts dedicated to the BH-galaxy coevolution problem in the last three decades, the answers to the root questions above remain disturbingly elusive so far, and BHs can still play the awkward role of the "stone-guest at the galaxy evolution supper".
The difficulty of converging toward a solution for the BHgalaxy coevolution problem probably arises because galaxies are complex systems that exchange matter, energy, and entropy with the environment through a network of interactions on many physical and temporal scales.In this respect, they are surprisingly reminiscent of other familiar complex systems: biological cells.This remarkable similarity prompted us to reconsider the problem of BH-galaxy coevolution from a different perspective: from the perspective of complex dynamical systems, which is one of the most successful approaches in evolutionary biology, neurology, and condensed matter.Although our evolving Universe is a clear example of the emergence of complexity from an initially high temperature and energy state, this approach has only been attempted sporadically so far, mainly because that it requires high-quality data and models, both of which have only recently begun to be available.
In complex dynamical systems, a rich dynamical behavior is generated from the interactions between a large number of subunits.These interactions generate properties that cannot be reduced to the behavior of the single subunit.Or, as in the quote by Phillip Anderson, "more is different".Complex systems share the following main characteristics: first, nonlinearity and sensitivity to initial conditions, the so-called butterfly effect.Second, Universal emergent properties (e.g., Laughlin 2005) that stem from the interaction of very many degrees of freedom and the naturally developing scale-invariance and power-law correlations.The same statistical theory can describe different systems, regardless of the system details.Third, self-organization is attained through feedback processes, which occur so that the microlevel interactions generate a pattern in the macrolevel that back-reacts onto the subunits.This is called coevolution in evolutionary biology, and it describes how organisms create their environment, which in turn influences the same organisms.This situation intriguingly resembles the BH-galaxy coevolution.Forth, the emergence of a singular behavior is reminiscent of a critical phenomenon, namely a second-order phase transition that is observed in systems that continuously transition between a more ordered and a more disordered phase by the tuning of a single parameter, usually the temperature.This transition occurs at specific values of the thermodynamic variables and exhibits features that are very different from more common first-order transitions that are observed in a more extended region in phase space.At the critical point, the range of spatiotemporal correlations diverges, leading to scale invariance and self-similarity.In recent years, it was observed that macroscopic systems can exhibit critical-like features without a tunable external parameter.In this case, long-range correlations emerge dynamically from the self-organization of local interactions.In this context, a purely Poisson dynamical system can reach a state without a characteristic scale as the critical state when a balance is reached from energy and mass inputs and outputs.It is more interesting when the criticality is self-organized, which is the result of correlated processes that are induced by feedback processes.Self-organized criticality (SOC) is a property found in extremely different physical system, from earthquakes to solar flares, magnetic materials, and the brain (e.g., Bak et al. 1987;Bak 1996;Beggs & Plenz 2003;de Arcangelis et al. 2006;Friedman et al. 2012;Fontenele et al. 2019).Dynamical systems exhibiting SOC show a highly stochastic behavior that fully reflects their complexity in a phenomenology that escapes any treatment with standard statistical tools.
We would like to reconsider micro-to macroscale observations and models of BH and galaxy growth by studying them from the perspective of complex dynamical systems.This paper is dedicated to the disk-wind system and therefore deals with microscales from several Schwarzschild radii r S1 (the innermost wind-launch radii) up to the outer edge of the accretion disk and the BH sphere of influence (hundreds of thousand r S ).We therefore searched in observations and models for three main characteristics of dynamical systems: 1) An absence of a characteristic scale and criticality, 2) self-organization through feedback processes, and 3) universality.We investigated whether these behaviors might be triggered by BH accretion and winds.Evidence of the same complex dynamics and universality class triggered by BH feedback from micro-to macroscales would highlight the existence of unified principles and point toward the BH necessity scenario.Conversely, if micro-to macroscale scaling relations are mainly produced by Poisson processes, this would point toward a BH coincidence scenario.
In standard geometrically thin optically thick Shakura & Sunyaev (1973) disks, deterministic steady-state solutions exist.The α viscosity that gave the name to these disk models includes both collisional and turbulent components and is the parameter governing the angular momentum transport.Although turbulence is an intrinsically chaotic process, an α disk retains steady-state solutions because of the key assumption that α is constant or varies slowly at most with respect to the other relevant disk timescales.A natural form of turbulence is produced by magnetic fields (e.g., Balbus 2003), and accretion disks are generally magnetized.Magnetohydrodynamic disk winds (MHDW) thus acquire particular importance in this context.MHDW models are nonlinear.In particular, the Grad-Shafranov equation is highly nonlinear.It expresses the poloidal projection of the momentum balance.The solution of this equation provides the detailed radial structure of the winds, but this is unfortunately an extremely complex equation without a general analytic solution.Simplified solutions based on self-similarity have been found by Blandford & Payne (1982) and Contopoulos & Lovelace (1994; also see Koenigl & Pudritz 2000 and Cui & Yuan 2020 for reviews).Several 2-dimensional and 3-dimensional numerical simulations of MHDWs have been published so far.In the AGN context, we recall the work of Fukumura et al. (2010) and that by Sadowski et al. (2016), Sadowski & Narayan (2016), Sadowski &Gaspari (2017, hereafter SG17), andGaspari &Sadowski (2017, hereafter GS17); who simulated a General Relativity MHDW system within a few 100r S .MHDWs are a natural context for a search for chaotic or complex system behaviors.Several studies have been carried out on theoretical (e.g., Winters et al. 2003) and observational grounds (e.g., Sukova et al. 2016).
Radiation-driven AGN winds (RDW) were discussed in a number of pioneering works during the 1990s (see, e.g., Murray et al. 1995;Proga et al. 2000;Elvis 2000 and references therein) to explain different observations in the framework of AGN atmosphere unification schemes.Hydrodynamical models were developed by various authors (see, e.g., Proga 2007;Kurosawa & Proga 2009;SG17;and references therein).Radiation driving is at the core of the wind feedback model developed by Silk & Rees (1998), Fabian (1999), King (2003) and Zubovas & King (2012, see King & Pounds 2015 for a review).We discuss both MHDW and RDW in the following sections.
A key process in SMBH feeding is chaotic cold accretion (CCA; Gaspari et al. 2013): The turbulent hot atmosphere recurrently condenses into cold clumps that feed the SMBH from the macro-to the microscale through recurrent/fractal inelastic collisions.These collisions create a self-similar time variability with a flicker-noise 1/ f power spectrum (Gaspari et al. 2017).CCA is thus a key example of a chaotic and emergent process in feeding/feedback astrophysics, which is important for understanding some of the results and assumptions in this work.Furthermore, Mineshige & Negoro (1999) tried to interpret the apparently chaotic X-ray variability of Galactic BH candidates in an SOC framework using a simple cellular automaton to describe the accretion toward the BH.
Winds are a necessary element for dissipating angular momentum and through this, for allowing efficient BH gas accretion.Nuclear (launching radii 10-100r S ) semirelativistic winds are frequently observed in AGN (e.g., Tombesi et al. 2010Tombesi et al. , 2011Tombesi et al. , 2012Tombesi et al. , 2013;;Serafinelli et al. 2019;Chartas et al. 2021; see the full reference list in Table A.1).These winds are likely magnetically driven in many cases because a radiation drive can only be efficient for sources that emit at or above their Eddington luminosity (Reynolds 2012).At larger scales (from a fraction of a parsec to a few parsec, i.e., from thousands to tens of thousand r S ; this is also known as the mesoscale; see GS17), VLT and ALMA interferometric observations have unveiled nuclear dust/molecular disks, rings, and tori in local Seyfert galaxies (Jaffe et al. 2004;Gallimore et al. 2004Gallimore et al. , 2016;;Garcia-Burrillo et al. 2016;Combes et al. 2019;GRAVITY Collaboration 2020a,b, 2021).In the few objects in which jets and nuclear disks/rings were observed simultaneously, these structures are not aligned.Because jets are thought to be perpendicular to the inner accretion disks, this suggests that the wide angle relativistic winds emerging from the inner accretion regions can impact the structures around the BH sphere of influence at least in some cases (e.g., tori, rings, or clouds).This can give rise to some feedback between disk winds and gas accretion, which we explore below.
Key open questions still remain unsolved: whether micromesoscale feedback is relevant and whether these feedback and nonlinear processes are clues for a complex dynamical system in action.To answer these questions, we first collected observations of semirelativistic winds (or ultrafast outflows, UFOs; Sect.2) because they simultaneously are the most common and most energetic nuclear outflows.We then interpret them in the framework of MHDW and radiation driven wind scenarios in Sects.3 and 4. We study the statistical properties of UFOs and derive the distribution functions of the ratio ω of the mass-outflow and inflow rates in Sect.5, where we also study the links between ω and the Eddington ratio λ = L bol L Edd .We then introduce in Sect.6 a simple cellular automaton to investigate how the dynamical properties of an idealized disk-wind system change following the introduction of simple feedback rules.We show that if feedback is present, the system can be driven toward an SOC state.We finally discuss our results in Sect.7 and present our conclusions in Sect.8.

Sample of observed ultrafast outflows
We compiled UFO observations from the literature that were fit with an homogeneous model (the photoionization model XSTAR within XSPEC).This model provides the gas outflow velocity v w , the column density N H , and the ionization parameter ξ.Most of the observations were taken from Tombesi et al. (2010Tombesi et al. ( , 2011Tombesi et al. ( , 2012Tombesi et al. ( , 2013) ) and Gofford et al. (2013Gofford et al. ( , 2015)), and they were furthermore incremented with results on single sources.The final sample includes 54 system observations in 40 AGN (see the appendix for details).
To calculate the mass-outflow rate Ṁw , we used the approach of Krongold et al. (2007), and in particular, their formula in Appendix A2 for the mass-outflow rate of a disk wind in a conical geometry, where f is a function of the angle of the line of sight to the central source and the accretion disk plane and the angle formed by the wind with the accretion disk.r is the distance of the illuminated face of the wind from the central luminosity source.The numerical estimate on the right-hand side applies when f = 1.5, corresponding to a roughly vertical disk wind and an average line-of-sight angle of 30 • .
To account for the scaling of the mass-outflow rate with source activity, we studied the mass-outflow rate normalized by the mass-accretion rate ω = Ṁw Ṁaccr , where Ṁaccr = L bol c 2 , and is the radiative efficiency.We show in the next section that ω has a profound significance in the context of accretion models.
While N H and v w in Eq. ( 1) are results of fitting the data with a photoionization model, r cannot be directly estimated from the X-ray data in most cases.It is usually assumed in the literature either as the radius at which the observed velocity equals the escape velocity, r min = 2GM BH v 2 w , or from the definition of the ionization parameter ξ = L ion n e r 2 , where L ion is the luminosity above the hydrogen ionization threshold.It follows from this that r ion ≈ L ion 1.23ξN H ∆r r because n e ≈ 1.23n H for a fully ionized gas with solar abundances, and N H = n H ∆r, where the thickness of the absorber is ∆r r.If r = r min , Ṁw (r min ) = 0.88 N H 10 24 and ω(r min ) = 5.03 N H 10 24 In this formulation, Ṁw (r min ) depends linearly on N H and M BH and inversely on v w , while ω(r min ) depends linearly on N H and inversely on v w and λ.The range covered by the measured N H is between 10 22 cm −2 and 10 24 cm −2 , that is, 2 dex, with a typical uncertainty of about 0.2-0.5 dex.The range covered by v w is 0.05-0.5c,about 1 dex, with a typical uncertainty of 0.1-0.2dex.The range of λ = L bol L Edd spans 2 dex (see the appendix for details).The uncertainty associated with the unknown geometry of the wind is about 0.2-0.3dex (Krongold et al. 2007).While the relative statistical and systematic uncertainties on Ṁw (r min ) and ω(r min ) are generally small, Ṁw (r min ) can by construction be regarded as a lower limit to the real mass-outflow rate.
If r = r ion , Ṁw (r ion ) = 0.24 L ion 10 44 and ω(r ion ) = 1.37 assuming that L ion ≈ 0.1L bol (see, e.g., Lusso et al. 2012).In this formulation, Ṁw (r ion ) depends linearly on both L ion and v w and inversely on ξ.It does not depend on N H . ω(r ion ) depends linearly on v w and inversely on ξ.The range covered by ξ is about 3 dex, and its uncertainty is about 0.1-0.2dex.∆r/r is generally difficult to estimate from the data.The few attempts have been made for warm absorbers rather than UFOs (see, e.g., Krongold et al. 2007).∆r/r is thought to be about 0.1 or lower in conical outflows.Critically, different flows can have quite different ∆r/r, and therefore, the scatter in ∆r/r can dominate the scatter on ω(r ion ).Therefore, assuming a fixed value may artificially reduce the intrinsic scatter of ω(r ion ) in a sample of flows.In the end, the real normalized mass-outflow rate should be somewhere in between ω(r min ) and ω(r ion ), and the scatter on ω(r min ) in a sample of flows is probably a lower limit to the real scatter.In our sample of 54 flows, the median values of r min and r ion are about 40r S and 300 ∆r/r 0.1 r S , and the average of the difference between ω(r min ) and ω(r ion ) is 0.8 dex.We plot in Fig. 1 the average between these two extremes, computed as log( ω) = log( ω(r ion )) − log( ω(r min )) = log( ω(r min )) + 0.4, as a function of the AGN luminosity in Eddington units λ.Ṁaccr was estimated from the bolometric luminosity as Ṁaccr = L bol c 2 assuming a radiative efficiency = 0.1.The wind mass-outflow rates in Fig. 1 were corrected for relativistic effects following Luminari et al. (2020).The corrections were usually small, however, and the results obtained with the uncorrected data are quite similar to the results presented below.
Two features stand out in Fig. 1.The first feature is the expected correlation between ω and λ (see Eq. ( 6)).The Spearman-rank correlation coefficient of the ω − λ correlation is −0.45 for 52 degrees of freedom.The best-fit slope of the logarithmic ω − λ correlation in Fig. 1 that is, ω scales with the ratio of the wind Eddington rate and the Eddington rate.The second feature is that the observed ω spans more than 3 dex (more than 2 dex at all λ).This range is wider than the statistical or systematic uncertainty on Ṁw , and it is also robust against the presence of outliers.Some of the systems with the highest ω in Fig. 1 are lensed systems, in which an usually large uncertainty on the lens magnification adds to the other systematic and statistical errors.Other high ω systems are lowluminosity highly obscured Seyfert 2 galaxies, where the photoelectric cutoff extends to the energies of the UFO absorption features, which means that they are more complex to characterize.Even when these problematic sources are excluded, the picture remains substantially unchanged.As discussed above, the scatter in Fig. 1 is most likely a lower limit to the real scatter.In the following section, we compare these findings with the prediction of disk-wind models.

Magnetohydrodynamic disk-wind scenario
Nuclear winds are intimately related to mass accretion through the conservation of angular momentum, and they are likely magnetically driven in most cases (e.g., Reynolds 2012).More recently, Wang et al. (2022) performed MHD simulations of winds from thin accretion disks.The authors suggested that events resembling observed UFOs can certainly be generated within about 600r S .Beyond this radius, winds have a smaller ionization parameter than observed UFOs.In our sample of 54 UFOs, all but two have r min < 600r S and 60% have r ion < 600r S .We also note that the ionization of the wind is governed by how radiative feedback is produced/injected in the wind (in addition to ∆r/r).The source geometry of the ionization radiation is definitely poorly constrained.It can be an extended corona (which produces high ionization at relatively large radii) or a lamp-post region that is confined in the innermost region of the accretion disk (which mainly produces ionization at small radii).
The MHDWs are therefore a plausible way to describe the observed UFOs.We do not discuss the fate of these winds and whether they will produce energy-driven or momentum-driven outflows in the galaxy.
In MHDW systems, the conservation of momentum rate reads (see the reviews of Koenigl &Pudritz 2000 andPudritz et al. 2007) where Jd and Jw are the disk and wind momentum rates, respectively, Ṁin and Ṁout are the mass-accretion and -outflow rates, r 0 and r A are the wind launch radius and the Alfvén radius, and Ω 0 is the Keplerian angular velocity at the launch radius (every annulus of the disk rotates close to its breakup speed).Mass accretion and mass outflows are thus intimately connected.
No accretion is possible without the launch of massive winds from the system, and the Alfvén radius is the lever arm of the outflow.Because r A r 0 > 1, relatively small mass outflows can cause relatively large inflows: The removal of angular momentum by the winds allows matter to accrete inward.The Alfvén radius marks the extent of the region in which magnetic stresses dominate the flow.this follows the poloidal magnetic field lines anchored to the disk.Therefore, beyond this radius, the outflow terminates and corotates with the disk, reaching a termi- , where v 0 is the wind launching velocity at r 0 .The acceleration from the disk occurs through a A36, page 4 of 19 centrifugal effect, and therefore, these models are called magneto-centrifugal disk wind (MCDW).
In all disk-wind models, the accretion and outflow rates Ṁin and Ṁout depend on the radius.The outflow rate probed by an X-ray UFO observation probably concerns the phase when the wind emerges from the system at tens or hundreds of Schwarzschild radii r S , well before the entrainment phase, in which the outflow rate can be significantly loaded (see, e.g., Serafinelli et al. 2019), and therefore, it is quite justifyable to assume Ṁw = Ṁout .The observed X-ray and bolometric luminosities are produced in regions from a few to some dozen r S , and therefore, we identify Ṁin with Ṁaccr .Therefore, Ṁout Ṁin ≈ ω.In MCDW models, ω = ( r 0 r A ) 2 ∈ [0.01, 1] for typical r 0 r A ratios, and r A r 0 .The yellow area in Fig. 1 marks the expectations of the MCDW model.ω depends on the intensity of the poloidal magnetic field at the radius r 0 : The energy density of the magnetic field should equal the kinetic energy density in the wind, or where V A is the poloidal Alfvén velocity, v p is the poloidal velocity , B p is the poloidal magnetic field, and ρ is the wind density.From Eq. ( 8), it follows that While it is useful to emphasize the key role of winds in accretion systems in a direct and simple way, the MCDW models described above are highly idealized and are based on a simplified analysis.They do not include important physics such as thermal and radiation pressure and general relativity effects.Sadowski et al. (2016), Sadowski & Narayan (2016), SG17, and GS17 developed more realistic MHDGR models and simulated the accretion/outflows in a box of ∼100r S for systems with a wide range of λ.They found a generally high kinetic efficiency w = L w Ṁ• c 2 ≈ 0.02−0.04 that was roughly constant with λ from highly sub-Eddington to super-Eddington regimes.
Here, Ṁ• is the accretion rate through the event horizon, and it thus contributes to the BH growth.Due to GR effects, Ṁ• can be lower than Ṁaccr (the accretion rate at larger radii, where most of the radiative luminosity is produced), and thus, w = The shaded green area in Fig. 1 marks the expectation of this model for v w /c ∈ [0.05, 0.5], the velocity range of most observed UFOs. Figure 1 shows that all observed UFO ωs are consistent with the predictions of the MHDGR model considering the statistical and systematic uncertainties.

Radiation-driven disk-wind scenario
Winds can be efficiently driven by radiation when the luminosity of the central source is sufficiently high and when the gas is sufficiently optically thick, as illustrated by the Eddington equation, where M BH is the black hole mass, m p is the mass of the proton, and σ is an absorption cross section.For completely ionized gas and no dust, σ = σ T , where σ T is the Thompson cross section, while when the gas is not completely ionized and/or dusty σ σ T , because the line and dust opacities can be orders of magnitude higher than the Thomson-scattering opacity.However, in the innermost nuclear regions, below tens or hundreds of gravitational radii, the gas is likely fully ionized and dust free, and therefore, an efficient wind-driving through radiation pressure requires huge luminosities and dense gas clouds.When the gas is sufficiently dense (Compton thick, N H 10 24 cm −2 ), each photon would scatter at least once before it escapes, yielding all its moment.Therefore, Ṁw The identification of UFOs with these continuous winds is problematic, however, because Compton-thick intervening gas would strongly depress the nuclear X-ray flux, making the detection of absorption lines difficult if not impossible.A solution for this apparent conundrum is that UFOs are episodic, as in the scenario suggested by King & Pounds (2015).The wind is launched when the central engine runs at L bol ∼ L Edd .When this phase is short, the wind is not continuous, a shell moves away from the nucleus, and by expanding, it gradually decreases its column density.At some point, absorption lines become observable.In this scenario, the wind mass-outflow rate at the launch radius is (King & Pounds 2015).Then, Ṁw should be nearly constant with distance from the nucleus and time because Ṁw ∝ N H × r × v w (see Eq. ( 1)).Finally, the scatter of the observed ω in this scenario depends on η and v w /c, as in the case of continuous winds.
The fast wind expands until it collides with the ambient gas, forming a reverse shock that slows the nuclear wind down, and a forward shock expands in the interstellar gas of the BH host galaxy.The AGN radiation can efficiently cool the shocked gas, deciding the nature and the fate of the shock.When the cooling is strong, the shocked gas looses most of its energy and only transmits its ram pressure Ṁv w to the post-shock gas (momentum-driven flows).Conversely, when cooling is negligible, the flow maintains most of its energy and expands adiabatically in the interstellar gas (energy-driven flows).Since cooling is mostly driven by the AGN radiation field, the cooling time near the BH will be shorter than the flow time and the flow will be driven by momentum.The cooling time increases as R 2 , but the flow timescale increases linearly with distance.This means that beyond a critical distance, the cooling time can exceed the flow time, and the flow becomes energy driven (Zubovas & King 2012;King & Pounds 2015).UFOs probe nuclear flows, and they therefore most likely belong to the momentum-driven regime.It is thus expected that most of the sources in compilations of UFOs and galaxy-scale outflows are compatible with momentum-driven flows (e.g., Bonanomi et al. 2023).Because nuclear wind is likely short lived, the observed UFOs may well not be causally related with the more long-lived galaxy-scale winds.

Population properties versus single-system evolution
The main feature highlighted by Fig. 1 is the broad range of ω that is covered by existing observations.The range is broader than the statistical errors on single determinations and the systematic errors that are certainly present in a biased and inhomogeneous sample such as we used in this analysis.Therefore, the wide ω range is likely intrinsic to the wind systems, suggesting that a broad distribution of this key parameter is produced by physical disk-wind dynamic systems.As discussed in the A36, page 5 of 19 previous sections, a broad range of ω can be easily explained by either the MHDW or the RDW scenarions, in which different physical drivers cause the scatter.To further investigate this crucial point, we calculated the cumulative and differential distribution functions of the number of systems in which ω exceeds a limit value as a function of this limit value.These functions are the equivalent of the classic LogN-LogS functions in astrophysics.The distribution functions are plotted in Fig. 2 for the full sample (left panel) and for the sources at z < 0.3, that is, high-redshift luminous quasars and lensed sources are excluded.Luminous quasars can be powered by accretion flows quite different than low-redshift AGN, while the uncertainty on the wind and accretion parameters in the lensed sources are increased because it is difficult to model the amplification at different spatial scales.From now on, we consider the sample of AGN at z < 0.3, which is both more homogeneous and robust.
Below ω ∼ 0.05, the cumulative distribution function in the right panel of Fig. 2 flattens considerably, most probably because selection effects limit the completeness of the UFO sample.Therefore, to be representative of the global AGN population, the functions in Fig. 2 must be corrected for these observational selection effects.Winds with a low mass-outflow rate can only be detected in spectra with the highest signal-to-noise ratio, which are rarer.This produces a flattening of the functions toward low ω values.A correction for this selection effect would increase the number of systems with small ω.We evaluated a correction for incompleteness as described in the appendix.Figure 3 shows the normalized cumulative and differential (per decade) distribution functions of ω for the sample of UFO systems at z < 0.3, corrected for incompleteness.The distribution functions were normalized by dividing by the total number of sources (34) in which 43 independent UFOs were observed.
The power-law slopes of the ω cumulative and differential (per decade) distribution functions above ω = 0.2 are −0.86 −0.40 +0.27 and −0.57−0.39 +0.29 , respectively (90% confidence error).The powerlaw slopes of the λ w cumulative and differential distribution functions above λ w = 0.1 are −1.04 −0.61 +0.42 and −1.00 −0.76 +0.44 , respectively.In the ω and λ w ranges 0.007-0.14and 0.004-0.04,the differential distributions are flat, with a slope consistent with 0. The application in the next years of this type of analysis to less biased and complete AGN samples, such as the SUBWAYS sample (Mehdipour et al. 2022;Matzeu et al. 2023) is crucial to reduce the statistical and systematic uncertainties.
When the wind-activity timescale is shorter than the AGN timescale, as is likely in many cases (e.g., Nardini & Zubovas 2018), and when the population functions in Fig. 3 are drawn from independent snapshots, we can argue that the wind-activity history in each single source can be statistically described by these population functions.If this is the case, the message suggested by the analysis in the previous sections is that UFOs can be seen as series of sporadically launched quasi-spherical shells, with several or many relatively small wind events for each main wind event in each AGN activity cycle.This scenario is consistent with both MHDW and RDW.Expanding on this idea, we suggest that wind events may have a fractal behavior, analogous to what was suggested by Gaspari et al. (2017) for CCA.Since λ should follow a power-law distribution down to an inflection point as well.λ w and ω are both independent of the BH mass.Thus, the distribution of λ can be obtained from the distribution functions of ω and λ w , both corrected for selection effects.We therefore assumed broken power-law distributions for ω and λ w (with the parameters given above), and randomly generated 10 000 ω and λ w values following these distributions, mimicking 10 000 wind and accretion events in the active lifetime of a single source.The λ distribution function obtained in this way is plotted in Fig. 4 (red histogram) and resembles a log-normal distribution that peaks at λ ≈ 0.5.To compare this distribution function with the observed λ distribution functions from AGN surveys, again assuming that the population properties reflect the evolution of single systems, we need to assess the time spent by the system on each λ value.To estimate this timescale, we considered that the relativistic winds carry momentum and energy and that they are likely at a wide angle (e.g., Nardini et al. 2015;Luminari et al. 2018).The wind propagates in the direction of least contrast, that is, nearly perpendicularly to the accretion disk, until it encounters matter with a density that is high enough to shock (e.g., Menci et    references therein).This is likely to occur at the limit of the BH sphere of influence.Beyond this region, the geometry of the matter is not governed by the BH gravity and is more complex, and it is in many case not aligned with the accretion disk.VLT interferometric observations and ALMA observations have unveiled nuclear dust/molecular disks/rings/tori on scales down to a fraction of one parsec in local Seyfert galaxies (Jaffe et al. 2004;Garcia-Burrillo et al. 2016;Gallimore et al. 2016;Combes et al. 2019;GRAVITY Collaboration 2020a,b, 2021), that is, a scale close to the BH sphere of influence.This can be parameterized as the BH Bondi radius , where c ∞ is the sound speed well beyond R B , and the latter equality gives R B in units of the Schwarzschild radius r S .c ∞ = ( γkT µm p ) 1/2 ≈ 400 km s −1 when the Compton temperature of the gas irradiated by the AGN is T C = 2 × 10 7 K, while µ is the mean atomic mass per proton, which we fixed to 1.2 (see, e.g., Gofford et al. 2015).Thus, R B = 2.8 × 10 5 r S , that is, 2.7 pc for M BH = 10 8 M .Therefore, it is plausible to assume that the relativistic wind propagates nearly freely until it impacts the matter around the BH sphere of influence, where the wind thrust may efficiently contrast gas accretion within the Bondi radius.If this is the case, where P b is the thermal pressure in the shocked-wind bubble, and f ( ω) is the fraction of the solid angle covered by the accreting matter that is intercepted by the wind.The thermal pressure in the wind bubble is related to the thermal energy injected by the wind, E b = 3/2P b V b , where V b = 4/3πr 3 is the bubble volume, and it can be expressed as E b ≈ 1/2 L w ∆t (Weaver 1977).Following the standard Eddington argument, the gas mass that can be just supported by the wind thrust is the BH massaccretion rate times a characteristic activity time ∆t, that is, M gas = Ṁaccr × ∆t.Therefore, Eq. ( 12) can be rewritten as follows: Ṁw can be removed from the equation, which can be solved for ∆t, When the radius at which the wind impacts with the accreting matter is the Bondi radius, and when we recall that in MHDW A36, page 7 of 19 systems, v w = √ 2v 0 √ ω , the equation reads The second term of the denominator is negligible unless f ( ω) 1, thus ∆t ≈ 140 yr × √ ω( v 0 /c 0.1 ) −1 if M BH = 10 8 M and R B = 8.3 × 10 18 cm.This is similar to the timescale derived by King & Pounds (2015), which again suggests that many wind/accretion episodes can easily occur during a typical AGN lifetime of several millions years.
∆t depends on √ ω.We recall that λ = λ w / ω and that ω < 1.This therefore suggests that high λ states are shorter than low λ states.When a source spends a shorter time in high λ state than in lower λ states, the probability of finding this source in the former state is lower than the probability of finding it in a state of low λ.Consequently, the "observed" λ distribution in a blind survey is skewed toward low λ values (the most likely states).Thus, the observed λ distribution should be compared with the expected distribution, weighted by the different times the sources remain in any given λ states, that is, λ × ∆t ∝ λ w √ ω .The black curve in Fig. 4 represents the distribution of λ × ∆t = λ w √ ω .The slope of the λ × ∆t distribution above 0.05 after correcting current UFO measurements for selection biases is −0.9 ± 0.5, which is consistent with both the Bongiorno et al. (2012) determination, −0.96 in the redshift range 0.3-0.8 and Aird et al. (2012Aird et al. ( , 2018) ) determination, −0.65.Ananna et al. (2022) estimated the λ Edd distribution in a sample of local AGN and reported a broken power law with a break λ Edd similar to our λ w threshold, but with a a slope above the break of ≥−2, which is steeper than that of the λ × ∆t distribution and that of previous determinations.Kauffmann & Heckman (2009) reported log-normal distribution functions due to the presence of low-luminosity AGNs in their samples.
We finally computed the total BH growth obtained by a large number (10 000-100 000) of accretion events with λ following the distribution function in Fig. 4 (the black histogram shows the constant time spent in each accretion episode).This is equivalent to the BH growth obtained at a constant λ ∼ 0.1 during the entire AGN activity time.This number is very sensitive to details of the λ w and Ω distributions: their minimum value, the inflection points, and the distribution slopes.Therefore, the above figure should not be considered a firm limit for real accreting systems.We limit ourselves to noting that if the three main assumptions of the above analysis are correct (1. a link between accretion and wind emission, i.e., Ṁaccr = 1/ ω Ṁw ; 2. an Eddington-like limit on the mass-accretion rate because the wind thrust contrasts the accretion flow at the Bondi radius; and 3. a random distribution of accretion/wind events), then systems can hardly accrete close to their Eddington rate for the full AGN time cycle.High accretion rate events, even when they exceed several Eddington units, are likely, but their number will be small in comparison to low accretion rate episodes and their duration will be shorter.

Is the disk-wind system a self-organized critical system?
Self-organized criticality (SOC) was first introduced by Bak et al. (1987,1996).It is a characteristic of some com-plex dynamical system.These systems reach the critical state not by tuning an external parameter, but because of the selforganization of internal interactions, among which feedback is the most important.SOC systems develop power-law scaling laws and are therefore scale invariant.Because of these intriguing similarities with disk-wind systems (natural development of feedback and power-law scaling laws), it is reasonable to wonder whether disk-wind systems are SOC systems.The first pioneering work on SOC in BH astrophysics (Mineshige et al. 1994a,b;Mineshige & Negoro 1999) followed the introduction of SOC systems in the literature by just a few years.These works were focused on interpreting the chaotic X-ray variability of stellar BH and AGN in the framework of SOC models.Mineshige and collaborators were successful in generating the typical 1/ flike fluctuations that are often observed in these sources using a simple cellular automaton based on simple SOC rules.Their conclusions were later questioned by Uttley et al. (2005) and Done et al. (2007) on the basis that SOC models naturally produce power-law distributions of the temporal fluctuations, while log-normal distributions are typically observed in stellar BH systems.In AGN, the situation is not clear cut because in the beststudied case, that of the narrow-line Seyfert 1 IRAS 13224-3809 (Alston et al. 2019), there are strong signs of nonstationarity and strong departures from a log-normal distribution of fluxes.
Even the presence of power-law distributions is not sufficient to assess that the system can be described by the theory of critical phenomena, however, because many mechanisms can be at work to generate a power-law distribution.In addition to the power spectrum of the time series of the events, SOC systems must show the typical relation between some scaling exponents (Kuntz & Sethna 2000;Sethna et al. 2001).At least three other fundamental scaling laws must be analyzed: the histograms of the size S and duration T of the burst f (S ) ∼ S −τ , f (T ) ∼ T −τ t , and the scaling of the average size of avalanches with a given duration S (T ) ∼ T α 1 .The parameters τ, τ t , and α 1 are critical exponents of the system, regardless of the microscopic details of the system.Scaling theory (Kuntz & Sethna 2000;Sethna et al. 2001) predicts exponent relations, and in particular, when the system is at criticality, α = α 1 .All critical exponents can be measured both in data streams and in models, and the above relation allows us to assess whether the system is consistent with an SOC state.We first reproduced the Mineshige et al. (1994a,b) cellular automaton according to the following steps: 1.We divided the disk plane into cells along the radial and azimuthal coordinates.We chose N φ = 128 azimuthal segments and N r = 64 radial segments (see the discussion below for a justification of this choice).2. We randomly chose one cell in the outermost ring and placed an accreting gas particle of mass m at each time step.M i, j is the mass of a cell at radius r i and azimuthal coordinate φ j .3. We verified whether the mass accumulated at any location exceeded a given threshold, that is, whether M i, j > M crit .
The results do not depend on the actual value of this threshold.We started the simulation by assigning a mass to all disk cells.We verified that the results did not change for different starting conditions, for example, adopting the condition of an initial empty disk, adopting a random distribution of cell masses, or starting by assigning a mass to all disk cells of M crit .Every time M i, j > M crit , three mass particles fell from the cell into adjacent cells at r i+1 , that is, an accretion event A36, page 8 of 19 occurred in the disk, Given the circular boundary condition in 1), if j = 128 j+1 = 1, and if j = 1 j − 1 = 128 in the above expressions.This accretion of mass may cause that the mass in some of the inner r i+1 cells might also exceed the critical threshold, and an avalanche thus propagates inward at the following time steps.We call an avalanche the propagation of an accretion event through the disk.4. We repeated the procedure n times to obtain a statistical population from which temporal and spatial distributions can be drawn.Mineshige & Negoro (1999) introduced an additional role to account for inward diffusion, in order to better match observed variability light curves.We did not consider diffusion here because we are not interested in reproducing the fast disk variability, but rather the long-term disk evolution.
In order to describe the statistical properties of the system, two main quantities were calculated for each avalanche: (a) the lifetime of the avalanche T , that is, the number of time steps of single avalanche propagation, and (b) the size of the avalanche S , that is, the total number of accretion episodes during the avalanche evolution.At each time step, several avalanches can be active on the disk.We finally followed Mineshige et al. (1994a,b) in calculating the total energy that is liberated at each time step (adding the contributions from each active avalanche).
For our setup, the maximum avalanche size was limited to 4096 ( 64 i (i * 2) − 1).This largest avalanche covered all 128 cells in the innermost ring.We therefore set the number of azimuthal elements (N r × 2).The original setup of Mineshige et al. (1994a,b), Mineshige & Negoro (1999) used N φ = N r = 64, implying that avalanches propagating inside r = 32 spread multiple masses on the same cell.This additional accumulation of mass on some cells with respect to the simple rule n. 3 above was counterbalanced by assuming an M crit proportional to r.Using the setup in step (1) above, we dropped this additional complexity and used a constant M crit .We were able to reproduce the results of Mineshige et al. (1994a,b) using both our simple setup and their original setup with M crit proportional to r.We verified that the results did not depend on the actual value of M crit provided that M crit ≥ 3m.
We ran this simple cellular automaton to build the distributions f (T ) and f (S ).To reduce the noise due to the limited number of avalanches with large sizes, which hinders an assessment of the real shape of the distribution f (S ) for high S values, we ran the cellular automaton to produce >10 avalanches of size S > 1000.We produced a total of 210 000 avalanches and discarded the first 10 000 to remove possible initial transients.
The upper panels in Fig. 5 show the distributions of f (T ), f (S ), and the average size for a given duration S (T ).Table 1 lists the slope τ t of the f (T ) distribution, the slope τ of the f (S ) distribution, the slope α 1 of the S versus T , and α = (τ t −1) (τ−1) .τ t and τ were obtained by fitting the f (T ) and f (S ) distributions with a power law plus an exponential-cutoff model.This cutoff is due to the limited size of the system and moved to higher S when systems with N r > 64 and N φ > 128 were considered.The point at T = 65 in the f (T ) distribution represents the avalanches exiting the system from the inner radius.The best-fit τ t and τ are close to the values of −3/2 and −4/3 that are expected for variants in the original Bak et al. (1987) model, introducing a preferential direction (Ben Hur et al. 1996).The errors on τ t and τ are dominated by systematic uncertainties, that is, the wiggle in the f (T ) distributions and the excess at high S in the f (S ) distribution for setup 1.These wiggles suggest that the system is supercritical and therefore out of SOC.This is further supported by the mismatch between α = (τ t −1) (τ−1) and α 1 (see Table 1).We recall that setup 1 is conceptually similar the original sand-pile model of Bak et al. (1987), with the introduction of a preferential direction.In the experiment of Bak et al. (1987), SOC is achieved through self-regulation between the mass entering the system and the mass exiting the system boundaries.In our implementation, the possibilities for the mass to exit the system are severely reduced with respect to the original experiment because the only egress is via the central compact object.This produces an excess in the number of avalanches close to the maximum size.
Setup 1 does not include feedback processes, and we considered it as a benchmark against which to compare results arising from more complex scenarios.To gain further insights, we added three additional ingredients to the baseline cellular automaton: a link between gas accretion and wind emission, as in MCDW models (setup 2), an Eddington-like limit on the mass-accretion rate due to the wind thrust that contrasts the accretion flow at the outer radius (setup 3), and a modification of setup 1 to account for the Eddington limit to the gas accretion due to radiation pressure (setup 4) as in RDW.
In setup 2, at each location of the disk where M i, j > M crit , we randomly extracted a poloidal Alfvén velocity in units of the Keplerian velocity at the magnetic wind footpoint v K0 from 0.01 to 100.This is a proxy for the poloidal magnetic field.We then used the relations in Fig. 3 of Cui & Yuan (2020) to estimate r A /r 0 and thus ω = (r A /r 0 ) −2 .The fiducial model of Cui & Yuan (2020) is a hot-accretion flow, but the authors discussed deviations with respect to this baseline.They discussed that winds emerging from thin disks may have different properties.In particular, thin disks are cooler and thus have lower sound speeds C s0 than thick disks.For this reason, the authors calculated models not only for C s0 = 0.5v K0 , as appropriate for hot disks, but also for C s0 = 0.1v K0 , as appropriate for cold disks.We used both models from Fig. 3 of Cui & Yuan (2020) in our cellular automaton, and the differences were negligible.
In setup 2, the mass transferred to M i+1, j , M i+1, j+1 and M i+1, j−1 was set to m in = m/(1 + ω), while a mass m out = m ω/(1+ ω) left the system in the form of a wind.The system has a preferential direction, similar to setup 1, but mass can leave the system at each radius, so that not all masses that accrete at the outer radius are forced to reach the inner radius.This is different from both the original Bak et al. (1987) setup and the directional setup of Ben Hur et al. (1996).τ t in Table 1, is −1.63, steeper than the value expected for directional avalanches.This indicates a different universality class.We find that setup (2) approaches criticality because α 1 ∼ α (see Table 1 and Fig. 5).However, it must be noted that the scaling of S versus T does not follow a pure power-law shape.The best-fit slope for T < 10 is 1.49, but it is 1.82 for T > 10.
Setup 3 is based on setup 2 and further includes an Eddington-like limit on the mass-accretion rate due to the wind thrust that contrasts the accretion flow at the outer radius.At each time step, the mass that accretes on the outer radius is reduced by a factor ω , that is, m in = m × (1 − ω ) at this location.If ω is small, m in = m, while if ω approaches one, there is A36, page 9 of 19 Fig. 5. Left panels: distribution of the avalanche duration f (T ).Central panels: distribution of the avalanche size f (S ).The gray points are the original distribution function, and the black points mark the distribution logarithmically binned.Right panels: average avalanche size S vs. avalanche duration T .Upper panels: setup 1, similar to the original simulation of Mineshige et al. (1994a,b).Central upper panels: setup 2, including a link between mass accretion and wind mass outflow as in MCDW models.Central lower panels: setup 3, including a link between mass accretion and wind mass outflow as in MCDW models and an Eddington-like limit on the mass-accretion rate due to the wind thrust that contrasts the accretion flow at the first radius.Bottom panels: setup 4.This is the same as setup 1, but includes a threshold on the disk luminosity above which the accretion is halted (similar to an Eddington limit).The red lines in the left and central panels are best-fit models.The red lines in the right panels bridge the interval of allowed α.

A36, page 10 of 19
Table 1.Avalanche exponents for the different setups.
negligible accretion at the outer radius.This is a rather rough feedback recipe, and it should not be regarded as a realistic diskwind model.We studied this setup to understand how and in which direction the introduction of this further feedback process modifies the dynamical status of the system.The inclusion of this further feedback brings the system even closer to SOC.In this setup, the scaling of S versus T does follow a pure power law.
Finally, we considered a different modification of setup 1 by introducing a threshold in disk luminosity above which the accretion is abruptly truncated.The disk luminosity was calculated following Mineshige et al. (1994a,b) and assuming that the kth gas particle falling from the ring i to the ring j liberates a potential energy proportional to 1 If L disk > L threshold , the gas particle is pushed outside the system by radiation pressure and the accretion is blocked.This causes the luminosity to decrease, resuming the accretion and starting a new cycle.The inclusion of this feedback regularizes the flow, and the system is much more closed to criticality than for setup 1 because α 1 ∼ α (see Table 1 and Fig. 5).However, as in setup 2, the scaling of S versus T does not follow a pure power law in this case either.The best-fit slope for T < 10 is 1.27, but it is 1.98 for T > 10.
Another test to assess whether a system is close to SOC is the analysis of the shape of the avalanche profiles (e.g., Nandi et al. 2022 and references therein).At criticality, the avalanche evolution with time should follow a universal shape (a parabola for Barkhausen noise) for all durations.More precisely, the shapes corresponding to different avalanche durations should collapse on top of each other when the average size at each given time t is normalized by T α coll −1 .To visualize the situation, we show in Fig. 6 three avalanches for each setup.The left panels show an avalanche that did not reach the innermost radius, the central panel shows a complete avalanche, and the right panel shows the (complete) avalanche with the largest size for that setup.The avalanche with the largest size for setup 1 fully covers the innermost radius, that is, it does not decrease in size during the evolution, but continues to expand down to the innermost radius.The exponent α coll is an independent estimate of the critical avalanche exponents and corresponds to the best collapse, that is, the collapse with the least scatter between the points for different avalanche durations.Figure 7 shows S (t, T ) /T α coll −1 for avalanches with durations of T = 20, T = 30, T = 40, T = 50, and T = 60 as a function of t/T .In setup 1 (upper left panel), the avalanche evolution does not follow a universal shape, but exhibits a pronounced increase in size toward the end of the avalanche for events with a long duration.In setups 2 and 3, the avalanche profile approaches a parabolic shape (Fig. 7, upper right and lower left panels).In setup 2, α coll is similar to α, but the slope of the S (T ) correlation is not constant, but changes from 1.49 at low T to 1.82 at high T .In setup 3, all three critical exponents α, α 1 , and α coll are consistent with each other.It is remarkable that simple recipes mimicking the conservation of angular momentum in real disk-wind systems and feedback between the nuclear wind and accretion flow are able to robustly drive the system toward an SOC state.In setup 4, the avalanche profile is more regular than in setup 1 and closer to a parabolic shape, confirming that the inclusion of a threshold for gas accretion is enough to regularize the flow.The scatter of the points is much higher than in setups 2 and 3, however, and the avalanche profile does not fully collapse onto a universal function, indicating that the system still has not reached an SOC state.We wish to stress that the SOC state does not require the symmetry of the avalanche profile, but just the collapse onto a universal function.
The avalanche profiles in setups 2, 3, and 4 approach a parabolic shape (Fig. 7), with a slight rightward asymmetry.This asymmetry in the literature has been detected in random processes with short and long temporal correlations in the processes (Baldassarri et al. 2003;Colaiori et al. 2004).Following Nandi et al. (2022), to quantify the level of asymmetry, we fit the average avalanche shapes in setups 2 and 3 with the asymmetric function where a is the measure of asymmetry.We found very good fits in both cases, with a best fit a = 0.24 and 0.28 for setups 2 and 3, respectively.These values are similar to those reported in Nandi et al. (2022) for a completely different system (a neuronal network).A higher value (a ∼ 0.5) is obtained for setup 4.
To further investigate the role of feedback in bringing the system toward an SOC state, we analyzed the distributions of the waiting times between one avalanche and the following one δt.The shape of this distribution can distinguish between simple Poisson processes (as in setup 1) and correlated processes (as in setup 3), as discussed in de Arcangelis et al. (2006) and references therein.For pure Poisson processes, the distribution of δt follows an exponential law, while for correlated processes, more complex forms can be obtained, including power laws, power laws with exponential tails, and even nonmonotonic functions.Figure 8 shows the residuals after subtracting from the δt distributions of setups 1, 2, and 3 their best-fit exponential function.The distribution of the waiting times of setup 4 is identical to setup 1.While the δt distributions obtained from setups 1 and 2 are consistent with an exponential form and display roughly regular residuals, the δt distribution from setup 3 clearly has a more complex shape, confirming that feedback can drive the system A36, page 11 of 19  toward an SOC state, where temporal correlation are present, even with such a minimal modeling approach.

Discussion
The emergence of power laws, the intrinsic nonlinearity of the system, and the presence of feedback at several scales suggest that disk-wind systems could be further investigated as complex dynamical systems.We therefore developed a cellular automaton implementing the main characteristics of a disk-wind system.We studied four setups.Setup 1 is similar to the original cellular automaton introduced by Mineshige et al. (1994a,b), Mineshige & Negoro (1999).When analyzed with statistical tools inspired by critical phenomena, it is found to be supercritical, that is, with an excess of large avalanches, flattening the f (S ) distribution (Fig. 5) and producing an average avalanche profile far from the expected universal scenario (Fig. 7 left panel).The reason is that avalanches cannot exit efficiently from the system in our geometry and that feedback is lacking.The only side from which an avalanche can exit the system is the innermost radius of the disk.We introduced in setup 2 a prescription to calculate the mass that leaves the disk in a wind, inspired by MHDW systems.This reduced the accumulation of mass in each given cell, and thus reduced the avalanches size and duration, steepening the f (S ) and f (T ) distributions.The system approached an SOC state, in which avalanches can hardly extend over the full size of the system.The average avalanche profile now approached a parabolic shape (Fig. 7, central panel).However, the scaling between S and T is not a pure power law, and the critical exponents α and α 1 in Table 1 are not consistent with α Coll .
In setup 3, we further introduced feedback between the outflowing wind and the mass that accreted at the outermost radius.This further reduced the inflow rate and thus the accumulation of mass in each cell.A lower mass accumulation implies that the mass of a given cell may not be brought to a value greater than M crit , so that this cell can give rise to accretion at a later time, thus restarting some avalanche after it has stopped.This increased the duration of these avalanches with respect to setup 2. The net effect was that this form of feedback brought the system toward an SOC state.The scaling between S and T is now a pure power law, and the critical exponents α, α 1 and α Coll in Table 1 are now more consistent with each other.
We introduced in setup 4 a prescription to account for the mass that left the disk in a wind similar to the Eddington luminosity threshold and thus inspired by RDW.This was again able to regularize the flow, steepening the f (S ) and f (T ) distributions.The system is closer to SOC than that in setup 1, but worse than in setup 2 and 3 (Fig. 7).
Some type of feeding regulated by kinetic feedback (whether magnetic or radiatively driven seems less relevant) is truly relevant in this system in general.The more intermittent and fractallike the process (e.g., CCA), the closer it approaches SOC characteristics.
Our analysis has two main limitations.The first limitation is the rather qualitative framework in which observations are interpreted.A better approach requires the comparison of the observations to a robust physical model.We plan to extend the WINE model (Luminari et al. 2018(Luminari et al. , 2020(Luminari et al. , 2021;;Laurenti et al. 2021) by including the three main assumptions discussed above and a full MHD structure.The second limitation is that the sample of AGNs with a UFO detection is naturally biased toward detections.These limitations may be circumvented by using the well-defined and unbiased SUBWAYS sample of 24 AGN at z = 0.1-0.5.XRISM2 (Kitayama et al. 2014), and on a longer timescale, Athena3 (Nandra et al. 2013;Ettori et al. 2013) highresolution microcalorimeter observations will revolutionize our understanding of the inner AGN scale with respect to currently available data, which mostly offer only CCD resolution.Studying the distribution of ω and λ w including upper limits will allow us to estimate the bias in using samples of published UFOs and correct for it, which will allow us to broaden the sample by including both SUBWAYS and previously published UFOs.This will reduce the statistical uncertainty.This is the first of a series of papers that studies the broad problem of AGN-galaxy coevolution from a different perspective, the perspective of complex dynamical systems.Figure 9 inserts popular BH disk-wind models and AGN-galaxy coevolution models into a complexity -spatial scale plane.Following GS17 (updated in Gaspari et al. 2020), we divided the spatial scale into three broad regions, the microscale, spanning roughly from the BH event horizon to well within the BH sphere of influence radius, the mesoscale, around the BH sphere of influence from 10 3 r S to 10 6 r S , and the macroscale, outward up to the edge of the Galaxy and its circumgalactic medium.This paper is dedicated to disk-wind systems and therefore studies microscales.On scales encompassing the BH sphere of influence, between the outer accretion disk and the nucleus of the galaxy (i.e., 100-100 000r S , or between a fraction of a parsec to tens of hundred parsecs) several wind models have been discussed, the RDWs, momentum-and energy-driven winds (e.g., Fauchere-Giguere & Quataert 2012; Zubovas & King 2012), and the entrainment wind phase of GS17.
On macroscales, three main broad scenarios have been discussed in the literature at length.The first scenario, the Eddington-limited BH growth, occupies the lower right part of the diagram in Fig. 9.This scenario dates back to the work of Silk & Rees (1998), Fabian (1999), King (2003), andGranato et al. (2004; see the reviews of Fabian 2012; King &Pounds 2015 andMarsden et al. 2020).Briefly, BH growth can only proceed when AGN winds are not strong enough to displace the accreting gas (negative feedback).On one hand, the wind thrust is proportional to L AGN /c, and L AGN ≤ L Edd , and on the other hand, the weight of the accreting gas is proportional to σ 4 (where σ is the bulge velocity dispersion), assuming an isothermal gas distribution.Thus, for a given σ, the BH can only grow up to a critical mass, giving rise to the M BH − σ relation or to a ridge in the distribution of M BH for each σ, the BH adjustment scenario of Volonteri (2012), and the M BH -X-ray scaling relations (Gaspari et al. 2019), which have been observed to be even tighter than the above optical relations (e.g., see the M BH -X-ray temperature relation).This broad scenario introduces the crucial concept of BH wind feedback that limits the growth of the same BH, but this simplified view has little room for random processes.
The second promising scenario, which we will also study in more depth with a follow-up study, is the CCA introduced by Gaspari et al. (2013Gaspari et al. ( , 2017Gaspari et al. ( , 2019) ) and later by Voit et al. (2015), Voit (2018).Briefly, gas in the BCG halo cools down and condenses into cold clouds and filaments.Random collisions between clouds and filaments can efficiently dissipate their angular momentum, thus enabling efficient inward flows toward the galaxy nucleus and its supermassive BH, starting both star formation in filaments close to the BCG central regions and an AGN cycle.The triggered AGN feedback then transports gas mass and energy back to the macroscale, stopping the star formation, reheating the bulk of the ICM to a higher entropy level, virtually ending also the process, regardless of complexity, that caused the mass deposition rate, and thus completing a single AGN cycle.Thereafter, the gas in the halo, no longer heated, can again condense, giving rise to a new star formation -AGN cycle.This is reflected in the correlations between nuclear activity, star formation rate, and entropy of the ICM in cluster cores, but with large scatter.Correlations of this type have been found: the power to excavate cavities in the ICM is proportional to the X-ray luminosity of the cool core and to the AGN radio luminosity.Similarly, the tight correlations between the SMBH mass and X-ray properties (Gaspari et al. 2019) seem to arise from the CCA raining cycles integrated over several billion years.Only the galaxies in clusters/groups with a low inner entropy (short cooling time) have an active SMBH and are actively forming stars (Cavagnolo et al. 2008(Cavagnolo et al. , 2009)).In this model, BH feedback is again used to regulate the BH growth, and random processes are introduced to both modulate gas accretion toward the galaxy nucleus and its feedback to the gas cooling.In these models, random collisions create a self-similar time variability with a flicker-noise 1/ f power spectrum (Gaspari et al. 2017), similar to SOC systems.CCA is thus an excellent example of a chaotic and emergent process in astrophysics, in which the collective behavior of a large number of microscopic constituents (multiphase clouds) drives a completely novel macroscopic physical process (flickering BH accretion and rain).Other forms of clumpy accretion have been proposed, such as the chaotic accretion through a sequence of randomly oriented merger episodes (King & Pringle 2006;King et al. 2008).Ishibashi & Courvoisier (2009, 2012) also suggested that shocks between accreting clouds can result in X-ray characteristic timescales/spectra with a composition of flicker and red-noise power laws.In the future, it will be interesting to compare the different flavors of chaotic/clumpy accretion models using the tools adopted in Sect. 5.
The third broad scenario includes cosmological simulations using both semi-analytic models (SAM; Fontanot et al. 2020;Menci et al. 2014Menci et al. , 2023 and references therein) and hydrodynamic numerical models (e.g., Di Matteo et al. 2005;Costa et al. 2014;Sijiacki et al. 2015).These simulations follow the growth of density perturbations of the dark matter field and incorporate baryon physics at some level of detail to describe star formation, galaxy and BH growth, and stellar and BH feedback processes.These models by construction include random processes and also include feedback on both galaxy and BH growth from star formation and BH winds and jets.The nuclear accretion in these models is driven by galaxy interactions and by internal galaxy instabilities, such as bars and clumps.These are again random processes.At the heart of star formation and nuclear BH accretion processes are instabilities in the gas disks that can grow and produce strong inward gas flows.Disk instabilities and gravitational disk fragmentation are suppressed or reduced when the fraction of gas to total mass is low (e.g., Lilly et al. 2013).This may occur because AGN feedback removes cold gas from the disk, but other mechanism can be at work as well, such as supernova feedback (especially in low-mass systems), inefficient gas accretion above a threshold galaxy mass (Dekel & Birnboim 2006), and mechanisms impeding inward gas transport (Genzel et al. 2014).
We plan to rethink these scenarios in an attempt to interpret micro-and macroscale observations in the framework of complex dynamical systems by searching for evidence (or absence of evidence) of nonlinearity, universality, and the development of scaling laws, self-organization, criticality, and/or other complex behaviors on scales from the accretion disk in the BH sphere of influence, the bulge, the galaxy disk, and the circumgalactic medium.Finding evidence of the same complex dynamical behavior and universality class (similar critical A36, page 15 of 19 exponents), triggered by BHs on systems from the micro-to the macroscale would highlight the existence of unified principles and would strongly point towards the BH necessity scenario.

Conclusions
We have compiled a large sample of previously published UFO observations and interpreted them in the framework of a scenario including three main assumptions: 1) links between accretion and wind emission as in MHDW and RDW, that is, Ṁaccr = (1/ ω) Ṁw in the former case, or a threshold in disk luminosity above which the accretion is blocked in the latter case; 2) an Eddington-like limit on the mass-accretion rate due to the wind thrust that contrasts the accretion flow at the Bondi radius; and 3) a random distribution of accretion/wind events.We find that the distribution functions of ω = Ṁw / Ṁaccr and λ w = Ṁw / ṀEdd can be described by power laws above some thresholds, suggesting scale-invariance and universal emergent properties in disk-wind systems.Under this assumption, the wind activity timescale is short, much shorter than the AGN timescale, which agrees with other determinations (e.g., Nardini & Zubovas 2018;King & Pounds 2015), and that it is proportional to √ ω.When the population properties ( ω and λ w = Ṁw / ṀEdd ) are drawn from independent snapshots, it might be argued that the wind activity history in each single source can be statistically described by the population functions.If this is the case, the emerging scenario is that of many relatively small wind events for each main wind event in each AGN activity cycle.Because λ = ṀBH / ṀEdd = λ w / ω, λ should also follow a power-law distributions down to an inflection point.The BH growth obtained during each AGN activity cycle is the sum of the growth at each λ, weighted for the time spent at that λ.When this time is similar to the wind activity timescale, then λ∆t ∝ λ w / √ Ω.The slope of the λ∆t distribution is similar to that of the λ distributions estimated by Bongiorno et al. (2012) and Aird et al. (2012Aird et al. ( , 2018) ) for X-ray selected AGN samples, but it is somewhat flatter than that found by Ananna et al. (2022) for a sample of local AGN.
Furthermore, we developed a simple theoretical/numerical model based on a cellular automaton with the goal of testing whether by using simple feedback rules, it is possible to build a dynamical system based on BH physics that approaches the state of SOC, beyond the purely Poisson baseline.In all our simulations, feedback seems to be a necessity to drive the system toward an SOC state.This is encouraging because current models of SMBH feeding, such as CCA, show that fractal distributions and flickering variability are crucial for achieving a consistent self-regulation (Gaspari et al. 2020, for a review).
The main question we started with was whether BH feedback is a coincidence or a necessity for the evolution of galaxies, groups, and clusters of galaxies.Our results consolidate the view that feedback is a necessity for the disk-wind system.Finding similar universality class and statistical properties (scale-invariance or feedback-driven SOC) on large scales as well would point toward the BH necessity scenario.

Fig. 1 .
Fig. 1. ω = Ṁw Ṁaccr vs. Eddington ratios for a sample of 54 UFO systems in 40 AGN.The filled circles represent z < 0.3 AGN, the open triangles show z > 0.3 AGN, and the open circles show lensed AGN.The shaded green area is the expectation of the Sadowski & Narayan (2016), SG17 and GS17 models, and the yellow area is the expectation of the analytic magneto-centrifugal model (e.g., Cui & Yuan 2020).

Fig. 2 .
Fig. 2. Cumulative distribution function of the number of systems in which ω exceeds a limit value.Left panel: 54 individual systems in 40 AGN.Right panel: 43 individual system in 34 sources at z < 0.3.

3.
Normalized, cumulative (red), and differential (per decade; blue) binned distribution functions of ω = Ṁw Ṁaccr (left panel) and λ w = Ṁw ṀEdd (right panel) for 43 individual UFO systems in 34 AGN at z < 0.3 (both corrected for selection effects; see the appendix).The black curve shows the uncorrected ω and λ w cumulative distribution functions.The black curve in the left panel is the same as the black curve in the right panel of Fig. 2.

Fig. 4 .
Fig. 4. Expected probability distribution function of λ (red histogram) and λ × ∆t based on the observed distributions of ω and λ w in a sample of UFO at z < 0.3.

Fig. 6 .Fig. 7 .
Fig.6.Trajectories of three avalanches for each setup (from top to bottom, setups 1, 2, 3, and 4), propagating through the disk.The left panels show avalanches that do not reach the innermost radius, the central panel shows complete avalanches, and the right panel shows the (complete) avalanches with the largest size for that given setup.The central black dot, located at the innermost radius, corresponds to the sink through which the avalanches exit the system.The color-coding indicates the time since the onset of the avalanche.The legend is reported at the top.

Fig. 8 .
Fig.8.Residuals after subtracting the best-fit exponential model from the cumulative distribution function of the waiting times δt between one avalanche and the next.The errors are calculated using Poisson statistics.From left to right, we show setups 1, 2, and 3.The residuals for setup 3 are quite large, indicating that in this case, the distribution function of δt is more complex than an exponential law.This indicates a correlated process.

Fig. 9 .
Fig. 9. Complexity vs. spatial scale plane for several disk-wind models and BH-galaxy coevolution models.