u/u/erica/BEPaper: BEsphere.3.tex

\documentclass[10pt,preprint]{aastex}
%\documentclass[12pt]{report}
\usepackage{natbib}
\usepackage{amsmath}
%\usepackage{fullpage}
\bibliographystyle{apj}

\begin{document}
\newcommand{\crc}{cells$/r_c$}
\newcommand{\rr}[1]{$R_{#1}$}


\shorttitle{The effect of a massive ambient medium on the collapse of a Bonnor-Ebert sphere: deviation from the canonical outside-in collapse}
\shortauthors{Kaminski et al.}

\title{Untitled - BE Crushing Solution} 
%\title{The effect of a massive ambient medium on the collapse of a Bonnor Ebert sphere: deviation from the canonical outside-in collapse}
%\title{A more physically accurate approach to the Bonnor Ebert sphere collapse problem: taking the ambient medium into account}
%\title{Taking the ambient medium into account when modeling the collapse of a BE sphere: a more physically motivated scenario}
%\title{Collapse of a Bonnor Ebert sphere in a massive ambient medium: deviations from the canonical outside-in collapse}
%\title{Considering the effects of a massive ambient medium on the collapse of a Bonnor Ebert sphere}
%\title{A new solution for the Bonnor Ebert collapse problem: the role of the ambient medium on modifying outside-in collapse}




\author{Erica Kaminski\altaffilmark{1}, Adam Frank\altaffilmark{1}}
\altaffiltext{1}{Department of Physics and Astronomy, University of Rochester, \
Rochester, NY 14620 \\Email contact: erica@pas.rochester.edu}


%{\it clumps} 

\begin{abstract}
A series of 3D hydrodynamic simulations were carried out to explore the effects of a massive ambient medium on the collapse of a marginally stable Bonnor-Ebert (BE) sphere. In particular, we sought to discern whether a massive ambient medium would be sufficient to drive the collapse of the sphere, whether this collapse would be triggered by the ram pressure of infalling material gravitationally accelerated by the BE sphere, and if this collapse would differ from the canonical {\it outside-in} collapse previous models have shown occurs for BE spheres. To these aims, various uniformly dense ambient environments were initialized to be in pressure equilibrium with the BE sphere. These ambient densities ranged from the conventional 'light' density ($\rho=0.01\rho(Rbe)$, where Rbe is the BE sphere's truncation radius), to a more physically realisitic ambient density equal to that of the sphere at its outer most edge ($\rho=\rho(Rbe)$). In each of the runs, ram pressure at the sphere's outer edge was found to exceed the critical threshold of external pressure on the BE sphere, although the resulting collapse varied. The collapse features were found to be sensitive to the density in the surrounding ambient medium. In the matched case, the sphere collapsed under a crushing wave. . . In contrast, the sphere in the light ambient medium remained dynamically stable, oscillating about its equilibrium values for $\sim$5 crossing times. These results reveal the importance of taking the ambient medium under consideration when modeling the collapse of a Bonnor-Ebert sphere. 
\end{abstract}

\keywords{a list to choose from}
\section{Introduction}

Various properties of the ``Bonnor-Ebert'' (BE) sphere \citep{bonnor1956, ebert1955}, a hydrostatic sphere in pressure equilibrium with its ambient environment, make it a good candidate for numerical modeling of protostellar collapse. First, as a candidate star forming structure is envisaged as gravitationally bound and unstable, it is easy to imagine a protostar evolving from an initially hydrostatic configuration. Indeed, spherical clumps have been observed in or near hydrostatic equilibrium, such as the Bok Globule B68 (check citation). Second, the stability criterion against gravitational collapse has been worked out analytically. Third, pushing the sphere out of the stability regime with various physical perturbations illuminate collapse characteristics. Such features of the collapse may help advance single star formation theory as well as provide clues to observational astronomers in identifying potential star forming sites.  

While the collapse of a BE sphere has been studied extensively over the years, the literature reveals studies of the BE sphere in precarious and unphysical situations; the BE sphere has largely been modeled as residing in artificially low density ambient mediums, seemingly to isolate the collapse of the sphere from the ambient environment \citep{foster1993, ogino1999, banerjee2004}. However, any actual collapsing cloud would not be discontinuous from the ambient medium. Additionally, ad hoc perturbations such as overall density increase within the simulation box, have been employed to force the sphere into collapse. In these ways, previous models have been physically unrealistic. Despite these certain drawbacks though, it is important to note that modeling the BE sphere in this way happened to reveal an important and unexpected feature of the collapse. 

The collapse of a BE sphere was originally anticipated to be an inside-out process, a characteristic feature of Frank Shu's similarity solution for the highly unstable singular isothermal sphere (SIS) \citep{shu1977}. In this well-known paper, \cite{shu1977} speculated that the collapse of {\it any} hydrostatic isothermal sphere, including Bonnor-Ebert spheres, would approach the SIS through a subsonic adjustment to a $1/r^2$ density distribution. The classic {\it inside-out} collapse of the SIS was then, as Shu proposed, a general feature of collapse, applicable to any hydrostatic isothermal sphere, unstable or not.   

Simulations of stable, flat-topped BE spheres in low density ambient environments, however, proved the contrary. \cite{foster1993} explored the collapse of BE spheres of varying truncation radii embedded in ambient backgrounds of uniform $\rho=0.01\rho(Rbe)$, where Rbe is the BE sphere's truncation radius. Their set ups showed that despite initial perturbation methods, as well as whether the sphere was initially in a stable or unstable hydrostatic regime, the collapse proceeded much differently than that of the SIS. Instead of an {\it inside-out} collapse, the collapse of a BE sphere was {\it outside-in}. Further, the collapse was not subsonic, rather supersonic. This directly refuted Shu's proposal. Studies of the collapse with more advanced fluid dynamic codes continued to provide support for outside-in collapse of the BE sphere {\citep[see e.g.][hereafter BPH04]{banerjee2004}. 

Now, although these early models illuminated unique and unexpected features of the BE collapse, they have been physically unrealistic given the discontinuous jump in density across the BE sphere/ambient boundary. More sophisticated simulatations have taken this under consideration. For instance, in calculating protostellar masses due to infall and dispersal of ambient gas, \cite{myers2008} modeled pre-stellar cores \citep{difrancesco2007, ward-thompson2007} as BE spheres seperated from a uniform ambient medium by a continuous, smooth density transition. \cite{hennebelle2003} looked at the effects of . 

Thus, there appeared to be a missing link in the studies of collapsing isothermal spheres between those embedded in highly unrealistic environments and those in more complex systems. This provided the motivation for the present work. We sought to study how the classic outside-in collapse of a BE sphere would be modified when it was placed in a non-negligible and continuous ambient medium. Further we did not artificially perturb the BE sphere into collapse, rather studied the role the ambient alone had in inducing collapse.

This paper is broken into the following sections. In Section 2 we will discuss our model and methods. In Section 3 we will present our findings. Namely, we found that the collapse resembled a crushing wave solution (Hannebelle, Whitmore, Di Francesco) early on in the presence of a massive ambient medium. As this wave propagated inward with time, the radial density and velocity profiles asymotically approached the outside-in results. In less massive ambient backgrounds, the sphere attempts to reequilibrate, but ultimately is pushed over the stability limit. Collapse resembles an intermediate state between the crushing wave and outside-in solution. In Section 4 we will present discussion and conclusions of these results.   


%"~" is $\sim$


% \cite{kleinwoods1998}

% \citep[see e.g.][]{truelove1998}


%$\propto n^2\Lambda(T)$

%\S~\ref{problem}


\section{Methods}\label{methods}
\subsection{Bonnor-Ebert sphere definitions}\label{BEdefs}
The equations that specify the parameters of the Bonnor-Ebert sphere are derived from the equations of hydrodynamics. Combining the condition for hydrostatic equilibrium in spherical coordinates, Poisson's equation for gravity, and the isothermal equation of state, 

\begin{equation}\label{HSE}
-\frac{\nabla P}{\rho} - \nabla \phi_{g} = 0
\end{equation}

\begin{equation}\label{Poisson}
-\nabla ^2 \phi_{g} = 4 \pi G \rho
\end{equation}

\begin{equation}\label{isoT}
P = \rho C_{s}^2
\end{equation}
where $\phi_{g}$ is the gravitational potential, P is the thermal pressure, ${\rho}$ is the density, and $C_{s}$ is the isothermal sound speed yields \citep[cf.][]{stahler2005},

\begin{equation}\label{rho(r)}
\rho(r) = \rho_{0} exp(-\phi_{g}/C_{s}^2)
\end{equation}

This is an equation for the radial density function of the hydrostatic isothermal sphere in terms of some central density, $\rho_{0}$, and the gravitational potential. The goal from here is to find $\phi_{g}$ so that the functional form of $\rho (r)$ is determined for a specified $\rho_{0}$. The following equation found by combining \ref{rho(r)} and \ref{Poisson} provides $\phi_{g}$,

\begin{equation}\label{phi}
\frac{1}{r^2} \frac{d}{dr}(r^2 \frac{d \phi_{g}}{dr}) = 4 \pi G \rho_{0} exp(-\phi_{g}/C_{s}^2)
\end{equation}
and hence the density profile of the BE sphere with central density $\rho_{0}$. It is common to cast \ref{phi} into a scale-free form by making the following variable substitutions, 

\begin{equation}\label{psi}
\psi = \frac{\phi_{g}}{C_{s}^2}
\end{equation}

\begin{equation}\label{xi}
\xi = (\frac{4 \pi G \rho_{0}}{C_{s}^2})^{1/2} r
\end{equation}
which leads to the famous {\it Lane-Emden} equation for an isothermal sphere, 

\begin{equation}\label{xi}
\frac{1}{\xi ^2} \frac{d}{d \xi}(\xi ^2 \frac{d \psi}{d \xi}) = exp(- \psi)
\end{equation}

There are two classes of solution to this equation. The first is the asymptotic solution, known as the {\it singular isothermal sphere} (SIS). It represents a sphere of infinite central density given by,

\begin{equation}\label{SIS}
\rho = \frac{C_{s}^2}{2 \pi G} r^{-2}
\end{equation}
This is the solution that was studied by Frank Shu in his analytical derivation of inside-out collapse \citep{shu1977}. Of greater physical relevance is the second solution, or {\it family} of  solutions, the non-singular Bonnor-Ebert spheres. These solutions adhere to the boundary conditions $\psi(0) = 0$ and $\psi'(0)=0$. Each  BE sphere is specified by a different truncation radius $R_{BE}$ and central density $\rho_{0}$, and is bound by an external pressure $P_{ext}$. They are characterized by a flat-topped density profile near their centers with outer densities that decline monotonically with radius (Fig.~\ref{fig_BE}). 


The stability of the BE sphere is determined by the value of $\xi$, also known as the BE sphere's {\it non-dimensional radius}. BE spheres are inherently unstable for configurations with $\xi > 6.451$, which is known as the critical radius, $\xi_{crit}$. The critical values of external pressure and radius of such a marginally stable BE sphere of mass $M_{BE}$ and with isothermal sound speed $C_{s}$ are given by (Spitzer),

\begin{equation}\label{Pcrit}
r_{crit} = 0.41 \frac{G M_{BE}}{C_{s}^2}
\end{equation}

\begin{equation}\label{Rcrit}
P_{crit} = 1.40 \frac{C_{s}^8}{G^3 M_{BE}^2}
\end{equation}

\begin{figure}[htbp]
\centering
\epsscale{.60}\plotone{Methods.eps}
\caption{{\small The density profile of a critical Bonnor-Ebert sphere as a function of $\xi$ in log-log space. The y-axis is in scaled units, normalized to the central density $\rho_{0}$ of the BE sphere. Given the scaled nature of this curve, it represents a family of solutions, each BE sphere given by a different $\rho_{0}$ and truncation radius.}}
\label{fig_BE}
\end{figure}

To initialize the BE sphere in our model we used an approximate analytic solution of the Lane-Emden equation \citep{liu1996}, which provided the density contrast of the BE sphere given its a) dimensional radius, b) non-dimensional radius ({\it$\xi$}), and c) central density \footnotemark[2].

\footnotetext[2]{documentation: https://clover.pas.rochester.edu/trac/astrobear/wiki/u/BonnorEbertModule}

\subsection{Model}
We modeled the collapse of the marginally stable BE sphere ($\xi = 6.5$) using AstroBEAR2.0\footnotemark[3], a highly parallelized, multidimensional, adaptive mesh refinement (AMR) code that solves the equations of hydrodynamics and magnetohydrodynamics \citep{andy2009,carroll2011}. AstroBEAR2.0 has a library of multiphysics tools, from heat conduction and resistivity, to radiative cooling and self-gravity. In the present work, AstroBEAR2.0 solved the Euler equations with a self-gravity source term. In conservative form, these equations are:

\footnotetext[3]{documentation: https://clover.pas.rochester.edu/trac/astrobear/}

(cite equations)


The self-gravity is handled by the elliptic solver, HYPRE\footnotemark[4], and is described in Appendix. 

\footnotetext[4]{documentation: https://computation.llnl.gov/casc/hypre/software.html}

The problem domain was a cube of side L $\approx 50$ pc, at the center of which we initialized the critical BE sphere. This BE sphere had the same properties as that used in BPH04. Namely, the radius of the sphere was r $\approx 1.62 pc$. The gas in the box was assumed to be atomic hydrogen ($\mu = 1$), and had a ratio of specific heats ($gamma = 1.0001$). Using the ideal gas equation of state, this set the simulation to be approximately isothermal. A ($gamma = 1.0$) could not be used, as to do so would introduce a division by zero in the solver. The central density was initialized to be $\rho_{0} = 2004 cm^-3$. This led to a temperature within the BE sphere of T $\approx$ 20 K.   


Given the symmetry of the problem it was possible to minimize computational costs by only simulating the 1st (positive) octant of the domain. That is, the simulated box actually extended from [0,25pc] in x-, y-, and z- directions, thereby encompassing only the 1st octant of the BE sphere centered at (0,0,0) (Figure 2a). Boundary conditions on the box were set to extrapolating (outflow) on all outside-facing sides, and reflecting on the faces that sliced through the BE sphere. The coarse grid was made up of $16^3$ cells, initialized with 5 levels of AMR within the BE sphere (Figure2b). Each level of AMR increased the effective number of cells in each direction by a factor of 2. Thus, at the start of the simulation the effective resolution was $\triangle x_{min} \approx 0.05$ pc, corresponding to $\approx 34 cells/R_{BE}$. Two additional levels of refinement were triggered during the BE sphere collapse when the ''Truelove condition'' \citep{truelove1997},

\begin{equation}\label{Truelove}
\lambda_{J} < 4 \triangle x
\end{equation}
was violated, where $\lambda_{J}$ is the local Jeans length associated with each cell center, and $\triangle x$ is the width of a cell on any given level. Essentially this is an AMR criterion that maintains adequate resolution of the Jeans length. Once the maximum refinement is reached, a sink particle can form inside one of the finest cells if a sequence of conditions are met that insure gravitational instability. These conditions are the same as those outlined in \cite{fedderrath2010}. The formation of a sink in those runs where collapse occurred determined the end of the simulation. By those points the density was not great enough to demand a switch to an adiabatic EOS (cite); an isothermal approximation was justified given the orders of density increase. Further, as this work was designed to track collapse properties of the BE sphere, the formation of a protostellar-like sink particle was a sufficient time to halt the evolution.

\subsection{Simulations}

We performed a set of 6 simulations, listed in Table~\ref{tab_runs}. In each run the marginally stable ($\xi = 6.5$) BE sphere was initialized to be in pressure equilibrium with an ambient medium of a different density. To test the stability of the BE sphere, we chose an ambient density of $\rho_{amb} = 0.01\rho(R_{BE})$, where $R_{BE}$ is the truncation radius of the BE sphere. This we named the ''sparse'' ambient case. We ran the simulation out for 5 crossing times $t_{sc}$, 

\begin{equation}
t_{sc} = r / C_{s} 
\end{equation}
and checked whether the sphere remained upright. Next we added a 10$\%$ overall density enhancement to the BE sphere in the light ambient case. This was the perturbation BPH04 and others used to initiate outside-in collapse of the BE sphere, and so we named it the ''classical'' case. Lastly, in the spirit of testing new conditions on the BE sphere, we did a series of runs that placed the sphere in ambient environments of varying density. Note these runs did not consist of any applied perturbations. These ambients ranged from the sparse ambient, to the denser ambient given by $\rho_{amb} = \rho(R_{BE})$. The latter run we named the ''matched'' case.  

\begin{table}[htbp]\centering\small
\label{tab_runs}
\begin{tabular}{rrc}
\hline\hline
%Run     Effective Resolution    
Run & Density of Ambient & Perturbation? \\
\hline
%$cells/r_c$ & $cells\ in\ z\times r$  & & $Approx.$\\
1 - Sparse  &    $\rho_{amb} = 0.01 \rho(R_{BE})$  &  No        \\
2 - Classic &    $\rho_{amb} = 1/100 \rho(R_{BE})$  &   Yes  \\
3a - Matched  &    $\rho_{amb} = \rho(R_{BE})$ & No      \\
3b  &    $\rho_{amb} \approx 1/3 \rho(R_{BE}) $  & No       \\
3c &  $\rho_{amb} = 1/10 \rho(R_{BE})$  &    No   \\
3d &  $\rho_{amb} \approx 1/31 \rho(R_{BE}) $  &   No       \\
\hline
\end{tabular}
\caption{\small Description of the 6 simulations of varying background densities. The first column lists the name of the run. The second column gives the ambient density of the simulation, and the third column states whether a perturbation was applied to initiate collapse. In only the ''classical'' run was a perturbation applied. This was a $10 \%$ overall density enhancement to the box, similar to \cite{banerjee2004}.}
\end{table}



\section{Results}\label{results}

\begin{figure}[htbp]
\centering
\epsscale{1}
\plottwo{rhoLight.eps}{vradLight.eps}
\caption{{\small Sparse case - Ambient is 1/100 $\rho (Rbe)$ - 5 crossing times}}
\label{light_case}
\end{figure}

\begin{figure}[htbp]
\centering
\plottwo{BPrho.eps}{vradBP.eps}
\caption{{\small Classic case - }}
\label{BP_case}
\end{figure}

\begin{figure}[htbp]
\centering
\plottwo{MatchedRho.eps}{MatchedVrad.eps}
\caption{{\small Matched Case}}
\label{Matched_case}
\end{figure}

\begin{figure}[htbp]
\centering
\plottwo{rhow3.eps}{vradw3.eps}
\caption{{\small Ambient is 1/3 $\rho(Rbe)$}}
\label{w3_case}
\end{figure}

\begin{figure}[htbp]
\centering
\plottwo{rhow10.eps}{vradw10.eps}
\caption{{\small Ambient is 1/10 $\rho(Rbe)$}}
\label{w10_case}
\end{figure}

\begin{figure}[htbp]
\centering
\plottwo{rhow31.eps}{vradw31.eps}
\caption{{\small  Ambient is 1/30 $\rho(Rbe)$}}
\label{w30_case}
\end{figure}




%The evolution of an adiabatic shocked clump progresses through several stages. The external shock will sweep across or ``pass over'' the clump radius in a time $t_{sp}=r_c/v_s$, where $r_c$ is th%e radius of the clump and $v_s$ the shock velocity. When the shock hits the clump, a new shock wave is transmitted into the clump. This transmitted shock crosses the clump in a time known as the %{\it cloud-crushing time} $t_{cc}$ (KMC94),

%\begin{equation}\label{eq_tcc}
%  t_{cc} = \frac{r_c}{v_{s,c}}\simeq\frac{\chi^{1/2}r_c}{v_s}
%\end{equation}

%\subsection{Radiative cooling}\label{prob_cooling}
% \cite{dysonwilliams1997, zeldovich2002}

%($10^5$--$10^7$ K)


%\begin{table}[htbp]\centering\small
%\begin{tabular}{rrrr}
%\hline\hline
% Run     Effective Resolution    
%Run & Effective Resolution & $\Delta x/L_{cool}$& $\Delta x/L_{cool,inspec.}$\\
%\hline
%$cells/r_c$ & $cells\ in\ z\times r$  & & $Approx.$\\
%12  &    $ 96\times    24$  &   -- & --       \\
%24  &    $192\times    48$  &    1 & $\sim 1$ \\
%48  &    $384\times    96$  &    2 & 1-2      \\
%96  &    $768\times   192$  &    4 &  3       \\
%192 &  $1,536\times   384$  &    7 &  5       \\
%384 &  $3,072\times   768$  &   15 &  7       \\
%768 &  $6,144\times 1,536$  &   29 & 11       \\
%1536& $12,288\times 3,072$  &   58 & 14       \\
%\hline
%\end{tabular}
%\caption{Resolutions of the 8 simulations, from 24 to 1,536 cells per clump radius $r_c$. The third column gives the number of cells per cooling length $L_{cool}$ from Eq.~\ref{eq_lcoolclump}, and the fourth column via visual inspection (see \S~\ref{analysis}). At \rr{12} the cooling length was unresolved.\label{tab_resolutions}}
%\end{table}


%\begin{figure}[htbp]
%\plotone{diagram_rotated.eps}
%\caption{{\small Important features of the flow are given in a synthetic Schlieren image of run \rr{1536} at $t\sim 0.5 t_{cc}$. The image has been reflected about the axis of symmetry, with the reflection showing the location of AMR refined regions. The bow shock wraps tightly around the clump from the strong cooling. The clump surface is ablated by its interaction with the postshock flow. A slip stream forms behind the clump. Transmitted shocks propagate internally through the clump. The external shock is susceptible to the cooling instability, and a conical reflected shock forms off of the axis of symmetry which reengages the flow.}}
%\label{fig_diagram}
%\end{figure}


%\begin{figure}[htbp]
%\centering
%\plotone{schlieren-montage-small.eps}

% Figure~\ref{fig_schlieren}

%\section{Discussion and Conclusion}\label{discussion}

%``Convergence'' has a clear meaning 



%\rr{100}--\rr{200}

%\acknowledgements


\bibliography{erica}
\end{document}