© ESO, 2011

1. Introduction

Dense cores within molecular clouds are birth places of stars as some of them are gravitationally bound and will subsequently collapse and form a star. Observations of starless, prestellar cores promote our understanding of the initial stages of star formation. Molecular line observations provide information about e.g. gas density, cloud structure, and kinematics of the gas. Furthermore, they can be utilised to probe gas kinetic temperature. Radiative transfer calculations suggest, that the external heating by the interstellar radiation field is attenuated by dust in dense cores thus causing the gas temperature to fall towards the centre of these cores (see, e.g. Zucconi et al. 2001; Evans et al. 2001; Keto & Field 2005). The temperature gradient is most pronounced within a radius of a few thousand AU from the core centre. Far infrared- and submillimeter dust continuum studies have shown dust temperature to decrease towards the centre of cores (see, e.g. Ward-Thompson et al. 2002; Juvela et al. 2002; Pagani et al. 2003, 2004; Schnee & Goodman 2005). High angular resolution observations showing this behaviour for gas have not been extensively conducted (one of the few is Crapsi et al. 2007).

An isolated, roundish starless core with high central density (n_c > 10⁵ cm^-3) is a preferred choice to look for the predicted temperature gradient. One nearby core that fulfills these criteria is Ophiuchus D (L1696A), north from the ρ-Ophiuchus molecular cloud. We adopt the new distance estimate of 120 pc to the ρ-Ophiuchus cloud for this core (Lombardi et al. 2008). There is no embedded protostar in southern part of Oph D (Kirk et al. 2007). Our reference position RA = 16^h28^m289, Dec =  − 24°19 ′09″(J2000), coincides with the dust emission maximum of the 850 μm SCUBA map of Kirk et al. (2005), the dust absorption maximum of 7 μm ISOCAM map of Bacmann et al. (2000), and the N₂D⁺ line emission peak according to the observation by Crapsi et al. (2005). Direct evidence of very cold gas in molecular clouds was obtained by Harju et al. (2008) when observing the ortho-H₂D⁺ line towards Oph D with APEX, using the same antenna pointing position as our reference point. They reported the kinetic temperature of the core derived from H₂D⁺ linewidth to be 6.0  ±  1.4 K.

Nitrogen bearing molecular species, such as NH₃, are found in high density cores (n_c > 10⁵ cm^-3) where most carbon bearing species have already depleted from gas phase (e.g., Willacy et al. 1998; Aikawa et al. 2001; Tafalla et al. 2002; Jørgensen et al. 2004). In addition to tracing the densest gas NH₃ is also a good thermometer of the gas (Ho et al. 1979), and therefore appropriate choice for probing the temperature gradient.

2. Observations and data reduction

Maps of the para-NH₃ (J, K) = (1, 1) and (2, 2) transitions at frequencies 23 694.4955 MHz and 23 722.6336 MHz were made with the Very Large Array (VLA) in five observing slots, each six hours in duration, between 9th and 15th June 2008 in the array’s DnC configuration. The synthesised beam of VLA in this hybrid compact configuration was 372  ×  235 (position angle PA = 372 east of north) with natural weighting applied. The FWHM of the primary beam is 2′ at this frequency. In order to improve the signal-to-noise ratio of the (2, 2) map we used a circular clean beam with FWHM of 10″ to construct the clean image. For consistency we also cleaned the (1, 1) map with the same beam to ensure identical resolution of the maps for the analysis. The transitions were observed simultaneously in 2IF mode with 128 channels and a bandwidth of 1.5625 MHz in both spectral windows (resulting in a channel width of 12.2 kHz or 0.154 km s^-1). The bandwidth was enough to include the main hyperfine group and the inner F₁ = 1−2 and F₁ = 2−1 satellite groups of the (1, 1) line.

Temporal variations in gain and phase were calibrated every 15 min using quasar 1626-298 (F_1.3   cm = 2.57 Jy¹). Bandpass- and absolute flux calibration were carried out in the beginning and the end of each 6 h period using quasar 1331+305 (F_1.3   cm = 2.59 Jy). The uncertainty of flux calibration was estimated to be of the order of 15%. At the time of the observation some of the antennas of the array were replaced by new Expanded VLA antennas, for which reason the online Doppler tracking was not available. In order to avoid abrupt phase jumps the sky frequency was calculated by hand for each phase calibrator – source – phase calibrator triplet.

The visibilities were calibrated and data reduction and imaging was carried out with Common Astronomy Software Applications (CASA)². The noise rms of the data cube cleaned with the 10″ beam is 1 mJy/beam, which corresponds to about 22 mK. We were unable to estimate the missing flux since no single dish data was available (the maximum in the full resolution clean map is 0.11 K km s^-1 whereas the negative contours reach a level of  − 0.06 K km s^-1). This should not affect the uncertainties of physical quantities derived from line intensity ratios, e.g. optical thickness and kinetic temperature, unless there is differential flux loss between (1, 1) and (2, 2) spectral lines.

We fitted the hyperfine structure of NH₃(1, 1)-line using our χ²-minimisation routine to obtain estimates of line total optical thickness τ₁₁, velocity 3_LSR, and linewidth (FWHM) Δ3. The hyperfine structure of (2, 2)-line could not be detected, so the main component was fitted with a single Gaussian line profile.

3. Results

3.1. Integrated intensity maps

The integrated intensity maps are shown in Figs. 1 and 2. The full angular resolution map of the ammonia (1, 1) transition is shown in Fig. 1. The signal-to-noise ratio of the (2, 2) transition was very small when cleaned to the full angular resolution, and therefore its map is not shown. In Fig. 2 the integrated intensity maps of (1, 1) and (2, 2) transitions smoothed to 10″ are presented. The emission of (1, 1) line was integrated over all available hyperfine components (main group and the inner satellite groups) and (2, 2) line over the main component. The RA and Dec offsets in these maps are relative to our reference position.

The integrated intensity maps show two prominent local emission maxima within the mapped area, especially in NH₃(1, 1) maps. Hereafter, we will refer to these as the northern and the southern peak. Such structure has not been seen in previous submillimeter or infrared observations nor in molecular line mappings of this core, with the exception of Motte et al. (1998) who report two cores in the southern end of Oph D in their 1.3 mm continuum map (D-MM1 and D-MM2, see their Fig. 2). Typically the southern end of Oph D has been characterised as a single prestellar core. The structure will be discussed in more detail in Sect. 4.

Fig. 1 Integrated intensity map of Oph D in (1, 1) transition of NH3. The VLA synthesised beam, 372  ×  235, is indicated in the bottom right corner. The contours are from 20% to 90% with steps of 10% of the map peak value 0.11 K km s-1. The reference position is RA = 16h28m289, Dec =  − 24°19 ′09″(J2000).

Fig. 2 Integrated intensity maps of Oph D in (1, 1) (left) and (2, 2) (right) transitions of NH3 smoothed to 10″ to improve the signal-to-noise ratio of (2, 2) line. The 10″ Gaussian clean beam is indicated in the bottom right corner. The contours are from 10% to 90% with steps of 10% of the map peak values 3.5 K km s-1 and 0.11 K km s-1 for (1, 1) and (2, 2) lines, respectively.

Fig. 3 Channel map of Oph D in the (1, 1) transition of NH3. Channel width is 0.154 km s-1.

3.2. Kinematics of the gas

The kinematics of the gas was studied with the aid of Gaussian fits to the hyperfine structure of the NH₃(1, 1) line which provide accurate estimates of radial velocity and line widths. The resulting 1σ uncertainties of these parameters are typically better than 0.01 km s^-1. The hyperfine fits were used to form another image cube with a single Gaussian line profile normalised with the integrated intensity of the full NH₃(1, 1) line. This image cube describes the radial velocity distribution of the ammonia. The velocity channel map based on this image cube is shown in Fig. 3. The channel width corresponds to the original resolution of the VLA correlator. The map covers velocities from 3.29 to 3.76 km s^-1 in four channels. The emission of the northern peak moves from northeast towards southwest as the velocity increases, whereas the southern peak emission appears to move slowly westwards along a roughly east-west path.

We created a local velocity gradient map for both transitions following the method presented in Goodman et al. (1993). Each velocity gradient is determined using a nine pixel grid centred at the pixel under consideration. Furthermore, we required that at least seven of the surrounding pixels have a usable velocity estimate available for the least squares fit. The results are presented in Fig. 4, where the velocity gradient vectors of (1, 1) transition are superimposed on the map of integrated intensity.

In the velocity gradient maps one can see two distinctly different velocity fields on top of the southern and the northern peaks, which seems to suggest dissimilar rotation of gas near these peaks. The effect can also be seen in the map of the (2, 2) transition, although not as prominently. The velocity gradients in the vicinity of both peaks seem homogeneous enough to apply the solid-body rotation fit also in a larger scale. The overall total gradient is determined for the two locally uniform areas separately instead of giving one value to the whole clump. The magnitude and direction of total gradients for northern and southern peaks are given in Table 1.

The magnitudes of the total velocity gradients, and  km s^-1 pc^-1, of the northern and the southern peak are comparable to the values typically detected in other cores using nitrogen compound tracers (Caselli et al. 2002a,b; Crapsi et al. 2004, 2007), but the agreement with previous observations of Oph D is poor (Goodman et al. 1993; Crapsi et al. 2005). This is probably due to significant difference in the angular resolution of the data, being  ~ 80″ (Benson & Myers 1989, see their Fig. 21) for NH₃ data used in Goodman et al. (1993) and 26″ for the N₂H⁺(1 − 0) data used in Crapsi et al. (2005). With these beamwidths the velocity field is averaged considerably, and more importantly the reported magnitude of the total gradient of the whole core will not be a good figure of rotation for this object, since, as we have shown, there exists two distinct velocity regions with almost opposite gradient directions.

Fig. 4 Velocity gradient map of NH3(1, 1) transition of Oph D. The arrows depict the direction (of increasing radial velocity and the magnitude of the velocity gradients. An arrow with a length of 10″ represents a gradient of 14.3 km s-1 pc-1.

Fig. 5 Spectra (gray bars) and chi-square fits (black solid line) of NH3(1, 1) and (2, 2) from the centre of the northern core (from the 10″ resolution map). The (2, 2) spectrum is offset by  −1 K for clarity. The splitting of the (1, 1) components is due to hyperfine structure, and is not to be taken as a sign of self-absorption.

3.3. Ammonia column densities and the kinetic temperatures using the LTE assumption

We fitted the hyperfine structure of NH₃(1, 1) line using a χ² routine, and the (2, 2) line with a simple Gaussian profile. Fits were performed for pixels having signal larger than 3σ_noise level. We made the usual assumption, that the different hyperfine states are populated according to LTE. We derived the excitation temperature T_ex, rotational temperature T₁₂ and the total (ortho + para) ammonia column density N(NH₃) following the procedure described by Harju et al. (1993). In the calculation of total ammonia column density we assumed the NH₃ ortho-to-para ratio to follow the LTE conditions. The kinetic temperature of the gas was calculated from the rotational temperature using the semi-empirically derived analytic expression from Tafalla et al. (2004)${\mathit{T}}_{\mathrm{kin}}\mathrm{=}\frac{{\mathit{T}}_{\mathrm{12}}}{\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{T}}_{\mathrm{12}}\mathit{/}\mathrm{42}\mathrm{)}\ln \left[\mathrm{1}\mathrm{+}\mathrm{1.1}\exp \mathrm{(}\mathrm{-}\mathrm{16}\mathit{/}{\mathit{T}}_{\mathrm{12}}\mathrm{)}\right]}\mathit{,}$(1)which is based on Monte Carlo simulations using NH₃  −  H₂ collision coefficients from Danby et al. (1988) and is found to be suitable for dense cores having temperatures less than 20 K. According to Maret et al. (2009), who present new collisional rates for NH₃  −  H₂, their new and Danby’s coefficients agree within a factor of two. They also found, that non-LTE Monte Carlo radiative transfer simulations using their new and Danby’s coefficients produce practically equal results for a parameter set typical of a dense and cold prestellar core. Thus, Eq. (1) characterises the relation between rotational and kinetic temperature of the gas robustly.

The results of the hyperfine fits and the physical parameters derived for the northern peak (RA = 16^h28^m293, Dec =  − 24°19′03″) and the southern peak (RA = 16^h28^m28 8, Dec =  − 24°19′23″) peaks are presented in Table 2. The maps of NH₃ optical thickness τ and column density N(NH₃) are shown in Fig. 6.

3.4. Mass and fractional abundance

In the absence of continuum data of our own, we use the 1.3 mm data of Motte et al. (1998) to estimate masses and fractional abundance of ammonia of our sources. Despite the difference of 5″ between the coordinates of our core centres and the D-MM1 and D-MM2 reported by Motte et al. (1998), we consider them the same cores.

Fig. 6 The map of NH3 optical thickness τ (left) and total column density N(NH3) [cm-2] (right) of Oph D.

If we assume the dust is isothermal (having temperature T_d) and the background emission has been subtracted from the observed intensity, the intensity of dust emission ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{dust}}$ can be calculated using the equation of radiative transfer: ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{dust}}\mathrm{=}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}}\mathrm{)}\mathit{.}$(2)Thermal dust emission in the mm/submm region is optically thin, so Eq. (2) becomes ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{dust}}\mathrm{\approx }{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}{\mathit{\tau }}_{\mathit{\nu }}$. The optical thickness τ_ν is given by the product of total surface density Σ, the mass absorption coefficient per unit mass of dust κ_d and the dust-to-gas mass ratio R_d. The surface density can be expressed with the aid of column density Σ = N(H₂)μ_H2m_H, where μ_H2 is the mean molecular weight per H₂ molecule and m_H is the mass of the hydrogen atom. The column density as a function of dust intensity and temperature becomes $\mathit{N}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{I}}_{\mathit{\nu }}^{\mathrm{dust}}}{{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}{\mathit{\mu }}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{m}}_{\mathrm{H}}{\mathit{\kappa }}_{\mathrm{d}}{\mathit{R}}_{\mathrm{d}}}\mathrm{·}$(3)The 1.3 mm flux density in a 15″ beam reported by Motte et al. (1998) is 45  ±  5 and 45  ±  15 mJy beam^-1 for southern and northern core, respectively. We assume that the gas and dust temperatures are equal, and use an average value of T_d = 8 K for the dust temperature for both cores. Adopting values of κ_d = 0.5 cm² g^-1 (Ossenkopf & Henning 1994) and R_d = 0.01, the molecular hydrogen column density for the cores becomes $\begin{array}{ccc}\mathit{N}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\mathrm{5.3}\mathrm{\pm }\mathrm{2.0}\mathrm{\times }{\mathrm{10}}^{\mathrm{22}}{\mathrm{cm}}^{-2}\mathrm{forthenorthernand}& & \end{array}$$\begin{array}{ccc}\mathit{N}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\mathrm{5.3}\mathrm{\pm }\mathrm{0.6}\mathrm{\times }{\mathrm{10}}^{\mathrm{22}}{\mathrm{cm}}^{-2}\mathrm{forthesoutherncore}\mathit{.}& & \end{array}$The fractional abundance of ammonia can be estimated using the calculated value of molecular hydrogen column density and our ammonia column density (Table 2), $\begin{array}{ccc}\mathit{X}\mathrm{\left(}{\mathrm{NH}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathit{N}\mathrm{\left(}{\mathrm{NH}}_{\mathrm{3}}\mathrm{\right)}}{\mathit{N}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}}\mathrm{\approx }\mathrm{1.8}\mathrm{\times }{\mathrm{10}}^{-8}\mathit{.}& & \end{array}$The uncertainty of factor  ~ 2 is typically associated with the dust opacity coefficient. This affects directly the column density, mass, and abundance estimates.

Motte et al. (1998) give mass estimates for both cores based on integrated flux densities. The reported masses of the cores are 0.15 M_⊙ and 0.16 M_⊙ for the southern and northern cores. Based on peak flux densities, the mass within the 15″ beam can be estimated with $\mathit{M}\mathrm{=}\frac{{\mathit{S}}_{\mathit{\nu }}{\mathit{d}}^{\mathrm{2}}}{{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}{\mathit{\kappa }}_{\mathrm{d}}{\mathit{R}}_{\mathrm{d}}}\mathit{,}$(4)where d = 120 pc is the distance to the source. This gives a mass estimate of M = 0.10   M_⊙ for both cores. The 15″ beam, however, encloses the cores only partially, and therefore the mass estimate is probably too small to describe the whole cores.

4. Modelling

4.1. Description of the model: BE spheres

The southern tip of Oph D cloud has in general been considered a single prestellar core. Our VLA data, however, indicates the existence of two cores instead of one. This hypothesis is backed up by the velocity gradient map (Fig. 4) as well as by the two separate local maxima of optical thickness and column density seen in Fig. 6. The velocity gradients seem to imply, that the direction and magnitude of rotation of the gas within the area of these two hypothetical cores is markedly different from each other. The secondary peaks of optical thickness and column density south from the most dominant peak also point towards the presence of two cores. Steinacker et al. (2005) modelled Oph D cloud with a complex 3D substructure based on ISOCAM 7 and 15 μm maps and a 1.3 mm map from IRAM 30 m telescope. In their model the southern part of the cloud also contains two dominant density concentrations. Also Motte et al. (1998) report two cores (D-MM1 and D-MM2 in their Fig. 2 and Table 2) in the southern part of Oph D.

We attempt to model the southern part of Oph D with two cores, each in hydrostatic equilibrium and with a radial temperature gradient, the so called modified Bonnor-Ebert spheres (Evans et al. 2001; Zucconi et al. 2001) with radii 2200 AU (northern core) and 2100 AU (southern core). The central positions of the core models were chosen to coincide with the local kinetic temperature minima derived from the observation. The calculation of the gas density and temperature profiles for both core models was done iteratively until satisfactory convergence was achieved. The initial core model was calculated by Sipilä et al. (2010, see Sect. 2.1. therein), who starting from an isothermal sphere derived the dust density and temperature profile for a core with a radius of 2400 AU. The southern core density is high enough (n_c > 10⁶ cm^-3) to assume the dust and gas temperatures to be equal (Burke & Hollenbach 1983), but for the northern core (n_c ~ 10⁵ cm^-3) the situation may not be as clear. Despite this uncertainty we used Sipilä’s dust temperature and density profile, scaled to the radii of our core models, as a first estimate of our cores.

In order to derive the gas density profile the following modification to the original substitutions introduced by Bonnor (1956) was made to include the non-isothermal temperature distribution T(r) of the core in the calculations (Sipilä et al. 2010, 2011): $\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{\lambda }\left(\frac{{\mathit{T}}_{\mathrm{c}}}{\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\right)\mathrm{e}\mathrm{-}\mathit{\psi }\mathit{,}$(5)$\mathit{r}\mathrm{=}\mathit{\xi }\sqrt{\frac{\mathit{k}{\mathit{T}}_{\mathrm{c}}}{\mathrm{4}\mathit{\pi Gm\lambda }}}\mathit{,}$(6)where T_c and λ are the temperature and density at the centre. The resulting equation of hydrostatic equilibrium was solved numerically using boundary conditions presented in Sipilä et al. (2010), and the obtained density and temperature profiles were fed to the radiative transfer Monte Carlo code (Juvela 1997) to simulate the spectra. We use molecular data provided by the Leiden Atomic and Molecular Database (Schöier et al. 2005). The method described in Keto et al. (2004) is used for handling overlapping hyperfine emission lines. For the simulation we constructed a 3D model which consists of these two cores with a projected separation of 2600 AU in the plane of the sky. In addition to the density and temperature, the fractional abundance of ammonia X(NH₃) in our model can be adjusted, as well as the turbulent velocity within the core. The non-thermal velocity dispersion component Δ3_NT is calculated from the observed linewidths Δ3_obs with an equation: $\mathrm{\Delta }{3}_{\mathrm{NT}}\mathrm{=}\sqrt{\mathrm{(}\mathrm{\Delta }{3}_{\mathrm{obs}}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{8}\ln \mathrm{2}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{kin}}}{{\mathit{\mu }}_{\mathrm{NH}\mathrm{3}}{\mathit{m}}_{\mathrm{H}}}}\mathrm{·}$(7)The prevailing non-thermal velocity dispersion in the northern and southern cores is between 40 − 50 m s^-1 (the isothermal sound speed for a core with average gas kinetic temperature of 8 K is c_s  ~  170 m s^-1). The turbulent velocities in our simulation model are therefore also kept at a low value.

Fig. 7 Line-of-sight average profiles of column density, kinetic temperature and optical thickness (τtot) towards three radial directions (E of N) from northern core centre outwards. The values derived from the observation are marked with small gray circles and from the model with a solid line.

4.1.1. Comparison between the model and observation

In dense prestellar cores the fractional abundance of NH₃ may increase towards the core centre (see e.g. Hotzel et al. 2004, and references therein). In more evolved cores, however, when the central density increases above n_c ~ 10⁶ cm^-3, also ammonia starts to vanish from the gas phase (Aikawa et al. 2005; Flower et al. 2006; Miettinen et al. 2010). A variety of fractional abundance values was tested in core models, using the estimate X(NH₃) ≈ 1.5  ×  10^-8 discussed in Sect. 3.4 as a starting point. In each simulation the abundance radial profile shape followed the current density profile shape, except in the case of depletion which approximated as a steep descending slope.

The best-fit parameters of the northern core turned it into a nearly isothermal sphere; the temperature at the centre is 8.8 K and increases only 0.4 K towards the outer edge. The optical thickness and column density could not be fitted simultaneously to match the observed values at the centre of the cores, and a relative error of  ~15% remains. The central density of the model is n_c = 2  ×  10⁵ cm^-3 and the fractional abundance gradually increases towards the centre peaking at the value X(NH₃) = 2  ×  10^-8. Mass of the northern core model is only 0.045 M_⊙, and could not be increased while maintaining a reasonable agreement with the observation.

In order to model the low temperatures observed in the southern core the model central density had to be increased to 4  ×  10⁶ cm^-3, similar to the results of Harju et al. (2008). The best-fit temperature profile has a temperature decrease from 10 K to 6.1 K from the surface to the centre. The fractional abundance of ammonia was then used to control the optical thickness and column density. It was quickly noticed that a swift decrease in abundance within the innermost 700 AU was required to obtain a match with the observed values. The abundance increases to a value of X(NH₃) = 2.5  ×  10^-9 before the depletion starts to take effect. The mass of the southern core model is 0.16 M_⊙.

The structure of the southern part of Oph D is certainly more complex than what can be modelled using two spherical cores in hydrostatic equilibrium. Nevertheless, a relatively close agreement with the observation was achieved with this simple model. The spectra produced by the radiative transfer calculations were analysed in the same way as the observed spectra. The results are presented in Figs. 7 and 8. In both figures the column density, kinetic temperature, and optical thickness are plotted in three radial directions from the core centre towards the edge. Each subplot contains the observed values in a sector  ± 20° from the mean direction plotted with small gray circles together with their 1σ uncertainties. The corresponding simulation result is shown with a solid line. The temperature and density of northern and southern core models are presented in Fig. 9.

Fig. 8 Line-of-sight average profiles of column density, kinetic temperature and optical thickness (τtot) towards three radial directions (E of N) from southern core centre outwards. The values derived from the observation are marked with small gray circles and from the model with a solid line.

The 3D model was also used to reproduce the H₂D⁺ line observed by Harju et al. (2008) with APEX. The fractional abundance of H₂D⁺ was taken from Sipilä et al. (2010) and applied for both cores. Otherwise the model was exactly the same as in NH₃ simulations. A single pointing spectrum at our and Harju’s reference position using APEX antenna beam FWHM of 17″ was simulated, and a linewidth of Δ3 = 0.35 km s^-1 was obtained. This is exactly the same linewidth produced by a single hydrostatic core model of Harju et al. (2008) and Sipilä et al. (2010) in their respective simulations of this observation. Although the observed narrow linewidth Δ3 = 0.26 km s^-1 could not be reproduced with our two core model, we can at least see that the simulation is in line with the previous single core models, and did not introduce further deviation from the observed value.

Fig. 9 The density (solid line) and temperature (dash line) profiles of the southern (left) and northern core (right) BE model.

4.2. Description of an alternative model: spheroid cores

The modified Bonnor-Ebert sphere models reproduce the observed spectra reasonably well, and is also an apparent choice for the roundish southern core. The shape of the northern core, however, appears more elongated in the VLA map. An oblate spheroid, whose axis of rotation is on a plane perpendicular to our line of sight, would describe the shape and velocity field of the northern core perhaps better than a sphere. To test what effect the rotation of the cloud has on its shape, we run simulations with the Pencil Code³ using a 3D-grid. Pencil Code is a modular solver of magnetohydrodynamic equations using a high order finite difference scheme. The solution method enables a straightforward combination of various “physics modules”, like self-gravity, and various equations of state to the general equations describing the time-evolution of a hydrodynamic flow. As an initial condition, we used an isothermal Bonnor-Ebert sphere, where an initial rotation rate corresponded to the angular velocity of the northern core Ω = 2.2  ×  10^-13 rad s^-1. No extra rotational energy was put into the system during the time development of the simulation, and magnetic fields were neglected. The rotation rate t₀Ω = 0.27, where the dynamical timescale for BE sphere is t₀ = (4πGρ_c) ^− 1/2 (Banerjee et al. 2004), is quite fast, but still not enough to play alone a significant role in supporting the core against gravitational collapse. The dynamics of the model over time were observed, and with this angular velocity, according to calculations, the core seems to settle in an spheroidical equilibrium state.

The experiment indicated that the Bonnor-Ebert sphere transforms into a spheroid (with an axial ratio of the projected ellipse less than 2:1). The convergence in the simulation is slow, however, and still after 1 Myr there are oscillations. At this point the central density has decreased to one fifth of the original value, and the density profile is somewhat constant along the longer axis until a rather sudden drop occurs at the edge of the core. In the direction of the axis of rotation the shape of the original hydrostatic equilibrium density distribution remains about the same, just rescaled.

In this alternative model a slightly spheroidical shape was adopted for the southern core as well. The density and temperature profiles of the spheroids are, as was indicated by the simulation conducted with the rotating sphere, radially scaled versions of those of modified Bonnor-Ebert spheres. In this model the central densities were n_c = 1.0  ×  10⁶ cm^-3 and n_c = 1.4  ×  10⁶ cm^-3 for the northern and southern cores. The fractional abundance in both cores has a decrease in the centre, i.e. ammonia is depleted. The northern core abundance reaches a peak value of X(NH₃) = 1.1  ×  10^-8 1000 AU from the centre and the southern core abundance X(NH₃) = 7  ×  10^-9 700 AU from the centre.

The best-fit temperature profile of the northern core is steeper than in the previously discussed Bonnor-Ebert sphere model. The model gas temperature decreases to 7.4 K in the centre of the core, and rises to 10.2 K at the surface. The southern core temperature profile has a minimum at the centre of 6.4 K and rises to 8.0 K towards the surface. Masses of the models are 0.11 M_⊙ and 0.07 M_⊙ for northern and southern core, respectively.

4.3. Comparison between the two models

The simulation results of the two spheroid model show that the temperature profile of the northern core is not in good agreement with the observed values. The column density and optical thickness, on the other hand, are slightly better matched than with the modified Bonnor-Ebert sphere model. The central depletion of ammonia and higher gas density seem to provide more promising match than the less dense core with incresing abundance towards the centre.

Both models of the southern core have problems with optical thickness and temperature profiles. The modified Bonnor-Ebert sphere model seems to do better with the temperature whereas with the optical thickness profile both models fail to explain the observed values consistently.

We estimated the goodness of fit by calculating the reduced χ² value between the model profile and the data: ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathit{D}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{(}{\mathit{O}}_{\mathit{i}}\mathrm{-}{\mathit{E}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}\mathit{,}$(8)where D is the number of degrees of freedom, O the observed value, E the modelled value and σ² the variance of observation. The χ²-values are presented in Table 3. The comparison between the two models reveals no clear winner; both models have strong and weak properties. For the southern core, however, the use of a modified Bonnor-Ebert sphere leads to better results.

It should be noted, however, that models with only one core, be it a spheroid or modified Bonnor-Ebert sphere, did not produce acceptable match with the observed values. Also, if we dismiss the temperature gradient of the southern core and make it isothermal, the temperature estimated from the simulated ammonia spectra is severely different from the observation. This indicates that a temperature gradient is necessary in the physical model of the southern core. The northern core is almost isothermal, and its best-fit radial temperature profile has only 0.4 K difference between centre and surface.

5. Discussion

As was seen in the previous section, our distinctly different models both produce a reasonable match with the observation, but none of the two can be clearly pointed superior. The model with two hydrostatic cores, modified Bonnor-Ebert spheres, is, however, on more solid ground theoretically and we will commit ourselves to that model in the following discussion.

A Bonnor-Ebert sphere satisfies the following virial equilibrium equation $\mathrm{2}\mathit{T}\mathrm{+}\mathit{U}\mathrm{=}\mathrm{3}\mathit{V}{\mathit{P}}_{\mathrm{ext}}\mathit{,}$(9)where P_ext is the external pressure, and T and U the kinetic and gravitational potential energies. We neglect the effects of possible magnetic field, and calculate the kinetic energy $\mathit{T}\mathrm{=}{\mathit{E}}_{\mathrm{therm}}\mathrm{+}{\mathit{E}}_{\mathrm{rot}}\mathrm{+}{\mathit{E}}_{\mathrm{turb}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\mathit{M}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{kin}}}{\mathit{\mu }{\mathit{m}}_{\mathrm{H}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{I}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\mathit{M}\frac{\mathrm{(}\mathrm{\Delta }{3}_{\mathrm{NT}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{8}\ln \mathrm{2}}$(10)and gravitational potential energy $\mathit{U}\mathrm{=}{\mathit{E}}_{\mathrm{grav}}\mathrm{=}\frac{\mathrm{-}\mathit{G}{\mathit{M}}^{\mathrm{2}}}{\mathit{R}}\mathit{f}\mathrm{\left(}{\mathit{\xi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{-}\mathit{G}{\mathit{M}}^{\mathrm{2}}}{\mathit{R}}\left[\frac{\mathrm{3}}{{\mathit{\xi }}_{\mathrm{0}}\frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}}\mathrm{-}\frac{{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}}{{\left(\frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}\right)}^{\mathrm{2}}}\right]$(11)of the cores to examine their relations. From the relation (9) we can estimate the external pressure assuming that the cores are Bonnor-Ebert spheres. The centrally peaked density profile was taken into account individually for both cores in the numerical calculation of moment of inertia, I, in the equation of E_rot. The function f(ξ₀) takes the density profile into account in gravitational potential energy for an isothermal core. We will use it although our cores have a temperature gradient (the slight difference between the isothermal and non-isothermal BE spheres is illustrated in Appendix B of Sipilä et al. 2011). The value of f(ξ₀) for a homogeneous density core is 0.6. The southern one of our core models is slightly supercritical ξ₀ = 6.75 (f(ξ₀) = 0.74) and the northern one has a markedly lower value ξ₀ = 1.49 (f(ξ₀) = 0.61). The energy ratios are presented in Table 4.

The ratio of thermal energy to gravitational energy of the northern core is 3.89. This would suggest, that the external pressure inhibits the core from dispersing. On the other hand, when the core mass is estimated using data partly from the observation and partly from the model the energy ratio is  ~ 0.5. The core mass in this case is estimated using the column densities within the area of projected ellipse in the interferometric map. When the average fractional abundance is estimated from the model Bonnor-Ebert sphere (X(NH₃) = 1.8  ×  10^-8 for the northern and X(NH₃) = 5  ×  10^-9 for the southern core), the total mass of the core can be calculated. The difference between this mass and the model mass for the northern core suggests that the spheroidal ellipsoid into which an initially rotating BE sphere evolves to, might be a more suitable core model. However, the mass of the spheroid model was also very low, only 0.11 M_⊙.

The modified BE sphere mass of the southern core is almost equal to the mass derived from the observations as described above. The southern core seems to be gravitationally bound, which can be seen from the ratio of kinetic and gravitational energy in Table 4. The mass of the core is very close to the estimation given by Motte et al. (1998), 0.15 M_⊙. Also the peak flux density of Motte’s 1.3 mm continuum data was equal for both cores. Due to their vicinity the external pressure for both cores should be roughly the same. The values are close when the northern core mass is the one derived from the observation. Using the model mass the difference is significant (see Table 4). These facts further support the conclusion, that the northern core model is not correct, despite the fact that the match with the observed quantities is satisfactory.

6. Summary and conclusions

We have presented an interferometric observation of the southern part of Ophiuchus D prestellar core. We mapped the core with VLA using the NH₃(1, 1) and (2, 2) inversion lines. The synthetic beam of VLA allowed spatial resolution in the order of 400 AU, but we were forced to decrease it to about 1000 AU in order to improve the (2, 2) line signal-to-noise ratio. From the observed spectra we derived the gas kinetic temperature distribution within the mapped area. The total column density of ammonia, the optical thickness of both lines, and the radial velocity distributions were also obtained.

Our observation shows that the southern part of Oph D comprises of two potentially prestellar cores. In the southern core the line-of-sight average gas temperature derived from the observed ammonia spectra drops from the outer edge towards the centre, reaching a minimum of T_kin = 7.4 K. Towards the centre of the northern core the line-of-sight average gas temperature from ammonia is T_kin = 8.9 K. The northern core is almost isothermal with only 0.4 K increase in temperature from centre to the edge. The observed masses of the cores, 0.21 M_⊙ for the northern and 0.17 M_⊙ for the southern core, are quite small. Their potential collapse could lead to the formation of brown dwarfs or low mass stars. The two cores were previously detected by Motte et al. (1998).

Our observation shows that there exists two rather homogeneous velocity gradient fields within the mapped area. They coincide well with the positions of the northern and southern cores. If we interpret the velocity gradients as rotational motion, we see that the two cores are rotating in almost opposite directions. The angular velocity of the northern core is Ω = 2.2  ×  10^-13 rad s^-1 and the southern Ω = 6.9  ×  10^-12 rad s^-1. Steinacker et al. (2005) proposed that Oph D would have been formed by turbulent fragmentation. The directions and magnitudes of the gas velocity field derived from our observation support this view.

A modified Bonnor-Ebert sphere approximates the southern core reasonably well. The density and temperature profiles are reproduced satisfactorily when the central number density of the model is n_c = 4  ×  10⁶ cm^-3 and the model temperature drops from 10 K at the surface to 6.1 K in the centre. Heavy ammonia depletion within 700 AU from the centre is required to succesfully match the southern core model with the observation.

The northern core is more challenging to model, and we did not find a single good solution. Two different parameter value combinations succeed in producing an acceptable match with the observation. The other one is a modified Bonnor-Ebert sphere with central density of n_c = 2 × 10⁵ cm^-3 and having ammonia abundance increase inwards, reaching X(NH₃) = 2 × 10^-8 at the centre. An almost equally good match with the observation can be achieved with a more dense spheroidical core with n_c = 1 × 10⁶ cm^-3 and heavy ammonia depletion within the innermost 1000 AU. Gas kinetic temperatures at the centres of the two core models are 8.8 K and 7.4 K, respectively.

The location of the cold spot in the southern part of Oph D according to Harju et al. (2008) is almost in between our two cores. They observed the ground state (1₁₀ → 1₁₁) rotational transition of ortho-H₂D⁺, which may peak at a different position than ammonia owing to chemistry. Nevertheless, the previous H₂D⁺ result is surprising as one might expect to find the coldest gas to concentrate in the core centres.

The minimum line-of-sight average gas kinetic temperature derived from the interferometric NH₃ observations is T_kin = 7.4  ±  1.1 K. The difference between this result and the temperature reported by Harju et al. (2008) (T_kin = 6.0  ±  1.4 K) is not statistically significant, and can furthermore be understood by considering the following points. First of all, H₂D⁺ probably traces denser, and thus colder gas than ammonia. Furthermore, our temperature represents the line-of-sight average, which inevitably increases the temperature estimate as ammonia thrives also in warmer envelope of the core. The fact that we could not utilise the full resolution interferometric map, but were forced to use a map cleaned with a 10″ restoring beam in our analysis, also influenced our results.

The cores seem to be gravitationally bound. The observed angular velocities of both cores, especially the northern one, are significant. If there is enough time for the cores to evolve, the rotation may lead to a formation of a disk. The southern core simulation results suggest that the core may be eventually collapsing, as the dimensionless radius of the Bonnor-Ebert sphere is ξ = 6.75, slightly above the critical value. There is no observational data to aid in the discussion regarding the role of magnetic field in the support of these cores. Some mechanisms, like a slight increase in the external pressure and, in the case of magnetic support, ambipolar diffusion is needed to trigger the collapse process.

Acknowledgments

The study has been funded by the Academy of Finland through grants 132291 and 127015. We thank all the VLA staff, especially Gustav van Moorsel, for their hospitality and support during our stay at the VLA. We are grateful to the anonymous referee for helpful comments. We also thank the editor, Malcolm Walmsley, for his useful comments. This research has made use of NASA’s Astrophysics Data System and the NASA/IPAC Infrared Science Archive, which is operated by the JPL, California Institute of Technology, under contract with the NASA.

References

All Tables

All Figures

	Fig. 1 Integrated intensity map of Oph D in (1, 1) transition of NH3. The VLA synthesised beam, 372  ×  235, is indicated in the bottom right corner. The contours are from 20% to 90% with steps of 10% of the map peak value 0.11 K km s-1. The reference position is RA = 16h28m289, Dec =  − 24°19 ′09″(J2000).
In the text

	Fig. 2 Integrated intensity maps of Oph D in (1, 1) (left) and (2, 2) (right) transitions of NH3 smoothed to 10″ to improve the signal-to-noise ratio of (2, 2) line. The 10″ Gaussian clean beam is indicated in the bottom right corner. The contours are from 10% to 90% with steps of 10% of the map peak values 3.5 K km s-1 and 0.11 K km s-1 for (1, 1) and (2, 2) lines, respectively.
In the text

	Fig. 3 Channel map of Oph D in the (1, 1) transition of NH3. Channel width is 0.154 km s-1.
In the text

	Fig. 4 Velocity gradient map of NH3(1, 1) transition of Oph D. The arrows depict the direction (of increasing radial velocity and the magnitude of the velocity gradients. An arrow with a length of 10″ represents a gradient of 14.3 km s-1 pc-1.
In the text

	Fig. 5 Spectra (gray bars) and chi-square fits (black solid line) of NH3(1, 1) and (2, 2) from the centre of the northern core (from the 10″ resolution map). The (2, 2) spectrum is offset by  −1 K for clarity. The splitting of the (1, 1) components is due to hyperfine structure, and is not to be taken as a sign of self-absorption.
In the text

	Fig. 6 The map of NH3 optical thickness τ (left) and total column density N(NH3) [cm-2] (right) of Oph D.
In the text

	Fig. 7 Line-of-sight average profiles of column density, kinetic temperature and optical thickness (τtot) towards three radial directions (E of N) from northern core centre outwards. The values derived from the observation are marked with small gray circles and from the model with a solid line.
In the text

	Fig. 8 Line-of-sight average profiles of column density, kinetic temperature and optical thickness (τtot) towards three radial directions (E of N) from southern core centre outwards. The values derived from the observation are marked with small gray circles and from the model with a solid line.
In the text

	Fig. 9 The density (solid line) and temperature (dash line) profiles of the southern (left) and northern core (right) BE model.
In the text