numerical-starformation.tex

\subsection{Star formation and feedback}
\label{sec.ov.star}

\subsubsection{Overview}

Due to the computationally unfeasible number of stars in a galaxy
($10^{11}$) and the lack of detailed understanding of star formation, a
number of phenomenological star formation models are included in
\enzo.
Broadly speaking, these methods all work in similar
ways: at a specified time interval, all grid cells that are at the
highest local level of refinement (i.e., that have no child cells) are
examined to see if they meet a set of criteria for star formation.
This may simply be a baryon density threshold, but can also include more
complex tests, such as an examination of local cooling and dynamical
time scales, molecular hydrogen fraction, metallicity, and converging
gas flows.  If the star formation criteria are met, some mass of gas
is taken away from the cell in question and a ``star particle'' with
the same mass is placed in the center of that cell with the same
velocity as the removed gas.  This star particle is then
allowed to inject mass, momentum, thermal energy, metals, and possibly
magnetic fields and/or cosmic ray populations into its local
environment.  In general, the particle is treated as an ensemble of
stars, with feedback properties occuring over time according to the
assumed initial mass function of the stellar population.

In the following sections, we describe some of the more widely-used
star formation and feedback methods used in \enzo.  We note that
similar methods have been employed in many other codes used for galaxy
formation, with comparable implementations for both star formation and
feedback in other grid-based codes.  With regards to particle-based
codes, star formation is broadly similar in implementation, though
feedback is typically implemented in a very different way due to the
Lagrangian nature of the method \citep[see,
e.g.,][]{sh03a,sh03b,hs03}.

% Specific example: Cen & Ostriker
\subsubsection{Cen \& Ostriker}
\label{sec:starform_cen}

The Cen \& Ostriker method is a heuristic model of star formation on
galactic scales.  This method, first described by \citet{CO1992},
assumes that stars form in substantially overdense, converging, and
gravitationally unstable gas.  Algorithmically, each cell at the
locally highest level of refinement is examined at each timestep to
see if it meets the critiera for star formation.  Star particles
are allowed to form in a cell if the following criteria are met:
\begin{eqnarray}
\rho_b/\bar{\rho}_b & \geq & \eta,  \\[2mm]
\label{cendens}
\div \myvec{v}_b & < & 0, \\
\label{cencont}
t_{\rm cool} & < & t_{\rm dyn} \equiv \sqrt{3 \pi / 32G \rho_{\rm tot}}, \\
m_{b} & > & m_{J} \equiv G^{-3/2} \rho_{b}^{-1/2}c^{3} 
\left[ 1 + \frac{\delta\rho_{d}}{\delta\rho_{b}} \right]^{-3/2}
\end{eqnarray}
where 
$\eta$ is the user-defined
overdensity threshold, 
and $m_{b}$ and $m_{j}$ are the baryonic mass in the
cell and the Jeans mass of the cell, and c is the isothermal sound speed
in the cell.  If all of these criteria are met, the mass of a star
particle is calculated as \(m_{*} = m_{b} \frac{ \Delta t}{ t_{\rm
    dyn} } f_{\rm *eff} \), where $f_{\rm *eff}$ is the star formation
efficiency parameter.

If $m_{*}$ is greater than a minimum star mass $m_{\rm *min}$, a particle
is created and given several attributes: mass, a unique index number,
the time of formation $t_{\rm form}$, the local dynamical free-fall time
$t_{\rm dyn}$ and the metallicity fraction of the baryon gas in the cell
$f_{\rm Zb}$.  There is a user-defined minimum dynamical time
$T_{\rm dyn,min}$, which is observationally motivated and affects the
feedback rates (see below).  The particle is placed in the center of
the cell and given the same peculiar velocity as the gas in the cell,
and is then treated in the same manner as the dark matter particles.
An amount of gas corresponding to the new particle's mass is
removed from the cell.

The star formation algorithm creates each star particle
instantaneously.  However, feedback should take place over a
significant timescale, as all of the stars contained within the ``star
particle'' would in reality form (and massive stars would die) over a
substantial period of time.  Therefore, we assume that for the
purposes of feedback that the mass of stars formed at a time $t$ with
timestep $\Delta t$ is:
\begin{equation}
\Delta m_{\rm sf} = \int_t^{t+\Delta t} \frac{dM}{dt} dt =
 \int_{\tau_0}^{\tau_1} m_* \tau e^{-\tau} d\tau =
 m_* \left[(1+\tau_0) e^{-\tau_0} - (1 + \tau_1) e^{-\tau_1} \right]
 \end{equation}
where $\tau_0 = (t-t_{\rm form})/t_{\rm dyn}$ and $\tau_1 = (t+\Delta t-t_{\rm form})/t_{\rm dyn}$.

During this timestep, the star particle returns
metal-enriched gas and thermal energy from supernovae and from stellar
winds.  Since massive stars have very short lifetimes, we assume that
there is an immediate feedback of some fraction $f_{\rm SN}$ of the rest
energy from the stars into the gas, such
that $E_{\rm add} = f_{\rm SN} \Delta m_{\rm sf} c^2$, where c is the speed of
light.  In addition, a fraction $f_{Z*}$ of the stellar mass is fed back 
in the form of metals.  Finally, a fraction of the mass $f_{m*}$
is added back into the gas
along with momentum in order to simulate the mass ejection from
all stars (not just supernovae).

There are six user-defined parameters in this algorithm: three deal
with star formation ($\eta$, $m_{\rm *min}$ and $t_{\rm
  dyn,min}$), and three deal with feedback ($f_{\rm SN}$, $f_{Z*}$ and
$f_{m*}$).  Some of these parameters are completely free, while others
can be guided by observation or theory.  For example, the supernova
feedback parameter, $f_{\rm SN}$, can be constrained by assuming that, for
every $200 M_\odot$ of stars created, one supernova occurs, and this
event feeds back approximately $10^{51}$ ergs of thermal energy,
giving:

\begin{equation}
f_{\rm SN} = \frac{10^{51}\, \mathrm{erg}}{200~M_\odot\, c^2} \simeq 3 \times 10^{-6}
\end{equation}

The metal yield $f_{Z*}$, defined as the mass in metals produced per
unit mass of stars created, can be constrained by, e.g., the
theoretical model of \citet{1995ApJS..101..181W}.  This model suggests
that $f_{Z*} = 0.02$ is an appropriate number.  The minimum dynamical
time is set to be $t_{\rm dyn,min} = 10^7$~years to agree with the SN timescales
seen in nearby OB associations.

The other parameters, such as the overdensity threshold $\eta$,
minimum star mass $m_{\rm *min}$, and mass ejection fraction $f_{m*}$ are
not well constrained either theoretically or observationally.  Indeed,
$m_{\rm *min}$ is a purely numerical parameter designed to keep the code
from producing too many star particles, and thus has no observational
or theoretical counterpart.  The $\eta$ parameter, on the other hand,
nominally has physical meaning (i.e., the density above which star
formation must occur on a very short timescale in a self-gravitating
cloud); however, in the vast majority of simulations the densities
that are reachable are nowhere near the densities of protostellar
clouds, and thus this parameter becomes a rough proxy for finding the
densest environments in a given simulation.  

% Specific example: Kravtsov
\subsubsection{Schmidt-Law method}
\label{sec:starform_kravtsov}

This method of star particle creation is designed to reproduce
the global Schmidt law of star formation~\citep{2003ApJ...590L...1K,
  1959ApJ...129..243S}.  This algorithm is deliberately minimal, and
is explicitly geared towards modeling star formation in a
phenomenological way on kiloparsec scales.  Stars are assumed to form
with a characteristic gas timescale $\tau_*$ such that $\dot{\rho}_* =
\rho_{\rm gas}/\tau_*$. 
This ``constant
efficiency'' model on the scale of star formation regions is well
motivated observationally
\citep{1996AJ....112.1903Y,2002ApJ...569..157W}.  Star formation is
only allowed to take place in very dense regions with $\rho_{\rm gas}
\geq \rho_{\rm SF}$, where $\rho_{\rm SF}$ is a constant proper (as
opposed to comoving) density threshold.  No other criteria are
imposed.  Typical choices for $\tau_*$ and $\rho_{\rm SF}$ are $\tau_*
= 4$ Gyr and $\rho_{\rm SF} = 1.64~$M$_\odot$~pc$^{-3}$~$(n_{\rm H}
\sim 50$~cm$^{-3})$.  The adopted timescale is derived from the
observationally-determined normalization of the Schmidt law, and the
density threshold is determined by observations of star forming
regions on $\sim 100$ pc scales.

Algorithmically, the star formation events in the Schmidt-law algorithm are
assumed to occur once every global timestep (with the constraint $\Delta t_0 \leq
10^7$ years).
In cells where star formation is determined to occur (i.e. $\rho_{\rm
  gas} \geq \rho_{\rm SF}$), star particles with a mass of $m_* =
\dot{\rho}_* V_{\rm cell} \Delta t_0$ (where $V_{\rm cell}$ is the
volume of the mesh cell) are assumed to form instantaneously in a
manner similar to that described in Section~\ref{sec:starform_cen}.
The \enzo\ implementation of this algorithm is similar, except that
instead of forming stars only at the root grid timestep, we allow
stars to form at the timestep of the highest level of resolution at
any particular point in space.  As can be seen from the equation for
$m_*$ above, this can result in very small stellar masses.  To avoid
memory and processor time issues related to having very large numbers
of star particles we impose a threshold mass $M_{\rm *,min}$ such that
a star particle only forms if $m_* \geq M_{\rm *,min}$.  An
appropriate choice of $M_{\rm *,min}$ does not significantly change
the star overall star formation history of a simulation, though it may
delay the onset of star formation in a given cell relative to a
simulation without a particle mass threshold.

Each ``star particle'' is assumed to represent an ensemble of stars
and is treated as a single-age stellar population (as in the previous
section).  Kravtsov assumes that the stellar initial mass function (IMF) is described by a Miller \&
Scalo functional form with stellar masses between $0.1$ and
$100~M_\odot$ \citep{1979ApJS...41..513M}.  All stars in this IMF with
$M_* > 8~M_\odot$ deposit $2 \times 10^{51}$ erg of thermal energy
and a mass $f_z M_*$ of metals into the cell in which they form
without delay, with $f_z \equiv \min(0.2, 0.01~M_*-0.06)$
(i.e. instantaneous deposition of metals).  The definition of $f_z$ is
a rough approximation of the results of \citet{1995ApJS..101..181W}.

% Specific example: H2-regulated
\subsubsection{\HH-regulated method}
\label{sec:starform_H2reg}

The methods described in the previous two sections are generally used
in simulations that have relatively poor resolution, $\Delta x \ga 1$~kpc.
At this physical scale, individual star forming regions are not
resolved, so all uncertainty about the behavior of molecular clouds is
folded into a density threshold for star formation.  In calculations
with much higher resolution, however (on the order of a few pc), individual
molecular clouds can be resolved, thus rendering these approximations
invalid.  To that end, \citet{2012ApJ...749...36K} have implemented
a star formation algorithm that is specifically geared to
high-resolution cosmological simulations of galaxy formation, where
the formation of molecular hydrogen is followed directly and stars
are allowed to form at the highest level of refinement when the local
H$_2$ fraction exceeds a pre-determined threshold.  Cells are examined
every root grid timestep, $\Delta t_0$, and cells at the highest level
that exceed the H$_2$ limit form star particles with masses
proportional to the inferred mass of molecular hydrogen in the
star-forming region
\citep{2008ApJ...689..865K,2009ApJ...693..216K,2010ApJ...709..308M}.
The mass of the particle is calculated as:
\begin{equation} 
m_p = \epsilon \rho_{\rm gas} (\Delta x_{m})^3 \frac{\Delta t_0}{t_*}
\end{equation}
where $\Delta x_{m}$ is the resolution of the maximum level of
refinement, $\epsilon$ is an efficiency parameter with a standard
value of 0.01 \citep[as motivated by][]{2007ApJ...654..304K}, and
t$_*$ is the local free-fall timescale.

The feedback method used in this method is identical to that described
in Section~\ref{sec:starform_cen}.

% Specific example: Pop III
\subsubsection{Population III star formation}
\label{sec:starform_pop3}

Unlike in the previous sections, under some circumstances it is both
possible and desirable to simulate stars individually, rather than
treating particles as ensembles.  One particular example of this is
Population III star formation
\citep{ABN02,2007ApJ...654...66O,2008ApJ...685...40W,2009Sci...325..601T},
where a given halo may only form one star of primordial composition.
To accommodate this, \enzo\ contains a star formation algorithm that
forms individual Population III stars directly \citep{2007ApJ...659L..87A,
  2008ApJ...685...40W, 2012MNRAS.427..311W}.  Using criteria similar
to~\citet{CO1992}, a star particle forms when a cell meets all of the
following conditions:

\begin{enumerate}
\item A baryon overdensity of $5 \times 10^5$ (corresponding to a
  hydrogen number density of 
  roughly $10^3$ cm$^{-3}$ at $z=10$),

\item A converging gas flow ($\div \myvec{v} < 0$), and

\item A molecular hydrogen mass fraction f$_{\rm H2} > 5 \times 10^{-4}$,
\end{enumerate}

These are comparable to the conditions typical of a collapsing metal-free molecular cloud roughly
10 million years before the birth of a Pop III main-sequence star.  If
multiple neighboring cells are flagged as being able to form stars, a
single star is created instead.  This star has a mass that is randomly
sampled from a stellar IMF with a functional form of

\begin{equation}
f(\log M) dM = M^{\alpha} \exp \Big[ -\Big( \frac{M_{\rm char}}{M}
\Big)^{\beta} \Big]
\end{equation}

Feedback from Population III stars created using this algorithm comes in multiple forms.
Radiative feedback using the Moray radiation transport algorithm
\citep{Wise11_Moray} is available, using the mass-dependent hydrogen
ionizing and Lyman-Werner photon luminosities and lifetimes of the
Population III stars from \citet{2002A&A...382...28S}.  At the end of
their main-sequence lifetimes, explosion energies, ejected gas mass,
and ejected metal are calculated using a variety of sources that
depend on the mass of the star, and are described in detail in Section
3.2.1 of \citet{2012MNRAS.427..311W}.