A Self-consistent Treatment of the Line-driving Radiation Force for Active Galactic Nuclei Outflows: New Prescriptions for Simulations

Flows driven by photons have been studied for almost a century, and a quantitative description of the radiative forces on atoms and ions is important for understanding a wide variety of systems with outflows and accretion disks, such as active galactic nuclei (AGN). Quantifying the associated forces is crucial to determining how these outflows enable interactive mechanisms within these environments, such as AGN feedback. The total number of spectral lines in any given ion of the outflow material must be tabulated in order to give a complete characterization of this force. Here, we provide calculations of the dimensionless line force multiplier for AGN environments. For a wide array of representative AGN sources, we explicitly calculate the photoionization balance at the proposed wind-launching region above the accretion disk, compute the strength of the line-driving force on the gas, and revisit and formalize the role of the commonly used ionization parameter ξ in ultimately determining the line-driving force. We perform these computations and analyses for a variety of AGN central source properties, such as black hole mass, initial wind velocity, and number density. We find that, while useful, the ionization parameter provides an incomplete description of the overall ionization state of the outflow material. We use these findings to provide an updated method for calculating the strength of the radiative line-driving using both the X-ray spectral index Γ X and the ionization parameter.


INTRODUCTION
Active Galactic Nuclei (AGN) are believed to be powered by accretion onto a supermassive black hole (SMBH).Quantifying and understanding the associated mass and energy flows are crucial to understanding the observed properties of quasars, blazars, Seyfert, and radio galaxies, as well as how these objects evolve and interact with their surroundings.These powerful outflows are suspected of being capable of influencing star formation rates far beyond the small-scale size of the AGN itself.This long-range "feedback" is likely responsible for producing several known correlations between the SMBH mass and a range of large-scale properties of the host galaxy and intergalactic medium (Magorrian et al. 1998;Silk & Rees 1998;Cavaliere et al. 2002;Ostriker et al. 2010;Fabian 2012;King & Pounds 2015;Hopkins et al. 2016;Harrison et al. 2018).
Corresponding author: Aylecia S. Lattimer As seen in Figure 1, various outflow components have been proposed to be present in many AGN (see also Crenshaw et al. 2003;Beckmann & Shrader 2012;Laha et al. 2021).The large number of observed outflow characteristics has resulted in a proliferation of taxonomic AGN types.Along with the development of these AGN categories came many attempts to build unified models that explain the wind types not as products of fundamental structural differences, but as a result of, e.g., viewing angle or evolutionary stage (Antonucci 1993;Elvis 2000;Netzer 2015;Padovani et al. 2017).For example, features associated with the broad-line region (BLR) are believed to sit at lower latitudes than the larger-scale distribution of features associated with the narrow-line region (NLR).Energetic outflows in those regions are observed in UV lines as broad absorption line (BAL) or narrow absorption line (NAL) clouds, with features of an intermediate nature called "mini-BAL."X-ray spectra indicate the presence of a dynamic coronal line region (CLR), some outflow regions that are highly ionized and escaping at nearly relativistic speeds (e.g., ultra-fast outflows; hereafter UFOs), and some with lower velocities and a range of ionization states (i.e., the elusive "warm absorber").The exhibited overlap in the measured properties of these features suggest that it is possible some could arise from the same highly structured and stratified structures being viewed with different spectral diagnostics.
Multiple explanations for AGN wind acceleration have been proposed.For example, thermal forces, magnetocentrifugal, and other magnetohydrodynamic (MHD) effects have all shown promise in being able to explain various key observations.It is likely that some combinations of effects may be responsible for accelerating some features.Although the physical processes responsible for accelerating AGN winds have not been conclusively identified, another possible mechanism driving for these outflows is radiative driving (Giustini & Proga 2019).This is often referred to as radiation pressure, where the force of radiation acts on the spectral lines in the material of the wind or outflow.The number of spectral lines in any given ion have a dominant effect on this force, which can act on the plasma as a whole, or on subsets of ions (Castor 1974;Castor et al. 1975, hereafter CAK).The subsequent absorption and re-emission of photons in a given spectral line results in a net radial transfer of momentum, giving rise to a line-driven wind.
We now know that photon-driven outflows may play a role in various other astrophysical environments, such as massive OB stars, as well as other accreting objects such as protostars and cataclysmic variables.This type of driving is also often thought to occur in the BAL region of AGN, where the blueshift of the spectral lines suggest that outflow velocities can reach up to 0.2c.Line-driven disk winds are a promising hydrodynamical scenario for these outflows, due in part to the large number of strong resonance lines present in the spectra of many broad absorption line quasars (BALQSOs; Murray et al. 1995;Proga et al. 2000;Proga 2007;Risaliti & Elvis 2010;Higginbottom et al. 2014;Nomura et al. 2016;Zhu et al. 2022).
However, in order for line-driving to launch and maintain an outflow from the accretion disk, there must exist a sufficient number of bound atoms to produce those lines.Thus, the gas must be relatively lowly ionized.This presents a challenge to the efficiency of linedriving for these winds, since the material is question is in close proximity to the central source and can easily become over-ionized (Higginbottom et al. 2014).A proposed solution to this issue is a so-called "failed wind" near the interior of the accretion disk, which acts to shield the outer material from becoming over-ionized (see, e.g., Proga et al. 2000;Proga & Kallman 2004;Quera-Bofarull et al. 2020).
The main parameter used to characterize these outflows is the line-force multiplier (discussed in detail in Section 4.3), which describes both the overall strength of the line-driving as well as how it varies in different acceleration regimes of an outflow.This work aims to produce improved models of these line-driving forces by building upon previous force multiplier calculations for AGN and other high-energy systems containing pho-toionized outflows (Stevens & Kallman 1990;Stevens 1991;Arav & Li 1994;Dannen et al. 2019;Chelouche & Netzer 2003;Everett 2005;Chartas et al. 2009;Saez & Chartas 2011).In this paper, we adapt the process outlined for massive stellar sources in Lattimer & Cranmer (2021) to calculate the force multiplier for winds launched from AGN disks, using an even larger database of atomic line-strength parameters.This paper also computes the line-driving force for a large grid of background AGN parameters (i.e., black-hole mass, accretion rate, and the radial location, speed, and acceleration of the outflow), whereas many earlier models limit these calculations to a smaller number of representative cases.Additionally, we explore the importance of the previously mentioned failed wind in launching and maintaining a line-driven disk wind.
In Section 2 we describe the assumed geometry of the system and the initial parameters of the models.Section 3 describes the physics taken into account in the model and the calculation of the thermal equilibrium, ionization balance, and luminosity along the radial line of sight (LOS) from the central source.In Section 4 we discuss the calculation of the weighted strengths of the spectral lines.We end with a discussion of our results and our conclusions in Sections 5 and 6 respectively.

SED of the Central AGN
The purpose of this paper is to describe the variety of possible forms the line-driving force multiplier M (t) may take in realistic AGN environments.To do this, we stop short of attempting to produce self-consistent dynamical models of the AGN outflows themselves.We include only the properties of the gas required to compute the ionization states of the elements producing spectral lines that affect M (t).Thus, the densities, velocities, and temperatures that we specify at various distances from the AGN central source are meant to be representative examples of the conditions likely to be found in these regions, not complete solutions to fluid or MHD conservation equations.
For our models, we consider an AGN central source embedded in an accretion disk.The initial SED of the central source is computed using a modified version of the python module QSOSED (Quera-Bofarull 2019).This module recreates the qsosed model of XSPEC (Arnaud 1996), following the method outlined in Kubota & Done (2018).These model SEDs are made up of three characteristic regions: the outer cool disk, the warm Comptonising inner disk, and the hot inner Comptonising region (i.e. the X-ray corona).For the warm Comptonising region, QSOSED follows the passive disk sce- nario as discussed in Petrucci et al. (2018); Kubota & Done (2018).The general geometry followed by both QSOSED and our models is shown in Figure 2.
The two main free parameters used for computing the SED are the mass of the black hole M BH and the black hole mass accretion rate in Eddington units ṁ (Beckmann & Shrader 2012), given by ṁ =

ṀEdd
. (1) Additionally, for all models we assume a black hole dimensionless spin absolute value a * of zero, meaning that all black holes in our models are non-rotating.This also determines that the radius r isco of the innermost stable circular orbit (ISCO) will be constant across all models when given in units of gravitational radii R g , defined as For a given set of AGN properties, we compute the line force multiplier M (t) along a radial ray of points that starts at the radius of the corona R cor and extends to some outermost distance r out .This ray allows us to compute the line force at many possible locations where there may be a line-driven wind.This is also used to compute the absorption of the SED due to material interior to each radial distance r along it.For the remainder of this work, we define r out as 3 kpc for each model.This was chosen to ensure that the outermost radial step of the model exceeds the largest distances reported for observed AGN outflow features (see, e.g.Laha et al. 2021).For completeness, we additionally define the outer edge of the disk to be the self-gravity radius R self as given in Laor & Netzer (1989).Figure 2 shows a generalized diagram of the modeled geometry.We note that although Figure 2 includes the curvature of the flux-tube for the sake of completeness, this curvature is not taken into account by the mass conservation equation used in this work.
The luminosity of the central source at the first radial step (R cor ), is calculated from the SED flux F ν that is produced by QSOSED: Likewise, the initial bolometric luminosity is also defined by QSOSED at R cor .Figure 3 shows the calculated luminosity as produced by QSOSED at R cor for a range of black hole masses and accretion rates.We see a linear correlation between the luminosity L ν and M BH .Thus, as seen in Kubota & Done (2018), we can set the initial SED using solely the black hole mass M BH and the accretion rate ṁ.We also see that the SED at X-ray frequencies depends strongly on ṁ, flattening significantly with decreasing accretion rate.Similarly, Kubota & Done (2018) found that the hot Comptonization region must truncate at the corona (R cor ) for low accretion rates ( ṁ ≲ 0.2).This is due to lower limits on the X-ray spectral index Γ X (Haardt & Maraschi 1991;Malzac et al. 2005;Stern et al. 1995), which describes the slope of the X-ray portion of the SED (see Section 6 for further discussion of Γ X ).Thus, we assume for all models that the disk does not extend inward beyond the radius of the hot corona such that the disk begins at R cor and extends to r out .This subsequently allows us to use R cor as the first radial point along the primary LOS.
Figure 4 shows the dependence of both R self and R cor on the mass accretion rate onto the black hole ṁ.From Equation (18) of Laor & Netzer (1989), we expect that the self-gravity of the disk will be lower for large M BH , due to the strong gravity of the central source.Conversely, for large values of ṁ, the larger amount of gas in the disk should increase the self-gravity, pushing R self outward.This is confirmed in Figure 4, where we see that R cor starts to exceed R self for Eddington accretion rates below a few times 10 −2 , putting a lower limit on the values of ṁ that we use in these models.

Gas Outflow Properties & Monte Carlo Process
To create our grid of models, we construct ranges of reasonable values for our initial parameters.The first of these are the black hole mass and the accretion rate as described above, which set the properties of the source SED.Additional input parameters include launch radius of the wind r L , number density exponent β, a shielding factor S, the wind fraction f wind , and the outflow velocity u wind .Grids of possible input values for each initial parameter are given in Table 1, and we describe each of these parameters in more detail below.
We must also define the radial grid that represents the line of sight between the parcel of gas and the central AGN.Beginning at R cor , each model steps radially through the line of sight of the gas parcel to the final observing radius r out from the central source.We found that 500 points along this line of sight were sufficient to achieve convergence in our final calculations of the force multiplier.Taking into account the truncation of the disk at R cor discussed above, and the fact that the calculation of the initial SED includes the SED of the corona itself, we therefore use the corona radius as the first radial step along the LOS.This precludes the ability to begin the radial grid inside the corona, at, for example, r isco .The rapid attenuation of the SED with increasing radius requires that the distance between each radial step be smaller at distances close to R cor .Because the radial grid must by necessity span several orders of magnitude, we space our steps by first creating a condensed grid, defined by bounds x 1 and x 2 , given by (5) We then space steps evenly in this condensed grid.To construct the final set of radial points to be used in our models, we convert these linearly spaced steps to fit the original bounds of our radial range (i.e.R cor and r out ).For example, for a given value x 1 < x < x 2 in the condensed grid, the radius in the final radial grid along the primary LOS would be given by log(r) = x 2 .This method preferentially clusters the radial points of the final grid near R cor , allowing us to reach consistent results with far fewer radial steps.For this work, we determined the number of radial points necessary to reach converged results by calculating results for a fixed representative subset of models using 100, 250, 500, 1000, and 1500 radial points.Using the method described above, we reached consistent results, where results from a given number of radial points matched those from the highest number of points, using only 500 radial steps.This drastically reduces the computation time required for a single model.Therefore, for each model we used a grid of 500 radial steps, beginning at R cor and ending at r out .Motivated by the assumption that AGN outflows are launched at the surface of the accretion disk and often become reoriented in the radial direction (see, e.g.Proga  & Kallman 2004;Quera-Bofarull et al. 2023), we define a so-called launch radius r L that falls somewhere between 100 and 1,000 R g for each model.We can then specify the hydrogen number density and outflow speed differently inside and outside this radius using Here, S is a dimensionless shielding enhancement factor introduced to vary the number density before the launch radius, in order to probe the importance of a proposed inner failed wind (see, e.g.Murray et al. 1995;Chelouche & Netzer 2003).For this work, log(S) is sampled uniformly between 0 and 2. Note that at radii less than the launch radius, we assume a more-or-less static atmosphere with a large enough scale height to safely assume a constant density.At radii greater than the launch radius, our radial grid intercepts a presumed outflow flux-tube that has a representative radial speed u wind and a density that decreases with increasing distance.The exponent β is varied between values of 1.01 and 4. For a time-steady spherically symmetric outflow with a constant outflow speed, mass-flux conservation would give β = 2. Larger values would be consistent with a radially increasing velocity, so β can be considered as a parameterization of the local degree of acceleration in the modeled outflow.Conversely, smaller values tend toward the commonlyused assumption of a constant representative density at different distances (see, e.g.Dannen et al. 2019).However, we do not use values β ≤ 1 as these values would lead to an infinitely large column density (see Equation ( 11)). Figure 5 shows representative versions of n H and u r with increasing radius.
The number density at the launch radius n L is given by mass-flux conservation: where m H is the mass of a hydrogen atom, and Ṁwind is a spherically-symmetric approximation of the wind's mass-loss rate.This is given as some fraction f wind of the black hole's mass accretion rate Ṁacc (see also Daly 2021): We constrain this f wind "efficiency factor" to fall within the range of 10 −2 and 10 (see, e.g.King 2010; King & Pounds 2015;Mou et al. 2017;Yi et al. 2017;Leighly et al. 2018).
In Equation ( 9), u wind is the outflow velocity, which we set at a constant value for all points along the line of sight, beginning at the launch radius, as in Equation ( 8).The wind velocity for each model is chosen from a range of 100 to 70,000 km s −1 .This range ensures that our models encompass the full range of possible outflow velocities.The lower limit of 100 km s −1 is consistent with low-velocity NAL outflows, while the upper limit of 70, 000 km s −1 represents high velocity BALs and mini-BALs (Giustini & Proga 2019).Because we dictate that the launch radius of the wind fall between 100 and 1,000 R g , we do not consider velocities that are typical of UFOs (∼ 0.4c), as these are often launched below this range, at ∼ 10 R g .
From the hydrogen number density n H (Equation ( 7)), the hydrogen column density N H at each radial step can also be found: Note that in Equation ( 11) we adopt the convention of calculating the column density as the intervening column density between the parcel of gas at its radially dependent position and an observer at r = ∞.See Section 3.3 for an in-depth description of how the SED absorption is computed along the ray.Additionally, we calculate the effective hydrogenionizing luminosity, L X , of the current SED, given by  7)) and radial flow speed ur (Equation ( 8)).The radius of the hot corona and launch radius of the wind are denoted by Rcor and rL respectively.The number density at the launch radius, nL, is also marked.
Here, J ν is the mean intensity for a given frequency.It is useful to define L X in terms of J ν , as we eventually use the mean intensity to account for the full radiation field seen by the parcel of gas at a given radius (see Section 3.3).The frequency limits ν min and ν max are defined as the frequencies corresponding to 1 and 1,000 Rydbergs, respectively (Fukumura et al. 2022).It is worth noting that L X is sometimes given as the integral over all energies above 1 Rydberg (13.6 eV).However, in this work we find the difference between this and the chosen limits of 1 to 1,000 Rydbergs to have a negligible difference on the values of L X .
From L X we can calculate the ionization parameter ξ at each radial step, which is defined as as originally given by Tarter et al. (1969).The ionization parameter is often used to describe the full ionization balance of the wind.However, it may fail to take into account the full ionization and recombination processes that may be present (see, e.g Shields 1974;Stevens 1991;Osterbrock & Ferland 2006;Krtička et al. 2022).

Thermal Equilibrium
In order to compute the ionization balance of each relevant element, we require values of the local temperature T at each point along the radial ray.To compute T , we assume that the parcel is in time-steady thermal equilibrium.Dyda et al. (2017) found that this equilibrium consists essentially of radiative heating/cooling and adiabatic cooling.This can be specified as The radiative term is specified by an analytic fit that contains terms for bremsstrahlung, Compton, X-ray, and UV-line processes.We define The quantity Λ − Γ is given by Equation 8of Dyda et al. (2017) as a function of T and the X-ray temperature T X , which we compute from the radially dependent SED.The X-ray temperature T X can be calculated from the mean photon energy: The integration limits ν min and ν max used here are the frequencies that correspond to the traditionally accepted energy bounds on what are considered X-rays, at 0.1 and 100 keV respectively.These are chosen to represent the full range of wavelengths that fall in the X-ray portion of the spectrum when calculating the X-ray temperature.
We note here that these limits are used solely in the calculation of T X .The entire frequency range of the source SED is used when calculating the other quantities (e.g., ξ) that are used in determining the temperature used at each radial step.Dyda et al. (2017) used a fixed X-ray temperature of T X = 1.2 × 10 8 K.We see from Figure 6, which shows the calculated values of T X along each step of the radial grid for our set of ∼ 550 models, that many of our models at some point along the primary LOS fall near this value.The number and criteria of what we consider complete models are discussed in further detail later in this section.We see a clear anti-correlation between ṁ and T X at the coronal base radius R cor .This is a consequence of the steepening in the SED at X-ray frequencies seen in Figure 3(b), with steeper spectra (i.e., correspondingly lower values of T X ) occurring for larger accretion rates.We also see that as radial distance increases, T X frequently increases.The changes in SED shape that are responsible for this variation occur due to frequency-dependent absorption, and are discussed further below in Section 3.3.
Rearranging Equation ( 8) of Dyda et al. (2017), we can define Γ − Λ as As in Dyda et al. (2017), the terms of Equation ( 17) represent the heating from Compton processes (A c ) and X-rays (A x ) and cooling from bremsstrahlung (A b ) and spectral line emission (A l ).Taking all quantities to be in cgs units, we set A c = 8.9 × 10 −36 , A x = 1.5 × 10 −21 , A b = 1.3 × 10 −26 , A l = 1.7 × 10 −18 , and the line temperature T l = 1.3 × 10 5 K (see Blondin 1994;Dyda et al. 2017).The adiabatic cooling term, which is present only when there is a nonzero outflow speed, was specified by Cranmer et al. (2007) for flux-tube streamlines with cross-sectional area A(r) as where u is the outflow speed along the flux tube and P is the gas pressure.If we assume power-law radial dependencies for these quantities, we can express everything in terms of the local quantities and not their derivatives.Thus, we assume T ∝ r −δ and ρ ∝ r −β .Combining this with mass-flux conservation along the expanding tube (i.e., ρuA = constant), we get We set δ = 0 for our model outflows.Setting δ = 0 assumes the temperature varies more slowly with distance than the density, and it also ensures that Q ad is always negative, making it a true cooling term.Finally, for an ionized, hydrogen-dominated plasma we use the ideal gas law to specify P ≈ 2n H k B T .
Combining Equations ( 15), ( 17), and ( 19), solving for T becomes a root-finding exercise, which for this work was done using the optimize module available in SciPy (Virtanen et al. 2020).We also implement a temperature floor T floor , which we set at T floor = 1000 K.In the cases where the algorithm fails to find a root, the model is terminated and marked as a failed run.In order to achieve an acceptable level of statistical significance, we initiated 576 models.After accounting for these failed runs, likely resulting from randomly occurring unphysical combinations of initial parameters, we were left with a total of 556 successful modeled outflows.All analyses in the remainder of this work use this total number of model runs.
Figure 7 shows the calculated temperature for each completed model as a function of the ionization parameter ξ.For all successful runs, we hereafter assume that the temperature T found by the thermal equilibrium calculation described above is equal to the electron temperature of the material, such that T e = T .Note that if all of our models had used identical values of T X and ignored adiabatic cooling, the results in Figure 7 would all fall along a single curve, similar to the results of Dannen et al. (2019) that are also displayed.Further quantitative comparison to previous work is given detailed in Appendix B.
However, our models show additional variability that depends on different combinations of the input parameters.Specifically, the variability at small radial distances (corresponding mostly to ξ ≳ 10 6 ) is mostly due to variations in T X at the coronal base, which in turn are driven by different values of ṁ (see Figure 6).The variability among our models at larger radii (r ≳ 10 18 cm, or ξ ≲ 10 5 ) is mostly due to different models exhibiting differing relative amounts of adiabatic cooling.At inner radii, corresponding to high ionization regions beginning at ∼ ξ ≳ 10 3 , our results differ from the two Dannen et al. (2019) curves due to the adiabatic cooling included in our thermal equilibrium calculation (see also Section 5.2).We note also that models with ξ ≲ 10 −6 tend to correspond to plasmas with equilibrium temperatures at or below the imposed floor of 1000 K.In subsequent sections we often neglect these models as not corresponding to "realistic" AGN conditions.

Ionization/Recombination Balance
After obtaining the local temperature T , the ionization/recombination equilibrium of the gas parcel at the current radius can be computed.There are five processes that need to be accounted for in this calculation: collisional ionization (CI i ), autoionization (AI i ), photoionization (P I i ), radiative recombination (RR i+1 ), and dielectric recombination (DI i+1 ).We note that for this work we do not include three-body (collisional) recombination (3BR i+1 ), as this effect generally does not represent a significant contribution to the ionization balance except in the high-density limit, which the AGN environments considered here traditionally do not experience (see, e.g., Rees et al. 1989).The five rates considered are given in units of inverse seconds (s −1 ).For a time steady case, we can write the ionization balance as Here the subscript ion indicates processes that ionize out of a given ionization state ion and into ion + 1, and the subscript oin+1 in the denominator terms indicated processes that recombine from state ion + 1 down to ion.Hence the numerator accounts for ionization processes, and the denominator accounts for recombination.The collisional ionization (CI i ) rate is given by where n e is the electron number density and q c,i (T ) is a temperature dependent collision rate coefficient, with units of [cm 3 s −1 ].Fitting formulas for q c,ion (T ), for all elements up to Ni (Z = 28), were provided by Voronov (1997) and are computed by a routine distributed by Dima Verner of the University of Kentucky.We estimate the coefficients for Z = 29 and Z = 30 according to the method developed in Cranmer (2000).
We can write the autoionization (AI i ) rate as AI ion = n e q a,ion (T ).( 22) Data for the rate coefficient fits to q a,i (T ) were provided in Landini & Monsignori Fossi (1990) for 13 of the most abundant elements.Where available, we use these fits to calculate the autoionization contribution to the overall ionization balance.Where the rate coefficients were not available, i.e. for less abundant elements, the autoionization contribution was set to zero.
The photoionization rate (P I ion ) can be described by where σ ν is the photoionization cross-section and J ν is the mean intensity, as discussed further in Section 3.3.The rate for radiative recombination can be written as Note that here the α rate coefficient is labeled with the stage that is being recombined into.For elements up to Zn (Z = 30), these rates were tabulated by Verner & Ferland (1996) We describe dielectric recombination as We calculate the dielectric recombination rates for elements up to Z = 28 from Mazzotta et al. (1998).However, we neglect dielectric recombination for elements with Z > 28.
Once the overall ionization balance is calculated from Equation (20) for each ion of each element considered, we can then compute the number density n ion for each ion.This is done following the iterative process described in Section 2.1 of Lattimer & Cranmer (2021).At each radial step the initial estimate of n e was given by n e = 0.8n H , where n H is the hydrogen number density of the wind (Equation ( 7)).The initial estimate was then refined using an undercorrection technique until convergence to a final value.This was done at the end of each iteration over the ionization balance by tabulating a new estimate of n e from the calculated ionization balance, which was then multiplied by the previous estimate.The square root of this product was then used as the estimate of n e for the next iteration.We found that 35 iterations were sufficient to reach a stable and self-consistent value.At each radial step we also calculate the fraction n ion /n el , as well as the mean ionization number I, which we define using the ion number i: This provides a uniform method of showing the relative distance between an element being completely neutral (I = 0) and being fully ionized (I = 1).We also provide ionic fractions n ion /n el for selected elements in Appendix B, as well as a comparison of these values to those obtained from other photoionization codes.
Figure 8 shows the mean ionization number I, varying with the ionization parameter ξ for each model.Traditionally, ξ has been used to describe the total ionization state of the outflow material.However, we see from Figure 8 that multiple values of I are possible for the same value of the ionization parameter, suggesting that ξ does not provide an adequately detailed and thus unique description of the outflow's ionization balance (see also, e.g.Kallman & McCray 1982;Stevens 1991;Devereux & Heaton 2013;Dannen et al. 2019;Kallman et al. 2021).We discuss the implications of this further in Section 6.

Luminosity & Flux Calculation
At each radial step along the line of sight, the local temperature of the gas T (Section 3.1), the ionization/recombination balance (Section 3.2), and the luminosity L ν of the source are recalculated to account for the attenuation of the SED due to both the increased distance from the central source and absorption from the intervening material.When calculating the large-scale absorption of the SED, we take into account continuum and bound-free sources of opacity.We also adopt the convention of assuming that the line opacity is involved only in the calculation of the line-driving force (see, e.g., CAK).The frequency-dependent opacity of the material is then given by where σ T is the Thomson scattering cross section.The cross-section for photoelectric absorption σ ν is given here as the photoionization cross section, and is calculated for each ion according to Equation (1) of Verner & Yakovlev (1995).The necessary analytic fitting coefficients are taken from Verner & Yakovlev (1994).σ ν is often equal to zero below a threshold ionization frequency.These threshold frequencies for each ion are also taken from Verner & Yakovlev (1994, 1995).The number density of each ion n ion is taken from the calculated photoionzation balance at each radial step along the line of sight (see Section 3.2).The above sources of opacity give us the optical depth ∆τ of the material between two radial steps ∆r: ∆τ = χ ν ∆r (28) The attenuated luminosity to be used in the next radial step's calculations can be found using the optical depth ∆τ and the current SED L ν .The current radial step is denoted by the subscript n, and the next step along the radial grid is denoted by the subscript n + 1. Putting together Equations ( 27), ( 28), and the luminosity at the current L ν,n step we have The incident luminosity (described above) is only one component of the local radiation field at an arbitrary distance r from the central AGN source.Below, we describe how we model the entire radiation field as an effective flux F eff that arises from two distinct sources of radiation.
We consider two lines of sight: the primary LOS from the central source, and the secondary LOS that is a proxy for photons that undergo multiple scattering.The primary LOS is calculated according to Equation ( 29) and a standard hydrogen density model (see Section 2.2, Equation ( 7)).Along the primary LOS, the AGN's SED undergoes attenuation due to both Thomson scattering and photoelectric absorption.
In addition to the incident luminosity from the primary LOS, ionizing photons that contribute to the total flux due to scattering and reprocessing must also be accounted for in order to calculate the total effective flux F eff at the observer (Higginbottom et al. 2014).This would ideally be done with a fully self-consistent hydrodynamic and radiative transfer model; however, for simplicity we do this by implementing a luminosity "floor" at each radial step.At the first radial step (n = 0), the total flux is approximated as simply the flux due to the luminosity of the source, since at this point we assume there has been no attenuation of the SED.At subsequent radial steps, we implement an approximate luminosity floor L f loor on the attenuation of the SED.
We use the secondary LOS to calculate L f loor , which differs from the primary LOS in two ways: 1) we assume no shielding (S = 1), and 2) it is only affected by Thomson scattering.The lack of shielding is consistent with our interpretation of this path as circumnavigating the failed-wind region due to scattering and reprocessing.The assumed lack of photoelectric absorption allows us to set the attenuated secondary LOS SED as a simple "gray" fraction of the unattenuated SED.
We have assumed that the secondary LOS is a proxy for a path that involves multiple scattering events.Thus, only a fraction of the luminosity L f loor will reach the observer's location in the form of a diffuse radiation field.The probability for multiple scattering is lower than that for single scattering, but computing this ac- curately is beyond the scope of this paper (see, e.g.Higginbottom et al. 2014).Thus, we account for it using an arbitrary fraction γ, which is also randomly selected for each model from a range of reasonable values (see Table 1), such that the value of L f loor,0 is given by where L ν,0 is taken from Equation (3).
At each additional step along the radial grid, we follow a process similar to that used for the primary LOS to calculate the value of L f loor , which is given by where n e is the electron number density as calculated from the ionization balance in Section 3.2, σ T is the Thomson scattering cross-section, and ∆r is the size of the radial step.We can subsequently define the floor flux F f loor at the next radial step as Finally, we compute the total effective flux F eff that reaches the observer at radial distance r n+1 from the central AGN.The flux F LOS from the primary radial LOS (i.e., the maximally attenuated SED) is similarly calculated from Equation (29) as We include the impact of the secondary LOS only when it exceeds the most attenuated parts of the primary SED, such that the total flux at the observing point for a given frequency ν at the next radial step r n+1 is given by Figure 9 shows the evolution with radii of the effective flux F eff for an example model.For this model, the floor imposed on the flux by the secondary LOS takes effect at a radius of ∼ 10 14 cm.This newly calculated flux F eff,n+1 is then used to determine the mean intensity at each wavelength J ν for the next radial step: J ν is then used to recalculate the local gas temperature, X-ray temperature T X (Equation ( 16)), hydrogenionizing luminosity L X (Equation ( 12)), and ionization parameter ξ (Equation ( 13)) for the next radial step.We then repeat the calculations detailed in Sections 3.1 and 3.2 to find the ionization balance of the new radial step, and so on.This process is continued until the final observational radius r out has been reached.

PUTTING THE LINES IN LINE-DRIVING
Next, we describe the computation of the radiativedriving line-force multiplier M (t) for all 500 radial grid zones in each of the successful Monte Carlo model outflows.Thus, we aim to provide realistic calculations of M (t) for a comprehensively broad range of AGNrelevant physical environments.

Spectral Line Database Selection
It has been shown historically that the inclusion of as many spectral lines as possible is critical in order to achieve a an accurate characterization of the linedriving acceleration (see, e.g., Gormaz-Matamala et al. 2019).For example, when CAK performed their analysis of the radiative driving using 900 C III multiplets, they predicted mass loss rates greater than 100 times those initially found by Lucy & Solomon (1970), who used 12 doublets of C, N, Si, and S.
Since a complete tabulation of the line force requires the inclusion of as many lines as possible, the atomic line data used in this work was retrieved from a variety of sources.As in Lattimer & Cranmer (2021), we used multiple atomic databases to fill in any gaps in the available data wherever possible.Databases used include the National Institute of Standards and Technology (NIST) (Kramida et al. 2018), version 10.0 of the CHIANTI database (Dere et al. 1997;Del Zanna et al. 2021), the database of lines used by the radiative transfer code CMFGEN1 (Hillier 1990;Hillier & Miller 1998;Hillier & Lanz 2001), the Opacity Project's TOPbase (Cunto & Mendoza 1992;Cunto et al. 1993), and the XSTAR database (Mendoza et al. 2021).For the AGN central sources considered here, the more highly ionized states of each element become increasingly important.Therefore we include, where available, all ionization states for all elements up to Zn (Z = 30).The CHIANTI and CMFGEN databases additionally include both theoretical and observed spectral line data.For the sake of completeness, our line list is comprised of both theoretical (where applicable) and observed transitions for the remainder of this work.The final line list contains 5,671,602 lines, an improvement of approximately 25% over the list used in Lattimer & Cranmer (2021).An overview of the database selection process is provided in Appendix A, and specific line counts and the database used for each ion are detailed in Table A1.The compiled line list is also available for public use online2 (Lattimer 2024a).

Weighted Line Strengths
As in Lattimer & Cranmer (2021), we follow Gayley (1995) in characterizing the distribution of spectral line strengths as a set of dimensionless ratios q i .This represents the ratio of the radiative force due to a single spectral line to the radiative force due to Thomson scattering.We also calculate the dimensionless weighting factors W i , which weight the line strength (q i ) according to the illumination of the atom from the local SED.The product q i W i thus gives the full ratio of radiative acceleration due to a specific line i to the acceleration on free electrons.The dimensionless line strength parameter q i can be expressed by and the dimensionless weighting factor W i by In Equation 36, n i is the number density of ions of a given stage that are also in the lower bound level of a given line.In Equation ( 37), L ν (ν 0 ) denotes the luminosity L ν in erg s −1 Hz −1 at the specific rest-frame frequency ν 0 of the line being considered.
As in Gayley (1995), we also define the sum of the line strengths Q as For simplicity, we assume Boltzmann excitation equilibrium for atomic level populations.We then use the total number of particles in a given ionization state to express the number of particles in the lower transition level n i as where the partition function U ion (T ) is found using a version of the Cardona et al. ( 2010) method described in Section 3.1 of Lattimer & Cranmer (2021).The use of Boltzmann excitation equilibrium is not rigorously valid for an AGN environment, and in future work this assumption will be revisited in order to explore how significantly atomic level populations could be affected by non-LTE processes.
We are able to determine n ion /n e for each ion using the full ionization balance, which is calculated as described in Section 3.2.To do this, we follow the procedure outlined in Section 2.2 of Lattimer & Cranmer (2021).We can then rewrite Equation (36) in terms of known quantities: (40) Using Equations ( 37) and ( 40) respectively, we can calculate the weighting function W i and the dimensionless line strength q i for any line in our list from the calculated ionization balance and flux at any given radial point along the primary LOS.
Figure 10 shows the distribution of the weighted line strengths q i W i for an example model at various radial distances spanning from R cor to r out .For many models, such as the one shown here, the high ionization of the wind material at the inner radii close to R cor leads to a prevalence of weak lines, as even the strongest lines (i.e.those occurring at the highest values of q i W i ) are pushed to extremely low values of q i W i (∼ 10 −6 ).However, as we move to larger radii, we see that the strongest lines shift to occur at much higher strengths (q i W i ≳ 10 1 ).This is due to the decrease in the overall ionization of the wind material via the absorption of the ionizing X-rays by the material at inner radii.This suggests that the shielding by the inner, highly ionized region does indeed play an important role in allowing the outer regions of the wind to become less ionized and reach higher line strengths, thus allowing radiative pressure to become a viable driving mechanism.

The Line-Force Multiplier
We define the line acceleration as the radiative acceleration due to electron scattering multiplied by the line force multiplier M (t).Here, t is a dimensionless optical depth parameter that that does not depend on the line strength, given by t = κ e ρv th dv dr for expanding atmospheres (CAK).In this case, t is less than the electron scattering optical depth, while for static atmospheres t is equivalent to this depth.
For calculating the force multiplier, we assume that the outflow is a) supersonic, b) contains no overlapping lines, and c) that we can apply the Sobolev approximation due to large velocity gradients in the outflow.This approximation holds in rapidly accelerating outflows where |dv/dr| is much larger than v th /H, the thermal speed divided by the outflow's radial scale height (Sobolev 1957(Sobolev , 1960)).
For an expanding wind, the full form of M (t) depends on various radiative transfer effects, such as line "self-shadowing", which arises due to differences in the Doppler-shifted local reference frames (Gayley 1995).Using our compiled line list, we therefore write the full calculation of the force multiplier M (t) as where the geometric finite-disk factor η is the same ratio F given by Equation ( 21) of Gayley (1995), defined as the ratio of the true line force to that derived in the limit of purely radial photons.Because this effect can be modeled independently of the atomic physics, we use a nominal value of η = 1 for the remainder of this work.Following the standard convention, we use the proton thermal velocity v 2 th = 2k B T /m p to calculate the dimensionless optical depth τ i , which depends on the CAK t parameter: CAK proposed a power-law form of M (t), expressed as where k and α are the fit constants of the power-law.
Given the computationally intensive manner of calculating the full form of M (t) from compiled line lists, this form has since been widely used.
Although t can be calculated at any point in a radiatively-driven outflow, for this work we evaluate M (t) over a grid of 100 evenly spaced values of t ranging from t = 10 −25 to t = 10 15 .In the limit of t ≫ 1, the corresponding regions in an AGN environment are not likely to be relevant to the line-driven outflows considered here because they are either deep in the optically thick corona or in the vicinity of the ISCO.We therefore include values of t > 1 in our calculations only for the purposes of evaluating any trends in the form of the force multiplier, as t = 1 roughly corresponds to the "photosphere" of the source, and thus these high values represent a realistically unobservable region of the accretion disk.

Plasma Properties
In order to capture as complete a picture of the outflow as possible for each model, we calculate the properties of the wind at each radial step, rather than only at r out .The evolution of the ionization parameter (Equation ( 13)) as the model steps radially outward in the wind is particularly important.In general, ionization decreases with increasing radius from the central source; at small radii very close to the black hole, the wind material will be highly ionized due to the high incident luminosity at these radii.
In ballistic wind codes such as QWIND (Quera-Bofarull et al. 2023), the dense inner region close to the SMBH is often incapable of successfully launching an outflow, due to the subsequent over-ionization of the wind material.However, this resulting failed wind region is thought to play a crucial role in shielding the material at larger radii from becoming similarly overionized.In our models, we simulated this effect by introducing a shielding parameter S, as discussed in Sec-  We use our grid of t to calculate the force multiplier at every radial step along the primary LOS that exhibits a value of N H that falls within the range of reasonable values for the wide variety of observed AGN outflow features.To ensure that any edge behavior is included, we define this range as 17 < log(N H ) < 25, which is slightly wider than the range given in Laha et al. (2021).Approximately 71.8% of our modeled points fell within this region, with ∼ 28% corresponding to what we consider unobservable values.
Figure 12 shows the calculated hydrogen column density in terms of radius and ξ for all models, with the points neglected as falling outside the observable range shaded in grey.Figure 13 additionally highlights the ranges of N H and ξ for various types of AGN, taken from Laha et al. (2021), and the fraction of our modeled points that fall within those regions.The NAL, BAL, UFO, and Warm Absorber N H regions encompass approximately 4%, 27%, 23%, and 15% respectively of our modeled points, and the observable ranges of ξ account for approximately 5%, 11%, 15%, and 15% respectively.In total, ∼ 40% of our modeled points fall within the N H regions and ∼ 30% fall within ξ ranges that correspond to at least one of these four AGN types.In Figures 12-13, the regions shaded in grey, which we designate as unphysical, are the result of our random selection of initial parameters.While each initial parameter is individually selected from ranges of reasonable values, it is likely that not all model runs resulted in physically feasible combinations; some unphysical combinations of parameters will inevitably arise.This is also discussed in Section 3.1 in regards to the thermal equilibrium calculation.We therefore use the ranges of calculated column density discussed above to reduce our total set of models (556 completed runs out of 576 initiated models) to only those that fall within observed ranges of column density after all the models have been completed.

Line-Driving Forces
Accounting for the calculation of the force multiplier at each radial step, each model produces up to 500 M (t) curves.After discarding the radii that do not fall within the reasonable range of N H (78,373 points), we are left with 199,627 valid radial steps, and thus the same number of M (t) curves, with each curve having 100 optical depth points.Figure 14 shows a randomly selected subset of these curves.As in Lattimer & Cranmer (2021), we find that M (t) takes the form of a saturated powerlaw, with a generic power-law component at generally higher values of t that then flattens into a constant value as t → 0. This is described by where C, k, α, and s are fit constants that set the shape of the curve.Here, C is a constant that sets the magnitude of the flat, saturated portion as t → 0. In our M (t) curves, we easily see that C = Q.The parameter α sets the slope of the power-law region at high values of t, and the parameters k and s set the position and sharpness of the transition from the high-to low-t regions respectively.We note from Figure 14 that this turnover to a constant value of M (t) = Q occurs at varying values of t for different models and radii; this is sometimes as low as t ≈ 10 −12 .At the smallest radii, M (t) is essentially constant and the turnover to the final saturation value can occur at optical depths as high as t ≈ 10 15 .In general, the value of t that marks the transition from the power-law to constant regimes decreases as radius increases.
Likewise, the saturation value Q varies across models and radii.Figure 15 shows Q (Equation ( 38)) with respect to the ionization parameter ξ, and Figure 16 shows the evolution of Q as the temperature T varies along the radial LOS.As anticipated, the values of Q follow a general upward trend with increasing radius; i.e, Q (and thus the strength of the radiative outflow) increases as overall ionization decreases.In addition to this general trend, we find that all models exhibit a steep decrease in Q as ξ → ∞, starting at ξ ≈ 10 5 .
We also compare our values of Q to those taken from Stevens & Kallman (1990) and Dannen et al. (2019).Our models generally agree with these values within the applicable range of ξ, with the notable exception of 10 2.5 ≲ ξ ≲ 10 5 .The disagreement in this region is likely due to our inclusion of adiabatic cooling in the thermal equilibrium calculations, which was not considered in Stevens & Kallman (1990).Dannen et al. (2019) computed the cooling rate as a function of ξ using XSTAR (Mendoza et al. 2021;Dyda et al. 2017).These differences are likely responsible for the discrepancies shown in Figure 15.Our increased cooling estimates consequently lower the calculated temperatures along the radial LOS (see Figure 7), which in turn decreases the overall ionization of the material.As M (t), and thus Q, is directly dependent on the ionization of the material, this results in an increased final value of Q.
We see from Figures 14, 15, and 16 that there is great diversity in the possible M (t) curves.It is therefore useful to normalize the values of M (t) by each curve's corresponding value Q.Additionally, we designate the optical depth where the M (t) curve transitions from the saturated Q portion to the power-law portion as t * .We define this turnover point as the value of the optical depth t at which the slope of the M (t) curve α differs significantly from 0. Here, we use a conservative condition of α > 0.25, where α is given by The variation of the slope α with both the original optical depth t and the normalized optical depth t/t * is shown in Figure 17.We see here that α increases rapidly between 0 ≲ log(t/t * ) ≲ 10.Many of the models achieve a final slope of α ≈ 1 at high t/t * .The asymptotic value of α → 1 as t → ∞ is understandable from Equation (42) as the remnant behavior of the strongest lines in the limit of infinite optical depth.Note that the widely used CAK-type values of α ≈ 0.5-0.7 occur within the slow transition over many orders of magnitude of t between these two limiting behaviors.
Figure 18 shows the same subset of force multiplier curves as in Figure 14, normalized by the corresponding Q and t * values.Unsurprisingly, we find that all variation in the magnitude of M (t) at log(t/t * ) ≲ 0 is dependent solely upon Q, as all models normalize to a single curve at this point.Conversely, for log(t/t * ) ≳ 0 the normalization of M (t) and t by Q and t * respectively does not remove all variation in the power-law segment of the curves.From Section 5 we see that there are many factors that determine the final form of the force multiplier for any given outflow at a given radius.We have performed an explicit, radially dependent calculation of the ionization balance and the resulting strength of the line-driving force.While this type of explicit process may be ideal for characterizing the line-driving force, such calculations are not always computationally feasible.Therefore, we provide here a method for approximating the force multiplier using both our results and commonly used observational AGN properties.As Q sets the maximum value of the force multiplier for a given model and radius, its careful and accurate treatment is of the utmost importance.Although each initial input parameter (see Table 1) contributes to the final shape of the M (t) curve and the value of Q, we find that Q is most strongly dependent on the accretion rate ṁ and the calculated X-ray temperature T X .Recall from Figure 3 and Equation ( 16) that T X itself is strongly correlated with ṁ, as shown in Figure 6.However, T X depends not only on ṁ, but also on the degree of attenuation undergone by the SED between the coronal base at R cor and the radial point being considered (Stevens 1991).Therefore, when quantifying Q it is necessary to include the dependencies on both T X and ṁ.
For the sake of accessibility and maintaining generality, it is useful to reparameterize Q in terms of its dependence on the ionization parameter ξ and the Xray spectral index Γ X , which is more commonly used to characterize AGN spectra than X-ray temperature.This parameter describes the shape of the SED in the observed X-ray range, defined here as falling between 0.5 and 7 keV: To find the value of Γ X for each of our models, we fit a simple power-law to the SED between this range at each radial step along the primary LOS.The relationship between Γ X and T X is illustrated in Figure 19.In log-space, the relationship between Γ X and the basal value of T X at R cor is well defined by a simple cubic polynomial, given by We note that this formula differs slightly from the analytic solution that is found if the SED is assumed to follow a strict power-law across the entire X-ray frequency range, which is not the case for our attenuated SED models.This analytic solution is given by where ν 1 = 0.1 keV and ν 2 = 100 keV, as in Equation ( 16).The curve produced by this analytic solution is shown in panel (b) of Figure 19.51) and used in Algorithm 1 is shown in purple (see Section 6.2).The form of this fit is given by Equation ( 45), with fit parameters α = 0.9, s = 0.3, and k = 500.
For some models and radii, the specific combination of attenuation along the LOS produces an SED that does not follow a power-law in this range of frequency.These often result in anomalously high or low values of Γ X , so for the remainder of this work we consider only SEDs with 0.5 ≤ Γ X ≤ 3.5 (see, e.g., Pounds et al. 1990;Nandra & Pounds 1994;Page et al. 2005;Brightman et al. 2013).Approximately 94% of the radial points in our set of successful models fall within this range.Figure 20 shows the dependence of Q on ξ, T X , and Γ X .Note that the value of Q is determined more or less uniquely by its location in the two-dimensional plane defined by ξ and Γ X .However, there are combinations of parameters where there seem to exist multiple possible values of Q, such as the vertical strip in the vicinity of log ξ ≈ 2.5.Unfortunately, the fitting formulae that we provide are not able to replicate the full complexity of the models in these regions.
We now turn our attention to characterizing Q in terms of ξ and Γ X .To do this, we use a k-nearest neighbors (KNN) algorithm (Goldberger et al. 2004;van der Maaten & Hinton 2008), which is trained and tested using our calculated values of these three parameters.As mentioned in Section 3, we neglect points where ξ < 10 −6 , and the corresponding values of Q and Γ X , when training the KNN model.
To construct this model, we used the KNN regression module for python available from sci-kit learn.In this method, the target data point is predicted by a local interpolation of the nearest points in the training set (Pedregosa et al. 2011).For constructing this model, each point in the local parameter-space made up by Q, ξ, and Γ X contributed uniformly to the determi-Figure 19.The evolution of the X-ray spectral index Γx with X-ray temperature TX and radius r.Panel (a) shows ΓX for the entire set of modeled radial points, with the points corresponding to Rcor marked in black.Panel (b) shows ΓX at Rcor.The analytic relationship corresponding to a pure power-law X-ray SED is marked with the grey dashed line, and the fit derived from our modeled points (Equation ( 48)) is denoted in red.
nation of a given query point.The initial training set was created by randomly splitting our model calculated values of Q and the corresponding Γ X and ξ points.Approximately 80% of the modeled points were allocated to the training set, with the remaining 20% reserved to judge the efficacy of the model.After initial fitting and validation using the split sets of points, we reallocated the test points to the training set in order to provide as refined a fit as possible.
To validate the KNN model after training, we compare the KNN predictions Q KN N with our calculated values Q. Figure 20 illustrates the dependence of Q and Q KN N on ξ and Γ X .The final panel (c) of Figure 20 gives Q KN N as calculated from generic grids of ξ and Γ X , where −6 < log(ξ) < 11 and 0.5 < Γ X < 3.5, with 250 points in each grid.Figure 21 gives Q calculated from the full LOS ionization balance model, overlaid by the KNN-predicted values.From Figure 21, we see the accuracy with which the KNN model is able to reproduce the values Q of the testing set when given the same input parameters.To further quantify the validity of these results, we also plot Q against Q KNN (Figure 22).The KNN model results in an r-value of r = 0.95; thus, we conclude that it provides an acceptable prediction of the value of Q for a given input ξ and Γ X .

Normalized Optical Depth t *
It is also useful to calculate the normalized optical depth t * for a given Q.We do this by fitting a second order polynomial to t * as a function of Q, as shown in Figure 23.While there is some variation in Q for any given value of t * , we find that the relationship between the two is well described by log(t * ) = −0.01log 2 (Q) − 1.18 log(Q) − 1.97.( 50) By estimating the saturation value Q of M (t) at low-t, and thus the turnover point t * , we are able to deduce the normalized force multiplier curve.By fitting Equation (45) to the normalized force multiplier curves of our modeled outflows, we determine the constants α, k and s.We find that the power-law slope of the high-t region for our models is well described by α = 0.9, and that the curves are overall reasonably characterized by setting the remaining parameters to k = 500 and s = 0.3 (shown in Figure 18).The general shape of a given normalized M (t) curve is therefore given by M norm = 500 (500 0.3 + (t/t * ) 0.27 ) 3.33 . (51) The normalized M (t) curve in turn allows the actual force multiplier curve to be determined for a given input of ξ and Γ X .A self-contained function for this calculation, as outlined above and summarized in Algorithm 1, is available online for public use 3 (Lattimer 2024b).However, we stress that this is an estimation, and where possible the full radiative transfer effects should be considered.

Conclusions & Future Work
Algorithm 1: Process for Estimating M (t) Input: ξ, ΓX , t Compute Q(ξ, ΓX ) from KNN model, as in Section 6.1 Compute t * (Q) using Equation (50) from Q(ξ, ΓX ) Compute t/t * Compute M (t)/Q(ξ, ΓX ) from Equation ( 51) → Compute M (t) using Q, t * , and Equation ( 51) In this work we have provided new calculations of the radiative line-driving force in physical environments rel- evant to AGN outflows.We carried out a comprehensive and self-consistent treatment of the ionization balance in the outflow material and calculated the resulting force multiplier M (t) for a variety of AGN and wind parameters.Because these parameters varied over a large multi-dimensional space of reasonable values, we found a much broader variety in the magnitudes of the linedriving force multiplier than have previously been seen in the results of other work.We found that the ionization parameter ξ, while often used to describe the ionization state of the material, is insufficient to completely describe the radiation force.Therefore, we provided an updated method of estimating the force multiplier us-   ing the results of the models that are presented in this paper.
This result has important implications for calculations of outflow properties that currently depend on ξ, and may shed light on the nature of BAL outflows in general.It has been suggested by several BAL outflow studies that these outflows may exist at much larger radii than initially thought (see, e.g. de Kool et al. 2001;Moe et al. 2009;Dunn et al. 2010;Borguet et al. 2013).For example, based on a sample of 42 outflows, Arav et al. (2018) found that up to 50% of quasar outflows are located at distances greater than 100 pc.Given these results, it may naturally appear difficult to reconcile these kiloparsec scales with the much smaller launch radii (0.01-0.1 pc) typically exhibited by simulations of line-driven disk winds (Matthews et al. 2023;Borguet et al. 2013).This would pose questions such as whether BAL outflows should be considered "disk winds" in the traditional sense at all, or if these outflows are due to some combination of driving mechanisms (Choi et al. 2022;Matthews et al. 2023).However, the models presented here are applicable out to R out = 3000 pc, and our results are consistent with the large spatial scales found by these previous studies.Our models suggest that line-driving remains a viable driving mechanism out to these large radii, with the inner shielding region of over-ionized material allowing these outer regions to reach high line strength and force multiplier values.
By estimating the saturation value Q of M (t) at low-t, and thus the turnover point t * , in terms of commonly used AGN parameters such as the ionization parameter ξ and the X-ray spectral index Γ X , we provide a new method of estimating the force multiplier M (t) at any given point in the wind.Γ X is dependent on the initial conditions that at the base of the wind that set the initial SED.For example, the accretion rate ṁ plays a prominent role in determining the initial SED, also sets Γ X at each radial point.The attenuation of the SED by the intervening gas also strongly affects the value of Γ X .Because of these implicit dependencies, we choose to parameterize Q in terms of ξ and Γ X , rather than column density N H as suggested in Stevens (1991).
For future work, we plan to incorporate the force multiplier results and the model framework presented here into QWIND, a publicly available ballistic AGN-outflow code (Risaliti & Elvis 2010;Quera-Bofarull et al. 2023).While these types of models neglect forces due to gas pressure gradients and MHD effects, they provide an efficient way to estimate the speeds, densities, and geometric streamline trajectories in the limiting case of supersonic radiation-driven flows.Currently, QWIND calculates the force multiplier as a function of both the optical depth parameter t and the ionization parameter ξ.By replacing the current force multiplier with the M (t) values calculated using the explicit ionization balance, we plan to explore to what extent the successful launching and subsequent strength of the outflow depends on the individual components of ξ that are explored in this work.In addition, a full dynamical model would allow us to replace our parameterized expressions for adiabatic cooling, which may have overestimated the magnitude of this effect, with a more self-consistent treatment of the cooling processes.
It will also be useful to refine the assumed atomic level excitation balance used for calculating the ionization balance of the outflow.The ionization and recombination rates described in Section 3.2 should be updated with the most recent results from laboratory and theoretical studies (e.g., Bryans et al. 2006;Badnell 2007;Dere et al. 2009;Kallman et al. 2021).Also, computing non-LTE level populations presents a consistent challenge in the calculation of the force multiplier, as doing so requires complete atomic data that is often not available for the requisite number of spectral lines.As discussed, the models presented here likewise assume a Boltzmann excitation equilibrium for atomic level populations.However, some work suggests that this may result in the erroneous inclusion of many spectral lines, which in turn could increase the calculated values of M (t) (see, e.g.Dannen et al. 2019).Therefore, replacing the current level populations with a non-LTE excitation equilibrium is an important next step that will shed light on the magnitude of this effect and consequently any influence it may have on the calculated properties of the modeled outflows.
Additionally, we note that although we implement an imposed floor on the SED, as discussed in Section 3.3, our absorption model is most accurately representative of a monolithic slab of wind material.In future work, it will be important to explicitly consider the stochastic or fractal structures of wind, which would attenuate the source SED differently to how we have calculated it here (Jin et al. 2022;Locatelli et al. 2022).
Lastly, in some cases thermal (Proga & Kallman 2002) or magnetic (Miller et al. 2006) driving appears to be a more promising wind-driving scenario for AGN disk winds than radiative driving.For example, it has been suggested that the outflow features observed from the low-mass X-ray binary GRO J1655-40, composed of a stellar mass black hole (M BH = 5.4 ± 0.3M ⊙ ) and a companion (M = 1.45 ± 0.35M ⊙ ), could be explained by magnetic driving (Miller et al. 2006(Miller et al. , 2008)).However, recent results from radiation hydrodynamic simulations have suggested that a thermal-radiative model may better match the observations (Tomaru et al. 2023).Furthermore, magnetic driving may not accurately reproduce the overall ion populations for the outflow material around GRO J1655-40 without including additional free parameters, implying that MHD driving is a less efficient method for launching and maintaining this outflow (Tomaru et al. 2023).Possible interactions between magnetic and radiative driving will be an interesting avenue to explore in future work.
Our models here provide a careful treatment of the radiative transfer processes, however, they do not currently incorporate the effects of MHD or thermal driving.In order to develop more comprehensive simulations of disk outflows, including determining which driving mechanisms are dominant, future work should include a thorough treatment of the MHD and thermal effects, as well as full time-dependence of the photoionization (Proga & Kallman 2004;Davis & Tchekhovskoy 2020;Nomura et al. 2021;Waters et al. 2021;Rogantini et al. 2022;Higginbottom et al. 2024).As these complex interactions are often simplified in modeling AGN disk winds, the development of a model that treats these self-consistently will be of utmost importance in determining how these effects contribute to the strength and sustainability of these outflows.
This work was supported by the Future Investigators in NASA Earth and Space Science and Technology (FINESST) program, through NASA grant 80NSSC22K1748.It was also supported by start-up funds from the Department of Astrophysical and Planetary Sciences at the University of Colorado-Boulder.CHIANTI is a collaborative project involving George Mason University, the University of Michigan (USA), University of Cambridge (UK) and NASA Goddard Space Flight Center (USA).This research made use of NASA's Astrophysics Data System (ADS).The authors would additionally like to thank the anonymous reviewer for their helpful insights.To compile the line list used in this work, we follow the general process outlined in Lattimer & Cranmer (2021).In cases where more than one database listed atomic data for a given ion, the database with the largest number of available transitions was used for each ionization state of each element.For several ions, the atomic data necessary to calculate the line strengths and weights were not available from the databases used.For example, the available sources do not provide data for a majority of the ions of both Cu and Zn.After compiling the line list from the total available data, any duplicate transitions were discarded, such that each line is only represented once in the final line list.Additionally, we exclude any listed lines produced by transitions where i > j and those with anomalously large gf values, which we define as gf > 10 2 .The CHIANTI and CMFGEN databases used here both include theoretical transitions which have not been observed experimentally.For the sake of completeness, in this work we do not differentiate between these two subsets of lines and include both observed and theoretical transitions in our final line list.Table A1 shows a breakdown of the final line counts n used by ion, with the source database for each also listed.For the initial source SED, we calculated T along the primary LOS.For this calculation, the temperature depends only on the ionization parameter and the shape of the input SED.Thus, we represent the LOS using a grid of 500 values of ξ, evenly distributed in logspace.We set the limits of ξ, corresponding to the first and last radial zones, as log(ξ in ) = 6 and log(ξ out ) = −2, in accordance with the ionization parameter range used in Mehdipour et al. (2016).Using this even spacing in ξ, we then computed the temperature using the method described in Section 3.1.We performed this calculation using both AGN1 and AGN2 as the initial input SED. Figure B2 shows the temperature curves corresponding to those produced by CLOUDY, SPEX, and XSTAR for AGN1 and AGN2.Comparative curves produced by our code are also plotted for several values of the X-ray temperature T X .The various energy limits used in determining T X are shown in Table B1.
We see from Figure B2 that our curves produce a closer fit to the CLOUDY, SPEX, and XSTAR curves at higher values of T X .When neglecting adiabatic cooling, our thermal equilibrium method depends solely on T X , which in turn depends only on the chosen X-ray frequency limits and the shape of the initial SED.When using our default frequency limits ν min and ν max , corresponding to energies E min = 100 eV and E max = 100 keV, we find log(T X ) ≈ 8.38 for AGN1 and log(T X ) ≈ 8.49 for AGN2.These fall near the middle of the curves shown in Figure B2.As the value of T X increases, the agreement between the calculated temperatures and those found by the three comparison codes improves.Generally, these higher values of T X correspond to the use of wider frequency ranges in calculating the X-ray temperature when using Equation ( 16).
The discrepancies between our calculated temperature values and those computed by CLOUDY, SPEX, and XSTAR are likely due to differences in the treatment of the considered heating and cooling processes, as we see similar differences to the values calculated by Dannen et al. (2019) (see Figure 7), who used XSTAR to calculate their temperature values.Additionally, we note that each heating and cooling process incorporated into Equation ( 17) has its own unique dependence on the shape of the SED.Although Equation ( 17) generalizes these dependencies to a single value of T X for the sake of simplicity, we consider any inaccuracies that may arise as a result to be negligible once the effects of adiabatic cooling are included (as in Figure 7(a)), since this has a much greater effect on the thermal equilibrium.

B.2. Ionic Fractions
We also present the ionic fractions, defined as n ion /n el , for select important elements.We note that this is a related but separate quantity from the mean ionization number I, shown in Figure 8.The ionic fractions detailed here represent the fraction of ions in a given ionization state when compared to the total number density of the given element.We compare these to the values taken from Mehdipour et al. (2016).When calculating the ionization balance from the unobscured AGN1 SED for comparison to the results from CLOUDY, SPEX, and XSTAR we neglect the attenuation of the SED due to the intervening medium.While this attenuation is included in our main set of models, we make this modification in order to maintain similar conditions within the wind material as those in Mehdipour et al. (2016).
The inner and outer edges of the radial grid along the line of sight are given by respectively, where log(ξ in ) = 6 and log(ξ out ) = −2, as defined in Section B.1 from Mehdipour et al. (2016).We also adopt their convention of using a constant value of the hydrogen density, with n H = 10 8 cm −3 , and set the ionizing luminosity L X = 1.17 × 10 44 erg s −1 for the unobscured AGN1 (Mehdipour et al. 2015(Mehdipour et al. , 2016)).Figures B3-B9 show the calculated ion fractions, alongside those taken from Mehdipour et al. (2016) for elements O, Si, and Fe.These elements were chosen to be consistent with Figure 8 wherever possible.
We see from Figures B3-B9 that our ionization balance calculations tend to produce results in good agreement with the three comparison codes for lower-Z elements.This trend is also shown for low ionization states.While the agreement between the four plotted curves begins to decrease for high-Z and high ionization states, such as high Fe ionization states, this is not unexpected given the deviations found between the different thermal equilibrium curves in Section B.1.We note also that the three comparison codes tend to show higher degrees of mutual disagreement for these same ions (see, for example, Fe VII).
The discrepancies between our ionic fractions and those calculated by CLOUDY, SPEX, and XSTAR are likely due to a combination of propagating effects from the thermal equilibrium calculation as well as differences in the treatment of the cloud of wind material.For example, Mehdipour et al. (2016) assumes a geometry of a central AGN with a single SED (AGN1) obscured by a surrounding cloud, which in turn produces a separate SED (AGN2).However, our models treat these two components together, making a direct comparison between the two difficult.Additionally, Mehdipour et al. (2016) used different set of elemental abundances (those from Lodders et al. (2009)) in their calculation of the photoionization equilibria resulting from the AGN1 and AGN2 SEDs.It is worth noting that these abundances, given in Table 1 of Mehdipour et al. (2015), do not include a number of elements that are included in Asplund et al. (2021), which we use in our calculations.This may also result in inconsistencies in the resulting ionization balances.Given these factors, we consider the agreement between the calculated ionization balance and those found by the three comparison codes to be within acceptable limits for the purposes of this work.

Figure 1 .
Figure 1.Idealized example illustration of AGN disk geometry showing the variety of possible observed outflow features and regions.See text for acronym definitions.

Figure 2 .
Figure 2. Modeled radial line-of-sight geometry at an arbitrary angle θ above the disk.Black dashed line shows the modeled line of sight and gray box represents the gas parcel.The radius of the hot corona and the edge of the disk are denoted by Rcor and R self respectively.A representative version of the flux tube, shown in green, is included for completeness.Distances and sizes are not to scale.

Figure 3 .
Figure 3. Luminosity Lν as calculated by QSOSED at the first radial step Rcor for (a) varying black hole mass MBH and (b) varying accretion rate in Eddington units ṁ.The curves corresponding to the lowest value of MBH and ṁ respectively are marked with a dotted black line in each panel.The portion of the X-ray region (0.5-7 keV) used to define X-ray spectral index ΓX is shaded in gray.

Figure 4 .
Figure 4. Dependence of the self-gravity radius R self and the corona radius Rcor on the mass accretion rate onto the black hole in Eddington units, ṁ.Rcor shows little dependence upon the mass of the central black hole (MBH).

Figure 5 .
Figure 5. Representative versions of the hydrogen number density nH (Equation (7)) and radial flow speed ur (Equation (8)).The radius of the hot corona and launch radius of the wind are denoted by Rcor and rL respectively.The number density at the launch radius, nL, is also marked.

Figure 6 .
Figure 6.Calculated X-ray temperature TX , plotted against the accretion rate ṁ as it varies along the primary LOS.Black crosses indicate the basal value of TX at the first radial step Rcor.

Figure 7 .
Figure 7.The calculated temperature T , varying with ionization parameter ξ and compared to Dannen et al. (2019).Panel (a) shows T as calculated from the full form of Equation (14).Panel (b) gives T as calculated neglecting the adiabatic cooling term in Equation (14).

Figure 8 .
Figure 8. Mean ionization number I, given by Equation (26) for (a) Z = 1 (b) Z = 8 (c) Z = 14 and (d) Z = 26 with respect to ionization parameter ξ.We see that high ionization parameter generally corresponds to radial points closer to the central source, i.e. high ξ → low r.

Figure 11 .
Figure 11.The hydrogen number density nH varying with radius for a representative subset of models, colored according to the input accretion rate ṁ.
tion 2.2, to vary the number density of this inner region prior to the launch radius.Figure 11 shows how the hydrogen number density n H varies along the radial line of sight for each model, according to the input parameters S and β.The step-function drops in n H that are seen here clearly indicate the point at which the shielded failed-wind region of each model transitions to the corresponding outflowing region.

Figure 12 .
Figure 12.Values of the hydrogen column density NH along the primary LOS for each model run as calculated from Equation (11).Dashed lines show the limits of what we consider a reasonable values of NH .Colored points indicate radii for which the force multiplier M (t) was calculated, i.e. those that we consider to fall within the observable range.

Figure 13 .
Figure 13.Colored points indicate the model-calculated values of NH that fall within ranges of NH and ξ (shaded) corresponding to (a) NAL, (b) BAL, (c) UFO, and (d) Warm Absorber AGN observations.The complete set of points of NH as calculated from Equation (11) are plotted in grey for each panel.

Figure 14 .
Figure 14.Calculated values of M (t) for a representative subset of models and radii, colored according to radius.

Figure 15 .
Figure 15.The evolution of Q with radius, plotted against the ionization parameter ξ for each model.Solid black line indicates values taken from Stevens & Kallman (1990), and solid green lines correspond to the values found in Dannen et al. (2019) for the two AGN considered .

Figure 16 .
Figure 16.The evolution of Q with radius, plotted against the calculated temperature T for each model. .6.1.Determining Q

Figure 17 .
Figure 17.M (t) slope α as defined in Equation (46), varying with (a) log(t) and (b) log(t/t * ).Solid black line in panel (b) corresponds to the calculated slope α for the best fit curve shown in Figure 18.

Figure 18 .
Figure18.Normalized force multiplier M (t/t * )/Q for the set of successful models, shown here in black.The best fit curve given by Equation (51) and used in Algorithm 1 is shown in purple (see Section 6.2).The form of this fit is given by Equation (45), with fit parameters α = 0.9, s = 0.3, and k = 500.

Figure 20 .
Figure 20.Dependence of Q on (a) ionization parameter ξ and X-ray temperature TX , and (b) ξ and X-ray spectral index ΓX .Panel (c) shows the dependence of Q KN N , as calculated via the KNN model described in Section 6.1, for general grids of ξ and ΓX .

Figure 21 .
Figure 21.Comparison of the calculated and and KNNpredicted values of Q for the testing subset of points (∼ 20% of the total points).Q calc points are plotted in black, with Q KNN overlaid in purple.

Figure 22 .
Figure 22.Correlation of Q calc and Q KNN , plotted in black.Points are plotted with 50% opacity, such that areas with high concentrations of points appear darker, and those with fewer points appear lighter.The line corresponding to Q calc = Q KNN is marked in purple.

Figure 23 .
Figure 23.The normalized optical depth t * , plotted against the saturation value Q for each model.The fit described by Equation (50) is shown in purple.

Figure B1 .
Figure B1.The source SEDs of NGC 5548, adapted from Mehdipour et al. (2015, 2016).The solid line (AGN1) depicts the unobscured AGN, and the dashed line (AGN2) represents the obscured SED of the same source.The range of frequency used in Section 3.1 to calculate TX , corresponding to 0.1 keV-100 keV, is shaded in gray.

Figure B2 .
Figure B2.Temperature as a function of ionization parameter ξ, calculated for the source SEDs given in Figure B1.Panel (a) shows the values calculated when using the AGN1 initial SED, and panel (b) shows those found using AGN2.Dashed, dotdashed, and dotted lines indicate the temperature equilibrium curves produced by CLOUDY, SPEX, and XSTAR respectively (Mehdipour et al. 2016).Solid colored curves show the thermal equilibrium temperatures calculated via the method described above and in Section 3.1.The blue shaded region indicates temperature values for the full set of models discussed in the body of this work.

Figure B3 .
Figure B3.Ionic fractions nion/nO for oxygen (Z = 8) ions.Colored dashed, dot-dashed, and dotted lines represent values calculated by CLOUDY, SPEX, and XSTAR respectively.Solid black lines indicated values found from the ionization balance described in Section 3.2.

Figure B4 .
Figure B4.Ionic fractions nion/nSi for silicon (Z = 14) ions.Colored dashed, dot-dashed, and dotted lines represent values calculated by CLOUDY, SPEX, and XSTAR respectively.Solid black lines indicated values found from the ionization balance described in Section 3.2.

Figure B6 .
Figure B6.Ionic fractions nion/nFe for iron (Z = 26) ions.Colored dashed, dot-dashed, and dotted lines represent values calculated by CLOUDY, SPEX, and XSTAR respectively.Solid black lines indicated values found from the ionization balance described in Section 3.2.

Table 1 .
The ranges of values and the associated units used to select the initial input parameters for each model run.Since these ranges by necessity cover many orders of magnitude for the majority of these parameters, the random sampling of each is done in logspace, with the exception of the density exponent β and the launch radius, rL.

Table A1 .
Database and number of lines (n) for each ion.A dash (-) indicates that no data were available.CMFGEN, NIST, CHIANTI, and TOPbase are abbreviated as CM, N, CH, and T respectively.

Table B1 .
Energy limits and the resulting values of TX used in calculating the temperature curves shown in FigureB2.