© ESO, 2016

1. Introduction

The winds and atmospheres of stars have been proposed to play an important role in the propagation, matter content, stability, and potential disruption of the jets of active galactic nuclei (AGN; e.g. Komissarov 1994; Bowman et al. 1996; Hubbard & Blackman 2006; Bosch-Ramon et al. 2012b; Perucho et al. 2014). In addition, the interaction of jets with stellar atmospheres or winds have also been suggested to lead to non-thermal emission that may be detectable from Earth, in both blazar and non-blazar AGN, and in the form of both transient and persistent radiation (e.g. Bednarek & Protheroe 1997; Barkov et al. 2010, 2012a,b; Khangulyan et al. 2013; Bosch-Ramon et al. 2012b; Araudo et al. 2013; Bosch-Ramon 2015; Bednarek & Banasiński 2015). It is noteworthy that there might already be direct observational evidence of such an interaction (e.g. Hardcastle et al. 2003; Müller et al. 2014), and that high-energy phenomena observed in some AGN might be interpreted in the context of jet-star interactions (e.g. Barkov et al. 2010, 2012a,b; Khangulyan et al. 2013). The actual extent of the dynamical and radiative impact of these interactions is still unknown, however.

In the jet-star interaction scenario, stars are expected to cross AGN jets at different distances from the jet base, producing shocks that can transfer kinetic energy to non-thermal particles. In the case of highly magnetized jets, these would present different mechanisms to transfer jet energy, now in the magnetic field instead of being in kinetic form, to non-thermal particles (see, e.g. Bosch-Ramon 2012, and references therein). Electric fields produced by jet-stretched magnetic field lines around an obstacle would also accelerate particles (e.g. Jones et al. 1996).

In the inner-most jet regions, stars are expected to move fast and the jet is narrower, which means that only a few stars can interact with the jet at a time. Here the emission might be released during relatively short events that are triggered by high-inertia targets, such as the external weakly bound layers of evolved giants (e.g. Barkov et al. 2010), or stars with very high mass-loss (e.g. Araudo et al. 2013). Persistent emission might also take place far from the jet base, as the jet propagates through the inner-most kpc regions of the galaxy, and even farther out (Araudo et al. 2013; Bednarek & Banasiński 2015). In all these situations, several ingredients are required to accurately estimate the interaction duration, rate, effect on the jet properties, and related radiation: a proper characterization of the stellar populations and their spatial distribution, both galaxy-type dependent, and a detailed description of the physics of the jet-star interaction and the associated non-thermal processes. An approximate study combining both hydrodynamics and radiation estimates was carried out, the result of which was that the emission from individual jet-star interactions may be significantly higher than previously thought (Bosch-Ramon 2015).

A proper understanding of the AGN jet radiation and its underlying physics needs an accurate characterization of the emission produced in stars interacting with AGN jets. To proceed in this direction, following previous works, we here combine relativistic hydrodynamical (RHD) simulations and radiation calculations, assuming no dynamical feedback from non-thermal processes (a rather good assumption in the cases studied here), to characterize the emission produced in the shocked jet region that forms when the jet flow is stopped by a stellar wind. First we simulate the interaction of a stellar wind and a relativistic jet until steady-state is reached. Then, the obtained hydrodynamical information is used to characterize the injection of non-thermal particles, their propagation, and emission. As matter and radiation densities are generally low unless interactions occur close to the jet base (see Barkov et al. 2012a; Khangulyan et al. 2013), hadronic processes are not expected to be relevant. Thus, the emitting particles are assumed to be leptons here, which might be electrons for a proton-dominated jet, or electrons and positrons (e^±) for a e^±-dominated jet. We consider a situation in which a star of relatively high mass-loss and luminosity, such as a red giant or a moderately early star, interacts with a jet of intermediate power. The results can be scaled, which allows deriving broader conclusions. Throughout this paper, primed quantities are in the fluid reference frame (FF).

2. Hydrodynamics

2.1. Simulated cases and streamline preparation

We performed axisymmetric RHD simulations of the interaction between a relativistic jet and a stellar wind within the jet. The simulations were conducted using a finite-difference code based on a high-resolution shock-capturing scheme that solves the equations of RHD in two dimensions (2D) in a conservation form (Martí et al. 1997). The code is parallelized using open message passing (OpenMP; Perucho et al. 2005). The simulations were run in a workstation with two Intel(R) Xeon(R) CPU E5-2643 processors (3.30 GHz, 4 × 2 cores, with two threads for each core) and four modules of 4096 MB of memory (DDR3 at 1600 MHz). The obtention of streamlines from the simulation results allowed characterizing the regions of interest as many 1D hydrodynamical structures, which is suitable for radiation calculations.

2.1.1. Jet-star interaction

We assumed a collisionless adiabatic and relativistic ideal gas with a dynamically negligible magnetic field. For simplicity, the gas has one particle species with a constant relativistic adiabatic index of 4/3 for both the jet and the stellar wind material. It is not required for our purposes at this point to determine whether the jet is made of protons (nuclei) plus electrons, or e^±-pairs. The physical size of the domain is r ∈ [0,l_r] with l_r = 2 × 10¹⁵ cm, and z ∈ [0,l_z] with l_z = 1.5 × 10¹⁵ cm. The total number of cells is 400 and 300 in the radial and axial directions, respectively. This resolution was chosen to have enough numerical dissipation to avoid a growth of instabilities that is fast enough to prevent the formation of a quasi-steady state (see Sect. 2.2). The star was located at (r₀,z₀) = (0,0.3 × 10¹⁵) cm, and its spherical wind was injected through a region with radius r_in = 7 × 10¹³ cm (14 cells), small enough not to be affected by the shock terminating this wind. The jet was injected at the bottom boundary of the grid. The jet streamlines were approximated as parallel instead of radially extending from the jet origin because the scales of the simulation were much smaller than the height of the jet at which the interaction takes place. The upper and right boundaries of the grid were set to outflow, while the left boundary was set to reflection. All these parameters are summarized in Table 1.

The physical parameters of the jet, which is in an inflow condition at the bottom of the computational grid, are the total jet power within the grid, L₀ ≈ 4 × 10³⁷ erg s^-1 (~10⁴⁴ erg s^-1 for a 1 pc jet radius), the Lorentz factor , with v₀ = v_z = 0.995 c and v_r = 0, and the specific internal energy ϵ₀ = 9 × 10¹⁸ erg g^-1. The jet power is computed as (1)where ρ₀ = 1.24 × 10^-27 g cm^-3 is the jet density, its pressure, and its specific enthalpy.

The stellar wind is a spherical inflow condition imposed at 7 × 10¹³ cm from the star centre. The wind physical parameters at injection are the mass-loss rate Ṁ = 10^-9M_⊙ yr^-1, the radial velocity v_sw = 2 × 10⁸ cm s^-1, and the specific internal energy ϵ_sw = 9 × 10¹³ erg g^-1. The derived stellar wind density at injection is ρ_sw ≈ 5.2 × 10^-21 g cm^-3. The stellar wind is taken to be homogeneous, meaning that it is not clumpy, and since it is supersonic, its density profile from the injection radius up to its termination is ∝1 /R², with R being the distance to the star centre. The wind properties correspond to those of a high-mass star with a modest mass-loss rate. The thrust of this wind, Ṁv_sw ≈ 1.3 × 10²⁵ g cm s^-2, would correspond to that of a red giant, although if a red giant wind had been simulated, the velocity would have been an order of magnitude lower (with Ṁ scaling accordingly), making the simulation much longer. However, since the relation between the jet momentum flux and the wind thrust determines to first order the shape of the interaction region, a lighter wind of equal thrust can be used to reduce the computational costs.

The star was assumed to be at rest. This is a reasonable assumption as long as the stellar velocity is much lower than v_sw. Otherwise, the jet-wind interaction geometry will strongly depart from axisymmetry, making the results obtained here less realistic. The Keplerian velocity for a 10⁸M_⊙ central black hole at a distance of 10 pc is v_K ≈ 2 × 10⁷ cm s^-1, a 10% of the adopted v_sw-value. We considered this v_K-value low enough at this stage, but a caveat must be made: for a more realistic red giant wind, v_K would become of the order of v_sw, and the star motion would then have to be taken into account. This requires 3D simulations, however, which are much more computationally expensive than the axisymmetric simulations carried out here, which are meant as a first step in coupling radiation and hydrodynamics studying jet-star interactions (we note, nevertheless, that work is being carried out in this direction). Finally we note that in general the jet-crossing time will be much longer than the simulated times (given in Sect. 2.2).

The grid was initially filled with the jet properties except within the stellar wind injection region, which was filled with the stellar wind properties. The jet and stellar wind physical parameters are summarized in Table 2.

The jet momentum flux (thermal plus kinetic pressure) at the bottom of the computational grid is (2)and the stellar wind momentum flux at a distance R is(3)where S_sw = 4πR², and is the pressure of the stellar wind. Equations (2) and (3) allow us to set the point where the jet and the stellar wind momentum fluxes are equal, at R_CD = 0.1 × 10¹⁵ cm from the star centre, locating the contact discontinuity (CD) on the simulation axis. This corresponds to a position z_CD = z₀−R_CD = 0.2 × 10¹⁵ cm with respect to the bottom of the grid (z = 0).

2.1.2. Streamline calculations

Streamlines of the jet flow were characterized in the numerical solution as described in Appendix A. An important parameter determining the streamline magnetic field is χ_B, the initial ratio of the Poynting-to-matter energy flux. The streamlines were divided into segments or cells. Each cell represents an annular element from the point of view of their volume and the determination of the non-thermal energy budget because the simulation is axisymmetric. On the other hand, the cells represent point-like non-thermal emitters in the radiation calculations below. To provide the 3D information of the structure of the whole non-thermal emitter, an azimuthal angle ψ was assigned to each cell, as shown in Fig.1, and the value of this angle was varied from 0 to 2π to cover the whole emitting volume. Figure 1 provides a sketch of the cell geometry in both contexts, the simulation results, and the radiation calculations.

2.2. Results

2.2.1. Steady-state

Figure 2 shows the density map and the trajectories of the computed streamlines for the jet-stellar wind simulation when (quasi-) steady-state is reached at t = 3.3 × 10⁷ s. The dynamical timescale of the problem is determined by the stellar wind, which is ~l_z/v_sw = 7.5 × 10⁶ s, as it carries more mass and v_sw ≪ v₀. Thus, the t-value reached seems appropriate. The two shocks formed are bow-shaped towards the star, which has a much lower thrust than the jet within the grid. The CD in the simulation is at cm, as obtained analytically in Sect. 2.1.1. Figure 3 illustrates the re-acceleration of the shocked jet material as it is advected upwards. This feature has a strong effect on the radiation because Doppler-boosting effects cannot be neglected. This is also shown in Fig. 4, where the Doppler-boosting enhancement of the emission in the jet direction is presented.

2.2.2. Instability growth

The simulation reaches a quasi-stationary numerical solution, with the shocked flow structure in a metastable state. There are recurrent perturbations coming from the numerical, spatial, and temporal discretization that grow as a result of the developing instabilities. Although the perturbations are of numerical origin, they can be considered to mimic the irregularities expected in real flows because they are hardly completely smooth or laminar. The double-shock structure presents variations in time caused by irregularities originated in the CD that grow as they are advected with the flow. The growth of these irregularities is mostly linked to the Kelvin-Helmholtz instability (KHI) because the velocity of the shocked jet flow along the CD grows very quickly from the simulation axis, which leads to a strong velocity difference with respect to the shocked wind; if the velocity of the CD perpendicular to itself were not zero, the Rayleigh-Taylor instability (RTI) would develop as well (Chandrasekhar 1961). In our case, the average velocity of the CD was zero (although there are small fluctuations in its position). The KHI is thus the dominant source of perturbation growth in the CD close and far from the axis.

The development of instabilities has imposed a limitation on the simulation parameters. At the onset of the simulation, instabilities develop in the interacting flows, with the related growing irregularities being advected out of the grid and leaving the quasi-steady interaction structure described above. For lighter jets, heavier winds, or a higher resolution, the development and growth of the perturbations are enhanced and produce smaller and denser fragments of stellar wind able to penetrate deeper into the jet flow and generate shocks in the entire computational domain. This is caused by instability growth, which causes small wind structures to quickly develop and propagate within the grid, triggering shocks in the jet flow. This can lead to a grid that is partially occupied by hot material that reaches the grid boundaries with subsonic velocities. This renders the simulation results unrealistic because waves develop in the grid boundaries and bounce back, which affects the flow dynamics inside the grid. In addition, since the flow is subsonic at the grid boundaries, the subsonic flow is evacuated too slowly and accumulates and can potentially end up filling the whole grid. This could be avoided by a larger grid, although if a higher resolution were used, the disruptive effect of instabilities (see e.g. Perucho et al. 2004) would be enhanced, pushing the grid size requirements even further.

Strong sensitivity to resolution and flow density contrast was already faced in previous similar axisymmetric simulations. The fast growth of perturbations close to the axis was seen for instance in Paredes-Fortuny et al. (2015) in the context of a relativistic pulsar wind interacting with a non-relativistic stellar wind. It was noted that this effect might have a partially numerical origin in the coordinate singularity plus a reflective boundary at r = 0. However, a similarly fast perturbation growth also appeared in relativistic 2D simulations for that scenario (Bosch-Ramon et al. 2012a; Lamberts et al. 2013) in planar geometry, which shows that the perturbations were not only due to an artefact of the conditions at r = 0. For the jet-stellar wind scenario, we ran low-resolution 3D simulations (not presented here) focusing on the region close to the star. We found that the shocked structure was also prone to develop instabilities, although the small grid size prevented a deeper analysis.

2.2.3. Perturbed state

Judging from the discussion above, this growth seems to be a physical effect to a large extent, although the fast growth of the instabilities may be partially linked to the axisymmetry of the simulations. It is thus worth considering some instance of the interaction region when it is affected by a strong perturbation before reaching steady-state. Such an instance is shown in Fig. 5, which presents the density map and trajectories of the computed streamlines at t = 3.7 × 10⁶ s. As seen in the figure, the perturbations deeply penetrate the jet, increasing the size of the interaction region and strongly modifying the streamline trajectories. This is an example of how important physical instability growth might be for the effective size of the stellar target the jet meets.

In the corner of the map in Fig. 5 where the top and right boundaries join, a reflection shock is visible in a small region. This reflection shock is an example of the presence of subsonic flow at the boundaries. However, this shock only affects a minor region of the simulation, even accounting for the cylindrical symmetry, so we have allowed for the presence of this small artefact. At our level of resolution, our simulation did reach a physically realistic steady-state (with the mentioned small fluctuations in the CD). Although the hydrodynamical solution shown in Fig. 5 does not correspond to steady-state because radiation is computed from the shocked jet material with a dynamical timescale ~l_z/c ≪ l_z/v_sw, we can assume that the flow is in a pseudo-steady state for emission computation purposes.

When characterizing the non-thermal emitters through streamlines, we discarded the lines that presented mixing between jet and stellar wind fluids that exceeded 50% in at least one cell of the streamline for the jet-star interactions in steady-state. The surface of any discarded streamline was added to the line immediate adjacent along the radial coordinate. This compensates for the lack of lines close to the jet axis because the hydrodynamical conditions are relatively similar there. On the other hand, the streamlines of the jet-star interaction in the perturbed state were computed regardless of the level of mixing because many lines were numerically affected. For this case, we just removed the segments of the streamlines that were affected by a mixing level >50%.

3. Radiation

3.1. Non-thermal emitter

After obtaining the streamlines from the RHD simulations, we used the hydrodynamic information to compute the injection, cooling, and radiation of the non-thermal population. The computational domain was divided into cells as described in Sect. 2.1.2. Each cell was considered homogeneous and characterized by its position and velocity vector information, pressure (P), density (ρ), section (S), magnetic field (B), and the flow velocity divergence (∇(Γv)) to compute adiabatic losses. A parameter accounting for wind mixing, going from 0 (100% jet material) to 1 (100% stellar wind), was also taken into account through the use of a tracer variable that is included in the RHD code.

As described in Appendix B, the code computes for each streamline (i) which cells have non-thermal particle injection as well as the luminosity injected into these non-thermal particles; and (ii) the steady particle energy distribution in each cell, N(E). The non-thermal particle injection was assumed to occur when the internal energy density and entropy grows, which is the case for shocks in our hydrodynamical approach (in the case of resistive magnetohydrodynamical simulations, magnetic reconnection would also lead to an increase of the internal energy and entropy). Another important parameter is the fraction of energy in non-thermal particles, χ_NT, which we fixed here to a modest 0.1, which was taken as a reference value. The radiation levels need to be scaled proportionally when higher or lower acceleration efficiencies are adopted. Note, however, that high values of χ_NT ≲ 1 would not be consistent with our hydrodynamical calculations if radiation cooling dominated (although in most of the instances computed in this work this is not the case).

After obtaining the particle energy distribution in each cell, as explained in Sect. 3.1.1, the code derives the synchrotron and inverse Compton (IC) radiation as seen from the observer. When the radiation coming from each line is known, the total spectral energy distribution (SED) can be obtained and emission maps derived. The code includes the possibility of accounting for time delays in the emission from non-stationary flows, but here this is not required as the emitting flow can be considered steady.

The hydrodynamic simulations are axisymmetric and 2D, but the actual emitter is a 3D structure. To properly compute Doppler boosting or obtain radiation maps for non-trivial geometries, the 1D emitters were therefore equally distributed azimuthally with respect to the symmetry axis, conserving the radial and vertical position and velocity components. The higher the number of 1D emitters, the better the sampling of the 3D emitter. It is enough to equally distribute only the emitting cells to calculate Doppler boosting (but not maps), which somewhat reduces the computational costs. This transformation only affects the radiative part of the code; the particle energy distributions in each cell or streamline are unaffected.

3.1.1. Emission

The luminosity per frequency unit (spectral luminosity) and per solid angle of the synchrotron emission in the FF for each cell and in the optically thin case was computed following (Pacholczyk 1970) (4)where N(E^∗) is the cell particle energy distribution, E^∗ the particle energy, , the critical energy, , , h the Planck constant, and θ′ the angle between B′ and the direction towards the observer. The function F(x) is defined through an integral of the K_5/3(z) Bessel function, but we adopted the following approximation (e.g. Aharonian 2004), which is valid in the interval 0.1 <x< 10: (5)The orientation of B′ is not fully defined in our approach, and for simplicity we made the approximation . As particles move isotropically in the FF, the SED in that frame can be computed simply through (6)where L′(ϵ′) is the spectral luminosity of each cell in the FF.

The IC radiation was computed in the FF using the kernel for the emission rate of electrons of energy E′ interacting with monodirectional target photons of energy with an angle θ′ (e.g. Aharonian & Atoyan 1981; Aharonian 2004; Khangulyan et al. 2014): (7)where , and z = ϵ′/E′. Convolving Eq. (7) with N′(E′) and the target photon energy distribution density in the FF, , the IC SED can be obtained: (8)A black body, centred at the star location, (0,0.3 × 10¹⁵) cm, was adopted to obtain n₀(ϵ₀), which in this case , with a temperature and luminosity of T_∗ = 3 × 10³ K and L_∗ = 3×10³⁶ erg s^-1, respectively, typical for a red giant. Other photon fields such as the cosmic microwave background, the radiation from an accretion disk, the starlight from the galaxy, or radiation from the jet itself, were not taken into account because we focus here on the dominant photon field on the spatial scales of the calculations: the stellar photon field. If the non-thermal particles were followed up to cover larger spatial scales, however, other photon fields might be dominant and therefore should included in the calculations (see Sect. 5).

The derived SEDs were obtained for each cell in the FF. Owing to relativistic and light retardation effects, the photon energy and SED as seen by the observer have to be transformed as ϵ = δϵ′ and ϵL(ϵ) = δ⁴ϵ′L′(ϵ′), where δ = 1/Γ(1−βcosθ_obs), with θ_obs being the angle between the flow and the observer direction in the laboratory frame (LF). For the IC, the angle θ′ is also to be transformed from the LF. Further details on the transformation of the relevant angles can be found in Appendix C. For a dense, external target photon field, the gamma rays produced might be absorbed through pair creation (functionality included in the code). Electromagnetic cascading in the emitter environment may be also relevant in some situations. In the present scenario opacities are well below 1 in the considered cases, therefore gamma-ray absorption is negligible.

3.2. Results

We applied the radiative code described in Sect. 3.1 in conjunction with the RHD results of Sect. 2 to compute the radiation in the two different scenarios: the star-jet interaction in the steady and the perturbed states. We let three parameters vary: the angle between the jet axis and the line of sight, φ; the initial ratio of Poynting-to-matter energy flux (roughly the ratio of magnetic-to-thermal pressure just downstream of the jet shock), χ_B; and the height in the jet, with respect to the jet base, where the jet-star interaction takes place: z_int. The values considered for these parameters are summarized in Table 3. We recall that the simulations were computed with a purely RHD code, meaning that the magnetic field is not dynamically relevant and only high enough to allow the plasma to behave as a fluid. However, we adopted χ_B = 10^-1 in some cases, which is still small enough to avoid violating the low-B condition, to have an example of a relatively high B-case.

3.2.1. Scalability of the results

Throughout the paper, the jet power, L₀, opening angle (taken ~1/Γ₀, as in Bosch-Ramon 2015, see references therein), and wind thrust, Ṁv_sw, were fixed to the following reference values: ≈10⁴⁴ erg s^-1, 0.1 rad, and 1.3 × 10²⁵ dyn, respectively. The star temperature and luminosity were also fixed to those of a red giant (see Sect. 3.1.1). However, the obtained results for the different values of φ, χ_B, and z_int explored can be easily generalized if Γ₀ and T_∗ are fixed, and L_∗ is approximated as ~Ṁv_swc (Bosch-Ramon 2015). An additional assumption to perform the generalization is that the radiation on the particle energy distribution (e.g. through synchrotron self-Compton) or the radiation itself (e.g. through internal pair creation) do not significantly interact. Under these conditions, the SEDs and mapped quantities presented below can be scaled as follows:

• The spectrum does not change if the escape-to-radiativetimescale ratio, is constant, where and t_esc ∝ R_CD. The quantity t_rad typically depends on the particle energy, which means that the timescale ratio has to be computed by fixing this energy to some particular value. The quantity t_esc would correspond to the typical timescale required for particles to escape the emitting region (see Bosch-Ramon 2015).

• The SED normalization changes as (Bosch-Ramon 2015) where R_j is the jet radius. In the adiabatic regime of the emitter, the SED normalization is also ∝t_esc/t_rad.

• Thus, in the adiabatic regime, the non-thermal luminosity, either synchrotron or IC, can be simply scaled as which is the product of the dependences stated in (ii).

These relations are more refined versions than, but based on, those presented in Eqs. (6) and (9) in Bosch-Ramon (2015).

For the jet Lorentz factor Γ₀, Bosch-Ramon (2015) indicated that the normalization of the observer luminosity should approximately scale as . However, this is strictly valid in the ultra-relativistic regime, and far from the shock and/or under fast expansion of the streamlines; the scaling is less sensitive to Γ₀ for relatively slow expansion and in the simulated region.

3.2.2. Spectral energy distributions and radiation maps

Figures 8 and 9 show the observer synchrotron and IC SEDs for the star-jet interaction in the steady-state, for χ_B = 10^-4 and 10^-1, φ = 0°−135°, and z_int = 10 pc. For high magnetization, synchrotron emission strongly dominates and reaches much higher photon energies. However, even in this case, most of the particle energy distribution is dominated by advection escape, and therefore the synchrotron and IC SED shapes do not depend significantly on χ_B. The advection escape dominance also implies that non-thermal particles behave as an adiabatic flow and that the ratio of non-thermal-to-thermal energy will approximately keep constant along the streamlines. As seen in the figures, the emission is significantly boosted for the jet on axis. Interestingly, for χ_B = 10^-1 the synchrotron SEDs for high φ-values present softer spectra because particles coming from the more Doppler-boosted outer lines (thus beamed away for high φ) have higher maximum photon energies. This effect is not as clearly seen in the IC SED because the softening of the SED is already strong as a result of the Klein-Nishina (KN) effect in the cross section.

In Fig. 10 the observer synchrotron and the IC SED are shown for different values of z_int from 1 to 1000 pc, and for φ = 0°, χ_B = 10^-4, and φ = 0°. It is clear from the figure that the closer to the jet base, the higher the ratio of radiative versus advection escape, which as mentioned is ∝1 /z_int.

Figure 11 presents the distribution of the bolometric luminosity per cell in the rz-plane for both synchrotron and IC, taking χ_B = 10^-4, z_int = 10 pc, and φ = 0°. We note that because of the azimuthal geometry, the computational cells correspond to annular physical regions. The maps show that the synchrotron emission is more widely distributed than IC emission. This is because the IC target photon field is concentrated towards the star.

Finally, in Figs. 12 and 13 we show the contribution of the different lines for z_int = 10 pc, χ_B = 10^-4, and φ = 0°, to the observer synchrotron and IC SED, and total energy distribution in the LF, respectively. The contribution to the emission varies substantially between different streamlines, depending on the non-thermal particle content, the role of adiabatic cooling or heating along the lines, the flow velocity and direction, the local magnetic field, and the relative position with respect to the source of target photons.

Figure 14 shows the observer synchrotron and IC SEDs for the star-jet interaction in the perturbed state, for χ_B = 10^-4 and 10^-1, φ = 0°−135°, and z_int = 10 pc. Figure 15 presents the distribution in the rz-plane of the bolometric luminosity per cell for both synchrotron and IC, taking χ_B = 10^-4, z_int = 10 pc, and φ = 0°.

Figure 16 shows the synchrotron and IC SEDs for the two cases studied in this work, in the low- and high-magnetization case and for φ = 0° and z_int = 10 pc. The figure shows that the synchrotron emission is significantly higher for the star-jet interaction in the perturbed state because a larger section of the jet is affected by the stellar wind, which is apparent from comparing Figs. 2 and 5 in Sect. 2. This implies that more jet energy available for radiation. On the other hand, the IC radiation levels change very little when compared to the synchrotron levels because the target photon density significantly drops with distance from the star. In the high-magnetization case, the differences between the two studied cases are smaller than for a low magnetization, and the synchrotron SED of the perturbed case shows a softer spectrum at the highest energies. This is most likely related to the high magnetic field, which through severe synchrotron cooling prevents the most energetic electrons from reaching regions of higher Doppler-boosting. Otherwise, the radiation of these particles would have led to more flux at the higher synchrotron energies and thus to a harder spectrum.

4. Conclusions

The results of the radiation calculations carried out using hydrodynamical information are in the line of those presented in Bosch-Ramon (2015). In that work, it was noted that the effective size of the obstacle is much larger than the distance at which this and the jet flow collided (i). It was also pointed out that Doppler-boosting plays a significant role (ii). All this is confirmed here:

• (i)

  For instance, for the jet-star interaction in steady-state withφ = 0°, z_int = 10 pc, and χ_B = 10^-4, the total observer luminosity is ≈3.5 × 10³¹ and 5 × 10³³ erg s^-1 for the synchrotron and IC components, respectively, whereas the synchrotron and IC luminosities from a region r<R_CD are ≈2 × 10²⁹ and 8 × 10³¹ erg s^-1, respectively, which is a factor ~100 smaller than for the whole grid.

• (ii)

  In addition, the total synchrotron and IC luminosity in the FF are ≈1.3 × 10³⁰ and 2 × 10³² erg s^-1, respectively, which is a factor ~30 smaller than the observer values when Doppler-boosting is not accounted for. It is worth noting that jets on kpc and larger scales may be less relativistic, as considered in Bednarek & Banasiński (2015), but closer to the galaxy centre, they are likely to have much higher Lorentz factors. Farther out, jet-boundary instabilities, shear-layer development, recollimation shocks, and mass-loss from stellar winds work together to decelerate the jets (see e.g. Perucho 2014, 2016).

As expected from previous work (e.g. Araudo et al. 2013; Bednarek & Banasiński 2015; Bosch-Ramon 2015), we find that advection escape dominates radiation losses in the star-jet interactions studied, that is, for moderately powerful jets and stellar winds. Radiation peaks from X-rays to MeV energies depending on χ_B for synchrotron emission, and in the 100–1000 GeV range for IC, with the cooler, either older or less massive, stars yielding the higher SED peak energy (~1 TeV, as in the present work). A study of collective interactions of many stars with the jet (see Araudo et al. 2013, for massive stars) is under way to also account for other stellar populations such as evolved stars and for the increased effective section of the obstacles and relativistic effects. However, Bosch-Ramon (2015) predicted significant collective emission, for which they took the stellar populations of the inner regions of the radio galaxy M 87 as a reference. These predictions receive further support here because we obtain similar IC fluxes and general behaviour for individual interactions, although we note that a more detailed prescription of the stars and their winds in the inner regions of AGN is required.

Instabilities might play an important role, and the shocked flow structure might be relatively unstable. This may enhance the emission for some individual jet-star interactions, although the collective emission would likely be at about steady-state radiation levels. Otherwise, for interactions taking place in the innermost jet regions and bright enough to be detectable by themselves, instability growth could temporary increase the emission levels and induce variability on scales ≳R_CD/v_sw in addition to other types of variability, such as that associated with jet crossing. It is worth noting that a longitudinal magnetic field in the jet would reduce the instability growth.

5. Final remarks

The question whether electrons or e^±-pairs are indeed accelerated in the jet shock is a key point. Bosch-Ramon (2015) showed that for luminosities per interaction similar to those obtained by us and under the same interaction conditions, collective jet-star interactions could yield detectable levels of emission in M 87, for instance. Therefore, our results together with those from Bosch-Ramon (2015) indicate that for acceleration efficiencies of χ_NT ≳ 0.1, collective star-jet interactions may be detectable in AGN in gamma rays unless (i) an AGN jet is a stronger emitter through a different mechanism; and/or (ii) the source is far and the jet seen off-axis; and/or (iii) the jet power is low. It is worth noting that the need of a relatively high χ_NT-value applies particularly to the low-B scenario because for higher magnetizations the synchrotron component could overcome the IC in the 100 MeV region, which would weaken this requirement.

The limited size of the grid makes accounting for all the emission produced in the interaction region difficult, in particular for the synchrotron radiation. Bosch-Ramon (2015) noted that non-negligible kinetic to internal energy transformation takes place far from the simulation axis, although with a weak r-dependence. However, Doppler-boosting effects are expected to be strong far from the obstacle, and synchrotron emission can come from farther regions than IC because the latter is affected by stellar field dilution with distance from the star. Larger grid simulations or jet-scale semi-analytic calculations are required to determine the emission contribution of these farther regions. As a reference, we point out that for the case with φ = 0°, z_int = 10 pc, and χ_B = 10^-4, ~10% of the synchrotron and 30% of the IC total luminosity come from a distance 3 × R_CD from the simulation axis. Therefore, the trend of the IC emission indicates convergence, but the synchrotron emission may be still far from that. This is linked to the fact that our simulations yield very low adiabatic cooling rates. The dominant source by far of non-thermal energy loss is advection escape from the grid. This is consistent with the modest density decrease with z in the shocked jet flow, as the density maps in Sect. 2.2 illustrate. Our results are complementary to those in Bednarek & Banasiński (2015), where jet emission from particles accelerated in jet-obstacle interactions was computed by accounting for the jet B and photon fields that are relevant on larger spatial scales, such as the galaxy bulge or the CMB. A larger grid would also help to study more unstable configurations, such as adopting a higher density wind-jet contrast or increasing the resolution for a more realistic setup, although the computation time required would severely increase. Finally, 3D simulations that account for the star motion also need to be carried out.

---

Acknowledgments

We acknowledge support by the Spanish Ministerio de Economía y Competitividad (MINECO) under grants AYA2013-47447-C3-1-P, and MDM-2014-0369 of ICCUB (Unidad de Excelencia “María de Maeztu”). This research has been supported by the Marie Curie Career Integration Grant 321520. V.B.-R. also acknowledges financial support from MINECO and European Social Funds through a Ramón y Cajal fellowship. X.P.-F. also acknowledges financial support from Universitat de Barcelona and Generalitat de Catalunya under grants APIF and FI (2015FI_B1 00153), respectively. D.K. acknowledges financial support by a grant-in-aid for Scientific Research (KAKENHI, No. 24105007-1) from the Ministry of Education, Culture, Sports, Science and Technology of Japan (MEXT). M.P. is a member of the working team of projects AYA2013-40979-P and AYA2013-48226-C3-2-P, funded by MINECO.

References

Appendix A: Streamline calculations

The simulations are axisymmetric, and by assuming that the flow is stationary, which is approximately valid for most of the jet material within the grid, we can therefore introduce a stream function, Φ, as (A.1)The entire flow can be divided into a set of regions, V_i, based on the value of the stream function: Φ_i−1< Φ(r,z) ≤ Φ_i, where 1 ≤ i ≤ N and Φ₀ = 0. Physically, each of these regions fulfils the condition that the same amount of matter flows through their cross section. Since at the bottom boundary of the computational domain the quantities v_r and Γρv_z are constant, with v_r = 0, the stream function can be computed on this surface from dΦ = 2πrΓρ(v_zdr−v_rdz) = Γ₀ρ₀v₀dS₀, where dS₀ = 2πrdr is a surface element. This means that at the bottom of the computational domain the cross sections of the regions V_i correspond to a sequence of annular rings: r_i−1<r ≤ r_i, with r₀ = 0, r_N = l_r, and .

These regions V_i form layers with axial symmetry in the 3D space. If their thickness is small enough, they can be considered homogeneous in the directions perpendicular to the fluid motion, and the properties of the non-thermal particles in each layer will depend only on the time spent by the fluid element of interest since its injection into the computational domain. In particular, the number density and energy distribution of the non-thermal particles responsible for the synchrotron and IC emission of the jet in a region V_i depend only on this time. To compute the energy distribution numerically, we selected a streamline for each volume element V_i and solved the particle transport equation along this line using the momentarily snapshot of the plasma properties obtained with our hydrodynamic simulations (for details see Appendix B).

Each streamline was divided in segments by an equidistant time step Δt: (A.2)where Δt is the smallest cell size divided by the highest speed in the simulation. The initial point for the streamline in the region V_i was selected as r_s,0 = (r_i−1 + r_i)/2 and z_s,0 = 0. This division of the streamlines leads to the splitting of the regions V_i into annular cells. The cross sections transversal to the flow velocity of the cells can be computed assuming conservation of the energy crossing those sections per time unit, (A.3)where h and v are the specific enthalpy and the modulus of the three-velocity, respectively. The subscript 0 denotes the conditions at the bottom boundary of the region V_i, that is, in particular, (A.4)After the next point along the streamline, (r_s,j,z_s,j), was obtained, the hydrodynamic parameters (v_r,j, v_z,j, ρ_j, and p_j) were computed through a bilinear interpolation of the corresponding values from the neighbour cells (Press et al. 1992). This was done iteratively until the fluid element followed along the streamline left the simulation grid.

The magnetic field at the inlet of the computational domain in the FF, , was calculated assuming that the Poynting flux is equal to a fraction χ_B of the matter energy flux, which yields (A.5)Then, from Eq. (A.5), the magnetic field at each cell of the streamline can be derived assuming that it is frozen into the plasma and is perpendicular to its direction of motion: (A.6)Finally a downsampling of each streamline through interpolation was performed to reduce the computation time of the radiative code, resulting in streamlines sampled with up to 200 equidistant cells.

Appendix B: Non-thermal particles

Appendix B.1: Injection of particles

Shocks are assumed here to be the sites of particle acceleration. Using streamline data, we can compute the energy rate of accelerated particles (or non-thermal luminosity) injected in cells were shocks occur. Injection is thus assumed to take place when there is an increment in the internal energy and a decrement in the flow velocity between two consecutive cells to approximately capture shocks. The rate of energy injection in the form of accelerated particles in the cell volume in the FF, , is taken to be a fraction χ_NT ≤ 1 of the generated internal energy per second in the cell k (), that is, the time derivative of the 00-component of the energy-momentum tensor times volume in the FF: (B.1)where v_+/− and Γ_+/− are the velocities and Lorentz factors in the cell right and left boundaries, respectively. Equation (B.1) is valid for a steady-state fluid.

The normalization of the source function for non-thermal particles, Q′, should satisfy the following condition: (B.2)At this stage, the expression we assume for Q′ is a power-law of index −2, typical for shock acceleration in different contexts (see below), with two cutoffs at high and low energies: (B.3)The steep exponential cutoff towards low energies is adopted to avoid numerical artefacts at the lowest energies, and also to make the number of transrelativistic particles small. The low-energy cutoff energy is arbitrarily fixed here to MeV.

In addition to the fraction of energy going to non-thermal particles χ_NT, it is also required to set the timescale at which particles gain energy when accelerated to determine the highest particle energy, since the high-energy cutoff is fixed to the energy at which the acceleration timescale is equal to the shortest cooling or escape timescale. For the acceleration timescale, we assume that the process is linked to the jet shock and adopt a phenomenological prescription with η = 2π(c/v)², which tends to ~10 when v → c. This is done for simplicity because particle acceleration in relativistic shocks is at present a complex and far from settled matter, and the simple prescription is enough for our purposes. The prescription is loosely based on the acceleration timescale for non-relativistic strong shocks (e.g. Drury 1983), although in our prescription v refers to the flow velocity in the laboratory reference frame (LF), and not to the shock velocity in the reference frame of the medium upstream of the shock (as in particle acceleration theory). The cooling or escape timescales are given by the synchrotron timescale (dominant at the highest energies), (a_s = 1.6 × 10^-3, cgs units¹); and the Bohm diffusion timescale, , with being the typical size of the emitter, and D′ = cE′/ 3qB′ being the Bohm diffusion coefficient. In this work, is fixed to 3 R_CD. The high-energy cutoff can be obtained by combining the different mentioned timescales: (B.4)When the internal energy increases, we set ∇(Γv) to zero because in the injected luminosity we already take the adiabatic heating that would take place into account. Formally, non-thermal particle injection should take place only at shocks, and adiabatic cooling or heating would operate otherwise. We choose the adopted approach for simplicity, however, as defining a shock accurately in streamlines is numerically more complex and demanding.

The emission results depend relatively weakly on the particle injection assumptions. This is true in particular for the IC luminosity in gamma rays, which is not significantly affected by E_c,high as long as it is ≳100 GeV, which is not a very demanding value. The relevant parameter is χ_NT, which is directly proportional to the normalization of the injected particle distribution and is not constrained.

Appendix B.2: Particle energy distribution

To compute the non-thermal particle evolution, we considered the hydrodynamical information constant over a time , much shorter than the dynamical timescale of the simulations. During this time, we let the particles evolve until they reach a stationary solution for their energy distribution. This calculation is made in many short time steps: (B.5)where the subindex i relates to the time step at time in the LF, k to the cell, and min [Δx_k,i/v_k,i] is the shortest cell crossing time at , with Δx_k,i being the cell length. To compute the particle energy distribution at each (LF time) and cell, the code, first order in space, time, and energy, follows three stages:

• First, the energy distribution after a of the particles injected at in the cell k, in the FF, is obtained from (dropping the indexes i and k): (B.6)where is the energy that particles had before advancing the FF time a , and is given implicitly by (B.7)Equation (B.6) is the solution of the equation adapted from Ginzburg & Syrovatskii (1964) to compute the evolution of particles in an homogeneous region, (B.8)where Ė′ is the particle energy-loss rate including all the relevant losses in the FF: IC (e.g. Khangulyan et al. 2014); synchrotron; and adiabatic cooling. The adiabatic losses are computed using the divergence of the spatial part of the flow four-velocity, (B.9)

• Second, even if no particles are accelerated in a cell, there are still particles that had come at from the previous cell following an energy distribution . These particles evolve as (B.10)under steady cooling conditions, .

• After particles are evolved in energy, flow advection is taken into account to include the contribution from those evolved particles arriving from the previous cell, and remove the locally evolved particles that flow to the next one. In the cell k the advection effect is (B.11)where (B.12)and and are the fraction of particles of a given cell, k−1 and k, respectively, which left that cell.

The computing tool described here can deal with non-linear processes, such as synchrotron self-Compton (SSC) because when the synchrotron radiation is computed, the particle distribution can be recalculated using a new ambient photon field that includes this component. Internal pair creation could be also taken into account as an additional source of injected particles, but this functionality has not been implemented yet. At present, nevertheless, these processes are not relevant in the studied scenario: with respect to IC on locally produced radiation, as SSC, IC on the external photon field of the star dominates these other IC channels. On the other hand, pair creation through photon absorption in external or local photon fields has a very low probability and can be neglected.

Appendix C: Transformation of angles to the fluid co-moving reference frame

Several sources of target photons can play an important role in inverse Compton scattering. However, given the small size of the computational domain, we here only considered the contribution from the star, which has been approximated as a point-like emitter. Therefore, IC scattering proceeds in the anisotropic regime (see e.g. Bogovalov & Aharonian 2000), and we have to account for several effects that involve relativistic transformations of angles. In particular, particle cooling and scattering occur in the FF, which means that to compute these processes, we need to transform the photon field and scattering angle to the FF.

We consider a 2D hydrodynamic flow and a coordinate system with the z-axis coinciding with the symmetry axis. The source of target photons is located at r_∗ = (0,0,z_∗) and the observer is in the rz-plane, that is, the line of sight is parallel to (C.1)A fluid element located at r = (rcosψ,rsinψ,z) (with r, ψ and z being cylindrical coordinates) moves with four-velocity (C.2)along the direction represented by (C.3)Here, is the fluid element Lorenz factor. Particles in the flow are illuminated by photons with velocity directed along (C.4)For the non-thermal emission calculations, the following angle-related parameters are required:

• (i)

  The Doppler factor is given by the following expression:(C.5)where .

• (ii)

  To compute the IC losses, we need to transform the photon field to the FF. This transformation is determined by the Doppler factor of the stellar photons (Khangulyan et al. 2014): (C.6)

• (iii)

  Gamma-gamma absorption on a field provided by a point-like source is determined by the distance to the source of target photons and by the angle between n_obs and n_ph (Khangulyan et al. 2008): (C.7)

• (iv)

  Finally, to compute anisotropic IC, the scattering angle in the FF is needed. To derive , let k = (k,kn_obs) be the momentum of a photon propagating towards the observer, and ω = (ω,ωn_ph) the momentum of a photon emitted by the star towards the fluid element. Since the scalar products (ku) and (ωu) are Lorentz invariant and the four-velocity of the fluid element in the FF is u′ = (1,0), we obtain that (C.8)and (C.9)Since (kω) is an invariant, (C.10)and consequently, (C.11)

All Tables

All Figures