Testing the disk-corona interplay in radiatively-efficient broad-line AGN

The correlation observed between monochromatic X-ray and UV luminosities in radiatively-efficient active galactic nuclei (AGN) lacks a clear theoretical explanation despite being used for many applications. Such a correlation, with its small intrinsic scatter and its slope that is smaller than unity in log space, represents the compelling evidence that a mechanism regulating the energetic interaction between the accretion disk and the X-ray corona must be in place. This ensures that going from fainter to brighter sources the coronal emission increases less than the disk emission. We discuss here a self-consistently coupled disk-corona model that can identify this regulating mechanism in terms of modified viscosity prescriptions in the accretion disk. The model predicts a lower fraction of accretion power dissipated in the corona for higher accretion states. We then present a quantitative observational test of the model using a reference sample of broad-line AGN and modeling the disk-corona emission for each source in the $L_X-L_{UV}$ plane. We used the slope, normalization, and scatter of the observed relation to constrain the parameters of the theoretical model. For non-spinning black holes and static coronae, we find that the accretion prescriptions that match the observed slope of the $L_X-L_{UV}$ relation produce X-rays that are too weak with respect to the normalization of the observed relation. Instead, considering moderately-outflowing Comptonizing coronae and/or a more realistic high-spinning black hole population significantly relax the tension between the strength of the observed and modeled X-ray emission, while also predicting very low intrinsic scatter in the $L_X-L_{UV}$ relation. In particular, this latter scenario traces a known selection effect of flux-limited samples that preferentially select high-spinning, hence brighter, sources.


Introduction
The development of an in-depth understanding of accretion physics in active galactic nuclei (AGN) has tended to lag behind in comparison to other accreting objects (e.g., X-ray binaries, cataclysmic variables, and protoplanetary disks), for which many more observational constraints are available. While it seems that the standard thin-disk model (Shakura & Sunyaev 1973, hereafter SS73) is not able to fully explain the plethora of accreting sources that we observe (e.g., Koratkar & Blaes 1999;Blaes 2007;Antonucci 2015), it is still unclear to what extent this simple but effective prescription has to be improved (Kishimoto et al. 2008;Capellupo et al. 2015Capellupo et al. , 2016. Since the first AGN X-ray spectral surveys were performed (e.g. Elvis et al. 1978;Turner & Pounds 1989), the need for an additional spectral component to extend the cold-disk's keV temperatures was evident. This so-called X-ray "corona" (e.g. Liang & Price 1977;Galeev et al. 1979) is now almost universally considered as a hot (∼10 9 K), optically thin (τ 1) plasma up-scattering the disk photons via thermal Comptonization (Haardt & Maraschi 1991, 1993Haardt et al. 1994;Stern et al. 1995), although an additional warm component is sometimes needed to fit the softest X-rays (Petrucci et al. 2018;Kubota & Done 2018, and references therein). The proximity of the corona to the central black hole was immediately suggested by its strong and fast variability (e.g. McHardy 1989) and by the reflection signatures (Lightman & White 1988;Pounds et al. 1990;Nandra et al. 1991;Williams et al. 1992;Tanaka et al. 1995), but in-depth information regarding its geometry and formation mechanism is still lacking.
The geometry of the corona can be constrained via the observation of X-ray reverberation lags (Fabian et al. , 2017De Marco et al. 2013;Uttley et al. 2014), that seem to show a, possibly non-static, corona extending vertically and radially over the underlying disk for a few and a few tens of gravitational radii, respectively (Wilkins et al. 2016). The compactness of the corona and the origin of the X-rays close to the black hole also appear to be confirmed by micro-lensing results (e.g. Mosquera et al. 2013;Reis & Miller 2013).
As far as the formation of the corona is concerned from the theoretical point of view, the most likely explanation for it is that it is magnetically-dominated with an efficient saturation of the magnetic field that is amplified via the magneto-rotational instability (MRI, Chandrasekhar 1960; and extending buoyantly upward (and downward) from the denser parts of the disk (Galeev et al. 1979;Stella & Rosner 1984;Di Matteo 1998;Merloni & Fabian 2002;Blackman & Pessah 2009). Magnetic reconnection can then keep the corona hot (e.g., Liu et al. 2002;Uzdensky & Goodman 2008;Uzdensky 2016;Beloborodov 2017;Werner et al. 2019;Ripperda et al. 2019). This scenario seems to be supported by magneto-hydrodynamic (MHD) simulations (Miller & Stone 2000;Uzdensky 2013;Bai & Stone 2013;Jiang et al. 2014;Salvesen et al. 2016;Kadowaki et al. 2018), although only qualitative comparisons with observations have been made so far (however, see Schnittman et al. 2013). Much effort has, nonetheless, been put into trying to shed light on the physics of the disk-corona system (see Blaes 2014) and this will continue with global 3D radiation-MHD simulations (e.g. Jiang et al. 2017a), that are now approaching sub-Eddington flows as well (Jiang et al. 2019).
Observationally, the increase in quality and quantity of available AGN X-ray-to-UV data from large samples can provide insightful, and more easily approachable, diagnostics. The smoking gun of the disk-corona interplay in radiatively efficient AGN is given by the non linear correlation observed between the 2 keV and 2500 Å monochromatic luminosities (e.g., Vignali et al. 2003;Strateva et al. 2005;Steffen et al. 2006;Young et al. 2009;Lusso et al. 2010;Lusso & Risaliti 2016, and references therein), that persists throughout the common observed X-ray and optical-UV bands . Despite the possible differences arising from different sample selections and regression techniques, most observations point towards a log L X − log L UV correlation with a slope ≈0.6, a dispersion that can be as small as σ ≈ 0.2 dex (Lusso & Risaliti 2016;Chiaraluce et al. 2018), and no apparent redshift dependency. Such a tight correlation paved the way for quasars to provide an alternative standard candle for cosmographic studies (Risaliti & Lusso 2015. The slope, being smaller than unity, indicates that from lowly to highly accreting AGN, the disk emission increases more than the corona emission (e.g., Kelly et al. 2008) with crucial implications for the physics governing the coupled disk-corona system. However, a solid and conclusive theoretical explanation, for what is one of the most studied multi-wavelength observables in AGN, is still lacking.
The goal of this paper is indeed to test a self-consistently coupled disk-corona analytic model against the observed L X −L UV . Given the existing gap between simulations and observations, we argue that the use of simplified (but motivated) prescriptions still represents a powerful tool to explain observed disk-corona scaling relations, as it was done with the X-ray photon index (or the X-ray bolometric correction) correlation with the Eddington ratio (Wang et al. 2004(Wang et al. , 2019Cao 2009;Liu & Liu 2009;You et al. 2012;Liu et al. 2012Liu et al. , 2016a, or with the log L X − log L UV itself (Lusso & Risaliti 2017, hereafter LR17;Kubota & Done 2018). We here rely uniquely on the log L X − log L UV relation, since monochromatic L X and L UV values can be directly obtained from spectral fits. Forward modeling monochromatic luminosities circumvents difficulties and issues typical of model comparisons with accretion rate, Eddington ratio or bolometric luminosity estimates (e.g., Richards et al. 2006;Davis & Laor 2011;Slone & Netzer 2012;Krawczyk et al. 2013;Capellupo et al. 2015Capellupo et al. , 2016Kilerci Eser & Vestergaard 2018).
We describe our disk-corona model in Sect. 2 (and Appendix A) and we briefly show its qualitative predictions in Sect. 3. Then, we outline the observational test that we put forward to thoroughly understand the disk-corona interplay in Sect. 4 and we show the results in Sect. 5. Throughout this work, we quote median values with 16th and 84th percentiles unless otherwise stated.

The disk-corona model
The disk-corona model adopted in this work is largely based on the prescriptions put forward by Merloni (2003, hereafter M03; see also Merloni & Fabian 2002, in which the standard conservation equations of a geometrically-thin and opticallythick accretion disk (SS73; Pringle 1981) are self-consistently coupled with the X-ray corona, indicated as the fraction f (e.g., Haardt & Maraschi 1991;Svensson & Zdziarski 1994) of accretion power (per unit area, Q + ) that is dissipated away from the cold disk (e.g., Stella & Rosner 1984;Di Matteo 1998): with Q cor = v D P mag and Q + = 3 2 c s τ rφ , where v D is the vertical drift velocity (taken proportional to the Alfvén speed via an order-unity constant b), P mag = B 2 /8π is the magnetic pressure, c s is the sound speed and τ rφ is the vertically-averaged stress tensor.
The stress tensor can be assumed to be dominated by Maxwell stresses (e.g., Hawley et al. 1995;Sano et al. 2004;Minoshima et al. 2015), from which we can write τ rφ = k 0 P mag , with k 0 being a constant of order unity (Hawley et al. 1995).
To build a self-consistent solution to the accretion problem, we need to relate the stress tensor (via the magnetic pressure) with local quantities that standard analytic models are familiar with. As a matter of fact, the α−prescription is not the only educated guess that is adopted to dodge our ignorance of the physical mechanism producing the disk viscosity. Within the same theoretical framework, fundamental modifications to the viscosity law can be introduced depending on whether the viscous stress is assumed to scale proportionally with the total (P tot , gas plus radiation) pressure (SS73), with the gas pressure alone (Lightman & Eardley 1974;Sakimoto & Coroniti 1981;Meyer & Meyer-Hofmeister 1982;Stella & Rosner 1984) or with the geometric mean of the two (M03; Ichimaru 1977;Taam & Lin 1984;Burm 1985). It was soon discovered that the first prescription leads to thermally and viscously unstable disks in the radiation-pressure dominated regions, with the first instability acting on shorter timescales (Lightman & Eardley 1974;Shakura & Sunyaev 1976;Pringle 1976). This encouraged many authors (Hoshi 1985;Szuszkiewicz 1990;Merloni & Nayakshin 2006;Grzȩdzielski et al. 2017a) to generalize the viscosity law. Recent simulations (albeit of gas-pressure dominated disks only) indeed seem to show a power-law stress-pressure relation (Sano et al. 2004;Minoshima et al. 2015;Ross et al. 2016;Shadmehri et al. 2018), with an index varying from zero to one according to the different assumptions.
Here, we address this issue generalising the model reported in M03 with: where α 0 is a constant, generally not equal to α SS73 = P mag /P tot . This behavior is physically motivated by the MRI prescriptions, as its growth rate was shown to depend on the P rad -to-P gas ratio (Blaes & Socrates 2001;Turner et al. 2002) influencing the level of the magnetic field saturation. Equations (1) and (2) provide the closure equation of the disk-corona system: where k 1 = 3k 0 /2b gathers the model's uncertainties in an order unity factor (M03). Its exact value only affects f at its maximum f max = 2α 0 /k 2 1 and not the nature of what is described throughout this paper. A135, page 2 of 16 The model is then completed with the equation of state: and with a density-and temperature-dependent opacity κ = κ (ρ, T ). We compute the opacity value self-consistently with the density and temperature at each radius with an iterative process, using as reference stellar opacity tables (at solar metallicity) from the Opacity Project (Seaton et al. 1994;Seaton 1995). This is important since the density and temperature regimes relevant for AGN disks imply opacities that can be significantly different from the electron scattering value (e.g., see Jiang et al. 2016;Czerny et al. 2016;Grzȩdzielski et al. 2017b). Further, we assume a downward component of the X-ray emission (η) and a disk albedo (a disk ), which modify the disk equations from the usual (1 − f ) factor (M03; Svensson & Zdziarski 1994) to: We here adopt η = 0.55 and a disk = 0.1, respectively (e.g., Haardt & Maraschi 1993). These are typical values for anisotropic Comptonization in plane-parallel geometry, although more generally the product η (1 − a disk ) can be a function of the photon index Γ (Beloborodov 1999;Malzac et al. 2001) and of the disk's vertical structure. For simplicity, we adopt dimensionless units for the black hole mass, the accretion rate, the radial distance and the vertical scale-height: The equations for h, mid-plane ρ (g cm −3 ), P (dyn cm −2 ) and T (K), with the closure equation for f , are reported in Appendix A in Newtonian approximation (however, see Merloni & Fabian 2003, for a relativistic derivation of the µ = 0.5 case), along with the related radial profiles (Fig. A.1). We note that a constant efficiency of 0 = 0.057, typical of nonrotating black holes, and a no-torque inner boundary condition (J(r) = 1 − √ r 0 /r, with r 0 = 3 and r out = 2000) are initially adopted.
Once m,ṁ, α 0 , µ and f max are fixed, one can numerically solve the closure equations for f at each radius (see the last rows of Eqs. (A.1) and (A.3), respectively). The left-hand side is equal to P rad /P gas and we can infer the correct regime and compute the main physical quantities at the mid-plane (ρ, P, T , κ). Then, the effective temperature at the surface is computed: where we take τ(r) = h(r) ρ(r) κ(r). Monochromatic optical-UV luminosities can be then easily computed in the multi-color blackbody approximation: where π B ν (T eff ) is the black-body flux at the frequency ν and temperature T eff (r). In this framework, the energy per unit area dissipated in the corona at each radius is Q cor (r) = f (r)Q + (r), although only a fraction (1−η) will contribute to what is observed as X-ray emission: Here, we did not include the component reflected by the disk (given by the fraction f η a disk ), so that we could easily extrapolate at each radius a monochromatic value, for instance L 2 keV , assuming a simple power-law spectrum within ν i = 0.1 and ν f = 100 keV: The model relies on the assumption that a plane-parallel geometry holds for bright radiatively-efficient sources, lying in a sweet spot of accretion rate (ṁ approximately from a few percent to Eddington). Hence, the accretion disk extends down to the innermost stable circular orbit (ISCO) and no advection is included. A scripted version of the model outlined in this section will be made publicly available online 1 .

Radial profiles for f and L 2 keV
In Fig. 1 we show as an example radial profiles of f and L 2 keV , obtained by solving Eqs. (A.1), (A.3), (9) and (10). For simplicity, we fixed α 0 = 0.02 and f max = 0.5 and used the three values of µ corresponding to the most-used viscosity laws, namely µ = 0, 0.5 and 1, for P mag proportional to P tot , P gas P tot and P gas , respectively. Other values of µ would support the same picture with analogous intermediate profiles. A range of typical m, m and X-ray spectral slopes was chosen, following the distribution of objects observed in the survey field adopted for the observational test (see Sect. 4 below), namely with median values (and related 16th and 84th percentiles) of log m = 8.7 +0.4 −0.5 , m = 0.2 +0.5 −0.1 and Γ = 2.1 ± 0.1. The solid (or solid-dashed) lines represent the median profiles, with the corresponding shaded areas showing the 16th and 84th percentiles of the distribution. The top panel of Fig. 1 shows how the standard µ = 0 law (e.g., SS73) results in f = f max at all radii (e.g., Svensson & Zdziarski 1994), whereas alternative viscosity laws (e.g., µ = 0.5 and µ = 1) show non-constant radial profiles for f : in the latter cases, the fraction of power dissipated in the corona is smaller in the regions strongly dominated by P rad . As it was shown in M03 in particular for the µ = 0.5 scaling, the higher suppression of the growth rate in P rad -dominated regions of the disk leads to such damped f -profiles. This directly influences the strength of the corona emission, as L 2 keV is proportional (through L X,tot ) to the product of f and Q + : Q + peaks at small radii in a very similar way across all models, therefore the ones with deeper f -profiles show flatter L 2 keV radial profiles and, hence, weaker coronae (bottom panel of Fig. 1). The exact shape of f (r) also affects the strength of the disk emission since the two are self-consistently coupled (see Eq. (5)).
We can also define the mean value of each f (r) profile (i.e. for each combination of m,ṁ and Γ): Colors are coded according to the choice of the viscosity law: stress proportional to P tot (P gas + P rad , µ = 0, black), to P gas (µ = 1, blue), or the geometric mean of the two (µ = 0.5, red). The continuous solid, or solid-dashed, lines represent the median profiles, with the related shaded areas showing the 16th and 84th percentiles, the scatter due to the range of sources (i.e. m andṁ) modeled. L 2 keV is proportional to the product of f and Q + (the accretion power per unit area). As Q + has very similar profiles across all models, those systems for which f is smaller produce weaker coronae in the central part. that is also a function of µ, α 0 and f max . Then, the mean value can be computed for the median f profiles in the examples in the top panel of Fig. 1: f median = 0.5, 0.13 or 0.05, for µ = 0, 0.5 and 1, respectively. Of course, within such a model the exact value of f median depends on its normalization f max , that is a free parameter in the model only bound to be <1. Nonetheless, simply from looking at f median as a function of µ (and from Fig. 1) we can see how, for the same set of inputs (e.g., m,ṁ), the different accretion prescriptions relate to the output corona luminosities: in a nutshell, going from µ = 0 to µ = 1 produces lower f median , thus weaker coronae. Changing µ also affects the logarithmic scatter in the radial profiles, from being absent in µ = 0 to increase with µ for µ 0 (see Fig. 1). The spread on a given f (r) is due to the scatter in m,ṁ and Γ, where the major role is played by the accretion rate (e.g., see Fig. 1 in M03). Crucially, f decreases with increasingṁ for all µ 0 models. This points in the same direction as the evidence of an X-ray bolometric correction (that is proportional to the inverse of f ) increasing with the accretion rate (e.g. Wang et al. 2004;Vasudevan & Fabian 2007, 2009Lusso et al. 2010;Young et al. 2010). This relation between f andṁ has also crucial implications for what our models predictions on the physical mechanisms behind the observed L X −L UV (see Sect. 3).

Local thermal stability
Before proceeding to a detailed observational test of the model, we briefly discuss here the stability issue for the various adopted viscosity laws. Jiang et al. (2016) showed that the presence of the iron bump in the opacity at ∼2 × 10 5 K stabilizes the flow in the disk regions around that temperature, where the cooling term has a different dependency and thermal runaway is avoided (Grzȩdzielski et al. 2017b). In the top panel of Fig. 1, a solid median line for the f -profile represents (thermally) stable regions of the median test source, whereas a dashed line highlights the instability regions. The vertical dot-dashed lines show instead where the median transition radius, from P radto P gas -dominated regions, lies. This highlights that, for the median test source, the stability region extends also well within P rad -dominated regions of the disk, confirming previous results (Jiang et al. 2016;Grzȩdzielski et al. 2017b). More quantitatively, we computed the thermal stability balance (e.g. Pringle 1976) for each test source (m,ṁ) at all radii with varying viscosity laws. The new stability regions in the inner P rad -dominated portions of µ = 0 and µ = 0.5 disks are ubiquitous, but they appear at different radii according to where the disk reaches the temperatures around the iron bump in κ (see also Fig. A.1). The µ = 1 case, as it is well known (e.g. Lightman & Eardley 1974), is stable throughout.

Predictions of the model on the L X −L UV relation
In this section we aim to test our disk-corona model (presented in Sect. 2) against the observed L X −L UV , a robust observable linked to the disk-corona physics. Before performing a more quantitative observational test (Sect. 4), we here outline the predictions of our model concerning the disk-corona energetics and the expected impact of our accretion prescription on the L X −L UV .
The schematic illustration in Fig. 2 summarizes the qualitative take-home messages of this work. The observed L X −L UV states that going from a lower to a higher accretion regime, the luminosity of the corona increases less than the disk luminosity, resulting in a slope smaller than one in the log-space. In our model for the disk-corona system, the luminosity outputs are directly modified by the viscosity prescription in the flow, determined by the parameter µ, and by the fraction of accretion power going into the corona, f (see Table 1 for a summary on the model parameters). Among all scenarios spanned by these two main unknowns, the qualitative behavior of the accretion disk-corona system, along its radial extent, is similar: higher accretion states have a more powerful disks and coronae, but wider P rad -dominated inner region and, only for modified viscosity prescriptions (i.e. µ 0), lower relative contribution of the corona to the total luminosity (see the upper diagram in Fig. 2). Thus, our model can provide a simple explanation for the observed slope of the L X -L UV relation, bridging in a simple but effective way the gap between the observed X-to-UV energetics and some aspects of MRI simulations. Changing µ not only affects the disk thermodynamics, but also changes the amount of power carried away by the corona (see Fig. 1). A constant radial profile for f (e.g., Svensson & Zdziarski 1994; here µ = 0) would naturally result in a L X −L UV close to a one-to-one relation. On the contrary, the alternative viscosity prescriptions, that we identify with µ 0, inherently result in a different diskcorona energetic coupling for varying accretion rates: in particular, higherṁ yield more damped f −profiles (see also M03). In this scenario, the outcome would be a slope of the L X −L UV that is smaller than one (see the lower diagram in Fig. 2).
It is worth stressing that it is only the relative fraction f that is more suppressed in the inner regions of systems with a modified viscosity and not the X-ray emission per se. Regardless the underlying assumption of a plane-parallel geometry for the diskcorona system (see Eq. (5)), the X-ray emission peaks in the innermost radii (e.g., see the µ = 0.5 case in Fig. 1). This will be addressed in details in Sect. 5.1. Figure 3 shows an example of a mock L X −L UV from the model realizations, with f max = 0.5 and different values of µ as in Fig. 1, with the same color coding. For the connected solid points, a single typical mass (log m = 8.7) is adopted, with increasingṁ = 0.03, 0.07, 0.17, 0.42, 1 (from left to right in L UV ). The relation is linear for a given mass, with the dashed lines indicating a slope of one to guide the eye. As qualitatively shown in the illustrative Fig. 2, models where f is constant in radius and in accretion state (black points) yield a slope close to one (even higher for the single mass); instead, alternative viscosity prescriptions (red and blue), that change the disk-corona energetic interplay via f (r), show a flatter slope. The underlying transparent points show the mock L X −L UV for a distribution of log m = 8.7 +0.4 −0.5 ,ṁ = 0.2 +0.5 −0.1 and Γ = 2.1 ± 0.1, that follows the typically observed objects (see Sect. 4).

Observational test: modeling the L X −L UV
In the previous section we qualitatively outlined what are the physical mechanisms identified as the origin of the observed slope smaller than one. A more quantitative test is needed to thoroughly investigate all aspects of the observed L X −L UV , including its normalization and intrinsic scatter. The exact value of the slope given by the models not only depends on the unknowns µ, f max and α 0 , but also on the details of the distributions of m,ṁ, Γ that are adopted for the calculations. Nonetheless, not all combinations of these three parameters are observed, because they do not exist in nature or we are biased against their detection (for example, the black masses of AGN follow a given mass distribution, some some mass ranges are more probable than others in any observed sample). That is why we select a reference sample of radiatively-efficient broad-line AGN (Sect. 4.1) and model the most likely values for m,ṁ, Γ for each source individually, based on the available data.

The reference sample of broad-line AGN
We built our reference sample starting from the 1787 AGN within the XMM-XXL north survey (Pierre et al. 2016) identified as broad-line AGN (BLAGN) by the Baryon Oscillation Spectroscopic Survey (BOSS) follow-up (Menzel et al. 2016;Liu et al. 2016b, hereafter L16). The X-ray spectral analysis on these sources was performed in L16 with the Bayesian X-ray Analysis software (BXA, Buchner et al. 2014), providing N H values, photon indexes (Γ) and the rest-frame 2−10 keV intrinsic luminosities (L 2−10 keV ). Furthermore, single-epoch virial black hole masses (M BH , e.g. Shen et al. 2008) and continuum luminosities (at 1350, 1700, 3000 and 5100 Å) were obtained on the BOSS spectroscopy with a fitting pipeline (for details, we refer to L16; Shen & Liu 2012). Luminosities were computed in L16 assuming H 0 = 70 km s −1 Mpc −1 , Ω m = 0.27 and Ω Λ = 0.73 2 .
Then, we applied some cleaning criteria to avoid, as much as possible, imprecise estimates for the intrinsic (accretionpowered) L X and L UV and to remain consistent with what is computed by the model. Firstly, among the monochromatic luminosity values available in the optical-UV from L16 we adopted L 3000 Å , obviously inducing a redshift cut in the sample (see Fig. B.1). Model wise there is no difference in computing L 3000 Å or the more standard L 2500 Å , and Jin et al. (2012) showed compatible correlations between the X-ray luminosity and each wavelength of the optical spectrum, although their coverage starts from 3700 Å. We verified a posteriori that this choice does not affect significantly the slope of the L X −L UV or the conclusions of our work.
Secondly, despite being defined as BLAGN, Liu et al. (2018) found that a fraction of these sources shows continuum reddening probably due to intervening dust along the line of sight, not accounted for by our model. Liu et al. (2018) defined a slope parameter α for the optical-UV continuum, to discern between the reddened sources and the bulk of blue BLAGN at each redshift. The contamination from extinction at L 3000 Å was minimized by conservatively selecting sources with α < −0.5 (see Liu et al. 2018, their Fig. 2).
Moreover, the XMM-XXL survey has a typical exposure time of ∼10 ks per pointing (L16; Pierre et al. 2016). Here, the analysis was restricted to sources with at least 10 counts in the EPIC-pn (Strüder et al. 2001) and EPIC-MOS (Turner et al. 2001) cameras on board XMM-Newton (Jansen et al. 2001), to exclude sources with extremely low-quality X-ray spectra. Table 1. Summary of recurring model parameters.

Parameter Definition Comments
Regulates the scaling between the magnetic stress µ = 0 → P mag ∝ P tot µ and the thermal pressure: Proportionality constant of the viscosity law (see above) Then, to exclude data contaminated by X-ray absorption, not accounted for in our modeling, we conservatively selected only sources in which the 84th percentile of the N H posterior distribution was smaller than 10 21.5 cm −2 , a value typically adopted to distinguish X-ray obscured and un-obscured sources (Merloni et al. 2014;see also Della Ceca et al. 2008).
Finally, we take L 2 keV as reference for the corona emission in the L X −L UV . In the model, we computed mock L X,tot with no reflection, so that we could easily extrapolate L 2 keV assuming a simple power-law spectrum. However, in L16 the reflection component was also included in the calculation of L 2−10 keV , as it is usually observed both in low-z (e.g. Nandra et al. 2007) and high-z (e.g. Baronchelli et al. 2018) spectra (see the average-AGN model in Buchner et al. 2014). Therefore, we consistently excluded from the analysis all the sources with a significant reflection component: given the high errors of the typical log R fit in L16 3 , we included only sources in which the 16th percentile was <−0.2 and the 84th was <0.5. We note that L16 included in the fit also a scattering contribution from ionized material inside the angle of the torus (see Buchner et al. 2014), although the fit normalizations are on the order of 10 −4 with respect to the main power-law component.
The final cleaned subsample, to which we will refer as XMM-XXL, consists of 379 sources with observed m (with 3 R is the ratio of the normalization of the reflection component with respect to the power-law component. median log m = 8.7 +0.4 −0.5 ), L 3000 Å , Γ (with median Γ = 2.1 ± 0.1) and L 2 keV . In Fig. 4 we show the log L X = α + β log L UV relation, with the best-fit linear regression given by: . For the main scope of this paper, it is sufficient to have a reference sample cleaned in accord with the physics described within the model.

Methodology of the observational test
In our modelṁ = λ edd = L bol /L edd , although we do not take as reference alsoṀ or λ edd from L16: the former is interpolated from the mass and a monochromatic optical luminosity (Davis & Laor 2011), while the latter depends on a A135, page 6 of 16 disk-luminosity estimate via L bol . Both approaches are based on standard-disk assumptions or calculations and using those values within our non-standard disk models would be an inconsistency. One can also estimate L bol applying bolometric corrections (BC) to the observed monochromatic optical-UV luminosities (Richards et al. 2006;Runnoe et al. 2012), although the many uncertainties in play (Krawczyk et al. 2013;Kilerci Eser & Vestergaard 2018) and the high scatter in the BCs (Richards et al. 2006;Lusso et al. 2012) discouraged us in relying on this approach. Then, for every source we iteratively obtain theṁ value yielding a model L 3000 Å consistent with the observed one within its errors (typically ∼0.01 dex). This approach is similar to the interpolation method put forward by Davis & Laor (2011), although we do it consistently for each different model, which is given by a choice of µ, α 0 and f max . The methodology then consists in fixing µ, α 0 and f max (see Table 1 for a summary on the model's parameters), which will be referred to as the model choice, within a discrete 3D grid in µ = [0, 0.2, 0.4, 0.5, 0.6, 0.8, 1], α 0 = [0.02, 0.2] and f max = [0.1, 0.2, 0.3, 0.5, 0.7, 0.9, 0.99]. Then, we take as input m, Γ and L 3000 Å from the observed data, allowing us to solve the equations of the model for each source and computeṁ and L 2 keV values (see Sect. 2). For each observed source of the reference sample, every model in the 3D grid can provide a mock entry for the L X −L UV . A proper comparison requires uncertainties to be assigned on the mocks, as the observed m, Γ and L 3000 Å come with their own measurement and systematic errors, where obviously the ∼0.4−0.5 dex systematics in the mass estimates (e.g. Shen 2013, and references therein) play the dominant role. As it is mentioned above, mock L 3000 Å values converge to the related observed quantities within their errors, hence we conservatively fixed the mock δL 3000 Å at the 90th percentile of the uncertainty distribution in the observed L 3000 Å (i.e. ∼0.03 dex). In order to compute uncertainties forṁ and L 2 keV , we ran each model 200 times on the same source, extracting the input values (m, Γ and L 3000 Å ) from a normal distribution with mean and standard deviation taken from the observed quantities and their errors. Then, the uncertainty onṁ and L 2 keV is taken from the dispersion of the 200 runs.

Results of the observational test
For all the models on the discrete 3D grid in the µ, α 0 and f max parameter space (see Sect. 4.2 and Table 1), we fit the L X -L UV distribution with a log-linear relation log L X = α + β log L UV . Three examples are shown in Fig. 5 for µ corresponding to the known analytic viscosity prescriptions (see Sect. 2). Ideally, a model should reproduce the observed L X −L UV in both normalization and slope. However, we can start decomposing the problem in two parts: a good match in the normalization ( α) would state that globally, for a given optical-UV luminosity distribution, the modeled corona emission was strong enough (see Sect. 5.1); instead, if the slope ( β) is matched, then the model accurately describes how the coronal strength varies from lowlyto highly-accreting sources (see Sect. 5.2). Moreover, as it can be seen from the examples in Fig. 5, our models come with their one intrinsic scatter, given by different m and Γ at a fixedṁ. This provides precious insights on the nature of the total observed scatter (see Sect. 5.3). define a score for the goodness of match: The r 2 score is computed drawing 1000 random samples from the observed log L 2 keV within their errors (i.e. y data,i ), and 1000 random regression lines from emcee's chains on the mock (i.e. y mock i ). Then, the median and the 84th-16th inter-quantile range are quoted from the resulting distribution of 1000 r 2 i scores. Negative scores indicate the data are poorly reproduced by the model; an r 2 = 0 would be obtained by a constant value corresponding to the mean of the observed log L 2 keV distribution. We can put a quality threshold and keep all the models that yield a positive score.
The r 2 score as a function of f max is shown in Fig. 6, where the choice of µ is color coded and the additional dependency on α 0 is represented with varying line-types (it is minor or absent, as in µ = 0). We show for simplicity only values of µ corresponding to the known analytic viscosity prescriptions (see Sect. 2). The other values used would accordingly show intermediate results.
For each viscosity law there is a preferred f max , that fixes the maximum coronal strength in a model. Models with higher µ need higher normalization f max , since they have a comparably weaker X-ray emission, in accord with their lower f median (see Fig. 1 and Sect. 2.1). Nonetheless, the law correspondent to µ = 1 does not produce adequately strong coronae even with f max = 0.99 and can be ruled out (see right image in Fig. 5). Furthermore, µ = 1 produces a radially flatter X-ray emission profile (see bottom panel of Fig. 1), in contrast with observations that hint for coronae peaking in the inner radii (e.g. Mosquera et al. 2013;Reis & Miller 2013;Wilkins et al. 2016). We explore this behavior more quantitatively in the top panel of Fig. 7 showing how the radius of the annulus at which the 2 keV emission peaks (r peak , or at which it is 90% of the total, r 90 ) varies with µ: as µ increases, most of the corona emission comes from annuli placed at larger and larger radii.
We also verified that our results do not depend on the sample adopted as reference. We performed the same analysis with the RM-QSO sources (Liu et al., in prep.;Shen et al. 2019), on which a similar analysis was performed and on which we applied compatible cleaning criteria and methodology, as described in Sects. 4.1 and 4.2. The results are shown in the r 2 score plot ( Fig. 6) with shaded areas, color coded for µ in the same way and using only α 0 = 0.02. There is generally a good agreement between the two samples, suggesting that our results are not dependent from the different data used.
It is worth stressing that the f max value at which each µ (possibly) matches the observed normalization is degenerate with the assumptions on the accretion efficiency and on the product η (1 − a disk ). Namely, higher accretion efficiencies and/or a higher fraction of the coronal emission beamed away from the disk would increase the normalization of the L X -L UV relation, and shift all curves of Fig. 6 to the left. This will be further examined in Sects. 5.6 and 5.7.
5.2. The slope of the L X −L UV Figure 6 allows us to track the models (i.e. combinations µ, α 0 and f max , see Table 1) that broadly reproduce the normalization α of the observed L X −L UV . Nonetheless, obtaining the correct normalization is simply a weighting exercise of the energetic outputs of the disk and the corona. It is the slope that carries the exact information on how the disk-corona interplay changes across the different accretion regimes of bright radiatively-efficient AGN (see Sect. 3). This would require a precise knowledge of the true slope of the L X −L UV . The observations suggest a value around ≈0.6 (Lusso & Risaliti 2016; LR17; our Appendix B) and we can acknowledge this value as reference. Our methodology, however, can be regarded as dataindependent, and it would applicable even if future works will update the current knowledge on the exact value of the slope.
In the middle panel of Fig. 7 we show how the modeled slope of the of the L X −L UV gets closer to the observed one for increasing µ (i.e. for more damped radial f -profiles), for a fixed α 0 = 0.02 and using only the f max corresponding to the highest r 2 -score. This is because models with increasing µ have higher logarithmic scatter in f (r), meaning that going from lowly-to highly-accreting sources the span in f i is larger, with high-ṁ objects having comparably weaker X-ray emission with respect to low-ṁ companions (see Sect. 3). We show this for µ = 0, 0.5 and 1, respectively 4 : log f = log f max log f = (−1.12 ± 0.24) − (0.15 ± 0.02) logṁ + (0.05 ± 0.03) log m log f = (−1.82 ± 0.36) − (0.27 ± 0.03) logṁ + (0.07 ± 0.04) log m where the steepest dependency fromṁ is obtained for larger µ. This test points in the same direction as the evidence of an X-ray bolometric correction increasing with the accretion rate (e.g. Wang et al. 2004;Vasudevan & Fabian 2007, 2009Lusso et al. 2010;Young et al. 2010), although we refrain to compare this observable with our regressions (e.g. Wang et al. 2004;Cao 2009;Liu & Liu 2009;You et al. 2012;Liu et al. 2012Liu et al. , 2016a, due to the many more uncertainties in play when deriving bolometric luminosities in comparison to the quantities entering in the L X −L UV (see the discussion in Sect. 4.2).

The scatter of the L X −L UV
The observed scatter of the L X −L UV for the sample used in this work is σ intr = 0.27 ± 0.01 (Sect. 4.1). As a matter of The green shaded area qualitatively shows the inner radii, where the bulk of X-ray emission is supposed to come from according to X-ray reverberation and micro-lensing. For increasing µ, the X-ray emission profile peaks at larger radii. Middle panel: for increasing µ the models obtain a slope of the L X −L UV closer to the observed one. The dark-green area represents the reference slope of the cleanest XXM-XXL (Appendix B), while the light-green refers to the slope quoted in LR17. Bottom panel: intrinsic scatter of the mock L X −L UV relations as a function of µ. The green area represents a tentative upper limit of the true scatter (Lusso & Risaliti 2016;Chiaraluce et al. 2018), that is only due to the physical properties of AGN. For simplicity, all panels show only the results obtained with a single f max , corresponding to the highest r 2 -score (e.g., Fig. 6), and fixed α 0 = 0.02. fact, this value represents an upper limit to the intrinsic dispersion inherent to the physics of the system, as the observed scatter is affected by a combination of instrumental and calibration issues, UV and X-ray variability, non-simultaneity of the multi-wavelength observations. A lot of effort has been put into trying to quantify as accurately as possible all these contaminants (e.g. Vagnetti et al. 2013;Lusso 2019, and references therein), with claims that the intrinsic scatter in the L X -L UV relation is smaller than 0.18−0.20 (Lusso & Risaliti 2016;Chiaraluce et al. 2018). Any successful model should be able to reproduce such a low scatter.
From the examples of mock L X −L UV relations plotted in Fig. 5, it can already be seen that our models come with their one intrinsic scatter. In our methodology (Sect. 4.2), the modeled  [0, 0.2, 0.4, 0.5, 0.6, 0.8, 1], as shown in the legend. For simplicity, we report for each µ only results obtained with a single f max , corresponding to the highest r 2 -score, and fixed α 0 = 0.02. Models that reproduce the observed slope α are also the ones that show weaker coronae (lower normalization α).
m was tuned to the observed L 3000 Å , hence the intrinsic scatter of the mock L X −L UV relations is simply the dispersion of the modeled L 2 keV , at a givenṁ, due to different m and Γ. We show this more quantitatively in the bottom panel of Fig. 7. The models dispersion varies with µ because changing the viscosity law induces a different logarithmic scatter in f (r) (see Fig. 1) and it also affects the distance (in gravitational radii) from which the bulk of the L 2keV is coming (see top panel of Fig. 7). The resulting σ intr of the models is likely a complex combination of these (and possible more) factors. All the models, with the exception of µ = 0, lie below the available observational constraints (Lusso & Risaliti 2016;Chiaraluce et al. 2018) of 0.18−0.20. This is another successful prediction of our model (see Sect. 3).

5.4.
A complete picture: the slope-normalization plane of the L X −L UV In the previous sections, we decomposed the match in either normalization or slope to have a better understanding on how our disk-corona models can relate to the observed L X −L UV . However, the goal would be to have a model that can fully encompass these observables. Hence, in Fig. 8 we display 1-, 2-and 3-sigma contours in the slope-normalization plane ( β − α) of the L X −L UV for both data and models. All regressions were performed with emcee normalizing both L X and L UV to the median value of XMM-XXL. The data contours are related to the cleanest XMM-XXL version (Appendix B) and to the RM-QSO sources 5 . Model contours are shown for µ = [0, 0.2, 0.4, 0.5, 0.6, 0.8, 1] using a single f max , corresponding to the highest r 2 -score (e.g., Fig. 6) for each µ, and a fixed α 0 = 0.02, for simplicity. Figure 8 shows that models reproducing the observed slope, namely the ones with higher µ (as in middle panel of Fig. 7), are also the ones that show weaker coronae (lower normalization α) and overly extended L 2 keV -emission (i.e. higher r peak and r 90 , top panel of Fig. 7).
5.5. The 3D plane: L X vs L UV vs m As shown by LR17, the L X −L UV relation for AGN is rather a three-dimensional problem, with the mass (or its proxy given by the full-width half-maximum of broad emission lines) playing a significant role as well. The observed L X −L UV − m plane from XMM-XXL can be fit by: log L 2 keV − 25 = (−2.28 ± 0.08) + (0.57 ± 0.02) (log L 3000 Å − 25) The comparison in the 3D plane states that the exact dependency is not obtained by any of the models, with µ = 1 being the closest in qualitatively retrieving the coefficients for L 3000 Å and m. We note that the mass is taken from the observations, thus this mismatch states that the luminosities in the model do not depend on the mass in the correct way.

The impact of the accretion efficiency
Throughout this work we adopted an efficiency 0 = 0.057, typical of non-rotating black holes (e.g. Shapiro 2005), for simplicity. Nonetheless, a high spin seems to be preferred to model the blurred relativistic iron line, detected both in the local Universe (Nandra et al. 2007;Reynolds 2013) and up to z ∼ 4 (e.g. Baronchelli et al. 2018). Moreover, flux-limited samples are known to be biased in preferentially detecting highspinning black holes (Brenneman et al. 2011;Vasudevan et al. 2016), simply because they are brighter than their non-rotating analogous (see Reynolds 2019). Then, we tested the model using maximally-spinning black holes, with radiative efficiency 0.3 and ISCO down to r 0 = 1.24r g (Thorne 1974). This has a major impact on the normalization axis of the L X −L UV . Everything else in the source being equal, in a spinning black hole matter can be accreted down to smaller distances with respect to their non-rotating companions, thus the accretion power in the system is much higher. As a matter of fact, changing the radiative efficiency has an impact on the numerical equation that regulates f (r): for the same m andṁ and r > 3 the values of f is higher, and the transition radius between P rad -and P gas -dominated regions moves at lower radii. This self-consistently affects the disk equations via the (1 −f ) factor (see Appendix A), hence the surface temperature is decreased at higher radii, where most of the disk emission at 3000 Å comes from. Then, the modeledṁ value needed to match the observed L 3000 Å is higher (see Sect. 4.2) and, consequently, L 2 keV ∝ f Q + is higher. In Fig. 9 we show the model contours in the correlation slope-normalization plane computed for both low and high radiative efficiency, for µ = 0.4, 0.5 and 0.6 only.
Interestingly, maximally-spinning sources yield a better match with the data contours, in particular for the viscosity law µ = 0.5, with f max = 0.9. For instance, Fig. 10 shows how the data and this high-spin model compare in the L X −L UV plane. We want to stress that using only a maximum spin for all sources is an extreme measure, but since the (unknown) observed spin distribution is likely dominated by high-spin values (Reynolds 2019), model contours of a more realistic diverse population of high-spinning sources would be closer to the high-efficiency ones in Fig. 9 rather than to the spin-zero case. We also note that, even if the modeled coronae would be somewhat weaker using a realistic spin distribution, with respect to the maximum-spin case, the model with µ = 0.5 can still be realized with a higher f max = 0.99. Thus, we speculate that the new empty red contours in Fig. 9 consist in a fair approximation of a realistic high-spin population model. The tension with the observed L X −L UV would be significantly relaxed.

The impact of the downward scattering component
The results shown in Fig. 6 are also degenerate with the assumptions on the value of the product η (1 − a disk ), that is on the assumed downward component of the X-ray emission (η) and on the disk albedo. The adopted value of η = 0.55 is typical for anisotropic Comptonization in a plane-parallel corona (Haardt & Maraschi 1993), although it is unclear how much it would change in different geometries or prescriptions. In a patchy corona (Haardt et al. 1994) η would unlikely part significantly from the one in the slab case. The only major difference would rather involve the transmission or absorption by the corona of the radiation reflected by the disk. However, we conservatively excluded from the reference sample adopted in the observational test all the sources with a non-negligible reflection component detected (see Sect. 4.1), allowing us to avoid its complicated modeling. In an outflowing corona (e.g. Beloborodov 1999; Malzac et al. 2001), the ratio between the downward and A135, page 10 of 16  Fig. 10. L X −L UV relation for the high-spin model with µ = 0.5, α 0 = 0.02 and f max = 0.9 (empty red points, corresponding to the empty red contour in Fig. 9), with the red line showing best fit slope from emcee. The connected filled points (dark red) show the single-mass trend (log m = 8.7) for varying accretion rate (0.03, 0.07, 0.17, 0.42, 1). For a comparison, the black contour shows where the data lie in the plane, with the related best-fit slope (black line).
the upward flux decreases with the bulk velocity of the corona (e.g. Janiuk et al. 2000). We tried to quantify possible offsets in the β-α plane due to different values of η, ranging from 0.6 (slightly enhanced downward scattering) to 0.4 (reduced downward scattering, roughly approximating an outflowing corona with β bulk ≈ 0.1 − 0.2, e.g. Janiuk et al. 2000). We show this in Fig. 9 for the µ = 0.5 case, with dark-red density spots (η = 0.4, 0.5 and 0.6 from higher to lower α, respectively). Changing the downward component by ∆η ∼ 0.1 induces a significant offset of ≈0.1 dex in α and a minor change in β.

Discussion
The L X −L UV relation has been studied for decades (starting with the better-known α OX parameter, Tananbaum et al. 1979), its robustness used for bolometric estimates (e.g. Marconi et al. 2004;Hopkins et al. 2007;Lusso et al. 2010) and recently even for cosmology (Risaliti & Lusso 2015. Nonetheless, there is currently no solid and exhaustive physical explanation for it. In Sect. 3 we outlined the qualitative predictions of our model and in Sect. 5 we obtained that concordance with current data can be obtained with a modified viscosity prescription in the accretion flow (µ = 0.5), provided the spin of the sources is high. Here, we briefly discuss whether other competing analytic disk-corona models succeed or not and then we try to investigate the impact of the assumptions in our model on the results.

Comparison with other models
LR17 tried to explain this relation with a very simplified, but effective, toy-model. Most of their assumptions are in common with our work (see Sect. 6.2), although our treatment is more complete and physically motivated. The assumption of the MRI amplifying the magnetic field to a lesser extent in P rad -dominated regions (Blaes & Socrates 2001;Turner et al. 2002;M03) is taken to the extreme with a step function for the f -profile: all the accretion power is emitted by the disk in P rad -dominated regions (i.e. f (r rad ) = 0), whereas it is equally distributed between disk and corona in P gas -dominated regions (i.e. f (r gas ) = 0.5). The resulting predicted slope and normalization of the L X −L UV are claimed to be consistent with the observations. The former can be confirmed by our analysis, as their f (r) step-function is nothing but an extremely damped f (r) beyond µ = 1, whose mock slope of the L X −L UV was the closest to the observed one. In the latter case, their match in normalization might be an involuntary artifact: with respect to the power transferred to the corona f , the observed luminosity is roughly halved if a downward scattering component is included (i.e. f (1 − η), with η ≈ 0.5). We verified this running our model with µ = 0 and α 0 = 0.02, forcing f = 0 in the P rad -dominated region and fixing both f = 0.50 and f = 0.99 in P gas -dominated radii. In Fig. 11 we show the related contours in the β− α plane along with our results of Fig. 8. This confirms that their step f -profile results in a slope consistent with the observed value, albeit producing overly weak coronae (too low normalization in the β− α). Hence, their toy-model does not reproduce the L X − L UV . Moreover, the X-ray emission from their toy-model inevitably peaks at the transition radius between P rad -and P gas -dominated regions. Indeed, their model with f gas = 0.99 yields r peak = 142 438 51 and r 90 = 790 1490 445 (i.e. produces extremely extended coronae). Kubota & Done (2018) coupled an outer standard disk with an inner warm Componising region, that produces the soft X-ray excess, and an innermost hot corona for the hard X-ray continuum. Their model fits remarkably well the broadband continua of three sources spanning a wide range of accretion rates. They also claim to reproduce the observed L X −L UV , using both the regression line and data points from LR17, although only displaying all the possible sources modeled within a grid of m = 10 6 −10 10 andṁ = 0.03 − 1 (their Figs. 7 and 8). Nonetheless, first-order normalization matches, even with m andṁ spanning within typical values, can be misleading. A more conclusive test would be, as we do, to match mock and data sources one by one.

Further assumptions and theoretical uncertainties
Only models with µ 0.4 are able to reproduce the observed normalization within the range of possible f max values, whereas for µ 0.5 they are off by 0.1−0.2 dex along the normalization. In Sects. 5.6 and 5.7 we showed how a higher accretion efficiency and/or a different downward scattering component may affect our results in the slope-normalization plane. Their impact would be significant and can possibly ease the tension between data and models: high-spin black holes and/or moderately outflowing coronae would be consistent with the observations. We now try to investigate some other simplifications of our model, likely to have a minor or less quantifiable effect on our conclusions.

Soft X-ray excess and thermal instability
The XMM-XXL L 2 keV value was interpolated from the L 2−10 keV fit in L16 after excluding sources with high reflection fraction (Sect. 4.1). The impact of the soft X-ray excess component can be considered negligible in that energy range, thus data points in the L X −L UV are likely not contaminated. However, our models do not include a soft X-ray excess generation mechanism, the monochromatic L 2 keV being extracted from a powerlaw spectrum within 0.1−100 keV. If a significant fraction of the power dissipated in the corona is actually used by a different mechanism producing the observed soft-excess, namely from a warm corona (e.g. Petrucci et al. 2018;Kubota & Done 2018;Middei et al. 2019), the mock L 2 keV would be overestimated to an unclear extent. Nonetheless, if the soft X-ray excess is produced by blurred relativistic reflection (e.g. Crummy et al. 2006;Garcia et al. 2019), the influence of this component on our analysis would have been excluded with our selection criteria (Sect. 4.1).
In Sect. 2, we briefly addressed the disk-instability problem (see Fig. 1, top panel) and despite the local stabilizing effect of the iron bump in the opacities, disks with µ = 0 (e.g., SS73) and 0.5 (e.g., M03) are globally unstable in P rad -dominated regions. An intriguing question may be whether the unstable regions in the disk are responsible for generating the soft-excess, possibly within inhomogeneous flows (e.g. . As a matter of fact, the higherṁ the wider the region where P rad dominates and the higher the soft-excess strength (e.g. Boissay et al. 2016). Nonetheless, a more thorough investigation of this scenario is beyond the reach of this paper.

Magnetically-dominated disks
In our model the stress tensor is dominated by Maxwell stresses as confirmed by simulations (e.g., Hawley et al. 1995;Sano et al. 2004;Minoshima et al. 2015), although the magnetic pressure is bound to be only a fraction of the product P µ gas P 1−µ tot via α 0 at the mid-plane. However, there are theories postulating disks that are P mag -dominated also in the denser regions (e.g. Begelman & Silk 2017, and references therein) and not only in the upper layers (e.g. Miller & Stone 2000), possibly solving a few long-standing issues of the standard accretions disk theory (Dexter & Begelman 2019). Simulations indeed showed that P mag can become an important competitor in supporting the disk vertically (Bai & Stone 2013;Salvesen et al. 2016), although heavily depending on the strength of the net vertical magnetic field, the origin of which is not fully understood, yet. If this imposed net vertical field is small (if β 0 = P tot /P mag >> 1), the buoyant escape of the toroidal component, amplified by MRI, is faster than its creation and a disk-corona system consistent with our model is formed. However, the evidence of disks that are magnetically-dominated even at the mid-plane is supported by Jiang et al. (2019), that recently performed a global 3D radiation-MHD simulation of two sub-Eddingtion flows. The structure of their simulated disks is significantly different from the standard SS73 model and reaches a complexity that our simplified prescriptions are not able to grasp. On the other hand, these simulations could not produce spectra and luminosities, yet. We here rely on the assumption that the energetics of P magdominated disks are not significantly different from standard thin disks at radii larger than ∼10r g (e.g., see Sdowski 2016).

Winds and outflows
In order to see if any known broad absorption line (BAL) quasars were present in our sample, we cross-matched the XMM-XXL catalog (L16; Menzel et al. 2016) with SDSS-DR12 (Pâris et al. 2017), that flagged 29580 BAL QSO after visual inspection.
Only two sources among the 379 used in our analysis were flagged, although they were both assigned zero indexes in the common metrics used for a more quantitative measurement of the BAL properties (Pâris et al. 2017). Hence, our sample has no contaminations from known BALs, although we can investigate the possible impact of un-modeled wind-dominated objects on our work. For instance, Nomura et al. (2018) recently developed a disk model compensating for the mass-loss rates of UV-driven winds, while consistently adjusting the temperature and emission of the underlying disk. They referred to a future work for a more complete modeling of the inner radii and the hard X-ray emission, but the influence on L 3000 Å values seems already significant, providedṁ 0.5. Since winds appear to act only from moderate to Eddingtonṁ, neglecting their presence would have an impact on the modeled L X −L UV slope. The wind carries away kinetic energy reducing the disk emission accordingly, thus for a given observed high L 3000 Å , our no-wind model would underestimateṁ for the possible outflowing sources contaminating our sample.

The larger-than-predicted disk argument
One of the most studied issues of the standard SS73 disk model is that observed sizes appear to be larger than expected at optical-UV wavelengths, using both microlensing effects (e.g. Morgan et al. 2010;Blackburne et al. 2011;Jiménez-Vicente et al. 2012) and flux variability lags across multiple bands in the so-called reprocessing scenario (e.g. Edelson et al. 2015;Fausnaugh et al. 2016Fausnaugh et al. , 2018Jiang et al. 2017b;Cackett et al. 2018;McHardy et al. 2018), in which often a compact X-ray emitting region (e.g., a lamppost corona) irradiates the disk inducing light-travel lags in the UV-optical bands. However, even combining all these results discordant with the theoretical predictions is not trivial (Kokubo 2018), particularly if different techniques are used (see Moreno et al. 2018;Vio & Andreani 2018). What is more, there are also numerous studies finding consistency with the sizes predicted by the standard SS73 theory (e.g. McHardy et al. 2016;Mudd et al. 2018;Yu et al. 2018;Edelson et al. 2019;Homayouni et al. 2019), thus we do not consider necessary to use the larger-thanpredicted argument to abandon all the standard prescriptions, yet.

No-torque inner boundary
For convenience, we adopted the no-torque condition with the stress vanishing at the inner edge. However, the presence of magnetic torques (Gammie 1999; would increase the disk effective temperature and the Q + emissivity in the innermost radii Dezen & Flores 2018) and, if applied to the disk only, it would cause instead a drop in the fraction f (Merloni & Fabian 2003). Without a proper MHD treatment, it is unclear how the modeled L 2 keV ∝ f Q + would be affected, and consequently the L X −L UV slope.

The vertical structure
Our model does not properly treat the vertical structure of the disk. The effective temperature is obtained from T eff ∝ T mid /τ 1/4 , where τ = h ρ κ assumes constant ρ and κ along the scale-height. Even keeping the approximation of a constant ρ, κ should change self-consistently with the decrease in temperature. A more thorough modeling of the disk vertical structure in supermassive A135, page 12 of 16 black holes was presented by Hubeny and collaborators, taking into account both scattering processes and free-free and boundfree opacities (Hubeny et al. 2000(Hubeny et al. , 2001. Their model also share some of our limits (e.g., stationary disk, α-prescription, no-torque boundary, vertical support from thermal pressure only), validating the comparison. The overall SED has lower (higher) fluxes at low (high) frequencies with respect to standard calculations, with the most significant impact on the modeling of the soft-excess . The computation of L 3000 Å should be affected in a minor way, with a small overestimation on the order of a color correction (e.g. Done et al. 2012), that is either roughly constant or weakly depending on m andṁ (e.g. Davis & El-Abd 2019). Our conclusions should not be significantly affected, although this would need to be improved for a proper SED modeling and time-lags predictions.

Conclusions
The gap between simulations and observations in AGN needs to be bridged and simplified, but motivated, analytic prescriptions still represent a powerful tool to explain the observed multi-wavelength scaling relations. For instance, the clear correlation observed between monochromatic logarithmic L X and L UV luminosities has been used for decades (in the shape of the more-known α OX parameter, Tananbaum et al. 1979) in many applications (even for cosmology, e.g. Risaliti & Lusso 2015. Despite this, a conclusive theoretical explanation for the observed correlation is still lacking. Being smaller than one, the observed slope indicates that, going from low-to high-accretion rate AGN, the X-ray emission increases less than the optical-UV emission. Any viable disk-corona model must be able to explain this. In this work, we tested a self-consistent disk-corona model (Sect. 2, see also M03) against the L X −L UV relation. We were able to identify the possible mechanism regulating the diskcorona energetic interplay, in terms of viscosity prescriptions (e.g., µ = 0.5) that naturally lead to an X-ray emission increasing less than the disk emission going to higher accretion rates (see Sect. 3).
We also put forward a quantitative observational test (Sect. 4), using a reference sample of AGN (Sect. 4.1) observed both in the (rest-frame) UV and in X-rays: taking from each source the observationally determined m,ṁ and Γ we were able to model an analogous mock object (Sect. 4.2) producing a set of mock L X −L UV . This allowed us to reach a deep understanding of the physics driving the slope, normalization and scatter of the L X −L UV (see Sect. 5).
We find that if the black-hole population is assumed to be non-spinning, results from this test are inconclusive: the viscosity prescriptions reproducing the slope of the observed L X −L UV relation, also produce overly weak coronae. Interestingly enough, the tension between the strength of the observed and modeled X-ray emission (i.e. in the normalization of the L X −L UV ) can be significantly relaxed adopting a more realistic high-spinning black-hole population and/or with moderatelyoutflowing coronae. We tested the former case adopting the efficiency (and the inner orbit) of maximally-spinning black holes, in which matter is able to accrete further into the potential well, resulting in a much higher accretion power and, consequently, in much stronger coronae (Sect. 5.6). Moreover, if the spin is high the X-ray emission profile peaks closer to the black hole, in even better agreement with X-ray reverberation and microlensing studies (Mosquera et al. 2013;Reis & Miller 2013;Wilkins et al. 2016). In particular, the disk-corona model testing maximally-spinning black holes with µ = 0.5 (i.e. magnetic stress proportional to the geometric mean of P gas and P tot , e.g. see M03), f max = 0.9, α 0 = 0.02 (see Table 1) provides the best match with the observations (Fig. 10), although the modeled slope is still somewhat larger than the observed one (Fig. 9). Going beyond this type of exercises, only 3D global radiation-MHD simulations will be able to better disclose the disk-corona physics (e.g. Jiang et al. 2017aJiang et al. , 2019, provided a clearer way of approaching the observations will be reached. We report the equations for h, mid-plane ρ (g cm −3 ), P (dyn cm −2 ), T (K, at the mid-plane) and the closure equation for The value of ξ can be obtained by studying the continuity of all the above quantities at the boundary between the radiation pressure-to the gas pressure-dominated regions. It corresponds to ξ 1.00k −1/3 0 . In Fig. A.1 we report examples of radial profiles for ρ, P tot , κ, h/r, T mid and T eff . Similar examples for f and L 2 keV are shown in Fig. 1. Once f max is fixed, the dominant variance among the models is given by the choice of the viscosity law (µ), while α 0 plays a minor role. This is shown in Fig. A.2, where profiles for f and L 2 keV show little difference in varying α 0 from 0.02 to 0.2.

Appendix B: The reference AGN sample
For the source-by-source modeling of XMM-XXL we used the 379 sources obtained following the methodology outlined in Sect. 4.2. In Fig. B.1 we show the distribution of L 2 keV (top panel) and L 3000 Å (bottom panel) in the luminosity-redshift plane of the 379 sources (red and blue respectively), with respect to the parent sample of BLAGN from L16 (black). The L X −L UV slope of this reference sample is 0.54 ± 0.02, from Eq. (12). This is incompatibly flatter than the values quoted in the recent literature, namely 0.64 ± 0.02 (Lusso & Risaliti 2016) or 0.63 ± 0.02 (LR17). The cleaning criteria applied in Sect. 4.1 were aimed to exclude low-quality data and to be consistent with the model, while in the above-mentioned literature the possible biases of flux-limited samples were treated carefully in order to reliably use quasars for cosmology (Risaliti & Lusso 2019).
We investigated whether this inconsistency in the slope would be bridged restricting the analysis to the brightest objects at all redshifts with a very crude and conservative selection. From the sensitivity curve of the XXL-N survey in the 0.5−10 keV band at half of the survey area (L16, their Fig. 3) we obtained the flux limit in that energy band. Then, we interpolated the flux limit at 2 keV using the mean photon index of the sample, obtaining the sensitivity curve shown in red in the top panel of Fig. B.1. In Menzel et al. (2016) a cut at r < 22.5 mag was applied. We converted this magnitude limit in a luminosity sensitivity only within 0.80 z 1.27, for which 3000 Å was actually detected in the r band. For different redshifts, we first computed a redshift dependent color correction for the other bands (u, g, i and z) performing a linear regression on the difference with the r-band magnitude. This provided a magnitude limit for L 3000 Å at all redshifts, consistently with the band in which that wavelength was actually detected, from which we obtained the related sensitivity line in the bottom panel of Fig. B.1. We then divided XMM-XXL in six redshift bins, making sure to have at least 30 counts per bin. For each bin, we excluded all the sources below the limits given by the sensitivity curves on both axis, evaluated at the maximum z of the bin to be conservative (Fig. B.2). The resulting cleanest subsample reaches accordance with the recent literature of the L X −L UV , with a slope of 0.59 ± 0.03.   The sensitivity surfaces at the minimum, median and maximum redshift of the bin are represented in red with a full area, a dashed line and a shaded area respectively. These surfaces are obtained from the sensitivity lines in Fig. B.1 at the above-mentioned redshifts. The sources above the shaded sensitivity area in each z-bin give the cleanest XMM-XXL sample.