杨光 G Yang

M Boquien

V Buat

D Burgarella

L Ciesla

F Duras

M Stalevski

W N Brandt

C Papovich

X-CIGALE: Fitting AGN/galaxy SEDs from X-ray to infrared

Keywords: methods, data analysis -methods, observational -galaxies, nuclei -quasars, general -X-rays, general 1 INTRODUCTION

CIGALE is a powerful multiwavelength spectral energy distribution (SED) fitting code for extragalactic studies. However, the current version of CIGALE is not able to fit X-ray data, which often provide unique insights into AGN intrinsic power. We develop a new X-ray module for CIGALE, allowing it to fit SEDs from the X-ray to infrared (IR). We also improve the AGN fitting of CIGALE from UV-to-IR wavelengths. We implement a modern clumpy two-phase torus model, SKIRTOR. To account for moderately extincted type 1 AGNs, we implement polar-dust extinction. We publicly release the source code (named "X-CIGALE"). We test X-CIGALE with X-ray detected AGNs in SDSS, COSMOS, and AKARI-NEP. The fitting quality (as indicated by reduced χ 2 ) is good in general, indicating that X-CIGALE is capable of modelling the observed SED from X-ray to IR. We discuss constrainability and degeneracy of model parameters in the fitting of AKARI-NEP, for which excellent mid-IR photometric coverage is available. We also test fitting a sample of AKARI-NEP galaxies for which only X-ray upper limits are available from Chandra observations, and find that the upper limit can effectively constrain the AGN SED contribution for some systems. Finally, using X-CIGALE, we assess the ability of Athena to constrain the AGN activity in future extragalactic studies.

The dust is thus heated and reemits the energy as infrared radiation. Considering the tight link between AGN multiwavelength SEDs and these physical processes, it is feasible to infer source properties from modelling the observed photometric data. On the other hand, the observed SED is often complicated, involving factors such as host-galaxy contributions and dust extinction. Misinterpretation of the SED could lead to unrealistic physical properties. Therefore, it is critical to decipher the observed data appropriately with a powerful and reliable SED fitting code.

The Code Investigating GALaxy Emission (CIGALE) is a state-of-the-art python code for SED fitting of extragalactic sources (Boquien et al. 2019). It employs physical AGN and galaxy models, and allows flexible combination between them. The current version of CIGALE can simultaneously fit the observed SED from UV to far-IR (FIR) and extract source physical properties such as AGN luminosity and host stellar mass (M ). However, the current CIGALE is not able to model X-ray fluxes, which often provide a unique view of AGNs.

X-ray observations have many advantages in AGN studies (see Brandt & Alexander 2015 for a review). Strong X-ray emission is nearly a universal property of the AGN phenomenon. X-rays are generated from the immediate vicinity of the BH, directly revealing the intrinsic AGN power. Therefore, X-ray fluxes are widely used as a tracer of BH accretion rate (e.g. Yang et al. 2018aYang et al. , 2019)). Thanks to their great penetrating power, X-rays are only mildly affected by obscuration in general. Also, AGNs are much more efficient in generating X-rays than their host galaxies. Therefore, the observed X-ray fluxes are often dominated by AGNs and have negligible galaxy contribution. Considering these advantageous properties, X-ray observations are widely used to select AGNs, especially in the distant universe, (e.g. Luo et al. 2017;Chen et al. 2018). These selections are often more complete and reliable than the selections at other wavelengths such as optical and IR.

Besides the lack of X-ray fitting capability, CIGALE's current AGN model [START_REF] Fritz | [END_REF], which covers the UV to IR, also has some other disadvantages. The model assumes that the central engine is surrounded by a dusty torus (i.e. the AGN unified model; Antonucci 1993;Urry & Padovani 1995;[START_REF] Netzer | The Physics and Evolution of Active Galactic Nuclei Netzer H[END_REF]Zou et al. 2019). The torus absorbs a fraction of the UV and optical emission from the central engine and reemits the energy as IR photons. When viewing from the equatorial direction, the central engine is obscured and only reemitted IR radiation can be observed (type 2 AGN). When viewing from the polar direction, the central engine is directly visible (type 1 AGN).

One disadvantage of the AGN model is that it assumes the dusty torus is a smooth structure. However, such smooth models for the torus are disfavored on physical grounds (e.g. Tanimoto et al. 2019). To reach a scale height consistent with observations, the dust grains in a smooth torus would have random velocities ∼ 100 km s -1 , corresponding to a temperature of ∼ 10 6 K. This high temperature far exceeds the dust-sublimation temperature (∼ 10 3 K). Another disadvantage of the AGN model is that the disk emission is assumed to be absolutely unextincted for the case of type 1. However, recent observations indicate that a non-negligible amount of extinction exists for some type 1 AGNs (e.g. Bongiorno et al. 2012;Elvis et al. 2012;Lusso et al. 2012), which can be attributed to the dust existing along polar directions (e.g. Stalevski et al. 2017Stalevski et al. , 2019;;Lyu & Rieke 2018). The current CIGALE cannot model the SEDs of these type 1 AGNs.

In this paper, we further develop CIGALE and enable it to fit X-ray data. The new development allows CIGALE to model AGN SED from X-ray to IR simultaneously and extract source properties such as AGN intrinsic luminosity and host-galaxy stellar mass (M ). Besides developing the X-ray part, we also improve CIGALE's capability in fitting the UV-to-IR SED of AGNs. We implement the latest version of SKIRTOR, a clumpy two-phase torus model derived from a modern radiative-transfer method (Stalevski et al. 2012(Stalevski et al. , 2016)). In addition, we introduce polar-dust extinction to account for the possible extinction in type 1 AGNs. We name the new version of CIGALE as "X-CIGALE".

The structure of this paper is as follows. In §2, we outline the scheme of our new code development. In §3, we test X-CIGALE on AGNs with X-ray detections from different surveys. We test fitting galaxies with only X-ray upper limits in §4. We summarize our results and discuss future prospects in §5.

Throughout this paper, we assume a flat ΛCDM cosmology with H 0 = 69.3 km s -1 Mpc -1 and Ω M = 0.286 (WMAP 9year results; Hinshaw et al. 2013). Quoted uncertainties are at the 1σ (68%) confidence level, unless otherwise stated. Quoted optical/infrared magnitudes are AB magnitudes.

THE CODE

We briefly summarize the mechanisms and features of CIGALE in §2.1. In §2.2, §2.3, and §2.4, we detail our new development of X-CIGALE, i.e. the X-ray fitting, SKIRTOR, and the polar-dust extinction. The new inputs/outputs introduced in X-CIGALE are listed in Appendix A.

A brief introduction of CIGALE

CIGALE is an efficient SED-fitting code which has been developed for more than a decade (Burgarella et al. 2005;Noll et al. 2009;Serra et al. 2011;Boquien et al. 2019). CIGALE is written in Python. X-CIGALE is built upon CIGALE, and the fitting algorithm of X-CIGALE is the same as that of CIGALE. Here, we only briefly introduce the algorithm, and interested readers should refer to Boquien et al. (2019) for a detailed description.

CIGALE allows the user to input a set of model parameters. The code then realizes the model SED for each possible combination of the model parameters, and convolves the model SED with the filters to derive model fluxes. By comparing the model fluxes with the observed fluxes, the code computes likelihood as L = exp(-χ 2 /2) for each model. CIGALE supports two types of analyses, i.e. maximum likelihood (minimum χ 2 ) and Bayesianlike. In the maximum-likelihood analyses, CIGALE picks out the model with the largest L value, and calculates physical properties such as M and star formation rate (SFR) from this single model. In the Bayesian-like analyses, for each physical property, CIGALE calculates the marginalized probability distribution function (PDF) based on the L values of all models. Finally, from this PDF, CIGALE derives the probability-weighted mean and standard deviation, and outputs them as the estimated value and uncertainty.

Among the above processes, one key step is the realization of model SEDs from input parameters. This procedure relies on a set of modules, and each module is responsible for a function that shapes the SED. For example, the "nebular emission" module adds the nebular-emission components to the SED, and the "dust attenuation" module extincts the SED. Our new development of X-CIGALE follows this module-based structure. We enable CIGALE to fit X-ray data by developing a new X-ray module ( §2.2); we implement SKIR-TOR templates and polar-dust extinction in a new SKIRTOR module ( §2.3 and §2.4).

The new X-ray module

In this section, we develop a new X-ray module to enable X-CIGALE to fit X-ray data. In §2.2.1, we detail the basic settings of this new module. In §2.2.2, we present the adopted X-ray SED for AGN and galaxy components. In §2.2.3, we present the relation that we used to link AGN X-ray with other wavelengths. We note that our new developments are for the majority AGN population in optical/X-ray surveys, and thus X-CIGALE may not be applicable to some minor populations such as radio-loud and broad absorption line (BAL) objects (e.g. Brandt et al. 2000;Miller et al. 2011). We leave the treatment of these particular AGNs to future works.

Basic settings

As presented in §1, the X-ray band has many advantages in studying AGNs. Therefore, we implement an X-ray module for X-CIGALE. The main goal of this module is to connect X-ray with other wavelengths, rather than to obtain detailed X-ray spectral properties (e.g.

X-CIGALE 3

photon index and hydrogen column density) by performing detailed X-ray spectral analyses. This is because the latter has already been well realized by many specialized X-ray codes such as XSPEC (Arnaud 1996) and Sherpa [START_REF] Freeman | Proc. SPIEVol. 4477, Astronomical Data Analysis[END_REF], and there is no need for X-CIGALE to perform similar analyses. Also, it is technically difficult to fit the X-ray spectra within the framework of X-CIGALE. X-CIGALE assumes that a sample of sources are observed with a single "filter transmission", as is the case in UV-to-IR data. However, at X-ray wavelengths, the transmission curve varies from source-to-source, as it might depend on many factors such as position on the detectors and observation date. For example, the soft-band transmission of Chandra has been continuously declining since its launch (e.g. [START_REF] O'dell | Proc. SPIE. p. 103970C[END_REF]. In fact, each source is associated with a unique transmission curve and the curve is taken into account when fitting the X-ray spectra with, e.g. XSPEC and Sherpa.

Therefore, the X-ray module of X-CIGALE is designed to work on the high-level X-ray data products, i.e. intrinsic X-ray fluxes in a given band. We require the X-ray fluxes to be corrected for telescope transmission. Fortunately, this correction is embedded in routine X-ray data processing and has already been applied in X-ray photometric catalogs (e.g. Yang et al. 2016;Luo et al. 2017). Since the transmission has already been considered, X-CIGALE only needs to adopt a uniform-sensitivity (i.e. boxcar-shaped) "filter". We have already included a few typical boxcar X-ray filters, e.g. 0.5-2 keV and 2-7 keV for convenience, while the user can easily generate the filters for any X-ray band.

In addition, we require the input X-ray fluxes to be "absorption-corrected". The absorption might be from the source itself, our Galaxy, and/or the intergalactic medium (IGM; e.g. Starling et al. 2013). However, we do not differentiate these types of absorption, as it is often infeasible to separate them in X-ray data analyses. The absorption correction can be obtained from routine X-ray data processing, e.g., spectral analyses via XSPEC/Sherpa or band-ratio analyses (e.g. Xue et al. 2016;Luo et al. 2017). The user may also choose to use hard X-ray bands where absorption corrections are generally small (e.g. Yang et al. 2018a,b).

In X-ray catalogs (e.g. Xue et al. 2016;Luo et al. 2017), X-ray fluxes ( f X ) are conventionally given in the cgs units of erg s -1 cm -2 , but X-CIGALE requires the input fluxes ( f cigale ) to be given in the units of mJy. Therefore, the user needs to convert the flux units with

f cigale = f x × 4.136 × 10 8 E up -E lo (1)
where E lo and E up refer to the lower and upper limits of the energy band in units of keV. The X-ray module covers rest-frame 10 -3 -5 nm, corresponding to ≈ 0.25-1200 keV. Such an energy range is sufficient for practical purposes: current X-ray instruments cannot observe energies significantly below rest-frame 0.5(1 + z) keV in general; the AGN flux is typically non-detectable above ≈ 1000 keV due to the existence of the cut-off energy in AGN X-ray spectra (see §2.2.2).

X-ray SED

To first-order approximation, the intrinsic AGN X-ray spectrum is typically a power law with a high-energy exponential cutoff, i.e.

f ν ∝ E -Γ+1 exp(-E/E cut ) (2)
where Γ is the so-called "photon index", widely adopted in X-ray astronomy, and E cut is the cutoff energy. We adopt this spectral shape in X-CIGALE. Detailed X-ray spectral fitting in the literature finds Γ ≈ 1.8 (e.g. Yang et al. 2016;Liu et al. 2017). We allow the user to set Γ in X-CIGALE. We set E cut = 300 keV, the typical value from the observations of Seyferts (e.g. Dadina 2008;Ricci et al. 2017). Note that since E cut is above the highest observable energy of most X-ray observatories (e.g. Chandra and XMM-Newton), the exact choice of E cut has practically negligible effects on the fitting with X-CIGALE for most cases. The adopted AGN X-ray SED is displayed in Fig. 1. Besides AGNs, galaxies can also emit X-rays, although the emission from galaxies is often much weaker than that from AGNs for X-ray detected sources. There are three main origins of galaxy X-ray emission: low-mass X-ray binaries (LMXB), high-mass X-ray binaries (HMXB), and hot gas. The strengths of these components can be modeled as a function of galaxy properties such as M and SFR. We adopt the recipe from Mezcua et al. (2018), where a Chabrier (2003)

/SFR) = 40.3 -62Z + 569Z 2 -1834Z 3 + 1968Z 4 (3)
where M and SFR are in solar units; t denotes stellar age in units of Gyr; Z denotes metallicity (mass fraction). The hot-gas luminosity (in units of erg s -1 ) is described as log(L hotgas 0.5-2keV /SFR) = 38.9 (4)

Similarly as for AGN, we also employ the SED shape in Eq. 2 for all three components, with E cut fixed at 100 keV (LMXB and HMXB; e.g. Zhang 1997; Motta et al. 2009) and 1 keV (hot gas; e.g. Mathews & Brighenti 2003). We allow the user to set Γ values for the LMXB and HMXB components. In our test fitting in §3, we set Γ to 1.56 and 2.0 for LMXB and HMXB, respectively (e.g. Fabbiano 2006;Sazonov & Khabibullin 2017). Adjusting these Γ does not affect the fitting results significantly, as the observed X-ray fluxes are often dominantly contributed by AGNs rather than galaxies. The X-ray continuum from hot gas can be modelled as free-free and free-bound emission from optically thin plasma (Γ = 1; e.g. Mewe et al. 1986). Therefore, we fix Γ = 1 for the hot-gas component in X-CIGALE. Fig. 1 shows the adopted X-ray SEDs of the three components. We add all three components for the total X-ray SED from galaxies.

The α ox -L 2500Å relation

As in §2.2.1, the main goal of X-CIGALE is to fit X-ray and other wavelengths simultaneously. Some known connections between X-ray and other wavelengths must be applied; otherwise, the fitting would be practically useless. We adopt the well-studied "α ox -L 2500Å " relation (e.g. Steffen et al. 2006;Just et al. 2007;Lusso & Risaliti 2017), where L 2500Å is AGN intrinsic (de-reddened) luminosity per frequency at 2500 Å and α ox is the SED slope between UV (2500 Å) and X-ray (2 keV), i.e.

α ox = -0.3838 log(L 2500Å /L 2keV ). (5)
The observed α ox -L 2500Å relation (Just et al. 2007) is written as where L 2500Å is in units of erg s -1 Hz -1 . The 1σ intrinsic dispersion of this α ox -L 2500Å relation is ∆α ox ≈ 0.1 (see Table 8 of Just et al. 2007). Here, ∆α ox is the α ox deviation from that expected from the α ox -L 2500Å relation, i.e.

α ox = -0.137 log(L 2500Å ) + 2.638 (
∆α ox = α ox -α ox (L 2500Å ). (7)
Observations have found that the α ox -L 2500Å relation does not have significant redshift evolution, indicating that the relation originates from fundamental accretion physics (Steffen et al. 2006;Just et al. 2007;Lusso & Risaliti 2017). We allow the user to set the maximum |∆α ox | allowed (|∆α ox | max ). Internally, X-CIGALE calculates all models with α ox from -1.9 to -1.1 with a step of 0.1. (Just et al. 2007).

Note that the α ox -L 2500Å relation above is derived from observations, assuming that the unobscured AGN emission is isotropic at both UV/optical and X-ray wavelengths. However, the UV/optical emission is unlikely isotropic, because it is from the accretion disk and the effects of projected area and limb darkening affect the angular distribution of the radiative energy. After considering these effects, the disk luminosity can be approximated as L(θ) ∝ cos θ(1 + 2 cos θ), where θ is the angle from the AGN axis (e.g. Netzer 1987). This angular dependence of disk emission is adopted in SKIRTOR, the UV-to-IR AGN module adopted in X-CIGALE (see §2.3). The X-ray emission should likely be less anisotropic than the UV/optical emission, because the X-rays originate from re-processed UV/optical photons via inverse Compton scattering. However, the exact relation between X-ray flux and viewing angle depends on model details, such as corona shape and opacity, which are poorly known (e.g. Liu et al. 2014;Xu 2015). For simplicity, we assume that the X-ray emission is isotropic.

Our assumption of anisotropic UV/optical emission and isotropic X-ray emission leads to a dependence of α ox -L 2500Å on 0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 log (nm, rest-frame) The "breaks" at 5 nm are caused by the wavelength limit of the X-ray module, which ends at 5 nm. As explained in §2.2.3, such breaks are generally not problematic for practical purposes, as most X-ray instruments only cover wavelengths 10 0.4 nm ( 0.5 keV, as marked).

viewing angle. We further assume that the observed α ox -L 2500Å relation for type 1 AGNs reflects the intrinsic α ox -L 2500Å relation for all AGNs at a "typical" viewing angle of θ = 30 • . This value approximates the probability-weighted viewing angle for type 1 AGNs, i.e.

θ ≈ ∫ 90 • -∆ 0 θ sin θdθ ∫ 90 • -∆ 0 sin θdθ ≈ 30 • , (8)
where ∆ denotes the angle between the equatorial plane and edge of the torus, i.e., half opening angle. The typical value is ∆ ≈ 40 • from observations (e.g. Stalevski et al. 2016). Although ∆ is a free parameter in X-CIGALE, we do not recommend the user choose other values than 40 • , as this value is favored by observations and is consistently adopted throughout the build-up of the X-CIGALE code. The weight sin θ is proportional to the probability for the viewing angle being θ.

We note that our SED fitting results ( §3 and §4) are not sensitive to the assumed typical θ, and will not change significantly if adjusting θ within the range of ≈ 10 • -50 • . In the X-CIGALE output (Appendix A), the α ox and L 2500Å always refer to the value at θ = 30 • , regardless of the actual viewing angle in the model. This α ox and L 2500Å design is to reflect AGN essential properties, independent of the viewing angle. By changing the integral ranges in Eq. 8 to (90 • -∆, 90 • ), we can derive the probability-weighted θ for type 2 AGNs, i.e. θ ≈ 70 • . These typical θ values (type 1: ≈ 30 • , type 2:≈ 70 • ) are used in our SED fitting ( §3 and §4).

SKIRTOR

The previous CIGALE AGN model responsible for the UV-to-IR SED is from [START_REF] Fritz | [END_REF]. This model assumes that the dusty torus is a smooth structure. However, more recent theoretical and observational works find that the torus is mainly made of dusty clumps (e.g. Nikutta et al. 2009;Ichikawa et al. 2012;Stalevski et al. 2012;Tanimoto et al. 2019). SKIRTOR is a clumpy twophase torus model (Stalevski et al. 2012(Stalevski et al. , 2016)), based on the 3D radiative-transfer code, SKIRT (Baes et al. 2011;Camps & Baes 2015). 2 In SKIRTOR, most (mass fraction = 97%) of the dust is in the form of high-density clumps, while the rest is smoothly distributed. In addition, SKIRTOR considers the anisotropy of the power source, AGN disk emission (see §2.2.3), while Fritz's model simply assumes isotropic disk emission. Therefore, we implement SKIRTOR within X-CIGALE. We recommend using SKIRTOR as the UV-to-IR SED model of AGNs, although X-CIGALE allows the user to choose between SKIRTOR and Fritz's model.

SKIRTOR adopts a disk SED that has a higher fraction of far-UV luminosity ( 100 nm) compared to observations (see §3.2.1 of Duras et al. 2017). Following Duras et al. (2017), we update SKIRTOR with a new disk SED (Feltre et al. 2012) that is supported by observations, i.e.

λL λ ∝              λ 2 8 ≤ λ < 50 [nm] λ 0.8 50 ≤ λ < 125 [nm] λ -0.5 125 ≤ λ < 10 4 [nm] λ -3 λ > 10 4 [nm]. (9)
We modify the disk SED with the following method. We denote the old and new intrinsic disk SEDs as L old λ,normed and L new λ,normed , respectively, where the subscript "normed" indicates the total power of these SEDs has been normalized to unity. Then the new observed disk SED component (which might be obscured) can be converted from the old one by multiplying by the factor, L new λ,normed /L old λ,normed . The new scatter component can be obtained in the same way; the dust reemitted component remains unchanged. The method above keeps energy balance. This method is also described on the SKIRTOR official webpage.3

Polar Dust

The extinction of type 1 AGN

In SKIRTOR (also in Fritz's model), the extinction of UV and optical radiation for type 1 AGN is assumed to be negligible. This assumption holds for most optically selected blue quasars. For example, Richards et al. (2003) found only ≈ 6% of their SDSS quasars are extincted. However, the assumption might not be true for, e.g. X-ray selected AGNs. For example, in the COSMOS AGN catalogs selected by XMM-Newton (Bongiorno et al. 2012), the fraction of extincted sources (E(B -V) ≥ 0.1) among broad-line AGNs is ≈ 40%.

To check the extinction of type 1 AGNs, we compare the median UV-optical SEDs of spectroscopically classified type 1 AGNs in SDSS and COSMOS (see §3 for details). These median SEDs are derived from the photometric data in §3. For each source in a sample, we interpolate the observed photometry to obtain F ν as a function of observed-frame wavelength. We then shift this interpolated SED to rest-frame wavelength and normalize F ν at 250 nm. Finally, at each wavelength, we obtain the median F ν of all the sources in the sample. Fig. 3 (left) shows the results. The SDSS median SED is similar to the typical unobscured quasar SED of F ν ∝ λ 0.5 (see Eq. 9). In contrast, the COSMOS median SED is significantly redder than F ν ∝ λ 0.5 (e.g. Elvis et al. 2012). We note that this difference in SED shape is observationally driven by selection effects. The SDSS sample consists of optically selected, and is thus biased toward blue and optically bright objects. The COSMOS sample consists of X-ray selected objects and does not suffer from significant bias in the UV/optical (e.g. Brandt & Alexander 2015). Although driven by selection effects, Fig. 3 (left) at least indicates that reddened AGN SEDs indeed exist, and we discuss the physical cause of the SED reddening below.

The red SED shape might be physically caused by the aforementioned dust extinction. However, another potential physical cause is host-galaxy contribution to the SED. Since the UV/optical SEDs of galaxies are generally redder than those of unobscured AGNs (e.g. Fig. 3 of Salvato et al. 2009), AGN-galaxy mixed SEDs tend to be redder than pure AGN SEDs. To investigate the cause of SED reddening, we can compare the magnitudes that sample restframe UV wavelengths for COSMOS and SDSS, as galaxy contributions to the photometry should be small at UV wavelengths. 4 In Fig. 3 (right), we show the r-band (rest-frame 2000 Å) magnitude distributions at z = 2 -2.5 and f ν,2keV = 3-10 (10 -7 mJy, the results are similar for other redshift/X-ray flux bins). The control of redshift and X-ray flux is to force the compared samples to have similar X-ray luminosities and thereby bolometric luminosities (assuming the X-ray bolometric correction factors are similar for the sources). The COSMOS AGNs are systematically fainter than the SDSS AGNs in UV/optical. Therefore, Fig. 3 (right) indicates that, at a given AGN bolometric luminosity, the rest-frame UV AGN luminosities of COSMOS are typically lower than those of SDSS, supporting the existence of dust extinction. We conclude that dust extinction is at least one of the physical causes of the SED reddening (Fig. 3 left), although galaxy SED contributions might enhance the reddening (e.g. Bongiorno et al. 2012).

The polar-dust model

From §2.4.1, it is necessary to account for dust extinction of type 1 AGNs in X-CIGALE. The geometry of the obscuring materials is sketched in Fig. 4, where the materials responsible for type 1 AGN obscuration are called "polar dust" (e.g. Lyu & Rieke 2018). The existence of polar dust has been proved by high-resolution mid-IR (MIR) imaging of local Seyfert galaxies (e.g. López-Gonzaga et al. 2014;Stalevski et al. 2017Stalevski et al. , 2019;;Asmus 2019). However, the physical properties of the polar dust could be complicated and vary for different objects. For example, it might be close to the dust-sublimation radius (∼ pc scale; e.g. Lyu & Rieke 2018) or on galactic scales (∼ kpc; e.g. Zou et al. 2019).

Considering these complexities, we do not build a grid of physical models and perform radiation-transfer simulations. Instead, we employ several empirical extinction curves, including those from Calzetti et al. (2000, nearby star-forming galaxies), Gaskell et al. (2004, large dust grains), and Prevot et al. (1984, Small Magellanic Cloud, SMC), and the user can choose among these curves. In our tests of X-CIGALE ( §3 and §4), we adopt the SMC extinction curve, which is preferred from AGN observations (e.g. Hopkins et al. 2004;Salvato et al. 2009;Bongiorno et al. 2012; but also see, e.g. Gaskell et al. 2004). The extinction amplitude (parameterized as E(B -V)) is a free parameter set by the user, and setting E(B -V) = 0 returns to the original torus.

Since the scheme of X-CIGALE maintains energy conservation, we need to implement dust emission to account for the radia-2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 log (nm, rest-frame) f , 2keV = 3 10 (10 7 mJy), z=2.0-2.5 SDSS COSMOS Figure 3. Left: Median SEDs of spectroscopic type 1 AGNs in SDSS (blue, optically selected) and COSMOS (orange, X-ray selected). The green curve shows a typical unobscured quasar SED. All the SEDs are normalized at 250 nm. The SDSS SED is similar to the typical quasar SED. The COSMOS SED is significantly redder than the SDSS SED, indicating the presence of dust extinction and/or host-galaxy contribution for the COSMOS sources. Right: The observed r-band magnitude distributions for SDSS and COSMOS type 1 AGNs. Both AGN samples are in the bins of z = 2-2.5 and f ν,2keV = 3-10 (10 -7 mJy). The SDSS AGNs are systematically brighter than the COSMOS AGNs, and this qualitative result also holds for other redshift/X-ray flux bins. This result indicates that dust extinction should be at least one of the causes for the red SED shapes of the COSMOS AGNs, although galaxy SED contributions might enhance the reddening. This figure does not compare the AKARI-NEP sample in §3.3, because spectroscopic identifications of type 1 AGNs are not available for AKARI-NEP ( §3.3.1).

tive energy absorbed by the dust. We assume the dust reemission is isotropic. For simplicity, we adopt the "grey body" model (e.g. Casey 2012), i.e.

L ν (λ) ∝ 1 -e -(λ 0 /λ) β c λ 3 e hc/λkT -1 , ( 10
)
where λ 0 is fixed at 200 µm, emissivity (β) and temperature (T) are free parameters set by the user. L ν in Eq. 10 is normalized so that total energy is conserved, i.e.

L emit total = ∫ 90 • -∆ 0 L extinct total (θ) sin θdθ (11)
where L emit total is the dust reemitted luminosity (angle-independent) and L extinct total is the luminosity loss caused by polar-dust extinction (angle-dependent). Note that the integral on the right-hand-side of Eq. 11 is to account for the fact that the polar dust only accounts for the obscuration in the polar directions while the polar-dust reemission is in all directions (see Fig. 4; e.g. Eq. 8). Fig. 5 shows the model SEDs for different extinction levels, where T = 100 K, β = 1.6, and ∆ = 40 • .

Our model above follows the AGN-unification scheme, i.e. AGN type is determined solely by the viewing angle, which is a free parameter in X-CIGALE. When the viewing angle is within the polar directions (type 1), the observed AGN disk emission suffers moderate (or none if E(B -V) = 0) extinction from the polar dust. When the viewing angle is within the equatorial directions (type 2), the observed AGN disk emission is strongly obscured by the torus. If the AGN type is known (e.g. from spectroscopy), the user can limit the viewing angle to the polar or equatorial direction ( §3.1 and §3.2). Otherwise, the user can adopt multiple viewing angles including both polar and equatorial directions, and let X-CIGALE choose freely between them ( §3.3). For example, if the user set "viewing angles = 30 • , 70 • ; E(B -V) = 0.1", then CIGALE will build two model SEDs. For the 30 • model (type 1), the UV/optical SED is reddened by the E(B-V) = 0.1 polar dust whose reemission also contributed to the IR SED. For the 70 • model (type 2), the polar dust does not affect the UV/optical SED (already obscured by torus), but its reemission still contributes to the IR SED.

Our polar-dust model above provides one possible scenario for the reddened type 1 AGNs, i.e. the viewing angle is small and the extinction is caused by dust along the polar directions. An alternative scenario is that the line-of-sight (LOS) intercepts the torus, but the extinction is only moderate by chance due to the inhomogeneity of torus. However, this scenario has not been well investigated with physical torus models in the literature, to our knowledge. Therefore, we focus on the polar-dust model in the current version of X-CIGALE, and future versions of X-CIGALE may include this alternative scenario when its SED templates are available.

TESTS ON X-RAY DETECTED SOURCES

In this section, we test X-CIGALE with three samples of AGNs, i.e. SDSS ( §3.1), COSMOS ( §3.2), and AKARI-NEP ( §3.3). The basic properties of these samples are summarized in Table 1. These three samples have different characteristics. The SDSS sample is optically bright type 1 quasars. The COSMOS sample is X-ray selected AGNs with broad multiwavelength coverage from u to Herschel/PACS 160 µm. This sample also has spectroscopic AGN classifications. The AKARI-NEP sample is small but has excellent MIR observations from AKARI.

SDSS

The sample and the models

The SDSS sample is optically selected from the DR14 quasar catalog [START_REF] Pâris | [END_REF]. All the sources are spectroscopically confirmed type 1 AGNs. In addition to the SDSS ugriz bands, the Pâris et al. ( 2018) catalog provides X-ray data from XMM-Newton archival observations when available (the 3XMM (3) r-band AB magnitude range (20%-80% percentile). ( 4) 2-10 keV X-ray flux range (20%-80% percentile) in units of 10 -14 erg s -1 cm -2 . (5) Redshift range (20%-80% percentile). ( 6) Types of AGNs included in the survey.

catalog; Rosen et al. 2016). We require the sources to be detected in the 2-12 keV band at > 3σ significance levels. Here, the choice of the hard X-ray band (2-12 keV) is to minimize the effects of X-ray obscuration (see §2.2.1). 5 We do not include the IR photometry compiled in the DR14 catalog, because our main goal here is to test X-CIGALE on the simple cases, i.e. the quasar-dominated SED. In the X-ray to optical wavelengths, the AGN component is dominant; but in the IR wavelengths, the galaxy component may be non-negligible. We discuss the cases of AGN-galaxy mixed SEDs 5 Alternatively, one can perform X-ray spectral fitting to obtain the absorption-corrected X-ray fluxes. However, extracting and analyzing the X-ray spectra from the public XMM-Newton archival data are beyond the scope of this work. in §3.2 and §3.3. We require the sources to have Galactic extinctions estimated in the DR14 quasar catalog, because we need to correct for Galactic extinction before providing the photometry to X-CIGALE. These criteria lead to a final sample of 1986 AGNs (Table 1).

For the SDSS sample, we can neglect the galaxy SED component, because the sources are optically bright quasars which often dominate the observed UV/optical SEDs. The AGN-dominant (X-ray to IR) models in X-CIGALE can be achieved by setting frac AGN to a value close to unity (e.g. 0.999). 6 The adopted AGN model parameters are listed in Table 2. The only free parameter in our fitting is polar-dust E(B -V), which affects the UV/optical SED shape. We further justify that it is necessary to have E(B -V) as a free parameter in §3.1.2. Other SKIRTOR parameters are fixed, because they only affect the IR SED shape where there is no band coverage for the SDSS sample (see §3.3.2 for the assessment of these parameters).

For the X-ray module, we adopt Γ = 1.8 for AGN (the dominant component in X-rays), the typical intrinsic photon index constrained by observations ( §2.2.2). Adopting other AGN Γ values (e.g. 1.4 or 2.0) do not affect our fitting results significantly. Our adopted Γ = 1.8 is slightly different from that assumed in the 3XMM catalog (Γ = 1.7; Rosen et al. 2016). Therefore, we scale the 2-12 keV fluxes by a factor of 0.96 to correct the effects of different Γ, and this correction factor is obtained using the PIMMS website. 7 For the LMXB and HMXB components, we set Γ = 1.56 and 2.0, respectively (see §2.2.2). We adopt |∆α ox | max = 0.2, and this |∆α ox | max value is ≈ 2σ scatter of the α ox -L 2500Å relation ( §2.2.3). Note that although the X-ray module has both parameters fixed, X-CIGALE internally calculates 9 models of different α ox values and selects |∆α ox | ≤ |∆α ox | max (see §2.2.3).

Fitting results

We run X-CIGALE with the model settings in §3.1.1 for the SDSS sample. The median reduced χ 2 ( χ 2 red ) and degrees of freedom (dof) are 1.4 and 5, respectively. These χ 2 red and dof values correspond to a p-value of 23%, well above the conventional 2σ (5%) or 3σ (0.3%) values. This result indicates that X-CIGALE is able to model the observed photometry of the SDSS quasars. Fig. 6 shows a random example of the SED fitting. Fig. 7 displays the E(B -V) distribution from the fitting. As expected (see §2.4.1), most (75%) SDSS AGNs have weak or no extinction with E(B -V) ≤ 0.1.

To evaluate the effects of the new X-ray module, we re-run X-CIGALE but without this module. We compare the AGN intrinsic L 2500Å between the fitting with X-ray (L 2500Å,X ) vs. without X-ray (L 2500Å,noX ) in Fig. 8. The L 2500Å,X and L 2500Å,noX are similar, and this similarity is as expected. SDSS sources are mostly unobscured type 1 AGNs due to their selection method ( §2.4.1), and thus the intrinsic AGN emission is directly observable at UV/optical wavelengths. Therefore, adding the X-ray module does not significantly change L 2500Å estimation for SDSS sources in general.

In Table 2, the only free parameter is polar-dust E(B -V), because this parameter affects the UV/optical SED which is covered by the SDSS bands ( §3.1.1). In other words, we consider that E(B -V) can be constrained by the photometric data. The constrainability of a model parameter can be evaluated by the "mock analysis" of X-CIGALE, which already exists in the previous version of CIGALE (see §4.3 of Boquien et al. 2019 for details). Briefly, after fitting the observed data, X-CIGALE simulates a mock catalog based on the best-fit model for each object. The photometric uncertainties are considered when simulating the mock data. X-CIGALE then performs SED fitting to the mock catalog and obtains Bayesian-like estimated values (and their errors) of the parameter. By comparing these estimated values and those used to generate the mock catalog (i.e. the "true" values), one can assess whether the parameter can be reliably constrained. The mock analysis serves as a sanity check to assess whether a physical parameter can be retrieved in a self-consistent way. This mock analysis can be invoked by setting "mock_flag=True" in the X-CIGALE configurations.

We run the mock analysis to test if polar-dust E(B -V) can be constrained. We compare the estimated and true values in Fig. 9 (left). The estimated and true values are well correlated, indicating that E(B -V) can be self-consistently constrained. In Fig. 10, we show the PDF of E(B -V) for the source in Fig. 6. Fig. 10 indicates that the E(B -V) is indeed well constrained in the Bayesian-like analysis.

In our fitting, aside from the model normalization (automatically determined by X-CIGALE; see §4.3 of Boquien et al. 2019), E(B -V) is the only free model parameter (Table 2). We also test freeing other parameters such as viewing angle and torus optical depth, and the mock-analysis results of E(B -V) are similar. For example, Fig. 9 (right) shows the result after setting the viewing angle to 0-90 • with a step of 10 • (i.e. all allowed values). The estimated and true values are still well correlated, indicating that E(B -V) and viewing angle are not strongly degenerate. This non-degeneracy is understandable, because, in our polar-dust model, E(B -V) is the only parameter responsible for modelling the observed UV/optical SED shapes of type 1 AGNs like the SDSS objects ( §2.4).

COSMOS

The sample and the models

The COSMOS sample is X-ray selected (> 3σ, 2-10 keV band) from the COSMOS-Legacy survey performed by Chan- Torus optical depth at 9.7 microns τ 9.7 7.0 Torus density radial parameter p (ρ ∝ r -p e -q | cos(θ )| ) 1.0 Torus density angular parameter q (ρ ∝ r -p e -q | cos(θ )| ) 1.0 Angle between the equatorial plane and edge of the torus ∆ 40 • Ratio of the maximum to minimum radii of the torus 20 Viewing angle θ (face on: θ = 0 • , edge on: θ = 90 dra (Civano et al. 2016). The COSMOS-Legacy catalog assumes Γ = 1.4. Similarly as in §3.1.1, we apply a correction factor of 0.87 (calculated with PIMMS) to the 2-10 keV fluxes to make them consistent with our adopted Γ = 1.8 for AGN. Marchesi et al. (2016) matched these X-ray sources with the optical/NIR counterparts in the COSMOS2015 catalog (Laigle et al. 2016) and compiled their spectroscopic information when available. We select the sources with spectroscopic classifications of AGN types. We adopt the photometric data in the COSMOS2015 catalog, including 14 broad bands from u to IRAC 8.0 µm. In addition, when available, we also include photometric data from Spitzer/MIPS (24 µm) and Herschel/PACS (100 µm and 160 µm), from the PEP survey (Lutz et al. 2011). We adopt the redshift measurements from Marchesi et al. (2016). These redshifts are either secure spectroscopic redshifts or high-quality photometric redshifts. There are a total of 590 objects in COSMOS (Table 1). Among these 590 objects, 206 and 384 are type 1 and type 2 AGNs, respectively.

The X-CIGALE model parameters for COSMOS are listed in Table 3. For the SKIRTOR and X-ray modules, the parameter setting is the same as in §3.1.1 except for frac AGN and the viewing angle. Here, we allow frac AGN to vary among 0.01, 0.1-0.9 (step 0.1), and 0.99, because, unlike in the case of SDSS, the AGN contribution to the observed SED may not be generally not dominant for the COS-MOS sample, especially in the IR bands. We set the viewing angle to 30 • and 70 • for the spectroscopic type 1 and type 2 AGNs, respectively. These values are approximately the probability-weighted θ for type 1 and type 2 AGNs, respectively, given a torus of ∆ = 40 • (see 2.2.3).

For the galaxy component, we adopt the model setting similar to that in Ciesla et al. (2015). Specifically, we adopt a delayed star-formation history (SFH), because it can characterize the SEDs of both early-type and late-type galaxies reliably (e.g. Ciesla et al. 2015;Boquien et al. 2019). Also, the delayed SFH only has relatively small parameter space (only two free parameters) and thereby high fitting efficiency. We adopt a Chabrier ( 2003) IMF with metallicity (Z) fixed to the solar value of 0.02. For the galactic dust attenuation, we adopt the dustatt_calzleit module in X-CIGALE (Calzetti et al. 2000;Leitherer et al. 2002). The allowed E(B -V) values for young stars are 0.1, 0.2, 0.3, 0.4, 0.5, 0.7, and 0.9. The E(B -V) ratio between the old and young stars is fixed to 0.44. The amplitude of the 217.5 nm UV bump on the extinction curve is set . The E(B -V ) PDF for the example source in Fig. 6. The PDF is normalized such that its integral is unity.

to 0 (SMC) and 3 (Milky Way). We adopt the Dale et al. (2014) model for galactic dust reemission. There is only one free parameter in this model, i.e. the α slope in dM dust ∝ U -α dU, where M dust and U are dust mass and radiation-field intensity, respectively. The α values are set to 1.5, 2.0, and 2.5.

Fitting results

We run X-CIGALE with the parameter settings in §3.2.1. The median χ 2 red values are 1.4 and 0.9, for type 1 and type 2 AGNs, respectively, while the median dof are 15 for both types. These median χ 2 red and dof corresponding to p-values of 0.12 and 0.59. These relatively large p-values for both type 1 and type 2 indicates that our models ( §2) are able to model AGN SEDs of different types. This result supports the AGN-unification scheme ( §1), on which our models are based. Fig. 11 displays two examples of the SED fitting in COSMOS. In Fig. 7, we compare the polar-dust E(B -V) of type 1 AGNs in COSMOS vs. SDSS. The COSMOS type 1 AGNs tend to have higher E(B -V) than SDSS type 1 AGNs, consistent with the diagnostic in §2.4.1. Fig. 12 compares L 2500Å,X and L 2500Å,noX for the COSMOS sample. The differences between L 2500Å,X and L 2500Å,noX are larger compared to those in SDSS (Fig. 8). This is because, for the SDSS sources, the observed optical fluxes are dominated by the AGN component, and thus AGN power can be effectively constrained even without X-ray data. In contrast, for COSMOS, the observed opticalto-IR fluxes are often not dominated by AGN, and X-CIGALE needs to decompose the fluxes into galaxy and AGN components. This SED decomposition process may be sometimes difficult, given that different models could result in similar model fluxes in opticalto-IR SED. Therefore, the X-ray data, which is often dominated by AGN, can be helpful in constraining the AGN power. Fig. 13 shows an example type 2 AGN SED fitted with vs. without X-ray. For this source, the observed UV-to-IR fluxes are dominated by the galaxy component, and thus AGN power cannot be effectively constrained without X-ray data.

In our fitting (Table 3), we set viewing angles at 30 • and 70 • for type 1 and type 2 AGNs, respectively. This model setting can be done for our COSMOS sample, where spectroscopic classification is available ( §3.2.1). However, for many photometric surveys (e.g. §3.3.1), the AGN spectra are not available, and thus spectrumbased AGN-type classification cannot be performed. In this case, the X-CIGALE user can set both 30 • and 70 • for the viewing angle, and allow X-CIGALE to freely choose between them. We test this configuration with our entire COSMOS sample, including spectroscopic type 1 and type 2 AGNs. The other parameters are the same as in Table 3. For spectroscopic type 1 (type 2) AGNs, 70% (28%) and 30% (72%) sources have the best-fit viewing angles of 30 • and 70 • , respectively. This means that, if one uses the best-fit viewing angle to perform AGN-type classification (i.e. the SED-based classification), the correct rate will be roughly ≈ 70% for both type 1 and type 2. E(B -V ) of starlight for the young population 0.1, 0.2, 0.3, 0.4, 0.5, 0.7, 0.9 E(B -V ) ratio between the old and young populations 0.44 Galactic dust emission: Dale et al. (2014) α slope in dM dust ∝ U -α dU 1.5, 2.0, 2.5

AGN (UV-to-IR): SKIRTOR

Torus optical depth at 9.7 microns τ 9.7 7.0 Torus density radial parameter p (ρ ∝ r -p e -q | cos(θ )| ) 1.0 Torus density angular parameter q (ρ ∝ r -p e -q | cos(θ )| ) 1.0 Angle between the equatorial plan and edge of the torus 40 • Ratio of the maximum to minimum radii of the torus 20 Viewing angle θ (face on: θ = 0 • , edge on: θ = 90 • ) a 30 • (type 1), 70 • (type 2) AGN fraction in total IR luminosity frac AGN 0.01, 0.1-0.9 (step 0.1), 0.99 Extinction law of polar dust SMC E(B -V ) of polar dust 0, 0.05, 0.1, 0.15, 0.2, 0.3 Temperature of polar dust (K) 100 Emissivity of polar dust 1.6 X-ray: This work

AGN photon index Γ 1.8 Maximum deviation from the α ox -L 2500Å relation 0.2 LMXB photon index 1.56 HMXB photon index 2.0 N . -(a)
For COSMOS, the viewing angles are set to 30 • and 70 • for the spectroscopic type 1 and type 2 AGNs, respectively. For AKARI-NEP, we allow X-CIGALE to choose between 30 • and 70 • for the entire sample, since spectroscopic classifications are not available.

AKARI-NEP

The sample and the models

The AKARI-NEP sample is also X-ray selected (> 3σ, 2-7 keV band) based on Chandra observations of the AKARI-NEP field (Krumpe et al. 2015). The Krumpe et al. (2015) catalog assumes Γ = 1.4. Similarly as in §3.1.1 and §3.2.1, we apply a correction factor of 0.94 (calculated with PIMMS) to the 2-7 keV fluxes to make them consistent with our adopted Γ = 1.8 (AGN). We match the X-ray with the multiwavelength catalog compiled by Buat et al. (2015) using a 1 matching radius. This multiwavelength catalog has 19 bands from u to Herschel/PACS (100µm and 160µm). Notably, these bands include an excellent set of 9-band MIR data from the AKARI telescope, allowing us to test X-CIGALE on the MIR wavelengths.

The final sample has 74 sources (Table 1).

For AKARI-NEP, we adopt the same X-CIGALE fitting parameters as for COSMOS (see Table 3) except the viewing angle. For the viewing angle, we allow X-CIGALE to choose freely between 30 • (type 1) and 70 • (type 2), since spectroscopic classifications of AGN type are not available for the AKARI-NEP sample.

Fitting results

We run X-CIGALE with model parameters in §3.3.1. The median χ 2 red is 1.2 (median dof= 17). The resulting p-value is 0.27, indicating overall good fitting quality. Fig. 14 shows two example fitted SEDs in the AKARI-NEP sample. Note that the MIR data can be well fitted with our model. Fig. 15 compares L 2500Å,X and L 2500Å,noX for AKARI-NEP. The differences between L 2500Å,X and L 2500Å,noX are larger compared to those in SDSS (Fig. 8). The reason is similar as discussed in §3.2.2, i.e. SED decomposition is needed for AKARI-NEP and such decomposition might be ambiguous without X-ray.

Parameter constrainability and degeneracy

The AKARI-NEP sample is suitable for investigating the constrainability and degeneracy of AGN model parameters, thanks to its excellent coverage at MIR wavelengths ( §3.3.1) where AGNdust emission peaks (Fig. 5). In this section, we present discussion on model-parameter constrainability and degeneracy using the AKARI-NEP sample. Since the population in this sample is AGN-galaxy mixed systems in general (e.g. Fig. 14), the results below might not be applicable to some particular sources, for which the observed MIR emission is dominated by AGNs (e.g. hot dustobscured galaxies, hot DOGs; Dey et al. 2008;Vito et al. 2018). For these AGN-dominated sources, the AGN parameters might be easier to constrain, as host-galaxy contributions to the MIR emission are negligible.

There are three parameters that have multiple values in our fitting (Table 3), i.e. viewing angle, frac AGN , and polar-dust E(B -V). The viewing angle determines AGN types (see §2.4.2), and we show that the spectroscopic AGN types can be recovered with ≈ 70% accuracy in §3.2.2. For frac AGN and E(B -V), we run the mock analysis as a sanity check of their constrainability ( §3.1.2), and the results are displayed in Fig. 16. For frac AGN , the estimated and true values are generally correlated (median errors = 0.13), although some sources have relatively large uncertainties. Therefore, the relative IR emission strength between AGN and galaxy (as measured by frac AGN ) can be effectively constrained. In contrast, unlike the case of SDSS ( §3.1.2), the estimated E(B -V) is relatively flat as a function of true E(B -V), indicating that E(B -V) cannot be well constrained in general. This result is understandable. For SDSS, the SED is dominated by type 1 AGNs, and the E(B -V) is directly related to the observed UV/optical SED shape. However, for AKARI-NEP, the SEDs are generally produced by both AGN and galaxy components, and the E(B -V) cannot be determined directly from the UV/optical SED shape (or other SED features).

In our fitting (Table 3), most of the torus and polar-dust parameters such as τ 9.7 and polar-dust temperature are fixed at single values. This is because these parameters are related to the MIR SED shape. Considering that model degeneracy is likely strong in the MIR, broad-band photometry data like AKARI-NEP might not be able to effectively constrain these parameters. Now, we test whether torus τ 9.7 can be well constrained. The X-CIGALE configuration is the same as in Table 3 except that τ 9.7 is allowed to vary among 3, 5, 7, 9, and 11 (all allowed values). The mock-analysis results are presented in Fig. 17 in Table 3 such as polar-dust temperature and torus opening angle, and found they cannot be effectively constrained either. Now, we analyze the reasons why τ 9.7 has generally large uncertainties, using the two sources in Fig. 14 as illustrative examples. These reasons also generally explain the large uncertainties of other unconstrained parameters. In Fig. 18, we show the τ 9.7 PDFs and 2D probability density maps (τ 9.7 vs. frac AGN ).

For J175535.47+660959.0, the frac AGN can be constrained to 0.46 ± 0.13. From the density map, at high frac AGN (≈ 0.6), the probability peaks at τ 9.7 = 11. However, at lower frac AGN , the peak shifts to lower τ 9.7 . Therefore, the τ 9.7 parameter is degenerate with frac AGN , i.e. the probability distribution of τ 9.7 depends on the value of frac AGN . This degeneracy makes the marginalized τ 9.7 PDF relatively flat, leading to the large uncertainty of this parameter. From the density map of J175520.20+660949.1, frac AGN is well constrained at a low level, but the τ 9.7 PDF is flat. This is because when frac AGN is low, the observed MIR SED is dominated by the galaxy component (Fig. 14 bottom). In this case, the models with different τ 9.7 are degenerate, in the sense that they have similar MIR SED shapes due to dominant galaxy contributions. The results above indicate that model degeneracy is responsible for the large uncertainties in the τ 9.7 estimation. X-ray data. For this type 2 source, the best-fit intrinsic AGN L 2500Å is significantly different in the two cases as labeled. The observed UV-to-IR fluxes are dominated by the galaxy component. Therefore, the intrinsic AGN power cannot be effectively constrained without the X-ray data.

FITTING WITH X-RAY UPPER LIMITS

In the previous section, we apply X-CIGALE to sources with X-ray detections. The X-ray sources are often only a small fraction ( 10%) of the entire sample in optical/IR surveys. However, many extragalactic studies need to constrain AGN emission for the X-ray undetected majority galaxy population (e.g. Buat et al. 2015;Vito et al. 2016;Bowman et al. 2019). This task can be achieved by X-CIGALE, based on X-ray flux upper limits. The detailed fitting algorithm is presented in §4.3 of Boquien et al. (2019), and we do not repeat it here. Below, we test this usage with 100 randomly selected AKARI-NEP galaxies (Buat et al. 2015), all of which are X-ray undetected. Here, we choose AKARI-NEP rather than SDSS or COSMOS, since AKARI-NEP has the best multiwavelength coverage among the three surveys ( §3). We focus on a relatively small sample (100), because the upper-limit analysis in X-CIGALE is time-consuming due to its complicated mathematical form of χ 2 (see §4. vary from source to source. The derivation of the precise values requires intensive simulations (e.g. Xue et al. 2016;Luo et al. 2017), which are beyond the scope of this work. Instead, we adopt a single conservative value, 5.3×10 -15 erg s -1 cm -2 (4.4×10 -7 mJy, 2-7 keV), for all the sources. For ≈ 80% of the survey area, the actual Chandra flux limits should be lower than this value (see Fig. 10 of Krumpe et al. 2015). We run X-CIGALE based on this upper limit. The parameter settings are the same as in §3.3.1 except that we allow frac AGN = 0, since it is possible that the AGN component does not exist for these upper-limit sources. For comparison, we also re-run X-CIGALE without the X-ray upper limit. Fig. 19 compares the Bayesian-like estimate of frac AGN of these two runs (frac AGN,noX vs. frac AGN,Xup ). As expected, frac AGN,Xup is generally lower than frac AGN,noX , since the X-ray upper limit can constrain AGN power. Notably, for ≈ 10% of the sources, the frac AGN,Xup is much lower (∆frac AGN,Xup > 0.2) than frac AGN,noX . This result indicates that, even when excellent MIR data are present ( §3.3.1), the AGN-galaxy decomposition might still be inaccurate/ambiguous without X-ray data. On the other hand, there are still ≈ 30% of sources that have non-negligible frac AGN (> 0.1) when the X-ray upper limit is used in the fitting. It is possible that a non-negligible IR flux is contributed by the AGN. However, another possibility is that the current X-ray upper limit is too high to effectively constrain AGN emission. In the future, Athena may clarify this problem with its great sensitivity. The total Chandra exposure time on the AKARI-NEP field is ≈ 300 ks. Given this amount of exposure time, Athena can reach its confusion-limited sensitivity of ∼ 1 × 10 -16 erg s -1 cm -2 (8 × 10 -9 mJy, 2-7 keV) for the entire AKARI-NEP field. 8 Assuming this X-ray flux limit, we re-run X-CIGALE. The resulting frac AGN,Xup is below 0.01 for all sources (see Fig. 19), indicating that the AGN SED contribution will be negligible if a source is undetected by Athena. Therefore, X-CIGALE, with future Athena observations, will have great power in unambiguously determining the presence of AGN. This feature will be extremely helpful for future extragalactic studies.

8 Chandra can, in principal, reach deeper sensitivity than Athena thanks to its superior angular resolution. However, reaching Athena-like (or deeper) flux limits will practically need large amounts of exposure time of Chandra. This has only been achieved in two small Chandra fields (only ∼ 500 arcmin 2 each), i.e. 7 Ms CDF-S and 2 Ms CDF-N (Xue et al. 2016;Luo et al. 2017).

SUMMARY AND FUTURE PROSPECTS

We have developed and tested X-CIGALE, a new version of the galaxy SED fitting code, CIGALE. Our development and test results are summarized below.

(i) We have developed a new X-ray module ( §2.2). The module is mainly designed to connect the intrinsic X-ray emission with other wavelengths, and X-ray obscuration and transmission should be corrected before providing the X-ray data to X-CIGALE. The X-ray module includes the X-ray emission from both galaxy and AGN. The galaxy component includes the emission from HMXB, LMXB, and hot gas. The AGN's X-ray SED is connected to its UV-to-IR SED using the well-known α ox -L 2500Å relation.

(ii) We have implemented a modern torus model, SKIRTOR, to fit AGN UV-to-IR SEDs ( §2.3). SKIRTOR adopts a clumpy two-phase torus, which is responsible for obscuring the UV/optical emission from the AGN disk. SKIRTOR is developed from a 3D radiative-transfer method, and thus obeys the energy-conservation law. However, SKIRTOR assumes that the AGN disk emission is absolutely unextincted when viewed from the polar direction. Therefore, SKIRTOR cannot model the SEDs of slightly extincted type 1 AGNs. To overcome this disadvantage, we introduce extinction from polar dust ( §2.4). The extinction amplitude, E(B-V), is a free model parameter set by the X-CIGALE user.

(iii) We have tested X-CIGALE on the AGNs with X-ray detections in SDSS, COSMOS, and AKARI-NEP §3. The three samples have distinctive characteristics in terms of AGN properties and available data. The fitting quality is good in general, with typical χ 2 ∼ 1 for all the samples. This result indicates that X-CIGALE is capable in modelling observed AGN SEDs under different circumstances. We also compare the fittings results with vs. without X-ray data. We find that the resulting AGN power is sometimes different in the two cases, when both AGN and galaxy components are present. Therefore, the AGN-galaxy SED decomposition may be unphysical without the constraints from X-ray data. We discuss constrainability and degeneracy of model parameters in the fitting of AKARI-NEP, for which excellent mid-IR photometric coverage is available.

(iv) We also test X-CIGALE on a random sample of galaxies with only Chandra X-ray upper limits in the AKARI-NEP field, where excellent MIR data are available ( §4). We compare the fitting results with and without the X-ray upper limits. After using the X-ray upper limits, frac AGN sometimes becomes lower, indicating that the current Chandra upper limit can effectively constrain AGN emission, as least for some systems. We also evaluate the potential of the future Athena mission by replacing the Chandra upper limit with the expected Athena value for a similar exposure time (≈ 300 ks). The resulting frac AGN is constrained to a negligible level (< 1%) for all the sources, indicating that Athena can robustly constrain AGN emission in general with a moderate amount of exposure time ( 300 ks).

We publicly release X-CIGALE on the official website of CIGALE. 9 As for the previous versions of CIGALE, X-CIGALE is open-source, allowing the user to modify the source code freely. In the future, we will further develop X-CIGALE and enable it to address special AGNs such as radio-loud and BAL objects ( §2.2). Besides the three surveys tested in this work (SDSS, COSMOS, and AKARI-NEP), the user can apply X-CIGALE to the existing multiwavelength surveys such as CDF-S (Luo et al. 2017) The intrinsic AGN L ν at 2500 Å at viewing angle = 30 • W Hz -1 agn.accretion_power

The intrinsic AGN disk luminosity averaged over all directions W

Figure 1 .

 1 Figure 1. An example X-ray SED model for a typical source with AGN L 2-10 keV = 10 43 erg s -1 , M = 10 11 M , SFR = 10 M yr -1 , T = 1 Gyr, and Z = 0.02. Different colors indicate different components. For this source, the X-ray luminosity is dominantly contributed by the AGN.

Figure 2 .

 2 Figure 2. The model SEDs for an unobscured AGN with log L 2500Å = 30 (cgs). Different colors indicate different α ox .The "breaks" at 5 nm are caused by the wavelength limit of the X-ray module, which ends at 5 nm. As explained in §2.2.3, such breaks are generally not problematic for practical purposes, as most X-ray instruments only cover wavelengths 10 0.4 nm ( 0.5 keV, as marked).

Figure 4 .

 4 Figure 4. Schematic plot (at the meridional plane) of the AGN obscuration system adopted by X-CIGALE (not to scale). The model in original CIGALE only includes torus obscuration. We add the obscuration of polar dust to account for type 1 AGN extinction.

Figure 5 .

 5 Figure 5. Same format as in Fig. 2, but for AGNs with different polar-dust E(B -V ). The top and bottom panels are for type 1 and type 2 AGNs, respectively. The SEDs are normalized at 10 µm. Except for E(B -V ), the other model parameters are the same (e.g. polar-dust temperature = 100 K and emissivity = 1.6) when generating these model SEDs.

Figure 6 .Figure 7 .

 67 Figure 6. A random example of fitted SED from SDSS. The fitting residuals (observed flux/model flux) are shown at the bottom. Some source properties are labelled on the plot. Since the fitting only uses the AGN component, the total SED is actually the same as the AGN SED.

Figure 9 .

 9 Figure 9. Left: Estimated vs. true values in the mock analyses of the SDSS sample for polar-dust E(B -V ). The black line indicates the 1:1 relation between estimated and true values. For clarity, the data points are from 100 randomly selected sources in the SDSS sample. The values on the x-axis are shifted slightly for display purposes only. The residuals (defined as estimated minus true) are displayed on the bottom panel, with the black line indicating a zero residual. For E(B -V ), the estimated and true values are well correlated, indicating that E(B -V ) model parameter can be effectively constrained by the data. Same format as left but from the fitting with free viewing angle.

Figure 11 .

 11 Figure 11. Same format as in Fig. 6 but for COSMOS type 1 (top) and type 2 (bottom). The orange and blue curves indicate galaxy and AGN model components.

Figure 12 .

 12 Figure12. Same format as Fig.8but for the COSMOS sample. The blue and orange colors indicate type 1 and type 2 sources, respectively. The differences between L 2500Å,X and L 2500Å,noX are larger compared to those in SDSS (Fig.8).

Figure 13 .

 13 Figure 13.An example COSMOS SED fitted with (top) vs. without (bottom) X-ray data. For this type 2 source, the best-fit intrinsic AGN L 2500Å is significantly different in the two cases as labeled. The observed UV-to-IR fluxes are dominated by the galaxy component. Therefore, the intrinsic AGN power cannot be effectively constrained without the X-ray data.

Figure 14 .

 14 Figure 14. Same format as in Fig. 11 but for two sources in AKARI-NEP. Notably, the MIR data are well fitted with our model.

Figure 15 .

 15 Figure15. Same format as Fig.8but for the AKARI-NEP sample. The differences between L 2500Å,X and L 2500Å,noX are larger compared to those in SDSS (Fig.8).

Figure 16 .Figure 17 .

 1617 Figure 16. Same format as in Fig. 9 but for frac AGN (left) and E(B -V ) (right) of the AKARI-NEP AGNs.

Figure 18 .

 18 Figure18. Top: the τ 9.7 PDF for the example source in Fig.14. Bottom: the 2D probability density map of τ 9.7 vs. frac AGN . The density map is normalized such that the 2D integral is unity.

  1 X-CIGALE then calculates |∆α ox | and discards the models with |∆α ox | > |∆α ox | max . In our test fitting ( §3 and §4), we adopt |∆α ox | max = 0.2, corresponding to the ≈ 2σ scatter of the α ox -L 2500Å relation

Table 1 .

 1 Sample properties

	Name	N	m r	f 2-10 keV	Redshift Type AGN
	SDSS	1986 19.2-20.9 2.2-10.5	0.6-1.9	1
	COSMOS	590	21.3-23.7	0.4-1.7	0.6-1.8	1 & 2
	AKARI-NEP	74	20.5-23.4	0.6-2.1	0.5-1.4	1 & 2
	N . -(1) Survey name. (2) Number of AGNs.

Table 2 .

 2 SDSS fitting parameters

	Module	Parameter	Values
	AGN (UV-to-IR):
	SKIRTOR

Table 3 .

 3 COSMOS and AKARI-NEP Fitting Parameters

	Module	Parameter	Values
	Star formation history:	e-folding time, τ (Gyr)	0.1, 0.5, 1, 5
	delayed model, SFR ∝ t exp(-t/τ)	Stellar Age, t (Gyr)	0.5, 1, 3, 5, 7
	Simple stellar population:	Initial mass function	Chabrier (2003)
	Bruzual & Charlot (2003)	Metallicity (Z)	0.02
	Galactic dust attenuation:
	Calzetti et al. (2000) & Leitherer et al. (2002)

Table A1 .

 A1 Output SED components for the new X-ray and SKIRTOR modules In X-CIGALE output, all the SED components are in the format of L λ in units of W nm -1 .

	Module	Component	Explanation
	X-ray ( §2.2.2)	xray.agn	The AGN corona
	10 -3 -5 nm	xray.galaxy	The total SED of HMXB, LMXB, and hot gas
	SKIRTOR ( §2.3)	agn.SKIRTOR2016_disk The AGN disk
	8-10 6 nm	agn.SKIRTOR2016_dust	The dust reemission
	N . -

Table A2 .

 A2 Additional output physical parameters for the new X-ray and SKIRTOR modules

	Module	Parameters	Explanation	Units
		xray.agn_Lnu_2keV	The AGN L ν at 2 keV	W Hz -1
		xray.agn_Lx_2to10keV	The AGN 2-10 keV luminosity	W
		xray.agn_Lx_total	The AGN total (0.25-1200 keV) X-ray luminosity	W
	X-ray	xray.alpha_ox	The AGN α ox	-
		xray.lmxb_Lx_2to10keV	The 2-10 keV LMXB luminosity	W
		xray.hmxb_Lx_2to10keV	The 2-10 keV HMXB luminosity	W
		xray.hotgas_Lx_0p5to2keV The 0.5-2 keV hot-gas luminosity	W
		agn.disk_luminosity	The observed AGN disk luminosity (might be extincted)	W
		agn.dust_luminosity	The observed AGN dust reemitted luminosity	W
	SKIRTOR	agn.luminosity	The sum of agn.disk_luminosity and agn.dust_luminosity	W
		agn.intrin_Lnu_2500A

© 2018 The Authors

MNRAS 000,1-17 (2018)

This α ox range corresponds to

2-10 keV X-ray bolometric corrections ranging from ≈ 10 to ≈ 500.

http://www.skirt.ugent.be/root/_landing.html

https://sites.google.com/site/skirtorus/sed-library

This statement breaks if the AGN host galaxies are highly star-forming in general. However, the AGN hosts tend to have normal levels of starformation activity, as shown by previous studies (e.g.Harrison et al. 2012;Stanley et al. 2015).MNRAS 000,1-17 (2018)

Due to a technical reason, this value cannot equal to 1.MNRAS 000,1-17 (2018)

https://cigale.lam.fr/ MNRAS 000,1-17 (2018)

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000,1-17 (2018)

ACKNOWLEDGMENTS

We thank the referee for helpful feedback that improved this work. GY acknowledges support from the French Space Agency CNES. VB received funding from Excellence Initiative of Aix-Marseille University -AMIDEX, a French "Investissements d'Avenir" programme. MS acknowledges support by the Ministry of Education, Science, and Technological Development of the Republic of Serbia through the projects Astrophysical Spectroscopy of Extragalactic Objects (176001) and Gravitation and the Large Scale Structure of the Universe (176003). WNB acknowledges support from NASA ADP grant 80NSSC18K0878 and the V.M. Willaman Endowment. This project uses Astropy (a Python package; see Astropy Collaboration et al. 2018).

with X-ray upper limits for 100 random galaxies in the AKARI-NEP field. The blue and orange points represent frac AGN,Xup obtained from Chandra and Athena upper limits, respectively. The solid black lines indicate the 1:1 relation. As expected, frac AGN,Xup is generally lower than frac AGN,noX , since the X-ray upper limit can constrain AGN power. The current Chandra upper limit can effectively constrain AGN power for some sources. The Athena upper limit strongly suppresses frac AGN,Xup to < 1% for all sources.

APPENDIX A: NEW INPUTS AND OUTPUTS IN X-CIGALE

The input parameters for the new X-ray and SKIRTOR modules can be found in Table 3. After fitting, X-CIGALE can output the best-fit model SEDs of different components, and the SED components for the new X-ray and SKIRTOR modules are summarized in Table A1. Besides the best-fit SEDs, X-CIGALE can also yield the maximum-likelihood and Bayesian-like values of source physical properties. These physical properties include not only the model parameters (Table 3), but also some additional quantities as listed in Table A2. New quantities can be added in the future as requested by the user of X-CIGALE. We remind that the quantities, "agn.intrin_Lnu_2500A" and "xray.alpha_ox" refers to the values as measured at a viewing angle of 30 • (see §2.2.3). The quantity "agn.accretion_power" refers to the intrinsic (unextincted) AGN disk luminosity averaged over all directions (weighted by sin θ; see §2.2.3). This quantity, paired with an assumed radiative efficiency, can be used to estimate BH accretion rate (e.g. Yang et al. 2017Yang et al. , 2019)).