A&A 467, 657-664 (2007)
DOI: 10.1051/0004-6361:20067042
K. Motoyama1 - T. Umemoto2 - H. Shang3
1 - Theoretical Institute for Advanced Research in Astrophysics, Dept. of Physics, National Tsing Hua University, 101, Sec. 2, Kuang-Fu Rd., Hsin-Chu, 30013, Taiwan
2 -
Nobeyama Radio Observatory, Nobeyama, Minamimaki, Minamisaku, Nagano 384-1305, Japan
3 -
Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 106, Taiwan
Received 29 December 2006 / Accepted 7 February 2007
Abstract
Context. Molecular clouds near H II regions tend to harbor more luminous protostars.
Aims. We investigate whether a radiation-driven implosion mechanism enhances the luminosity of protostars near regions of high ionizing fluxes.
Methods. We performed numerical simulations to model collapse of cores exposed to UV radiation from O stars. We investigated the dependence of mass loss rates on the initial density profiles of cores and variation of UV fluxes. We derived simple analytic estimates of accretion rates and final masses of protostars.
Results. The radiation-driven implosion mechanism can increase accretion rates of protostars by 1-2 orders of magnitude. On the other hand, mass loss due to photo-evaporation is not high enough to have a significant impact on the luminosity. The increase in accretion rate results in luminosity 1-2 orders of magnitude higher than those of protostars that form without external triggering.
Conclusions. Radiation-driven implosion can help explain the observed higher luminosity of protostars in molecular clouds near H II regions.
Key words: stars: formation - ISM: HII regions - methods: numerical
Massive stars play a very important role in star formation in giant molecular clouds. Strong radiation from massive stars impacts on the surrounding environment, and may trigger star formation nearby. Two scenarios of triggered star formation near H II regions have been suggested. One of them is "collect and collapse'' (Elmegreen & Lada 1977; Whitworth et al. 1994). As the H II region expands, it pushes ambient gas into a shell. Subsequent fragmentation and collapse of the swept-up shell may form the next generation of stars. For example, a core within one fragment of the molecular ring surrounding the Galactic H II region Sh104 contains a young cluster (Deharveng et al. 2003). Radiation-driven implosion (hereafter, RDI) is the other scenario (Bertoldi & McKee 1990; Bertoldi 1989). When an expanding H II region engulfs a cloud core, the strong UV radiation will compress the core. Bright-rimmed clouds of cometary shapes found at the edges of relatively old H II regions are explained by Lefloch & Lazareff (1994) with a RDI model. To understand star formation in giant molecular clouds, we study the effects of the RDI.
Observational results indicate that the H II regions can strongly affect star formation nearby. Dobashi et al. (2001) investigated the maximum luminosity of protostars as a function of the parent cloud mass. They found that molecular clouds near H II regions contain more luminous protostars than others. Moreover, Sugitani et al. (1989) showed that the ratios of luminosity of the protostar to core mass of bright-rimmed globules are much higher than those in dark globules. These results indicate that the RDI affects the luminosity of protostars near H II regions.
Recent studies suggest that star formation triggered by external
compression can increase accretion rates and luminosity of
protostars. Motoyama & Yoshida (2003) investigated collapses of
cores compressed by external shock waves and found that accretion
rates were enhanced. Some observations support this result.
Belloche et al. (2006) observed the Class 0 protostar IRAS 4A
in the star formation cloud NGC1333, whose star formation is suspected
to be triggered by powerful molecular outflows
(Sandell & Knee 2001). They showed a numerical model of
collapse triggered by external pressure, which reproduces the
observed density and velocity profiles. Moreover, they deduced a
high accretion rate of (0.7-2)
yr-1. Their result indicates that external compression increases the
accretion rate of a protostar, which produces high accretion
luminosity. In this paper, we attribute the origin of higher
luminosity in protostars to star formation activities induced by the
RDI.
Some studies have been made on RDI. Bertoldi (1989) and Bertoldi & McKee (1990) developed an approximate analytical solution for the evolution of a cloud exposed to ionizing UV radiation. Lefloch & Lazareff (1994) investigated the RDI with numerical simulations. However, their calculations did not include self-gravity of the gas. Recently, SPH simulations including self-gravity have been performed. Kessel-Deynet & Burkert (2003) investigated the effects of initial density perturbations on the dynamics of ionizing clouds, and Miao et al. (2006) focused on cloud morphology. Despite the many studies on RDI, little attention has been given to the effects of RDI on accretion rates.
In this paper, we explore whether the RDI can enhance accretion rates by orders of magnitude compared to those of protostars that form without an external trigger. RDI is expected to enhance the accretion rates as a result of strong compression, and it would increase the luminosity of protostars. On the other hand, photo-evaporation of the parent core due to RDI decreases the mass of the protostar, and would subsequently decrease the luminosity of the new protostar. We investigate which effect dominates.
The paper is organized as follows. In Sect. 2, we describe our numerical method. In Sect. 3, we present our results. In Sect. 4, we discuss luminosities of protostars based on our numerical results. In Sect. 5, we summarize our main conclusions.
We simulate evolution of cores exposed to diffuse UV radiation. We
assumed that a core of neutral gas is immersed within an H II
region and is exposed to diffuse radiation from the surrounding
ionized gas. In an actual H II region, the core is not only
exposed to diffuse radiation but also to direct radiation from the
ionizing star. The diffuse flux of Lyman continuum photons is ![]()
of the direct flux of Lyman continuum photons
(Canto et al. 1998; Whitworth & Zinnecker 2004). The spherical core
has an effective cross-section to diffuse radiation of
and
to direct UV radiation of
,
respectively. The ratio of total
number of diffuse photons to total number of direct photons is
0.6. Our assumption underestimates the total amount of ionizing UV
photons by a factor of
2.6. Since the purpose of this paper
is to investigate the effects of RDI on the accretion rate of a
protostar, we neglected direct radiation and considered only diffuse
radiation for simplicity. Effects of direct radiation and morphology
of the core are beyond the purpose of the present paper, and will
be discussed in our next paper using two dimensional simulations.
The ionizing UV photons enter from the outer boundary of the
computational domain. The number flux of these UV photons is
expressed as
We assumed that the sound speed is a function of the ionization
fraction x. The hot ionized gas has a temperature of
104 K in an H II region. To imitate conditions of the
hot ionized gas (x=1), the sound speed of the ionized gas
was
set to 13 km s-1. The neutral gas (x=0) has two
states: warm gas, if the density is lower than a critical density
g
cm-3, where
is the mass of hydrogen atom, and cold gas,
if the density is higher than the critical density
.
We
introduce this artificial condition for the initial pressure balance
between the core and the rarefied warm gas surrounding it. The
rarefied warm gas is ionized in a very short time and will hardly
affect the time evolution of the core. The sound speed of neutral
gas
is expressed as
![]() |
(2) |
![]() |
(3) |
We assume that the initial core is a Bonnor-Ebert sphere, whose mass
M0 is set to 3
or 30
.
A Bonnor-Ebert sphere
is a self-gravitating isothermal sphere in hydrostatic equilibrium
with a surrounding pressurized medium (Bonnor 1956).
Gravitational stability of a Bonnor-Ebert sphere is characterized by
a non-dimensional physical truncation radius
,
where
is the truncation radius,
G is the gravitational constant, and
is the central
density. If
is smaller than the critical value
,
the Bonnor-Ebert sphere is stable. All
Bonnor-Ebert spheres that we use in the calculations are
gravitationally stable (
).
It is reasonable to assume that the initial density distribution of
the core has spherical symmetry. The density distribution of the core is
assumed to remain unchanged while the core is immersed within an
expanding H II region. The sound crossing-time of the
Bonnor-Ebert sphere in model C is
yr, and the
expansion velocity of the H II region is approximated well
with
at
,
where
and
are the radii of the H II region and the
Strömgren sphere. If we assume
pc and an ambient
density of
,
the timescale on which the
H II region engulfs the core,
,
is
yr. This timescale is shorter than the sound-crossing
time, and the adoption of spherical symmetry is acceptable for
simplicity.
We have computed seven models of collapse of cores exposed to UV
radiation from an O star. We summarize the parameters of our
simulations in Table 1. For the three models
(A, B, and C) we assume that the mass of the Bonnor-Ebert sphere is 3
and the ionizing star is an O8V star. These three models
differ in the non-dimensional physical truncation radii of the
Bonnor-Ebert spheres
.
Bonnor-Ebert spheres with larger
have higher central densities and smaller radii. For
other three models (D, E, and F) we assume that the ionizing star is
an O4V star. In model G, we computed the case of a more
massive 30
core, exposed to UV radiation from an O8V star.
We perform simulations in one spatial dimension using the Godunov
method. Our hydrodynamical code has second-order accuracy in both
space and time. The hydrodynamic equations are
![]() |
(5) |
We solved the radiative transfer equation for the diffuse UV
radiation using a two-stream approximation.
Canto et al. (1998) first treated the diffuse UV radiation
using an iterative two-stream approximation.
Pavlakis et al. (2001) simplified the method and showed that a
non-iterative version serves as a good approximation for
centrally-peaked density distributions. We adopted the latter for
our Bonnor-Ebert spheres. The equations of evolution for the
ionization fraction and propagation of diffuse UV radiation are
given by
We verified our procedure above by a test problem described in Lefloch & Lazareff (1994). We simulated the time evolution of uniform neutral gas illuminated by an ionizing flux increasing linearly with time in one dimensional Cartesian coordinates. An ionization front propagates with constant velocity, and neutral gas accumulates ahead of the ionization front. We adopted the same parameters as Lefloch & Lazareff (1994). Deviations between numerical results and the analytic solutions were less than 1% and 5% for the density and the velocity, respectively. Moreover, we checked the accuracy of our code in spherical coordinates by comparing with the self-similar solution of Shu (1977).
The sink cell method (e.g. Boss & Black 1982) was used to
avoid the long computational times encountered near the center. As
the gravitational collapse proceeds, the infall velocity increases
at the central region and the time step that satisfies the
Courant-Friedrichs-Levy (CFL) condition, which is necessary for
numerical stability, reduces. If the central density exceeds the
reference density
,
the central 10 grids within 10 AU will be treated as sink cells. The mass that
enters the sink cell is treated as a point mass located at r=0.
The boundary condition at the outer boundary of the sink cell was given
by extrapolation from the neighboring cells. Since the infall
velocity near the sink cell is supersonic, this method does not affect
the flow of outer gas.
We used non-uniform grids in order to obtain high resolution at the
central region. Sizes of the ith grid from the center were given as
,
where
is 0.02, and the size of the most inner grid is 1 AU. The outer boundary is located at
.
For example, we allocate 1473 grids in model C. For the core itself, 1200 grids were allocated. Reflecting and
free boundary conditions were imposed at the inner and outer
boundaries, respectively.
Table 1: Parameters adopted in the simulations.
![]() |
Figure 1:
Time evolution of a) density profile, b) velocity profile,
c) ionization fraction and d) pressure profile in the model with
|
| Open with DEXTER | |
To illustrate common features of the time evolution, we describe
results in model C with
listed in Table 1. Figure 1 shows the time
evolution of density, velocity, ionization fraction, and pressure.
Initially, the core was in hydrostatic equilibrium with the rarefied
warm neutral gas. At t=0, UV radiation, whose UV photon flux is
,
started and radiated from the outer boundary. The outer
layer of the core was ionized by UV radiation. Temperature and
pressure at this outer layer increased drastically. The hot ionized
gas then expanded outwards, accelerated to
40 km s-1, and drove a shockwave into the core. At
yr, the shockwave reached the center of the core. The
dashed line in Fig. 1a shows that the core is
compressed into a small region of radius
at this time. The density becomes
103 times higher than the
initial density and the core becomes gravitationally unstable.
The core collapses subsequently. We see from Fig. 1a and b that density and velocity of the compressed core
increase rapidly after
due to
self-gravity. After the central density exceeds the reference value
,
the sink-cell treatment turns on, at which point a
protostar forms at the center of the core. At this phase, inner
region of the core accretes onto the central protostar, and the
ionized outer layer expands outward.
Velocities of the shock front and the ionization front remain nearly
constant, on the other hand. Figure 2 shows locations of
the ionization and the shock fronts as functions of time. The velocity
of the shock front is a constant value of
.
The analytic estimation described in Sect. 4.1 shows a similar velocity of
.
For convenience, we define the ionization front at the location where
the ionization fraction is 0.9. How the ionization front is defined
hardly affects the results, since, as shown in Fig. 1c, the transition layer from neutral to ionized gas is very
thin. The average velocity of the ionization front is
,
and analytic estimation shows a similar
velocity of
.
Some fraction of the incident UV photons is absorbed in the layer of
the photo-evaporated gas. In this layer, ionized hydrogen atoms
recombine with electrons at the rate of
per unit
volume per unit time. Some UV photons are used for ionization of
these recombined hydrogen atoms. We define Q, the ratio of UV photon
flux at the ionization front to the incident UV photon flux, as the
photo-evaporation rate of the core. Figure 3 shows Q as a
function of radius. Spitzer (1978) derived an analytic
estimation of Q as
As a reference, we simulate the spontaneous collapse of the core
without UV radiation. Since the thermal pressure balances the
gravitational force in the Bonnor-Ebert sphere, we slightly enhance
the density to
,
where
is the
density of a marginally stable Bonnor-Ebert sphere (
)
to initiate collapse. This core whose mass is 3
collapses by self-gravity, and the average value of the accretion rate
is
.
The RDI indeed increases accretion rates. Figure 4
shows the accretion rate as a function of time in model C with
.
The accretion rates are measured at the radius of 10 AU. The horizontal line indicates the average accretion obtained
from the reference model. At
,
the
central density reaches the threshold density
,
and the
sink cell method turns on. The accretion rate increases rapidly
afterwards, and lasts for
104 yr.
Figure 5 shows the average values of accretion
rate as a function of distance from the ionizing star. The
accretion rates increase with decreasing distance from the ionizing
star, and with increasing UV photon flux. If the cores are
located within 10 pc from an ionizing star, the accretion
rates can go 1-2 orders higher than models without UV radiation.
The accretion rates not only depend on the UV photon fluxes but also
on
and M0. Cores with larger values of
have smaller
radii and higher densities. Shock velocities and densities in
post-shock regions are larger in the models with small
.
Gas of higher density has a shorter free-fall time, and reaches an
accretion rate higher in the models with smaller
.
We will
compare this with analytic estimates. In the
model G, not only a small
free-fall time but also the large mass of compressed core is responsible
for the high
accretion rate.
Photo-evaporation can reduce the actual mass that accretes onto the
protostar. In the outer layer, hot ionized gas expands outward at a
velocity up to
.
Photo-evaporation
of gas takes place at
,
except for the case of
very small cores. The escape velocity of the Bonnor-Ebert sphere in
model C is
,
for example, is much
smaller than that of the photo-evaporated flow. Photo-evaporation
prevents a large amount of gas from accreting onto the protostar, and
the final mass of the protostar can be small due to mass loss by photo-evaporation.
The final masses of protostars decrease with decreasing distance from
the ionizing star. Figure 6 shows the final masses of protostars as functions of distance from the ionizing star. The mass
of a protostar is defined as mass that has entered sink cells, and
calculations continue until accretion rates become sufficiently
small (
yr-1). Mass loss due to
photo-evaporation at the surface of the core is proportional to the
incident UV photon flux and final masses of protostars are smaller
near the ionizing star. For an O8V or an O4V star, mass loss due to
photo-evaporation varies from several percent to several tenths of
a percent of the initial core mass. The more massive the ionizing star
becomes, the smaller the final mass.
![]() |
Figure 2:
Location of the ionization front and shock front.
The thick (black) and thin (red) solid lines label the
position of shock front and ionization front in model C with
|
| Open with DEXTER | |
![]() |
Figure 3:
The ratio of UV photon flux at the ionization front to
incident UV photon flux as a function of radius in the model C with
|
| Open with DEXTER | |
The final masses not only depend on the UV photon fluxes but also on the sizes of cores. Figure 6 shows that mass loss due to photo-evaporation
is significant in the model G, in which the initial radius of the core is 10 times larger than that in the model C. The larger core has
a larger photo-evaporation rate and longer time to be exposed to UV radiation
because of the larger surface area and longer crossing time of the shock wave. This is the reason why the core loses large amounts of
mass in model G. For the same reason, final masses are small in the
models with small
.
In this section, we derive analytic estimates of final mass and
accretion rates of protostars under the influence of RDI. For
simplicity, we assume that a core initially has a uniform density n0 and a
radius R0. The core is exposed to UV radiation whose
number flux of UV photons is
.
Spherical symmetry is assumed.
![]() |
Figure 4:
Accretion rate as a function of time in the model C with
|
| Open with DEXTER | |
![]() |
Figure 5: The average accretion rates as functions of distance from the ionizing star. Open diamonds, open squares, open circles label results from the model A, B, and C, respectively. Filled diamonds, filled squares, filled circles label models D, E, and F, respectively. The crosses are for the model G. The horizontal line indicates the average accretion rate in the reference model. The dotted line, dash-dotted line and dashed line indicate analytic estimation for models A, B, and C, respectively. |
| Open with DEXTER | |
The final mass of protostar
is expressed as
![]() |
Figure 6: Ratio of final mass of protostar to initial core mass as a function of distance from the ionizing star. The open diamonds, the open squares, the open circles indicate results in the models A, B, and C, respectively. The filled diamonds, the filled squares, the filled circles indicate results in the models D, E, and F, respectively. The crosses indicate results in the model G. The dotted line, the dash-dotted line and the dashed line indicate analytic estimations for model A, B, and C, respectively. |
| Open with DEXTER | |
The mass loss flux
due to photo-evaporation is expressed as
The shock velocity
is derived following
Bertoldi (1989) to estimate the shock crossing time
.
Assuming that ram pressure of the flow equals thermal pressure of
the post-shock region, we obtain
![]() |
(19) |
The accretion rate, over a free-fall time, is roughly estimated by
![]() |
(23) |
Table 2: Star formation regions of potential triggered origins.
The luminosity of a protostar is mainly produced by accretion, because
nuclear burning does not influence much of the luminosity during the
protostellar collapse phase. The accretion luminosity L is
expressed as
Enhanced luminosity around bright-rimmed clouds are seen in
observations. Sugitani et al. (1989) observed three
bright-rimmed globules associated with cold IRAS sources in
13CO(J=1-0) emission. They compared the physical parameters with four
isolated dark globules and found that the ratios of the luminosity
to globule mass, 3-13
,
are 1-2 orders higher
than the values of isolated dark globules, 0.03-0.3
,
as shown by our calculations. These objects are
good candidates for influence by the RDI.
Protostars near H II regions may spend most of their lifetime
in the state of high luminosity. If the duration of the high accretion
rate is much shorter than the lifetime of the protostar, it is
relatively difficult to catch the protostar in its high-luminosity
state by observations. In the simulations, the protostellar phase
begins from the time when the sink-cell method is turned on and
lasts until 99% of the final mass has accreted within the sink
cells. The lifetime of the protostar in model C is
yr with d=9 pc. As shown in Fig.
4, the accretion rate is enhanced for a duration
of
,
comparable to the lifetime of the
protostar. It is likely that protostars near H II regions are
observed in the state of higher luminosity.
Evidence of triggered star formation has been suggested at the edges of many H II regions. Some observations favor the scenario of sequential star formation around a massive star. Lee et al. (2005) selected candidates of pre-main sequence stars based on 2MASS colors in the Orion region, and revealed the spatial distribution of classical T Tauri stars using their spectroscopic observations. The pre-main sequence stars form between the ionizing star and the bright-rimmed clouds, but not deeply inside the cloud. The youngest stars are found near the cloud surfaces, forming an apparent age gradient. As the H II region expands, it might be able to trigger star formation sequentially from near to far. Several such examples are summarized in Table 2.
RDI may contribute to the triggering mechanism around a massive
star. Lee et al. (2005) found that only bright-rimmed
clouds with strong IRAS 100
and H
emission show
signs of ongoing star formation. This suggests that star formation
may have been triggered by strong compression from the ionization
shock, because these emissions are tracers of dense gas and
ionization. Small-scale age gradients of YSOs near the ionization
fronts are seen in some molecular clouds.
Fukuda et al. (2002) observed heads of molecular pillars
and
in the Eagle Nebula (M 16). They found that
protostars and starless cores are ordered in age sequence. The more
evolved objects are closer to the cloud heads (and the ionizing
star) and the younger objects are farther from the exciting star.
Similar age gradients are also seen in bright-rimmed clouds
(Sugitani et al. 1995) and the cometary globule IC1396N
(Getman et al. 2007). These age gradients near the ionization front
seem to suggest that star formation may be sequentially triggered by
ionization shocks.
Stars formed facing the H II regions appear to be more luminous. Yamaguchi et al. (1999) carried out 13CO(J=1-0) observations toward 23 southern H II regions and identified 95 molecular clouds associated with H II regions. They found that IRAS point sources tend to be more luminous on the side facing an H II region. Their result is consistent with enhanced luminosity as a result of the RDI.
Shock speeds derived from observations are similar.
White et al. (1999) and Thompson et al. (2004) derived
shock velocities from pressure in pre- and post-shock regions.
Getman et al. (2007) derived shock velocity from the age
gradients of YSOs. Their estimated values of a few
are consistent with our results. For example,
Thompson et al. (2004) estimated that the shock velocity in
bright-rimmed cloud SFO11 is
.
SFO11
is illuminated by an O6V star whose projected distance is 11 pc.
These conditions are close to those in model C with d=9, whose
ionizing star is an O8V star. As described in section
3.1, the averaged shock velocity is
in this model.
Lee & Chen (2006) found that classical T Tauri stars exist near the surfaces of bright-rimmed clouds and intermediate mass stars exist relatively farside clouds. One possible difference in stellar mass may be the density or mass in the parent cores. Photo-evaporation may reduce the mass that would otherwise accrete onto the stars. In our simulations a few percent to several tenths of a percent of the initial mass may be lost by photo-evaporation. The mass of young stars may be smaller near the surface of bright-rimmed clouds compared to inner regions.
We simulate the evolution of cores exposed to strong UV radiation around H II regions. We investigate possible effects of radiation-driven implosion on the accretion luminosity of forming protostars. We also developed analytic estimates of accretion rates and final masses of protostars, which are compared to the numerical results.
The main findings are:
Acknowledgements
The authors would like to thank Wen-Ping Chen, Huei-Ru Chen, Jennifer Karr, and Ronald Taam for helpful discussions. This work is supported by the Theoretical Institute for Advanced Research in Astrophysics (TIARA) operated under Academia Sinica, and the National Science Council Excellence Projects program in Taiwan administered through grant number NSC95-2752-M-001-008-PAE, and NSC95-2752-M-001-001. Numerical computations were in part carried out on VPP5000 at the Astronomical Data Analysis Center, ADAC, of the National Astronomical Observatory of Japan.