{'fraction_non_alphanumeric': 0.04276167519156525, 'fraction_numerical': 0.043075234679679386, 'mean_word_length': 4.195804907850524, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 1, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discovery of extended infrared emission around the neutron star RX J0806.4-4123

Following up on a faint detection of a near-infrared (NIR) source at the position of the X-ray thermal isolated neutron star RX J0806.4-4123, we present new Hubble Space Telescope observations in the H-band. The NIR source is unambiguously detected with a Vega magnitude of 23.7 +/- 0.2 (flux density of 0.40 +/- 0.06 microJy at lambda =1.54 microm. The source position is coincident with the neutron star position, and the implied NIR flux is strongly in excess of what one would expect from an extrapolation of the optical-UV spectrum of RX J0806.4-4123. The NIR source is extended with a size of at least 0.8arcsec and shows some asymmetry. The conservative upper limit on the flux contribution of a point source is 50%. Emission from gas and dust in the ambient diffuse interstellar medium can be excluded as cause for the extended emission. The source parameters are consistent with an interpretation as either the first NIR-only detected pulsar wind nebula or the first resolved disk around an isolated neutron star.


INTRODUCTION
The X-ray pulsar RX J0806.4-4123 belongs to a group of seven nearby (< 1 kpc), X-ray thermal isolated neutron stars (XTINSs, also dubbed the "Magnificent Seven") whose defining properties are: the lack of detected non-thermal X-ray emission, radio-quietness, relatively large pulse periods (3 − 11 s) and inferred dipole magnetic fields on the order of ∼ 10 13 G. The observational properties of the XTINSs (e.g., Haberl 2013; Kaplan & van Kerkwijk 2009) place this neutron star population between those of the rotation-powered pulsars and the magnetars. Since their X-ray luminosities exceed their spin-down energies by a factor of at least 10, an explanation of the XTINS properties requires an additional energy source other than the energy loss of a rotating dipole. Currently, it is commonly assumed that the most promising model for the heating of magnetars and XTINSs is the decay of their large magnetic fields. Viganò et al. (2013) presented simulations of a magnetothermal evolution model capable of unifying the diverse neutron star populations. In this model, magnetars and XTINSs are hot because they are heated by the decay of their large magnetic fields. In an alternative model, the observational properties of the XTINSs and magnetars are explained by formation of and influence by supernova fallback disks (e.g., Ertan et al. 2014). However, the fallback disk model alone cannot explain the giant outbursts of magnetars. Acceptance of this model is also hampered by the lack of fallback disk detections. So far, there is only one example of a possible (passive) fallback disk around the 3.9 kpc distant anomalous X-ray pulsar 4U 0142+61, detected with Spitzer IRAC by Wang et al. (2006). Ertan et al. (2007) reported that an active irradiated disk model can also reproduce the observed optical and infrared fluxes of 4U 0142+61. RX J0806.4-4123 has the second-longest period among the XTINSs, P = 11.37 s. Its spin-down power isĖ = 1.6 × 10 30 erg s −1 , its inferred dipole magnetic field is B = 2.5 × 10 13 G, and its total X-ray luminosity at a distance of 250 pc is L X = 2.6 × 10 31 erg s −1 (Kaplan & van Kerkwijk 2009;Haberl 2007;Posselt et al. 2007;Haberl et al. 2004). Previously, we used the Very Large Telescope (VLT) in the near-infrared (NIR) to carry out a direct-imaging search for substellar companions around isolated XTINSs (Posselt et al. 2009), and we noticed a slight H-band flux enhancement at the location of RX J0806.4-4123. Motivated by an additional Herschel 160 µm detection very close to the neutron star position (Posselt et al. 2014), we followed up on the NIR source with the Gemini tele-  Kaplan et al. (2011), and the red cross marks the 2002 Chandra position that has a positional uncertainty of 0. 6 (Haberl et al. 2004). scope. We indeed found indications of the VLT source in another NIR band. The two independent faint flux enhancements together, H VLT ≈ 23.3 ± 0.5(1σ) and J Gemini ≈ 24.8 ± 0.5(1σ), resulted in a combined NIR detection (H and J bands) significance of 3.1σ (Posselt & Luhman 2016). From the UV-optical spectral slope (measured by Kaplan et al. 2011), one would expect the XTINS to have a NIR magnitude on the order of 28, much fainter than the NIR detection. Possible interpretations of the NIR flux enhancements could be a substellar companion, a fallback disk or even a very unusual pulsar wind nebula. In this paper, we report on the results of recently obtained Hubble Space Telescope (HST ) NIR imaging follow-up observations that were carried out to confirm the previous weak NIR detection.

Hubble Space Telescope data
We obtained NIR images (program GO-14745) of RX J0806.4-4123 with the HST Wide Field Camera 3 (WFC3; Kimble et al. 2008) on 2016 September 11 (MJD 57642.98), employing the F160W filter with a pivot wavelength of 1.537µm. The total exposure time was 5417 s. We used the WIDE-6 POS-TARGs dither pattern (optimizes the subsampling of the pixels; HST Instrument Science Report 2016-14 1 ) for our 6 individual exposures in MULTIACCUM mode (linear SPARS timing sequence). We stacked and processed our images The target, STmag=27.4, R maxmin = 0.63 s4, STmag= 27.5, R maxmin = 0.17 Figure 2. The comparison of images before (only sky subtraction) and after PSF (and sky) subtraction. We show the target (left) and a similarly faint source (right). The images employ the same limits, logarithmic scale and color map. Rmaxmin is defined as the ratio of the maximum−minimum in the PSF-subtracted image (i.e., the residual range; right subpanels) to the maximum−minimum in the image without PSF subtraction (left subpanels) and was obtained from image statistics in the inner 0. 7 × 0. 7 regions. For the target, Rmaxmin is about a factor 4 larger than it is for s4, underlining that there is substantial "residual flux" after PSF subtraction in case of the target. within PyRAF using AstroDrizzle (version 2.1.8) of the DrizzlePAC software (Gonzaga et al. 2012;Fruchter et al. 2010). Since our target source is very faint, we used the inverse variance map (IVM) weighting scheme for the final combination of the data (Gonzaga et al. 2012). We experimented with the drizzle parameters for a balance of spatial resolution and sampling noise. We chose a final pixel scale of 0. 0642 pixel −1 and a pixel fraction of 0.7. We registered our astrometry to the Gaia data release 1 (Gaia Collaboration et al. 2016a,b); employing the common 117 reference sources in the field of view, the rms of our astrometry fit is 0. 019. RX J0806.4-4123 was observed in X-rays in 2002, 2005, and 2015 with Chandra. From these observations the proper motion was found to be rather small, < 50 mas/yr, (Haberl et al. 2004;Motch et al. 2009;Posselt et al., in prep.) RX J0806.4-4123 was also observed and detected with HST in the UV (ACS/SBC F140LP; STmag= 23.61 ± 0.11) and in the optical (ACS/WFC F475W, STmag= 27.92 ± 0.22) in (Kaplan et al. 2011. The neutron star's position in X-rays and in particular in the optical F475W band are consistent with an unambiguous HST detection in the F160W band, see Figure 1. To measure the source flux we used aperture measurements with a sky annulus for the background and noise estimate. We used the daophot and photcal IRAF packages. In order to account for systematic errors and an apparent small extension, we measured values for four different aperture radii (r = 0. 4, 0. 5, 0. 6, 0. 7), six different sky annuli (inner radius from r = 0. 96 to 1. 6 with different widths ranging from r = 0. 6 to 1. 5 where the latter is limited by a nearby diffraction spike). We also checked 10 slightly different aperture center positions (maximal difference ∼ 0. 3). The largest systematic effect on the measured magnitude came from increasing the aperture radius from 0. 4 to 0. 6 (up to ∼ 0.3 brighter magnitude), indicating that the source may have an extension of up to r = 1. 2. The effect of changing the aperture center position reached ∼ 0.1 in magnitude. The magnitude effect of different sky annuli was much smaller (∼ 0.03). Considering these different measurements, we obtained the following median values and 1σ uncertainties (including systematic errors), a Vega magnitude of 23.7 ± 0.2, an ST magnitude of 27.1 ± 0.2, corresponding to a flux density of 0.40 ± 0.06 µJy.
To test the hypothesis that the source is extended, we subtracted the point spread function (PSF) from the emission. In a drizzled image, pixels are correlated due to the drizzling process, and different PSFs are expected for individual sources due to the different sampling and signal-to-noise properties of each object (Anderson 2014). Therefore, we apply the "bundle" approach by Anderson (2014) which uses the flat-fielded flt images derived from each of the 6 exposures. Accurate (4 times supersampled) "effective" PSFs were constructed for the flt domain and are available at the Space Telescope Science Institute's WFC3 webpages 2 (Anderson 2016). We chose 12 apparent isolated point sources as reference sources. Most of them are similarly faint as the target, but we also included a few moderately bright objects not afflicted by diffraction spikes. Using a pixel raster of 41 × 41, we derived bundles for the target as well as for the 12 reference sources. We created stack images of each source (i) after sky subtraction but without PSF subtraction, and (ii) after PSF and sky subtraction. Comparing these 'before' and 'after' images allowed us to assess the goodness of the PSF subtraction. In general, the flt-based PSF subtraction works well. In the resulting images, residuals are visually negligible for objects as faint as the target. For brighter objects, the relative residuals (with respect to unsubtracted image) are negligible as well. We refer to the Appendix for the study of several reference sources. There are clear extended-emission residuals in the case of the target, see Figure 2 in which we also show the PSF-subtraction result for a similarly faint source for comparison.
Our results led us to conclude that the NIR source at the location of RX J0806.4-4123 is definitely an extended one. The maximum extension of the target source is reached at a position angle of ∼ 20 • East of North. Along the line of maximum extension, the angular size of the full major axis is at least 0. 8. At r = 0. 4 from the center, i.e., at a distance equal to the semimajor axis, the flux has decreased to 23% of the peak flux value (for comparison, a similar flux decrease for a reference point source, e.g., s4 in Figure 7, results in an "extension" of 0. 3). In the perpendicular direction, the full minor axis was measured to be ≈ 0. 5. If one assumes that there is indeed an additional point source, the flux of that point source represents ∼ 20% of the total F160W flux according to our PSF subtraction. In order to estimate a firm upper limit of the contribution of a possible point source to the total flux, we slightly varied the position of the subtracted point source. Using 18 positions within r = 0. 15 of the maximum NIR flux, we estimated a standard deviation of 2.5% in flux. In addition, we manually increased the subtracted point source flux from the fit value to larger ones until the residual "hole" was inconsistent with background fluctuations of the source-free region. This criterion was met when the point source constitutes 42% of the total flux. Therefore, the estimated conservative upper limit on the contribution of a point source is 50%. Turnaround observing mode. The program ID was GS-2016B-FT-23. We obtained 492, 12, and 226 dithered exposures with integration times of 10 s, 12 s, and 15 s, respectively, resulting in a total exposure time of 2.3 h. The field of view of the Hawaii-2 array has a diameter of 6. 1 and a pixel scale of 0. 18 per pixel at the f/16 telescope focal surface. The airmass of the individual exposure ranged from 1.02 to 1.32, the average airmass was 1.08 For the data reduction, we used the f2 package within the Gemini IRAF Package (version 1.12) with IRAF (version 2.16). We applied the usual corrections (dark current subtraction, bad pixel mask correction, flatfielding). After constructing sky frames with nisky, the flat-divided sky frames were subtracted from the respective flat-divided science images. Taking into account the different exposure times, we combined the data from individual nights as well as the data from all nights. Using the Gaia data release 1 (Gaia Collaboration et al. 2016a,b) and the Graphical Astronomy and Image Analysis Tool (GAIA) of the JAC Starlink Project (Draper et al. 2016), we calibrated the astrometry of the reduced images to an absolute astrometric accuracy of 0. 13(1σ).
We selected 16 stars (magnitudes K s = 11.92 − 14.62) from the 2MASS point source catalog (Cutri et al. 2003) with quality flag AAA and no obvious nearby neighbors in the combined image to establish a photometry scale. We applied the daophot and photcal IRAF packages and used an aperture radius of 7 pix for the 2MASS stars (the FWHM is about 3 pix). We determined a zero point of 25.50 ± 0.02. With this calibration, all but the three brightest sources (brighter than 12.35) have measured flux values within 2.2σ of their 2MASS magnitudes. The reason for the deviation of the brightest reference source is likely a slightly nonlinear behavior of the detector in some of the longer exposures due to imperfect screening. Based on the other 13 calibration stars, reliable K s magnitudes are 13.0 and fainter.
There is no significant K s source at the position of RX J0806.4-4123, though there may be a very faint (less than 1σ) flux enhancement (see Figure 3). Formally, we measured the flux in an aperture with radius of 7 pix and obtained a magnitude 22.4±0.4 (1σ), resulting in an upper 3σ limit of K s < 21.2 (corresponding to < 2.2µJy) at the position of the neutron star. Such an upper limit is consistent with ∼ 3σ detections in the surrounding area (radius r ∼ 30 ), where such faint sources are found to have magnitudes of 20.8 to 21.4, depending on the varying background emission.

DISCUSSION
The NIR source at the location of RX J0806.4-4123 is unambiguously confirmed with our HST F160W observations. The measured F160W flux, 0.40 ± 0.06 µJy, is clearly in excess of the expected neutron star flux considering the extrapolation of its UV-optical flux (see Figure 4). In addition, the object has a size of at least 0. 8 (Section 2) and an inhomogeneous brightness distribution, consisting of a brighter 'core' and elongated 'halo' emission ( Figure 2). In order to estimate the likelihood of an unassociated background source, we measure a source density of 0.07 sources per square arcsec in the F160W image which includes all sources with STmag≤ 27.4 (conservative limit since the target has STmag= 27.1 ± 0.2). The separation of the flux peak of the extended emission from the neutron star position in the 2010 HST F475W image is 0. 15. Recent Chandra observations indicate a proper motion of RX J0806.4-4123 of ≈ 20 mas yr −1 (Posselt et al., in prep.; the previous limit was µ < 86 mas yr −1 (2σ), Motch et al. 2009). Considering a (conservative) total position uncertainty of 0. 4 (including proper motion 2σ < 50 mas yr −1 , F475W position error < 0. 1), and N = 1 trials, we estimate a by-chance probability of 3.7% 3 . Therefore, we regard an unassociated back-ground source as unlikely.
Even though the by-chance probability for an unassociated background source is low, a background galaxy might explain the detected NIR emission. The optical and UV data obtained with HST ( Figure 2 by Kaplan et al. 2011) do not show extended emission, and the fluxes are consistent with those measured for the other six of the Magnificent Seven, i.e., it is reasonable to assume that there is no significant flux contribution from a galaxy to the F475W and L140P flux measurements of RX J0806.4−4123. The colors of galaxies indicate that any putative galaxy counterpart must be either very red (e.g. Daddi et al. 2000;Franx et al. 2003;van Dokkum et al. 2006;Guo et al. 2012), or, considering the lack of F475W flux, be at redshift z 3.5 (e.g. Giavalisco 2002). The same result follows from simulations using the mock galaxy catalog by Williams et al. (2018). These simulations also show that the probability to detect galaxies with consistently red observed colors (either star-forming or non-starforming) in the field is low -only ∼ 0.2 per square arcmin are expected to have F160W fluxes within the measured flux uncertainties while also being sufficiently faint in F475W. The analysis of actual, sufficiently deep observations results in similar numbers. Analysis of the data on the Hubble Ultra Deep Field (HUDF; Rafelski et al. 2015) yield 0.39+/−0.17 galaxies per square arcmin if similarly red and faint sources as our target are considered. The uncertainty results from the large Poisson error in the very small 12.8 arcmin 2 field of view. 4 Furthermore, we note that high-redshift galaxies are in general more compact, having smaller sizes than low-redshift galaxies (Shibuya et al. 2015;van der Wel et al. 2014). The detected ∼ 1 NIR extension would indicate an unusually large high-redshift galaxy. In the mock galaxy simulations, for example, < 5% of galaxies with consistent observed colors are large enough to also be consistent with the detected NIR extension. Thus, we conclude that a highz galaxy is an unlikely counterpart of the extended NIR emission at the position of RX J0806.4−4123. The detected extended emission could also be due to coincidence with emission of the interstellar medium (ISM). Scattered stellar radiation, hot dust, or line emission from the ISM gas are potential origins. We discuss these possibilities in more detail in Section 3.1.
If the NIR emission is associated with RX J0806.4−4123, then L F160W ∼ 6 × 10 27 erg s −1 , assuming a distance of 250 pc. Using the X-ray flux and spin-down power of the NS, the NIR efficiency is η F160W = L F160W /Ė = 4 × 10 −3 , and the luminosity ratio L F160W /L X = 3 × 10 −4 . The F160W flux value, 0.40 ± 0.06 µJy, and a size of 0. 5 × 0. 8 translate into an average brightness temperature 5 T B = 282 K. If our source were a blackbody with this temperature, then one would expect a very high IR flux, e.g., F 4.5µm = 49 Jy which is clearly above the current flux limits by Spitzer. However, thermal emission with a temperature, T = 950 K and a filling factor 6 of 7.3 × 10 −11 would be consistent with the F160W and J-band flux measurements as well as the upper limits obtained for K S and 4.5µm. It follows from these estimates that thermal emission is, in principle, consistent with the (N)IR measurements if the emitting area is a very small fraction of the apparent area of the detected extended source. Thermal emission could come from warm gas or X-ray irradiation and subsequent cooling of solid bodies such as dust or a substellar companion. At a distance of 250 pc, an extension of ∼ 0. 8 corresponds to 3 × 10 15 cm or 200 AU. This rules out a substellar companion as the sole origin of the NIR excess. However, the existence of a warm substellar companion cannot be entirely ruled out since the conservative upper limit on the contribution of a point source to the total flux is 50%. Such flux limit translates into a substellar companion mass limit of < 10 Jupiter masses (using a bolometric correction of 3.5 (Golimowski et al. 2004), the calculations presented by Burrows et al. (2001), assuming a distance of 250 pc and an age of 10 Myr). Warm gas and heated dust, either in the ISM or cirumstellar matter, are discussed in Sections 3.1 and 3.2.
Alternatively, the F160W measurement can also be explained by non-thermal emission. Assuming a power law F ν ∝ ν α for non-thermal emission, the current spectral constraints imply either a rather steep slope α < −2.5, or an exponential cutoff between the F160W and F475W bands. Charged particles producing curvature and/or synchrotron emission would be the likely source. With respect to an association to a neutron star, the most obvious origins are magnetospheric emission (for the possible flux contribution from a point source) 5 The brightness temperature T B is defined from the equation where Fν is the flux density measured in the solid angle Ω subtended by the source, and Bν (T B ) is the Planck function. 6 The filling factor is the ratio of projected emitting area and projected source area. or a pulsar wind and the shocks it produces. We discuss such a pulsar wind nebula in Section 3.3.

Emission from the ISM ?
In principle, an ISM density enhancement along the line of sight could scatter NIR emission from ambient light and produce the observed extended NIR emission. The so-called cloudshine is a large-scale (typically, several arcminutes) phenomenon that is caused by ambient star light around large and dense molecular clouds (Foster & Goodman 2006;Lehtinen & Mattila 1996). There are no other extended NIR features and no indication of a molecular cloud (e.g., Dobashi et al. 2005;Dame et al. 2001) in the field of RX J0806.4-4123 which one would expect if the NIR emission were due to cloudshine. Instead of the interstellar radiation field and a large, dense molecular cloud, one could postulate that the NIR emission could originate from local point source emission that is scattered by a small (i.e., unaccounted for) ISM enhancement at a slightly larger distance than the one of RX J0806.4-4123 7 . However, in this case the expected thermal NIR flux from the XTINS surface alone is too low to produce the large observed NIR excess by scattering. An additional NIR emitter would be required as well as an unlikely "wall" in the ISM.
If the NIR excess flux is due to genuine emission from the ISM, it could originate either from ISM gas or dust in the diffuse ISM or in an ISM density enhancement. There are a few emission lines of ISM gas in the H and the other NIR-bands. Strongly forbidden rotational-vibrational emission lines of H 2 are observed, for example, in reflection nebulae or shocked regions in star forming regions and around AGB stars (e.g., Shull & Beckwith 1982). These H 2 emission lines are due to collisional excitation or UV fluorescence, and in accordance with their detection in dense environments they require high densities and significant molecular fractions (Black & van Dishoeck 1987;Hollenbach & McKee 1989). The low N H value of RX J0806.4-4123 implies not only a low density, but also a very low 7 The absorbing hydrogen column density inferred from the Xray spectrum of RX J0806.4-4123 is between N H = (0.4 ± 0.1) × 10 20 cm −2 (Haberl et al. 2004) and N H = (1.7 ± 0.2) × 10 20 cm −2 (Kaplan & van Kerkwijk 2009), depending on the spectral model. Thus, using the relation between N H and A V by Foight et al. (2016), the extinction of the neutron star, A V < 0.1, rules out a residence in dense ISM. The color excess along the line of sight toward the neutron star is only E(B − V ) = 0.14 ± 0.08 (1σ) even at a distance of 1.5 kpc (well beyond the likely distance of the XTINS) according to the a recent 3D map of the ISM by Capitanio et al. (2017) (their spatial resolution along the line of sight is ∼ 200 pc at 1 − 1.5 kpc). molecular fraction, 2N H2 (2N H2 + N HI ) −1 < 10 −4 , because inefficient self-shielding from cosmic rays leads to dissociation of H 2 (Burgh et al. 2007;Savage et al. 1977). Therefore, the NIR excess flux around the neutron star is unlikely to be produced by H 2 emission lines.
Interstellar dust emission peaks in the far infrared. It can also emit NIR photons if it is hot enough. The required high temperatures (1000 −2000 K) in thermal equilibrium could only be reached very close to the XTINS (unresolved by HST ) by typical (size ∼ 0.1µm) dust grains. A non-equilibrium process, the so-called stochastic heating of very small grains 10Å or large polycyclic aromatic hydrogen (PAH) molecules by UV photons, was found to produce the large-scale NIR emission of reflection nebulae (e.g., Draine & Li 2001;Sellgren 1984). There are two reasons why very small grains/PAH molecules of the ISM cannot be responsible for the extended emission around RX J0806.4-4123. First, there is not a sufficient amount of very small grains in the ISM surrounding the XTINS to produce the observed NIR flux. Using the constraints on the XTINS spectrum in UV (Kaplan et al. 2011), and our flux measurement at 1.6µm, one can estimate the required number density, n vsg , of very small grains (∼ 10Å) which are stochastically heated by ∼ 10 eV UV photons (Sellgren 1984). The value, n vsg ∼ 10 −3 cm −3 , is at least four orders of magnitude larger than what one can expect for the typical ISM with a hydrogen number density n H ∼ 1 cm −3 , dust-to-gas ratio of 6%, and typical grain size and composition distribution (all from, e.g., Zubko et al. 2004; we used their BARE-GR-S dust model which matches observations of the ISM in the solar neighborhood). The lack of known "NIR halos" around other strong UV sources such as white dwarfs supports our estimate that the diffuse ISM does not provide enough small grains for substantially enlarged NIR flux close to these objects. The second reason why very small ISM grains are an unlikely origin for the observed NIR emission is the strong X-ray emission of the XTINS. Voit (1991) estimated that X-ray photons efficiently evaporate grains of sizes 10Å. Overall, we conclude that the observed NIR emission around RX J0806.4-4123 is unlikely to come from either gas or dust of the photonheated ISM. We discuss the interaction of potential pulsar wind particles with the ISM in Section 3.3.

Emission from circumstellar material?
If the NIR emission comes from the immediate vicinity of RX J0806.4-4123, it could be due to a surrounding disk or torus structure, similar to what has been inferred for the magnetar 4U 0142+61 from the Spitzer detection by Wang et al. (2006) or for white dwarfs from NIR detections (e.g., Melis et al. 2010). The stochastic heating of very small grains, which was discussed in the previous section for the ISM, is a possibility for circumstellar matter, too. However, close to a neutron star, one expects quick destruction (evaporation by X-rays) and removal (by Poynting-Robertson drag; e.g., Figure  10 by Posselt et al. 2014) of such small grains. Hence, regular replenishment, e.g., by collisions of asteroids or large grains, would be required to keep such a 'debris disk' bright in the NIR.
Considering alternatively an active fallback disk model, Ertan et al. (2017Ertan et al. ( , 2014 showed that such a disk is able to produce the X-ray properties and the optical excess of the XTINSs. Here, we use their model to check whether such a fallback disk can reproduce the measured NIR excess. The model employs dissipative heating which is most relevant in the active inner part of the disk, for RX J0806.4-4123, from the inner disk radius, r in ∼ 10 9 cm to 5 × 10 11 cm. The model also considers X-ray irradiation which is most relevant for the passive outer part of the disk. For details on the model and the model parameters for RX J0806.4-4123, we refer to Ertan et al. (2017Ertan et al. ( , 2014. The outer radius of the passive disk is not well known and depends on the initial disk conditions and the Xray luminosity in the early phase of evolution. However, the local irradiation temperature is very small for large radii -for RX J0806.4-4123, it decreases to < 30 K for r > 10 13 cm. For studies of the flux at infrared and shorter wavelengths, it is sufficient to use an outer radius r out ∼ 10 13 cm. The extension of the NIR source at the distance of RX J0806.4-4123 is r ∼ 10 15 cm (or ∼ 100 AU). This raises the question of whether there is any possibility for the disk model to explain not only the observed NIR flux, but also the observed source extension. Drawing on the findings in protoplanetary disk studies, favorable properties of the disk and viewing geometry can perhaps provide such an explanation. For protoplanetary disks, it was shown that in the case of a flared disk NIR emission can be scattered and produce extended emission (e.g., Anthonioz et al. 2015;Mulders et al. 2013). First, we check whether the disk model can reproduce the observed spectral energy distribution. Then we outline a potential explanation for the observed source extension.
Using r in = 10 9 cm and r out = 10 13 cm, we obtain the model spectrum given in Figure 4. The disk model curve in Figure 4 is within the error range of the data, illustrating that the fallback disk model can reproduce the observed NIR flux. Within this disk model, the work done by the magnetic torque heats the inner rim of the disk which then produces the optical and UV fluxes. In contrast, the NIR emission of the model spectrum is produced from the surface of the irradiated inner disk within a few times 10 9 cm of the XTINS. One implication is that a completely passive disk could also produce the NIR emission, but not the optical-UV emission.
One possibility to explain the observed source extension is based on the assumption of disk flaring and scattering off the outer passive disk surface. The thickness of a disk can be described with the pressure scale height, h c s /Ω K , where c s is the sound speed and Ω K is the local Keplerian angular velocity of the disk. In flared disks the flaring is described by the aspect ratio, h/r. The outer radii of protoplanetary disks around T Tauri stars reach several hundred AU and the h/r values can approach unity (see, e.g., Dullemond & Monnier 2010; Armitage 2010; Chiang & Goldreich 1997 for reviews). Chiang & Goldreich (1997) (CG97) calculated self-consistently the aspect ratio and the temperature profiles of these disks. We follow the result of CG97 for our estimates, but adapt the model parameters (luminosity, mass) to the properties of RX J0806.4-4123. CG97 found h/r ∝ r 2/7 for an optically thick inner disk region. Outside this region there is a radially isothermal optically thin region along which the aspect ratio increases even faster with increasing r as h/r ∝ r 1/2 . The disk photosphere, defined as the disk surface layer corresponding to the optical depth of unity to the irradiating flux, is estimated by CG97 to be H 4h. In the luminosity regime of RX J0806.4-4123, L * ∼ L X = 2 × 10 31 erg s −1 , using h/r ∝ L 1/7 * M −4/7 (CG97) with M = 1.4M , we estimate the aspect ratio of the disk for RX J0806.4-4123 as H/r 2.7×10 −2 r 2/7 AU where r AU is the radial distance in AU. For r AU = 100, we find H/r 0.1. This could be taken as a lower limit, since the outermost part of the passive disk could be optically thin to its own emission, along which H increases faster with increasing r, as r 1/2 . In CG97, the radially isothermal disk region has a temperature ∼ 20 K. For RXJ0806.4-4123, such a temperature is reached at r ∼ 8 × 10 12 cm. Thus, a large fraction of the outer passive disk of RXJ0806.4-4123 could be radially isothermal with H/r ∝ r 1/2 . In this case, the result is H/r ≈ 0.3 at 100 AU. Hence, depending on the extension of the radially isothermal region, H/r could be between 0.1 and 0.3 at the outermost disk. For our estimates below we take H/r = 0.2 at r out 100 AU.
The NIR luminosity illuminating and being scattered off the outer disk surface, L scat , is a small fraction of the total NIR luminosity emitted from the inner disk, L in . As seen from the inner disk, the outer passive disk subtends a solid angle Ω d from the mid-plane of the disk (θ = 90 • ) to the surface of the outer disk (θ = θ d ).
For H/r 0.2 at r 100 AU, this angle is θ d 78 • . Assuming an albedo of β 0.9, we estimate that up to 4 % of the total NIR emission from the inner disk could illuminate and be scattered by the outer disk surface, L scat /L in 0.04.
The observed ratio of the NIR flux scattered off the outer disk to the NIR flux from the inner disk (r ∼ 10 9 cm) depends on the details of the viewing and the scattering geometry, actual h/r profile, temperature and the extension of the currently passive disk. For instance, the scattering geometry matters because the scattered emission is unlikely to be purely isotropic. The angular distribution of photons scattered by a dust grain depends on the grain composition and grain size (Mulders et al. 2013;Draine 2003). Large grains in particular show very anisotropic scattering, where one scattering produces a conical beam with the small solid angle Ω sc . The NIR light emitted from the inner disk is estimated to be scattered from a thin layer of the disk surface. In order to estimate the contribution of the scattered component to the observed NIR flux, we use a simplified picture and represent the total scattered emission by two narrow fan beams opening from the two surfaces of the disk. In the case of single scatterings, the solid angle of one fan beam can be estimated by integration over the azimuthal angle as Ω beam ≈ 2π 1/2 Ω 1/2 sc . Inside this solid angle, the ratio of the observed scattered flux to the observed inner disk flux, F scat /F in , is related to the ratio of the luminosities as where i is the inclination angle between the line of sight and the normal to the mid-plane of the disk. If we assume, for example, Ω beam ∼ 1 sr, and cos i = 0.25 (i.e., i 75 • ), we obtain F scat /F in ≈ 1 and F in /F tot ≈ 0.5 which are consistent with the observation. For this scenario, half of the H-band emission is produced at the (unresolved) inner disk, while the remaining half is the contribution from the scattered light. In this example, the inclination of the observer and the inclination of the scattering disk surface are very similar, and the observer remains inside the beaming solid angle of the near-side disk surface. Instead of the simple fan-beam geometry for the scattered NIR radiation, the NIR light could actually be scattered in all directions. Our explanation remains valid if the scattered flux density is still con-centrated near the disk surface. If the disk is observed at small inclination angles (close to face-on view), the contribution of the scattered light to the observed flux decreases depending on the efficiency of the "beaming", while the flux from the inner disk increases. In other words, the observed manifestations of any large disks around XTINSs could vary significantly depending on the optical properties of the disks, the inclination cos i and the observed wavelength.
Equation 1 shows that for an observer inside the solid angle of the beamed scattering, Ω beam , the observed F scat /F in could be significantly greater than L scat /L in . In our simplistic model of an assumed scattering anisotropy, the observer sees not the entire disk but only a fraction of the disk surface close to the line of sight, which, at a given outer radius, reduces the size of the observed extended emission in comparison to the actual size of the disk. However, the observed size increases with increasing contribution from multiple scatterings. A given light ray can go through multiple scatterings, both local and at different radii, which could significantly change its azimuthal direction. The details of the extended emission depend on the actual geometry and the scattering properties of the disk and require detailed modeling to quantify the "beaming" of the NIR emission and its observable azimuthal distribution. We also note that we extrapolated the results from protoplanetary disks to our case of a potential neutron star disk. The composition of such a disk is very likely quite different from that of a protoplanetary disk. For example, one could expect a higher metallicity and possibly a higher dust-to-gas ratio in supernova fallback material.
In addition to thermal continuum, such a disk could also show line emission. Modeling of different disk geometries, densities and compositions is beyond the scope of this paper, in particular because we currently have a good measurement in only one photometric band of the (N)IR.
Alternatively to a disk with favorable geometry, scattering of the emission from the inner disk by the dust grains in an additional outer cold dusty belt or torus structure around RX J0806.4-4123 could also explain the extended NIR emission. Our previous Herschel 160µm detection for the region of RX J0806.4-4123 had resulted in an estimate of dust grain locations at r d = 2.3 × 10 16 a −1/2 µm cm to 2.2 × 10 15 a −1/2 µm cm for a reasonable dust temperature range of 10 to 22 K and in dependence of the grain size, a µm , in µm (Posselt et al. 2014). We had speculated about the existence of a belt structure to explain the large ratio of the lumi-nosity at 160µm to the XTINS' X-ray luminosity (for details, see Posselt et al. 2014). The association of the Herschel 160µm detection with RX J0806.4-4123 is not as firm as the Hubble NIR detection because of the worse spatial resolution of Herschel. Nevertheless, it is intriguing that the spatial dimensions of the potential torus are of about the same order as the NIR extension. A dusty torus alone cannot explain the NIR emission because the scattered NIR emission from RX J0806.4-4123 would be expected to be much fainter. The implied two-component structure is not uncommon for debris disks around main sequence stars (e.g., Su et al. 2013) and may also be a possibility around a neutron star.
Overall, a disk could be the source of the detected NIR flux, while the observed source extension could be due to scattering by cold dust grains located at larger separations from the neutron star, either in a disk with favorable properties or in a surrounding dusty torus.

A Pulsar Wind Nebula?
Extended emission around a pulsar could also be a pulsar wind nebula (PWN). Most PWNe are known from high-energy and radio observations (e.g., Reynolds et al. 2017), about a dozen have been detected in H α (Brownsberger & Romani 2014), two have recently been discovered in the ultraviolet (e.g., Rangelov et al. 2017). A few high-energy PWNe were also detected in the (near)-infrared, e.g., the PWNe of Vela (Shibanov et al. 2003) and PSR B0540-69 (Mignani et al. 2012), but there is only one serendipitous discovery of an infraredonly PWN with Spitzer at wavelengths longer than 8µm (Wang et al. 2013). The NIR emission around RX J0806.4-4123 could come either from line emission in the shocked ambient medium or from synchrotron emission of the shocked pulsar wind. In the case of shocked ambient medium, possible origins are emission lines of H 2 and [Feii]. These shock-excited lines are frequently observed in the dense environments of star-forming regions (e.g., Herbig-Haro objects), evolved stars with eruptive mass loss, and supernova remnants (e.g., Lee et al. 2014;Reipurth & Bally 2001). Such a dense environment, however, is unlikely for RX J0806.4-4123 based on the X-ray absorption in the XTINS spectrum and 3D models of the ISM (Sec. 3.1). Hence, in the case of the PWN scenario, the emission is probably due to synchrotron emission from the shocked pulsar wind.
At lowĖ, particles cannot be accelerated to sufficiently high energies to produce X-ray emission, but emission at lower energies could be possible. According to Reynolds et al. (2017), the maximum synchrotron photon energy E max of the accelerated electrons in the shocked pulsar wind can be estimated as where ζ ∼ 1 is a numerical factor,Ė 30 is the spin-down power in units of 10 30 erg s −1 , B −6 is the magnetic field at the location of the shocked pulsar wind in units of µG, and σ is the magnetization which is defined as the ratio of the Poynting flux to the particle enthalpy flux. Expected PWN magnetic field values are of the order of 10 − 100µG depending on distance, while σ is unknown but usually assumed to be 1 or even σ 1 (Reynolds et al. 2017;Kennel & Coroniti 1984). The spin-down power of RX J0806.4-4123 isĖ 30 = 1.6 (Kaplan & van Kerkwijk 2009). Thus, E max ∼ 1 eV (corresponding to λ ≈ 1.2µm) is possible for σ 1.
The X-ray efficiencies, η X = L X /Ė of observed PWNe vary greatly, reaching values between 10 −6 and 10 −2 (Reynolds et al. 2017). The similarly defined NIR efficiency of RX J0806.4-4123, η F160W = 4 × 10 −3 , is within the PWN efficiency range seen at other wavelengths. For a pulsar moving supersonically with the total velocity v PSR through the ambient medium with density ρ a , the stand-off radius R s of the bow shock apex can be described by R 2 s =Ė PW f (4πc ρ a v 2 PSR ) −1 , wherė E PW = ξ wĖ Ė is the fraction of spin-down power carried away by the wind, and the factor f takes into account a possible anisotropy of the pulsar wind (f = 1 for an isotropic wind). Using an observed source extension radius r ∼ 0. 4 which corresponds to R s ∼ 1.5 × 10 15 cm (Section 2), ξ w f 1, and ρ a ≈ n a m p (m p is the mass of a proton), we obtain where v ⊥ is the transverse velocity component. The proper motion of RX J0806.4-4123 is unusually low, µ < 86 mas yr −1 (2σ), with recent Chandra observations indicating ≈ 20 mas yr −1 (Posselt et al., in prep.). At a distance of 250 pc, these values correspond to v ⊥ < 100 km s −1 (2σ) and 24 km s −1 , respectively. This does not significantly exceed the speed of sound in the typical diffuse ISM (c s ∼ 10 km s −1 for T ∼ 8000 K). The constraint on the velocity and number density in Equation 3 depends on the actual R s , but it is uncomfortably low, particularly considering that an anisotropic wind and a typically assumed fraction of the spin-down power emitted with the wind (ξ w ∼ 0.5) will lower the value further. R s could, however, be smaller since it is difficult to establish the exact location of the pulsar in the extended NIR emission. Hence, a PWN can currently not be ruled out as a possible origin of this emission.
If RX J0806.4-4123 were a usual rotation-powered pulsar, the PWN would be an obvious explanation. Such an explanation is, however, somewhat surprising for a member of the Magnificent Seven which show only thermal X-ray and no radio or γ-ray emission. A part of the detected NIR emission (the point-source contribution that can be at maximum 50%) could in principle come from the magnetosphere of RX J0806.4-4123. However, the optical-UV spectrum of the neutron star, F ν ∝ ν α with α = 1.62 ± 0.20(1σ) (Kaplan et al. 2011), is consistent with, and thus likely dominated by, a Rayleigh-Jeans tail of thermal emission (for which α = 2). If there is a PWN, then a lack of detectable non-thermal emission from the neutron star at higher energies than NIR is unprecedented. It is noteworthy that the XTINS spectra of other members of the Magnificent Seven, in particular RBS 1774 8 , show much stronger deviation from a Rayleigh-Jeans tail, α = 0.53 ± 0.08(1σ) (Kaplan et al. 2011). However, as discussed by Kaplan et al. (2011), a magnetospheric origin of the deviations from Rayleigh-Jeans is not the only possible explanation, in particular considering the high optical luminosities of these XTINSs. The current uncertainties in the models of strongly magnetized neutron star atmospheres at long wavelengths could be, for example, another explanation for the unusual spectra of the Magnificent Seven. If the extended NIR emission around RX J0806.4-4123 is indeed a PWN, it could indicate an interesting new avenue to study the properties of pulsar winds without very energetic particles originating from the (inner) magnetosphere. It is interesting that in comparison to typical spectral slopes of X-ray detected PWNe, α X ≈ −0.5 (Kargaltsev & Pavlov 2008), the slope of the NIR PWN, α PWN −2.5, is rather steep according to Perhaps it also has a viewing geometry favorable for the detection of the fallback disk. Could these properties be the reason why we have seen extended NIR emission only for this neutron star so far? Looking at the Magnificent Seven with comparably deep NIR limits, we compare the two brightest and best studied members with RX J0806.4-4123. RX J1856.6-3754 at a parallactic distance of 123 +11 −15 pc, has a transverse velocity of 192 +17 −28 km s −1 andĖ ∼ 3 × 10 30 erg s −1 , and its H-band limit is 21.5 mag (Tetzlaff et al. 2011;Walter et al. 2010;Posselt et al. 2009;van Kerkwijk & Kaplan 2008). It also has a rather puzzling H α nebula (Brownsberger & Romani 2014;van Kerkwijk & Kaplan 2008;van Kerkwijk & Kulkarni 2001). RX J0720.4-3125 at a parallactic distance of 280 +210 −85 pc, has a transverse velocity of 143 +108 −44 km s −1 andĖ = 4.7 × 10 30 erg s −1 , and its H-band limit is 23.1 mag (Tetzlaff et al. 2011;Kaplan & van Kerkwijk 2005;Posselt et al. 2009). In comparison to the Vega magnitude of 23.7 ± 0.2 measured for RX J0806.4-4123 with WFC3/F160W, the H-band flux limits are shallower for RX J0720.4-3125 (factor 2) and RX J1856.6-3754 (factor 8). The smaller distance of RX J1856.6-3754 is not small enough to counterbalance the effect of its shallower observation, while for RX J0720.4-3125 the large distance uncertainty makes any NIR luminosity limit highly uncertain. Thus, there is currently no indication that the NIR emission around RX J0806.4-4123 is exceptional among the Magnificent Seven.
Reports of NIR-detected PWNe (or disks) are rare for other isolated neutron stars. Rotation-powered pulsars and magnetars are detected in the NIR, but there is neither such a steep slope between the optical and the NIR fluxes nor are the PWNe only detected at these wavelengths (e.g., Danilenko et al. 2011;Durant & van Kerkwijk 2005;Durant et al. 2004). Magnetars are on average at several kpc, hence any small-sized extended emission would be difficult to establish. An example of a nearby rotation-powered pulsar is Geminga which is similarly distant (d = 250 +230 −80 pc) as RX J0806.4-4123, but has a transverse velocity of ≈ 211 km s −1 and a much higher spin-down powerĖ = 3.3 × 10 34 erg s −1 (Verbiest et al. 2012;Faherty et al. 2007;Bertsch et al. 1992). As Figure 6 by Danilenko et al. (2011) shows, Geminga's optical-IR spectrum is relatively flat in F ν though interestingly there may be some indication of a rise toward the mid-infrared. Geminga has a prominent X-ray PWN, and Shibanov et al. (2006) reported a detection of the PWN bow shock in the I band, but not in the HST NICMOS F160W image. From the higherĖ of Geminga, one would expect PWN emission at shorter wavelengths (equation 2). The situation could be differ-ent for rotation-powered pulsars withĖ as low as those of the Magnificent Seven. The old PSR J0108−1431 haṡ E = 5.8 × 10 30 erg s −1 , only slightly larger than the value of RX J0806.4-4123. At a parallactic distance of 240 +124 −61 pc, its transverse velocity is 194 +104 −51 km s −1 and its H-band limit is 21.4 (Deller et al. 2009;Posselt et al. 2009). Again, the relatively shallow observations do not allow a meaningful comparison with RX J0806.4-4123.

CONCLUSIONS
The extended NIR emission around the isolated neutron star RX J0806.4-4123 is unlikely to come from an unrelated field object. It can be interpreted as coming either from a disk with favorable viewing geometry or a PWN created by shocked pulsar wind particles of relatively low energy. A flux contribution of up to 50% could come from a point-like source which could represent the unresolved inner disk or the pulsar magnetosphere in the respective interpretations. The spectra and the details of the spatial shape are expected to be quite different for the two explanations and can be probed with future high-resolution observations with the James Webb Space Telescope.
Most pulsars do not have deep NIR observations, mostly because the emission of the neutron stars at these wavelengths was expected to be just an unexciting extension of the optical-UV spectral slope. RX J0806.4-4123 is a good example that neutron stars keep surprising us. Whether the detected extended emission around RX J0806.4-4123 is a lucky coincidence of low pulsar speed and favorable viewing geometry or this is a common property among (some types of) pulsars needs to be checked with systematic deep NIR surveys.
We thank the referee, Marten van Kerkwijk, for his comments and suggestions which helped to improve the quality of this manuscript. We are very thankful to Jay Anderson for his helpful support in the WFC3/F160W data reduction, in particular his crucial help with respect to an optimal PSF subtraction. We also thank Mario Gennaro and Howard Bond for helping us to optimize our observation setup.
B.P. acknowledges support for this work (program 14745) by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.
U.E. acknowledges research support from TÜBİTAK (The Scientific and Technological Research Council of Turkey) through grants 117F144 and and from Sabancı University.
C.C.W acknowledges support from the National Science Foundation Astronomy and Astrophysics Fellowship grant AST-1701546.
The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium.
This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC,https:// www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.  Figure 6. The goodness of the PSF-subtraction for point sources using the Bundle software on the FLT-images. Rmaxmin is defined as the ratio of the maximum−minimum in the PSF-subtracted image (right subpanels) to the maximum−minimum in the sky-subtracted image (left subpanels). The statistics was obtained in the central region using a square with 0. 7 length. Reference sources 1 to 4 and 12 are plotted using the same limits and the same color maps as for the target. For the other reference sources, we adapted the limits and color maps to better show the larger dynamical ranges.
reference source 9, STmag= 26.5, R maxmin = 0.07 reference source 10, STmag= 26.5, R maxmin = 0.05 reference source 11, STmag= 22.6, R maxmin = 0.01 reference source 12, STmag= 28.7, R maxmin = 0.62 Figure 7. Continuation of Figure 6. Figure 8. Comparison of image statistics in the inner 0. 7 × 0. 7 regions for the (bundle-stack) images with PSF subtraction ("psf+sky sub") and without PSF subtraction ("sky sub") for the reference sources (boxes and triangles) and the target (cross and asterisk symbols). The sizes of the symbols scale with the F160W magnitude. The brightest reference source (smallest symbol) has STmag=19.0, the faintest reference source has STmag=28.7. The latter is the object which is deviating from the other reference sources. Its PSF-subtracted image reveals a potential nearby faint neighbor.