REVISITING THE PARALLAX OF THE ISOLATED NEUTRON STAR RX J185635−3754 USING HST/ACS IMAGING

The initial data processing consists of cleaning the images of cosmic rays and correcting the pixel positions for distortions in the detector.

Cosmic rays are a problem with all CCD detectors. While largely cosmetic, occasionally a cosmic ray does hit within the PSF of a source, and will affect both the astrometry and the photometry. We begin with the flat-fielded _flt.fits images. We reject all pixels with a data quality flag set to 2048 or greater (those identified as saturated pixels and cosmic rays) by setting the data value equal to the median of the surrounding eight pixels. About 2% of the pixels in each image are so affected. This is mostly done to simplify the data reduction code, where large data values can draw off a mean or median. But in the end this is of little consequence. We save this cosmetically corrected image as a fits file.

The ACS is a radial bay instrument on the HST, hence there is significant and asymmetric distortion in the images (Anderson & King 2004). We use Anderson's img2xym_HRC Fortran code to find the stars in the cosmetically corrected images and to correct for the instrumental distortion. This code determines positions by doing a PSF fit using the filter-specific PSF. According to Anderson & King (2004), the distortion correction corrects to better than 0.01 pixel in each coordinate for sufficiently bright stars. The output of this code is a list of raw and corrected X and Y positions in the instrumental frame, along with an instrumental magnitude. The code does not return any estimate of the uncertainty in the position. We adopt the mean pixel scale of 0.02827 arcsec pixel⁻¹.

Using the thresholds we selected (HMIN = 5, FMIN = 150), we identify 21 stars in the field that are common to at least five of the visits, in addition to the neutron star. Thirteen of these stars are recovered in all visits, and five are seen in seven of the visits. The others lie close enough to edge of the detector that instrumental rotation leaves them outside the field of view on occasion. Only the fifth visit contains all 21 stars.

We estimate the measurement uncertainties σ at each iteration of the various fits. We take the measurement error on a particular star to be the standard deviation of the position measurements in the four images obtained during a visit after transforming to a common frame. As the transformation changes, so do these uncertainties. Where a positional measurement disagrees by more than 1 pixel from the mean of the other three, we flag that measurement as bad and ignore it. Most such instances are due to cosmic rays affecting the initial PSF fit. The uncertainties σ vary from star to star and for each visit, but are identical for the four measurements of a star in a particular visit. We set a minimum σ = 0.01 pixel, the accuracy to which Anderson & King (2004) claim to be able to measure the positions. As expected, σ increases with instrumental magnitude (Figure 1). In Figure 1, we also show the formal uncertainties on the positions as extracted using the IDL StarFinder procedure (Diolaiti et al. 2000). The two are in excellent agreement.

Zoom In Zoom Out Reset image size

The initial transformations are accomplished by first shifting the measured positions within a given visit by the nominal dither, and then rotating the (X,  Y) coordinates of the stars to nominal north using the ORIENTAT keyword in the fits images. We wrote three independent χ² minimization astrometric fitting codes, and wrote another incorporating the IDL function MPFITFUN.³ Our fitting codes differ in small details, including how they solve for the parallax ellipse, how they reject outliers, and how they estimate the measurement uncertainties. We refer to these models as A1 through A4 in Table 3. As a test, we also fit the data using a publicly available astrometric package, ATPa (model T). Finally, we generated an independent Bayesian estimate of the parallax and proper motion (referred to as model B). The various models yielded consistent results.

This is a general description of astrometric fitting via χ² minimization. In detail this corresponds to model A1. Differences between the models are described below. Our notation is such that i refers to a star, k refers to an epoch, and j refers to one of the four images at that epoch. There are a total of JK = 32 images and I = 21 field stars.

We register all images by solving for their centers and relative rotations. The plate scale is not necessarily the same for all images. We considered the case where the scale varied independently in the X- and Y-directions, but found no need for this. We limited our fitting to a scalar scale change. The stars have proper motion and a correction for parallax might need to be included.

The measured positions form the data sets x_i,j,k and y_i,j,k. If the positions of all field stars appeared on all frames there would be a total of 672 (x, y) points, but the positions of some stars on some frames are invalid, either because they fall off the frame, or because they exceed the 1 pixel error tolerance we adopted.

The function to be minimized is

$$\begin{equation} \chi ^2=\sum _{i,j\ne 0,k\ne 0}\left[{(X_{i,j,k}-x_{i,j,k})^2 \over \sigma _{x,i,j,0}^2+\sigma _{x,i,j,k}^2} + {(Y_{i,j,k}-y_{i,j,k})^2\over \sigma _{y,i,j,0}^2+\sigma _{y,i,j,k}^2}\right], \end{equation} \tag{ 1 }$$

where the sum is over the valid indices. The model predictions for the positions, relative to a reference position (x_i,0,0, y_i,0,0), are

$$\begin{eqnarray} X_{i,j,k} &=& M_{j,k}[\cos (\theta _{j,k})x_{i,0,0}- \sin (\theta _{j,k})y_{i,0,0}\nonumber\\ &&+\Delta x_{j,k}+\mu _{x,i}t_k+P_{x,k}\Pi _i],\nonumber \\ Y_{i,j,k} &=& M_{j,k}[\cos (\theta _{j,k})y_{i,0,0}+ \sin (\theta _{j,k})x_{i,0,0}\nonumber\\ &&+\Delta y_{j,k}+\mu _{y,i}t_k+P_{y,k}\Pi _i]. \end{eqnarray} \tag{ 2 }$$

The reference position is taken to be the first frame (j = 0) of either the first or fifth epochs (k = 0 or 4) except when the stellar data for that frame is invalid, in which case the next frame of the epoch is used. In all cases, we use the notation 0, 0 to refer to the reference image. M refers to the plate's relative scale factor, or magnification, θ to the plate's relative rotation, μ_x and μ_y to the star's proper motion, t to the difference in epoch between an image and the reference image, and Π = 1/D, where D is the stellar distance. P_x and P_y are the relative parallactic displacements, which depend on the time of year. If the parallax Π is in units of pc⁻¹, then P_x and P_y have units of pixels pc. With these definitions, M_0,0 = 1, Δx_0,0 = Δy_0,0 = θ_0,0 = t₀ = P_x,0 = P_y,0 = 0. It follows that X_i,0,0 = x_i,0,0 and Y_i,0,0 = y_i,0,0. Initially, it is assumed μ_x,i = μ_y,i = Π_i = 0 for all i.

The measurement errors, σ_x,i,k and σ_y,i,k, are the standard deviations among the measured positions for each star at each epoch, which normally numbers J_i,k = 4.

$$\begin{equation} \begin{array}{ll} \sigma _{x,i,k}=\sqrt{\sum _j(X_{i,j,k}- x_{i,j,k})^2/J_{i,k}},\\ \sigma _{y,i,k}=\sqrt{\sum _j(Y_{i,j,k}- y_{i,j,k})^2/J_{i,k}}. \end{array} \end{equation} \tag{ 3 }$$

These errors are same for all values of j for each value of k. As discussed above, we set a floor of 0.01 pixel as the minimum value for σ.

The plate registration is done in the standard way by solving the 4(JK − 1) = 124 simultaneous equations for the minimization of χ²:

$$\begin{equation} \begin{array}{llll} {\partial \chi ^2/\partial \Delta x_{j,k}}=0,\\ {\partial \chi ^2/\partial \Delta y_{j,k}}=0,\\ {\partial \chi ^2/\partial \theta _{j,k}}=0,\\ {\partial \chi ^2/\partial M_{j,k}}=0. \end{array} \end{equation} \tag{ 4 }$$

Note that in these equations, (j, k) ≠ (0, 0).

We use a Newton–Raphson approach to solving the set of equations f_i(x_j) = 0. Assume an initial guess x_0j and make the Taylor expansion to determine the corrections dx_0j:

$$\begin{equation} f_i(x_{0j})+{\partial f_i\over \partial x_j}dx_{0j}=f_{0i}+A_{ij}dx_{0j}=0. \end{equation} \tag{ 5 }$$

Solving this system gives dx_0j = −(A_ij)⁻¹f_0i. The new guess becomes x_1j = x_0j + dx_0j, a new Taylor expansion is made, and a new correction dx_1j is estimated and so forth until convergence is achieved.

Next, the stellar parameters, their proper motions and distances, are determined from the 3I = 63 simultaneous equations

$$\begin{equation} \begin{array}{lll} {\partial \chi ^2/\partial \mu _{x,i}}=0,\\ {\partial \chi ^2/\partial \mu _{y,i}}=0,\\ {\partial \chi ^2/\partial \Pi _{i}}=0. \end{array} \end{equation} \tag{ 6 }$$

Now, using the non-zero values for μ_x,i, μ_y,i, and Π_i, the positional errors and plate registrations are repeated. The scheme converges after about four iterations.

The target star parameters are established once convergence is achieved by minimizing

$$\begin{equation} \chi ^2_{\rm NS}=\sum _{j,k}\left[{(X_{{\rm NS},j,k}-x_{{\rm NS},j,k})^2 \over \sigma _{x,{\rm NS},j,0}^2+\sigma _{x,{\rm NS},j,k}^2} + {(Y_{{\rm NS},j,k}-y_{{\rm NS},j,k})^2 \over \sigma _{y,{\rm NS},j,0}^2+\sigma _{y,{\rm NS},j,k}^2}\right], \end{equation} \tag{ 7 }$$

utilizing the above expressions for $i=\rm NS$ and solving the equations for i = NS, where NS is the index corresponding to the neutron star.

The transformations are based on between 18 and 21 stars in each visit. The median residual error in the transformation (the difference between the modeled stellar positions and the reference positions) is 0.05 pixels; the median measurement error is 0.07 pixels. In principle, the accuracy of the transform depends on the instrumental magnitude of the stars used (Figure 1). Examining the spread of errors, it is clear that the 12 stars with instrumental magnitudes <−9 exhibited less scatter in the residuals. Including only these bright stars gives a median residual error of 0.04 pixels with a median measurement error of 0.02 pixels. However, no systematic uncertainties are introduced by using all the stars in the field, and neither the results nor their uncertainties depend on the choice of reference stars.

We determine the parallactic ellipse using standard formulae (e.g., in the Astronomical Almanac), for the nominal coordinates of the target, α = 18^h56^m3511 and δ = −37°54'305.

In model A2, we use only the nine bright reference stars that appeared in all 32 images to determine the image transformations. This model uses an iterative approach, wherein the plate transformation is fit with all stellar motions fixed at zero. The proper motion is then estimated as a linear trend in the residuals. This is fit, subtracted, and the plate is re-transformed. This process is repeated until the solution converges. The neutron star is not included in the fit. This approach implicitly assumes that all parallaxes are zero. In the second step, the parallax and proper motion are fit to the transformed positions of the neutron star. Finally, this process is repeated using all the stars in the image. No significant parallaxes were found, validating the first assumption.

In model A3, we estimate the individual uncertainties empirically in the _flt images using the IDL xstarfinder procedure (Diolaiti et al. 2000). An empirically determined PSF is fitted to each star in each image, resulting in astrometric uncertainties ranging from ∼ 0.01 pixels for bright objects to ∼0.09 pixel for the faintest objects, including the target. We use these uncertainties along with the distortion-corrected positions from the img2xym_HRC output. The first image in the first visit is used as a reference frame. We fit the 19 stars visible in the first image (including the neutron star) and performed a χ² minimization as described above, fitting the proper motions and parallaxes simultaneously.

Model A4 is based on MPFIT. This model can be run using either some (we selected the brightest 12 stars, those with instrumental magnitudes <−9) or all the stars in the field for determining the transformations. An iterative approach is used, similar to that used in model A2, except that the parallaxes and proper motions are fit simultaneously. Parallaxes can be forced to be zero.

We transformed the positions of all stars in the images to the reference frame (the first image of the first visit) allowing for five free parameters (translations and scale factors in both in X and Y and orientation). We rejected stars which failed to transform (including the target), leaving a subsample of 12 stars to define the reference frame. We estimated the measurement uncertainties by fitting a straight line to the mean errors of these 12 stars as a function of instrumental magnitude. This differs from the scheme used in the A models.

Thus, having positional measurements and corresponding uncertainties of the target in the common reference of frame, we computed posterior probabilities of the parallax and proper motions of the target (all other parameters are treated as nuisance parameters, which are not of immediate interest but which must be accounted for in the analysis of those parameters which are of interest).

The expected geocentric right ascension and declination α_GCRS(t), δ_GCRS(t) of the target at a given epoch t can be expressed in the form

$$\begin{equation} \left\lbrace \begin{array}{@{}ll@{}} \alpha _{\rm GCRS}(t) = \alpha _{\rm BCRS}(t_0) + \Delta \alpha _1(t,t_0,\alpha,\delta,\ mu_{\alpha },\mu _{\delta },\Pi,v_r)\\ \qquad + \Delta \alpha _2(t,\alpha,\delta,\Pi), \\ \delta _{\rm GCRS}(t) = \delta _{\rm BCRS}(t_0) + \Delta \delta _1(t,t_0,\alpha,\delta,\ mu_{\alpha },\mu _{\delta },\Pi,v_r)\\ \qquad+ \Delta \delta _2(t,\alpha,\delta,\Pi), \\ \end{array} \right. \end{equation} \tag{ 8 }$$

where α_BCRS(t₀), δ_BCRS(t₀) are the barycentric coordinates at a reference time t₀, and the pairs (Δα₁, Δδ₁) and (Δα₂, Δδ₂) are the corresponding displacements of the star's position in the right ascension and declination owing to the spatial motion and parallactic motion of the Earth. (μ_α, μ_δ) are the right ascension and declination components of the proper motion, Π is the annual parallax, and v_r is the radial velocity.⁴

In principle, those expected positions and actual measured ones in different epochs may differ by some uncertainty _i owing to the measurements errors plus any real signal in the data that cannot be explained by the model:

$$\begin{equation} \left\lbrace \begin{array}{@{}ll} \alpha _{\rm GCRS}({\rm expected})=\alpha _{\rm GCRS}({\rm measured}) +\epsilon _{\alpha }(t), \\ \delta _{\rm GCRS}({\rm expected})=\delta _{\rm GCRS}({\rm measured}) +\epsilon _{\delta }(t). \end{array} \right. \end{equation} \tag{ 9 }$$

In the absence of detailed knowledge of the effective noise distribution, the most conservative choice of the distribution of positional differences (_i) is a Gaussian and therefore for the likelihood, we may write

$$\begin{eqnarray} &&p(D|M(\theta),I) = \prod _{i=1}^{n_{\rm obs}}\frac{1}{2\pi \epsilon _i}\nonumber\\ &&\quad\times\exp \left[-\sum _{i=1}^{n_{\rm obs }}(d_i-d_i({\rm model})))^2/\big(2\epsilon _i^2\big)\right], \end{eqnarray} \tag{ 10 }$$

where by D(d_i, i = 1, n_obs) we denote either the right ascension or the declination of the target at a given epoch.

Using Bayes theorem, we may estimate posterior probabilities of parameters:

$$\begin{equation} p(M(\theta)|D,I)=\frac{p(D|M(\theta),I)p(M(\theta),I)}{p(D|I)}, \end{equation} \tag{ 11 }$$

where M(θ) is a model with parameters θ(α_BCRS(t₀), δ_BCRS(t₀), Π, μ_α, μ_δ, V_r), and p(M(θ), I) the prior probabilities of the considered parameters. p(D|I) serves as the normalization.

We used the flat priors for the proper motions and the uninteresting parameters. For the parallax, we used both so-called flat priors and the Jeffreys priors. The flat prior is a uniform probability distribution within some reasonable interval. In the Jeffreys prior, formulated as

$$\begin{equation} p(\Pi |I) = 1/(\Pi *\rm Log[\Pi _{\rm hi}/\Pi _{\rm lo}]), \end{equation} \tag{ 12 }$$

the probability density is uniform in the logarithm of the parallax (Π_hi and Π_lo are the upper and lower bounds of the interval considered for the parallax). Use of this prior is motivated by the invariance of conclusions with respect to scale changes in time (Jaynes 1998). This prior is also form invariant with respect to re-parameterization in terms of distance. That is, an investigator working in terms of distance and using a 1/distance prior will reach the same conclusions, as an investigator working in terms of Parallax and using a 1/parallax prior. This would not be true for a prior of another form, for example, a uniform prior. To the extent that parameterization in terms of distance and parallax are both equally "natural," this form of invariance is desirable. Of course, if the prior range of a parameter is small when expressed as a fraction of the center value, then the conclusions will not differ significantly from those arrived at by assuming a uniform prior.

Since the measurement errors play a crucial role in the estimation of the unknown parameters, we considered the following three different cases.

• 1.  
  Positional uncertainties are known for each epoch (from scaling the uncertainties as a function of instrumental magnitude, as described above).
• 2.  
  Positional uncertainties are unknown but are the same at all epochs ($\epsilon _i=\rm const\times \sigma$).
• 3.  
  Positional uncertainties are unknown but are different for different epochs ($\epsilon _i=\rm const\times \sigma _i$).

We checked the plate scale in the reference image by using five stars that are in common with the Two Micron All Sky Survey catalog. The mean plate scale is 28.38 ± 0.53 mas pixel⁻¹, where the uncertainty is the standard deviation of the 10 independent measurements. There are four stars in common with the NOMAD-1 catalog (Zacharias et al. 2005). The mean plate scale is 29.55 ± 1.3 mas pixel⁻¹. These are both consistent with the mean plate scale from Anderson & King (2004).

To test our sensitivity to the small expected parallax, we ran a series of tests wherein we replaced the target with a test star of known parallax and proper motion after transforming the frames. We set the brightness to the mean brightness of the target. We selected a large proper motion similar to that expected from the target. We applied a Gaussian-distributed error to the position of the test star, and then fit the data in an attempt to recover the known motions. Results are shown in Table 2.

The eight visits to the field occurred at four distinct times of the year, within tolerances of 10 days. On a given date, the parallactic factor is the same for all the stars, so we can assume that any difference in position is attributable solely to proper motion. We determined a mean proper motion for each star from the average of these four measurements for each star in the field. Within the uncertainties, the proper motion of the target is fully consistent with those from the global minimizations.

The mean proper motions for the 21 field stars are μ_α = −2.4 ± 4.7, μ_δ = 0.4 ± 3.0 milliarcsec yr^-1. There is no evidence for systemic motions of the field stars.

The full fit solves for four global parameters (zero-point translations, rotation, and plate scale). The largest change in plate scale relative to the reference image is 6×10⁻⁵; the largest change in orientation from that of the reference image is 0.0055° (20 arcsec). Neither of these are significant. We also ran the fitting routines for the cases where the rotation or the scale factor were assumed not to change. This did not affect the results significantly.

We also solved for the case where the parallax of the reference stars is fixed at zero. This also had no significant effect on the results.

Finally, we ran the positions through a publicly available astrometry code, ATPa (G. Anglada Escudé 2010, in preparation; see also Boss et al. 2009). This code was developed for use in the Carnegie Astrometric Planet Search.

We de-rotated the positions reported by img2xym_HRC using the nominal orientation, so that the (x,y) axes correspond to (−R.A., decl.), respectively. We wrote out the coordinates in the correct format to input into the ATPa astrometric solution tool. We used the nominal HRC plate scale. Since the geometric distortions have already been removed, we used calibration rank 1 in ATPa, which corrects for linear offsets, rotations, and magnification changes.

Because this code is still under development, and we ran it as a black box, we use it only as a check that our other codes are producing reasonable results. Results are stated in Table 3 as model T. We do not include the ATPa results in our final parallax determination.