Open Access
Issue
A&A
Volume 669, January 2023
Article Number A21
Number of page(s) 10
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202244531
Published online 23 December 2022

© The Authors 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1 Introduction

High-mass gamma-ray binaries are among the most luminous high-energy sources in the Galaxy. These systems host a massive star and a compact object, and their radiation extends from radio waves to high-energy (HE; 0.1–100 GeV) and very-high-energy (VHE; >100 GeV) gamma rays, with the nonthermal emission being dominated by energies > 1 MeV (see, e.g., Dubus 2013; Paredes et al. 2013; Paredes & Bordas 2019, for a review). Observations of these systems also show gamma-ray orbital modulation, indicating that the related emitting region cannot be located too far from the binary (e.g., for LS 5039 and LS I +61 303; Aharonian et al. 2006; Albert et al. 2006, respectively). Approximately 10 systems pertaining to this class have been detected so far, but the power origin has not been firmly identified yet in most of them, with the important exceptions of PSR B1259-63 (Johnston et al. 1992), PSR J2032+4127 (Lyne et al. 2015), and probably LS I +61 303 (Weng et al. 2022), in which the recent detection of radio pulsations strongly points at the presence of a pulsar (see also Yoneda et al. 2020 and Volkov et al. 2021 for the contested detection of X-ray pulses from LS 5039). The uncertainty in the power origin in gamma-ray binaries has fueled a debate between two competing scenarios: the accretion-powered and the pulsar-wind-powered scenarios (e.g., Maraschi & Treves 1981; Paredes et al. 2000; Martocchia et al. 2005; Dubus 2006; Romero et al. 2007; Bosch-Ramon & Khangulyan 2009; Massi et al. 2017; see however, e.g., Papitto et al. 2012 for a pulsar model with regime transitions).

In the accretion-powered scenario, particle acceleration is typically thought to occur in the relativistic jets of a microquasar, which are fed by the accretion of stellar material onto the compact object. Jet propagation through the stellar wind medium on binary scales has been numerically studied in the relativis-tic (e.g., Perucho & Bosch-Ramon 2008; Perucho et al. 2010; Charlet et al. 2022) and nonrelativistic regimes (e.g., Yoon & Heinz 2015; Yoon et al. 2016). Leptonic and hadronic models have been adopted for the broadband emission from the relativistic jets of high-mass microquasars (e.g., Romero 2008). The effects of orbital motion on the jet kinematics and the broadband nonthermal emission have also been considered in several (semi-)analytical and numerical works (see, e.g., Bosch-Ramon & Barkov 2016; Molina & Bosch-Ramon 2018; Molina et al. 2019; Barkov & Bosch-Ramon 2021b, for some recent studies including orbital effects).

In the standard pulsar-powered scenario, a significant fraction of the pulsar spin-down power is dissipated at the interaction of a relativistic wind from the pulsar and the stellar wind (e.g., Maraschi & Treves 1981; Tavani & Arons 1997). At the shock front, particles can be efficiently accelerated and emit X-rays via synchrotron and gamma rays via inverse Compton (IC) scattering on the stellar photon field (e.g., Tavani et al. 1994; Kirk et al. 1999; Dubus 2006; Khangulyan et al. 2007); although, hadronic scenarios have also been proposed (e.g., Neronov & Chernyakova 2007). Relativistic hydrodynamic simulations show that fast reacceleration of the shocked pulsar flow should occur beyond the sonic point (at the outskirts of the binary; see, e.g., Bogovalov et al. 2008), which would influence the shocked pulsar wind emission by affecting adiabatic losses, the magnetic field, and Doppler boosting (e.g., Kong et al. 2012; Khangulyan et al. 2014b; Dubus et al. 2015; Molina & Bosch-Ramon 2020). The magnetization (and anisotropy) of the pulsar wind can have only a moderate effect on the dynamics and overall structure of the flow (Bogovalov et al. 2012), unless the pulsar wind magnetic field becomes dynamically dominant (Bogovalov et al. 2019). Instability growth and turbulence can impact the structure of the shocks and produce shocked two-wind mixing, and couple with orbit-related Coriolis forces that turn the shocked flows into an unstable spiral structure (e.g., Bosch-Ramon & Barkov 2011; Bosch-Ramon et al. 2012, 2015; Lamberts et al. 2012, 2013). Detailed modeling of the shocked-flow dynamics and radiation in the case of LS 5039, including orbital motion, has been performed by Molina & Bosch-Ramon (2020; semi-analytically), taking the emission at the Coriolis turnover into account, and by Huber et al. (2021a,b; numerically), accounting for the complex fluid dynamics through coupled relativistic hydrodynamics (RHD)-nonthermal particle calculations. The role of eccentricity in the large-scale evolution of the shocked flows has also been explored, using RHD simulation data to estimate the associated nonthermal emission (e.g., Barkov & Bosch-Ramon 2018, 2021a).

Smooth flows are generally considered in the pulsar-wind-powered scenario, but small-scale inhomogeneities can inherently arise from radiation instabilities in the acceleration phase of hot-star winds (e.g., Lucy & Solomon 1970; Runacres & Owocki 2002; Puls et al. 2008). Large-scale structures may also be present and possibly originate from the circumstellar disks of fast-rotating stars (see, e.g., Okazaki et al. 2011; Chernyakova et al. 2014, in the context of gamma-ray binaries), corotating interaction regions between different velocity streams, or magnetically confined regions at the wind base (e.g., Cranmer & Owocki 1996; Lobel & Blomme 2008). In nonaccreting systems, the radiative consequences of the presence of clumps were analytically explored, for instance, by Zdziarski et al. (2010) in the context of LS I +61 303. A general study, also of analytical nature, of the dynamical and radiative effects of a clumpy stellar wind interacting with a pulsar wind was carried out by Bosch-Ramon (2013). Later on, axisymmetric RHD simulations by Paredes-Fortuny et al. (2015) and nonthermal emission computations based on RHD simulation data by de la Cita et al. (2017) were performed to numerically study the effects of the arrival of a single clump in the two-wind interaction region. These works concluded that, depending on its size, a clump can noticeably perturb the interaction region and largely affect the output of the associated nonthermal radiation. For instance, in the mentioned works and in Chernyakova et al. (2014), it was argued that the collision of a large clump of stellar material with the pulsar wind offers a plausible explanation for the observed gamma-ray flares (e.g., PSR B1259-63) and quick X-ray variability (e.g., LS 5039, LS I +61 303) in high-mass gamma-ray binaries. Indeed, short-term variability of high-energy emission seems to be a common feature in gamma-ray binaries (e.g., Bosch-Ramon et al. 2005; Smith et al. 2009; Takahashi et al. 2009; Rea & Torres 2011; An et al. 2015; Tam et al. 2018). Studies of the clumpy-wind presence have also been performed, for instance, in the context of massive-star binaries (Pittard 2007), accreting X-ray sources (e.g., Oskinova et al. 2012), and microquasars (e.g., Araudo et al. 2009; Owocki et al. 2009; Perucho & Bosch-Ramon 2012; de la Cita et al. 2017; López-Miralles et al. 2022).

As explained above, so far only simple analytical treatments or costly numerical calculations have been used to model the emitting flow under a clumpy stellar wind in the pulsar-powered scenario. Here, we present a novel approach to modeling the effects of a clumpy stellar wind on the emitting region in that scenario; this approach aims at balancing realistic modeling and computational efficiency. The unshocked and shocked clump propagation and dynamics are modeled adopting semi-analytical Monte Carlo and hydrodynamic calculations, allowing for a general exploration of the effects of multiple clumps on the geometry evolution of the shocked two-wind region. The associated nonthermal emission and its variability induced by the clumps are analytically modeled treating the emitting particles either in an adiabatic or a radiative regime, with a particle power-law energy distribution. Gamma-ray absorption is taken into account. This method is very fast to implement but still provides the most relevant information; although, formally it can only be applied consistently to the region where the shocked pulsar wind is subsonic. The paper is structured as follows: the physical scenario and its modeling are explained in Sect. 2; the results are presented in Sect. 3; and a summary and discussion are provided in Sect. 4.

2 Physical system and modeling

2.1 Physical system

Although this study applies to high-mass gamma-ray binaries hosting young pulsars in general, we consider a system with orbit and distance similar to those of LS 5039 (e.g., Casares et al. 2005), a powerful and well-studied source. The adopted components for the binary system are a main-sequence O-type star with luminosity L* = 1039 erg s−1 and temperature T = 4 × 104 K, and a nonaccreting pulsar with spin-down power set to Psd = 3 × 1036 erg s−1. The orbital parameters are the system eccentricity e = 0.35, orbital period P = 3.91 days, and semimajor axis α = 2.1 × 1012 cm. The orientation of the system with respect to the observer is such that superior and inferior conjunction occur at orbital phases ϕ = 0.058 and ϕ = 0.716, respectively. We assume a distance of dobs = 3 kpc and an inclination with respect to the line of sight of i = 60° (favored in LS 5039 if it harbors a neutron star; Casares et al. 2005). The actual role of i is however secondary for us, as it mostly affects the orbital spectral energy distributions (SEDs) and light curves; for suborbital variability, i mainly influences the normalization of the curves.

The massive star powers a supersonic wind that is driven by line radiation pressure. This wind is assumed to be inhomogeneous, with most of its mass carried by clumps. We adopt a wind mass-loss rate of = 10−7 M yr−1 and a wind velocity approximated to a constant uw = 2 × 108 cm s−1. The pulsar drives an ultrarelativistic wind that is assumed to be cold, isotropic, and weakly magnetized (see however Bogovalov et al. 2012, Bogovalov et al. 2019; Derishev & Aharonian 2012; Bosch-Ramon 2021, and references therein). Typical bulk Lorentz factors are Γw = 105 (Khangulyan et al. 2012; Aharonian et al. 2012), but our results will not strongly depend on this value as the shocked clumps do not accelerate up to relativistic velocities; Γw can also affect the energy range of the flow emitting particles (e.g., Kennel & Coroniti 1984), but this has little effect on our conclusions and is neglected. The system parameters are summarized in Table 1, and a schematic of the system is shown in Fig. 1.

Table 1

System parameters.

2.2 Modeling

2.2.1 Smooth-wind contact discontinuity

The stellar and pulsar winds collide to form a double-bow-shock structure made of shocked wind material. The shape of the contact discontinuity (CD) separating the shocked stellar and pulsar winds can be described, in the smooth-wind case, by a simple approximation of the boundary of equilibrium between the two wind ram pressures (Cantó et al. 1996): (1)

where ϑ and ϑ1 correspond to θ and θ1 in Fig. 1 from Cantó et al. (1996). Equation (1) gives the shape of a one-dimensional CD cut, but the latter axisymmetry in the smooth-wind case allows us to derive the actual geometry of the two-dimensional CD by rotating this one-dimensional cut around the CD symmetry axis.

The ratio of the pulsar and the stellar wind momentum rates, (2)

is ≈0.08 for the selected parameters, a usual value in the literature. Such an η value means that the stellar wind momentum rate dominates and the CD bends over the pulsar, but most of the luminosity arriving at the CD still comes from the pulsar wind. The distance of the stagnation point from the pulsar is R0 = dη1/2/(1 + η1/2), where d is the orbital separation distance. The calculations are done in the rotating frame associated with the pulsar orbital motion, with ωorb being the angular velocity of the orbit described above. In this frame, due to the orbit-induced Coriolis force, the CD is slightly tilted with respect to the star–pulsar direction. Here, the tilting of the CD is introduced by rotating the entire structure counter-wise with respect to the orbital motion by an angle ~orb/uw. Given that most of the clumps are expected to get destroyed in the shocked flow (Pittard 2007; Bosch-Ramon 2013; Paredes-Fortuny et al. 2015), the actual CD shape is approximated here as a smooth-case CD plus distortions produced by the largest clumps (see below).

We treat the pulsar wind after the termination shock as subsonic because it strongly simplifies the emitting flow hydrodynamics; although, far enough from the shock the flow can become supersonic again. The region we consider, of size ≲d, is somewhat larger than the corresponding region in the smooth-wind case (Bogovalov et al. 2008), but flow irregularities induced by the clump presence may justify this choice. Further away from the star, beyond the sonic point, shocked clumps are expected to have already dissolved in the shocked two-wind structure, while the incoming stellar wind tends to be more homogeneous (see, e.g., Rubio-Díez et al. 2022, and references therein). These regions however would present flow reacceleration, local instabilities, wind mixing, and orbital effects through the Coriolis force, all likely to produce additional variability. How these factors interact with a clumpy stellar wind is, however, beyond the scope of this work.

thumbnail Fig. 1

Three-dimensional schematic (top) and two-dimensional cut (bottom) of the physical scenario comprising a massive star with a clumpy stellar wind and a pulsar. The dashed lines in the bottom panel show the trajectories of clumps that have crossed the smooth-wind contact discontinuity.

2.2.2 Clump effects on the CD

To account for the nonlinear nature of the formation of clumps, we characterize their mass distribution through a power-law function: . The index k and the volume filling factor f, which is the average-to-clump density ratio and links the clump size and mass, are treated as free parameters. The initial minimum and maximum clump radii adopted here are Rc,min = 0.01 R* and Rc,max = 0.1 R*, respectively (see Sect. 3.3 in Bosch-Ramon 2013, for a discussion). Table 2 summarizes the parameters selected for the four simulated cases with different degrees of clumpiness of the stellar wind. We note that k = 2 yields a top-heavy clump distribution, whereas k = 3 yields a bottom-heavy one. The computational method for the clump dynamics is detailed below and presented in Fig. 2.

Clumps, assumed to be spherically symmetric and of constant density, are isotropically launched from the stellar surface with their direction, initial mass and radius, and injection time being determined by a random number generator. The clump injection rate is chosen according to . Clumps move away from the stellar surface with a clump radial velocity component of uw. As the calculation is done in the rotating frame, we add an azimuthal clump velocity component of the form −orb (against orbital motion), where r is the distance from the star, to account for the effect of the orbit. Stellar rotation is neglected at this stage. During this time, clumps are assumed to expand such that f is constant, so their radii evolve as: Rc(r) = Rc(R*)(r/R*)2/3.

Depending on their orientation, some clumps reach the CD, while others leave the simulation grid without interacting with the shocked two-wind region. Clumps that reach the CD may or may not be able to cross the entire shocked two-wind region, depending on their radii. The critical radius for penetration is (Bosch-Ramon 2013) (3)

where ∆R is the thickness of the shocked two-wind region at each location. Hydrodynamic simulations for a smooth stellar wind and η values similar to ours (Bogovalov et al. 2008; Bosch-Ramon et al. 2015) indicate that ∆R is ~1/3 of the distance to the pulsar (rp) there (which in turn depends on ϑ, ϑ1, η, and d). Clumps with sizes smaller than that given in Eq. (3) get destroyed and dissolve in the shocked two-wind medium, not affecting the CD geometry substantially. Clumps with larger sizes can penetrate into the (until then) unshocked pulsar wind.

Once the clumps reach the pulsar-side boundary of the shocked two-wind region (i.e., the smooth pulsar wind termination shock), the ram pressure of the pulsar wind starts to impact the clumps. Adapting a semi-analytical hydrodynamical model from Barkov et al. (2012) for the evolution of a clump under the impact of a relativistic radial flow, we solve the equation of motion of the clumps in the pulsar wind zone. The radial component (outward from the pulsar) of the clump Lorentz factor evolves according to: (4)

where is the pulsar wind ram pressure, Γc is the clump Lorentz factor, and Mc is the clump mass. The motion of the clumps is initially decelerated toward the pulsar and then, after the pulsar wind ram pressure overpowers that of the stellar wind, accelerated away from the pulsar. Gravitational forces are not dynamically relevant compared with the wind and clump thrusts in the region of interest and can be neglected.

Within the shocked two-wind region, clumps cannot expand as they are confined by the thermal pressure of the shocked winds, but in the (previously) unshocked pulsar wind zone quick clump expansion can occur at the shocked-clump sound speed in all directions except against the pulsar wind (see, e.g., del Palacio et al. 2019). For simplicity, at this stage we assume that the shocked expanding clumps are still uniform and spherical (valid at least for the initial stages of clump expansion; Paredes-Fortuny et al. 2015), with sound speed: (5)

where is the internal specific enthalpy, and ρc is the clump mass density. As clumps do not have time to become relativistic, the adiabatic index can be fixed to 5/3, and both Γc and h are always barely above 1.

Expansion is switched off if the radius of an expanding clump reaches half the distance to the pulsar when the undisturbed CD is initially crossed (measured from the clump center). This accounts for the fact that the approximation used for the dynamical evolution of the clumps that are shocked by the pulsar wind (computed taking the clump centers as their reference positions) is valid until an expanding clump affects much of the originally unshocked pulsar wind zone (equivalently, the clump size becomes significant with respect to rp).

At that stage, the clumps quickly get deflected away from the pulsar (as ) and re-approach the original CD. Once the center of an expanded clump has intercepted that surface again, we assume that it quickly mixes with the shocked two-wind region.

We note that other Rc(r) prescriptions may be also possible because the actual clump geometry and evolution are very complex (e.g., Sundqvist et al. 2018; El Mellah et al. 2020). In the case of linear growth, for instance, the relevant quantities at the CD would be within a factor of 2 with respect to those in the adopted case, which does not qualitatively affect the conclusions of this work. The actual geometry of clumps when they fully cross the shocked two-wind region is not important, as their later expansion once shocked by the pulsar wind is fast in all directions except toward the pulsar; thus it can be considered to be roughly isotropic.

Using the above approach, which albeit crude attempts to capture the main clump evolution features found by Paredes-Fortuny et al. (2015), we can determine the location and size of the penetrating clumps at any given time, which allows us to characterize the corresponding shape and temporal evolution of the entire distorted CD. This is done by tracing radial lines from the pulsar in all directions. If one of these lines intersects the surface of a shocked clump before the line reaches the smooth-wind CD, then a distorted-CD point is placed at the intersection location. If the line reaches the smooth-wind CD, then the CD is not distorted at this location. When this is done for all relevant directions from the pulsar, one obtains a new grid of points representing the distorted CD (see Fig. 3 for a sketch of the computation process). Complementary to Fig. 1, Fig. 3 illustrates the CD between the pulsar wind and the stellar wind once it is affected by the clump presence. Clump dynamics is assumed to be affected only by the pulsar wind and not by the presence of other clumps, as the latter would be a higher-order effect and is neglected at this stage.

Table 2

Degrees of inhomogeneity of the stellar wind.

thumbnail Fig. 2

Flowchart of the clump dynamics computational process.

thumbnail Fig. 3

Two-dimensional snapshots of a CD section at 30 min intervals showing the derivation of the distorted CD shape. Radial lines (light gray lines) are traced from the pulsar (red point; not in scale) in all directions. If a line intersects the surface of a shocked clump (light blue disks) before reaching the smooth-wind CD (black dashed line), then a distorted-CD point is placed at this location; otherwise, the intersection point of the radial line with the smooth-wind CD is added to the grid. The new grid of points (dark blue points on the blue line) representing the distorted CD is fitted by a spline (blue solid line).

2.2.3 Emitter characterization

For simplicity, we take the CD surface as the location of the emission and divide this surface into N surface elements that are treated as point-like emitters. The magnetic field and particle distribution, which are determined by the local conditions on each specific point on the CD, are treated as uniform within each individual emitter. For convenience, we express the magnetic field (B) energy density for each emitter as a fraction ηB of the shocked pulsar wind pressure at the CD: uB = ηBPW. The fraction ηB is fixed throughout this work to 0.1, which corresponds to a magnetic-to-stellar photon energy density ratio of ~ 10−3 at the CD close to the pulsar location. This ηB value lies between the weakly and strongly magnetized pulsar wind cases and is consistent with a hydrodynamics approach to characterize the CD. Taking different ηB values would not significantly change our conclusions, as radiation losses would still be dominated by IC and because our main focus is on the relative flux change in specific energy bands (although broadband spectral changes would be expected). Only electron/positron cooling processes are considered at this stage, that is, IC and synchrotron emission, because they are the most efficient radiation processes in compact binaries (see, e.g., Bosch-Ramon & Khangulyan 2009).

The complexity of the emitter structure is very high, as shown by previous works dealing only with some of the aspects of the scenario studied here (see, e.g., Dubus et al. 2015; de la Cita et al. 2017; Huber et al. 2021a, where RHD simulations were combined with radiative calculations). However, as we are interested in the typical clump-induced emission time variability at certain energies and not in deriving a detailed SED, we simplify the N point-like emitters as being in one of two extreme cooling regimes: radiative or adiabatic.

In the fully radiative regime, the energy of the emitting particles is radiated away right after reaching the CD. In this case, the particle energy distribution is determined by the acceleration (injection) energy distribution and the relevant radiative losses only (see below). In the adiabatic regime, the total energy density in the post-shock region, ≈3 Pw times the relevant volume (known only approximately; see below), determines the nonthermal population there. In this case, the particle energy distribution coincides with the acceleration one. Intermediate cases are more realistic, with an expected transition between regimes from the higher particle energies (radiative; more important close to the pulsar) to the lower particle energies (adiabatic; more important far from the pulsar). This transition would manifest itself in features of the IC and synchrotron SED that are hinted when comparing our results in both regimes. Albeit simple, our approach allows for the derivation of general conclusions on the radiation behavior because results for both cooling regimes are provided. The normalizations of the particle distributions in both regimes are just approximately consistent with each other due to our simplified treatment of the hydrodynamics in the adiabatic case.

In our model, relativistic particles are injected following a power-law distribution in energy, E, with an exponential cutoff of the form (6)

where we adopt p = 2 as a fiducial value. Significantly softer (harder) injection energy distributions would yield substantially lower (higher) fluxes in the higher-energy regions of the synchrotron and IC components. Further, injection energy distributions more complex than a power-law function, such as an additional narrow component peaking around Γw (see, e.g., Dubus et al. 2015), would lead to more complex SEDs. Again, given our focus on the clump-induced relative variability within specific energy bands, the adopted injection distribution is enough for our purposes.

The maximum electron energy, characterized by Emax in Eq. (6), is mainly constrained by acceleration, radiative losses, and escape. As particles with multi-TeV energies exist in gamma-ray binaries, we assume a sufficient acceleration rate, with efficiency ηacc = 0.1 and timescale tacc = E/(ηaccqcB), where q is the electron charge. As noted, the relevant cooling channels in the shocked pulsar wind are synchrotron emission and IC. The energy-loss rate for synchrotron (Ėsyn) is derived from Blumenthal & Gould (1970) and for IC (ĖIC) from Khangulyan et al. (2014a) assuming an anisotropic distribution of blackbody photons. Equating the acceleration timescale tacc with the total radiation cooling timescale, tcoo1 = E/(|Ėsyn| + |ĖIC|), yields . We also consider the Hillas criterion (Hillas 1984), which is similar to the constraint from advection escape and gives an upper limit of the maximum energy: , where the spatial scale rs represents the size of the source, taken here simply as ~d. Then, , which is typically ~10 TeV.

The minimum electron distribution energy has been fixed to Emin = 100 MeV. This parameter strongly affects the lowest-energy part of the IC SED and weakly the normalization of the distribution of particles when their total energy is fixed. The former affects the contribution of IC emission below X-rays, whereas the latter effect can be “absorbed” by changing the total energy in the power-law distribution of particles. This nonthermal energy, which affects the normalization of the SED linearly, is taken here to be the entirety of the shocked pulsar wind energy, but it can be lower (e.g., most of the energy may be in a narrow particle energy distribution).

In the radiative regime, the electron injection distribution is normalized by the pulsar wind luminosity within the solid angle of each emitter, ∆LΩ: (7)

The corresponding cooled particle population can be described as: (8)

On the other hand, in the adiabatic regime the electron distribution ∆N(E) follows the energy dependence in Eq. (6) and is normalized by the internal energy ∆E stored within an emitter of approximate volume dV = dSsphhs, where dSsph = rpdΩ, and hs is the thickness of the shocked pulsar wind shell. To determine hs, we take hs = 0.2 rp based on relativistic hydrodynamical simulations (e.g., Bosch-Ramon et al. 2015). Thus, one has: (9)

Once the particle energy distribution for each point-like emitter has been obtained, the associated synchrotron and IC SEDs are computed using the formulae given in Blumenthal & Gould (1970) and Khangulyan et al. (2014a), respectively. Given that the flow is only weakly relativistic in the subsonic region of the shocked pulsar wind and the orientation distribution of the flow velocity is broad in the CD (even broader in the case of a distorted CD; Paredes-Fortuny et al. 2015), Doppler boosting has not been taken into account. Gamma-ray absorption in the stellar photon field through pair creation has been calculated using the cross section from Gould & Schréder (1967). The radiation of the created pairs may affect the X-ray and gamma-ray emission (e.g., Bednarek 2007; Bosch-Ramon et al. 2008), but we neglect this contribution at this stage. Due to the clump presence and orbital motion, changes in the interaction angles of IC and gamma-ray absorption and in the CD geometry (also affecting synchrotron radiation through B) must be taken into account.

3 Results

3.1 Smooth case

We have computed the SEDs and light curves in the case of a smooth stellar wind along the entire orbit including regions of the CD up to 1.5d from the star; cases with different emitting sizes mostly affect the normalization and are not presented here. Figure 4 shows the synchrotron and IC SEDs for the two cooling regimes (adiabatic and radiative) and two orbital phases: inferior (ϕ = 0.716 -INF-) and superior conjunction (ϕ = 0.058 -SUP-).

The results obtained are similar to those found in previous works (see Sect. 1). In the adiabatic regime, synchrotron and IC radiation are rather hard, with the synchrotron component peaking at soft gamma rays, and the IC component softening above ~1–10 GeV due to the Klein-Nishina (KN) effect in the cross section and peaking around ~10–100 GeV. In the radiative regime, the synchrotron peak broadens, but the SED is still hard down to ~10 eV due to IC losses in the KN limit, where it becomes flat due to dominant IC losses in the Thomson regime. The (unabsorbed) IC SED is still moderately hard at ~0.1–100 GeV due to IC KN losses, softens above that energy due to synchrotron becoming dominant plus the KN effect, and flattens below ~0.1 GeV due to IC Thomson losses. Below ~1 MeV, corresponding to (cooled) electrons with E < Emin, the IC spectrum hardens again. The effect of gamma-ray absorption in the SED is the strongest and starts at the lowest gamma-ray energy (~30 GeV) at superior conjunction when the photonphoton interaction angle is the largest. Absorption becomes less severe the farther the pulsar is from that orbital phase, becoming a minor effect around inferior conjunction, mostly due to the interaction angle decreasing (an effect smoothed further by the extended emitter; e.g., Khangulyan et al. 2008). Pulsar-star distance changes are a secondary factor for small-to-moderate eccentricities that enhances absorption toward periastron. For simplicity, soft X-ray absorption due to the photoelectric effect has not been considered at this stage. We did not study radio wavelengths either, as this emission is expected to be severely free-free absorbed so close from the star (see, e.g., Molina & Bosch-Ramon 2020, and references therein).

Figure 5 shows the orbital light curves in the case of a smooth stellar wind for the radiative and adiabatic regimes and three energy bands: X-rays (1–100 keV; synchrotron), HE gamma rays (0.1–10 GeV; IC), and VHE gamma rays (>100 GeV; IC). For the X-ray light curves, we only show the synchrotron components (adiabatic and radiative), because IC X-rays are produced deep in the adiabatic regime (see Fig. 4), and the corresponding curve is lower than any synchrotron curve. The HE and VHE gamma-ray light curves are by contrast shown only for IC (both regimes), which is the dominant emission at those energies. The light curves are mostly explained by changes in the emitter spatial volume and angular size as seen from the pulsar (affecting particle normalization), in the distance from the emitter to the pulsar (affecting B and synchrotron in the X-rays) and the star (more mildly affecting IC target photons), and by angular effects (affecting IC in the HE and VHE and gamma-ray absorption in the VHE). The fact that the emitter is extended with size ~d and has a complex geometry changing due to eccentricity introduces some variations, but overall the results are as expected from previous work, with X-rays peaking around apastron, and the HE and the VHE gamma-rays peaking around superior and inferior conjunction, respectively. The VHE light curve presents a secondary peak before periastron, but similar features were already found in previous work (e.g., Khangulyan et al. 2008; Dubus et al. 2008; Takahashi et al. 2009; Molina & Bosch-Ramon 2020).

thumbnail Fig. 4

SEDs for the smooth-wind case for the (a) adiabatic (green) and (b) radiative (blue) regimes considering regions up to 1.5d from the star at inferior conjunction (INF; ϕ = 0.716) and superior conjunction (SUP; ϕ = 0.058). The synchrotron (colored dashed) and IC (colored solid) components are plotted, and the unabsorbed spectra (colored dotted) are also shown for comparison. As a reference, the shaded areas indicate the synchrotron (~10–2 eV) and IC (~107 eV) photon energies corresponding to the (adiabatic-to-radiative) transition particle energy for all individual emitters in the CD.

thumbnail Fig. 5

Light curves for a smooth stellar wind for the adiabatic (green dashed) and radiative (blue solid) regimes including regions of the CD up to 1.5d from the star for synchrotron emission in the 1–100 keV band (top) and IC emission in the 0.1–10 GeV (middle) and >100 GeV (bottom) energy bands. The vertical dashed lines labeled SUP and INF mark the superior and inferior conjunction, respectively. We show two orbital cycles for clarity.

thumbnail Fig. 6

Two-dimensional cut on the orbital plane of the smooth CD (solid gray line) at an intermediate orbital phase (ϕ = 0.28) showing the trajectories of multiple clumps. The light-gray shaded area gives a rough estimate of the thickness of the shocked two-wind region. The blue and yellow lines correspond to high-density and low-density clumps for stellar wind filling factors of f = 0.01 and f = 0.1, respectively. The clump radii when launched from the stellar surface are 0.1 R* (solid) and 0.06 R* (dashed). The colored shaded areas show the regions along the pulsar wind shock front within which the largest high-density (blue shaded area) and low-density (yellow shaded area) clumps can penetrate.

thumbnail Fig. 7

Evolution over time of the radii of the clumps whose trajectories are shown in Fig. 6. Line patterns are the same as in Fig. 6.

thumbnail Fig. 8

Two-dimensional snapshots of the CD on the orbital plane at 1 h intervals within a period of 3 h for (a) a top-heavy (k = 2) high-density (f = 0.01), (b) a bottom-heavy (k = 3) high-density (f = 0.01), (c) a top-heavy (k = 2) low-density (f = 0.1), and (d) a bottom-heavy (k = 3) low-density (f = 0.1) clump distribution. The black dashed lines correspond to a smooth stellar wind. The horizontal gray lines mark the CD size used to compute the light curves in Fig. 9.

3.2 Clumpy-wind case

To study the short-term variability induced by a clumpy stellar wind, we focus on an intermediate orbital phase (ϕ = 0.28), at which the line of sight is perpendicular to the pulsar-star line, and the three mentioned energy bands (X-rays and HE and VHE gamma rays) and cooling regimes (adiabatic and radiative). In Fig. 6, we show the geometry on the orbital plane of the smooth CD and the region filled with the shocked two-wind flow. We also show the trajectories of clumps of different initial masses/radii and densities that manage to penetrate the shocked two-wind region close to the stagnation point, where the minimum initial clump radii for penetration are ~0.02R* for f = 0.01 and ~0.05 R* for f = 0.1.

Although the largest clumps with f = 0.01 can reach the unshocked pulsar wind region almost everywhere in the explored CD region (blue shaded area in Fig. 6), for the largest clumps with f = 0.1 the CD crossing is restricted to the vicinity of the stagnation point (yellow shaded area). Figure 7 shows the time evolution of the radii of these clumps from the time of entering the shocked two-wind region until the time clump expansion is halted. All these clumps reach the same final size (half rp at the crossing point) because they share the same penetration location; although, smaller or lighter clumps tend to evolve faster. Clumps entering farther away from the smooth-CD symmetry axis may not reach again the smooth-CD location within the simulated region.

In Fig. 8, we plot cuts of the shocked two-wind region on the orbital plane, for stellar winds of different degrees of inhomogeneity over a period of 3 h at 1 h intervals. The arrival of clumps to the CD deforms the overall interaction structure, generating relatively quick and global variations in the geometry of the emitting shocked pulsar wind. In the most extreme case (top-left panel), clumps frequently reach quite close to the pulsar, radically distorting the CD geometry. Otherwise, for a bottom-heavy low-density clump distribution (bottom-right panel), the results in the explored narrow phase interval essentially coincide with a smooth-wind scenario mainly because of the small number of clumps that are large enough to cross the shocked two-wind region. In intermediate cases (bottom-left and top-right panels), moderate but significant CD variations may still be seen.

Figure 9 shows the light curves for different stellar wind scenarios in both cooling regimes over a period of 3 h. The light curves are computed considering an emitter size extending up to a distance of 1.5 d from the star. Referring to the adiabatic X-ray light curves (top left panel) as an example, one can see that the light curves of the two most inhomogeneous stellar winds (blue solid and dashed lines) exhibit large and rapid variability. For a relatively low degree of clumpiness (green dashed–dotted line), large flaring events still occur, albeit more sparsely, while for the least clumpy case (yellow dotted line) variability subsides almost completely. The same trends in variability are observed in all panels. The light curves for the smooth-wind case are not shown for simplicity but largely coincide with the bottom-heavy low-density light curves.

Using Fig. 4, one can evaluate which light curve, either radiative or adiabatic, is more realistic. In all the studied energy bands, the emission in the region of interest is produced in the radiative regime, and thus, despite the simplicity of the approach, this regime provides to first order a reasonably realistic account of the behavior of the emission. In this regime, the relative changes in flux are moderate in the explored narrow phase interval, of ~20% (top-heavy and bottom-heavy dense clumps), 10–100% (top-heavy light clumps; note the large X-ray spike in the top right panel), and ≲10% (bottom-heavy light clumps), with variability timescales of ~0.1–1 h.

It is worth noting that the adiabatic light curves present significantly higher variability than the radiative ones in the dense-clump cases. Thus, although radiative losses are generally dominant on the scales we explored, including supersonic flow regions (not considered here) could lead to higher variability, as adiabatic losses can become stronger there (Khangulyan et al. 2014b). For the emission produced more deeply in the radiative regime, the impact of hydrodynamics is expected to be small.

We also assessed whether there is a reduction in variability due to emission coming from CD regions potentially more or less affected by clumps. We compared our results with light curves for emitting regions extending closer (up to 1.2 d) or farther from the star (up to 1.8 d). Larger sizes of the emitting region do not smooth out the relative changes in the flux, mostly affecting just the normalization, which indicates that the clump presence still affects the emitter relatively far from the pulsar.

thumbnail Fig. 9

Light curves over a period of 3 h for the adiabatic (left column) and radiative (right column) regimes at energies 1–100 keV (synchrotron only; top rows) and 0.1–10 GeV and >100 GeV (IC only; middle and bottom rows, respectively) for a top-heavy high-density ((k = 2, f = 0.01; solid), a bottom-heavy high-density (k = 3, f = 0.01; dashed), a top-heavy low-density (k = 2, f = 0.1; dashed-dotted), and a bottom-heavy low-density (k = 3, f = 0.1; dotted) clump distribution.

4 Discussion

In general, the predicted luminosities in the explored energy bands are loosely similar to those of LS 5039 at a few-kpc distance. For the corresponding fluxes, present instrumentation may be able to trace the short-term relative variations in X-rays and VHE gamma rays (see typical rate uncertainties in, e.g., Martocchia et al. 2005; Aharonian et al. 2006, for X-rays and VHE gamma rays, respectively), and the future Cherenkov Telescope Array would reach significantly finer levels of variability at VHE (e.g., Paredes et al. 2013). At HE gamma rays, several orbital cycles must be folded in phase to get statistically meaningful orbital light curves, which prevents studying (nonextreme) variability on suborbital timescales.

Clump interaction rate

The frequency of interaction of the largest clumps with the shocked pulsar wind can be estimated from the mass distribution as: (10)

where the solid angle of interaction ∆Ω = π(dp/d)2 depends on the radius dp of the region that can be penetrated by clumps. For a moderately clumpy stellar wind (i.e., top heavy, low density), dp ≈ 0.5d (yellow shaded area in Fig. 7), which yields an interaction rate of s−1 or ~21 clumps per orbit and an average time between consecutive interactions of s. The duration of each individual interaction is on the order of the clump residence time in the region τc = dp/u ~ 6 × 103 s. This implies that approximately 40% of the time one large clump is interacting with the CD, and ≈ 14% of the time two large clumps may simultaneously interact with the CD, meaning several times per orbit. Even in the least clumpy case explored here, a few of the largest clumps would still be expected to (individually) interact with the CD per orbit.

Effect of η

Although our study focused on η = 0.08, we have run a few trials in low resolution for different η values for completeness. The increase in the shocked two-wind shell thickness of the termination shock for larger η values shrinks the penetrable region to the very close vicinity of the pulsar. Larger-mass clumps that can penetrate tend to grow significantly and cause large distortions of the CD, but the overall effect is a reduced variability due to a low rate of interactions. Otherwise, for smaller η values, the termination shock wraps closer to the pulsar, extending the penetrable region to almost the entire CD length, effectively increasing the interaction rate. Clumps entering relatively close to the stagnation point will have long residence times (although expansion will halt soon after crossing), while those entering from the outer regions will evolve and exit more quickly. In this second case, overall, variability becomes faster and more intense than in the fiducial calculations.

5 Summary

We have studied the effects of different clumpy stellar winds interacting with a relativistic pulsar wind in a high-mass binary, adopting different filling factors (f = 0.1 and 0.01) and levels of clumpiness (, top heavy; and , bottom heavy). Our results for the SED and orbital light curves in the smooth-wind case are similar to those presented in the literature. The presence of clumps does not seem in general to strongly modify the emitter, and thus smooth-wind models seem enough to capture the most prominent behavior of the source, but there is clear clumpy-wind-induced variability in the X-ray and HE and VHE gamma-ray light curves that can be detectable by high-energy instrumentation.

Regarding the different cases explored, one finds that the bottom-heavy light-clump distribution can be, most of the time, indistinguishable from a smooth wind in the wind-interaction region explored in this work. In the remaining cases (denser and/or top-heavy clump distributions), the presence of clumps may be traced in the X-ray light curve and less prominently in the VHE gamma-ray light curve, with predicted flux variations of ~ 10–20% in the radiative regime. These variations could reach ~100% in some rare events, that is, the arrival to the CD of particularly large clumps or even of several large clumps simultaneously. These events, despite being more sporadic, would be very noticeable and can occur even in the least clumpy case explored in this work. Adding extra factors such as clump-induced effects further downstream of the emitting flow may reveal a stronger impact of clumps.

Acknowledgements

The authors want to thank the referee, Achim Felmeier, for constructive and useful comments that helped to improve the manuscript. The authors wish to acknowledge Edgar Molina for helpful discussions. This work has received financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grant PID2019-105510GB-C31 and through the “Unit of Excellence Maria de Maeztu 2020-2023” award to the Institute of Cosmos Sciences (CEX2019-000918-M). V.B-R. is Correspondent Researcher of CONICET, Argentina, at the IAR.

References

  1. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 460, 743 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Aharonian, F. A., Bogovalov, S. V., & Khangulyan, D. 2012, Nature, 482, 507 [NASA ADS] [CrossRef] [Google Scholar]
  3. Albert, J., Aliu, E., Anderhub, H., et al. 2006, Science, 312, 1771 [NASA ADS] [CrossRef] [Google Scholar]
  4. An, H., Bellm, E., Bhalerao, V., et al. 2015, ApJ, 806, 166 [NASA ADS] [CrossRef] [Google Scholar]
  5. Araudo, A. T., Bosch-Ramon, V., & Romero, G. E. 2009, A&A, 503, 673 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Barkov, M. V., & Bosch-Ramon, V. 2018, MNRAS, 479, 1320 [NASA ADS] [CrossRef] [Google Scholar]
  7. Barkov, M. V., & Bosch-Ramon, V. 2021a, Universe, 7, 277 [NASA ADS] [CrossRef] [Google Scholar]
  8. Barkov, M. V., & Bosch-Ramon, V. 2021b, MNRAS, 510, 3479 [Google Scholar]
  9. Barkov, M. V., Aharonian, F. A., Bogovalov, S. V., Kelner, S. R., & Khangulyan, D. 2012, ApJ, 749, 119 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bednarek, W. 2007, A&A, 464, 259 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237 [Google Scholar]
  12. Bogovalov, S. V., Khangulyan, D. V., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2008, MNRAS, 387, 63 [Google Scholar]
  13. Bogovalov, S. V., Khangulyan, D., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2012, MNRAS, 419, 3426 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bogovalov, S. V., Khangulyan, D., Koldoba, A., Ustyugova, G. V., & Aharonian, F. 2019, MNRAS, 490, 3601 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bosch-Ramon, V. 2013, A&A, 560, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Bosch-Ramon, V. 2021, A&A, 645, A86 [EDP Sciences] [Google Scholar]
  17. Bosch-Ramon, V., & Barkov, M. V. 2011, A&A, 535, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Bosch-Ramon, V., & Barkov, M. V. 2016, A&A, 590, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Bosch-Ramon, V., & Khangulyan, D. 2009, Int. J. Mod. Phys. D, 18, 347 [Google Scholar]
  20. Bosch-Ramon, V., Paredes, J. M., Ribó, M., et al. 2005, ApJ, 628, 388 [NASA ADS] [CrossRef] [Google Scholar]
  21. Bosch-Ramon, V., Khangulyan, D., & Aharonian, F. A. 2008, A&A, 482, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Bosch-Ramon, V., Barkov, M. V., Khangulyan, D., & Perucho, M. 2012, A&A, 544, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Bosch-Ramon, V., Barkov, M. V., & Perucho, M. 2015, A&A, 577, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Cantó, J., Raga, A. C., & Wilkin, F. P. 1996, ApJ, 469, 729 [Google Scholar]
  25. Casares, J., Ribó, M., Ribas, I., et al. 2005, MNRAS, 364, 899 [Google Scholar]
  26. Charlet, A., Walder, R., Marcowith, A., et al. 2022, A&A, 658, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Chernyakova, M., et al. 2014, MNRAS, 439, 432 [NASA ADS] [CrossRef] [Google Scholar]
  28. Cranmer, S. R., & Owocki, S. P. 1996, ApJ, 462, 469 [NASA ADS] [CrossRef] [Google Scholar]
  29. de la Cita, V. M., del Palacio, S., Bosch-Ramon, V., et al. 2017, A&A, 604, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. del Palacio, S., Bosch-Ramon, V., & Romero, G. E. 2019, A&A, 623, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Derishev, E. V., & Aharonian, F. A. 2012, in American Institute of Physics Conference Series, High Energy Gamma-Ray Astronomy: 5th International Meeting on High Energy Gamma-Ray Astronomy, eds. F. A. Aharonian, W. Hofmann, & F. M. Rieger, 1505, 402 [NASA ADS] [Google Scholar]
  32. Dubus, G. 2006, A&A, 451, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Dubus, G. 2013, A&AR, 21, 64 [NASA ADS] [CrossRef] [Google Scholar]
  34. Dubus, G., Cerutti, B., & Henri, G. 2008, A&A, 477, 691 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Dubus, G., Lamberts, A., & Fromang, S. 2015, A&A, 581, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. El Mellah, I., Grinberg, V., Sundqvist, J. O., Driessen, F. A., & Leutenegger, M. A. 2020, A&A, 643, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Gould, R. J., & Schréder, G. P. 1967, Phys. Rev., 155, 1404 [Google Scholar]
  38. Hillas, A. M. 1984, ARA&A, 22, 425 [Google Scholar]
  39. Huber, D., Kissmann, R., Reimer, A., & Reimer, O. 2021a, A&A, 646, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Huber, D., Kissmann, R., & Reimer, O. 2021b, A&A, 649, A71 [EDP Sciences] [Google Scholar]
  41. Johnston, S., Manchester, R. N., Lyne, A. G., et al. 1992, ApJ, 387, L37 [Google Scholar]
  42. Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 710 [Google Scholar]
  43. Khangulyan, D., Hnatic, S., Aharonian, F., & Bogovalov, S. 2007, MNRAS, 380, 320 [Google Scholar]
  44. Khangulyan, D., Aharonian, F., & Bosch-Ramon, V. 2008, MNRAS, 383, 467 [Google Scholar]
  45. Khangulyan, D., Aharonian, F. A., Bogovalov, S. V., & Ribó, M. 2012, ApJ, 752, L17 [NASA ADS] [CrossRef] [Google Scholar]
  46. Khangulyan, D., Aharonian, F. A., & Kelner, S. R. 2014a, ApJ, 783, 100 [Google Scholar]
  47. Khangulyan, D., Bogovalov, S. V., & Aharonian, F. A. 2014b, in International Journal of Modern Physics Conference Series, 28, 1460169 [NASA ADS] [CrossRef] [Google Scholar]
  48. Kirk, J. G., Ball, L., & Skjæraasen, O. 1999, Astropart. Phys., 10, 31 [NASA ADS] [CrossRef] [Google Scholar]
  49. Kong, S. W., Cheng, K. S., & Huang, Y. F. 2012, ApJ, 753, 127 [NASA ADS] [CrossRef] [Google Scholar]
  50. Lamberts, A., Dubus, G., Lesur, G., & Fromang, S. 2012, A&A, 546, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Lamberts, A., Fromang, S., Dubus, G., & Teyssier, R. 2013, A&A, 560, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Lobel, A., & Blomme, R. 2008, ApJ, 678, 408 [NASA ADS] [CrossRef] [Google Scholar]
  53. López-Miralles, J., Perucho, M., Martí, J. M., Migliari, S., & Bosch-Ramon, V. 2022, A&A, 661, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Lucy, L. B., & Solomon, P. M. 1970, ApJ, 159, 879 [Google Scholar]
  55. Lyne, A. G., Stappers, B. W., Keith, M. J., et al. 2015, MNRAS, 451, 581 [Google Scholar]
  56. Maraschi, L., & Treves, A. 1981, MNRAS, 194, 1P [Google Scholar]
  57. Martocchia, A., Motch, C., & Negueruela, I. 2005, A&A, 430, 245 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Massi, M., Migliari, S., & Chernyakova, M. 2017, MNRAS, 468, 3689 [Google Scholar]
  59. Molina, E., & Bosch-Ramon, V. 2018, A&A, 618, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Molina, E., & Bosch-Ramon, V. 2020, A&A, 641, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Molina, E., del Palacio, S., & Bosch-Ramon, V. 2019, A&A, 629, A129 [EDP Sciences] [Google Scholar]
  62. Neronov, A., & Chernyakova, M. 2007, Ap&SS, 309, 253 [NASA ADS] [CrossRef] [Google Scholar]
  63. Okazaki, A. T., Nagataki, S., Naito, T., et al. 2011, PASJ, 63, 893 [NASA ADS] [Google Scholar]
  64. Oskinova, L. M., Feldmeier, A., & Kretschmar, P. 2012, MNRAS, 421, 2820 [NASA ADS] [CrossRef] [Google Scholar]
  65. Owocki, S. P., Romero, G. E., Townsend, R. H. D., & Araudo, A. T. 2009, ApJ, 696, 690 [NASA ADS] [CrossRef] [Google Scholar]
  66. Papitto, A., Torres, D. F., & Rea, N. 2012, ApJ, 756, 188 [NASA ADS] [CrossRef] [Google Scholar]
  67. Paredes, J. M., & Bordas, P. 2019, in Frontier Research in Astrophysics - III [Google Scholar]
  68. Paredes, J. M., Martí, J., Ribó, M., & Massi, M. 2000, Science, 288, 2340 [Google Scholar]
  69. Paredes, J. M., Bednarek, W., Bordas, P., et al. 2013, Astropart. Phys., 43, 301 [NASA ADS] [CrossRef] [Google Scholar]
  70. Paredes-Fortuny, X., Bosch-Ramon, V., Perucho, M., & Ribó, M. 2015, A&A, 574, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Perucho, M., & Bosch-Ramon, V. 2008, A&A, 482, 917 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Perucho, M., & Bosch-Ramon, V. 2012, A&A, 539, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Perucho, M., Bosch-Ramon, V., & Khangulyan, D. 2010, A&A, 512, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Pittard, J. M. 2007, ApJ, 660, L141 [NASA ADS] [CrossRef] [Google Scholar]
  75. Puls, J., Vink, J. S., & Najarro, F. 2008, A&AR, 16, 209 [Google Scholar]
  76. Rea, N., & Torres, D. F. 2011, ApJ, 737, L12 [NASA ADS] [CrossRef] [Google Scholar]
  77. Romero, G. E. 2008, in Revista Mexicana de Astronomia y Astrofisica Conference Series, 33, 82 [NASA ADS] [Google Scholar]
  78. Romero, G. E., Okazaki, A. T., Orellana, M., & Owocki, S. P. 2007, A&A, 474, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Rubio-Díez, M. M., Sundqvist, J. O., Najarro, F., et al. 2022, A&A, 658, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Runacres, M. C., & Owocki, S. P. 2002, A&A, 381, 1015 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Smith, A., Kaaret, P., Holder, J., et al. 2009, ApJ, 693, 1621 [NASA ADS] [CrossRef] [Google Scholar]
  82. Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Takahashi, T., Kishishita, T., Uchiyama, Y., et al. 2009, ApJ, 697, 592 [Google Scholar]
  84. Tam, P. H. T., He, X. B., Pal, P. S., & Cui, Y. 2018, ApJ, 862, 165 [NASA ADS] [CrossRef] [Google Scholar]
  85. Tavani, M., & Arons, J. 1997, ApJ, 477, 439 [NASA ADS] [CrossRef] [Google Scholar]
  86. Tavani, M., Arons, J., & Kaspi, V. M. 1994, ApJ, 433, L37 [NASA ADS] [CrossRef] [Google Scholar]
  87. Volkov, I., Kargaltsev, O., Younes, G., Hare, J., & Pavlov, G. 2021, ApJ, 915, 61 [Google Scholar]
  88. Weng, S.-S., Qian, L., Wang, B.-J., et al. 2022, Nat. Astron., 6, 698 [NASA ADS] [CrossRef] [Google Scholar]
  89. Yoneda, H., Makishima, K., Enoto, T., et al. 2020, Phys. Rev. Lett., 125, 111103 [Google Scholar]
  90. Yoon, D., & Heinz, S. 2015, ApJ, 801, 55 [NASA ADS] [CrossRef] [Google Scholar]
  91. Yoon, D., Zdziarski, A. A., & Heinz, S. 2016, MNRAS, 456, 3638 [NASA ADS] [CrossRef] [Google Scholar]
  92. Zdziarski, A. A., Neronov, A., & Chernyakova, M. 2010, MNRAS, 403, 1873 [Google Scholar]

All Tables

Table 1

System parameters.

Table 2

Degrees of inhomogeneity of the stellar wind.

All Figures

thumbnail Fig. 1

Three-dimensional schematic (top) and two-dimensional cut (bottom) of the physical scenario comprising a massive star with a clumpy stellar wind and a pulsar. The dashed lines in the bottom panel show the trajectories of clumps that have crossed the smooth-wind contact discontinuity.

In the text
thumbnail Fig. 2

Flowchart of the clump dynamics computational process.

In the text
thumbnail Fig. 3

Two-dimensional snapshots of a CD section at 30 min intervals showing the derivation of the distorted CD shape. Radial lines (light gray lines) are traced from the pulsar (red point; not in scale) in all directions. If a line intersects the surface of a shocked clump (light blue disks) before reaching the smooth-wind CD (black dashed line), then a distorted-CD point is placed at this location; otherwise, the intersection point of the radial line with the smooth-wind CD is added to the grid. The new grid of points (dark blue points on the blue line) representing the distorted CD is fitted by a spline (blue solid line).

In the text
thumbnail Fig. 4

SEDs for the smooth-wind case for the (a) adiabatic (green) and (b) radiative (blue) regimes considering regions up to 1.5d from the star at inferior conjunction (INF; ϕ = 0.716) and superior conjunction (SUP; ϕ = 0.058). The synchrotron (colored dashed) and IC (colored solid) components are plotted, and the unabsorbed spectra (colored dotted) are also shown for comparison. As a reference, the shaded areas indicate the synchrotron (~10–2 eV) and IC (~107 eV) photon energies corresponding to the (adiabatic-to-radiative) transition particle energy for all individual emitters in the CD.

In the text
thumbnail Fig. 5

Light curves for a smooth stellar wind for the adiabatic (green dashed) and radiative (blue solid) regimes including regions of the CD up to 1.5d from the star for synchrotron emission in the 1–100 keV band (top) and IC emission in the 0.1–10 GeV (middle) and >100 GeV (bottom) energy bands. The vertical dashed lines labeled SUP and INF mark the superior and inferior conjunction, respectively. We show two orbital cycles for clarity.

In the text
thumbnail Fig. 6

Two-dimensional cut on the orbital plane of the smooth CD (solid gray line) at an intermediate orbital phase (ϕ = 0.28) showing the trajectories of multiple clumps. The light-gray shaded area gives a rough estimate of the thickness of the shocked two-wind region. The blue and yellow lines correspond to high-density and low-density clumps for stellar wind filling factors of f = 0.01 and f = 0.1, respectively. The clump radii when launched from the stellar surface are 0.1 R* (solid) and 0.06 R* (dashed). The colored shaded areas show the regions along the pulsar wind shock front within which the largest high-density (blue shaded area) and low-density (yellow shaded area) clumps can penetrate.

In the text
thumbnail Fig. 7

Evolution over time of the radii of the clumps whose trajectories are shown in Fig. 6. Line patterns are the same as in Fig. 6.

In the text
thumbnail Fig. 8

Two-dimensional snapshots of the CD on the orbital plane at 1 h intervals within a period of 3 h for (a) a top-heavy (k = 2) high-density (f = 0.01), (b) a bottom-heavy (k = 3) high-density (f = 0.01), (c) a top-heavy (k = 2) low-density (f = 0.1), and (d) a bottom-heavy (k = 3) low-density (f = 0.1) clump distribution. The black dashed lines correspond to a smooth stellar wind. The horizontal gray lines mark the CD size used to compute the light curves in Fig. 9.

In the text
thumbnail Fig. 9

Light curves over a period of 3 h for the adiabatic (left column) and radiative (right column) regimes at energies 1–100 keV (synchrotron only; top rows) and 0.1–10 GeV and >100 GeV (IC only; middle and bottom rows, respectively) for a top-heavy high-density ((k = 2, f = 0.01; solid), a bottom-heavy high-density (k = 3, f = 0.01; dashed), a top-heavy low-density (k = 2, f = 0.1; dashed-dotted), and a bottom-heavy low-density (k = 3, f = 0.1; dotted) clump distribution.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.