A NuSTAR SURVEY OF NEARBY ULTRALUMINOUS INFRARED GALAXIES

During its two-year baseline mission, NuSTAR observed a sample of nine of the nearest (z < 0.078) ULIRGs, out of the total sample of ∼25 ULIRGs within that volume, based on the selection of Sanders et al. (2003). The sample was split into a "deep" sample of five ULIRGs observed for >50 ks, and a "shallow" sample of four objects observed for ∼25 ks. The deep sample consists of the four ULIRGs at z < 0.078 that are brightest at 60 μm. This selection by low redshift and high infrared brightness was done to maximize the likelihood of obtaining high-quality NuSTAR data. The Superantennae was added to the deep sample because previous observations at >10 keV suggested the presence of a Compton-thick AGN. All five deep survey sources have previous observations above 10 keV by BeppoSAX or Suzaku PIN; three have reported detections. In addition, the deep sample targets have simultaneous soft X-ray coverage from either Chandra (Mrk 231; Teng et al. 2014) or XMM-Newton (the other four). Table 1 lists the targets and exposure times for the entire sample. The shallow sample targeted nearby ULIRGs showing AGN signatures in their optical spectra. Thus, the shallow sample, like the deep sample, is biased toward high detection rates at X-ray energies.

The targets in our sample have been very well studied across the electromagnetic spectrum. In this section, we note the information at other wavelengths and past X-ray observations that are relevant to our study.

2.2.1. The Deep Survey

2.2.2. The Shallow Survey

The XMM-Newton data were obtained using the EPIC array in full window imaging mode. For these observations, the medium optical blocking filter was applied. The data were reduced using XMM-Newton Science analysis Software version 13.5.0. The most up-to-date calibration files, as of 2014 May, were used for the reduction. We followed the XMM-Newton ABC guide²³ to extract images and spectra for our Deep Survey targets. In particular, during the data reduction process, portions of data with high background flares were removed. The Mrk 273 data were highly affected by these flares; thus, the calibrated data only contain ∼17% of the original total exposure. All of the Deep Survey targets appear to be point sources. In order to ensure that we are probing the same spatial scale as the NuSTAR data, we used the same source extraction regions as NuSTAR data (see below). The background spectra are extracted from nearby source-free areas on the same chip using extraction regions of the same size. The extracted spectra were binned to 15–50 counts per bin, depending on the source count rate, such that χ² statistics may be used.

For IRAS 13120–5453, archived Chandra data were used to extract a low-energy X-ray spectrum (PI: Sanders). These data were reduced using CIAO 4.5 with CALDB version 4.5.6. The standard calibration procedures were followed for reducing ACIS-S data in VFAINT mode using the chandra_repro script.²⁴ As with the XMM-Newton data, the same NuSTAR source region was used to extract the source spectrum. The background spectrum was extracted using a region of the same size in a nearby source-free area. The spectrum was binned to at least 15 counts per bin so that χ² statistics can be applied.

The NuSTAR observations were reduced using the NuSTAR Data Analysis Software (NuSTARDAS) that is part of HEASoft version 6.15.1, with NuSTAR calibration database version 20131223. The script nupipeline was used to produce calibrated event files for each of the two focal plane modules (FPMA and FPMB; Harrison et al. 2013). The good time intervals of these events are listed in Table 1. All spectra were binned such that χ² statistics can be used.

The E < 20 keV NuSTAR background is spatially non-uniform over FPMA and FPMB, as a result of stray light being incompletely blocked by the aperture stop. To correct for this aperture background and the instrumental background, we followed the procedure in Wik et al. (2014) to create simulated total background events for each source. The simulated backgrounds were scaled for exposure time, size of the extraction region, and the response of individual chips for each observation. The simulated backgrounds were used to create background-subtracted images for photometry in Section 4 and background spectra for the spectral analysis of the faint sources (IRAS 13120–5453, the Superantennae, and Arp 220) in Section 5. We conservatively estimate that the broadband systematic uncertainties in the derived background spectra are ≲5%.

With the exception of Mrk 273, source spectra were extracted using circular apertures with 1' radii. For Mrk 273, due to the projected vicinity of a background source (Mrk 273X; Xia et al. 2002), the source spectrum was extracted using a circle with 08 radius. Spectral analysis was performed using HEASoft version 6.15.1. For objects with multiple observations, the spectra were co-added using the FTOOL addascaspec. When modeling the spectra, an additional constant factor, typically of the order of a few per cent, is applied to account for the cross-normalization between FPMA and FPMB, and between FPMA and XMM-Newton EPIC-pn. These cross-normalization constants were allowed to vary for IRAS 05189–252 and Mrk 273 since these are the brightest sources, where more degrees of fit are possible; for the fainter IRAS 13120–5453, Arp 220, and the Superantennae they were held fixed. The cross-normalization values used were current at the time of the modeling. Subsequently, cross-normalization values have been published by Madsen et al. (2015) in the NuSTAR calibration paper. The values we used differ by only a few per cent from those of Madsen et al. (2015), and the difference has negligible impact on our results.

We assumed the abundances of Wilms et al. (2000) and the photoelectric cross sections of Verner et al. (1996) in our spectral modeling with XSPEC. The assumed column densities due to Galactic absorption (${N}_{{\rm{H,Gal}}}$) are given in Table 1. All errors quoted in this paper are at the 90% confidence level (Δχ² = 2.706 for a single parameter).

IRAS 05189–2524 is the X-ray brightest ULIRG in our sample. Observed in two NuSTAR pointings separated by about eight months, the second of which was divided into two data sets, IRAS 05189–2524 has shown minor variability between these epochs. In Figure 2, we show the light curve for our NuSTAR observations. The average count rate in the second and third observations changed by ∼20% relative to the average count rate (∼0.1 counts per second) in the first observation. However, we note that this variation is smaller than the standard deviation (∼0.03 counts per second or ∼30%) derived when data points from all three observations are combined. Therefore, the variability is not statistically significant. Emission lines at 6.4 and 6.8 keV were detected, with a Δχ² of 65 for four degrees of freedom.

Zoom In Zoom Out Reset image size

Although the light curve shows apparent minor variability, there is no obvious variation in the spectral shape between NuSTAR observations. In particular, we do not see the drop in the 2–10 keV emission measured by Teng et al. (2009) from their Suzaku data; the source appears to have reverted back to its previous "high" state. Therefore, we combined the spectra from all three NuSTAR data sets for our broadband modeling.

Modeling multiple epochs of historic X-ray data, Teng et al. (2009) found that the best-fit model to explain the sudden change in spectral shape is an increase in the line-of-sight column from two partial covering absorbers. Following this result, we fit the XMM-Newton plus NuSTAR broadband spectrum with a double partial covering model for the AGN component. Following the best-fit model for Mrk 231 in Teng et al. (2014), we also included a MEKAL and a cutoff power law (Γ = 1.1 with cutoff energy at 10 keV) for the non-thermal emission from high-mass X-ray binaries (HMXBs) to account for the ∼80 M_⊙ yr⁻¹ starburst, estimated from the infrared luminosity. The photon index for the AGN, Γ = 2.51 ± 0.02, is very steep, but is not far from the top of the range (e.g., 1.5 < Γ < 2.2; Nandra & Pounds 1994; Reeves & Turner 2000) observed in other AGNs. The two partial covering absorbers have N_H of 5.2 ± 0.2 × 10²² cm⁻² and ${9.3}_{-0.7}^{+1.0}\times {10}^{22}$ cm⁻² with 98 ± 0.2% and 74${}_{-1.6}^{+1.2}$% covering fractions, respectively. This model is shown (${\chi }_{\nu }^{2}\;\sim \;$1.07), along with the spectrum, in Figure 3. The derived Γ is steeper than Γ_eff estimated from the photometry because the Γ_eff calculation did not account for the flat power-law contribution from the HMXBs (Γ_HMXB = 1.1).

Zoom In Zoom Out Reset image size

Given the high 2–10 keV flux levels, we do not expect this source to be highly obscured. Both the MYTorus and BNTorus models poorly describe the broadband spectrum (${\chi }_{\nu }^{2}\sim 1.8$ for both).

From our modeling, we find that the AGN in IRAS 05189–2524 is currently in a Compton-thin state. Its intrinsic 2–10 keV luminosity is $3.7\times {10}^{43}$ erg s⁻¹, about 0.8% of the bolometric luminosity. Our present result is consistent with those from the multiple-epoch fitting by Teng et al. (2009). The large drop in the observed 2–10 keV flux in the 2006 Suzaku X-ray Imaging Spectrometer (XIS) data is likely due to an intervening absorber, since a change by a factor of 30 in intrinsic flux is rare in AGNs (Gibson & Brandt 2012). The incomplete nature of the time series data makes it impossible to precisely determine the timescale of variability. If it is of the order of years, as is consistent with the data, then the Compton-thick absorber responsible for the Suzaku variability must be within a few parsecs of the nucleus. However, a more extreme case of intrinsic variability (a factor of ∼260) was observed in the narrow-line Seyfert 1 galaxy PHL 1092 (Miniutti et al. 2012) so we cannot absolutely rule out the possibility of strongly varying intrinsic flux. Indeed, our best-fit model, even in the high state, should not have been detectable by Suzaku PIN in the 2006 observation, which speaks to the much greater sensitivity above 10 keV of NuSTAR compared to Suzaku PIN.

We first fit the broadband spectrum of IRAS 13120–5453 with a simple power law modified only by Galactic absorption and a MEKAL model for the starburst component. Because of the relatively high Galactic column density (N_H,Gal ∼ 4.6 × 10²¹ cm⁻²), it is difficult for the MEKAL model to constrain the line components below 2 keV. The best-fit power law requires Γ ∼ −2 to fit the shape of the >2 keV spectrum, implying a highly obscured AGN continuum. Additionally, the spectrum shows three strong emission lines at 1.86, 3.40, and 6.78 keV which correspond to Si xiii, Ar xviii, and Fe xxv, respectively. The Δχ² for the iron line is 8.35 for two degrees of freedom. The strong lines of the α-elements Si (EW ∼ 0.20 keV) and Ar (EW ∼ 0.53 keV) may be an indication of a strong starburst. Since the lines are stronger than the predictions of the MEKAL plasma model, if they are starburst in origin then a more complex starburst model is needed, for example with multiple temperatures, or highly non-solar abundances.

To constrain the intrinsic absorption of the AGN continuum, we added an absorption component and fixed Γ at the canonical value of 1.8. This new model requires that the intrinsic absorber have a column density of ∼4 × 10²⁴ cm⁻², implying the AGN is Compton-thick. For both torus models, we fixed the torus inclination angle at 85°, since the optical data suggest the AGN is Type 2. Since we do not have an independent measure of the SFR in IRAS 13120–5453 as we did for Mrk 231 and IRAS 05189–2524, we allowed the hot gas temperature and the normalizations of the HMXB model components to vary. We also included Gaussian components to model the lines at 1.9, 3.4, and 6.8 keV.

The MYTorus model seems to fit the spectrum well (${\chi }_{\nu }^{2}\sim 0.80$); however, the model cannot constrain the error of Γ within the bounds of the MYTorus model (1.4 < Γ < 2.5). For the best-fit model, we fixed Γ at 1.8, resulting in a column density of ${3.1}_{-1.3}^{+1.2}\times {10}^{24}$ cm⁻², consistent with the assertion by Iwasawa et al. (2011) that IRAS 13120–5453 is Compton-thick based on the strength of the Fe line. The 2–10 keV absorption-corrected luminosity for IRAS 13120–5453 is 1.25 × 10⁴³ erg s⁻¹. The starburst component is absorbed by a column of ${5.6}_{-5.6}^{+14.3}\times {10}^{21}$ cm⁻². The thermal component has a temperature of 0.56${}_{-0.31}^{+0.16}$ keV and a 0.5–2 keV luminosity of $1.8\times {10}^{41}$ erg s⁻¹. The non-thermal HMXB component has a 0.5–8 keV luminosity of 4.5 × 10⁴¹ erg s⁻¹. The luminosities of both the thermal and non-thermal components are consistent with a SFR of ∼170 M_⊙ yr⁻¹, based on the relations of Mineo et al. (2012a, 2012b). Although high, this SFR is within the range observed for ULIRGs. For comparison, Mrk 231 has a SFR of ∼140 M_⊙ yr⁻¹ (Rupke & Veilleux 2011). The absorption-corrected 2–10 keV luminosity of IRAS 13120–5454 is ∼0.67% of its AGN bolometric luminosity. The best-fit MYTorus spectrum is shown in Figure 4. Given the data quality, the BNTorus model places poor constraints on the torus opening angle and Γ. However, the values for the column density, intrinsic X-ray luminosity, and the starburst components are consistent with the results from MYTorus.

Zoom In Zoom Out Reset image size

Mrk 273 is the only source in a past Suzaku survey of ULIRGs that was detected above 10 keV (Teng et al. 2009) by PIN. With only a marginal detection (1.8σ), the Suzaku data required a double partial covering model and implied that at least one of the partial covering absorbers is Compton-thick (N_H ∼ 1.6 ×  10²⁴ cm⁻²). A variability analysis using the Suzaku and historic X-ray data suggests that the variations in spectral shape between 2 and 10 keV are due to changes in the column density or the fraction of the partial covering.

Mrk 273 is well detected by NuSTAR up to ∼30 keV, above which the background dominates. The MYTorus and BNTorus models fit the data nearly equally well. As with the Mrk 231 analysis (Teng et al. 2014), we included model components that account for HMXB and thermal emission for a starburst that is forming stars at a rate of ∼160 M_⊙ yr⁻¹ (Veilleux et al. 2009a). With the torus inclination angle fixed at 85°, the best-fit MYTorus model (${\chi }_{\nu }^{2}=0.91$) suggests that the direct intrinsic AGN emission (or the zeroth order emission) has ${\rm{\Gamma }}={1.43}_{-0.03}^{+0.17},$ and the global N_H is (4.4 ± 0.1) × 10²³ cm⁻². The leaky fraction of the absorber is only ${3.1}_{-1.8}^{+2.4}$%. This best-fit model and the data are shown in Figure 5.

Zoom In Zoom Out Reset image size

The BNTorus model gives similar results (${\chi }_{\nu }^{2}=0.90$). With the torus inclination angle fixed at 87°, Γ is constrained to be ${1.29}_{-0.17}^{+0.20}$ and the line-of-sight N_H is ${3.5}_{-0.8}^{+1.1}\times {10}^{23}$ cm⁻². The opening angle of the torus is at least ∼43°. The BNTorus model infers a leaky fraction of ${4.5}_{-2.7}^{+4.3}$%. These parameters are consistent with those derived using the MYTorus model.

Teng et al. (2009) concluded that the spectral variability seen in Mrk 273 was due to changes in the column density. By fitting the multiple-epoch data together and assuming a common model, they showed that the older measurements of the column density with Chandra and XMM-Newton were a factor of 2–4 lower than the Suzaku measurement. Our NuSTAR analysis is consistent with this result, as the column density we derived using the torus models is a factor of ∼4 lower than the Suzaku measurement. Similar to our Mrk 231 results (Teng et al. 2014), the torus models favor a relatively flat power-law intrinsic photon index for an AGN. The intrinsic 2–10 keV luminosity of Mrk 273 is ∼8.6 $\times {10}^{42}$ erg s⁻¹, representing ∼0.3% of the AGN bolometric luminosity.

Although it is the nearest ULIRG, Arp 220 is the faintest source at NuSTAR energies in our Deep Survey. This source was observed by Suzaku in 2006, but was undetected above 10 keV (Teng et al. 2009). The Suzaku spectrum, with a lack of detection above 10 keV and an ionized Fe line with large EW (∼0.42 keV), suggests that the AGN is very heavily obscured. The direct emission from the AGN is completely obscured by a high column density, leaving behind a purely reflected spectrum. Since the NuSTAR flux between 15 and 40 keV is 35 times lower than the upper limit derived from the Suzaku observations, our NuSTAR observations can more tightly constrain the X-ray properties of Arp 220.

When modeling the new broadband X-ray spectrum, we first revisited the ionized reflection model favored by Teng et al. (2009), which did not include an HMXB component. The MEKAL plus ionized reflection model (reflionx; Ross & Fabian 2005) is well fit to the broadband data (${\chi }_{\nu }^{2}\sim 1.25$). All model parameters are consistent with those derived from the Suzaku XIS data alone. To better model the shape of the spectrum below 2 keV, we added a second MEKAL component. The best-fit hot gas temperatures are ${0.10}_{-0.10}^{+0.02}$ and ${0.50}_{-0.25}^{+0.20}$ keV. Both these temperatures are consistent with those observed in ULIRGs (e.g., Franceschini et al. 2003; Ptak et al. 2003; Teng et al. 2005, 2009; Teng & Veilleux 2010). The underlying reflected power law has ${\rm{\Gamma }}={1.76}_{-0.32}^{+0.22},$ assuming the input ionization parameter ξ is 10³ erg cm s⁻¹. The reflected 2–10 keV luminosity is ∼9.0 × 10⁴⁰ erg s⁻¹. As the detailed modeling of Mrk 231 by Teng et al. (2014) has shown, it is important to constrain the HMXB component that also contributes to the X-ray spectrum. When an HMXB component (SFR = 200 M_⊙ yr⁻¹ assuming the starburst infrared luminosity derived by Veilleux et al. 2009a) is added to the model, the reflection component is no longer required.

We attempted to use the BNTorus model to constrain the properties of the obscuring torus. However, due to the poor photon statistics above 10 keV, the model does a poor job in deriving robust values and errors for parameters that characterize the torus. The model suggests a line-of-sight column density of at least $5.3\times {10}^{24}$ cm⁻² with the best-fit value tending toward $\gt {10}^{25}$ cm⁻². Similarly, the MYTorus model has difficulty constraining the properties of the torus. Using this model, the column density is at least 1.2 × 10²⁴ cm⁻², with the nominal value hitting the upper bound of the model at 10²⁵ cm⁻². In this case, it is not possible to measure the intrinsic X-ray luminosity of the AGN. Therefore, if an AGN is present in Arp 220, it is highly Compton-thick and the N_H cannot be constrained with the NuSTAR >10 keV data. This result is consistent with measurements by two groups that find that the western nucleus of Arp 220 is embedded in a column of 1.3 × 10²⁵ cm⁻² (Downes & Eckart 2007) to 1.5 × 10²⁵ cm⁻² (Scoville et al. 2014).

It is possible that the observed X-ray spectrum of Arp 220 does not have an AGN component at all. If the column density is $\gt {10}^{25}$ cm⁻², then no direct emission can be detected if the absorber is not "leaky." What the observed global X-ray spectrum represents is simply the thermal and non-thermal emission from the major starburst. Therefore, we also modeled the Arp 220 spectrum without an AGN component. The data are well fit by a two-temperature MEKAL plus a strong ionized Fe K line. The iron line likely originates from the strong starburst. A strong bremsstrahlung component would be expected to accompany the 6.7 keV emission. Thus, we also included a redshifted bremsstrahlung component in the two-temperature MEKAL model. The MEKAL luminosities were fixed such that they are consistent with the SFR based on the relation of Mineo et al. (2012b). This best-fit model (${\chi }_{\nu }^{2}=1.22$) is shown in Figure 6.

Zoom In Zoom Out Reset image size

To summarize, the X-ray emission from Arp 220 appears to be consistent with only a starburst. However, there is the possibility that a very deeply buried AGN is present in this source.

Braito et al. (2009) reported that the Superantennae was detected by Suzaku PIN above 10 keV. Although Braito et al. (2009) claim a signal-to-noise ratio of ∼10 in their 15–30 keV detection, the source spectrum is only ∼5.5% above the PIN background, which has a systematic uncertainty of ∼1.5%. They describe the source as a Compton-thick AGN (${N}_{{\rm{H}}}$ ∼ 3–4 × 10²⁴ cm⁻²) shining with an intrinsic 2–10 keV luminosity of a few times 10⁴⁴ erg s⁻¹, at the level of a luminous quasar. This galaxy was observed twice by NuSTAR, with a temporal separation of about four months. In both epochs, the Superantennae is weakly detected, with a 15–30 keV flux that is 30 times lower than that measured by Braito et al. (2009) with Suzaku. There appears to be no significant variability between the two sets of observations in terms of the spectral shape or the strengths of the Fe lines. Due to the lack of discernable variability, we have co-added the two sets of NuSTAR spectra in order to improve the signal-to-noise ratio of the overall spectrum. Emission lines at 6.5 and 6.9 keV were detected with Δχ² = 30 for four degrees of freedom.

5.5.1. Broadband Fitting

We have applied three different models to the broadband spectrum of the Superantennae. All three include thermal and power-law components for the starburst in addition to the typical AGN component. We first tested whether a simple absorbed power law can explain the spectral shape. A point in favor of the Compton-thick AGN scenario is that the power-law spectrum inferred from only <10 keV data is relatively flat (Γ ∼ 1.3; Braito et al. 2003). With the broader energy coverage of the XMM-Newton plus NuSTAR data, we find that the spectrum, after accounting for the starburst contribution, can be well fit with a standard power-law model for the AGN component (${\chi }_{\nu }^{2}=1.30$). The best-fit result requires only Galactic absorption and Γ of 1.54 ± 0.13, consistent with that measured from typical AGNs. We do not detect the presence of a strong neutral Fe line. The apparent Fe emission can be described by two narrow Gaussians with central energies at 6.54${}_{-0.07}^{+0.16}$ and ${6.87}_{-0.10}^{+0.37}$ keV. These lines have EWs of ${288}_{-94}^{+370}$ and ${296}_{-163}^{+521}$ eV, respectively. This fit is shown in Figure 7.

Zoom In Zoom Out Reset image size

For completeness, we also applied the torus models to the spectrum to test whether the AGN can also be Compton-thick as suggested by Braito et al. (2009). Both the MYTorus and BNTorus models give similar parameter values; however, the BNTorus models cannot constrain the opening angle of the torus. The MYTorus model fits the data very well (${\chi }_{\nu }^{2}=1.27$). The best-fit model implies that the underlying nuclear spectrum, with ${\rm{\Gamma }}={1.54}_{-0.14}^{+0.17},$ is obscured by a column of ${N}_{{\rm{H}}}={4.2}_{-3.1}^{+5.8}\times {10}^{24}$ cm⁻². A nominally small fraction, 13${}_{-9}^{+87}$%, of the direct emission is leaked through the absorber, but the large error bars clearly indicate a poorly constrained parameter.

The Compton-thin and Compton-thick models are statistically equivalent. Although technically the Compton-thick model has a smaller reduced χ², it is a more complicated model that does a poor job of constraining the leaked emission component. If the leaked fraction parameter is left completely free, the parameter errors reach an unphysical value. If the leaked component is removed, MYTorus cannot account for most of the 3–9 keV emission in the spectrum. Therefore, we favor the Compton-thin interpretation of the Superantennae spectrum. With this model, the intrinsic 2–10 keV luminosity of the AGN is 1.7 × 10⁴² erg s⁻¹, several hundred times lower than that reported by Braito et al. (2009). The 2–10 keV to bolometric luminosity ratio for this AGN is 0.08%.

In these models, the 15–30 keV flux for the Superantennae is ∼1.7 × 10⁻¹³ erg s⁻¹ cm⁻², a factor of 30 lower than the flux measured in this band by Braito et al. (2009) with Suzaku. There does not appear to have been any notable variability below 10 keV for the Superantennae. Braito et al. (2009) noted two other sources within the Suzaku PIN non-imaging field of view that are also AGNs and have similar fluxes to the Superantennae. While they used XMM-Newton data to constrain the spectral properties of these background sources and concluded that they can be described by unabsorbed power laws, it is possible that their contributions above 10 keV were not fully accounted for. Other field sources could also have contaminated the Suzaku >10 keV measurement. These include a field source within the NuSTAR field of view in the first observation of the Superantennae. The field source is a point source ∼85 from the Superantennae (R.A.: 19:32:48.3, decl.: –72:33:52.0). Although fainter than the Superantennae, the count rate of this source is 36% and 45% of those of the Superantennae in the 3–10 keV and 10–20 keV bands, respectively. These numbers suggest that the field source can harden the apparent 15–30 keV Suzaku PIN spectrum of the Superantennae at >10 keV, leading to the previous conclusion that the source is Compton-thick. There may be other field sources within the Suzaku field of view that are outside the NuSTAR field of view with similar properties. Therefore, we conclude that the Suzaku data were contaminated and that the Superantennae most likely hosts a Compton-thin AGN.

5.5.2. Iron Line Variability

IRAS 08572+3915, IRAS 10565+2448, and IRAS 14378−3651 were undetected by NuSTAR in 25 ks. Given the NuSTAR sensitivity limits, this implies that the intrinsic 2–10 keV luminosities of these sources are below ∼5 × 10⁴² erg s⁻¹ for the typical redshifts of our sources of ∼0.05, assuming the standard canonical AGN power-law model. Otherwise, the strong X-ray continua of these sources would have been detected above 10 keV.

Using the observed count rates from our NuSTAR observations at the locations of our targets, we determined upper limits to the observed 2–10 keV luminosities of these sources, under the assumption that the obscuring column density is not high. In the derivation, we assumed only a power-law component with Γ = 1.8 and Galactic absorption. No additional column density was assumed. Using the 3–10 keV count rates for each source, we used WebPIMMS²⁵ to estimate the unabsorbed 2–10 keV luminosity. The total 3–10 keV count rates extracted from circular regions with 1' radii for IRAS 08572+3915, IRAS 10565+2448, and IRAS 14378−3651 are 3 × 10⁻³, 2 × 10⁻³, and 3 × 10⁻³ counts per second, respectively. Assuming no intrinsic obscuration, these rates correspond to upper limits to the 2–10 keV luminosity of ∼6 × 10⁴¹ and 7 × 10⁴¹ erg s⁻¹ for IRAS 08572+3915 and IRAS 10565+2448, respectively, and an intrinsic 2–10 keV luminosity of 1 × 10⁴² erg s⁻¹ for IRAS 14378–3651. These limits and measurement are consistent with those previously measured by Chandra (e.g., Teng et al. 2009; Teng & Veilleux 2010; Iwasawa et al. 2011), and suggest that the intrinsic 2–10 keV to bolometric luminosity ratios in these sources are ≲0.01, ≲0.09, and 0.08%, respectively, assuming no source obscuration. However, as we saw in the case of Arp 220 (Section 5.4), it is not possible to rule out very high obscuring column densities.

For several decades, what powers the enormous infrared luminosities of ULIRGs has remained an unanswered question. Many studies concluded that the lack of strong X-ray detections in ULIRGs implied that these sources are highly obscured (e.g., Franceschini et al. 2003; Ptak et al. 2003; Teng et al. 2005; Iwasawa et al. 2011). However, these studies lacked sensitive detections at energies above 10 keV, which are necessary to disentangle the effects of obscuration in order to robustly measure the intrinsic X-ray luminosities of the AGNs in ULIRGs. With the launch of NuSTAR, which is ∼100 times more sensitive than Suzaku PIN in the 10–40 keV energy band, broadband (0.5–30 keV) X-ray spectroscopy has allowed us to estimate the intrinsic X-ray luminosity of five of the nine ULIRGs in our sample and place constraints on the remaining four.

Our observations reveal that the ULIRGs in our sample have surprisingly low observed fluxes in high-energy (>10 keV) X-rays. Of the nine ULIRGs in our NuSTAR sample, six were detected well enough to enable detailed spectral modeling of their broadband X-ray spectra. Of these six, only one, IRAS 13120–5453, has a spectrum consistent with a Compton-thick AGN. We cannot rule out the possibility that a second ULIRG in the sample, Arp 220, is highly Compton-thick (N_H > 10²⁵ cm⁻²). Thus, by the strictest definition (N_H > 1.5 × 10²⁴ cm⁻²), these NuSTAR data show that most of the ULIRGs in our sample are not Compton-thick.

However, detailed analysis of the NuSTAR data on IRAS 05189–2524 and Mrk 273 shows that the hard X-ray fluxes of these sources have varied compared to similar Suzaku observations in 2006 (Teng et al. 2009). The observed variability in both sources can be explained by a change in the absorbing column. The column densities for both IRAS 05189–2524 and Mrk 273 have reduced by a factor of a few since the 2006 observations. These changes are sufficient to alter the classification of these AGNs from being Compton-thick to Compton-thin.

In our sample of nine ULIRGs, three (IRAS 05189–2524, Mrk 231, and Mrk 273) have strong hard X-ray continua above 10 keV after correcting for obscuration. The remaining targets have low count rates. Only two ULIRGs (IRAS 05189–2524 and IRAS 13120–5453) have intrinsic 2–10 keV AGN luminosities above 10⁴³ erg s⁻¹. Although the AGNs in many of these ULIRGs dominate the spectral energy distributions at other wavelengths (e.g., Veilleux et al. 2009a), these active black holes do not appear to produce as much X-ray emission as would be expected for typical AGNs. When compared to their bolometric luminosities, the AGNs in our sample of ULIRGs are not emitting as much X-ray power as Seyfert galaxies. The intrinsic 2–10 keV to bolometric luminosity ratio is in the range 0.03%–0.81%. For comparison, this value ranges from 2% to 15% for radio-quiet quasars and Seyfert galaxies (Elvis et al. 1994). However, there is evidence that objects with Eddington ratios near or above unity have smaller 2–10 keV to bolometric luminosity ratios (0.3%–0.7%; Vasudevan & Fabian 2009; Lusso et al. 2010, 2012). The low 2–10 keV to bolometric luminosity ratios for our sample could imply that these sources have high accretion rates. This is unsurprising, as ULIRGs are mergers that may be rapidly growing their central black holes. The high accretion rates of ULIRGs are consistent with the detection of ionized Fe K lines in five of the six ULIRGs with spectra in our sample (IRAS 05189–2524, IRAS 13120–5453, Mrk 231, Arp 220, and the Superantennae). Iwasawa et al. (2012) suggested that there is a link between the presence of ionized lines in the COSMOS sample and high accretion rate. The exception is Mrk 273, which has a low accretion rate (see below).

For the three sources in our sample that have both dynamically measured black hole masses (e.g., Veilleux et al. 2009b) and spectrally derived AGN fluxes, we calculated Eddington luminosities, Eddington ratios, and 2–10 keV to Eddington luminosity ratios, tabulated in Table 4. These parameters distinguish the brightest objects in our sample each into a unique category. IRAS 05189–2524 is the X-ray brightest and accreting at a super-Eddington rate. Mrk 231 is the X-ray faintest and also accreting at a super-Eddington rate. Mrk 273 has the most massive black hole and the lowest Eddington ratio. Many studies have found that Γ becomes softer with increasing λ_Edd (e.g., Shemmer et al. 2005, 2008; Brightman et al. 2013). Comparing the model-derived Γ of our three brightest sources with their λ_Edd values, this Γ–λ_Edd correlation appears to hold with the exception of Mrk 231. More recently, Yang et al. (2014) found a correlation between Γ and the dimensionless ratio of 2–10 keV and Eddington luminosities, which holds for black hole accretion systems including both black hole binaries and AGNs. The authors found that Γ decreases with increasing l_X up to ∼10⁻³ but steepens again above that l_X value. (See also Constantin et al. 2009.) Our limited data points from three sources are consistent with this phenomenological model. If both the Γ–λ_Edd and Γ–l_X relations are true and indicative of the accretion processes in most AGNs, then Mrk 231 is an outlier. This implies that the 2–10 keV to bolometric luminosity correction is different for Mrk 231 than for most AGNs, as suggested by Teng et al. (2014).

Table 3.  Best-fit Parameters for the ULIRGs in the Deep Survey

Model Parameter	IRAS 05189–2524	IRAS 13120–5453	Mrk 273	Arp 220	Superantennae	Comment on Parameter
(1)	(2)	(3)	(4)	(5)	(6)	(7)
B/A	1.02${}_{-0.03}^{+0.03}$	1.05 (f)	0.99${}_{-0.06}^{+0.06}$	1.05 (f)	1.05 (f)	FPMB-A cross-normalization
XMM/A or CXO/A	0.86${}_{-0.02}^{+0.02}$	1.20* (f)	0.88${}_{-0.14}^{+0.15}$	0.90 (f)	0.85 (f)	XMM-Newton-FPMA or Chandra (*)-FPMA cross-normalization
kT (keV)	0.16${}_{-0.01}^{+0.01}$	0.56${}_{-0.06}^{+0.05}$	0.64${}_{-0.15}^{+0.15}$	0.16${}_{-0.02}^{+0.02},$ 0.69${}_{-0.10}^{+0.10}$	0.51${}_{-0.17}^{+0.20}$	MEKAL gas temperature from the starburst
Brems. kT (keV)	...	...	...	9.13${}_{-2.04}^{+3.12}$	...	bremsstrahlung temperature from the hot gas
Abs. 1 (1022 cm−2)	5.19${}_{-0.18}^{+0.20}$	315.7${}_{-129.4}^{+232.5}$	43.8${}_{-5.7}^{+9.5}$	0.38${}_{-0.05}^{+0.04}$	...	neutral absorber 1
cf 1	0.98${}_{-0.01}^{+0.01}$	1 (f)	1 (f)	1(f)	...	covering factor 1
Abs. 2 (1022 cm−2)	9.32${}_{-0.68}^{+0.95}$	...	...	...	...	neutral absorber 2
cf 2	0.74${}_{-0.02}^{+0.01}$	...	...	...	...	covering factor 2
Abs. HMXB (1022 cm−2)	$\gt 2.66$	0.58${}_{-0.58}^{+1.32}$	...	...	...	neutral absorber applied to the HMXB power-law component
ΓHMXB	1.1 (f)	1.1 (f)	1.1 (f)	...	1.1 (f)	HMXB cutoff power-law index with cutoff energy at 10 keV
ΓAGN	2.51${}_{-0.02}^{+0.02}$	1.8 (f)	1.43${}_{...}^{+0.17}$	...	1.54${}_{-0.07}^{+0.07}$	AGN power-law index (MYTorus lower limit fixed at 1.4)
Inc (deg)	...	85 (f)	85 (f)	...	...	inclination angle
Eline(keV)	6.43${}_{-0.05}^{+0.05}$	1.86${}_{-0.05}^{+0.11}$	...	6.78${}_{-0.07}^{+0.07}$	6.53${}_{-0.11}^{+0.16}$	line 1
EWline(keV)	0.074${}_{-0.035}^{+0.033}$	0.104 (uc)	...	0.899${}_{-0.398}^{+0.469}$	0.296${}_{-0.107}^{+0.447}$	EW of line 1
E${}_{{\rm{line}}}$ 2 (keV)	6.80${}_{-0.04}^{+0.02}$	3.40${}_{-0.06}^{+0.11}$	...	...	6.88${}_{-0.08}^{+0.37}$	line 2
EWline 2 (keV)	0.117${}_{-0.040}^{+0.028}$	0.460${}_{-0.279}^{+0.348}$	...	...	${0.330}_{-0.148}^{+0.538}$	EW of line 2
Eline 3 (keV)	...	6.86${}_{-0.12}^{+0.28}$	...	...	...	line 3
EWline 3 (keV)	...	0.848${}_{-0.562}^{+0.674}$	...	...	...	EW of line 3
Const. (C-thin)	...	...	0.03${}_{-0.02}^{+0.02}$	...	...	Compton-thin fraction
f0.5–2 (10−13 erg s−1 cm−2)	${1.50}_{-0.50}^{+0.48}$	0.50${}_{-0.24}^{+0.08}$	${0.98}_{-0.14}^{+0.11}$	0.98${}_{-0.04}^{+0.04}$	0.75${}_{-0.05}^{+0.04}$	observed 0.5–2 keV flux
f2–10 (10−12 erg s−1 cm−2)	${3.87}_{-0.13}^{+0.11}$	0.26${}_{-0.33}^{+0.30}$	${0.76}_{-0.25}^{+0.03}$	0.12${}_{-0.01}^{+0.01}$	0.23${}_{-0.01}^{+0.01}$	observed 2–10 keV flux
f10–30 (10−12 erg s−1 cm−2)	${2.74}_{-0.16}^{+0.11}$	0.89${}_{-0.78}^{+0.29}$	${2.97}_{-1.40}^{+0.09}$	0.05${}_{-0.02}^{+0.02}$	0.25${}_{-0.03}^{+0.03}$	observed 10–30 keV flux
LMEKAL (erg s−1)	1.40 × 1041	2.03 × 1041	1.47 × 1041	2.80 × 1041	1.77 × 1041	intrinsic MEKAL 0.5–30 keV luminosity
LBrems (erg s−1)	...	...	...	1.66 × 1040	...	intrinsic bremsstrahlung 0.5–30 keV luminosity
LHMXB (erg s−1)	4.70 × 1041	7.62 × 1041	7.09 × 1041	...	2.79 × 1041	intrinsic HMXB 0.5–30 keV luminosity
L0 (erg s−1)	5.58 × 1043	3.15 × 1043	2.40 × 1043	...	4.57 × 1042	intrinsic AGN 0.5–30 keV luminosity
L0 (2–10 keV) (erg s−1)	3.69 × 1043	1.25 × 1043	8.55 × 1042	...	1.70 × 1042	intrinsic AGN 2–10 keV luminosity
LAGN/Lbol,AGN (%)	0.81	0.67	0.28	...	0.08	2–10 keV X-ray-to-bolometric luminosity ratio for the AGN
χ2/dof	714.2/666	36.3/48	189.5/209	122.9/101	146.22/112	goodness of fit

Note. (f) denotes a fixed parameter and (uc) detnotes an unconstrained parameter. Col. (1): model parameter for each fit. Col. (2): best-fit model: const. × N${}_{{\rm{H,Galactic}}}$ (MEKAL + Abs${}_{\text{nuclear HMXB}}\;\times$ cutoffPL${}_{\text{nuclear HMXB}}$ + Abs${}_{{\rm{1}}}\;\times$ Abs${}_{{\rm{2}}}\;\times$ PL_AGN + Line(6.4 keV) + Line(6.7 keV)). Col. (3): best-fit model: const. × N${}_{{\rm{H,Galactic}}}$ (MEKAL + Abs${}_{\text{nuclear HMXB}}\;\times$ cutoffPL${}_{\text{nuclear HMXB}}$ + Line(6.7 keV) + MYTorus × PL_AGN). Col. (4): best-fit model: const. × N${}_{{\rm{H,Galactic}}}$ (MEKAL + cutoffPL${}_{\text{nuclear HMXB}}$ + MYTorus × PL_AGN + Const.${}_{{\rm{C}}-\mathrm{thin}}\;\times$ PL${}_{{\rm{AGN}}}$). Col. (5): best-fit model: const. × N${}_{{\rm{H,Galactic}}}\times$ Abs₁ ×(MEKAL₁ + MEKAL₂ + zbremss + Line(6.7 keV)). Col. (6): best-fit model: const. × N${}_{{\rm{H,Galactic}}}$ (MEKAL + cutoffPL${}_{\text{nuclear HMXB}}$ + Line(6.7 keV) + Line(6.9 keV) + PL_AGN). Col. (7): comments on model parameter. In the second row, the Chandra-FPMA cross-normalization is marked with an asterisk; the others are for XMM-Newton-FPMA.

Download table as:  ASCIITypeset image

In Figure 8, we compare the intrinsic 2–10 keV to bolometric luminosity ratio, also known as the X-ray bolometric correction, to the fraction of the bolometric luminosity attributed to the AGN from infrared measurements. This latter ratio, L_bol,AGN, is the average AGN fraction calculated from six independent methods that include fine-structure line ratios, mid-infrared continuum ratios, and the EWs of the aromatic features (Veilleux et al. 2009a). Included in the figure are also two lines that help guide the eye for where a pure AGN contributing 100% of the luminosity in the infrared with 1% and 10% X-ray bolometric correction would lie for a given AGN fraction.

Zoom In Zoom Out Reset image size

X-ray bolometric corrections of 2%–15% are typical for Seyferts and radio-quiet quasars (Elvis et al. 1994), with much lower bolometric corrections of ∼0.3%–0.7% seen for objects accreting at close to the Eddington rate (e.g., Vasudevan & Fabian 2009; Lusso et al. 2010, 2012). While two of the eight ULIRGs lie within the range of bolometric corrections that is typical for typical Seyferts, the majority appear surprisingly faint in the X-rays. Their X-ray bolometric corrections are more in line with those of objects accreting at close to the Eddington rate.

This disagreement between infrared and X-ray diagnostics is particularly large for Mrk 231 and IRAS 08572+3915: while Spitzer diagnostics find that these two sources are heavily dominated by AGNs, both are very underluminous in the X-rays. The ratio of intrinsic 2–10 keV luminosity to bolometric luminosity for Mrk 231 is only 0.03%, and the ratio for IRAS 08572+3915 is even lower (assuming no intrinsic absorption). Teng et al. (2014) established that Mrk 231 is intrinsically X-ray weak, likely related to the powerful wind detected in this broad absorption line (BAL) quasar. There is growing evidence emerging that some AGNs with strong outflows, such as some BAL quasars, are intrinsically X-ray weak (e.g., Luo et al. 2013, 2014; Teng et al. 2014), suggesting that intrinsic X-ray weakness and strong winds may be linked. IRAS 08572+3915, known to have a strong outflow on kiloparsec scales (Rupke & Veilleux 2013), may be another example of an intrinsically X-ray weak AGN with powerful winds. In fact, IRAS 08572+3915 is a more extreme case of intrinsic X-ray weakness than Mrk 231. Efstathiou et al. (2014) identified IRAS 08572+3915 as the most infrared-luminous galaxy, with its AGN contributing ∼90% of its total power output.

Ultrafast outflows have recently been discovered in two nearby ULIRGs, IRAS F11119+3257 (Tombesi et al. 2015) and Mrk 231 (Feruglio et al. 2015), that also host large-scale molecular outflows. Although large-scale outflows are quite different from the disk winds directly observed in BAL quasars, the presence of both kinds of outflows in Mrk 231 suggests that the two phenomena may be related.

In Figure 8 we considered the fraction of bolometric luminosities attributable to the AGN. Now we consider diagnostics of the intrinsic AGN luminosities. One such method is the [O iv] 26 μm luminosity to 2–10 keV luminosity relation (e.g., Melendez et al. 2008; Diamond-Stanic et al. 2009; Rigby et al. 2009; Weaver et al. 2010). [O iv] 26 μm is far more robust to extinction than is the optical emission line [O iii], and thus a more suitable diagnostic for highly obscured objects like ULIRGs, for which star formation also contributes significantly to their overall power.

We first investigate the possibility that the [O iv] line is contaminated by star formation. ULIRGs have significant circumnuclear star formation (see Teng et al. 2014 for a detailed analysis on Mrk 231), and the hottest stars are capable of creating some of the 55 eV photons required to triply ionize oxygen. We investigate this possibility by comparing the [O iv] fluxes to [Ne v]. Even the hottest stars should not generate significant amounts of the 97 eV photons required to quadruply ionize neon within a galaxy; accordingly, [Ne v] is not detected in the mid-infrared spectra of star-forming galaxies. Goulding & Alexander (2009) find a tight relation between the [O iv] 26 μm and [Ne v] 14.3 μm emission in AGNs, shown in their Figure 10, and argue that this is a method of determining whether an AGN dominates the [O iv] luminosity. Four ULIRGs in our sample have both [O iv] and [Ne v] detected; in the other objects neither line is detected. These four objects lie close to the relation of Goulding & Alexander (2009), at the upper right corner with 41 < log (L[O iv]/erg s⁻¹) < 42.2. If anything, these objects have somewhat higher [Ne v]/[O iv] ratios than the relation. This argues that the [O iv] emission, like the [Ne v] emission, is powered by the AGN, and that there is no significant [O iv] contamination from extremely hot stars.

We therefore proceed to compare the measured intrinsic 2–10 keV luminosities of our ULIRGs with the expected intrinsic 2–10 keV luminosities predicted by the measured [O iv] luminosity and the relation of Rigby et al. (2009). Figure 9 shows that the [O iv] relation, which was calibrated from empirical measurements of Seyfert 1 galaxies, overpredicts the intrinsic 2–10 keV luminosity for all sources in our sample. We now discuss possible reasons for this discrepancy.

Zoom In Zoom Out Reset image size

First, it is worth noting that [O iv] was only detected in half the sample. We attribute this to the difficulty of detecting a weak line over a very strong continuum in the Spitzer IRS spectra that have moderate signal-to-noise ratio. It is possible, though implausible, that the actual [O iv] fluxes for the non-detections are very low, enough to pull down these four sources to match the dashed line in Figure 9.

The second possibility is that the AGNs in some ULIRGs have different UV-to-X-ray ratios than Seyfert galaxies from which the relation is derived. The [O iv] and [Ne v] emission arises in the narrow line region, where 55 and 97 eV (respectively) photons that originate in the accretion disk are able to triply ionize oxygen and quadruply ionize neon. Thus, these high ionization fine-structure lines are measures of the intrinsic UV luminosity. The X-ray luminosity is also tied to the accretion disk, as it is produced by Compton upscattering of UV photons produced in the accretion disk by hot electrons in the corona. Thus, Figure 9 suggests that the AGNs in ULIRGs may have different X-ray-producing efficiencies than do typical AGNs in Seyfert galaxies, which we speculate might be due to differences in the structure or geometry of the corona.

To contextualize, the long effort to determine the AGN content of ULIRGs has involved two parallel paths: extinction-robust hard X-ray photons and extinction-robust mid-infrared signatures. Both approaches assume that the AGN within a ULIRG has the properties of a typical Seyfert AGN, and is merely highly obscured. Our work reveals that these assumptions may break down in certain regimes. Objects like Mrk 231 and IRAS 08572+3915 have infrared signatures of AGN dominance, but are severely underluminous in the X-rays, and the cause is not high extinction. The rest of the sample appears X-ray underluminous as well (see Figures 8 and 9). Imanish & Terashima (2004) reached similar conclusions for IRAS 08572+3915 and three other ULIRGs by comparing Hα to Chandra and XMM-Newton spectra. The ULIRGs in our sample have systematically low intrinsic 2–10 keV luminosities for their measured [O iv] 26 μm luminosities, compared to a relation calibrated with Seyfert 1 AGNs. Both lines of evidence argue that the AGNs within some, perhaps many, ULIRGs may not be typical Seyfert 1 nuclei. We suggest that the violent conditions within some ULIRGs—high accretion rates, strong winds—may be feeding back and affecting the X-ray properties of their AGN.