Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system
Xavier Paredes-Fortuny, Valentí Bosch-Ramon, Manel Perucho, Marc Ribó
AAstronomy & Astrophysics manuscript no. clump c (cid:13)
ESO 2018October 5, 2018
Simulations of an inhomogeneous stellar wind interacting with apulsar wind in a binary system
X. Paredes-Fortuny , V. Bosch-Ramon , M. Perucho , and M. Ribó ,(cid:63) Departament d’Astronomia i Meteorologia, Institut de Ciènces del Cosmos, Universitat de Barcelona, IEEC-UB, Martí i Franquès1, E-08028 Barcelona, Spain, e-mail: [email protected] Departament d’Astronomia i Astrofísica, Universitat de València, Av. Vicent Andrés Estellés s / n, 46100 Burjassot (València), SpainReceived - / Accepted -
ABSTRACT
Context.
Binary systems containing a massive star and a non-accreting pulsar present strong interaction between the stellar and thepulsar winds. The properties of this interaction, which largely determine the non-thermal radiation in these systems, strongly dependon the structure of the stellar wind, which can be clumpy or strongly anisotropic, as in Be stars.
Aims.
We study numerically the influence of inhomogeneities in the stellar wind on the structure of the two-wind interaction region.
Methods.
We carried out for the first time axisymmetric, relativistic hydrodynamical simulations, with Lorentz factors of ∼ Results.
For typical system parameters, and adopting a stellar wind inhomogeneity with a density contrast χ (cid:38)
10, clumps withradii of a few percent of the binary size can significantly perturb the two-wind interaction region, pushing the two-wind interface to (cid:46)
40% of the initial distance to the pulsar. After it is shocked, the inhomogeneity quickly expands and is disrupted when it reaches thesmallest distance to the pulsar. It eventually fragments, being advected away from the binary system. The whole interaction region isquite unstable, and the shocked pulsar wind can strongly change under small perturbations.
Conclusions.
We confirm the sensitive nature of the two-wind interaction structure to perturbations, in particular when the stellar windis inhomogeneous. For realistic over-dense regions of the stellar wind, the interaction region may shrink by a factor of a few, with theshocked flow presenting a complex spatial and temporal pattern. This can lead to strong variations in the non-thermal radiation.
Key words. hydrodynamics – X-rays: binaries – stars: winds, outflows – radiation mechanisms: non-thermal – gamma rays: stars
1. Introduction
In binaries consisting of a massive star and a young non-accreting pulsar, the relativistic wind of the pulsar interacts withthe non-relativistic wind of the stellar companion. This can re-sult in e ffi cient particle acceleration and in the production ofnon-thermal radiation, from radio to gamma rays. At present,PSR B1259 −
63 has been confirmed as a member of this class ofobjects, and there are several other high-mass binaries hosting acompact object whose nature is still unknown and may belong tothis class as well (e.g., Dubus 2013; Paredes et al. 2013).The properties of the non-thermal radiation in high-mass bi-naries hosting a non-accreting pulsar are determined to a largeextent by the dynamics of the two-wind interaction structure. Tostudy di ff erent aspects of this structure and its evolution alongthe orbit, heavy and detailed numerical simulations have beenperformed. Romero et al. (2007) carried out non-relativistic,smoothed particle hydrodynamic (SPH) simulations in three di-mensions (3D) to study the orbital evolution of a hypothetic two-wind interaction region in LS I +
61 303. Bogovalov et al. (2008)studied numerically the collision of the pulsar and the stellarwind in PSR B1259 −
63, treating the flows as laminar fluids byrelativistic hydrodynamical simulations in two dimensions (2D).For the same system, Okazaki et al. (2011) and Takata et al.(2012) performed 3D SPH non-relativistic simulations to studythe tidal and wind interactions between the pulsar and the decre- (cid:63)
Serra Húnter Fellow. tion disc of the Be star, and Bogovalov et al. (2012) studied thecollision of a magnetized anisotropic pulsar wind and the stellarwind adopting the same geometry and using a similar treatmentas adopted in Bogovalov et al. (2008). Lamberts et al. (2012a)performed non-relativistic hydrodynamical simulations in planargeometry of the two-wind interaction structure for several or-bits on large scales. Lamberts et al. (2012b) and Lamberts et al.(2013) performed 2D relativistic hydrodynamical (RHD) simu-lations in planar geometry of the interaction between the pulsarand the stellar wind on the scales of the binary system. Finally,Bosch-Ramon et al. (2012) performed 2D RHD simulations inplanar geometry aimed at studying the interaction between thestellar and pulsar winds on scales at which the orbital motion isimportant.Some relevant conclusions of the numerical work done upto now are the importance of instabilities, in particular theRayleigh-Taylor (RT) and Kelvin-Helmholtz (KH) instabili-ties, in the spatial and temporal properties of the two-wind in-teraction structure on both small and large scales, the occur-rence of shocked-flow re-acceleration, the impact of orbital mo-tion through lateral, strong shocks, the e ff ective mixing of theshocked winds downstream of the flow, the strong e ff ects of thepulsar wind ram pressure on Be discs, and the relatively mi- RT: instability developed in the contact surface between two fluids ofdi ff erent densities, with the lightest fluid exerting a force on the densestone. KH: instability developed in the contact surface between two fluidsof di ff erent parallel velocities (see Chandrasekhar 1961).Article number, page 1 of 32 a r X i v : . [ a s t r o - ph . H E ] F e b & A proofs: manuscript no. clump nor role of pulsar wind anisotropies and magnetic fields. How-ever, all the simulations of binaries hosting a pulsar consideredsmooth winds, while the stellar wind is expected to be inhomo-geneous, in particular for earlier spectral types (Lucy & Solomon1970). The wind inhomogeneities (clumps hereafter) are ex-pected to play an important role in the structure of the two-windinteraction region, and therefore in the non-thermal emission.The origin of these clumps can be explained by di ff erent mech-anisms depending on their scale, either the stellar surface at thewind formation region, or a circumstellar disc in the Be system.For instance, the impact of a fragment of a Be star disc on thetwo-wind interaction region has been proposed to cause the GeVflare observed in the gamma-ray binary PSR B1259 −
63 (Abdoet al. 2011; Chernyakova et al. 2014). In addition, some sort ofwind clumping could explain short flares on scales of secondsto hours found in the X-ray light curves of LS 5039 and LS I +
61 303 (e.g., Bosch-Ramon et al. 2005; Paredes et al. 2007;Smith et al. 2009; Li et al. 2011). The clump-pulsar wind in-teraction problem has previously been studied analytically (seeBosch-Ramon 2013), but numerical calculations of the interac-tion between an inhomogeneous stellar wind and a relativisticpulsar wind have not been conducted yet.In this work, we present 2D axisymmetric RHD simulationsof the interaction between a wind clump and a relativistic pulsarwind. Orbital motion has not been accounted for because of theaxisymmetric nature of the simulations; given the small spatialand temporal scales considered, this is a reasonable approxima-tion. We simulated clumps with di ff erent sizes and densities tostudy the global variations induced in the two-wind interactionstructure, as well as on the evolution of the clumps under theimpact of the pulsar wind.We note that this is the first time that relativistic simulationsare carried out in a more realistic set-up, with axisymmetry alongthe two-star line, without suppressing the development of insta-bilities. The numerical results are compared with the analyticalresults found in Bosch-Ramon (2013). The paper is structured asfollows: in Sect. 2 we describe the physical framework: the two-wind interaction region, the origin of the clumps in the windsof massive stars, and the analytical equations; in Sect. 3, we de-scribe the numerical simulations: the numerical set-up and theobtained results; finally, a discussion and the conclusions can befound in Sects. 4 and 5.
2. Physical framework
The collision between the stellar wind and the wind from theyoung pulsar creates a shock structure bow-shaped towards thewind with the lower momentum flux. The contact discontinuityis located where the ram pressure of the two winds balances. Asdiscussed in Bosch-Ramon (2013), the presence of clumps cansignificantly distort the overall interaction structure, althoughtheir dynamical impact is determined by their size and mass (see,e.g., Pittard 2007, for a similar non-relativistic scenario). Windsformed by small and light clumps will behave as uniform winds,while massive clumps will tend to come less frequently and bemore damaging for the global stability of the interaction region(see Perucho & Bosch-Ramon 2012 for the e ff ect of an inhomo-geneous wind on a high-mass microquasar jet). A sketch of thephysical scenario is shown in Fig. 1.The instability of the line-driving mechanism in the innerwind is thought to be an important source of clumps in the windstructure of hot stars of spectral types O and B (Lucy & Solomon1970). Runacres & Owocki (2002) suggested that the outer evo-lution of the inhomogeneous wind can be approximated as a pure Fig. 1.
Cartoon representing a density map of the physical scenario: aninhomogeneous stellar wind interacts with the relativistic pulsar wind.The shocked stellar wind and the shocked pulsar wind are separated bya contact discontinuity. The distances are not to scale. The simulationspresented in this paper only considered a stellar wind with one clumpcentered at the axis joinning the two stars. gas dynamical problem, and the stellar wind clumps initiatedclose to the star can survive up to long distances ( d ∼ (cid:12) ).Mo ff at (2008) suggested that all hot stellar winds are likely to beinhomogeneous because of radiative instabilities, with a multi-scale distribution of masses and sizes. Therefore, consideringthe outer wind evolution as a gas dynamical problem, the clumpscould expand with the wind flow and large clumps could reachthe interaction region between the stellar and the pulsar wind.However, the relatively small initial sizes ( ∼ . R ∗ ), and lim-itations on the clump growth (which can be parameterized as R c = R c ( d / R ∗ ) α , with α ≤ R c being the initial clump ra-dius) could prevent the formation of clumps as large as, R c ∼ R ∗ ,for instance, which is required for the strongest variations of thetwo-wind interaction structure (see Bosch-Ramon 2013, and ref-erences there in). Typical massive star radii are R ∗ ∼ R (cid:12) .Large clumps di ff erent from those related to radiative insta-bilities may be found in the stellar wind. For instance, early-type stars known as Be stars present a decretion disc formedby material ejected from the stellar equator by rapid rotation(Hanuschik 1996). The truncation of this disc, either caused bytidal forces or by direct pulsar wind ram pressure (see Okazakiet al. 2011), could explain the presence of large clumps in thestellar wind formed by chunks of disc. This possibility has beenconsidered for instance to explain the GeV flare detected fromthe gamma-ray binary PSR B1259 −
63 by Chernyakova et al.(2014) (see Sect. 4.2.1). Other types of large-scale structures inthe stellar wind, of size of the order of the stellar radius havealso been inferred from the observed discrete absorption compo-nents (DACs) in the UV (e.g., Kaper et al. 1997). The appearanceof DACs can be interpreted as co-rotating dense structures pro-
Article number, page 2 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system duced at the stellar surface and extending several tens of stellarradii (e.g., Lobel 2008), and can be attributed to the dynami-cal e ff ects of rotation, magnetic fields, or non-radial pulsations(Cranmer & Owocki 1996). The arrival of any of those structuresat the two-wind interaction region probably modify the latter inspace and time.Bosch-Ramon (2013), developed an analytical model for thetwo-wind interaction dynamics accounting for the lifetime ofclumps under the pulsar-wind impact. It was concluded that fora clump radius R c (cid:28) χ − / ∆ R , being ∆ R the thickness of thetwo-wind interaction region and χ = ρ c /ρ w the density con-trast between the clump density and the wind density, the clumpis destroyed and deflected by the shocked wind medium beforecrossing the two-wind interaction region. Otherwise, for R c ap-proaching χ − / ∆ R or larger, the clump will penetrate into theunshocked pulsar wind. When R and R (cid:48) are the initial and finaldistances between the pulsar and the contact discontinuity, fol-lowing Bosch-Ramon (2013), this approximate relation applies R (cid:48) ∼ R − χ / R c , (1)which for R c ∼ R (cid:48) implies R c ∼ R (cid:48) ∼ R (cid:0) + χ / (cid:1) (where R (cid:48) > . (2)The simulations presented here are intended to investigatein detail the dynamical consequences of the presence of clumpsin the stellar wind. Below we compare analytical estimates withnumerical results.
3. Numerical simulations
The simulations were performed using a finite-di ff erence codebased on a high-resolution shock-capturing scheme that solvesthe equations of relativistic hydrodynamics in two dimensionsin a conservation form (Martí et al. 1997). The code is paral-lelized using open message passing (OMP) (Perucho et al. 2005).The simulations were run in a workstation with two Intel(R)Xeon(R) CPU E5-2643 processors (3.30 GHz, 4 × γ = . r ∈ [0 , a ] and z ∈ [0 , a ], where a = × cm. Theadopted resolution is 150 ×
250 cells. The pulsar is located at( r , z ) = (0 , a ), and the star at ( r , z ) = (0 , a ), outside the simu-lated grid. This yields a two-star separation of d = . × cm.The typical orbital separation distances in gamma-ray binariesare ∼ cm (see, e.g., Dubus 2013, and references therein).The initial conditions of the simulation were computed inspherical coordinates assuming adiabatic gas radial propagationand solving the Bernoulli equation for the pulsar and stellarwinds. The initial two-wind separation point is derived as theapproximate location where the on-axis wind ram pressures be-come equal ( z = a ).The chosen physical parameters of the pulsar wind are thetotal luminosity L sd = erg s − , the Lorentz factor Γ =
5, andthe specific internal energy (cid:15) pw = . × erg g − . As a resultof resolution limitations, the adopted Lorentz factor is smaller Table 1.
Simulation parameters.
Parameter Value γ . ρ . × − g cm − a × cm l r a (2 . × cm) l z a (4 . × cm) n r n z Notes.
Polytropic index γ , density unit ρ , distance unit a , physical r grid size l r , physical z grid size l z , number of cells in the r axis n r , andnumber of cells in the z axis n z . than expected in pulsar winds, Γ ∼ –10 (see Khangulyanet al. 2012; Aharonian et al. 2012, and references therein), buthigh enough to capture important relativistic e ff ects (see the dis-cussion in Bosch-Ramon et al. 2012). For the stellar wind, thephysical parameters are a mass-loss rate ˙ M = − M (cid:12) yr − , aradial velocity v sw = − , and a specific internal en-ergy (cid:15) sw = . × erg g − . The derived pulsar wind and stel-lar wind densities are ρ pw = . × − g cm − and ρ sw = . × − g cm − , respectively. All the previous quantities,except for L sd and ˙ M , are given at a distance r = a with respectto the pulsar and star centres. The values of L sd , ˙ M and v sw werechosen as representative values of gamma-ray binaries hosting apulsar because L sd is to be high enough to power the gamma-rayemission, and the stellar wind properties correspond to those ofan OB star. All these parameters are summarized in Tables 1 and2. The previous physical values lead to a pulsar-to-stellar windthrust ratio of η = F pw S pw F sw S sw = (cid:16) ρ pw Γ v h pw + p pw (cid:17) S pw ˙ M v sw + p sw S sw ≈ L sd ˙ M v sw c ∼ . , (3)where F pw ( F sw ) is the momentum flux of the pulsar wind (stellarwind), S pw ( S sw ) the spherical surface at a distance r = a withrespect to the pulsar (star), h pw the specific enthalpy of the pulsarwind given by h pw = + (cid:15) pw c + p pw ρ pw c , and p pw / sw the pressure ofthe pulsar / stellar wind given by p pw / sw = ( γ − ρ pw / sw (cid:15) pw / sw .The pulsar wind is introduced by defining an injector withthe mentioned properties and radius 3 a (2 . × cm; 15 cells)as a boundary condition, and the stellar wind is injected accord-ing to its characterization at the upper boundary of the grid. Thelower and right boundaries of the grid are set to outflow, whilethe left boundary of the grid is set to reflection.The pulsar wind is injected with a Lorentz factor of 5, butbecause of adiabatic propagation, it reaches Lorentz factors of ∼ r , z ) = (0 , a )after the steady state is reached. The clumps are characterizedby their radius R c and their density contrast χ with respect to theaverage stellar wind value at their location. We present here theresults for four di ff erent cases: χ =
10 and R c = a , 2 . a , 5 a and χ =
30 and R c = a , corresponding to di ff erent degrees ofwind inhomogeneity, from modest χ =
10 and R c = a to ratherhigh χ =
10 (30) and R c = a .The axisymmetric nature of the simulations, in particular thereflection boundary conditions, the presence of a coordinate sin-gularity, plus the concentration of fluxes, all at r =
0, lead to thegeneration of perturbations that quickly grow, posing di ffi culties Article number, page 3 of 32 & A proofs: manuscript no. clump
Table 2.
Pulsar and stellar parameters in code and CGS units.
Parameter Pulsar wind Stellar wind ρ . ρ (1 . × − g cm − ) 1 . × ρ (2 . × − g cm − ) (cid:15) . c (9 × erg g − ) 2 . × − c (1 . × erg g − ) v . c (2 . × cm s − ) 0 . c (3 × cm s − )( r , z ) (0 , a ); (0 , × cm) (0 , a ); (0 , . × cm) Notes.
Density ρ , specific interal energy (cid:15) , and wind velocity v at a distance r = a with respect to the pulsar and star centres located at ( r , z ). to achieve a steady state. This e ff ect makes the simulation alsoquite sensitive to the initial conditions, and the initial two-windcontact discontinuity has to be carefully chosen so as not to en-hance the axis perturbations further. In this context, the modestresolution of most of the simulations and the pulsar wind Lorentzfactor adopted were chosen to stabilize the solution of the cal-culations through a larger numerical dissipation and moderatewind density contrast. Given the fast growth of any perturbation(numerical or physical) in the simulated problem (as explainedbelow), higher resolutions will enhance (as also shown below)the instability of the simulated structures through the penetra-tion of small fragments of stellar wind material into the pulsarwind region. On the other hand, higher Lorentz factors wouldhave enhanced the wind density contrast and thereby the insta-bility growth rate. Finally, the grid size was also limited to pre-vent the instabilities from fully disrupting the two-wind interac-tion structure within the simulated region. With our resolution,larger grid sizes would have given enough time to the instabil-ities to develop and grow towards the downstream edge of thegrid. The resulting structures would be similar to those found inthe higher resolution simulations also presented in this work, andwould also lead to the collapse of the two-wind interaction re-gion. Albeit small, the grid size does not significantly a ff ect thehydrodynamics within as the flow is supersonic at the outflowboundaries. We remark that because the clump is the strongestdynamical factor, its impact can be studied neglecting the role ofnumerical perturbations on the axis after a (quasi-)steady statehas been reached. Figure 2 shows the density map of the simulated flow in the(quasi-)steady state at t = ∼ a from the pulsar. The thickness of the shocked two-wind region is ∼ a on the simulation axis. After being shocked,winds are accelerated side-wards because of the strong pressuregradient in the downstream region, as pointed out in Bogov-alov et al. (2008). This is seen in the velocity maps presented inthis section, for which the two-wind interaction structure has notyet been a ff ected by the clump arrival. Because of this flow re-acceleration, the ratio of the momentum-flux to pressure showsthat the flow becomes supersonic at the boundaries.These axisymmetric simulations are quite sensitive to the ini-tial conditions, and in particular, reaching steady state, is highlydependent on the initial set-up of the simulation. It is still anopen question, to be answered with 3D simulations, to whichextent this sensitivity to the initial set-up is caused by the ge-ometry of the grid. Nevertheless, it is expected that instabilitiesquickly grow if perturbations are present, as shown by planarcoordinate relativistic simulations (Bosch-Ramon et al. 2012;Lamberts et al. 2013). In the present simulations, the boundaryconditions at the axis induce perturbations there that grow un- r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β Fig. 2.
Density distribution by colour at time t = ff erent locations ( β = v/ c ). Theaxes units are a = × cm. The density units are ρ = . × − g cm − . The pulsar and the star are located at ( r , z ) = (0 , a ) and( r , z ) = (0 , a ), respectively. The same units are used in all the figures. der the RT instability. These perturbations are then amplified bythe KH instability as the flow propagates along the contact dis-continuity. For grid sizes larger than 30 a in the r direction, thestructures resulting from instability growth can reach deep intothe pulsar wind, eventually disrupting the two-wind interactionregion and filling the pulsar vicinity with shocked material. Forsmall grids, on the other hand, these structures are advected outof the grid before it happens. The RT and KH instabilities growfaster for faster (and lighter) pulsar winds, which for a fixed η implies a stronger wind density contrast and thus a more un-stable contact discontinuity. A Lorentz factor of Γ ∼ r ∼ − a in the two-windinteraction region in Fig. 2, the opening angle and width of theshocked wind zone are similar to those found in Bogovalov et al. Article number, page 4 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system (2008, 2012). Those simulations were similar to ours given theiraxisymmetry, although the simulated flows were laminar.Figures 3, 4, 5, and 6 show the density map evolution forthe clumps characterized by χ =
10 and R c = a , 2 . a , 5 a ,and χ =
30 and R c = a , respectively. The times of the mapsnapshots were chosen such that the images are illustrative ofthe structure evolution. The bigger clumps, characterized by χ =
10 and R c = . a and R c = a , and the smallest but densestone, with χ =
30 and R c = a , strongly perturb the two-windinteraction region, pushing the contact discontinuity to less thana half its initial distance (see Fig. 2) to the pulsar (see Figs. 4, 5,and 6, respectively).We now focus on a clump with χ =
10 and R c = . a , asit is illustrative of a strong impact by a clump of intermediatesize. For this case, we also show Figs. 7, 8, 9, 10, and 11, whichdisplay the evolution of the clump in the density map, the tracermap (1 for the clump material and 0 otherwise), the pressuremap, the map of the momentum flux over the pressure (greaterthan a few for a super-sonic flow), and the map of β . These fig-ures show how the two-wind interaction region is pushed by theclump, until the contact discontinuity reaches a minimum dis-tance to the pulsar at R (cid:48) num ≈ . a ( R (cid:48) num ≈ . a for the termina-tion shock). After that, the clump starts to be pushed backwardsby the pulsar wind, is shocked, and decelerates, eventually dis-rupting, with its fragments driven away from the simulation axisby the shocked flow. All this is clearly seen in the density maps(Figs. 4 and 7), and in the tracer (Fig. 8). In addition, Figs. 9,10 and 11 provide information on the presence of shocks, ap-parent as sudden increases in pressure or drops in the ratio ofmomentum-flux to pressure, and relevant for particle accelera-tion; flow re-acceleration, important for non-radiative cooling,flow relativistic motion, and flow-to-sound-speed relation; andflow speed and direction downstream of the shocks, determin-ing Doppler e ff ects on the flow emission. Equivalent images toFigs. 7, 8, 9, 10, and 11 for the other three simulated clumpcases are available in the on-line material. In all the simulations,a t = χ =
10 and R c = a , the system recovers steady stateafter the clump has been advected. In contrast, for the cases with χ =
10 and R c = . a , and χ =
30 and R c = a , the instabil-ities eventually lead to the collapse of the two-wind interactionregion, filling the region close to the pulsar with shocked pul-sar wind. The clump with R c = a and χ =
10 only slightlyperturbs the two-wind interaction region and pushes the contactdiscontinuity to ∼ / R c = χ = ff erent combinations ofthese parameters yield results in line with those shown.To show the impact of increasing resolution, we present theresults of two simulations with the same set-up but, a resolu-tion 2 and 1.5 times (i.e. 300 ×
500 and 225 ×
375 cells, respec-tively) higher than the resolution adopted for most of the sim-ulations in the paper. For the highest resolution case, we wereunable to reach steady state (without clump) because the higher r [ a ]01020304050 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β Fig. 12.
Density distribution by colour at time t = ×
500 cells. The remaining plot propertiesare the same as those of Fig. 2. resolution allowed the rapid development of instabilities alreadywithin the grid, leading to the collapse of the two-wind inter-action region. A temporal state of this simulation taken shortlybefore the collapse of the interaction region is shown in Fig. 12.For the intermediate-resolution case, we reached (quasi-)steadystate, but compared with the lower resolution simulations, it pre-sented a more unstable pattern in the two-wind interaction regionaway from the simulation axis and a more compressed two-windinteraction region close to the axis singularity. A sequence of im-ages showing the evolution of a clump characterized by χ = R c = . a for the intermediate resolution simulation isshown in Fig.13. The higher resolution allowed the developmentof denser small-scale structures that pushed the pulsar wind ter-mination shock closer to the neutron star. However, the generalbehaviour of the higher and lower resolution clump simulationsis similar, which suggests that the main features resulting fromthe lowest resolution simulations, albeit smoother, are reliable.
4. Discussion
The interaction of a relativistic pulsar wind with the non-relativistic wind of a massive star was simulated through rela-tivistic axisymmetric hydrodynamic calculations. The two-windinteraction structure reaches a (quasi-)steady state, although thesolution is metastable and quite sensitive to the initial set-upparameters, such as the grid size or the wind density con-trast. Therefore, even though the collision between the twowinds forms a structure similar to those previously found in ax-isymmetric relativistic simulations, it shows some irregularity(Fig. 2). In the dynamical problem solved here, the axisymmet-ric geometry may introduce numerical perturbations to the phys-
Article number, page 5 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 10 − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β Fig. 3.
Density distribution by colour for the case with clump parameters χ =
10 and R c = a for the times shown at the top of each plot. Theremaining plot properties are the same as those of Fig. 2.Article number, page 6 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 10 − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β Fig. 4.
Density distribution by colour for the case with clump parameters χ =
10 and R c = . a for the times shown at the top of each plot. Theremaining plot properties are the same as those of Fig. 2. Article number, page 7 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 10 − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β Fig. 5.
Density distribution by colour for the case with clump parameters χ =
10 and R c = a for the times shown at the top of each plot. Theremaining plot properties are the same as those of Fig. 2.Article number, page 8 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 10 − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β Fig. 6.
Density distribution by colour for the case with clump parameters χ =
30 and R c = a for the times shown at the top of each plot. Theremaining plot properties are the same as those of Fig. 2. Article number, page 9 of 32 & A proofs: manuscript no. clump r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = Fig. 7.
Zoom-in of the density distribution by colour for the case with clump parameters χ =
10 and R c = . a for the times shown at the top ofeach plot. The remaining plot properties are the same as those of Fig. 4. The pulsar and the star are located at ( r , z ) = (0 , a ) and ( r , z ) = (0 , a ),respectively. ical solution that grow a ff ected by the RT instability close to theaxis, and by the KH instability farther away. However, similarirregularities have been also found in relativistic 2D simulationsin planar geometry that have been attributed to the developmentof KH instabilities (Bosch-Ramon et al. 2012; Lamberts et al.2013), which implies that the irregularities found in this workare not only related to the axisymmetric geometry adopted.The growth of instabilities strongly a ff ects the shockedpulsar-wind region in its (quasi-)steady state. Given the veryhigh sound speed, this region can change within time intervalssignificantly shorter than the dynamical time of the simulation, ∼ d /v sw , dominated by the slow stellar wind. Nevertheless, wefind that the interaction region globally follows the expected ge-ometry both for the contact discontinuity and for the width of thetwo-wind collision region (Bogovalov et al. 2008, 2012).The arrival of clumps can have a very strong impact on thewhole interaction structure, overcoming the e ff ect of any pos-sible numerical perturbation related to the geometry of the cal-culations. The clumps trigger violent RT and KH instabilitiesbecause the structure is prone to su ff er them, and thus the inter-action between the clump and the two-wind collision region ishighly non-linear, leading not only to a secular modification ofthe global geometry, but also to quick changes of the shockedpulsar-wind region. This is apparent in all the map time se-quences (in particular in the zoomed-in sequence: Fig. 7), as thelargest variations a ff ect the shocked pulsar wind.Below we discuss a few relevant points: the comparison be-tween the analytical and the numerical approximation; the clumpe ff ect on the global structure and radiation, with a mention ofthe GeV flare of PSR B1259 −
63 as a possible instance of matterclump-perturbation of the two-wind interaction region, and workunder development.
An analytical estimate of the minimum distance between thetwo-wind interaction region and the pulsar, after the clumphas penetrated the shocked pulsar wind, has been presented inSect. 2. For comparison, the analytical and numerical values ofthe minimum distances from the contact discontinuity, and thetermination shock, to the pulsar, are shown in Table 3 for all thesimulated clumps.The final distance between the pulsar and the contact dis-continuity, for a given density contrast and a clump radius, com-puted using Eq. (1), agree well with the numerical results ex-cept for the clump with radius R c = a (see Table 3). Thereason is that such a large clump is already outside of the ap-plication range of Eq. (1), which strictly applies only to clumpswith R c (cid:28) R (cid:48) . For R c → R (cid:48) or bigger, the clump behaves moreas an homogeneous wind than as a discrete obstacle, but χ timesdenser than the average stellar wind. In this case, the minimumdistance from the contact discontinuity to the pulsar can be com-puted from Eq. (3) and assuming momentum flux equality as R (cid:48) = η / χ d / (1 + η / χ ) (cid:39) η / χ d, with η χ = χ − η . For χ =
10 itgives R (cid:48) an ∼ . a , which agrees well with the numerical resultfor the clump with R c = a ( R (cid:48) num ≈ . a ).The time required for the clump with χ =
10 and R c = . a to become shocked is ∼ ∼ ∼ χ / R c /v sw ≈ Article number, page 10 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 0 . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β Fig. 8.
Tracer distribution by colour for the case with clump parameters χ =
10 and R c = . a for the times shown at the top of each plot. Thetracer value ranges from 0 (pulsar and stellar wind) to 1 (clump). The remaining plot properties are the same as those of Fig. 2.Article number, page 11 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 10 − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β Fig. 9.
Pressure distribution in units of ρ c by colour for the case with clump parameters χ =
10 and R c = . a for the times shown at the top ofeach plot. The remaining plot properties are the same as those of Fig. 2.Article number, page 12 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 10 F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β Fig. 10.
Momentum flux over pressure distribution by colour for the case with clump parameters χ =
10 and R c = . a for the times shown at thetop of each plot. The momentum flux is given by F p = ρ Γ v (1 + (cid:15) + P /ρ ) + P . The remaining plot properties are the same as those of Fig. 2.Article number, page 13 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 0 . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β Fig. 11. β distribution by colour for the case with clump parameters χ =
10 and R c = . a for the times shown at the top of each plot. Theremaining plot properties are the same as those of Fig. 2.Article number, page 14 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]01020304050 z [ a ] t =
270 s 10 − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]01020304050 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]01020304050 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]01020304050 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]01020304050 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β r [ a ]01020304050 z [ a ] t = − D e n s it y [ ρ ] . . . . . . . . . β Fig. 13.
Density distribution by colour for the case with clump parameters χ =
10 and R c = . a and a resolution of 225 ×
375 cells for the timesshown at the top of each plot. The remaining plot properties are the same as those of Fig. 2. Article number, page 15 of 32 & A proofs: manuscript no. clump
Table 3.
Clump impact on the size of the two-wind interaction region. χ R c (a) R (cid:48) num (a) R (cid:48) an (a) R (cid:48) (TS)num ( a )10 1 13 14.8 1110 2.5 8.5 10 7.510 5 7.5 2.2 630 1 10 12.5 8 Notes.
Density contrast χ , clump radius R c , numerical minimum dis-tance with respect to the contact discontinuity R (cid:48) num (analytical: R (cid:48) an ),and numerical minimum distance with respect to the termination shock R (cid:48) (TS)num . On the scales simulated in this work, stellar winds with a mod-est inhomogeneity degree ( χ =
10 and R c = a ) presentnon-negligible variations of the interaction structure, and evensmaller / lighter clumps can provide continuous perturbations forthe development of instabilities in the contact discontinuity. Inaddition, medium-sized or denser clumps ( R c = .
5, 5 a and χ = R c = a and χ =
30) can lead to strong variationsin the size of the two-wind interaction structure. Both small andmedium-sized / denser clumps can generate quick and global vari-ations in the shocked pulsar wind. This a ff ects the location of thepulsar wind termination shock (i) and also introduces seeds forsmall-scale relativistic and transonic turbulence that would growdownstream of this shock (ii).(i) The consequences of large variations in the terminationshock location have previously been discussed in Bosch-Ramon(2013). In short, they can induce variations in the cooling, ra-diative as well as non-radiative, channels, and non-linear radi-ation processes, through synchrotron self-Compton or internalpair creation, for large reductions of the emitter size.(ii) The relativistic flow variations on small spatial and tem-poral scales downstream of the pulsar wind shock, already ap-parent despite the modest resolution in Fig. 11, would lead toa complex radiative pattern in time and direction. This complexradiative pattern is caused by Doppler boosting because of thecomplex orientation of the fluid lines and the relativistic speedsachieved through re-acceleration of the shocked pulsar wind (seeKhangulyan et al. 2014, for a study of the impact of Dopplerboosting on radiation).Finally, weak shocks are present in the shocked pulsar wind(e.g., Fig. 10), which suggests that further particle acceleration,additional to that occurring in the pulsar wind termination shock,could take place already deep inside the binary system, wellbefore the postshock flow has been a ff ected by the orbital mo-tion (see, e.g., Bosch-Ramon et al. 2012, for a discussion of thelarger-scale evolution of the shocked structure). − The flare in PSR B1259 −
63, observed by Fermi about 30 daysafter periastron passage (Abdo et al. 2011), has a potential con-nection with the Be disc through the disruption and fragmenta-tion of the latter. This may have led to the impact of a densepiece of disc on the two-wind interaction structure, largely re-ducing the size of the pulsar wind termination shock. This re-duction could allow e ffi cient Compton scattering by a populationof GeV electrons on local X-ray photons, as proposed in Dubus& Cerutti (2013) (see also Khangulyan et al. 2012, for a simi-lar proposal involving infrared photons) as a result of the strongenhancement of the X-ray photon density. This idea is worth to be studied in more detail, although it remains unclear why sucha modification of the GeV emitter has no clear e ff ects at otherwavelengths (Chernyakova et al. 2014).It is worth estimating what fraction of the mass of a Be stellardisc the simulated clumps would represent. The mass of a typicalBe disc can be derived from the disc mass-loss rate, ∼ − –10 − M (cid:12) yr − , times the disc extension, say ∼ ∼ − (Okazaki2001). This yields a disc mass of ∼ –10 g. In particular, forPSR B1259 −
63, Chernyakova et al. (2014) estimated a disc massof 2 × g, with a disc radius of 0.42 AU, roughly similar tothe values just mentioned. The masses of the simulated clumpsare ∼ × g, ∼ g, ∼ g and ∼ × g for theclumps characterized by χ =
10 and R c = a , 2 . a , 5 a and χ =
30 and R c = a , respectively. For comparison, the mass ofthe clump with R c = . a would be ∼ . To distinguish the importance of numerical artefacts in our re-sults and study a more realistic case, a 3D version of the simu-lations presented here is under way to more accurately charac-terize the instabilities that a ff ect the two-wind interaction region.In addition, we are planning to carry out calculations of the ra-diation outcome expected from the two-wind interaction region,and most importantly, from the clump interaction with this struc-ture, making full use of the dynamical information provided bythese simulations.
5. Conclusions
We presented, for the first time, 2D axisymmetric RHD simula-tions of the interaction between an inhomogeneous stellar windand a relativistic pulsar wind, focusing on the region inside thebinary system. We simulated clumps with di ff erent sizes anddensities to study di ff erent degrees of the stellar wind inhomo-geneity. The presence of the clumps results in significant varia-tions of the interaction region, which are expected to stronglya ff ect the non-thermal radiation as well. Therefore, we con-firm the sensitive nature of two-wind interaction structure underthe presence of the stellar wind inhomogeneities. The shockedflow presents a complex spatial and temporal pattern, with fastchanges in the shocked pulsar wind. This can lead to strongshort time-scale flux variability in the non-thermal radiation ofgamma-ray binaries, which might be observed for instance ingamma rays with the future Cherenkov Telescope Array (CTA;Acharya et al. 2013; Paredes et al. 2013). Acknowledgements.
We thank the anonymous referee for his / her constructiveand useful comments. We acknowledge support by the Spanish Ministerio deEconomía y Competitividad (MINECO) under grants AYA2010-21782-C03-01, AYA2010-21322-C03-01, AYA2010-21097-C03-01, AYA2013-47447-C3-1-P, FPA2010-22056-C06-02 and FPA2013-48381-C6-6-P. We also acknowledgesupport by the “Generalitat Valenciana” grant “PROMETEO-2009-103”. Thisresearch has been supported by the Marie Curie Career Integration Grant 321520.X.P.-F. also acknowledges financial support from Universitat de Barcelona andGeneralitat de Catalunya under grants APIF and FI, respectively. V.B-R. alsoacknowledges financial support from MINECO and European Social Fundsthrough a Ramón y Cajal fellowship. Article number, page 16 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system
Appendix A: Physical quantity maps for differentclump parameters
The density zoom, tracer, pressure, ratio of momentum-flux topressure, and velocity maps are presented in Figs. A.1–A.15 fordi ff erent clump parameters: χ =
10 and R c = a ; χ =
10 and R c = a ; and χ =
30 and R c = a . References
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJ, 736, L11Acharya, B. S., Actis, M., Aghajani, T., et al. 2013, Astroparticle Physics, 43, 3Aharonian, F. A., Bogovalov, S. V., & Khangulyan, D. 2012, Nature, 482, 507Bogovalov, S. V., Khangulyan, D., Koldoba, A. V., Ustyugova, G. V., & Aharo-nian, F. A. 2012, MNRAS, 419, 3426Bogovalov, S. V., Khangulyan, D. V., Koldoba, A. V., Ustyugova, G. V., & Aha-ronian, F. A. 2008, MNRAS, 387, 63Bosch-Ramon, V. 2013, A&A, 560, A32Bosch-Ramon, V., Barkov, M. V., Khangulyan, D., & Perucho, M. 2012, A&A,544, A59Bosch-Ramon, V., Paredes, J. M., Ribó, M., et al. 2005, ApJ, 628, 388Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stabilityChernyakova, M., Abdo, A. A., Neronov, A., et al. 2014, MNRAS, 439, 432Cranmer, S. R. & Owocki, S. P. 1996, ApJ, 462, 469Dubus, G. 2013, A&A Rev., 21, 64Dubus, G. & Cerutti, B. 2013, A&A, 557, A127Hanuschik, R. W. 1996, A&A, 308, 170Kaper, L., Henrichs, H. F., Fullerton, A. W., et al. 1997, A&A, 327, 281Khangulyan, D., Aharonian, F. A., Bogovalov, S. V., & Ribó, M. 2012, ApJ, 752,L17Khangulyan, D., Bogovalov, S. V., & Aharonian, F. A. 2014, International Jour-nal of Modern Physics Conference Series, 28, 60169Lamberts, A., Dubus, G., Lesur, G., & Fromang, S. 2012a, A&A, 546, A60Lamberts, A., Dubus, G., Fromang, S., & Lesur, G. 2012b, in American Instituteof Physics Conference Series, Vol. 1505, American Institute of Physics Con-ference Series, ed. F. A. Aharonian, W. Hofmann, & F. M. Rieger, 406–409Lamberts, A., Fromang, S., Dubus, G., & Teyssier, R. 2013, A&A, 560, A79Li, J., Torres, D. F., Zhang, S., et al. 2011, ApJ, 733, 89Lobel, A. 2008, in Clumping in Hot-Star Winds, ed. W.-R. Hamann, A. Feld-meier, & L. M. Oskinova, 81Lucy, L. B. & Solomon, P. M. 1970, ApJ, 159, 879Martí, J. M., Müller, E., Font, J. A., Ibáñez, J. M., & Marquina, A. 1997, ApJ,479, 151Mo ff at, A. F. J. 2008, in Clumping in Hot-Star Winds, ed. W.-R. Hamann,A. Feldmeier, & L. M. Oskinova, 17Okazaki, A. T. 2001, PASJ, 53, 119Okazaki, A. T., Nagataki, S., Naito, T., et al. 2011, PASJ, 63, 893Paredes, J. M., Bednarek, W., Bordas, P., et al. 2013, Astroparticle Physics, 43,301Paredes, J. M., Ribó, M., Bosch-Ramon, V., et al. 2007, ApJ, 664, L39Perucho, M. & Bosch-Ramon, V. 2012, A&A, 539, A57Perucho, M., Martí, J. M., & Hanasz, M. 2005, A&A, 443, 863Pittard, J. M. 2007, ApJ, 660, L141Romero, G. E., Okazaki, A. T., Orellana, M., & Owocki, S. P. 2007, A&A, 474,15Runacres, M. C. & Owocki, S. P. 2002, A&A, 381, 1015Smith, A., Kaaret, P., Holder, J., et al. 2009, ApJ, 693, 1621Takata, J., Okazaki, A. T., Nagataki, S., et al. 2012, ApJ, 750, 70 Article number, page 17 of 32 & A proofs: manuscript no. clump r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = Fig. A.1.
Zoom-in of the density distribution by colour for the case with clump parameters χ =
10 and R c = a for the times shown at the top ofeach plot. The coloured arrows represent the three-velocity at di ff erent locations ( β = v/ c ). The axes units are a = × cm. The pulsar and thestar are located at ( r , z ) = (0 , a ) and ( r , z ) = (0 , a ), respectively.Article number, page 18 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 0 . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β Fig. A.2.
Tracer distribution by colour for the case with clump parameters χ =
10 and R c = a for the times shown at the top of each plot. Thetracer value ranges from 0 (pulsar and stellar wind) to 1 (clump). The remaining plot properties are the same as those of Fig. A.1.Article number, page 19 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 10 − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β Fig. A.3.
Pressure distribution in units of ρ c by colour for the case with clump parameters χ =
10 and R c = a for the times shown at the top ofeach plot. The remaining plot properties are the same as those of Fig. A.1.Article number, page 20 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 10 F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β Fig. A.4.
Momentum flux over pressure distribution by colour for the case with clump parameters χ =
10 and R c = a for the times shown at thetop of each plot. The momentum flux is given by F p = ρ Γ v (1 + (cid:15) + P /ρ ) + P . The remaining plot properties are the same as those of Fig. A.1.Article number, page 21 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 0 . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β Fig. A.5. β distribution by colour for the case with clump parameters χ =
10 and R c = a for the times shown at the top of each plot. Theremaining plot properties are the same as those of Fig. A.1.Article number, page 22 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = Fig. A.6.
Zoom-in of the density distribution by colour for the case with clump parameters χ =
10 and R c = a for the times shown at the top ofeach plot. The remaining plot properties are the same as those of Fig. A.1. Article number, page 23 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 0 . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β Fig. A.7.
Tracer distribution by colour for the case with clump parameters χ =
10 and R c = a for the times shown at the top of each plot. Thetracer value ranges from 0 (pulsar and stellar wind) to 1 (clump). The remaining plot properties are the same as those of Fig. A.1.Article number, page 24 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 10 − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β Fig. A.8.
Pressure distribution in units of ρ c by colour for the case with clump parameters χ =
10 and R c = a for the times shown at the top ofeach plot. The remaining plot properties are the same as those of Fig. A.1. Article number, page 25 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 10 F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β Fig. A.9.
Momentum flux over pressure distribution by colour for the case with clump parameters χ =
10 and R c = a for the times shown at thetop of each plot. The momentum flux is given by F p = ρ Γ v (1 + (cid:15) + P /ρ ) + P . The remaining plot properties are the same as those of Fig. A.1.Article number, page 26 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 0 . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β Fig. A.10. β distribution by colour for the case with clump parameters χ =
10 and R c = aa
10 and R c = aa for the times shown at the top of each plot. Theremaining plot properties are the same as those of Fig. A.1. Article number, page 27 of 32 & A proofs: manuscript no. clump r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = r [ a ]101520253035 z [ a ] t = Fig. A.11.
Zoom-in of the density distribution by colour for the case with clump parameters χ =
30 and R c = a for the times shown at the topof each plot. The remaining plot properties are the same as those of Fig. A.1.Article number, page 28 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 0 . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . T r ace r . . . . . . . . . β Fig. A.12.
Tracer distribution by colour for the case with clump parameters χ =
30 and R c = a for the times shown at the top of each plot. Thetracer value ranges from 0 (pulsar and stellar wind) to 1 (clump). The remaining plot properties are the same as those of Fig. A.1.Article number, page 29 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 10 − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β r [ a ]010203040 z [ a ] t = − − − − P r e ss u r e [ ρ c ] . . . . . . . . . β Fig. A.13.
Pressure distribution in units of ρ c by colour for the case with clump parameters χ =
30 and R c = a for the times shown at the topof each plot. The remaining plot properties are the same as those of Fig. A.1.Article number, page 30 of 32. Paredes-Fortuny et al.: Simulations of an inhomogeneous stellar wind interacting with a pulsar wind in a binary system r [ a ]010203040 z [ a ] t =
670 s 10 F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β r [ a ]010203040 z [ a ] t = F p / P . . . . . . . . . β Fig. A.14.
Momentum flux over pressure distribution by colour for the case with clump parameters χ =
30 and R c = a for the times shown at thetop of each plot. The momentum flux is given by F p = ρ Γ v (1 + (cid:15) + P /ρ ) + P . The remaining plot properties are the same as those of Fig. A.1.Article number, page 31 of 32 & A proofs: manuscript no. clump r [ a ]010203040 z [ a ] t =
670 s 0 . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β r [ a ]010203040 z [ a ] t = . . . . . . . . . . . β Fig. A.15. β distribution by colour for the case with clump parameters χ =
30 and R c = aa