© ESO, 2011

1. Introduction

The ultraviolet radiation from a massive star in an environment of low mass star formation produces several interesting phenomena, among them the appearance of the so-called proplyds (an acronym that stands for PROto-PLanetarY Disk; see O’Dell et al. 1993; and O’Dell & Wen 1994). Proplyds are young stellar objects (YSO) with a cometary-shaped ionization front that contain a central, low mass young star surrounded by an accretion disk, a photoevaporated wind, and, in some cases, a monopolar or a bipolar outflow. The proplyds in the Orion nebula (M 42) were first reported as unresolved optical emission features by Laques & Vidal (1979) and then as radio sources by Churchwell et al. (1987). Using the Hubble Space Telescope (HST), O’Dell & Wen (1994) and O’Dell & Wong (1996) reported the discovery of tens of proplyds around the Trapezium cluster.

The proplyds present strong optical emission lines such as Hα, [O iii] 5007 Å, [O i] 6300 Å  (Bally et al. 2000) and also some ultraviolet (Henney et al. 2002) and IR lines (Takami et al. 2002). With HST high spatial resolution images and spectra, it was possible to see that the emission comes from different parts of the object. Proplyds present a cometary shaped feature, very strong in Hα. Several of them show another bow shaped feature, detached from the cometary structure itself, visible in Hα, [O iii] 5007 Å and in IR lines, lying further from the low mass star than the ionization front (see Fig. 1). In some objects, the disks are seen in emission in [O i] 6300 Å. It is also possible to see that the optical lines are double or even triple peaked. For example, in the central parts of the proplyd 167 − 317 (LV2), Hα has three components: a low velocity component, peaking at a radial velocity of  ~30 km s^-1, and two high velocity components, centered at  ~100 km s^-1 and  ~−75 km s^-1 (Vasconcelos et al. 2005).

Johnstone et al. (1998) explained the observations with a model in which the radiative field and the wind from θ¹ Ori C interact with a low mass YSO, generating a photoevaporated wind, an ionization front (IF) and a bow shock. They also proposed the existence of a photodissociated region (PDR), in which the photoevaporated wind is neutral (the H is neutral). They separate their models in two classes that depend on the amount of EUV and FUV flux arriving at the proplyd’s disk surface. The EUV model, valid for a YSO near to or very far from θ¹ Ori C (d < 10¹⁷ cm or d > 10¹⁸ cm), presents a very thin PDR. In this case, the EUV flux determines the photoevaporated mass loss rate. On the other hand, for the FUV model, the PDR is thick and can contain a shock front within it. Here, the FUV flux determines the mass loss rate. Störzer & Hollenbach (1999) further improved the models by Johnstone et al. (1998) including results from PDR codes.

Some axisymmetric numerical simulations were also done in order to reproduce the characteristics of proplyds. García-Arredondo et al. (2001) carried out axisymmetric simulations of the interaction of the stellar wind from θ¹ Ori C with the photoevaporated proplyd flow. They took available observational data to constrain stellar parameters such as the wind density and velocity and the ionization front parameters. Their numerical calculations also assume a hemispherical proplyd head and a cylindrical tail. Studying the interaction between these two winds, they have successfully reproduced the arc emission for the proplyds near θ¹ Ori C.

Richling & Yorke (2000) performed 2D, axisymmetric hydrodynamical simulations of photoevaporating disks including both ionizing (hν ≥ 13.6 eV) and dissociating (6   eV < hν < 13.6 eV) radiation. In their models, disk structures are formed through collapse simulations of 1 and 2 M_⊙ molecular clumps (Yorke & Bodenheimer 1999), which are then exposed to the radiation field by switching on the external UV radiation field in the calculation. They studied the effects of distance from the UV photon source on emission line maps and the effect of the presence of a spherical wind from the proplyd star, which promotes the appearance of collimated, bipolar microjets. Instead of calculating the dissociation of the H₂ molecule, they followed the ionization of C. They took into account photoelectric heating and cooling by fine-structure lines such as [C ii] 158 μm, [O i] 63 μm and [O i] 145 μm. An important point in their work is the inclusion of the diffuse radiation field, which is responsible for the tails of the proplyds (see Fig. 1).

In this work we present the first fully three-dimensional numerical simulations of disks exposed to FUV and EUV radiation fields. As discussed in Johnstone et al. (1998), the diffuse ionizing and dissociating radiation field is important for determining the shape of the ionization front behind the (directly) illuminated disk surface. Our present simulations do not include the diffuse field, resulting in the tails that do not have the correct morphology (see Cerqueira et al. 2006b). We concentrate on obtaining a description of the head of the proplyd flow, and study the effect of different orientations of the disks with respect to the impinging UV photon field (and the stellar wind).

The present paper is organized as follows. In Sect. 2, we describe the code and the simulations. In Sect. 3, we present our results. Finally, in Sect. 4, we draw our main conclusions.

2. The numerical method, the numerical setup, and the simulated models

The 3D numerical simulations have been carried out with the YGUAZÚ-A adaptive grid code (Raga et al. 2000, 2002) using a 5-level binary adaptive grid. The YGUAZÚ-A code integrates the gasdynamic equations employing the flux vector splitting scheme of van Leer (1982) together with a system of rate equations for atomic/ionic species. In our simulations, we consider 4 species: H i, H ii, C i and C ii. This code has been extensively employed for simulating different astrophysical flows such as jets (Masciadri et al. 2002; Cerqueira et al. 2006a), interacting winds (González et al. 2004), photoevaporating clumps (Cerqueira et al. 2006b) and supernova remnants (Velázquez et al. 2001a, 2004). It was also tested with laser generated plasma laboratory experiments (Sobral et al. 2000; Raga et al. 2001; Velázquez et al. 2001b).

In this work, instead of solving an energy equation, we prescribe a temperature law, given by $\mathit{T}\mathrm{=}\mathrm{(}{\mathit{T}}_{\mathrm{1}}\mathrm{-}{\mathit{T}}_{\mathrm{2}}\mathrm{)}\mathrm{·}{\mathit{x}}_{\mathrm{HII}}\mathrm{+}{\mathit{T}}_{\mathrm{2}}\mathrm{·}{\mathit{x}}_{\mathrm{CII}}\mathrm{+}{\mathit{T}}_{\mathrm{3}}\mathrm{·}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{x}}_{\mathrm{CII}}\mathrm{)}\mathit{,}$(1)where T₁ = 10   000 K, is the characteristic temperature of an H ii region, T₂ = 1000 K, the typical PDR temperature, T₃ = 10 K, the temperature of the molecular gas, x_H   II is the hydrogen ionization fraction and x_CII is the carbon ionization fraction (x_C   II = 1 when all the C   I is C   II). This prescription is justified if the thermal equilibrium time scale is much smaller than the dynamical time scale (Lefloch & Lazareff 1994), which is the case here, for both the ionized and the molecular gas. It is clear from Eq. (1), that we have possible temperatures ranging from 10 K to 10⁴ K. In order to study the dependence of the PDR geometry with temperature, we also compute models in which we set T₂ = 3000 K.

Following Richling & Yorke (2000) we do not treat the dissociation of H₂, but consider that the dissociation front and the carbon ionization front coincide. We then solve rate equations for hydrogen and carbon taking into account photoionization, collisional ionization, radiative and dielectronic recombination.

We model the O star as a point source located outside the computational domain. This source emits EUV photons (with hν > 13.6 eV) at a rate S_EUV = 7.2 × 10⁴⁸ s^-1 and FUV photons (with 6   eV < hν < 13.6 eV) at a rate S_FUV = 1.78 × 10⁴⁹ s^-1 (or zero, in order to study the effects of the EUV radiation field only; see Sect. 3 below).

For the transfer of FUV photons, we consider the optical depth due to the ionization edge of the C I photoionization cross section, and to the absorption of photons by dust: ${\mathit{\tau }}_{\mathrm{FUV}}\mathrm{=}\mathrm{1.22}\mathrm{\times }{\mathrm{10}}^{-17} \mathrm{c}{\mathrm{m}}^{\mathrm{2}}{\mathit{N}}_{\mathrm{CI}}\mathrm{+}\mathrm{1.15}\mathrm{\times }{\mathrm{10}}^{-21} \mathrm{c}{\mathrm{m}}^{\mathrm{2}}{\mathit{N}}_{\mathrm{HI}}\mathit{,}$(2)where N_CI is the C I column density integrated along the path of the ray. The second term (representing the FUV dust absorption) is calculated assuming that the optical extinction is given by A_V = N_HI / 2 × 10²¹      cm^-2 (appropriate for a standard grain distribution). We then multiply this number by 2.4 to obtain the FUV extinction. This is consistent with empirically derived optical/UV extinction curves (see, e. g., Cardelli et al. 1988). We are also implicitly assuming that there is no dust in the region where H is fully ionized.

Our simulations include a flat, neutral disk, with a constant thickness H = 15.6 AU and an outer radius of 50 AU. The disk is positioned at distances of 0.02, 0.1 and 0.2 pc from the ionizing photon source. With this particular choice of parameters, we can in principle ensure that the photoevaporated flow of the nearest disk would be EUV-dominated, while for the more distant one, FUV-dominated. Using previous analytical calculations of Johnstone et al. (1998) and the disk radius given above, we can roughly estimate the minimum distance from the ionizing source in order to have an FUV-dominated flow as d_min ~ 0.1 pc. The disk radius that would support an FUV-dominated flow is, for d = 0.1 pc, r_d ≈ 30 AU, and r_d ≈ 120 AU for d = 0.2 pc (see Eqs. (15) and (16) in Johnstone et al. 1998). We choose to simulate disks at distances from the ionizing photon source that are i) at a distance that is suitable for an EUV-dominated flow (d = 0.02 pc), ii) close to the limiting cases for which we should expect an FUV dominated flow to occur (d = 0.1 pc, r_d = 50 AU), and iii) FUV-dominated flow (d = 0.2 pc). In this last case, the disk radius r_d = 50 AU is found to obey the r_d ≲ 122 AU condition, and the system would develop an FUV-dominated flow. The only caveat is that, as pointed out by Johnstone et al. (1998), there are several sources of uncertainties in evaluating these (analytical) numbers.

The density of the disk is assumed to be constant in both the vertical and radial directions. We start the simulations assuming pressure equilibrium between the PDR region and the disk itself: ${\mathit{n}}_{\mathrm{d}}\mathrm{=}{\mathit{n}}_{\mathrm{0}}\mathrm{(}{\mathit{T}}_{\mathrm{2}}\mathit{/}{\mathit{T}}_{\mathrm{3}}\mathrm{)}\mathit{,}$(3)where T₂ and T₃ are the PDR and disk temperatures (see Eq. (1)) and n₀ is the number density at the base of the photoevaporated flow, which can be estimated as: ${\mathit{n}}_{\mathrm{0}}\mathrm{\simeq }{\mathit{N}}_{\mathrm{D}}\mathit{/}{\mathit{r}}_{\mathrm{d}}\mathit{,}$(4)for a FUV dominated flow. Here, N_D is the column density of the PDR attained for τ_FUV ~ 1 − 3, that is of the order of 10²¹ cm^-2 (Johnstone et al. 1998). Using r_d = 50 AU, we obtain n₀ ~ 10⁶ cm^-3. We therefore set n_d = 10⁷ cm^-3. For an EUV dominated flow, the density is found by equating the EUV flux to the recombinations in the ionized flow: $\frac{{\mathit{S}}_{\mathrm{EUV}}}{\mathrm{4}\mathit{\pi }{\mathit{d}}^{\mathrm{2}}}\mathrm{=}\mathrm{3}{\mathit{\alpha }}_{\mathrm{r}}{\mathit{n}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{d}}\mathit{,}$(5)where α_r = 2.6 × 10^-13 cm³ s^-1 is the recombination coefficient for Hydrogen at 10⁴ K, d is the distance from the disk to the ionizing photon source, and S_EUV = 7.2 × 10⁴⁸ s^-1 and r_d = 50 AU in all of the models presented in this paper. We should note that the factor of 3 on the right hand side of Eq. (5) is a result of assuming a spherical divergence for the photoevaporated wind.

We then calculate the densities n₀ for different distances and find n₀ = (0.98 − 5.0) × 10⁵ cm^-3 for d = (0.1 − 0.02) pc. For all cases, we treat the disk as an infinite reservoir, since at each time step we reset its initial conditions, as described above.

We assume the presence of a 0.2 M_⊙ star at the disk center and consider the gravitational force that it exerts on the flow (by including the appropriate source terms in the momentum and energy equations). A Keplerian rotation velocity law (consistent with the 0.2 M_⊙ central star) is imposed on the disk material. The gravitational force is important as it inhibits the formation of a photoevaporated wind in the central regions of the disk. We expect a photoevaporated flow to start beyond the gravitational radius, which can be estimated as r_g = GM _⋆  / a², where a is the local sound speed. For FUV dominated flows, a ≃ 3 km s^-1, so that r_g = 2.9 × 10¹⁴ cm, or  ≈ 20 AU. The inhibition of the photoevaporated wind in the central regions of the disk can be seen in our FUV dominated flow simulations (which have r_g ≈ 20 AU), but not in our EUV dominated flow simulations (which have r_g ≈ 2 AU, smaller than our numerical resolution). The size of the remnant disk, which is known to be typically smaller than the gravitational radius, will not be investigated in this work, since we are treating the disk as an infinite reservoir (Johnstone et al. 1998).

For some models, we add a plane-parallel wind (flowing parallel to the ionizing/dissociating external radiation field, with velocity v_w = 100 km s^-1 and with number density n_w = 500 cm^-3), which represents the wind from the O-type star. The wind velocity of θ¹ Ori C is estimated to be of the order of 500 to 1650 km s^-1 (Bally et al. 1998), or even larger (up to about 3600 km s^-1, Walborn & Nichols 1994), which is by far larger than the adopted velocity here. On the other hand, the ram pressure estimated for both the stellar wind and the photoevaporated flow for proplyds in the Trapezium region is  ~ 3 × 10^-8 dyn cm^-2 to 5 × 10^-7 dyn cm^-2 for the proplyds closer to θ¹ Ori C (at distances  ≃ 0.01 pc; see, e.g., Garcia-Arredondo et al. 2001). Therefore, our particular choice of parameters for the stellar wind give us a ram pressure comparable to that estimated for a group of proplyds in the Trapezium cluster. In particular, with the adopted parameters, we have P_sw = 8.3 × 10^-8 dyn cm^-2.

The initial setup of the simulations is shown in Fig. 2. The flat, dense structure with which we model the circumstellar disk is placed on the yz-plane, at x = 1   203 AU. The disk axis, which we define here as the axis perpendicular to the disk plane, and that contains the disk center at (x,y,z) = (1203, 334, 334) AU, coincides with the x-axis for an inclination angle θ = 0°. For larger inclination angles, the disk axis rotates clockwise, with respect to the x-axis, as depicted in Fig. 2. The remaining volume of the computational domain is filled with an ionized environment of uniform density n_a = 500 cm^-3. All models except M1a have an initial wind of v_w = 100 km s^-1 and density n_w = 500 cm^-3, entering the computational domain through the x = 0 boundary. The ionizing/dissociating photon source is placed along the x-axis, outside the computational domain to the left.

In Table 1 we present all the relevant parameters used in our models: the dissociating photon rate, S_FUV (in units of 10⁴⁹ s^-1; second column), the distance from the disk to the external star (in pc; third column), the inclination angle θ between the disk and the x axis (fourth column), the velocity of the wind from the external star (in km s^-1; fifth column) and the temperature of the PDR region (in K; sixth column). All the models have an ionizing photon flux S_EUV = 7.2 × 10⁴⁸ s^-1 (similar to the values used by Richling & Yorke 2000) and r_d = 50 AU.

The simulations have been computed in a 5-level, binary adaptive grid with a maximum resolution (along the three axes) of 5.2 AU for all models. The computational domain has a size (x_max,y_max,z_max) = (1337, 668, 668) AU.

3. Results and discussion

In Fig. 3, we show the xz-midplane density stratifications (left) and Hα maps (right) for model M1a (top) and model M1b (bottom) at t = 2000 years. The models present only an EUV radiation field from the ionizing source (at a 0.1 pc distance from the disk). The difference between the two models is the presence of a stellar wind in model M1b, as described in Table 1. The disk axis has no inclination with respect to the direction of the impinging photons (i. e., θ = 0° see Fig. 2 and also Table 1).

The Hα emission in model M1a, is restricted to a region near the disk surface (and facing the ionizing source), and the edge of the shadow on the non-illuminated side of the disk. For model M1b, we have Hα emission close to the disk (with morphology and intensity similar to those seen in model M1a), and also fainter emission regions distributed along the interaction region between the stellar wind and the photoevaporated flow. It is known that some proplyds in Orion show faint Hα arcs facing θ¹ Ori C (Bally et al. 1998), which correspond to the stellar wind/photoevaporated flow interaction emission found in our models.

In Fig. 4 we show xz-midplane density maps for model M1b at t = 60, 80, 100 and 120 yr, in order to illustrate the initial evolution of such an interaction. At t = 60 yr, the stellar wind and the photoevaporated flow start to interact at x ≃ 668 AU. At later times, the interaction region moves towards the disk, and stabilizes around the position of the point of ram pressure balance between the stellar wind and the photoevaporated flow. We can also see this effect in Fig. 5, where we show the photoevaporated wind mass loss rate (Ṁ; dashed line), the distance from the disk to the interaction region i) as measured directly from the xz-midplane density maps (d_d − IR,m; dotted-line) and ii) as calculated (d_d − IR,c; full-line) by balancing the ram pressures from the stellar wind (P_sw) and from the photoevaporated flow (P_pef). The expression for d_d − IR,c, that takes into account the logarithmic outward acceleration expected for the photoevaporated wind is given by (see Henney et al. 1996):

${\mathit{d}}_{\mathrm{d}\mathrm{-}\mathrm{IR}\mathit{,}\mathrm{c}}\mathrm{\simeq }{\mathit{R}}_{\mathrm{0}}\mathrm{·}[\mathrm{1}\mathrm{+}\mathrm{4}\ln (\frac{{\mathit{R}}_{\mathrm{0}}}{{\mathit{r}}_{\mathrm{d}}}){]}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}$(6)where r_d is the disk radius, and R₀: ${\mathit{R}}_{\mathrm{0}}\mathrm{=}\sqrt{\frac{\mathit{Ṁ}{\mathit{a}}_{\mathrm{II}}}{\mathrm{4}\mathit{\pi }{\mathit{P}}_{\mathrm{sw}}}}\mathrm{·}$(7)In this equation, a_II is the sound speed of the photoionized flow, and P_sw is the ram pressure of the stellar wind, which is P_sw = 8.3 × 10^-8 dyn cm^-2 for our chosen parameters (see Sect. 2). In Fig. 5, we see that both the measured and the calculated distances are comparable.

In Fig. 6, we show the xz-midplane density stratifications (left) and Hα maps (right) for models M2a (top), M2b (middle) and M2c (bottom), at t = 2000 years. The models have both EUV and FUV radiation fields (from the stellar source). The distance from the disk to the ionizing source is 0.2 pc, 0.1 pc and 0.02 pc for models M2a, M2b and M2c, respectively; see Table 1. The density of the disk is determined from Eqs. (3) and (4) for the FUV dominated models M2a and M2b, and Eqs. (3) and (5) for the EUV dominated model M2c. The disk axis is aligned with the x-direction (θ = 0° see Fig. 2 and also Table 1), and the stellar wind is present in the three models.

We see that a high density region shrouds the disk in models M2a and M2b. This region is related to the photodissociation region expected from the analytical models of Johnstone et al. (1998) and Störzer & Hollenbach (1999). It has a temperature T₂ ≃ 1000 K (see Eq. (1) and Table 1). This region is composed of C ii and H i. We also note that the environment (the left part of the computational domain) and the photoevaporated wind outside the PDR have a temperature of about 10⁴ K, indicating that they are photoionized. The photodissociation region separates the disk from the ionized photoevaporated wind. In model M2c (in which the disk is at 0.02 pc from the ionizing photon source) the photodissociation region is limited to a thin shell in contact with the surface of the disk.

The ionization front radius (as defined in Fig. 2) is found to grow from r_IF = 2.4r_d = 120 AU at t = 200 yr to r_IF = 3.8r_d = 190 AU at t = 2000 yr for model M2a, while for model M2b, r_IF = 1.6r_d = 80 AU (constant throughout the simulation). In models M2a and M2b, the outer radius of the photodissociation region (r_w, defined in Fig. 2) grows with time. This lateral expansion occurs at a velocity of  ≈ 0.1 km s^-1, at  ~ 1 / 10 of the  ~ 1 km s^-1 outwards flow velocity of the gas within the photodissociation region. This indicates that the ionization front in these models is pressure confined by the base of the photoevaporated flow. The ratio between the cross section radius r_w and the axial standoff distance r_IF (see Fig.  2) has values r_w / r_IF = 1.4 and 1.9 for models M2a and M2b (at t = 2000 yr), respectively. Interestingly, Johnstone et al. (1998) predict r_w = 1.3r_IF from their simplified, analytic model. Models M2a and M2b have detached, almost hemispherical ionization fronts. On the other hand, in model M2c the ionization front lies very close to the surface of the disk (with just one pixel of separation between them). It is interesting to note that observations of proplyd 170 − 337 (Johnstone et al. 1998), that is situated at 0.1 pc from θ¹ Ori, show that its IF radius is approximately 80 AU. From the analytical calculations of Johnstone et al. (1998) for FUV dominated flows, ${\mathit{r}}_{\mathrm{IF}\mathrm{14}}\mathrm{\simeq }\mathrm{0.1}\mathrm{(}{\mathit{ϵ}}^{\mathrm{2}}\mathit{/}{\mathit{f}}_{\mathrm{r}}{\mathit{S}}_{\mathrm{49}}\mathrm{)}{\mathit{d}}_{\mathrm{17}}^{\mathrm{2}}$, where r_IF14 is the radius of the IF in units of 10¹⁴ cm, d₁₇ is the distance from the source in units of 10¹⁷ cm, S₄₉ is the ionizing photon rate in units of 10⁴⁹ s^-1, f_r is the fraction of EUV photons absorbed by recombination in the ionized region of the flow, and ϵ is a dimensionless parameter of order unity (see Johnstone et al. 1998). Taking (ϵ² / f_rS₄₉) = 16 (Johnstone et al. 1998), and d₁₇ = 3.08, we obtain r_IF ≃ 100 AU, which is consistent with the values obtained from models M2a and M2b.

In Fig. 7 we show the velocity field in a region close to the accretion disk superimposed on density (left panels) and Mach number maps (right panels), M = v / a, where a is the local sound speed and v is the local velocity of the gas, for models M2a (top panels) and M2b (bottom panels), at t = 200 yr. From Fig. 7, we see that the photoevaporated flow has four different dynamical regimes, indicated in the top-right panel of Fig. 7 by the labels 1 to 4, which refer to: 1) a supersonic neutral outflow; 2) a subsonic, neutral outflow, both of them inside the PDR region and separated by a shock; 3) a supersonic, ionized outflow outside the PDR region and 4) a subsonic, neutral inflow, localized around the disk center (at r ≲ r_g).

As we have stated before in Sect. 2, the disk radius adopted in our models is larger than the gravitational radius. For instance, for the parameters of the M2 models (M2a, M2b and M2c), we have r_g ≃ 20 AU. We can see from the velocity field maps that between disk radial distances r = 0 (the disk axis) and r ≃  ± 3Δz (where Δz = 5.2 AU is the spacing of the highest resolution grid, see Sect. 2), we do not have an outflow. For r ≳ 20 AU, the outflow starts supersonically from the disk surface, with Mach numbers  ≈ 1.1 (region 1 in Fig. 7), and reach an internal shock (before the IF; see Fig. 7) with Mach numbers  ≈ 1.8 (note that the PDR has a sound velocity a_I = 3.2 km s^-1). Between this internal shock and the outer boundary of the PDR, the outflow is subsonic, with Mach numbers  ≈ 0.4 − 0.6 (region 2 in Fig. 7). Outside the PDR and before reaching the interaction region between the two winds, we have a supersonic outflow with Mach numbers  ≈ 1.2 − 1.6 (with a sound velocity a_II = 14.2 km s^-1).

The structure obtained from models M2a and M2b is similar to the scenario proposed by Johnstone et al. (1998) for the FUV dominated flow. From this we conclude that both the M2a and M2b models produce photoevaporated flows in the FUV dominated regime.

In Fig. 8 we show the xz-midplane density stratifications (left) and Hα maps (right) for models M3a (top) and M3b (bottom panels) at t = 2000 yr. The models have an EUV+FUV radiation field from the ionizing source, that is located at 0.2 pc (M3a) and 0.1 pc (M3b) distance from the disk. In both cases, the disks have a 45° of inclination with respect to the direction of the impinging photons (θ = 45° see Fig. 2 and Table 1). In Fig. 9 we present the same results, but for models M4a and M4b (with θ = 75°). As can be seen in Table 1, models M2a, M3a and M4a differ only by the values of the disk inclination. The same is true for the set of models M2b, M3b and M4b.

In the inclined disk models (models M3a, M3b, M4a and M4b, see Figs. 8, 9) the internal structure of the PDR region is the same as in the M2a and M2b models (with θ = 0°, see Fig. 7 and Table 1). There is a neutral supersonic outflow (with M ≈ 2 in all models) near the illuminated surface of the disk. For inner disk radial distances (from the disk axis), there is a gas inflow. In Fig. 10 we show velocity vectors superimposed on density maps for models M3a (top-left), M3b (top-right), M4a (bottom-left) and M4b (bottom-right). As in the M2a and M2b models, we again have an internal shock that separates the supersonic from the subsonic neutral flows, inside the PDR. Outside this internal shock we have the outer edge of the PDR region, that separates the subsonic, neutral flow from the external, ionized supersonic flow. These structures are characteristic of an FUV dominated flow, as discussed previously.

The Hα emission in the region close to the ionization front of models M3 and M4 (Figs. 8 and 9) has a general morphological resemblance to the Hα morphology of the “face-on” models (M2, see Fig. 6). However, the Hα maps of the inclined disk models show clear side-to-side asymmetries, which are more pronounced for the models with θ = 45° (M3a, b, see Table 1 and Fig. 8) than for the ones with θ = 75° (M4a, b, see Fig. 9).

From the t = 2000 yr time-frames of models M3 and M4 (see Figs. 8 and 9), we have calculated the axial stand-off distance (r_IF) and width of the cross section presented by the ionization front (r_w, see Fig. 2). We obtain the following results:

• –

  M3a: r_IF = 3.6r_d = 180 AU; r_w = 4.4r_d = 220 AU;

• M3b: r_IF = 2.0r_d = 100 AU; r_w = 2.4r_d = 120 AU;

• M4a: r_IF = 2.8r_d = 140 AU; r_w = 3.7r_d = 185 AU;

• M4b: r_IF = 1.9r_d = 95 AU; r_w = 2.1r_d = 105 AU,

with r_w / r_IF ≲ 1.3 for all models. These results illustrate the reduction in the size of the photodissociation region that is obtained for increasing values of the disk inclination angle θ.

In Fig. 11 we present the mass loss rate as a function of time, for all the simulated models except model M2d. The mass loss rate for models M1a and M1b (without an FUV radiation field; see Table 1) shows a small time variability (as mentioned above), but is almost independent of the presence of the stellar wind. Also, the mass loss rate has approximately the same value even if the distance from the disk to the ionizing photon source changes by a factor of two, even for the are FUV dominated models, as is the case for those with a given inclination angle (with respect to the ionizing photon source) as M2a and M2b (θ = 0°), M3a and M3b (θ = 45°), and M4a and M4b (θ = 75°). This is due to our simplified model, in which the PDR temperature is assumed to be position independent and with the same value for all models. Looking at a series of models with increasing values of the inclination angle θ (e.g., models M2a, M3a and M4a, which otherwise have identical paramenters, see Table 1), we see that the mass loss rate as a monotonically decreasing function of θ. This result is to be expected, given the lower effective cross sections (to the impinging UEV/FUV radiation field) presented by the disk for increasing inclination angles. In particular, if we take the mean value of Ṁ for models M2a and M2b (with θ = 0), for models M3a and M3b (with θ = 45°) and for models M4a and M4b (with θ = 75°), we find Ṁ_M2 = 1.4 × 10^-7   M_⊙ yr^-1, Ṁ_M3 = 6.2 × 10^-8   M_⊙ yr^-1Ṁ_M4 = 4.8 × 10^-8   M_⊙ yr^-1, which implies Ṁ_M3 ≃ 0.4Ṁ_M2 and Ṁ_M4 ≃ 0.3Ṁ_M2.

In Fig. 12, we show the xz-midplane density stratifications (left) and Hα maps (right) for model M2d at t = 200 yr. This model has the same parameters as the M2b model (θ = 0°, d = 0.1 pc; see Table 1) but a different PDR temperature T₂ = 3000 K (see Eq. (1)). The ionization front radius is larger, by a factor of  ≃ 3, for the model M2d when compared with those with T₂ = 1000 K. The mass loss rate is  ≈ 5 × 10^-7   M_⊙ yr^-1 for M2d (about 3.5Ṁ_M2). We have run equivalent simulations for models M3b and M4b, that is, a model with θ = 45° and d = 0.1 pc, but with T₂ = 3000 K, and a model with θ = 75° and d = 0.1 pc, but with T₂ = 3000 K. For these models (not shown here), we found the same trend of larger mass loss rates and ionization front standoff distances for larger values of T₂.

4. Conclusions

This work describes the first fully three-dimensional numerical gasdynamic simulations of the photoevaporation of neutral disks, subject to an interaction with the FUV/EUV photon flux and the wind from an external O star.

With this initial setup, we study the time-evolution of the photoevaporated flow. We have calculated models in the EUV dominated regime, in which one obtains an ionization front (IF) close to the surface of the disk, and a photodissociation region (PDR) not resolved by our numerical scheme between the surface of the disk and the IF, and in the FUV dominated regime, in which a PDR (well resolved in our numerical simulations) is formed between the surface of the disk and the IF. Both EUV and FUV dominated models develop an interaction region between the photoevaporated wind and the wind from the hot star. In some cases, this interaction region lies within our computational domain. For these cases, the interaction region lies close to the position of stellar wind/photoevaporated flow ram pressure balance.

We have simulated models of disks with different inclinations with respect to the direction of the impinging EUV and FUV photon fluxes, and at different distances from the ionizing photon sources ranging from 0.02 pc to 0.2 pc. The model at a distance of 0.02 pc from the O-star is clearly EUV dominated, since the IF is just above the disk surface. Models at distances from 0.1 pc to 0.2 pc show detached PDRs. These PDRs show structures, such as the presence of a neutral supersonic flow near the disk and, farther out but before reaching the IF, a subsonic neutral flow, that are consistent with the main features predicted by Johnstone et al. (1998) for a FUV-dominated flow. We note that these structures occur in our simulations even for the inclined disk models. The typical size of the IF was found to be in agreement with analytical results. In particular, the limits predicted by Johnstone et al. (1998) for the ionization front radius, 2.5r_d ≲ r_IF ≲ 4r_d, were recovered in most of the models, with or without disk inclination.

We found that, for a given inclination of the disk, the mass loss rate seems to be independent of the distance to the FUV/EUV photon source (for the FUV-dominated flows). In our simulations this is due to the constancy of the PDR temperature. For higher disk inclinations, we obtain smaller mass loss rates, which is a direct result of the smaller effective cross section presented by the disk to the impinging FUV/EUV photon flux. Taking the mean values of the mass loss rate from our simulations, we obtain Ṁ_M2 = 1.4 × 10^-7   M_⊙ yr^-1, Ṁ_M3 = 6.2 × 10^-8   M_⊙ yr^-1Ṁ_M4 = 4.8 × 10^-8   M_⊙ yr^-1, for models with θ = 0°, θ = 45° and θ = 75°, respectively. These values are surprisingly similar, in spite of the very different inclination angles used in the simulations.

Our models predict Hα emission maps with moderate side-to-side asymmetries, which directly result from the misalignment of the axis of the accretion disk with respect to the direction to the photoionizing source. Asymmetries are seen both in the emission associated with the IF and in the emission of the photoevaporated wind/expanding H II region interaction working surface. It is clear that asymmetries in these two structures are indeed observed in many of the Orion proplyds (Bally et al. 2000). A detailed attempt to model a specific, asymmetrical proplyd with a 3D photoionizing disk simulation will be presented in a future paper.

Acknowledgments

We thanks the anonymous referee for his/her comments and suggestions, which contributes to improve the quality of the manuscript. M.J.V. would like to thanks FAPESB for partial financial support (PPP project 7606/2006). A.H.C. wish to thank the Brazilian agency CNPq for partial financial support (projects 307036/2009-0 and 471254/2008-8), and also PROPP-UESC (projects 635 and 802). A.R. acknowledges support from the CONACyT grants 61547, 101356, and 101975. We thank A. Esquivel for sharing with us a set of IDL routines for plotting the YGUAZU-a output.

References

All Tables

All Figures