Fitting AGN/Galaxy X-Ray-to-radio SEDs with CIGALE and Improvement of the Code

It is generally believed that the observed X-rays are a result of a "disk-corona" structure. The disk emits UV/optical photons, and a fraction of them are up-scattered to X-ray wavelengths by the high-energy electrons in the corona (i.e., inverse Compton scattering). The angular dependence of the X-ray emission is related to the detailed physical properties of the corona, such as shape and optical depth (e.g., Sunyaev & Titarchuk 1985; Xu 2015).

The observations of Liu et al. (2014) found that type 2 AGNs tend to systematically have lower intrinsic X-ray luminosity (L_X) than type 1 AGNs at a given [o iv] 25.89 μm luminosity. Assuming that the [o iv] emission from the narrow-line region (NLR) is isotropic, they interpreted their result as an indicator of AGN X-ray anisotropy, because type 2 AGNs have larger viewing angles (as measured from the AGN axis) than type 1 AGNs under the scheme of AGN unification (e.g., Antonucci 1993; Urry & Padovani 1995; Netzer 2015). Also, based on the X-ray and high-resolution mid-IR observations of nearby AGNs, Asmus et al. (2015) suggested that AGN X-ray emission might be anisotropic (see their Section 5.4 for details).

However, x-cigale assumes that AGN intrinsic X-ray emission is isotropic, and does not allow anisotropic modeling. This assumption could overestimate the X-ray emission for type 2 viewing angles.

In x-cigale, the AGN X-ray emission is modeled using the α_ox–L_ν,2500Å relation from Just et al. (2007), where L_ν,2500Å is the intrinsic-disk emission at a viewing angle of 30° ¹⁴ and α_ox is the AGN SED slope connecting L_ν,2500Å and L_ν,2keV . For type 1 AGNs, whose viewing angles are near 30° (Yang et al. 2020), the SEDs are similar for isotropic and anisotropic X-ray models in the framework of x-cigale. However, for type 2 AGNs, whose viewing angles are much larger than 30°, the isotropic and anisotropic models will predict significantly different X-ray emissions, at a given AGN power.

To test the effectiveness of x-cigale (isotropic X-ray emission), we use a spectroscopic type 2 AGN sample from the Chandra COSMOS-Legacy survey (Civano et al. 2016; Marchesi et al. 2016). We require these AGNs to have signal-to-noise ratio (S/N) > 3 in the hard band of Chandra (2–7 keV), and we apply absorption corrections to the hard-band fluxes based on the correction factors from Marchesi et al. (2016), because x-cigale requires that the input X-ray fluxes are intrinsic (Yang et al. 2020). ¹⁵ The absorption corrections from Marchesi et al. (2016) are based on a standard hardness-ratio analyses.

We remove sources with L_X < 10^42.5 erg s⁻¹, for which the X-ray emission might originate from normal galaxies rather than AGNs (e.g., Aird et al. 2017). We adopt the 14 broadband photometric data (u to Infrared Array Camera (IRAC) 8 μm) from the COSMOS2015 catalog (Laigle et al. 2016). We also use the Multiband Imaging Photometer for Spitzer (MIPS) 24 μm, Photodetector Array Camera and Spectrometer (PACS) 100/160 μm, and Spectral and Photometric Imaging Receiver (SPIRE) 250/350/500 μm photometry from the "super-deblended" catalog of Jin et al. (2018). There are a total of 296 type 2 AGNs, spanning a redshift range of 0.3–1.6 (10th–90th percentiles).

Our fitting parameters are listed in Table 1. For the star formation history (SFH), we adopt a delayed τ SFH model and a Bruzual & Charlot (2003) simple stellar population model. We adopt the Calzetti et al. (2000) galaxy attenuation law and the Dale et al. (2014) dust IR spectral templates. For AGN IR emission, we adopt the SKIRTOR clumpy torus model (Stalevski et al. 2012, 2016). We fix the torus half-opening angle to the default 40°, which is observationally preferred (e.g., Stalevski et al. 2016). Under this setting, there are four type-2 viewing angles (60°, 70°, 80°, and 90°) available in SKIRTOR, and we allow all of these values in our fits (see Table 1). The full SKIRTOR model set has another five parameters such as 9.7 μm optical depth and ratio of outer to inner radius. These parameters generally have minor effects on the broadband SED shapes (e.g., Yang et al. 2020), and thus we leave them at the default values to reduce the needed computing resources. In summary, we employ four templates (corresponding to different viewing angles) out of the total 19200 SKIRTOR models. In x-cigale, AGN X-ray and UV/optical emissions are related with the α_ox–L_ν,2500Å relation of Just et al. (2007), where α_ox is the UV/X-ray slope calculated at the typical AGN type-1 viewing angle of θ = 30° (Yang et al. 2020), i.e.,

$$\begin{eqnarray}&&{\alpha }_{\mathrm{ox}}=-0.3838\mathrm{log}\displaystyle \frac{{L}_{\nu ,2500{\rm{\mathring{\rm A} }}}(30^\circ )}{{L}_{\nu ,2{keV}}(30^\circ )},\end{eqnarray} \tag{ 1 }$$

where L_ν,2500Å and L_ν,2keV are the monochromatic AGN luminosities per frequency at rest frame 2500 Å and 2 keV, respectively. Although the α_ox–L_ν,2500Å relation is reasonably tight, it has a non-negligible intrinsic scatter of ≈0.1 in terms of α_ox (e.g., Steffen et al. 2006; Just et al. 2007). x-cigale considers the scatter by constructing different models around the α_ox–L_ν,2500Å relation, and the user can set the maximum deviation from the relation, $| {\rm{\Delta }}{\alpha }_{\mathrm{ox}}{| }_{\max }$ (see Section 2.2.3 of Yang et al. 2020 for details). In our fits, we set $| {\rm{\Delta }}{\alpha }_{\mathrm{ox}}{| }_{\max }=0.2$ (Table 1), about 2σ of the intrinsic scatter (e.g., Just et al. 2007). We set the photon index Γ = 1.8, a typical value for distant X-ray AGNs (e.g., Yang et al. 2016; Liu et al. 2017).

Table 1. Model Parameters for the Type 2 AGNs in COSMOS

Module	Parameter	Symbol	Values
Star formation history $\mathrm{SFR}\propto t\exp (-t/\tau )$	Stellar e-folding time	τstar	0.1, 0.5, 1, 5 Gyr
	Stellar age	tstar	0.5, 1, 3, 5, 7 Gyr
Simple stellar population; Bruzual & Charlot (2003)	Initial mass function	−	Chabrier (2003)
	Metallicity	Z	0.02
Dust attenuation; Calzetti et al. (2000)	Color excess	E(B − V)	0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.7, 0.9 mag
Galactic dust emission; Dale et al. (2014)	Slope in dMdust ∝ U−α dU	α	2
AGN (UV-to-IR) SKIRTOR	AGN contribution to IR luminosity	fracAGN	0–0.99 (step 0.1)
	Viewing angle	θ	60, 70, 80, 90°
	Polar-dust color excess	E(B − V)PD	0, 0.2, 0.4
X-ray	AGN photon index	Γ	1.8
	Maximum deviation from the αox–Lν,2500Å relation	$| {\rm{\Delta }}{\alpha }_{\mathrm{ox}}{| }_{\max }$	0.2
	AGN X-ray angle coefficients	(a1, a2)	(0, 0) / (0.5, 0) / (1, 0) / (0.33, 0.67)a

Note. For parameters not listed here, we use the default values. Bold font indicates new parameters in cigale v2022.0 introduced in this work. (a) Each set of angle coefficients is for one x-cigale run. (a₁, a₂) = (0, 0) indicates the x-cigale (isotropic) run.

Download table as:  ASCIITypeset image

Figure 1 (top left) shows an example SED fit for one of the COSMOS type 2 AGNs. Under the scheme of AGN unification (e.g., Antonucci 1993; Urry & Padovani 1995; Netzer 2015), our type 2 AGNs should also follow the α_ox–L_ν,2500Å intrinsic relation (Equation (1)). To test this point, we plot the Δα_ox = α_{ox, fitted} − α_{ox, expected} (i.e., deviation from the α_ox–L_ν,2500Å relation) distribution from our x-cigale run in Figure 2. The Δα_ox values tend to be systematically negative, with a median value of −0.093 (corresponding to a factor of 1.75 lower in terms of L_ν,2keV /L_ν,2500Å) for x-cigale. This result suggests that, with the assumption of isotropic AGN X-ray emission in x-cigale, the observed X-ray fluxes of our type 2 AGNs tend to systematically lie below the expectations from the α_ox–L_ν,2500Å relation. One natural solution to this issue is allowing intrinsic X-ray anisotropy, so that an AGN viewed at type 2 angles will have lower X-ray fluxes than viewed at type 1 angles. We perform this code-implementation task in Section 2.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Considering the evidence for X-ray anisotropy in Sections 2.1 and 2.2, we modify the code so that the user can model L_X as a second-order polynomial function of the cosine of the viewing angle (e.g., Netzer 1987):

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{{\rm{X}}}(\theta )}{{L}_{{\rm{X}}}(0)}={a}_{1}\cos \theta +{a}_{2}{\cos }^{2}\theta +1-{a}_{1}-{a}_{2},\end{eqnarray} \tag{ 2 }$$

where the coefficients (a₁ and a₂) are free parameters set by the user, and θ is the viewing angle (face-on = 0, edge-on = 90°). The constant term in Equation (2), 1 − a₁ − a₂, guarantees that the right-hand side equals the left-hand side when θ = 0. (a₁, a₂) = (0, 0) means isotropic L_X. In cigale v2022.0, α_ox is still calculated using Equation (1).

We repeat the fitting of type 2 AGNs (Section 2.1) but using cigale v2022.0. We perform three runs each with (a₁, a₂) set to (0.5, 0), (1, 0), and (0.33, 0.67), respectively, while the other parameters are the same as the run in Section 2.2 (see Table 1). (a₁, a₂) = (0.33, 0.67) means the same angular dependence as that of the AGN disk emission (L_disk) in the UV-to-IR AGN module (the SKIRTOR model; Stalevski et al. 2012, 2016). (a₁, a₂) = (1, 0) is equivalent to a thin disk geometry with an angle-independent X-ray intensity. The angle dependence of (a₁, a₂) = (0.5, 0) is between (a₁, a₂) = (1, 0) and the isotropic case. Figure 3 displays the viewing angular dependence under these (a₁, a₂) settings. The angular dependence is stronger on the order of (0.5, 0), (1, 0), and (0.33, 0.67).

Zoom In Zoom Out Reset image size

The SED fit for an example source is displayed in Figure 1. For this source, the anisotropic models have better fitting quality than the isotropic model (as indicated by the reduced χ² labeled in Figure 1). Figure 2 displays the Δα_ox distributions from the cigale v2022.0 runs. The settings of (a₁, a₂) = (1, 0) and (a₁, a₂) = (0.33, 0.67) lead to systematically positive Δα_ox, with median values of 0.069 and 0.120, respectively. This result indicates that the angular dependence defined by these two parameter sets is overly strong. In contrast, the Δα_ox from (a₁, a₂) = (0.5, 0) is more evenly distributed around zero compared to those from (a₁, a₂) = (1, 0), (0.33, 0.67), and the x-cigale (isotropic) result. The Δα_ox median (−0.027) is the smallest among all four models, indicating that (a₁, a₂) = (0.5, 0) is likely the most physical model among the tested ones.

To assess the overall fitting quality, we calculate the difference of the Akaike information criterion (ΔAIC; Akaike 1974) between the cigale v2022.0 (anisotropic) and the x-cigale (isotropic) fits. This quantity is defined as ΔAIC = 2Δk + Δχ², where Δk is the difference in the number of free parameters. Δk is zero for our case here. A lower ΔAIC indicates a stronger probability of anisotropic models. For example, ΔAIC < −4 means the anisotropic model is more than ≈7 (e^−ΔAIC/2) times more probable than the isotropic model, indicating a strong support for the former (e.g., Burnham & Anderson 2002).

The example source in Figure 1 has ΔAIC < −4, indicating that the anisotropic models are preferred over the isotropic model. When inspecting the residuals in Figure 1, one could be puzzled that the main difference between the x-cigale versus cigale v2022.0 fits is in the IR rather than X-ray. We note that the root of this difference is not related to the IR AGN emission model, as all fits are based on the same IR AGN models (Table 1). Instead, the actual cause is X-ray angle dependence, which is the only different setting among the fits. We briefly explain this cause below.

For our COSMOS type 2 sample, the AGNs have emission mostly in X-ray and IR as their UV/optical radiation is obscured. The X-ray/IR ratio is an observable quantity closely related to the X-ray angle dependence (e.g., Asmus et al. 2015). To see this point, we can write the X-ray/IR ratio as

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{{\rm{X}}}(\theta )}{{L}_{\mathrm{IR}}(\theta )}=\displaystyle \frac{{L}_{{\rm{X}}}(\theta )}{{L}_{{\rm{X}}}(30^\circ )}\times \displaystyle \frac{{L}_{{\rm{X}}}(30^\circ )}{{L}_{\mathrm{UV}/{\rm{O}}}(30^\circ )}\times \displaystyle \frac{{L}_{\mathrm{UV}/{\rm{O}}}(30^\circ )}{{L}_{\mathrm{IR}}(\theta )},\end{eqnarray} \tag{ 3 }$$

where L_X, L_IR, and L_UV/O are AGN X-ray, IR, and intrinsic UV/optical luminosities, respectively. In Equation (3), the second factor $\tfrac{{L}_{{\rm{X}}}(30^\circ )}{{L}_{\mathrm{UV}/{\rm{O}}}(30^\circ )}$ is roughly a constant value, as cigale directly links AGN X-ray and UV/optical emission by the α_ox–L_ν,2500Å relation at a 30° viewing angle (Yang et al. 2020). The third factor $\tfrac{{L}_{\mathrm{UV}/{\rm{O}}}(30^\circ )}{{L}_{\mathrm{IR}}(\theta )}$ is also about a constant (depending on the dust-model details), because the IR emission originates from the UV/optical photons absorbed by dust, and cigale strictly keeps energy conservation (Boquien et al. 2019; Yang et al. 2020). Therefore, the X-ray/IR ratio is approximately proportional to the first factor $\tfrac{{L}_{{\rm{X}}}(\theta )}{{L}_{{\rm{X}}}(30^\circ )}$, which is the X-ray angle dependence. For type 2 viewing angles, compared to the X-ray isotropic model, our tested anisotropic models have lower $\tfrac{{L}_{{\rm{X}}}(\theta )}{{L}_{{\rm{X}}}(30^\circ )}$ values (Figure 3) and thereby systematically lower X-ray/IR ratios. For the source in Figure 1, the observed IR/X-ray ratio is more similar to the anisotropic model values than the isotropic one. This is the reason why the anisotropic configurations model the observed data better than the isotropic one.

This source in Figure 1 is also a representative example demonstrating that different bands are not modeled independently. cigale templates are rigid across all bands, finding a solution that minimizes the "global" χ², although such a solution might not minimize residuals in some bands.

Figure 4 shows the distribution of ΔAIC. For the model of (a₁, a₂) = (0.5, 0), 13% of the sources have ΔAIC < −4, while only 2% sources have ΔAIC > 4. The overall distribution is toward the negative sign (median ΔAIC = −0.11). This result indicates that the fitting quality has an overall improvement from the isotropic model to the (a₁, a₂) = (0.5, 0) anisotropic model. However, from Figure 4, the models of (a₁, a₂) = (1, 0) and (0.33, 0.67) have similar or even worse fitting quality than that of the isotropic model.

Zoom In Zoom Out Reset image size

From the Δα_ox and the AIC analyses above, an isotropic AGN X-ray model is disfavored compared to the anisotropic model with (a₁, a₂) = (0.5, 0). Therefore, AGN X-ray emission is likely weaker toward larger viewing angles, qualitatively consistent with the observations of Liu et al. (2014) and Asmus et al. (2015; see Section 2.1). On the other hand, the amplitude of this viewing-angle dependence is moderate, since (a₁, a₂) = (0.5, 0) results in better fitting quality than (a₁, a₂) = (1, 0) and (0.33, 0.67), which have stronger angular dependence (see Figure 3). The conclusion that AGN X-rays have weaker angular dependence than UV/optical [(a₁, a₂) = (0.33, 0.67)] is understandable. The X-ray photons result from the inverse Compton scattering of the UV/optical seed photons, and the strength of anisotropy is suppressed by this scattering process (e.g., Xu 2015; Yang et al. 2020).

We caution that our conclusion of X-ray angle dependence is for the overall type 2 AGN population rather than individual sources, as our analyses above are based on the statistical analyses of the entire type 2 sample. It might be possible that individual AGNs have different angular dependence, because the structure of the AGN corona, which produces the X-ray photons (Section 2.1), could vary among individual sources (e.g., Ricci et al. 2018; Tortosa et al. 2018).

We set Γ = 1.8 in our runs (Table 1), but the actual power-law photon index for an X-ray AGN may range from ≈1.6 to ≈2.2 (e.g., Yang et al. 2016; Liu et al. 2017). The photon-index parameter can affect the model-predicted X-ray fluxes. To assess this effect, we repeat our runs allowing Γ to vary between 1.6 and 2.0. The resulting Δα_ox and ΔAIC distributions are similar to those in Figures 2 and 4, and the (a₁, a₂) = (0.5, 0) configuration is still the most favored model. In Table 1, we adopt large viewing angles (≥60°) assuming the classic unification model, i.e., type 2 AGNs are obscured by the torus. However, some recent observations suggest that type 2 AGNs might also have small viewing angles and be obscured by polar dust (e.g., Mountrichas et al. 2021a; Ramos Padilla et al. 2022). To consider this possibility, we test new cigale runs allowing all available viewing-angle values in SKIRTOR (0–90° with a step of 10°). The result still favors the (a₁, a₂) = (0.5, 0) anisotropic model, consistent with our original result. Based on the tests above, we consider our main conclusion not to be critically dependent on the adopted parameters of the photon index and viewing angle in Table 1.

In the code of cigale v2022.0, we set the default (a₁, a₂) to (0.5, 0) based on our results above. For general purposes of AGN modeling, the user does not need to change these default values. For the specific purposes of studying AGN X-ray anisotropy, the user can test different (a₁, a₂) values in different runs and select the best parameters, like in our approach above. This method allows further studies of X-ray angular dependence for different AGN samples (e.g., high-accretion rates versus low-accretion rates), and thereby can provide insight into the properties of AGN coronae.

Both AGNs and normal galaxies can emit X-ray photons. Normal-galaxy X-rays originate primarily from point sources of X-ray binaries and diffuse hot gas. AGNs tend to be more luminous than normal galaxies at X-ray wavelengths. As a consequence, most of the X-ray detected sources in extragalactic surveys are AGNs. However, normal galaxies become increasingly important as survey depth improves. The 7 Ms Chandra Deep Field-South (CDF-S) survey has ≈30% of the X-ray detections classified as normal galaxies, and such sources dominate the faintest detections (Luo et al. 2017). It is thereby expected that many more normal galaxies will be detected in deep surveys by future X-ray telescopes with large collecting areas such as Athena and Lynx. Therefore, it is critical to have realistic recipes for normal-galaxy X-ray modeling.

x-cigale has both AGN and galaxy X-ray components (Yang et al. 2020). The latter includes the emission from high-mass X-ray binaries (HMXBs), low-mass X-ray binaries (LMXBs), and hot gas. The AGN component has been well tested (e.g., Yang et al. 2020; Zou et al. 2020; Mountrichas et al. 2021b), but this is not the case for the galactic component. Below, we test and improve the modeling of galaxy X-ray emission.

We test the galaxy X-ray modeling in x-cigale using the 7 Ms CDF-S survey, which is the deepest X-ray survey to date (Luo et al. 2017). We take advantage of this unique data set to study the X-ray emission from normal galaxies in the distant universe, as galaxies' X-ray power is typically low (L_X ≲ 10^42.5 erg s⁻¹) and beyond the sensitivity of most X-ray surveys.

We first select sources classified as "galaxy" instead of "AGN" or "star" by Luo et al. (2017). The classification is based on X-ray and other multiwavelength data. We then restrict the sample to only sources within the GOODS-S field (Guo et al. 2013), where deep multiwavelength coverage is available from UV to far-IR (FIR). We compile the UV-to-IRAC4 data from Guo et al. (2013) and Spitzer/Herschel mid-IR to FIR data from the ASTRODEEP team (T. Wang 2020, private communication; GOODS-S Herschel catalog). We discard sources with MIPS 24 μm S/N < 3, as reliable IR data are essential in constraining star formation rate (SFR; which scales with ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$) and possible low-level AGN activity. There are a total of 39 X-ray detected galaxies in our sample. We adopt the redshift measurements compiled by Luo et al. (2017), which are secure spectroscopic redshifts or high-quality photometric redshifts. The redshifts cover a range of z = 0.10–1.06 (10th–90th percentiles), with a median of z = 0.44.

We first perform SED modeling of these galaxies using x-cigale. The fitting parameters are summarized in Table 2. The galaxy settings are similar to those in Section 2, except that we allow two metallicity values of Z = 0.004 and 0.02, as the ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$–SFR scaling relation depends on metallicity (e.g., Fragos et al. 2013a). We still allow a moderate AGN component in the fitting (frac_AGN ≤ 0.2). Although the sources are classified as galaxies by Luo et al. (2017), some could possibly be low-luminosity AGNs (e.g., Young et al. 2012; Ding et al. 2018).

Table 2. Model Parameters for the CDF-S Normal Galaxies

Module	Parameter	Symbol	Values
Star formation history $\mathrm{SFR}\propto t\exp (-t/\tau )$	Stellar e-folding time	τstar	0.1, 0.5, 1, 5 Gyr
	Stellar age	tstar	0.5, 1, 3, 5, 7 Gyr
Simple stellar population; Bruzual & Charlot (2003)	Initial mass function	−	Chabrier (2003)
	Metallicity	Z	0.004, 0.02
Dust attenuation; Calzetti et al. (2000)	Color excess of the nebular lines	E(B − V)	0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 mag
Galactic dust emission; Dale et al. (2014)	Slope in dMdust ∝ U−α dU	α	2
AGN (UV-to-IR) SKIRTOR	AGN contribution to IR luminosity	fracAGN	0, 0.01, 0.03, 0.1, 0.2
	Viewing angle	θ	30°, 70°
	Polar-dust color excess	E(B − V)PD	0
X-ray	Deviation from the expected $\mathrm{log}{L}_{{\rm{X}}}^{\mathrm{HMXB}}$	δHMXB	− 0.5 to 0.5 (step 0.1) dex
	Deviation from the expected $\mathrm{log}{L}_{{\rm{X}}}^{\mathrm{LMXB}}$	δLMXB	− 0.5 to 0.5 (step 0.1) dex

Note. For parameters not listed here, we use the default values. Bold font indicates new parameters in cigale v2022.0 introduced in this work.

Download table as:  ASCIITypeset image

We compare the model versus observed X-ray fluxes in Figure 5 for x-cigale (left panel). For many (21%) of our sources, the offsets between the model and observed fluxes are more than 0.5 dex. We show an example SED fit with such an issue in Figure 6 (left panel). Therefore, x-cigale is not able to model well all of the observed X-ray fluxes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

x-cigale assumes that galaxy X-ray emission from HMXBs and LMXBs can be calculated from the scaling relations of ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$–SFR and ${L}_{{\rm{X}}}^{\mathrm{LMXB}}$–M_⋆ (Fragos et al. 2013a). However, this is an oversimplified assumption, because these relations are just an approximation for the overall galactic population, and scatters around them exist. For example, the content of globular clusters at a given M_⋆, which is not modeled in x-cigale, can significantly affect ${L}_{{\rm{X}}}^{\mathrm{LMXB}}$ (e.g., Lehmer et al. 2020). Also, since the HMXB and LMXB emissions are from discrete point sources, ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$ and ${L}_{{\rm{X}}}^{\mathrm{LMXB}}$ inevitably suffer from statistical fluctuations that are especially strong in low-SFR and/or low-M_⋆ galaxies (e.g., Lehmer et al. 2019, 2021).

To model the ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$ and ${L}_{{\rm{X}}}^{\mathrm{LMXB}}$ dispersions of individual galaxies in greater detail, we introduce two new free parameters, δ_HMXB and δ_LMXB, to account for the scatters of the ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$–SFR and ${L}_{{\rm{X}}}^{\mathrm{LMXB}}$–M_⋆ scaling relations, i.e.,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\left(\displaystyle \frac{{L}_{2-10\ \mathrm{keV}}^{\mathrm{HMXB}}}{\mathrm{SFR}}\right) & = & 40.3-62Z+569{Z}^{2}-1834{Z}^{3}\\ & & +\,1968{Z}^{4}+{\delta }_{\mathrm{HMXB}}\\ \mathrm{log}\left(\displaystyle \frac{{L}_{2-10\ \mathrm{keV}}^{\mathrm{LMXB}}}{{M}_{\star }}\right) & = & 40.3-1.5\mathrm{log}t-0.42{(\mathrm{log}t)}^{2}\\ & & +\,0.43{(\mathrm{log}t)}^{3}+0.14{(\mathrm{log}t)}^{4}+{\delta }_{\mathrm{LMXB}},\end{array}\end{eqnarray} \tag{ 4 }$$

where M_⋆ and SFR are in solar units; t denotes stellar age in gigayears; Z denotes metallicity (mass fraction). The parameters of δ_HMXB and δ_LMXB are logarithmic deviations from the scaling relations, with positive/negative values meaning higher/lower HMXB and LMXB luminosities, respectively. The user can set multiple values for each parameter to enable a more flexible XRB prescription.

In addition to introducing δ_HMXB and δ_LMXB, we also implement another update of the code. The code provides three SFR parameters: the instantaneous SFR, the average SFR over 10 Myr, and the average SFR over 100 Myr. While x-cigale adopted the instantaneous SFR when calculating ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$, we adopt the average SFR over 100 Myr in cigale v2022.0. This change is because the HMXB emission has strong variability on ∼10 Myr timescales (e.g., Linden et al. 2010; Garofali et al. 2018; Antoniou et al. 2019), but we do not have well-informed calibrations for how the ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$ varies on such short timescales. On the other hand, the ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$ dependence on longer timescales of ∼100 Myr has been carefully characterized (e.g., Lehmer et al. 2019, 2021). Although the instantaneous SFR and the 100-Myr averaged SFR are similar for a smooth SFH, they can differ significantly if a recent burst/quenching is present in the SFH.

With cigale v2022.0, we re-fit our sample of CDF-S normal galaxies. We set δ_HMXB from −0.5 to 0.5 dex with a step of 0.1 dex, while keeping the other parameters unchanged (Table 2). This parameter range is chosen because Lehmer et al. (2021) found the 2σ scatter of ${L}_{{\rm{X}}}^{\mathrm{HMXB}}/$ SFR is ≈0.5 dex at SFR ≈4 M_⊙ yr⁻¹, which is the median SFR of our sample. We also set δ_LMXB to the same values as δ_HMXB.

Figure 5 (right) compares the cigale v2022.0 resulting model versus observed X-ray fluxes. The model fluxes agree much better with the observed fluxes compared to those fitted by x-cigale (Figure 5, left). The offsets between the model and observed are all within 0.5 dex. Figure 6 (right) shows an example SED fit with cigale v2022.0. For this example, compared to x-cigale, cigale v2022.0 has a better fit not only to the X-ray data but also to the UV data. This is because x-cigale is forced to use a stellar population model that corresponds to a relatively high X-ray emission, although this model does not well fit the observed UV fluxes. This example highlights the importance of introducing δ_HMXB and δ_LMXB, without which inappropriate stellar models might be selected.

Figure 7 displays the distributions of the fitted δ_HMXB and δ_LMXB, respectively. Both distributions have slightly positive median values, i.e., 0.04 dex (HMXB) and 0.12 dex (LMXB). These near-zero medians indicate that the ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$ and ${L}_{{\rm{X}}}^{\mathrm{LMXB}}$ scaling relations (Fragos et al. 2013a) are good approximations for the overall galactic population detected in deep X-ray surveys at z ≲ 1. The slightly positive trend of the distributions suggests that the scaling relations might have systematic offsets. But the positive trend is expected due to a selection effect, because our X-ray data are flux-limited and thus tend to select higher ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$ and ${L}_{{\rm{X}}}^{\mathrm{LMXB}}$ sources. A larger normal-galaxy X-ray sample, from, e.g., eROSITA (e.g., Vulic et al. 2021), is needed to investigate the nature of the positive trend of δ_HMXB and δ_LMXB. Both of the δ_HMXB and δ_LMXB distributions have substantial scatters (standard deviations ≈0.2 dex; Figure 7). These scatters are likely caused by, e.g., globular-cluster contents and statistical fluctuations (see Section 3.3).

Zoom In Zoom Out Reset image size

We set frac_AGN ≤ 0.2 in our runs (Table 2) since the sources were classified as normal galaxies by Luo et al. (2017), but this quantitative choice is rather arbitrary. One might worry that our fitting results could heavily depend on the assumption of frac_AGN ≤ 0.2, as AGNs can also contribute to the observed X-ray fluxes. To assess the effects of this assumption, we perform a new run including higher frac_AGN values, i.e., frac_AGN = 0, 0.01, 0.03, 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6, while keeping the other parameters the same. The resulting δ_HMXB and δ_LMXB distributions are identical to those in Figure 7, indicating that these two parameters are not strongly degenerate with frac_AGN. Figure 8 displays example 2D and 1D probability distributions of frac_AGN and δ_LMXB from the new fit with frac_AGN = 0–0.6. The frac_AGN is tightly constrained at a low level (≲0.01), and the situations are similar for all of our sources. The tight AGN constraints are understandable as the AGN emission is constrained not only by the X-ray data but also by the UV-to-IR data. In summary, we conclude that our fitting results based on the parameters in Table 2 do not depend on the assumption of frac_AGN ≤ 0.2.

Zoom In Zoom Out Reset image size

We set the default values of δ_HMXB and δ_LMXB both to 0, corresponding to the standard Fragos et al. (2013b) scaling relations. For luminous X-ray sources (e.g., L_X ≳ 10^42.5 erg s⁻¹), the observed X-ray fluxes are likely dominated by AGNs, and thus, the user can just keep the default values of δ_HMXB = δ_LMXB = 0. For less luminous sources, galactic X-ray emission could dominate the observed fluxes, and the user can adopt different values of δ_HMXB and δ_LMXB (e.g., Table 2) to allow for more flexible XRB modeling. For specific galactic populations, the user could allow multiple values for only one of δ_HMXB and δ_LMXB to save memory and reduce computation time. For example, for quiescent galaxies, ${L}_{{\rm{X}}}^{\mathrm{HMXB}}$ should be negligible compared to ${L}_{{\rm{X}}}^{\mathrm{LMXB}}$. In this case, the user can adopt multiple values for δ_LMXB while keeping δ_HMXB = 0.

x-cigale adopts a single fixed SED shape of an AGN accretion disk from Schartmann et al. (2005), which is a broken power law. Although the Schartmann et al. (2005) recipe is a good approximation for the overall disk SED shape, it might not be sufficiently accurate for individual sources. This is because the observed UV/optical slopes of type 1 quasars have non-negligible intrinsic dispersions (e.g., Elvis et al. 1994), possibly due to different black hole (BH) masses, accretion rates, and spins (e.g., Koratkar & Blaes 1999).

To test whether this single spectral shape is sufficient to account for the observed SEDs, we use the SDSS DR14 type 1 quasar sample (Pâris et al. 2018) that has XMM-Newton X-ray detections (see Section 3.1.1 of Yang et al. 2020 for details). We further apply a magnitude cut (r_AB < 21.8 mag) and a redshift cut (z > 1) to ensure that the observed SEDs are dominated by AGNs, as SDSS normal galaxies with r_AB < 21.8 mag are always below z = 1 (e.g., Sheldon et al. 2012). Therefore, these cuts allow us to model SEDs with pure-AGN templates as below, avoiding potential degeneracy issues (see Section 4.4).

The final sample has 1080 sources. We run x-cigale on the SDSS ugriz and the 2–10 keV fluxes. The inclusion of X-ray photometry is to better constrain the AGN intrinsic emission. The fitting parameters are listed in Table 3. We set frac_AGN (fractional AGN IR luminosity) to a value close to unity (0.9999), so that the observed UV/optical SED is totally AGN dominated, which is the case for the SDSS quasars after our magnitude and redshift cuts. We also allow different levels of polar-dust extinction.

Table 3. Model Parameters for the SDSS Quasars

Module	Parameter	Symbol	Values
AGN (UV-to-IR) SKIRTOR	AGN contribution to IR luminosity	fracAGN	0.9999
	Viewing angle	θ	30°
	Polar-dust color excess	E(B − V)PD	0., 0.01, 0.02, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6 mag
	Intrinsic-disk type	—	Schartmann et al. (2005)
	Deviation from the default UV/optical slope	δAGN	− 1 to 1 (step 0.1)
X-ray	AGN photon index	Γ	1.8
	Maximum deviation from the αox–Lν,2500Å relation	$| {\rm{\Delta }}{\alpha }_{\mathrm{ox}}{| }_{\max }$	0.2

Note. For parameters not listed here, we use the default values. Bold font indicates new parameters in cigale v2022.0 introduced in this work.

Download table as:  ASCIITypeset image

We compare the resulting x-cigale model versus observed u − z colors in Figure 9 (left). In this figure, we only consider sources having both u and z S/Ns above 5. One notable issue is that a "plateau" exists at model u − z ≈ 0.5. This is because the model u − z cannot be bluer than the intrinsic-disk color (u − z = 0.5), but a significant fraction (32%) of sources have observed u − z < 0.5. Due to this issue, the offsets between the model and observed u − z have a positive median value of 0.15. Figure 10 (left) shows an example SED. The observed SED is even bluer than the zero-extinction model SED. To address this issue, it is necessary to allow a flexible SED shape for AGN intrinsic-disk emission.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To allow for deviations from the default Schartmann et al. (2005) optical spectral slope as suggested by the observed quasars (Section 4.2), we introduce a free parameter, δ_AGN, i.e.,

$$\begin{eqnarray}\lambda {L}_{\lambda }\propto \left\{\begin{array}{lr}{\lambda }^{2} & 0.008\leqslant \lambda \leqslant 0.05\ [\mu {\rm{m}}]\\ {\lambda }^{0.8} & 0.05\lt \lambda \leqslant 0.125\ [\mu {\rm{m}}]\\ {\lambda }^{-0.5+{\delta }_{\mathrm{AGN}}} & 0.125\lt \lambda \leqslant 10\ [\mu {\rm{m}}]\\ {\lambda }^{-3} & 10\lt \lambda \leqslant 1000\ [\mu {\rm{m}}]\end{array}\right..\end{eqnarray} \tag{ 5 }$$

In cigale v2022.0, we also allow the user to choose the disk continuum from the SKIRTOR model (Stalevski et al. 2012, 2016), i.e.,

$$\begin{eqnarray}\lambda {L}_{\lambda }\propto \left\{\begin{array}{lr}{\lambda }^{1.2} & 0.008\leqslant \lambda \leqslant 0.01\ [\mu {\rm{m}}]\\ {\lambda }^{0} & 0.01\lt \lambda \leqslant 0.1\ [\mu {\rm{m}}]\\ {\lambda }^{-0.5+{\delta }_{\mathrm{AGN}}} & 0.1\lt \lambda \leqslant 5\ [\mu {\rm{m}}]\\ {\lambda }^{-3} & 5\lt \lambda \leqslant 1000\ [\mu {\rm{m}}]\end{array}\right..\end{eqnarray} \tag{ 6 }$$

From Equations (5) and (6), the major differences between the Schartmann et al. (2005) and SKIRTOR disk continuum are the wavelength boundaries and power-law indices at far-UV (λ < 125 nm), where observational constraints are weaker compared to those at longer wavelengths (e.g., Stevans et al. 2014; Lusso et al. 2015).

We re-fit the SDSS quasar SEDs using cigale v2022.0. We still adopt the Schartmann et al. (2005) disk SED shape (Equation (5)), since it agrees better with observations than the SKIRTOR disk model (Equation (6); e.g., Duras et al. 2017). We allow δ_AGN to vary from −1 to 1 with a step of 0.1 (see Table 3). This δ_AGN range covers nearly all of the observed SED slope variations in the SDSS quasar sample of Davis et al. (2007; see their Figure 2).

We compare the model versus observed u − z colors in Figure 9 (right). The model colors agree better with the observed colors than those from the x-cigale fits. The median offset of the model and observed u − z is close to zero (0.05), while the x-cigale fits have a median of 0.15. Also, the scatter has been reduced from 0.24 (x-cigale fits) to 0.16 (cigale v2022.0 fits). Figure 10 (right) shows an example SED fit. The observed SED can be well fitted with a reduced ${\chi }_{r}^{2}=0.6$ (the value is 3.4 for the x-cigale fit). Figure 11 displays the 2D and 1D probability distributions of δ_AGN and E(B − V)_PD of this example fit. There is an anticorrelation between δ_AGN and E(B − V)_PD in the 2D probability distribution. This anticorrelation is expected, because, e.g., a higher δ_AGN (redder intrinsic slope) and a lower E(B − V)_PD (weaker polar-dust extinction) can roughly cancel out the effects of one another. Therefore, there is a natural degeneracy between these two parameters.

Zoom In Zoom Out Reset image size

Figure 12 displays the distribution of the fitted δ_AGN and E(B − V)_PD. The median of the δ_AGN distribution is −0.36, with a significant scatter of 0.34. This negative median value could be intrinsic, as the default Schartmann et al. (2005) spectral shape is based on the observed quasar SEDs, without considering polar-dust extinction. Indeed, the polar-dust extinctions are non-negligible [median E(B − V)_PD = 0.03], although heavy extinctions of E(B − V)_PD > 0.1 are rare (6%). However, we note that the degeneracy between δ_AGN and E(B − V)_PD (e.g., Figure 11) might also contribute to the negative trend of δ_AGN. We caution that both of the δ_AGN and E(B − V)_PD distributions quantitatively depend on the adopted specific extinction law. We adopt the default SMC law here, and refer to Buat et al. (2021) for a detailed discussion on the effects of different laws.

Zoom In Zoom Out Reset image size

Given the degeneracy between δ_AGN and E(B − V)_PD, one might think of adapting the original disk SED shapes (Schartmann and SKIRTOR; Section 4.3) to the bluest variation in our sample and attributing all observed SED variations to E(B − V)_PD. This idea has the benefit of simplicity, but we do not adopt it for three reasons. First, it is unlikely that all AGNs share the same intrinsic SED shape over wide ranges of BH masses and accretion rates (e.g., Whiting et al. 2001; Richards et al. 2003). Second, the current approach is more flexible than adapting the template, and the philosophy of cigale highlights flexibility rather than simplicity. Third, although E(B − V)_PD and δ_AGN are degenerate over UV/optical wavelengths, they can be better differentiated given excellent IR coverage, because dust-reddened AGNs have polar-dust IR re-emission (included in cigale), but intrinsically red AGNs do not. In cigale, polar dust is an obscuration structure with an optical depth much smaller than the torus (see Section 2.4 of Yang et al. 2020 for details). cigale assumes that the polar-dust IR emission follows a "graybody" model (e.g., Casey 2012) with temperature and emissivity as free parameters. Future IR missions, e.g., the James Webb Space Telescope (JWST) and Origins, will be able to detect (tightly constrain) the polar-dust IR re-emission and thereby differentiate the two cases of dust-reddened versus intrinsically red SEDs.

We set the default δ_AGN and E(B − V)_PD to the median values of our fits, i.e., δ_AGN = −0.36 and E(B − V)_PD = 0.03. When fitting type 1 AGNs, the user is recommended to adopt multiple values of these two parameters (such as those in Table 3), because both parameters have significant scatters based on our fits (see Figure 12). When fitting type 2 AGNs, for which the AGN disk emission is almost entirely obscured by the dusty torus, the user can keep the default δ_AGN to save memory and reduce computational time. However, it is still recommended to adopt multiple values of E(B − V)_PD when mid-IR to FIR coverage is available for the type 2 sources, because E(B − V)_PD sets the strength of polar-dust re-emission that could contribute significantly at IR wavelengths.

Finally, we stress the importance of our redshift and magnitude cuts (Section 4.2). These cuts guarantee the observed SEDs are dominated by AGNs, allowing us to model the data with pure-AGN templates. Actually, we have tested fits with AGN+galaxy mixed templates by freeing frac_AGN, and found our parameters of interest (δ_AGN and E(B − V)_PD) could be significantly affected. This is because a blue observed SED can be either explained by a low-dust star-forming galactic component or an AGN component. Figure 13 displays such an AGN–galaxy degeneracy for an SDSS quasar, for which δ_AGN and E(B − V)_PD are strongly affected when a galactic component is allowed. Therefore, our redshift/magnitude cuts and pure-AGN approach are crucial for our investigation of δ_AGN and E(B − V)_PD. The user who is interested in these parameters should be cautious of the degeneracy effects when modeling SEDs that have non-negligible galactic components. On the other hand, some studies indicate that some other source properties such as frac_AGN, AGN bolometric luminosity (L_AGN), and SFR are not strongly affected by the degeneracy issue, especially when good multiwavelength coverage is available (e.g., Mountrichas et al. 2021c; Yang et al. 2021; Thorne et al. 2022). This is understandable, considering that those properties can be constrained by multiwavelength data simultaneously, e.g., L_AGN is related to X-ray, UV/optical, and IR wavelengths. However, in contrast, δ_AGN and E(B − V)_PD are very sensitive to the detailed SED-shape modeling at UV/optical wavelengths, and thus they are more strongly affected by the AGN–galaxy degeneracy than properties like L_AGN.

Zoom In Zoom Out Reset image size

x-cigale can only account for radio emission from SF (Boquien et al. 2019). This SF radio emission has two components: one is a thermal component contributed by the "nebular" module; the other is a synchrotron component contributed by the "radio" module. The latter is often dominant, and is calculated in x-cigale using the radio-IR correlation parameter q_IR (e.g., Helou et al. 1985), i.e.,

$$\begin{eqnarray}&&{q}_{\mathrm{IR}}=\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{SF},\ \mathrm{IR}}}{{L}_{\nu ,21\mathrm{cm}}\times 3.75\times {10}^{12}{Hz}}\right)\end{eqnarray} \tag{ 7 }$$

where L_{SF, IR} is the total star-forming IR luminosity (mostly in FIR) and L_ν,21cm is the corresponding radio synchrotron luminosity at 21 cm (1.4 GHz). The default value of q_IR is 2.58 in x-cigale. Other than q_IR, which sets the normalization at 21 cm, there is another free parameter (α_SF) that controls the power-law slope of the SF synchrotron emission, i.e.,

$$\begin{eqnarray}&&{L}_{\nu ,\mathrm{SF}}\propto {\nu }^{-{\alpha }_{\mathrm{SF}}}.\end{eqnarray} \tag{ 8 }$$

The default value of α_SF is 0.8 in both x-cigale and cigale v2022.0.

x-cigale does not have AGN emission at radio wavelengths. However, AGNs may have powerful jets that emit strong radio radiation (i.e., radio-loud AGNs), and the jets can play an important role in AGN–galaxy coevolution (e.g., Fabian 2012). The physical origin of AGN radio jets is still controversial. One popular theory is the Blandford-Znajek (BZ) process (e.g., Blandford & Znajek 1977; Blandford et al. 2019). The BZ mechanism considers that the jet is powered by the rotational energy of the BH through the magnetic field threading the horizon (e.g., Davis & Tchekhovskoy 2020). Recent observations suggest that the magnetic flux/topology close to the BH instead of the BH spin could be the determining factor of the jet-launching process (e.g., Zhu et al. 2020). Other than the jets, other processes such as AGN winds, coronae, and shocks can also emit at radio wavelengths (e.g., Panessa et al. 2019).

Similar to the procedures of the previous sections (Sections 2, 3, and 4), we first compile a proper radio-selected sample and then perform SED modeling with x-cigale in this section. We adopt all of the >5σ radio detections in the VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017a, 2017b). We also collect the VLA 1.4 GHz fluxes when available from Schinnerer et al. (2010). Delvecchio et al. (2017) matched the radio sources to the COSMOS2015 catalog (Laigle et al. 2016). We adopt these matching results and obtain the UV-to-IRAC4 broadband photometry (14 bands) from COSMOS2015. We discard the radio sources without COSMOS2015 counterparts, as the UV-to-IRAC4 data are necessary to model the stellar population. The sample contains 6497 sources in total. We also include Spitzer/MIPS (24 μm), Herschel/PACS (100 μm and 160 μm), and Herschel/SPIRE (250 μm, 350 μm, and 500 μm) photometry from the "super-deblended" catalog of Jin et al. (2018). We do not include X-ray fluxes here due to the reason presented in Section 5.4, i.e., we want to keep our SED fits and subsequent source classifications independent from the X-ray information. We adopt the redshift measurements from Delvecchio et al. (2017), which are spec-z (if available) or photo-z. The median redshift is 1.18, and the 10th–90th percentile range is z = 0.42–2.56.

We first fit the photometric data above with x-cigale. The model parameters are listed in Table 4. We set q_IR to a range of 2.4–2.7 (step 0.1), based on the observations of Delvecchio et al. (2021). Figures 14 (left) and 15 (left) display the resulting model fluxes versus the observed values for 3 GHz and 1.4 GHz, respectively. The model fluxes are systematically lower than the observed ones, e.g., 28% (11%) of sources have observed 3 GHz fluxes more than three (10) times higher than the model fluxes. In contrast, no sources have model fluxes greater than three times higher than the observed values. This result of "radio excess" strongly indicates that an AGN radio component is needed to explain the observed radio fluxes for many sources (Azadi et al. 2020). Figure 16 (left) shows an example SED fit with significant radio excess.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Model Parameters for the COSMOS Radio Sources

Module	Parameter	Symbol	Values
Star formation history $\mathrm{SFR}\propto t\exp (-t/\tau )$	Stellar e-folding time	τstar	0.1, 0.5, 1, 5 Gyr
	Stellar age	tstar	0.5, 1, 3, 5, 7 Gyr
Simple stellar population; Bruzual & Charlot (2003)	Initial mass function	−	Chabrier (2003)
	Metallicity	Z	0.02
Dust attenuation; Calzetti et al. (2000)	Color excess of the nebular lines	E(B − V)	0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 mag
Galactic dust emission; Dale et al. (2014)	Slope in dMdust ∝ U−α dU	α	2
AGN (UV-to-IR) SKIRTOR	AGN contribution to IR luminosity	fracAGN	0–0.99 (step 0.1)
	Viewing angle	θ	30°, 70°
	Polar-dust color excess	E(B − V)PD	0, 0.2, 0.4 mag
Radio	SF radio-IR correlation parameter	qIR	2.4, 2.5, 2.6, 2.7
	SF power-law slope	αSF	0.8
	Radio-loudness parameter	RAGN	0.01, 0.02, 0.05, 0.1, 0.2, 0.5,..., 1000, 2000, 5000, 10000
	AGN power-law slope	αAGN	0.7

Note. For parameters not listed here, we use the default values. Bold font indicates new parameters in cigale v2022.0 introduced in this work.

Download table as:  ASCIITypeset image

We add a new AGN component to the radio module of x-cigale. To quantitatively model AGN radio emission, we employ the radio-loudness parameter, R, defined as (e.g., Ballo et al. 2012),

$$\begin{eqnarray}&&{R}_{\mathrm{AGN}}=\displaystyle \frac{{L}_{\nu ,5\mathrm{GHz}}}{{L}_{\nu ,2500{\rm{\mathring{\rm A} }}}},\end{eqnarray} \tag{ 9 }$$

where L_ν,5GHz and L_ν,2500Å are the monochromatic AGN luminosities per frequency at rest frame 5 GHz and 2500 Å, respectively. R_AGN is a free parameter that allows for any values ≥0 (R_AGN = 0 means no AGN radio emission).

Here, we adopt L_ν,2500Å as the intrinsic (polar-dust absorption corrected) luminosity observed at a viewing angle of 30°, and this quantity is available for x-cigale models (see Section 2). This definition ensures that R is a physical quantity inherent to the AGN itself and does not depend on the viewing angle (e.g., Padovani 2016; Padovani et al. 2017). Therefore, R works consistently for both type 1 and type 2 AGNs. Currently, we assume L_ν,5GHz is isotropic in x-cigale. In the future, we will model the radio anisotropy, which can be important for, e.g., blazars and BL Lac objects.

We assume a power-law AGN SED over the wavelength range of 0.1–1000 mm, i.e.,

$$\begin{eqnarray}&&{L}_{\nu ,\mathrm{AGN}}\propto {\nu }^{-{\alpha }_{\mathrm{AGN}}},\end{eqnarray} \tag{ 10 }$$

where we allow the user to freely set the α_AGN slope. We set the default value as α_AGN = 0.7 (e.g., Randall et al. 2012; Tiwari 2019). We caution that the power-law shape is an overall simplistic assumption, as the real AGN radio SEDs might be more complicated. The formula in Equation (10) mainly serves as a correction for the AGN contribution to radio fluxes, especially for the cases where only one or two radio bands are available like our COSMOS radio sample. In the future, we will explore more realistic and complicated radio models based on multiband radio data.

Using cigale v2022.0, we re-fit the photometric data of the COSMOS radio sources (Section 5.1). We set R_AGN to a wide logarithmically spaced grid from 0.01 to 10000 (see Table 4) based on the observations of quasars (e.g., Zhu et al. 2020). We fix the radio slope in Equation (10) at the default value of α_AGN = 0.7, as most (73%) of our sources only have one radio band (3 GHz) available.

Unlike the results from x-cigale, the model radio fluxes agree well with the observed fluxes (see Figures 14 and 15). The offsets between the model and observed radio fluxes are mostly (99.91% for 3 GHz and 98.6% for 1.4 GHz) within 0.5 dex. Therefore, we conclude that our implementation of the AGN radio component is indeed useful in explaining the observed radio flux, although the good fit of radio data is not surprising given only one (or two) radio band(s) is available for each source in our sample. Interestingly, the model 1.4 GHz fluxes tend to be slightly lower than the observed ones (median offset = 0.08 dex), suggesting that the typical AGN radio SED in the COSMOS sample is steeper than our assumed α_AGN = 0.7 (Table 4). However, this systematic offset might also be a selection effect, as the relatively shallow 1.4 GHz data may miss AGNs with flatter radio SEDs and therefore lower 1.4 GHz fluxes.

Figure 16 compares example SED fits from x-cigale and cigale v2022.0. The observed radio fluxes are dominated by the AGN component from the cigale v2022.0 fit. The x-cigale fit is not able to explain the radio fluxes due to the lack of AGN radio emission. Compared to the cigale v2022.0 fit, the x-cigale fit has a stronger galactic IR component. This is because x-cigale only has galactic radio emission, which is related to galactic IR emission through the radio-IR correlation (Section 5.1). The high observed radio flux forcibly elevates not only galactic radio emission but also its IR emission as a consequence. In the cigale v2022.0 fit, the radio flux can be mostly explained by the AGN component, and thus the strong requirement of a galactic component is relaxed.

cigale v2022.0 can calculate AGN rest-frame 1.4 GHz luminosity (P_AGN,1.4GHz; e.g., Padovani 2016) as a measure of AGN radio strength. In this work, we consider the sources with P_AGN,1.4GHz/δ P_AGN,1.4GHz > 2 (where δ P_AGN,1.4GHz is the P_AGN,1.4GHz uncertainty from cigale v2022.0) as radio AGNs and the rest as radio SF galaxies. This definition guarantees that the AGN radio component is statistically significant (>2σ) for the classified radio AGNs. We note that our definition of AGN/SF is based on the radio-band decomposition, because our focus here is radio emission. For example, if a source has AGN features at other wavelengths (e.g., X-ray; see below) but its AGN radio emission is insignificant, it will be classified as a radio SF galaxy here. There are a total of 3221 radio AGNs, 50% of the sample. This high fraction indicates that radio AGNs are common among the sources selected by deep radio surveys. The radio-AGN fraction depends on radio fluxes. The fractions are 47% and 85% for sources with 3 GHz fluxes below and above 0.2 mJy, respectively. This significant radio-flux dependence is also found by Smolčić et al. (2017b), who used empirical criteria to classify AGNs and SF galaxies.

Figure 17 displays R_AGN versus P_AGN,1.4GHz and their distributions for our radio AGNs (P_AGN,1.4GHz/δ P_AGN,1.4GHz > 2). The red dashed line marks the conventional threshold (i.e., R_AGN =10; Kellermann et al. 1989) for radio-quiet (RQ) versus radio-loud (RL) AGN classifications. The numbers of RQ and RL AGNs are 823 (26%) and 2398 (74%), respectively. We remind the reader that this RQ/RL classification demonstrates an advantage of cigale v2022.0, which simultaneously models multiwavelength data in a consistent way (Section 5.3). Such a task is challenging for empirical approaches, because the AGN UV/optical emission is often heavily obscured and not directly observable (e.g., Figure 16).

Zoom In Zoom Out Reset image size

X-ray emission is a good tracer of the BH-accretion process (e.g., Brandt & Alexander 2015; Brandt & Yang 2021). It is intriguing to investigate the X-ray emission of our classified radio types. We adopt the Chandra COSMOS-Legacy survey (Civano et al. 2016; Marchesi et al. 2016). A total of 801 of the radio-selected sources in our radio sample are detected in X-ray. Figure 18 displays the fractions of X-ray detected sources among different radio types. The error bars represent binomial uncertainties calculated using astropy.stats.binom_conf_interval. The uncertainties are negligible compared to the differences across different radio types, thanks to our relatively large sample sizes.

Zoom In Zoom Out Reset image size

The X-ray fraction of the radio AGNs is 1.9 times higher than that of the radio SF. This result indicates that there is a positive link between AGN radio and X-ray emission, broadly consistent with the findings in the literature (e.g., Merloni et al. 2003; Laor & Behar 2008). Among the radio AGNs, the RQ population has a higher X-ray detected fraction than the RL population (Figure 18). This is expected, because RQs should have higher L_ν,2500Å than RLs at a given P_AGN,1.4GHz (see Equation (9)), and L_ν,2500Å is strongly correlated with AGN L_X due to the α_ox–L_ν,2500Å relation (e.g., Steffen et al. 2006; Just et al. 2007). However, the X-ray fraction of the RL AGNs is still significantly higher than that of the radio SF population (13% versus 9%). Assuming that the radio emission in RL AGNs is mainly from jets (Section 5), the elevated X-ray fraction of the RL population suggests a connection between jets and X-ray emission. This connection suggests that AGN jets could actively produce X-rays (e.g., Harris & Krawczynski 2006), or that there is a positive link between jets and the X-ray emitting coronae (e.g., Zhu et al. 2020).

We note that, since we do not include the X-ray data in our x-cigale run (Section 5.2), the X-ray detection fraction is independent of our SF, RQ, and RL classifications. Therefore, the X-ray fraction dependence on the radio type should be intrinsic, not a bias due to our SED-fitting procedure.

We set the default α_AGN = 0.7 (see Section 5.3) and R_AGN = 10 (i.e., the boundary between RL and RQ AGNs). For α_AGN, when there are multifrequency radio data spanning a large wavelength range, the user can adopt multiple α_AGN values to better model the observed radio fluxes. When only one or two radio bands within a narrow wavelength range (like our case) are available, the user can just keep the default α_AGN to save memory and reduce computational time. For the parameter of R_AGN, we recommend that the user adopt multiple values based on our fits (e.g., Figure 17). The user can narrow the range of R_AGN in some cases. For example, if the sources have spatially extended radio structures (strong evidence for radio-loud AGNs), then R_AGN can be set to >10 values.