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

Observations of APM 08279+5255 were carried out with the IRS on board 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 μm) was used, covering the rest-frame 6.7 μm 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."⁸

The Basic Calibrated Data 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 subtraction.⁹ Sky frames at each map position are 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 Δλ =0.169 μm wavelength bin (observed frame).

Observations of APM 08279+5255 with SOFIA/HAWC+ were performed on 2017 October 18 at an altitude of ∼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 μm (Band D). Total intensity observing in On-the-fly Mapping mode was used, producing a continuous telescope motion with 1.5 minutes Lissajous scans, covering (regions slightly smaller than) the standard HAWC+ fields-of-view of 27 × 17 (Band A), 42 × 26 (Band C), and 73 × 45 (Band D). Five scans corresponding to a total integration time of 333 s were obtained in Band A, 12 scans totaling 870 s were obtained in Band C, and 12 scans totaling 936 s were obtained 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° and 37°. No direct measure of precipitable water vapor 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 overburdened at the zenith. The assumed zenith opacities are τ_z = 0.11, 0.21, 0.24 for Bands A, C, and D, respectively.

The data were reduced and processed using the 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 71, 110, and 195 at 53, 89, and 154 μm, respectively. After reduction, the pixel scales are 100, 155, and 275 for Bands A, C, and D, respectively.

In the analysis, we include photometry published in the literature, as listed in Table 1. 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, which are available through the NASA/IPAC Infrared Science Archive. The photometry in all IRAC bands reported in Table 1 was extracted using a 38 diameter aperture, with aperture corrections applied. The MIPS photometry at 24 μm 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, tabulated in the Herschel/PACS point-source catalog (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 μm (Table 1) 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.

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₀) constrained through large velocity gradient 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_filling = 0.75 and 0.25 for the warmer and cooler components, respectively, and a dust absorption coefficient of

$$\begin{eqnarray}&&{\kappa }_{\nu }={\kappa }_{0}{\left(\nu /{\nu }_{0}\right)}^{\beta },\end{eqnarray} \tag{ 1 }$$

where κ₀ = 0.4 cm² g⁻¹, ν₀ = 250 GHz, and the emissivity spectral index β is fixed to 2. They find that the warmer component contributes ∼85% to the FIR luminosity (${L}_{\mathrm{FIR}}$ =2 × 10¹⁴ μ_L⁻¹ L_⊙) 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 ∼100 μm, which has a ∼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 outlier¹¹ at 3.5 μm of S_3.5 = 27.3 ± 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 μm–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:

$$\begin{eqnarray}&&{S}_{\nu }(\nu )={S}_{\nu ,d}(\nu )+\displaystyle \sum _{i=\mathrm{flat},\mathrm{steep}}{c}_{i}\times {\nu }^{{\alpha }_{i}},\end{eqnarray} \tag{ 2 }$$

where ν is the frequency in GHz, c_i is a coefficient, α_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

$$\begin{eqnarray}&&{S}_{\nu }=[{B}_{\nu }({T}_{d})-{B}_{\nu }({T}_{\mathrm{CMB}})]\displaystyle \frac{(1-{\exp }^{-{\tau }_{v}})}{{\left(1+z\right)}^{3}}{{\rm{\Omega }}}_{\mathrm{app}},\end{eqnarray} \tag{ 3 }$$

where Ω_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_dust) and area filling factor (F_filling):

$$\begin{eqnarray}&&{\tau }_{\nu }=\displaystyle \frac{{\kappa }_{\nu }\,{M}_{\mathrm{dust}}}{(\pi \,{r}_{0}^{2})\,{F}_{\mathrm{filling}}}.\end{eqnarray} \tag{ 4 }$$

The area filling factor clearly takes a value in the interval [0, 1]. We define $x\equiv {M}_{\mathrm{dust}}/{r}_{0}^{2}$ to reduce the extra dimension and degeneracies between M_dust and r₀. For the second term of Equation (2), we adopt the best-fitting results from Riechers et al. (2009) because 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,cold < T_d,warm, T_d,warm ≤ 500 K, ${F}_{\mathrm{Filling}}^{\mathrm{cold}}\,\in [0,1]$, ${F}_{\mathrm{Filling}}^{\mathrm{warm}}\in [0,1]$, 10⁷ M_⊙ pc⁻² < x_cold < 10¹⁰ M_⊙ pc⁻², and 10⁴ M_⊙ pc⁻² < x_warm < 10¹⁰ M_⊙ pc⁻². We also impose the criterion that x_warm < x_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 (PDFs) 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,cold = 110${}_{-3}^{+3}$ K, T_d,warm = 296${}_{-13}^{+17}$ K, ${F}_{\mathrm{Filling}}^{\mathrm{cold}}$ = 0.9${}_{-0.1}^{+0.1}$, x_warm = 2.5${}_{-1.0}^{+1.2}$ $\,\times {10}^{7}$ M_⊙ pc⁻², and x_cold = 1.4${}_{-0.1}^{+0.1}$ $\ \times \ {10}^{9}$ M_⊙ pc⁻². The best-fit values and their uncertainties are quoted based on the 16th, 50th, and 84th percentiles of the PDFs (see the Appendix). The best-fit SED is shown in Figure 3, where we overplot the model presented by Riechers et al. (2009)¹² (hereafter, R09 fit) for comparison.

Zoom In Zoom Out Reset image size

Previous works suggested that APM 08279+5255 is a typical example of a coexisting 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 × (4/μ_L) M_⊙ yr⁻¹. In such AGN host galaxies, hot dust emission with a characteristic temperature of ∼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 μm, whereas warm dust emission from a dusty torus near the central AGN with a characteristic temperature of ∼400 K (Schartmann et al. 2005; Burtscher et al. 2015) should peak between 3 and 40 μm. 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 (because its IR SED clearly peaks around rest frame 20 μm). 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 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 β 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, θ (describing the vertical distribution of clouds); the average number of clouds along an equatorial line of sight, N₀; the total V-band optical depth, τ_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₀ through

$$\begin{eqnarray}&&{N}_{\mathrm{tot}}\propto \int {N}_{0}\,\displaystyle \frac{\eta (r,z,\phi )}{{R}_{\mathrm{cl}}^{2}}\,{dV},\end{eqnarray} \tag{ 5 }$$

where η (r, z, ϕ) describes the cloud distribution in the r-, z-, and ϕ-space; R_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_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:

$$\begin{eqnarray}&&\left(\displaystyle \frac{{R}_{\mathrm{cl}}}{0.5\,\mathrm{pc}}\right)\propto {\left(\displaystyle \frac{{L}_{\mathrm{scale}}}{{10}^{46}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{\displaystyle \frac{1}{2}}{\left(\displaystyle \frac{{T}_{\mathrm{sub}}}{1500{\rm{K}}}\right)}^{-2.8},\end{eqnarray} \tag{ 6 }$$

where T_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_scale based on the apparent X-ray luminosity (2–10 keV; not lensing-corrected) of APM 08279+5255, which is $\mathrm{log}\left({L}_{{\rm{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, 10], a ∈ [−2.5, 0.5], θ ∈ [0, 65]°, τ_V ∈ [15, 90], and i ∈ [0, 90]°.

We place priors on the model parameters to ensure they are physically sensible. We place flat priors of ${M}_{d}^{\mathrm{cold}}\in [{10}^{7},{10}^{10}]$ M_⊙ on the dust mass and ${T}_{d}^{\mathrm{cold}}\in [30,300]$ K on the dust temperature. The luminosity scaling factor is allowed to vary between $\mathrm{log}\left({L}_{\mathrm{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 ∼1 μm. 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 μ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 overpredicts 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_hot ∈ [10, 1300] pc on its emitting radius, ${M}_{d}^{\mathrm{hot}}\in [10,{10}^{5}]$ M_⊙ on its dust mass, and ${T}_{d}^{\mathrm{hot}}\in [1000,1800]$ K on its temperature.

Zoom In Zoom Out Reset image size

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 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 μm), MIR (rest frame 3–40 μm), and FIR (rest frame 42.5–122.5 μm) in the bottom section of Table 2.

Table 2.  SED Model Parameters for APM 08279+5255

Properties	Units	Values
Analytical: Two-temperature MBB+Radio Power Lawsa
Td,warm	(K)	296${}_{-15}^{+17}$
Td,cold	(K)	110${}_{-3}^{+3}$
${F}_{\mathrm{Filling}}^{\mathrm{cold}}$		0.9${}_{-0.1}^{+0.1}$
xwarmb	(107 M⊙ pc−2)	2.5${}_{-1.0}^{+1.2}$
xcoldb	(109 M⊙ pc−2)	1.4${}_{-0.1}^{+0.1}$
LIR	(1015μ−1 L⊙)	1.8${}_{-0.6}^{+0.7}$
LFIR	(1014μ−1 L⊙)	1.8${}_{-0.2}^{+0.1}$
Semi-analytical: Including Clumpy Torusc
rhot	(pc)	655${}_{-409}^{+418}$
${M}_{d}^{\mathrm{hot}}$	(μ−1 M⊙)	11${}_{-1}^{+5}$
${M}_{d}^{\mathrm{cold}}$	(109μ−1 M⊙)	1.4${}_{-0.2}^{+0.3}$
${T}_{d}^{\mathrm{hot}}$	(K)	1385${}_{-49}^{+10}$
${T}_{d}^{\mathrm{cold}}$	(K)	100${}_{-1}^{+2}$
N0		8.3${}_{-2.2}^{+1.3}$
a		−1.0${}_{-0.1}^{+0.1}$
θ	(°)	55${}_{-11}^{+4}$
τV		42${}_{-5}^{+5}$
i	(°)	15${}_{-8}^{+8}$
log(Lscale/L⊙)		48.8${}_{-0.1}^{+0.1}$
LIR	(1015μ−1 L⊙)	3.5${}_{-0.1}^{+0.1}$
LNIR	(1014μ−1 L⊙)	5.9${}_{-0.6}^{+0.4}$
LMIR	(1014μ−1 L⊙)	2.7${}_{-0.7}^{+0.3}$
LFIR	(1014μ−1 L⊙)	1.7${}_{-0.1}^{+0.5}$
${L}_{\mathrm{FIR}}^{\mathrm{Torus}}$	(1013μ−1 L⊙)	4.7${}_{-0.5}^{+0.3}$
${L}_{\mathrm{FIR}}^{\mathrm{Cold}}$	(1014μ−1 L⊙)	1.2${}_{-0.7}^{+0.5}$

Notes. The different IR luminosities are obtained by integrating over the rest frames: 1–1000 μm for L_IR; 1–3 μm for L_NIR; 3–40 μm for L_MIR; and 42.5–122.5 μm for L_FIR. Uncertainties quoted are derived based on the 16th and 84th percentiles of the posterior PDFs (see Figures 5 and 6).

^aSee Section 4.1. ^b $x\equiv {M}_{d}/{r}_{0}^{2}$; see Equation (4). ^cSee Section 4.2.

Download table as:  ASCIITypeset image

We find an FIR luminosity of ${L}_{\mathrm{FIR}}^{\mathrm{Torus}}$ = 4.7${}_{-0.5}^{+0.3}$ $\ \times \ {10}^{13}$ μ_L⁻¹ L_⊙ for the torus component. This corresponds to ∼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 ×(4/μ_L) M_⊙ yr⁻¹ (assuming the Chabrier 2003 initial mass function). This is consistent with that derived from the two-temperature MBB model, which is SFR = 2970 ×(4/μ_L) M_⊙ yr⁻¹. We note that the IR luminosity derived from this model is ∼2 times higher than that yielded by the fully analytic two-temperature dust model (Section 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 Figures 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 μm), we find that dust emission due to the central AGN accounts for ∼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 μ_L ∼ 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 ≲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 ∼1 μm 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 μm for this reason. If we assume that the residual flux density at 1 μm unaccounted for by the torus model results from stellar emission in the host galaxy, we infer a stellar mass of μ_LM_* = (3.6–5.3) × 10¹³ M_⊙. 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_⊙ ${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 μ_LM_* ≃ 4 × 10¹³ M_⊙, or 10¹³ M_⊙ assuming a lensing magnification factor of μ_L = 4. Such a stellar mass would place APM 08279+5255 among the most massive galaxies at z ∼ 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 μ_L ≈ 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.