© ESO, 2011

1. Introduction

High-mass X-ray binaries (HMXB) are the parent population of high-mass microquasars (HMMQ), which show the peculiarity of presenting radio jet activity (e.g. Mirabel & Rodrígues 1999; Ribó 2005). Jets of HMMQs propagate in a medium characterized by the stellar wind of the massive primary star of the system. On the scales of the binary, the stellar wind can significantly affect the dynamics of the jet, trigger shocks suitable for particle acceleration, and under some conditions, disrupt the jet by a combination of asymmetric recollimation shocks and hydrodynamical instabilities (see Perucho & Bosch-Ramon 2008; and Perucho et al. 2010).

When launched, powerful enough jets can traverse the companion wind, and once at distances from the injection point that are much longer than the separation distance (z ≫ d_orb, where z is the axial coordinate and distance to the compact object), they do not suffer significant interference by the supersonic wind, which comes from behind the jet. At this stage, the jet head propagates with a speed  ~$\sqrt{{\mathit{n}}_{\mathrm{j}}\mathrm{\left(}{\mathit{d}}_{\mathrm{orb}}\mathrm{\right)}\mathit{/}{\mathit{n}}_{\mathrm{w}}\mathrm{\left(}{\mathit{d}}_{\mathrm{orb}}\mathrm{\right)}}\hspace{0.17em}{\mathit{v}}_{\mathrm{j}}\mathrm{~}\mathrm{0.1}\hspace{0.17em}{\mathit{v}}_{\mathrm{j}}$ assuming free expansion, where v_j is the initial jet velocity. In this propagation phase, a strong shock is present at the tip of the jet, and particle acceleration may tap a significant fraction of the jet luminosity there.

When the stellar wind is dynamically negligible, the jet can propagate as far and fast as the inertia of the swept ISM permits and, eventually, the jet head is shocked, thereby inflating a wide cocoon that pushes a shell of shocked ISM (Bordas et al. 2009), as in FRII galaxies (e.g. Fanaroff & Riley 1974; Kaiser & Alexander 1997). However, for systems hosting OB stars, the dynamics of the jet may be severely affected before reaching the unperturbed ISM.

To better understand the hydrodynamical jet evolution at different stages of an HMXB, we performed 2D simulations of HMMQ jets for three different scenarios: soon after the SNR explosion (case 1), the phase of wind/ISM interaction (case 2), and the phase in which the wind/ISM interaction is affected by the system proper motion (case 3). We also estimated the broadband nonthermal emission that may be produced in cases 1 and 2, and predicted their main observational features.

2. The physical picture

For a young HMXB, say of age t_src ~ 3 × 10⁴ yr, with typical ISM densities n_ISM ~ 1 cm^-1 and moderate peculiar velocities v_p ≲ 10⁷ cm s^-1, the system may still be embedded inside an adiabatically expanding SNR (case 1). In such a situation, if a jet is launched, after  ~10 yr it will encounter a region of shocked stellar wind confined by the shocked SNR ejecta, which in turn drive a strong shock in the ISM. Once in the high-pressure region of the shocked stellar wind, the jet head will be slowed down to medium subsonic velocities (~1000 km s^-1 for typical jet and wind properties, see below). At farther distances, the jet head will encounter the shocked SNR ejecta, eventually reaching the shocked ISM. Finally, the jet will enter the normal ISM having lost significant amounts of thrust, producing a shock in the ISM and possibly suffering disruption. Intermittent jet activity with a low duty cycle will not clean the way for later ejections, since the sound speed in the shocked wind and the SNR ejecta is high. High enough duty cycles ( ≳ 10%) for the jet to keep the way clean can be understood in terms of an averaged persistent activity. The case of an HMMQ persistent jet interacting with an SNR has a clear example in SS 433 (see Velázquez & Raga 2000; Zavala et al. 2008).

Older HMXBs have already left the SNR, and the stellar winds can directly interact with the ISM (case 2). In this case, the stellar wind inflates a bubble filled with shocked wind material, and this hot bubble drives a shock in the ISM that becomes radiative even for moderate ages of the source (t_src ~ 10⁵ yr; see Castor et al. 1975). Therefore, a jet launched in this HMXB phase will interact first with the shocked wind bubble, later reaching a thin, dense, and cold shell of shocked ISM. Qualitatively, the situation is similar to the previous case, but now the ISM shell is much denser because of radiative cooling. Eventually, the jet can reach the unperturbed ISM, shocking it. As in case 1, the jet may eventually be disrupted. Cygnus X-1 could correspond to the situation represented by case 2 (Martí et al. 1996; Gallo et al. 2005; Russell et al. 2007). Since there are no clear signatures of a reverse shock (Bordas et al. 2009) in this source and only the shocked ISM thermal radiation is detected, it may be that the Cygnus X-1 jet head is fully disrupted. A shock driven in the ISM can still develop in such a conditions since the disrupted jet region has much higher pressure than the environment, as simulations of jets of FRI radio galaxies show (Perucho & Martí 2007; see also Bordas et al. 2011).

If the HMXB age is old enough, say t_src ~ 10⁶ yr, the wind-driven ISM shock can be slower than v_p. Then, a bow-shaped shock forms, driven by the system proper motion and pointing in its direction (case 3). In this scenario, the shocked ISM material is slowly evacuated in the opposite direction to the system motion, and the shocked wind material is also advected in the same direction but at much faster velocities, similar to the shocked wind sound speed (~10³ km s^-1). If the jet is launched in this stage, for jet directions that are misaligned enough with the binary motion the shocked wind material will impact the jet roughly from the side. At farther distances, the jet will encounter the shocked ISM, roughly as dense and cold as in case 2. The dense shocked ISM layer plus lateral impact of the shocked wind will stop the advance of the jet and lead to total jet destruction before reaching the unshocked ISM.

Sketches of the scenarios corresponding to cases 1, 2, and 3 are presented in the top, middle, and bottom panels of Fig. 1.

Fig. 1 Sketches of a jet propagating through different media depending on the age of the system: the shocked stellar wind, the shocked SNR ejecta, and the shocked and unshocked ISM (top); the shocked stellar wind, and the shocked and unshocked ISM (middle); through the unshocked and shocked stellar wind, and interacting with the shocked ISM under a significant peculiar velocity (bottom).

3. Studying the jet termination

3.1. Medium characteristics

We have adopted an analytical approach for characterizing the medium in which the jet propagates. In the three cases considered here, the ISM density has been fixed to n_ISM = 1 cm^-3, with a temperature T_ISM = 10⁴ K, as expected in regions close to massive stars. The stellar wind velocity and the stellar mass-loss rate have been taken to be v_w = 2 × 10⁸ cm s^-1 and Ṁ_w = 10^-6   M_⊙/yr, respectively (typical of an O star; see, e.g., Puls et al. 2009), which corresponds to a kinetic luminosity of the wind ${\mathit{L}}_{\mathrm{w}}\mathrm{=}\mathit{Ṁ}{}_{\mathrm{w}}\hspace{0.17em}{\mathit{v}}_{\mathrm{w}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{\approx }\mathrm{1.3}\mathrm{\times }{\mathrm{10}}^{\mathrm{36}}$ erg s^-1.

In case 1, the system is assumed to be at the center of an SNR. The ejected mass of the SNR has been taken to be 2 M_⊙, and the explosion energy 5 × 10⁴⁹ erg. These moderate values have been adopted accounting for the fact that the energetics of supernovae leading to the formation of some HMXBs can be quite low (a relatively extreme case may be Cygnus X-1, see Mirabel & Rodrigues 2003). We have calculated the properties of the shocked ejecta of the SNR and the shocked ISM by adopting the solution of a spherical explosion in its adiabatic phase (e.g. Blandford & McKee 1976) for t_src ~ 3 × 10⁴ yr. The radius of the whole region affected by the SNR is ${\mathit{z}}_{\mathrm{bs}}\mathrm{\approx }\mathrm{(}{\mathit{E}}_{\mathrm{SNR}}\mathit{/}{\mathit{m}}_{\mathrm{H}}{\mathit{n}}_{\mathrm{ISM}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{5}}\hspace{0.17em}{\mathit{t}}_{\mathrm{src}}^{\mathrm{2}\mathit{/}\mathrm{5}}\mathrm{\approx }\mathrm{3}\mathrm{\times }{\mathrm{10}}^{\mathrm{19}}$ cm for the adopted parameter values, where m_H is the hydrogen mass. The density of the shocked ISM is the one for a strong adiabatic shock, n_ISMs = 4   n_ISM, and the temperature, T_ISMs ≈ 4 × 10⁵ K. From z_bs and E_SNR, we derive a characteristic pressure within that radius, P_SNR ~ 2   E_SNR/3   V_bs ≈ 3 × 10^-10 erg cm^-3 (T_SNR ≈ 10⁸ K; ${\mathit{v}}_{\mathrm{bs}}\mathrm{=}\mathrm{4}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{z}}_{\mathrm{bs}}^{\mathrm{3}}\mathit{/}\mathrm{3}$), as well as the locations of the contact discontinuity between the shocked ISM and the shocked ejecta, z_c ≈ 0.9   z_bs. The location of the contact discontinuity between the shocked stellar wind and the SNR ejecta is computed by taking the same pressure on both sides, obtaining z_ws ≈ 8 × 10¹⁸ cm. The shocked wind density is n_ws ≈ 0.02 m^-3, and the temperature T_ws ≈ 10⁸ K. We simulate the jet once it is already inside the shocked wind region, since the dynamical impact of the medium on the jet is negligible at closer distances.

In case 2, we have adopted the solution given by Castor et al. (1975) that describes the interaction of a stellar wind with the ISM at the stage when the shock in the latter is already radiative. We adopted an age of 10⁵ yr in this case, and the radius of the whole interaction region is z_bs = 0.76   (L_w/   m_H   n_ISM)^1/5   t^3/5 ≈ 2 × 10¹⁹ cm. The pressure in the shocked wind region is ${\mathit{P}}_{\mathrm{ws}}\mathrm{\approx }{\mathit{n}}_{\mathrm{ISM}}\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}\hspace{0.17em}\mathit{ż}\begin{array}{}\mathrm{2}\\ \mathrm{bs}\end{array}\mathrm{\approx }\mathrm{3}\mathrm{\times }{\mathrm{10}}^{-11}$ erg cm^-3 (T_ws ≈ 6 × 10⁷ K), and the density n_ws = 4 × 10^-3 cm^-3. The shocked ISM region is a very thin, cool, and dense layer. For simplicity, the width of this region has been taken equal to one simulation cell: 10¹⁷ cm, with n_ISMs ≈ 10² cm^-3 and T_ISMs = 2 × 10³ K. As in case 1, the jet injection point is already in the wind shocked region.

In case 3, taking v_p = 10⁷ cm s^-1 and t_src = 10⁶ yr, one gets ż_bs ~ 10⁶   cm   s^-1 ≪ v_p, which means that the wind-driven shock in the ISM becomes strongly asymmetric, with the shocked wind material advected backwards with respect to this shock. The wind driven shock in the ISM is now closer to the HMMQ than in case 2, at z_ws = 10¹⁹ cm. The shock in the stellar wind is at z_w ≈ 8 × 10¹⁸ cm, in which the (pre-shock) stellar wind density and temperature are n_w = 2.4 × 10^-4 cm^-3 and T_w = 10³ K. The density, pressure, and temperature in the shocked wind are n_ws = 10^-3 cm^-3, P_ws = 1.4 × 10^-11 erg cm^-3, and T_ws = 10⁸ K, respectively, and a lateral shocked wind velocity of v_ws ⊥  = 10⁸ cm s^-1 has been assumed. On the top of the shocked wind region, a thin shocked ISM shell of 10¹⁷ cm width, n_ISMs = 100 cm^-3, and T_ISMs = 10³ K is located. In this simulation, the jet is launched in the unshocked wind region, because it is more extended than in previous cases, and the unshocked wind is taken as coming from behind the jet. The values of the main ambient medium parameters, for the three cases considered here, are summarized in Table 1.

3.2. Simulations

We performed 2D numerical simulations of jets evolving in the three different scenarios described above. In cases 1 and 2 the jets are simulated using an axisymmetric grid in cylindrical coordinates, whereas the transversal motion of the shocked wind in case 3 requires a slab jet in planar coordinates. We expect from our previous experience (see Perucho & Bosch-Ramon 2008; Perucho et al. 2010) that the results of the slab simulation will be qualitatively correct. The jets are injected with power L_j = 3 × 10³⁶ erg s^-1 at z₀ = 3 × 10¹⁸   cm for cases 1 and 2, and z₀ = 1.25 × 10¹⁸   cm for case 3. In case 3, the power corresponds to a jet square section (i.e. the power crossing a surface equal to the jet width square). In cases 1 and 2, the jet is assumed to have a radius at the injection point R_j0 = 0.1   z₀. In case 3, the jet half-side length has been taken as 1.5 × 10¹⁷ cm. We have fixed the jet velocity at v_j = 10¹⁰ cm s${-1}^{}$, the Mach number, M_j ≈ 18, and the jet temperature, T_j = 1.2 × 10⁹   K at injection in all three cases. The jets are composed by protons and electrons and have densities of ρ_j = 1.8 × 10^-29 g cm^-3 for cases 1 and 2, and ρ_j = 8.3 × 10^-29 g cm^-3 for case 3. The adiabatic exponent for the jets, with the given temperature and composition is Γ ≈ 1.58. The ambient medium in the different simulations follows the description given in Sect. 3.1. We note that the shock driven by the jet head in the medium (wind, SNR, and ISM) is called here bow shock (which is more a sound wave rather than a real shock in the hot shocked wind and SNR ejecta), and the jet shocked region is called cocoon, following Bordas et al. (2009).

The simulations were carried out with the finite-volume code Ratpenat, which solves the equations of relativistic hydrodynamics in conservation form using high-resolution-shock-capturing methods. Ratpenat was parallelized with a hybrid scheme with both parallel processes (MPI) and parallel threads (OpenMP) inside each process (see Perucho et al. 2010). The simulations were performed using the computational facility Tirant, at the “Servei d’Informàtica” of the “Universitat de València”, with 12, 8, and 4 processors for cases 1, 2, and 3, respectively. In the axisymmetric simulations, the numerical grid box expands transversally to 100 R_j0, with an extended grid of 100 R_j0 added to the uniform grid. In the axial direction, the grid extends to 150   R_j0 in case 1, and to 100   R_j0 in case 2. The numerical resolution of these simulation is eight cells per initial jet radius. The computational grid, including the extended region is 1200 × 1000 cells. In the case of the simulation of the slab jet, the grid expands 100   R_j0 to each side of the axis, plus 100   R_j0 extra on each side and 100   R_j0 in the axial direction. The resolution used here is four cells per jet initial half-side length. The dimensions of the grid are thus 1200 × 400 cells. The boundary conditions adopted in the Cartesian and cylindrical cases are reflection at the jet bottom side, injection right at the jet base, and outflow in the rest. Reflection in the jet axis has been adopted in the cylindrical simulations.

Fig. 2 Axial cuts of density (top, left), pressure (top, right), temperature (down, left), and velocity (down, right) at the last snapshot of the simulation in case 1 (tsrc = 8900 yr). Units are cgs.

Fig. 3 From top to bottom, maps of density, pressure, temperature, and axial velocity at the last snapshot of the simulation in case 1 (tsrc = 8900 yr). Units are cgs.

Fig. 4 Evolution of pressure (top), density (middle), and temperature (bottom) in the bow shock (left) and cocoon (right) regions for case 1. Units are cgs.

Fig. 5 Position of the jet head (left) and mean radius of the bow-shock region (right) for case 1. The bow shock has already left the shown grid by the end of the simulation. Units are cgs.

3.3. Results

The simulation of case 1 has been run during 8900 yr of evolution. During the first 4500 yr, the jet propagates up to z ≈ 3 × 10¹⁹   cm, with an advance velocity v_h ≈ 2 × 10⁸   cm   s^-1. The jet is initially underpressured with respect to the hot ambient medium (shocked wind). This generates a strong conical recollimation shock in the jet soon after injection, followed by pinching. At z ≈ 3 × 10¹⁹   cm, the jet head reaches the shocked ISM, and its advance is stopped. The speed of the head is reduced to 2 × 10⁷   cm   s^-1, and the jet is only able to propagate up to z ≈ 3.2 × 10¹⁹   cm by the end of the simulation. The expansion of the bow shock during this phase is mainly radial (perpendicular to z), with an expansion to a distance of 4.3 × 10¹⁹   cm from the jet axis (outside the frame of the image). The reverse shock becomes stronger after the collision with the shocked ISM, and pinching is enhanced along the jet. By the end of the simulation, the SNR shock in the ISM has expanded farther, relaxing the density and pressure conditions that the jet finds at arrival, which makes it easier for the jet to carve its way across this region. At the given velocity, the jet would need about 3000 yr more to complete its path through the shocked ISM.

In case 1, the pressure in the cocoon and shocked medium (bow shock) regions is basically constant from the early stages of the simulation: P_s ≈ 5 × 10^-10   erg/cm³. The shocked medium density starts increasing when the jet head reaches the ISM shocked by the SNR from 3 × 10^-26   g/cm³ to 5 × 10^-26   g/cm³ during the period 4500–8900 yr. The temperature in the shocked ambient medium falls slowly during the whole run from 2 × 10⁸ to 1.2 × 10⁸   K. The cocoon density increases from 3 × 10^-28   g/cm³ to 7 × 10^-28   g/cm³ during the first 2000 yr owing to the jet low pressure with respect to the ambient medium, and remains fairly constant afterwards because of the recollimation and pinching shocks. For the same reason, the temperature in the cocoon decreases from 3 × 10¹⁰ to 10¹⁰ K during these first 2000 yr, and becomes constant from that moment on. Turbulence and mixing of jet and external material takes place in the cocoon. Turbulence follows, down to the scales allowed by the simulation resolution, a spectrum similar to Kolmogorov, and the entrained external material actually dominates the mass content of the cocoon. The axial cuts of density, pressure, temperature, and velocity in the last frame, presented in Fig. 2, show the rich structure in the jet after it has been pinched by the ambient pressure, with plenty of shocks, up to the point where the jet material reaches the strong reverse shock in the interaction with the ISM shocked by the SNR. The recollimation and pinching shocks appear early in the simulation given the low initial lateral pressure of the jet with respect to the medium, with the pinching structures getting enhanced when the jet head reaches the shocked ISM. The maps for the density, pressure, temperature, and axial velocity at the last snapshot (8900 yr) of the evolution are shown in Fig. 3. The temporal evolution of the average cocoon and bow-shock density, pressure, temperature, and velocity are shown in Fig. 4. The temporal evolution of the jet head position and the mean radius of the bow shock are shown in Fig. 5.

The simulation of case 2 was run during 9800 yr of evolution. The jet crosses the shocked wind region from injection at z ≈ 3 × 10¹⁸   cm to the contact discontinuity with the shocked ISM at z ≈ 2 × 10¹⁹   cm, at 1500 yr after injection, at a mean velocity of v_h ≈ 3.5 × 10⁸   cm   s^-1. As in case 1, the jet is initially underpressured and collimated by the ambient medium, which produces pinching oscillations along the jet body. After reaching the shocked ISM, the jet head gets flattened by the high-density wall. At that point, the pressure in the bow-shock and cocoon regions starts to decrease more slowly with time. The same happens with the density in the cocoon, whereas the density of the shocked ambient starts to increase, since the head of the jet interacts with a denser medium. The backflow expands sideways, filling the simulated region, whereas the head of the jet acts as a drill on the shocked ISM region owing to its much higher pressure and crosses it at t ≈ 4300 yr after injection. Thus, it needs 3000 yr to cross this small shell. As in case 1, the internal structure of the jet becomes very rich thanks to the interaction between the previous pinching patterns and the backward moving waves produced at the shock of the jet flow with the stagnation point, at z ≈ 2.4 × 10¹⁹   cm. The axial cuts of the jet, presented in Fig. 6, show this effect and the increase in the temperature and decrease in the velocity after each collimation shock, in which strong pressure jumps are produced. After the jet head crosses the shocked ISM layer, it propagates in the ISM with a supersonic speed of  ≈ 2.4 × 10⁷   cm   s^-1, an order of magnitude smaller than the one in the shocked wind region. The maps for the density, pressure, temperature, and axial velocity at the last snapshot (9800 yr) of the evolution are shown in Fig. 7. The temporal evolution of the average cocoon and bow-shock density, pressure, temperature, and velocity are shown in Fig. 8. The temporal evolution of the jet head position and the mean radius of the bow shock are shown in Fig. 9. Overall, the jet evolution of case 2 is similar to case 1, with the cocoon also affected by turbulence and strong mixing.

Fig. 6 Same as in Fig. 2 but for case 2 (tsrc = 9800 yr).

Fig. 7 Same as in Fig. 3 but for case 2 (tsrc = 9800 yr).

Fig. 8 Same as in Fig. 4 but for case 2.

Fig. 9 Same as in Fig. 5 but for case 2.

Fig. 10 Maps of density (top) and pressure (bottom) at the last snapshot of the simulation for case 3 (tsrc ≈ 2000 yr). The shocked wind moves from top to bottom.

Fig. 11 Maps of temperature (top) and axial velocity (bottom) at the last snapshot of the simulation for case 3 (tsrc ≈ 2000 yr). The shocked wind moves from top to bottom.

Fig. 12 Same as in Fig. 10 but for a source age of 300 yr.

Fig. 13 Same as in Fig. 11 but for a source age of 300 yr.

Simulation 3 was run during  ≈ 2000 yr. The jet was injected in the grid in the wind region at z = 1.25 × 10¹⁸   cm. When the jet reaches the shocked wind region, at z = 8 × 10¹⁸   cm, it gets the impact from the shocked wind that has been deflected by the propagation of the binary system in the ISM (see Fig. 1, bottom). This wind is enough in the slab simulation to distort the bow shock of the jet and disrupt the latter as follows. After crossing the wind/shocked wind discontinuity, the bow shock is deviated to the right from injection in the images by the shocked wind coming from top to bottom. The backflow is also deviated in this direction dy the generated pressure gradient. This generates an excess of backflow material in the unshocked wind region to the right of the jet, resulting in a deviation of the jet to the left. The combination of this process with the shocks generated at the different discontinuities destroys the jet even inside the wind region. When the head of the jet reaches the shocked ISM region after 300 yr, at z = 10¹⁹   cm, the jet does not have enough thrust to penetrate this layer, and all the injected flow is thus converted into a backflow that fills both the unshocked and the shocked wind regions. The combination of disruption and advection does not allow the jet to reach the normal ISM regardless of the time elapsed, and thus jets with the given power will be trapped inside the wind/shocked wind region and advected backwards with respect to the system motion. The maps for the density, pressure, temperature, and axial velocity, in the last snapshot (2000 yr) and after 300 yr of evolution, are shown in Figs. 10–13.

4. Nonthermal activity

To estimate what may be the potential nonthermal luminosity of the radiation produced when the jet head propagates within the SNR or the shocked stellar wind, we used a model for the cocoon emission similar to the one presented in Bordas et al. (2009). In short, the model assumes the cocoon as an homogeneous emitter, which is true to first order provided that the flow is roughly sonic. Our simulations show turbulent motion and jet/medium mixing in the cocoon, and the flow becomes trans-sonic far from the jet tip, but we neglect these effects at this stage. In the radiation calculations, following Bordas et al. (2009), the magnetic field (B) energy density (u_B) was assumed to be a 10% of the matter energy density, and the nonthermal injection luminosity, 1% of the jet power. For an optically thin emitter, the radiation luminosity will just be proportional to this fraction. Postshock magnetic field values close to, or not far below, equipartition are expected if the jet keeps a non-negligible fraction of its total energy as magnetic energy, which is a reasonable assumption provided that jets are magnetically launched (e.g. Komissarov et al. 2007; Spruit 2010). Unlike in Bordas et al. (2009), here only the cocoon has been considered as a powerful emitter, since the recollimation and reverse shocks suffered by the jet are strong, but the bow shock is too weak to accelerate particles. The results for the cocoon radiation are presented in Sect. 4.1, in which we focus only on electrons, given the low cocoon densities that prevent efficient gamma-ray production from proton-proton interactions. These low densities also imply that relativistic Bremsstrahlung will be inefficient there.

When the jet head impacts the shocked ISM thin shell, the highest energy particles accelerated in the cocoon region may reach this shell, which in case 2 reaches densities n_ISMs ~ 100 cm^-3. Given the relativistic Bremsstrahlung (electrons) and proton-proton interaction (if protons are also accelerated) timescales and the source age, about 1–10% of the luminosity in relativistic particles reaching that shell could be radiated at very high energies. We note that denser media would increase the efficiency, which is simply  ∝ n_ISMs (if not in the saturation regime).

For the case when the jet has reached the undisturbed ISM, if it is not disrupted, the nonthermal behavior can be described as in Bordas et al. (2009). This would apply for instance to the emission produced during the second half of case 2 simulation. We have not consider particle acceleration in the turbulent and sheared flow in the recollimated jet and cocoon regions, but stochastic and shear acceleration may be relevant there (see Rieger et al. 2007, and references therein). Some thermal radiation may be produced in the shocked ISM. This thermal component would tend to be quite extended and, given the shock velocities, peak in the UV.

We do not consider case 3 here since jet full disruption makes conditions seem less suitable for particle acceleration. It might be that stochastic and shear acceleration may still take place in the jet-disrupted medium, although estimating the nonthermal emission for this case requires an in-depth study that is beyond the scope of this work.

4.1. Cocoon emission: cases 1 and 2

The simulation results are used to characterize the physical properties of the shocked regions. The position and velocity of the bow shock, the internal energy density (u_th; and hence u_B = 0.1   u_th), the temperature, and the rest mass density are used at any source age to account for the evolution of the particle energy distribution, which is affected by both adiabatic and radiative losses: relativistic Bremsstrahlung (negligible), inverse Compton (minor) and synchrotron (dominant). For the IC emission, the radiation field provided by the companion star dominates the CMB. The star is taken to be a bright O, with luminosity 10³⁹ erg s^-1 and temperature 4 × 10⁴ K. At a distance of  ~10 pc, the energy density ratio is u_∗/u_CMB ~ 7, increasing towards the central engine like  ∝ z^-2.

We assumed that the recollimation and reverse shocks accelerate particles, simplifying them as just one accelerator with average properties. The shock velocity was taken to be v_s ~ v_j, which yields an acceleration efficiency Ė ~ (1/2π)   (v_s/c)²   q   B   c ≈ 0.02 × q   B   c (e.g. Drury 1983). As mentioned before, the assumed nonthermal injection luminosity is 1% that of the jet, i.e. 3 × 10³⁴ erg s^-1, and the magnetic-to-internal energy density ratio is 0.1. This gives a B-strength of about  ~4–6 × 10^-5 G depending on t_src. These B-values are one order of magnitude higher than those inferred by Tudose et al. (2006) and Safi-Harb & Petre (1999) for the jet termination regions in Circinus X-1 (a low-mass microquasar) and SS 433, respectively. On the other hand, the ages of these sources are also about one order of magnitude older than those considered here, and thus a lower magnetic field is naturally expected in the former cases. With the adopted magnetic field, synchrotron losses are dominant in the cases studied here. For the conditions assumed, maximum energies are typically  ~10³ TeV. We notice that the adopted nonthermal fraction is quite conservative, and there is room to increase it. As a reference, the values of the radiation model parameters are listed in Table 2 for cases 1 and 2 when the source age is  ~10⁴ yr.

Fig. 14 Synchrotron and IC SEDs for cases 1 (black/thick lines) and 2 (blue/thin lines) computed for tsrc = 103 (left), tsrc = 3 × 103 (middle), and tsrc = 9 × 103 yr (right).

The spectral energy distributions (SED) of the nonthermal synchrotron and IC emission for cases 1 and 2, for evolution times t_src = 10³,3 × 10³, and 9 × 10³ yr, are shown in Fig. 14. Depending on the age of the source and the case, and given that u_B ≫ u _∗ , the synchrotron luminosity can be between  ~10–1000 times higher than the IC one. This difference is more pronounced in case 1. This is a result of the smaller amount of energy pumped by the stellar wind compared to that of an SNR, which leads to lower pressure and therefore lower B in case 2. The IC dominance may be only possible if the B-strength were well below equipartition. At the end of the simulations (t_src = 9 × 10³ yr), the synchrotron fluxes are, for a 3 kpc distance source, of  ~ 20 mJy (case 2) and 100 mJy (case 1), and  ~4 × 10^-12 (case 2) and 7 × 10^-12 erg cm^-2 s^-1 (case 1), in radio and X-rays, respectively. The IC emission yields fluxes of  ~6 × 10^-14 (case 1) and 10^-13 erg cm^-2 s^-1 (case 2), and  ~5 × 10^-13 (case 1) and 3 × 10^-13 erg cm^-2 s^-1 (case 2), in the ranges 0.1–100 GeV and  >100 GeV, respectively. These fluxes might be in principle significantly higher, since the nonthermal fraction is not much constrained.

The cocoon is in pressure equilibrium with the medium. Since the dynamical timescales of the SNR/ISM and wind/ISM interaction structures are longer than the structures generated by the jet, the magnetic field is likely to be approximately constant with time as long as the jet remains inside the SNR/ or the wind/ISM structure. In this context, the radio emission grows with the source age due to particle accumulation, but the X-ray fluxes do not change much, because X-ray synchrotron emitting electrons have time to cool down from quite early times. Concerning IC, the increase in the available cooling time owing to longer t_src compensates for the dilution of the stellar photon field (and subsequent reduction of the emission efficiency) produced by the increase in cocoon size. This renders fairly constant total IC fluxes, although the spectrum below  ~1 TeV becomes softer with time. In Fig. 15, we present the evolution of the specific flux at 5 GHz, and the bolometric fluxes between 1–10 keV, 0.1–100 GeV, and  >100 GeV up to 10⁴ yr. In general, all the lightcurves show a smooth increase for the considered times.

From the morphological point of view, the sources discussed here would appear extended on scales of  ~10′   (d/3   kpc)^-1 in radio, soft X-rays, and high-energy gamma-rays. At hard X-rays and  ~1 TeV, the emission would be rather concentrated close to the jet reverse shock and possibly the jet recollimation shocks, since the emitting electrons, with lifetimes  ≈ 1.3 × 10²   (B/10^-4   G)^-2   (E/10   TeV)^-1 yr, may not have time to fill the cocoon.

Fig. 15 Evolution of the computed fluxes at 5 GHz (top, left), 1–10 keV (top, right), 0.1–100 GeV (bottom, left), and  > 100 GeV (bottom right) up to 104 yr, for cases 1 (black/thick lines) and 2 (blue/thin lines).

5. Discussion

In the light of our results, we briefly discuss here the two HMMQs with very clear evidence of jet interaction with the ISM: SS 433 (case 1) and Cygnus X-1 (case 2). In the case of SS 433, the jet is strongly interacting with its supernova remnant (e.g. Zealey et al. 1980). This source is therefore a clear example of a jet/SNR interaction, but the luminosity of the jet,  ~10⁴⁰ erg s^-1, is much higher than the value adopted here. Nevertheless, there is evidence of nonthermal radio and X-ray emission in the SNR/jet interaction regions of SS 433, which implies that particle acceleration is taking place there (e.g. Safi-Harb & Petre 1999), and possibly also gamma-rays, although at present only upper-limits are available (see Bordas et al. 2010, and references therein). It is worth noting that the interaction of an SS 433-like jet with its external medium has been detected even in nearby galaxies (e.g. Pakull et al. 2010).

Concerning Cygnus X-1, there is a thin emission region that seems to trace the bow shock formed by the jet in an inhomogeneous ISM (the counter-jet is not detected; see, e.g., Gallo et al. 2005). An SNR/jet interaction can indeed be discarded, because the SNR apparently was very weak, as reported by Mirabel & Rodrigues (2003). The structure is elongated, so a shocked wind/jet origin seems unlikely from observations, although the wind energetics (~10³⁶ − 10³⁷ erg ^-1) plus medium inhomogeneity might be a possibility, albeit remote. On the other hand, the deprojected distance of the bow shock to the system is about 10 pc (Gallo et al. 2005), far enough for the jet head to have already left the shocked wind/ISM region. As noted above, the non-detection of emission from the cocoon region may be related to jet disruption. In such a case, radiation may still be generated from particles accelerated via stochastic or shear acceleration or the jet recollimation shock(s), which is expected even if the jet disrupts at farther distances (see Bordas et al. 2011, for a similar situation in extragalactic sources). The bow shock in Cygnus X-1 does not show clear signatures of the system proper motion, which can be explained by the low peculiar velocity of the source with respect to its birth association, Cygnus OB3 (Mirabel & Rodrigues 2003), plus a relatively young jet age (Gallo et al. 2005).

Other high-mass systems may be good candidates for studying the interaction of the jet with its surroundings. Cygnus X-3, for instance, has a powerful jet, and radiation has been detected in the surroundings of this source (Heindl et al. 2003; Sánchez-Sutil et al. 2008), although the nonthermal nature of this emission and its relation to Cygnus X-3 still need to be confirmed. Hints of extended X-ray emission have been also found in LS I +61 303 (Paredes et al. 2007; see however Rea et al. 2010), but the microquasar nature of this source is disputed (see, e.g., Bosch-Ramon & Khangulyan 2009, and references therein).

The number of powerful HMMQs interacting with the environment may be on the order of 10 in our Galaxy¹ (see, e.g., Perucho et al. 2010). The small number of sources and their preferred location in the inner regions of the Galaxy (Bosch-Ramon et al. 2005) make their detection difficult because of strong local/foreground absorption of radio and soft X-ray emission. Diffuse hard X-rays, clustered in the acceleration regions (recollimation and reverse shocks, and possibly the bow shock generated in the jet-ISM interaction), may be better probes for unveiling young HMMQs. Very high-energy radiation could be detected from jet/medium interactions by forthcoming Cherenkov instrumentation (e.g. CTA), which does not suffer from strong gamma-ray background contamination. It could be, however, that most HMMQs would be old and fast enough to have a well-developed bow shock driven by the proper motion of the system. In such a case, the jet would be destroyed and mixed with material from the system motion-driven bow shock, and its large-scale radiation, if any, could be hard to detect.

6. Summary

In the present work, we have explored different possibilities for propagating of a HMMQ jet, assuming different initial times for the jet activity to start with respect to the formation of the HMXB: (i) shortly (~40 000 yr) after the SNR explosion of the progenitor of the compact object; (ii) at epochs (~10⁵ yr) at which only a spherical wind/ISM shock is present; and (iii) at the stages when the system peculiar motion bow shapes the wind/ISM shock to create strong shocked wind convection backwards (>10⁵ yr). For the two first cases, the jet needs about 10⁴ yr to reach the unperturbed ISM, a period much shorter than the potential duration of an HMMQ, which may be up to the lifetime of the remaining massive star, of  ~10⁶ yr. For the third case, the jet gets disrupted and mixed with the shocked stellar wind, and advected away from the stellar bow shock, in just  ~10³ yr. For the cases in which the jet reaches the normal ISM, its propagation can be described as in Bordas et al. (2009), although the impact of crossing the shocked wind/SNR ejecta/shocked ISM can lead to jet disruption, in a similar way to FRI jets

(e.g. Perucho & Martí 2007). When the jet is embedded in the hot shocked wind/SNR ejecta, significant nonthermal emission can be generated in the cocoon region, mainly in radio and X-rays, but also gamma rays. Hard X-rays and very high-energy photons may be the best probes for unveiling the presence of microquasar jets interacting with their environment.

Acknowledgments

The authors thank the referee, Philip Hughes, for his constructive and very encouraging comments. The research leading to these results has received funding from the European Union Seventh Framework Program (FP7/2007-2013) under grant agreement PIEF-GA-2009-252463 V.B.-R. and P.B. acknowledge support by the Ministerio de Educación y Ciencia (Spain) under grant AYA 2007-68034-C03-01, FEDER funds. M.P. acknowledges support from a “Juan de la Cierva” contract of the Spanish “Ministerio de Ciencia y Tecnología” and from the Spanish “Ministerio de Educación y Ciencia” through grants AYA2007-67627-C03-01, CSD2007-00050, and AYA2007-67752-C03-02. P.B. acknowledges support from the German Federal Ministry of Economics and Technology through DLR grant 50 OG 1001. P.B. also acknowledges the excellent working conditions at the INTEGRAL Science Data Center. The authors acknowledge the “Servei d’Informàtica” of the “Universitat de València” for the computing time allocated for this project in “Tirant”.

References

All Tables

All Figures

	Fig. 1 Sketches of a jet propagating through different media depending on the age of the system: the shocked stellar wind, the shocked SNR ejecta, and the shocked and unshocked ISM (top); the shocked stellar wind, and the shocked and unshocked ISM (middle); through the unshocked and shocked stellar wind, and interacting with the shocked ISM under a significant peculiar velocity (bottom).
In the text

	Fig. 2 Axial cuts of density (top, left), pressure (top, right), temperature (down, left), and velocity (down, right) at the last snapshot of the simulation in case 1 (tsrc = 8900 yr). Units are cgs.
In the text

	Fig. 3 From top to bottom, maps of density, pressure, temperature, and axial velocity at the last snapshot of the simulation in case 1 (tsrc = 8900 yr). Units are cgs.
In the text

	Fig. 4 Evolution of pressure (top), density (middle), and temperature (bottom) in the bow shock (left) and cocoon (right) regions for case 1. Units are cgs.
In the text

	Fig. 5 Position of the jet head (left) and mean radius of the bow-shock region (right) for case 1. The bow shock has already left the shown grid by the end of the simulation. Units are cgs.
In the text

	Fig. 6 Same as in Fig. 2 but for case 2 (tsrc = 9800 yr).
In the text

	Fig. 7 Same as in Fig. 3 but for case 2 (tsrc = 9800 yr).
In the text

	Fig. 8 Same as in Fig. 4 but for case 2.
In the text

	Fig. 9 Same as in Fig. 5 but for case 2.
In the text

	Fig. 10 Maps of density (top) and pressure (bottom) at the last snapshot of the simulation for case 3 (tsrc ≈ 2000 yr). The shocked wind moves from top to bottom.
In the text

	Fig. 11 Maps of temperature (top) and axial velocity (bottom) at the last snapshot of the simulation for case 3 (tsrc ≈ 2000 yr). The shocked wind moves from top to bottom.
In the text

	Fig. 12 Same as in Fig. 10 but for a source age of 300 yr.
In the text

	Fig. 13 Same as in Fig. 11 but for a source age of 300 yr.
In the text

	Fig. 14 Synchrotron and IC SEDs for cases 1 (black/thick lines) and 2 (blue/thin lines) computed for tsrc = 103 (left), tsrc = 3 × 103 (middle), and tsrc = 9 × 103 yr (right).
In the text

	Fig. 15 Evolution of the computed fluxes at 5 GHz (top, left), 1–10 keV (top, right), 0.1–100 GeV (bottom, left), and  > 100 GeV (bottom right) up to 104 yr, for cases 1 (black/thick lines) and 2 (blue/thin lines).
In the text