This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Cannonball or Bowling Ball: A Proper Motion and Parallax for PSR J0002+6216

S. Bruzewski Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM 87131, USA F.K. Schinzel An Adjunct Professor at the University of New Mexico. National Radio Astronomy Observatory, P.O. Box O, Socorro, NM 87801, USA G.B. Taylor Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM 87131, USA P. Demorest National Radio Astronomy Observatory, P.O. Box O, Socorro, NM 87801, USA D. A. Frail National Radio Astronomy Observatory, P.O. Box O, Socorro, NM 87801, USA M. Kerr Space Science Division, US Naval Research Laboratory, Washington, DC 20375, USA P. Kumar Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM 87131, USA S. Bruzewski bruzewskis@unm.edu
(Received 2023-09-25; Revised 2023-10-25; Accepted 2023-10-30)
Abstract

We report the results of careful astrometric measurements of the Cannonball pulsar J0002+6216 carried out over three years using the High Sensitivity Array (HSA). We significantly refine the proper motion to μ=35.3±0.6\mu=35.3\pm 0.6 mas yr-1 and place new constraints on the distance, with the overall effect of lowering the velocity and increasing the inferred age to 47.60±0.8047.60\pm 0.80 kyr. Although the pulsar is brought more in line with the standard natal kick distribution, this new velocity has implications for the morphology of the pulsar wind nebula that surrounds it, the density of the interstellar medium through which it travels, and the age of the supernova remnant (CTB 1) from which it originates.

††journal: ApJ††software: Astropy (Astropy Collaboration et al., 2022), CASA (CASA Team et al., 2022), corner.py (Foreman-Mackey, 2016), Difmap (Shepherd, 1997), emcee (Foreman-Mackey et al., 2013) matplotlib (Hunter, 2007), Numpy (Harris et al., 2020), Scipy (Virtanen et al., 2020)

1 Introduction

Identifying compact objects associated with supernova remnants provides a unique window into the diverse outcomes of core-collapse supernovae. A veritable zoo of such young neutron stars have been found (Popov, 2023) including central compact objects (CCOs) and various types of magnetar, such as soft gamma-ray repeaters (SGRs) and anomalous X-ray pulsars (AXPs) (Gaensler et al., 2001). Young radio pulsars comprise the majority of neutron stars found in supernova remnant (SNR) associations (Green, 2019; Ferrand & Safi-Harb, 2012). The study of these systems has helped constrain initial pulsar periods, magnetic fields, and beaming fractions, in addition to their kick velocity at birth (Frail & Moffett, 1993; Hansen & Phinney, 1997; Igoshev et al., 2022; Kapil et al., 2023; Johnston et al., 2005).

Schinzel et al. (2019, henceforth Paper I) noted that the 115 ms γ\gamma-ray and radio pulsar PSR J0002+6216 was found at the head of a bow-shock pulsar wind nebula (PWN). Follow-up radio and X-ray observations resolved the PWN and showed that it was consistent with a high Mach shock formed from the wind from the energetic high velocity PSR J0002+6216 coming into ram pressure balance with its surrounding medium (Kumar et al., 2023). Remarkably, the tail of the PWN extends at least 7-10 arcminutes from PSR J0002+6216, pointing backward toward the geometric center of the Galactic SNR CTB-1 (G116.9+0.2) 28 arcminutes away. The authors argued that CTB 1 is the remnant of the supernova that produced PSR J0002+6216, indicating that PSR J0002+6216 was born with an unusually high velocity (VPSR>1000V_{\rm PSR}>1000 km s-1) that allowed it to escape the parent SNR.

Not all of the young pulsars found in or near SNRs are real associations. Often the case is made on the basis of positional coincidence and some agreement between ages and distances. The gold standard for making such claims is accurate pulsar proper motion and parallax measurements. Often both the magnitude and direction of proper motion can reinforce the claimed association by showing that the pulsar likely originated near the geometric center of the SNR (Deller et al., 2012; Van Etten et al., 2012; Shternin et al., 2019; Johnston & Lower, 2021; Long et al., 2022). In other cases, proper motion measurements either exclude an association or require more complex scenarios to remain tenable (Brisken et al., 2006; de Vries et al., 2021; Espinoza et al., 2022). Pulsars (or PWN) with large offsets from the geometric center of the SNR, including pulsars outside the SNR (e.g, Ng et al., 2012; Motta et al., 2023), are especially important for constraining the high birth velocity tail of PSRs (Frail et al., 1994; Hobbs et al., 2005; Verbunt et al., 2017). The burden of proof is high, since in these cases the probability of a chance association is greatly increased (Gaensler & Johnston, 1995). Proper motion measurements of such associations have given mixed results at best (Zeiger et al., 2008; Hales et al., 2009). In the case of PSR J0002+6216/CTB-1, Paper I bolstered their case by using nearly a decade of data from the Fermi Gamma-Ray Space Telescope to infer a proper motion of μ=115±33\mu=115\pm 33 mas yr-1 and θμ=121∘±13\theta_{\mu}=121^{\circ}\pm 13. This value, although of modest significance, agrees in magnitude and direction with the PSR-SNR angular offset, assuming an age for the system of 10 kyrs.

In this work, we have undertaken VLBI observations over the course of several years towards providing an accurate measure for the proper motion and parallax of PSR J0002+6216. Having such values allows us to provide a more definitive association between PSR J0002+6216 and CTB-1, as well as to probe the implications of the SNR age implied from the more accurately measured velocity of the pulsar. Section 2 will discuss our observational setup, data processing, and the fitting of parallax and proper motion to our data. Section 3 describes the best-fit parameters, and Section 4 discusses their implications. Finally Section 5 concludes and summarizes this work.

2 Data

2.1 Observations

Observations were performed under project code BS278 (VLBA/19B-048) with the High Sensitivity Array (HSA111https://science.nrao.edu/facilities/vlba/HSA), comprised of NRAO’s Very Long Baseline Array (VLBA) and the phased Karl G. Jansky Very Large Array (VLA), as well as the Effelsberg 100 m radio telescope. The observations spanned a period of just over two years (2020 February 10 – 2022 July 30). The antennas unavailable during a specific observation, together with additional summarizing characteristics, are noted in Table 1.

Initially, we performed a test observation (BS278Z) to check for suitable phase reference calibrators, both for VLBA-only in-beam calibration (selected from the VLA L-band image presented in Paper I) and nearby bright calibrator sources selected from the VLBA calibrator catalog222https://obs.vlba.nrao.edu/cst. This test did not yield a suitable in-beam calibrator and only identified one suitable phase calibrator, J0003+6307, separated by 0.87∘ with a flux density of 46.5±\pm1.3 mJy and a peak of 10.48±\pm0.24 mJy beam-1. Given the weak nature of this calibrator, we included a secondary phase calibrator J2339+6010, which is at a distance of 3.53∘ from our target with a flux density of about 70 mJy and a peak of about 36 mJy beam-1. In addition, we observed J0014+6117 at a distance of 1.71∘ (about 30 mJy / 15 mJy beam-1 peak) as a VLA phasing calibrator. J0014+6177 can also be used as cross-check calibrator. However, it seems to be resolved out on the longest baselines. J0319+4130 and J0137+3309 served as fringe finders, bandpass calibrators, and flux references. The calibrators and phase centers are summarized in Table 2.1.

Table 1: HSA Observation Summary
Date of Obs. S A Missing Issues
2019-07-20aaNon-detections Z 11 EF,Y –
2019-10-09aaNon-detections A 12 – NL,EF gaps
2020-02-10aaNon-detections B 11 MK EB,Y gaps
2020-05-04 C 11 SC KP,FD,OV gaps
2020-08-08 D 12 – SC,KP,NL,LA gaps
2020-10-05 E 11 HN KP,SC gaps
2020-12-14 F 12 – OV high winds
2021-03-28aaNon-detections G 11 OV SC,PT,MK 1h lost
2021-06-18 H 11 SC Y late, FD timing
2021-09-20 I 12 – –
2022-02-14 J 11 HN MK missing data
2022-07-30 K 10 KP,MK –

Note. — The S column denotes the segment of BS278, and the A column counts the number of participating antennas. Issues: Gaps are noted where a significant amount of data are missing from correlation.

The Z and A segments of BS278 were observed using the 18 cm receiver and a frequency range of 1.568–1.824 GHz and the RDBE personality of the VLBA backend with a recording data rate of 2 Gbps. This resulted in non-detection of our target, primarily due to strong radio interference at some of the antenna sites. Subsequent observations (BS278C onward) thus used a lower frequency range centered on a mostly interference-free part of the spectrum, 1.268–1.524 GHz in the 21 cm wavelength band. For Effelsberg we utilize the center horn of the 7-beam 21 cm receiver. In addition, the phased-VLA was used to record the full allowable bandwidth of 1 GHz (starting with segment C) to determine the pulsar ephemerides using VLA’s WIDAR correlator for later binned correlation of the HSA observations. The phased-VLA observations yielded detections in all cases with an average flux density of ∼30\sim 30 uJy. The pulsar is linearly polarized and shows a rotation measure of −178.5±2.5-178.5\pm 2.5 rad m-2.

The HSA correlation was performed using the VLBA DiFX correlator (Deller et al., 2011) using 256 correlation channels per one of the eight dual-polarization spectral windows with 32 MHz bandwidth each and a 0.5 s dump time. This allows full coverage of 256 MHz in bandwidth centered on 1.396 GHz. As mentioned before, for each phased VLA observation, pulsar timing data was recorded to obtain pulsar ephemerides of our target pulsar, which were then supplied to DiFX for binned correlation, where 10 bins in pulsar spin phase were formed to increase the signal-to-noise ratio of the pulsar emission, with the first bin (“ON”) centered on the main pulse of the pulsar. The VLBA DiFX correlator is set such that the flux density remains the same as the unbinned equivalent. In addition, multiple phase centers were placed on locations of known bright sources that were identified from the wide-field VLA L-band observation mentioned above. In addition, we updated the target phase center four times throughout the campaign to account for the pulsar proper motion and keep the pulsar close to the correlation phase center. Correlation phase centers are provided in Table 2.1.

Table 2: Calibrators and Phase Centers
Name R.A. Dec. Dist. Description
h:m:s d:m:s deg.
  J0001+6051 00:01:07.099891 +60:51:22.79829 1.43 Position Check Calibrator
J0003+6307 00:03:36.511479 +63:07:55.87126 0.87 Primary Phase Calibrator
J0014+6117 00:14:48.792109 +61:17:43.54262 1.71 VLA phasing calibrator
J0137+3309 01:37:41.299543 +33:09:35.13377 3.99 Position Check Calibrator
J0319+4130 03:19:48.160114 +41:30:42.10568 – Fringe Check / Bandpass Calibrator
J2236+2828 22:36:22.470849 +28:28:57.41328 – Flux Check / Fringe Check
J2339+6010 23:39:21.125227 +60:10:11.84951 3.53 Secondary Phase Calibrator
J0003+6219 00:02:53.527 +62:19:16.940 – Correlation phase center
J0002+6205 00:01:44.120 +62:04:30.690 – Correlation phase center
J0001+6228 00:00:43.195377 +62:27:59.722930 – Correlation phase center
J0000+6222 00:00:17.959930 +62:22:20.034252 – Correlation phase center
J0000+6215 00:00:24.397047 +62:15:21.717124 – Correlation phase center
J0000+622A 00:00:17.753 +62:22:14.354 – Correlation phase center
TARGET 00:02:58.17 +62:16:09.4 – Correlation phase center for segment B
TARGET 00:02:58.205 +62:16:09.510 – Correlation phase center for segments D–E
TARGET 00:02:58.207 +62:16:09.500 – Correlation phase center for segments F–H
TARGET 00:02:58.210 +62:16:09.4958 – Correlation phase center for segments I–K

Note. — List of Calibrators and correlation phase centers for BS278C-K, including updated target positions. Coordinates are in J2000.

2.2 Calibration & Model-Fitting

The correlated data were calibrated using the Astronomical Image Processing System (AIPS; Greisen, 2003) using the included VLBAUTIL procedures. Standard astrometric calibration was performed applying corrections for Earth’s orientation, the ionosphere, for digital sampling, and delay. After this, bandpass, sampler, gain, and parallactic angle corrections were applied, followed by determination and application of phase corrections using the calibrator, J0003+6307. For applying calibration solutions to the different phase centers and bins the VLBAMPHC procedure was used. After this, the calibration corrections were applied to the correlated data and split by observed source and written to separate UVFITS formatted files for further analysis.

Then we used the Difmap package (Shepherd, 1997) for additional manual flagging of radio frequency interference, primarily from Global Navigation Satellite System interference, such as GPS, Galileo, or GLONASS in the 21 cm wavelength band. We then used uu,vv-plane model-fitting of a delta function, as a good representation for an unresolved point source, to determine the position of our target and calibrators for each epoch. The detections were typically at the 5σ\sigma level, which did not allow for robust fitting of a Gaussian function to the visibilities. We also used Difmap to produce restored clean images for the calibrators and target.

One issue arising from this approach is that, while Difmap’s model fitting excels in this particular case of faint sources, it does not provide much in the way of uncertainties for the calculated best-fit coordinates. Ideally we would like to have an estimate at each epoch for how well the position of a source of a certain flux density can be localized amid some level of background noise. We do this by making a copy of our original visibility data, then using CASA (CASA Team et al., 2022) to replace all data-points with an artificial point source of known position and flux density (set to approximately match the properties of PSR J0002+6216). Random simple noise is then added to the data, which we re-scale using the existing weights. Difmap is then provided this artificial data, as well as a perturbation of the position to use as an initial starting model, and made to perform model fitting for 50 steps, after which the final model is recorded to disk.

Conveniently, none of the algorithms involved in this simulation are all that intensive, and so the entire pipeline can be run using a single script on a single core. Given that the only difference between two simulations will be the seed of the random noise, the problem is embarrassingly parallelizable, and so we perform 1000 such simulations at each epoch. Once the final models of all simulations have been collected and tabulated, we simply need to characterize the scatter of the models to establish our fitting uncertainty. We note that the scatter in Right Ascension and Declination seem to be highly correlated in all epochs, likely due to our uu,vv-coverage, such that we cannot disregard the off-diagonal terms of the covariance matrix at this stage. Instead we opt to record the covariance matrix at each epoch, measured over the coordinates Δ​α∗=(α−α¯)⋅cos⁡δ¯\Delta\alpha^{*}=(\alpha-\bar{\alpha})\cdot\cos\bar{\delta} and Δ​δ=(δ−δ¯)\Delta\delta=(\delta-\bar{\delta})333Note that this cosine correction will appear fairly regularly throughout this work.

2.3 Atmospheric Effects

For our observations, the positional uncertainties on a given date will be dominated by the effects of the ionosphere. This can be corrected for at some level by ionospheric corrections performed during calibration, but for corrections at the milliarcsecond level, typically further refinements are required. Similar work (Deller et al., 2016, 2019; Ding et al., 2023) has typically made use of one or more relatively bright in-beam calibration sources, which (if truly extra-galactic and stationary) can be simultaneously phase-calibrated, effectively removing atmospheric effects by providing an absolute position offset between the two sources at the sub-milliarcsecond level.

Unfortunately, such calibration is not possible in the case of J0002+6216. The area immediately surrounding the pulsar is devoid of any other sources out to the edges of the VLBA field of view. Furthermore, as this source is quite faint (∼\sim22–30 uJy in single dish/phased VLA observations), we require the additional sensitivity from the inclusion of Effelsberg and the phased VLA, which have an even smaller field of view, further ruling out this option. The typical detection in the “ON” bin of our observations, resulted in a flux density of ∼\sim14 uJy, suggesting that a significant fraction is lost due to residual ionospheric phase errors that dominate in this frequency range.

Lacking a suitable source in the field of view, we opt to leverage a nearby but out-of-beam source. This secondary source, dubbed SRC2, has a peak flux density near 700 uJy beam-1, meaning its own positional uncertainties (as measured via simulations described in Section 2.2) are more than an order of magnitude better constrained than those of our target. Its position seems to experience ionospheric variations on the same scale as other calibrators, and has no apparent proper motion across our observations. As SRC2 is fairly close to our target, we can assume that it lies behind roughly the same region of the ionosphere (as described approximately by Taylor et al., 1999), and should have its position perturbed similarly. With this in mind we use the displacement of SRC2 from its median position to apply a correction to the position of J0002+6216 at each epoch. This should very roughly approximate the results of the aforementioned simultaneous phase-calibration, albeit with higher uncertainties. These updated positions are reported in Table 3.

Table 3: Target Positions
Date Position (J2000) σα∗\sigma_{\alpha}^{*} σδ\sigma_{\delta}
MJD mas\mathrm{mas} mas\mathrm{mas}
58973.44 00h02m58.20495s +62d16m09.50991s 1.29 1.11
59069.70 00h02m58.20652s +62d16m09.50328s 1.18 1.07
59127.89 00h02m58.20686s +62d16m09.50229s 1.19 1.09
59197.80 00h02m58.20773s +62d16m09.49950s 1.21 1.22
59383.63 00h02m58.21090s +62d16m09.49124s 1.28 1.12
59477.80 00h02m58.21150s +62d16m09.49155s 1.17 1.10
59624.17 00h02m58.21284s +62d16m09.48623s 1.13 1.05
59790.60 00h02m58.21578s +62d16m09.47532s 1.23 1.12

Note. — Position of J0002+6216 at each epoch it was detected. Uncertainties are derived from simulations of position fitting to our UV-data, as well as systematic uncertainty from the ionosphere. Here we report only the diagonal terms of the covariance matrix as our uncertainties, as the correlation is fairly small.

In order to fully treat the uncertainty, we need a probe of the systematic uncertainty produced by the ionosphere. As a proxy for this, we choose to measure the scatter of each of our calibrators relative to each of their mean positions across all epochs. We then measure the covariance of this distribution and add it to the simulated positional covariance at each epoch. The scale of this systematic uncertainty is approximately 1 mas (found from the diagonal elements in the covariance matrix), and is only very mildly correlated in R.A. and Declination. The uncertainties reported in Table 3 include this systematic atmospheric term.

2.4 Parallax and Proper Motion Fitting

For fitting parallax and proper motion we adopt the model used by Gaia as described in Lindegren et al. (2016) Equation 2, similarly assuming that the effect of radial proper motion is small enough as to be undetectable in our observations. Using that equation, a position can be found as a function of time, given the set of variables {α0\alpha_{0}, δ0\delta_{0}, ω\omega, μα∗\mu_{\alpha}^{*}, μδ\mu_{\delta}} and choosing a reference epoch. For our analysis we choose a reference epoch of MJD 58849, which corresponds to January 1st 2020, slightly prior to our first observation.

We implement a Gaussian prior on the total proper motion of 115±33115\pm 33 mas yr-1, as reported in Paper I, though we note that the final parameter estimates do not seem to be strongly effected by this prior. Because our position uncertainties are covariant, the equation for log-likelihood takes on the form

ln⁡ℒ=−12​∑i=0n[2​ln⁡2​π+ln⁡(detVi)+xiT​Vi−1​xi],\ln{\mathcal{L}}=-\frac{1}{2}\sum_{i=0}^{n}\left[2\ln{2\pi}+\ln{(\det{V_{i}})}+x_{i}^{T}\ V_{i}^{-1}x_{i}\right]\mathrm{,} (1)

where at a given epoch ii, xi=(mi−oi)x_{i}=(m_{i}-o_{i}) represents the difference between the model position mi=(αm,δm)im_{i}=(\alpha_{m},\delta_{m})_{i} and observed position oi=(αo,δo)io_{i}=(\alpha_{o},\delta_{o})_{i}, and ViV_{i} represents the observed/estimated covariance. We then pass the data, uncertainties, prior, and log-likelihood functions to the emcee Python package. We use 200 walkers initialized in a small space around initial approximate estimates provided by least-squares fitting of a straight line, and the scatter relative to that. The walkers progress for 20,000 steps, after which we calculate the maximum auto-correlation time τ\tau (typically of order 50 steps) and use that to discard the first 3​τ3\tau steps (prior to burn-in) and thin by a factor of τ/2\tau/2. These sample chains are flattened, and then can be evaluated to find best fit values and uncertainties for each of our parameters as shown in Table 4.

Table 4: Estimated and Calculated Values
Parameter Estimate
α0\alpha_{0} 00h02m58.20346s ± 0.89\pm\ 0.89 mas
δ0\delta_{0} +62d16m09.51294s ± 0.81\pm\ 0.81 mas
ω\omega 0.63±0.450.63\pm 0.45 mas
μα∗\mu_{\alpha}^{*} 32.52±0.5932.52\pm 0.59 mas yr-1
μδ\mu_{\delta} −13.71±0.53-13.71\pm 0.53 mas yr-1
μtot\mu_{\text{tot}} 35.30±0.6035.30\pm 0.60 mas yr-1
θPA\theta_{\text{PA}} 112.86±0.83112.86\pm 0.83 deg
Age 47.60±0.8047.60\pm 0.80 kyr
v2kpcv_{\text{2kpc}}aaThese velocities assume distances with no uncertainty. 334.90±5.66334.90\pm 5.66 km s-1
v3kpcv_{\text{3kpc}}aaThese velocities assume distances with no uncertainty. 502.35±8.4502.35\pm 8.4 km s-1

Note. — Values provided are the median of all samples for each parameter. Uncertainties are reported at the 1-σ\sigma level.

3 Results

Refer to caption
Figure 1: The datapoints along with the median fit produced by our program. Here, color corresponds to time, such that purple (the beginning of the viridis scale) represents the beginning of the year. The closer a data point is to a similar color in the fit, the better the fit. Ellipses surrounding each data point represent the 1, 2, and 3σ\sigma positional uncertainties.

Our best-fit model is illustrated in Figure 1. Here, we mainly focus on the estimates of parallax and proper motion, as those have the largest physical implications. From the best-fit values for proper motion along each axis, we can calculate the total proper motion magnitude μtot=μα∗2+μδ2=35.30±0.60\mu_{\mathrm{tot}}=\sqrt{\mu_{\alpha}^{*2}+\mu_{\delta}^{2}}=35.30\pm 0.60 mas yr-1 and position angle θPA=arctan⁡(μα∗/μδ)=112.86∘±0.83∘\theta_{\mathrm{PA}}=\arctan{\left({\mu_{\alpha}^{*}}/{\mu_{\delta}}\right)}=112.86^{\circ}\pm 0.83^{\circ}. As the total proper motion is generally of more interest than its individual components, we illustrate the distribution in this parameter against parallax in Figure 2. As the full five-parameter corner plot is quite large and does not show any correlations of note, we opt to make it available upon request.

Refer to caption
Figure 2: A 2D histogram of the relevant parameter distributions. Here we show the total proper motion, calculated using its individual components.

If we compare our total proper motion and position angle with the initial estimates of these parameters found from Fermi data in Paper I (μtot=115±33\mu_{\mathrm{tot}}=115\pm 33 mas yr-1 and θPA=121∘±13∘\theta_{\mathrm{PA}}=121^{\circ}\pm 13^{\circ}), we find differences of 2.37​σ2.37\sigma and 0.59​σ0.59\sigma, respectively. The agreement in position angle is expected, as the prior estimate of that value also matched well with the direction of the orientation of the pulsar wind nebula as described in Paper I and later in Kumar et al. (2023), with the later finding a position angle of θPA=111.13∘±0.52∘\theta_{\mathrm{PA}}=111.13^{\circ}\pm 0.52^{\circ}.

Next we examine the parallax of the source, for which a value of ω=0.63\omega=0.63 can be naively inverted to find a distance and velocity of 1.59 kpc and 244 km s-1 respectively. That said, such a naive inversion will suffer from the Lutz-Kelker bias, and can produce meaningless or aphysical values in cases where the parallax is not incredibly well constrained (Verbiest et al., 2010). Ultimately a proper measurement of the distance to J0002+6216 via parallax is thwarted by the effects of the ionosphere. The scale of its perturbation on our observations is above the scale we expect to see the parallax at, and in fact Paper I showed that we can largely rule out parallaxes larger than 1 mas or so.

Perhaps the best estimate for distance currently possible comes from the association with the remnant CTB-1. The new estimates of proper motion trace back a path directly through the geometric center of the SNR (as seen in Figure 3 and discussed more in Section 4), making a chance association quite unlikely. This geometric argument has been sufficient for other systems such as the recently discovered ”Mini Mouse” (Motta et al., 2023). If we assume this association to be true, then our best constraints for the distance to the pulsar are once again approximately 1.5-4 kpc (Landecker et al., 1982; Hailey & Craig, 1994; Fesen et al., 1997; Yar-Uyaniker et al., 2004), plus an upper limit d≲7d\lesssim 7 kpc derived from X-ray and γ\gamma-ray efficiencies (Zyuzin et al., 2018; Wu et al., 2018).

4 Discussion

Refer to caption
Figure 3: The possible history of the pulsar, accounting for uncertainties in position and proper motion. While PSR J0002+6216 is likely to have originated at or near the geometric center of CTB 1 (G116.9+0.2 in this Figure), here we illustrate the possible positions if we assume that it originated somewhere beyond the SNR, out to approximately the characteristic age of the pulsar. The black circles represent that possible 1​σ1\sigma, 3​σ3\sigma, and 5​σ5\sigma positional scatter at each epoch, whereas the shaded regions represent the paths between those epoch contours. Colored circles represent remnants found in Green (2019).

Given our significantly improved measurement of the proper motion, we can begin to discuss the implications of the pulsar’s space-velocity. This velocity can now be written as

v=334.90​(d2​kpc)​km​s−1​,v=334.90\left(\frac{d}{2\ \mathrm{kpc}}\right)\mathrm{km\ s}^{-1}\text{,} (2)

where we have chosen to use 2 kpc as our reference point to match the assumptions in Paper I. Our new estimate has the effect of decreasing the speed by a factor of more than three from the original estimate of over 1000 km s-1. For the range 1.5–4 kpc we find velocities 251–670 km s-1, and more specifically for the value of 3.1±0.43.1\pm 0.4 kpc obtained by Hailey & Craig (1994) we find 519.09±67.56519.09\pm 67.56 km s-1. While any of these values would still place PSR J0002+6216 among a rare category of high-velocity pulsars compared to the overall distribution (Hansen & Phinney, 1997; Verbunt et al., 2017), it could perhaps now be considered for demotion from the ranks of other Cannonball pulsars (see Hobbs et al., 2005; Shternin et al., 2019) to something more akin to a vigorously thrown bowling ball.

Next, we examine the implications of our proper motion result for PSR J0002+6216, using three different measures of the pulsar age. This approach does not require estimates of other (uncertain) parameters such as the distance, and thus should yield secure results. In the first instance, we assume that the age of PSR J0002+6216 is equal to its characteristic age (τc\tau_{c}=306 kyrs). In this case, our measured proper motion places the pulsar birthplace more than 2∘2^{\circ} away, well outside the boundaries of CTB 1. In this instance, there is no physical PSR-SNR association, just a line-of-sight coincidence. There are at least two weaknesses with this hypothesis. Ng et al. (2012) have shown that long PWN tails that point back toward SNRs have a very low posterior probability of occurrence. If anything, the chance probability is even smaller for this system since both our proper motion result and the PWN tail of PSR J0002+6216 can be accurately traced back to the geometric center of CTB 1.

Another objection could be making the (common) assumption that τc\tau_{c} can be used as a proxy for the actual age of the system. This has been shown to have questionable accuracy, at least for young pulsars in which τc\tau_{c} is seen to be substantially lower than or greater than the quoted age of the associated SNR (Blazek et al., 2006; Popov & Turolla, 2012; Johnston & Lower, 2021). This discrepancy with characteristic age implies a relatively long birth spin period, close to the present period of 115 ms, and has marked implications for kick models which predict a close association between the direction of proper motion and the rotation axis (e.g. Ng & Romani, 2007; Coleman & Burrows, 2022). Examining this alignment would require more sensitive and full Stokes measurement of the pulsar profile to improve existing measurements of both the rotation and magnetic axes relative to the line of sight (Wu et al., 2018). Figure 3 illustrates the path of the pulsar over time, out to the distance that would be covered under this assumption of age.

Motivated by these findings, we consider the implications if the pulsar age is equal to the age of the SNR of 10±\pm2 kyrs, derived in Paper I from several sources in the literature. Using the measured proper motion, we find that the pulsar will have traveled only 6′ from current location, placing it outside the southeastern edge of CTB 1 and thus is not associated. However, the difficulty in adopting this age is that we would need to explain how PSR J0002+6216 was born inside the tail of its own PWN. In our original discovery image (Paper I) the tail is at least 7′ in length, and in the lower resolution Canadian Galactic Plane Survey images (CGPS, Kothes et al., 2006) it can be traced 11′ back to the edge of CTB 1. It seems more likely that the true age of PSR J0002+6216 is greater than 10 kyr.

Finally, we can examine the kinematic age of the system. If we assume that the pulsar originated from the geometric center of the remnant, then we can use the 28±′1′{}^{\prime}\pm{1}^{\prime} offset between the pulsar and the geometric center to find a kinematic age for a system of approximately 47.6±0.847.6\pm 0.8 kyr. This has the effect of implying CTB 1 is significantly older than the estimates discussed above. However, this is not completely unexpected, as there is substantial evidence that CTB 1 is in its radiative phase, and more specifically belongs to a class of mixed morphology remnants, as evidenced by its ring-like radio emission and centrally located X-ray emission (Lazendic & Slane, 2006). This means that typical Sedov-Taylor scaling relationships (Sedov, 1959; Taylor, 1950) one might use for younger remnants may not produce accurate results from X-ray properties and could vary significantly depending on various assumptions.

To illustrate the effect of this increased age on the parameters of the system, we can model the approximate evolution of CTB 1 over time, following the procedure roughly described in Vink (2012). The remnant initially follows the Sedov-Taylor equations up until radiative losses become important, and further evolution is directed by momentum conservation. Adopting the equation for the latter phase from Toledo-Roy et al. (2009), our model can be written as

R​(t)≈{5.03​(E51/n0)1/5​t2/5​pc,t<tr​a​drr​a​d​(1.58​ttr​a​d−0.58)1/4,t≥tr​a​dR(t)\approx\begin{cases}5.03\ ({E_{51}}/{n_{0}})^{1/5}\ t^{2/5}\ \mathrm{pc},&t<t_{rad}\\ r_{rad}\left(1.58\ \frac{t}{t_{rad}}-0.58\right)^{1/4},&t\geq t_{rad}\end{cases} (3a)
tr​a​d≈44.6​(E51/n0)1/3​kyrt_{rad}\approx 44.6\ ({E_{51}}/{n_{0}})^{1/3}\ \mathrm{kyr} (3b)
rr​a​d≈23​(E51/n0)1/3​pcr_{rad}\approx 23\ ({E_{51}}/{n_{0}})^{1/3}\ \mathrm{pc} (3c)

where tr​a​dt_{rad} and rr​a​dr_{rad} are the time and radius at which the remnant transitions into the radiative phase. Note that this evolution is dependent on the term (E51/n0)({E_{51}}/{n_{0}}), the ratio of the supernova energy (in units of 105110^{51} ergs) and the pre-shock hydrogen density (in units of cm-3). This means that for a given choice of energy, density and distance, one can accurately describe the angular radius as a function of time, allowing constraints on those parameters knowing the current angular size of CTB 1 (∼17.5′\sim 17.5^{\prime}) and its presumed age (∼48\sim 48 kyr). From this we can see that from reasonable choices of distance, (E51/n0)≈({E_{51}}/{n_{0}})\approx 0.01–0.5, which implies a high-density medium, a weak supernova, or some combination of the two.

Furthermore, for mixed morphology remnants Cox et al. (1999) derived that the time/radius at which the radio shell will begin to form can be written as

ts​h​e​l​l≈53​E513/14​n0−4/7​kyrt_{shell}\approx 53\ E_{51}^{3/14}n_{0}^{-4/7}\ \mathrm{kyr} (4)
rs​h​e​l​l≈24.9​E512/7​n0−3/7​pcr_{shell}\approx 24.9\ E_{51}^{2/7}n_{0}^{-3/7}\ \mathrm{pc} (5)

allowing one to place constraints similarly using the inner edge of the radio emission ring (which has an angular size of ∼13′\sim 13^{\prime}). Figure 4 shows such constraints for the inner and outer radii, assuming a distance of 3 kpc so as to probe the supernova energy and local pre-shock density. Most choices of distance seem to require fairly typical energies, and the high value for density could be expected given that most mixed-morphology remnants seem to form in regions of denser than typical ISM. At the 3 kpc distance, the best match to the constraints comes when E51≈0.7E_{51}\approx 0.7 and n0≈4n_{0}\approx 4 cm-3.

Refer to caption
Figure 4: Constraints in the density/energy space derived from the observed outer and inner radii of CTB 1, assuming a distance of 3 kpc (near that found in Hailey & Craig, 1994). The red region represents outer radii which are within 1′1^{\prime} of 17.5′17.5^{\prime}, while the blue region represents inner/shell radii which are within 2′2^{\prime} of 13′13^{\prime}. The space where both regions overlap is shown for 3 kpc as well as a few other distances.

Interestingly, this result of a higher than expected density also has implications for the bow shock of PSR J0002+6216. As described in Kumar et al. (2023), the pulsar is moving rapidly through the local ISM, generating a bow shock around it and trailing a long and highly collimated tail of emission for several arcminutes. Paper I estimated a Mach number ℳ≈200\mathcal{M}\approx 200 from initial measurements and comparisons with systems with similar morphology (Ng et al., 2012; Kargaltsev & Pavlov, 2008). Assuming the distance and density used as an example above (and a 10% helium mixture as in Paper I), we estimate that the pulsar would have a mach number ℳ≈110\mathcal{M}\approx 110, still broadly consistent with the observed collimation of the tail. Without an increase in density, the decrease in velocity due to the updated proper motion would lower the overall Mach number, resulting in a less striking collimation. Change in various parameters also has the effect of lowering the standoff distance of the shock, derived by equating the wind pressure and the ram pressure of the pulsar (Frail et al., 1996) to

rS​O=0.0015​(n03.5​cm−3)−1/2​(d3​kpc)−1​pc​,r_{SO}=0.0015\ \left(\frac{n_{0}}{3.5\ \mathrm{cm}^{-3}}\right)^{-1/2}\left(\frac{d}{3\ \mathrm{kpc}}\right)^{-1}\mathrm{pc}\text{,} (6)

which would be a few times smaller than the best estimates from Kumar et al. (2023), and among the smallest standoffs known.

5 Conclusions

Using several years worth of HSA data, we have significantly refined the estimate of the proper motion of the “Cannonball” pulsar J0002+6216, to a new value of μtot=35.30±0.60\mu_{\mathrm{tot}}=35.30\pm 0.60 mas yr-1. This new estimate is more than three times smaller than the previous best estimates made in Schinzel et al. (2019) using Fermi data, and as such would act to lower the spatial velocity perpendicular to the line of sight. This new proper motion also provides for a precise distance-independent kinematic age measurement for the system of 47.60±0.8047.60\pm 0.80 kyr, significantly older than previously estimated from the pulsar/remnant, but also still significantly younger than the pulsar’s characteristic age.

We also attempt to provide a more accurate measure of the distance to the source, as that will also play a role in determining the kinematic properties. While a value for the parallax of ω=0.63±0.45\omega=0.63\pm 0.45 mas was obtained, the high uncertainties make a definitive distance difficult to determine. These uncertainties are dominated by scattering in the ionosphere, but we also see uncertainty simply from trying to fit a position to a fairly weak source. Refinement of this value would require both increased sensitivity and more complex observing strategies, which could be used to isolate and remove these ionospheric effects. We note that this is likely to be an excellent target for the ngVLA in the coming decade, as it will have baselines similar to the present-day VLBA, but with significantly more antennas allowing for a dramatic increase in sensitivity.

As an example, in one hour the HSA444http://old.evlbi.org/cgi-bin/EVNcalc.pl can be expected to reach a thermal noise level of approximately 13 uJy at 1.4 GHz, which then necessitates longer and more difficult to schedule observations to provide enough signal to noise for a 14 uJy source such as PSR J0002+6216. In comparison, a one hour observation on the ngVLA555https://ngvla.nrao.edu/page/performance at 2.4 GHz should reach a thermal noise of 0.24 uJy. Even if we assume that the pulsar will dim at this higher frequency by a factor of about 2, following the standard pulsar spectral index of α≈−1.4\alpha\approx-1.4 (Bates et al., 2013), we should still expect an approximate signal-to-noise ratio of 27 after just one hour. Combine this with the possibility of more advanced observation and calibration techniques (such as simultaneous subarray observations of a calibrator source), and it is likely that ngVLA could provide an estimate of the parallax at significantly reduced uncertainty.

Beyond the parallax, further questions also remain for this system. The proper motion, now reduced in magnitude but with substantially less uncertainty, combined with the expansion of CTB 1  seem to be indicative of a fairly large ISM density. The implied value would be beyond the low density (∼\sim0.1–1.0 cm-3) which has typically been inferred from the remnant’s relatively weak Hα\alpha emission (Landecker et al., 1982) or from the non-detection of tracers such as CO (Zhou et al., 2023). While this can be somewhat mitigated if one assumes a weaker supernova energy (E51<1E_{51}<1), the collimation of the pulsar’s bow shock is decoupled from this term and still requires a higher density given this lower velocity. Further work will be necessary to properly constrain the excess density and identify its origin.

6 Acknowledgements

SB, FKS, and GBT acknowledge support from NASA Fermi Guest Investigator Program, grants 80NSSC18K1725, 80NSSC19K1508, 80NSSC22K1643. Work at NRL is supported by NASA. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Support for this work was provided by the NSF through the Grote Reber Fellowship Program administered by Associated Universities, Inc. National Radio Astronomy Observatory and the Student Observing Support Program. We thank Tom Maccarone, Adam Deller, and David Green for insightful discussions toward the goals of this work, as well as Jason Wu for his supporting of early HSA observations from Effelsberg.

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.

This research has made use of NASA’s Astrophysics Data System and has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has made use of data, software and/or web tools obtained from NASA’s High Energy Astrophysics Science Archive Research Center (HEASARC), a service of Goddard Space Flight Center and the Smithsonian Astrophysical Observatory, of the SIMBAD database, operated at CDS, Strasbourg, France.

References

  • Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  • Bates et al. (2013) Bates, S. D., Lorimer, D. R., & Verbiest, J. P. W. 2013, MNRAS, 431, 1352, doi: 10.1093/mnras/stt257
  • Blazek et al. (2006) Blazek, J. A., Gaensler, B. M., Chatterjee, S., et al. 2006, ApJ, 652, 1523, doi: 10.1086/507938
  • Brisken et al. (2006) Brisken, W. F., Carrillo-Barragán, M., Kurtz, S., & Finley, J. P. 2006, ApJ, 652, 554, doi: 10.1086/507765
  • CASA Team et al. (2022) CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501, doi: 10.1088/1538-3873/ac9642
  • Coleman & Burrows (2022) Coleman, M. S. B., & Burrows, A. 2022, MNRAS, 517, 3938, doi: 10.1093/mnras/stac2573
  • Cox et al. (1999) Cox, D. P., Shelton, R. L., Maciejewski, W., et al. 1999, ApJ, 524, 179, doi: 10.1086/307781
  • de Vries et al. (2021) de Vries, M., Romani, R. W., Kargaltsev, O., et al. 2021, ApJ, 908, 50, doi: 10.3847/1538-4357/abcebe
  • Deller et al. (2012) Deller, A. T., Camilo, F., Reynolds, J. E., & Halpern, J. P. 2012, ApJ, 748, L1, doi: 10.1088/2041-8205/748/1/L1
  • Deller et al. (2011) Deller, A. T., Brisken, W. F., Phillips, C. J., et al. 2011, PASP, 123, 275, doi: 10.1086/658907
  • Deller et al. (2016) Deller, A. T., Vigeland, S. J., Kaplan, D. L., et al. 2016, ApJ, 828, 8, doi: 10.3847/0004-637X/828/1/8
  • Deller et al. (2019) Deller, A. T., Goss, W. M., Brisken, W. F., et al. 2019, ApJ, 875, 100, doi: 10.3847/1538-4357/ab11c7
  • Ding et al. (2023) Ding, H., Deller, A. T., Stappers, B. W., et al. 2023, MNRAS, 519, 4982, doi: 10.1093/mnras/stac3725
  • Espinoza et al. (2022) Espinoza, C. M., Vidal-Navarro, M., Ho, W. C. G., Deller, A., & Chatterjee, S. 2022, A&A, 659, A41, doi: 10.1051/0004-6361/202142598
  • Ferrand & Safi-Harb (2012) Ferrand, G., & Safi-Harb, S. 2012, Advances in Space Research, 49, 1313, doi: 10.1016/j.asr.2012.02.004
  • Fesen et al. (1997) Fesen, R. A., Winkler, F., Rathore, Y., et al. 1997, AJ, 113, 767, doi: 10.1086/118297
  • Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  • Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  • Frail et al. (1996) Frail, D. A., Giacani, E. B., Goss, W. M., & Dubner, G. 1996, ApJ, 464, L165, doi: 10.1086/310103
  • Frail et al. (1994) Frail, D. A., Goss, W. M., & Whiteoak, J. B. Z. 1994, ApJ, 437, 781, doi: 10.1086/175038
  • Frail & Moffett (1993) Frail, D. A., & Moffett, D. A. 1993, ApJ, 408, 637, doi: 10.1086/172622
  • Gaensler & Johnston (1995) Gaensler, B. M., & Johnston, S. 1995, PASA, 12, 76, doi: 10.1017/S1323358000020075
  • Gaensler et al. (2001) Gaensler, B. M., Slane, P. O., Gotthelf, E. V., & Vasisht, G. 2001, ApJ, 559, 963, doi: 10.1086/322358
  • Green (2019) Green, D. A. 2019, Journal of Astrophysics and Astronomy, 40, 36, doi: 10.1007/s12036-019-9601-6
  • Greisen (2003) Greisen, E. W. 2003, in Astrophysics and Space Science Library, Vol. 285, Information Handling in Astronomy - Historical Vistas, ed. A. Heck, 109, doi: 10.1007/0-306-48080-8_7
  • Hailey & Craig (1994) Hailey, C. J., & Craig, W. W. 1994, ApJ, 434, 635, doi: 10.1086/174765
  • Hales et al. (2009) Hales, C. A., Gaensler, B. M., Chatterjee, S., van der Swaluw, E., & Camilo, F. 2009, ApJ, 706, 1316, doi: 10.1088/0004-637X/706/2/1316
  • Hansen & Phinney (1997) Hansen, B. M. S., & Phinney, E. S. 1997, MNRAS, 291, 569, doi: 10.1093/mnras/291.3.569
  • Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  • Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974, doi: 10.1111/j.1365-2966.2005.09087.x
  • Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  • Igoshev et al. (2022) Igoshev, A. P., Frantsuzova, A., Gourgouliatos, K. N., et al. 2022, MNRAS, 514, 4606, doi: 10.1093/mnras/stac1648
  • Johnston et al. (2005) Johnston, S., Hobbs, G., Vigeland, S., et al. 2005, MNRAS, 364, 1397, doi: 10.1111/j.1365-2966.2005.09669.x
  • Johnston & Lower (2021) Johnston, S., & Lower, M. E. 2021, MNRAS, 507, L41, doi: 10.1093/mnrasl/slab092
  • Kapil et al. (2023) Kapil, V., Mandel, I., Berti, E., & Müller, B. 2023, MNRAS, 519, 5893, doi: 10.1093/mnras/stad019
  • Kargaltsev & Pavlov (2008) Kargaltsev, O., & Pavlov, G. G. 2008, in American Institute of Physics Conference Series, Vol. 983, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and More, ed. C. Bassa, Z. Wang, A. Cumming, & V. M. Kaspi, 171–185, doi: 10.1063/1.2900138
  • Kothes et al. (2006) Kothes, R., Fedotov, K., Foster, T. J., & Uyanıker, B. 2006, A&A, 457, 1081, doi: 10.1051/0004-6361:20065062
  • Kumar et al. (2023) Kumar, P., Schinzel, F. K., Taylor, G. B., et al. 2023, ApJ, 945, 129, doi: 10.3847/1538-4357/acba93
  • Landecker et al. (1982) Landecker, T. L., Roger, R. S., & Dewdney, P. E. 1982, AJ, 87, 1379, doi: 10.1086/113227
  • Lazendic & Slane (2006) Lazendic, J. S., & Slane, P. O. 2006, ApJ, 647, 350, doi: 10.1086/505380
  • Lindegren et al. (2016) Lindegren, L., Lammers, U., Bastian, U., et al. 2016, A&A, 595, A4, doi: 10.1051/0004-6361/201628714
  • Long et al. (2022) Long, X., Patnaude, D. J., Plucinsky, P. P., & Gaetz, T. J. 2022, ApJ, 932, 117, doi: 10.3847/1538-4357/ac704b
  • Motta et al. (2023) Motta, S. E., Turner, J. D., Stappers, B., et al. 2023, MNRAS, 523, 2850, doi: 10.1093/mnras/stad1438
  • Ng et al. (2012) Ng, C. Y., Bucciantini, N., Gaensler, B. M., et al. 2012, ApJ, 746, 105, doi: 10.1088/0004-637X/746/1/105
  • Ng & Romani (2007) Ng, C. Y., & Romani, R. W. 2007, ApJ, 660, 1357, doi: 10.1086/513597
  • Popov (2023) Popov, S. B. 2023, Universe, 9, 273, doi: 10.3390/universe9060273
  • Popov & Turolla (2012) Popov, S. B., & Turolla, R. 2012, Ap&SS, 341, 457, doi: 10.1007/s10509-012-1100-z
  • Schinzel et al. (2019) Schinzel, F. K., Kerr, M., Rau, U., Bhatnagar, S., & Frail, D. A. 2019, ApJ, 876, L17, doi: 10.3847/2041-8213/ab18f7
  • Sedov (1959) Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (New York: Academic Press)
  • Shepherd (1997) Shepherd, M. C. 1997, in Astronomical Society of the Pacific Conference Series, Vol. 125, Astronomical Data Analysis Software and Systems VI, ed. G. Hunt & H. Payne, 77
  • Shternin et al. (2019) Shternin, P., Kirichenko, A., Zyuzin, D., et al. 2019, ApJ, 877, 78, doi: 10.3847/1538-4357/ab1905
  • Taylor (1950) Taylor, G. 1950, Proceedings of the Royal Society of London Series A, 201, 159, doi: 10.1098/rspa.1950.0049
  • Taylor et al. (1999) Taylor, G. B., Carilli, C. L., & Perley, R. A. 1999, Synthesis Imaging in Radio Astronomy II (Astronomical Society of the Pacific)
  • Toledo-Roy et al. (2009) Toledo-Roy, J. C., Velázquez, P. F., de Colle, F., et al. 2009, MNRAS, 395, 351, doi: 10.1111/j.1365-2966.2009.14517.x
  • Van Etten et al. (2012) Van Etten, A., Romani, R. W., & Ng, C. Y. 2012, ApJ, 755, 151, doi: 10.1088/0004-637X/755/2/151
  • Verbiest et al. (2010) Verbiest, J. P. W., Lorimer, D. R., & McLaughlin, M. A. 2010, MNRAS, 405, 564, doi: 10.1111/j.1365-2966.2010.16488.x
  • Verbunt et al. (2017) Verbunt, F., Igoshev, A., & Cator, E. 2017, A&A, 608, A57, doi: 10.1051/0004-6361/201731518
  • Vink (2012) Vink, J. 2012, A&A Rev., 20, 49, doi: 10.1007/s00159-011-0049-1
  • Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  • Wu et al. (2018) Wu, J., Clark, C. J., Pletsch, H. J., et al. 2018, ApJ, 854, 99, doi: 10.3847/1538-4357/aaa411
  • Yar-Uyaniker et al. (2004) Yar-Uyaniker, A., Uyaniker, B., & Kothes, R. 2004, ApJ, 616, 247, doi: 10.1086/424794
  • Zeiger et al. (2008) Zeiger, B. R., Brisken, W. F., Chatterjee, S., & Goss, W. M. 2008, ApJ, 674, 271, doi: 10.1086/525276
  • Zhou et al. (2023) Zhou, X., Su, Y., Yang, J., et al. 2023, ApJS, 268, 61, doi: 10.3847/1538-4365/acee7f
  • Zyuzin et al. (2018) Zyuzin, D. A., Karpova, A. V., & Shibanov, Y. A. 2018, MNRAS, 476, 2177, doi: 10.1093/mnras/sty359