The semicentennial binary system PSR J2032+4127 at periastron: X-ray photometry, optical spectroscopy and SPH modelling

TL;DR

This paper studies the PSR J2032+4127 binary system during its 2017 periastron passage, combining X-ray and optical observations with SPH modelling to understand wind interactions and disk truncation.

Contribution

It provides the first detailed multi-wavelength observational analysis of PSR J2032+4127's periastron and employs SPH simulations to model the system's geometry and wind interactions.

Findings

X-ray flux increased approaching periastron and dipped around it.

Optical H? line strength declined, indicating disk truncation.

SPH models successfully reproduce observed behaviors and geometry.

Abstract

X-ray photometry and optical spectra are presented covering the periastron passage of the highly-eccentric, ~50 year binary system PSR J2032+4127 in November 2017. This system consists of a 143 ms pulsar in orbit around a massive OB star, MT 91-213. The data show dramatic changes during the encounter as the pulsar wind collided with the stellar wind. The X-ray flux rose on the approach to periastron, then underwent a major dip in the few days around periastron, and then gradually declined over the next few weeks. The optical spectroscopy revealed a steady decline in the H? line strength on the approach to periastron (from an Equivalent Width of -15A to -7A) implying a truncation of the OB star's circumstellar disk by the approaching neutron star. Smooth Particle Hydrodynamic (SPH) modelling is used here to model the system within the context of the observed behaviour and predict the geometrical configuration of the circumstellar disk with respect to the pulsar's orbit.

Tables

Table 1. Table 1: Journal of parameters derived from the H α 𝛼 \alpha line. EW raw is a standard equivalent width measurement derived from the raw data without fitting. The subsequent columns list the properties of three Gaussian fits corresponding to the blue (1) and red (2) narrow components and the wide (3) component. The columns list central wavelength ( λ 𝜆 \lambda , Å), standard deviation ( σ 𝜎 \sigma , Å), peak height (A) and integrated area under the fits, corresponding to -EW, (Å) for each of the three components. All properties are measured with respect to a normalized continuum.

MJD	EWraw	${\lambda }_1$	${\lambda }_2$	${\lambda }_3$	${\sigma }_1$	${\sigma }_2$	${\sigma }_3$	$A_1$	$A_2$	$A_3$	-EW1	-EW2	-EW3
57504.115590	-4.70	6558.32	6566.58	6563.05	1.49	1.24	7.44	0.16	0.14	0.23	0.58	0.42	4.21
57529.199445	-6.89	6558.14	6566.36	6562.97	1.49	1.21	7.65	0.18	0.18	0.34	0.68	0.54	6.44
57559.154624	-8.17	6558.22	6566.35	6563.05	1.63	1.39	7.32	0.24	0.21	0.39	0.99	0.73	7.25
57587.112859	-8.69	6558.10	6566.58	6563.44	1.73	1.46	6.60	0.32	0.29	0.41	1.40	1.05	6.84
57592.130870	-8.57	6558.32	6566.57	6563.54	1.79	1.56	6.69	0.35	0.30	0.39	1.56	1.17	6.47
57646.018603	-6.97	6558.90	6566.77	6563.41	1.55	1.41	6.02	0.33	0.30	0.33	1.27	1.07	5.03
57675.963940	-5.33	6559.67	6566.59	6566.67	2.12	1.97	12.81	0.43	0.38	0.06	2.29	1.85	1.99
57887.118891	-14.58	6557.39	6564.36	6563.50	2.80	1.79	6.92	0.48	1.31	0.36	3.39	5.87	6.18
57901.179391	-15.01	6557.59	6564.30	6563.26	2.72	1.74	6.90	0.61	1.48	0.30	4.17	6.42	5.19
57908.149722	-15.32	6556.64	6564.32	6561.72	1.82	1.47	5.24	0.30	1.10	0.79	1.37	4.05	10.40
57913.168435	-15.10	6557.52	6563.99	6562.65	2.64	1.89	7.86	0.64	1.32	0.28	4.24	6.24	5.60
57918.187500	-14.57	6558.10	6564.28	6562.97	2.73	1.70	6.96	0.66	1.41	0.27	4.55	5.97	4.77
57926.054738	-14.33	6558.10	6564.37	6563.72	2.84	1.75	6.49	0.70	1.27	0.27	5.01	5.54	4.42
57935.133659	-14.19	6557.50	6564.56	6561.74	1.62	1.38	4.66	0.43	1.16	0.75	1.75	4.01	8.77
57998.943637	-11.81	6559.59	6565.15	6562.67	1.97	1.76	7.05	0.85	0.91	0.23	4.21	4.01	4.06
58019.913414	-11.05	6559.34	6566.02	6563.46	1.43	0.81	4.19	0.56	0.33	0.81	2.01	0.66	8.52
58025.919956	-9.96	6559.56	6566.00	6563.94	1.68	0.99	4.01	0.59	0.29	0.69	2.50	0.71	6.95
58030.955914	-9.72	6559.62	6566.04	6563.94	1.55	0.97	4.08	0.56	0.26	0.69	2.18	0.62	7.07
58051.819740	-7.50	6559.99	6566.64	6563.53	1.34	1.30	3.88	0.53	0.22	0.53	1.77	0.73	5.13
58056.825851	-7.45	6560.00	6566.53	6563.48	1.33	1.27	3.95	0.48	0.20	0.54	1.59	0.63	5.33
58058.916801	-7.07	6560.03	6566.64	6563.35	1.37	1.49	3.71	0.52	0.23	0.49	1.77	0.86	4.52
58059.935327	-7.24	6560.05	6566.48	6563.48	1.30	1.24	4.01	0.54	0.23	0.48	1.75	0.72	4.88
58061.835944	-7.21	6560.03	6566.69	6563.34	1.23	1.28	3.86	0.50	0.20	0.53	1.55	0.63	5.10
58067.823851	-7.04	6560.11	6566.75	6563.31	1.27	1.35	3.86	0.53	0.20	0.49	1.70	0.69	4.72
58069.811054	-6.90	6560.03	6566.90	6563.35	1.17	1.11	3.80	0.51	0.17	0.53	1.50	0.46	5.05
58071.832268	-6.87	6560.13	6567.08	6563.43	1.27	1.12	3.94	0.48	0.17	0.50	1.51	0.47	4.96
58072.898812	-6.81	6560.19	6566.66	6563.46	1.22	1.15	3.80	0.52	0.16	0.51	1.59	0.47	4.85
58074.870135	-6.89	6560.23	6566.89	6563.44	1.20	1.24	3.83	0.47	0.18	0.52	1.41	0.56	5.02
58076.930428	-6.22	6560.21	6567.05	6563.33	1.13	1.12	3.50	0.42	0.20	0.53	1.18	0.55	4.68
58115.828472	-11.35	6560.76	6567.17	6563.07	1.31	1.27	5.62	0.70	0.36	0.60	2.29	1.14	8.43
58212.246712	-9.32	6558.82	6565.29	6560.57	0.79	1.31	3.16	0.35	1.10	0.67	0.69	3.60	5.32
58237.167561	-7.84	6558.97	6565.25	6562.51	1.50	1.19	1.01	0.89	1.15	0.59	3.36	3.45	1.50

Table 2. Table 2: Model parameters. The inclination of the binary orbit plane to the circumstellar disk plane is characterized by the angle, β 𝛽 \beta , and the azimuth of tilt measured from the x 𝑥 x -axis (i.e., the apastron direction), is γ 𝛾 \gamma . The parameter r disk / a subscript 𝑟 disk 𝑎 r_{\mathrm{disk}}/a is the initial disk radius normalized by the periastron distance.

Model	Be disk					Stellar wind
	$\beta$	$\gamma$	$r_{\mathrm{disk}}/a$	${\rho }_0$		$\dot{M}$	$V_{\mathrm{wind}}$
				(g/cm3)		($M_{\odot }/\mathrm{yr}$)	(km/s)
Model 1	${45}^{\circ }$	$-{45}^{\circ }$	0.02	${10}^{-11}$		${10}^{-8}$	${10}^3$
Model 2	${45}^{\circ }$	$+{30}^{\circ }$	0.02	${10}^{-11}$		${10}^{-8}$	$2\cdot {10}^3$
Model 3	${45}^{\circ }$	$+{30}^{\circ }$	0.1	${10}^{-11}$		${10}^{-8}$	$2\cdot {10}^3$

Equations

$$\eta$$

$${\alpha }_{\mathrm{ff}}=\frac{2^{5/2}{\pi }^{1/2}{e}^6}{(3mkT{)}^{3/2}c}\frac{Z^2{n}_{\mathrm{i}}{n}_{\mathrm{e}}}{{\nu }^2}\bar{g}$$

$$\bar{g}=\frac{3^{1/2}}{\pi }\left[\ln \frac{(2kT{)}^{3/2}}{\pi e^{2}Zm\nu }-\frac{5{\gamma }_{\mathrm{E}}}{2}\right]$$

$${\tau }_{\mathrm{ff}}=\int {\alpha }_{\mathrm{ff}}ds,$$

$$\frac{dE}{dt}=\frac{2}{3c^2}{M}_B^2{\Omega }^{4}si{n}^2\alpha$$

Full text

††pubyear: 2018††pagerange: The semicentennial binary system PSR J2032+4127 at periastron: X-ray photometry, optical spectroscopy and SPH modelling.–References

1 Introduction

The binary system PSR J2032+4127 is now well established as a 143 ms pulsar in partnership with the massive B0Ve star MT 91-213 (Massey & Thompson, 1991; Camilo et al., 2009) in the Cygnus OB2 association. The distance to MT91-213 is estimated at 1.38 kpc from GAIA data (Jennings et al., 2018). The system has a binary period of $\sim$48 years and an eccentricity of 0.978 (Lyne et al., 2015; Ho et al., 2017). The pulsar was originally detected in a blind search through FERMI $\gamma$-ray data (Abdo et al., 2009), and subsequently the same period was detected in the radio (Camilo et al., 2009). It is also associated with a MAGIC, VERITAS andHEGRA TeV source (Aharonian et al., 2002, 2005; Aliu et al., 2014; Abeysekara et al., 2018).The only other known similar system is PSR B1259-63, which hosts a 48 ms radio pulsar and an O9.5Ve star in a binary partnership with a period of 3.4 years (see, for example, Dubus (2013) for a review of the many papers written on this system).

The most recent periastron passage occurred on 13 November 2017 (=JD 2458070) and, given the long binary period, this offered the first opportunity to study this event with an array of telescopes at different wavelengths (see, for example, (Li et al., 2018; Kolka et al., 2017; Bednarek, Banasiński, & Sitarek, 2018)). In this paper we present the X-ray measurements obtained by the Neil Gehrels Swift Observatory plus ground-based optical spectroscopy from the Liverpool Telescope. Some of the X-ray data have already been presented in Li et al. (2017); Petropoulou et al. (2018), but the entire data set is combined here with new optical data and then used as the basis for Smooth Particle Hydrodynamical (SPH) simulations to help understand the system behaviour.

X-ray spectroscopic studies of the system using XMM-Newton and NuStar from the same multi-wavelength campaign are presented in a companion paper by (2019). In addition, detailed radio observations throughout the periastron passage may be found in (2019).

2 X-ray observations

The source was monitored over the 0.3 - 10 keV range throughout the periastron passage by the Neil Gehrels Swift observatory (Gehrels et al., 2004). The observations cover over a year and typical exposure times used were in the range 1-3 ks depending upon the spacecraft visibility window at that time. The XRT lightcurve was produced following the instructions described in the Swift data analysis guide (http://www.swift.ac.uk/analysis/xrt/). This resulted in collecting up to $\sim$100 background-subtracted counts per observation at the peak brightness just before the time of periastron, and then a flux very close to zero at the lowest post-periastron point - see Figure 1. Though versions of the X-ray lightcurve have been published by Li et al. (2018) and Abeysekara et al. (2018), we reproduce the light curve again here as we wish to use it to make clear comparisons with our models later in this paper.

In addition to [MT 91] 213, the double supergiant binary system BD +40 4220 (=Schulte 5) was seen in almost every Neil Gehrels Swift Observatory observation (except when it occasionally dropped off the edge of the field of view). See Figure 2 .

The X-ray measurements reveal a steady increase in the X-ray flux of PSR J2032+4127 on the approach to periastron, followed by a sharp dip around periastron. Then a subsequent sharp recovery and slow general decline after periastron. A possible second, even sharper dip is seen about 25d after periastron at JD 2458095.

It is interesting to compare the X-ray flux before and after periastron in the context of the separation of the two objects - see Figure 3. In this plot the X-ray measurements are shown as a function of the separation in light-seconds obtained using a Keplerian model of the orbit, and assuming an orbital inclination to the line-of sight of 60 degrees and a stellar mass ratio of 11.1 (, 2019). Globally the X-ray brightness on periastron approach is consistently much higher than post periastron for the same binary separation. The only exception to this pattern of behaviour is a short period immediately after periastron when the separation is 400-600 light-seconds.

3 Optical Spectroscopy

3.1 Observations

Medium resolution optical spectroscopy of the counterpart of PSR J2032+4127 (MT91-213) was carried out using the fibre fed integral field FRODOSpec Spectrograph (Morales-Rueda et al., 2004) on the 2.0m Liverpool Telescope (Steele et al., 2004) over the period 2016 April 25 – 2018 April 28. Observations were carried out using the high resolution VPH gratings in both the red and blue arms. These provide wavelength coverage from 3900–5100 and 5900-8000Å with a resolution $R\sim 5000$. A total of 33 epochs of data were taken spread over the 2 year observing programme. Data taking was intensified around the time of periastron, although poor weather prevented some observations at that time. Total exposure time per epoch was 1800-s. This was split into $3\times 600$-s integrations to provide protection against cosmic rays.

Data reduction was carried out using the observatory supplied pipeline (Barnsley, Smith, & Steele, 2012). For each integration this involves the processes of bias removal, flat-fielding, distortion correction, fibre tracing, spectral extraction, sky subtraction and wavelength calibration. No absolute flux calibration was carried out. The three integrations per epoch were combined using median filtering to remove cosmic rays from the final spectra.

3.2 Spectral Type Reassessment

The original classification of MT91-213 was made from a low resolution spectrum presented by Massey & Thompson (1991) in their Figure 11. They obtained a spectral type of B0V and noted the line broadening due to rapid rotation. Since our spectra are of higher resolution, we decided to re-visit this classification using the blue arm spectra. Unfortunately the signal-to-noise ratio of the individual blue spectra is poor due to the low throughput of the instrument. So we therefore selected the four spectra corresponding to the time of weakest H$\alpha$ emission (Section 3.3) and combined them together to make a 2 hour equivalent exposure that best represents the spectrum of the underlying B star. This is presented in Figure 4 along with a spectrum of the B0.2V star $\tau$ Sco taken from the data set of Steele, Negueruela, & Clark (1999). Comparisons were made with these and other spectra (Walborn & Fitzpatrick, 1990). We note that the lack of HeII 4541 Å rules out an earlier classification (e.g. O9.5V). In our classification we did not use the strength of the Balmer series lines due to contamination from circumstellar disk emission. This visibly affects the H$\alpha$ and H$\beta$ lines and also contributes to some infilling of at least the H$\gamma$ line. Overall we confirm that a spectral type in the range B0V – B0.2V remains appropriate for MT91-213.

We also used our spectrum to evaluate the rotation speed of the object using the the full-width half maximum (FWHM) – rotation speed relationship of Steele, Negueruela, & Clark (1999). From the FWHM of the Hei 4026, 4387 and 4471 Å lines we derived $v\sin i=230\pm 15$ km/s, where the uncertainty reflects the spread in the values derived from the 3 lines. If MT91-213 is rotating at the mass-shedding limit, then the escape velocity off the surface = 850 km/s. Thus we can estimate that $\sin i=0.27$ or i $\geq$ $16^{\circ}$.

3.3 Emission Lines

We confine our emission line analysis to the higher signal-to-noise, red-arm spectra. In this region we see obvious emission from H$\alpha$, HeI 6678 Å and HeI 7065 Å. Line profiles for all of these lines are presented in Figure 5. Visual inspection shows how the initially symmetrical double peaked line profiles of all three lines became strongly asymmetric in the year before periastron, before starting to return to a symmetrical state post periastron.

Equivalent widths of all three lines were measured with respect to the nearby continuum. In Figure 6 we plot the equivalent widths of all three species as a function of time, and also against each other. All three are well correlated (and so therefore likely formed in similar regions of the circumstellar disk). The maximum equivalent width of all three species is around MJD 57910 and the minimum at MJD 58073 which is very shortly after periastron.

As already described, it is clear that all of the line profiles are double peaked - as expected from a disk emission. The H$\alpha$ profile also shows a very broad component at all times. We found the H$\alpha$ profile was well fit by the sum of three Gaussian functions. We fitted all of the H$\alpha$ line profiles by such a combination, using least squares minimization to derive all fit parameters for a single epoch in parallel. Unfortunately the HeI series lines were too weak to fit robustly in most cases. Our analysis therefore will concentrate only on H$\alpha$. The results of the fitting are presented in Table 1. The three components consist of two narrow components and one broad component. The narrow components (“blue” and “red”) have mean wavelengths of 6559.1 and 6566.0 Å. The mean line widths of the individual components are 1.6 (blue) and 1.4 Å (red). These standard deviations correspond to full width half maxima of 3.8 and 3.3 Å respectively.

In the upper panel of Figure 7 we plot the ratio of the peak height of the fitted blue and red components, whereas this quantity was stable in 2016, in 2017-18 significant B/R variability is apparent. The rise to a maximum of this ratio at the time of periastron is striking indicating significant circumstellar disk distortion.

The mean separation between the two narrow components is 6.9 Å which corresponds to a velocity of 315 km/s. Over the course of the entire data set this quantity decreases from a mean of $\sim 8$ Å in 2016 to $\sim 6$ Å in 2018 (Figure 7 - centre panel).

The wide component has mean wavelength 6563.2 Å and a mean full width half maximum of 12.6 Å. This width is dominated by electron scattering and should not be interpreted kinematically (Hanuschik, 1989).

In Figure 7 (lower panel) we plot how the mean wavelength of the two narrow Gaussian components varies with time. In comparison with the 2016 data (MJD<57700) the 2017 data shows considerable variation of these quantities. At the start of the observing season (MJD$\sim 57900$) the blue and red components are observed at wavelength considerably bluer than their 2016 mean. These components both then move towards the red end of the spectrum over the course of the year until by periastron they reach a maximum. By 2018 the values appear to have returned closer to their 2016 mean. In contrast the wavelength of the wide component remains stable over the whole time period.

3.4 Circumstellar disk size

For a Keplerian disk a simple estimate of the radius of the circumstellar disk ${R_{CD}}$ may be made. The peak separation of the H$\alpha$ component is 315 km/s which, if we interpret that as a Doppler shift, gives a Keplerian disk velocity of $\sim 160$ km/s. Assuming $\sin i=1$ then this corresponds to a disk radius of $R_{CD}\sim 15R_{OB}\sim 0.7AU$.

Assuming an orbital period of 48 years and an eccentricity of 0.978 then the periastron distance will be 0.73 AU. So the pulsar’s closest approach to the OB star is close to the edge of the circumstellar disk revealed by the H$\alpha$ measurements. Plus the occurrence during 2017 of a strong $B/R$ peak cycle, the equivalent width minimum at periastron and the correlated wavelength shift of the two components over 2017 all point to the influence of the neutron star passage on the disk (Reig, 2011), with tidal forces initiating a disk asymmetry (Hanuschik, 1994; Okazaki, 1997).

4 SPH Simulations

Takata et al. (2012) showed that the temporal change in the X-ray emission from PSR B1259$-$63/LS 2883 around periastron is reproduced by a model where the hydrodynamic simulations of the pulsar wind interacting with the Be-star disk and wind are used to calculate the synchrotron emission from the shocked pulsar wind. It is interesting to study whether the same approach works for PSR J2032+4127/MT 91 213 and, if it works, to probe Be-disk parameters, such as the orientation, the density distribution, and the outer radius.

4.1 Numerical Setup

The simulations are performed with the same SPH code as in Okazaki et al. (2011) and Takata et al. (2012). It was developed originally by Bate et al. (1995) and then modified so as to work for Be decretion disks and stellar winds (Okazaki et al., 2002, 2008, 2011). In the code, the Be circumstellar disk, the Be-star wind, and the pulsar wind are modeled by ensembles of gas particles of different particle masses with negligible self-gravity, whereas the Be star and the pulsar are represented by sink particles with the appropriate gravitational masses. Gas particles which fall within a specified accretion radius are accreted by the sink particle. The radiative cooling is taken into account with the cooling function generated by CLOUDY 90.01 for an optically thin plasma with solar abundances (Ferland, 1996). We set the binary orbit as the $x$-$y$ plane with the apastron and periastron in the $+x$- and $-x$-directions, respectively.

In each simulation, we first set up a Be disk around a spherical Be star, which is inclined to the binary orbit plane with the tilt angle, $\beta$, and the azimuth of tilt measured from the $x$-axis (i.e., the apastron direction), $\gamma$ (see Fig. 8 for model geometry). The initial disk has a typical density distribution in the form, $\rho_{0}(r/R_{\mathrm{OB}})^{-3.5}$, in the equatorial plane of the Be star and is in hydrostatic equilibrium in the direction perpendicular to it. Here, $\rho_{0}$ is the quantity called the base density and $R_{\mathrm{OB}}$ is the radius of the Be star. The disk temperature is assumed to be isothermal at the temperature of $0.6\,T_{\mathrm{eff}}$ (Carciofi & Bjorkman, 2006). The initial disk consists of $10^{5}$ particles.

Then, we start a simulation at $t=-0.01P_{\mathrm{orb}}$ ($\sim 175$ days prior to periastron). The mass ejection from the star, the Be wind, and the pulsar wind are turned on at that moment. We emulate the mass ejection from the Be star by injecting particles with the local Keplerian speed at a constant rate of $10^{-8}M_{\odot}\,\mathrm{yr}^{-1}$ at a radius just outside the Be star’s equator. For simplicity, we assume that the winds coast without any net external force, assuming in effect that gravitational forces are either negligible (i.e. for the pulsar wind) or are cancelled by radiative driving terms (i.e. for the Be wind). The relativistic pulsar wind is emulated by a non-relativistic flow with a velocity of $10^{4}\,{\rm km\;s}^{-1}$ and an adjusted mass-loss rate so as to provide the same momentum flux as a relativistic flow with the same assumed energy. We assume that all the spin down energy $L_{\rm sd}=1.5\times 10^{35}\,{\rm erg~{}s}^{-1}$ goes into the kinetic energy of a spherically symmetric pulsar wind. We also assume the Be wind to be spherically symmetric with the mass loss rate $\dot{M}$ and a constant velocity $V_{\mathrm{wind}}$, both of which are parameters given in each simulation. We run simulations until $t=0.01P_{\mathrm{orb}}$ ($\sim 175$ days after periastron).

We study three models, of which the parameters are summarized in Table 2. These models have different sets of the initial azimuthal direction, $\gamma$, and radius, $r_{\mathrm{disk}}$, of the Be disk and the Be-wind velocity, $V_{\mathrm{wind}}$: In Models 1 and 2, the initial disk radius is the same, but the disk azimuth and wind velocity are different. In Models 2 and 3, on the other hand, only the initial disk radius is different. In all models, the other Be-star parameters are fixed, as well as those of the pulsar.

When the global form of the pulsar wind interacting with the stellar wind is considered, a key parameter is the ratio of wind momentum fluxes, $\eta$, given by

[TABLE]

Using the parameters adopted above, we have $\eta=0.08$ in Model 1 and $\eta=0.04$ in Models 2 and 3. Consequently, the collisional shock between the pulsar and stellar winds is wrapped around the pulsar in all the three models studied below.

4.2 Hydrodynamic Interaction between the Pulsar Wind and the Disk and Wind of the Be Star

Figures 9–11 present several snapshots taken from the simulations for Models 1–3, respectively. In each figure, the upper panels show the column density along the binary orbit axis, whereas the lower panels display the density in the binary orbit plane. The snapshots are taken at $t\sim-50\,\mathrm{d}$, $t\sim-10\,\mathrm{d}$, $t\sim+10\,\mathrm{d}$, and $t\sim+50\,\mathrm{d}$ from left to right, respectively. In each panel, the white circle at the center is the Be star and the small white dot moving around it is the pulsar. The Be disk, the Be-star’s wind, and the pulsar wind can be easily distinguished by eye, because of large differences in their densities.

In Model 1 (Fig. 9), where the Be disk has $\beta=45^{\circ}$ and $\gamma=-45^{\circ}$, the pulsar is initially in the region where the wind of the Be star cannot reach owing to the shadowing presence of the Be disk, so that no head-on collision between the pulsar wind and the Be wind occurs (first snapshot). The pulsar stays in the disk shadow until $t\sim-30\,\mathrm{d}$. At this stage, the interaction surface between the pulsar wind and the disk outer radius covers only a small solid angle around the pulsar, implying that only a small fraction of the pulsar wind power is available for particle acceleration during this stage. At $t\sim-30\,\mathrm{d}$, however, the pulsar exits the disk shadow. After this the pulsar wind directly collides with the Be wind, forming an interaction surface covering a much larger solid angle around the pulsar than in the previous disk-shadow stage (second snapshot). This colliding-wind stage lasts until a few days after periastron, when the pulsar enters the shadow of the Be disk again. This second disk-shadow stage lasts only briefly. At $t\sim+10\,\mathrm{d}$ (third snapshot), the pulsar exits the disk shadow, so that the pulsar wind resumes direct collision with the Be wind, which lasts until the end of the simulation (last snapshot).

In contrast, in Model 2 (Fig. 10), the pulsar initially moves in the stellar wind (first snapshot). The ram pressure of the pulsar wind in the wind-wind collision region increases with decreasing distance between the pulsar and the Be star. This lasts until $t\sim-20\,\mathrm{d}$ when the pulsar enters the shadow region of the Be disk. As mentioned above, the interaction surface between the pulsar wind and the Be disk is small in this stage (second snapshot). Then, a few days before periastron, the pulsar exits the shadow region and resumes direct collision with the Be wind (third snapshot). This strong interaction phase continues until $t\sim+20\,\mathrm{d}$ when the pulsar enters the disk shadow again and stays there for the next $\sim 100\,\mathrm{d}$, during which the interaction between the pulsar wind and the Be disk is very weak (last snapshot).

In Model 3 (Fig.11), the initial disk has the same disk parameters as in Model 2, except that the disk radius is five times larger. The interaction between the pulsar and the Be star in Model 3 is, at first glance, very similar to that in Model 2: The pulsar wind directly collides with the Be wind for $t\lesssim-20\,\mathrm{d}$, $0\,\mathrm{d}\lesssim t\lesssim+20\,\mathrm{d}$, and $t\gtrsim+120\,\mathrm{d}$, whereas the pulsar is in the shadow region of the Be disk for $-20\,\mathrm{d}\lesssim t\lesssim 0\,\mathrm{d}$ and $+20\,\mathrm{d}\lesssim t\lesssim+120\,\mathrm{d}$. A closer inspection, however, exposes a subtle difference. In Model 3, when the pulsar is in the disk shadow, the pulsar wind termination shock, where the non-thermal X-ray emission originates, forms at a position much closer to the pulsar, and so has a larger solid angle around the pulsar, than in Model 2. We will see how this difference in the shock geometry influences the non-thermal X-ray flux in the next section.

4.3 Emission Model for Shocked Pulsar Wind

We adopt the same emission model as Takata et al. (2012) to calculate the non-thermal X-ray emission at each orbital phase, using data from the above hydrodynamic simulations. First, the simulation volume is divided into $201^{3}$ uniform grids, and at each grid point the pulsar wind pressure is calculated. Then, from the 3-D distribution of the pulsar wind pressure, synchrotron emission is locally evaluated by using an assumption for the local magnetic fields (for the calculation scheme, see (Takata et al., 2012)). In these calculations, we implicitly assume that particles are accelerated through the first-order Fermi mechanism at the pulsar wind termination shock, and have a power-law energy distribution as $f(\Gamma)\sim\Gamma^{-p}$, where $\Gamma$ is the Lorentz factor and $f(\Gamma)$ is the number density of shocked particles per Lorentz factor at each grid point. To calculate the energy distribution of the accelerated particles we adopt $p=2$ as the canonical value, and assume that 50% of the total pulsar wind pressure is due to the magnetic fields. The total emission is obtained by integrating the local emissions over the whole simulation volume.

The resulting simulated X-ray light curves for three models are compared in Fig. 12. From the difference between the Model 1 (top panel) and Model 2 (middle panel) light curves, we can easily see how the azimuthal direction of the Be circumstellar disk, $\gamma$, and the Be-wind velocity, $V_{\mathrm{wind}}$, affect the non-thermal X-ray flux. In these models with typical disk density and size, the X-ray flux is much stronger when the pulsar wind collides with the stellar wind than when it does with the Be disk. Note that the wind velocity affects only the former flux, increasing the shocked pulsar wind pressure. As a result, in Model 1, the X-ray flux is low until the pulsar exits the disk shadow and starts interacting with the Be wind at $t\sim-30\,\mathrm{d}$. Then, an outburst occurs, of which the X-ray flux has a broad peak around periastron and gradually declines afterwards, with a short break due to the second transit of the pulsar through the disk shadow region.

In contrast, in Model 2, the continuing collision between the pulsar wind and the Be wind causes gradual increase in the X-ray flux towards the peak at periastron, which is interrupted for $\sim 20\,\mathrm{d}$ just before periatron, when the pulsar moves through the disk shadow. After reaching the peak around periastron, the X-ray flux starts to go down rapidly at $t\sim 20\,\mathrm{d}$, when the pulsar enters the disk shadow, to a level much lower than the pre-perastron level. This low state continues until the end of the simulation. Note that the characteristic features of the Model 2 light curve, particularly those prior to periastron, globally agree with the observed X-ray variations, except that the model predicts the occurrence of the X-ray dip slightly earlier than the observation. This discrepancy in the timings of the simulated and observed X-ray dip may be resolved by a small change in the disk geometry adopted for Model 2, i.e., with the azimuthal direction $\gamma$ slightly larger than $30^{\circ}$ and/or the tilt angle $\beta$ slightly smaller than $45^{\circ}$.

On the other hand, the effect of the initial disk radius on the X-ray flux is rather subtle. At first glance, the X-ray light curves for Model 2 ($r_{\mathrm{disk}}=0.02a\sim 16R_{\mathrm{OB}}$) and Model 3 ($r_{\mathrm{disk}}=0.1a\sim 81R_{\mathrm{OB}}$) are very similar: a gradual increase of flux towards the peak around periastron, interrupted by a brief dip just before periastron, and a subsequent, long-lasting low state. A closer look, however, reveals that there are differences at epochs when the X-rays mainly originate from the collision of the pulsar wind with the disk outer part. In Model 3, the dip feature prior to periastron is less conspicuous than in Model 2, with small flares that partly fill in the dip. Moreover, the post-periatron decline is much slower in Model 3 than in Model 2. These differences arise from the fact that the timescale for the pulsar to cross the disk is comparable to that to strip the disk gas around the orbit. Since our simulation is non-relativistic, where the pulsar wind has the same momentum flux as, but much smaller inertia than, the relativistic wind does (e.g., Bosch-Ramon et al., 2012), the differences between the simulated X-ray light curves for Model 2 and Model 3 disk parameters could turn out to be much smaller in relativistic simulations.

In summary, from the comparison between the observed X-ray light curve and those calculated on the basis of SPH simulations, it is found that the gradual increase of flux before $t\sim-30\,\mathrm{d}$ and the peak around periastron are due to the collision between the pulsar wind and the wind of the Be star, whereas the dip just before periastron and the low X-ray state after periastron occur when the pulsar is in the shadow of the Be disk, where no direct wind-wind collision takes place. The disk orientation and density are also roughly constrained to be $\beta\sim 45^{\circ}$ and $\gamma\sim 30^{\circ}$ (see Fig. 8 for the definition of these angles) and $\rho_{0}\sim 10^{-11}\,\mathrm{g\;cm}^{-3}$, respectively. The disk size, however, is not well constrained from our non-relativistic simulations.

4.4 Free-Free Absorption in the Radio

In (2019), the monitoring observations of radio pulses at 1520 MHz show that the radio pulsation started weakening just before periastron and then were gone for 2$-$3 weeks from several days after periastron. It is of interest to see whether either model discussed above, particularly Model 2, is consistent with this radio behavior. Below we calculate the free-free absorption at 1520 MHz along the line of sight to the pulsar, using the SPH simulation data.

The free-free absorption coefficient at a radio frequency $\nu$ is given by

[TABLE]

(e.g., Rybicki & Lightman, 1979), where $m$ is the electron mass, $Z$ is the atomic charge, $n_{\mathrm{i}}$ and $n_{\mathrm{e}}$ are the ion and electron densities, respectively, $T$ is the temperature of the circumstellar gas, and $\bar{g}$ is the average Gaunt factor, which is given by

[TABLE]

for $\nu\gg\nu_{\mathrm{pe}}$, where $\nu_{\mathrm{pe}}\sim 9\times 10^{3}n_{\mathrm{e}}^{1/2}\,\mathrm{Hz}$ is the plasma frequency and $\gamma_{\mathrm{E}}\simeq 0.5772$ is Euler’s constant (Spitzer, 1978). Then, the optical depth, $\tau_{\mathrm{ff}}$, is calculated by numerically integrating $\alpha_{\mathrm{ff}}$, using the simulation data, along the line of sight to the pulsar as

[TABLE]

with $ds$ being the line element along the line of sight. We adopt that the longtitude of periastron is $40^{\circ}$ (Ho et al., 2017). The inclination angle of the binary orbit is assumed to be $70^{\circ}$.

Figure 13 shows the variations of the optical depth for the free-free absorption at 1520 MHz for Models 1-3, from which a few interesting features are observed. First, in all models the optical depth becomes very large for a while after periastron. This large optical depth is due to the occultation of the pulsar by the Be disk. Specifically, the condition of $\tau_{\mathrm{ff}}\gg 1$ necessary for completely obscuring pulsed emission is met from $t\sim-5\,\mathrm{d}$ to $t\sim+20\,\mathrm{d}$ in Model 1 (dotted line), from $t\sim-3\,\mathrm{d}$ to $t\sim+15\,\mathrm{d}$ in Model 2 (solid line), and from $t\sim-3\,\mathrm{d}$ to $t\sim+28\,\mathrm{d}$ in Model 3 (dashed line). Second, the transition between the optically thin and thick absorption states occurs more gradually in Model 1 than in Models 2 and 3, owing to the difference in disk orientation. Finally, the duration of the optically thick absorption state is longer for a larger disk size. Particularly, Model 3 predicts too-long optically thick state to be consistent with the radio monitoring observation.

In conclusion, Model 2, of which the simulated X-ray light curve is in best agreement with the X-ray observations, also provides the variation pattern of the free-free absorption that is most consistent with the radio observations.

5 Discussion

5.1 Pulsar wind contributions

We can compare the long-term X-ray flux measured by Neil Gehrels Swift Observatory for the system presented in this paper with the only other known similar system, PSR B1259-63 during a periastron passage. PSR J2032+4127 is at a distance of 1.3 kpc, whereas PSR B1259-63 lies 2.3 kpc away. So if we want to directly compare the X-ray luminosities we need to divide all the values for PSR J2032+4127 by a factor of 3.1 to correct for the relative distances. Neil Gehrels Swift Observatory measurements for PSR B1259-63 were extracted from the Neil Gehrels Swift Observatory archive for the periastron passage that took place in Dec 2010. The binary separation values for the partners in PSR B1259-63 were obtained using the published orbital elements and assuming an orbital inclination of 23 degrees and stellar mass ratio of 14.2 (Shannon, Johnston, & Manchester, 2014). The result of the scaled comparison is shown in Figure 14, from which it can be immediately seen that PSR B1259-62 is significantly more X-ray bright than PSR J2032+4127.

The luminosity from a pulsar wind is believed to stem from the loss of rotational energy from a classic dipole field which is given by Equation 5:

[TABLE]

where $M_{B}$ is the magnetic dipole moment, $\Omega$ is the spin frequency and $\alpha$ is the angle between the rotation and magnetic axis. Assuming that $M_{B}$ and $\alpha$ are the same for both systems, this implies that the energy in the pulsar wind should simply scale by the factor of $\Omega^{4}$. So putting in spin periods of 48ms for PSR B1259-63 and 143ms for PSR J2032+4127 implies a flux ratio of $\sim$76. Observationally the scaled luminosity difference between the two systems is a factor of $\sim$40. So assuming that the magnetic dipole moment does not vary much between neutron stars the explanation for this factor of $\sim$2 difference may lie with the parameter $\alpha$.

5.2 Coronal X-rays

Coronal X-ray emission is expected from isolated OB stars in the range $10^{30}-10^{32}$ erg/s (Cassinelli et al., 1994; Nazé et al., 2011). So there may well be such a contribution in addition to any energy from the colliding winds. The peak observed Swift count rate of $\sim$0.04 cts/s corresponds to a 2 - 10 keV X-ray luminosity of $\sim 4\times 10^{32}$ erg/s assuming $N_{H}=7.7\times 10^{-21}cm^{-2}$ cm and $\Gamma=2.0$ (values from Ho et al. (2017)) and a distance of 1.38 kpc. The baseline upper limit on the luminosity, well before periastron, was estimated assuming an XRT count rate limit of $\leq$0.004 cts/s. This corresponds to an limit on the luminosity of $\sim 5\times 10^{31}$ erg/s. The latter value would be consistent with coronal X-ray emission from a star of this type, B0Ve, so the increased X-ray luminosity due to the interaction between the two stars in the system peaks at the order of $\sim 3.5\times 10^{32}$ erg/s,

Incidentally the optical object in PSR B1259-63 is very similar with a classification of O9.5Ve, so similar coronal X-ray emission would be expected.

6 Conclusions

We have presented detailed optical spectroscopy and X-ray photometry from the first observed periastron passage of the semicentennial binary system PSR J2032+4127. This rare event has allowed us to study the detailed interaction of a neutron star and a circumstellar disk plus stellar wind around an OB star in a setting where normally this doesn’t occur. Our detailed SPH models have allowed us to have a better understanding of the probable orbital configuration of the system and gain a better insight into the behaviour of this colliding wind system.

Acknowledgements

We are grateful to the Neil Gehrels Swift Observatory for the valuable X-ray coverage of this system. The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council. We are also grateful to Jumpei Takata for giving us his code for calculating synchrotron emission using SPH data, and we thank Tsuguya Naito for allowing us to adapt his original figure to create Fig. 8 in this work. Fermi work at NRL is supported by NASA.