Context. Dense stellar winds may mass-load the jets of active galactic nuclei, although it is unclear on what time and spatial scales the mixing takes place.

Aims. Our aim is to study the first steps of the interaction between jets and stellar winds, and also the scales on which the stellar wind mixes with the jet and mass-loads it.

Methods. We present a detailed 2D simulation – including thermal cooling – of a bubble formed by the wind of a star designed to study the initial stages of jet-star interaction. We also study the first interaction of the wind bubble with the jet using a 3D simulation in which the star enters the jet. Stability analysis is carried out for the shocked wind structure to evaluate the distances over which the jet-dragged wind, which forms a tail, can propagate without mixing with the jet flow.

Results.The 2D simulations point to quick wind bubble expansion and fragmentation after about one bubble shock crossing time. Three-dimensional simulations and stability analysis point to local mixing in the case of strong perturbations and relatively low density ratios between the jet and the jet dragged-wind, and to a possibly more stable shocked wind structure at the phase of maximum tail mass flux. Analytical estimates also indicate that very early stages of the star jet-penetration time may be also relevant for mass-loading. The combination of these and previous results from the literature suggests highly unstable interaction structures and efficient wind-jet flow mixing on the scale of the jet interaction height.

Conclusions. The winds of stars with strong mass loss can efficiently mix with jets from active galactic nuclei. In addition, the initial wind bubble shocked by the jet leads to a transient, large interaction surface. The interaction between jets and stars can produce strong inhomogeneities within the jet. As mixing is expected to be effective on large scales, even individual asymptotic giant branch stars can significantly contribute to the mass-load of the jet and thus affect its dynamics. Shear layer mass-entrainment could be important. The interaction structure can be a source of significant non-thermal emission.

1. Introduction

Jets of active galactic nuclei (AGN) are collimated outflows that originate at the core of galaxies, very likely in the vicinity of a supermassive black hole (SMBH) that accretes matter from its environment (e.g. Begelman et al. 1984). Among the AGN population, the well-known morphological dichotomy between FRI and FRII radio galaxies (Fanaroff & Riley 1974) is thought to be related to jet power (e.g. Rawlings & Saunders 1991). On the one hand, the extragalactic FRII jets, with powers typically over 1044 erg s-1, show an edge-brightened structure in radio. In this case, the jets show a remarkable stability and reach the interaction site with the ambient still collimated, generating a strong forward bow shock, and a reverse shock, known as the hot spot (e.g. Cyg A, Carilli & Barthel 1996; McKean et al. 2016). On the other hand, the edge-darkened jets in radio are related to lower powers and, although they may show a collimated morphology on parsec scales, there is a transition to large opening angles or decollimation (the flaring zone) along the first kiloparsecs that leads to the absence of a hot spot at the ambient interaction regions (e.g. 3C 31, Laing & Bridle 2002a). The suggested reason for the observed jet deceleration has been typically thought to be the entrainment process of cold, slow gas that is incorporated and mixed with the jet flow (Bicknell 1984; Laing 1996).

Mass entrainment seems a natural outcome of jet propagation. Within their first kiloparsecs, jets evolve inside their host galaxies, which contain large amounts of gas, dust, and stars, mainly close to the nucleus (e.g. Burbidge 1970), so there is plenty of material the jet could interact with. It seems unavoidable therefore that stars and medium inhomogeneities will frequently interact with, and in some cases penetrate into, the AGN innermost jet regions. In the case of FRI jets, two main processes have been invoked to explain entrainment: mixing in a turbulent shear layer between the jet and the ambient (e.g. De Young 1986, 1993; Bicknell 1994; Wang et al. 2009), and injection from stellar mass loss (e.g. Phinney 1983; Komissarov 1994; Bowman et al. 1996; Hubbard & Blackman 2006).

Laing & Bridle (2002b) constructed a model of the jet in 3C 31 using the basic conservation laws and the velocity field inferred by Laing & Bridle (2002a), and concluded that the continuous deceleration in the jet requires a monotonic increase in the entrainment rate at large distances, which is incompatible with mass-load from stars, whose density falls rapidly with distance, so entrainment from the galactic atmosphere across the boundary layer of the jet was favoured far from the nucleus. Closer to the active nucleus, mass-load by stars could still be relevant. Concerning the process of entrainment through a turbulent shear layer, Perucho & Martí (2007) showed via a 2D axisymmetric simulation that a recollimation shock in a light jet formed in reaction to steep interstellar density and pressure gradients may trigger large-scale non-linear perturbations that lead to jet disruption and mixing with the external medium.

Recent work by Laing & Bridle (2014) has shown that continuous deceleration of FRI jets by small-scale instabilities at the shear layer could better explain the jet properties inferred from observations and modelling downstream of the flaring region, where these jets increase their brightness and also start to decollimate. This is in agreement with Bicknell (1984) and Wang et al. (2009), who developed models based on this idea. Numerical simulations of the development of Kelvin-Helmholtz instabilities in 2D and 3D also support the slow and continuous entrainment and deceleration due to the growth of short-wavelength Kelvin-Helmholtz instability modes (Perucho et al. 2005, 2010). The relative irrelevance of mass-load by stellar winds from an old population of stars (as expected in giant elliptical galaxies typically hosting jets) on large scales has been confirmed by numerical simulations (Perucho et al. 2014).

On the one hand, mass-load could still represent an efficient deceleration mechanism in the case of low-power jets (Lj ≤ 1042 erg s-1, Komissarov 1994; Perucho et al. 2014). Following the conclusions of Laing & Bridle (2014), it remains possible that the weaker sources in the sample of ten FRI jets in Laing & Bridle (2014), namely M 84 and NGC 193, are decelerated primarily by stellar mass-loading. On the other hand, even a single star with large mass-loss rates may significantly affect the dynamics of relatively weak jets. Hubbard & Blackman (2006) concluded that a Wolf-Rayet star can temporally quench the whole jet flow for jet powers Lj ≤ 1042 erg s-1. Huarte-Espinosa et al. (2013) have shown that a star with a powerful wind (Ṁ = 10-4M⊙ yr-1) can produce observable structures in jets with relatively low powers (Lj = 2.4 × 1043 erg s-1).

Regarding the powerful jets of FRII AGN, observations reveal a decrease in the jet flow velocity from parsec scales, where Lorentz factors are of the order of tens (see e.g. Lister et al. 2013), to kiloparsec scales, where Lorentz factors are closer to 1 (Mullin & Hardcastle 2009). This deceleration could also be partly caused by stellar wind mass-load, although it is more probable that it is due to mass-load at the shear layer and to the conversion of kinetic energy into internal energy at conical shocks along the jet.

As noted above, however, jet interaction with clouds and stars is still expected given the richness of the AGN environment. Even if these objects play a minor role in the jet evolution on kpc scales, this does not prevent a significant dynamical and radiative impact of their interactions with the jet in other contexts. For instance, stellar wind and cloud mass-load in the first parsec may explain the transition from pair- to proton-dominated jets (see Bosch-Ramon et al. 2012, Khangulyan et al. 2013, and references therein). In addition, it has been claimed that interaction with stars and clouds can explain the presence of knots in M 87 (Blandford & Koenigl 1979; Coleman & Bicknell 1985) and in the jet of Centaurus A; knots detected at tens to hundreds of parsecs from the source could correspond to this kind of interaction (Worrall et al. 2008; Goodger et al. 2010). More recently, VLBI observations of the jet in Centaurus A have also revealed a jet structure interpreted as a possible jet- star interaction at sub-parsec scales (Müller et al. 2014). Also, mass-loading by stellar winds has recently been claimed to explain the pressure imbalance between the lobes and the ambient medium in Centaurus A (Wykes et al. 2013; see also Wykes et al. 2015 and references therein). Finally, the radiative counterpart of the interaction of relativistic flows with stellar winds or clouds may be important as well, due to the strong shock produced and the subsequent particle acceleration, which could be – at least in some cases – responsible for the production of steady and variable gamma-ray emission in AGN jets (Bednarek & Protheroe 1997; Barkov et al. 2010; Barkov et al. 2012b; Bosch-Ramon et al. 2012; Khangulyan et al. 2013; Araudo et al. 2013; Bosch-Ramon 2015; Bednarek & Banasinski 2015; de la Cita et al. 2016; Banasinski et al. 2016; Aharonian et al. 2017). Therefore, in short, entrainment by stellar winds cannot be neglected for the study of jet content, long-term jet and lobe evolution, or even to explain emission features at high energies.

The global impact of stars or clouds on the jet propagation and content has been studied analytically or numerically in the past by different authors (e.g. Komissarov 1994; Bowman et al. 1996; Steffen et al. 1997; Choi et al. 2005; Hubbard & Blackman 2006; Sutherland et al. 2007; Jeyakumar 2009). More specific analytical calculations of the dynamical interaction between one cloud or red giant (RG) with an AGN jet have been carried out in the context of radiation studies (e.g. Araudo et al. 2010; Barkov et al. 2010, 2012a; Khangulyan et al. 2013), whereas Bosch-Ramon et al. (2012) performed relativistic hydrodynamic simulations of the interaction of a RG in the innermost region of the jet. In the latter, the star was simulated as a high-density core, surrounded by an atmosphere with a radial profile in density, and in pressure equilibrium. The main results were the formation of a powerful shock; the formation of a potential site of non-thermal processes; and the expansion, disruption, and advection of the material ablated by the jet. The density profile of the wind was shown to play an important role: homogeneous clouds impacted by the jet quickly expand, while a r-2-density profile leads to a smooth atmosphere ablation (with some dependence of the mass extraction rate on the numerical resolution; Bosch-Ramon et al. 2012). In jet regions located farther downstream, it is the wind of the star and not its atmosphere that interacts with the jet, and the shocked wind density profile will determine the properties of the jet-wind interaction at different stages of the process: (i) the star carries a wind bubble in (ram) pressure equilibrium with its environment, simulated here in 3D and treated analytically in the discussion; (ii) right after the star has entered the jet there is a transitory stage corresponding to the interaction between the jet and the stellar wind bubble that has survived jet penetration (this paper concentrates on phases 1 and 2); and (iii) the steady wind-jet interaction (phase 3, not studied here; see de la Cita et al. 2016). Although they may not be the most relevant regarding mass-load, the first two stages imply a sudden quick release of mass into the jet, and a potential large target for jet interaction, energy dissipation, and non-thermal emission, hence the importance of its characterization.

In this paper, we focus on the first stage of the jet-wind interaction, and present relativistic hydrodynamical simulations of a RG or asymptotic giant branch (AGB) star surrounded by its wind interacting with a relativistic jet with moderate power at 100 pc from the galactic nucleus. We assume that the RG/AGB star has just entered into the jet and is releasing the wind bubble formed prior to its entrance in the jet. The initial stages of the process involve the interaction of the jet with the wind bubble. The high density and compression of the gas in the bubble anticipate that cooling must be included, and it represents an important point of this work. Taking into account that the formation of a cometary tail is a natural outcome of the interaction in parallel with the dispersion of the wind bubble, we also study the typical disruption scales of the tail using linear Kelvin-Helmholtz instability analysis. In other words, we study the wind bubble evolution inside the jet and the formation and properties of the comet-like tail, which is also present during the interaction. This tail is eventually incorporated into the jet flow via the growth of instabilities, giving an approximate distance at which this gas could be completely mixed with the jet flow for the case of the formation of a stable channel of shocked stellar wind. The aim is to estimate the scales at which all these processes occur, extending and complementing previous analytic (e.g. Araudo et al. 2010; Barkov et al. 2010, 2012a; Khangulyan et al. 2013) and numerical (Bosch-Ramon et al. 2012; Perucho et al. 2014; Bosch-Ramon 2015; de la Cita et al. 2016) studies of this scenario performed by the authors (see also Hubbard & Blackman 2006).

The points described in the previous paragraph were investigated using numerical simulations. The high resolution and simulation time required to follow the evolution of this thin layer of dense cold gas makes running this simulation in 3D unfeasible. Therefore, we ran a 2D axisymmetric simulation to study the physical scenario in as much detail as possible. This decision has a clear negative effect in terms of the evolution of the cometary tail and mixing, which is limited in the 2D case by the lack of development of helical instability modes and turbulence. Nevertheless, this caveat does not affect the amount of gas that will eventually be mass-loaded, as the latter is the amount of gas in the stellar wind envelope, or the basic dynamics of the bubble. In addition, we performed a short 3D simulation for the initial stages of the interaction of a star with the jet during its first interaction with the jet boundary. At this stage, we neglected the role of the magnetic field in the jet, accounting only for its ram pressure. Relativistic magnetohydrodynamic simulations of this scenario will be performed in our future work. We emphasize that we do not investigate here the absolute jet mass-load by stellar winds, but rather the evolution of the shocked wind in a single jet-star interaction, and the possible distance scales of tail-jet mixing.

The paper is structured as follows: in Sect. 2 we introduce the physical and the numerical set-up of the simulations; in Sect. 3 we present the results of our simulations; Sect. 4 is devoted to the discussion of our results; and in Sect. 5 we summarize our work.

2. Simulations

For the simulations, we used the finite-volume code Ratpenat, which solves the equations of relativistic hydrodynamics in conservation form using high-resolution shock-capturing methods. Ratpenat is a hybrid parallel code: MPI + OpenMP (e.g. Perucho et al. 2010). For visualization we used IDL software and LLNL VisIt (Childs et al. 2012).

The conservation equations for a relativistic flow in 2D cylindrical coordinates (R,z), assuming axisymmetry and using units in which c = 1, are (1)with the vector of variables (2)the vector of fluxes and the source terms (5)with uμ (μ = 0,1,2,3) being the natural velocities. The five unknowns D,Dl,SR,Sz and τ, stand for the densities of the total and leptonic rest masses, the radial and axial components of the momentum, and the energy (excluding the rest mass energy), respectively. All five variables are defined in the laboratory frame, and are related to the quantities in the local rest frame of the fluid (primitive variables) according to where ρ and ρl are the total and the leptonic rest-mass densities, respectively; vR,z are the components of the velocity of the fluid; W is the Lorentz factor (W = (1 − vivi)− 1 / 2, with summation over repeated indices); and h is the specific enthalpy defined as (10)where ε is the specific internal energy and p is the pressure. The system is closed by the Synge equation of state (Synge 1957), as described in Appendix A of Perucho & Martí (2007). This equation of state accounts for a mixture of relativistic Boltzmann gases (in our case, electrons, positrons, and protons). The code also integrates an equation for the jet mass fraction, f. This quantity, set to 1 for the injected jet material and 0 otherwise, is used as a tracer of the jet material through the grid. Cooling has been introduced as a source term, Λ, in the energy equation. The term is defined (in cgs) units as in the approximation given by Myasnikov et al. (1998): Here, we take ne = ρe/me and nZ = np = ρp/mp (Z = 1, we assume hydrogen for simplicity). The equation of state that we use allows us to account for the number of electrons or pairs in a given cell, with ρe and ρp the electron and proton densities, respectively. The cooling fraction is thus computed for each cell using the described procedure. We note in this respect that the particles in each cell are assumed to be in thermodynamic equilibrium and contribute to the local temperature as described in Synge (1957). We consider that the thermal radiative losses are only relevant for the denser, non-relativistic wind in our simulation, where v ≪ c (the shocked wind velocity is vw = 10-4c, and the shocked wind velocity at the dense regions is even smaller, see below). Therefore, we consider that these losses are isotropically emitted and that u0 = 1, ui = 0 in Eq. (5).

The 3D simulation was performed solving the relativistic hydrodynamics equations in Cartesian coordinates. In this case, no cooling source terms were included, and an ideal gas equation of state was used, i.e. a single population of particles was considered.

2.1. Two-dimensional, radiative, axisymmetric simulation (SW2c)

In our 2D simulation we focus on the initial, transitory phase of interaction between a stellar wind envelope (bubble) and a jet at ~100 pc from the central black hole. In particular, we study the effects of cooling and the presence of a dense shell of shocked wind gas surrounding the stellar wind bubble. The unit distance of the simulation is the size of the injector, R0 = 1.1 × 1016 cm, and the unit density is the density given to the wind at the injection region, ρw,0 = 1.5 × 10-19 g cm-3. The grid is formed by a homogeneous resolution region involving a domain in cylindrical coordinates (r,z) of [ 0, 100 ] R0 × [ 0, 200 ] R0 (1.1 × 1018 cm×2.2 × 1018 cm) with a resolution of 6.4 cells per R0, for a total of 640 × 1280 cells.

Snapshot of rest-mass density (bottom half) and modulus of velocity (top half) of simulation SW2c at t = 7.4 × 106 s, in which we take advantage of the axisymmetric nature of the simulation. The plot only includes the central, homogeneous grid to show the set-up of the simulation.

Figure 1 shows a snapshot of the rest-mass density during the initial stages of simulation SW2c to illustrate the set-up of the simulation. The numerical grid is filled by the jet flow, but for a circular region occupied by the wind injector, the bubble, and a dense shell of shocked wind gas expected for a supersonic wind. The injector at the star location is simulated as the cells within 1 R0, where the properties of the flow are fixed as a boundary condition. It is located at position (0,20 R0), i.e. on axis, at 20 R0 (z = 2.2 × 1017 cm) from the simulation boundary in which the jet flow is injected. The wind is injected with radial velocity vw = 10-4c. The stellar wind parameters result in a mass-loss rate Ṁw = 10-5M⊙ yr-1. The boundary conditions are reflection on the jet axis (r = 0 R0), (jet) injection at z = 0 R0 (left boundary as seen in Fig. 1), and outflow at r = 100 R0 = 1.1 × 1018 cm and z = 200 R0 = 2.2 × 1018 cm (top and right in the figure, respectively).

The set-up of the 2D axisymmetric simulation requires that the wind injector be located at the symmetry axis. This is a simplification of the scenario because we locate the wind bubble inside the jet from the beginning of the simulation, i.e. we do not follow the process of entrance of the star into the jet. This also implies that the dense shell of shocked wind gas that surrounds the jet remains basically untouched until the star is fully inside the jet. This condition is relaxed in our 3D simulation (see next section), which was designed to study the interaction of a large bubble of stellar wind formed outside the jet, with the jet boundary.

The bubble region is defined around the injector and is initially 10 R0 in size. This initial bubble size implies a penetration velocity 10 times higher than vw. The bubble region is filled with a wind with decreasing density ρw(r) = ρw,0(R0/r)2 and an initial (constant) temperature Tw = 103 K, which is dynamically irrelevant. A dense shell of shocked wind material is set surrounding the wind bubble: the density of the bubble at its outermost 10% radial region is taken as 10 times higher than that obtained from a pure r-2 density profile to mimic the presence of cooled shocked wind at its termination from the beginning of the simulation.

The presence of such a shell, i.e. formed by shocked wind gas that surrounds the wind bubble, is probably an important ingredient for the mass-load of the jet during the first interaction phase. This shell can be dense enough to change the mass-load pattern in the long-term along this phase. Furthermore, thermal cooling of the wind should be considered in this case because it can be relevant for the evolution of the shocked shell.

We used a version of the code that includes the Synge equation of state (Synge 1957) with two populations of particles, namely leptons (electrons and positrons) and baryons (protons). This allowed us to set up the jet as composed of e+/− pairs, i.e. thermally hot, and the simulated stellar wind as a proton-electron flow. Furthermore, we used thermal cooling terms, following the approximation used in Myasnikov et al. (1998) to account for the cooling of the dense lower-temperature wind and shocked-wind gas. The jet conditions were chosen to obtain a jet power Lj ~ 1044 erg s-1 and the jet radius Rj = 10 pc, i.e. we consider the interaction to take place at ~100 pc from the central black hole for a jet with opening angle of ξ = 0.1 radians. The rest-mass density of the pairs that form the jet is ρj = 6.3 × 10-29 g cm-3, the jet velocity vj = 0.9798 c (corresponding to a Lorentz factor γj = 5), the temperature Tj = 109 K, and the pressure Pj = 10-8 dyn cm-2.

The simulation was run in Tirant at the University of València, in 128 cores, taking approximately 8 × 105 computing hours, mainly owing to the necessary cooling-time check at the time-step calculation routine. The dynamical timescales of the whole simulated jet-bubble interaction process is given by the scale of the bubble divided by the wind velocity, i.e. ~1010 s. At the resolution used, each simulation time-step represents ≃0.045 R0/c × 3.67 × 105s/ (R0/c) 1.65 × 104 s. Therefore, the simulation requires ≃106 code steps to be completed, which is very demanding in terms of computing time, taking into account that each time-step took ≈25 s using 128 cores (which means that it would take around 1 h for a single core to make one step). Altogether, the complex thermodynamics used and the timescale of the evolution of the system convert this simulation into a challenging approach to the scenario, which is currently impossible to reproduce in 3D with the same resolution and similar simulation runtimes. This would multiply the time needed to run the simulation by a factor equal to the number of cells used in the third dimension, i.e. a factor larger than one hundred. In summary, this simulation offers a good approach to the physics of the system because, even if the set-up and geometry miss 3D effects, the computing effort and its emphasis are focused on the (thermo-)dynamics of the interaction.

It is relevant to note that after the first stage is passed where this bubble is shocked and dragged away by the jet, only a central region of the injector remains. In a real scenario, the size of this region is determined by jet-wind ram pressure equilibrium as long as the star is still within the jet cross-section. This means that once the wind bubble is gone, our simulation does not properly resolve the shocked jet-wind structure in the second steady stage of the interaction because this stage occurs within our injection boundary condition, and is therefore not considered here.

2.2. Three-dimensional simulation (SW3)

As a complement to our 2D simulation, we ran a short 3D simulation of a RG/AGB star wind bubble penetrating a relativistic jet (simulation SW3), focusing on the entrance of the star into the jet. This simulation was performed at Mare Nostrum, at the Barcelona Supercomputing Centre, using 256 cores over 5 × 105 computing hours.

In 3D, the injector is simulated as a group of cells within 1 R0 where the properties of the flow are fixed as a boundary condition at a given instant. In contrast, in the 2D simulation the injector propagates across the grid and enters the jet flow from its surrounding medium.

The grid is formed by the domain, in the x,y,z Cartesian coordinates, [ − 10, 10 ] R0 × [ 0, 240 ] R0 × [ 0, 160 ] R0, involving a total of 64 × 768 × 512 cells. The grid is divided into two regions. In the region where y = [ 0, 2.7 × 1017 ] cm (i.e. [ 0, 25 ] R0), we set up the ambient medium (see Figs. 5 and 6), in which the initial bubble and injector are embedded. The jet occupies the region y> 2.7 × 1017 cm, and is set to flow in the positive z direction (upwards in Figs. 5 and 6). The separation between the jet and the ambient medium is thus located at y = 25 R0, along the x- and z-axis. In order to avoid numerical noise, a shear layer has been introduced between the ambient and the jet (see e.g. Perucho et al. 2004). The shear layer has a width of ≃12 R0 (extending 6 R0 around y = 25 R0). The wind injector is placed at (0, 13, 20) in R0 units at t = 0, i.e. at a distance of ≃12 R0 from the shear layer. The wind region occupies a spherical region in the grid with 64 cells of radius within this region.

The sizes of the bubble and the whole simulated region are small enough with respect to the size of the jet (2 × 1017 cm − 1.5 × 1018 cm vs. 6 × 1019 cm) that we can approximate the simulated jet region as a slab, thus omitting the curvature of the jet surface. The boundary conditions are (jet) injection at z = 0 R0 and outflow at the other five boundaries of the simulation box. The time-step of the simulation is ≃0.075 R0/c = 2.75 × 105 s.

The set-up of this simulation is simplified with respect to the simulation SW2c. We use a single ideal gas with adiabatic exponent of 5 / 3 for the whole grid. This implicitly imposes a jet dominated by a non-relativistic proton gas, but this does not significantly affect the results as we focus on the non-relativistic bubble, and largely simplifies and speeds up our calculations. In addition, we note that the jet is strongly supersonic, and so the nature of its content does not significantly affect the results obtained from this set-up. The jet conditions have also been chosen to obtain a jet power Lj ~ 1044 erg s-1 with a jet radius of Rj = 10 pc, as for simulation SW2c. Therefore, the rest-mass density of the jet gas is ρj = 6.3 × 10-29 g cm-3, the jet velocity vj = 0.9798 c (corresponding to a Lorentz factor γj = 5), the temperature Tj = 109 K, and the pressure Pj = 10-8 dyn cm-2.

The medium pressure is set in equilibrium with the jet thermal pressure in order to avoid transversal motions, resulting in a temperature T ≃ 103 K and a density ρa = 4.53 × 10-23 g cm-3. An initial region with size 10 R0, the bubble, is defined around the injector, with a wind with decreasing density ρw(r) = ρw,0(R0/r)2 and an initial (constant) temperature Tw = 103 K, which is dynamically irrelevant. The initial size of the bubble is chosen to be in equilibrium with the thermal pressure of the ambient medium, which is one one-hundredth of the jet ram pressure. These parameters were chosen to emphasize the bubble being shocked and disrupted by the jet ram pressure. In this case, the injector approaches the jet with a velocity vk = 1.1 × 107 cm s-1, which would be the Keplerian velocity for a black hole mass of MBH ≃ 3 × 108M⊙ at a distance of 100 pc. The initial distance between the injector and the shear layer and this velocity give a time ≃1010 s for the beginning of the interaction.

3. Results

3.1. SW2c

Figure 1 shows that the jet-wind interaction evolves as expected during the initial stages of simulation: a bow-shaped shock forms in the jet flow within a short time. Ideally, the jet and wind properties give a stagnation point at the equilibrium position given by (11)from the star (i.e. Rs ≃ R0). In reality, the shock propagates through the shell and the initial wind bubble1 and stops not far from the injector after a time t ≃ 1.5 × 1010 s.

Figure 2 shows the axial position of the shock produced in the jet with time. It stabilizes at z = 1.8 − 1.9 × 1017 cm; i.e. r = 3 − 4 × 1016 cm from the star. This is about three times Rs, significantly farther than it should be. The reason for this is partially that the steady state jet-wind interaction region is not resolved, and the propagation of the stagnation point towards the star stops when it approaches the injector boundary condition due to limited resolution: at r = 3 − 4 × 1016 cm from the star, the number of cells between this position and the injector is only ≃20. We note that in a jet-wind interaction, Rs is approximately the distance from the star to the wind termination shock, and not to the shock in the jet, but this could hardly explain a jet shock position at ~3 Rs. Therefore, the simulation fails to properly describe the jet-wind interaction close to the star for t ≳ 1.5 × 1010 s. Nevertheless, the evolution of the shocked initial wind bubble is still properly described up to later times far from the star.

Figure 3 shows snapshots of simulation SW2c at different times along its evolution. We observe an initial phase (phase 1) where the shock crosses the shell and the wind region and the mass flux slowly increases (see also Fig. 4), and a second phase (phase 2) after the shock has crossed the shell, the bubble expands and fragments, and a large amount of material is ablated and accelerated downstream. In this series of images, we also observe during phase 1 that a cometary tail with fairly homogeneous properties is generated after the bow shock forms. This tail is close to pressure equilibrium with the shocked jet gas surrounding it and expands slowly within the shocked jet gas.

Figure 4 shows the evolution of the mean density in the wake of the interaction within the simulated region (z = 1.1 × 1018 cm from the left boundary). The value at each cell is weighted with the cell volume, which increases with radius in our axisymmetric representation. The plot shows that a single star produces a continuous increase of the mean local jet rest-mass density during phase 1 (this is expected qualitatively beyond numerical effects), and a rapid increase up to more than three orders of magnitude with respect to the initial jet density in the region at the end of phase 2. Clearly, the presence of a dense shell and the inclusion of cooling favours the local increase in density in the shocked regions and therefore the increase in the density of the material dragged downstream. During phase 1, the jet is mass-loaded, but this mass is confined to very thin, heavy, and slow jets within the jet itself. In this regard, we note that the large number of stars expected to fill the jet (e.g. Araudo et al. 2013; Wykes et al. 2013, 2015; Müller et al. 2014; Bosch-Ramon 2015) would lead to a highly inhomogeneous structure with a substantial amount of its mass in the form of dense, slow, and narrow jets. Such a configuration of jets does not imply a stable configuration in the long run. Regarding numerical effects, Klein et al. (1994) claim that convergence in non-radiative shock-cloud interaction simulations is reached for ~100 cells per cloud radius. In our simulation, we used 64 cells per cloud radius, i.e. ≃2 / 3 of the suggested resolution, and thus we are not far below the convergence resolution. Therefore, we do not expect numerical diffusion to play a crucial role. In addition, although our simulation is radiative, the cooling time is much longer than the calculation time-step, meaning that cooling should not significantly affect numerical diffusion either.

Mixing is relevant for mass-load, as it determines how the new matter is integrated in the jet flow. The axisymmetric nature of the simulation prevents the growth of the disruptive helical modes of Kelvin-Helmholtz instability (see e.g. Perucho et al. 2005, and compare with the tail structure of simulation SW3 shown in next section). As a result, the degree of mixing of the wind gas and the jet flow is small until the end of phase 2, when the initial wind bubble is mostly carried away by the jet. Although we cannot follow mixing during phase 1 in 2D simulations, we can state that the amount of gas ablated from the bubble is eventually mixed and loaded within the jet flow and use the properties of the cometary tail to study the scales in which this must happen (see Sect. 4). On the contrary, the disruptive nature of phase 2 favours rapid mixing within our grid. It is also relevant to mention that the jet radius of 10 pc and the star orbital velocity of 1.1 × 107 cm s-1 mean that the crossing time of the star would be of ~6 × 1012 s. Our simulation thus shows that the initial wind bubble is completely driven away by the jet well before the star has crossed the jet.

The final stage of the simulation is more like the case of a homogeneous cloud than of the inhomogeneous cloud in Bosch-Ramon et al. (2012). This is a consequence of the generation of an arc-shaped, cold, dense homogeneous region previous to the disruption phase (see right column in Fig. 3). The region is then destroyed like the leftovers of the homogeneous cloud in that previous paper, after the shock had crossed the cloud. Cooling is an important qualitative difference between this stage in the simulation presented here and that of the inhomogeneous cloud in Bosch-Ramon et al. (2012). This process increases the density of the shocked wind region, so it plays an important role in the details of the shocked bubble evolution, also enhancing the fragmentation of the disrupted cloud.

3.2. SW3

Figures 5 and 6 show cuts of the 3D grid through the centre of the injector at different times (t = 1.28 × 1010 s, and t = 1.8 × 1010 s) for density and tracer, respectively (the tracer is defined as zero for the ambient medium and the bubble, and 1 for the jet flow; it takes values between 0 and 1 in cells where mixing takes place). As the star advances towards the jet and the wind bubble touches the shear layer, it triggers a small wavelength (≪Rj) perturbation and a conical shock that are advected downstream along the shear layer and into the jet, respectively (the extension of the shock is not shown beyond y = 1018 cm in the right panel of Fig. 5). The perturbation is shown as a double wave being advected downstream in the left panel of Fig. 5. The double wave is caused by the extension of the shear layer. This short wavelength perturbation at the jet surface is not followed in our simulation, but its development can be interesting as it may trigger the growth of an instability at the jet shear layer. The wind bubble is distorted at the interaction site, showing elongation and erosion of its outer layers. The tracer plots show that strong mixing takes place at the boundary between the bubble and the jet because of the large velocity gradient there. Another interesting feature of this initial interaction is the wave that propagates upstream along the shear layer. The axial velocity within the shear layer is small enough to allow for this upstream wave motion, unlike the case within the jet which advects all the generated waves. The wave increases the internal energy of the region, expands the shear layer, and favours upstream motion (vz< 0) in the turbulent mixing region of shocked bubble and jet gas. This situation generated problems in the simulation, as the wave reached the bottom boundary of the grid and the inflow condition there was altered from that point on. This compelled us to stop the simulation. This is an issue that will be fixed in future planned simulations, and it does not affect the results presented here in the intended range of applicability.

Normalized mean rest-mass density across the whole simulated jet cross-section for simulation SW2c, at z = 1.1 × 1018 cm (coinciding with half the grid). The dashed line indicates the approximate time at which the transition to phase 2 takes place in the simulation.

Regarding the flow downstream of the interaction region, the conical shock enters farther into the jet flow (towards the right-hand side in Figs. 5 and 6), and the shocked wind and jet material form a complex flow structure where turbulent mixing is already observed within the small grid region. The faster and hotter jet material is mixed with slower and colder shear-layer flow. A short cometary tail forms downstream of the star from wind bubble material eroded by the shear layer and jet flow. This tail is rapidly disrupted, thus favouring the mixing and acceleration of the wind gas; the mixed material visible in blue in the right panel of Fig. 6 propagates downstream with velocity ≃0.4 c. The tail is perturbed by the turbulent shocked jet flow and the transversal velocity of the star, which already gives transversal velocity to the tail gas.

Figures 7 and 8 show a projected 3D image of the rest-mass density (normalized to the wind density at the injector) at times t = 1.45 × 1010 s. The colour scales indicate the flow speed (at the top left of the images) and its rest-mass density (top right). In these plots, the star enters the jet from the right of the images and propagates to the left into the jet flow. The jet, shown in blue, propagates from the top to the bottom, on the left-hand side of the images. The current lines of the jet show how the jet flow deviates from its original direction and slightly decelerates at the site of the shock triggered by the entrance of the star. The lines starting at the wind region show the advection and tangling of the stellar wind gas downstream of the interaction region.

The tail can be seen in those plots as a well-defined, light blue region downstream of the position of the star. The images also show that the cometary tail can be easily disrupted, thus triggering efficient mixing with shocked jet material within a distance ≤7 × 1016 cm. Longer 3D simulations up to later stages of the interaction are planned in order to follow the star penetration into the jet and study the beginning of phase 1 and the eventual formation of a steady tail.

In summary, this simulation provides the following results: 1) The entrance of an obstacle within the jet may trigger the development of small-scale instabilities at the jet boundary and favour shear mixing when these instabilities grow to non-linear amplitudes (e.g. Perucho et al. 2010); 2) a conical shock is advected by the jet and propagates through its cross section, which can also trigger transitory instabilities with longer wavelengths (in the sense that the star entrance is a single event, as opposed to periodical perturbations); 3) a shock wave propagates upstream along the shear layer, which can lead to some energy dissipation upstream in the jet; and 4) mixing between stellar-wind tail and jet gas is very fast in 3D owing to strong non-symmetric perturbations. The mass flux of the tail is actually affected by numerical factors such as resolution (affecting the rate at which the gas is extracted from the bubble), but we would also expect such a turbulent tail to form, due to the growth of small-scale instabilities at the contact discontinuity. A numerical experiment of this (simplified) scenario in 3D is feasible, but will require future simulations with an improved set-up to avoid the shear layer upstream wave, and a much longer running time for completion of the penetration stage.

Two-dimensional cuts (at x = 0 R0) of rest-mass density at different times, for simulation SW3, in cgs units. The cuts were made across the centre of the injector. Left panel: t = 1.3 × 1010 s. Right panel: t = 1.8 × 1010 s. We have cut the plotted region at y = 108 cm because of the lack of relevant structure farther into the jet. We note the changes in the colour scales.

Two-dimensional cuts (at x = 0 R0) of jet mass fraction (tracer) at different times, for simulation SW3, in cgs units. The cuts were made across the centre of the injector. The tracer is zero for the ambient medium, 1 for the jet flow, and it takes values between 0 and 1 in cells where mixing takes place. Left panel: t = 1.3 × 1010 s. Right panel: t = 1.8 × 1010 s. We have cut the plotted region at y = 108 cm because of the lack of relevant structure farther into the jet.

Three-dimensional representations of rest-mass density (simulation SW3) at t = 1.4 × 1010 s including one cut at constant x-coordinate (at half grid, left panel), and the same cut plus a second one at constant z-coordinate (at the position of the wind injector, right panel).

Three-dimensional representation of rest-mass density (simulation SW3) at t = 1.4 × 1010 s including three cuts at half the grid size in the x-coordinate and at the position of the wind injector in the z- and y-coordinates.

4. Discussion

4.1. Mass-loading and mixing

The simulations presented here show that mixing and acceleration are efficient and occur within small scales compared to the jet interaction height (~1 pc versus 100 pc in these simulations). Simulation SW2c shows that radiation cooling and the presence of a shocked wind shell surrounding the bubble lead to explosive disruptive bubble-jet mixing in phase 2. Although it is short, simulation SW3 already shows that even within the small numerical box, i.e. at distances ≤1018 cm, mixing is hardly avoidable even in phase 1 (see Figs. 6 and 8). The lateral motion of the star, together with the complex interplay of waves as the star/wind bubble system enters the jet, produce non-linear perturbations of the tail and its rapid destruction. We expect that 3D simulations including radiative losses could generate denser but still rapidly disrupting tails, able to efficiently mass-load the jet locally.

These results, together with previous studies, and in particular the unstable nature of phase 3 suggested in de la Cita et al. (2016), indicate that the spread in the jet of the loaded matter via turbulent mixing, plus collective mass-load by different stars/clouds, should globally increase the jet density and decelerate it (Wykes et al. 2013, 2015; Perucho et al. 2014, see also Sect. 4.3). The actual mass loaded in the jet, i.e. the actual stellar mass-loss rates and the collective effects of many stars, will determine whether this process is important for the jet content (e.g. Perucho et al. 2014). Although efficient mixing is predicted, the mass-loading process may still lead to strong jet inhomogeneities in density and composition if the loaded material is not rapidly spread across the jet cross-section.

We have computed the amount of mass dragged by the jet after ~1010 s ≃ 103 yr, the dynamical time of the bubble, and it is ≈2 × 10-2M⊙, which is compatible with the adopted total mass of the bubble plus the dynamical time and plus mass-loss rate of the star. Making the raw assumption that the loaded mass is homogeneous in time (although it is actually concentrated towards the end of phase 2, see Fig. 4), we obtain a mean input of ~10-5M⊙/ yr for a single AGB star during phases 1 and 2 (and in phase 3 as well, as wind keeps being injected in the jet). Such a mass-load rate is of the order of the mass flux expected in FRI jets (see e.g. Laing & Bridle 2002b; Perucho & Martí 2007). This is therefore a critical stage in the process of mass-loading by stars and should be taken into account in future numerical simulations of jet deceleration (Perucho et al. 2014). Interestingly, ~0.01% of ~1 M⊙-stars are at present AGB stars. Thus, the jet of a radio galaxy such as M87 may contain ~103 AGB stars up to kpc scales (see the discussion in Bosch-Ramon 2015). We note nevertheless that this conclusion applies to mass-load when the bubble is well inside the jet, but the mass-load may also be relevant during the early jet-penetration stage of the bubble, when it is in direct interaction with the jet boundary.

4.1.1. Mass-load in the shear layer

In this work, we have assumed that once the bubble is inside the jet, the bubble radius (Rb) is determined by equating the time needed by the bubble to enter the jet to the time needed for the jet to shock the bubble. However, the radius of the bubble before entering the jet is determined by equating the wind and the medium (ram) pressures. This yields a value for Rb outside the jet that is larger by a factor of ~max( (Pjet/PISM)1 / 2 (vw/vorb), 1) than inside the jet, with (where vorb is the star velocity). The consequence of this is that the mass-load in the jet shear layer will be also larger by the same factor. Given that SW3 shows that the matter loaded in the shear layer efficiently spreads inside the jet, it seems reasonable to assume that all the mass accumulated in the bubble before jet penetration is eventually entrained by the jet but when Rb ≳ Rj(z), i.e. close enough to the black hole. Upstream of this point, the specific bubble-jet contact geometry should be taken into account.

In fact, the bubble radius outside the jet may be so large that the mass injected by the bubble to the jet shear layer surpassed that injected in phase 3. This would happen for Rb ≳ Rj (vw/vorb), which will occur at (12)where n0 is the number density of the ambient medium (assuming hydrogen), ξ is the jet opening angle (taken to be 0.1 radians in the paper, see Sect. 2), and we adopt the quantity format Ax = A/ 10x, with A in the corresponding units.

We thus conclude from these basic estimates that first, shear-layer entrainment that may dominate mass-load occurs during the early jet penetration stage. Then, this work shows that phases 1 and 2 lead to an explosive release of mass, smaller than in the penetration phase and well inside the jet. And finally, phase 3 mass-load occurs more smoothly than phase 2 mass-load, but still in an unstable manner all through the jet (de la Cita et al. 2016).

4.1.2. Scalability of the results

The star mass-loss rate adopted here is high, more typical of an AGB than a red giant. Nevertheless, we note that for the 3D adiabatic simulation SW3, the results can be scaled keeping Ṁ/Lj constant, where Lj is the jet kinetic power. In addition, the general discussion on stability should qualitatively apply as well to lower mass-loss rate stars for the same Lj, although a quantitative assessment requires specific simulations, and the same applies to faster winds, such as those of massive stars. Simulation SW2c is not scalable as it includes more detailed physics that require fixing units. Lighter winds could still cool down effectively, given the dynamical timescales of the simulations, as the shell cooling timescales are short enough; this is safely the case for a range of Ṁ spanning few orders of magnitude below the adopted value. This can be seen by estimating the cooling effect on the region of the stellar wind shock. The post-shock temperature should be ≈4 × 104 K, which implies a very fast cooling regime. The cooling rate is Q = λn2, where λ = 10-22 ergcm3s-1. The wind density, without compression, is 105 cm-3, so the cooling time would already be tcool ≈ 106 s; accounting for shock compression, tcool ≈ 2 × 105 s. These timescales should be compared with the adiabatic evolution time, tadd ~ Rs/vw ≈ 3 × 109 s, and to the simulation time-step, ≃1.6 × 104 s. Thus, for SW2c, one has tadd ≫ tcool. Although wind cooling is thus of particular importance for RG/AGB stars, it must be taken into account even for the case of lighter stellar winds. Regarding massive stars, their winds are less dense and much faster, and therefore the cooling will be significant only for much closer shock radii, i.e. for very powerful jets, or much closer to the jet base.

4.1.3. Tail stability and mixing scales

We focus now on the stability properties of wind-tails, a relevant aspect to mixing. In 3D, the motion of the star across the jet already provides the tail with a transversal velocity, inducing the growth of helical structures. This perturbation is non-linear and favours mixing locally or, at least, close to the interaction region comparing to the jet scale. On the contrary, the tail of shocked wind material in 2D appears to generate a stable structure as the shock crosses the bubble (phase 1). The axisymmetry and the absence of any perturbation in a region that is in pressure equilibrium with its surroundings avoids the trigger of disruptive instability modes. Despite these simplifications, the properties of the tail obtained in the simulation are of the same order as those for the 3D case: the tail shows densities of 103 − 104 times the shocked jet density and typical velocities ≃0.05 − 0.1 c during this transitory phase that lasts for ~1010 s. Assuming that the mass flux is indeed mimicking an actual physical process despite numerical contamination (see Sect. 3.2), we took these typical values to run a stability analysis of the tail/shocked jet system.

We calculated the solutions to the stability problem of a sheared flow (see e.g. Perucho et al. 2007). The formation of a shear layer is expected in situations in which there is some mixing-dissipation at the contact discontinuity. The formation of the tail from the erosion of the external layers of the bubble should lead to such a structure. We adopted the spatial approach in our stability calculation, so we obtained the real and imaginary parts of the wavenumber component along the z-direction, kz,r and kz,i, respectively, in terms of a real value of the frequency ω for solutions of the type δ(r,θ,z,t) = A0e− i(krr + nθ + (kz,r + ikz,i)z + ωt) = A0ekz,ize− i(krr + nθ + kz,rz + ωt), where kr is the radial wavenumber (a function of kz and ω), and n = 0,1 for the pinching and helical modes, respectively (see e.g. Perucho et al. 2005). Neglecting the imaginary part of the radial wavenumber, kr, with respect to the imaginary part of the axial wavenumber, kz, we can write the axial dependence of the amplitude as A(z) = A0ekz,iz. The units are given by the speed of light and the radius of the tail (Rt ≃ 1017 cm).

Solutions of the linear perturbation equation for typical profiles of rest-mass density and axial velocity obtained from the simulations (left panel: pinching, symmetric mode; right panel: helical mode). The tail density is taken as ρt/ρj = 103; the tail and jet velocities are vt = 0.05 c and vj = 0.95 c, respectively; and the sound speeds are cs,t = 0.04 c and cs,j = 0.54 c. The dashed and dotted lines show the line corresponding to phase velocities, w/k, equal to the relativistic composition of the velocity and sound speed for each of the flows, i.e. 0.98 c for the jet and 0.09 c for the tail. These velocities separate the solution plane in different allowed and forbidden regions (see e.g. Payne & Cohn 1985, for the non-relativistic case).

The solutions to the stability problem are given in Fig. 9. The left panel shows the solutions for the pinching mode (n = 0), and the right panel shows the solutions for the helical mode (n = 1), for the stability equation using ρt/ρj = 103 (ρt/ρj = 104 has also been calculated, but is not shown), vt = 0.05 c, and vj = 0.95 c, and specific internal energies of ϵt = 2.5 × 10-3c2 and ϵj = c2, which are typical values given by our simulations. The solutions have been derived for the case of a transition with a thickness of the same order as the tail radius Rt, as also given by the simulations. Each of the curves observed in the panels displayed in Fig. 9 stand for a particular unstable mode that can be triggered with different frequencies and wavelengths: the left panel shows the solutions for the pinching symmetric mode, whereas the right panel shows the solutions for the helical mode. The upper curves within each panel stand for kz,r, whereas the lower ones represent the corresponding kz,i. The inverse of the latter is known as the growth length, which defines the distance at which the amplitude of the perturbation e-folds. The solutions were obtained point to point using the shooting method (a combination of a Runge-Kutta integrator with variable step and a root-finding Müller method, Roy Choudhury & Lovelace 1984; Perucho et al. 2007), which is the reason for the discontinuous appearance of the curves (which is irrelevant for our analysis, though). The denser tail (ρt/ρj = 104) results in slightly longer, but similar, growth lengths (smaller values of kz,i). Increasing the size of the layer would increase the growth distances, i.e. increase the stability, of short wavelength modes (λ ≤ Rsh, with Rsh the size of the transition), whereas the properties of long wavelength modes would not change significantly (e.g. Perucho et al. 2005).

In the case of the pinching mode, the solutions of the fundamental mode (a mode with no nodes between the tail axis and the shear layer of the tail) appear only at small phase speeds (k ≫ w) and show small relative growth rates. The helical surface mode is, on the contrary, the most unstable for low frequencies (long wavelengths). Body modes are characterized by their different radial structure and they appear as different curves in the solution, higher order modes showing up as kr grows. The minimum growth lengths obtained from the solutions of the linear problem correspond to maximum values of . Considering a perturbation of ten per cent of the background value, A0 = 0.1 (the amplitude is normalized to the background value) for variable X, the perturbation will become fully non-linear when A ~ 1. Taking into account that A(z) = A0ekiz, and using the , A becomes ~1 at z ≃ 20 Rt ≃ 1018 − 1019 cm, for Rt = 5 × 1016 − 1017 cm. This means that the tail material is expected to mix with the jet within distances of ~1 pc if an initially stable tail is perturbed. From previous stability studies (e.g. Perucho et al. 2005) we can claim that these distances are reduced when the relative density of the tail with respect to the jet decreases (i.e. closer to the active nucleus), and that they are longer when this relative density increases (e.g. towards the end of phase 2). The increase in the mass-flux at the end of phase 2 could result in an important increase in the growth lengths and result in the formation of longer tails. Nevertheless, if our estimates are correct, one would expect to find jets with an inhomogeneous structure caused by tails ~1 pc long downstream of jet/star interaction sites. Such a scenario was recently suggested for VLBI observations of the jet in Centaurus A by Müller et al. (2014). We can conclude that tails are, in principle, stable transient structures due to their inertia (Perucho et al. 2005), but that the strong perturbations expected – and observed in simulation SW3 – may force efficient mixing relatively close to the interaction site (≃1 pc).

It is also relevant to mention that the explosive nature of the disruption of the wind bubble will contribute to efficient mixing on the scales of the interaction region. Such explosive events may be bright enough to be observed in close sources. If we consider a mean jet radius of Rj ~ 10 pc along its trajectory within the host galaxy and a mean orbital velocity of vo ~ 107 cm s-1, the mean crossing time of the stars is tc ~ 6 × 1012 s. Taking into account that the disruption phase lasts for ~1010 s (see Fig. 4), we expect that about one one-thousandth of the red giant stars in the jet should be at this stage and possibly produce bright spots within the jet, as do those observed in Centaurus A (Worrall et al. 2008; Goodger et al. 2010).

4.2. Non-thermal activity

If particles are accelerated at the jet shock, recent numerical studies of cloud/stellar wind-jet interactions and of their related non-thermal emission show that the interaction region will produce significant amounts of detectable radiation (Bosch-Ramon 2015; de la Cita et al. 2016; Vieyro et al. 2017). The main factors behind the high apparent radiation efficiency are jet energy dissipation, which actually occurs on scales much larger than the jet shock stand-off distance from the obstacle (Bosch-Ramon 2015), and Doppler boosting effects caused by a post-shock flow that is also moving at relativistic speeds (see also the analytical studies of Barkov et al. 2012b; Khangulyan et al. 2013, where the obstacle itself moves at relativistic velocities). Confirming and extending the results of Bosch-Ramon et al. (2012) for a homogeneous obstacle, the present work suggests that the initial wind bubble can also trigger substantial non-thermal activity. Given its size, much larger than the more or less steady interaction structure simulated by e.g. Bosch-Ramon (2015) and de la Cita et al. (2016), the bubble represents a temporal target for the jet that may eventually cover a significant fraction of the jet section (as in Khangulyan et al. 2013). The details of the subsequent radiation will strongly depend on the height of the interaction region, in particular the balance between radiative and non-radiative losses (Bosch-Ramon 2015), and in certain cases the bubble penetration into the jet may lead to detectable transient emission. A quantitative study of this emission is beyond the scope of this work, but it is worth noting that the star penetration stage should be also accounted for when assessing the role of stars in jet high-energy emission. We also note that such a study would also strongly benefit from the inclusion of the magnetic field in the simulations.

4.3. Stellar wind induction

If non-thermal particles are produced, their presence can increase the wind mass-loss rate if the jet shock can reach close to the stellar surface. This will happen if the jet has a very high ram pressure, or when the interaction takes place close to the jet base (Khangulyan et al. 2013). This mass-loss rate amplification is close to the well-known effect of wind induction by X-ray radiation in close X-ray binaries and in AGNs (Basko & Sunyaev 1973; Dorodnitsyn et al. 2008a,b). We estimate in what follows the wind mass-loss rate following Khangulyan et al. (2013).

One can estimate the heating rate due to non-thermal particles interacting with the stellar photosphere as , where χ is the efficiency of particle acceleration. The excited mass flux can be estimated as g/s cm2, where α-12 is the efficiency of wind induction, and R∗ and M∗ the stellar radius and mass, respectively. The parameter α-12 is ≃1 in the case of X-ray heating; we take it here as a fiducial value, but it should be derived for the specific case of non-thermal heating. The total mass-loss rate in the wind can be estimated as (13)To compare the intensity of the induced wind and the normal mass-loss rate, we recall that the mass-loss rate of RG/AGB star is assumed to be Ṁ = 10-5M⊙/yr or 0.7 × 1021 g/s. Combining this value with Eq. (13), we get the condition under which the induced stellar wind starts to dominate over the normal one: (14)where .

A jet can be significantly decelerated by the wind if its power is smaller than the critical value erg s-1. In the induced wind case, the jet deceleration limit does not depend on jet power, and the maximum jet radius for which the jet will be decelerated by the wind can be estimated as (15)where N∗ is number of stars located in the jet at a given distance from the SMBH. In conclusion, the induced stellar wind could be a dynamically important factor for a star with a big radius where the jet is passing close to its base (not simulated here) and independently of its power.

4.4. Potential effects of magnetic field

A further step of this study should be performed including magnetic fields in order to derive more conclusive results, taking into account that the magnetic field can be a relevant factor at the interaction region. Moreover, a gas cloud embedded in a jet, like the scenario studied in this paper, can trigger the growth of Rayleigh-Taylor instability (Chandrasekhar 1961). In this respect, the magnetic field can also play an important role (Imshennik 1972). A brief summary of the role of magnetic field is in order here.

Following seminal papers by Klein et al. (1994) and Poludnenko et al. (2002) on the interaction between clouds and shocks in a hydrodynamical context, several works have studied the influence of magnetic fields on the interaction between gas clouds and shocks in the context of supernova remnants or dense clouds interacting with the interstellar medium (see e.g. Mac Low et al. 1993; Jones et al. 1996; Miniati et al. 1999; Gregori et al. 1999). In Mac Low et al. (1993) and Jones et al. (1996) the authors study the effect of parallel and perpendicular fields on clumps of gas interacting with winds using axisymmetric simulations mainly. These works indicate that the orientation of the field has relevant consequences for the interaction. On the one hand, when the magnetic field is parallel to the shock surface, it is amplified by the interaction and, when stretched around the cloud, it prevents its disruption in that direction, but can amplify cloud disruption in the axial direction, thus leading to asymmetrical cloud disruption. On the other hand, when the field is perpendicular to the shock (axial), the lines are stretched around the cloud and show reversals that may induce reconnection in this region. The magnetic tension of the axial field prevents the formation of vortexes, also reducing the fragmentation of the cloud. For oblique fields, the evolution resembles one case or the other more closely depending on the angle between the lines and the direction of motion (Jones et al. 1996; Miniati et al. 1999). Gregori et al. (1999) presented 3D simulations in which the magnetic field is perpendicular to the propagation, deriving similar results. To our knowledge, there are still no works that study the evolution of the tail, although in Miniati et al. (1999) the asymmetric nature of the simulations allowed oscillations in the tail that could couple to a kink mode of the current-driven instability or a helical KH mode.

5. Summary

We have performed 2D and 3D simulations of the first stages of the interaction between a star surrounded by its wind bubble and an AGN jet. At the beginning of the interaction, the wind bubble is shocked while being slowly eroded by the shocked jet flow. This leads to a transient continuous comet-like tail of wind material, which likely gets disrupted locally, as 3D simulations and analytical stability calculations show. This phase is followed by the eventual expansion, disruption, and spread of the material initially present in the wind bubble. For this to occur, the presence of a previous shocked wind shell and radiation cooling are important, showing that the case of a r-2 density profile wind actually resembles that of a homogeneous cloud. Although steady jet-wind interaction when the two flows balance the other’s pressure has not been simulated here, complementary results from de la Cita et al. (2016) show that the interaction structure is also likely unstable, which would favour jet-wind mixing. Finally, although the jet-wind interaction dynamics tends to distribute the stellar wind mass on large jet regions, some degree of inhomogeneity may be expected in extragalactic jets. The mass-loss rate adopted, 10-5M⊙/ yr, which is also the mean mass-load in the jet, represents a significant amount of the expected mass fluxes in FRI jets. In fact, the mass entrained through the shear layer could be a substantial fraction of, or even dominate, the total mass loaded by a star into the jet. If a moderate number of AGB stars are present in the centre of galaxies, as expected, then AGN stars could easily be important channels for mass-loading in AGN jets.

The transitory wind bubble penetration of the jet should be studied as a potential source of non-thermal emission in addition to other stages of the jet-star interaction, and the associated energetic particles may induce additional jet mass-load for jet-star interactions close to the stellar atmosphere. The magnetic field must be included in order to more accurately characterize jet-wind mixing, and the scenario-related non-thermal processes.

Acknowledgments

The authors acknowledge Sarka Wykes for the comments on the manuscript and the discussion. M.P. acknowledges support by the Spanish “Ministerio de Economía y Competitividad” grants AYA2013-40979-P and AYA2013-48226-C3-2-P. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement PIEF-GA-2009-252463. V.B.-R. acknowledges support by the Spanish Ministerio de Economía y Competitividad (MINECO/FEDER, UE) under grants AYA2013-47447-C3-1-P and AYA2016-76012-C3-1-P with partial support by the European Regional Development Fund (ERDF/FEDER), MDM-2014-0369 of ICCUB (Unidad de Excelencia “María de Maeztu”), and the Catalan DEC grant 2014 SGR 86. V.B.R. also acknowledges financial support from MINECO and European Social Funds through a Ramón y Cajal fellowship. This research has been supported by the Marie Curie Career Integration Grant 321520. B.M.V. acknowledge the partial support by the NSF grant AST-1306672 and DoE grant DE-SC0016369. The authors thankfully acknowledge the computer resources, technical expertise, and assistance provided by the “Centre de Càlcul de la Universitat de València” through the use of Tirant, the local node of the Spanish Supercomputation Network, and that also provided by the Barcelona Supercomputing Center “Centro Nacional de Supercomputación”.

All Figures

Snapshot of rest-mass density (bottom half) and modulus of velocity (top half) of simulation SW2c at t = 7.4 × 106 s, in which we take advantage of the axisymmetric nature of the simulation. The plot only includes the central, homogeneous grid to show the set-up of the simulation.

Normalized mean rest-mass density across the whole simulated jet cross-section for simulation SW2c, at z = 1.1 × 1018 cm (coinciding with half the grid). The dashed line indicates the approximate time at which the transition to phase 2 takes place in the simulation.

Two-dimensional cuts (at x = 0 R0) of rest-mass density at different times, for simulation SW3, in cgs units. The cuts were made across the centre of the injector. Left panel: t = 1.3 × 1010 s. Right panel: t = 1.8 × 1010 s. We have cut the plotted region at y = 108 cm because of the lack of relevant structure farther into the jet. We note the changes in the colour scales.

Two-dimensional cuts (at x = 0 R0) of jet mass fraction (tracer) at different times, for simulation SW3, in cgs units. The cuts were made across the centre of the injector. The tracer is zero for the ambient medium, 1 for the jet flow, and it takes values between 0 and 1 in cells where mixing takes place. Left panel: t = 1.3 × 1010 s. Right panel: t = 1.8 × 1010 s. We have cut the plotted region at y = 108 cm because of the lack of relevant structure farther into the jet.

Three-dimensional representations of rest-mass density (simulation SW3) at t = 1.4 × 1010 s including one cut at constant x-coordinate (at half grid, left panel), and the same cut plus a second one at constant z-coordinate (at the position of the wind injector, right panel).

Three-dimensional representation of rest-mass density (simulation SW3) at t = 1.4 × 1010 s including three cuts at half the grid size in the x-coordinate and at the position of the wind injector in the z- and y-coordinates.

Solutions of the linear perturbation equation for typical profiles of rest-mass density and axial velocity obtained from the simulations (left panel: pinching, symmetric mode; right panel: helical mode). The tail density is taken as ρt/ρj = 103; the tail and jet velocities are vt = 0.05 c and vj = 0.95 c, respectively; and the sound speeds are cs,t = 0.04 c and cs,j = 0.54 c. The dashed and dotted lines show the line corresponding to phase velocities, w/k, equal to the relativistic composition of the velocity and sound speed for each of the flows, i.e. 0.98 c for the jet and 0.09 c for the tail. These velocities separate the solution plane in different allowed and forbidden regions (see e.g. Payne & Cohn 1985, for the non-relativistic case).

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

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