Multiphase AGN winds from X-ray irradiated disk atmospheres
DDraft version January 25, 2021
Typeset using L A TEX twocolumn style in AASTeX62
Multiphase AGN winds from X-ray irradiated disk atmospheres
Tim Waters, Daniel Proga, and Randall Dannen Department of Physics & AstronomyUniversity of Nevada, Las Vegas4505 S. Maryland PkwyLas Vegas, NV, 89154-4002, USA
Submitted to ApJABSTRACTThe mechanism of thermal driving for launching accretion disk winds is interconnected with classicalthermal instability (TI). Indeed, the effective scale height of irradiated accretion disk atmospheres isdetermined by the extent of the cold branch of the S-curve demarcating thermally unstable zones. Ina recent paper, we demonstrated that as a result of this interconnectedness, radial wind solutions ofX-ray heated flows are prone to becoming clumpy. In over two decades of numerical work, however,only smooth thermally driven disk wind solutions that approach a steady state have been found. Inthis paper, we show that the Bernoulli function determines whether or not the entropy mode cangrow due to TI in dynamical flows. Based on this finding, we identify a critical radius, R u , beyondwhich TI should accompany thermal driving, resulting in clumpy disk wind solutions. Our numericalsimulations reveal that clumpiness is a consequence of buoyancy disrupting the stratified structure ofsmooth solutions. Namely, instead of a thin transition layer separating the highly ionized disk windfrom the cold phase atmosphere below, TI seeds the formation of hot spots below the transition layerthat rise up, fragmenting the atmosphere. We find that these hot spots first appear within large scalevortices that form within R u . This results in the continual production of characteristic cold phasestructures that we refer to as irradiated atmospheric fragments (IAFs). These IAFs resemble tsunamisupon interacting with the disk wind, as crests develop as they are advected outwards. The multiphasecharacter of these solutions results from the subsequent disintegration of the IAFs, which takes placewithin a turbulent wake that reaches high enough elevations in the wind so as to be observable fromsightlines as high as 45 ◦ . We discuss the properties of the absorption measure distribution (AMD) indetail, showing that dips in the AMD are not expected within TI zones. We also compute syntheticX-ray absorption lines to show that this complicated dynamics should result in spectral signaturessuch as a less sudden desaturation of O viii Ly α and multiple absorption troughs in Fe xxv K α . Keywords: galaxies: active - methods: numerical - hydrodynamics - radiation: dynamics INTRODUCTIONThe luminosity associated with disk-mode accretionin active galactic nuclei (AGNs) may power a wide vari-ety of the AGN outflows that have been observed (e.g.,Giustini & Proga 2019). In order of decreasing kineticpower, these include extremely high-velocity outflows,broad absorption line quasar outflows, ultrafast out-
Corresponding author: Tim [email protected] flows, several other classes of broad and narrow absorp-tion line winds, X-ray obscurers, and warm absorbers(Rodr´ıguez Hidalgo et al. 2020; Laha et al. 2020, andreferences therein). If the radiation primarily transfersits momentum to the gas, the resulting wind is referredto as being line-driven, dust-driven, or radiation-driven,depending on if the dominant source of opacity is, re-spectively, scattering opacity due to spectral lines, dustopacity, or electron scattering opacity. It is common toemploy a simple treatment of the gas thermodynamics(e.g., an isothermal approximation) when modeling such a r X i v : . [ a s t r o - ph . GA ] J a n Waters, Proga, & Dannen winds, which precludes the possibility that the gas canbecome multiphase due to thermal instability (TI; Field1965). To explain X-ray obscurers and warm absorbers,which are multiphase components within a more highlyionized outflow (e.g., Kaastra et al. 2002, 2014), a morecompelling wind launching mechanism results from ra-diation primarily transferring its energy rather than itsmomentum to the gas, i.e. thermal driving. The appealof this explanation is due to the fundamental connectionbetween thermal driving and TI, our focus in this paper.This connection was first utilized in a pioneering studyby Begelman et al. (1983; hereafter BMS83) to revealhow the irradiation of optically thin gas results in a ther-mally driven wind beyond the Compton radius, R IC = GM bh ¯ mk T C = 5 . × − µM T C / K pc , (1)the distance where the sound speed of gas heated tothe Compton temperature T C exceeds the escape speedfrom the disk. (Here M bh is the black hole mass, M = M bh / M (cid:12) , and ¯ m = µm p is the mean particlemass of the plasma.) Compared to T C , the thermalizedgas in the disk itself is much colder at T ∼ K. Priorto the work of BMS83, it had been shown by Krolik et al.(1981) that the process of heating gas from ∼ K to T C (hereafter referred to as ‘Comptonization’) is highlyprone to TI. In particular, Krolik et al. (1981), Kallman& McCray (1982) and Lepp et al. (1985) calculated avariety of so-called S-curves in the AGN environment,revealing that they commonly have thermally unstableregions (‘TI zones’) at temperatures where recombina-tion cooling, bremsstrahlung, and Compton heating arethe dominant heating and cooling processes.As BMS83 commented, the effective scale height ofan irradiated accretion disk is set not by the temper-ature profile inside the disk but rather by the higheststable temperature on the cold branch of the S-curve(of order 10 K). Above this temperature, TI leads torunaway heating and an accompanying steep drop indensity, thereby facilitating Comptonization by rapidlyremoving any cold gas. The magnetohydrodynamical(MHD) turbulent nature of the disk is not expected toalter this picture because the Compton radius is at dis-tances where the local dissipation energy (and hence themagnetic energy) is far less than that due to radiativeheating. While self-gravity and dust physics may be im-portant within the accretion disks at these radii, theyshould be dynamically unimportant within the irradi-ated disk atmosphere and disk wind, which have com-paratively low density and are hotter than dust subli-mation temperatures.The spectral energy distribution (SED) determinesseveral properties of the irradiated gas, in particular the Compton temperature and the net radiative cool-ing rates relevant for assessing the number and temper-ature ranges of TI zones. How the existence and prop-erties of TI zones depend on the underlying SED wassystematically explored in series of papers by Chakra-vorty et al. (2008, 2009, 2012) with further dependenciesassessed by Ferland et al. (2013), Lee et al. (2013) andR´o˙za´nska et al. (2014). SEDs for individual systems areinferred by fitting the observed spectra to theoreticalmodels for the intrinsic components of an AGN spec-trum (e.g., Done et al. 2012; Jin et al. 2012; Ferlandet al. 2020). In this work, we present calculations foran S-curve obtained in this way for the well studied sys-tem NGC 5548 (see Mehdipour et al. 2015). While thevelocity, density, and temperature structure of thermaldisk winds can be sensitive to the SED, our main resultsregarding TI are expected to hold so long as the S-curvefeatures at least one TI zone.SED fitting calculations for NGC 5548 performed bydifferent groups using photoionization modeling are inagreement that the corresponding S-curves have promi-nent TI zones (e.g., Arav et al. 2015; Mehdipour et al.2015). The SED found by Mehdipour et al. (2015),which has been implemented into
CLOUDY (Ferland et al.2017), yields a Compton temperature of T C ≈ K.The Compton radius for NGC 5548 is therefore about R IC = 0 . §
2) generally have velocities and col-umn densities in the right range (300 − / s and10 − cm − ) to account for the sample of warmabsorbers collected by Laha et al. (2014). However, theionization parameters were found to be systematicallytoo large, so they concluded that there must be an ad-ditional process such as condensation due to TI that canexplain the presence of lower ionization gas at high ve-locity. Note that most previous studies considering therole of TI in explaining the warm absorber phenomenonhave done so in the context of a thermal wind beingdriven from the torus (e.g., Krolik & Kriss 2001; Kall-man & Dorodnitsyn 2019), not from the disk. ultiphase AGN disk winds §
2, we presentthe framework that is now widely used to model ther-mally driven winds from irradiated accretion disk atmo-spheres, and we describe the dynamics of steady statesolutions. In §
3, we present our results. We first connectour study of TI in radial winds to disk winds and derive a characteristic radius beyond which thermally driven out-flows are expected to be clumpy. Then we present ournumerical simulations of clumpy disk winds. We focusour discussion in § § THE STRUCTURE OF IRRADIATEDATMOSPHERES & THERMAL DISK WINDSIn this section, we summarize the flow structure ofsteady state thermal wind solutions obtained using a hy-drodynamical modeling framework in which the effectsof radiation are treated in the optically thin limit. Wefirst review this framework because the overall takeawayfrom this work is that these solutions serve as initialconditions for multiphase disk wind solutions. Addi-tionally, these thermal wind models are complementaryto the latest efforts to model MHD winds together withthe internal dynamics of the disk in full 3D, using eitherMRI turbulent disks (Zhu & Stone 2018; Jacquemin-Ide et al. 2020) or strongly magnetized disks, where amagnetic diffusivity is included to permit accretion (e.g.,Sheikhnezami & Fendt 2018). While the thermodynam-ics used in the above works is highly simplified comparedto that employed here, it may be feasible to combineboth approaches (see e.g., Waters & Proga 2018).2.1.
Modeling framework
The ‘classical disk wind solution’ approach taken hereis to not explicitly treat the disk physics. Rather, aboundary condition (BC) along the midplane is applied,providing a reservoir of mass and angular momentum.A self-consistent choice of parameters and domain sizewill prevent the wind from becoming Compton thick,so that it becomes possible to model how irradiationlaunches a thermal wind using an optically thin heatingand cooling function. As discussed below, given an SED,these solutions are governed by just three dimensionlessparameters. Formulating the problem in terms of theseparameters enables a comparison with our prior resultson radial winds, as well as with an already large bodyof published results.In X-ray irradiated environments, the density BC can-not be chosen arbitrarily; it is constrained by the valueof the pressure ionization parameter, Ξ ≡ (4 πJ/c ) /p ,where J is the mean intensity of the ionizing radia-tion (integrated over frequency and solid angle) localto a given parcel of gas with pressure p . For parsec-scale gas irradiated by an effective point source, 4 πJ Waters, Proga, & Dannen is the radiation flux F X , making Ξ the ratio of the ra-diation pressure p rad = F X /c and gas pressure. Tak-ing the central source to have a total ionizing luminos-ity L X [defined as energies between 1 and 10 Ry; seeKrolik et al. (1981)], the ionizing flux at distance r is F X = L X e − τ ( r ) / πr , where τ ( r ) = (cid:82) r n ( r (cid:48) ) σ ( r (cid:48) ) dr (cid:48) isthe optical depth along the line of sight to the source.Hence, without invoking an optically thin approxima-tion, the density is poorly constrained due to the highlynonlinear dependence on optical depth. AGN warm ab-sorbers are located at parsec scale distances where thisapproximation may apply, thereby reducing the ioniza-tion parameter to Ξ = L X / (4 πr c n k T ).The density in these models is therefore determinedrelative to the midplane BC on the ionization parame-ter, Ξ = L X / (4 πr c n k T ), where r is the radius atwhich n = n and T is the equilibrium temperature onthe S-curve corresponding to Ξ ; n thus evaluates to n = 2 . × ( L X / erg s − )Ξ ( T / K) ( r / pc) cm − . (2)This normalization density is considered to be typicalof the density at heights in the irradiated layers of theaccretion disk where the plasma first becomes Comptonthin. As indicated above, in the classic disk wind solu-tion approach, this value is assigned to the midplane ofthe computational domain at ( r, θ ) = ( r , ◦ ). The in-ner and outer boundaries of the computational domainare set relative to r , and the midplane density profile isgiven by n ( r ) = n ( r /r ) to keep Ξ( r ) = Ξ along themidplane.The fiducial radius r can be expressed in terms of thehydrodynamic escape parameter (HEP), defined asHEP = v c s = GM bh (1 − Γ) rc s . (3)Here, c s = (cid:112) γkT / ¯ m is the adiabatic sound speed, v esc the effective escape speed from the disk, and Γ ≡ L/L
Edd is the Eddington parameter for a system withtotal luminosity L . The √ − Γ reduction to the es-cape speed is due to radiation pressure from the centralengine. Using Eq. (1), we have r = 1 γ HEP T C T R IC (1 − Γ) , (4)where HEP is the HEP assigned to T ; this is the mainparameter controlling the strength of the thermal wind. In Waters & Proga (2018), we describe a modified approachthat introduces a shielded adiabatic accretion disk into the setup,which in 3D would allow combining these irradiated disk solutionswith MRI turbulent disks.
If we insist on having a strong wind beyond r , we re-quire r ≥ R IC (1 − Γ) [Proga & Kallman (2002); thiscriterion is just that of BMS83 for Γ = 0], giving theupper boundHEP < γ T C T (requirement for a strong wind) . (5)To complete the specification of a thermal wind model,we must determine the heating and cooling rates corre-sponding to the radiation field. In recent years, we havedeveloped methods to compute these from the obser-vationally inferred intrinsic SED as self-consistently asis currently possible (Dyda et al. 2017; Dannen et al.2019) using the photoionization code XSTAR (Bautista& Kallman 2001). These calculations require introduc-ing the density ionization parameter, ξ = L X /n H r ,where n H = ( µ H /µ ) n is the hydrogen number density(with µm p ≡ ρ/n and µ H m p ≡ ρ/n H for mass den-sity ρ ). In this work, we use rates derived from theunobscured SED for NGC 5548 (see Mehdipour et al.2015), which has f X ≡ L X /L = 0 .
36. The S-curvecorresponding to this SED has a Compton temperature T C = 1 . × K. Notice that for a fixed SED andelemental abundances (needed to determine the heatingand cooling rates in photoionization calculations), thethree dimensionless parameters governing thermal diskwind solutions are Γ, Ξ , and HEP .2.1.1. Characteristic timescales
It is useful to arrive at characteristic timescales interms of these dimensionless parameters. A fiducialdynamical time is the orbital period at r , t orb ( r ) =2 πr /v φ ( r ), where v orb ( r ) = (cid:112) GM bh (1 − Γ) /r is themidplane orbital speed (and also the escape speed); us-ing Eq. (4), it evaluates to v orb ( r ) = 3 . × (cid:113) ( T , /µ ) (HEP / − , (6)where T , ≡ T / K. The orbital period is therefore t orb ( r ) = 5 . × µ / M (1 − Γ) T / , (HEP / / yrs . (7)The other relevant timescale is the cooling time, t cool ≡ E Λ = 6 . T n C yrs . (8)Here E = c v T is the gas internal energy density,Λ = 10 − ( n/ ¯ m ) C is the total cooling rate inunits of erg s − g − with C the rate in units of10 − erg cm s − , and n = n/ cm − . Using Eq. (2)and Eq. (4), the cooling time at r evaluates to t cool ( r ) = 2 . µ M (1 − Γ) f X Γ(HEP / Ξ C hrs , (9) ultiphase AGN disk winds Figure 1.
An example of the steady state density struc-ture of thermally driven disk wind solutions near the Comp-ton radius R IC . The solution consists of a tenuous outflow(brownish gas) above a somewhat denser corona (yellowishgas) that lies above a dense and nearly hydrostatic cold phaseatmosphere (blue and green gas). A set of streamlines areoverlaid and color-coded according to where they originate:along the inner spherical boundary (gray), parallel to atmo-sphere/disk wind boundary (red), and parallel to the mid-plane at z = 0 . R IC (white). As the red streamlines show,the base of the disk wind is in the upper layers of the at-mosphere where poloidal velocities are less than 10 km / s.The lowest black contour denotes 500 km / s and the next two1000 km / s and 1500 km / s. The white dotted line marks thesonic surface. The circulation apparent in the cold phase at-mosphere increases with radius and permits the formation ofa vortex for r > R IC (see Fig. 5), which in turn triggers amultiphase outflow at even larger radii. where we have eliminated L X in favor of f X using L X = f X Γ L Edd . Comparing Eq. (7) and Eq. (9), we find that t cool (cid:28) t orb for the densest gas along the midplane.An important property of thermal wind solutions isthat they are scale-free for a fixed ionizing radiationfield. That is, they apply to any mass black hole sys-tem if the SED is assumed to be unchanged. The onlyterm that can potentially break the scale-free nature ofthe hydrodynamic equations is the non-adiabatic sourceterm. When non-dimensionalized, it enters the equa- The implicit dependence of the cooling rate ( C in Eq. (8))on density and temperature can make these solutions sensitive tothe S-curve and hence to the SED, but because AGN SEDs areoverall similar yet uncertain for any given system, a good firstapproximation can be found by simply adopting a representativeS-curve and applying it to different systems. Figure 2.
Phase diagram showing the S-curve correspond-ing to the SED of NGC 5548 (solid black curve). Blue dotsand vertical lines mark Ξ (the midplane BC on Ξ) for eachof the solutions in Fig. 3, log(Ξ ) = 0 .
5, log(Ξ ) = 0 .
92, andlog(Ξ ) = 1 .
05. The Balbus contour is plotted as a dashedgray curve; all points to the left of it are thermally unsta-ble and to the right thermally stable according to Balbus’criterion for TI (Balbus 1986). The two unstable regions en-countered while tracing the S-curve starting from low tem-peratures are referred to as the upper and lower TI zones.The base of the lower TI zone is conventionally denoted asΞ c , max and is marked. Overplotted on this ( T, Ξ)-plane arethe phase-‘tracks’ for each of the streamlines in Fig. 1. (Thewhite streamlines are plotted in cyan here.) Notice that thedisk wind streamlines pass through both TI zones. tions multiplied by t /t cool ( r ), where t is a charac-teristic dynamical time such as t orb ( r ) above or R g /c ( R g = GM bh /c ). Because both t and t cool ( r ) areproportional to M bh , the equations have no mass depen-dence. This is ultimately due the density scaling as n ∝ M − in the ionization parameter framework (see Eq. (2)and Eq. (4)). Notice also that t cool ( r ) /t orb ( r ) ∝ HEP − / , implying that the entropy equation will be-come increasingly stiff in a stronger gravitational field.2.2. Basic flow structure of irradiated disks
As shown in Fig. 1, the density structure that resultsfrom this modeling framework consists of three distinctregions: a dense disk atmosphere (blue to green regions),a less dense corona (light yellow to dark brown regions),and a tenuous disk wind (light brown to whitish regions).That the solution consists of these basic componentswas already apparent from the early study by Woodset al. (1996; see also Proga & Kallman 2002), but thefocus subsequently shifted to simulating only the corona
Waters, Proga, & Dannen and disk wind regions to determine if the velocities ofthese models can reach those observed in low-mass X-ray binaries (e.g., Luketic et al. 2010; Higginbottom &Proga 2015; Higginbottom et al. 2017). As discussedin § just belowΞ c , max on the S-curve. In Fig. 2, we mark the locationof Ξ c , max on the S-curve used in this work. This S-curve was calculated from the SED of NGC 5548 (seeMehdipour et al. 2015) by Dyda et al. (2017). We alsoshow three values of Ξ ; the midplane BC of the solutionshown in Fig. 1 corresponds to the middle blue line atlog(Ξ ) = 0 .
92 ( ξ = 50).The velocity structure can also be grouped into threesets of streamlines (although these groupings do notcorrespond to that of the density field): those origi-nating along the inner boundary of the computationaldomain, those with footpoints nearly parallel to theatmosphere/disk wind boundary, and those along themidplane. We show this streamline geometry in theFig. 1, and we plot the temperature along these stream-lines on the ( T, Ξ)-plane in Fig. 2. As the gas becomesless bound, the velocity field transitions from stream-lines with appreciable poloidal divergence for 1 R IC (cid:46) r (cid:46) R IC , to nearly parallel streamlines for r > R IC .Also, contours of constant velocity (black lines de-note 500 km / s, 1000 km / s, and 1500 km / s) become morewidely spaced. The disk atmosphere is characterized bystreamlines that turn back into the midplane or by seg-ments of streamlines near the footpoints that curve be-fore transitioning into the wind, indicative of bound gaswith some amount of circulation. The caption of Fig. 1provides further details on the velocity dynamics.The highest temperatures are reached near the innerboundary, corresponding to the gray streamlines closestto the Compton branch of the S-curve in Fig. 2. Thishot gas is responsible for launching the fastest outflow athigh latitudes. We will refer to this as the radial wind re-gion because the streamlines do not originate from thedisk. We note that the midplane BC is the source ofmatter for gas along the inner boundary, as the pres-sure gradients above the atmosphere establish the den-sity profile along r in , but this gas then becomes the baseof a disconnected wind. Temperature decreases with dis-tance along these wind streamlines because of the nearspherical expansion, which permits the timescale for adi-abatic cooling to be shorter than the timescale for radia-tive heating. The radiative heating rate scales similarlyto t − . Outside the atmosphere, t cool is shortest in thedisk corona, explaining why gas becomes nearly Comp-tonized only in this region. Radiative heating timescalesare shortest here due to the nearly 1 /r scaling of den-sity in the radial wind region. The streamlines originating near the atmosphere/diskwind boundary show similar behavior on the ( T, Ξ)-plane to the radial wind region streamlines only at largedistances. The temperature peaks occuring between ∼ K and 5 × K are at points beyond the sonicsurface (dotted white line in Fig. 2), and the decreasingtemperature thereafter is again due to adiabatic cool-ing, which is less efficient for non-spherical expansion.The footpoints of these streamlines are located on the S-curve slightly below Ξ c , max , and the temperature sharplyrises with distance upon entering the region of runawayheating. All streamlines pass through both the lowerand upper TI zones (the regions falling to the left of theBalbus contour, the dashed line in Fig. 2).2.3. Scale height dependence of irradiated diskatmospheres on Ξ Recent work applying these solutions to explain theobserved changes in wind velocities accompanying statechanges in low-mass X-ray binaries has pointed out howthe extent of the atmosphere is sensitive to Ξ (see theappendix of Tomaru et al. 2019). This dependence isillustrated in Fig. 3, where we compare temperaturemaps for solutions with ξ = 100, ξ = 50, and ξ = 5(Ξ = 11 .
32, Ξ = 8 .
24, and Ξ = 3 . somewhat below Ξ c , max . Including thecold phase atmosphere incurs much greater computa-tional expense because the cooling time becomes ordersof magnitude smaller, as shown in the bottom panels.Once the atmosphere is present, taking a lower Ξ re-sults in it extending to larger heights above the diskmidplane before the transition into an outflowing coronaand disk wind occurs at Ξ > Ξ c , max .A comparison of these solutions reveals that the verti-cal heights of accretion disks and their irradiated atmo-spheres scale oppositely with temperature. The blackdashed lines on the temperature maps in Fig. 3 mark,for the same midplane temperature, the classical scaleheight of a disk, h = (cid:15) r , where (cid:15) = v esc /c iso =1 / √ γ HEP (with c iso = (cid:112) kT / ¯ m the isothermal soundspeed). A cooler disk implies a thinner, less pressur-ized disk, but a cooler irradiated atmosphere requires itto extend to higher latitudes and be of higher pressure.The latter property follows immediately from the defi-nition of the pressure ionization parameter, Ξ ≡ p rad /p . ultiphase AGN disk winds Figure 3.
Top panels : Temperature maps showing how the height of the disk atmosphere varies with the value of the ionizationparameter assigned to the midplane of the compuational domain ( ξ = 100 .
0, 50.0, and 5.0 correspond to Ξ = 11 . right below Ξ c , max ,which results in a solution like that shown in the left panel. However, the presence of a cold phase atmosphere is necessaryto capture TI in the simulations. In this work, we present results for ξ = 50 solutions. Even smaller values of ξ results ina lower temperature atmosphere that extends to greater heights above the midplane (compare middle and right panels). This T -dependence of the height of the disk atmosphere is opposite that of an underlying disk, as seen by comparing the classicaldisk scale height h = (cid:15)r , marked with the black dashed lines, with the height of the atmosphere. The temperature and velocitystructure of the disk wind, meanwhile, are comparatively insensitive to Ξ . To show this clearly, we again plot the sonic surfaceas well as the velocity contours from Fig. 1. Bottom panels : Maps of the cooling time (in years), calculated for NGC 5548( M bh = 5 . × M (cid:12) ). The velocity structure is again overlaid and the orange contour marks where t cool /t ff = 0 . For a fixed radiation flux, i.e. at a given radius in the at-mosphere, gas in thermal equilibrium at a lower Ξ has ahigher gas pressure. The greater extent of the cold phasefollows because the atmosphere is nearly in hydrostaticequilibrium, so a fixed pressure gradient is required tobalance the gravitational and centrifugal forces. Whenstarting from a higher pressure, therefore, cold phasegas must extend to larger heights to ‘climb’ the S-curveat this particular pressure gradient.It is the runaway heating process that occurs whengas reaches Ξ c , max on the S-curve that is responsible forthe atmosphere transitioning into a disk wind (BMS83).Thus, while the height of the atmosphere is sensitiveto the value of Ξ , the temperature at the base of thedisk wind is the same in each of the panels of Fig. 3. Itis therefore expected that the run of temperature withdistance in the wind is not sensitive to Ξ . Because the velocity structure of the disk wind is set by the temper-ature at Ξ c , max , it also is nearly the same from panel topanel. To see this clearly, we overplot the sonic surface(white contour) and three constant velocity surfaces.2.4. Stability of the disk wind vs. atmosphere to TI
To emphasize how TI can operate differently in thecold phase atmosphere versus in the disk wind, it isuseful to draw a comparison with studies of the cir-cumgalactic medium (CGM), where TI is invoked asa condensation process only. Focus has been givento the importance of the ratio t cool /t ff (with t ff = r √ / (cid:112) GM bh /r the free-fall timescale) in assessing thestability of CGM atmospheres to TI (e.g., Sharma et al.2012; Gaspari et al. 2015). We computed this ratio forthe solutions here, because the relevant dynamical timeis the orbital period, which is just t dyn = π √ t ff . We Waters, Proga, & Dannen find that t cool /t ff < .
01 below the orange contourplotted in the bottom panels of Fig. 3. This ratio is aqualitative diagnostic of whether cooling times are shortenough to permit condensation out of the hot phase,which they are. The fundamental criterion for trigger-ing TI, however, is that gas occupy the TI zone longenough for perturbations to grow. As explained in § t cool /t ff (cid:28) t heat /t ff (cid:28)
1. That is, the only wayfor TI to operate in the atmosphere (at least initially)is for perturbations to undergo runaway heating to formheated layers or pockets of warmer gas. We identifiedthis as the dominant process triggering clumpiness inthe radial wind solutions presented in D20. Because theatmosphere is cooler than the temperature of the TIzone, perturbations can only grow near the disk windinterface where T reaches T (Ξ c , max ). RESULTSThe solutions shown in Fig. 3 only reach a steady stateif the outer boundary is located within about 15 R IC .With larger domain runs, the innermost flow region canstill reach a steady state but the outer regions will not.This turns out to be because the inner flow region isjust within the radius where TI can first be triggereddue to dynamical changes in the atmosphere. In § R u , beyond whichIAFs are expected to form due to TI and enter the diskwind. This radius typically exceeds 100 R IC . In § From clumpy radial winds to clumpy disk winds
It is beyond the scope of this work to develop a com-prehensive theory of ‘dynamical TI’. It is clear, how-ever, that the list of qualitative requirements identifiedby D20 for triggering TI in radial winds holds also fordisk winds. Therefore, we first recall our list from D20.3.1.1.
Necessary conditions for triggering TI in outflows • For radiation pressure due to a central source, thegas pressure along a streamline must decrease asor slower than r − upon entering a TI zone toallow Ξ ∝ / ( r p ) to remain constant or decrease.Otherwise gas will quickly exit the TI zone, as is clear from Fig. 2. In more complicated radiationfields, the gas pressure gradient along streamlineswill still likely need to permit Ξ = 4 πJ/p to varysimilarly with distance in the TI zone. • The flow must not be so fast that perturbationsdo not have time to become nonlinear during theirpassage through a TI zone. • The stretching of perturbations due to accelera-tion terms (from the nonlinear term v · ∇ v ) mustbe less efficient than the amplification of pertur-bations from TI.Bullets (i) and (ii) are due to pressure gradients in theflow field and (iii) to velocity gradients, none of whichare accounted for in the local theory of TI. (For priorapplications of dynamical TI in global flows see Balbus& Soker 1989 and Mo´scibrodzka & Proga 2013; see alsothe 3D simulations of Kurosawa & Proga 2009.)3.1.2. Connection between TI and the Bernoulli function
In D20, we presented four qualitatively different typesof flow solutions (Models A, B, C, & D) in 1D. Themodel A solution had the highest HEP and is suscep-tible to TI only during a transient phase of evolution,while models B, C, & D are clumpy wind solutions thatnever reached a steady state. After improving our nu-merical methods (see Appendix A), all of these solutionsreach a steady state, which allowed us to perform a fol-lowup analysis based on the Bernoulli theorem. Thisanalysis is presented in Fig. 4 and amounts to the fol-lowing empirical result: • The above conditions for an outflow to be unstableto TI cannot be met if the flow enters a TI zonewith a negative value of the Bernoulli function, Be = v c s γ − − GM bh r (1 − Γ) . (10)To understand this result, we consider what must hap-pen to trigger TI starting from a steady state flow field.First recall that isobaric TI only amplifies the entropymode [unless the S-curve has a negative slope on the( ξ, T )-plane; see Waters & Proga 2019], and that theentropy mode, unlike sound waves, is advected with theflow. Therefore, only those streamlines that already passthrough a TI zone can amplify perturbations. Next re-call the close connection between Bernoulli’s theorem While the local dispersion relation for TI is unchanged ina uniform velocity field (due to Galiliean invariance), bullet (ii)arises because TI zones may span only a small range of pressureson the S-curve, and the flow pressure typically drops off rapidly. ultiphase AGN disk winds Figure 4.
Top panel:
Steady state versions of our 1Dradial wind solutions from D20 plotted on the ( T, Ξ)-plane.Plus signs denote sonic points and a black dot marks Ξ c , max .Model A is stable to TI and exits the TI zone, whereas mod-els B and C are unstable and settle onto the S-curve insidethe TI zone. Only models B and C will turn into clumpywind solutions if they are perturbed (by applying linear per-turbations at the base, for example). Bottom panel:
TheBernoulli function (normalized to its value at r out ) versusradius. Plus signs mark sonic points. Black dots mark thelocations where Ξ = Ξ c , max , the entry to the lower TI zone.As discussed in § Be >
0. From thedynamical criterion that a sign change in Be at the locationwhere Ξ = Ξ c , max implies stability, we identified the char-acteristic ‘unbound’ radius R u beyond which a transition toinstability should occur (see § R IC , R u is a prop-erty of the S-curve; the vertical line marks R u = 98 . R IC . and the steady state entropy equation. The heat equa-tion can be expressed as T Ds/Dt = −L , which in asteady state becomes v · ∇ s = − L T , (11)where s is the entropy per unit mass. Meanwhile,Bernoulli’s theorem is v · ∇ Be = −L . (12) By eliminating L , we arrive at Crocco’s theorem( ∇ Be = T ∇ s along streamlines) for steady state so-lutions, but it is clear from the derivation that theadvected value of the Bernoulli function closely tracksthat of entropy. An abrupt jump in the value of Be therefore implies a similarly large change in the steadystate entropy profile.The key insight revealed by this Bernoulli analysis isthe effect this profile has on an advected entropy modethat can become unstable to TI. Only for an entropyprofile that does not suffer a dramatic jump at the en-try to a TI zone can the 3 bullet points from § Be profiles to be smooth at the TI zone.3.1.3. Characteristic radius at which thermally drivenoutflows become clumpy
The sign of Be is negative in the disk atmosphere be-cause this cold phase gas is virialized. It becomes pos-itive in the disk wind because runaway heating makesthe enthalpy large and accelerates the flow, i.e., the firsttwo terms in Eq. (10) dominate the third. For an en-tropy mode to avoid an abrupt change as it enters a TIzone, Be must not change signs, hence our previous bul-let point. The radius at which Be first becomes positivejust as T reaches T c , max ≡ T (Ξ c , max ) from below de-fines the characteristic ‘unbound radius’ beyond whichthermally driven winds are expected to be clumpy. Theexact radius is only implicitly defined by the Bernoullifunction, but as shown in Appendix B, a lower bound is R u = 2 γ γ − γ + 1 T C T c , max R IC (1 − Γ) . (13)Eliminating R IC using Eq. (1) shows R u to be pro-portional to ‘the other’ characteristic distance setby the radiation field for a given black hole mass, GM bh ¯ m/kT c , max . Like R IC , therefore, it is a prop-erty of the S-curve for any particular system; the ra-tio T C /T c , max = 4 . × for our S-curve, giving R u = 140 R IC (1 − Γ) for γ = 5 /
3. The vertical line inFig. 4 marks this radius for our radial wind solutions(which have Γ = 0 . R u is simply the ap-proximate radius at which gas entering the TI zone canescape the system, it characterizes both radial wind anddisk wind solutions.3.1.4. The role of a cold phase vortex
Well within the distance R u , the gas dynamics in thecold phase atmosphere changes significantly, as we show0 Waters, Proga, & Dannen
Figure 5.
The same ξ = 50 solution from Figs. 1 and 3 butshowing the full domain, which extends to 40 R IC . In addi-tion to the streamlines overlaid in Figs. 1, we add a set paral-lel to the midplane with footpoints z = 2 . R IC spaced 1 R IC apart along x (black). Beyond 25 R IC , these streamlines areall outflowing and do not enter the disk wind, showing thatthe disk wind becomes dynamically disconnected from theinterior of the outflowing atmosphere. Within 25 R IC , someblack streamlines flow into a vortex centered at r ≈ . R IC .This vortex spans about 9 R IC , corresponding to ≈ . M bh = 5 . × M (cid:12) ) and 0 .
03 pc for an AGNwith M bh = 10 M (cid:12) . In the lower panel, we zoom-in some-what to better show the vortex, which turns out to be criticalfor the onset of TI. in Fig. 5. For r (cid:38) R IC , the atmosphere transitionsinto a slow outflow. Prior to that, a large scale scalevortex appears in the flow at 16 R IC (cid:46) r (cid:46) R IC ,and this dynamics is associated with a steepening ofthe slope of the atmosphere/disk wind boundary at15 R IC (cid:46) r (cid:46) R IC . We are not aware of any previousauthors reporting on the presence of such a counter-clockwise vortex in the atmosphere, although we notethat (Woods et al. 1996) identified clockwise vortices at10 × smaller radii being driven by the stair-step natureof their computational grid. In our simulations, thisdynamics is driven by the pressure gradients along thesteepest bend in the S-curve at about log( T ) = 5 . Values for HEP = 36, ξ = 50 (Ξ = 8 . . M bh scaling scaling r . × ld HEP − M bh n . × cm − HEP M − v orb ( r ) 3 . × km s − HEP / - t orb ( r ) 4 . × yr HEP − / M bh t cool ( r ) 2 . × yr HEP − M bh Table 1.
Characteristic quantities from § § §
4. These valuesare for NGC 5548 ( M bh = 5 . × M (cid:12) ). Using the HEP -scalings, the values for our HEP = 225 run from Fig. 1 canbe obtained from the HEP = 36 values; e.g., r decreasesby a factor of 225 /
36 = 6 .
25. Using the M bh -scalings, oursolutions can be applied to other black hole masses. Forexample, to rescale to M bh = 10 M (cid:12) (52 × smaller), divideall times by a factor of 52 and multiply the densities by 52.The velocities, temperatures, and ionization parameters areunchanged for the same Ξ and Γ. sarily normal to the atmosphere’s surface, whereas indisk winds they are nearly parallel to it beyond 25 R IC .Within 15 R IC , atmosphere streamlines can connect tothe disk wind, as Fig. 1 shows. However, as dis-cussed above, these streamlines are stable to TI. Beyond25 R IC , the base of the disk wind is disconnected fromthe atmosphere’s interior. The temperature at the baseis that of the lower TI zone, T c , max , but pressure fallsoff with distance, meaning Ξ increases outward and gasdoes not enter the TI zone. Thus, starting from a steadystate at least, the only path for matter to enter the windfrom the atmosphere at radii r > R IC is through thevortex. We sum this up with a final bullet point: • In the case of thermally driven disk winds, pres-sure gradients along the S-curve must drive a vor-tex in the atmosphere, as this allows perturbationsto enter the TI zone and become nonlinear.3.2.
Numerical simulations of clumpy disk winds
Thus far we have examined solutions for a few differ-ent Ξ , all with HEP = 225 and Γ = 0 .
1, correspond-ing to r ≈ / M bh = 5 . × M (cid:12) , the massof NGC 5548 (e.g., Kriss et al. 2019). These domainshave r in = r / . R IC and r out /r in = 40. Given thetheoretical expectations from § § R IC but shouldat least include the transition from bound to outflow-ing streamlines in the atmosphere, which occurs around20 R IC for HEP = 225. The outer boundary should ex-tend beyond R u = 140 R IC (1 − Γ), where clumpy wind ultiphase AGN disk winds Figure 6.
Density maps for our mid-res run capturing the onset of TI and the subsequent formation of an IAF. Velocityarrows and contours at 500 km s − and 10 km s − are overlaid as in previous plots. The Be = 0 contour is shown in red and itcoincides with the contour T = T c , max marking the position of the lower TI zone. To visualize the slow velocity field within thedisk atmosphere, we overplot a second set of velocity arrows (gray) that have 10 × smaller magnitude. This reveals a vortex at( x, z ) ≈ (75 , R IC , the formation site of a hot spot (first visible at t = 1049). This spot forms due to TI (see Fig. 7). dynamics is expected to occur. To have maximal nu-merical resolution at minimal cost, it is useful to mini-mize the dynamic range. Based on these considerations,here we present solutions for HEP = 36, which gives r = 13 . R IC for Γ = 0 .
1. Our radial domain is thenassigned as r in = r with a dynamic range r out /r in = 22,giving r out = 1 . R u .This fiducial setup was arrived at after analyzing morethan a dozen simulations with different domain sizesand governing parameters (HEP , Ξ , and Γ). Theprocesses at play are robust and hence captured by asingle simulation. Our numerical methods are detailedin Appendix A; in summary, we solve the equations ofnon-adiabatic gas dynamics using the Athena ++ code.Our simulations are done in spherical coordinates as-suming axisymmetry, and we employ static mesh refine- ment (SMR). We present results at two resolutions: our mid-res ( hi-res ) run has three (four) levels of SMR.In Table 1, we quote the relevant model parametersdescribed in § from 225 to 36 mostly preserves the dynamics in the diskwind, i.e., our HEP = 36 solutions are nearly equiva-lent to HEP = 225 solutions at locations in the windbeyond where HEP <
36 along the midplane. The loca-tion of vortices in the atmosphere, meanwhile, is some-what sensitive to HEP because the dynamics of gas inthe cold phase depends more closely (due to the higherdensities) on the value of t cool /t dyn ∝ / √ HEP .To better enable a direct comparison with prior stud-ies, in our figures we quote times in units of the Keple-rian orbital period at R IC (2 πR IC / (cid:112) GM bh /R IC ): t kep ( R IC ) = 3 . × µ / M ( T C / K) yrs . (14)2 Waters, Proga, & Dannen
Figure 7.
Analysis of hot spot/bubble dynamics.
Top Panels:
Colormaps of the number density as in Fig. 6 but zooming inon a 50 R IC × R IC region centered on the hot spot. Three streamlines are overplotted (blue, orange, and green). In the leftpanel, the footpoints are located within the vortex and the orange and green streamlines are wrapped around the vortex. Inthe right panel, they are located at smaller radii to illustrate that as the hot spot becomes a hot bubble, it separates from thevortex. Bottom Panels:
Phase diagrams similar to Fig. 2 showing the tracks for the streamlines in the panels above. For thehot spot, the streamlines enter the lower TI zone only (left diagram). For the hot bubble, the streamlines enter both the lowerand upper TI zones (right diagram).
With µ = 0 . t kep ( R IC ) = 844 yrs for M bh = 5 . × M (cid:12) , while t kep ( R IC ) = 16 . M bh = 10 M (cid:12) (the black hole mass assumed in D20). We continue toexpress distances in units of R IC , while number densitiesare computed for M bh = 5 . × M (cid:12) and given in cgsunits. In the caption of Table 1, we provide an exampleof rescaling densities to different black hole masses. 3.2.1. The formation of IAFs
Referring to Fig. 6, the first appearance of unsteadybehavior begins after 10 t kep ( R IC ). Shown are fivestages capturing how the first IAF in this run forms. Thered contour here is the surface where Be = 0, which co-incides with the surface where Ξ = Ξ c , max . This impliesthat streamlines enter the TI zone with Be <
0, and ultiphase AGN disk winds Figure 8.
As in Fig. 6 but for the hi-res run. Rather than through the growth of an embedded hot spot on the bottom edge of avortex, the IAF here initially forms by hot gas from the wind region flowing into the atmosphere along the top edge of a vortex.As the bottom panels show, embedded hot spots also form, indicating that both IAF formation channels operate simultaneously.Additionally, in this hi-res run a turbulent wake forms in the wind downstream of the IAF and a vortex shedding process canbe seen taking place. so TI cannot easily amplify entropy modes (see § x, z ) ≈ (75 , R IC in the 1-st panel gets com-pressed enough to reach the temperature of the TI zone.Because it is circulating, the flow passes through the TIzone multiple times, giving more time for TI to operate.It is apparent from the 2-nd panel that a ‘hot spot’ isforming due to the Be = 0 contour appearing within theatmosphere (as a red ‘dot’). This spot becomes a hotbubble by the 3rd panel. As the bubble expands fur-ther in the next 3 frames, a thin layer of the atmospherepeels off, hence our referring to this as an irradiatedatmosphere fragment.To show that the dynamics just described is attributedto TI, in Fig. 7, we analyze the flow behavior on the( T, Ξ)-plane as in Fig. 2. In the top panels, three stream- lines are plotted on zoomed-in versions of the 2-nd and3-rd frames shown in Fig. 6. All three pass through thehot spot, which is located near the bottom of the vor-tex where the streamlines bunch together, indicative ofcompressional heating. The bottom left panel of Fig. 7shows that these streamlines take several passes throughthe lower TI zone. The bottom right panel of Fig. 7shows that the transition from hot spot to hot bubblecoincides with the flow also passing through the upperTI zone. The right-most tracks on the ( T, Ξ)-plane thatavoid the upper TI zone correspond to the streamlinesexiting the vortex and entering the disk wind.A second formation channel for an IAF occurs in Fig. 6also: rather than an embedded hot spot forming, the fi-nal two panels shows hot gas from the disk wind carvingits way into the atmosphere. This channel is the firstto occur in the hi-res run, as shown in Fig. 8. Notice4
Waters, Proga, & Dannen
Figure 9.
Density map snapshots as in Fig. 6 but showing the subsequent evolution of the flow. The full domain is now plotted.The white circular line marks the ‘unbound radius’ R u that approximates the location in the atmosphere where the contours Be = 0 and T = T c , max no longer coincide. As shown by the red contour, Be > R u ; the flow at these distancesbecomes highly susceptible to TI. The brown ‘spots’ in the atmosphere at t ≥ the IAF begins forming at an earlier time, which is notunexpected because higher resolution can permit fastergrowing entropy modes. The depression at r ≈ R IC starts out as just a ripple on the surface of the atmo- sphere. The ripples, in turn, are due to the small levelof amplification that entropy modes undergo when gaswith Be < § ultiphase AGN disk winds The evolution of IAFs
Notice from both Fig. 6 and Fig. 8 that once an IAF isexposed to the disk wind, it begins to disintegrate, send-ing multiphase gas into the wind. The classic turbulentflow patterns associated with vortex shedding accom-pany this process in the case of the second formationchannel, when a small filament appears at the top ofthe more distant surface of the depression in the atmo-sphere. This budding IAF acts as an obstruction to thewind, resulting in what clearly looks to be a K´arm´anvortex street in the downstream flow. This is visible inthe 3-rd, 4-th, and 5-th panels of Fig. 8.The disintegration of the filaments/crests of IAFscombined with the downstream mixing and gradualheating of the shedded multiphase gas within their wakeswill together be referred to as ‘evaporation’. We usethis term to cast equal emphasis on the wake dynamics,which is what truly distinguishes smooth and clumpyevolution in terms of global wind diagnostics (see § t = 1200, the final snapshot fromFig. 6. The IAF formation dynamics described abovehappens continuously, and the small scale IAFs shownup close in Fig. 6 gradually disrupt the flow at largerradii. As seen in the t = 1600 panel, a larger scaleIAF forms just within R u . As shown in the final threepanels, this IAF becomes increasingly filamentary andseparates from the atmosphere, resulting in a smallerscale clump becoming entrained in the wind. This clumpevaporates on a timescale of about ∆ t = 100, corre- sponding to 10 yrs for M bh = 10 M (cid:12) and 1 . M bh = 5 . × M (cid:12) (see Eq. (14)).This ‘early’ evolution, while useful for examining theformation/disintegration dynamics of IAFs, should beconsidered ‘transient’ behavior that eventually leads to afully developed clumpy disk wind solution starting fromsmooth initial conditions. As shown below, the fullydeveloped solution features long lived hot bubbles withinthe disk atmosphere at the base of the large scale IAFs.In Fig. 9, these bubbles are noticeable in the final fivesnapshots (see the brown gas to the lower left of theIAFs between 100 R IC and 200 R IC ).3.3. Smooth versus clumpy solutions
We now quantify the difference between smooth andclumpy solutions by comparing our runs at late timeswith an earlier phase when the flow resembles a steadystate solution. In Fig. 10, we present a streamline analy-sis similar to that of Fig. 7 but now focusing on the flowat large radii. The red streamline in both temperaturemaps traverses the atmosphere, and as seen on the ac-companying phase diagrams, traces only points on thecold stable branch. For the smooth flow, neither theorange nor blue streamline enter a TI zone. Only thegreen streamline that enters the wind through the atmo-sphere traverses the lower TI zone, and it does so withΞ increasing. This is indicative of the flow still beingmarginally stable to TI at this stage despite this regionbeing beyond R u .For the clumpy flow, the question arises, does TI op-erate only at the formation sites of IAFs or also withinthe IAFs as they are advected downstream? In Fig. 10,the green streamline in the late time snapshot showsthat TI is indeed operational in the distant flow. Thisstreamline passes through a hot bubble and is seen tooccupy both the upper and lower TI zones. Similarly,there are two separate tracks of the orange streamlinein the upper TI zone, one reaching a lower tempera-ture than the other. This is due to TI also taking placein the less hot bubble at r ≈ R IC . Finally, the bluestreamline passing through the crest of a fully developedIAF actually starts off in the upper TI zone, likely dueto its footpoint being in hot gas from the nearby bub-ble. It passes isobarically to the cold phase as the IAFis traversed, and then pressure decreases downstreamfrom there as it enters the yellow gas of the lower diskwind. Thus, while Fig. 6 shows that the dynamics ofhot bubbles and IAFs are coupled, here we see that thebubbles themselves stay heated due to their being ther-mally unstable. This continued heating of the hot gasfurther expands the bubble, forcing more of the cool gasenveloping it to enter the TI zone.6 Waters, Proga, & Dannen
Figure 10.
Comparison of the temperature structure and TI dynamics of smooth and clumpy wind stages.
Top panels :Colormaps of the temperature at 10 orbital periods at R IC when the solution is still smooth (left) and after evolving thesolution for 1 orbital period at r out when we end the simulation (right). For each snapshot, we plot four streamlines anchored atvarious heights at x = 150 R IC . Bottom panels : Phase diagrams showing the tracks for the streamlines in the panels above. Inboth cases, the red streamline only occupies the cold branch of the S-curve. In the clumpy wind snapshot, the green streamlinepasses through the inner hot bubble, while the orange streamline passes through both the inner and outer bubbles, visible astwo separate tracks through the upper TI zone. Both bubbles undergo runaway heating associated with TI.
One can ask a related but more observationally rel-evant question: does a line of sight analysis reveal thesame behavior, namely a paucity of gas at ξ ranges cor-responding to TI zones in smooth flows compared toclumpy ones? To answer this question, we switch toanalyzing the high-resolution run because it yields themost accurate calculation of the absorption measure dis-tribution (AMD). Fig. 11 shows maps of ξ , again for asmooth stage of evolution and also at two later, clumpy stages that capture how the crest of an IAF can breakoff, leaving an isolated, outflowing clump. The uppersightline in the bottom panel passes through this clump,whereas the same sight line at t = 2141 passes just abovethe crest. Thus, there is a much broader range of ξ at t = 2235 compared to t = 2141, as is also apparentfrom the corresponding phase diagrams, which are nowplotted using the ( T, ξ )-plane. There are again multiplepasses through both TI zones for sightlines that inter- ultiphase AGN disk winds Figure 11.
Comparison of the ionization structure and TI dynamics of smooth and clumpy wind stages along two differentsightlines for the hi-res run.
Left panels : colormaps of the ionization parameter ξ for an early smooth phase (top panel) and attwo later stages when IAFs become fully developed. Velocity arrows and contours are again overlaid. The radius R u is shown inwhite. Two sightlines at 54 ◦ and 63 ◦ from the polar axis are drawn; the red portions of these lines trace the radial range r > R u used in subsequent figures. Right panels : Corresponding phase diagrams of gas along the 54 ◦ (left plot) and 63 ◦ sightline (rightplot). The solid gray line is the S-curve. The shaded gray area shows thermally unstable regions according to Balbus’ criterionfor TI; the boundary of this region defines the Balbus contour shown in previous phase diagrams. The upper right panel showsthat sightlines passing through both the disk atmosphere and the base of the disk wind span a large range of ionization statesand only skim the instability region. The smooth flow results in a continuous track on the phase diagram but would still beconsidered ‘multiphase’ if observed along the 63 ◦ sightline. The subsequent clumpy flow regime shows an even larger range ofionization states. Depending on the sightline and evolutionary stage, passing through an IAF can entail traversing brownish,yellowish, greenish, and bluish gas both upon entry and exit (corresponding to multiple passes through the instability region). Waters, Proga, & Dannen
Figure 12.
Comparison of wind properties along θ at r out for smooth and clumpy stages of evolution. These two statescorrespond to the t = 588 and t = 2235 snapshots shownin Fig. 11. The solid and dashed lines in the second panelcorrespond to the total velocity v and its radial component v r , showing that the velocity field of the disk wind region isalmost purely radial. The next two panels show the mass lossrate and kinetic power. Their integrated values are comparedin Table 2. Consistent with the table, we exclude the region85 ◦ < θ < ◦ to allow for an underlying accretion disk. Thedotted vertical line divides the remaining domain equally. sect cold phase gas, implying higher column densities Quantity Smooth Clumpy < . ◦ > . ◦ < . ◦ > . ◦ (cid:104) v (cid:105) [km s − ] 857 123 835 1054 π ˙ m/ ˙ M (cid:12) π ˙ K/L [10 − ] 1.04 0.11 1.04 0.07 Table 2.
Global wind diagnostics of smooth versus clumpywind stages corresponding to the t = 588 and t = 2235snapshots in Fig. 11. The top entries are the mean outflowspeeds at r out , (cid:104) v (cid:105) = (cid:82) v r ( r out ) sin θ dθ/ (cid:82) sin θ dθ . The nextentries are the integrated profiles shown in the lower panelsof Fig. 12: 4 π ˙ m/ ˙ M (cid:12) is the mass loss rate in units of M (cid:12) yr − and 4 π ˙ K/L is the kinetic luminosity in units of L = 6 . × erg s − (the luminosity for M bh = 5 × M (cid:12) assumingΓ = 0 . > . ◦ values exclude a 5 ◦ region above themidplane to allow for an underlying disk. The entries for ˙ m and ˙ K account for midplane symmetry, so their sum gives˙ M and P K , the total values over a 4 π solid angle. of gas in the temperature range of the lower TI zone, T ≈ − × K. For sightlines through the smoothsolution, however, gas only occupies this temperaturerange if the upper atmosphere is skimmed. This con-trast in what a distant observer will see indicates thatsmooth and clumpy wind models are very different, withthe clumpy model defying the standard narrative thataccompanies discussions of the AMD (see § M = 4 π (cid:82) π/ ˙ m sin θ dθ and P K = 4 π (cid:82) π/ ˙ K sin θ dθ , where ˙ m = ρv r r and˙ K = (1 / ρv v r r . These calculations require v r inaddition to ρ and v = (cid:113) v r + v θ + v φ , which is shown bythe dashed curve in the 2-nd from top panel of Fig. 12.Notice that the profiles of v r and v overlap above thedisk atmosphere, indicating that the disk wind regionsare radially outflowing. Within the atmosphere, v φ isdominant, explaining the higher velocity for θ > ◦ . InTable 2, for both smooth and clumpy states, we computevalues to the left and right of the vertical line in Fig. 12,which approximately separates the subsonic, multiphaseflow from the highly ionized, supersonic flow above.That the solutions are so similar in terms of bulkwind properties is merely a reflection of the underly- ultiphase AGN disk winds DISCUSSIONWe found here that the standard framework for mod-eling thermally driven disk winds allows steady statesolutions to be obtained only when the computationaldomain excludes the cold phase atmosphere and the ac-companying large scale vortices that form at r (cid:38) R IC .Clumpy disk wind solutions result otherwise, and weshowed these solutions to be analogs of the clumpy ra-dial wind solutions obtained by D20.TI plays two related but distinct roles: (i) it firstcauses the formation of small hot spots in the upperatmosphere that subsequently expand, becoming largerhot bubbles; and (ii) it facilitates continuous heating ofnew gas entering the bubbles. As a result of (ii), thebubbles further expand and rise, causing a fragmentedlayer of the atmosphere to be lifted into the disk windabove. The resulting structure, which we refer to asan IAF, forms crests and filaments that can break off,sending smaller clumps into the wind. The evaporationof the IAF and any resulting clumps supplies the diskwind with a substantially higher column of lower ioniza-tion gas than is present in the smooth solutions at thesame viewing angle.IAFs are large-scale structures and the dynamics justdescribed occurs on timescales spanning 10 − years.According to these models, the system is frozen in astate like those shown in the bottom two panels ofFig. 11 on timescales of years. In § § § Absorption diagnostics
Ionized AGN outflows can show absorption from ionsthat are widely separated in ionization energy. Fromthe distribution of ionic column densities, it is possible to reconstruct the AMD (Holczer et al. 2007), which isthe distribution of hydrogen column density N H as afunction of the density ionization parameter, ξ . Con-versely, from a model of the gas distribution along theline of sight, the AMD can be computed as the slope ofthe column density N H when plotted against log( ξ ),AMD ≡ dN H d log( ξ ) . (15)Because N H is a monotonically increasing function bydefinition, the AMD is a positive quantity that shouldbe overall large in high density regions characterized bysmaller ionization parameters, such as inside an AGNcloud. The generic shape of the AMD across regionswith sharp density gradients, such as the surface ofclouds (IAFs in these models), is unclear. We there-fore present an in depth analysis of the AMD for thetwo sightlines discussed already in § I , I ν = I e − τ ν ( θ ) , (16)where τ ν ( θ ) = (cid:82) κ ν ( r, θ ) ρ ( r, θ ) dr is the optical depthat r out for a fixed inclination θ (the viewing angle for adistant observer). These calculations therefore neglectcontributions due to the wind’s intrinsic emission alongthe line of sight, as well as scattering into the line ofsight. Following the methods of Waters et al. (2017), weevaluate κ ν ( r ) = κ ( r ) φ ( ν ) /φ ( ν D ) by post-processingour simulation results to Doppler broaden the lines ac-cording to the radial velocity and temperature profiles; φ ( ν ) is the (radially dependent) Gaussian profile witha width ∆ ν = ν v th /c set by the thermal velocity v th ( r ) = (cid:112) kT ( r ) /m ion and with the line center fre-quency Doppler-shifted to ν D ( r ) = ν [1 + v r ( r ) /c ]. Weuse lookup tables generated from XSTAR output to ob-tain the line-center opacity κ ( r ) = κ [ ξ ( r ) , T ( r )] foreach ion. These opacity tables are generated from thesame grid of photoionization calculations used to tabu-late our net cooling function L that defines our S-curve.0 Waters, Proga, & Dannen
Figure 13.
AMD analysis for the two sightlines shown in Fig. 11. The top, middle, and bottom subpanels show the columndensity, the AMD, and the column density weighted line of sight velocity (as defined in the text). Calculations for the fullsightline are shown in the left panels, while the right panels show a partial calculation only along the red portion of the sightlineshown in Fig. 11, covering the range r > R u . Black, blue, and orange curves correspond to the three snapshots shown in Fig. 11.The ξ -ranges of the upper and lower TI zones are marked by pairs of solid and dashed vertical lines in the AMD plots. ultiphase AGN disk winds Absorption measure distribution (AMD)
The middle panels in the left column of Fig. 13 showsthe AMD computed along the two sightlines plotted inFig. 11. To aid our analysis, in the right middle panelswe also show the AMD over just the outer portion ofthe sightline extending from r > R u (the red portionin Fig. 11). Additionally, in the upper panels we plotthe cumulative column density to show that the shapeof the AMD is consistent with the slope of these curves.Finally, in the lower panels we plot the column den-sity weighted radial velocity v r ≡ (cid:82) ∆ N v r dξ/ (cid:82) ∆ N dξ ,where ∆ = 0 . ξ ) overwhich the integrals were evaluated; v r is indicative ofthe expected blueshifts of absorption lines.Focusing first on the smooth state (upper panel inFig. 11), notice that starting from the outer bound-ary, the two sightlines trace gas with increasing ioniza-tion as r decreases. These radial profiles are thereforein direct correspondence with the AMD. For the 63 ◦ sightline, the cold phase atmosphere gas (blue color) isall at the same low ionization and there is a large col-umn of it compared to the less dense gas beyond (greenthrough brown colors), explaining the more than two or-der of magnitude dip in the AMD. The AMD of the 54 ◦ sightline spans a much smaller range of ξ and is overallsmall (a few 10 cm − ) as only low density gas is beingprobed. There is a dip upon exiting the upper TI zonethe ξ -range of (marked with vertical dashed lines) be-cause the higher ionization gas beyond has a smaller risein column density, as shown in the top panels of Fig. 13.Before discussing the AMD for the clumpy states,we assess how the analysis just given fits in with themost common discussion of the AMD in the literature(‘the standard narrative’ hereafter), in which referenceis made to the theory of local TI. Namely, dips in theAMD are interpreted as indicating thermally unstable ξ -ranges, the outcome of local TI being to replace theseranges by disparate regions of the S-curve that are con-nected isobarically. Below we will discuss in more detaillocal versus ‘dynamical’ TI, but here we emphasize thatthe dip in the AMD is not due to dynamical TI becausethe flow is still smooth at this stage. While the pro-nounced dip in the AMD for the 63 ◦ sightline is merelya column density effect, it is not coincidental that it oc-curs near ξ c , max , the entry to the lower TI zone. Theflow structure along this sightline actually epitomizesthe intrinsic connection between thermal driving and TIdiscussed in §
2. To recap, the presence of a TI zone is Technically, in the theory of local TI, gas will not typicallyreach the upper stable branch of the S-curve due to the effects ofthermal conduction (see Proga & Waters 2015). responsible for the sharp transition from a bound atmo-sphere to a thermally driven wind at ξ = ξ c , max . Gaswith T > T c , max is much more tenuous in a steady equi-librium state. Hence, for sightlines through the atmo-sphere, the AMD will naturally exhibit a dramatic falloffonce ξ c , max is reached. Compared to the standard nar-rative, therefore, pronounced dips in the AMD here arecaused by the full structure (see §
2) of thermally drivenwind solutions being probed.The most basic result in Fig. 13 is that the AMD forclumpy states span a larger range of ξ . The AMD alwaysextends to lower ξ for clumpy compared to smooth statesbecause IAFs undergo compression. It can extend tolarger ξ when gas within bubbles are probed. We showedbubbles to be the primary sites were TI takes place,hence they are subjected to continuous runaway heatingthat increases ξ .Notice that beyond ξ c , max , the AMDs for both clumpystates undergo dips or remain relatively flat between thepairs of vertical lines marking the ‘unstable’ ξ -rangesthat correspond to the upper and lower TI zones. Dipscan also occur beyond the upper TI zone. For smoothstates, the AMD rises slightly within the upper TI zone.These features are inconsistent with the standard nar-rative. Rather than featuring a dip within the unstable ξ -ranges, the AMD shape for smooth, thermally drivenwind solutions is simpler still: gas continues to occupythe ‘unstable’ ξ -range of the S-curve because this samerange is actually stable just slightly off the S-curve. Asshown in the top panels of Fig. 11, gas does not actuallypass through the upper TI zone; hence for t = 588, theflow is formally stable in the range between the verticaldashed lines in Fig. 13, and it lies below the equilibriumcurve in a region of net radiative heating. A steady statedynamical equilibrium is maintained by this heating be-ing balanced by adiabatic cooling due to expansion. Forthe 63 ◦ sightline, gas also passes through the lower TIzone, but it has not yet become unstable due to the pre-dominance of the stabilizing effects discussed in § ξ beyond the upper TIzone are simply due to the overheated gas within bub-bles having low column densities.To understand the AMD of the clumpy flow states indetail, it is helpful to contrast the regimes of local TIand dynamical TI. The depletion of gas within TI zones— the principle outcome of local TI — still takes placewithin global outflow solutions once they cease to besmooth. In local TI, there is finite supply of the initialunstable gas reservoir, while in global flows, the TI zonecan be continually replenished with new gas. In papersthat reference the theory of local TI when interpretingthe AMD (e.g., Behar 2009; Detmers et al. 2011; Ad-2 Waters, Proga, & Dannen hikari et al. 2019), the implicit assumption is thereforebeing made that thermal timescales are much shorterthan dynamical ones to allow depletion to dominate re-plenishment. This would entail, for example, that a hotbubble in pressure equilibrium with the cold atmosphereforms rapidly (compared to orbital timescales) at the on-set of TI and then no longer evolves thermally, therebykeeping the TI zone devoid of gas. Such a scenario isinconsistent with the force equation, however, becausehot bubbles tend to buoyantly rise. The dynamics ac-companying the bubbles appearing in our simulationsindeed involves the effects of buoyancy, as this is whatleads to the formation of IAFs (see § outside of TI zones at even larger ξ . Thisdip probes the small column of gas within a hot bub-ble. Referring to the lower (63 ◦ ) sightline in the bottompanel of Fig. 11, the corresponding partial AMD (or-ange curve in bottom right AMD panel of Fig. 13) showsprominent dips at ξ ≈ . × and ξ ≈ . × .The lower ionization dip is from the entrained bubbleat r ≈ R IC , while the higher ionization one is thedarker brown gas between R u and the first IAF encoun-tered along this portion of the sightline. This latter gasis seen to originate from another bubble located justbelow the sightline that is also subjected to TI. How-ever, both of these dips become mostly filled when weinclude the full sightline because the gas at small radiicovers the same ionization parameter range and is alsodenser. This analysis therefore casts serious doubt onthe whole notion that the AMD can be used to drawany conclusions about TI. In principle it can, but thereis degeneracy with ‘normally’ heated gas now that wehave shown that TI should not lead to dips within the ξ -ranges of TI zones.The AMD of the same sightline but for the previoussnapshot (blue 63 ◦ line) does show a wide dip in theTI zone for the partial AMD calculation. This dip alsogets filled in upon including the full sightline, but whatis interesting is the red portion of the sightline tracesonly one region of yellowish/brownish gas — that of theevaporating flow of the IAF just within r = R u . This isthe most highly ionized gas probed and is seen to be atthe ξ marking the location of the dip. Furthermore, theevaporative layer of the other much thicker IAF along this sightline contains the less ionized (light green col-ored) gas with ξ ≈
800 that is within the same TI zone.The primary dip in the AMD near ξ c,max is unambigu-ously the dark green evaporation layer. Thus, we havearrived at the expected AMD shape of a single IAF: anoverall peak for ξ < ξ c,max , a large dip near ξ = ξ c,max (that corresponds to ‘exiting the atmosphere’ as in theAMD for the smooth phase) followed by secondary dipsat larger ξ that often coincide with TI zones but areattributable to a sharp falloff in column density in theIAF’s wake rather than to TI.The orange AMD curves in the top panels of Fig. 13provide another example of this AMD shape. This isthe upper sightline in the bottom panel of Fig. 11 thatpasses through an isolated clump. Such clumps arisewhen the crest of an IAF becomes detached after beingablated by the wind. This is clearly the most low ioniza-tion gas along this sightline, showing that the primarydip in the AMD beginning near ξ c,max is the signatureof this clump. There would be a secondary dip in theupper TI zone but it is filled in by the larger column ofbrownish gas at smaller radii. Notice that this evapo-rative structure is reminiscent of the cometary tails dis-cussed in the context of occulation events (e.g., Bianchiet al. 2012). The evaporating wakes here are depicted asmoving transverse to the line of sight rather than alongit, but the ionization structure they inferred is never-theless consistent with the dynamics taking place in oursimulations.In summary, the AMD for clumpy multiphase diskwind models can be understood as a combination of thatexpected for a smooth outflow with that for one or moreembedded IAFs (which are equivalent to discrete cloudswhen probed only along the line of sight). Specifically,the AMD shape has several variants depending on thesightline. For sightlines intersecting cold phase gas, in-dicative of either passage through the atmosphere or anIAF, we find • a nearly flat distribution for ξ < ξ c,max corre-sponding to compressed cold phase gas; • a pronounced dip just above ξ c,max , implyingsomewhat larger blueshifts because IAFs are sub-jected to ram pressure by the disk wind.For sightlines passing through hot bubbles, where TI isactively taking place, • a dip should occur at ξ larger than that of any TIzones (for log( ξ ) > • the low- ξ range of the AMD does not extend to ξ c,max and blueshifts are systematically higher. ultiphase AGN disk winds Figure 14.
Synthetic X-ray absorption line profiles calculated for the 54 ◦ (top panels) and 63 ◦ (bottom panels) sightlinefor each snapshot in Fig. 11. Shown are profiles for the K α resonance lines of Fe xxvi and Fe xxv followed by the doubletsSi viii Ly β and O viii Ly α ; the individual components of these blended doublets are shown as blue and red colors. Bottomsubpanels are scatter plots showing ξ as a function of the line of sight velocity (the x-axis is v los = − v r ); red dots correspondto the red portions of the sightlines in Fig. 11. Hashed areas mark the ξ -ranges of the upper/lower TI zones. Waters, Proga, & Dannen
The last point was not discussed above, but it is evidentfrom the blue curve in the upper panels of Fig. 11 andis relevant to the synthetic line profile calculations thatwe now present.4.1.2.
Synthetic absorption line profiles
We frame our discussion of X-ray absorption line pro-files around the four bullet points above. These AMDproperties correspond, respectively, to the following ex-pectations: • Deeper absorption should occur for unsaturatedlines in clumpy compared to smooth states for ionsprobing the cold phase gas. • There should be a spectral signature correspond-ing to the prominent dip near ξ c , max : for ions withpeak abundance at ξ > ξ c , max , a gradual transitionfor core to wing as the line desaturates; for ionswith peak abundance at ξ < ξ c , max , a sharp desat-uration should occur right at ξ c , max . For the 54 ◦ sightline here, this occurs at v ≈ −
200 km s − . • Deeper absorption should occur at the lower veloc-ities characterizing gas in the atmosphere for ionsdirectly probing gas in the hot bubbles (Fe xxv K α and Fe xxvi K α for our S-curve). For ions probingthe cold phase gas, less absorption should accom-pany bubbles due to the reduction in atmosphericcolumn density caused by their presence. • Multiple absorption troughs should mimic theshapes of absorption line doublets for ions prob-ing the evaporative wakes of the IAFs, due to theAMD showing dips at different velocities.In Fig. 14, line profile calculations are presented forions with peak abundance in the cold phase at ξ < ξ c , max (O viii ), in the evaporative wakes of IAFs at ξ > ξ c , max (Si viii ), and in the hot bubbles or the highly ionized gasat small radii (Fe xxv , Fe xxvi ). The third and fourthbulleted expectations are readily seen to be realized forFe xxv K α and Fe xxvi K α . The abundance of Fe xxv peaks in what is essentially the brightest yellow gas inFig. 11, while for Fe xxvi the peak is for the brown(higher ionization) gas, which has higher velocity. Thisexplains the overall line shapes for the 63 ◦ sightline(bottom panels of Fig. 14). Comparing the amountsof brown gas along the red portion of this sightline inFig. 11, the t = 2235 snapshot should show deeper ab-sorption in Fe xxvi K α at velocities corresponding tothe bubbles ( v ≈
125 km s − ) compared to t = 588.This is clearly the case, and it occurs more prominentlyeven for Fe xxv K α , although for Fe xxvi K α the effectserves to make it appear as though this singlet is blendedwith another line of comparable opacity. For the 54 ◦ sightline, the shape of Fe xxv K α at t = 2235 is such that it masquerades as a blended 2:1doublet (where the blue component has an oscillatorstrength twice that of the red component). Si viii Ly β is an actual 2:1 doublet line with blended components.The effects of clumpiness are mostly concealed for it be-cause this line mainly traces the green gas in Fig. 11,which is confined to the outer portion of the sightlineand hence rather localized in velocity space. However,comparing the shapes of the Si viii Ly β profiles for the63 ◦ sightline at t = 588 and t = 2235, the effects ofevaporation dynamics in IAFs can be inferred: a moregradual transition from core to wing implies increasedion opacity across a broader velocity range. This reflectsthe underlying IAF dynamics, where the warmer gas inthe wake is being accelerated to higher velocity.This spectral signature for evaporation dynamics iseven more prevalent in the line profiles for O viii Ly α ,which is another prominent 2:1 X-ray doublet line. It isvery optically thick in the cold phase gas (blue colors inFig. 11), while it is marginally optically thin to evaporat-ing gas (green and yellow colors). For the 54 ◦ sightlineat t = 588 and t = 2141, Fig. 14 therefore shows thatthis line is unsaturated because there is no cold phasegas present along the line of sight. For t = 2235 in theupper set of plots, the line becomes saturated due to thedense clump in the wind (see Fig. 11). The wing is nolonger Guassian-like as in the profile at t = 588. Thisextended wing feature in clumpy compared to smoothsnapshots is seen most clearly for the profiles of the 63 ◦ sightline. Note the difference with the extended wings ofthe Fe xxvi K α profiles. Fe xxvi is not a direct tracerof evaporation dynamics, as its K α line forms in thehottest gas far downstream in an IAF’s wake.The O viii Ly α line profiles for the 54 ◦ sightline ex-hibit the expected behavior given in the first bulletabove. There is no cold gas along this sightline forthe smooth state and progressively more in subsequentstates. We note that the Si viii Ly β profiles are also con-sistent with our expectations despite showing oppositetrends in moving from t = 588 to t = 2235 when com-paring both sightlines. Because this ion’s abundancepeaks in the evaporative wakes (green gas) rather thanin the cold phase, it deepens in sightlines passing aboverather than through the atmosphere.As for the second bullet, here we have only consideredions with peak abundances above ξ c , max , so the core towing transition is gradual rather than sharp. In otherwords, a prominent dip in the AMD results in a sharploss of opacity only if the ion abundance is not signif-icantly increasing beyond the dip. For high ionizationUV lines, such as the doublet C iv λλ , ultiphase AGN disk winds Dynamical TI in the presence of other forces
We previously showed that radiation forces alter thebasic outcome of local TI in an anisotropic radiationfield (Proga & Waters 2015). Line driving at the dis-tances where IAFs form appears to be too weak (dueto over-ionization and large optical depth in lines) tosignificantly affect these multiphase disk wind solutions(Dannen et al., in preparation). It is useful nonethe-less to draw a comparison with these previous cloud ac-celeration simulations to address the conflicting claimsregarding the role of TI in both the broad line region(BLR; e.g., Beltrametti 1981; Mathews & Ferland 1987;Mathews & Doane 1990; Elvis 2017; Wang et al. 2017;Matthews et al. 2020) and intermediate to narrow lineregions of AGNs (e.g., Krolik & Vrtilek 1984; Stern et al.2014; Goosmann et al. 2016; Bianchi et al. 2019; Borkaret al. 2021). In Proga & Waters (2015), we were track-ing a single cloud as it formed and accelerated (mainlyby line driving), and we envisioned that globally, theenvironment of this cloud would be similar to that sur-rounding the isolated clump shown in the bottom panelof Fig. 11. Rather than the result of condensation outof a warm unstable phase, this clump is an ejected pieceof the crest of an IAF. ‘Cloud formation due to TI’has therefore taken on a new meaning in this paper,as clumps ejected from IAFs ultimately owe their exis-tence to TI. It is thus necessary to distinguish betweenlocal TI, in which clouds form via condensation, anddynamical TI. In this work, the latter is only directlyresponsible for forming hot bubbles, which then go onto fragment the cold phase gas in the upper layers of thedisk atmosphere.At much smaller radii where line-driving becomes veryimportant, it is first of all unclear if TI is even rele-vant. In the original discovery of clumpy disk winds byProga et al. (1998), for example, the wind was taken tobe isothermal, thereby showing that clumps can arisedue to radiation forces and gravity alone. However, ithas been proposed that BLR clouds are the result ofcondensation out of the hot phase of line-driven winds(Shlosman et al. 1985), i.e. local TI was envisioned.Let us suppose, therefore, that TI and line driving areinterdependent in the BLR. It then becomes unclear ifthe above distinction between dynamical TI and localTI still needs to be drawn. Waters & Li (2019) ar-gue that the local TI approximation is valid in the BLRand furthermore show that if the flow dynamics can bemodeled using multiphase turbulence simulations, then the resulting ionization properties are consistent withthose of local optimally emitting cloud (LOC; see e.g.,Leighly & Casebeer 2007) models. The properties ofmultiphase turbulence can resolve some but not all ofthe shortcomings of the original two phase model of Kro-lik et al. (1981) identified by Mathews & Ferland (1987)and Mathews & Doane (1990), who concluded that BLRclouds are unlikely to be formed by TI. The remainingissues depend on the validity of the local TI approxi-mation, which will typically break down wherever gra-dients in the bulk flow become steep, e.g., within thetransition layers expected between disks, atmospheres,coronae, and winds. If the line emitting gas is producedwithin such layers, then these issues should be revisitedafter the much richer multiphase gas dynamics accompa-nying dynamical TI under strong radiation and/or mag-netic forces is better understood.Our finding that the sign of the Bernoulli functionat the entrance to a TI zone determines whether or notthe flow is stable to TI is an important development thatwill hopefully lead to a more complete understanding ofdynamical TI. The need for a sign change in Be mayturn out to be unique to thermally driven flows at largedistances where gravity is small. More fundamental isthe accompanying interpretation that individual entropymodes propagating along streamlines get disrupted uponencountering a large gradient in Be . By consideringBalbus’ criterion for TI, (cid:18) ∂ L /T∂T (cid:19) p < , (17)and the fact that in a steady state, L /T = − T − v · ∇ Be by Bernoulli’s theorem (see Eq. (10)), we see directlythat ∇ Be along streamlines is relevant to the onset ofTI. A formal inquiry into this connection requires a La-grangian perturbation analysis, however.Since a generalized Bernoulli function is central to thetheory of MHD outflows (e.g., Spruit 1996), it will beinteresting to determine if a jump in Be caused by mag-netic pressure alone can inhibit TI. In extending thesesolutions to MHD, the disk atmospheres are expected tobe unstable to the MRI. It remains to be seen if the re-sulting MHD turbulence can upend the important roleof the atmospheric vortices in this work, which resultin the episodic production of IAFs. This is a perti-nent issue to address, as these vortices can become muchweaker upon including the pressure gradient necessaryfor a radial force balance in the midplane BC (see Ap-pendix A). While this term is only a small correctionfor a cold disk and is therefore typically omitted (e.g.,Tomaru et al. 2018; Higginbottom et al. 2018), our testruns showed that IAFs develop on significantly longer6 Waters, Proga, & Dannen timescales upon including it. An MRI disk setup willnot fully alleviate the ambiguity regarding this term. We point out that by extending the current model-ing framework to include an MRI disk setup, it becomespossible to investigate one of the most plausible theoriesenvisioned for the origin of BLR clouds — dense clumpsbeing ejected from disk atmospheres by large scale mag-netic fields (Emmering et al. 1992; de Kool & Begelman1995). IAFs, if they can somehow still form at thesesmall radii, can no longer be advected outward as grav-ity is too strong (i.e. the bound in Eq. (5) is violated).It is entirely conceivable, however, that they could fa-cilitate the mass loading of cold phase atmospheric gasat the base of a magnetic flux tube, perhaps in concertwith a combination of line driving and magnetic pressureforces. However, solving the non-adiabatic MHD equa-tions in the BLR is a very computationally demandingproblem, even in 2D. The HEP of the BLR is found byrearranging Eq. (4),HEP = 3 . × µ (1 − Γ)( T / K) ( r / − . (18)For M bh = 5 . × M (cid:12) and a typical BLR radiusof 10 ld (characteristic of NGC 5548; see e.g., Krisset al. 2019), HEP ≈ . Recalling from § t cool /t orb ∝ HEP − / , applying the same methods usedhere would incur a timestep at least 10 × smaller thanour current one, and our current one is already abouttwo orders of magnitude smaller than the timestep nec-essary to perform global MRI disk simulations usingideal MHD. Nevertheless, we aim to investigate such arole for dynamical TI in the near future. SUMMARY AND CONCLUSIONSPrior work has shown that within ∼ R IC , the pres-ence of a TI zone on the S-curve simply aids Comp-tonization by providing an efficient heating mechanism(see § § Ambiguity arises because the midplane BC serves to approx-imate the angular momentum supply of the underlying disk butnot its temperature, which is significantly colder than the actualtemperature assigned to θ = π/
2. Optical depth effects in thetransition from disk to atmosphere will dictate the value of thispressure gradient, so radiation MHD is ultimately needed. the large scale vortices that are present in the atmo-spheres of our solutions within 80 R IC , streamlines inthe atmosphere transition from being connected to thedisk wind to being separately outflowing. On a phase di-agram, this atmospheric outflow occupies the cold phasebranch of the S-curve below the temperature of the TIzone. The only other pathway into the TI zone is alongstreamlines passing through the vortices themselves, andthe circulation permits multiple passes through the zone,giving TI more time to amplify entropy modes. Thesestreamlines are confined to a small range of radii rela-tive to the domain size but nevertheless occupy a regionspanning > R IC (see Fig. 7).The resulting nonlinear perturbations produce hotspots that then seed the formation of an IAF. Namely,IAFs are what we call the large, filamentary and outwardpropagating vortical structures that arise as hot spotsbecome hot bubbles and disrupt the transition layer be-tween the atmosphere and disk wind. They are liftedinto the disk wind by the expansion and buoyancy ofthe bubbles, which undergo continuous runaway heatingas a consequence remaining thermally unstable. Onceexposed to the wind, the filamentary layers of the IAFevaporate and supply the disk wind with multiphase gasover a large solid angle. Occasionally these layers breakinto smaller clumps that also evaporate before enteringthe supersonic flow. Importantly, however, this evap-oration dynamics occurs on timescales of thousands tohundreds of thousands of years in AGNs, so IAFs canprovide a physical explanation for the fact that warmabsorber type outflows are very common.Our simulations show that the IAFs transition frommere bulges and depressions in the atmosphere to largescale filaments right around the ‘unbound radius’, R u ≈ . × T C / K T c , max / K R IC (1 − Γ) , (19)which was identified as the distance beyond which ourradial wind solutions from D20 can become clumpy. Wenote that R IC (1 − Γ) is the generalized Compton radiusupon accounting for radiation pressure (Proga & Kall-man 2002) and that R u (like R IC ) is a property of theS-curve. Therefore, the theory of multiphase thermallydriven disk winds is predictive: this radius is sensitiveto both the luminosity and the SED, making it possi-ble to compare R u with the inferred locations to warmabsorbers across a population of AGNs with differentvalues of Γ and T C /T c , max .In conclusion, here we have uncovered what are likelythe simplest type of multiphase disk wind solutions. Itis also noteworthy that this year marks 25 years sincethe pioneering numerical work by Woods et al. (1996), ultiphase AGN disk winds ACKNOWLEDGEMENTSWe thank Sergei Dyda for regular discussions overthe course of this project. TW thanks Hui Li and theTheoretical Division at Los Alamos National Labora-tory (LANL) for allowing him to retain access to the in-stitutional computing (IC) clusters. These calculationswere performed under the LANL IC allocation award w19_rhdccasims . Support for this work was providedby the National Aeronautics Space Administration un-der ATP grant NNX14AK44G and through ChandraAward Number TM0-21003X issued by the Chandra X-ray Observatory Center, which is operated by the Smith-sonian Astrophysical Observatory for and on behalf ofthe National Aeronautics Space Administration undercontract NAS8-03060.REFERENCES
Adhikari, T. P., R´o˙za´nska, A., Hryniewicz, K., Czerny, B.,& Behar, E. 2019, ApJ, 881, 78,doi: 10.3847/1538-4357/ab2dfcArav, N., Chamberlain, C., Kriss, G. A., et al. 2015, A&A,577, A37, doi: 10.1051/0004-6361/201425302Balbus, S. A. 1986, ApJL, 303, L79, doi: 10.1086/184657Balbus, S. A., & Soker, N. 1989, ApJ, 341, 611,doi: 10.1086/167521Banda-Barrag´an, W. E., Zertuche, F. J., Federrath, C.,et al. 2019, MNRAS, 486, 4526,doi: 10.1093/mnras/stz1040Bautista, M. A., & Kallman, T. R. 2001, ApJS, 134, 139,doi: 10.1086/320363Begelman, M. C., McKee, C. F., & Shields, G. A. 1983,ApJ, 271, 70, doi: 10.1086/161178Behar, E. 2009, ApJ, 703, 1346,doi: 10.1088/0004-637X/703/2/1346Beltrametti, M. 1981, ApJ, 250, 18, doi: 10.1086/159344Bianchi, S., Guainazzi, M., Laor, A., Stern, J., & Behar, E.2019, MNRAS, 485, 416, doi: 10.1093/mnras/stz430Bianchi, S., Maiolino, R., & Risaliti, G. 2012, Advances inAstronomy, 2012, 782030, doi: 10.1155/2012/782030Borkar, A., Adhikari, T. P., R´o˙za´nska, A., et al. 2021,MNRAS, 500, 3536, doi: 10.1093/mnras/staa3515Chakravorty, S., Kembhavi, A. K., Elvis, M., & Ferland, G.2009, MNRAS, 393, 83,doi: 10.1111/j.1365-2966.2008.14249.xChakravorty, S., Kembhavi, A. K., Elvis, M., Ferland, G.,& Badnell, N. R. 2008, MNRAS, 384, L24,doi: 10.1111/j.1745-3933.2007.00414.x Chakravorty, S., Misra, R., Elvis, M., Kembhavi, A. K., &Ferland, G. 2012, MNRAS, 422, 637,doi: 10.1111/j.1365-2966.2012.20641.xCowie, L. L., & McKee, C. F. 1977, ApJ, 211, 135,doi: 10.1086/154911Crenshaw, D. M., Kraemer, S. B., & George, I. M. 2003,ARA&A, 41, 117,doi: 10.1146/annurev.astro.41.082801.100328Dannen, R. C., Proga, D., Kallman, T. R., & Waters, T.2019, ApJ, 882, 99, doi: 10.3847/1538-4357/ab340bDannen, R. C., Proga, D., Waters, T., & Dyda, S. 2020,ApJL, 893, L34, doi: 10.3847/2041-8213/ab87a5de Kool, M., & Begelman, M. C. 1995, ApJ, 455, 448,doi: 10.1086/176594Detmers, R. G., Kaastra, J. S., Steenbrugge, K. C., et al.2011, A&A, 534, A38, doi: 10.1051/0004-6361/201116899Done, C., Davis, S. W., Jin, C., Blaes, O., & Ward, M.2012, MNRAS, 420, 1848,doi: 10.1111/j.1365-2966.2011.19779.xDyda, S., Dannen, R., Waters, T., & Proga, D. 2017,MNRAS, 467, 4161, doi: 10.1093/mnras/stx406Ebrero, J., Kaastra, J. S., Kriss, G. A., et al. 2016, A&A,587, A129, doi: 10.1051/0004-6361/201527808Elvis, M. 2017, ApJ, 847, 56,doi: 10.3847/1538-4357/aa82b6Emmering, R. T., Blandford, R. D., & Shlosman, I. 1992,ApJ, 385, 460, doi: 10.1086/170955Ferland, G. J., Done, C., Jin, C., Landt, H., & Ward, M. J.2020, MNRAS, 494, 5917, doi: 10.1093/mnras/staa1207Ferland, G. J., Kisielius, R., Keenan, F. P., et al. 2013,ApJ, 767, 123, doi: 10.1088/0004-637X/767/2/123 Waters, Proga, & Dannen
Ferland, G. J., Chatzikos, M., Guzm´an, F., et al. 2017,RMxAA, 53, 385. https://arxiv.org/abs/1705.10877Field, G. B. 1965, ApJ, 142, 531, doi: 10.1086/148317Gaspari, M., Brighenti, F., & Temi, P. 2015, A&A, 579,A62, doi: 10.1051/0004-6361/201526151Giustini, M., & Proga, D. 2019, A&A, 630, A94,doi: 10.1051/0004-6361/201833810Goosmann, R. W., Holczer, T., Mouchet, M., et al. 2016,A&A, 589, A76, doi: 10.1051/0004-6361/201425199Grafton-Waters, S., Branduardi-Raymont, G., Mehdipour,M., et al. 2020, A&A, 633, A62,doi: 10.1051/0004-6361/201935815Higginbottom, N., Knigge, C., Long, K. S., et al. 2018,MNRAS, 479, 3651, doi: 10.1093/mnras/sty1599Higginbottom, N., & Proga, D. 2015, ApJ, 807, 107,doi: 10.1088/0004-637X/807/1/107Higginbottom, N., Proga, D., Knigge, C., & Long, K. S.2017, ApJ, 836, 42, doi: 10.3847/1538-4357/836/1/42Holczer, T., Behar, E., & Kaspi, S. 2007, ApJ, 663, 799,doi: 10.1086/518416Jacquemin-Ide, J., Lesur, G., & Ferreira, J. 2020, arXive-prints, arXiv:2011.14782.https://arxiv.org/abs/2011.14782Jin, C., Ward, M., & Done, C. 2012, MNRAS, 422, 3268,doi: 10.1111/j.1365-2966.2012.20847.xKaastra, J. S., Steenbrugge, K. C., Raassen, A. J. J., et al.2002, A&A, 386, 427, doi: 10.1051/0004-6361:20020235Kaastra, J. S., Kriss, G. A., Cappi, M., et al. 2014, Science,345, 64, doi: 10.1126/science.1253787Kallman, T., & Dorodnitsyn, A. 2019, ApJ, 884, 111,doi: 10.3847/1538-4357/ab40aaKallman, T. R. 2010, SSRv, 157, 177,doi: 10.1007/s11214-010-9711-6Kallman, T. R., & McCray, R. 1982, ApJS, 50, 263,doi: 10.1086/190828Kriss, G. A., De Rosa, G., Ely, J., et al. 2019, ApJ, 881,153, doi: 10.3847/1538-4357/ab3049Krolik, J. H., & Kriss, G. A. 2001, ApJ, 561, 684,doi: 10.1086/323442Krolik, J. H., McKee, C. F., & Tarter, C. B. 1981, ApJ,249, 422, doi: 10.1086/159303Krolik, J. H., & Vrtilek, J. M. 1984, ApJ, 279, 521,doi: 10.1086/161916Kurosawa, R., & Proga, D. 2009, ApJ, 693, 1929,doi: 10.1088/0004-637X/693/2/1929Laha, S., Guainazzi, M., Dewangan, G. C., Chakravorty, S.,& Kembhavi, A. K. 2014, MNRAS, 441, 2613,doi: 10.1093/mnras/stu669Laha, S., Reynolds, C. S., Reeves, J., et al. 2020, NatureAstronomy, doi: 10.1038/s41550-020-01255-2 Lee, J. C., Kriss, G. A., Chakravorty, S., et al. 2013,MNRAS, 430, 2650, doi: 10.1093/mnras/stt050Leighly, K. M., & Casebeer, D. 2007, in AstronomicalSociety of the Pacific Conference Series, Vol. 373, TheCentral Engine of Active Galactic Nuclei, ed. L. C. Ho &J. W. Wang, 365Lepp, S., McCray, R., Shull, J. M., Woods, D. T., &Kallman, T. 1985, ApJ, 288, 58, doi: 10.1086/162763Lodders, K., Palme, H., & Gail, H. P. 2009, LandoltBörnstein, 4B, 712,doi: 10.1007/978-3-540-88055-4 34Luketic, S., Proga, D., Kallman, T. R., Raymond, J. C., &Miller, J. M. 2010, ApJ, 719, 515,doi: 10.1088/0004-637X/719/1/515Mathews, W. G., & Doane, J. S. 1990, ApJ, 352, 423,doi: 10.1086/168548Mathews, W. G., & Ferland, G. J. 1987, ApJ, 323, 456,doi: 10.1086/165843Matthews, J. H., Knigge, C., Higginbottom, N., et al. 2020,MNRAS, 492, 5540, doi: 10.1093/mnras/staa136Mehdipour, M., Kaastra, J. S., Kriss, G. A., et al. 2015,A&A, 575, A22, doi: 10.1051/0004-6361/201425373Mizumoto, M., Done, C., Tomaru, R., & Edwards, I. 2019,MNRAS, 489, 1152, doi: 10.1093/mnras/stz2225Mo´scibrodzka, M., & Proga, D. 2013, ApJ, 767, 156,doi: 10.1088/0004-637X/767/2/156Nakatani, R., & Yoshida, N. 2019, ApJ, 883, 127,doi: 10.3847/1538-4357/ab380aPonti, G., Fender, R. P., Begelman, M. C., et al. 2012,MNRAS, 422, L11, doi: 10.1111/j.1745-3933.2012.01224.xProga, D., & Kallman, T. R. 2002, ApJ, 565, 455,doi: 10.1086/324534Proga, D., Stone, J. M., & Drew, J. E. 1998, MNRAS, 295,595, doi: 10.1046/j.1365-8711.1998.01337.xProga, D., & Waters, T. 2015, ApJ, 804, 137,doi: 10.1088/0004-637X/804/2/137Rodr´ıguez Hidalgo, P., Khatri, A. M., Hall, P. B., et al.2020, ApJ, 896, 151, doi: 10.3847/1538-4357/ab9198R´o˙za´nska, A., Czerny, B., Kunneriath, D., et al. 2014,MNRAS, 445, 4385, doi: 10.1093/mnras/stu2066Sharma, P., McCourt, M., Quataert, E., & Parrish, I. J.2012, MNRAS, 420, 3174,doi: 10.1111/j.1365-2966.2011.20246.xSheikhnezami, S., & Fendt, C. 2018, ApJ, 861, 11,doi: 10.3847/1538-4357/aac5dcShlosman, I., Vitello, P. A., & Shaviv, G. 1985, ApJ, 294,96, doi: 10.1086/163278 ultiphase AGN disk winds Spruit, H. C. 1996, in NATO Advanced Study Institute(ASI) Series C, Vol. 477, Evolutionary Processes inBinary Stars, ed. R. A. M. J. Wijers, M. B. Davies, &C. A. Tout, 249–286Stern, J., Behar, E., Laor, A., Baskin, A., & Holczer, T.2014, MNRAS, 445, 3011, doi: 10.1093/mnras/stu1960Stone, J. M., Tomida, K., White, C. J., & Felker, K. G.2020, ApJS, 249, 4, doi: 10.3847/1538-4365/ab929bTomaru, R., Done, C., Odaka, H., Watanabe, S., &Takahashi, T. 2018, MNRAS, 476, 1776,doi: 10.1093/mnras/sty336Tomaru, R., Done, C., Ohsuga, K., Nomura, M., &Takahashi, T. 2019, MNRAS, 490, 3098,doi: 10.1093/mnras/stz2738Townsend, R. H. D. 2009, ApJS, 181, 391,doi: 10.1088/0067-0049/181/2/391 Wang, J.-M., Du, P., Brotherton, M. S., et al. 2017, NatureAstronomy, 1, 775, doi: 10.1038/s41550-017-0264-4Waters, T., & Li, H. 2019, arXiv e-prints,arXiv:1912.03382. https://arxiv.org/abs/1912.03382Waters, T., & Proga, D. 2018, MNRAS, 481, 2628,doi: 10.1093/mnras/sty2398—. 2019, ApJ, 875, 158, doi: 10.3847/1538-4357/ab10e1Waters, T., Proga, D., Dannen, R., & Kallman, T. R. 2017,MNRAS, 467, 3160, doi: 10.1093/mnras/stx238Woods, D. T., Klein, R. I., Castor, J. I., McKee, C. F., &Bell, J. B. 1996, ApJ, 461, 767, doi: 10.1086/177101Zhu, Z., & Stone, J. M. 2018, ApJ, 857, 34,doi: 10.3847/1538-4357/aaafc9 Waters, Proga, & Dannen
APPENDIX A. NUMERICAL METHODSUsing
Athena ++ (Stone et al. 2020), we solve the equations of non-adiabatic gas dynamics, accounting for the forcesof gravity and radiation pressure due to electron scattering opacity: DρDt = − ρ ∇ · v , (A1) ρ D v Dt = −∇ p − G M bh ρr (1 − Γ) , (A2) ρ D E Dt = − p ∇ · v − ρ L . (A3)For an ideal gas, p = ( γ − ρ E , and we set γ = 5 /
3. Our calculations are performed in spherical polar coordinates,( r, θ, φ ), assuming axial symmetry about a rotational axis at θ = 0 ◦ . The S-curve of our net cooling function L isshown in Fig. 2 and was obtained by running a large grid of photoionization calculations using XSTAR (see Dyda et al.2017) for the SED of NGC 5548 obtained by Mehdipour et al. (2015). The abundances are those of Lodders et al.(2009), but for simplicity, in this work we set µ H = 1, keeping only µ = 0 . ρ L = n ( C − H ), and we interpolate from tables of the total cooling ( C )and total heating ( H ) rates given in units of erg cm s − ; these tables were made publicly available by Dannen et al.(2020). To obtain the clumpy outflow solutions in D20, we used bilinear interpolation to access values from the tables(for details, see Dyda et al. 2017). We suspected that this basic algorithm could be somewhat inaccurate near theS-curve where rates are close to zero. In the testing phase of this work, we verified this suspicion, confirming thatinterpolation errors introduce perturbations into the flow that are continually amplified by TI, allowing even our 1Dsolutions to remain clumpy indefinitely. Specifically, the %-difference in rates evaluated near the S-curve where C ≈ H could exceed 1000% based on a comparison of looking up rates tabulated from an analytic function (using a similarsampling as our actual tables).We therefore implemented a GSL bicubic interpolation routine that we found reduces the percent differences forrates near the S-curve to less than 1% and returns rates that are on average about 10 × more accurate. As discussedin § Initial conditions, boundary conditions, and computational domain
Our initial conditions are ρ ( r, θ max ) = ρ ( r /r ) (with ρ = ¯ m n ), v ( r, θ max ) = v φ ( r ) ˆ φ , and p ( r, θ max ) = n k T ( r /r ) , where θ max denotes zones with cell edges at π/ je zones in Athena ++ notation). All other zones have ρ ( r, θ ) = ρ ( T C /T )( r /r ) , v ( r, θ ) = (1 . v (cid:112) . − r in /r + 10 − v )ˆ r , and p ( r ) = p ( r, θ max ), where v = (cid:112) GM bh /r .Reflecting BCs are applied at θ = π/
2. In our fiducial runs, the rotation profile along this boundary is enforced to besub-Keplerian to account for the radiation force, v φ ( r ) = (cid:112) GM bh (1 − Γ) /r at θ = π/
2. We also ran tests accountingfor the contribution of the pressure gradient in the radial force balance along the midplane, given by v φ r = GM bh r (1 − Γ) + 1 ρ dpdr . (A4)For the isothermal midplane BC that we use, dp/dr = ( − /r ) p ( r ) at θ = π/
2. This midplane BC is applied insideour heating and cooling routine (implemented using
EnrollUserExplicitSourceFunction() ). We reset the midplanedensity, velocity, and pressure on every call (2 × per timestep for a 2nd order integrator) to the above values. Onlythe density and velocity are required to be reset, as Athena ++ will find the same solution without specifying thepressure also, but only if a timestep 2 × smaller is used. This stability constraint arises because the midplane containsthe densest gas and thus gas with the shortest cooling times in the domain — see Fig. 3.We apply a ‘constant gradient’ BC at r in , in which all primitive variables are linearly extrapolated from the firstactive zone into the ghost zones. At r out , we apply modified outflow BCs that prevent inflows from developing at theboundary. Specifically, v r is set to 0 in the ghost zones if it is found to be less than zero (in Athena ++ notation, ultiphase AGN disk winds prim ( IVX , k , j , ie + i ) = if prim ( IVX , k , j , ie + i ) < ). Without this latter BC, a transient disturbance willenter from the outer boundary after hundreds of orbital times at R IC (but before IAFs develop).Our calculations utilize static mesh refinement (SMR) to give an effective grid size, N r × N θ ( nx1 × nx2 in Athena ++ notation), of 1056 ×
576 for our mid-res run and 2112 × ×
72, and SMRlevels are set at 15 ◦ < θ < ◦ , 30 ◦ < θ < ◦ , and 35 ◦ < θ < ◦ for our mid-res run. For our hi-res run, a fourthSMR level is applied at 35 ◦ < θ < ◦ . Our radial domain size is r in = 13 . R IC and r out = 285 R IC .A.2. Timestep constraint
In Waters & Proga (2018), we used a semi-implicit solver to include the heating and cooling source term as describedby Dyda et al. (2017). For this work, we had to abandon this routine in favor of a simpler explicit solve after findingthat the semi-implicit solve causes a small fraction of the coldest zones within the atmosphere to ‘jump off’ the S-curve at every timestep when these zones are visualized on the ( T, Ξ)-plane. This behavior indicates that numerically,multiple temperature solutions can exist for the same timestep, ∆ t . This is a generic property of implicit routines(Townsend 2009), with the simple remedy being to resort to an explicit solve. Note that Townsend’s exact integrationroutine requires L to be a function of temperature alone, whereas here, L = L ( T, ξ ).An explicit solve requires implementing a custom timestep (using