Testing the stability of supersonic ionized Bondi accretion flows with radiation hydrodynamics

TL;DR

This paper analyzes the stability of ionized Bondi accretion flows onto young stellar objects, revealing marginal stability and recurring collapse of HII regions through analytic models and 1D radiation hydrodynamics simulations.

Contribution

It provides a new analytic steady-state solution for ionized accretion flows and demonstrates their marginal stability and instability via simulations.

Findings

Steady-state two-temperature solution predicts compact HII regions.

Ionization self-consistently treated leads to marginal stability.

HII regions tend to recur and collapse over time.

Abstract

We investigate the general stability of 1D spherically symmetric ionized Bondi accretion onto a massive object in the specific context of accretion onto a young stellar object. We first derive a new analytic expression for a steady state two temperature solution that predicts the existence of compact and hypercompact HII regions. We then show that this solution is only marginally stable if ionization is treated self-consistently. This leads to a recurring collapse of the HII region over time. We derive a semi-analytic model to explain this instability, and test it using spatially converged 1D radiation hydrodynamical simulations. We discuss the implications of the 1D instability on 3D radiation hydrodynamics simulations of supersonic accreting flows.

Equations

$${\dot{M}}_a=4\pi r^2\rho (r)v(r),$$

$$\frac{v^2(r)}{2C_s^2}+\ln (\rho )-\frac{2R_B}{r}=\mathrm{constant},$$

$$r_B(r)=\frac{GM}{2c_s^2(r)}$$

$$\approx 110.9\left(\frac{M}{{\mathrm{M}}_{\odot }}\right){\left(\frac{C_s}{2{\mathrm{km}\mathrm{s}}^{-1}}\right)}^{-2}\mathrm{AU},$$

$$\frac{v^2(r)}{2C_s^2}+\ln (\rho )-\frac{2R_B}{r}=\ln (\rho {}_{\infty }),$$

$$\rho (r)=-\frac{\rho {}_B{R}_B^2{C}_s}{r^{2}v(r)}.$$

$${\left(\frac{v(r)}{C_s}\right)}^2-\ln \left[{\left(\frac{v(r)}{C_s}\right)}^2\right]-4\ln \left(\frac{r}{R_B}\right)+3-4\frac{R_B}{r}=0,$$

$$v(r)=-C_s\sqrt{-W\left(-{\left(\frac{R_B}{r}\right)}^4\exp \left(3-4\frac{R_B}{r}\right)\right)},$$

$$v(r)=\{\begin{array}{ll}-C_s\sqrt{-W_{-1}\left(-\omega (r)\right)}& r\le R_B,\\ -C_s\sqrt{-W_0\left(-\omega (r)\right)}& r>R_B,\end{array}$$

$$\omega (r)={\left(\frac{R_B}{r}\right)}^4\exp \left(3-4\frac{R_B}{r}\right).$$

$$\rho {}_i(R_I)v_i(R_I)$$

$$\rho {}_i(R_I)\left(v_i^2(R_I)+C_{s,i}^2\right)$$

$$\Gamma =\frac{\rho {}_i(R_I)}{\rho {}_n(R_I)}=\frac{v_n(R_I)}{v_i(R_I)},$$

$$C_{s,i}^2\Gamma {}^2-\left(v_n^2(R_I)+C_{s,n}^2\right)\Gamma +v_n^2(R_I)=0,$$

$$\begin{gathered} \Gamma =\frac{1}{2C_{s,i}^2}(v_n^2(R_I)+C_{s,n}^2 \\ \pm \sqrt{{\left(v_n^2(R_I)+C_{s,n}^2\right)}^2-4C_{s,i}^2{v}_n^2(R_I)}). \end{gathered}$$

$$v_R$$

$$v_D$$

$$\rho {}_i(r)=\frac{\rho {}_i(R_I)R_I^2{v}_i(R_I)}{r^2{v}_i(r)}.$$

$$\begin{gathered} \frac{v_i^2(r)}{2C_{s,i}^2}+\ln (\rho {}_i(r))-\frac{GM}{C_{s,i}^{2}r}= \\ \frac{v_i^2(R_I)}{2C_{s,i}^2}+\ln (\rho {}_i(R_I))-\frac{GM}{R_I{C}_{s,i}^2}, \end{gathered}$$

$$\begin{gathered} {\left(\frac{v_i(r)}{C_s}\right)}^2-\ln \left[{\left(\frac{v_i(r)}{C_s}\right)}^2\right]-4\ln \left(\frac{r}{R_I}\right)+\ln \left(\frac{v_i^2(R_I)}{C_{s,i}^2}\right) \\ -4\frac{R_{B,i}}{r}-\frac{v_i^2(R_I)}{C_{s,i}^2}+4\frac{R_{B,i}}{R_I}=0, \end{gathered}$$

$$v_i(r)=-C_{s,i}\sqrt{-W(-\omega {}_i(r))},$$

$$\begin{gathered} \omega {}_i(r)={\left(\frac{R_I}{r}\right)}^4{\left(\frac{v_i(R_I)}{C_{s,i}}\right)}^2 \\ \exp \left(4\frac{R_{B,i}}{R_I}-4\frac{R_{B,i}}{r}-\frac{v_i^2(R_I)}{C_{s,i}^2}\right). \end{gathered}$$

$$\rho (r)$$

$$v(r)$$

$$\rho {}_n(r)$$

$$\rho {}_i(r)$$

$$v_n(r)$$

$$v_i(r)$$

$$\omega {}_n(r)=$$

$$\omega {}_i(r)=$$

$$\exp \left(4\frac{R_{B,i}}{R_I}-4\frac{R_{B,i}}{r}-\frac{v_i^2(R_I)}{C_{s,i}^2}\right).$$

$$\rho {}_i(R_I)$$

$$v_i(R_I)$$

$$\Gamma =\frac{1}{2}\left(\frac{v_n^2(R_I)}{C_{s,i}^2}+\frac{1}{P_c}-\sqrt{{\left(\frac{v_n^2(R_I)}{C_{s,i}^2}+\frac{1}{P_c}\right)}^2-4\frac{v_n^2(R_I)}{C_{s,i}^2}}\right),$$

$$Q(t)=4\pi {\int }_{R_c}^{R_I(t)}{n}_H^2(r,t)\alpha {}_B(T)r^2\mathrm{d}r,$$

Full text

††pubyear: 2019††pagerange: Testing the stability of supersonic ionized Bondi accretion flows with radiation hydrodynamics–A

1 Introduction

The problem of a spherically symmetric accretion flow onto a massive star including a rigorous treatment of the hydrodynamics was first introduced by Bondi (1952), and is consequently referred to as the Bondi problem. Mestel (1954) studied the effect of ionizing radiation on the original Bondi problem, using an approximation whereby both the ionized and the neutral region are assumed to be isothermal, with different isothermal sound speeds and a sharp transition from ionized to neutral regime at the ionization front. Keto (2002, 2003) showed that this work explains the existence of a steady state solution in which an ionized region is trapped inside the accreting flow. This explains the existence of trapped compact and ultracompact Hii regions around massive accreting stars. These ultracompact Hii regions are indeed observed around massive protostars, with observed size estimates varying from $\sim{}10^{3}$ AU (Churchwell, 2002) to as small as $\sim{}10$ AU (Ilee et al., 2016). He also showed how accretion can continue while the protostar is already emitting UV radiation, which has implications for the final mass of a young stellar object.

Recently, this work was extended to 3D using a Monte Carlo radiation hydrodynamics code (Lund et al., submitted to MNRAS). They were able to numerically simulate a trapped Hii region, and also showed the change in regime from a trapped R type ionization front into an expanding D type ionization front under a changing source luminosity, with the subsequent shut down of accretion onto the central star. However, their work also casted doubt about the stability of the steady state solution of Keto (2002), as their trapped Hii regions showed signs of periodic collapse, whereby the ionization front periodically collapsed onto the central star and moved outwards again.

In this work, we will analyse this instability in more detail, using a high resolution 1D spherically symmetric radiation hydrodynamics code, and a semi-analytic stability analysis.

Note that massive stars are expected to form through an accretion disk (Beltrán & de Wit, 2016; Kuiper & Hosokawa, 2018), which leads to more extended ultracompact Hii regions, and stabilises the ionization front (Sartorio et al, submitted to MNRAS). This renders the assumption of spherical symmetry an unlikely scenario for accretion onto a protostar. Nevertheless, this work is still relevant for our understanding of trapped Hii regions and can serve as a useful benchmark test for models that combine hydrodynamics, photoionization and gravitational accretion.

2 Analytics

In this section, we present the analytical solution to the original Bondi problem, and the extended Bondi problem with ionizing radiation. This work is a summary of the original work of Bondi (1952) on steady state accretion, and extends the work of Mestel (1954) and Keto (2002) with a full analytic expression for two temperature steady state accretion. We will summarize the main equations and present the full analytic steady state solution, before addressing the stability of this analytic solution under self-consistent ionization.

2.1 Bondi accretion

We assume spherically symmetric flow onto a massive object of constant mass $M$ with a constant accretion rate (Bondi, 1952)

[TABLE]

with $\rho{}(r)$ and $v(r)$ the mass density and fluid velocity as a function of radius $r$.

If we furthermore neglect the gravitational force due to the mass of the accreting fluid and assume the potential is dominated by $M$, and assume an isothermal equation of state with a constant sound speed $c_{s}(r)=\sqrt{{\rm{}d}P(r)/{\rm{}d}\rho{}(r)}=C_{s}$, the constancy of the accretion rate leads to a specific form of Bernoulli’s equation:

[TABLE]

where we introduced the (constant) Bondi radius

[TABLE]

with $G=6.674\times{}10^{-11}~{}{\rm{}m}^{3}~{}{\rm{}kg}^{-1}~{}{\rm{}s}^{-2}$ Newton’s constant.

Assuming that the velocity vanishes at infinity, we find the value of the constant:

[TABLE]

with $\rho{}_{\infty}{}$ the rest density of the accreting material. By evaluating this equation at $r=R_{B}$, we can find the density $\rho{}_{B}$ at the Bondi radius: $\rho{}_{B}=e^{3/2}\rho{}_{\infty}{}\approx{}4.48\rho{}_{\infty}{}$.

Using (1), we find an expression for the density as a function of radius:

[TABLE]

Substituting this in (5), we find

[TABLE]

which has general solution (Cranmer, 2004)

[TABLE]

where $W(x)$ is the Lambert W function, and the negative root was chosen since we want material to be accreting onto the central object.

The Lambert W function is a complex valued function with an infinite number of branches, but two of these branches, $W_{0}(x)$ and $W_{-1}(x)$, are real valued in the range $[-{1}/{e},0]$, with $W_{0}\left(-{1}/{e}\right)=W_{-1}\left(-{1}/{e}\right)=-1$, $W_{0}(0)=0$, and $W_{-1}(x)\rightarrow{}-\infty{}$ for $x\rightarrow{}0$. Since we want the velocity to be a strictly increasing function of radius, we see that we start on the $W_{0}$ branch for $r>R_{B}$, and switch to the $W_{-1}$ branch at $r=R_{B}$. For $r\rightarrow{}0$, the velocity diverges.

We end up with the following analytic expressions for the velocity for steady state accretion onto a massive object:

[TABLE]

with

[TABLE]

2.2 Ionizing radiation

To include ionization from a central isotropic point source, we follow Mestel (1954) and assume a sharp transtion from ionized to neutral at a well defined radius $R_{I}$. Both inside and outside the ionized region the gas obeys an isothermal equation of state, with the constant sound speed given by $C_{s,i}$ and $C_{s,n}$ respectively. We will characterize the change in equation of state in terms of the pressure contrast $P_{c}={C_{s,i}^{2}}/{C_{s,n}^{2}}$. Since the temperature in the ionised region will be higher than that in the neutral region, we always have $C_{s,i}>C_{s,n}$, and hence $P_{c}>1$.

Since (6) and (9) are completely determined by the boundary condition at $r\rightarrow{}+\infty{}$, this solution is still valid for $r\geq{}R_{I}$, with $C_{s}=C_{s,n}$, and $R_{B}=R_{B,n}={GM}/(2C_{s,n}^{2})$. To find the solution for $r<R_{I}$, we impose conservation of mass and momentum across the ionization front to find the following Rankine-Hugoniot conditions:

[TABLE]

where $\rho{}_{i}(r)$ and $v_{i}(r)$, and $\rho{}_{n}(r)$ and $v_{n}(r)$ are the density and velocity in the ionized and neutral regions respectively.

If we introduce the jump factor

[TABLE]

we find the following equation for the jump as a function of the known neutral velocity $v_{n}(R_{I})$:

[TABLE]

with solutions

[TABLE]

$\Gamma{}$ only has real solutions for $|v_{n}(R_{I})|\geq{}v_{R}$ or $|v_{n}(R_{I})|\leq{}v_{D}$ (note that $v_{n}(r)$ is always negative because we are in an accretion regime), with

[TABLE]

which are commonly referred to as the R type and D type solutions. We will limit ourselves to R type solutions. In the region between the R type and D type solutions there is no solution with a single shock at the position of the ionization front, but more complex double front solutions could be possible; we do not study these in this work.

It is obvious that $|v_{n}(R_{I})|>C_{s,i}$, so that the R type solution implies supersonic neutral gas at the ionization front. Furthermore, $\Gamma{}>1$ for $|v_{n}(R_{I})|\geq{}v_{R}$, which implies $|v_{i}(R_{I})|<|v_{n}(R_{I})|$.

$\Gamma{}$ has two solutions in the R type region, corresponding to a large and a small jump (corresponding to the positive and negative sign in (15)). We call these the weak and strong R type solution. If we assume a density structure that initially evolved through Bondi accretion without ionization, then the solution will naturally evolve into the weak solution with the smallest jump, and this is the solution we will assume below. In this case, the velocity in the ionized region will also be supersonic (see Appendix A for a detailed analysis of weak and strong R type solutions and imposed velocity restrictions).

Inside the ionized region, equation (1) still gives us the density:

[TABLE]

The Bernoulli equation (2) now becomes

[TABLE]

since we now need to apply boundary conditions at $r=R_{I}$ instead of $r\rightarrow{}+\infty{}$.

Similarly as before, we can combine this equation with (6) to get

[TABLE]

where we introduced the Bondi radius for the ionized region $R_{B,i}={GM}/(2C_{s,i}^{2})$. This equation has the general solution

[TABLE]

with

[TABLE]

Note that as in the neutral case, this solution will switch from the $W_{0}$ to the $W_{-1}$ branch at a critical radius $R_{c}$, for which $v_{i}(R_{c})=-C_{s,i}$. However, we know that this cannot happen for a weak R front, as $v_{i}(r)<-C_{s,i}$. We hence only need to consider the $W_{-1}$ branch.

The full solution for spherically symmetric accretion with ionization (in the case $|v_{n}(R_{I})|\geq{}v_{R}$) is:

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

The parameters $\rho{}_{i}(R_{I})$ and $v_{i}(R_{I})$ are given by the jump conditions:

[TABLE]

with

[TABLE]

so that the entire solution is parametrised in terms of $P_{c}$, $R_{I}$, $\rho{}_{B,n}$, $C_{s,n}$ and $M$.

The density and velocity profile of the two temperature solution is plotted in Fig. 1 for three different values of the ionization radius $R_{I}$. It can be seen that higher values of $R_{I}$ lead to higher densities and lower velocities inside the ionised region; the solution in the neutral region is the same for all models.

2.3 Self-consistent ionization

The derivation above did not make any assumptions about how the ionization happens, and just assumed ionization leads to a constant ionization radius $R_{I}$. However, in general, the ionization radius $R_{I}(t)$ will vary with time, and the time evolution will be linked to the hydrodynamical quantities, $\rho{}(r,t)$ and $v(r,t)$. To find out if the steady-state solution obtained above is actually stable, we have to include this effect.

For simplicity, we will assume that the dynamical time scale $t_{\rm{}dyn}$ over which the hydrodynamical quantities evolve is much larger than the recombination time scale, $t_{\rm{}rec}=\left(n_{\rm{}H}\alpha{}_{B}(T)\right)^{-1}$, where $n_{\rm{}H}$ is the number density of hydrogen atoms and $\alpha{}_{B}(T)$ is the collisional recombination rate of ionized hydrogen to all excited levels of neutral hydrogen, but excluding ionizations to the ground state (the so-called on-the-spot approximation), which depends on the temperature $T$. Under this assumption, the ionization structure of the system instanteneously adapts to any change in density, and we can model ionization by simply changing the equation of state.

To get the ionization radius $R_{I}$, we use a Strömgren approximation, and assume that all material inside the ionization region is completely ionized. At the ionization radius, we have a strong discontinuity, and outside the ionization radius the material is completely neutral. In this case, the ionization balance equation is given by (Osterbrock & Ferland, 2006)

[TABLE]

where $Q(t)$ is the total ionizing luminosity of the star. Since the steady state density profile is very centrally peaked, and since we do not realistically expect the profile to keep up until $r=0$, we also introduced an inner cut off radius $R_{c}$, below which we assume no more ionizations to occur.

Under the two temperature approximation, the recombination rate $\alpha{}_{B}(T_{i})=\alpha{}_{i}$ is constant. If we further assume our gas to be composed of hydrogen only, we can rewrite the ionization balance equation as

[TABLE]

where $m_{\rm{}H}=1.674\times{}10^{-27}~{}{\rm{}kg}$ is the hydrogen atomic mass. Taking the time derivative of this equation, we find the following equation for the time evolution of $R_{I}(t)$:

[TABLE]

where we assumed the ionization radius is a continuous function of time, and used the general form of Leibniz’ rule:

[TABLE]

2.4 Stability analysis

For what follows, we will assume a constant source luminosity $Q(t)$, and assume we already have a steady state solution as derived in 2.2. We will now perturb this solution by introducing a small density perturbation outside the ionization radius. As shown by (37), this will not affect the ionization radius, and the perturbation will accrete until it reaches the ionization front radius. Once it reaches the ionization front and crosses it by a small distance $\varepsilon{}$, the ionization front radius will change according to

[TABLE]

since this is the only part of the integral that is non trivial. This expression does not lend itself to a rigorous linear perturbation analysis because of the presence of an integral in a differential equation, so that we will restrict ourselves to a semi-analytic analysis below.

What happens next depends on the sign of ${\partial{}\rho{}(r,t)}/{\partial{}t}$. If this sign is positive (a positive perturbation), the value of the integral will be positive and the ionization front will move inwards. For a negative perturbation, the ionization front will move outwards.

In either case, the perturbation does not cause a restoring force that balances out the perturbation and restores the initial ionization front radius, as the perturbation will be further accreted inside the ionized region, and will keep contributing to (39). Depending on the size of the perturbation, the ionization radius will either settle into a new dynamic equilibrium value, or will keep moving in the same direction. In the limit of a very large positive perturbation, the ionization front will move with the accreting bump and collapse entirely. This qualitative behaviour of the ionization front radius constitutes a first important result: the change in ionization front caused by a perturbation cannot be undone while that perturbation is still present inside the ionized region. The rate at which the ionization radius changes depends on the size and the shape of the perturbation.

Once the initial perturbation leaves the ionized region by crossing the cut off radius $R_{c}$, the situation is nominally reset to the original two temperature situation. However, since the profile behind the perturbation adjusted to a new ionization radius while the perturbation was still present, the amount of available material to ionize will now no longer match the initial, stable two temperature solution. For a positive perturbation, the ionization front will have shrunk and the density will now be too low; for a negative perturbation the reverse will have happened. We now end up in the opposite scenario of what we started out with: the initially positive perturbation now results in a profile with a negative perturbation, while the initially negative perturbation is now a positive perturbation. This is a second important result: once a perturbation leaves the ionized region, a new perturbation of opposite sign is created.

The new perturbation of opposite sign has to ensure the ionizing luminosity $Q(t)$ is exactly balanced by recombinations (see (36)). Since the new perturbation will be created at a larger radius than the original perturbation, and since the Bondi accretion profile for the density scales as $\sim{}r^{-1.4}$, a larger density perturbation is needed to achieve this. This leads to the final important result: the newly created perturbation will always be larger than the original perturbation.

In conclusion, any initially small perturbation will lead to an oscillation of the ionization front radius on a time scale of the free fall time (the time a perturbation created at the ionization front needs to reach the inner cut off radius) and with a growing amplitude, until the amplitude is large enough to cause the ionization front to move at the same speed as the perturbation. Physically, there are four interesting scenarios, depending on the sign and the size of the initial perturbation:

A large positive perturbation will cause a collapse of the ionization front whereby the ionization front gets trapped inside the perturbation and collapses onto the central cut off radius. 2. 2.

A small positive perturbation will initially cause a collapse of the ionization front, but the ionization front will move at a slower pace than the accretion velocity of the perturbation. When the perturbation reaches the central cut off radius, the ionization front will go through a phase of expansion, followed by one or more subsequent oscillations with growing ionization front speed, until the ionization front gets trapped inside a perturbation. 3. 3.

A large negative perturbation will initially cause an expansion of the ionization front until the perturbation reaches the central cut off radius. At this point the ionization front will collapse as in the case of a large positive perturbation. 4. 4.

A small negative perturbation will initially cause an expansion of the ionization front, followed by multiple oscillations with growing ionization front speed, until the ionizatio front gets trapped inside a perturbation.

In the next section, we will test our semi-analytic stability analysis by numerically seeding known perturbations of a given sign and size that correspond to these four scenarios.

3 Numerical simulations

3.1 1D code

3.1.1 Algorithm

We use a textbook (Toro, 2009) implementation of a 1D spherically symmetric second order finite volume method111Available from https://github.com/bwvdnbro/HydroCodeSpherical1D. We use a very conservative version of the slope limiter of Hopkins (2015). Fluxes are estimated using an HLLC Riemann solver (Toro, 2009) with a polytropic index $\gamma{}=1.001$ (we tested our results against an exact Riemann solver), and are limited using the per-face slope limiter of Hopkins (2015). After every time step the pressure is reset using the isothermal equation of state, effectively mimicking the use of an isothermal Riemann solver. The spherical source terms needed to make the method spherically symmetric are added using a second order Runge-Kutta method, as suggested by Toro (2009).

Gravity is added as an extra source term in the momentum equation, using an operator splitting method that uses a leapfrog scheme (Springel, 2010). We assume that the central mass is much larger than the total mass of the fluid, such that we only consider the mass of the central object, as an external force.

We know from Section 2 that the central density is divergent, which will cause our numerical method to break down for small radii. We therefore set up an inner boundary radius with outflow boundary conditions, inside which we do not follow the hydrodynamics. At a much larger radius, we set up an inflow boundary where we impose the constant accretion rate.

To include ionization, we assume the Strömgren approximation introduced above, so that the ionization radius is determined by finding the root of the following function:

[TABLE]

where $r_{j}$ is the central radius of cell $j$; $\rho{}_{j}$ is its density. $\Delta{}r$ is the radial size of a single cell, while

[TABLE]

is the ionizing luminosity that makes it out of the inner mask (we will treat $Q^{\prime}$ as a simulation parameter).

Once the ionization radius is found, we set the pressure of the region based on the following adapted isothermal equation of state:

[TABLE]

where $P_{c}$ is the pressure contrast between the ionized and neutral region that was introduced in Section 2.

3.1.2 Simulation parameters

We want to set up physically plausible initial conditions for which our assumptions hold, i.e. we want the density to be high enough to guarantee $t_{\rm{}rec}<<t_{\rm{}dyn}$, while at the same time having a density that is low enough to actually be ionized by a physically plausible input luminosity $Q^{\prime}$. For a central stellar mass of $18~{}{\rm{}M}_{\odot}{}$, we expect an ionizing luminosity $Q_{i}\sim{}10^{48}~{}{\rm{}s}^{-1}$ (Keto, 2003; Sternberg et al., 2003). For the steady state solution, the ionization balance equation is given by

[TABLE]

We have to numerically integrate this equation to find the ionizing luminosity that corresponds to a given set of two temperature Bondi parameters.

Fig. 2 shows the ionizing luminosity $Q_{i}$ as a function of the neutral Bondi density $\rho{}_{B,n}$, for a central mass $M=18~{}{\rm{}M}_{\odot}{}$, neutral sound speed $C_{s,n}=2.031~{}{\rm{}km~{}s}^{-1}$ (corresponding to a hydrogen only gas with $T_{n}=500~{}{\rm{}K}$), pressure contrast $P_{c}=32$ (corresponding to an ionized temperature $T_{i}=8000~{}{\rm{}K}$), and assuming an ionization front radius $R_{I}=30~{}{\rm{}AU}$ (similar to Lund et al., submitted to MNRAS) and inner cut off radius $R_{c}=10~{}{\rm{}AU}$ (due to numerical precision issues, we cannot compute the argument to the Lambert W function for radii $<\sim{}10~{}{\rm{}AU}$). Also shown is the corresponding recombination time scale at the ionization front, computed as $T_{d}={m_{\rm{}H}}/{(\rho{}_{n}(R_{I})\alpha{}_{i})}$, with $\alpha{}_{i}=4\times{}10^{-13}~{}{\rm{}cm}^{3}~{}{\rm{}s}^{-1}$. For realistic ionizing luminosities $\sim{}10^{48}~{}{\rm{}s}^{-1}$, we find a neutral Bondi density $\sim{}10^{-19}~{}{\rm{}g~{}cm}^{-3}$, and a corresponding recombination time scale $\sim{}10^{-2}~{}{\rm{}yr}$. At the ionization radius, the free fall time is $T_{ff}=\sqrt{{2R_{I}^{3}}/{(GM)}}=8.723~{}{\rm{}yr}$ (the free fall time taking into account the Bondi velocity as initial velocity is of the same order of magnitude), so the recombination time is more than a factor 100 smaller than the typical dynamical time in the entire ionized region.

We also tested this assumption using an expensive 1D time dependent Monte Carlo radiation transfer algorithm (Tootill, private comm.) and found an excellent agreement between the time dependent results and the results obtained by assuming instanteneous ionization.

In what follows, we will hence use a neutral Bondi density $\rho_{B,n}=10^{-19}~{}{\rm{}g~{}cm}^{-3}$. We will use a simulation box that spans from $R_{c}=10~{}{\rm{}AU}$ to $R_{m}=100~{}{\rm{}AU}$, subdivided into 2700 equal sized radial cells, such that our desired ionization radius is on a cell boundary.

3.1.3 Stable solution

As a first test of the 1D code, we will try to obtain the stable two temperature solution derived in Section 2, starting from a constant density and velocity, and setting a constant ionization front radius $R_{I}=30~{}{\rm{}AU}$. The analytic solution is entirely fixed by specifying the neutral Bondi parameters, ionization front radius and pressure contrast; similarly the numerical solution should be entirely determined by the inflow boundary conditions at $R_{m}$, the mass of the external point mass, and the ionization front position and pressure contrast if we do not restrict the solution at $R_{c}$ (which we can effectively do by using open boundaries). As initial conditions, we hence set $\rho{}(r)=\rho{}_{n}(R_{m})$ and $v(r)=v_{n}(R_{m})$ everywhere.

Fig. 3 shows the density and velocity at 4 different times during the simulation. By $t=20~{}{\rm{}yr}$, the simulation has settled into a steady state solution. Fig. 4 shows the relative difference between this steady state solution and the analytic steady state solution found in Section 2, at time $t=40~{}{\rm{}yr}$. Both solutions are in excellent agreement. We conclude that our 1D code is capable of resolving the necessary physics to study this problem.

3.1.4 Stability test

As a second test, we test the stability of the steady state solution when the ionization front radius is computed self-consistently. To this end, we use the steady state output of the previous test, and compute the ionization radius based on equation (40).

Fig. 5 shows the resulting density and velocity profiles at four different times during the simulation. Initially, the steady state solution is stable. However, after $\sim{}10~{}{\rm{}yr}$, the density around the ionization front starts to break up into small peaks that travel inwards with the accretion flow. These small peaks cause tiny fluctuations in the ionization front radius, which gradually grow. After $\sim{}20~{}{\rm{}yr}$, the ionization front radius is no longer stable, and starts to oscillate, as can been seen from Fig. 6.

It is important to note that we did not seed the small oscillations that lead to the collapse of the steady state solution; these oscillations are seeded numerically by accumulated round off error. We conclude that the steady state solution is only marginally stable, as even very small scale noise can cause it to become unstable.

3.1.5 Convergence

The result for the stability test discussed above depends on numerically seeded noise, which depends on resolution. This makes it impossible to perform a standard convergence test to check whether the time evolution of the ionization front corresponds to a physical effect or is simply caused by numerical issues. To show that our results do converge, we have to modify our set up so that we can obtain results that are independent of the spatial resolution for a sufficiently high spatial resolution. We identify two issues that need to be addressed: (a) we need to seed the initial collapse of the ionization front with a known perturbation that is the same for all resolutions, (b) we need to make sure that newly created instabilities that cause the recurring collapse of the ionization front are consistent between resolutions. Since these depend on the shape of the ionization front, we need to make sure we resolve the ionization front consistently between resolutions.

To obtain a spatially resolved ionization front, we convert the strong jump into a smooth linear transition. As before, we compute the ionization front radius $R_{I}$ from equation (40). However, now the neutral fraction $x_{\rm{}H}(r)$ of the gas does not just jump from [math] to $1$ across the ionization front, but is given by

[TABLE]

where $S$ is the maximal slope of the transition, and the transition width is $W={3}/{(2S)}$. The transition itself is a third order polynomial:

[TABLE]

This particular transition shape was chosen because it allows us to control the slope of the transition, while ensuring continuous first derivatives in the transitions points $r=R_{I}-W/2$ and $r=R_{I}+W/2$.

The adapted equation of state now becomes

[TABLE]

Fig. 7 shows the stable two temperature solution for a range of simulations with different values of the transition width $W$ and for different spatial resolutions. As expected, the transition region smooths out the jump from the neutral to the ionized region, but outside the transition region the density and velocity still follow the analytic solution.

Fig. 8 shows the ionization front radius as a function of time for the same simulations when we self-consistently update the ionization state. All simulations still become unstable after some time, with the instability still being seeded by numerical noise. Since the noise depends on the resolution, the instability is seeded differently for different resolutions, and we do not obtain a converged result.

To seed the instability, we use the output of the set up runs shown in Fig. 7 and add a small gaussian-like density perturbation filter $\phi{}(r)$ to the density, such that $\rho{}^{\prime}(r)=\phi{}(r)\rho{}(r)$:

[TABLE]

with $u={|r-r_{0}|}/{h}$, where $r_{0}$ is the center of the gaussian-like bump and $h$ its width. $\rho{}_{p}$ is the amplitude of the perturbation, and the shape of the perturbation is based on the cubic spline kernel commonly used in SPH codes (see e.g. Springel, 2005). We will use $r_{0}=65~{}{\rm{}AU}$ and $h=5~{}{\rm{}AU}$. We will choose 4 different amplitudes for the perturbation to test the various scenarios derived above: $\rho{}_{p}=0.1$, $\rho{}_{p}=2$, $\rho{}_{p}=-0.1$ and $\rho{}_{p}=-0.01$. These correspond respectively to a small positive seed that should cause an ionization front movement slower than the accretion velocity, a large positive seed that should cause the immediate collapse of the ionization front, a negative seed that should initially cause an expansion of the ionization front followed by a collapse, and a negative seed that should cause a number of oscillations.

Fig. 9 shows the time evolution of a simulation with a small positive seed perturbation. The perturbation initially moves inwards with the accretion flow until it reaches the ionization front. Once it arrives there, it causes an increase in the density inside the ionization region, which causes the ionization front radius to shrink according to equation (37). The density enhancement continues to move inwards with the accretion flow, and for a while the ionization front moves inwards as well. The region outside the ionization front radius in the meantime adjusts to the lower pressure and settles into a neutral Bondi accretion, with an overall lower density.

When the initial perturbation reaches the inner outflow boundary, the mass inside the perturbation is effectively lost from the system, causing the density inside the ionized region to go down again, and subsequently the ionization front radius to expand again. However, since the density inside the original ionization front radius partially settled into a neutral Bondi profile, the overall density inside the ionized region will be lower than in the steady state two temperature solution, and the new ionization front radius will expand beyond the original ionization front radius. This is not a steady state solution anymore, so the ionization front radius will consequently shrink again as the inner density adapts to the ionization again. During that process, a large peak in density builds up around the ionization front radius, which in combination with the ionization front evolution (equation (37) leads to a runaway collapse of the ionization front towards the mask. Once this peak reaches the inner outflow boundary, the situation is reset again, and the whole process repeats itself.

Fig. 10 shows the evolution of the ionization front radius as a function of time for the simulations with a negative seed perturbation with $\rho{}_{p}=-0.1$. We clearly see how the repeated runaway collapse of the ionization front leads to a periodic oscillation of the ionization front radius with a period of $\sim{}10~{}{\rm{}yr}$, roughly consistent with the free fall time of the system. It is clear that the combination of a smooth transition that damps small scale noise and the seeding of the instability with a well defined perturbation allow us to achieve numerical convergence for the high resolution simulations. This in principle should allow us to use this test as a benchmark for this type of instability. Note that the convergence is better for initial conditions with a larger transition region from fully ionized to fully neutral. For what follows, we will hence limit ourselves to simulations with an ionization transition of $5~{}{\rm{}AU}$ and a resolution of 2700 cells.

Note that the expanding ionization front at the end of each oscillation moves relatively fast, so that we cannot guarantee that the dynamical time for this expansion is significantly larger than the recombination time; a main assumption for our ionization treatment. However, the overall density inside the central region is low enough to guarantee the build up of a new ionization front at the expected radius in a time that is short enough not to significantly alter our dynamics, and we can hence still assume our results are qualitatively correct, albeit with a slightly longer oscillation period. We confirmed this using a time dependent Monte Carlo radiative transfer simulation. Similarly, moving the inner outflow boundary to smaller radii will lead to a small increase in oscillation period, but will not change the qualitative behaviour, as we still expect accreted material at small enough radii to be lost from our system of interest.

Fig. 11 shows the time evolution of the ionization radius for the 4 different scenarios discussed in 2.4. In the case of a positive density perturbation, the ionization front initially collapses. If the bump is small, the ionization front does not move at the same speed as the bump, and the collapse halts when the bump moves accross the inner outflow boundary. If the bump is large enough, the ionization front moves along with the bump and both collapse onto the inner outflow boundary.

For a negative density bump, the ionization front initially expands until the density bump is accreted accross the inner boundary. At this point, the expanded ionization region will start to contract again. For a large negative bump, this immediately leads to a collapse of the ionization front. For smaller bumps, the collapse is only partial and leads to an oscillation of the ionization front with an amplitude that increases over time.

In all cases, the initial perturbation eventually evolves into a periodic oscillation of the ionization front, i.e. all scenarios naturally evolve into the case where the ionization front movement becomes coupled to the movement of the density perturbation.

These simulations hence confirm the semi-analytic model we derived before.

3.1.6 TORUS comparison

To investigate the impact of our approximations on the result of the 1D models, we reran our 1D test with the radiation hydrodynamics code torus (Harries, 2000), using the “detailed” model as described in Haworth et al. (2015). This includes polychromatic radiation transport, diffuse field radiation (i.e. not the on-the-spot approximation) as well as helium, metals and the associated line cooling in the thermal balance. In these calculations a hydrodynamics-only calculation is run until a steady state is reached. The medium is the fully ionized, which lowers the optical depth and hence speeds up the convergence of the first photoionization equilibrium calculation. Subsequently, photoionization and hydrodynamics steps are performed iteratively to follow the RHD evolution. The ionizing luminosity of the source was tuned to approximately reproduce the same ionization front radius as in the models above.

The resulting time evolution of the ionization front is shown in Fig. 12. As in the models above, the ionization front periodically shrinks and then reappears. This shows that even with a more detailed treatment of the photoionization, perturbations of the two-temperature solution are enhanced by the spherical symmetry of the system and never dampen out. It is interesting to note that in this case, the ionisation front does not collapse all the way up to the mask. The reason for this is that the ionization front in this case is located upstream from the density peak that causes the collapse, so that the density peak is absorbed by the mask before the ionization front reaches the mask radius.

3.2 3D simulation

The presence of a 1D instability has important implications for 3D RHD simulations of spherically symmetric Bondi accretion. If spherical symmetry is not explicitly imposed on the radiation field or on the coupling between radiation field and hydrodynamics, then the instability will be seeded differently in different directions. Once initial asymetries have been seeded, the problem is no longer spherically symmetric and true 3D effects will govern the hydrodynamical evolution.

We explore these effects by running the same Bondi set up described above in a full 3D RHD simulation using the Monte Carlo RHD code CMacIonize (Vandenbroucke & Wood, 2018). Since we are not interested in the detailed results, we limit ourselves to low resolution simulations and do not explicitly check for convergence. More realistic simulations will be subject of future studies. We run two simulations with different grids used to discretize the fluid. Both simulations follow a fluid in a box of $100~{}{\rm{}AU}\times{}100~{}{\rm{}AU}\times{}100~{}{\rm{}AU}$ with inflow boundary conditions and a point mass (and ionization source) at the centre. Just as for the 1D simulations, a central sphere with radius $10~{}{\rm{}AU}$ is masked out. The mass, initial ionization radius, neutral Bondi density, neutral gas temperature and pressure contrast are set as before.

A first set up uses a regular Cartesian grid of $128\times{}128\times{}128$ cells. To compute the ionization front position, we use $10^{7}$ Monte Carlo photon packets and 10 iterations. Initial instabilities will be seeded by two mechanisms: Poisson noise inherent to the Monte Carlo method, and discretization error due to the low resolution grid used. The latter will have a strong dependence on specific grid directions (e.g. the coordinate axis directions and diagonals), and from the top panel of Fig. 13 it is clear that they constitute the dominant seeding mechanism, resulting in a very symmetrical, yet no longer spherically symmetric ionization region.

The higher pressure in the directions with a larger ionization radius will cause a tangential force that evacuates these regions and pushes material into directions that are still neutral, reenforcing the differences between both regions. While some directions with larger ionization radii still periodically collapse due to the radial instability, the directions with a smaller ionization radius actually remain stable. At later times, the ionization region has a chaotic, asymmetric shape, while the density field shows clear signs of grid symmetries.

To eliminate the effect of Cartesian grid symmetries, we run the same simulation using an unstructured Voronoi grid consisting of 100,000 cells, using $10^{6}$ photon packets for 10 iterations. This is the same grid structure used to run CMacIonize as a moving-mesh code, but without actually moving the mesh in between time steps (since this would be very complicated for a 3D problem with spherical accretion). In this case, the results look less artificial (see bottom panel of Fig. 13), but the qualitative interpretation is the same: the combination of the 1D instability and 3D effects will cause the ionization front to break up into a non-spherically symmetric structure, and the ionization front will stabilize in some directions.

We conclude that unless ionization is treated in a perfectly spherically symmetric manner, 3D simulations will never result in a perfectly symmetric ionization region, irrespective of any additional 3D effects that could be included in these simulations. This is valuable background knowledge for the interpretation of future 3D studies of Bondi accretion that will include rotation, and a general limitation for RHD codes in studying problems that involve both accretion and photoionization.

4 Conclusion

In this work, we studied the stability of spherically symmetric ionized Bondi accretion flows around young stellar objects. We derived an anlytic expression for a two temperature steady state solution in which an R type ionization front with a size of the order of $10~{}{\rm{}AU}$ is trapped in the highly supersonic accreting flow, consistent with the predictions of Keto (2002, 2003). We were able to accurately reproduce this analytic profile in a constrained 1D RHD simulation.

However, we also showed that this steady state solution is only marginally stable if ionization is treated self-consistently, leading to a collapse of the ionization front under small perturbations in density or ionizing luminosity. We were able to confirm our semi-analytic model of this instability using 1D RHD simulations, and were able to obtain spatially converged results by introducing a noise suppressing transition layer and a seeded instability.

Both the constrained simulations that lead to the two temperature steady state solution and the converged simulations of the instability could serve as benchmark problems for 1D RHD codes.

Finally, we showed that the existence of a 1D instability has profound implications for 3D simulations of ionized accretion flows when the radiation field is not treated in a perfectly spherically symmetric manner. The instability makes it impossible to retain spherical symmetry, and will amplify grid symmetries.

Acknowledgements

We want to thank the anonymous referee for constructive and insightful comments regarding the convergence of our simulations. BV and KW acknowledge support from STFC grant ST/M001296/1. NS would like to thank CAPES for graduate research funding. KL acknowledges support from the Carnegie Trust. DFG thanks the Brazilian agencies FAPESP (no. 2013/10559-3) and CNPq (no. 311128/2017-3) for financial support. TJH is funded by an Imperial College London Junior Research Fellowship.

Appendix A Analysis of strong and weak R type solutions

Here we perform a detailed analysis of equation (15) for the case of an R type front, with $v_{n}(R_{I})<-v_{R}$. We begin by rewriting the equation in terms of the dimensionless variables $x=-{v_{n}(R_{I})}/{C_{s,i}}$ and $A={C_{s,n}}/{C_{s,i}}$:

[TABLE]

Note that $A=\sqrt{{1}/{P_{c}}}<1$. To get a real solution, we require

[TABLE]

and hence the neutral region is supersonic. For the critical value $x_{c}=1+\sqrt{1-A^{2}}$,

[TABLE]

which implies $v_{i}(R_{I})=-C_{s,i}>v_{n}(R_{I})$. We now explore how the value of $\Gamma{}$ changes for $x>x_{c}$ in the case of a strong and a weak R type front.

The first derivatives of $\Gamma{}(x,A)$ are

[TABLE]

For $x>x_{c}$ and $0<A<1$, these derivatives have no real roots, which means $\Gamma{}(x,A)$ is monotonous in these regions. Furthermore, we have that ${\partial{}\Gamma{}(x,A)}/{\partial{}x}\rightarrow{}\pm{}\infty{}$ for $x\rightarrow{}x_{c}$, which means that $\Gamma{}>\Gamma{}_{c}$ for a strong R front, and $\Gamma{}<\Gamma{}_{c}$ for a weak front.

For the weak front, we find (note the signs):

[TABLE]

or

[TABLE]

which means that the ionized region is always supersonic. For a strong R type front the condition $\Gamma{}>\Gamma{}_{c}$ does not impose similar restrictions, so that a strong front can be subsonic or supersonic.

Note that there is some confusion in literature about strong and weak R type fronts and whether or not they are supersonic or subsonic. From our analysis, it follows that a weak R type front, defined as the front with the smaller jump in density and velocity across the ionization front, is always supersonic, and is the only front expected to exist in nature. The strong front, defined as the front with a larger jump in density and velocity across the ionization front, can be either supersonic or subsonic. For the critical case $v_{n}(R_{I})=-v_{R}$, the strong and weak front are the same, and the ionized velocity is the ionized sound speed.