Probing the spectrum of the magnetar 4U 0142+61 with JWST

The Near Infrared Camera (NIRCam; Rieke et al. 2023) was used to image the magnetar and measure its flux in two spectral bands. The NIRCam observations were carried out on 2022 September 21 from 00:34:19 to 00:43:48 UTC. The F250M and F140M filters were employed in the long and short-wavelength NIRCam channels, using the NRCBLONG and NRCB1 detectors, respectively, in the SUB400P subarray configuration. The pivot wavelengths and bandwidths for these filters are listed in Appendix A and Figure 1 shows the wavelength dependencies of the filter throughputs.

We used the RAPID readout pattern, having seven total dithered integrations and 10 groups per integration. The effective scientific exposure time in each of the F250M and F140M filters is 115.9 s. We use the pipeline-processed data (calibration software version 1.11.4; Bushouse et al. 2023). In the drizzle process, a pixel shrinking of 1.0, pixel scale ratio of 1.0, and inverse variance map (IVM) weighting scheme were used for the final resampling of the data. The nominal pixel sizes are 62.7 mas and 30.7 mas for the F250M and F140M images, respectively.

Both the 250M and F140M images show a strongly nonuniform background, dominated by the so-called 1/f noise patterns (Schlawin et al., 2020). The details of the source photometry for these noisy images are described in Appendix A. The net source flux densities in the F250M and F140M filters are $f_{\nu}=14.5\pm 0.4$ $\mu$Jy and $5.1\pm 0.3$ $\mu$Jy, respectively. The statistical uncertainties here and throughout the paper are reported at the $1\sigma$ level unless noted otherwise.

For the F250M filter, which has a larger field of view, we compared the radial profile of 4U 0142+61 with those of other (likely stellar) sources in the field and found them to be consistent with each other. Hence, we conclude that there is no indication of any extended NIR emission around the magnetar. The field of view for the F140M filter was smaller and no other sources were detected, leaving nothing to directly compare the radial profile to. The empirical PSF FWHM for the F250M (F140M) filter is 85 (48) mas, which constrains the angular size of the source to be $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}0.1^{\prime\prime}$. At a distance of 3.6 kpc, the F250M limit corresponds to a physical limit on the size of a disk or IR wind-nebula of $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}10^{16}$ cm or $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}2\times 10^{5}$ $R_{\odot}$.

To examine a contemporaneous spectral energy distribution (SED) from IR to hard X-ray, we carried out nearly simultaneous observations of the magnetar with the NuSTAR and Swift X-ray observatories.

The broadband X-ray spectrum of 4U 0142+61 has been studied in detail in many previous works (see e.g., Rea et al. 2007; Gonzalez et al. 2010; Enoto et al. 2011; Tendulkar et al. 2015). It has been shown in these works that the soft part of the X-ray spectrum, $E\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}$ 10 keV, requires two spectral components to get an acceptable fit, while a third component (namely a hard PL) is required if energies $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}$ 10 keV are included in the fit. We jointly fit the Swift and NuSTAR spectrum of 4U 0142+61 to determine whether or not the hard PL tail connects smoothly to the IR spectrum, thus may be produced by the same particle population. To fit these spectra, we use an empirical absorbed BB+PL+PL model (model 1 hereafter), multiplied by a constant (frozen to 1 for Swift-XRT), to account for any cross-calibration differences between Swift and NuSTAR FPMA/B (i.e., const$\times$tbabs$\times$(bbodyrad+pow+pow) in XSPEC). We also fit an empirical BB+BB+PL (model 2 hereafter) also multiplied by a constant (i.e., const$\times$tbabs$\times$(bbodyrad+bbodyrad+pow)), as the previous work by Tendulkar et al. (2015) showed that it can fit the spectra equally well, but provides a significantly different photon index for the hard power-law tail.

We report the best-fit values for both models in Table 1. Statistically acceptable fits were achieved for both models, although the BB+PL+PL model fits the data better than the BB+BB+PL model. Overall, the best-fit parameters for model 1 from our dataset agree with the values found by Tendulkar et al. (2015). However, we find higher temperatures for both thermal components in model 2 and a harder photon index (i.e., $\Gamma_{h}$=0.83 for our fits versus 1.03 from Tendulkar et al. 2015). These differences may be due to changes in the NuSTAR calibration between the NuSTAR observations and analyses. In any case, we are primarily focused on the hard PL component of the spectrum, as it is the most likely to connect to the IR. We find that the photon index of the hard power-law tail is $\Gamma_{h}=0.35\pm 0.07$ ($\alpha_{h}=-0.65\pm 0.07$) or $\Gamma_{h}=0.83\pm 0.06$ ($\alpha_{h}=-0.17\pm 0.06$) for model 1 and model 2, respectively. Figure 3 shows the best-fit for model 1, the fit residuals, and the unfolded spectrum and individual model components.

Wang et al. (2006) reported on a Spitzer Infrared Array Camera (IRAC) observation of 4U 0142+61, taken on 2005 January 17, however, Wang & Kaspi (2008) re-analyzed these data and reported more accurate spectral fluxes (based on PSF fitting photometry instead of aperture photometry) $f_{\nu}=32.1\pm 1.2$ and $48.8\pm 7.6$ $\mu$Jy at 4.5 and 8 $\mu$m, respectively. Spitzer also observed the source four times after it underwent an X-ray outburst on 2007 February 7. The IRAC observations spanned 2007 February 14-21 and measured average spectral fluxes, $f_{\nu}=32.1\pm 2.0$ and $59.8\pm 8.5$ $\mu$Jy at 4.5 and 8 $\mu$m, respectively (Wang & Kaspi, 2008), which are consistent with the observation from 2005. These points agree well with the JWST NIRCam photometry and MIRI spectrum, obtained almost 18 years later, with the largest offset between all of the Spitzer IRAC photometry and our PL fit being only about 1.5$\sigma$ for the 4.5$\mu$m flux measured from the 2005 observation (see Figures 4, 5, and 6).

Spitzer also observed the source with the Infrared Spectrograph (IRS) on 2006 January 22. The analysis of these data was reported by Wang et al. (2008), who discuss the possible detection of a silicate emission feature at 9.7 $\mu$m. Unfortunately, the Spitzer IRS spectrum is too noisy for a meaningful comparison with the MIRI LRS spectrum. However, we found no evidence of the 9.7 $\mu$m emission feature in the LRS spectrum.

From the Spitzer observation with the Multiband Imaging Photometer (MIPS) on 2006 February 19, about a year after the Spitzer IRAC observations, Wang et al. (2008) reported a $3\sigma$ upper-limit of 38 $\mu$Jy at 24 $\mu$m. Posselt et al. (2011) re-analyzed this MIPS observation and found a possible source at the position of 4U 0142+61 in two independent, re-processed AORs (Astronomical Observation Requests), with a flux of $41$ $\mu$Jy, but it was only a $2\sigma$ significance detection. We use the latter measurement for this paper, since the analysis was performed using an updated Spitzer data reduction pipeline, and because it is close to the 3$\sigma$ upper-limit derived by Wang et al. (2008)⁴⁴4Note that the $3\sigma$ upper-limit from Posselt et al. (2011) would be about 60$\mu$Jy. This upper limit (or possible detection) from the MIPS observation strongly disagrees with the PL fit from the JWST data (see Figure 4). We find that a spectral break of about $\Delta\alpha\simeq 1.5$ (or $\Delta\alpha\simeq 1.0$, assuming a 60$\mu$Jy upper-limit) would be necessary to be consistent with the MIPS upper-limit assuming the break occurs at the long-wavelength end of the MIRI LRS spectrum (i.e., at 11$\mu$m).

There have been numerous optical and NIR observations of 4U 0142+61 (e.g., Hulleman et al. 2000, 2004; Dhillon et al. 2005; Morii et al. 2005; Durant & van Kerkwijk 2006; Muñoz-Darias et al. 2016; Chrimes et al. 2022). Wang et al. (2006) used published results from Hulleman et al. (2004) and Israel et al. (2004), which included four optical ($B,V,R,I$) and three NIR ($J,H,K_{s}$) points (see Figure 4). We also have gathered the results using different observations from those used in Wang et al. (2006) and plotted them in Figure 5 (to avoid overcrowding the plots due to many data points measured in similar bands), together with the best PL fit of the JWST data. In Figure 5, the 2007 Spitzer IRAC fluxes from Wang & Kaspi (2008) are plotted. The $J,H$, and $K_{s}$ points, shown by filled triangles, were taken from the Gemini observations of 2004 November 2, reported by Durant & van Kerkwijk (2006). We added two square points that span the minimum and maximum $K$-band values reported by Durant & van Kerkwijk (2006). The F125W and F160W fluxes, in bands centered at 1.25 and 1.6 $\mu$m, were obtained by Chrimes et al. (2022) from HST WFC3/IR observations of 2018 January 1, and the optical GTC/OSIRIS fluxes in the SDSS $z,r,i,g$ bands (which were all observed on the same night) were adopted from the 2013 observations reported by Muñoz-Darias et al. (2016). We chose these datasets as they consist of observations across several filters which were obtained in the same observing run (i.e., nearly simultaneous for each individual dataset).

As shown in Figures 4 and 5, the results of these different observing campaigns exhibit offsets of different sizes from the PL fit to the MIRI LRS data. Notably, the points in Figure 4 show a strong discrepancy with the high-frequency extension of the LRS PL fit, as well as with the results of more recent NIR-optical observations. However, the strong discrepancies between the NIR data and LRS PL model still remain even in the more recent observations. The discrepancies, defined as $\Delta\chi=$(data-model)/$\sigma$ between the broadband data points from Figure 5 and the PL fit to the LRS spectrum, are shown in Figure 6.

The largest discrepancy between the model and photometry ($\approx 15\sigma-25\sigma$) occurs in the HST bands, primarily due to the small uncertainties on the HST points. The second largest discrepancy occurs in the $K_{s}$ band, where the model overpredicts the $K_{s}$ flux by a factor of about 1.7 (or $\sim$11$\sigma$). Additionally, the $H$ and $J$ band points lie below the model but are less discrepant than the $K_{s}$ band point (at about 10$\sigma$ and 5$\sigma$, respectively). The two HST points, which lie above the nearby $J$ and $H$ band points, are still about a factor of 1.2 smaller than flux predicted by the best-fit model to the LRS data. Surprisingly, the NIRCam points agree remarkably well with the best-fit PL model, even though these points were not included when fitting the model. One of these points lies a bit redward of the $K_{s}$ band point, while the other is between the $H$ and $J$ bands. JWST is supposed to have an absolute flux calibration that is much better than $10\%$ (typically better than $2\%$), making calibration an unlikely cause for the discrepancies (Gordon et al., 2022). The optical points from Muñoz-Darias et al. (2016) also agree well with the extrapolated PL spectrum (as does their best-fit spectral index, $\alpha=0.8\pm 0.5$, albeit with large uncertainties), but these points are the most reliant on the best-fit $A_{V}$, which has changed throughout the various calibration updates (see Section 2.2.2). We discuss potential causes for these discrepancies, such as variability and/or different emission sites for the IR and optical emission, in the following subsections.

Since the IR-optical observations of the magnetar were taken at different epochs, an obvious explanation for the discrepancies between the MIRI LRS spectrum (and NIRCam photometry) and previous broadband photometry is source variability (see Figure 6). Originally, Hulleman et al. (2004) reported variability at NIR wavelengths, particularly in the $K$ band, but noted that the source did not appear to be variable at optical wavelengths. $K$-band variability of greater than one magnitude, over a timescale of a few days, was reported by Durant & van Kerkwijk (2006). These authors also found that the source was possibly variable at optical wavelengths, by about half a magnitude in the $I$ band. On the other hand, Muñoz-Darias et al. (2016) reported a lack of variability in the optical $z$, $i$, $r$ and $g$ bands, but they did not comment on the possibility of NIR variability. Chrimes et al. (2022) found that the NIR HST fluxes were consistent with the source being variable when compared to the NIR results of Hulleman et al. (2004) and Durant & van Kerkwijk (2006). Additionally, using AKARI observations, Kohmura et al. (2013) found a 64% lower flux at 4 $\mu$m compared to the Spitzer 4.5 $\mu$m flux and suggested that this discrepancy was due to variability.

There is also the MIPS 24 $\mu$m upper limit, which is about a factor of three lower than the expected flux extrapolated from the best PL fit to the MIRI data. The factor of three discrepancy in flux is consistent with the extent of the variability reported by Durant & van Kerkwijk (2006) in the NIR (demonstrated in Figure 5 for the $K_{s}$ band). Therefore, we cannot exclude the possibility that the mid-IR and NIR fluxes vary by about a factor of three, giving rise to the observed discrepancies. However, the good agreement of the Spitzer IRAC fluxes, from both Wang et al. 2006 and Wang & Kaspi 2008, with the JWST LRS + NIRCam spectrum, observed 17+ years later, as well as the lack of (or small scale) variability observed at optical wavelengths challenge this scenario. Additionally, Table 1 in Wang & Kaspi (2008) shows that the IRAC flux of 4U 0142+61 did not vary over the week long observing campaign.

It could also be the case that the IR and optical emission are produced by different mechanisms and/or at different sites near the magnetar, leading to variability in the IR that is not observed in the optical. However, there is a fair agreement between the optical observations and extrapolated best-fit PL model, suggesting that both emission components are produced by the same mechanism. Thus, the question of the origin and true extent of the possibly wavelength-dependent variability remains open.

Another explanation of the possibly different IR and optical spectra of 4U 0142+61 is the presence of a fallback disk whose contribution to the observed emission may be different in different wavelength ranges. The discovery of optical pulsations (Kern & Martin, 2002; Dhillon et al., 2005), with a pulsed fraction of about 25%–30%, significantly higher than $\sim 4\%$ in X-rays, has led to perception that the optical emission comes from the magnetar’s magnetosphere. If true, it means that the hypothetical fallback disk would be more easily detected in the IR. Wang et al. (2006) noticed that the previously reported $V,R,I$ fluxes are consistent with a PL model with $\alpha=-0.3$, at an assumed $A_{V}=3.5$. They used this PL in combination with a disk model to fit the IR + optical data (see Figure 4). That model assumes that the optically thick disk is heated by the X-ray emission from the magnetar, the disk’s local effective temperature decreases with radial distance from the magnetar as $T\propto r^{-3/7}$, and each point on the disk surface emits BB radiation with this temperature. Wang et al. (2006) found that the best-fit disk model had inner and outer disk temperatures $T_{\rm in}=1200$ K and $T_{\rm out}=715$ K, and inner and outer disk radii $r_{\rm in}=2.9R_{\odot}$ and $r_{\rm out}=9.7R_{\odot}$, respectively, at $d=3.6$ kpc and $\cos i=0.5$, where $i$ is the disk inclination. Because of the small difference between the inner and outer temperatures, the disk spectrum resembles a BB spectrum with a temperature of 920 K.

We reproduced the spectrum from Wang et al. (2006) using a simple disk spectrum model given by Equation (1). Figure 4 shows the model disk+PL spectrum together with the results of the optical/IR photometry (adopted from Wang et al. 2006 with updated IRAC fluxes from Wang et al. 2008) and the JWST data. It is immediately apparent that the JWST data are strongly inconsistent with this model. The MIRI spectrum does not have the anticipated curvature, and it overshoots the model at the long wavelength end. Additionally, the NIRCam F140M point greatly exceeds the anticipated model flux. Moreover, the most recent optical points from Muñoz-Darias et al. (2016) do not follow a PL with $\alpha=-0.3$, at least for the best-fit $A_{V}=3.9$ found from the PL fit of the MIRI LRS data, as we see from Figure 5.

Although the PL+disk model used by Wang et al. (2006) is inconsistent with the JWST data, it does not mean that another disk model, perhaps in combination with a PL component, cannot describe the IR or IR+optical data. For instance, in the approximation of a multi-temperature BB flat disk, the flux density spectrum is

If the radial dependence of the local effective temperature of a disk obeys a PL, $T(r)=T_{\rm in}(r/r_{\rm in})^{-\beta}=T_{\rm out}(r/r_{\rm out})^{-\beta}$, then the disk spectrum is close to a PL with the slope $\alpha=2/\beta-3$, i.e.,

Thus, the spectrum of a disk with a sufficiently large ratio of the inner and outer temperatures, $T_{\rm in}/T_{\rm out}\gg 2/\beta-1$ (i.e., the large ratio of the outer and inner radii, $r_{\rm out}/r_{\rm in}\gg(2/\beta-1)^{1/\beta}$) has a broad PL part. If we assume that the JWST PL spectrum is due to such a disk, then $\beta=2/(\alpha+3)=0.505$ for $\alpha=0.96$ (or $\beta=0.5$ for $\alpha=1$).

In a self-consistent disk model, the radial dependence of temperature should be derived from a balance between heating and emitted energies. In original models for a passive protostellar flat disk, illuminated by a central star, the dependence $T\propto r^{-3/4}$ was derived, which corresponds to $\alpha=-1/3$ (see Adams et al. 1987, and references therein). However, most protostellar disks exhibit flattish SEDs in the IR range, with $\alpha\approx 0.25$–1 (e.g., Beckwith et al. 1990), which corresponds to $\beta\approx 0.5$–0.6, similar to the slope we found for 4U 0142+61. The difference between the observed and model temperature dependencies is likely due to some unrealistic assumptions in the original models. In particular, realistic disks should be flared rather than flat, and the local emission spectra are not BBs because the temperature of an irradiated passive disk should decrease inward (toward the disk mid-plane), and the disks become optically thin at large wavelengths (Chiang & Goldreich, 1997). Therefore, the spectrum given by Equation (1) may be a crude (albeit useful) approximation.

The entire spectrum of an optically thick, multi-temperature disk has three distinct regions. For the flat disk spectrum given by Equation (1), these are a Rayleigh-Jeans tail at $h\nu\ll kT_{\rm out}$, where $f_{\nu}\propto\nu^{2}$, a PL part given by Equation (2), and a high-frequency tail resembling a Wien spectrum, $f_{\nu}\propto\nu^{2}\exp(-h\nu/kT_{\rm in})$ at $h\nu\gg(2/\beta-1)kT_{\rm in}$. One could assume that the high-frequency part might correspond to the cutoff of the optical spectrum between the $V$ and $B$ bands (Hulleman et al., 2004), seen in Figure 4. However, such an interpretation requires a rather high temperature $T_{\rm in}\sim 7000$ K, which is well above the dust sublimation temperature for grains of any chemical composition (see e.g., Koch-Miramond et al. 2002; Temple et al. 2021). This means that it would be a gaseous disk. Ertan et al. (2007) argue that the Spitzer NIR data, together with the data from preceding NIR-optical observations, can be interpreted as emission of a gaseous fallback disk with temperatures $T_{\rm in}\approx 6500$ K, $T_{\rm out}\approx 360$ K, and radii $R_{\rm in}\approx 7\times 10^{9}$ cm, $R_{\rm out}\approx 1.9\times 10^{12}$ cm. The $T(r)$ dependence provided in their Table 1 can be approximated by a PL, $T(r)\propto r^{-0.51}$, which leads to the spectral slope $\alpha=0.92$, virtually coinciding with our measurement. Ertan & Cheng (2004) and Ertan et al. (2007) discuss how such a disk could be formed and how its emission could be strongly pulsed.

Looking at Figure 5, one might assume that the high-frequency cutoff of the disk spectrum corresponds to the $K$ band, which would not require such high values of $T_{\rm in}$. However, such an assumption does not look plausible because the $K$ band is in between the simultaneously observed NIRCam 2.5 $\mu$m and 1.4 $\mu$m points that lie on the extrapolation of the LRS spectrum towards NIR. Although it is somewhat suspicious that all the $K$-band measurements are below the LRS spectrum, it is likely due to a chance coincidence.

One of appealing features of the disk interpretation of the (N)IR spectrum is that it could successfully capture the MIPS upper limit. Therefore, we also attempted to fit the disk model given by Equation (1) to the MIRI LRS spectrum, including the $2\sigma$ MIPS upper limit. We find, however, that this disk model is not able to account for the MIPS upper limit while also matching the slope of the LRS spectrum, particularly the points at wavelengths longer than 8$\mu$m. Specifically, the rollover from the Rayleigh-Jeans tail to the flat part of the disk spectrum is not sharp enough to capture the longer wavelength portion of the LRS spectrum, so it undershoots the LRS data. Additionally, the slope of the flat part of this disk model is too soft and under-predicts the NIRCam points. If the MIPS point is ignored in the fit, then the best-fit disk model can successfully capture the MIRI LRS and NIRCam data, but it overpredicts the MIPS upper-limit by a factor of $\sim 2$–3.

It is important to note, however, that at long wavelengths the disk may become optically thin, so its flux at a given $\lambda$ becomes smaller by a factor of $1-\exp({-\tau_{\lambda}})$, where the disk’s optical thickness $\tau_{\lambda}$ decreases with increasing wavelength (e.g., Draine & Lee 1984). Therefore, we cannot exclude the possibility that a fallback disk with $T_{\rm in}/T_{\rm out}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}15$ (hence, $r_{\rm out}/r_{\rm in}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}200$) significantly contributes to at least the IR spectrum of 4U 0142+61.

Overall, variability of nonthermal (magnetospheric) emission seems to be a more natural explanation of the MIPS upper limit discrepancy with the MIRI + NIRCam PL spectrum. However, the true cause of this discrepancy and the origin of the 4U 0142+61’s IR-optical emission can hardly be firmly established without nearly simultaneous observation with JWST and HST.

As we show in Section 2.2.2, the NIRCam and MIRI LRS data in the $\lambda=1.4$–11 $\mu$m band are well fit by an absorbed PL model with $\alpha\simeq 0.96$. Since PL spectra are typical for emission of relativistic particles, this might suggest that the observed IR spectrum comes from the magnetosphere of 4U 0142+61 and can be connected with a nonthermal component of the X-ray spectrum, particularly, with the hard X-ray component, as suggested by Muñoz-Darias et al. (2016). However, the slopes of the IR and hard X-ray spectra are quite different, i.e., they cannot be connected without a large spectral break, $\Delta\alpha\simeq 1.6$ or $\Delta\alpha\simeq 1.1$ for models 1 and 2, respectively (see Figure 7). This suggests that even if the IR originates from the magnetosphere, the IR and hard X-ray emissions either come from two separate particle populations or are produced by different emission mechanisms. Future NIRCam timing observations can help to solidify the location of these IR/NIR emitting particles, the emission mechanism, and, depending on the IR spectrum at longer wavelengths, the strength of IR/NIR pulsations. Additionally, measuring the phase-shifts between the (N)IR and X-ray emission could help to determine the relative locations of the emission sites.