© The Authors 2025

1. Introduction

LS 5039 is a binary system hosting an O6.5 V star and a compact object of unclear origin (Casares et al. 2005), and a powerful gamma-ray source with a nonthermal radiation luminosity, L_NT ∼ 10³⁶ erg s⁻¹, peaking around ∼10 − 100 MeV (e.g., Collmar & Zhang 2014). As was shown in Casares et al. (2005), a black hole in the system is likely only if the latter is pseudo-synchronized, although hints of P ≈ 9 s period pulsations found in X-rays would favor a pulsar if confirmed (Yoneda et al. 2020). Despite evidence of X-ray pulsations being under debate (Kargaltsev et al. 2023; Makishima et al. 2023), they would strongly indicate the presence of a slowly rotating, magnetar-like neutron star (NS), make the standard scenario of a pulsar in the ejector regime less likely due to a low wind power, and favor a magnetospheric origin for the energy fueling the nonthermal emission (Yoneda et al. 2020). It is worth noting that LS I +61 303, a similar high-mass binary and gamma-ray emitter, has been found to host a relatively slowly rotating NS that may also present magnetar-like behavior (e.g., Weng et al. 2022, and references therein). Given all this, a model for energizing the nonthermal emitter with the energy stored in the NS magnetosphere may be required in these sources.

Yoneda et al. (2020) explored analytically the possibility of feeding the nonthermal emission in LS 5039 with energy from the magnetosphere, released due to its interaction with the stellar wind, but the specifics of this interaction were not provided. Moreover, the presence of a magnetar in the system should be reconciled with the absence of evidence for a very young supernova remnant (SNR) near LS 5039 (Moldón et al. 2012). Magnetar flares have been discussed in the literature (see, e.g., Beloborodov 2020; Khangulyan et al. 2022, and references therein) as the sources of fast radio bursts (FRBs; Lorimer et al. 2007), particularly in the context of flares of a magnetar in a high-mass binary (Bozzo et al. 2008) when there is periodicity (see, e.g., Lyutikov et al. 2020b; Barkov & Popov 2022), so LS 5039 hosting a magnetar would fit in this framework. Thus, in this work we propose a mechanism for LS 5039 based on recurrent periods of magnetar-like flaring activity. These periods start with a weak NS (or pulsar) wind that cannot balance the ram pressure of the wind from the companion star, and the shocked pulsar wind can reach or even enter the light cylinder of the NS, directly touching its magnetosphere (as is proposed in Yoneda et al. 2020). This can trigger magnetic perturbations in the magnetosphere that can quickly grow, and ultimately lead to magnetic reconnection and the release of large amounts of energy into the otherwise weak pulsar wind beyond the light cylinder (Barkov et al. 2022; Sharma et al. 2023). The injection of magnetospheric energy lasts as long as the pulsar wind termination shock touches the magnetosphere, and during this period the shocked pulsar wind mediates the energy transfer between the latter (i.e., the energy source) and the region resulting from the star-pulsar wind interaction (i.e., where -or downstream of which- the nonthermal emitter is located). Eventually, though, the enhanced pulsar wind drives the whole shocked flow structure beyond the light cylinder, stopping the injection of energy. This severely reduces the pulsar wind ram pressure and the shocked pulsar wind quickly expands inward until it again touches the magnetosphere, moment at which the cycle can repeat. This process leads to recurrent periods of power-enhanced pulsar wind that feed the overall nonthermal emission in the system. To show how the shocked structure evolves, we carried out relativistic hydrodynamics (RHD) simulations that mimic the effect of the magnetic active periods on the shocked winds, including the role of orbital motion. We also explored, using analytical tools, the processes near the NS relevant for the nonthermal emission.

Regular pulsar winds (i.e., not enhanced by magnetar activity) can already be significantly magnetized and nonspherically symmetric (e.g., Bogovalov et al. 2019). Moreover, the part of the wind enhanced by magnetar activity can also cover a relatively narrow solid angle of a few steradians, which increases the wind anisotropy (Sharma et al. 2023). In addition, three-dimensional (3D) relativistic magnetohydrodynamics (RMHD) simulations (Porth et al. 2014; Barkov et al. 2024) showed that the shocked pulsar wind structure is very sensitive to changes in the environment. Thus, the magnetic field is indeed a major factor in the problem, and should eventually be included in simulations of the discussed scenario. However, the 3D RMHD solution turned out to be on average similar to that for a spherical RHD wind, which also shows strong instabilities (see sect. 3 in Barkov et al. 2024). Furthermore, the simulations presented here took about 100 000 CPU hours; since an RMHD simulation requires at least a 10 times higher spatial resolution, the computational time should be ∼10⁹ CPU hours, which is extremely hard to achieve. For all this, at this stage RHD simulations were employed as they provide qualitatively relevant results, while keeping the computing costs at a reasonable level and allowing for an exploratory study.

The article is structured as follows: In Sect. 2 the simulations are described, and their results are shown in Sect. 3. Then, in Sect. 4, a discussion is carried out about the regime of interaction between the NS and the surrounding medium, the origin of the nonthermal energy and the age of the system, and the nonthermal emitter.

2. Simulation setup

The RHD simulations were performed in 2D and 3D geometry in Cartesian (2D and 3D) and cylindrical (axisymmetric 2D) coordinates using the PLUTO code¹ (Mignone et al. 2007). Spatial linear interpolation, a third-order Runge-Kutta approximation in time, and an HLLC Riemann solver were adopted (Li 2005). PLUTO is a modular Godunov-type code entirely written in C and intended mainly for astrophysical applications and high-Mach-number flows in multiple spatial dimensions. The simulations were performed on the CFCA XC50 cluster of the National Astronomical Observatory of Japan (NAOJ) and RIKEN HOKUSAI Bigwaterfall. The flow was approximated as an ideal, relativistic adiabatic gas, with one particle species and a polytropic index of 4/3. The size of the Cartesian domain is x and y ∈ [ − 32, 32] in 2D and 3D, and z ∈ [ − 20, 20] in 3D. The resolution is nonuniform in the entire computational domain, with a total number of cells: N_X = N_Y = 780 (2D and 3D), and N_Z = 650 (3D; see details in Table 1). The size of the axisymmetric domain is r ∈ [0, 0.7], and z ∈ [0.7, 2.0], with a total number of cells: N_r = 260, and N_Z = 360. The higher-resolution central part features N_r = 130, and N_Z = 260, with size r ∈ [0, 0.25], and z ∈ [0.7, 1.2].

In this paper, the unit of length a = 2.2 × 10¹² cm corresponds to the semimajor axis for an orbital period of T_orb ≈ 3.9 days and a binary total mass of ≈30 M_⊙ (Casares et al. 2005). This work was intended as a first exploration of the phenomenon under study. Thus, the results of the Cartesian 2D simulations were meant to be illustrative of the large-scale evolution (Bosch-Ramon et al. 2012), but the difficulty in capturing in detail the magnetospheric scale (∼10⁻² times smaller than a; see Sect. 4) precluded in that case a detailed characterization of the shocked pulsar wind evolution right between the star and the NS. The latter limitation also affects the Cartesian 3D simulations. Nevertheless, the tendencies in the flow evolution on the smallest captured scales hinted at the role of the shocked pulsar wind when enhanced by strong magnetospheric activity. In the Cartesian simulations, a corresponds to ≈64 cells. On the other hand, the axisymmetric 2D simulations allowed us to study the shocked two-wind region closest to the magnetosphere, with a corresponding now to ≈520 cells, but without including orbital motion, the latter playing a minor role on those scales. Finally, the equation of state adopted, albeit simple, is already enough to derive meaningful conclusions from the simulations (see the explorations carried out in Bosch-Ramon et al. 2012, 2015; Barkov et al. 2024).

2.1. Setup of the winds

The wind of the massive star is assumed to be cold and move radially with constant velocity v_w = 2400 km/s and Mach number M_w = 7. The pulsar wind is also cold and moves radially with Lorentz factor Γ_pw = 1.9 and Mach number M_pw = 17, and the star-to-pulsar wind momentum rate ratio is η = 0.05 in the weak wind phase in the Cartesian simulations and η = 4 × 10⁻⁴ in the axisymmetric simulations, where

$$\begin{array}{c}\hfill \eta =\frac{{\dot{E}}_{\mathit{sd}}}{c\dot{M}{v}_w},\end{array}$$(1)

with Ė_sd being the spin-down power of the pulsar and Ṁ the massive star mass-loss rate. The winds were modeled similarly to those in Bosch-Ramon et al. (2015). During flaring periods, in the Cartesian case η grows such that its orbit-averaged value reaches 0.16, but in one case it reaches 0.39 (in the axisymmetric case only one flare is presented; see the next section).

2.2. Setup of the magnetar flaring periods

The flaring periods were approximated via a function representing recurrent additional injections of power, Ė_pw ∝ Ė_fl(1+δρ) ∝ ρ_pw(1+δρ), which has the following time profile:

$$\begin{array}{c}\hfill \delta \rho =A_{\mathit{fl}}\frac{5}{\pi }\frac{e^{-\overline{y}}}{{\mathrm{\Delta }}^{1/2}},\mathrm{with}\overline{y}\equiv 15.917\frac{1-cos\overline{x}}{\mathrm{\Delta }},\end{array}$$(2)

where $\overline{x}=2\pi t{N}_{\mathit{fl}}/T_{\mathrm{orb}}$ is a dimensionless time variable normalized to the orbital period, T_orb, ρ_pw is the weak pulsar wind density, N_fl the number of flaring periods per orbit, Δ the relative width of the flaring period, and A_fl its amplitude. The last two terms quantify how much energy per orbit is released by the flaring periods compared to the one from the weak pulsar wind alone, T_orb < Ė_fl > and T_orb < Ė_sd >, respectively². The time profile of the flaring periods was chosen to be sharp, while avoiding numerical artifacts in the simulations. A profile with Δ = 0.1 is presented in Fig. 1 for better visualization, but in the simulations a narrower profile with Δ = 0.05 was used. In the Cartesian simulations, the energy injection averaged over the orbit including both the flaring and weak wind periods releases 3.3 (≡ < η >  = 0.16) and 7.8 (≡ < η >  = 0.39) times more energy than the weak pulsar wind alone (η = 0.05). The maximum enhancement of injected power within the flaring periods is higher: 72 (max(η_fl) = 3.6) and 214 (max(η_fl) = 10.7) times larger than in the weak wind phase. The parameter values used are presented in Table 2. In the axisymmetric simulations, we were exploring a different spatial and temporal scale, as was allowed by the more focused setup, and the flare was made significantly shorter to mimic the response of the shocked structure on the scales of the magnetospheric radius. Thus we took Δ = 0.003, and the adopted maximum enhancement of injected power during the flaring period was higher, by a factor of 870, which corresponds to max(η_fl) = 0.35 from η = 4 × 10⁻⁴.

Fig. 1. Example of the time profile of a flaring period; δρ is the density on top of the regular pulsar wind density during flaring periods, x = 2πtNfl/Torb the dimensionless time (see Sect. 2.1 for details), and Δ = 0.1 the width of the pulse.

3. Results

We carried out simulations for four 2D models and one 3D model in Cartesian coordinates including orbital motion, and one simulation for one 2D axisymmetric model. In the 2D Cartesian models, we varied the number of flaring periods during the orbit (N_fl = 10, 30, 90) for an orbit-averaged power of the flaring component alone of < Ė_fl > = 2.3 < Ė_sd >, and for N_fl = 30, we varied the flaring period power ×3 (i.e., < Ė_fl > = 6.8 < Ė_sd >). For higher flaring period rates, the interaction structure becomes smoother. In the case 2Dn90e10, with N_fl = 90, the leading edge of the shock is the most uniform (see Fig. 2, left top panel). The result of multiple flaring periods can be seen most clearly in the radial shock fronts that form on the trailing side of the interaction structure. Non-flare solutions are also shown, nfeta0.05 and nfeta0.5, which illustrate two cases with a strong difference in η, 0.05 (weak wind) and 0.5 (enhanced wind), but without the effect of the flaring periods (left and right bottom panels in Fig. 2). The main differences between the flare and non-flare solutions are: 1) the free pulsar wind zone (FPWZ), which is significantly smaller in the cases with flares, and 2) the spiral tail is significantly more loaded by the stellar wind in those cases. Point 2 is best illustrated by comparing the density of the shocked pulsar wind in the cases with flares and those without. Due to heavy mass-loading, the shocked pulsar wind in the spiral arm is slower in the cases with flares.

Fig. 2. Colored density maps on the orbital plane, with colored arrows representing the modulus of the 3-velocity vector for the 2D models. The flaring period rates are Nf = 90 (top left, 2Dn90e10), 30 (top right, 2Dn30e10), and 10 (middle left, 2Dn10e10) flaring periods per orbit, keeping the total energy budget of the system constant. In the model 2Dn30e30, we increased the power of the flaring period by a factor of 3, while keeping Nf = 30 (middle right panel). The non-flaring cases are presented at the bottom left (nfeta0.05) and right (nfeta0.5).

For the cases 2Dn30e10 and 2Dn10e10, the individual flaring periods are more powerful and more separated in time compared to 2Dn90e10 (top right, middle left, and middle right panels in Fig. 2, respectively). As a result, the impacts of individual flaring periods become clearly visible in the density structure of both the leading and the trailing sides of the interaction structure. Increasing the orbit-averaged flaring period power, as seen comparing the 2Dn30e10 and the 2Dn30e30 models (top right and middle right panels in Fig. 2, respectively), leads to a larger two-wind interaction structure, while the FPWZ size decreases when increasing this power.

In the 3D simulations (model 3Dn30e10; Fig. 3), we switched on the flaring periods for short periods of time in a series of five flaring periods to mimic the flaring period impact within a fraction of the orbit (they occur for one sixth of the orbit), as the higher computational demands in 3D prevented us from probing the effect of the recurrent flaring periods for as long a time as in 2D. As we know from previous work (see, e.g., Bosch-Ramon et al. 2015), the results in 3D are more turbulent than in 2D in general, so the flares enhance turbulence further. The impact of the flaring periods on the overall two-wind interaction structure can be better seen in Fig. 4, although the limitations of the simulations prevented us from reproducing in detail the evolution of the pulsar wind termination shock; this is exemplified by a zoom-in of the interaction structure for the 3D simulations shown in Fig. 5, where the pulsar wind termination shock does not actually reach the magnetospheric scales. On the other hand, in the 2D axisymmetric simulations presented in Fig. 6, which do not account for orbital motion but probe more realistically the spatial and temporal scales close to the NS, the cycle of weak-enhanced-weak pulsar wind was better captured. The series of events that the simulations show can be described as follows. Initially, the magnetospheric flares are not active, and the weak pulsar wind termination shock can reach close to the NS, mimicking the situation when this shock is about to touch the magnetosphere (Fig. 6, top left³). When this happens, the flaring period resumes (Fig. 6, top right). Due to its high density, the shocked stellar wind reacts slowly to the power-enhanced pulsar wind and the shocked pulsar wind stays in contact with the magnetosphere, still exciting the magnetospheric activity for some time (Fig. 6, middle left). Eventually, as the shocked stellar wind is pushed outward, the pulsar wind termination shock moves away from the magnetosphere, and the extra energy injection stops (Fig. 6, middle right). At that point, the weakening of the pulsar wind allows the shocked pulsar wind to quickly propagate inward, toward the magnetosphere (Fig. 6, bottom panels), the system again becoming ready to restart the cycle. Using 2D axisymmetric simulations on small scales more systematically to better inform simulations on large scales (which here deal in a simplified manner with those small scales) is to be tackled in future studies.

Fig. 3. Colored density maps with colored arrows representing the modulus of the 3-velocity vector for the 3D model (3Dn30e10) for various cuts: XY (top; orbital plane), XZ (middle), and ZY (bottom).

Fig. 4. Colored density maps with colored arrows representing the modulus of the 3-velocity vector for the 3D model (3Dn30e10) for cuts on the orbital plane, at orbital phases ϕ = 0.5435 (top left), ϕ = 0.575 (top right), ϕ = 0.5812 (middle left), ϕ = 0.6033 (middle right), ϕ = 0.6467 (bottom left), and ϕ = 0.6576 (bottom right).

Fig. 5. Zoom-in of three snapshots of the cycle followed by the interaction structure: weak pulsar wind (left), enhanced pulsar wind (middle), and the end of energy injection (right).

Fig. 6. Colored density maps obtained from an axisymmetric simulation showing qualitatively the evolution of a flare, from top left to bottom right: 0 (weak wind), 48 (flare starts), 88 (flare continuation), 141 (end of flare), 166 (pulsar wind approaches magnetosphere), and 186 s (flare about to resume) from the beginning of the flare. The grid size is ≈0.65 a; the magnetosphere radius would be ≈0.02 a in this case.

4. Summary and discussion

In the previous section, we described how magnetospheric energy can be recurrently injected in the two-wind interaction region, and we illustrated this process with RHD simulations. Some caveats were raised in Sect. 1, as the pulsar wind enhanced by magnetar flares is expected to be magnetized and not spherically symmetric (e.g., Sharma et al. 2023). However, these differences from a spherical RHD wind are alleviated by the fact that other common factors (stellar wind, the nonsteady nature of the shocked flows, orbital motion, etc.) also play a strong role in the evolution of the interaction structure. Moreover, strong relativistic waves are expected between the magnetosphere and the stellar wind during magnetar flares, which can lead to a more spherical shocked pulsar wind, and to a reduction in the wind magnetization. Assessing all this requires detailed future analysis, and better connecting processes taking place on different spatial and temporal scales. Putting these caveats aside, in what follows some aspects of the proposed scenario are discussed; namely, the regime in which the stellar wind interacts with the NS surroundings, the nonthermal energy origin and the age of the source, the latter being relevant as no recent SNR has been found near LS 5039. The article finishes with a few remarks regarding the nonthermal emitter. In all the expressions below, we adopted the convention A_ξ = A/10^ξ, with A being in cgs units but for Ṁ, which is in units of solar masses per year.

4.1. NS-medium interaction regime

A NS interacting with a stellar wind can have two regimes in LS 5039: (1) ejector and (2) georotator (or propeller; see below; Illarionov & Syunyayev 1975; Lipunov 1992), which are discussed below. The state of the system is characterized by several radii: The distance at which pulsar and stellar winds would be in ram pressure balance⁴:

$$\begin{array}{c}\hfill R_s=a\frac{{\eta }^{1/2}}{1+{\eta }^{1/2}}\approx 3.5\times {10}^{10}{a}_{12.3}{\eta }_{-3.5}^{1/2}\mathrm{cm}(\mathrm{assuming}\eta \ll 1);\end{array}$$(3)

the light cylinder radius,

$$\begin{array}{c}\hfill R_{\mathit{LC}}=cP/2\pi \approx 3.8\times {10}^{10}{P}_{0.9}\mathrm{cm};\end{array}$$(4)

the Alfvenic radius at the balance location between the stellar wind ram pressure and the dipolar magnetic pressure of the NS magnetosphere:

$$\begin{array}{c}\hfill R_A={\left(\frac{a^2{\mu }^2}{2\dot{M}{v}_w}\right)}^{1/6}\approx 2.7\times {10}^{10}\frac{a_{12.3}^{1/3}{\mu }_{33}^{1/3}}{{\dot{M}}_{-6.4}^{1/6}{v}_{w,8.3}^{1/6}}\mathrm{cm},\end{array}$$(5)

where μ = B_dpR_NS³ is the dipolar NS magnetic moment, and B_dp the dipolar magnetic field at the NS surface; and the gravitational capture radius:

$$\begin{array}{c}\hfill R_g=\frac{2G{M}_{\mathit{NS}}}{v_w^2}\approx 9.3\times {10}^9{v}_{w,8.3}^{-2}\mathrm{cm},\end{array}$$(6)

for M_NS = 1.4 M_⊙ and a wind speed normalization close to, but a bit lower than, its velocity at infinity of ≈2.4 × 10⁸ cm s⁻¹ (Casares et al. 2005), as we took into account the compactness of the system. For B_dp ≳ 4 × 10¹³ G (and taking R_NS = 10⁶ cm), one has R_g < R_A, which is the reason we focus on the georotator regime hereafter. Nevertheless, the propeller regime could still be compatible with the triggering of magnetic flares. We did not account for a multipolar magnetic field (see below) when deriving R_A as such a field drops much faster than B_dp.

The induced magnetospheric flares could be possible either in a marginal ejector regime, or in the georotator regime. Let us discuss first the ejector regime.

4.1.1. Ejector regime

The ejector regime takes place if R_s > R_LC. From the dipole formula, the spin-down power becomes

$$\begin{array}{c}\hfill {\dot{E}}_{\mathit{sd}}=\frac{{\mu }^2{\mathrm{\Omega }}^4}{c^3}(1+sin{(\alpha )}^2)\sim 1.5\times {10}^{34}\frac{{\mu }_{33}^2}{P_{0.9}^4}{\mathrm{erg}\mathrm{s}}^{-1},\end{array}$$(7)

where Ω = 2π/P. For simplicity, in Eq. (7) we neglected the dependence on the angle between between the magnetic momentum and the rotation axis, which is ∼(1 + sin(α)²) (see, for numerical force-free simulations, e.g., Bogovalov 1999; Spitkovsky 2006, and for the analytical expression (2/3)sin(α)², e.g., the review by Beskin 2010), and we simplified it to ∼1, but it does not change our conclusions significantly. The marginal ejector regime would correspond then to R_s temporarily reaching R_LC. Taking R_s = R_LC, one can derive a minimum dipolar magnetic field for the ejector regime:

$$\begin{array}{c}\hfill B_{\mathit{cr}}\approx \frac{c^3}{8{\pi }^3}\frac{P^3\sqrt{\dot{M}{v}_w}}{R_{\mathit{NS}}^{3}a}\approx 2\times {10}^{15}\frac{P_{0.9}^3{\dot{M}}_{-6.4}^{1/2}{v}_{w,8.3}^{1/2}}{R_{NS,6}^3{a}_{12.3}}\mathrm{G},\end{array}$$(8)

where we assumed η ≪ 1. In the marginal ejector regime, the magnetospheric flares would occur when B_dp ∼ B_cr. If B_dp = B_cr, the dipole spin-down time would be

$$\begin{array}{c}\hfill T_{\mathit{sd}}\approx \frac{1}{24}\frac{I{P}^3{c}^3}{{\pi }^2{\mu }^2}\approx \frac{8{\pi }^4}{3c^3}\frac{a^{2}I}{P^3\dot{M}{v}_w}\approx 500\frac{I_{45}{a}_{12.3}^2}{P_{0.9}^3{\dot{M}}_{-6.4}{v}_{w,8.3}}\mathrm{yr},\end{array}$$(9)

where I ≈ 10⁴⁵ g cm² is the NS moment of inertia.

The T_sd value obtained assuming B_dp = B_cr is similar to the τ ∼ 500 yr inferred from the X-ray data by Yoneda et al. (2020) and Makishima et al. (2023) for the pulsar spin-down. Under the same assumption, the spin-down power, Ė_sd, is ∼7 × 10³⁴ erg s⁻¹, which is also close to the value inferred by Yoneda et al. (2020), and significantly below L_NT ∼ 10³⁶ erg s⁻¹ in LS 5039. This spin-down power would correspond to the weak pulsar wind phase in the simulations, that is, between flaring periods. We recall that in the Cartesian simulations including orbital motion, η = 0.05 was adopted due to resolution limitations, making the simulation results qualitative, and in reality η would be much smaller, as in the axisymmetric simulation. For Ė_sd ∼ 7 × 10³⁴ erg s⁻¹, Ṁ ∼ 4 × 10⁻⁷ M_⊙ yr⁻¹, and v_w ∼ 2 × 10⁸ cm s⁻¹, one obtains η ∼ 5 × 10⁻⁴.

4.1.2. Georotator regime

For B_dp < B_cr, which implies that R_s < R_LC, the system is in the georotator regime and the stellar wind directly interacts with the NS magnetosphere, and the maximum spin-down power associated with this interaction can be derived as follows. First, the mass rate of the stellar wind intercepted by the NS magnetosphere is obtained:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathit{NS}}\approx \dot{M}\frac{R_A^2}{4a^2}={\left(\frac{{\mu }^2{\dot{M}}^2}{2^7{v}_w{a}^4}\right)}^{1/3},\end{array}$$(10)

which, following (Lyutikov et al. 2020a), allows for the maximum power dissipated on the magnetosphere boundary to be derived, as the incoming flow is forced to corotate with the NS magnetic field⁵:

$$\begin{array}{c}\hfill {\dot{E}}_{\mathit{Gr}}\approx \frac{{\pi }^2}{2^{2/3}{P}^2}{\left(\frac{{\mu }^4\dot{M}}{v_w^2{a}^2}\right)}^{1/3}\approx 2.5\times {10}^{34}\frac{{\mu }_{32}^{4/3}{\dot{M}}_{-6.4}^{1/3}}{P_{0.9}^2{v}_{w,8.3}^{2/3}{a}_{12.3}^{2/3}}{\mathrm{erg}\mathrm{s}}^{-1}.\end{array}$$(11)

This value is similar to the Ė_sd value obtained above, so the small amount of energy released in this regime when no flares happen may also be considered equivalent to the weak pulsar wind phase in the simulations.

From the maximum power estimate given in Eq. (11), in the georotator regime the spin-down time is

$$\begin{array}{c}\hfill T_{sd,Gr}\gtrsim 2^{5/3}\frac{I{v}_w^{2/3}{a}^{2/3}}{{\mu }^{4/3}{\dot{M}}^{1/3}}\approx 400\frac{I_{45}{v}_{w,8.3}^{2/3}{a}_{12.3}^{2/3}}{{\mu }_{32}^{4/3}{\dot{M}}_{-6.4}^{1/3}}\mathrm{yr}.\end{array}$$(12)

Taking T_sd, Gr = τ ∼ 500 yr from Yoneda et al. (2020), one can estimate the minimum dipolar NS magnetic moment as (τ units are seconds: 500 yr ≈10^10.2 s)

$$\begin{array}{c}\hfill {\mu }_{min}\approx 8\times {10}^{31}\frac{I_{45}^{3/4}{v}_{w,8.3}^{1/2}{a}_{12.3}^{1/2}}{{\tau }_{10.2}^{3/4}{\dot{M}}_{-6.4}^{1/4}}{\mathrm{G}\mathrm{cm}}^3,\end{array}$$(13)

so the minimum dipolar NS magnetic field is

$$\begin{array}{c}\hfill B_{dp,min}\approx 8\times {10}^{13}\frac{I_{45}^{3/4}{v}_{w,8.3}^{1/2}{a}_{12.3}^{1/2}}{R_{NS,6}^3{\tau }_{10.2}^{3/4}{\dot{M}}_{-6.4}^{1/4}}\mathrm{G}.\end{array}$$(14)

Adopting B_dp = B_cr, R_g is ∼0.3 R_A for the source specific parameters, whereas taking B_dp = B_dp, min, R_g becomes ∼0.8R_A, still marginally in the georotator regime. Importantly, the lower value of B_dp in the latter case can lead to a much longer ejector phase; for μ = μ_min,

$$\begin{array}{c}\hfill T_{\mathit{sd}}({\mu }_{min})\approx 4\times {10}^3\frac{P_{0.3}^3{\tau }_{10.2}^{3/2}{\dot{M}}_{-6.4}^{1/2}}{v_{w,8.3}{a}_{12.3}{I}_{45}^{1/2}}\mathrm{yr},\end{array}$$(15)

where P has been normalized now to 2 s to account for pulsar rotation braking while still fulfilling R_s > R_LC (ejector). If one sets B_cr = B_dp, min, one can obtain the longest pulsar period for which this condition is still fulfilled:

$$\begin{array}{c}\hfill P_{max}\approx 2.8\frac{I_{45}^{1/4}{a}_{12.3}^{1/2}}{{\tau }_{10.2}^{1/4}{\dot{M}}_{-6.4}^{1/4}}\mathrm{s}.\end{array}$$(16)

Substituting it into Eq. (15) yields the duration of the ejector phase corresponding to P_max:

$$\begin{array}{c}\hfill T_{sd,max}\approx 1.1\times {10}^4\frac{{\tau }_{10.2}^{3/4}{a}_{12.3}^{1/2}{I}_{45}^{1/4}}{v_{w,8.3}{\dot{M}}_{-6.4}^{1/4}}\mathrm{yr}.\end{array}$$(17)

An important conclusion is reached here: if the source were in the georotator regime between flares, and not marginally in the ejector regime as first discussed, the previous ejector phase could have lasted for ∼10⁴ yr due to the much lower B_dp and an initial pulsar period of ∼P_max.

4.2. Nonthermal energy origin

The magnetic energy of the NS can be estimated as

$$\begin{array}{c}\hfill E_{\mathrm{mag}}\approx \chi R_{\mathit{NS}}^3{B}_{\mathit{NS}}^2\approx 1.7\times {10}^{47}{R}_{NS,6}^3{B}_{NS,15}^2\mathrm{erg},\end{array}$$(18)

with the constant χ depending on the field topology (1/6 for a dipole field, as adopted by us, and slightly smaller for a higher multipole field). Here, B_NS includes magnetic field components beyond the dipolar one, as magnetar-like activity requires of reconnection events that occur when the magnetic field has high multipoles. This multipolar NS magnetic field can be characterized through an amplification factor, x_B, with respect to B_dp, which can be significantly larger than 1: B_NS = x_BB_dp ≫ B_dp. The multipolar magnetic field can be generated through the dynamo mechanism; for instance, when the compact object formed (see, e.g., Landau & Lifshitz 1980; Thompson & Duncan 1993; Mösta et al. 2015; Raynaud et al. 2020; Sharma et al. 2023; Barrère et al. 2025; Igoshev et al. 2025, and references there in).

The dissipation of the NS magnetic field takes place on a timescale of

$$\begin{array}{c}\hfill {\tau }_{\mathrm{mag}}\lesssim \frac{E_{\mathrm{mag}}}{L_{\mathit{NT}}}\sim 5\times {10}^3\frac{R_{NS,6}^3{B}_{NS,15}^2}{L_{NT,36}}\mathrm{yr},\end{array}$$(19)

where we took into account that not all the energy released in the magnetic flares goes into nonthermal emission, as some of this energy goes to make work on the environment or heat the NS, or to nonthermal particles that do not radiate efficiently. In the case of B_NS ∼ B_dp, min and the georotator regime, the magnetic dissipation timescale would be ≲30 yr, which is completely implausible. To sustain L_NT for ∼10⁴ yr for instance, at least B_NS ∼ 1.4 × 10¹⁵ G is needed, meaning that x_B ≳ 20.

A marginal ejector case with B_NS ∼ B_dp ∼ B_cr, and thus x_B ∼ 1, could sustain L_NT for a few centuries until the georotator regime were reached, and the subsequent georotator regime, albeit powerful, would evolve very fast (in ∼10 yr, as B_dp should still be ∼B_cr to feed L_NT, see Eq. 12), all of this implying a very young SNR. Therefore, a present georotator regime with a large χ_B (i.e., with a relatively small B_dp) would be favored for LS 5039. In addition, in this case, even if the present georotator phase is just a few centuries old, and the NS slows down its rotation significantly during that timescale, this regime may last much longer with substantial nonthermal activity, say for ∼10⁴ yr. Thus, such a type of source can be an active accelerator for a relatively long time. A sketch of the scenario based on a present georotator regime for LS 5039 is shown in Fig. 7.

Fig. 7. Sketch of the georotator scenario evolution, showing the pulsar angular velocity evolution with time with the ejector and the georotator (present) phases and a dashed green line marking the transition.

4.3. Lack of SNR evidence

As has been shown, in the georotator regime, the total time elapsed since the NS was formed can be ≫τ ∼ 500 yr, the age from Yoneda et al. (2020). There is first the period in which the pulsar was still in the ejector regime, which as explained may have lasted for ∼10⁴ yr depending on the initial period and magnetic field. Then, in the georotator regime, the system may have evolved for some time until reaching the present stage, with a current decaying time of several centuries. Nevertheless, regardless of the georotator decaying time, since the multipolar magnetic field can be strong even when the NS rotation slows down, the actual magnetic field dissipation phase can last much longer (e.g., τ_mag ∼ 10⁴ yr). From all this, the SNR may have occurred long ago and become by now already inconspicuous and difficult to detect.

Another, less likely, possibility for the absence of observational evidence of a SNR associated with LS 5039 could be a very small ejecta mass (e.g., up to 0.1 M_⊙) in the case of a white dwarf (WD) turning into a NS (see eg. Woosley & Baron 1992; Kuroda et al. 2025). Such a scenario is unlikely as the rate of such a kind of SN explosion, particularly in high-mass binaries, is very low. Such a scenario may have occurred if the WD progenitor had transferred most of its mass to the secondary star, nowadays an O6.5 V star (Casares et al. 2005), preventing its initial direct collapse into a NS. In this quite improbable scenario, the NS might be quite young but no SNR would be visible.

4.4. The nonthermal emitter

The proposed scenario is not very different from the standard colliding wind scenario (e.g., Tavani & Arons 1997), as the former still features a relativistic outflow leaving the magnetosphere, but its power source is expected to be more irregular. Unfortunately, the simulations presented do not allow us properly connect the processes close to the NS with those on larger scales, and further work is required to properly characterize the development of an enhanced pulsar wind, and how it propagates and terminates, all of this in the context of the evolution of the whole interaction structure. It is known that the acceleration of electrons up to the required energies in LS 5039 is difficult (Khangulyan et al. 2008; Bosch-Ramon & Khangulyan 2025), and only an adequate characterization of that region can assess its ability to produce the particles with the required energies. Further downstream, the differences between the two scenarios are qualitatively smaller, although the more perturbed nature of the shocked flows due to magnetospheric flaring may suggest that particle acceleration on the scales of the spiral formation is stronger in that case. We emphasize that highly relativistic flows can still form outside the magnetosphere when strong flares occur, which can be a basic ingredient (as in the model proposed by Derishev & Aharonian 2012; Bosch-Ramon & Khangulyan 2025) for the production of ultra-high-energy photons such as those observed in LS 5039 (Alfaro et al. 2025). The overall nonthermal emitter is, however, likely to spread over a much larger region (e.g., Khangulyan et al. 2008; Bosch-Ramon et al. 2012). We finish by noting that the timescales inherent to the scenario explored here are diverse, from the light-crossing timescale of the magnetosphere (seconds) and the stellar wind crossing time of the interaction region (minutes to hours), to the orbital period (days).

Acknowledgments

The authors want to thank the referee for a useful report that helped to improve the clarity of the manuscript. VB-R acknowledges financial support from the Spanish Ministry of Science and Innovation under grant PID2022-136828NB-C41/AEI/10.13039/501100011033/ERDF/EU and through the María de Maeztu award to the ICCUB (CEX2024-001451-M), and from the Generalitat de Catalunya through grant 2021SGR00679. VB-R is Correspondent Researcher of CONICET, Argentina, at the IAR. Part of this work was supported by the Russian Science Foundation grant 23-22-00385 (sections 2 and 4) and BASIS foundation grant #24-1-2-25-1 (sections 1 and 3). The calculations were performed on CFCA XC50 cluster of the National Astronomical Observatory of Japan (NAOJ) and RIKEN HOKUSAI Bigwaterfall.

References

All Tables

All Figures

	Fig. 1. Example of the time profile of a flaring period; δρ is the density on top of the regular pulsar wind density during flaring periods, x = 2πtNfl/Torb the dimensionless time (see Sect. 2.1 for details), and Δ = 0.1 the width of the pulse.
In the text

Fig. 2. Colored density maps on the orbital plane, with colored arrows representing the modulus of the 3-velocity vector for the 2D models. The flaring period rates are Nf = 90 (top left, 2Dn90e10), 30 (top right, 2Dn30e10), and 10 (middle left, 2Dn10e10) flaring periods per orbit, keeping the total energy budget of the system constant. In the model 2Dn30e30, we increased the power of the flaring period by a factor of 3, while keeping Nf = 30 (middle right panel). The non-flaring cases are presented at the bottom left (nfeta0.05) and right (nfeta0.5).

In the text

	Fig. 3. Colored density maps with colored arrows representing the modulus of the 3-velocity vector for the 3D model (3Dn30e10) for various cuts: XY (top; orbital plane), XZ (middle), and ZY (bottom).
In the text

	Fig. 4. Colored density maps with colored arrows representing the modulus of the 3-velocity vector for the 3D model (3Dn30e10) for cuts on the orbital plane, at orbital phases ϕ = 0.5435 (top left), ϕ = 0.575 (top right), ϕ = 0.5812 (middle left), ϕ = 0.6033 (middle right), ϕ = 0.6467 (bottom left), and ϕ = 0.6576 (bottom right).
In the text

	Fig. 5. Zoom-in of three snapshots of the cycle followed by the interaction structure: weak pulsar wind (left), enhanced pulsar wind (middle), and the end of energy injection (right).
In the text

Fig. 6. Colored density maps obtained from an axisymmetric simulation showing qualitatively the evolution of a flare, from top left to bottom right: 0 (weak wind), 48 (flare starts), 88 (flare continuation), 141 (end of flare), 166 (pulsar wind approaches magnetosphere), and 186 s (flare about to resume) from the beginning of the flare. The grid size is ≈0.65 a; the magnetosphere radius would be ≈0.02 a in this case.

In the text

	Fig. 7. Sketch of the georotator scenario evolution, showing the pulsar angular velocity evolution with time with the ejector and the georotator (present) phases and a dashed green line marking the transition.
In the text