© The Authors 2024

1. Introduction

Among the pulsar population, millisecond pulsars (MSPs) are characterized by their weaker surface magnetic fields (B_psr ∼ 10⁸ − 10⁹ G), faster spin periods (P_spin ∼ 10⁻³ s), and smaller spin-down rates (̇P_spin ∼ 10⁻²⁰ s s⁻¹). The recycling scenario postulates that mass accretion from a binary companion (Alpar et al. 1982) can lead to their spin-up to millisecond spin periods. Indeed, this scenario is supported by the fact that MSPs are commonly found in a binary system (Manchester 2017).

In recent years, significant progress has been made in understanding binary systems known as “redbacks” (RBs) and “black widows” (BWs), which are also referred to as spider pulsars (Eichler & Levinson 1988; Roberts 2013). These systems consist of rotation-powered millisecond pulsars that reside in compact binary orbits. Indeed the orbital period of these systems is typically P_orb≤ 1 day (Chen et al. 2013; Breton et al. 2013), and the two components are separated by a_IB ∼ 10¹¹ cm. The companion is a degenerate or semi-degenerate star of approximately 0.01 to 0.05 solar masses (M_⊙) for BWs, and about 0.1 to 0.5 M_⊙ for RBs (Chen et al. 2013). The first black widow pulsar discovered was PSR B1957+20 (Fruchter et al. 1988), and the first redbacks discovered were PSR J0024–7204W (Bogdanov et al. 2005) and PSR J1023+0038 (Archibald et al. 2009). The term “spider systems” evokes the possible ablation of the companion star (Phinney et al. 1988; Ginzburg & Quataert 2020) due to the impact of the pulsar wind. Additionally, the notion that these systems could lead to the complete ablation of the companion star may serve as the missing link connecting binary pulsar systems to isolated millisecond pulsars.

The discovery of several dozen spider systems over the last decade has been made possible through follow-up observations of Fermi/LAT unidentified sources (e.g., Li et al. 2018; Strader et al. 2019; Clark et al. 2020). The effort involved in searching for spider systems has been motivated by indications that they may harbor massive neutron stars beyond 2 M_⊙ (e.g., Linares 2020), although it has been shown that these measurements can suffer from large systematic errors (e.g., Voisin et al. 2020a; Romani et al. 2021; Clark et al. 2023). These comprehensive multiwavelength observations have provided valuable insights into their properties and characteristics. The radio light curves show a pulsar eclipse. This could be caused by plasma ejected by the companion star (e.g., Polzin et al. 2018).

The recent emergence of a subclass of millisecond pulsars (MSPs), known as transitional millisecond pulsars (tMSPs) (PSR J1023+0038 Archibald et al. 2009, PSR J1824–2452I Papitto et al. 2013, PSR J1227–4853 Bassa et al. 2014; Roy et al. 2015) has posed significant challenges with respect to our understanding of the physical processes governing the evolution of these astrophysical objects. These systems undergo transitions between the radio-pulsar and accretion states on timescales much shorter than typical stellar evolution timescales.

The first discovered tMSP was PSR J1023+0038 (Archibald et al. 2009). Initially, optical observations in 2001 misclassified the previously known radio source as a cataclysmic variable with an accretion disk around the pulsar, based on the detection of double-peaked emission lines (Bond et al. 2002; Szkody et al. 2003). However, in 2004, the emission lines had disappeared, leading to the subsequent identification and categorization of the source as a radio pulsar (Thorstensen & Armstrong 2005; Archibald et al. 2009). PSR J1023+0038 (Papitto & Torres 2015) remained in the radio pulsar state until 2013 when the radio pulsations ceased, and there was a significant increase in optical/UV, X-ray, and γ-ray activity, indicating a reactivation of accretion on the pulsar. The transition took place within less than two weeks (Stappers et al. 2014) and the system has since remained in this state.

Multiwavelength observations from radio to gamma rays have the potential to constrain every component of the system. Pulsar timing provides accurate ephemeris of the orbital motion up to a degeneracy between masses and the inclination angle of the system. This degeneracy can be lifted thanks to the modeling of the optical light curve of the companion, leading to mass measurements which have shown that these neutron stars tend to be rather massive (e.g., Clark et al. 2023; Linares 2020; Van Kerkwijk et al. 2011). Pulsed X-ray emission was also detected during the accretion state of PSR J1023+0038, allowing for the evolution of the rotational phase during accretion to be tracked (Jaodand et al. 2016).

Unpulsed X-ray emissions revealed peaks of emissions near pulsar inferior or superior conjunctions, depending on the system, which have been interpreted as the signature of an intrabinary shock (IBS) wrapping around one or the other star depending on the balance between the two winds (Sanchez & Romani 2017; Wadiasingh et al. 2017; Gentile et al. 2014; Kandel et al. 2021; Papitto & Martino 2021). Indeed, the tangential discontinuity separating the companion and the pulsar wind is an efficient site of particle acceleration (Harding 1990; Arons & Tavani 1993; Cortés & Sironi 2022) as the shocked material from the pulsar wind can reach relativistic bulk velocities leading to a non-thermal emission. Synchrotron radiation (SR) and inverse compton (IC) are two processes responsible for the emission of these high-energy photons. In the accretion state of tMSPs, the disk is characterized by flaring light curves in optical and X-ray (Bogdanov et al. 2015; Linares et al. 2022). In the case of PSR J1023+0038, the onset of the last accretion state has been accompanied by a fivefold increase in unpulsed gamma-ray flux the origin of which is yet unclear (Stappers et al. 2014). The unique behavior of tMSPs together with this consistent set of multi-wavelength observations offers an opportunity to study wind-wind interactions and accretion physics in the same system, on top of providing strong evidence in support of the recycling scenario. The various states observed in spider systems are undoubtedly influenced by the balance between the ram pressure exerted by the pulsar wind on the matter lost by the companion. The interplay of gravitational and rotational forces significantly shapes the flow of stellar matter (Wadiasingh et al. 2018).

Two main scenarios have been studied: the case of an accreting neutron star and the case where the pulsar is active, generating a wind that interacts with the wind of the companion that results in a shock surrounding the companion star (Benvenuto et al. 2015; Sanchez & Romani 2017). The interplay between these forces drives the observed transitions and underlies the complex behavior exhibited by tMSPs. Until now, the numerical investigation of the interaction between the wind of a companion star and a pulsar wind in a binary system has predominantly focused on systems involving a massive star. These studies have employed non-relativistic smoothed particle hydrodynamic (SPH) simulations in 3D (Romero et al. 2007), as well as classical and relativistic hydrodynamic simulations in 2D (Bogovalov et al. 2012; Lamberts et al. 2013), where the orbital motion and considered axisymmetric cases were neglected (Paredes-Fortuny et al. 2015). Recent advances include 3D relativistic simulations, including the effects of orbital motion, conducted by Bosch-Ramon et al. (2015), Huber et al. (2021).

In this paper, we employ hydrodynamical numerical simulations to investigate the wind-wind interaction in 2D, at the equatorial plane between the pulsar wind and the wind from the low-mass companion star. Bearing in mind the transitional nature of tMSPs, we focus on the region of the parameter space surrounding the limit between accreting and non-accreting neutron star (NS) behavior. Additionally, we model the X-ray light curve in the IBS regime.

The paper is structured as follows. In Sect. 2 we describe the physics of the numerical model we used in our study. This includes the description of the numerical model of the pulsar and companion wind. In Sect. 3, we summarize the IBS modelization, as well as the model of SR emission used to compute the associated X-ray light curves. The simulation setup is described in Sect. 4. Finally, we discuss the outcomes of our study in Sect. 5 and present the application for the observations of the pulsar spider systems in Sect. 6.

2. Physical and numerical model

To construct a comprehensive model of spider systems, it is essential to explore various physical processes, phenomena and interactions involved in the system. These include the irradiation of the donor star’s atmosphere by the pulsar, the transfer of gas from the companion star and the subsequent formation of the accretion disc, as well as the wind-wind interaction between theflow from the companion and the pulsar wind. By investigating these aspects, we aim to develop a coherent understanding of the complex dynamics at play in the system (see Table 1).

2.1. Orbital configuration

In this study, we chose to investigate the case of a RB system prone to a transitional behavior, consisting of a pulsar with a mass of M_psr = 1.4 M_⊙ and a companion star with a mass of M_c = 0.4 M_⊙.

The separation distance between the pulsar and the companion star is a_IB = 10¹¹ cm, resulting in a revolution time of s (equivalent to 3.56 hours), where G is the gravitational constant. Additionally, the center of mass of the system is located at a distance of r_center ≈ 0.78a_IB from the companion star.

In our analysis, we assume a circular orbit, following the common circularization pattern observed in spider systems. It is worth noting that in such systems, the companion star is typically tidally locked (Hurley et al. 2002); however, the eccentricity is usually very small (e.g., Voisin et al. 2020b).

2.2. Winds

We consider the fact that the winds from the pulsar and the companion are non-magnetized, isotropic, supersonic, and originate from a surface surrounding each one of them. The relation between the wind mechanical luminosity L_w, mass flux ṁ_w, and speed v_w is:

(1)

At the first level, the interaction between the pulsar wind and the companion wind can be characterized by the ratio of their momentum fluxes (Canto et al. 1996; Parkin & Pittard 2008; Bogovalov et al. 2012; Sanchez & Romani 2017), given by:

(2)

where ṁ_c,w and v_c,  w are the mass flux and speed of the companion wind, and ṁ_psr,w and v_psr,  w are the mass flux and speed of the pulsar wind.

The influence of the orbital properties on the wind from the companion star and on the wind from the NS can be characterized by two dimensionless quantities. Both are given at the surface where the wind is set. The first expresses the binding of the flow to the system, given by the ratio between wind speed and the escape speed v_esc as follows:

(3)

The second expresses the influence of the rotation of the system, specifically the Coriolis force. It is determined by the ratio of the wind speed to the rotation speed v_rot ∼ 490 km/s:

(4)

2.2.1. Pulsar wind

The supersonic pulsar wind is emitted from the surface of a sphere with a radius of R_psr, w = 0.025a_IB, which surrounds the pulsar, as demonstrated by the model outlined in (Meyer et al. 2022).

Throughout this paper, we consider a pulsar spin-down luminosity of L_sd = 10³⁵ erg/s, which is typical of known tMSPs, where a fraction (1 − b) = 0.1 accounts for the radiated luminosity in γ-rays. The radiative luminosity in the hard X-ray band is assumed to be negligible compared to the γ-ray luminosity. The speed of the pulsar wind is set (as seen in Table 2) so that the physical hierarchy of the system is respected, that is, the effects of gravity and rotation on the pulsar wind are nearly negligible.

However it’s important to acknowledge that this speed is significantly lower than the realistic pulsar wind speeds, which can reach a Lorentz factor of 10⁶ (e.g., Kennel & Coroniti 1984). The choice of a reduced pulsar wind speed can lead to significant changes in its properties, affecting compression rates, velocities, and consequently the positioning of shocks. It can also affect the development of associated instabilities. Since the momentum ratio is expected to primarily determine the behavior of the interaction between the two winds (Wilkin 1996), the mass flux of the wind ṁ_psr,w is set such that momentum flux of a relativistic wind with a mechanical luminosity bL_sd is reproduced, following:

(5)

We consider that this setup already provides a valuable analog to the actual system at a significantly reduced computational cost compared to a relativistic version. The pulsar wind in our simuations is set with a mechanical luminosity following Eq. (1). In the 2D approach the mass flux at a distance r is: . To reproduce correct momentum flux at a distance r = a_IB, the initial mass flux of the wind, , setting L_psr, w = 5 × 10³¹ erg/s.

In addition we make the assumption that all Poynting flux has been converted into kinetic flux¹. However, it is important to note that a realistic pulsar wind is expected to be strongly asymmetric. This asymmetry arises from the angular distribution of the Poynting flux, which follows a sin²θ profile according to the analytical model (Michel 1973), and even a sin⁴θ profile based on numerical simulations (Tchekhovskoy et al. 2013), where θ represents the polar angle.

2.2.2. Companion wind

In this work, we explore the effect of varying the two parameters characterizing the companion wind. The mass-loss rate is constrained by stellar evolution considerations to remain below ṁ_c,w,max = 10⁻⁸ M_⊙/yr (Chen et al. 2013; Benvenuto et al. 2015). In order to reproduce the correct hierarchy of speeds between the two winds we set v_c, w < v_c, max = 300 km/s = 0.03 v_psr, w. This limit is not very constraining, as larger speeds would cause the matter to escape in all directions instead of being constrained to flow through L₁.

The wind of the companion is characterized by a Mach number ℳ_c, w such that the surface temperature of the companion is T_c ∼ 4480 K, typical value for RB systems. We explore distinct scenarios by varying the initial speed v_c,  w (with values between 60 and 200 km/s) and the mass flux of the wind ṁ_c (with values between 4.5 × 10⁻¹⁰ M_⊙/yr and 9 × 10⁻¹⁰ M_⊙/yr). This is summarized in Table 3.

The wind of the companion is set with sufficient specific kinetic energy to overcome the gravitational potential barrier to the Lagrange point L₁. This criterion is satisfied by considering the conservation of Bernoulli energy along the flow streamline, given by:

(6)

where γ represents the polytropic index, v_c, w denotes the wind speed of the companion in the rotating frame, and c_s, w = v_c, w/ℳ the sound speed of the companion wind. Lastly, Φ represents the effective gravitational energy in the rotating frame, accounting for the gravitational potential energy of the system, it is given in case of synchronous rotation by:

(7)

where r_psr, r_c and r_center are respectively the distance from pulsar, the companion star and the center of mass.

The fluid elements at the surface of the companion star that can escape the gravitational potential and cross the Lagrange point L₁ must have a minimal specific kinetic energy. Thus, the Bernoulli energy at the Lagrange point, L₁, is then ξ = Φ. As a consequence, the wind speed at the surface of the companion star must exceed the escape speed given by:

(8)

where Φ_L1 represents the gravitational potential at the Lagrange point L₁, ℳ the mach number of the companion wind at its surface, and Φ_c, surface the gravitational potential at the surface of the companion star.

In this paper, the wind emanates from the equipotential surface with an effective potential Φ_c, surface = Φ(x = R_c, w, y = 0), ((x, y) coordinates are defined in Fig. 1).

Fig. 1. Sketch of the grid setup used in simulations. The companion star is at the center of the cylindrical grid, and the pulsar at a distance aIB along X axis. Wind injection surfaces are shown in red, gravitational potential of Lagrange 1 point is shown in blue. Top left arrow indicates sens of rotation.

The influence of the gravity on the companion star flow is then given by the ratio between wind speed and the threshold speed necessary for the wind to fill the Roche Lobe and cross the Lagrange point L₁ (Eq. 8):

(9)

The wind is set with speed such that mass transfer occurs through Roche-lobe overflow (RLOF), so the ξ_c,  grav is set to be on the order of 1. As a result, the dimensionless rotational parameter ξ_c,  rot is about one-third, so the rotation is expected to play a significant role.

We define the Roche lobe filling ratio based on the position of the surface along the line of center, R_c, w = 0.3a_IB, which corresponds to about 80% of the Roche lobe along this direction. Indeed, as the companion star is typically irradiated by the pulsar wind, this leads to its expansion, resulting in the filling of a significant portion of its Roche lobe, exceeding 90% (e.g., Crawford et al. 2013) and also exceeding the radius of an isolated main-sequence star with a similar mass. This behavior is consistent with a quasi-RLOF state, facilitating mass and angular momentum transfer from the companion star to the pulsar (Benvenuto et al. 2015).

Furthermore, the companion star is tidally locked to the pulsar (Phinney et al. 1988). Combined with irradiation, this results in a potentially strong anisotropy of the surface temperature of the companion with a hotter “day” side. However, for the purposes of this initial study, our main focus is on the effect of wind-wind interaction. Therefore, we do not consider this effect in the cases presented in this paper.

Our current HD approach does not allow us to capture the dynamics associated with the potential companion magnetic field. Some studies have suggested that it may be an important element in high-energy emission scenarios (Sim et al. 2024) or in irradiation patterns (Sanchez & Romani 2017).

2.3. Heating of the companion star’s wind

The observed high effective surface temperature of the companion star may be influenced by various factors. These factors are likely to act in a complementary manner. One proposed mechanism is the irradiation by high-energy pulsar photons, considering that a fraction (1 − b) of the pulsar’s spin-down luminosity L_sd is deposited at the surface of the companion star (Benvenuto et al. 2015).

Another proposed mechanism is the irradiation feedback process, as suggested in Hameury & Ritter (1997), where a fraction of the released accretion luminosity is irradiated onto the photospheric layers of the companion star, primarily in the X-ray regime.

Additionally, a portion of the pulsar wind particles can undergo reprocessing in the intrabinary shock and travel along the magnetic field lines of the companion star, reaching its surface. This mechanism, proposed to explain the heating of low-mass companions in black widow systems (Sanchez & Romani 2017; Yap et al. 2019), involves the irradiation of the companion star. Furthermore, direct irradiation resulting from the shock region between the pulsar wind and the wind of the companion has been proposed. It should be noted, however, that this mechanism alone is unable to reproduce the observed modulation in the light curve (Bogdanov et al. 2011). These various mechanisms can collectively contribute to the heating and increase in the surface temperature of the companion star. In this work, we investigate the effect of heating the wind upstream of the companion by the pulsar wind by using the description given in Tavani & London (1993):

(10)

with L_sd the spin-down luminosity of the pulsar, σ_T the Thomson cross section, and m_p the proton mass. The heating parameters Υ is the measure of the quality of irradiating spectrum of X-rays and γ-rays, thus Υ_γ = (1 − b) = 0.1 and Υ_X/Υ_γ ≪ 1 as luminosity in the X-rays is assumed negligible. T_c ∼ 10⁸ K is the Compton temperature and T_x ∼ 5 × 10⁶ K is the temperature at which the relevant heavy elements are almost completely ionized. These temperatures act as thresholds below which we consider the heating to be negligible.

This simplified approach cannot capture all relevant physics. Notably, the transferred matter through L1 may exhibit different physical properties, such as being in an optically thick regime. While we cannot account for such physics, we consider this approximation to be sufficient for simulating the interaction between the two winds at a significantly reduced computational cost compared to more comprehensive methods.

We also introduce the power ratio, which relates the companion wind power to the fraction of pulsar spin-down power geometrically intercepted by the companion (assuming spherical symmetry of the pulsar wind):

(11)

where f = R_c/(2πa_IB) is the fraction of power intercepted by the cross section of the companion in the equatorial plane with R_c the companion star radius, setting:

(12)

2.4. The cooling

We also incorporated an optically thin radiative cooling function Λ(T) that depends on the local temperature T. We adopted the cooling function Λ(T) described in the study by Schure et al. (2009), which we calibrated for solar metallicity. Moreover, we maintained a floor temperature of 100 K.

3. Intra-binary shock and X-ray emission

3.1. IBS analytical model

The study of the interaction of two supersonics, isotropic stellar winds has been done in the past. Numerical simulations of the collision between two massive-star winds have been made using full gas-dynamics equations in the radiative regime in Stevens et al. (1992). In this regime, the shocked region is thin and its locus can be brought down to the surface of ram-pressure balance as done in Dyson et al. (1993).

Analytical expressions of the shock surface between two spherical winds in the thin-shell regime have been discussed in Canto et al. (1996), as an extension of Wilkin (1996). These neglect the effect of orbital motion. The derivation rests upon conservation of mass, linear and angular momenta through the surface defined by ram-pressure balance. In this section we summarize the results of Canto et al. (1996) that are relevant to our study.

The result is a surface of revolution around the intra-binary axis. The surface intersects the axis at the stagnation point located at a distance r_sh from the companion such that,

(13)

where β is the momentum flux ration defined in Eq. (2). The surface is asymptotically conical at large distances. Its opening angle θ_∞ with respect to the companion-to-pulsar direction is given by the solution of:

(14)

For example, when β → 0 (pulsar wind infinitely dominating) we have θ_∞ → π, corresponding to a bow shock wrapping around the companion and asymptotically cylindrical.

These models offer a simple solution, providing an asymptotic angle for approximating the solid angle of Doppler emission within the shock. However, these models consistently neglect the influence of gravity and orbital motion.

In the case of a close and thus fast-rotating binary system, such as a redback system, the Coriolis force is important as it breaks the cylindrical symmetry about the intra-binary axis, causing the stagnation surface to wrap and bend around the companion (Lemaster et al. 2007).

In this study, we perform hydrodynamical simulations of two isotropic winds, incorporating the effects of both gravity and orbital motion. It has been shown in Bogovalov et al. (2012) that a hydrodynamical approach provides a viable approximation for numerically modelling the interaction between a weakly magnetized pulsar wind and a stellar wind.

3.2. IBS X-ray emission

The intrabinary-shock is an efficient site of particle acceleration and X-ray emission (Harding 1990; Arons & Tavani 1993). The pulsar wind reaches relativistic bulk speeds along the shock, resulting in relativistic beaming of emissions by the shocked material. As a consequence, emissions from the shock are modulated with the orbital phase. Due to the beaming effect, the geometry of the shock surface and the inclination angle of the system largely determine the shape of the X-ray light curve. The X-ray flux is then showing a maximum at the inferior conjunction of the orbit for RBs (shock surrounding the pulsar), and a minimum for BWs (shock surrounding the companion). The single or double peak we can observe depends mainly on the opening angle of the intra-binary shock.

In this study, we model the X-ray emission within the shock. For this purpose we evaluate the SR of an electron population in the IBS and solve the radiative transfer equation along different lines of sight, allowing us to study the link between the geometry of the shock and the associated light curve in X-rays.

Taking physics variables from our simulations allows us to differ from previous analytical works. First of all the direction of particle propagation no longer strictly adheres to the shock control line, but instead aligns with the velocity of the flux within the IBS. Secondly the emission in the shock depending on the electron population is modeled using the density in the simulation.

We consider a relativistic electron population in each shock cell i, set with a power law as:

(15)

with γ_e, min < γ_e < γ_e, max. The constant p depends on the process of particle acceleration. X-ray spectra of tMSPs show a hard spectrum corresponding to a particle distribution index p ∼ 1.3 (Papitto & Martino 2021, and references therein). As already noted in (e.g., Kandel et al. 2019), a value close to 1 suggests that the dominant mechanism of particle acceleration is magnetic reconnection (Sironi & Spitkovsky 2014). However, determining a dominant acceleration process based solely on a hard spectrum is challenging, especially when considering the spectrum of injected particles from the pulsar wind (van der Merwe et al. 2020). As a result, the acceleration process is not considered in this work. For simplicity, we adopt the fiducial value p = 1.3. N₀ is the normalization coefficient, defined as (Gómez et al. 1995):

(16)

where e_th, e = ϵ_ee_th and n_e = ϵ_en. e_th is here the thermal energy, n the number density and ϵ_e = 0.1 represents the fraction of electrons contribution. We set γ_R = γ_e, max/γ_e, min, and γ_e is obtained via the observed frequency of synchrotron in X-ray band, ν_{min − max} = 2.4 × 10¹⁶ − 2.4 × 10¹⁸ Hz (corresponding to the range of energy 0.1 − 10 keV; Bogdanov et al. 2009; Huang et al. 2012):

(17)

If we consider the magnetic field along the shock is defined (assuming a split monopole geometry Thompson et al. 2004; Weber 1967) toroidal outside the light cylinder as:

(18)

with

(19)

Using the typical values for MSPs (P = 10⁻³ s, ̇P = 10⁻²⁰s s⁻¹, r_surf = 10 km), we obtain B_surf ∼ 10⁸ G. With r_LC = c/Ω = 50 km, we have B_LC ∼ 8 × 10⁵ G. For each simulation we compute emission assuming a magnetic field at the stagnation point of B_sh = 80 G, meaning γ_e ∼ 10⁴ − 10⁵.

We also consider the radiative losses to be instantaneous and the synchrotron emissivity occurring in the electron rest-frame as:

(20)

with

(21)

We assume that particles move in a circular motion around the magnetic field lines in the restframe, meaning that the pitch angle α is 90°. The total specific emission in each cell is then computed as:

(22)

In order to obtain the specific intensity received by an observer, the quantity of Eq. (22) is transformed,

(23)

where δ is the Doppler factor:

(24)

with Γ the bulk Lorentz factor of pulsar material, θ the angle between the flow direction and sky position from the considered line of sight. The flow direction is directly computed using the velocity of the matter in the simulation.

The orbital variability of the total emission is obtained by solving the radiative transfer equation along a line of sight for a different viewing angles during all the orbital period, following:

(25)

We consider the optical depth τ_ν ≪ 1.

Since we are in 2D, the stars always eclipse the line of sight. The radius of the stellar companion being far more predominant than the neutron star (NS), we consider the emission to be blocked only if the star is in the LOS. The associated X-ray light curve of IBS occurring in our simulations are computed to constrain our model.

4. Simulation setup

The simulations were conducted using the Message Passing Interface-Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC; Keppens et al. 2021). We incorporated the treatment of a rotating frame around the center of mass, which can be shifted from the grid center, into the hydrodynamic module of the AMRVAC code.

In this modified version of AMRVAC, we solved the conservation equations for mass, momentum, and energy. The effects of the gravity of both stars and the rotation of the reference frame were included as source terms in the momentum and energy equations. Additionally, heating and optically thin radiative cooling were accounted for through source terms in the energy equation. The following equations were solved:

(26)

(27)

(28)

where ρ and P represent the fluid density and pressure, respectively; v denotes the fluid velocity vector in the reference frame. The total energy in the co-rotating frame is given by , where e_th represents the thermal energy. In this paper, we assume a polytropic equation of state such that and γ = 5/3 the polytropic index. After conducting some tests, we did not include viscosity in the above equations because it did not have a significant effect over the simulation timescales.

The gravitational potential is defined as , where r_c represents the distance from the companion star and r_psr represents the distance from the pulsar. Regarding the rotating reference frame, it has an angular velocity Ω centered on the central mass, and within the grid, the distance from this center is denoted as r_center. Regarding the heating, our study primarily investigates the impact of X-ray and γ-ray irradiation from the pulsar on the companion wind. This is represented by the source term S_psr (Eq. 10). The cooling term Λ(T) is discussed in Sect. 2.4.

The simulations are conducted on a 2D polar grid, covering a radial extent of r = [3 × 10⁹, 2 × 10¹¹] cm and a toroidal extent of ϕ = [0, 2π] (see Fig. 1). A logarithmic grid is employed in the radial direction. Initially, the grid is refined to the first level with a resolution of (32 × 64) cells. Given the substantial scale difference between the shock scales, the pulsar wind zone, and the overall simulation box investigated in this project, the use of adaptive mesh refinement (AMR) is essential to accurately resolve the significant variations encountered in the spider system with a pulsar. In our simulations, the grid is allowed to be refined up to 9 additional levels, with a doubling of resolution at each new level of refinement. After studying the impact of resolution, we determined that this level of resolution is the minimum necessary to adequately capture the physics close to the pulsar. The refinement and coarsening processes utilize Lohner’s error estimator to assess second-order variations in density and radial velocity. Moreover, the refinement is only allowed in regions where the winds are encountered.

Regarding the internal and external boundary conditions, the companion star is placed at the grid center, and its wind is imposed within the equipotential surface with an effective potential Φ_c, surface. Concerning the boundary condition around the NS, during the accretion phase study, an inflow Neumann boundary condition projected on the radial direction in the pulsar’s frame is applied. The radial outer boundaries are treated as open boundaries, utilizing a Neumann boundary condition with vanishing derivatives for all quantities. During the active pulsar phase the wind is emitted from the surface of a sphere with a radius of R_psr, w. The equations are solved using the total variation diminishing (TVD) method, specifically the Harten-Lax-van Leer-Contact (HLLC) method (Harten 1997), combined with the Koren flux limiter (Vreugdenhil & Koren 1993).

5. Simulations of characteristics states and transition

In the context of transitional systems, there is an intricate progression from a LMXB state, characterized by discernible accretion onto the NS, to a pulsar state where radio emissions remain unobscured Stappers et al. (2014). Within this section, our principal aim is to replicate both of these states through the utilization of our computational model. Subsequently, we focus on the influence of the companion wind characteristics on the state of the system and the corresponding observables.

Here, we study different stellar wind configurations, each with a different effect on the behavior of the system. Our overarching objective is to delineate the critical threshold that marks the transition from a pulsar accreting state to a radio pulsar state.

5.1. Accretion stream without a pulsar wind

First, we set the parameters for an accreting neutron star and the wind from the companion star. Specifically, the wind of the NS is inactive and the companion wind is adjusted to ensure that mass transfer occurs by RLOF. Table 3 presents the matrix of cases we have studied. Rows are labeled from “0” through “3” in order of increasing companion wind speed v_c,  w, while columns labeled “a” and “b” are in order of increasing companion mass flux, ṁ_c,w.

Due to the effects of gravity and orbital motion, the supersonic accretion stream leaving L₁ does not directly collide with the surface of the pulsar but rather orbit it eccentrically. Then if the dynamics is dominated by the gravity of the neutron star, the angular momentum of the stream with respect to the neutron star l_stream gives the circularization radius following (Hendriks & Izzard 2023):

(29)

This leads to a resulting circularization radius of , with q = M_psr/M_c = 3.5 the mass ratio of the system and x₁ ∼ 0.37 the distance of L₁ from the companion star (Lu et al. 2023). In our simulations we obtain radii within a range R_circ, b0 ∼ 0.92R_circ and R_circ, b3 ∼ 1.35 R_circ (R_circ, b1 ∼ 0.98 R_circ and R_circ, b2 ∼ 1.1 R_circ). In fact, as the speed of the companion wind increases, the angular momentum of the accretion stream increases, so that it circulates at a greater distance from the neutron star.

Due to the rotation of the system, the returning stream then collides with itself (the incident stream). This interaction redirects the subsequent paths of the incident stream and the returning stream following the rotation direction as illustrated in Fig. 2. The impact does not disrupt the upstream trajectory of the incident stream since no disturbance can propagate upstream in supersonic flows. However, each new collision induces additional deflections that result in the circularization of the flow, following the dynamics explained in Lubow & Shu (1975). Within these interacting streams, a residual disk of material forms from earlier interactions, with the material slowly spiraling inward toward the NS.

Fig. 2. Comparison of the accretion stream in cases (b1) and (b3) without pulsar wind, that is for two different companion wind velocities and equal companion mass loss. Left is the companion and right is the NS. Purple dashed lines indicate boundary surfaces, while the blue dashed line indicates the isocontour of the effective gravitational potential through the L1 Lagrange point. The black arrow shows the direction of orbital rotation. The color scale gives matter density.

The variability in the mean density enveloping the pulsar underscores the interaction between the returning stream and the incident stream (Fig. 3, left column). Within a ring delimited by the internal boundary of the pulsar wind and the estimated circularization radius R_circ, two types of mass loss are observed. The first one, occurring twice per period, is associated with an angular momentum variation, indicating changes in the circularization radius over time as the stream approaches and then moves away from the control ring. The second one is manifesting approximately 10 times per orbital period which corresponds to the Keplerian period at the circularization radius. Notably, in the case of (b3), while the second mass loss persists, the first one is less pronounced, indicating a stationary accretion stream state. This can be attributed to the increased density within the flow, resulting in less perturbation of the accretion stream. This stationary behavior is particularly pertinent when examining the interaction with the pulsar wind.

Fig. 3. Variability of the accretion stream regime with and without pulsar wind. The panels show, from top to bottom, the temporal variation of mean density (in g cm−3) within a ring, extending from the internal boundary of the pulsar wind to the circularization radius Rcirc. Left column correspond to simulations without a pulsar wind (see Sect. 5.1), while right column to simulation with both winds. First line corresponds to scenario (b1), second line to scenario (b3).

The increase in companion wind speed leads to an increase in the specific kinetic energy and thus Bernoulli energy within the accretion stream. This results in a lateral expansion of the stream over a wider equipotential surface in the vicinity of L₁. Specifically, the width of the accretion stream is 3 times larger in case (b1) than in case (b0), and 3.2 times larger in case (b3) than in case (b0). Additionally, examining cases (a0) to (a3) reveals that the thickness remains consistent with cases (b0) to (b3), respectively (from 0.048 a_IB to 0.288 a_IB), and that this parameter is only influenced by the speed.

Another consequence of this phenomenon is the enhanced ram pressure within the accretion stream. The total ram pressure in the accretion stream shows a direct correlation with the initial speed of the stellar wind.

Near L₁, in the inertial frame, the kinetic energy of the accretion stream is comparable to the thermal energy, with Mach numbers ℳ_L1 < 1 (between ℳ_L1, b0 = 0.7 and ℳ_L1, b3 = 0.3). As the stream approaches the neutron star, it is accelerated, with the acceleration being more pronounced for scenario (b0) compared to scenario (b3). Consequently, the transition to the supersonic regime occurs earlier at lower speeds. As the flow is accelerated along the stream, the energy becomes mainly kinetic when approaching the NS, accounting for more than 90% of the energy for (b0) and 70% for (b3).

Within the accretion stream, the Bernoulli energy (Eq. 6) remains constant along the streamlines, indicating an adiabatic flow. However, the thickness of the stream decreases under tidal effect, while the mass flux stays constant, causing the overflowing stream to become more confined near the stream center.

Regarding mass loss in the system, apart from the matter not accreted by the NS and ejected in the close environment, we observe a condensed mass loss phenomenon at the L₂ point behind the companion star. This differs from the perspective proposed by Lu et al. (2023), which focuses on potential mass loss of an accretion disk at the L₃ point. However, it is noteworthy that the total mass lost remains negligible when compared to the mass efficiently captured by the NS gravitational field.

5.2. Pulsar wind influence on accretion stream

The pulsar wind is modeled as described in Sect. 2.2.1 and the settings are summarized in Table 2. Figure 4 shows the evolution of the scenario (b1), which serves as a representative example of all simulations involving the interaction between the accretion stream and the pulsar wind. These cases are represented by orange cells in Table 4.

Fig. 4. Visualization of the interaction between the accretion stream and pulsar wind for the case (b1). Left is the companion and right is the NS. Purple dashed lines indicate boundary surfaces, while the blue dashed line indicates the isocontour of the effective gravitational potential through the L1 Lagrange point. The color scale gives matter density. Green text top left indicates the time in unit of orbital period, P.

The fast pulsar wind reaches the Roche Lobe of the companion star and a shock forms around the companion star (see panel (1)). If the ram pressure within the stream is sufficiently high, the accretion stream emerges (panel (2)) and reaches the NS (panel (3)). The circularization radius remains the same across the different stellar parameters and is smaller compared to inactive pulsar wind scenarios (Sect. 5.1), with R_circ, b1 ∼ R_circ, b3 ∼ 0.8 R_circ. The influence of the pulsar wind leads to a decrease in the velocity of the companion wind, which in turn leads to a decrease in the specific angular momentum and a smaller circularization radius. The interaction with the pulsar wind also causes the stream center to be shifted from the L₁ point by approximately 0.05 a_IB in the direction orthogonal to the line of centers (visible in panel (7)). On top of that we observe that the transition to a supersonic regime is now occuring precisely at L₁, ℳ_L1 ∼ 1.

During the third phase, the fast pulsar wind interacts strongly with the accretion stream. It compresses it laterally, pushing it further away from the pulsar. The closest distance between the accretion stream and the pulsar varies between 0.1 − 0.2a_IB (panels (4) and (5)). During this variable phase, the accretion stream continues to wrap around the pulsar until the head of the accretion stream completes its tour around the pulsar and intersects with the side of the incoming stream at about 0.85P. At the same time, the stream rotates and converges toward the pulsar and confines the pulsar wind (panel (6)). The system then enters an accretion stream regime, akin to what was observed with an isolated companion wind, the pulsar wind now mixed in the accretion flow surrounding the NS (panel (7)).

Unlike the periodic behavior we used to observe considering an isolated companion wind (see Sect. 5.1), the accretion stream state, when interacting with a pulsar wind, exhibits increased instability. We now witness seemingly chaotic behaviour on timescales of a fraction of orbital period (Fig. 3, right column). Due to those instabilities, the ram pressure of the pulsar wind close to the NS on occasion surpass that of the accretion stream. In the case of the weak accretion stream (b1), around time t ∼ 1.4 P, the pulsar wind significantly disturbs the accretion arm, pushing it outward. This process returns the system to a state similar to the initial transition phase at time t ∼ 0.8 P (as seen in panels (8) and (9)). This initiates an unstable behavior in which the accretion stream is expelled several times until the returning stream successfully catches up with the incident one and the system returns to an accretion stream state akin to panel (7) (at about t ∼ 5 P), before undergoing disruption once again. This instability highlights the intricate interaction between the two winds.

However, scenarios (b2) and (b3) exhibit a more stable accretion stream state, remaining relatively unaffected by the pulsar wind, thereby maintaining the configuration shown in panel (7), Fig. 4. The system then enters into a quasi-stationary state, where instabilities quasi-periodically lead to the accretion stream being disrupted away before gradually rebuilding (Fig. 3, bottom right panel). However the matter surrounding the NS is not entirely expelled, maintaining the system in a accretion stream state. These instabilities occurring around the pulsar are intricately tied to the accumulation of matter that remains unaccreted. This accumulation leads to an increase in the total angular momentum, reaching a point where it becomes too high to sustain a stable flow. Consequently, mass expulsion occurs followed by renewed accumulation, effectively restoring the system to its previous density. Therefore the rate at which these instabilities manifest is directly related to the pace at which matter accumulates. Faster accumulation results in more frequent instabilities. In Sect. 5.1, we demonstrated that increasing the initial speed of the companion wind leads to an increase in density within the stream. This results in an increase in the frequency of instabilities. Indeed in our simulations, instabilities occur approximately every 0.8 P for scenario (b3) and approximately every 1 P for scenario (b2). It is noteworthy that each instability is distinct, even though the system exhibits a repeating pattern resembling a limit cycle.

A similar variability is also observed in PSR J1023+0038 while in the LMXB state, as it notably exhibits a distinctive bimodal behavior in X-ray luminosity (Papitto et al. 2015; Bogdanov et al. 2014; Archibald et al. 2015; De Martino et al. 2013), characterized by low and high modes. The high mode prevails approximately 80% of the time, characterized by a luminosity of approximately 3 × 10³³ erg/s. In contrast, the low mode occurs about 20% of the time, and is characterized by a luminosity around 5 × 10³² erg/s. X-ray flares, detected occasionally, exhibit luminosities reaching up to 5 × 10³⁴ erg/s but are observed in less than 2% of the time. However, distinct from our simulations where instabilities occur on the order of hours, transitions between the low and high states occur rapidly, with a timescale of approximately 10 seconds, and no periodicity seems to emerge (Tendulkar et al. 2014; Bogdanov et al. 2015).

The mechanism governing transitions between the high and low states remains elusive and seems reliant on inhomogeneities in the accretion flow. The discovery of this bimodal flux distribution (Bogdanov et al. 2015; Patruno et al. 2013) triggered discussions regarding possible accretion in the propeller regime with matter reaching the magnetic poles of the NS for the high mode (Papitto & Torres 2015; Campana et al. 2016; Romanova et al. 2005). The low mode on the other hand suggests shock emission from the interaction of the relativistic pulsar wind with the accretion disk (Linares 2014; Zelati et al. 2018), surviving just outside the light cylinder (Ekşİ & Alpar 2005) and advancing back and forth toward the corotation radius. However the corotating radius and light cylinder, which are in our case approximately R_co ≈ 16.7 km and R_LC ≈ 50 km, remains unresolved in our simulation, and we can only suggest a potential link between these two seemingly chaotic variabilities. Moreover, we cannot speculate on a possible IR/optical signature of the behavior observed in our model based on our current considerations.

5.3. Radio pulsar state and transition

Depending on the mechanical luminosity ratio between the companion wind and the pulsar wind, there are cases where the pulsar wind dominates, potentially preventing the formation of the accretion stream from the companion, as in cases (a0), (b0), and (a1) (Fig. 5). Conversely, there are certain cases where the companion wind is strong enough to give rise to an accretion stream, even in cases where its mechanical luminosity is lower than that of the pulsar, as observed in the above cases (panel (7) Fig. 4). The occurrence of these different regimes can be described by the mechanical luminosity ratio between the companion and pulsar winds, as listed in Table 4. The transition between these states is explained by examining the ram pressure in both the accretion stream and the pulsar wind.

Fig. 5. Stationary “PSR” state outcome for case (b0) obtained in case of a dominant pulsar wind. Left is the companion and right is the NS. Purple dashed lines indicate boundary surfaces, while the blue dashed line indicates the isocontour of the effective gravitational potential through the L1 Lagrange point. The color scale gives matter density.

Indeed as explained Sect. 3, wind-wind interaction is usually characterized by the momentum flux ratio β (Eq. 2), determining the shape of the stagnation surface between the two winds but also its precise location r_sh (Eq. 13) (Stevens et al. 1992; Johnstone et al. 2015). However, this ratio describes the entire competition between two isotropic winds when neglecting the effects of the Coriolis force and gravity. In our scenario, where gravity and orbital motion significantly influence wind interaction, it becomes crucial to analyze the ram pressure’s evolution within the accretion stream. Indeed in the context of the RLOF process, the kinetic energy, ram pressure and density within the stream passing through L₁ surpasses that of an isotropic wind. This distinction arises from the condensed and confined flow, channeled through a single point. As a result, there is no smooth transition from a pulsar-dominated to a companion-dominated wind as in the isotropic case.

The radius r_t is the stagnation position from the companion star, where the accretion stream and the pulsar wind are balanced, such that the ram pressures are equal, ρ_c,  w(r_t) v_c,  w²(r_t) = ρ_psr,  w(r_t) v_psr,  w²(r_t). This radius for each cases investigated is shown in Fig. 6 represented by the different intersections between the curve of the pulsar wind and accretion stream ram pressure.

Fig. 6. Evolution of the ram pressure as a function of the distance from the companion star. The colored curves correspond to the accretion stream (for simulations without a pulsar wind), while the black curve corresponds to the pulsar wind.The critical radius rt, crit is shown here with the vertical green line. The vertical gray line indicates the position of L1.

We have established the existence of a critical transition radius r_t, crit such that for r_t < r_t, crit an accretion stream is formed but not for r_t > r_t, crit. In the latter case, the system enters a radio pulsar state characterized by an intra-binary shock wrapping around the companion star. Its stagnation radius lies before the critical radius and an embryonic arm shields the companion from the pulsar wind.

In the isotropic-wind picture the transition between pulsar-dominated and companion dominated occurs at β = 1 with an IBS symetrically wrapping around one or the other component. Here, due to the effect of gravity, the settings required for the companion wind to become dominant and create an accretion stream are far less demanding in energy. We thus propose a more insightful approach by computing the ratio of ram pressure between the pulsar wind and the accretion stream (in simulations without a pulsar wind) at the critical radius r_t, crit. This ratio is defined as:

(30)

By construction, the transition occurs when η = 1. On the other hand whenever η < 1, indicating that the pulsar wind dominates the ram pressure at r_t, crit, the system is in a pulsar state. Conversely, a dominant companion wind at r_t, crit, η > 1, leads to the formation of an accretion stream. The transition is occurring for relatively low value of β compared to analytical works (see Table 4).

We observe a critical radius at r_t, crit ∼ 0.41a_IB, consistent with the instances where the accretion stream does not appear (scenarios (b0) and (a1) Fig. 6). The time required for the stream to fully emerge is prolonged when r_t → r_t, crit⁻.

The distinctive discontinuity in the dominance of the winds underscores the contrast to isotropic systems due to gravity, orbital motion, and the RLOF regime. The existence of this tipping point might provide a plausible explanation for the transitional phenomena observed in spider systems, suggesting that these systems may be approaching a configuration close to this critical limit. However, we note that the instabilities observed in our simulations occur over relatively short timescales (a few hours) compared to the state transitions observed in tMSPs, which occur over several years (Stappers et al. 2014).

6. Orbital variability of the X-ray flux in the pulsar state

In this section we apply the model of X-ray emissions exposed in Sec. 3.2 to the pulsar state in order to produce light curves. We set aside the accretion stream state since no stationary IBS forms and therefore the model does not apply.

6.1. Shock detection

We use the shock detection approach proposed by Lehmann et al. (2016). First, the candidate shock cells are identified using the norm of the relative density gradient,

(31)

and the convergence of velocity using:

(32)

Candidate cells are identified based on a density gradient or a velocity convergence surpassing a user-defined threshold specific to each simulation. This approach enables us to selectively focus on the termination front of the pulsar wind, typically corresponding to the region of the shock that is of particular interest in our study.

6.2. Analogous relativistic bulk velocity

The emission model requires a relativistic bulk velocity within the shock, particularly to account for relativistic beaming. We prescribe that the relativistic specific momentum be proportional to the simulated non-relativistic one, that is,

(33)

where v_R is the relativistic velocity, Γ the corresponding Lorentz factor, v_sim the simulated velocity, and v_sim, max the maximum detected velocity within the shock. In practice the latter corresponds to the asymptotic velocity along the arm. We follow Wadiasingh et al. (2017), Kandel et al. (2019) in order to set the asymptotic Lorentz factor to Γ_max = 1.8.

6.3. X-ray light curves

The shock and light curve in Fig. 7 correspond to case (b0) and are representative of the three cases in the pulsar state (a0, a1, b0). The light curve has a periodic double-peak feature around the companion star inferior conjunction, aligning with expectations for an IBS created due to a dominating pulsar wind. In analytical studies, the directions where the maximum emission occurs typically corresponds to the asymptotic direction of the shock, as explained in Sect. 3.

Fig. 7. Left panel: Detected shocks cells for the pulsar radio state (b0) (also see Fig. 5), with the color scale representing the intensity of the rest-frame synchrotron emission normalized to its maximum value. Arrows indicate the lines of sight corresponding to the vectical lines with the same colors on the right panel. Right panel: Associated X-ray light curve. The orange and blue dashed lines correspond respectively to the star and pulsar inferior conjunctions (at 0.25 and 0.75 of the orbital phase). The vertical magenta and green dashed lines correspond to the angles with the maximum intensity received from the system.

However in our case we note that it does not correspond exactly to this direction. In our simulation the flow crosses the shock surface at a non-negligible angle which mostly explains the observed difference. This allows the regions closer to the apex of the shock, where the rest-frame emission is strongest, to make a larger contribution to the total flux. Consequently, the maximum emission, resulting from Doppler boosting, is achieved closer to the star’s inferior conjunction. Additionally, an eclipse in the emission at precisely 0.25 of the orbital phase further accentuates the double-peak pattern, resulting from obstruction by the companion star itself.

Furthermore the double-peak feature is not perfectly centered at the star’s inferior conjunction (green and magenta arrows on left panel Fig. 7). Indeed, due to the orbital motion, the shock is not entirely symmetrical when surrounding the companion.

In addition, we observe an amplitude asymmetry of about 20% between the peaks: the leading peak has a lower luminosity compared to the trailing one (see right panel of Fig. 7). Indeed, part of the emission contributing to both peaks originates from the apex region of the shock, where a significant fraction is eclipsed by the companion star. As a result of orbital motion, the position of the apex region is shifted and the eclipse is asymmetric. This accounts for about half of the amplitude difference. Another effect is that particles are relatively more accelerated in the trailing arm of the shock, accounting for the remaining amplitude difference. These effects should be all the more pronounced that the orbit is viewed edge-on, and we therefore propose that such asymmetry is an indication of large inclination of the system.

The emission observed near the pulsar inferior conjunction corresponds to the apex of the shock surface where the plasma can be considered at rest and emits isotropically without being eclipsed. This represents approximately 10% of the highest emission.

Comparing our generated light curve with known systems, we can draw parallels with PSR B1957+20 (Huang et al. 2012; Kandel et al. 2021), which has a double-peak pattern around the star’s inferior conjunction. Despite the differences in stellar properties between our simulations and this system, we expect a similar light curve shape due to its edge-on inclination Clark et al. (2023). In Kandel et al. (2021), we can clearly see a pronounced eclipse between the two peaks, where the minimum flux even falls below the flux at pulsar inferior conjunction. This suggests that the shock may be closer to the companion star or more relativistically beamed than in our simulation. Additionally, the asymmetry between the two peaks of the light curve of PSR B1957+20 is consistent with that of Fig. 7, with the leading peak being smaller and the trailing one being larger. This supports the idea that such an asymmetry should be regarded as an indication of the system being nearly edge-on.

PSR J2129–0429 is an example of a RB system where the X-ray light curve is interpreted by an IBS wrapping around the pulsar instead of the companion (Roberts et al. 2015; Al Noori et al. 2018; Kong et al. 2018). Indeed, a double peak is observed around inferior conjunction of the pulsar. No asymmetry is observed between the amplitudes of the two peaks which is consistent with the fact that the pulsar cannot eclipse the shock. For the same reason, the dip between the two peaks can only be due to beaming in this model. A shock surrounding the pulsar is usually observed in RB cases and obtained in studies where the gravity and orbital motion are not factored in. In such studies, an isotropic companion wind with dominating ram pressure results in the bending of the shock toward the pulsar (Romani & Sanchez 2016; Wadiasingh et al. 2018; Kandel et al. 2019). As explained in Sect. 5.1 the configuration of the stellar wind in our simulation is set to create an accretion stream through RLOF. As a result the IBS created arises from the impact with the accretion stream rather than an isotropic companion wind. Instead, a situation similar to the aforementioned studies would require the companion wind speed to be much greater than the system’s gravitational escape velocity, thereby preventing RLOF.

We also note different effects of the momentum flux ratio β on the obtained light curve. The difference between the angles at which the maximum of emission occurs, corresponding to the two peaks, is linked to the opening angle of the shock. This angle diminishes as the ratio β decreases. This phenomenon arises from the fact that a powerful pulsar wind generates a shock that progressively surrounds the companion, resulting in a less distinct double-peak feature. The stagnation radius of the shock created is also strongly influenced by the momentum fluxes ratio. It decreases from r_sh = 0.4a_IB close to the critical radius, to r_sh ∼ 0.3a_IB as the shock will directly hit the surface of the companion for β ≪ 1.

7. Conclusions

We employed the AMRVAC 2.0 code to perform high-precision 2D HD simulations of the interaction between the winds of a stellar companion and a pulsar in a RB system. The goal of our study has been to create a model taking into account the effects of gravity, orbital motion, and heating from the pulsar wind.

To this end, we created a simplified model for both winds. Numerous approximations have been made in prior studies and this work is meant as a step toward a more self-consistent treatment of the complex phenomenology in these rich systems. The pulsar wind is non-relativistic, but settings were selected to reproduce an interaction with a typical MSP of L_sd = 10³⁵ erg/s, achieved by simulating a wind with an equivalent momentum flux. For the companion wind, we explored different values of initial speed and mass flux in order to characterize the process of accretion through RLOF in a spider binary.

This approach allowed us to shed light on the influence of the parameters of the companion wind on the shape and density of the resulting accretion stream. We replicated the two characteristic states observed in transitional systems: the accretion stream state and the radio pulsar state.

In contrast to analytical works, we observed that the transition between these two states is not gradual but instead undergoes a sharp shift at a specific tipping point. The transition radius, r_t, marks when the ram pressure of the accretion stream becomes dominant over the pulsar wind. We defined the critical radius, r_t, crit, such that for r_t < r_t, crit, the stream emerges and reaches the NS. Otherwise the pulsar wind dominates and the system enters a radio pulsar state characterized by an IBS encircling the companion star. We note that the transition occurs for momentum flux ratios β ∼ 0.14. This is much smaller than the value of 1 at which the system becomes dominated by the companion wind in analytical works that neglect gravity. This underscores the importance of incorporating gravity and orbital motion in our model.

Furthermore, we witnessed the instabilities induced by the pulsar wind while in the accretion stream state. In proximity to the tipping point, we identified an unstable behavior characterized by the intermittent creation and disruption of the accretion stream state. Although these alternations of pulsar and radio phases occur on timescales of a few orbital periods in our simulations (i.e., on a much shorter timescale than the few years observed for tMSPs), we suggest that this hints at a plausible explanation for the transitional phenomena observed. Indeed, these systems could be in a state close to a similar tipping point.

The IBS was chosen as the observable for comparison with known observations. We constructed the corresponding X-ray light curves modeling the SR emission occurring in the shock for the radio pulsar state. This state corresponds to a BW scenario where the pulsar wind dominates. The RB case, where the shock is surrounding the pulsar, was not stable with the physical parameters used in this work. The orbital motion shifts the position of the peaks with respect to the inferior conjunction of the companion. Notably, the leading peak is less intense than the trailing peak due to both the orbital motion and the occultation of the shock front apex by the companion. We suggest that such an asymmetry could be indicative of high inclination. Indeed, a similar asymmetry is observed in the light curve of the system PSR B1957+20, which is almost edge-on.

Our 2D approach confines us to the orbital plane. A 3D model would allow us to verify if the accretion flow remains within the orbital plane or if the expulsion of matter by the pulsar wind occurs in a preferred direction. Additionally, 3D modeling would enable examination of the impact of the system’s inclination on the observables. This would allow for observations of the shock morphology due to the interaction with the accretion flow in the direction perpendicular to the orbital plane.

The follow-up to this study will be to consider a ray tracing algorithm for the heating model. The ideas to consider a relativistic wind for the pulsar and doing MHD simulation to take into account the magnetic field impact on the transitions are also really important to constrain our system as much as possible.

Acknowledgments

We would like to express our gratitude to the referee for their insightful and pertinent questions, which significantly contributed to improving the quality of this paper.

References

All Tables

All Figures

	Fig. 1. Sketch of the grid setup used in simulations. The companion star is at the center of the cylindrical grid, and the pulsar at a distance aIB along X axis. Wind injection surfaces are shown in red, gravitational potential of Lagrange 1 point is shown in blue. Top left arrow indicates sens of rotation.
In the text

Fig. 2. Comparison of the accretion stream in cases (b1) and (b3) without pulsar wind, that is for two different companion wind velocities and equal companion mass loss. Left is the companion and right is the NS. Purple dashed lines indicate boundary surfaces, while the blue dashed line indicates the isocontour of the effective gravitational potential through the L1 Lagrange point. The black arrow shows the direction of orbital rotation. The color scale gives matter density.

In the text

Fig. 3. Variability of the accretion stream regime with and without pulsar wind. The panels show, from top to bottom, the temporal variation of mean density (in g cm−3) within a ring, extending from the internal boundary of the pulsar wind to the circularization radius Rcirc. Left column correspond to simulations without a pulsar wind (see Sect. 5.1), while right column to simulation with both winds. First line corresponds to scenario (b1), second line to scenario (b3).

In the text

Fig. 4. Visualization of the interaction between the accretion stream and pulsar wind for the case (b1). Left is the companion and right is the NS. Purple dashed lines indicate boundary surfaces, while the blue dashed line indicates the isocontour of the effective gravitational potential through the L1 Lagrange point. The color scale gives matter density. Green text top left indicates the time in unit of orbital period, P.

In the text

	Fig. 5. Stationary “PSR” state outcome for case (b0) obtained in case of a dominant pulsar wind. Left is the companion and right is the NS. Purple dashed lines indicate boundary surfaces, while the blue dashed line indicates the isocontour of the effective gravitational potential through the L1 Lagrange point. The color scale gives matter density.
In the text

	Fig. 6. Evolution of the ram pressure as a function of the distance from the companion star. The colored curves correspond to the accretion stream (for simulations without a pulsar wind), while the black curve corresponds to the pulsar wind.The critical radius rt, crit is shown here with the vertical green line. The vertical gray line indicates the position of L1.
In the text

Fig. 7. Left panel: Detected shocks cells for the pulsar radio state (b0) (also see Fig. 5), with the color scale representing the intensity of the rest-frame synchrotron emission normalized to its maximum value. Arrows indicate the lines of sight corresponding to the vectical lines with the same colors on the right panel. Right panel: Associated X-ray light curve. The orange and blue dashed lines correspond respectively to the star and pulsar inferior conjunctions (at 0.25 and 0.75 of the orbital phase). The vertical magenta and green dashed lines correspond to the angles with the maximum intensity received from the system.

In the text