Constraining the Active Galactic Nucleus and Starburst Properties of the IR-luminous Quasar Host Galaxy APM 08279+5255 at Redshift 4 with SOFIA

TL;DR

This study uses SOFIA and Spitzer data to analyze the IR properties of a high-redshift quasar-host galaxy, revealing detailed dust temperatures, star formation rates, and the quasar's contribution to IR luminosity with advanced modeling.

Contribution

It introduces a semi-analytic IR spectral decomposition method that accounts for a clumpy torus, providing more accurate estimates of dust properties and IR luminosity in high-z quasars.

Findings

Cold dust dominates FIR energy budget (66%)

Star formation rate corrected for lensing is approximately 3075 Msun/yr

Quasar contributes 30% to FIR luminosity, dominating total IR emission

Abstract

We present far-IR photometry and infrared spectrum of the z=3.9114 quasar/starburst composite system APM 08279+5255 obtained using the Stratospheric Observatory for Infrared Astronomy (SOFIA)/HAWC+ and the Spitzer Space Telescope Infrared Spectrograph (IRS). We decompose the IR-to-radio spectral energy distribution (SED), sampled in 51 bands, using (i) a model comprised of two-temperature modified blackbodies (MBB) and radio power-laws and (ii) a semi-analytic model, which also accounts for emission from a clumpy torus. The latter is more realistic but requires a well-sampled SED, as possible here. In the former model, we find temperatures of T_warm = 296^17_15 K and T_cold = 110^3_3 K for the warm and cold dust components, respectively. This model suggests that the cold dust component dominates the FIR energy budget (66%) but contributes only 17% to the total IR luminosity. Based on the torus models, we infer an inclination angle of i=15^8_8 degree and the presence of silicate emission, in accordance with the Type-1 active galactic nucleus nature of APM 08279+5255. Accounting for the torus' contribution to the FIR luminosity, we find a lensing-corrected star formation rate of SFR=3075x(4/mu_L) Msun yr^-1. We find that the central quasar contributes 30% to the FIR luminosity but dominates the total IR luminosity (93%). The 30% correction is in contrast to the 90% reported in previous work. In addition, the IR luminosity inferred from the torus model is a factor of two higher. These differences highlight the importance of adopting physically motivated models to properly account for IR emission in high-z quasars, which is now possible with SOFIA/HAWC+.

Figures

Tables

Table 1. Table 1 Photometry data for APM 08279 + + 5255.

Wavelength/Band	Frequency	Flux Density	Instrument/Band	Reference
($\mathrm{\mu m}$)	(GHz)	(mJy)
1.03	291060	5.00 $\pm$ 0.16	HST/NICMOS	1
1.25	239834	7.19 $\pm$ 0.20	Keck I/NIRC	2
1.55	193414	7.53 $\pm$ 0.19	HST/NICMOS	1
1.65	181692	8.97 $\pm$ 0.25	Palomar/P200 InSb	2
1.90	157786	7.81 $\pm$ 0.96	HST/NICMOS	1
2.15	139438	9.40 $\pm$ 0.26	Palomar/P200 InSb	2
3.4	88174	11.01 $\pm$ 0.23	WISE/W1	3†
3.5	85650	27.3 $\pm$ 1.0	Palomar/P200 InSb	2
3.8	79433	9.50 $\pm$ 1.22	Subaru/IRCS	4
3.6	83275	10.74 $\pm$ 0.0023	Spitzer/IRAC	3†
4.5	66620	7.88 $\pm$ 0.0021	Spitzer/IRAC	3†
4.7	63994	7.20 $\pm$ 1.98	Subaru/IRCS	4
4.6	65172	10.86 $\pm$ 0.20	WISE/W2	3†
5.8	51688	15.49 $\pm$ 0.007	Spitzer/IRAC	3†
12	24983	62.92 $\pm$ 0.927	WISE/W3	3†
12	24983	101	IRAS	5
12.5	24000	74 $\pm$ 11	Keck II/MIRLIN	2
17.9	16750	103 $\pm$ 50	Keck II/MIRLIN	2
22	13626	187.6 $\pm$ 4.84	WISE/W4	3†
24	12491	186.8 $\pm$ 0.779	Spitzer/MIPS	3†
25	12000	226.1 $\pm$ 16.2	IRAS	5
53	5656	567 $\pm$ 71	SOFIA/HAWC$+$	3
60	5000	511.1 $\pm$ 51.1	IRAS	5
70	4283	654.1 $\pm$ 8.745	Herschel/PACS	3
89	3368	682 $\pm$ 50	SOFIA/HAWC$+$	3
100	3000	951.1 $\pm$ 228	IRAS	5
154	1947	696 $\pm$ 43	SOFIA/HAWC+	3
160	1874	758.8 $\pm$ 13.2	Herschel/PACS	3
250	1200	627.2 $\pm$ 8.1	Herschel/SPIRE	3
350	857	431.4 $\pm$ 6.3	Herschel/SPIRE	3
350	850	386 $\pm$ 32	CSO/SHARC II	6
450	660	342 $\pm$ 26	CSO/SHARC II	6
450	660	285 $\pm$ 11	JCMT/SCUBA	7
450	660	211 $\pm$ 47	JCMT/SCUBA	8
500	600	249.6 $\pm$ 8.1	Herschel/SPIRE	3
850	350	84 $\pm$ 3	JCMT/SCUBA	7
850	350	75 $\pm$ 4	JCMT/SCUBA	8
1194	251	34 $\pm$ 0.55	PdBI	9
1219	245.905	31.4 $\pm$ 2.0	PdBI	10
1266	236.797	26.6 $\pm$ 1.3	PdBI	10
1350	225	24 $\pm$ 2	JCMT/SCUBA	8
1400	214.0	17.0 $\pm$ 0.5	PdBI	11
1400	211.1	16.9 $\pm$ 2.5	PdBI	12
1490	201.166	16.5 $\pm$ 0.8	PdBI	10
1958	153.132	5.4 $\pm$ 0.3	PdBI	10
2710	110.7	2.17 $\pm$ 0.19	CARMA	13
2760	108.6	2.08 $\pm$ 0.22	CARMA	13
2840	105.4	2.04 $\pm$ 0.08	CARMA	13
3200	94.3	1.2 $\pm$ 0.3	PdBI	11
3200	93.879	1.3 $\pm$ 0.2	PdBI	12
3200	93.856	1.1 $\pm$ 0.02	NOEMA	14
3300	90.2	1.3 $\pm$ 0.2	PdBI	12
3300	90.2	0.66 $\pm$ 0.18	PdBI	15
7000	43.3	0.90	VLA	16
13000	23.4649	0.30 $\pm$ 0.05	VLA	16
13000	23.4649	0.41 $\pm$ 0.07	VLA	17
13000	23.4649	0.38 $\pm$ 0.02	VLA	18
20000	14.9399	0.30 $\pm$ 0.09	VLA	18
36000	8.4601	0.45 $\pm$ 0.02	VLA	18
36000	8.4601	0.45 $\pm$ 0.03	VLA	1
60000	4.5276	0.55 $\pm$ 0.04	VLA	18
200000	1.400	1.16 $\pm$ 0.03	VLA	18
200000	1.400	0.92 $\pm$ 0.13	VLA	19 aaVLA FIRST Catalog Database version 2014dec17.
900000	0.334	5.10 $\pm$ 0.80	VLA	20

Table 2. Table 2 SED Model Parameters for APM 08279 + + 5255.

Properties	Units	Values
Analytical: Two-temperature MBB + Radio Power-laws aaSee §4.1
$T_{d,\mathrm{warm}}$	(K)	296${}_{-15}{}^{}{}_{}{}^{+17}$
$T_{d,\mathrm{cold}}$	(K)	110${}_{-3}{}^{}{}_{}{}^{+3}$
$F_{\mathrm{Filling}}^{\mathrm{cold}}$		0.9${}_{-0.1}{}^{}{}_{}{}^{+0.1}$
$x_{\mathrm{warm}}$ bb$x\equiv M_d/r_0^2$, see Equation 4.	(107 $M_{\odot }$  pc-2)	2.5${}_{-1.0}{}^{}{}_{}{}^{+1.2}$
$x_{\mathrm{cold}}$ bb$x\equiv M_d/r_0^2$, see Equation 4.	(109 $M_{\odot }$  pc-2)	1.4${}_{-0.1}{}^{}{}_{}{}^{+0.1}$
$L_{\mathrm{IR}}$	(10${}_{}{}^{15}\mu$-1  $L_{\odot }$)	1.8${}_{-0.6}{}^{}{}_{}{}^{+0.7}$
$L_{\mathrm{FIR}}$	(10${}_{}{}^{14}\mu$-1  $L_{\odot }$)	1.8${}_{-0.2}{}^{}{}_{}{}^{+0.1}$
Semi-analytical: Including Clumpy Torus ccSee §4.2
$r_{\mathrm{hot}}$	(pc)	655${}_{-409}{}^{}{}_{}{}^{+418}$
$M_d^{\mathrm{hot}}$	($\mu$-1  $M_{\odot }$)	11${}_{-1}{}^{}{}_{}{}^{+5}$
$M_d^{\mathrm{cold}}$	(10${}_{}{}^9\mu$-1  $M_{\odot }$)	1.4${}_{-0.2}{}^{}{}_{}{}^{+0.3}$
$T_d^{\mathrm{hot}}$	(K)	1385${}_{-49}{}^{}{}_{}{}^{+10}$
$T_d^{\mathrm{cold}}$	(K)	100${}_{-1}{}^{}{}_{}{}^{+2}$
$N_0$		8.3${}_{-2.2}{}^{}{}_{}{}^{+1.3}$
$a$		$-$1.0${}_{-0.1}{}^{}{}_{}{}^{+0.1}$
$\theta$	($°$)	55${}_{-11}{}^{}{}_{}{}^{+4}$
${\tau }_V$		42${}_{-5}{}^{}{}_{}{}^{+5}$
$i$	($°$)	15${}_{-8}{}^{}{}_{}{}^{+8}$
$\log \left(L_{\mathrm{scale}}/L_{\odot }\right)$		48.8${}_{-0.1}{}^{}{}_{}{}^{+0.1}$
$L_{\mathrm{IR}}$	(10${}_{}{}^{15}\mu$-1  $L_{\odot }$)	3.5${}_{-0.1}{}^{}{}_{}{}^{+0.1}$
$L_{\mathrm{NIR}}$	(10${}_{}{}^{14}\mu$-1  $L_{\odot }$)	5.9${}_{-0.6}{}^{}{}_{}{}^{+0.4}$
$L_{\mathrm{MIR}}$	(10${}_{}{}^{14}\mu$-1  $L_{\odot }$)	2.7${}_{-0.7}{}^{}{}_{}{}^{+0.3}$
$L_{\mathrm{FIR}}$	(10${}_{}{}^{14}\mu$-1  $L_{\odot }$)	1.7${}_{-0.1}{}^{}{}_{}{}^{+0.5}$
$L_{\mathrm{FIR}}^{\mathrm{Torus}}$	(10${}_{}{}^{13}\mu$-1  $L_{\odot }$)	4.7${}_{-0.5}{}^{}{}_{}{}^{+0.3}$
$L_{\mathrm{FIR}}^{\mathrm{Cold}}$	(10${}_{}{}^{14}\mu$-1  $L_{\odot }$)	1.2${}_{-0.7}{}^{}{}_{}{}^{+0.5}$

Equations

$${\kappa }_{\nu }={\kappa }_0(\nu /{\nu }_0{)}^{\beta },$$

$$S_{\nu }(\nu )=S_{\nu ,\mathrm{d}}(\nu )+\sum _{i=\mathrm{flat},\mathrm{steep}}{c}_i\times {\nu }^{{\alpha }_i},$$

$$S_{\nu }=[B_{\nu }(T_d)-B_{\nu }(T_{\mathrm{CMB}})]\frac{(1-{\exp }^{-{\tau }_v})}{(1+z{)}^3}{\Omega }_{\mathrm{app}},$$

$${\tau }_{\nu }=\frac{{\kappa }_{\nu }{M}_{\mathrm{dust}}}{(\pi r_0^2)F_{\mathrm{filling}}}.$$

$$N_{\mathrm{tot}}\propto \int N_0\frac{\eta (r,z,\varphi )}{R_{\mathrm{cl}}^2}dV,$$

$$\left(\frac{R_{\mathrm{cl}}}{0.5\text{pc}}\right)\propto {\left(\frac{L_{\mathrm{scale}}}{10^{46}{\text{erg s}}^{-1}}\right)}^{\frac{1}{2}}{\left(\frac{T_{\mathrm{sub}}}{1500\text{K}}\right)}^{-2.8},$$

Full text

††slugcomment: Accepted to ApJ55affiliationtext: NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

1. Introduction

Active galactic nucleus (AGN)-starburst galaxies at high redshifts are powerful probes of the early evolution of black holes (BH) and galaxy assembly, and they may be the progenitors of present-day massive spheroidal galaxies. Circumstantial evidence such as the tight BH mass-host properties relations observed in the nearby universe (e.g., $M_{\rm BH}-\sigma_{*}$ and $M_{\rm BH}-M_{\rm bulge}$) and the similar trends found between the cosmic star formation rate (SFR) and BH accretion rate histories suggest that the growth of supermassive BHs (SMBHs) and their host galaxies are tightly linked (and possibly causally linked; see reviews by Alexander & Hickox 2012; Heckman & Best 2014; Madau & Dickinson 2014, and, e.g., Sparre et al. 2015 for evidence from simulations). The leading SMBH-host co-evolution picture suggests that galaxies evolve via hierarchical galaxy mergers, in which massive galaxies evolve from a starburst-dominated phase to a quasar-dominated phase and finally settle as ‘red and dead’ spheroidal galaxies (e.g., Sanders et al. 1988; Di Matteo et al. 2005; Hopkins & Beacom 2006; Hopkins et al. 2008; Treister et al. 2010; Narayanan et al. 2010; Page et al. 2012). Recent studies report that many of the most luminous AGNs at high redshift appear to be offset from the local BH-host relations, where high-$z$ AGNs appear to have more massive black holes than expected given their host galaxy properties (e.g., Walter et al. 2004; Riechers et al. 2008a, b; Merloni et al. 2010; Wang et al. 2013). This offset suggests that AGNs and their host galaxies may not grow “in tandem”, especially during their peak phases of activity ($z$  = 1$-$4). Currently, we still lack a clear picture relating the interplay between AGN activity and star formation in their host galaxies at early cosmic epochs, which is essential for understanding their formation and evolution. In theory, useful insights can be gained by determining how the far-infrared (FIR) luminosity (which depends on the AGN accretion rate and SFR), the interstellar medium (ISM) conditions, and morphology of AGN host galaxies differ between the quasar phase ($z$ $\sim$ 1 $-$ 4) and the present-day epoch and between the different phases of the starburst-quasar evolutionary sequence in the proposed SMBH-host co-evolution picture.

Previous studies indicate that $\sim$30% of high-$z$ quasars exhibit strong emission at (sub-)millimeter wavelengths (e.g., Chapman et al. 2005; Alexander et al. 2008; Coppin et al. 2008), indicating large quantities of dust contributing to IR luminosities of $>10^{\rm 13}$ $L_{\odot}$. The large IR luminosities in these quasars are commonly attributed to cold dust heated by their star formation, corresponding to SFRs of $\gtrsim$1000 $M_{\odot}$  yr*-1* (e.g., Wang et al. 2008). One of the key challenges in improving our understanding of the AGN-host co-evolution picture stems from the limited number of comprehensive studies that have measurements covering the rest-frame mid-IR parts of the dust spectral energy distributions (SEDs) of high-$z$ AGN-starburst systems, with only a few detections attained with The Infrared Astronomical Satellite (IRAS), the Spitzer Space Telescope, and the Herschel Space Observatory (e.g., Rowan-Robinson et al. 1991, 1993; Barvainis et al. 1995; Polletta et al. 2008). Thus, the fraction of the cold dust heated by young stars (and thus the SFR) in high-$z$ quasar hosts remains uncertain (e.g., Omont et al. 2003; Carilli et al. 2005; Netzer et al. 2014). By probing the hot and warm dust heated by the AGN in high-$z$ quasars, one can better understand the co-eval growth of SMBHs and their host galaxies in their early stages of evolution. In addition, determining the geometry of the dusty torus together with the host properties, one can investigate whether dust-obscured and unobscured AGNs in the early Universe represent quasars at different stages of evolution (see the review by Netzer 2015) or whether they are simply a manifestation of different viewing geometries under the well-known AGN unification framework (i.e., Type-1 versus Type-2 AGN; e.g., Antonucci 1993; Urry & Padovani 1995).

In this paper, we demonstrate the new capabilities of the Stratospheric Observatory for Infrared Astronomy (SOFIA) enabled by the new sensitive FIR instrument, the High-resolution Airborne Wideband Camera$+$ (HAWC$+$; Dowell et al. 2013; Harper et al. 2018), together with ancillary IR photometry, to constrain the AGN contribution to the IR luminosity of APM 08279$+$5255. The target is among the most apparently IR-luminous, gas-rich, strongly lensed AGN-starburst systems known to date, at a redshift of $z$ $\sim$ 3.91 (RA, Dec J2000 = 08h31m41$\fs$7, $+$52$\arcdeg$45′17$\farcs$5; Egami et al. 2000; Weiß et al. 2007). Existing studies suggest that APM 08279$+$5255 is a rare starbursting quasar, in which both the AGN and host galaxy are undergoing extremely active phases, with an SFR of 5000 $\times$ $(4/\mu_{L})$ and an Eddington ratio of $\lambda_{\rm edd}$ = 0.4 (Saturni et al. 2018). Since APM 08279$+$5255 is magnified by $\mu_{\rm L}$ $\sim$ 4$-$100 (Egami et al. 2000; Riechers et al. 2009), it has an apparent bolometric luminosity of $L_{\rm bol}$ = 7$\times 10^{15}$ $L_{\odot}$ (Irwin et al. 1998) and IR luminosity of $L_{\rm IR}$ $=$ 1015 $\mu_{L}$-1 $L_{\odot}$, making it one of the most ideal high-$z$ galaxies for a demonstration study for SOFIA/HAWC$+$. Due to its extreme brightness, it is one of the best-studied high-$z$ objects and thus is currently one of the rare cases with a rich set of photometry. This enables more realistic SED modeling for this source for the first time, motivating our work presented here. Applying complex and realistic models to $z$ $>$ 0 AGNs has been done previously, but mainly out to $z$ $\lesssim$ 2 (Williams et al. 2018; Zhuang et al. 2018). While Leipski et al. (2014) and Lyu et al. (2016) use similar approaches to decompose the SEDs of $z$ $\sim$ 6 AGNs, these sources have less complete SED coverage than APM 08279$+$5255. In most cases, their IR SEDs are constrained by less than ten data points. By combining new SOFIA/HAWC+ photometry and a Spitzer/Infrared Spectrograph (IRS) spectrum with ancillary data covering multiple IR bands (Table 2.3), we improve constraints on the IR SED of APM 08279$+$5255 relative to previous work and model its SED with different AGN torus models to infer its torus properties.

This paper is structured as follows. In §2, we summarize the observations and procedures used to reduce the data. We also report ancillary data used in our analysis. In §3, we present the results. In §4, we detail the SED modeling analysis. Finally, we discuss the results and implication of our findings in §5, and we summarize the main results and present our conclusions in §6. Throughout this paper, we use a concordance $\Lambda$CDM cosmology with parameters from the WMAP9 results: $H_{0}$ = $69.32$ km s*-1* Mpc*-1*, $\Omega_{\rm M}$ = $0.29$, and $\Omega_{\Lambda}$ = $0.71$ (Hinshaw et al. 2013).

2. New Observations and Ancillary Data

2.1. New Data: Spitzer/IRS

Observations of APM 08279$+$5255 were carried out with the Infrared Spectrograph (IRS) onboard the Spitzer Space Telescope on 2008 December 4 for a total observing time of 4.5 hr (Program ID: 50784; PI: Riechers). The first order of the long wavelength, low resolution module (LL1: 19.5$-$38.0 $\micron$) was used, covering the rest-frame 6.7 $\micron$ polycyclic aromatic hydrocarbons (PAH) feature. These observations were obtained in spectral mapping mode, wherein the target was placed at six positions along the slit. The slit was positioned relative to a reference star using “high accuracy blue peak-up”111https://irsa.ipac.caltech.edu/data/SPITZER/docs/irs/irsinstrumenthandbook/85/.

The Basic Calibrated Data (BCD) files were produced by the Spitzer pipeline, which includes ramp fitting, dark sky subtraction, droop correction, linearity correction, flat-fielding, and wavelength calibration. Before spectral extraction, the two-dimensional dispersed frames were corrected for “rogue” (unstable) pixels using irsclean, latent charge removal and residual sky subtraction222https://irsa.ipac.caltech.edu/data/SPITZER/docs/dataanalysistools/cookbook/23/. Sky frames at each map position is created by combining the corrected frames separately. Uncertainty frames associated with the combined frames were produced based on the variance of each pixel as a function of time in the sky frames. One-dimensional spectra were extracted at the position of the continuum emission using the (spice) software. A final spectrum was obtained by averaging over all map positions. Based on the residual sky background, the median noise of the final spectrum is 5.0 mJy per $\Delta\lambda$ = 0.169 $\micron$ wavelength bin (observed-frame).

2.2. New Data: SOFIA/HAWC$+$

Observations of APM 08279$+$5255 with SOFIA/HAWC$+$ were performed on 2017 October 18 at an altitude of $\sim$43,000 ft during cycle 5 (Program ID: 05_0165; PI: Riechers). A single pointing was used to observe at 53 (Band A), 89 (Band C), and 154 $\micron$ (Band D). Total Intensity observing in On-the-fly Mapping mode was used, producing a continuous telescope motion with 1.5 min Lissajous scans, covering (regions slightly smaller than) the standard HAWC$+$ fields-of-view (FoVs) of 2$\farcm$7 $\times$1$\farcm$7 (Band A), 4$\farcm$2 $\times$2$\farcm$6 (Band C), and 7$\farcm$3 $\times$4$\farcm$5 (Band D). Five scans corresponding to a total integration time of 333 s were obtained in Band A, 12 scans totaling 870 s in Band C, and 12 scans totaling to 936 s in Band D. The scans were centered in the middle of the R0/T0 subarrays. As such, the maps are made from R0/T0 only.

The observations were performed at elevation angles between 30$-$37$\arcdeg$. No direct measure of Precipitable Water Vapor (PWV) was used in the calibration. Instead, the HAWC$+$ pipeline uses a standard atmospheric model (the ATRAN model; Lord 1992) to compute the predicted water vapor overburden at the zenith. The assumed zenith opacities are $\tau_{z}$ = 0.11, 0.21, 0.24 for Bands A, C, and D, respectively.

The data were reduced and processed using the crush333http://www.sigmyne.com/crush/ software, in which an atmospheric correction based on the ATRAN model was applied. For flux calibration, we applied a fixed instrumental conversion factor to convert readout counts to janskys per pixel. The conversion factor was determined by comparing the readout counts to the true fluxes of the primary and secondary flux calibrators, Uranus and Neptune, using the Bendo et al. (2013) temperature model. In the post-processed images, the PSF FWHMs are 7$\farcs$1, 11$\farcs$0, and 19$\farcs$5 at 53, 89, and 154 $\micron$, respectively. After reduction, the pixel scales are 1$\farcs$00, 1$\farcs$55, and 2$\farcs$75 for Bands A, C, and D, respectively.

2.3. Ancillary Data

In the analysis, we include photometry published in the literature, as listed in Table 2.3. These include photometry covering observed-frame NIR-to-radio wavebands.

We also report unpublished Spitzer Space Telescope/IRAC and MIPS photometry based on the Spitzer Enhanced Imaging Products (SEIP), which are available through the NASA/IPAC Infrared Science Archive (IRSA). The photometry in all IRAC bands reported in Table 2.3 was extracted using a 3$\farcs$8 diameter aperture, with aperture corrections applied. The MIPS photometry at 24 $\micron$ was extracted through PSF-fitting to the high-S/N detection. We retrieved the Spitzer images from the Spitzer Heritage Archive.

We also report Herschel/PACS photometry, which are tabulated in the Herschel/PACS Point source catalog (HPPSC; Marton et al. 2017). PACS photometric maps were generated using the JScanam task within the Herschel Interactive Processing Environment (HIPE version 13.0.0; Ott 2010). Sources were identified using sussextractor, and the flux densities were obtained by performing aperture photometry.

Extraction of the Herschel/SPIRE photometry at 250, 350, and 500 $\micron$ (Table 2.3) was performed using sussextractor within HIPE version 15.0.1 on Level 2 SPIRE maps obtained from the Herschel Science Archive. These maps were reduced and processed using the SPIRE pipeline version SPGv14.1.0 within HIPE, with calibration tree SPIRE_CAL_14_3. The sussextractor task estimates the flux density from an image convolved with a kernel derived from the SPIRE beam. The flux densities measured by sussextractor were confirmed using the Timeline Fitter, which performs photometry by fitting a 2D elliptical Gaussian to the Level 1 data at the source position given by the output of sussextractor. The fluxes obtained from both methods are consistent within the uncertainties.

3. Results

We detect emission in all three bands covered by SOFIA/HAWC$+$ at high significance (Figure 1). Point-source photometry were extracted using PSF-fitting and are listed in Table 2.3. Fitting a power-law to the IRS spectrum yields $S_{\nu}\propto\lambda^{0.01}$.

The PAH feature at 6.2 $\micron$ covered by our IRS observations remains undetected down to a $3\sigma$ limit of $<$3.4$\times 10^{-19}$ W cm*-2* (Figure 2). This limit was calculated by considering the residual flux density after removing the continuum and assuming a linewidth covering rest-frame wavelengths from 5.9$-$6.5 $\micron$. The continuum level was calculated by fitting a power-law over rest-frame wavelengths of 5.5$-$6.7 $\micron$.

4. Analysis: SED Modeling

4.1. Analytic: Two-temperature MBB + Radio Power-laws

The IR-to-radio SED of APM 08279$+$5255 has been modeled previously by, e.g., Rowan-Robinson (2000), Beelen et al. (2006), Weiß et al. (2007), and Riechers et al. (2009); some of these models have their dust-emitting sizes ($r_{0}$) constrained through large velocity gradient (LVG) modeling of the multi-$J$ CO and HCN line ratios and spatially resolved CO imaging. Existing models show that the observed dust SED can only be adequately described by a model with at least two dust components — a warmer and a cooler dust component.

Weiß et al. (2007) adopted a model assuming that the overall size of the dust-emitting region is similar to that of the CO line emission, with relative area filling factors of $F_{\rm filling}$ = 0.75 and 0.25 for the warmer and cooler components, respectively, and a dust absorption coefficient of

[TABLE]

where $\kappa_{0}$ = 0.4 cm2/g, $\nu_{0}$ = 250 GHz, and the emissivity spectral index $\beta$ is fixed to $2$. They find that the warmer component contributes $\sim$85% to the FIR luminosity ($L_{\rm FIR}$  = 2 $\times$ 1014 $\mu_{L}$-1  $L_{\odot}$) and suggest that this component is heated by the central AGN. However, in these previous models, the parameters describing the warm dust peak relied heavily on the constraints imposed by the IRAS measurement at $\sim$100 $\micron$, which has a $\sim$25% uncertainty due to the limited S/N and the very large IRAS beam and pixel sizes. Here, we update the SED by including our recently obtained SOFIA/HAWC$+$ data and additional photometry (e.g., Herschel) obtained after the Riechers et al. (2009) study, which extended the Weiß et al. (2007) study by also modeling the radio SED. We note an outlier444This outlier is unlikely to result from contamination due to bright rest-frame optical lines, such as H$\alpha$, [N ii], [S ii], and [Ar iii], given that the broadband WISE/W1 flux would also be contaminated. at 3.5 $\micron$ of $S_{\rm 3.5}$ = 27.3 $\pm$ 1.0 mJy, which was obtained with the Palomar Telescope (Egami et al. 2000).

Following Riechers et al. (2009), we first fit a simplified four-component model (two-temperature dust, free-free, and synchrotron) to the rest-frame FIR-to-radio photometry, covering observed frame 53 $\micron-90$ cm, to illustrate how well the previously employed functional form fits all of the data obtained to date. The functional form is expressed as follows:

[TABLE]

where $\nu$ is the frequency in GHz, $c_{i}$ is a coefficient, $\alpha_{i}$ is a spectral index, and $i$ represents either the flat (free-free) or steep (synchrotron) power-law components. The first term takes the functional form of a two-temperature MBB model, with each dust temperature component described by

[TABLE]

where $\Omega_{\rm app}$ is the apparent solid angle (the magnification factor is absorbed into this parameter). The optical depth of each component is parameterized through the dust mass ($M_{\rm dust}$) and area filling factor ($F_{\rm filling}$):

[TABLE]

The area filling factor clearly takes a value in the interval [0, 1]. We define $x\equiv M_{\rm dust}/r_{0}^{2}$ to reduce the extra dimension and degeneracies between $M_{\rm dust}$ and $r_{0}$. For the second term of Equation 2, we adopt the best-fitting results from Riechers et al. (2009) since no new constraints are added to the radio SED in this work.

We impose uniform priors on the model parameters as follows: 10 K$<T_{d,\rm cold}<T_{d,\rm warm}$, $T_{d,\rm warm}\leq 500$ K, F${}_{\rm Filling}^{\rm cold}\in[0,1]$, F${}_{\rm Filling}^{\rm warm}\in[0,1]$, 107 $M_{\odot}$  pc${}^{-2}<x_{\rm cold}<10^{10}$ $M_{\odot}$  pc*-2*, and 104 $M_{\odot}$  pc${}^{-2}<x_{\rm warm}<$1010 $M_{\odot}$  pc*-2*. We also impose the criterion that $x_{\rm warm}<x_{\rm cold}$. These conditions are physically motivated based on previous results of Beelen et al. (2006) and Weiß et al. (2007). We find the best-fit model by sampling the posterior probability distributions (PDF) using a Markov Chain Monte Carlo (MCMC) approach.

After discarding the first 2000 iterations in the burn-in phase, we thin the chains via 40 iterations based on the auto-correlation time to obtain independent samples from the posterior PDFs, yielding a best-fit model with $T_{d,\rm cold}$ = 110${}^{+3}_{-3}$ K, $T_{d,\rm warm}$ = 296${}^{+17}_{-13}$ K, $F_{\rm Filling}^{\rm cold}$ = 0.9${}^{+0.1}_{-0.1}$, $x_{\rm warm}$ = 2.5{}^{+1.2}_{-1.0}$$\times 10^{7} $M_{\odot}$  pc*-2*, and $x_{\rm cold}$ = 1.4{}^{+0.1}_{-0.1}$$\times 10^{9} $M_{\odot}$  pc*-2*. The best-fit values and their uncertainties are quoted based on the 16th, 50th, and 84th percentiles of the PDFs (See Appendix). The best-fit SED is shown in Figure 3, where we overplot the model presented by Riechers et al. (2009)555 The apparent difference in the R09 fit shown here compared to Riechers et al. (2009) results from the fact that we adopt a source size of $r_{0}$ = 1150 pc, following Weiß et al. (2007, cf. $r_{0}$ = 1300 pc). The R09 fit here therefore appears to undershoot the photometry at $\lambda_{\rm rest}$ $\sim$ 100 $\micron$.

(hereafter, R09 fit) for comparison.

4.2. Semi-analytic: Including Clumpy Torus Models

Previous works suggest that APM 08279$+$5255 is a typical example of a co-existing starburst and quasar in a galaxy, where its strong FIR luminosity due to cold dust emission has been attributed to the presence of a dust-obscured starburst, with a lensing-corrected SFR of 5000 $\times$ $(4/\mu_{L})$ $M_{\odot}$  yr*-1*. In such AGN host galaxies, hot dust emission with a characteristic temperature of $\sim 1500$ K (e.g., Barvainis 1987; Deo et al. 2011) due to reprocessed X-ray-to-optical AGN radiation is expected to peak in the NIR between 1 and 3 $\micron$, whereas warm dust emission from a dusty torus near the central AGN with characteristic temperature of $\sim 400$ K (Schartmann et al. 2005; Burtscher et al. 2015) should peak between 3 and 40 $\micron$. Based on the observed SED (Figure 3), hot/warm dust emission due to a putative torus in APM 08279$+$5255 is likely an important contributor to its IR emission (since its IR SED clearly peaks around rest-frame 20 $\micron$). We thus describe the SED of APM 08279$+$5255 (including the IRS spectrum) using a physically motivated model, with the UV-to-MIR part corresponding to emission from the accretion disk and and clumpy torus around it, the FIR part described by a single-temperature MBB function, and the radio part corresponding to free-free and synchrotron emission. The cold dust and radio components are described using the same parameterization as Equations 2 and 3. For the dust component, we assume an opacity defined using Equation 1. The spectral index $\beta$ is fixed to 2.0, following Priddey & McMahon (2001), to minimize the number of free parameters. For the AGN/torus component, we employ SED grids derived by performing radiative transfer calculations of three-dimensional clumpy tori (Hönig & Kishimoto 2010). In this model, the optical depths and sizes of the dust clouds depend on their distance from the central source. The clouds are randomly placed according to a spatial distribution function. The torus SED is then calculated by summing the direct and indirect emission of the clumps. More specifically, the model parameterizes the distribution of dust clouds in the tori using five parameters: the radial dust-cloud distribution power-law index, $a$; the half-opening angle, $\theta$ (describing the vertical distribution of clouds); the average number of clouds along an equatorial line-of-sight, $N_{0}$; the total V-band optical depth, $\tau_{V}$; and the inclination angle, $i$. The MIR SED is redder if more dust clouds are distributed at larger distances. Therefore, the power-law index $a$ is mostly constrained by the MIR photometry (a more negative $a$ would imply a more compact dust distribution; see Hönig et al. 2010 for details). The total number of clouds within the torus is thus related to $N_{0}$ through

[TABLE]

where $\eta(r,z,\phi)$ describes the cloud distribution in the $r$-, $z$-, and $\phi$-space; $R_{\rm cl}$ is the sublimation radius of a given cloud; and $V$ is the volume of the torus. The torus models are scaled based on a luminosity factor ($L_{\rm scale}$), which is the luminosity of the accretion disk. Physically, this determines the sublimation radius of the accretion disk/torus where dust clumps are distributed, according to the following scaling properties:

[TABLE]

where $T_{\rm sub}$ is the sublimation temperature (e.g., Barvainis 1987). The numerical values in the denominators are the physical properties assumed in the model. In the fitting process, we rescale the model flux to the distance of APM 08279$+$5255. We initialize $L_{\rm scale}$ based on the apparent X-ray luminosity (2$-$10 keV; not lensing-corrected) of APM 08279$+$5255, which is $\log\left(L_{X}/L_{\odot}\right)$ = 49.4 (Gofford et al. 2015).

Since the Hönig et al. (2010) templates are calculated based on discrete model parameters, we interpolate the SED grids such that each parameter can be sampled as a continuous variable. This is done in order to determine the PDFs of the parameters and thus obtain reliable uncertainties on the model parameters. The parameters are interpolated in six dimensions (including the dimension of wavelength) within the following boundaries: $N_{0}\in$[0, 10], $a\in[-2.5,0.5]$, $\theta\in[0,65]\arcdeg$, $\tau_{V}\in$[15, 90], and $i\in[0,90]\arcdeg$.

We place priors on the model parameters to ensure they are physically sensible. We place flat priors of $M_{d}^{\rm cold}\in[10^{7},10^{10}]$ $M_{\odot}$ on the dust mass and $T_{d}^{\rm cold}\in[30,300]$ K on the dust temperature. The luminosity scaling factor is allowed to vary between $\log\left(L_{\rm scale}/L_{\odot}\right)\in[47,49.4]$. In addition to unobscured torus emission, stellar emission from the host galaxy may be significant shortward of rest-frame $\sim 1\,\micron$. Modeling the stellar emission would require introducing multiple additional free parameters and properly accounting for potential differential magnification of each of the components. Thus, for simplicity, we fit the models to the photometry longward of 1 $\mu$m only.

The above model fails to yield a good fit to the NIR photometry — the torus component is skewed to the rest-frame NIR wavebands and over-predicts the fluxes observed in near- and mid-IR bands (top panel of Figure 4). Motivated by the work of Leipski et al. (2014), we add an additional hot dust component to the overall model. For this hot dust component, we place flat priors of $r_{\rm hot}\in[10,1300]$ pc on its emitting radius, $M_{d}^{\rm hot}\in[10,10^{5}]$ $M_{\odot}$ on its dust mass, and $T_{d}^{\rm hot}\in[1000,1800]$ K on its temperature.

After sampling the target posteriors using MCMC with 100 walkers and 5$\times 10^{5}$ iterations, and after discarding the first 2$\times 10^{5}$ iterations as the burn-in phase and thinning the chains by 50 iterations, we identify the best-fit parameters together with the their uncertainties based on the 16th, 50th, and 84th percentiles of the PDFs. The best-fit model is shown in the bottom panel of Figure 4. We summarize the best-fit parameters and the apparent luminosities integrated over the total IR, NIR (rest-frame 1$-$3 $\micron$), MIR (rest-frame 3$-$40 $\micron$), and FIR (rest-frame 42.5$-$122.5 $\micron$) in the bottom section of Table 2.

We find an FIR luminosity of $L_{\rm FIR}^{\rm Torus}$ = 4.7{}^{+0.3}_{-0.5}$$\times 10^{13} $\mu_{L}$-1  $L_{\odot}$ for the torus component. This corresponds to $\sim$30% of the total FIR emission, which is in contrast to the 90% reported in previous work (Weiß et al. 2007). Correcting for the torus contribution to the FIR emission, we find an SFR of $3075$ $\times$ $(4/\mu_{L})$ $M_{\odot}$  yr*-1* (assuming the Chabrier 2003 initial mass function). This is consistent with that derived from the two-temperature MBB model, which is SFR = 2970 $\times$ $(4/\mu_{L})$ $M_{\odot}$  yr*-1*. We note that the IR luminosity derived from this model is $\sim$2 times higher than that yielded by the fully analytic two-temperature dust model (§4.1). This discrepancy results mostly from including a hot dust and a dusty torus component in the latter model, which is physically motivated and better describes the IR SED of APM 08279$+$5255 (see Figure 3 and 4). That is, the latter model describes the rest-frame near- and mid-IR photometry and the IRS spectrum, whereas the former does not (by construction). This result highlights the importance of accounting for the NIR emission using physically motivated models.

Including the near- and far-IR emission (i.e., covering rest-frame 1$-$122.5 $\micron$), we find that dust emission due to the central AGN accounts for $\sim$93% of the total IR luminosity. This is similar to the large fraction (83%) inferred for the warm dust component in the fully analytic model. Based on the latest lens model (Riechers et al. 2009), the magnification factor is found to vary by only a factor of $\mu_{L}$ $\sim$ 2$-$3 across a spatial extent of 20 kpc. Thus, differential magnification may cause the above percentage to be uncertain to a factor of 2$-$3. However, if the cold dust is concentrated within the central 1 kpc (see Weiß et al. 2007 and Riechers et al. 2009), the difference in the magnification factor between the AGN and more-extended cold dust would be $\lesssim$25%. With this uncertainty in mind, dust heating due to the central quasar appears to be more significant than that from the star formation in the host galaxy of APM 08279$+$5255, and the fraction of the IR luminosity associated with the AGN is greater than those observed in other high-$z$ AGN/starburst composite systems studied to-date (e.g., Pope et al. 2008; Coppin et al. 2010; Lyu et al. 2016).

As noted above, the SED shortward of $\sim 1\,\micron$ potentially has contributions from both the AGN and stars in the host galaxy, and we have not attempted to fit the photometry shortward of 1 $\micron$ for this reason. If we assume that the residual flux density at 1 $\micron$ unaccounted for by the torus model results from stellar emission in the host galaxy, we infer a stellar mass of $\mu_{L}M_{*}$ = (3.6$-$5.3) $\times$ 1013 $M_{\odot}$. The range of stellar mass results from the range of $H$-band mass-to-light ratio of $L_{H}/M_{*}$ = $5.8-8.5$ ($L_{\odot}$ $M_{\odot}^{-1}$) adopted (e.g., Hainline et al. 2011), which corresponds to different population synthesis models (Bruzual & Charlot 2003 and Maraston 2005) and star formation histories (instantaneous star formation versus constant). Taken at face value, the stellar mass of APM 08279$+$5255 is $\mu_{L}M_{*}\simeq$ 4 $\times$ 1013$M_{\odot}$, or $10^{13}$ $M_{\odot}$ assuming a lensing magnification factor of $\mu_{L}$ = 4. Such a stellar mass would place APM 08279$+$5255 among the most massive galaxies at $z$ $\sim$ 4 known. Alternatively, the very large inferred stellar mass could suggest that the lensing magnification for the stellar component is (at least) an order of magnitude higher than the Riechers et al. (2009) value of $\mu_{L}\approx 4$ or that the residual flux is not due entirely to stellar emission, perhaps because the torus models employed do not adequately capture the geometry or/and dust properties of APM 08279$+$5255. Currently, the data at hand preclude us from favoring either scenario.

5. Discussion

In the case of the two-temperature dust SED model, we find that previous model consistently underestimates the more comprehensive set of IR photometry included in this work, especially in the rest-frame wavelength range covering $\lambda_{\rm rest}\in$[10, 100] $\mu$m (Figure 3), with a maximum residual of 53%. The differences is mainly a direct consequence of the higher dust temperatures found here (110 K vs. 65 K and 296 K vs. 220 K). In contrast to the results reported by Weiß et al. (2007), we find that the cold dust component dominates the FIR luminosity of APM 08279$+$5255 ($\sim$ 66%). Regardless of these differences, we find comparable apparent total IR luminosities (rest-frame 1$-$1000 $\micron$). The (F)IR luminosity is summarized in the top section of Table 2.

When fitting the SED with clumpy AGN torus models in addition to the cold dust and radio components, we find a torus inclination angle consistent with the Type-1 nature of APM 08279$+$5255 (under the AGN unification model). The best-fit model suggests the presence of silicate emission at rest-frame 10 and 18 $\micron$666The 10 $\micron$ feature is attributed to stretching of the Si-O bonds in silicates, whereas the 18 $\micron$ feature is attributed to the O-Si-O bending mode., as commonly observed in Type-1 AGNs with Spitzer (e.g., Hao et al. 2005; Siebenmorgen et al. 2005; Sturm et al. 2005). If these features are real, their presence is consistent with the AGN unification model, in which an AGN torus seen close to pole-on ($i$ = 0$\arcdeg$) should exhibit silicate emission features in the MIR (e.g., Hao et al. 2005; see also Stalevski et al. 2012), as dust on the surface of the torus is heated to temperatures sufficiently hot to produce silicate emission features. However, if these features are fictitious, this could be evidence that different torus models are needed to describe AGN/starburst systems at high redshifts.

As mentioned in §4.2, the best-fit semi-analytic model including the torus component without an extra hot dust component overestimates the near-to-mid-IR fluxes. The MIR slope is also inconsistent with the IRS spectrum. We follow the work by Leipski et al. (2014) for $z$ $\sim$6 quasars and include an extra hot dust component in the model. The extra hot dust component improves the fit and yields good agreement with the IR photometry. However, the physical motivation for this extra hot dust component is unclear, as the hottest dust emission should arise self-consistently in the clumpy AGN torus models computed via radiative transfer, and there should be no need to add an additional ad hoc hot dust component. This need for an additional hot dust component may suggest that these powerful quasars in the early Universe have very different dust compositions and distributions than the more-nearby systems used to motivate AGN torus models. If this issue persists for larger samples of high-$z$ AGN, further improvements in AGN torus models would be warranted.

6. Summary and Conclusions

We observe the IR-luminous, strongly-lensed AGN-starburst system APM 08279$+$5255 at $z$ = 3.9114 using the new sensitive FIR camera HAWC$+$ onboard SOFIA and the IRS onboard Spitzer. While APM 08279$+$5255 is a Type-1 broad absorption line AGN, it clearly shows strong mid-IR emission (even considering the viewing angle of APM 08279$+$5255 in the classical AGN unification model). Hence, SED modeling tools such as cigale are inappropriate777 Modeling tools such as cigale expect Type-1 AGNs to emit predominately in UV and optical wavebands and Type-2 AGNs to emit predominately in MIR wavelengths. Since a large amount of dust is still present in APM 08279$+$5255, the emission in the MIR waveband is significantly underestimated in cigale when one properly considers it as a Type-1 AGN. for APM 08279$+$5255. We combine the new data with ancillary data and decompose the IR SED into hot dust, torus, and cold dust components and the radio SED into thermal and non-thermal components (free-free and synchrotron). To facilitate comparison with previous studies, we also fit simplified two-temperature modified blackbody models.

We constrain the structure of the putative dusty torus around the central AGN by fitting 3D radiative transfer models of clumpy tori. Despite the uncertainties and degeneracies between the model parameters, some results are clear from the fit. For instance, the best-fit inclination angle and the silicate features predicted by the best-fit model are consistent with the Type-1 AGN nature of APM 08279$+$5255, which may be evidence supporting the unification model of AGN.

In the modified blackbody model presented in this work, we find that the warm dust component contributes only $34\%$ of the far-IR luminosity. Accounting for the torus’ contribution to the FIR luminosity, we find a lensing-corrected SFR of $3075$ $\times$ $(4/\mu_{L})$ $M_{\odot}$  yr*-1*. Differences compared to previous models are pointed out in §4.2.

Accounting for emission due to the torus at near-to-far-IR wavelengths (rest-frame 1$-$122.5 $\micron$), we find that dust heating from the quasar contributes about 93% to its total IR emission. Therefore, the AGN plays an important role in driving the total dust luminosity of APM 08279$+$5255, more so than in other high-$z$ AGN/starburst composite systems studied to date (e.g., Pope et al. 2008; Coppin et al. 2010; Lyu et al. 2016). Based on the residual flux density at 1 $\micron$ unaccounted for by the torus model, we find a stellar mass of $M_{*}$ $\simeq$ (4/\mu_{L})$$\times 10^{13} $M_{\odot}$. At face value, assuming the Riechers et al. (2009) value of $\mu_{L}\approx 4$, the resulting stellar mass would place APM 08279$+$5255 among the most evolved and massive galaxies at $z$ $\sim$ 4 known, but it may also indicate that the true magnification is significantly greater than 4 and/or that the best-fit AGN torus model is significantly underestimating the contribution of the AGN at $\sim$1 $\micron$.

Our results demonstrate the use of the new and sensitive FIR camera onboard SOFIA to constrain the physical properties of the dusty tori of quasars in the early Universe. The differences between the simplified model and the more complex but physically motivated model highlight the need for well-sampled IR SEDs of high-$z$ quasars in order to disentangle the effect of dust heating due to AGN versus star formation and thus to more accurately determine the SFRs in AGN/starburst composite systems. We note that the high IR luminosity of APM 08279$+$5255 may imply that our results are unlikely to be representative of high-$z$ AGNs. The Origins Space Telescope (OST) will have the capabilities to characterize the torus and dust components in high-$z$ AGNs with a much broader range of properties, enabling similar types of studies to be carried out for more typical and representative systems.

We thank the referee for the careful reading and for providing constructive comments on the manuscript. The authors thank Dominik Riechers for his contributions to the SOFIA and Spitzer proposals that resulted in observations used here, and for access to data and data products used in this work. The authors also thank Lee Armus, Chris Carilli, Aaron Evans, Patrick Ogle, and Fabian Walter for their role in this study, and Agatha Hodsman for providing the reduction of the Spitzer/IRS data that is used in this work. TKDL thanks Denis Burgarella for helpful discussion about the SED modeling code cigale. TKDL acknowledges support from NASA through the New York Space Grant, award NAS2-97001 from the SOFIA MSO, and support from the Simons Foundation. TKDL acknowledges the hospitality of NASA’s Goddard Space Flight Center and the Space Telescope Science Institute, where part of this manuscript was prepared. The Flatiron Institute is funded by the Simons Foundation. CMC thanks the University of Texas at Austin College of Natural Sciences, NSF grants AST-1714528 and AST-1814034, as well as support from a 2019 Cottrell Scholar Award from the Research Corporation for Science Advancement. This work is based in part on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. Financial support for this work was provided by NASA through award #SOF 05-0165 issued by USRA. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. This work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. The IRS was a collaborative venture between Cornell University and Ball Aerospace Corporation funded by NASA through the Jet Propulsion Laboratory and Ames Research Center. This research has made use of the NASA/ IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. Facilities: SOFIA (HAWC$+$), Spitzer (IRS, IRAC, MIPS), IRSA, Herschel (PACS, SPIRE)