AGN Outflow Shocks on Bonnor–Ebert Spheres

The role active galactic nuclei (AGNs) play in the evolution of their host galaxies remains an outstanding problem in theoretical astrophysics (Silk & Rees 1998). The discrepancy between the number of luminous galaxies predicted in a ΛCDM universe and that observed by astronomers persists and is mainly explained by negative feedback through AGN quenching (Weinmann et al. 2006). In many of today's cosmological simulations, this feedback is treated through the injection of thermal energy into the galaxy (e.g., Springel et al. 2005; Sijacki et al. 2007), others actually include jets and outflows (Kimm et al. 2012; Dubois et al. 2010), and both contribute to achieving the necessary quenching to match the observed luminosity functions.

Observers still investigate the different impacts of jets and outflows on the host galaxy. Recent observational evidence shows greater star formation rates (SFRs) in radio-loud quasars than in radio-quiet quasars. Kalfountzou et al. (2012) analyze approximately 20,000 quasars from the Sloan Digital Sky Survey (SDSS) and employ the [O ii] emission line to calculate SFRs. They find higher SFRs in the radio-loud AGNs, indicating AGN jet–triggered star formation. Kalfountzou et al. (2014) utilize Herschel-ATLAS (Eales et al. 2010) and find that the far-infrared (FIR) data are critical to calculating SFRs in AGNs, particularly in quasars where radiation from the accretion disk can contaminate calculations of the SFRs calculated from optical and ultraviolet (UV) data. Once the FIR data is incorporated into the SFR calculations, the group finds both that high-redshift radio-loud AGNs can have high SFRs and that these bursts of star formation happen simultaneously with the growth of supermassive black holes. These results are in tension with the idea of universal AGN quenching. Similarly, Zinn et al. (2013) study several hundred AGNs from the Chandra Deep Field South, compare the SFRs as determined from archived IR data. They find higher SFRs in the radio-loud AGNs, concluding that AGN jets trigger star formation. The differences in SFRs between radio-loud and radio-quiet quasars underscore the importance of understanding the feedback mechanisms on the scale of actual gas clouds rather than galactic scales.

However, with all the restrictions inherent to cosmological simulations on large physical scales, many of the critical physical processes in AGN feedback cannot be resolved due to density, velocity, and related timescales and dynamic ranges. Gaibler et al. (2012) simulate AGN jet feedback on a smaller, single-galaxy scale and show that mildly relativistic jets may induce star formation in their host galaxies. The discrepancy between the two approaches requires a closer look on smaller time and physical scales to treat the problem.

Wagner et al. (2016) also explore both positive and negative AGN feedback from a theoretical perspective and identify the challenges of simulations. One such problem is the huge range of physical scales important to AGN feedback, scaling from the pc size of the accretion disk, out to tens of kpc to the edges of the host galaxies, and extending even further to Mpc scales in interaction with the intergalactic medium. Covering these spatial scales with good resolution while simultaneously capturing the extreme velocities close to the speed of light is virtually impossible. Including self-gravity and detailed radiative processes would make this computationally prohibitive. It is important, then, to focus on specific aspects of AGN feedback with a given suite of simulations. As an example, Bieri et al. (2016) examine the feedback effect from the pressure confinement of a jet-driven bubble on a disk radio galaxy, avoiding the computational costs of an actual jet by artificially introducing overpressure. This frees up computational power to include self-gravity. They find that the pressure causes an increased fragmentation of dense clouds in the host galaxy and a subsequent increase in star formation (positive feedback).

Another specific aspect of AGN feedback that is not well understood is the impact of winds from expanding jet cocoons or outflow bubbles on dense clouds of gas on small scales (i.e., cloud scales), though there is some observational evidence of this. Russell et al. (2016) analyze observations from Atacama Millimeter/submillimeter Array (ALMA) and finds that AGN jets can blow giant molecular clouds (GMCs) away from the galactic nucleus with jet-driven bubbles before gravity draws them back, similar to mechanical pumps. These observations reveal star formation along the outer edges of the clouds that is potentially triggered by the expansion of the jet bubble. Employing observations from the Measuring Active Galactic Nuclei Under MUSE Microscope (MAGNUM) survey on NGC 5643, a radio-quiet AGN with outflows, Cresci et al. (2015b) find positive feedback. The observations show star formation in clumps within the double-sided ionization cones with high-velocity gas, likely caused by the compression of the clouds from the outflow. They also find a ring of star formation at 2.3 kpc, possibly triggered by the AGN outflow, a phenomenon also seen in Gaibler et al. (2012) and Dugan et al. (2014).

Simulations of AGN feedback on smaller scales provide a complementary view of the observations, as well as of large-scale simulations. The fate of the clouds and hence potential star formation is eventually decided on small scales. Mellema et al. (2002), Cooper et al. (2009), and Zubovas et al. (2014) have performed important studies to this end. Understanding the relationship between these winds and the clouds of gas available to form stars in a galaxy is critical to understanding the role AGN feedback has in galaxy evolution.

In this paper, we seek to build on past research of AGN feedback on clouds of gas by focusing on realistic AGN winds and gravitational collapse on the scale of a single cloud. We simulate high-velocity winds characteristic of AGN jet cocoons and outflow bubbles striking Bonnor–Ebert (BE) spheres (Ebert 1955; Bonnor 1956). We have better resolution than in the past, and we include gravity, heating, and cooling for the entire temperature range (including very low temperatures).

We organize this paper as follows: in Section 2, we discuss some previous literature; in Section 3, we describe our simulation setup; and, in Section 4, we give our results. In Section 5, we include our findings, and we conclude in Section 6.

A large number of cloud-crushing simulations have been performed in the past (of which we can only mention a selection), but few have employed values for the pressure, velocity, and density of the incoming shock wave that are realistic for AGN feedback, as well as some crucial physical processes. Adiabatic cloud-crushing simulations may be scaled to the velocities and densities relevant for AGN feedback. However, this is not possible when nonlinear processes, such as gravitational collapse or cooling, are included and have a significant impact. Many of the simulations in the literature have not included gravity, so exploring its stabilizing effect, subsequent collapse, and potential star formation has been impossible for them. We seek to include both realistic AGN wind parameters and gravity in order to probe whether or not a wind resulting from an AGN jet cocoon or an AGN outflow bubble can cause the collapse of the dense clouds of gas in the interstellar medium (ISM) and form stars.

Few studies in the past have specifically focused on the impact of radio jet cocoons on a cloud of gas. Mellema et al. (2002) investigate this topic through 2D simulations. They list three main phases of wind-cloud interaction, first described by Klein et al. (1994). In the first phase, the blast wave overtakes and interacts with the cloud. In the second phase, the cloud is compressed through ram pressure and possibly thermal pressure if it is underpressured with respect to the environment. This phase lasts for the cloud crushing time,

$$\begin{eqnarray}&&{t}_{\mathrm{cc}}={\chi }^{1/2}{R}_{{\rm{c}}}/{v}_{{\rm{w}}},\end{eqnarray} \tag{ 1 }$$

where ρ_c and ρ_w are the cloud and wind densities, respectively; $\chi ={\rho }_{{\rm{c}}}/{\rho }_{{\rm{w}}}$; R_c is the cloud radius; and v_w is the wind velocity. In the third phase, shocks from different sides of the cloud collide, forming a rarefaction wave back through the shocked gas, a process that typically takes a few cloud crushing times.

Mellema et al. (2002) use initial values for the cloud of ρ = 10 m_p cm⁻³, temperature 10⁴ K, M = 9.5 × 10⁵ M_⊙, and size 200 pc, and for the wind they choose v_w = 3500 km s⁻¹. Their cloud parameters are meant to resemble those of a GMC. They show that, with radiative cooling, the cooling time for the cloud is substantially shorter than both the cloud crushing time and the Kelvin–Helmholtz (KH) time,

$$\begin{eqnarray}&&{t}_{\mathrm{KH}}={R}_{{\rm{c}}}({\rho }_{{\rm{c}}}+{\rho }_{{\rm{w}}})/({v}_{{\rm{c}}}-{v}_{{\rm{w}}}){({\rho }_{{\rm{c}}}{\rho }_{{\rm{w}}})}^{1/2},\end{eqnarray} \tag{ 2 }$$

where v_c is the cloud velocity. They simulate uniform gas distributions in the shape of a circle and an ellipse and find that, at the conclusion of both simulations, less than 1% of the original cloud has mixed in with the ambient gas. This indicates the slow evaporation of the cloud under the conditions of a radio jet cocoon. They also find that radiative cooling occurs quickly enough that most of the cloud collapses into dense filaments that are likely to form stars, agreeing with the scenario of jet-induced star formation. However, self-gravity is not simulated, and star formation is not explicitly calculated.

Some simulations of shocks resulting from other astrophysical phenomena cover parameter spaces pertinent to AGN feedback. For example, Orlando et al. (2005) simulate a supernova shock wave passing over a small uniform spherical gas cloud and find that thermal conduction and radiative cooling are critical to cloud survival. They use two different scenarios, each with different velocities and temperatures: wind velocities of 250 and 430 km s⁻¹, densities of 0.4 m_p cm⁻³ in both, and temperatures of 1.6 × 10⁶ and 4.7 × 10⁶ K. The sphere has a density of 1 m_p cm⁻³ and a temperature of 10³ K, while the ISM has a density of 0.1 m_p cm⁻³ and a temperature of 10⁴ K. They show the importance of thermal conduction in suppressing destructive instabilities in the first case, though the cloud does expand and eventually evaporate. In the second case, they show that radiative cooling is pivotal to the formation of cold, dense clumps after the wind passes over the cloud. Gravity is not included, and no star formation is calculated.

Also relevant to AGN feedback, Cooper et al. (2009) explore shocks from starbursts on fractal clouds and uniform spheres alike and address the importance of radiative cooling. They use wind values of ρ = 0.1 m_p cm⁻³, v = 1.2 × 10⁵ cm s⁻¹, and T = 5 × 10⁶ K and sphere values of r = 5 pc, M = 523 M_⊙, $\bar{\rho }=91$ m_p cm⁻³, and T = 5 × 10³ K. The maximum resolution is 0.13 pc. The winds form dense filamentary structures like those observed in starburst-driven winds. They also find that radiative cooling is essential to the longer survival time of the clouds and filaments. The wind does not heat or accelerate radiative clouds as much as it does their adiabatic counterparts, decreasing the impact of KH instabilities and increasing the lifetime of the clouds and filaments. The geometry of the cloud does impact their results, and the wind destroys less-dense regions more quickly. While they indicate that self-gravity may cause the filaments to collapse, they do not include gravity or look at subsequent star formation in this study.

Pittard et al. (2010) perform 2D cloud-crushing simulations with velocities relevant to AGN feedback and find that KH instabilities dominate in the cloud's destruction. They assume values of r = 2 pc, central density ρ_c = 60 m_p cm⁻³, and temperature of 8000 K, and the outer edge of the cloud and the ambient gas are in pressure equilibrium with the surrounding gas of density ρ_a = 0.24 m_p cm⁻³ and temperature T_a = 8 × 10⁶ K. Cloud density follows the profile $\rho (r)={\rho }_{\mathrm{amb}}\left[{\rm{\Psi }}+(1-{\rm{\Psi }})\times \tfrac{1}{2}\left(1+\tfrac{\alpha -1}{\alpha +1}\right)\right]$, where ρ_amb is the ambient density, Ψ is a parameter selected to create a specific density contrast between the center and edge of the sphere, $\alpha =\exp [\min {(20,10((r/{r}_{{\rm{c}}})}^{2}-1))]$, and r_c is the cloud radius. They simulate velocities of 650, 1300, and 4300 km s⁻¹. They also find the formation of dense filaments in many of the simulations and note that the presence of turbulence in the post wind speeds up the destruction time of the cloud. However, they do not include gravity.

Zubovas et al. (2014) perform smoothed particle hydrodynamics (SPH) simulations aimed at determining the role of AGN outflow or jet-generated external pressure on turbulent, spherically symmetric molecular clouds, as well as fractal clouds. They begin with the initial conditions of a massive cloud of 10⁵ M_⊙ with a 10 pc radius for an average particle density of 380 cm⁻³. They simulate turbulent velocities of 3–10 km s⁻¹ and winds of up to 300 km s⁻¹ and include gravity. In the case of the shocked sphere, they employ a uniform density distribution for the cloud. The ISM, which comprises the wind, has a density of ρ ∼ 1 m_p cm⁻³ and a temperature of T = 10⁵ or 10⁷ K. They conclude that the pressure confinement of the shock wave increases the SFR and efficiency, i.e., that a higher fraction of the molecular gas is converted into stars. They also find that the pressure confinement decreases the sizes of the resulting star clusters, and that cloud rotation and shear against the ISM are not important unless the shear velocity is substantially greater than the sound speed of the constraining ISM. Their results are consistent with positive AGN feedback in gas-rich clouds with increased pressure (Gaibler et al. 2012; Bieri et al. 2016).

We summarize the parameters of these studies in Table 1.

We seek to improve past research on wind-cloud interactions by focusing on the high-velocity winds and shocks characteristic of AGN jet cocoons and outflow bubbles propagating through the host galaxy's ISM on a cloud-sized scale. We use BE spheres (Ebert 1955; Bonnor 1956) in place of uniform-density spherical clouds, and our spheres have higher peak densities, which are more representative of dense clumps within the ISM. We employ realistically high velocities in our winds to emulate those found from detailed models and simulations of AGN jets and winds. We include self-gravity, spatially correlated velocity perturbations, and a model of cooling and heating processes for all gas phases involved, in contrast to previously employed simplifications.

3.1. Numerics and Setup

3.2. Bonnor–Ebert Spheres

BE spheres are hydrostatic solutions of self-gravitating isothermal clouds embedded in a low-density environment. In particular, their nonuniform density profiles and pressure balance with the ambient medium make BE spheres a more realistic description for dense clumps of molecular gas than uniform spheres, which are not stable gravitationally. Barnard 68 provides an example of a hydrostatic blob of gas in pressure equilibrium with its environment and a density profile similar to that of a BE sphere (Alves et al. 2001). Because of their density profile and pressure balance, BE spheres are commonly employed as models to investigate low-mass star formation on small scales. For low ratios (<14.1) of the sphere's central density to its outer-edge density, BE spheres are stable against small perturbations. To assign values that would match those of a dense clump, we want to have a central and outer density a few orders of magnitude greater than the ambient medium. To make collapse more difficult, we select a ratio of central density to outer density to maximize the stability of the sphere, a value of 3. We have the ambient medium initially at rest and with a density of 1 m_p cm⁻³, comparable to the wind densities of many of the simulations. The BE sphere has an outer-edge density of 300 m_p cm⁻³ and a central density of 900 m_p cm⁻³, greater than those in other simulations (see Table 1). We deem these high densities necessary to simulate the high-density clumps within the ISM. We select a temperature of 19 K, which results in a sound speed of approximately 0.4 km s⁻¹. This choice of parameters results in a radius of 1.15 pc, an initial freefall time of approximately 3 Myr, and a cloud mass of 72 M_⊙.

We derive our density profile below, beginning with the equation of state, equation of gravitational stability, and spherically symmetric Poisson equation of gravity:

$$\begin{eqnarray}&&P=\displaystyle \frac{\rho {k}_{{\rm{B}}}T}{\mu {m}_{p}}=\rho {c}_{s}^{2},\end{eqnarray} \tag{ 3 }$$

 

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dr}}=-\displaystyle \frac{d{\rm{\Phi }}}{{dr}}\rho (r),\end{eqnarray} \tag{ 4 }$$

 

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}({r}^{2}{\rm{\Phi }})=4\pi G\rho .\end{eqnarray} \tag{ 5 }$$

We can combine the equation of hydrostatic equilibrium and the ideal gas law before substituting it into the above equation:

$$\begin{eqnarray}&&-{c}_{s}^{2}\displaystyle \frac{d}{{dr}}\mathrm{ln}\rho =\displaystyle \frac{d{\rm{\Phi }}}{{dr}}.\end{eqnarray} \tag{ 6 }$$

After substitution back into the spherically symmetric Poisson equation, we get a form of the Lane–Emden equation, a second-order ordinary differential equation (ODE) that essentially defines a BE sphere:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\displaystyle \frac{\partial \mathrm{ln}\rho }{\partial r}\right)=4\pi G\rho /{c}_{s}^{2}.\end{eqnarray} \tag{ 7 }$$

To compute the density profile, we define a central density, ρ₀, and then integrate numerically to the outer edge. We also apply the boundary condition that Φ(r = 0) = 0. BE spheres are well-discussed in the literature; for more information, see Kaminski et al. (2014) and references therein.

3.3. Parameter Space

We explore a parameter space that encompasses winds and shocks from starbursts to AGN jet cocoons and quasar outflows/winds. Most important, our wind velocities range from 300 to 3000 km s⁻¹ to properly simulate the impact AGN feedback could have on a BE sphere. Most studies in the literature do not explore winds that exceed velocities of 1200 km s⁻¹ despite galaxy-scale simulations and observations that show wind velocities greater than 1000 km s⁻¹, even when AGNs are in quasar mode as opposed to radio mode. While AGN jets and their cocoons may typically occupy the lower-density and higher-velocity regime, quasar winds are expected to be more mass-loaded.

For a spherically symmetric galaxy, the energy-conserving analytical models of a jet or outflow-driven bubbles depend only on the power of the jet or outflow (Wagner et al. 2012). The equations for the bubble's radius and resulting wind velocity are

$$\begin{eqnarray}&&{R}_{{\rm{b}}}={{At}}^{3/5},\end{eqnarray} \tag{ 8 }$$

 

$$\begin{eqnarray}&&{v}_{{\rm{w}}}=(3/5){{At}}^{-2/5,}\end{eqnarray} \tag{ 9 }$$

 where

$$\begin{eqnarray}&&A={\left(\displaystyle \frac{125{P}_{\mathrm{jet}}}{384\pi {\rho }_{{\rm{a}}}}\right)}^{1/5},\end{eqnarray} \tag{ 10 }$$

in which P_jet is the power of the jet or outflow and ρ_a is the density of the ambient medium. Using Equation (9), we can calculate the bubble velocity and radius on full-galaxy scales to demonstrate the range and locations of wind velocities from AGN feedback. We also vary the density of the ambient medium from 0.5 to 10 m_p cm⁻³ to further demonstrate the range of resulting values. We see that, in the first few Myr with that ambient density range, the bubble will have a radius less than 5 kpc and a velocity in the thousands of km s⁻¹. After 20 Myr, the bubble will have a radius between 7 and 14 kpc and a velocity in the hundreds of km s⁻¹. Therefore, to understand the variable impact of AGN feedback, one must explore the full parameter space of AGN winds and the unique combination of high velocities in the thousands of km s⁻¹ and densities resulting from black hole feedback.

Also unlike previous simulations in the literature, our winds are highly pressurized, meant to mimic the highly pressurized environments inside AGN jet cocoons and outflow bubbles. To be consistent with both theoretical and observational values for AGN winds, we simulate values within the ranges listed in Gaibler et al. (2012), Wagner et al. (2012, 2013) and Liu et al. (2013). We list the parameters for each of our simulations in Table 2. The run labels contain the density and velocity of the shock, so "d1-v3000" indicates wind density and velocity of 1 m_p cm⁻³ and 3000 km s⁻¹, respectively. We simulate winds with the combinations of densities of 0.5, 1, 3, and 10 m_p cm⁻³ with velocities of 300, 1000, and 3000 km s⁻¹. It should be noted that our model of a static BE sphere models star formation in low-mass or isolated sites rather than big cloud complexes or GMCs, where turbulence within the cloud plays a crucial stabilizing role. While for isolated or low-mass star formation sites, the AGN wind parameters are probably not affected much, massive molecular cloud complexes should strongly modify the propagation and properties of the AGN wind and may considerably shield the star-forming cores from the destructive processes. This scenario, however, is beyond the setup used in the present study.

Table 2.  Employed Values for the Simulations

Run	ρ	${v}_{{\rm{w}}}$ a	vw/cb	${P}_{\mathrm{ram}}$ c	Pd	Te	χf	${t}_{\mathrm{cc}}$ g	${t}_{\mathrm{KH}}$ h	${ \mathcal M }$ i
Label	(mp cm−3)	(km s−1)		(dyn cm−2)	(dyn cm−2)	(K)		(kyr)	(kyr)
d0-v300	0.5	300	0.001	7.52e–10	1e–09	8.96e+06	6e+02	95.84	166.09	1.07e+03
d1-v300	1.0	300	0.001	1.50e–09	1e–09	4.48e+06	3e+02	67.77	117.51	1.07e+03
d3-v300	3.0	300	0.001	4.51e–09	1e–09	1.49e+06	1e+02	39.13	68.00	1.07e+03
d0-v1000	0.5	1000	0.0033	8.35e–09	1e–09	8.96e+06	6e+02	28.75	49.83	3.57e+03
d10-v300	10.0	300	0.001	1.50e–08	1e–09	4.48e+05	3e+01	64.29	112.59	1.07e+03
d1-v1000	1.0	1000	0.0033	1.67e–08	1e–09	4.48e+06	3e+02	20.33	35.25	3.57e+03
d1-v1000-uni	1.0	1000	0.0033	1.67e–08	1e–09	4.48e+06	3e+02	20.33	24.54	3.57e+03
d3-v1000	3.0	1000	0.0033	5.01e–08	1e–09	1.49e+06	1e+02	11.74	20.40	3.57e+03
d0-v3000	0.5	3000	0.01	7.52e–08	1e–09	8.96e+06	6e+02	9.58	16.61	1.07e+04
d1-v3000	1.0	3000	0.01	1.50e–07	1e–09	4.48e+06	3e+02	6.78	11.75	1.07e+04
d10-v1000	10.0	1000	0.0033	1.67e–07	1e–09	4.48e+05	3e+01	6.43	11.26	3.57e+03
d3-v3000	3.0	3000	0.01	4.51e–07	1e–09	1.49e+06	1e+02	3.91	6.80	1.07e+04

Notes. Run d1-v1000-uni: uniform sphere with the same mass and radius as the BE sphere.

^aWind velocity. ^bWind velocity/c. ^cRam pressure of wind. ^dWind pressure. ^eWind temperature. ^f $\chi ={\rho }_{{\rm{c}}}/{\rho }_{{\rm{w}}}$. ^gCloud-crushing timescale. ^hKH timescale. ⁱMach number with respect to the sound speed of the sphere.

Download table as:  ASCIITypeset image

The lower-density value of 0.5 m_p cm⁻³ is included particularly because thermal instability of the dense gas starts at 1 m_p cm⁻³. Run d1-v1000-uni used a uniform sphere with the same mass and radius as the BE sphere. Additionally, the table includes the cloud crushing time and the KH timescale.

We can also derive the acceleration and therefore the velocity of a spherical cloud by a bubble-driven wind,

$$\begin{eqnarray}&&{m}_{{\rm{c}}}\displaystyle \frac{{{dv}}_{{\rm{c}}}}{{dt}}={C}_{D}{\rho }_{{\rm{w}}}{v}_{{\rm{w}}}^{2}\times \pi {R}_{{\rm{c}}}^{2},\end{eqnarray} \tag{ 11 }$$

 in which C_D is the drag coefficient, ${\rho }_{{\rm{w}}}$ is the density of the bubble-driven wind, and R_c is the radius of the cloud. We assume C_D ∼ 1. Employing Equation (11), we can calculate the theoretical cloud velocity as a function of time using the values for the winds we list in Table 2 and the average values for the BE sphere. This calculation results in expected velocities from 10 to 1000 km s⁻¹ on timescales less than 100 kyr. Yet this strongly depends on the (much-changing) cloud parameters.

3.4. Heating and Cooling

Although the multiphase nature of the ISM with a huge range of densities and temperatures is well-established, hydrodynamic simulations on galaxy scales often only describe the warm and hot phases properly. The cooler parts of the ISM are difficult to model and are widely mimicked by imposing a minimum pressure according to a polytropic equation of state or a minimum temperature. And, while resolution often makes this necessary—in particular with self-gravity, where the Jeans length needs to be well resolved to avoid numerical artifacts (Truelove et al. 1997)—this also means that the actual star-forming phase is never addressed and modeled explicitly, and all star formation–related results bear considerable uncertainty.

Since our study focuses on the pc scale, we necessarily need a thermodynamical model of the gas spanning the full range of temperatures, from ∼10 K in cold "molecular" clouds to >10⁷ K in the AGN wind. We handle this by adding a source term for heating and cooling to the energy equation with

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=\left(\displaystyle \frac{\rho }{{m}_{{\rm{p}}}}\right){\rm{\Gamma }}\,-{\left(\displaystyle \frac{\rho }{{m}_{{\rm{p}}}}\right)}^{2}{\rm{\Lambda }}(T/\mu ),\end{eqnarray} \tag{ 12 }$$

a heating parameter ${\rm{\Gamma }}=2\times {10}^{-26}$ erg s⁻¹, and a cooling function ${\rm{\Lambda }}(T/\mu )$ for the entire temperature regime (Figure 1, where μ is the mean molecular weight in m_p at the respective temperature: μ ≈ 2.3 for molecular gas and 0.62 for fully ionized gas). Numerically, the heating and cooling are computed with a linearized implicit method and substepping within the hydrodynamic time step. The high-temperature domain (T > 10⁴ K) is modeled after the collisional ionization equilibrium cooling function of Sutherland & Dopita (1993) with a solar metallicity, and the low-temperature domain (T < 10⁴ K) is modeled after Koyama & Inutsuka (2002) with their corrected fitting formula as given in Vázquez-Semadeni et al. (2007). For the transition between the two regimes, the high-temperature term of Koyama & Inutsuka (2002) is damped to continuously connect to and not exceed the Sutherland & Dopita (1993) function. This has the additional benefit of avoiding the overestimated Lyα cooling of the former and should bring us closer to the behavior found in detailed calculations of nonequilibrium chemistry and thermal balance (Micic et al. 2013) while still being able to stick to simple one-fluid hydrodynamics. The Koyama & Inutsuka (2002) fit includes many low-temperature processes, has often been employed for converging flow studies without full chemistry, and yields a two-phase behavior for the cold and warm neutral gas on the thermal equilibrium curve. Thermal instability can drive equilibrium gas of 1–240 m_p cm⁻³ (see Figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The cooling time,

$$\begin{eqnarray}&&{\tau }_{\mathrm{cool}}=\displaystyle \frac{3\,{k}_{{\rm{B}}}\,{m}_{{\rm{p}}}\,T}{2\,\mu \,\rho \,{\rm{\Lambda }}(T/\mu )},\end{eqnarray} \tag{ 13 }$$

for the cloud is initially 20 kyr at the outer edge; it decreases with rising density toward the center and decreases strongly during compression and collapse. Gas that is pushed out of thermal equilibrium (e.g., by compression) thus cools on short timescales. For comparison, we calculate the cloud ablation time,

$$\begin{eqnarray}&&{t}_{{\rm{abl}}}={R}_{{\rm{C}}}/{v}_{{\rm{w}}},\end{eqnarray} \tag{ 14 }$$

where R_C is the initial cloud radius, and find it to be in the range of 4 (for the 300 km s⁻¹ winds) to 0.4 kyr (3000 km s⁻¹). The corresponding cooling lengths for hydrodynamical processes within the cloud,

$$\begin{eqnarray}&&{L}_{\mathrm{cool}}={c}_{{\rm{s}}}\,{\tau }_{\mathrm{cool}},\end{eqnarray} \tag{ 15 }$$

 are initially 0.008 pc for the outer edge of the cloud and become much smaller after compression and collapse.

While ideally the cooling length should be resolved, this is simply impossible for our setup and the maximum resolution of 0.01 pc (cf. the detailed study on this by Yirak et al. 2010). When the cooling shells are unresolved, cooling at the neighboring cells is overestimated, but it is underestimated for the shell itself. As a result, the cloud could be more susceptible to KH and Raleigh–Jeans instabilities that ablate the clump. The tendency for gravitational collapse might be weaker than if the cooling shells were properly resolved. However, with additional resolution, other physical factors, such as thermal conduction or magnetic fields, become important and drastically complicate the required physics (see also Glover & Mac Low 2007). Predicting the impact of this lack of resolution is difficult and requires a much more sophisticated numerical model. However, these effects will mostly occur near the "surface" of the cloud and have less impact on the cloud cores and gravitational collapse.

3.5. Spatially Correlated Perturbations

While AGN blast waves and winds globally exhibit an approximately spherically symmetric expansion, they are not well-ordered flows and show considerable perturbation on smaller scales. Overpressured, expanding jet cocoons of AGN jets show the velocity components of the propagation of the cocoons and their leading bow shocks through a multiphase ISM (Figure 9 of Wagner et al. 2012), as well as velocity components of strong turbulence (Figure 15 of Gaibler et al. 2009). The two together generate stochastic flow perturbations on the scales that are relevant for interaction of such winds with cold clouds. To mimic this effect, we use the spatially correlated perturbations generated by a superposition of comoving plane waves on top of a constant inflow velocity ${{\boldsymbol{v}}}_{{\rm{w}}}$,

$$\begin{eqnarray}&&{\boldsymbol{v}}={{\boldsymbol{v}}}_{{\rm{w}}}+\displaystyle \sum _{\lambda }\,\displaystyle \sum _{i}^{{N}_{i}}{a}_{i}{{\boldsymbol{n}}}_{i}\sin ({k}_{\lambda }{{\boldsymbol{n}}}_{i}\cdot {{\boldsymbol{x}}}^{\prime }+{\phi }_{i}),\end{eqnarray} \tag{ 16 }$$

where k_λ represents the three discrete wave numbers considered (corresponding to wavelengths of 1/3, 1, and 3 pc); a_i, ${{\boldsymbol{n}}}_{i}$, and ${\phi }_{i}$ are the amplitude, random direction, and random phase of the N_i waves at each wave number, respectively; and ${{\boldsymbol{x}}}^{\prime }={\boldsymbol{x}}-{{\boldsymbol{v}}}_{{\rm{w}}}t$ is the location vector in the frame comoving with the AGN wind. The wave numbers are chosen to generate perturbations that are smaller than, similar to, and larger than the sphere itself, so that interactions on various scales can be excited. We choose N_i = 10 plane waves at each wave number, plane-wave amplitudes scaling $\propto {k}_{\lambda }^{-1/3}$, and the overall rms amplitude of the velocity perturbations normalized to 20% of the wind velocity ${{\boldsymbol{v}}}_{{\rm{w}}}$. Perturbations are only applied to the velocity variable on the boundary but clearly propagate into the other flow variables inside the computational domain. We note that these perturbations also avoid artifacts from an artificially symmetric setup (e.g., during the compression of the cloud). With respect to the interaction of the wind with the cloud, the choice of wave numbers excites perturbations on scales smaller, similar to, and greater than the truncation radius of the BE sphere.

3.6. Clump Analysis

To better investigate star formation, we employ a clump finder to look at the individual dense clumps that form after the BE sphere is shocked, rather than the sphere as a whole. This approach allows us to identify smaller pockets of gravitation instability and collapse that we would not see if we looked at the entire sphere as a single body, giving us a better idea of the trend of star formation. We use the RAMSES Clump Finder (Bleuler & Teyssier 2014) to find and track dense clumps resulting from the wind and collapse of the BE sphere in the various simulations. Typically, clump finders employ either a cell-based method, in which clumps and/or sink particles are formed based on conditions within an individual cell, or peak-based conditions, in which clumps and/or sink particles are based on the hydrodynamic state surrounding the density peaks. This new RAMSES Clump Finder uses a more sophisticated implementation of the peak-based method but with restrictions to exclude small density fluctuations, thus providing more realistic separations into individual clumps.

This algorithm works in five main steps. First, all cells with a density above a predetermined threshold are identified. Second, local density maxima are identified and given a peak identification number, and all the other cells are associated with one of these peaks based on steepest ascent. Third, the saddle points between the previously identified peaks are identified and determined to be significant or not based on the predetermined peak-to-saddle ratio. Fourth, all insignificant peak areas are merged with the closest significant ones, and the process is repeated until all the identified peaks are significant. Finally, the merger history for each clump is recorded and each cell within the clump is assigned the clump identification number for further investigation. For more detail, see Bleuler & Teyssier (2014). We modify the standard clump finder so that instead of analyzing data on individual clumps we can do so on each cell within a clump. Cell-by-cell analysis is particularly important for calculating the binding energy of each clump and subsequent collapse.

To probe gravitational instability in the resulting clumps, we experiment with two physical parameters taken as inputs for the clump finder: the density threshold and the peak-to-saddle ratio. The two values that best allow for the detection of gravitational instability are a threshold of 10^5.5 m_p cm⁻³ and a peak-to-saddle ratio of 1.2. These values allow for a focus on the ultradense center of the resulting clumps, precisely where collapse occurs and stars are formed. We define our gravitationally unstable regions as those that are Jeans unstable and cannot support themselves with their own internal pressure and that have binding energy in excess of their internal kinetic energies,

$$\begin{eqnarray}&&{M}_{{\rm{J}}}=5\langle R{\rangle }_{m}\langle P{\rangle }_{m}/(G\langle \rho {\rangle }_{m})\lt {M}_{\mathrm{clump}},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{i}\,{\phi }_{i}{m}_{i}\gt \displaystyle \sum _{i}\,{m}_{i}{v}_{{\mathrm{int}}_{i}}^{2}/2,\end{eqnarray} \tag{ 18 }$$

where ϕ is the gravitational potential, and we use the mass-weighted averages of radius, pressure, and density to calculate the Jeans mass. We calculate the resulting SFR,

$$\begin{eqnarray}&&\mathrm{SFR}=\displaystyle \sum _{\mathrm{clumps}}{m}_{\mathrm{clump}}/{t}_{\mathrm{ff}},\end{eqnarray} \tag{ 19 }$$

where t_ff is the local freefall time of the clump, calculated with the mass-weighted average density. Only clumps with masses >0.1 M_⊙ are included in the SFR calculation. We calculate the SFR with smaller, separate clumps rather than with larger, aggregated clumps to better probe gravitational instability. Because these are the most central, dense, and unstable parts of the clumps resulting from the wind and collapse, we calculate the SFR with complete efficiency. We integrate this SFR over the duration that the clumps are gravitationally unstable for the final stellar mass created, which we tabulate in Table 2. We do not remove mass from the simulation as time passes and stars are formed. As a result, the SFRs we calculate may just be indicative of a trend of star formation.

In all the simulations, we find that the wind drives shock waves into the sphere, followed by a very strong, general compression due to ram and thermal pressure of the wind and post-shock region. However, depending on the parameters of the simulation, specifically the ram pressure, we find that the cloud can be ablated before it has the chance to collapse. The simulations without collapse have KH times of about 20 kyr or shorter, meaning that the surface instabilities grow quickly and the cloud is stripped before gravity takes over. In the other simulations, however, the opposite is true: cloud density increases dramatically through the initial compression, followed by gravitational collapse before KH instabilities or the perturbations in the wind can ablate the sphere.

In all of the simulations with ram pressures smaller than that in d3-v1000, with the exception of d10-v300, the wind compresses the sphere to form dense filaments in the direction of the incoming wind. Figure 3 shows the column densities of all 12 simulations in ascending order of ram pressure. In the simulations with dense filaments, we see these filaments accelerated across the computational domain faster with increasing ram pressure. The highest-density parts of these filaments are toward the back end of the filament, away from the incoming wind and shielded by the slightly less dense front part of the filament. This is important for subsequent star formation because it prevents the shock from driving internal (small-scale) perturbations into the filament and clumps that would prevent gravitational collapse. The exception simulation, d10-v300, has an incoming wind density of 10 m_p cm⁻³; the spatially correlated perturbations cause this initial density to grow to over 100 m_p cm⁻³, and these dense perturbations ablate the sphere rapidly.

Zoom In Zoom Out Reset image size

Figure 4 shows the density projections for the four simulations with the highest ram pressures with smaller time intervals to resolve the rapid evolution of the sphere under those circumstances. Run d0-v3000 evolves in a manner similar to that of d0-v1000 and d1-v1000, forming a dense, peanut-shaped filament. By contrast, runs d1-v3000, d10-v1000, and d3-v3000 all show rapid cloud ablation without the formation of dense filaments stretched out in the direction of the wind velocity. In our analysis of star formation, we see that these dense filamentary structures are the result of gravitational collapse outpacing cloud ablation and are critical to forming stars.

Zoom In Zoom Out Reset image size

Generally, as the ram pressure and momenta of the wind increase, the cloud crushing time and the time that elapses until maximum density is achieved decreases. Figure 5 is a mass-weighted phase plot of density versus time for all 12 simulations. In the six simulations with some of the higher ram pressures, a vertical black line indicates the time when approximately 20% of the dense gas has been forced out of the computational domain. Runs d10-v300, d1-v1000-uni, d3-v1000, d10-v1000, d1-v3000, and d3-v3000 all show this effect. Runs d0-v1000 and d1-v1000, both with velocities of 1000 km s⁻¹, take 60 and 70 kyr to reach peak densities of 4 and nearly 5 orders of magnitude greater than the initial peak density, respectively. All the simulations with velocities of 3000 km s⁻¹—runs d0-v3000, d1-v3000, and d3-v3000—take 15, 25, and 35 kyr, respectively, to reach maximum densities roughly 3.5 orders of magnitude greater than the initial maximum density. However, runs d0-v300, d1-v300, and d3-v300, all of which have a wind velocity of 300 km s⁻¹, take roughly 150, 160, and 170 kyr, respectively, to reach peak density, climbing 4–5 orders of magnitude. In these three simulations with wind velocities of 300 km s⁻¹, as the ram pressure and momenta of the incoming winds increase, the time to peak density is slightly decreased, likely because of a small delay to the point where gravitational collapse can take over. Runs d0-v300, d1-v300, and d3-v300 show by far the most gravitational collapse of any of the simulations.

Zoom In Zoom Out Reset image size

The two simulations with wind densities of 10 m_p cm⁻³, runs d10-v300 and d10-v1000, reach peak densities of only 2 orders of magnitude greater than the initial value, and they do so quickly relative to the other simulations with the same velocity. This is because the spatially correlated perturbations with densities starting at 10 m_p cm⁻³ can quickly grow to 100–150 m_p cm⁻³ with gravity, and these very dense perturbations have an extremely destructive impact, quickly ablating the cloud before any collapse can occur.

Despite the similar density and morphological evolutions of the uniform sphere in run d1-v1000-uni and the BE sphere of the same mass in run d1-v1000, both with wind velocities of 1000 km s⁻¹ and densities of 1 m_p cm⁻³, the two runs have considerably different star formation, demonstrating the importance of the gas distribution within the cloud. The uniform sphere is already in an unstable state. We believe that the different density structures, specifically higher densities at the outer edges, better shield the inner core, allowing for earlier star formation.

The dense gas inside the sphere is accelerated to a variety of velocities depending on the simulation parameters, but, in all the runs, the entire sphere is initially accelerated to a velocity of at least 10 km s⁻¹. Figure 6 shows the evolution of the velocities of all gas with density greater than 150 m_p cm⁻³, basically limiting these plots to the gas in the sphere, dense filaments, and gas that has just been stripped. In each row, looking left to right, the velocity of the wind is increased by a factor of roughly 3, and the velocity of the sphere is likewise increased by a factor of roughly 3. In run d0-v3000, more than 50% of the dense gas has a velocity of ∼100 km s⁻¹ very early on. All the runs with a ram pressure greater than or equal to that in d0-v1000 show continuous or late acceleration of the dense gas.

Zoom In Zoom Out Reset image size

The differences between runs d1-v1000 and d3-v1000 and between runs d1-v3000 and d3-v3000 indicate how impactful the increase to a wind density of 3 m_p cm⁻³ is; all of the percentile lines show that the gas reaches much higher velocities, as much as 50 km s⁻¹ or more, far more quickly for the high wind density runs.

In runs d10-v1000 and d10-v300, the wind density is 10 m_p cm⁻³, and this density combined with the spatially correlated velocity perturbations in gravity results in lumps of gas with densities of around 100 m_p cm⁻³ traveling at high velocities. The high-density, high-velocity gas in the wind is represented by the horizontal line at a velocity 300 km s⁻¹ in the phase plot for run d10-v300 and the top left feature of velocity ∼1000 km s⁻¹ in run d10-v1000. All of this gas has a density of at least 150 m_p cm⁻³. This high-velocity, high-density gas in the wind quickly ablates the sphere, preventing it from reaching high densities or velocities, as well as preventing any collapse.

In runs d0-v300, d1-v300, and d0-v1000, the lowest ram pressure runs, the wind compresses the front of the sphere inward far more slowly than in the other simulations, allowing the gravitational potential to subsequently pull the back side of the sphere inward. Collapse-induced dynamics overcome wind-induced dynamics as the gas from the back side of the sphere collides with gas from the front side, slowing this initially high-velocity gas.

Here, Equation (11) provides an interesting point of comparison. Without ablation or collapse, uniform spheres would theoretically be accelerated to velocities similar to what we see in our simulations, with some differences. For example, in run d0-v300, with the lowest ram pressure, the analytical model and our simulations both show velocities of one to tens of km s⁻¹. However, in the analytical model, the velocity continually increases. In the simulation, the median velocity is initially high as the sphere is compressed, then it decreases as gravitational collapse takes over and the densities increase by several orders of magnitude. The same happens in d0-v1000 and d1-v300. On the higher end of the velocity spectrum, with wind velocities of 3000 km s⁻¹, the cloud never achieves the velocities in the simulations that the analytical models predict, not even the fastest 10% of the cloud. The median velocities of the clouds in runs d0-v3000, d1-v3000, and d3-v3000 are around 100 km s⁻¹, whereas the analytical models predict velocities closer to 1000 km s⁻¹. The reason for this difference is that, besides accelerating it, the wind also strongly changes the cloud's properties, particularly by decreasing its effective cross section by compressing it. Additionally, the cloud ablation time is much shorter than the time it would take to reach these velocities (∼100 kyr). When the values for the cloud radius in Equation (11) are decreased by factors of 2–10 and the density increases by 2–4 orders of magnitude as the cloud collapses into dense filaments, the expected values for the cloud velocity are likewise decreased to values under 100 km s⁻¹. Overall, the simulations show that the fluid dynamics of nonuniform spheres experiencing turbulent winds behave differently than simple estimates may predict because the cloud properties evolve with time.

As the dense filaments form and the gas is accelerated across the box, gravitational collapse can occur on scales as large as a clump and as small as our maximum resolution. Figure 7 shows the mass-weighted phase plot of the Jeans length of the individual cells versus time for each of the 12 simulations. We overplot a red line at 0.01 pc representing the maximum resolution of our simulation. With AMR codes, this also represents the maximum resolution of gravity, meaning that once the smallest cells obtain a particular density, the gravity within each cell can no longer be resolved. This means that gravitational collapse could occur inside the smallest cell without the code being able to resolve it. In turn, any gravitational collapse we calculate on a larger scale is therefore a lower limit on gravitational collapse. With better spatial resolution, we could likely detect more collapse, deeper potentials, and more subsequent star formation. Since the collapsing cores would continue to accrete material beyond the end of the simulation, the masses found in the stars would be expected to be larger than the masses found in our analysis. Yet feedback from the newly formed (proto-)stars may again limit the accretion of additional material. To address these issues, a much more detailed model of star formation than that attempted for the present study would be necessary. We note that a region of gas with the same mass and freefall time as our BE sphere but a different and unstable density distribution would, by arguments of freefall time, form roughly 5 M_⊙  over 250 kyr. Since star formation is generally found to show efficiencies well below unity, stellar masses beyond 2 M_⊙  (see Table 3) might already be understood as a sign of enhanced star formation.

Zoom In Zoom Out Reset image size

Table 3.  Stars Formed

Run	ρ	${v}_{{\rm{w}}}$ a	${P}_{\mathrm{ram}}$ b	${ \mathcal M }$ c	${t}_{\mathrm{cc}}$ d	${t}_{\mathrm{KH}}$ e	SFf	$\bar{v}$ g
Label	(mp cm−3)	(km s−1)	(dyn cm−2)		(kyr)	(kyr)	M⊙	(km s−1)
d0-v300	0.5	300	7.5e–10	730	96	166	37.4	3.2
d1-v300	1.0	300	1.5e–09	730	68	118	47.9	5.1
d3-v300	3.0	300	4.5e–09	730	39	68	29.4	17.4
d0-v1000	0.5	1000	8.4e–09	2400	29	50	3.9	9.4
d10-v300	10.0	300	1.5e–08	730	21	38	0.0	⋯
d1-v1000	1.0	1000	1.7e–08	2400	20	35	4.0	33.3
d1-v1000-uni	1.0	1000	1.7e–08	2400	20	25	5.7	14.1
d3-v1000	3.0	1000	5.0e–08	2400	12	20	0.1	⋯
d0-v3000	0.5	3000	7.5e–08	7300	10	17	0.0	⋯
d1-v3000	1.0	3000	1.5e–07	7300	7	12	0.0	⋯
d10-v1000	10.0	1000	1.7e–07	2400	6	11	0.0	⋯
d3-v3000	3.0	3000	4.5e–07	7300	4	7	0.0	⋯

Notes. Run d1-v1000-uni: uniform sphere with the same mass and radius as the BE sphere.

^aWind velocity. ^bRam pressure of wind. ^cMach number with respect to the sound speed of the sphere. ^dCloud-crushing timescale. ^eKH timescale. ^fCumulative stellar mass formed. Note that the stellar masses reported in this table are likely underestimated because real protostars would continue accreting mass from their envelopes beyond the durations of these simulations. ^gMass-weighted average velocity of star-forming clumps at the end of the simulation.

Download table as:  ASCIITypeset image

Figure 7 shows how much mass exists in the maximally resolved cells that are too dense to resolve the Jeans length. For the most part, this occurs because of the compression phase of each simulation. Runs d0-v1000, d1-v1000, d0-v3000, and d1-v3000 show this effect the most. In these simulations, after the clumps reach their peak densities, pushback from growing internal pressures overcomes the gravitational potential and moves the gas out of the regime where it would collapse inside a single cell. This figure shows, however, that for the most part, the potential within each cell is reasonably well resolved.

4.1. Clump Analysis

To investigate gravitational collapse and subsequent star formation, we evaluate the gravitational stability or instability of the individual clumps that form in the dense filaments. To best probe gravitational instability, we employ a threshold density of 10^5.5 m_p cm⁻³ and a peak-to-saddle ratio of 1.2 to identify the clumps most likely to collapse. However, to show that these are indeed the dense central cores, we also look at clumps with threshold densities of 10^4.5 m_p cm⁻³ so that we can properly analyze the more diffuse envelopes. Figure 8 shows the evolution of total clump mass versus time with the lower threshold of 10^4.5 m_p cm⁻³ for all 12 simulations. Figure 9 shows the evolution of total clump mass versus time with the higher threshold of 10^5.5 m_p cm⁻³. The discrepancy in mass but similarity in evolution of the two plots shows that the clumps with the higher density threshold are truly the dense cores of more massive, larger clumps and thus are shielded from the turbulent gas outside the larger clump.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows that only those simulations with velocities of 300 and 1000 km s⁻¹ and densities of 0.5, 1, and 3 m_p cm⁻³ form the dense clumps where stars form. The green lines show the sums of the mass of all of the clumps with masses >0.1 M_⊙, the blue lines show the total mass of the Jeans-unstable clumps, and the red lines show the total mass of the gravitationally bound and Jeans-unstable gas that will form stars.

None of the simulations with a wind velocity of 3000 km s⁻¹ form stars. None of the simulations with a wind density of 10 m_p cm⁻³ form stars because the velocity perturbations can cause the wind density to climb to up to 150 m_p cm⁻³, rapidly ablating the sphere. The simulation in which we have a uniform sphere rather than a BE sphere forms stars earlier as a result of the dense edges of the sphere, indicating that the distribution of gas does in fact make a difference in potential star formation.

The black vertical line in six of the simulations marks the time when approximately 20% of the dense gas (ρ > 150 m_p cm⁻³) has been forced out the far end of the box. These lines show that, in the simulations with the highest ram pressure, the clump analysis is valid because the clumps are stripped and destroyed before 20% of the dense gas has been forced out. Therefore, we expect that no star formation would occur in these simulations, even if we were able to track the simulations longer.

Furthermore, we see that nearly all of the gas in the dense clumps is Jeans-unstable, indicating that the internal (small-scale) velocities of the clump, rather than the thermal pressure, support it. This is significant because the wind and subsequent turbulence from the AGN feedback transfers a significant amount of kinetic energy to the clump, including to the internal velocities. However, as we focus on the higher density threshold clumps that comprise the cores of the lower density threshold clumps, we see the part of the clump that is shielded from the surrounding turbulence and lacks the high internal velocities to support the clump against gravitational collapse.

As the clumps collapse, we calculate the resulting cumulative stellar mass formed using the method described in Section 3.5. We find that, with increasing ram pressure, stars are formed earlier, and less stellar mass is formed. Figure 10 shows the cumulative stellar mass formed versus time in all the simulations in order of ram pressure. The simulations with wind velocities of 1000 km s⁻¹ begin forming stars earlier, around 50 kyr into the simulations, as the BE sphere is compressed to much higher densities. However, the SFR in these runs decreases with time because the clump masses decrease with time as they are stripped. The SFR increases with time in the simulations with wind velocities of 300 km s⁻¹ because the peak densities in these simulations tend to increase steadily through the simulation, and the resulting SFR mirrors the same pattern.

Zoom In Zoom Out Reset image size

An important result of the study is that we find that there is a decrease in the mass of stars formed with increased ram pressure in the wind, leading up to a threshold value of wind ram pressure of 1.67 × 10⁻⁸ dyn cm⁻². Table 3 shows the final mass of stars formed in each of the simulations. These simulations are ordered by the ram pressures of their respective winds. We note that d10-v300 does not form stars even though the listed ram pressure is lower than the threshold value; this is because perturbations cause the wind's peak densities to grow from 10 m_p cm⁻³ to over 150 m_p cm⁻³ in some places, creating ram pressures that are over an order of magnitude greater than the threshold.

Star clumps from the simulations with higher ram pressure clearly have higher velocities of up to 30 km s⁻¹, also shown in Table 3. These velocities are consistent with what Equation (11) predicts (tens of km s⁻¹). These results are also consistent with the positive radial velocity boost of stars formed in galaxies with AGNs from theoretical work by Gaibler et al. (2012), Silk et al. (2012), and Dugan et al. (2014); however, we find smaller velocities here than in those studies. We believe that a sphere larger than the one we simulate would survive longer and therefore have more time to be accelerated to higher velocities by the shock, reaching velocities closer to those calculated in the aforementioned studies.

These simulations cover the broad range of physical parameters for the winds caused by AGN feedback from jets or radio-quiet quasars. These winds can drive the peak densities of the cloud by up to 5 orders of magnitude through the formation of dense filaments and clumps, subsequently resulting in a drastic shortening of the freefall time. We find a threshold value for the ram pressure of the incoming wind below which the resulting clumps' binding energy overcomes the internal kinetic energy and pressure, resulting in gravitational collapse and star formation (shown in Table 3). The increased ram pressure of the incoming wind decreases the total mass of the resulting stars but increases the velocities of formed stars, allowing some to achieve velocities of up to 35 km s⁻¹ within the simulated time.

The results require some interpretation to be applicable to the debate between positive and negative feedback. We have taken a dense cloud of gas that would not collapse by itself and applied realistic physical circumstances expected for AGN feedback that induce star formation in half of the simulations. In this sense, our results are complementary to those of previous studies that shock uniform spheres that would collapse on their own, preventing the spheres from collapsing and showing negative feedback. However, nearly all of those do not include gravity and shock uniform spheres. In our suite of simulations, we include gravity and simulate BE spheres and find that the distribution of gas within the sphere can make a considerable difference in the results.

Our simulations show the formation of dense filaments reported in Mellema et al. (2002), Orlando et al. (2005), Cooper et al. (2009), Pittard et al. (2010), and Zubovas et al. (2014). Our results are in agreement with the latter three studies in that the geometry of the cloud impacts the results. Only Zubovas et al. (2014) includes gravity, and we both report star formation. However, we build on their results by identifying the velocity, density, and ram pressure range in which star formation is possible. We also include substantially higher velocities characteristic of AGN jet cocoons and outflow bubbles. With respect to the three phases of uniform cloud shocking introduced by Klein et al. (1994) and identified in Mellema et al. (2002), among others, Cooper et al. (2009) includes a fourth phase in which the cloud is ablated and destroyed. We find another possible phase in which gravitational collapse outpaces ablation.

Comparisons between cloud collapse time and ablation time have been made previously. Wagner & Bicknell (2011) argue that, for clouds less than 50 pc in size, the collapse time would be longer, resulting in negative feedback. However, Mellema et al. (2002) and Gaibler et al. (2012) speculate that only the outer, more diffuse envelopes of these clouds might be ablated, leaving behind the dense inner core to continue collapsing. Our results corroborate this idea. We find that a stratified density distribution of the sphere may change the manner in which the cloud is stripped, the mass that is lost, and how the remaining core is shielded and collapses, as in Nakamura et al. (2006). With a uniform sphere, shielding from the sphere's dense edge prevents the wind from injecting kinetic energy into the sphere, which leaves thermal pressure as the main mechanism of support. This is considerably different from the other simulations, in which the internal velocity dispersions are crucial to the clump's support. It is also worth noting that the extreme ram pressures resulting from AGN feedback quickly change the nature of the cloud being shocked, increasing the density and ablation time, decreasing the freefall time, and causing star formation. In the longer term, however, after the outer envelopes have been stripped and the cores have collapsed to form stars, the remaining gas budget for continued star formation may be depleted, causing long-term negative feedback after a burst of positive feedback. Strong fragmentation of the cloud may result in cloudlets too small to collapse or form stellar objects.

Within an entire galaxy, as a jet cocoon or outflow bubble expands through it, some clouds will be struck with velocities of 3000 km s⁻¹, while others, shielded by other clouds, will be struck with slower velocities of 1000 or even 300 km s⁻¹, depending on the location. Wagner et al. (2016) discusses the idea that positive and negative feedback can occur simultaneously within the same galaxy, depending on the size of the clouds being struck with shocks, and that the larger clouds will collapse under the pressure and form stars while the smaller clumps will be blown away. Since we do not systematically change the cloud properties in our study, we cannot assess the impact of variations in the typical sizes of clouds or cloud complexes, e.g., with redshift (Gaibler 2014). We only vary the winds striking the clouds, but we too find that simultaneous positive and negative feedback is possible in the same galaxy, provided that the winds depend on time and location within the galaxy. We also include gravity and have spatial resolutions that are much better than those achievable in galaxy-wide simulations.

We see that the star formation efficiency decreases with ram pressure. This is an interesting result that fits as a data point on the spectrum between negative and positive feedback: that AGN feedback can trigger star formation, but at a possibly lower efficiency than would otherwise occur, depending on the parameters of the feedback. There is, however, another intricacy that might be relevant: if star formation is not entirely prevented by high ram pressures, the timescale for collapse and the formation of stars within a cloud should be much shorter and happen more synchronously, rather than be spread over the (long) freefall time of the uncompressed cloud or cloud complex. Correspondingly, the impact of (negative) stellar feedback within the cloud by other young stars should be considerably smaller and thereby increase the efficiency of star formation.

Among other interesting conclusions, Cooper et al. (2009) find that the wind destroys regions of lower density faster than those of higher density. Combined with our results, in which the density of the core clumps can climb up to 5 orders of magnitude in ∼50 kyr, this would indicate that the clumps would survive for long periods of time. Cooper et al. (2009) also find that radiative cooling increases the survival time of the clouds. The clouds in Cooper et al. (2009) survive up to 1 Myr, achieving velocities of 150–400 km s⁻¹, and it is reasonable to believe that our dense clumps could survive that long if they were not already collapsing. Those velocities are higher than our clump velocities at the simulation's end, but that is to be expected because our simulations have maximum times of 250 kyr and our clumps are much denser.

The results from our simulations are similar to those of our previous studies, Gaibler et al. (2012) and Dugan et al. (2014) which show both that feedback in radio-loud galaxies may be positive and that stars formed during this period of feedback will have a positive velocity bias. Specifically, these stars are typically formed with greater radial and vertical (off the disk) velocity distributions. While we find velocities of up to 35 km s⁻¹ within the simulated time span for our small cloud, we note that larger clouds might survive longer and would be accelerated to higher velocities due to the larger effective cross section and thus form stars with higher peculiar velocities.

In turn, this also provides further computational support to a previously existing theory linking AGN feedback with hypervelocity stars (Silk et al. 2012). In this scenario, the cocoon or bow shock of a jet accelerates a pocket of gas to high velocities while compressing it to form stars. Silk et al. (2012) make analytical arguments based on energy and momentum, but here we provide computational evidence to the same end. Again, because our simulations are limited in the size of the cloud and the time we have to accelerate the cloud, we do not achieve stellar velocities as high as some of those observed in the Milky Way. However, what we show is that the phenomenon is possible on smaller scales.

It would be interesting in the future to extend the simulations to larger gas complexes, such as GMCs. For these, the assumption of having BE spheres is no longer suitable, and a turbulent GMC would need to be included. This may address the idea that the larger the cloud is, the more positive the feedback could be.

Our simulations show that the unique combination and range of ram and thermal pressure that AGN feedback provides can either trigger or prevent star formation in dense clumps within the ISM, depending mostly on the ram pressure. Our setup includes a broad range of parameters expected for winds originating from AGN jet cocoons or quasar winds. In the simulations that form stars, the wind compresses the cloud in the high-pressure post-shock region, which is then followed by the simultaneous stripping and collapse of the sphere. The initial compression from the shock and not gravity causes the peak density in the sphere to climb by as much as 5 orders of magnitude in as little as 50 kyr, reaching regimes where gravitational collapse overcomes internal pressure and kinetic energy. During this phase, the freefall time of the cloud drops from 3 Myr initially to as low as 10 kyr in subsequent clumps.

We analyze the clumps resulting from this collapse, as well as subsequent star formation within the clumps. We find indications of both positive and negative feedback, depending on the parameters of the wind. We show that the stratified density profile of BE spheres can have a considerable impact on gravitational collapse and star formation. We find that an increase in ram pressure and Mach number of the incoming wind decreases the total mass of stars formed but increases the velocities of the stars. Furthermore, star formation happens on a much shorter timescale and more synchronously. We identify threshold values for ram pressure of ∼2 × 10⁻⁸ dyn cm⁻², above which star formation is not expected to occur because the high internal velocities generated within the clumps support them against gravitational collapse long enough for the cloud to be destroyed.

Our results indicate that simultaneous positive and negative feedback will be possible in a single galaxy, as AGN wind parameters will vary with location within a galaxy. Similar phenomena are observed in AGNs (Cresci et al. 2015a, 2015b; Carniani et al. 2016).

Z.D. was supported by a Centre for Cosmological Studies Balzan Fellowship. V.G. was supported by the Sonderforschungsbereich SFB 881 "The Milky Way System" (subproject B4) of the German Research Foundation (DFG). The research of J.S. has been supported at IAP by the ERC project 267117 (DARK) hosted by Université Pierre et Marie Curie—Paris 6 and at JHU by NSF grant OIA-1124403. We would like to thank Alex Wagner for his input.