Relativistic jet feedback – III. Feedback on gas discs

ABSTRACT

We study the interactions of a relativistic jet with a dense turbulent gaseous disc of radius ∼2 kpc. We have performed a suite of simulations with different mean density, jet power, and orientation. Our results show that: (A) The relativistic jet couples strongly with the gas in the inner kpc, creating a cavity and launching outflows. (B) The high pressure bubble inflated by the jet and its back-flow compresses the disc at the outer edges, driving inflows within the disc. (C) Jets inclined towards the disc affect the disc more and launch sub-relativistic, wide-angled outflows along the minor axis. (D) Shocks driven directly by the jet and the jet-driven energy bubble raise the velocity dispersion throughout the disc by several times its initial value. (E) Compression by the jet-driven shocks can enhance the star formation rate in the disc, especially in a ring-like geometry close to the axis. However, enhanced turbulent dispersion in the disc also leads to quenching of star formation. Whether positive or negative feedback dominates depends on jet power, ISM density, jet orientation with respect to the disc, and the time-scale under consideration. Qualitatively, our simulations compare favourably with kinematic and morphological signatures of several observed galaxies such as NGC 1052, NGC 3079, 3C 326, and 3C 293.

hydrodynamics, methods: numerical, galaxies: evolution, galaxies: high-redshift, galaxies: ISM, galaxies: jets

1 INTRODUCTION

Outflows from active galactic nuclei (AGN) have long been identified to be one of the dominant mechanisms that influences stellar mass assembly in massive galaxies (Silk & Rees 1998; Di Matteo, Springel & Hernquist 2005; Bower et al. 2006; Croton et al. 2006). Feedback from the central supermassive black hole is thought to occur primarily in two ways. In the quasar mode (or the establishment mode, Wagner et al. 2016), outflows are driven by luminous AGN accreting at near Eddington rates, and disperse star-forming gas from the galaxy. Such winds are primarily considered to be wide-angled, driven by the radiation pressure from the luminous AGN (e.g. Proga 2000; Ostriker et al. 2010; Bieri et al. 2017). In the radio mode (or maintenance mode; Harrison 2014), it is envisaged that large-scale jets from the central black hole heat the IGM, thereby preventing cooling flows and mass build-up of galaxies (Fabian 2012). However, it is still unclear how the central AGN affects the gas in the galaxy, and over what time-scales (Schawinski et al. 2009, 2015). On-going star formation in AGN host galaxies (Harrison et al. 2012) also complicates the simple interpretation presented above, motivating the need for more detailed investigation of how the central black hole may affect the ISM of the galaxy.

Observational evidence for both forms of feedback is plentiful. Powerful radio jets extending to ∼100 kpc and heating the ambient gas, especially in clusters of galaxies (Fabian 2012), are an archetypal maintenance mode scenario. AGN powered galactic outflows have been detected in several systems (Nesvadba & Lehnert 2008; Rupke & Veilleux 2011; Tombesi et al. 2015; Collet et al. 2016; Nesvadba et al. 2017). Suppressed star formation in systems harbouring an AGN-driven outflow (Nesvadba et al. 2010; Wylezalek & Zakamska 2016) are an indicator of establishment mode feedback. The role of jets from supermassive black holes has largely been considered in the context of maintenance mode feedback (Gaspari, Brighenti & Temi 2012a; Gaspari, Ruszkowski & Sharma 2012b; Yang & Reynolds 2016). It is however not well understood whether galaxies with radio jets are also directly affected by the jet itself, before they evolve into large-scale structures.

Recent observational studies (Nesvadba et al. 2010, 2011; Alatalo et al. 2015) have shown that star formation may be significantly regulated or suppressed by turbulence and outflows induced by the radio jet. However, recent theoretical (Gaibler et al. 2012; Bieri et al. 2016; Fragile et al. 2017) and observational (Lacy et al. 2017; Salomé, Salomé & Combes 2015; Bicknell et al. 2000) papers have proposed a jet induced positive feedback scenario where jet-driven shocks may induce star formation in galaxies. Hence, the manner in which a relativistic jet, often well collimated at radio wavelengths, affects the large-scale interstellar medium (ISM) of the galaxy is not well understood. This is an open question, requiring both theoretical and observational investigation.

In earlier papers (Sutherland & Bicknell 2007; Wagner & Bicknell 2011; Wagner, Bicknell & Umemura 2012; Mukherjee et al. 2016), we have simulated the interaction of a relativistic jet with a turbulent ISM. We have shown that the primary driver of feedback in these cases is the evolution of an energy bubble inflated by the jet. In the initial phases when the jet is confined within the galaxy’s ISM, the energy bubble progresses nearly spherically, strongly interacting with the galaxy’s ISM, creating multiphase outflows. We have previously considered a spherical distribution of gas, which is appropriate for an elliptical system formed by a merger and triggering an AGN (Gatti et al. 2016).

However, in the high-redshift Universe (⁠|$z$| ∼ 2), a significant fraction of the galaxies have gaseous discs (Glazebrook 2013; Wisnioski et al. 2015). If an AGN is triggered in such galaxies, the jet (or nuclear wind) encounters a disc-like morphology of the ISM. Observations of feedback from radio-loud galaxies at high redshifts (Collet et al. 2016; Nesvadba et al. 2017) confirm this scenario. Since a disc-like geometry provides lower column depth along the path of the jet, it is not clear whether jet feedback can have a significant effect on the galaxy’s ISM.

Previously, Sutherland & Bicknell (2007), Gaibler, Khochfar & Krause (2011), Gaibler et al. (2012), Dugan, Gaibler & Silk (2017), and Cielo et al. (2018) carried out simulations of jets and/or winds from a central black hole affecting a gas disc in a galaxy. These results suggest that a jet can indeed affect the gas distribution as well as the star formation kilo-parsecs away from the jet axis. However, these papers do not address the effects of different jet power, ISM density, and jet orientation on the jet–ISM interaction in a systematic manner. Sometimes, essential physics such as relativistic effects, external gravity, or atomic cooling were ignored.

Adhering to the framework followed in our previous papers (Wagner & Bicknell 2011; Wagner et al. 2012; Mukherjee et al. 2016), we present here simulations of a relativistic jet interacting with a multiphase turbulent disc. Our primary motivation is to understand how the jet and its associated energy bubble couple with the disc on extended scales (≳ 1 kpc), the dependence on jet orientation with respect to the disc, the effect on the disc kinematics and turbulence, and the implications for the star formation rate (SFR; positive or negative). In our simulations we use the pluto hydrodynamic code and employ a relativistic Riemann solver, atomic cooling, external gravity, and establish a rotating disc with a turbulent velocity dispersion.

We structure the paper as follows. In Section 2, we describe the numerical set up of the simulation. We discuss the simulation code and the method used to set up the turbulent disc and the relativistic jet. In Section 3, we present the results of the simulations. We discuss the impact of jets launched both perpendicular to the disc (Section 3.1) and inclined towards it (Section 3.2). In Section 4, we discuss the possible impact of the jet feedback on the SFR. Finally, in Section 5, we summarize the results and discuss their implications. We have found qualitative similarities between the results our simulations and the observed kinematic and morphological structure of several galaxies such as 3C 326, 3C 293, NGC 1052, and NGC 3079, which we discuss in that section.

2 SIMULATION SET UP

2.1 Code and numerical schemes

We have performed 3D hydrodynamic simulations of a relativistic jet interacting with a turbulent gas disc before breaking out into the galactic halo. The simulations have been carried out using the relativistic hydrodynamic module of the publicly available pluto code (Mignone et al. 2007). The numerical techniques employed are identical to Mukherjee et al. (2016). We use the HLLC Riemann solver to solve the relativistic hydrodynamic equations (Mignone & Bodo 2005), with a piecewise parabolic reconstruction scheme (PPM method; Colella & Woodward 1984). A third-order Runge–Kutta integration scheme is used to perform the time integration.

The simulations were performed on a Cartesian (X–Y–Z) domain of physical dimensions: 4kpc × 4kpc × 8kpc, with a grid of 672 × 672 × 784 cells. The X and Y grids are uniform with resolution ∼6 pc. The Z grid is uniform over the central region (∼± 1.6kpc, 544 grid points), with the same resolution as the other two dimensions. Beyond ∼± 1.6kpc, a geometrically stretched grid was applied over 120 cells with a stretching ratio of r_s ∼ 1.017.¹ Thus, the central region with the gas disc is well resolved with a uniform resolution of ∼6pc³.

Energy losses due to non-equilibrium atomic cooling ([ρ/μ(T)]²Λ(T), where Λ is the cooling function in ergs⁻¹cm⁻³) are included using the results from the mappings v code (Sutherland & Dopita 2017). We tabulate the cooling function Λ as a function of p/ρ, and, during the simulation, we interpolate the cooling rate at every computational cell. The mean molecular weight, μ(T), also obtained from mappings v, monotonically varies from μ ∼ 1.24 at low temperatures to μ ∼ 0.6 when the gas is fully ionized. For temperatures beyond ∼10⁹ K, the cooling function is extended assuming bremsstrahlung emission, similar to Krause & Alexander (2007). A lower temperature threshold of ∼1000 K was enforced, below which cooling is turned off. Molecular cooling is not considered in this work. In the dense cores of the clouds, atomic cooling is strong enough to quickly drive the gas to the cooling floor. Thus, neglecting molecular cooling does not affect the dynamics of the gas in our simulations.

2.2 Two phase initial density

The initialization of the two-phase ISM closely follows that of Paper I with some differences. An external gravitational potential (ϕ(r)) is set up following the double isothermal prescription of Sutherland & Bicknell (e.g. 2007) and Mukherjee et al. (e.g. 2016) which accounts for contribution from baryonic (stellar) and dark components. For the assumed parameters (listed in Table 1), the total baryonic (stellar) mass contained within ∼10kpc is ∼1.1 × 10¹¹M_⊙. A hot isothermal halo (T ∼ 10⁷ K) in hydrostatic equilibrium is established. The density of the halo gas is determined by the following parameters: a central density n_h0 (R = 0), R being the spherical radius, mean molecular weight μ ∼ 0.6, and atomic mass unit m_a = 1.6605 × 10⁻²⁴gcm⁻³. The density is given by

\begin{eqnarray*} n_\text{h}(R)=n_{\text{h}0}\exp \left[\frac{-\phi ({R})}{k_\text{B} T/(\mu m_a)}\right], \end{eqnarray*}

(1)

where k_B is the Boltzmann constant.

Table 1.

Parameters of the ambient gas and gravitational potential common to all simulations.

Parameters		Value
Baryonic core radius	rB	1 kpc
Baryonic velocity dispersion	σB	250 km s−1
Ratio of DM to Baryonic	λ	5
core radius
Ratio of DM to Baryonic	κ	2
velocity dispersion
Halo temperature	Th	107 K
Halo density at r = 0	nh0	0.5 cm−3
Turbulent velocity dispersion of warm clouds	σt	250 km s−1

Notes. (a) For the assumed gravitational potential, total Baryonic (stellar) mass within a radius of 5 and 10 kpc are |$M_\text{B}(r\lt 5 \mbox{ kpc})= 9.52\times 10^{10}\, \text{M}_\odot$| and |$M_\text{B}(r\lt 10 \mbox{ kpc})= 1.04\times 10^{11}\, \text{M}_\odot$|⁠, respectively. (b) The total mass of dark matter within a radius of 5 and 10 kpc are M_DM(r < 5kpc) = 2.56 × 10¹¹M_⊙ and M_DM(r < 10kpc) = 9.1 × 10¹¹M_⊙, respectively.

Fig. 1 shows the variation of the halo density with radius for the simulations presented here. Beyond the core radius the halo density decreases (n_h ∼ 0.2 at r ∼ 1kpc) rapidly as a power law.

Figure 1.

Variation of the density (n_h(R)) of the hot ambient medium as a function of spherical radius (R).

We embed a turbulent gas disc of radius ∼2kpc into the halo. The density of the turbulent disc is described by a distribution of fractal clouds. The clouds follow a single point lognormal distribution and a two point power-law correlation with a Kolmogorov spectrum (Sutherland & Bicknell 2007; Wagner & Bicknell 2011; Mukherjee et al. 2016). The fractal density is created using our publicly available pyFC routines;² the cloud distribution is then interpolated into the pluto computational domain. The fractal cubes generated are characterized by the following parameters: physical size 1.8kpc × 1.8kpc × 1.8kpc on a grid of dimension 350 × 350 × 350 points, maximum wavelength λ_max = 1.8/6 = 0.3kpc,³ with a mean μ = 1 and dispersion of |$\sqrt{5}$| describing the lognormal distribution. The average size of the largest clouds is, thus, λ_max/2 = 150pc.

The fractal density distribution is then apodized with a mean warm phase density profile given by equation (2) below. The parameters defining the number density of the warm gas n_w(r, |$z$|⁠) as a function of cylindrical coordinates r and |$z$|⁠, are the central density n_w(0, 0), and the total dispersion σ_t. We neglect the contribution of the thermal component (⁠|$kT/(\mu \text{m})$|⁠) to the turbulent dispersion, which for a cold disc at T = 10³ K is only a few kms⁻¹, much smaller than the imposed velocity dispersion. The mean warm density is

\begin{eqnarray*} \frac{n_\text{w}(r,z)}{n_{\text{w}0}}\!\!=\exp \left\lbrace -\frac{1}{\sigma _\text{t}^2} \left[ \phi (r,z) - \epsilon ^2 \phi (r,0) -(1-\epsilon ^2) \phi (0,0) \right] \right\rbrace \nonumber\\ \end{eqnarray*}

(2)

(Sutherland & Bicknell 2007). In the above we assume the azimuthal velocity component to be a fraction ε of the Keplerian velocity:

\begin{eqnarray*} v_{\phi }=\epsilon \left[\frac{r \partial \phi (r,z)}{\partial r}\right]^{1/2}. \end{eqnarray*}

(3)

For this work, we assume ε = 0.93 for all simulations. The fractal cloud is set to be in pressure equilibrium with the ambient gas. Thus, initially the temperature of a cloud cell is T_cloud = n_hT_h/(n_w). For a cloud temperature greater than a critical value (T_cloud > 3 × 10⁴ K), the fractal density is replaced with halo gas, and this places a lower bound on the cloud density. Beyond the critical temperature the clouds are considered to be thermally unstable (Sutherland & Bicknell 2007).

In addition to the rotation, we impose a turbulent velocity on the distribution of warm clouds. Each Cartesian component of the velocity is drawn from an independent random Gaussian distribution with dispersion ∼100kms⁻¹. The turbulent velocity field is also created using the pyFC routine on a uniform grid of size same as the density field. The velocity field is initialized with a Kolmogorov power spectrum (pre-apodization) and λ_max = 0.6kpc, twice that of the density field. The λ_max = 0.6kpc for the velocity field is chosen to be larger than that of the density distribution to ensure that the clouds do not immediately blow up at the start of the simulation due to the imposed turbulent velocity. The turbulent velocity field thus created is then interpolated on to the computational domain, and apodized with the medium. Computational cells not deemed as the warm medium (T > T_cloud), are replaced by halo gas, with the velocity set to zero.

The fractal disc is then evolved without injecting the jet. The clouds interact and settle into a turbulent disc, similar to the results on cloud settling discussed in Paper I. The density after settling is shown in Fig. 2. The settled density thus forms a clumpy turbulent disc with average cloud sizes of ∼100 pc, similar to clumpy star-forming regions observed in high-redshift galaxies (Jones et al. 2010; Livermore et al. 2015). The power spectrum of the settled disc is close to the Kolmogorov slope (see more in Appendix  B). During the settling, the velocity dispersion decreases, typical of decaying turbulence without a driving force (also discussed in Paper I). The effective dispersion of the settled disc which the jet encounters is ∼80–100kms⁻¹, similar to the observed dispersion in high-redshift galaxies (Förster Schreiber et al. 2009; Wisnioski et al. 2015).

Figure 2.

Left: Density (log (n[cm⁻³])) of the turbulent disc in the X–Z plane after brief settling. Right: Density in the X–Y plane. The jet is injected into this settled turbulent disc structure.

2.3 Jet parameters

The jets are injected as a bi-conical outflow with an opening angle of 20°, as shown in the schematic in Fig. 3. The jet radius at the injection zone, is 30 pc. This ensures at least 10 computational cells across the diameter of the jet. The injection occurs from a spherical inner boundary inside the computational domain. The computational cells inside the inner boundary are set to have values characteristic of the halo gas and are not evolved with time. The radius of the spherical inner domain is 60 pc. For simulations A, B, F, G, and H (see Table 2), the jet axis is aligned with the Z-axis. For simulations C, D, and E, the jet axis is inclined from the Z-axis by an angle θ_inc. The simulation parameters are listed in Table 2.

Figure 3.

The geometry of the jet injection through an internal spherical pluto boundary. r_jet = 30 pc is the jet radius. r_sph = 60 pc is the radius of the internal spherical region. The origin of the outflow is the centre of the grid at (0,0,0). The outflow is directed radially in the shape of a cone whose virtual apex is below (0,0,0), with a half-angle of α = 10°.

Table 2.

Simulations.

Simulation label	Power (ergs−1)	n|$w$|0 (cm−3)a	|$\theta _{\rm inc}\, ^b$|	|$\chi \, ^c$|	|$\Gamma \, ^d$|	Gas mass (⁠|$10^9\, \text{M}_{\odot }$|⁠)
A	1045	100	0	10	5	2.05
B	1045	200	0	10	5	5.71
C	1045	200	20	10	5	5.71
D	1045	200	45	10	5	5.71
E	1045	200	70	10	5	5.71
F	1046	100	0	80	10	2.05
G	1046	200	0	80	10	5.71
H	1046	400	0	80	10	20.68

^a The density of the dense gas at r = |$z$| = 0.

^b Inclination of the jet axis with respect to the Z-axis.

^c The ration of rest mass energy to relativistic enthalpy given by equation (5).

^d The Lorentz factor. Γ = 10 corresponds to jet injection velocity of 0.995c and Γ = 5 corresponds to 0.979c.

The following parameters characterize the relativistic jets injected in the simulation domain. The Lorentz factor of the jet is |$\Gamma =1/\sqrt{1-\beta ^2}$|⁠, β = |$v$|/c being the jet velocity in units of the speed of light. We assume an ideal equation of state (Mignone & McKinney 2007), where the specific enthalpy (h) is related to the pressure and density by

\begin{eqnarray*} \rho h=\rho c^2+\frac{\gamma _{\rm ad}}{\gamma _{\rm ad}-1} p \end{eqnarray*}

(4)

The adiabatic index is assumed to be γ_ad = 5/3, which is appropriate for the thermal gas in the turbulent disc into which the energy from the jets is deposited (as discussed in Mukherjee et al. 2016). We define the parameter χ for a jet with adiabatic index γ_ad and jet pressure p_jet as

\begin{eqnarray*} \chi=\frac{\rho c^2}{\rho h -\rho c^2}=\left(\frac{\gamma _{\rm ad}-1}{\gamma _{\rm ad}} \right) \frac{\rho c^2}{p_{\rm jet}} \end{eqnarray*}

(5)

The kinetic power injected by the jet across a surface A_jet is given by

\begin{eqnarray*} P_{\rm j} &=\left(\Gamma ^2 \rho h - \Gamma \rho c^2 \right) c \beta A_{\rm jet} \nonumber \\ &=\frac{\gamma _{\rm ad}}{\gamma _{\rm ad} -1} c p_{\rm jet} \Gamma ^2 \beta A_{\rm jet} \left(1 + \frac{\Gamma -1}{\Gamma } \chi \right). \end{eqnarray*}

(6)

For a given value of jet power (P_j), Lorentz factor (Γ), and χ (as in Table 2), the pressure inside the jet is obtained from equation (6). The density inside the jet can then be computed fromequation (5).

3 RESULTS

We have performed a suite of simulations with different jet powers and orientations with respect to the disc (see Table. 2). In the following section we discuss the results of thesimulations.

3.1 Vertical jets

In the case of the simulations with jets perpendicular to the disc,⁴ although briefly frustrated by the turbulent ISM, the jets breaks out quickly, within a few 100 kyr. This is much shorter than the longer (∼1 Myr) confinement time-scales observed in simulations with spherically distributed gas and similar gas density and jet power (Mukherjee et al. 2016, 2017).

Fig. 4 shows the density on a logarithmic scale (log (n[cm⁻³])) and vertical velocity V_|$z$| normalized to 100kms⁻¹ in the X–Z plane for simulation B. The white and yellow dotted lines in the density image represent contours of Γ ∼ 2 and ∼1.05, respectively, corresponding to speeds of β ∼ 0.9c and ∼0.5c, respectively. The contours denote the region of the flow with bulk relativistic velocities, and trace the evolution of the jet.

Figure 4.

Top: Density (log (n[cm⁻³])) at different times in the X–Z plane for simulation B with P_j = 10⁴⁵ergs⁻¹, n_|$w$|0 = 200cm⁻³ (see Table 1), launched along the Z-axis. The white and magenta lines denote contours of β = 0.5c and β = 0.9c, respectively. Bottom: The vertical component of the velocity V_|$z$| normalized to 100kms⁻¹ for the dense gas (defined here as n > 0.1cm⁻³). The jet launches a central outflow, but the pressure from the energy bubble drives an inflow at the outer edges.

The figures show the evolution of the jets after break out from the disc. The jets create high pressure energy bubbles driving a strong shock through the ambient medium. The jets and the energy bubble expand rapidly into the ambient medium, sweeping up the gas in the forward shock and evacuating cavities above and below the disc. Initially, the swept up gas lies in a thin shell at the forward shock, which at later stages shows signature of Rayleigh–Taylor instabilities at the contact discontinuity. The low densities (∼0.1–0.2cm⁻³) and high temperatures (T ∼ 10⁷ K) of the swept up gas, imply cooling times of ∼1.4 × 10⁸ yr. Thus, the shell of halo gas is not expected to feed back into the galaxy for times greater than 100 Myr.

The gas inside the bubble is composed of a mixture of non-thermal plasma (whose mass fraction is given by the jet tracer), and outflowing thermal gas dredged up from the disc around the central nucleus. The pressure inside the bubble is nearly uniform, as shown in a snap shot in the left-hand panel of Fig. 5. The mean pressure inside the bubble decreases with time (see Fig. 6) as a result of the adiabatic expansion of the bubble as well as cooling losses. The time evolution of the pressure is a power law with a steeper slope at later stages. As discussed later in Appendix  A, the pressure evolution is close to that described by the self-similar solutions of an adiabatically expanding energy bubble. The departure from the ideal relations of the bubble pressure versus time is the result of (1) energy losses from cooling and (2) non-spherical growth of the bubbles.

Figure 5.

Left: Pressure in the X–Z plane normalized to |$p_0=\rho _0 c^2=9.2 \times 10^{-4}\, \mbox{dynes cm}^{-2}$|⁠, where ρ₀ = 0.6165m_a, m_a being the atomic mass unit. The pressure inside the bubble is nearly homogeneous. Along the jet axis multiple recollimation shocks can be seen up to a height of ∼1 kpc. Middle and right: Temperature in the X–Z plane at times corresponding to the middle and right-hand panel of Fig. 4. The white and cyan lines represent contours of β = 0.5c and β = 0.9c, respectively.

Figure 6.

The black lines denote the evolution of the mean pressure (volume averaged) inside the energy bubble, as a function of time. Simulations B (P_j = 10⁴⁵ergs⁻¹, θ = 0°), D (P_j = 10⁴⁵ergs⁻¹, θ = 45°), and G (P_j = 10⁴⁶ergs⁻¹, θ = 0°) are shown. The solid red lines are fits to the portions of the plots which can be well described by a power law. The best-fitting exponent of the time for the different sections are presented above the respective plots.

The interaction of the jet with the gas clumps and the large-scale energy bubble have a significant effect on the gas kinematics as we outline below. These simulations have the following features:

• Vertical outflow from the nuclear region:

  In the central region (r ≲ 500pc) the jets launch an outflow of shocked gas. Clouds accelerated by the interaction with the jet are ejected with speeds ≳ 200kms⁻¹. The outflowing clumps form cometary tails because of ablation by the surrounding flow. The diffuse ablated material is swept away at higher speeds exceeding ≳ 500kms⁻¹ (as also noted in Mukherjee et al. 2016). At late stages the outflowing material is seen to extend beyond 1kpc from the centre, forming diffuse filamentary structures. The total amount of gas ejected from the disc to heights greater than ≳ 1.5 kpc, is between ∼2 × 10⁵ and 4 × 10⁶M_⊙. The mass of the ejected gas increases with greater jet power (as shown in Table 3), as well as inclination which we discuss later. The amount of mass ejected is a small fraction of the total mass in the disc. This is because the cold dense rapidly cooling clumps are more resilient to feedback, as also discussed in earlier works (Nayakshin & Zubovas 2012; Mukherjee et al. 2016). It must be noted that inadequate spatial resolution results in numerical diffusion of the ablated material (Cooper et al. 2008). Hence, the velocities of the ablated material should be considered as an upper limit.⁵

• Outward radial flows in the disc plane:

  As the jet-driven energy bubble expands, its internal pressure drives strong outward radial flows into the plane of the disc. Fig 7 shows the density, temperature and radial velocity in the X–Y plane for simulation G, with a jet kinetic power of ∼10⁴⁶ergs⁻¹. The snapshots correspond to a time (∼0.81 Myr) when the jet has evolved significantly and the energy bubble has spread over the disc. The shocks from the jet-driven outflow engulf nearly all of the disc. The central region (r ∼ 1.5 kpc) shows strong radial outflow with velocities exceeding ∼500kms⁻¹.

  Fig. 8 shows the evolution of the mass-weighted mean radial velocity as a function of radius at different times for simulation G, and Fig. 9 shows the mass-weighted probability distribution function (PDF) of the radial velocity for simulations B, D, and G. At the beginning |$\bar{V}_\text{r}=0 \mbox{ km s}^{-1}$| as there is no net radial flow in the disc. This also results in a velocity PDF (shown in black in Fig. 9) symmetric about 0. At t ∼ 100 kyr, the central region has a radial velocity of ∼150kms⁻¹ up to ≲ 1kpc, which is the approximate extent of the energy bubble inside the disc at that time. At later stages, the radial outflows extend to the entire disc and the velocity PDF develops an extended tail of high radial velocity (≳ 1000kms⁻¹).

  The radially driven outflows evacuate the central region (radius ∼1kpc) of the initial dense volume filling gas. The evacuated region is composed of non-thermal plasma originating from the jet together with remnant filaments of gas clouds, and material ablated from them. A similar result has been discussed in Gaibler et al. (2012), who pointed out that the initial blast wave of the jet-driven bubble, while confined to the disc, causes the lateral evacuation of disc material. In our simulations, the vertical extent of the density at the jet launch sites is small (≲ 200pc). Hence, the jet confinement times are smaller than those in Gaibler et al. (2012). However, even after jet break out, which happens fairly quickly for our simulations (t ≲ 100 kyr), the cavity continues to expand due to the lateral expansion of the jet-driven bubble.

• Non-circular gas motions due to jet-driven outflows:

  The jet-driven outflows result in departure of the gas kinematics from the rotation initially established. In Fig. 10, we plot the radial variation of the anisotropy parameter which we define here as

  \begin{eqnarray*} \beta _{\rm aniso}= \frac{\langle V_\text{r}^2 \rangle + \langle V_z^2 \rangle }{2 \langle V_\phi ^2 \rangle }, \end{eqnarray*}

  (7)

  where |$\langle V_i^2 \rangle$| implies a mass-weighted average of the square of the ith velocity component. A rotating disc has β_aniso ∼ 0, whereas β_aniso ≳ 1 implies significant non-circular motions.⁶ The grey dashed line denotes β_aniso at t = 0 for simulation B. Although for a purely circular velocity field β_aniso is expected to be 0, the initial small non-zero value (β_anisp ≲ 0.5) results from the turbulent velocity field imposed on the rotating disc. The values of β_aniso for other simulations are similar at t = 0. At later times β_aniso significantly increases in the inner kpc, as a result of the outward radial and vertical outflows induced by the jet, as discussed in the earlier sections. Thus, the central region has significant non-circular motions of the gas clouds, which are being directly impacted by the jet.

• Inflows into the disc driven by the energy bubble: The evolving high pressure bubbles spread over the disc, driving shocks and inflows into the disc towards the mid-plane; this is apparent in the velocity plots of Fig. 4. This behaviour contrasts with the outflows launched by the jet from the central nucleus. From the temperature plots in the second and third panels of Fig. 5 and the corresponding velocity maps in Fig. 4, we see that the shocks driven by the high pressure bubble penetrate the disc with inflow speeds of ∼500kms⁻¹. Shocks also spread radially outwards from the centre due to the lateral expansion of the bubble. Consequently, the cold component of the disc is gradually engulfed with shocks driven by the bubble, which raises the temperature of the gas to ≳ 10⁵ K. Such a mechanism is usually expected to negatively affect star formation inside the disc. However, there may also be an enhancement in star formation as a result of the increased pressure compressing the disc (e.g. Bieri et al. 2016), and also from the density enhancement resulting from radiative shocks (Gaibler et al. 2012). We discuss this further in Section 4.

• Turbulence in the disc:

  The jet-induced flows, driven into the disc, shock the initially cold cloud clumps and induce turbulence. To probe the relative strength of compressible (⁠|$\nabla \times \mathbf {v}=0$|⁠) and incompressible motions (⁠|$\nabla \cdot \mathbf {v}=0$|⁠), we present in Fig. 11 the compression ratio defined as (e.g. see Iapichino et al. 2011): |$r_\text{c}=|\nabla \cdot \mathbf {v}|^2 /(|\nabla \cdot \mathbf {v}|^2 + |\nabla \times \mathbf {v}|^2)$|⁠. Similar to Bourne & Sijacki (2017), we find that the lobe is dominated by incompressible turbulence with low compression ratio (∼0.1–0.2). The outer bow shock has a higher compression ratio ∼0.8. Within the disc, there are certain regions with higher r_c where the shocks from the energy bubble impinge on the disc, or due to cloud-cloud collisions. However, overall, the disc is dominated by incompressible turbulence as well. The mean value of the compression ratio is r_c ∼ 0.17 in the jet lobes and r_c ∼ 0.2 in the disc.

  Compressible motions result from shocks generated by the expanding bubble and sound waves and weak shocks due to the turbulent motions in the cocoon. Incompressible flows result from turbulence and shear flows generated by the jet–ISM interaction. While compressible modes can decay, being carried by sound waves, incompressible modes can persist for a longer time (Reynolds, Balbus & Schekochihin 2015) contributing to the turbulent energy budget in the gas.

  Fig. 12 shows the increase in dispersion in the disc due to jet induced turbulence. We compute the mass-weighted mean velocity dispersion |$\tilde{\sigma }^2=\frac{1}{3}\left(\tilde{\sigma }_\text{r}^2+\tilde{\sigma }_\phi ^2+\tilde{\sigma }_z^2\right)$|⁠, where the dispersion for each component is computed as

  \begin{eqnarray*} \tilde{V}_i&=\frac{1}{M}\Sigma \rho V_i \nonumber \\ \tilde{\sigma }_i^2 &=\frac{1}{M} \Sigma \rho \left(V_i-\tilde{V}_i\right)^2 \end{eqnarray*}

  (8)

  where the index i corresponds to the (x, y, |$z$|⁠) components. The dispersion has been computed in annular rings as a function of cylindrical radius. This is strictly appropriate only for azimuthally symmetric discs. For jets pointed in the Z direction one expects that the mean jet-induced motions will on average be azimuthally independent. Nevertheless, for inclined jet as well, an azimuthal average conveys a good sense of the radial dependence of the jet-induced velocity dispersion.

  We compute the dispersions separately for gas with temperatures above and below T = 10⁴ K, to distinguish the shocked gas from the undisturbed gas, which is initially at a temperature corresponding to the cooling floor. This is done to investigate the effect on the kinematics of the different gas phases, as the jet-driven shocks create a multiphase environment (as discussed in Paper I). The top two panels of Fig. 12 show the radial distribution of the mean dispersion at different times for simulations G (P_j = 10⁴⁶ergs⁻¹) and B (P_j = 10⁴⁵ergs⁻¹), respectively. The dispersion of the shocked gas is represented in red, the cold undisturbed component in blue and the total dispersion including all gas phases is in green. We see that the dispersion in the shocked gas is increased by several times its initial value, indicating the strong turbulence induced by the shocks driven into the disc. For simulation B the dispersion is raised to ∼400kms⁻¹ throughout the disc. For simulation G with a higher jet power, the value is higher at ∼600kms⁻¹.

  The dispersion in the disc varies with time, depending on the evolution of the energy bubble. At t = 0.1 Myr, only the central region (≲ 1 kpc) has an enhanced dispersion, as it demarcates the location of the energy bubble. At later times the dispersion is the enhanced throughout the disc. The dispersion in the central region decreases with time. This occurs because (a) the turbulence naturally decays with time, (b) the mass from the central region is evacuated by the jet and the computed dispersion being mass-weighted diminishes. Overall, at late stages, after the jet-driven shocks have penetrated the disc, the dispersion increases on average. The cold undisturbed component remains steady at the initial value of ∼80kms⁻¹, with a slight decrease due to the inherent decay of turbulence, as discussed earlier in Section 2.2.

  The last two panels of Fig. 12 presents the radial distribution of the velocity dispersion (mass-weighted) for various simulations. The dispersion is computed at a time when the jets have significantly evolved, and the energy bubble has spread over the disc within the simulation domain. The simulations with jet power P_j = 10⁴⁶ergs⁻¹ (F, G, H) are presented in purple in the third panel, and those with P_j = 10⁴⁵ergs⁻¹ (A, B) are in brown in the last panel. The dispersions have been computed for all gas masses, irrespective of temperature or phase.

  We firstly notice that with a factor of 10 increase in jet power, the mean dispersion increases by a factor of 2–3 for the runs with similar mean density. Overall, the dispersion in the disc increases to several hundred kilometres per second. Secondly, the runs with lower mean density (n_|$w$|0 ∼ 100cm⁻³) have a higher dispersion for the same jet power. Gas clumps with lower density are more susceptible to fragmentation into smaller clumps as a result of hydrodynamic instabilities, which then further interact with the jet-driven flows, raising the overall dispersion of the disc. The dispersion in simulation H with n_|$w$|0 ∼ 400cm⁻³ is not very different from that in simulation G with n_|$w$|0 ∼ 200cm⁻³.

  The above is well demonstrated in Fig. 13, which shows the mass-weighted probability distribution function (PDF) of the density of simulations F, G, and H. A mass-weighted density PDF, P_m(ρ), highlights the regions with higher mass i.e. the gas disc in this case. Thus, it is a better tool to describe the nature of the density distribution in the gas disc than a volume-weighted PDF.⁷ The dashed lines represent the PDF at the start of the individual simulations. The solid lines show the density PDF at the times corresponding to those in Fig. 12.

  Beyond a critical density of ∼100–300cm⁻³, all the simulations show an increase in the PDF from its starting value. This results from the compression and enhancement of densities in the post-shock regions. The shocks driven into the clouds quickly become radiative resulting in a shell of high density surrounding the clouds, making them less prone to ablation. Similar results have also been discussed in Mukherjee et al. (2017) and Sutherland & Bicknell (2007). However, simulation F with n_|$w$|0 = 100cm⁻³, is more extended at lower densities (n ≲ 50cm⁻³), compared to that of G (n_|$w$|0 = 200cm⁻³) and H (n_|$w$|0 = 400cm⁻³). This implies that the lower density clouds in simulation H are more easily dispersed, compared to the denser ISM where clouds remain more intact. The enhanced cloud ablation also raises the velocity dispersion in the turbulent gas. The present simulations do not include self-gravity. The gas more prone to ablation below the critical density of ∼100–300cm⁻³ is not expected to be affected by its self-gravity, as the free-fall time-scales are much larger than the dynamical scales at such densities. However, self-gravity may enhance the PDF at the high-density end. We plan to explore this in a future work.

Figure 7.

Top left: The density in X–Z plane for simulation G with P_j = 10⁴⁶ergs⁻¹, n_|$w$|0 = 200cm⁻³. The image is zoomed on to the central region, focusing on the dense disc. The white line denotes contours of constant β = 0.9c, representing the relativistic component of the jet. Top right: The density in X–Y plane for simulation G showing the evacuation of the cavity. Bottom left: Temperature in X–Y plane for simulation G, showing the radial spread of the shocks driven by the energy bubble. Bottom right: Cylindrical radial velocity in the X–Y plane, normalized to 100kms⁻¹.

Figure 8.

Evolution of the mean cylindrical radial velocity (mass-weighted) as a function of radius, at different times for simulation G.

Figure 9.

The mass-weighted PDF of the radial velocity for simulations B, D, and G. The black solid line shows the PDF at t = 0. All the simulations start from an identical ISM. The PDF shows an extended tail at high velocities due to the radial outflows driven by the jet.

Figure 10.

The ratio β_aniso as defined in equation (7). The central regions (≲ 500pc) show significant non-circular motions (β_aniso ≳ 1). The plots have been generated at times similar to that in Fig. 12.

Figure 11.

Left: Vertical component of the vorticity (⁠|$\omega=\nabla \times \mathbf {v}$|⁠) at time t = 0.48 Myr in the X–Z plane. Right: Compression ratio |$r_\text{c}=|\nabla \cdot \mathbf {v}|^2 /(|\nabla \cdot \mathbf {v}|^2 + |\nabla \times \mathbf {v}|^2)$|⁠.

Figure 12.

Top: Mass-weighted velocity dispersion in the disc as a function of radius at different times for simulation G (P_j = 10⁴⁶ergs⁻¹). The red lines show the velocity dispersion for gas with T > 10⁴ K representing the shocked gas. The blue lines are for gas with T < 10⁴ K, representing the cold undisturbed component. The green lines show the dispersion for all mass, without phase separation. Second: Same as above, for simulation B (P_j = 10⁴⁵ergs⁻¹). Third and fourth: Total velocity dispersion for each simulation, following the spread of the energy bubble over the disc. The third panel shows simulations (F, G, H) with P_j = 10⁴⁶ergs⁻¹ but different mean densities (n_|$w$|0 = 100, 200, 400cm⁻³). The fourth panel is for P_j = 10⁴⁵ergs⁻¹ (A, B) with mean densities n_|$w$|0 = 100, 200cm⁻³, respectively.

Figure 13.

Mass-weighted density probability density function for simulations F, G, and H. These simulations have the same kinetic power (P_j = 10⁴⁶ergs⁻¹) but differ in density (n_|$w$|0: 100, 200, 400cm⁻³ for F, G, H, respectively). The PDF at t = 0 is presented as a dashed line. All simulations show an increase in the PDF at the higher densities ≳ 100cm⁻³ from its initial value. The lower density ISM in simulation F has a higher PDF at lower densities, due to greater cloud ablation.

Table 3.

Mass ejected beyond 1.5 kpc.

Simulation	Mass ejecteda
A	|$3.71\times 10^3\, \text{M}_\odot$|
B	|$2.17\times 10^5\, \text{M}_\odot$|
C	|$6.92\times 10^4\, \text{M}_\odot$|
D	|$4.11\times 10^6\, \text{M}_\odot$|
E	|$3.55\times 10^6\, \text{M}_\odot$|
F	|$2.23\times 10^6\, \text{M}_\odot$|
G	|$4.33\times 10^6\, \text{M}_\odot$|
H	|$1.98\times 10^6\, \text{M}_\odot$|

^a The amount of mass in |$z$| > 1.5 kpc, with cloud tracer >0.98. The values presented are at times corresponding to the end of the simulation, with the initial value at t = 0 subtracted.

3.2 Inclined jets

Recent papers (Gallimore et al. 2006; Combes 2017) have proposed that there is a lack of preferred orientation of jets with respect to the plane of the disc (or the galaxy’s major axis). There are a number of galaxies in which the jets are inclined to the disc normal (e.g. NGC 3079, Cecil et al. 2001; NGC 1052, Dopita et al. 2015), including galaxies in which the jet is in the plane of the disc (e.g. IC 5063, Morganti et al. 2015). In order to study such interactions we have performed simulations (simulations C, D, and E), for which the jet axis is inclined towards the disc. Fig. 14 shows the evolution of the density in the X–Z plane for simulation D where the jet is inclined at 45° to the plane of the disc. Zubovas et al. (2013) and Zubovas & Bourne (2017), The evolution of jets is markedly different in this case, compared to the cases of the perpendicular jets discussed earlier. Since the jets encounter a larger column depth of clouds along their path, they interact more strongly with the dense ISM. The jets do not immediately clear a channel along the inclined axis of launch. Rather, the dense clouds deflect the jets, which then decelerate and launch a sub-relativistic outflow along the minor axis following the path of least resistance (see Fig. 14).

Figure 14.

Left 3 panels: Density (log (n[cm⁻³])) in the X–Z plane at different times for simulation D with P_j = 10⁴⁵ergs⁻¹, θ_inc = 45° (Table 1). The contours denote regions with relativistic gas motions with white giving contours of β = 0.4c and β = 0.9c in magenta. The inclined jets decelerate, launching a spherically evolving bubble along the minor axis. Right: Vertical component of the velocity V_|$z$| in units of 100kms⁻¹ for gas with density n > 0.1cm⁻³ at t ∼ 2.14 Myr, the last density panel in the previous set. The inclined jet couples more with the ISM, resulting in greater outflow from the disc.

Fig. 15 shows the time evolution of the gas speed (β = |$v$|/c). The first panel shows the evolution of the jet for simulation B with the launch axis perpendicular to the disc. It can be clearly seen that the main jet stream remains relativistic with β ∼ 0.9c up to the hotspots, where it becomes turbulent and decelerates. The second and third panels show the evolution of the velocity for simulation D with the jet inclined at 45°. The jets remain relativistic in the central region (≲ 200pc), where they evacuate a cavity. Beyond the central cavity, as the jets interact with the clouds, they quickly decelerate and launch a sub-relativistic outflow (⁠|$v$| ∼ 0.01–0.1c) of non-thermal plasma.

Figure 15.

Left: Evolution of the velocity magnitude, β = |$v$|/c, for simulation B with θ_inc = 0°. Middle and right: Evolution of β at two different times for simulation D with θ_inc = 45°. The times correspond to the middle and right-hand panels of Fig. 14. The contours (white for β = 0.4c and cyan for β = 0.9c) denote regions with relativistic gas motions. The figures show sub-relativistic wider angled outflows in the inclined simulations, as compared to the vertical relativistic outflows for jets launched perpendicular to the disc.

The wind creates a spherical bubble which clears the ambient gas as it rises, similar to the vertical jets. However, the evolution of the bubble is more spherical in morphology, driven by the wide-angled wind, when compared to the ellipsoidal bubbles created by the vertical jets. The evolution of the bubble is also slower in the case of the inclined jets (as shown later in Fig. 16). At later stages, the northern jet approximately proceeds along the axis of launch to some extent. The southern jet however, is deflected due to strong interactions with a cloud and is bent towards the minor axis, along which it continues into the previously evacuated cavity.

Figure 16.

Left and middle: Density in the X–Z plane for simulations C and E with θ_inc = 20° and 70°, respectively. Right: A three panel plot showing the time evolution of the mass-weighted mean radial (spherical) velocity, the kinetic energy of the dense gas (n > 1cm⁻³), and volume of the energy bubble in units of kpc³, for simulations with different angles of inclination. Computational cells with a jet tracer value higher than 1e − 8 are considered to form the bubble.

The strong interaction of the jet with the dense gas in the disc results in dense gas being vertically lifted off the disc, as well as in the radially driven outflows described earlier. The rightmost panel in Fig. 14 shows the vertical velocity in units of 100kms⁻¹ at t ∼ 2.14 Myr. Clouds outflowing at ∼1000kms⁻¹ can be seen being driven to heights of ∼1.5–2 kpc, forming cometary tails of stripped gas. With the jets inclined, the interacting clouds experience an outward vertical component of the imparted momentum resulting in an outflow away from the disc. This is in contrast to the results of simulation B, where, despite the entrainment of clouds near the core of the jet, the dominant effect is the compression of the disc from the energy bubbles, resulting in inflows at the outer edges (≳ 500pc) of the disc.

The strength of the interaction increases with the angle of inclination as the jet drills through a larger gas column. Fig. 16 shows the density slices for simulations D and E, with jet inclinations θ_inc = 20° and 70°, respectively. Simulation C proceeds similarly to that of B, the inclination angle being ∼20°. In simulation E, the jet is frustrated inside the disc, similar to that of simulation C (with θ_inc = 45°). The vertical outflow of clouds is more widespread as the jet ploughs through the disc. The third plot shows the evolution of the mean spherical radial velocity (mass-weighted), the kinetic energy in the dense gas (for n ≳ 0.1cm⁻³) normalized to its initial value and the volume of the energy bubble for simulations with different angles of inclinations (simulations B, C, D, and E).

From the evolution of the mean radial velocity (mass-weighted) and kinetic energy, it is clear that jets aligned more closely to the disc interact more strongly. For simulations B and C with smaller inclination angles (measured from the disc-normal), the mean radial velocity increases until about t ∼ 0.5 Myr, decreasing thereafter as the jet–ISM interaction decreases in intensity. For simulations D and E, the mean radial velocity increases for a longer time (∼1 Myr), as the jets are frustrated inside the disc. The kinetic energy similarly shows an initial increase followed by a decrease. For simulations E and D, with jets more inclined to the disc, both the mean radial velocity and kinetic energy are significantly higher than in simulations B and C. Similarly, the evolution of the volume of the energy bubble presented in the third panel shows that simulations D and E closely follow each other, and evolve at a much slower rate than that of B and C. This is because of the lower velocity of the wide-angled sub-relativistic winds, which power the bubbles formed from the inclined jets.

The jet-driven outflows do not unbind a significant fraction of the ISM (by mass), as also discussed in Mukherjee et al. (2016). The mean radial velocity, though higher for simulations D and E, is much less than the velocity dispersion of the galaxy.⁸

Instead, the jet–ISM interactions induces local turbulence within the disc, as evidenced by the higher kinetic energy of simulation E, as compared to others. This also raises the velocity dispersion of the gas disc. In Fig. 17, we present the radial variation of the velocity dispersion for simulations B, C, D, and E. In the central regions, the resultant velocity dispersion of simulation E was found to be higher than that of simulation B by a factor of ∼3. To summarize, jets more closely inclined with the disc, couple more strongly with the ISM, over a longer time, inducing more localturbulence.

Figure 17.

The mass-weighted velocity dispersion as a function of radius for simulations B, C, D, and E. The times have been chosen such that the jet-driven energy bubble has spread over the disc.

3.3 Qualitative similarities with some observed galaxies

3.3.1 Shocked disc in NGC 1052

NGC 1052 is a nearby (⁠|$z$| = 0.0045) LINER⁹ galaxy with a prominent radio jet (Claussen et al. 1998; Wrobel 1984). Observational studies (Kadler et al. 2004a,b; Sugai et al. 2005; Dopita et al. 2015) have shown signatures of an inclined radio jet strongly interacting with its environment. We list below some of the key inferences drawn from these observational signatures.

• The inner sub-arcsecond region shows strong O iii emission (Sugai et al. 2005), characterized by two ridges extending east-northeast and west-south-southwest. The spectrum shows signatures of two components emitting in O iii, a narrow low-velocity component that is spatially extended and a less extended component with higher velocity, with values reaching ∼1000kms⁻¹ (see Section 3.2 of Sugai et al. 2005).

• Mapping the O iii in different emission channels reveals spatially distinct components with a central core where the outflow originates and two others which appear as working surfaces of a jet-driven shock (Fig. 4 of Sugai et al. 2005).

• The outflow axis (which follows the jet) is inclined to the rotation axis of the disc by ∼50° (Sugai et al. 2005; Dopita et al. 2015).

• Above and below the disc there exists outflowing |$\text{H}\, \alpha$| bubbles with diameters of ∼1.5 kpc, extending along the minor axis of the galaxy (Dopita et al. 2015). The bubbles have an internal dispersion of ∼150kms⁻¹.

• The gas in the disc in strongly disturbed with a high-velocity dispersion ∼300kms⁻¹.

• Kadler et al. (2004a) have found spatially extended X-ray emission which is thermal in origin (with T ∼ 0.4–0.5 keV). The extent of the X-ray emitting region is co-incident with radio emission at 1.5 GHz (Kadler et al. 2004a,b) and the turbulent gas disc, pointing to the presence of strongly shocked gas (Dopita et al. 2015).

Dopita et al. (2015) have put forward a model of a jet inclined to a gas disc to explain the above (see schematic in fig. 7 of Dopita et al. 2015). Our simulations of inclined jets, specifically simulation D with the jet inclined at 45°, qualitatively address several of these results and confirms the suggestion put forth in Dopita et al. (2015). As the jet impinges on the disc, it shocks and ablates the clouds. The resulting morphology is that of an arc near the head of the jets (see Fig. 18). Such an arc of shocked clouds, can indeed appear as a ridge of O iii emission, similar to the bow-shock shaped emission profile observed in the high-velocity channel in fig. 4 of Sugai et al. (2005).

Figure 18.

Temperature in X–Z plane for simulation D (P_j = 10⁴⁵ergs⁻¹, n_|$w$|0 = 200cm⁻³, θ_inc = 45°). The jet head is surrounded by an arc of dense clouds with shocked surface layers. Shocks from the jet-driven energy bubble spread through the disc. The bubble created on top of the disc is composed of high-temperature non-thermal plasma, as well as shocked thermal gas from ablated clouds. The contours in white (β = 0.4c) and magenta (β = 0.9c) denote regions with relativistic gas motions.

As the jet decelerates, it launches a slower moving bubble of hot gas that rises along the minor axis (see discussion in Section 3.2 and Fig. 15). Such a low-density hot bubble is similar to the observed |$\text{H}\, \alpha$| bubbles reported in Dopita et al. (2015). The two velocity components in the O iii emission may be explained by (a) shocks permeating the surfaces of slower moving dense clouds before they are accelerated and (b) emission from the fragments ablated from the clouds, which are swept outwards at speeds of ∼1000kms⁻¹, as shown in Fig. 14.

Although the jet head is confined within the inner few 100 pc (see Fig. 18), the energy bubble created by the jet proceeds farther. Thus, the region of shocked gas extends a few kpc beyond the collimated part of the jet. This implies that the observed high-velocity dispersion in the gas disc in NGC 1052 is driven by turbulence induced by the jet, as proposed by Dopita et al. (2015), and also demonstrated in Fig. 17.

Thus, although simulation D is not a precise model of NGC 1052, it shows remarkable qualitative similarities with the observed kinematics and morphology of the shocked gas in the galaxy. This strengthens the argument that wide-scale shocks driven by the jet and the expanding energy bubble may be dominant contributors to the observed emission from NGC 1052 (as suggested by spectral studies in Dopita et al. 2015) and to the kinematics over scales of a few kpc away from the nucleus.

3.3.2 The superbubble in NGC 3079

NGC 3079 is another galaxy which shows kpc scale bubbles arising from a nuclear outflow (Duric & Seaquist 1988; Irwin & Seaquist 1988; Veilleux et al. 1994; Cecil et al. 2001). One striking feature of this galaxy is the presence of ionized filaments extending to about a kpc from the disc, tracing the contours of the bubble (Cecil et al. 2001). Cecil et al. (2001) model the bubbles as arising from the wide-angled mechanically driven nuclear outflow. The source of the outflow has been debated in the literature. The central nucleus does indeed harbour star-forming molecular gas (Lawrence et al. 1985), pointing to the possibility of a starburst-driven wind (Veilleux et al. 1994; Sofue & Vogler 2001; Shafi et al. 2015). However, this is disputed in Hawarden et al. (1995), who favour an AGN¹⁰ driven wind.

Nevertheless, the galaxy harbours a jet whose orientation is misaligned with that of the large-scale outflow (bubbles) by ∼65° (Irwin & Seaquist 1988; Trotter et al. 1998). Irwin & Seaquist (1988) have proposed that such a jet itself can inflate a bubble to form the observed structures. This is similar to our results for inclined jets, where the jets remain relativistic in the central cavity, beyond which they decelerate and launch a slower outflow along the minor axis. As shown in Figs 14 and 16 the inclined jets launch a vertical outflow that creates a superbubble with a morphology usually expected from a nuclear starburst-driven bubble (Veilleux, Cecil & Bland-Hawthorn 2005; Cooper et al. 2008; Vijayan et al. 2018).

The column density plots in Fig. 19 show filamentary structures protruding from the disc, extending to heights of ∼1 kpc from the dense central disc. This is similar to the simulations of filamentary outflows from a purely starburst-driven wind (Cooper et al. 2008). We note that our results do not model the outflows in NGC 3079 precisely. However, our results do qualitatively demonstrate that, as an inclined jet interacts with the disc, it can also drive a wide-angled outflow with gas clouds/filaments being lifted from the disc. Thus, jets inclined to the disc can significantly affect the ISM over a radius of a few kpc.

Figure 19.

Column density along the Y-axis for simulation D (P_j = 10⁴⁵ergs⁻¹, n_|$w$|0 = 200cm⁻³, θ_inc = 45°). The top panel shows the disc at beginning before the jet is injected. The bottom panel is at t ∼ 2.14 Myr when the energy bubble has developed (see Fig. 14).

Similar results have also been explored recently in other simulations such as those (Dugan et al. 2017; Cielo et al. 2018). Jets inclined to the disc may in fact be more common, as the orientation of the jet with respect to the large-scale disc may not be related, and in fact suggested to be random in Gallimore et al. (2006); Combes (2017). Such interactions have also been observed in other systems such as NGC 4258 (Ogle, Lanz & Appleton 2014) and 3C 293 (Mahony et al. 2016a). Another system with a kpc scale bubble, possibly inflated by a jet, is NGC 6764 (Hota & Saikia 2006; Croston et al. 2008). Wide-angled radio bubbles extend along the minor axis of galaxy, similar to the morphology of bubbles insimulation D.

3.3.3 Enhanced turbulent dispersion in gas discs

An increase in disc turbulence as a result of the wide-spread effect of jet feedback has been observed in some sources where direct interactions of jets and a turbulent ISM has been detected. Such results are similar to what we find from our simulations, as we outline below.

• 3C 326: In 3C 326, the emission line kinematics of the multiphase ISM show very high local dispersion, with an FWHM of ∼600kms⁻¹ in the ionized gas (Nesvadba et al. 2010) and warm H2 (Nesvadba et al. 2011), and ∼350kms⁻¹ in the CO (1–0). The H2 (1–0) emission lines show complex line features requiring multiple Gaussian components to fit the profile. The authors conclude in these papers that the emission is primarily shock excited. Such broad complex line profiles indicate that the jets significantly couple with the entire disc, increasing its turbulence. This is analogous to the increased velocity dispersion in the shocked gas (∼300–600kms⁻¹) presented in Fig. 12. The shocked gas in our simulations is hot and ionized (with T ∼ 10⁵ K). Thus, our results are in agreement with the dispersion observed in the ionized components in 3C 326. Since our simulations lack molecular chemistry, we are unable to follow the cooling of the shocked gas to form molecules. However it is likely that the shocked components with high-velocity dispersion will show enhanced gas kinematics in the molecular phase when allowed to cool, as we also argue in Mukherjee et al. (2018).

• 3C 293: Similarly, 3C 293 hosts a powerful radio galaxy where fast outflows have been detected in ionized (Emonts et al. 2005; Mahony et al. 2016a) and neutral (Morganti et al. 2003; Mahony et al. 2013) gas, along with emission of thermal X-rays from shocked gas close to the jet (Lanz et al. 2015). Mahony et al. (2016b) find that the ionized gas exhibits a large width in the emission lines (FWHM ∼200–600kms⁻¹) along slits positioned both along the jet axis and also along slits displaced from the jet axis. Such enhanced kinematics away from the jet axis, is analogous to our results that a jet-driven energy bubble can have a wide-spread effect on the disc, at a few kpc away from the central axis.

4 IMPACT ON STAR FORMATION RATE

4.1 Positive feedback

Relativistic jets can both promote star formation via density enhancements from shocks (positive feedback) or quench it by mass expulsion or by creating turbulence. Shocks from jets increase the density in the radiative post-shock gas, facilitating collapse. Our current numerical set up lacks an explicit mechanism to follow star formation within the gas discs. However, one can approximately estimate the SFR within a computational volume as

\begin{eqnarray*} \mbox{SFR} =\epsilon _{\rm SFR} (\rho d^3x)/t_{\rm ff}. \end{eqnarray*}

(9)

Here, t_ff = (3π/(32Gρ))^1/2 is the free fall time. ε_SFR is an efficiency factor that defines the fraction of gas within a volume that is converted to stars. A similar definition of SFR has been widely used to estimate the SFR in computational simulations of galaxy evolution (Springel & Hernquist 2003; Dubois & Teyssier 2008) and also in recent studies of positive feedback (such as Gaibler et al. 2012; Bieri et al. 2015; Bieri et al. 2016; Dugan et al. 2017).

For the results presented here we compute the SFR for computational cells with density above a threshold n > 100cm⁻³ (similar to Gabor & Bournaud 2014), to select the dense regions with free-fall time-scales of a few million years, which will be prone to star formation. Such densities are typical of the mean density in giant molecular clouds (McKee 1999). The value of ε_SFR depends on the local physical conditions of the gas cloud such as internal velocity dispersion, magnetic fields, and gas dispersal from stellar outflows. Such processes are sub-grid for our simulations since we do not resolve the relevant scales. Typically a value of ε_SFR ∼ 0.01–0.1, used in numerical simulations of star formation in molecular clouds (Federrath & Klessen 2012, 2013), matches well with observations. Similar values have also been estimated from observations of SFRs in high-density star-forming regions in the Milky Way Galaxy (Krumholz & Tan 2007).

For our work, we find that assuming a value of ε_SFR = 0.015 yields beginning SFRs which match well with the Kennicutt–Schmidt relation (Kennicutt 1998). This is shown in Fig. 20 where we present the SFR densities of simulations A, G, and H. Simulations B, D, G have similar initial configurations (with n_|$w$|0 = 200cm⁻³), and hence similar SFR. Similarly simulations A and F with n_|$w$|0 = 100cm⁻³ share identical starting conditions in the gas disc. The initial SFRs at the start of different simulations, discussed later in this section, is presented in Table 4.

Figure 20.

The mean SFR surface density (total SFR/area of disc) compared with mean gas surface density (total mass/area of disc). The dotted line shows the Kennicutt–Schmidt relation (Kennicutt 1998). The initial SFRs (t = 0) for simulations A, G, and H are shown in partially filled symbols. The filled symbols for simulations G and H show the final SFRs at the end of the simulation. Other simulations have initial SFR similar to those of A, G, or H, being started from a similar initial set up, as discussed in the text.

Table 4.

Star formation rates at t = 0.

Simulation label	SFR (t = 0) (⁠|$\text{M}_\odot \, \mbox{yr}^{-1})^a$|
A	5.02
B	17.44
D	17.44
E	17.44
F	5.02
G	16.99
H	86.02

^a Assuming ε_SFR = 0.015 in equation (9). The values scale linearly with ε_SFR.

The above definition of SFR used in this work (and others such as Gaibler et al. 2012; Bieri et al. 2015; Bieri et al. 2016; Dugan et al. 2017) does not depend on the temperature of the gas. An accurate estimate of star formation, must however, account for local Jeans length, ensure a converging flow and that the gravitational potential energy is higher than the local turbulent kinetic energy (Federrath et al. 2010b). However, although the absolute magnitudes of SFRs may change on employing a more robust prescription for star formation, the simple star formation law assumed here, allows for ready comparison of the relative changes in the SFR between the different simulations, with varying mean densities, jet power and orientation, as presented in this work. For the assumed definition, the resultant SFR matches well with the Kennicutt–Schmidt relation (Kennicutt 1998). This ensures that the initial SFR, based on which the relative strength of positive feedback estimates are made, are realistic.

Following the above relation, we compute the evolution of the SFR in the gas discs in our simulations. In Fig. 21, we present the SFR of different simulations normalized to their initial values (see Table 4). We present the relative changes in the SFRs to better highlight the effect of the jet on the SFR in each simulation. We see that the evolution of the star formation strongly depends on the jet power, ISM density, and jet orientation.

Figure 21.

Evolution of the SFR (defined in equation 9) for different simulations. The values of the SFR are normalized to their initial values (presented in Table 4). Top: Simulations F, G, and H with P_j = 10⁴⁶ergs⁻¹ in black. In grey we present the evolution of the SFR without the jets for discs, with densities n_|$w$|0 = 100, 200, 400cm⁻³, respectively. Bottom: SFR for simulations A, B, D, and E. Simulations G and H show a substantial increase in SFR. Simulations A and F with n_|$w$|0 = 100cm⁻³ decline after an initial increase. Simulations D and E with jets inclined to the disc (θ_inc = 45°, 70°, respectively) show an increasing trend in the star formation.

Simulation B with jet power ∼10⁴⁵ergs⁻¹ shows only a modest increase in the overall SFR from its initial value. Simulation A with n_|$w$|0 ∼ 100cm⁻³ shows an initial decline before catching up to its original value and declining again. Simulations D and E with the jet inclined to the disc show a gradual increase by ∼40 per cent at late times. Simulation E with a higher inclination angle (θ_inc = 75°) shows a higher increase in SFR. Jets inclined to the disc couple directly with a larger column of gas, resulting in higher SFR estimates.

Jets with power 10⁴⁶ergs⁻¹ show a much sharper increase in the SFR. Simulation H with a higher density (n_|$w$|0 ∼ 400cm⁻³) shows an increase by a factor of 1.8 during the run time of the simulation, with an increasing trend. Simulation G (with n_|$w$|0 = 200cm⁻³) rises initially and flattens at ∼1.6. However, simulation F with n_|$w$|0 = 100cm⁻³ shows a decline in SFR after an initial increase. This is likely because of shearing of clouds as the jets engulf the disc.

Thus, it is evident that higher jet powers can, in fact, potentially enhance star formation as they drive stronger shocks into the ISM, resulting in higher density compression. In the top panel, we show the evolution of the SFR for the discs with the jet turned off (presented in grey). Without the jet, the SFR decreases with time. Thus, any enhancement in SFR, is purely due to the impact of the jets. These findings are consistent with the recent results on positive feedback such as Gaibler et al. (2012), Dugan et al. (2017), and Bieri et al. (2016). However, discs with lower ISM densities show a decline in star formation after an initial increase due to destruction of clouds by the jet-driven flows.

Although the jets can affect the disc as a whole, jet induced star formation is likely to occur closest to the jet axis, where the pressure from the jet is highest, leading to larger compression of shocked gas. To demonstrate this we present in Fig. 22 the SFR surface density (SFRD) in the X–Y plane, i.e. SFR per unit area, defined as: ∫(ρ/t_ff)d|$z$|⁠. The SFRD presented in Fig. 22 is normalized to a mean star formation rate density, defined as the total SFR at t = 0 divided by the area of the disc (=4π kpc², for a disc of radius 2 kpc). Thus, Fig. 22 shows the relative change in star formation at different locations in the disc to highlight the effect of jet activity on the SFR.

Figure 22.

Top: Star formation rate density (SFRD) for simulation H (P_j = 10⁴⁶ergs⁻¹, θ_inc = 0°), defined as SFR (equation 9) per unit area. The SFRD is normalized to the mean star formation rate density (SFRD0), defined as the SFR(t=0)/(π(2[kpc])²), where the radius of the disc is ∼2 kpc. Bottom: Same as above but for simulation D (P_j = 10⁴⁵ergs⁻¹, θ_inc = 45°) inclined at 45° to the disc. The figures show the relative change in SFRD under the influence of the jet. Simulation H shows an enhancement in SFRD in a circular ring around the central cavity. Simulation D shows patchy increase in SFRD with an elongated quenched region (south-west) tracking the channel carved by the jet.

For simulation H, we see a strong enhancement in the SFR in the regions surrounding the central cavity (≲ 1kpc). This would imply that a star-forming ring may be created around the jet axis, with a young stellar population. Similar results have also been presented in Gaibler et al. (2012). Although the jet has progressed considerably by this time, and the jet-driven bubble has spread over the disc, we do not see appreciable enhancement of SFR (or density) at the outer edges.

For simulation D, the central cavity is not symmetric, as the jet is inclined to the disc and drills through the sides before eventually breaking out. Such a disruption of the disc is clearly visible as an extended region of quenched star formation (cavity). Patches of enhanced star formation are found at the outer edges. Thus, our results show that a relativistic jet can both quench star formation in the central regions due to destruction of clouds, as well as enhance SFRs in the immediate surroundings due to compression from the jet-driven shocks.

Stars formed as a result of positive feedback will have motions which are significantly perturbed from the regular stellar kinematics (circular rotation in the present case). Similar results have been presented earlier by Dugan et al. (2014), Zubovas et al. (2013), and Zubovas & Bourne (2017), where stars formed as a result of positive feedback were found to have non-circular orbits with high radial velocities. In Fig. 23, we present the PDF of the cylindrical radial velocity of the gas in simulations G and H (which show a strong increase in SFR), weighted by the SFR given by equation (9), for gas with density n > 100cm⁻³. The grey lines denoting the PDF at t = 0 are symmetric about V − r = 0kms⁻¹. As the jet evolves (black lines), the velocity PDF develops an extended tail, up to Vr ≲ 500kms⁻¹. This is similar to Fig. 9, with smaller upper limit of the radial velocity due to the choice of the density threshold of n > 100cm⁻³. Although the simulations presented here lack explicit numerical schemes to track star formation in situ (as done in Gaibler et al. 2012; Dugan et al. 2014), the SFR-weighted velocity PDF demonstrates that positive feedback is expected to give rise to hyper-velocity stars.

Figure 23.

Probability density function (PDF) of the cylindrical radial velocity weighted by the SFR (given by equation 9, for n > 100cm⁻³), for simulations G and H. The grey line denotes the PDF at t = 0, which is symmetric about V_r = 0kms⁻¹. In black we show the PDF when the jet has evolved. The PDF is asymmetric, with a tail extending beyond V_r ≳ 500kms⁻¹. The PDF is an approximate indicator of the kinematics of the stars expected to form as a result of positive feedback.

4.2 Negative feedback

The above considerations of SFR are based purely on arguments of collapse induced by density enhancement. However, the SFR also significantly depends on the local turbulent velocity dispersion (e.g. see Federrath & Klessen 2012, for a review). The increased turbulence can provide additional support against collapse. True collapse will set in at spatial scales where turbulent velocity fluctuations are less than of the order of the local sound speed, such that global turbulent support does not affect the local dynamics (Vázquez-Semadeni, Ballesteros-Paredes & Klessen 2003; Federrath et al. 2010a; Federrath & Klessen 2012). Such scales are comparable to (or smaller than) the Jeans length¹¹ which our simulations do not resolve. However we note that the enhanced turbulence induced in the disc [as discussed earlier in Section 3.1, point (v)], will also add to turbulent support preventing collapse.

If we assume that the turbulence in the disc behaves according to the standard theory of a scale invariant cascade (e.g. Kolmogorov 1941), the velocity dispersion (σ_|$v$|) at smaller length scales l, depends on dispersion (σ_V) at larger scales L as: σ_|$v$| ∼ σ_V(l/L)^ζ (e.g. equation 12 of Federrath & Klessen 2012). Numerical simulations (Kritsuk et al. 2007; Schmidt et al. 2009) and measurements of velocity structure function in giant molecular clouds (Brunt 2003; Heyer & Brunt 2004) have shown that the exponent to vary is ζ ∼ 0.4–0.5. For length scales ≲ 1 pc, at which collapse may occur, the velocity dispersion would thus be σ_|$v$| ∼ 12kms⁻¹(σ_V/400kms⁻¹)(1pc/1kpc)^0.5. This is still much higher than the thermal sound speed expected in the cold molecular gas [c_s ∼ 0.2(T/10K)^1/2(μ/2)^−1/2kms⁻¹]. This implies that turbulent support may in fact prevent collapse at such scales resulting in negative feedback.

An often used indicator for ascertaining whether a gas disc in equilibrium will be prone to fragmentation and hence star formation is the Toomre parameter Q, defined as |$Q \sim \left(\frac{v_\phi }{r}\right) \sigma _v/\left(\pi G \Sigma \right)$|⁠. Here, σ_|$v$| is the velocity dispersion and Σ the density per unit area. We have approximated the epicyclic frequency as |$\kappa ^2=\frac{2 \Omega }{r}\frac{{\rm d}}{{\rm d}r}\left(r^2 \Omega \right) \sim \langle v_{\phi } \rangle /r$|⁠. A gas disc is will be gravitationally stable if Q > 1 (Toomre 1964; Wang & Silk 1994). If we consider the fractal disc in our simulations as a symmetric thick disc on average, the jet induced dispersion is expected to raise the Toomre Q, preventing further fragmentation. We present the cylindrically radial profile of Q for simulation H in Fig. 24. The starting velocity dispersion at the time of injection of the jet in simulation H was selected to be lower than in the other simulations (∼40kms⁻¹) so that initially, Q < 1. As the jet interacts with the disc, raising its dispersion, Q increases by a factor of 2 or more throughout the disc. Hence, on average, the increased velocity dispersion, even on the outer edges of the disc which are not directly affected by the central jet-driven outflow, may render the disc stable to fragmentation.

Figure 24.

The distribution of the Toomre parameter at different times for simulation H. Starting from an initial Q < 1, the jet driven increase in dispersion enhances the Q to values greater than 1.

Thus, to conclude, density enhancement near the jet axis can make the clouds more prone to gravitational collapse promoting star formation. However, overall, an increase in turbulent kinetic energy may lower the SFR. This indicates that both positive and negative feedback may operate in galaxies impacted by relativistic jets. Simulations with higher spatial resolution are needed to better assess the relative importance of these two forms of feedback.

5 SUMMARY AND DISCUSSION

In this paper, we have studied the interaction of relativistic jets and the energy bubbles they create with dense turbulent gas discs. We have carried out a suite of simulations with different jet powers, densities and jet inclinations. Our primary focus has been in understanding the evolution of the kinematics, morphology, and energetics of the multiphase ISM under the influence of the relativistic jet and its implications for the evolution of the galaxy properties such as SFR. Below we summarize the results and the their implications.

• Jet-driven outflows: The jets launch an energy bubble that drives shocks into the disc. The jet couples with the ISM driving fast multiphase outflows. Clouds are dredged up vertically along the axis of the jet with outflows velocities ∼500kms⁻¹. The dense clouds accelerated by the jets form cometary tails as they are ablated by the fast flows. Typically, a cloud is expected to survive a few cloud crushing times (Fragile et al. 2004; Scannapieco & Brüggen 2015), where the cloud crushing time is the time for the incumbent shock to pass the length of the cloud. For a shock of speed |$v$|_sb ∼ 1000kms⁻¹, and density contrast of χ = 10⁴ between the up-stream and down-stream fields (which would occur for a bubble density of n_b ∼ 0.01cm⁻³ and cloud density of n_c ∼ 100cm⁻³), the cloud crossing time (t_cc) is (Fragile et al. 2004)

  \begin{eqnarray*} v_\text{sc} &= \left(\frac{n_\text{b}}{n_\text{c}}\right)^{1/2} v_\text{sb} \end{eqnarray*}

  (10)

  \begin{eqnarray*} &= \chi ^{-1/2} v_\text{sb} \end{eqnarray*}

  (11)

  \begin{eqnarray*} t_\text{cc} &= \frac{L_\text{c}}{v_\text{sc}} \end{eqnarray*}

  (12)

  \begin{eqnarray*} &= \chi ^{1/2} \frac{L_\text{c}}{v_\text{sb}} \end{eqnarray*}

  (13)

  \begin{eqnarray*} & \sim 1 \mbox{ Myr } \left(\frac{\chi }{10^4}\right)\left(\frac{L_\text{c}}{100 \mbox{ pc}}\right) \left(\frac{v_\text{sb}}{1000 \mbox{ km s}^{-1}}\right). \end{eqnarray*}

  (14)

  This is comparable to the time-scales of the simulation. In our simulations, the gas ejected vertically from the disc, is mostly shock heated to high temperatures, and breaks up to form diffuse cometary tails. The ablated gas being lighter (n ∼ 0.1–10cm⁻³) is swept up and attains much higher velocities of |${\sim } 1000 \mbox{ km s}^{-1}$|⁠. However, we note that to attain numerical convergence of jet–cloud interaction, a much higher resolution is required, with the clouds resolved by at least 100 computational cells (Banda-Barragán et al. 2016). As the jet breaks out from the disc, the energy bubble spreads over the disc, driving shocks inward into the disc and compressing them.

• Non-circular motions in the disc: Besides vertical outflows, the lateral expansion of the jet-driven energy bubble also launches radial outflows into the plane of the disc. The jet-driven outflows result in a departure of the gas kinematics from the circular rotation initially established (as shown in Fig. 10). The strong radial motions of the gas has two observational signatures. First, lines of sight probing the disc would observe extreme gas kinematics as highly broadened emission lines. Significant distortion in the galaxy’s rotation curve in a position–velocity diagram would also be observed. Such features have been observed in detail for the galaxy IC 5063 (Dasyra et al. 2015, 2016; Morganti et al. 2015; Oosterloo et al. 2017). Detailed simulations targeted towards explaining the observed kinematics in IC 5063 have been presented in a separate publication (Mukherjee et al. 2018).

  Secondly, if the radially driven gas cools to form stars, non-circular motions would leave a population of stars with strongly elliptical orbits. Our present work lacks the numerical machinery to address the star formation properties of the gas clumps. However, similar suggestions have been made in recent papers Dugan et al. (2014), raising the possibility that hyper-velocity stars (or gas parcels) with non-circular radial orbits being are a possible signature of AGN activity.

• Turbulence inside the disc: The jets drive shocks and hydrodynamical flows that penetrate the disc, increasing its turbulence. The turbulence in the shocked gas increases by four to six times its initial value (as shown in Fig. 12). This increase occurs not just in the central nuclear region where the jet–ISM interaction is expected to be strongest, but all throughout the disc. This implies that the onset of jet activity produces large-scale effects over several kpc. For a disc with an identical initial ISM, jets with power 10⁴⁶ergs⁻¹ induce velocity dispersions higher by a factor ∼1.5 than a jet with power 10⁴⁵ergs⁻¹. Lower density clouds are more prone to ablation, resulting in higher velocity dispersion in ISM with lower mean density, for the same jet power. Our simulations support the kinematic observations of jet–ISM interactions in galaxies such as 3C 326 (Nesvadba et al. 2010, 2011) and 3C 293 (Mahony et al. 2016b).

• Inclined jets: We have carried out simulations of jets inclined to the disc at various angles of inclination. We find that an inclined jet behaves very differently from a jet launched perpendicular to the axis. Being inclined to the disc, the jets encounter a larger column depth of gas with which they interact before breaking out. The inclined jets decelerate considerably and launch a sub-relativistic wind along the galaxy’s minor axis following the path of least resistance. Such a wind is reminiscent of the wide-angled outflow often associated with AGN winds (Rupke & Veilleux 2011) or nuclear starburst (Veilleux et al. 2005).

  Similar results of AGN outflows inclined to the disc have recently been explored in simulations of Dugan et al. (2017). The starting outflow in Dugan et al. (2017) was assumed to be non-relativistic (though ultra-fast, akin to Wagner, Umemura & Bicknell 2013). However, as we demonstrate in our simulations, such wide-angled outflow may also arise from an inclined relativistic jet itself. The ensuing bubble evolves spherically, in contrast to the directional outflows seen in powerful radio jets. Often the difference in morphology of the observed winds is used to differentiate between a wind-driven feedback from a jet-driven feedback (Rupke & Veilleux 2011). However, our simulations show that even a relativistic jet can give rise to a morphology that that may be otherwise attributed to a quasar-driven wind.

  An inclined jet also perturbs the gaseous disc more strongly. The velocity map in Fig. 14 shows a gas distribution resembling a wide-angled outflow, compared to an inner outflow followed by compression at outer edges as seen in the case of a jet perpendicular to the disc. The stronger jet–ISM interaction results in an increased velocity dispersion and a higher mean radial velocity with a larger value of the inclination angle. The jet-driven bubble of an inclined jet expands more slowly, than that of a perpendicular jet, as the bubble for the inclined jet is powered by a sub-relativistic wind. The ensuing morphology of the shocked gas and bubbles have strong similarities with some observed systems such as NGC 1052 and NGC 3079.

• Positive versus negative feedback

  • Positive feedback: In Section 4, we have estimated the impact of the jet on the SFR in the galaxy. Assuming a Kennicutt–Schmidt star formation law (Kennicutt 1998), we find that regions close to the jet axis may experience an increase in SFR resulting from the enhanced density produced by strong radiative shocks. Simulations with P_j = 10⁴⁵ergs⁻¹ show only a modest increase in star formation (e.g. by ∼6 per cent for simulation B), whereas higher power jets (P_j = 10⁴⁶ergs⁻¹) show a significant overall increase in SFR (e.g. by 80 per cent for simulation H). The increase of SFR with higher jet power is due to the enhanced compression of the gas by the stronger shocks. At lower jet powers, the ablation by the outflow wins over shock compression. Similar results were also shown in Zubovas & Bourne (2017), where it was found that cloud fragmentation rate (a proxy for star formation) peaks at a critical AGN luminosity. Inclined jets, with stronger interactions with the ISM, are also seen to exhibit an increase in star formation.

    Jet-induced star formation is a viable mechanism to explain enhanced star formation in several observed sources with star-forming gas clumps located near a radio jet, e.g. the Minkowski object (Salomé et al. 2015; Lacy et al. 2017), 3C 285 (Salomé et al. 2015), 4C 41.17 (Bicknell et al. 2000). There is also evidence of enhanced SFR in galaxies hosting an AGN (Zinn et al. 2013; Bernhard et al. 2016). Several other theoretical works have also explored the possibility of AGN-induced star formation in galaxies, such as Nayakshin & Zubovas (2012), Wagner & Bicknell (2011), Gaibler et al. (2012), Dugan et al. (2017) etc., where compression of gas by the AGN-driven outflow can trigger star formation, similar to what we find.

    The enhancement in the SFR surface density is highest close to the jet axis, appearing as a circum-nuclear star-forming ring (≲ 1kpc). Similar results have also been presented by Gaibler et al. (2012). If such positive feedback is indeed at play, a star-forming ring with a young stellar population around the jet axis is a possible observational signature for such systems.

    Positive feedback can result in the formation of hyper-velocity stars, and more generally and stars with significant radial orbital perturbations, with kinematics significantly different from the rest of the stellar population (see Fig. 23). Recently, star formation (⁠|${\sim } 15\, \text{M}_\odot \, \mbox{yr}^{-1}$|⁠) has been detected inside a galactic outflow (Maiolino et al. 2017) with stellar velocities higher than that of the stars in the centre (by about ≲ 100kms⁻¹). Similar kinematic imprint of hyper-velocity stars resulting from star formation triggered by an AGN-driven outflow is likely to be a general feature of star-forming galaxies with a central active supermassive black hole (Silk et al. 2012; Zubovas et al. 2013; Zubovas & Bourne 2017; Wang & Loeb 2018). If a star-forming galaxy undergoes several episodes of AGN/jet activity (e.g. Rampadarath et al. 2018) which also result in positive feedback, such galaxies will have a population of stars whose formation has been triggered by positive feedback, with radial velocity asymmetries (as apparently found in stacked SDSS star-forming galaxies Cicone, Maiolino & Marconi 2016). Although supernova-driven outflows can also impact the stellar kinematics (e.g. El-Badry et al. 2017, for low-mass galaxies), the kinematic signature of star formation induced by an AGN-driven outflow will be different, in having a strong radial motion imprint of the newly formed stars. These dynamical signatures will accumulate and be reinforced over successive AGN episodes.

  • Negative feedback: In our simulations, we do not find any appreciable mass-loss due to outflows which would result in significant quenching. However, the increased local velocity dispersion in the ISM induced by the jet-driven shocks can provide additional turbulent support preventing formation of stars. Such an effect is not captured by the positive feedback estimates inferred from a simple Kennicutt–Schmidt star formation law that solely depends on the local density (e.g. equation 9, as well as in Gaibler et al. 2012; Bieri et al. 2016; Dugan et al. 2017). In our simulations we find that turbulence is strongly enhanced (∼400–600kms⁻¹) to several times the initial values. The dispersion computed here is on scales of the global disc (∼1 kpc).

    Observational support of turbulence-driven quenching by relativistic jets has been recently presented for some galaxies, e.g. 3C 326 (Nesvadba et al. 2010), a large fraction of the Molecular Hydrogen Emission Galaxy (Ogle et al. 2010), NGC 1266 (Alatalo et al. 2015), etc. Although high-velocity outflows have been detected in several of these galaxies (Nesvadba et al. 2010; Ogle et al. 2010; Guillard et al. 2012), they do not carry enough mass to unbind the gas reservoir completely. Such a scenario is also true for NGC 1266 with a jet of much lower power (Nyland et al. 2013) that does not drive outflows that escape the galaxy (Alatalo et al. 2011). This points to a quenching mechanism in which the jet can drive turbulence on scales comparable to the gas extent of the galaxy, as in our simulations.

  The current simulations represent only a few million years of the lifetime of a galaxy. Such time-scales are much smaller than the dynamical time of the galaxies. Hence, the estimates of positive or negative feedback presented here should be considered more as a singular event associated with one episode of jet activity. Long-term effects of jet-driven feedback on a disc pressurized by the jet cocoon have also been recently explored by Bieri et al. (2015, 2016). Bieri et al. (2016), who keep the disc artificially pressurized for ∼400 Myr, resulting in enhanced fragmentation of the disc and subsequent star formation. However, in a realistic system, the pressure in the energy bubble will decrease with time as a result of its expansion, as shown in Fig. 6. Thus, it is still uncertain whether there can be sustained positive feedback to significantly affect a galaxy’s stellar mass assembly.

  We conclude that an onset of jet activity may induce positive feedback in regions directly being impacted by the jet, resulting in localized regions of enhanced star formation during early phases of the jet evolution. Over larger spatial extents, enhanced turbulence from the jet-driven energy bubble may suppress star formation for the duration of jet activity. The above inferences are indirect, as our simulations currently lack the numerical tools to track potential star-forming clumps and lack the resolution to model turbulence in the star-forming clumps. In future we aim extend the numerical set-up to address these questions in greater detail.

ACKNOWLEDGEMENTS

We gratefully acknowledge Nicole Nesvadba, Christoph Federrath, Michael Dopita, Elizabeth Mahony, Sylvain Veilleux, and Luisa Ostorero for helpful discussions. We thank the referee for the thorough review and insightful comments. The simulations have been carried out in the supercomputing facilities at NCI and Pawsey, made available to researchers at ANU under the NCMAS and ANU partner share schemes. We thank the timely help from the NCI and Pawsey support staff. DM’s visit to JHU was supported by a Balzan grant from New College, Oxford. AYW has been supported in part by ERC Project No. 267117 (DARK) hosted by Université Pierre et Marie Curie (UPMC) – Paris 6, PI: J. Silk.

Footnotes

1

Resolution of successive cells in the Z direction are Δ|$z$|^{k + 1} = r_s × Δ|$z$|^k.

2

https://pypi.python.org/pypi/pyFC.

3

This corresponds to setting k_min = 6 as an input parameter to the pyFC cube generator.

4

Simulations A, B, F, G, and H, as in Table 2.

5

High-resolution simulations of wind-filament interaction (Banda-Barragán et al. 2016, 2018) show that numerical convergence requires 120 resolution elements along the filament.

6

The definition of the anisotropy parameter employed here is better suited to describe the non-radial motion in a gas disc. This is different from the standard definition of the anisotropy parameter, often used to quantify stellar kinematics in elliptical galaxies (e.g. Binney & Tabor 1995).

7

See Li, Klessen & Mac Low (2003) for a discussion on mass- versus volume-weighted PDFs.

8

It must be noted here that the mass-weighted mean radial velocity is averaged over the entire ISM. This includes the undisturbed dense cores and, thus, does not truly reflect the velocity of individual out-flowing filaments, which can be much higher, as can be seen in Figs 4, 14, and 15.

9

Low ionization nuclear emission line region (Heckman 1980).

10

NGC 3079 is classified as a Seyfert 2 (Ford et al. 1986) with a LINER nucleus.

11

The Jeans length is

\begin{equation*} \lambda _\text{J}=(\pi c_\text{s})/\sqrt{G \rho } \sim 1.38 \mbox{ pc } \left(\frac{c_\text{s}}{0.2 \mbox{ km s}^{-1}}\right)\left(\frac{n}{10^3 \mbox{ cm}^{-3}} \right)^{-1/2}, \end{equation*}

where |$c_\text{s}=\left(\frac{\gamma k_\text{B} T}{\mu m_a}\right)^{1/2}$| is the speed of sound. A value of μ ∼ 2 has been assumed.

12

During apodization computational cells with T > 3 × 10⁴ K are replaced with ambient gas. See Section 2.2 for details.

APPENDIX A: BUBBLE DYNAMICS

In this appendix, we discuss the self-similar expansion of a bubble driven by a nuclear wind expanding adiabatically into an ambient medium with a power-law density profile. The bubble has a radius R_B, volume V_B, and pressure p_B which is assumed to be uniform (as is typical of an overpressured bubble). The bubble is powered by a nuclear wind with energy flux F_E. The bubble is expanding into an ambient medium, whose density varies as a power law beyond a scale radius R_c (typically the core radius of the galaxy):

\begin{eqnarray*} \rho _a &= \rho _\text{c} \left(\frac{R}{R_\text{c}}\right)^{-\beta }. \end{eqnarray*}

(A1)

Thus, the total mass enclosed in a radius R_B (for halo with β ≤ 2) is

\begin{eqnarray*} \int _0^{R_\text{B}} \rho _a 4 \pi R^2 \text{d}R &= \int _0^{R_\text{B}} \rho _\text{c} \left(\frac{R}{R_\text{c}}\right)^{-\beta } 4 \pi R^2 \text{d}R \nonumber \\ &=\frac{4}{3-\beta } \pi R_\text{B}^3 \rho _\text{c} \left(\frac{R}{R_\text{c}}\right)^{-\beta }. \end{eqnarray*}

(A2)

We consider the bubble radius R_B to be defined by the forward shock. The velocity of the bubble radius is then |$v$|_B = dR_B/dt. Neglecting the mass-inflow by the wind such that the mass interior to the forward shock is comprised of swept up gas, the mass and momentum equations can be written as

\begin{eqnarray*} \frac{{\rm d}}{{\rm d}t} \int \rho _a v_\text{B} \text{d}^3x &= \int p \hat{r}\cdot \text{d}S, \end{eqnarray*}

(A3)

\begin{eqnarray*} \frac{{\rm d}}{{\rm d}t}\left(\frac{4}{3-\beta } \pi R_\text{B}^3 \rho _\text{c} \left(\frac{R_\text{B}}{R_\text{c}}\right)^{-\beta } \frac{{\rm d} R_\text{B}}{{\rm d}t} \right) &= 4 \pi R_\text{B}^2 p_\text{B}, \end{eqnarray*}

(A4)

where we have inserted equation (A2) in equation (A3). With an energy flux F_E into the bubble, the energy equation can be written as

\begin{eqnarray*} &\frac{{\rm d}}{{\rm d}t} \left[\frac{4 \pi }{3} \frac{p_\text{B}}{\gamma -1} R_\text{B}^3 \right] + 4 \pi R_\text{B}^2 p_\text{B} \frac{{\rm d} R_\text{B}}{{\rm d}t}=F_E. \end{eqnarray*}

(A5)

For a self-similar expansion, the bubble radius is a power law in time: R_B = At^α. Inserting this into equation (A4) we find the time evolution of the pressure as

\begin{eqnarray*} p_\text{B}=\frac{\rho _\text{c} \alpha }{(3-\beta) R_\text{c}^{-\beta }} A^{2-\beta } \left[(4-\beta)\alpha -1\right] t^{(2-\beta)\alpha -2}. \end{eqnarray*}

(A6)

Inserting equation (A6) in the energy equation (equation A5), and assuming the energy flux to be time-independent, we obtain

\begin{eqnarray*} \alpha=\frac{3}{5-\beta }, \end{eqnarray*}

(A7)

with the coefficient A given by

\begin{eqnarray*} A=\left[\frac{R_\text{c}^{-\beta } (3-\beta)}{6 \pi \rho _\text{c}}\frac{(5-\beta)^3}{(7-2\beta)(11-\beta)} F_E \right]^{1/(5-\beta)}. \end{eqnarray*}

(A8)

Thus, the expressions for evolution of radius, velocity of the bubble surface, and pressure are

\begin{eqnarray*} R_\text{B} &= \left[\frac{1}{6 \pi }\frac{(5-\beta)^3(3-\beta)}{(7-2\beta)(11-\beta)} F_E \right]^{1/(5-\beta)} \nonumber \\ & \times R_\text{c}^{-\beta /(5-\beta)} \, \rho _\text{c}^{-1/(5-\beta)} \, F_E^{1/(5-\beta)} \,\, t^{3/(5-\beta)}, \end{eqnarray*}

(A9)

\begin{eqnarray*} V_\text{B} &= \frac{3}{5-\beta } \left[\frac{1}{6 \pi }\frac{(5-\beta)^3(3-\beta)}{(7-2\beta)(11-\beta)} F_E \right]^{1/(5-\beta)} \nonumber \\ & \times R_\text{c}^{\beta /(5-\beta)} \, \rho _\text{c}^{-1/(5-\beta)} \, F_E^{1/(5-\beta)} \,\, t^{-(2-\beta)/(5-\beta)}, \end{eqnarray*}

(A10)

\begin{eqnarray*} p_\text{B} &= \frac{3(7-2\beta)}{(3-\beta)(5-\beta)^2} \left[\frac{(5-\beta)^3(3-\beta)}{6\pi (7-2\beta)(11-\beta)}\right]^{(2-\beta)/(5-\beta)} \nonumber \\ &\times R_\text{c}^{-3 \beta /(5-\beta)} \, \rho _\text{c}^{3/(5-\beta)} \, F_E^{(2-\beta)/(5-\beta)} \, t^{-(\beta + 4)/(5-\beta)}, \end{eqnarray*}

(A11)

\begin{eqnarray*} p_\text{B} &= \frac{3(7-2\beta)}{(3-\beta)(5-\beta)^2} \left[\frac{(3-\beta)(5-\beta)^3}{6 \pi (7-2\beta)(11-\beta)} \right]^{2/3} \nonumber \\ &\times \rho _\text{c}^{1/3} R_\text{c}^{\beta /3} F^{2/3} R_\text{B}^{-(\beta +4)/3}. \end{eqnarray*}

(A12)

For β = 0, the above equations reduce to the expressions derived in Weaver et al. (1977) for a uniform external medium. The exponent of the time for the evolution of the bubble pressure is similar to that of Kaiser & Alexander (1997), who also have studied the self-similar evolution of a jet.

Hence, for β ∼ 1.5–2, which is typical for the double isothermal potentials used in our simulations, the bubble pressure will thus be ∝t^−1.57 or ∝t⁻², respectively. For simulation D, which launches a spherical sub-relativistic wind, the pressure initially (t < 0.4 Myr) varies as ∼t^−1.21, with a later evolution of ∼t^−1.54, close to the predicted bubble solution. The initial slower decrease occurs because the bubble is trapped inside the disc and suffers radiative losses. Vertical jets show a faster and non-spherical expansion, resulting in a faster decrease in the bubble pressure (as shown in Fig. 6).

The strength of the shock driven into the ambient medium can be approximately estimated by equating the ram pressure of the shock to that of the pressure in the energy bubble.

\begin{eqnarray*} p_{\rm bubble} &= \rho _a v_{\rm shock}^2, \end{eqnarray*}

(A13)

\begin{eqnarray*} p_{\rm bubble} &\sim 8.95\times 10^{-8} \!-\! 8.95 \times 10^{-9} \mbox{dyne cm}^{-2} \end{eqnarray*}

(A14)

\begin{eqnarray*} & ({\mbox{for } {t}} \sim 0.1\!-\!0.5 \mbox{ Myr respectively, } \mbox{ for simulation B}) \nonumber \\ n_a &\sim 0.1 \mbox{ cm}^{-3}\mbox{(ambient density in 2-4 kpc)}, \end{eqnarray*}

(A15)

\begin{eqnarray*} v_{\rm shock} &\simeq 3\times 10^3 \!-\!10^4 \mbox{ km s}^{-1}, \end{eqnarray*}

(A16)

\begin{eqnarray*} \mathcal {M} &= v_{\rm shock}/\left(\frac{K T_a}{\mu m_a} \right)^{1/2} \simeq 8.1\!-\!26. \end{eqnarray*}

(A17)

Another approximate estimate ofthe shock velocity can be obtained from the rate of expansion of the volume of the bubble. For a spherically expanding bubble, the rate of the expansion of the bubble is |$\text{d}r_\text{b}/\text{d}t=1/(4 \pi r_\text{b}^2) \text{d}V_\text{b}/\text{d}t=1/(36 \pi V_\text{b}^2) \text{d}V_\text{b}/\text{d}t$|⁠, where |$V_\text{b}=(4 \pi r_\text{b}^3)/3$| is the volume of the bubble. In Fig. A1, we show the dr_b/dt for simulations B and F, plotted until the bubble reaches the edge of the simulation domain. As can be seen, for simulation B the rate of expansion up to ∼0.5 Myr is ∼0.012c = 3.6 × 10³kms⁻¹, close to the approximate analytical estimates derived above. The rate of expansion of the bubble in simulation F is slower, as explained in Section 3.2. Estimates of shock velocities obtained from observational modelling of emission lines observed in gas clouds embedded in the cocoon (Sutherland, Bicknell & Dopita 1993; Bicknell et al. 2000), yield values comparable to those derived above.

Figure A1.

Rate of expansion of the bubble for simulations B and F computed as |$\text{d}r_\text{b}/\text{d}t=1/(4 \pi r_\text{b}^2) \text{d}V_\text{b}/\text{d}t=1/(36 \pi V_\text{b}^2) \text{d}V_\text{b}/\text{d}t$|⁠, where |$V_\text{b}=(4 \pi r_\text{b}^3)/3$| is the volume of the bubble.

APPENDIX B: VELOCITY POWER SPECTRUM

In Fig. B1, we show the evolution of the velocity power spectrum computed for a domain of size 1kpc³, centred at coordinates (1,1,0). The domain was chosen to be away from the central axis of the jet which experiences strong outflows, and to be well within the height of the disc, such that a Fourier transform with periodic boundary conditions may be performed. For the settling disc (left-hand panel in Fig. B1), we note that the slope of the power spectrum at t = 0 is ∼k^−1.2. Although the turbulent velocity field is initially created with a slope of k^−5/3, the apodization with the warm medium when being interpolated on to the pluto grid¹² results in departure from a 5/3 slope.

Figure B1.

Velocity power spectrum for settling disc (left), simulation F (middle), and simulation H (right) at different times.

As the disc is settled, the slope of the power spectrum converges to a Kolmogorov slope (∼k^1.666) in the inertial range of 3 < k < 15 (corresponding to physical scales of ∼333–66pc, respectively). The jet when injected into the settled disc thus initially experiences a settled medium with a Kolmogorov power spectra in the inertial range. With the onset of the jet (middle and right-hand panels of Fig. B1), the power spectra flattens significantly in the range 3 < k < 15, with a mean slope of ∼k⁻¹. This implies that the said range of wavenumbers do not represent the inertial range of turbulent cascade anymore. This is because the jet and energy bubble forces turbulence on the gas at such scales which are comparable mean cloud sizes (∼100pc) . At higher wavenumbers, the power spectrum may be affected by numerical dissipation and hence the choice of resolution (Kritsuk et al. 2007). Higher resolution simulations need to be performed to investigate the nature of the cascade of turbulent energy driven by the jet to smaller length scales.