Abstract
We present the detection at 89 μm (observed frame) of the Herschel-selected gravitationally lensed starburst galaxy HATLAS J1429-0028 (also known as G15v2.19) in 15 minutes with the High-resolution Airborne Wideband Camera-plus (HAWC+) onboard the Stratospheric Observatory for Infrared Astronomy (SOFIA). The spectacular lensing system consists of an edge-on foreground disk galaxy at z = 0.22 and a nearly complete Einstein ring of an intrinsic ultra-luminous infrared (IR) galaxy at z = 1.03. Is this high IR luminosity powered by pure star formation (SF) or also an active galactic nucleus (AGN)? Previous nebular line diagnostics indicate that it is star formation dominated. We perform a 27-band multiwavelength spectral energy distribution (SED) modeling including the new SOFIA/HAWC+ data to constrain the fractional AGN contribution to the total IR luminosity. The AGN fraction in the IR turns out to be negligible. In addition, J1429-0028 serves as a testbed for comparing SED results from different models/templates and SED codes (magphys, sed3fit, and cigale). We stress that star formation history is the dominant source of uncertainty in the derived stellar mass (as high as a factor of ∼10) even in the case of extensive photometric coverage. Furthermore, the detection of a source at z ∼ 1 with SOFIA/HAWC+ demonstrates the potential of utilizing this facility for distant galaxy studies including the decomposition of SF/AGN components, which cannot be accomplished with other current facilities.
Export citation and abstract BibTeX RIS
1. Introduction
High-redshift dusty star-forming galaxies (DSFGs) are characterized by their very high intrinsic infrared (IR) luminosities, LIR ∼ 1012–1013 L⊙, and dust enshrouded intense star formation (SF) with star formation rates (SFRs) exceeding ∼100–1000 M⊙ yr−1 (for a review see Casey et al. 2014). Large area far-IR (FIR) to submillimeter/millimeter surveys such as the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010; Bussmann et al. 2013), the Herschel Multi-Tiered Extragalactic Survey (HerMES; Oliver et al. 2012; Bussmann et al. 2015), and the South Pole Telescope Sub-millimeter Galaxy (SMG) survey (Vieira et al. 2013; Spilker et al. 2016; Strandet et al. 2016) have discovered large samples of gravitationally lensed DSFGs whose properties are being studied in unprecedented detail thanks to the significant gain in both brightness and spatial resolution provided by gravitational lensing (Bussmann et al. 2013; Calanog et al. 2014; Spilker et al. 2016).
DSFGs are in a unique phase of galaxy formation and evolution and provide a laboratory for studying the co-evolution of galaxies and their super massive black holes (Alexander & Hickox 2012). X-ray is arguably the best indicator of active galactic nuclei (AGNs). X-ray observations of DSFGs show that the majority have bolometric luminosities that are dominated by star formation rather than AGN emission (e.g., Alexander et al. 2005a; Wang et al. 2013). However, these galaxies can be highly obscured by large columns of dust and gas and therefore are difficult to detect in the X-ray (e.g., Alexander et al. 2005b). Radio data provide another route to distinguishing star-forming galaxies from AGN. Star-forming galaxies follow a tight FIR-to-radio correlation over five orders of magnitude in galaxy luminosity, while radio-loud AGN produce excess radio emission above this relation (e.g., Yun et al. 2001). Alternatively, emission line diagnostics such as optical nebular emission lines can serve as indirect AGN indicators. However, rest-frame optical emission lines are often not detectable in these highly dust-obscured systems. None of these indicators provide a quantitative AGN fraction that can be used to correct for AGN contamination in measured LIR, SFR, and stellar mass, which are fundamental parameters that determine the nature of DSFGs.
Galaxy emission in the mid-IR contains rich information about the underlying sources within the galaxy, including polycyclic aromatic hydrocarbon (PAH) features (e.g., Pope et al. 2008; Riechers et al. 2014), which trace star formation, and additional hot dust emission from around any AGN. Combined with photometric data covering ultraviolet (UV), optical, near-IR (IR) and FIR multiwavelength spectral energy distribution (SED) modeling techniques are powerful tools to decompose star formation and AGN activity, and to quantitatively constrain the AGN fraction, in addition to providing self-consistent constraints on stellar masses, SFRs, stellar ages, dust properties, etc. (da Cunha et al. 2008; Noll et al. 2009; Kirkpatrick et al. 2012; Ciesla et al. 2015; Leja et al. 2017; Małek et al. 2017).
This paper focuses on a strongly lensed DSFG (magnification factor μ ∼ 10), HATLAS J142935.3-002836 (J1429-0028 hereafter; H1429-0028 in Messias et al. 2014; G15v2.19 in Calanog et al. 2014), standing out from the GAMA-15 field of Herschel-ATLAS (Eales et al. 2010) due to its extremely bright fluxes with S160 μm = 1.1 ± 0.1 Jy. This system consists of an edge-on foreground disk galaxy at z = 0.218 and a nearly complete Einstein ring of an intrinsic ultra-luminous IR galaxy at z = 1.027. The brightness enables detections from the optical to radio wavelengths, construction of lens models, and de-blending the source/lens photometry (Calanog et al. 2014; Messias et al. 2014). Messias et al. (2014) presented a detailed analysis of the gas and dust properties, dynamical information, and derive galaxy parameters from SED modeling of J1429-0028. However, a few critical questions remain unaddressed: (1) Is there an energetically important AGN in the system and what is the fractional contribution to the total IR luminosity, if any? (2) How did the galaxy build up the stellar mass, i.e., are we able to constrain the star formation history (SFH)? Given that we have obtained extensive photometric coverage from UV to radio, the best case among all DSFGs, J1429-0028 can serve as a laboratory to test systematic uncertainties in deriving physical properties via panchromatic SED modeling.
In this paper, we present new photometric observations from the High-resolution Airborne Wideband Camera-plus (HAWC+; Dowell et al. 2013; Smith et al. 2014) onboard the Stratospheric Observatory for Infrared Astronomy (SOFIA; Temi et al. 2014), which is the only current facility that covers the 50–100 μm wavelength range, a crucial regime for disentangling SF/AGN activity in galaxies at redshifts of z ∼ 0.5–2. In combination with our existing multiwavelength data, we will address the remaining questions about J1429-0028 in order to probe the nature of this DSFG in unprecedented detail. The paper is organized as follows. The SOFIA observations and data reduction are presented in Section 2. We perform multiwavelength SED modeling including AGN templates in Section 3. The results are discussed in Section 4 and the conclusions are summarized in Section 5. We adopt a Chabrier (2003) initial mass function (IMF) throughout the paper. We assume a ΛCDM cosmological model with H0 = 70 km s−1 Mpc−1, ΩM = 0.3, and ΩΛ = 0.7. The IR luminosity is integrated over a rest-frame of 8–1000 μm.
2. Observations
2.1. SOFIA/HAWC+
J1429-0028 was observed with the HAWC+ instrument (Dowell et al. 2013; Smith et al. 2014) onboard SOFIA on 2017 May 11 under the Cycle 5 program PID05_0087 (PI: Cooray). HAWC+ is a FIR imager and polarimeter with continuum bandpasses from 40 to 300 μm. We obtained the observation in band C at 89 μm (rest-frame 44 μm) with a bandwidth of ∼17 μm in the total-intensity OTF map configuration (Figure 1). The image has a field of view of 4
2 × 27 and a resolution (FWHM) of 7
8. The total effective on-source time is 939 s. The raw data was processed through the CRUSH pipeline v2.34-4 (Kovács 2008), and the final Level 3 data product was flux calibrated. The flux calibration error is about 10%. The resulting map has an rms noise level of ∼50 mJy beam−1. We extracted the photometry of the source using an 1-FWHM radius aperture and the resultant flux density is 748.4 ± 101.1 mJy. Using a lensing magnification factor of 9.7 ± 0.7 (Messias et al. 2014), the demagnified flux density at 89 μm is 77.2 ± 11.8 mJy. The minor contribution to the observed flux density from the foreground lensing galaxy is estimated to be at most 3%, which is within the quoted uncertainties.
Figure 1. Left: SOFIA/HAWC+ 89 μm detection of J1429-0028. The source is unresolved. The hatched circle shows the beam FWHM of 78. Right: three-color image of the gravitationally lensed system using HST F105W (blue), F160W (green), and Keck Ks (red) imaging data (Timmons et al. 2015).
Download figure:
Standard image High-resolution image2.2. Multiwavelength Photometry
Multiwavelength photometry from rest-frame UV to radio has been obtained from SDSS (u, g, r, i, z), VIKING (Z, Y, J, H, Ks), HST (F110W), Spitzer (3.6, 4.5 μm), WISE (3.4, 4.6, 12, 22 μm), Herschel (100, 160, 250, 350, 500 μm), IRAM-30 (1.2 mm), ALMA (1.28, 2.8 mm), and JVLA (7 GHz) (Messias et al. 2014). The de-blending of the foreground lens and background source is discussed in detail in Messias et al. (2014). We use the de-blended magnification-corrected photometry for the source to perform SED modeling (Table 1).
3. SED Models
Messias et al. (2014) performed SED fitting on J1429-0028 using the Multi-wavelength Analysis of Galaxy Physical Properties (magphys; da Cunha et al. 2008), which did not include AGN templates or radio emission. In this work, we employ two SED fitting codes, sed3fit (Berta et al. 2013) and the Code Investigating GALaxy Emission (cigale version v0.11; Noll et al. 2009; Roehlly et al. 2014), which include the additional components and rely on the energy balance technique, i.e., the energy of the absorbed starlight by dust is equal to the re-radiated energy by dust at infrared wavelengths.
sed3fit is a combination of the original magphys code and the Fritz et al. (2006) AGN library (updated by Feltre et al. 2012) which is effectively fit to data in a simultaneous three-component (stars, dust, AGN) model. Here we use the updated version of magphys for high-z IR models (da Cunha et al. 2015), which extends the parameter space to include properties that are more likely applicable to high-redshift SMGs. cigale includes up-to-date and customized SFH models, various stellar population synthesis models and IR dust re-emission templates, AGN templates, and radio emission, to constrain the AGN contribution in the IR via an energy balance manner similar to sed3fit. Both cigale and SED3FIT/MAGPHYS fit galaxy SEDs using a Bayesian approach to generate the posterior probability distribution function (PDF) for each parameter of interest, marginalized over all other parameters. The output values of the analyzed parameters in cigale are calculated as the weighted mean and standard deviation derived from the PDFs. SED3FIT/MAGPHYS outputs the 16%, 50% (median), 84% values of the PDFs. The two codes are the most widely used panchromatic SED fitting codes, and whether or not they derive consistent results is of interest and may have impact on how we interpret the data. It has been suggested that multiple SED fitting approaches should be used to cross-check the results (Hayward & Smith 2015).
In the subsections below, we describe the various models employed in sed3fit and cigale respectively. Since there are more than one SFH, dust attenuation law, and dust emission models that one can easily change in cigale, we will also test how different choices of models and parameters can affect the derived physical properties by varying one component at a time.
3.1. Stellar Component
We adopt the Chabrier (2003) IMF and the stellar population synthesis models of Bruzual & Charlot (2003, hereafter BC03) in both sed3fit and cigale. The original magphys code, which Messias et al. (2014) used, implements the updated version of BC03 models (CB07), where the new prescription by Marigo & Girardi (2007) for the thermally pulsating asymptotic giant branch evolution of low- and intermediate-mass stars is incorporated. CB07 models have been tested to result in 50%–80% lower stellar masses in HUDF galaxies than determined from the BC03 models (Bruzual 2007).
3.1.1. SED3FIT/MAGPHYS
The original magphys assumes an exponentially decreasing τ SFH, where τ is the star formation e-folding time. However, Messias et al. (2014) adopted the models in Rowlands et al. (2014), because the standard models have limitations to cover the physical parameter space for a DSFG. Messias et al. (2014) assumed both exponentially increasing and decreasing SFHs by distributing τ as a Gaussian between −1 and 1 Gyr−1. Bursts of star formation are superimposed at random times on the continuous SFH with a 50% probability of experiencing a burst in the last 2 Gyr. sed3fit and the high-z version of magphys assume an underlying continuous SFH with superimposed star formation bursts of random duration and amplitude. Each star formation burst can last between 30 and 300 Myr, and is set to occur at any random time in the previous 2 Gyr with a 75% probability. The priors on metallicity are a uniform distribution from 0.2 to 2 times solar.
3.1.2. CIGALE
In cigale, we first adopt the delayed-τ SFH model that rises in SFR up to a maximum, followed by an exponential decrease, which is motivated by high-redshift galaxies (Lee et al. 2010). Based on observations as well as galaxy evolution models for high-redshift galaxies (Lee et al. 2010; da Cunha et al. 2015), the SFHs overall should be initially rising with time before declining, although the real SFHs can be much more complicated (e.g., Behroozi et al. 2013; Pacifici et al. 2013; Simha et al. 2014; da Cunha et al. 2015). We also investigate the two-component τ-SFH model, which contains a young and an old stellar population (similar to the two components in sed3fit/magphys but in different forms), because real systems are likely to experience multiple episodes of star formation. There is evidence that the true stellar mass is likely better reproduced with the double SFHs (Michałowski et al. 2014), although the mass of the old stellar population is typically poorly constrained. Each stellar component is characterized by an exponentially declining SFH with two parameters: the stellar population age and the e-folding time. The two stellar populations are linked by the burst fraction (i.e., the mass fraction of the young stellar population). cigale also allows self-defined SFH as an input, but we do not attempt more complicated models. Possible values for metallicity are 0.0001, 0.0004, 0.004 (SMC), 0.008 (LMC), 0.02 (solar), and 0.05.
3.2. Dust Attenuation
3.2.1. SED3FIT/MAGPHYS
sed3fit/magphys uses the two-component dust attenuation model of Charlot & Fall (2000) with each component described by a power law. This model accounts for the fact that young stars in their birth clouds are more dust-attenuated than intermediate-age and old stars in the diffuse ISM.
3.2.2. CIGALE
Three dust attenuation models can be used in cigale: (1) a modified version of the Calzetti et al. (2000) attenuation law, (2) a single power-law model as defined in Charlot & Fall (2000) with a UV bump added, and (3) the two-component model of Charlot & Fall (2000) as in sed3fit/magphys. cigale enables both models (1) and (2) to allow differential attenuation in young and old stellar populations by using a reduction factor, which is defined as the ratio of the V-band attenuation of the old stellar population to the young population. We test all the dust attenuation laws and adopt the one that produces the best-fit SED.
3.3. Dust Emission
3.3.1. SED3FIT/MAGPHYS
The dust emission in sed3fit/magphys consists of four components: (1) mid-IR continuum from host dust, (2) warm dust in thermal equilibrium, (3) cold dust in thermal equilibrium, and (4) PAH empirical templates. The cold and warm dust components in thermal equilibrium emit as modified blackbodies with fixed emissivities. This dust model outputs the dust luminosity, dust mass, and dust temperatures for each component as well as a luminosity-weighted average dust temperature.
3.3.2. CIGALE
We first employ Dale et al. (2014) IR dust emission models (updated Dale & Helou 2002 models), which have been successfully applied to SMGs (Ma et al. 2015, 2016). The semi-empirical Dale et al. (2014) templates are parameterized by the power-law slope of the dust mass distribution over heating intensity α, dM/dU ∝ U−α. The dust emission templates are connected to the attenuated stellar population models by the dust luminosity, Ldust, which sets the basis of the energy balance.
We also test other physically motivated IR models. We use an updated version of the Draine & Li (2007) model that is implemented in cigale (Aniano et al. 2012; Ciesla et al. 2014). The Draine & Li (2007) model is motivated by the physical nature of dust (composition, geometry, and size distribution), and parameterized by the PAH mass fraction, minimum and maximum intensities of the interstellar radiation field, and the relative dust mass fraction heated by the diffuse radiation and the photodissociation regions. The normalization of the model to the photometric data yields the dust mass, Mdust, in addition to Ldust.
3.4. AGN Emission
The AGN emission in both SED codes is based on the AGN templates from Fritz et al. (2006), which consist of two components. One is the isotropic emission of the central AGN and the other is an improved model of the emission from the dusty torus heated by the central engine. Part of the central emission is absorbed by the dusty torus and re-emitted at longer wavelengths. sed3fit also includes the emission of the accretion disk in an updated library by Feltre et al. (2012). The full library comprises 2376 models with each one computed at 10 different line-of-sight angles. sed3fit provides a set of 10 AGN templates spanning the whole range of colors covered by the full library. A few parameters are specified to characterize the dust torus, including inner/outer radius, radial and angular dust distributions, angular opening angle, etc. Type 1 (unobscured), Type 2 (obscured), as well as intermediate-type AGN can be approximately described by those parameters. The most important parameter is the fractional contribution of AGN to the total IR luminosity (LIR,total = LStarburst + LAGN). We refer the reader to Berta et al. (2013) and Małek et al. (2017) for the detailed implementation and allowed parameter ranges in each code.
3.5. Radio Emission
The SED fit is also extended to radio wavelengths based on the well-established FIR-to-radio correlation, q, for star-forming galaxies. A tight FIR-to-radio correlation has been observed for star-forming galaxies over five orders of magnitude in galaxy luminosity (e.g., Yun et al. 2001). Radio-loud AGN, on the contrary, exhibit elevated radio emission that would be at least a factor of ∼2 higher (Rush et al. 1996; Morić et al. 2010). This radio emission model can be used to test whether or not the observed radio photometric data is consistent with the FIR-to-radio correlation.
sed3fit/high-z magphys computes the radio emission as the sum of a thermal (free–free emission) and a nonthermal component with fixed power-law slopes. A Gaussian prior distribution is assumed for the coefficient, q, with the mean at 2.34 (local value; Yun et al. 2001; Bell 2003) and 1σ value of 0.25 (e.g., Ivison et al. 2010). In cigale, a nonthermal (synchrotron emission from galaxies) component is added and the thermal emission is handled by the nebular module. The default coefficient of the FIR/radio correlation is 2.58 and the spectral index of the synchrotron power-law emission is 0.8 (same as in sed3fit). We test different prior values of 2.2–2.6 for q as it does not necessarily follow the local relation and may evolve with redshift (Delhaize et al. 2017).
Table 1. Input Photometry for SED Fitting with MAGPHYS, SED3FIT, and CIGALE
| Survey/Facility | Filters | Flux Density (mJy) |
|---|---|---|
| SDSS | u | <1.36 × 10−3 |
| g | <1.67 × 10−3 | |
| r | <2.72 × 10−3 | |
| i | <2.93 × 10−3 | |
| z | (2.77 ± 1.69) × 10−3 | |
| VIKING | Z | <3.12 × 10−3 |
| Y | (2.72 ± 1.51) × 10−3 | |
| J | (4.15 ± 1.62) × 10−3 | |
| H | (7.08 ± 1.60) × 10−3 | |
| Ks | (1.88 ± 0.37) × 10−2 | |
| HST | F110W | (2.69 ± 0.72) × 10−3 |
| Spitzer | 3.6 μm | (5.17 ± 1.18) × 10−2 |
| 4.5 μm | (5.69 ± 1.30) × 10−2 | |
| WISE | 3.4 μm | (4.57 ± 1.07) × 10−2 |
| 4.6 μm | (5.53 ± 1.28) × 10−2 | |
| 12 μm | (4.68 ± 1.24) × 10−1 | |
| 22 μm | (5.80 ± 1.48) × 10−1 | |
| SOFIA | 89 μm | 77.2 ± 11.8 |
| Herschel | 100 μm | 81.0 ± 15.9 |
| 160 μm | 114.8 ± 22.3 | |
| 250 μm | 76.7 ± 14.7 | |
| 350 μm | 46.3 ± 8.9 | |
| 500 μm | 22.6 ± 4.4 | |
| IRAM-30 | 1.2 mm | 0.95 ± 0.26 |
| ALMA | 1.28 mm | (5.43 ± 1.02) × 10−1 |
| 2.8 mm | (5.03 ± 0.97) × 10−2 | |
| JVLA | 7 GHz | (1.75 ± 0.23) × 10−1 |
Note. For the non-detections, the flux density upper limits are given at the 3σ level. We note that the background source flux densities listed in Table 4 of Messias et al. (2014) are the best-fit model flux densities derived from magphys rather than the input photometry.
Download table as: ASCIITypeset image
4. Results and Discussion
Table 2 lists the properties derived from the PDF analysis of sed3fit and cigale, including the instantaneous SFR, SFRs averaged over the last 10 and 100 Myr, M*, stellar population age (mass-weighted age), V-band extinction, IR luminosity, dust mass, dust temperature, and the AGN fraction in the IR luminosity. sed3fit also calculates the luminosity-averaged dust temperature as an output. Using the SFR and M*, we also derive the specific SFR (sSFR), which is a critical diagnostic for understanding the star formation mode of this galaxy.
Table 2. Physical Properties Derived from MAGPHYS, SED3FIT, and CIGALE
| Models and Properties | A. MAGPHYSa | B. High-z MAGPHYS | C. SED3FIT | D. CIGALE | E. CIGALE | F. CIGALE |
|---|---|---|---|---|---|---|
| Stellar population | CB07 | BC03 | BC03 | BC03 | BC03 | BC03 |
| SFH | modified τ SFH | delayed-τ SFH | delayed-τ SFH | delayed-τ SFH | double-τ SFH | double-τ SFH |
| + random bursts | + random bursts | + random bursts | (fburst = 0–1) | (fburst ≤ 0.1) | ||
| Dust attenuation | Charlot+2000 | Charlot+2000 | Charlot+2000 | power law | power law | power law |
| Dust emission | Rowlands+2014 | daCunha+2015 | daCunha+2015 | Draine+2014 | Draine+2014 | Draine2014 |
| AGN emission | ⋯ | ⋯ | Feltre+2012 | Fritz+2006 | Fritz+2006 | Fritz+2006 |
| SFR (M⊙ yr−1) | ⋯ | ⋯ |
|
714 ± 75 | 582 ± 29 | 516 ± 30 |
| SFR10 Myr (M⊙ yr−1) | ⋯ | ⋯ |
|
662 ± 56 | 582 ± 29 | 516 ± 30 |
| SFR100 Myr (M⊙ yr−1) |
|
|
|
361 ± 36 | 296 ± 21 | 266 ± 13 |
| M* (1010 M⊙) |
|
|
|
2.1 ± 0.5 | 3.4 ± 1.1 | 15.1 ± 1.2 |
| sSFR (Gyr−1) | 3.0 ± 1.6 | 0.51 ± 0.39 | 0.21 ± 0.17 | 17.2 ± 4.4 | 8.7 ± 2.9 | 1.8 ± 0.2 |
| ageM (Myr) |
|
|
|
<160 | <600 | 1991 ± 314 |
| AV | ⋯ |
|
|
4.3 ± 0.2 | 4.2 ± 0.2 | 4.1 ± 0.2 |
| Ldust (1012 L⊙) |
|
|
|
4.59 ± 0.23 | 4.41 ± 0.22 | 4.13 ± 0.21 |
| Mdust (108 M⊙) |
|
|
|
6.4 ± 0.4 | 6.3 ± 0.5 | 6.2 ± 0.6 |
| Tdust (K) | ⋯ |
|
|
⋯ | ⋯ | ⋯ |
 (%) |
0 | 0 |
|
<0.90 | <0.72 | <0.75 |
|
0.50 | 0.43 | 0.56 | 0.63 | 0.66 | 0.87 |
Notes. SFR is the current star formation rate. SFR10 Myr is the SFR averaged over the last 10 Myr. SFR100 Myr is the SFR averaged over the last 100 Myr. ageM is the mass-weighted stellar population age. sSFR is based on SFR100 Myr. Tdust is the weighted dust temperature of the birth cloud and ISM dust temperatures. 
from CIGALE is given at 3 σ. 
is the reduced χ2 from the SED fitting.
Download table as: ASCIITypeset image
We have checked the PDFs to make sure that the way sed3fit and cigale calculate the average values of parameters, either median (sed3fit) or weighted mean (cigale), does not make much of a difference (i.e., within 1σ uncertainty).
4.1. Star Formation History
A single stellar population with a delayed-τ SFH model from cigale (Model D) yields the highest SFR and the lowest M* among the model configurations. The stellar population age is poorly constrained and the upper limit is consistent with the timescale of a recent starburst during the SMG phase (e.g., Tacconi et al. 2006; Ivison et al. 2011).
The double SFH model from cigale (Model E) produces a higher stellar mass and a lower current SFR. The initial input burst fraction is set as a free parameter, i.e., 0–1. The SEDs favor an almost pure young burst (burst fraction ∼0.8) for ∼100 Myr, suggesting that the stellar mass in the entire system was built up rapidly in the last 100 Myr. The prodigious star formation is most likely triggered by major mergers, based on the source-plane morphology and dynamical mass estimates (Messias et al. 2014). A high burst fraction further supports this scenario. However, this burst fraction derived from SED modeling is significantly higher than those from merger simulations (<50%; Hopkins et al. 2013). Smith & Hayward (2015) use merger simulations to test how well magphys is able to recover SFHs, where the true SFH is known. They fail to obtain a good estimate of the SFH for the merger simulations despite being able to get a reasonable fit to the synthetic photometry. Likewise, we test cigale using merger simulations to check whether the burst fraction can be reasonably recovered. The comparison suggests that cigale tends to overestimate the burst fraction likely due to the similar effect of the "outshining" problem of the young stellar populations over older stellar populations (e.g., Reddy et al. 2012). If we restrict the burst fraction to be ≤0.1 (Model F, typical burst fraction suggested by merger simulations; Hopkins et al. 2013), the resultant stellar mass increases by a factor of ∼5 without much change in the SFR or the quality of the fit, and is in agreement with the results in Messias et al. (2014).
The sed3fit model (Model C) with delayed-τ SFH and random bursts yields the highest stellar mass and lowest SFR. The blue end of the best-fit SED from SED3FIT shown in Figure 2 is much higher than magphys and cigale best-fit SEDs due to the unconstraining upper limits, and is dominated by AGN emission, although there is no evidence of a Type 1 AGN in the system (Timmons et al. 2015). This significantly affects the stellar SED and the mass-to-light ratio, resulting in a much higher stellar mass and lower SFR.
Figure 2. Best-fit SED from cigale (black; Model E) compared to the best fits from sed3fit (green; Model C) and magphys (orange; Model A). Different cigale SED components are color coded. The unattenuated stellar emission from sed3fit (purple) and magphys (brown) is also displayed in comparison with the one from cigale (blue). The data points from left to right are SDSS ugri and VIKING Z 3σ upper limits, SDSS z, VIKING Y, HST/F110W, VIKING J, H, Ks, WISE1, Spitzer/IRAC 3.6, 4.5 μm, WISE2, WISE3, WISE4, SOFIA/HAWC+ 89 μm (filled red square), Herschel/PACS 100, 160 μm, Herschel/SPIRE 250, 350, 500 μm, IRAM-30 1.2 mm , ALMA 1.28, 2.8 mm, and JVLA 7 GHz. The inset shows the normalized probability distribution function of the AGN fraction in the total IR luminosity. The bottom panel shows the relative residual fluxes of the cigale SED fit.
Download figure:
Standard image High-resolution imageWardlow et al. (2011) test the effect of SFHs on stellar masses and ages for SMGs with 17-band photometry. They cannot reliably distinguish the different SFHs and estimate that this results in an additional factor of ∼5 uncertainty in the mass-to-light ratios, and therefore, the derived stellar masses. Even in the case of extensive photometric coverage such as J1429-0028, we cannot draw definitive conclusions about the SFH (stellar ages, burst fractions) solely based on the SED modeling.
4.2. Dusty System
All the SED fits indicate high dust extinction, AV ∼ 4–7, and ultra-luminous dust emission >4 × 1012 L⊙. The dust luminosity is the best constrained parameter with least uncertainty; the dust luminosities derived from different models and SED codes are all consistent with each other. The best-fit FIR SEDs from magphys, sed3fit, and cigale are similar while the MIR slopes are different due to lack of data. The derived dust masses from different models and SED codes are in the range of ∼3–6 × 108 M⊙. The cigale FIR dust emission is based on the updated Draine & Li (2007) model, and we adopt in sed3fit the high-z version of the dust emission model from magphys (da Cunha et al. 2015), because the original (low-z) library does not properly constrain the physical properties for this kind of galaxy. The average dust temperature is about 40 K, which is consistent with that of star formation dominated DSFGs.
4.3. Negligible AGN Contribution
Previous nebular line diagnostics show no indication of an unobscured Type 1 AGN in J1429-0028 (Timmons et al. 2015). However, an obscured Type 2 AGN could potentially contribute to the observed very high LIR. We set the input AGN type as a free parameter, i.e., Type 1, Type 2, and intermediate-type AGN. Regardless of the different combinations of component models and SED codes, the AGN contribution to the IR luminosity is negligible (3σ upper limit 
). We note, however, that the best-fit SED in the blue end from sed3fit (green) is significantly higher than the ones from magphys (gold) and cigale (black) due to an unexpected Type 1 AGN component dominating the UV spectrum and also affecting the optical-NIR part. This is simply caused by the unconstraining upper limits in the blue end. We cannot rely on the SED fitting alone to infer AGN activity in this case (more discussion in Section 4.4).
The JVLA radio data at 7 GHz is also consistent with the FIR-to-radio correlation for star-forming galaxies. Radio-loud AGN would push the radio emission at least a factor of 2 higher. We note that in this particular galaxy the new SOFIA data at 89 μm does not make much of a difference given the presence of PACS 100 μm. A 53 μm band is preferred for constraining the slope of the MIR SED. For sources that do not have PACS data, the 89 μm band can be critical.
In Figure 3 we compare the AGN fraction and IR luminosity of J1429-0028 with those of Chandra X-ray selected AGN in the Boötes legacy field (Jannuzi & Dey 1999; Murray et al. 2005) with MIR and FIR counterparts detected by Spitzer (Ashby et al. 2009) and Herschel (Oliver et al. 2012), which represents the largest AGN sample with a multiwavelength SED fitting analysis (Brown et al. 2018). Figure 3 shows the general trend of higher AGN fraction with increasing IR luminosity. The colored symbols denote the X-ray detected AGN in different redshift bins. The gray X's are the X-ray non-detections in the survey. J1429-0028 has an intrinsic IR luminosity that is higher than most of the sources at similar redshifts, comparable to the average LIR at higher redshifts. The average (median) AGN fraction in the X-ray selected sample in the redshift bin 0.8 < z < 1.2 is ∼10%, while the 3σ upper limit of J1429-0028 is significantly below the majority of the X-ray detected AGNs.
Figure 3. AGN fraction 
vs. LIR. The colored symbols are the X-ray detected AGN in different redshift bins by Brown et al. (2018), which represents the largest AGN sample with a multiwavelength SED fitting analysis. The gray X's are the X-ray non-detections. There is a general trend of higher AGN fractions with increasing LIR. The downward arrow denotes the 
upper limit of J1429-0028.
Download figure:
Standard image High-resolution image
4.4. Systematic Uncertainties in SED Modeling
Given that the AGN fraction is negligible in the IR, we now compare the other critical physical properties with and without the inclusion of the AGN component. The previous best effort on SED modeling of this source (Messias et al. 2014) was based on the original magphys, which employs the CB07 stellar population synthesis models rather than the BC03 models (used in the new high-z version of magphys). The inferred stellar masses from the BC03 models are about 1.25–1.5 times higher than those from the CB07 models (Bruzual 2007). The high-z version of magphys (Model B) results in a factor of ∼3 higher stellar mass for J1429-0028. The additional factor of ∼2 may be due to a combination of different choices of SFHs and IR dust emission models. We further compare the results from the high-z magphys (Model B) and sed3fit (Model C) for which the only difference is the inclusion of the AGN component. The physical parameters are consistent with each other within 1σ uncertainties. However, we stress that one should use multiple AGN indicators otherwise an unexpected AGN component can appear in the SED in the case of loose photometric data constraints in the blue end, which affects the derived stellar mass, SFR, etc. Deriving reliable stellar masses from SED fitting for DSFGs has been extremely challenging. We have demonstrated that even in the case of extensive photometric coverage from UV to radio for a single galaxy, the systematic uncertainties can be as high as a factor of ∼10. Therefore, one should explicitly state the models involved in SED modeling, especially the SFHs, and convert to similar assumptions if possible whenever we compare stellar mass measurements.
5. Summary and Conclusion
We have presented the SOFIA/HAWC+ 89 μm detection of the gravitationally lensed DSFG at z = 1.027. We conduct a detailed analysis of multiwavelength SED modeling from rest-frame UV to radio, including the new SOFIA/HAWC+ data, and test how different combinations of component models or parameters and SED codes affect the derived physical properties. A significant AGN contribution to the IR luminosity has been ruled out at a high confidence level regardless of model choices. A significant fraction of the total stellar mass ∼1011 M⊙ in the galaxy is built up in rapid starburst in the last ∼100 Myr with an instantaneous SFR of >500 M⊙ yr−1, likely triggered by a major merger. The inferred sSFR places J1429-0028 above the star-forming main sequence at z ∼ 1, further indicating the strong starbursting nature of the galaxy.
We stress that even in the case of extensive photometric coverage, the uncertainty in stellar mass of a DSFG can be as high as a factor of ∼10 and is dominated by the assumed SFH. Therefore it is of crucial importance to state explicitly the model assumptions associated with the derived stellar masses. One should also combine other AGN indicators with SED fitting including an AGN component to infer the AGN activity.
SOFIA/HAWC+ is currently the only facility that covers the wavelength regime 50–100 μm and such observations become critical for targets that do not have coverage in that spectral region. Although this particular galaxy does not contain an energetically important AGN, this paper demonstrates the potential of using SOFIA/HAWC+ to constrain the AGN fraction in Herschel-selected galaxies due to the lack of MIR data.
We thank the anonymous referee for the very constructive comments that have greatly improved the manuscript. We thank Christopher Hayward and Desika Narayanan for helpful discussions. This work is based 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 project PID05_0087. J.L.W. acknowledges the support of an STFC Ernest Rutherford Fellowship under grant number ST/P004784/1. D.R. acknowledges support from the National Science Foundation under grant number AST-1614213. This research has made use of NASA's Astrophysics Data System.