The central 1000 AU of a pre-stellar core revealed with ALMA. I. 1.3 mm continuum observations

TL;DR

This study uses ALMA observations to reveal the detailed structure of the central region of a pre-stellar core, showing a smooth, elongated kernel with a flat central zone, advancing understanding of early star formation processes.

Contribution

First high-resolution ALMA imaging of a pre-stellar core's central 1000 AU, revealing a smooth, elongated kernel consistent with magneto-hydrodynamical models.

Findings

Identified a compact, elongated kernel with a flat central zone in the pre-stellar core.

Demonstrated that the kernel's substructure can be explained by a smooth density profile.

Showed that lower-density cores are filtered out in high-resolution ALMA observations.

Abstract

Stars like our Sun form in self-gravitating dense and cold structures within interstellar clouds, called pre-stellar cores. Although much is known about the physical structure of dense clouds just before and soon after the switch-on of a protostar, the central few thousand astronomical units (au) of pre-stellar cores are unexplored. It is within these central regions that stellar systems assemble and fragmentation may take place, with the consequent formation of binaries and multiple systems. We present ALMA Band 6 observations (ACA and 12m array) of the dust continuum emission of the 8 Msun pre-stellar core L1544, with angular resolution of 2'' x 1.6'' (linear resolution 270 au x 216 au). Within the primary beam, a compact region of 0.1 Msun, which we call a "kernel", has been unveiled. The kernel is elongated, with a central flat zone with radius Rker ~ 10'' (~ 1400 au). The average number density within Rker is ~1 x 10^6 cm^{-3}, with possible local density enhancements. The region within Rker appears to have fragmented, but detailed analysis shows that similar substructure can be reproduced by synthetic interferometric observations of a smooth centrally concentrated dense core with a similar central flat zone. The presence of a smooth kernel within a dense core is in agreement with non-ideal magneto-hydro-dynamical simulations of a contracting cloud core with a peak number density of 1 x 10^7 cm^{-3}. Dense cores with lower central densities are completely filtered out when simulated 12m-array observations are carried out. These observations demonstrate that the kernel of dynamically evolved dense cores can be investigated at high angular resolution with ALMA.

Equations

$$M_{\ker }=\frac{D^2{S}_{1.3mm}}{{\kappa }_{1.3mm}{B}_{1.3mm}(T_d)},$$

$$L_J=\sqrt{\frac{\pi c_s^2}{G{\rho }_0}},$$

$$c_s=\sqrt{\frac{k_{\mathrm{B}}T}{{\mu }_{\mathrm{p}}{\mathrm{m}}_{\mathrm{H}}}},$$

$$F(r)=\frac{F_0}{1+(r/R_{\mathrm{flat}}{)}^{\alpha }},$$

$$r=\sqrt{x_{maj}^2+(y_{min}\times AR{)}^2},$$

$$n(r)=\frac{n_0}{1+(r/R_{\mathrm{flat}}{)}^{\alpha }},$$

Full text

stars: formation — ISM: clouds — techniques: image processing, interferometric — (magnetohydrodynamics:) MHD — methods: numerical

††journal: ApJ††facilities: ALMA††software: ZeusTW; (Krasnopolsky et al., 2010; Zhao et al., 2016), RADMC3D (Dullemond et al., 2012), CASA (McMullin et al., 2007)

1 Introduction

Dense cloud cores represent the initial conditions in the process of individual star formation (Shu et al., 1987). Prior to the formation of young stellar objects in their central regions, these objects are called starless cores. They accrete material from the surrounding less dense cloud and a fraction of them eventually become gravitationally unstable, entering the pre-stellar core phase (André et al., 2014). Pre-stellar cores contract under the pull of gravity, only counteracted by magnetic and thermal pressure, with small contribution from turbulence (Goodman et al., 1998). They can form single (Evans et al., 2015) or multiple (Pineda et al., 2015) stellar systems depending on the interplay between magnetic fields and the partially ionised material present within these clouds (Commerçon et al., 2011; Zhao & Li, 2013). Detailed dynamical models of pre-stellar core evolution have been built, showing that the contraction of rotating and magnetised cold cores is a complex process which does not necessarily lead to the formation of a protoplanetary disk (Li et al., 2014) or multiple stellar systems (Reipurth et al., 2014). To shed light on pre-stellar core evolution and the formation of multiple stellar systems, we ultimately need observational constraints for model predictions. In particular, to study the structure and gas motions within the central 1000 au of a pre-stellar core requires high sensitivity observations of dust continuum emission and molecular lines at size scales similar to those of protoplanetary disks ($\leq$200 au).

Several attempts have been made in the past to detect the central regions and possible substructure of dense cloud cores before star formation (i.e. starless cores, or, if self-gravitating, pre-stellar cores), but no positive results have been collected so far. No millimetre dust continuum emission was detected toward the pre-stellar core L1544 in the Taurus Molecular Cloud Complex, using the IRAM Plateau the Bure Interferometer (Crapsi et al., 2007). The conclusion was that, if an early-stage protostar is present in the center of L1544, its mass cannot exceed 0.2 Jupiter masses (MJ). No dust continuum detection was found within 11 starless cores in the Perseus Molecular Cloud, using CARMA, the Combined Array for Research in Millimeter-wave Astronomy (Schnee et al., 2010). The implication was that Perseus starless cores are characterised by shallow density profiles in the inner few thousand au. Other interferometric observations (with CARMA and SMA, the Submillimeter Array) have also been carried out to investigate the possible presence of substructure toward five bright starless cores in Perseus and Ophiuchus (Schnee et al., 2012). Once again, the non-detections implied that the observed cores have flat density profiles out to at least 1200 au, but no images of these central regions could be recovered.

Bourke et al. (2012) used the SMA to study the dense starless core N6 in the nearest cluster forming region, the $\rho$Oph A cloud in the Ophiuchus Molecular Cloud Complex, and found a poorly resolved compact source with size $\sim$1000 au and mass 0.005-0.01 M*⊙, not seen with single dish observations. A compact structure with mass $\leq$20 MJ was detected toward the starless core SM1N in the same cluster forming region $\rho$Oph A with the Atacama Large Millimetre and sub-millimetre Array (ALMA) in Cycle 1 (Friesen et al., 2014), but the sensitivity was not high enough to study its structure. More recent observations (Friesen et al., 2018) have shown that both SM1N and N6 may be at a later evolutionary stage, with presence of compact dust continuum emission and possible evidence of CO outflow. Another starless core, L1689N in Ophiuchus, has been imaged with the Atacama Compact Array (ACA, also known as Morita Array), by Lis et al. (2016), who found a compact dust continuum source with size similar to N6 but a mass of 0.2-0.4 M⊙*; however, the ACA angular resolution was not high enough to study the inner structure. Moreover, L1689N is located on the East side of the Class 0 binary source IRAS 16293-2422 and it is impacted by the outflow driven by one of the protostars, so its structure and kinematics have been probably affected by the nearby star formation activity. Still focusing on the actively star-forming Ophiuchus complex, Kirk et al. (2017) carried out ALMA Cycle 2 observations of 60 starless and protostellar cores and only detected one starless core, without spatially resolving it.

Cluster forming regions, such as Ophiuchus and/or active sites of star formation, tend to have complex structure, with clear presence of multiple dense cores within a higher density medium (e.g. Friesen et al., 2014; Wang et al., 2014; Henshaw et al., 2016; Ohashi et al., 2018); however, the study of these regions does not allow us to quantify how much of the observed structure is due to the original dynamical evolution of the pristine cloud, i.e. before star formation, and how much of such structure is instead a result of protostellar feedback. More quantitative work on the initial conditions of the process of star formation, including the detailed inner structure of isolated dense cores and their possible fragmentation, is better done in regions not affected by feedback, such as isolated dense cores in molecular clouds (e.g. the pre-stellar core L1544 in Taurus, studied here) or Bok Globules, although this is a difficult task. For example, recent work with (Cycle 1) ALMA toward starless cores in the Chamaeleon Molecular Cloud (Dunham et al., 2016) failed to detect any of the 56 starless cores in their sample. Besides the high sensitivity needed to observe these cold regions, which now is available with ALMA, another problem associated with the detection of nearby starless cores is that only those with high volume densities within compact central regions, or ”kernels” (extending the definition used by Myers, 1998, for cluster-forming molecular cloud cores, to low-mass pre-stellar cores) are expected to be revealed by interferometers; these centrally concentrated cores on the verge of star formation (the so-called pre-stellar cores) are dynamically evolved and have short life times (André et al., 2014). To summarise, pre-stellar cores are rare and difficult to find.

Here we present the ALMA (12 m array and ACA) 1.3 mm dust continuum emission map of the 8 M*⊙* pre-stellar core L1544 in Taurus at a distance of 135 pc (Schlafly et al., 2014, ; consistent with the recent Taurus distance measurements of Galli et al. (2018)), which for the first time resolves its central 1000 au. In past work, L1544 has been modelled as a Bonnor-Ebert (BE) sphere (Bonnor, 1956; Ebert, 1957) in quasi-static contraction (Keto & Caselli, 2010; Keto et al., 2014). The BE sphere has a central region, called the flat zone, with radius $R_{\rm flat}$, where the volume density is constant, and which can be theoretically determined by multiplying the sound speed with the free-fall time at the central density (Keto & Caselli, 2010).

The observations are described in Section 2, observational results are in Section 3, the comparison with smooth cloud models is in Section 4, hydro- and magneto-hydro simulations of contracting pre-stellar cores are presented in Section 5, discussion and conclusions are in Section 6.

2 Observations

We conducted ALMA observations of the L1544 pre-stellar core during Cycle 2 (ESO Project ID 2013.1.01195.S, PI Caselli). Both the 12-m (Main array) and the 7-m (ACA) arrays were employed. The 12-m array observations were carried out on 2014 Dec 27 and 29, while the 7-m array observations were carried out on 2014 June 15, July 20 and 29, and August 6 and 11. The single pointing observations used a correlator configuration where 1 of the 4 correlator spectral windows was used for continuum centred at 228.973 GHz and a 2 GHz bandwidth, while the other 3 spectral windows were used for N2D+, DCO+, and D13CO+. These data will be presented in a follow-up paper dedicated to the kinematics (Paper II).

The 12-m array observations used the quasar J0423-020 and J0510+1800 as bandpass and gain calibrators, respectively, while flux calibration was done using Uranus. The 7-m array observations used the quasar J0510+1800 as bandpass and gain calibrator, while flux calibration was done using J0510+1800, Ganymede, and Uranus observations. The 12-m array data were calibrated using the Common Astronomy Software Applications package (CASA; McMullin et al., 2007) version 4.2.2, while the 7-m array data were calibrated using CASA version 4.2.1. Before combining the data from both arrays we calculated the relative weights using the statwt() command. We perform a joint de-convolution of both datasets (ACA and 12-m array) in CASA (version 4.4.0) using multifrequency synthesis (mfs), a robust weighting of 0.5, and a taper of 1.5*′′, which results in a beam with major and minor axes of 2.05${}^{\prime\prime}\times 1.61^{\prime\prime}$ (277 au $\times$ 217 au at the distance of Taurus) with a position angle of 59.5∘. Since the continuum emission is extended, we use the multi-scale clean technique, with the multi-scale parameter of [0, 2, 6, 18] arcsec. The noise level was estimated using the root mean square over an emission free section of the image and found to be 36 $\mu$Jy beam-1*.

3 Results

Figure 1 shows the ACA-only dust continuum emission map of L1544 (left panel), the 12-m (Main) array map (central panel), and the ACA+Main array map (right panel), all without primary-beam correction to maintain a roughly constant noise across the image. Figure 2 presents a magnification of the combined ACA+Main Array map of the 1.3 mm dust continuum emission.

Figures 1 and 2 show that the high density kernel can now be viewed with unprecedented detail. The kernel radius ($R_{\rm ker}$), derived from the area within the 5 $\sigma$ contour of the combined map is 10.3″ (1446 au). Assuming a dust temperature $T_{d}$ = 6.5 K (Crapsi et al., 2007; Keto & Caselli, 2010) and a dust opacity at 1.3 mm $\kappa_{\rm 1.3mm}$ = 0.009 cm2 g*-1* (corresponding to dust grains with thick icy mantles and 106 yr of coagulation; Ossenkopf & Henning, 1994), consistent with recent findings (Chacón-Tanarro et al., 2019), the kernel mass can be found using:

[TABLE]

where D is the distance of the source (135 pc), $S_{\rm 1.3mm}$ is the primary-beam-corrected flux integrated within the 5 $\sigma$ contour and within the 50% primary beam response in the right panel of Figure 1 (42$\pm$2 mJy, taking into account a 5% calibration error), and $B_{\rm 1.3mm}(T_{d})$ is the Planck function at temperature $T_{d}$. Substituting the numerical values, we find $M_{\rm ker}$ = 0.100$\pm$0.005 M*⊙111The 1 $\sigma$ error of the mass only contains the 1 $\sigma$ error of the observed flux; however, uncertainties on the source distance and on the dust opacity effectively increase this error to about 50 $\%$ of the measured value.. The average H2 number density is then $n({\rm H_{2}})$ = 3 $M_{\rm ker}$/($4\pi\mu_{\rm H_{2}}m_{\rm H}R_{\rm ker}^{3}$) = (1.00$\pm$0.05)$\times$106* cm*-3*, where $\mu_{\rm H_{2}}$ is the molecular weight per hydrogen molecule (=2.8; Kauffmann et al., 2008) and $m_{\rm H}$ is the H-atom mass; this $n({\rm H_{2}})$ value is in agreement with previous work by Crapsi et al. (2007). It is interesting to note that the kernel is elongated close to North-South, unlike the surrounding dense core, whose major axis is tilted by about 45*∘*, as deduced by larger scale observations; we will come back to this misalignment in Section 5.

Fragments are apparent at a 3-$\sigma$ level in the Main array image (central panel in Fig. 1) and at a 2-$\sigma$ level in the combined (ACA+Main Array) map (Fig. 2). The fragments are not spatially resolved, i.e. their size is less than the ALMA synthesized beam of 2.05*′′* $\times$ 1.60*′′* (277 au $\times$ 216 au). Using the above temperature and dust opacity, the mass enclosed within the fragments (the 11 $\sigma$ contour in the central panel of Fig. 1) is 0.003 M*⊙* (3 MJ), i.e. each fragment is about 1 MJ. The fragments appear to be centrally concentrated structures, with peak number densities of around 5$\times$106 cm*-3*, assuming spherical symmetry and a radius of 92 au222The equivalent radius within the 11 $\sigma$ contour of Fig 1, central panel, is 276 au, so, assuming the presence of 3 fragments, we have simply divided this number by 3.. Therefore, the fragments could be local density enhancements, where material has currently assembled, but they cannot be close to gravitational collapse, as their size is significantly smaller that the size scale below which stable oscillations rather than gravitational collapse can occur. In fact, considering an isothermal cloud where gravity is only counteracted by thermal pressure, one can find the size scale below which stable oscillations rather than gravitational collapse will occur. This is the Jeans length $L_{J}$, defined by:

[TABLE]

where $G$ is the gravitational constant, $\rho_{0}$ is the mass density and $c_{s}$ is the sound speed:

[TABLE]

with $k_{\rm B}$ the Boltzmann constant, $T$ the gas temperature, $\mu_{\rm p}$ the mean molecular weight per free particle (2.37; Kauffmann et al., 2008). Assuming a temperature of 6.5 K and the mass density of the kernel of 4.5$\times$10*-18* g cm*-3* (derived from the measured mass, 0.1 M*⊙*, radius, 1446 au, assuming spherical symmetry and uniform density), the Jeans length of the kernel is about 3200 au, about 30 times larger than the size upper limit of the fragments. However, the fragments may also be imaging artefacts due to extended emission viewed with an interferometer and in the next session we will explore this point.

4 Comparison with smooth cloud models

4.1 The centrally concentrated Bonnor-Ebert sphere

We simulated similar ALMA (ACA and Main Array) observations to determine if the current best smooth model of the core could reproduce the observations. The model is generated using the physical parameters found by Keto & Caselli (2010) and refined by Keto et al. (2014), such that they reproduce the previous set of single dish observations of the dust continuum and molecular line emission. The structure deduced by Keto & Caselli (2010) and Keto et al. (2014) is a Bonnor-Ebert sphere (Bonnor, 1956; Ebert, 1957), with a temperature gradient which reproduces the one observed by Crapsi et al. (2007). The high-resolution continuum image is then processed using the task simobserve in CASA 4.6.0, which generates simulated data with a similar integration time and antenna configuration for both ACA and the Main array observations. The two datasets are then imaged using the same parameters as for the real observations, briggs weighting of 0.5, taper of 1.5*′′*, and multi-scale clean. The resulting image (Figure 3) shows a core that is clearly more centrally concentrated than deduced by the observations, with a peak flux about two times larger than that observed. The central zone of a Bonnor-Ebert sphere, where the volume density is constant, is called the flat zone, and its radius is the flat radius, $R_{\rm flat}$; the flat zone is expected to shrink during the dynamical evolution of the core toward the formation of a protostar (see e.g. Keto & Caselli, 2010). L1544 has a larger flat radius than predicted by Keto & Caselli (2010), thus suggesting that it is at an earlier evolutionary phase, or maybe that the L1544 dynamical evolution differs from that of a contracting Bonnor-Ebert sphere (see, e.g., the three-dimensional L1544 structure derived by Doty et al., 2005, using single dish dust continuum emission data), with possible formation of multiple fragments instead of a single peak.

4.2 The smooth elongated structure derived from ALMA observations

The comparison between the observations (Fig. 1, right panel) and the Bonnor-Ebert model (Fig. 3) indicates that L1544 is not as centrally concentrated as predicted. Here we consider a parametric model based on the smooth structure observed by ALMA, taking into account the elongated nature of the core to test further the existence of fragments within the core centre. For this, we fix the peak flux of the model to match the peak flux of the combined (ACA+Main Array) image; this clearly will underestimate the peak flux of the simulated combined image, as the more extended structure already filtered out by ALMA is not taken into account in the model. However, this is a compromise to allow a better fit to the central 3000 au, where most of the flux resides, compared to the outer regions which are better traced by single dish observations (see Section 4.3).

The following function is used to model the structure observed in Figure 1, right panel:

[TABLE]

where $F(r)$ is the flux density at radius $r$, defined as the elliptical distance to the core center:

[TABLE]

where $x_{maj}$ and $y_{min}$ are the coordinates of $r$ along the core major and minor axes, respectively, in a system of reference where the major axis is along the $x$-axis and the minor axis is along the $y$-axis, and $AR$ is the major to minor axis ratio of the whole modelled structure. The best-fit model is calculated by minimising the $\chi^{2}$ of the flux difference between the model (after applying the primary beam correction) and the combined ACA+Main Array data over the 50% primary beam response region. This includes core centre, orientation, axes ratio, and the radial profile parameters. The best-fit model has the following parameters: $F_{0}$=0.772 mJy beam*-1*, $R_{\rm flat}$=10.84*′′* (equivalent to 1463 au), $\alpha$=2.32, $AR$=1.83, position angle PA=166*∘*, measured due East from North. The peak value is fixed and the minimisation is carried out over the region with primary beam response higher than 50%.

The model is then passed through the ALMA simulator observation and processed in the same way as for the Bonnor-Ebert model. Figure 4 shows the results of this analysis. The left panel shows the simulated ACA+12m array map, the central panel shows the simulated 12m-array-only observations and the right panel shows the residuals obtained from the subtraction of the data in the central panel of Fig. 1 and the model in the central panel of Fig. 4. Although the simulated image does not reproduce the observed substructure, the residual image does not show fragments at a significant level, thus implying that the observed substructure is not significant. The next subsection focusses on the longer baselines to highlight the compact emission and further test the smooth-elongated structure model.

4.2.1 The smooth elongated structure observed with baselines longer than 13 k$\lambda$

We now show the comparison between data and simulated observations of the smooth elongated structure model when only baselines longer than 13 k$\lambda$ are used, to give more weight to the substructure. With the aim of confirming or discarding the presence of substructure within the central 1000 au of the pre-stellar core L1544, we generate synthetic ALMA observations using an array configuration similar to the one used during the real observations (Configuration 2 of Cycle 1 included in CASA) and the same on-source time using the simobserve task in CASA. We run 100 models with different starting hour-angles to generate different $uv$-coverage.

We image both the real data and the synthetic ones with the same CLEAN parameters, but since we want to determine how reliable is the identification of small scale structure identified in the image, we only use baselines longer than 13 k$\lambda$ and natural weighting. The resulting set of images have similar beams and a selection of those reproducing substructure similar to that observed is shown in Figure 5. They show that small substructures similar to those identified in the ALMA data of L1544 can be produced when a radial emission profile with a flat central section is used as input image, thus confirming that incomplete cancellation of Fourier components (sidelobes) could produce it.

4.3 Structure derived from single dish observations

As already seen in the previous sub-section, the presence of an extended flat region in the central 1000 au could produce artefacts when viewed with an interferometer. However, in the previous subsection, the peak flux of the simulated combined image was underestimated because the model was constructed from the ALMA data, which did not account for the more extended structure due to filtering. Recently, single-dish dust continuum emission observations of L1544 at 1.1 mm, with the AzTEC camera at the Large Millimeter Telescope (LMT), and at 3.3 mm, with the MUSTANG-2 camera at the Green Bank Telescope (GBT), have been carried out toward L1544 (Chacón-Tanarro et al., 2019). The beam sizes are 12.6*′′* (corresponding to a linear size of 1701 au) at 1.1 mm and 9.7*′′* ($\equiv$1310 au) at 3.3 mm. These observations provided an updated density profile of L1544, which has been found to be similar to the one already proposed by Crapsi et al. (2007), taking into account dust opacity variations, which are found to be consistent with the presence of bare grains in the outer part of the core and the presence of dust grains with thick icy mantles toward the central 2000 au (Chacón-Tanarro et al., 2019):

[TABLE]

where $n(r)$ is the number density of H2 at radius $r$, $n_{0}$ is the central number density. The parameters found by Chacón-Tanarro et al. (2019) are: $n_{0}$ = 1.6$\times$106 cm*-3*, $R_{\rm flat}$ = 17.3*′′* ($\equiv$2336 au), and $\alpha$ = 2.6 (we note that the density profile in Eq. 6 reproduces well the density profile of a Bonnor-Ebert sphere, when $\alpha$ = 2.5, as noticed by Tafalla et al., 2004). Although these single-dish observations cannot resolve the L1544 kernel, we use this structure to test if it can reproduce our ALMA observations. Figure 6 shows the ALMA simulated observations of the pre-stellar core with the spherically-symmetric physical structure derived by Chacón-Tanarro et al. (2019). From Figure 6, it is interesting to see that the core structure closely resembles that observed by ALMA (12m array + ACA), including the flux in the ACA + 12m and 12m-only images, as well as the apparent fragmentation in the central panel, although not at a significant level as in the central panel of Fig. 1. The large black area in the residual image in the right panel of Fig. 6 is due to the presence of extra emission in the model cloud in the East side, as the model is spherical, while the data clearly show elongated structure (with a peak in the West side of the primary beam).

5 Hydro- and Magneto-hydro simulations of contracting pre-stellar cores

Density perturbations within starless cores likely dissipate before becoming gravitationally unstable, as it has been shown by Toci et al. (2018) analytically and also with three-dimensional (3D) hydrodynamic simulations. These authors suggest that solenoidal modes may eventually promote fragmentation via formation of vortical structure, but vortices can be disrupted by magnetic fields due to magnetic stresses. To investigate further this point, we carry out 3D non-ideal magneto-hydrodynamic (MHD) turbulent simulations using ZeusTW code (Krasnopolsky et al., 2010) to model the evolution of the pre-stellar core L1544. We adopt the same chemical network as in Zhao et al. (2016) for obtaining the magnetic diffusivity of ambipolar diffusion (AD).

We first evolve under AD a uniform, spherical, magnetised core with total mass $M_{\rm c}=8.1~{}M_{\sun}$ (similar to the total mass of the L1544 core within the C18O(1–0) contours; Tafalla et al., 1998) and radius $R_{c}=5\times 10^{4}$ au, which corresponds to an initial number density for molecular hydrogen $n({\rm H_{2}})\approx 2.3\times 10^{4}$ cm*-3*. The free-fall time of the core is $t_{\rm ff}\approx 7\times 10^{5}$ yr. The initial core is threaded by a uniform magnetic field along the rotation axis with a constant strength $B_{0}\approx 14.8~{}\mu$G (consistent with Crutcher & Troland, 2000), which gives a dimensionless mass-to-flux ratio $\lambda$ ($\equiv{M_{\rm c}\over{\pi R_{\rm c}^{2}B_{0}}}2\pi\sqrt{G}$) of unity. The initial condition of the first stage is similar to that adopted by Ciolek & Basu (2000), therefore, AD is necessary to allow the magnetised core to collapse, instead of oscillating in the radial direction.

After the core evolves under AD for $\sim$0.7 Myr into a centrally concentrated structure with central number density of $5\times 10^{4}$ cm*-3*, we perturb the core with a $\mathcal{M}\approx 5$ turbulence of Kolmogorov type (with spectral power law exponent -11/3). A small amount of solid body rotation with $\omega=4\times 10^{-14}$ s*-1*, as measured by Ohashi et al. (1999), is also added to the core. Due to the strong magnetic support in the critical core, it is difficult for the turbulence motions to compress the flow compared to the frequent formation of shocks in non-magnetised clouds. As a reference, we also perform a non-magnetised model (HD) for the same core, in which clumps are easily formed by turbulent compression. The core in the HD case also flattens along to the angular momentum direction, due to angular momentum conservation in the rotating cloud; however, the thickness of the flattened region is larger than the MHD model, which results in higher column density than the MHD model for a similar central number density.

In Fig. 7, we show for both the non-ideal MHD and HD models three frames when the peak number density reaches 2$\times$106 cm*-3*, 5$\times$106 cm*-3* and 1$\times 10^{7}$ cm*-3*, respectively. At such stages, the root-mean-square velocities are subsonic for both models, with $\mathcal{M}$ about 0.8 for the HD model and 0.5 for the MHD model, respectively, mainly dominated by the bulk infall motions. The central $\sim$10000 AU of the core has the shape of an oblate ellipsoid. We need to tilt the core by $\theta$=45∘ to match the projected axial ratio (q$\sim$0.45–0.5) in the plane of sky (Crapsi et al., 2007) in the 5000 AU scale. This inclination angle is larger than the 16∘ derived in Ciolek & Basu (2000), because the core is not a perfect slab and the thinner edges of the oblate ellipsoid contribute less to the column density. No obvious substructures are present in the projected column density maps (obtained by direct line-of-sight integration of the number density) in Fig. 7. It is interesting to note that only simulated cores with peak number densities of 1$\times 10^{7}$ cm*-3* reach average number densities as high as those derived by Chacón-Tanarro et al. (2019) within 2336 au (the flat radius found in their work; see Section 4.3), thus indicating that single-dish observations are consistent with the presence of higher density regions, more compact than the telescope beam, within the flat radius.

In the 3D number density map, the central region in the HD model is clearly more clumpy and warped than its MHD counterpart, in agreement with previous work (e.g. Commerçon et al., 2012); the peak of the gravitational potential in the HD model is also shifted away from the cloud centre. These additional asymmetric structures are developed in the HD model because of the flow’s susceptibility to turbulent compression. In the projected maps of Fig. 7, it is difficult to visualise the structure of the core kernel, so we display in Figure 8 the 3D number density maps of the MHD and HD cores which reach peak number density of 1$\times 10^{7}$ cm*-3*. This clearly shows that the HD simulation produces substructure within the kernel, while the MHD simulation produces a smooth structure all the way down to the core centre.

Another noticeable feature of the kernel in our ALMA map is the misalignment between its major axis and that of the NH3(1,1) map observed by Crapsi et al. (2007) at a slightly larger scale (90″; see left panels of Fig.7). Our MHD models are able to reproduce such a misalignment, which is a natural outcome of rotating an elongated kernel (having a non-zero angle $\phi$ with the tilting axis) around the tilting axis of the oblate ellipsoid. The original orientation angle $\phi$ of the kernel can be obtained via the Rodrigues rotation formula, explained in the following text and with the help of Figure 9. Assume that vector ${\bf v}$ = (cos$\phi$, sin$\phi$, 0) represents the original orientation of the kernel in the plane of the oblate ellipsoid, and vector ${\bf k}$ = (1, 0, 0) represents the tilting axis around which the ellipsoid rotates. After inclining the ellipsoid by an angle $\theta$, ${\bf v}$ also rotates to a new vector ${\bf v}_{\rm rot}$ = (cos$\phi$, sin$\phi$ cos$\theta$, sin$\phi$ sin$\theta$), whose projection on the plane of sky ${\bf v}_{\rm proj}$ = (cos$\phi$, sin$\phi$ sin$\theta$) can be related to the misalignment angle $\delta$ between the major axis of the kernel and that of the nucleus. Therefore, tan$\phi$ = tan$\delta$/sin$\theta$, which yields $\phi\approx 45^{\circ}$ for the kernel of L1544. This implies that the kernel is by itself an elongated substructure that is well-distinguished from the larger oblate core.

5.1 Simulated ALMA observations

We first apply radiative transfer to the simulation models using RADMC-3D (Dullemond et al., 2012). The temperature distribution of the core, ranging from $\simeq$5 K to $\simeq$15 K, is computed taking into account the external illumination from the interstellar radiation field through RADMC-3D; the temperature in the central 3000 au is $\simeq$6 K, consistent with the measured values (Crapsi et al., 2007). Simulated ALMA Main Array observations of the HD and MHD cores in Fig. 7, following the same procedure as for previous simulated observations presented in Sections 4, can be found in Fig. 10, where the panels have been centred at the peak flux positions of the cores in Fig. 7. The order of the panels in Fig. 10 is the same as in Fig. 7, with the MHD cores in the left panels and HD cores in the right panels; the cores have central number densities increasing from top to bottom: 2$\times 10^{6}$ cm*-3*, 5$\times 10^{6}$ cm*-3* and 1$\times 10^{7}$ cm*-3*. It is interesting to note that the less centrally concentrated cores (those with peak number densities 2$\times 10^{6}$ cm*-3* and 5$\times 10^{6}$ cm*-3*) are completely filtered out by the Main Array because of their smooth structure, while the cores with peak number densities of 1$\times 10^{7}$ cm*-3* can be well detected, both in the MHD and HD cases. This could explain why detecting starless dense cores with interferometers has so far been unsuccessful (as reviewed in Sect. 1): only (short-lived) dynamically evolved cores on the verge of star formation (in our case, those with peak number densities larger than 5$\times$106 cm*-3*), can be detected, as also discussed by Dunham et al. (2016). In particular, we note that the 107 cm*-3* MHD core has a peak flux about 1.5 times lower than the observed core (see central panel of Fig. 1), while the 107 cm*-3* HD core has a peak flux about 1.5 higher than observed. This implies that with the present dust continuum emission observations it is not possible to distinguish between MHD and HD models and that high sensitivity and high spectral resolution molecular line observations are needed to quantify the role of magnetic fields in the dynamical evolution of pre-stellar cores. The fact that the observed peak flux is close to that of the simulated cores with central peak density 107 cm*-3* implies structural similarities with L1544, with the MHD simulation also reproducing the misalignment between kernel and dense core major axes. It is finally interesting to note that, on the one hand, the simulated observations of the MHD core with peak density $10^{7}$ cm*-3* presents fragmentation at the three-sigma level within the central 1000 au, despite a smooth model core (see Fig. 8), showing that the detection of fragments can be due to interferometric artefacts; on the other hand, the simulated observations of the HD core with peak density $10^{7}$ cm*-3* shows a smooth structure (see Fig. 10), despite a fragmented model core (see Fig.8), due to the small size of the fragments, well embedded within the kernel.

6 Discussion and conclusion

ALMA (Main Array + ACA) observations of the dust continuum emission at 1.3 mm of the 8 M*⊙* pre-stellar core L1544 have revealed a kernel (defined by the 5$\sigma$ contour) with radius of 1446 au, mass of 0.1 M*⊙* and average H2 number density of 1$\times 10^{6}$ cm*-3*. Within the kernel, two fragments at the 3-$\sigma$ level in the Main Array image, and three fragments at the 2-$\sigma$ level in the Main Array+ACA image, are found. The fragments, if real, are unresolved structures, with radii $\leq$92 au, mass $\leq$1 MJ, and peak volume densities $\simeq$5$\times 10^{6}$ cm*-3*; they cannot be close to gravitational collapse.

The comparison of our ALMA data with simulated observations of various smooth cloud models shows that: (1) the Keto et al. (2014) Bonnor-Ebert sphere is more centrally concentrated, with central fluxes about two times larger, than observations; (2) in the case of the smooth cloud model with structure deduced from our ALMA observations, substructure similar to the observed fragments appears in the simulated observations, confirming that incomplete cancellation of Fourier components can produce the observed fragments; (3) when the latest (spherically symmetric) physical structure of the L1544, deduced from single dish observations, is used as input in the simulated observations, a good match with the peak fluxes in the ACA+Main Array and Main Array-only images is obtained, although substructure at the observed $\sigma$ levels is not reproduced.

Non-ideal MHD and HD simulations of a dense core similar to L1544 have been run until the central density of the core reaches peak number densities 2$\times 10^{6}$ cm*-3*, 5$\times 10^{6}$ cm*-3* and 1$\times 10^{7}$ cm*-3*. The MHD simulated cores have a smooth structure and present a misalignment between the major axis of the dense core and that of the kernel, similar to that found with observations. The HD simulated cores present substructure and warped kernels. When ALMA Main Array simulated observations are carried out, only model cores with peak number densities of 1$\times 10^{7}$ cm*-3* can be detected. These cores have average number densities, within a flattening radius of $\simeq$2300 AU, close to those deduced by recent single dish observations ($\simeq$1.6$\times 10^{6}$ cm*-3*; Chacón-Tanarro et al., 2019). Thus, only short-lived dynamically evolved and centrally concentrated pre-stellar cores can be studied at high angular resolution; this explains why it has been so far difficult to detect starless cores with interferometers, including ALMA. The simulated observation of the high density MHD core shows substructure within the kernel at levels similar to the observed ones, despite the smooth model core; this again shows that interferometric artefacts can mimic fragmentation. The simulated observations of the high density HD core shows a smooth structure despite the model core being highly fragmented toward the central region; this is due to the small size of the fragments compared to the synthesized beam, suggesting that even higher angular resolution observations should be carried out to test MHD and HD predictions. In the future, we will compare the simulated and measured velocity fields to investigate more quantitatively the non-ideal MHD and HD predictions (Pineda et al., in prep.).

With the current data, we cannot claim evidence of fragmentation. Nevertheless, fragmentation can take place. The Taurus Molecular Cloud is known to host a large number of binary pairs: between 2/3 and 3/4 of the whole population of young stellar objects are multiple systems of two or more stars, while the rest appears to have formed as single stars (Kraus et al., 2011). The separation of the observed stellar multiples shows values around a few hundred au, thus indicating that fragmentation may happen within the flat central region of L1544. Basically, star formation does happen and it tends to produce multiple systems. Core fragmentation is one of the main expected modes of binary/multiple system formation (Reipurth et al., 2014), but the present data, toward one of the most massive and dynamically evolved starless dense cores, do not support this scenario, suggesting that either L1544 will form a single star or that fragmentation may happen at a later stage of evolution, maybe soon after the formation of the first protostar (see e.g. Pineda et al., 2015). Alternatively, fragmentation could indeed happen at pre-stellar stages, but higher sensitivity observations are needed to test this scenario. Our data are also pointing out that caution should be taken when claiming fragmentation within pre-stellar cores with ALMA.

We acknowledge the very helpful comments of our referee, Neal Evans, and Anika Schmiedeke for her help on the radiative transfer. P.C., J.E.P., B.Z. and A.C.-T. acknowledge the financial support of the European Research Council (ERC; Advanced Grant PALs 320620). MP acknowledges funding from the European Unions Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 664931.