© The Authors 2023

1 Introduction

Class 0 protostars are the youngest objects in the protostellar phase that lead to the formation of solar-type stars. They are embedded protostars, in which most of the mass is still in the form of a dense cold envelope surrounding the central developing protostellar embryo (Andre et al. 1993; Andre & Montmerle 1994). Because of their embedded nature, their spectral energy distribution peaks at submillimeter wavelengths, as most of the accretion luminosity produced in the center gets reprocessed by the dense envelope. The accretion of the envelope material is also accompanied by the ejection of material via a (bi)polar jet and/or outflow system. Both the radiative energy liberated at the accretion shock on the protostellar embryo and the mechanical energy transported by the jet and outflow contribute to opening and shaping outflow cavities.

The luminosity of a protostar is the sum of the luminosity of the central protostellar photosphere and the accretion luminosity, which corresponds to the fraction of energy carried by the infalling matter that is radiated away at the accretion shock. Generally, constraints on protostellar accretion come from the total luminosity of protostars, given that the dusty envelope absorbs the radiative energy coming from the center and re-emits it at (far-)infrared and submillimeter wavelengths. Several surveys (Kenyon et al. 1990, 1994; Evans et al. 2009; Dunham et al. 2010, 2013) supported the existence of a luminosity problem because the measured total luminosities are on average too low compared to what can be predicted by protostellar accretion models (Myers et al. 1998; Young & Evans 2005) that consider a typical protostellar lifetime (Dunham et al. 2015; Kristensen & Dunham 2018). While Offner & McKee (2011) and Myers (2011, 2014) constructed models that reasonably predict the observed luminosities without including outburst activity, there is evidence suggesting that protostellar luminosities can be variable on timescales of months or years (Fischer et al. 2019, 2022; Lee et al. 2021b; Park et al. 2021; Zakri et al. 2022), which is also supported by chemical studies of protostellar envelopes (Visser et al. 2012, 2015; Jørgensen et al. 2013; Anderl et al. 2016; Frimann et al. 2016) and the variability of emission line outflow tracers as a function of the distance from protostars (Arce & Goodman 2001; Plunkett et al. 2015). Such a variability could point to episodic accretion (Baraffe et al. 2009; Vorobyov & Basu 2010), episodic ejection (see, e.g., recent models from Commerçon et al. 2022), viscous and magnetic instabilities at the boundary between the stellar magnetosphere and the accretion disk (Kulkarni & Romanova 2008; D’Angelo & Spruit 2012; Takasao et al. 2019), binary interactions (Bonnell & Bastien 1992), or instabilities of the structure of the protostellar embryo itself (Vaytet et al. 2018).

The accretion shock converts the kinetic energy of the infalling circumstellar gas into a UV radiation field that is rapidly reprocessed by the high density gas and grains surrounding the protostellar embryo. Hence, the protostellar accretion rate and conditions responsible for reprocessing the accretion radiation in the first inner au are key to setting the temperature and alignment efficiency of dust grains, the emission of which is in turn widely used to probe both the mass of material in protostellar environments and the magnetic field topology through their polarized emission. Indeed, aspherical grains can align their minor axis parallel to the ambient magnetic field lines via the phenomenon of radiative torque alignment (RAT; Draine & Weingartner 1996, 1997; Lazarian & Hoang 2007; Andersson et al. 2015). Param-agnetic dust grains get internally aligned via, for example, the nuclear relaxation, the inelastic relaxation, and the Barnett internal relaxation processes (Purcell 1979; Lazarian & Draine 1999; Lazarian & Efroimsky 1999; Lazarian & Hoang 2007; Hoang & Lazarian 2009; Hoang et al. 2022), and interact with the external magnetic field such that the grain angular momentum performs a Larmor precession around the ambient magnetic field (Dolginov & Mitrofanov 1976). Then, the radiative torques, which result from the different scattering efficiencies of an impinging photon beam decomposed into right- and left-hand circular polarized light, drive the alignment between the grain angular momentum and the magnetic field. However, when embedded in a strong radiation field, grain precession can preferably happen around the anisotropic component of the radiation field rather than the magnetic field orientation (k-RAT; Lazarian & Hoang 2007; Tazaki et al. 2017).

The efficiency of the RAT grain alignment mechanism is a balance between the radiative torques’ efficiency (which is spinning dust grains up) and the collisional de-alignment of dust grains due to gas pressure (the rotational damping is also affected by the infrared emission of dust grains, but it is negligible in protostellar envelopes; Draine & Lazarian 1998; Hoang et al. 2021). In practice, this balance determines for which grains the local conditions enable their alignment with magnetic field lines via RAT. A dust grain is considered aligned if its alignment is sufficiently stable, that is, if the grain’s RAT-induced angular momentum, J_max(ψ) (where ψ is the angle between the radiation field vector and the magnetic field), is sufficiently large compared to the grain thermal angular velocity, J_th (i.e., for J_max(ψ)/J_th ≈ 3; Hoang & Lazarian 2008). This ratio of angular momentum defines the condition for grains to be aligned, which translates into a critical grain size, a_align. Dust grains larger than a_align can be considered aligned, and this condition is thus set by the environment parameters (i.e., radiation field, gas density, and temperature). Given the large mean wavelength| of the photons propagating in the envelope of protostars from the central object (i.e., the mean wavelength of the radiation field that has been reddened by the inner envelope material only, on the order of ~ 1–10 µm; see Hoang et al. 2021), it was found that dust grains, in order to be aligned, must be larger than the typical size of interstellar medium (ISM) dust grains (i.e., larger than ~10 µm; Le Gouellec et al. 2019; Valdivia et al. 2019; Hull et al. 2020). A homogeneous and high level of grain alignment efficiency was found in the diffuse ISM (Planck Collaboration XII 2020; Reissl et al. 2020). Looking at the inner regions of protostellar cores, Le Gouellec et al. (2020) find that the grain alignment efficiency in the envelope of Class 0 protostellar cores is on average relatively high, on the order of perfectly aligned grains (see also Kuffmeier et al. 2020). To reproduce the polarization fractions observed in the ISM and in dense environments, it was proposed that the efficiency of RAT can be greatly increased for grains containing super-paramagnetic inclusions, such as iron inclusions (Yang 2021; Chau Giang et al. 2023). Indeed, while the magnetic relaxation of ordinary paramagnetic grains has been found to inefficiently align thermally rotating grains with the magnetic field (mechanism introduced by Davis & Greenstein 1951, consisting of a dissipation of the grain rotational energy due to the rotating magnetization with respect to the grain main axis), super-paramagnetic relaxation can enhance the alignment degree induced by RAT (Lazarian & Hoang 2008; Hoang & Lazarian 2016b,a). Namely, the super-paramagneticity of dust grains increases the fraction of grains subject to RAT-induced supra-thermal rotation, up to which dust grains can be considered perfectly aligned.

Atacama Large Millimeter/submillimeter Array (ALMA) observations revealed polarized dust emission in several Class 0 circumstellar envelopes at ~50−1000 au scales (Hull et al. 2017, Hull et al. 2020; Cox et al. 2018; Maury et al. 2018; Sadavoy et al. 2018a,b, 2019; Kwon et al. 2019; Le Gouellec et al. 2019; Takahashi et al. 2019; Ko et al. 2020), and the prominence of dust polarization features following outflow cavity walls was noted. Le Gouellec et al. (2020) present regions of low polarization with organized magnetic fields, noting that the grain alignment efficiency does not seem to be totally homogeneous throughout the inner envelope and appears to depend on the local physical conditions, such as the irradiation field (see Pillai et al. 2020). In this work we focus on the specific role of the radiation field in the grain alignment in protostellar envelopes. This allows us to discuss the properties of the aligned dust grains. In Le Gouellec et al. (2023, Paper I hereafter), we compared ALMA dust polarization and radiation-sensitive molecular line observations for a handful of Class 0 protostellar objects, which had a variety of radiation field and polarized intensity enhancement morphologies. For Class 0 protostars, where the spatial distribution of the circumstellar material is isotropic and only one clear outflow emanates from the central object, we found in Paper I that observations suggest the local irradiation is potentially correlated to the regions of enhanced polarized dust emission.

The current paper is structured as follows. In Sect. 2 we present the magnetohydrodynamic (MHD) model of a proto-stellar core that we synthetically observe with the polarization radiative transfer code POLARIS (Reissl et al. 2016), which incorporates the RAT grain alignment mechanism. We analyze the synthetic dust polarization maps and compare the grain alignment efficiency resulting from different levels of irradiation with those found in Le Gouellec et al. (2020) in Sect. 3. Finally, in Sect. 4 we discuss the results of the radiative transfer calculations in light of the radiation field characteristics, dust evolution, and accretion variability in Class 0 protostars. We draw our conclusions in Sect. 5.

2 Model of a collapsing protostellar core

In this section we present the radiative transfer calculations we have performed onto a radiation nonideal MHD simulation. This simulation, representing a prototypical star-forming object, is then used a series of radiative transfer calculations, where a range of several parameters are explored, namely the protostar luminosity (L_★), the maximum dust grain size (a_max), and the fraction of grains aligned in suprathermal rotation (f_high−J).

2.1 RAMSES MHD simulation

We used a detailed model of the evolution of a protostellar object that includes many of the relevant physics and allows us to accurately test how the central luminosity gets reprocessed in realistic envelope conditions. Further details and consideration are presented in Appendix A. We used the RAMSES code (Teyssier 2002; Fromang et al. 2006; Commerçon et al. 2011b) with the implementation of sink particles (Krumholz et al. 2004; Bleuler & Teyssier 2014), ambipolar diffusion (Masson et al. 2012), and adaptive mesh refinement (AMR) to simulate the gravitational collapse of a protostellar core, with radiation nonideal MHD calculations. The reference scenario we adopted here follows the collapse of a magnetized (the mass-to-flux ratio µ is 5), intermediate-mass starless core of 30 M_⊙ initial mass, without initial turbulence, with an initial density profile of ρ ∝ 1/(1+r²). We used the hybrid method developed by Mignon-Risse et al. (2020), which consists of a hybrid radiative transfer method using the gray M1 closure relation (which implements the zeroth and first moments of the equation of radiative transfer; Levermore 1984; Rosdahl et al. 2013; Rosdahl & Teyssier 2015) for the radiation emanating from the protostellar embryo, and the gray flux-limited diffusion approach (Levermore & Pomraning 1981; Commerçon et al. 2011a, 2014) for photons emitted elsewhere in the simulation. After it has formed, the central protostellar embryo is represented by a sink particle. A jet is implemented by hand at the creation of the sink particle (the method used is the one of Verliat et al. 2022). The resolution of the highest level of the AMR grid is ~5 au.

We present in Fig. 1 the snapshot of this simulation that we use from now on. The simulation is 38.52 kyr old, and the central sink, which mimics the central protostellar embryo, is 14.8 kyr old, and has a mass of 1.2 M_⊙. As no initial turbulence is implemented, the collapse proceeds isotropically and the core exhibits a symmetric structure, which favors the accumulation of material toward the equatorial mid-plane. The jet rapidly clears out outflow cavities, which also causes the accumulation of material toward the outflow cavity walls, visible in the column density map and slice of mass density in Fig. 1, top and bottom panel, respectively. The line of sight we use from now on is an edge-on projection of the core.

Fig. 1 Snapshot of our MHD simulation. Top panel: system seen edge-on. The color scale and white contours represent the gas column density. The yellow streamlines trace the orientation of the density-weighted average magnetic field along the line of sight. The white arrows represent the velocity field in a slice of the simulation centered on the sink particle. Bottom panel: the contours and white arrows are the same as in the top panel. The color scale represents a slice of gas mass density, centered on the sink particle. In both panels, the little white circle is the full width half maximum of a circular Gaussian kernel we use to slightly smooth the highest resolution of the AMR grid, for visibility.

2.2 Radiative transfer of the dust emission with POLARIS

2.2.1 Radiative transfer parameters

We performed radiative transfer calculations on this simulation using the POLARIS code (Reissl et al. 2016), which calculates the local dust temperature and dust grain alignment efficiency of oblong dust grains with respect to the magnetic field orientation following the RAT theory developed in Lazarian & Hoang (2007) and Hoang & Lazarian (2014). The POLARIS code computes first the propagation of photons alongside the dust temperature via a Monte Carlo analysis. Then, using the density, radiation field, temperature, and dust grain properties, POLARIS computes the RAT-induced grain angular momentum and compares it to the thermal angular velocity to derive the a_align parameter in each cell of the grid. Finally, POLARIS solves the radiative transfer for the four Stokes parameters, at a given wavelength and line of sight. The only source of heating we include is the one emanating from the central source. We assumed that the impact of the external heating from the interstellar radiation field (ISRF) or the nearby young stars can be neglected in those embedded objects because of the small spatial scales targeted by these synthetic observations. The luminosity from the protostar photosphere and the accretion shock is modeled, for simplicity reasons, as a blackbody source located at the sink particle, whose temperature we vary in our different sets of calculations, and is derived as follows.

The parameters of the equivalent protostellar embryo were derived from the models developed by Kuiper & Yorke (2013; see also Hosokawa & Omukai 2009), which provide the radius and luminosity of the central protostellar object (R_★ = 1.23 R_⊙ and L_★ = 0.58 L_⊙ for our MHD model snapshot). Then the effective accretion luminosity needs to be derived. This is, however, a very degenerate quantity for several reasons: a fraction of the accretion energy can be radiated away through mechanical processes such as jets, the accretion is known to be episodic and can undergo bursting phases of accretion, and the accretion rate onto the sink may be different from the actual accretion rate of material onto the central protostellar embryo (because the accretion can be regulated by a disk, or happen directly form the inner envelope; see Lee et al. 2021a). Therefore, the determination of the exact value of the luminosity escaping from the sink particle is degenerate. In order to account for the variable accretion activity, we explore a range of luminosities, from 1 to 100 L_⊙, by varying the temperature of the central blackbody T_★, and choosing for R_★, the equivalent radius derived from the models cited above (where T_★ and R_★ are the parameters implemented in POLARIS that set the characteristics of the central blackbody, i.e., the source of heating). In addition, we also note that in our simulation the accretion rates onto the sink are large, with average values of Ṁ_acc ~ 10⁻⁵−10⁻⁴ M_⊙ yr⁻¹ (the accretion rate onto the sink in our simulation can vary by ~1−2 orders of magnitude), compared with the common bolometric luminosities of low- and intermediate-mass Class 0 star-forming objects (i.e., Ṁ_acc ~ 10⁻⁷−10⁻⁵ M_⊙ yr⁻¹; Evans et al. 2009). This may be due to the simulation being fairly massive and/or the spatial resolution, which do not allow us to resolve the disk below 5 au, where the accretion mechanisms are occurring. In order to reproduce different states of accretion (i.e., from steady or low accretion activity) to episodes of high accretion activity, we thus adopt different central luminosities (values of L_★ from 1 to 100 L_⊙) for the same MHD model snapshot that represent the fraction of accretion kinetic energy that is radiated away in the accretion luminosity.

We used a population of dust grains with typical ISM composition, with 62.5% astronomical silicates and 37.5% graphite grains (Mathis Rumpl Nordsieck, MRN dust; Mathis et al. 1977; optical properties from Weingartner & Draine 2001). This composition governs the ultimate maximum number of aligned grains, as silicates can be aligned with the magnetic field much more easily than graphite or carbonaceous grains (Andersson et al. 2015 and references therein; see, however, the recent work by Lazarian 2020). The dust grains are assumed oblate with an aspect ratio of 0.5 (Hildebrand & Dragovan 1995) and they follow a standard MRN-like distribution (Mathis et al. 1977) with cutoff sizes of a_min = 2 nm and various values of a_max. As the maximum grain size seems to be a key parameter that controls the resulting amount of polarized flux obtained in synthetic observations (Valdivia et al. 2019), we vary this parameter a_max from 0.2 to 50 µm (see Table 1). We chose these values in light of recent works that have hinted at the presence of grains larger than the typical ~ 0.5 µm ISM dust grain maximum size in Class 0/I envelopes (e.g., Miotello et al. 2014; Valdivia et al. 2019; Le Gouellec et al. 2019; Galametz et al. 2019; Agurto-Gangas et al. 2019; Hull et al. 2020; Nakatani et al. 2020; Ohashi et al. 2021). Indeed, the reprocessed radiation field impinging on the dust grains in the protostellar envelope needs to comprise photons whose wavelength are comparable to the size of dust grains in order to efficiently align the grains via RAT. Given the density of the inner envelope structures, only low-energy submillimeter photons can propagate, which in turn corresponds to grains with sizes >1 µm. We assume a gas-to-dust ratio of 100.

The radiative transfer computations carried out by the POLARIS code include several approximations to describe the effects of grains’ paramagneticity on the grain dynamics. Among the grains aligned via RAT (i.e., grains with a ≥ a_align), a given fraction of grains aligned at supra-thermal rotation is specified, that is to say, the fraction of aligned grains at the high-J attractor point, /high-J (Lazarian & Hoang 2007). Grains are considered imperfectly internally aligned at low-J, and perfectly internally aligned at high-J. For the 1 − f_high−J fraction of low-J aligned grains, the internal alignment is inefficient because of internal thermal fluctuations in the grains (Lazarian & Roberge 1997). To model the imperfect internal alignment of these grains, a Boltzmann distribution is used to describe the precession of grain’s major axis of inertia with its angular momentum. This version of POLARIS does not include the effect of wrong internal alignment of low-J aligned grains, which means we might overestimate the alignment efficiency of low-J aligned grains in our study. In addition, we consider that the high-J aligned grains are perfectly internally aligned. This is true for certain conditions or gas density, irradiation, and level of iron inclusions locked in the dust. Indeed, the amount of iron inclusions in dust increases significantly the grain magnetic susceptibility, which increases the rate of Barnett relaxation. Following Sect. 4 of Hoang et al. (2022), we see that for the typical conditions of our model − n_H ~ 10⁶ − 10⁸ cm⁻³ and u_rad/u_ISRF ~ 10² −10⁵ for the envelope ≲2000 au (see Sect. 3) − grains whose size a ≲ 10–50 µm should be efficiently internally aligned, depending on the level of iron inclusions locked in the dust. We do not consider maximum grain sizes a_max ≥ 50 µm. Grains are considered perfectly externally aligned in both the low- J and high- J cases in our radiative transfer models, which also represents an approximation, especially for the low- J aligned grains whose magnetic relaxation is less efficient. In our use of POLARIS, this makes the f_high−J parameter the parameter that describes the magnetic response of dust grains, because super-paramagnetic grains (i.e., grains with iron inclusions; Jones & Spitzer 1967) are mainly driven to high-J state due to the effect of radiative torques and increased efficiency of magnetic relaxation. In a given model, f_high−J is fixed for all grains (however, see Chau Giang et al. 2023, who varied f_high−J using a step function based on the local magnetic relaxation conditions). Table 1 presents all the sets of radiative transfer runs that we have performed, varying a_max, f_high−J, and the luminosity, L_★. Because the focus of the paper is the role of the radiation field in the dust polarization in protostars, we present the impact of a_max, and f_high−J, on the resulting dust polarization in Appendices B and C, respectively.

Fig. 2 Dust polarization radiative transfer results of the fiducial case (amax = 10 µm, fhigh−j = 1, and L★ = 20 L⊙). The plots are from one run of POLARIS, where our simulation has been synthetically observed in dust polarization at 0.87 mm. Top left panel: the color scale and black contours represent the dust continuum total intensity (Stokes I). The white arrows represent the velocity field in a slice of the simulation centered on the sink particle. Top right panel: the color scale represents the polarized intensity (). The yellow streamlines trace the orientation of the density-weighted average magnetic field along the line of sight. The line segments represent the polarization position angle orientations. Bottom left panel: the color scale represents the fractional polarization, Ƥfrac. Bottom right panel: the color scale is the optical depth computed during the ray-tracing of the radiative transfer. In the bottom panels, the white contours trace the total intensity.

2.2.2 Synthetic observations

We chose a line of sight where the model is seen edge-on to produce with POLARIS the output dust emission maps from the model, at 0.87, 1.3, and 3 mm, at a distance of 400 pc, in maps 8000 au in size with pixel sizes of 8 au. We present our fiducial case of the dust polarization maps calculated by one POLARIS radiative transfer run in Fig. 2. This fiducial case corresponds to the radiative transfer with a_max = 10 µm, f_high−J = 1, and L_★ = 20 L_⊙, performed at 0.87 mm. We note that the cavity walls exhibit enhanced emission in the dust continuum emission map, which is an important feature generally observed with ALMA toward star-forming objects (Maury et al. 2018; Kwon et al. 2019; Le Gouellec et al. 2019; Hull et al. 2020). This is likely due to the implementation of the jet, which has pushed aside material that accumulated in a shell-like structure, at the junction between the collapsing envelope and the outflowing gas. This enhancement in the dust continuum map is less visible when analyzing the column density map (see Fig. 1). The radiation field in the direction of the outflow heats the dust in the cavity walls, enhancing the total (Stokes I) thermal dust emission. Our radiative transfer results also reproduced the enhanced polarized dust emission toward the outflow cavity walls, another common signature observed in ALMA observations (Hull et al. 2017, 2020; Maury et al. 2018; Kwon et al. 2019; Le Gouellec et al. 2019; Ko et al. 2020). The lines of sight toward the base of the cavities (up to ~700 au from the center), and the inner mid-planes (up to ~400 au), exhibit low fractional polarization (~1%), which is explained by the highly disorganized magnetic field and high gas volume density in these two regions, respectively. In contrast, the outflow cavity walls, and outer mid-planes, exhibit enhanced signal in polarized dust emission, up to ~15%. Analyzing the optical depth map, we can assume the emission is optically thin in the envelope, except for the inner ~200 au, which exhibit an optical depth larger than unity.

The enhancement that we obtain toward the cavity walls in total intensity and polarized dust emission is explained by the propagation of the photons computed by POLARIS. Figure 3, which presents the 2D slices of temperature before and after the radiative transfer calculation for the fiducial case (left and right, respectively), shows how the material in outflow cavities is preferentially heated, compared with the equatorial mid-planes. In the radiative transfer results, the temperature reaches ~100K in the central 100, and at 1000 au from the protostar, reached ~40 K in the direction of the cavities but only ~20 K in the direction of equatorial mid-planes. The differences in the morphology of the temperature maps between before and after the radiative transfer are large toward the outflow cavities. During the nonideal MHD simulation, the central source of heating only corresponds to the luminosity of the central protostellar object (L_★ = 0.58 L_⊙ in the simulation; as the point is to discuss only the different heating mechanisms, the fiducial model with L_★ = 20 L_⊙ is used for the radiative transfer results shown in Fig. 3). However, RAMSES takes into account the dynamics of the gas to compute the local temperature. Given the velocity in the corresponding cells, the values of the ion-neutral gas friction parameter gets significantly high, and causes the encountered high gas temperature in a thin layer along the cavity walls. This physics is not accounted for in POLARIS, which only derive the local heating given the propagation of photons, and opacity of cells encountered in the photons propagation pathways. In other words, the dynamics of the outflow can represent a source of heating, and in turn a modification in the local physical conditions for grain alignment, that we do not have access to in our radiative transfer calculations (see Sect. 4.1.2).

We used the Common Astronomy Software Applications (CASA, version 5.8) simulator to interferometrically filter the synthetically observed maps, mimicking ALMA observations. For each simulation, we combine synthetic observations from ALMA configurations C-2, C-4, and C-6, with an exposure time of three hours per antenna configuration. After an appropriate slight smoothing, the resulting synthesized beams (resolution elements) of these filtered maps have an effective size of 100 au, when observed at 0.87 mm (and 140 au, when observed at 1.3 mm). We present in Fig. 4 an example of the resulting polarized dust emission maps after the ALMA synthetic observations. The large-scale emission (i.e., the emission corresponding to low spatial frequency signal) is removed and a standard sky noise model is applied to the data. We recover a significant part of the total intensity signal inside the central ~2000 au, while the detection of the polarized intensity suffers more from the resulting sensitivity (as noticed in Reissl et al. 2017; Kuffmeier et al. 2020; Le Gouellec et al. 2020, the intrinsic grain alignment efficiency is generally lower in radiative transfer models compared to observations; see Sect. 3.2). The main features in the polarized intensity map are the polarized outflow cavity walls and equatorial mid-plane. We note also that ~ 100–200 au above and below from the equatorial plane, we recover a horizontal apparent magnetic field, suggesting a significant toroidal magnetic field component (also seen in Figs. 1 and 2). As shown in (Le Gouellec et al. 2020), interferometric filtering has a strong impact on the polarization, as it artificially increases the polarization fraction values. Indeed, Stokes I and Stokes Q and U have different power spectrum (when analyzed as a function of spatial scales), due to the fact that Stokes Q and U are sensitive to the total intensity as well as the grain alignment conditions and magnetic field morphology. Therefore, to compare these synthetic observations with observations and discuss the grain alignment efficiency, this instrumental filtering is necessary (see Sect. 3.2).

Fig. 3 Temperature maps from RAMSES and POLARIS. The left panel presents the temperature map of the RAMSES simulation cube before the POLARIS radiative transfer, from the slice centered on the sink particle (in the simulation the only central source of heating is the luminosity of the central protostellar object, with L★ = 0.58 L⊙). The right panel shows the temperature resulting from the radiative transfer performed by POLARIS for the reference case (amax = 10 µm, fhigh−J = 1, and L★ = 20 L⊙). The point is not to compare the temperature values across the two maps, but to highlight that the heating mechanisms are different - in the RAMSES simulation the dynamics of the gas induce a significant heating of a thin layer surrounding the high velocity jet.

Fig. 4 Dust polarization radiative transfer results performed at λ = 0.87 µm with fhigh−J = 1.0, L★=100 L⊙, and amax = 10 µm, after spatial filtering with the CASA simulator, mimicking ALMA interferometric observations. Top panel: synthetically observed dust continuum total intensity (Stokes I, color scale) plotted with the velocity field of the central slice in the MHD simulation (white arrows) and the density-weighted magnetic field lines (yellow lines). Stokes I is shown when I >3 σI, where σI = 0.10 mJy beam−1. Bottom panel: the color scale is the total linearly polarized intensity, which is shown where Ρ>3 σP, where σΡ = 10 µJy beam−1. The line segments represent the magnetic field orientations inferred from the dust polarization map. They are plotted where Ρ > 5 σΡ. The white circle at the bottom-left corner is the resolution of the ALMA synthetic observations and is 100 au in size.

3 Analysis of the synthetic polarized dust emission maps

3.1 Dust polarization maps

We present in Fig. 5 set III of radiative transfer runs performed at 0.87 mm (while we vary L_★, f_high−J = 1 and a_max = 10 µm; see Table 1). For each value of L_★ (i.e., 1.0, 5.0, 20, 50, and 100 L_⊙), we show the total intensity, polarized intensity, polarization fraction, dust temperature, radiation field, and minimum size of aligned dust grains (a_align). For the radiation field, we plot u_rad/u_ISRF, where with u_λ is the spectral energy density at the wavelength λ, and u_ISRF = 8.64 × 10⁻¹³ erg cm⁻³ is the interstellar radiation field from Mathis et al. (1983).

Regarding the temperature and radiation field maps, Fig. 5 shows the progressive heating of the inner envelope as L_★ increases. In the most luminous case we consider (L_★ = 100 L_⊙), the temperature reaches ~ 70 K and the radiation field reaches ~ 10⁴ u_ISRF, at 1000 au in the direction of the outflow cavities. In the lowest luminous case (L_★ = 1 L_⊙), the polarized intensity does not reach observable values, while in the higher luminous cases (L_★ ≥ 20 L_⊙), the polarization fraction in the outflow cavity walls can reach values of ~15%, with substantial amount of polarized flux, allowing detection after filtering the emission from large scales, as expected when observing the model with an interferometer such as ALMA (see Sect. 3.2). As the luminosity (L_★) increases, the efficiency of the radiative torques increases as well, and more and more grains get aligned; in other words, a_align decreases with increasing irradiation field strength. The equatorial mid-planes remain poorly irradiated compared to the outflow cavities and cavity walls, and exhibit the highest values of a_align in the core. Indeed, given the high dust column density in the mid-planes, the radiation field is reprocessed very quickly toward lower energy photons. This is not totally hindering the production of polarized flux in the mid-planes: while the signal is significantly lower than in the cavity walls, the magnetic field lines appear very organized in the mid-planes, which allow the emission to remain significantly linearly polarized.

3.2 Dust grain alignment efficiency

To derive the average grain alignment efficiency in our models, we use the method presented in Le Gouellec et al. (2020), which consists in using the product 𝒮 × Ƥ_frac as an estimation of the dust grain alignment efficiency, where 𝒮 is the dispersion of polarization angles in the plane of the sky (computed over the 8 neighboring pixels at each location, with maps that have been re-gridded to 4 pixels per beam area) and Ƥ_frac is the polarization fraction (sensitive to the grain alignment efficiency, and disorganization of the magnetic field in the line of sight). To do so, we derive the distributions of 𝒮 × Ƥ_frac values as a function of column density in our models, before and after spatial filtering with the CASA simulator, and compare them with the 𝒮 × Ƥ_frac distributions obtained in Le Gouellec et al. (2020).

There are many different cell sizes in the RAMSES simulations due to the use of AMR, which degrades the spatial resolution in some regions of our radiative transfer maps. We made sure that the resolution achieved by the artificial ALMA observations on the models (i.e., 100 au when observed at 0.87 mm and 140 au at 1.3 mm) is larger than the size of the AMR cells in the 2000 au central region of the synthetic observations. Beyond this central region, the AMR cell sizes are larger, and the statistics cannot be computed. However, in those regions the emission is anyway not recovered after spatial filtering, as it is shown in Fig. 4. In the filtered maps, we apply the same criteria pixel selection as in Le Gouellec et al. (2020), which is based on a cutoff in Stokes I, whose value is obtained when the signal to noise ratio in polarized intensity is, on average, higher than 5.

We present the evolution of 𝒮 × Ƥ_frac values as a function of the normalized column density in our set of models III in Fig. 6, sets IV and V in Fig. 7, before and after filtering, for the two wavelengths of observation, 0.87 and 1.3mm. Alongside those results, we plot for reference, in gray, the distributions obtained in Le Gouellec et al. (2020) from several ALMA observations, where we split the 0.87 and 1.3mm observations. The initial distribution of Le Gouellec et al. (2020) merged 0.87, 1.3, and 3 mm observations. However, given the results obtained in Fig. B.2, splitting the 𝒮 × Ƥ_frac distribution by wavelength is necessary. Despite the fact that most of the ALMA observations were realized at 1.3 mm, four cores were observed at 0.87 mm (Serpens Emb 6, Emb 8, Emb 8(N), and NGC 1333 IRAS4A). We note that, on average, the dust grain alignment efficiency, as traced by 𝒮 × Ƥ_frac, is higher at 1.3 mm than at 0.87 mm in our radiative transfer modeling (comparing the models implementing the same parameters) and in the ALMA observations. This is most likely due to higher polarization fraction values at 1.3 mm (this increase in polarization fraction with wavelengths for our different models can be seen in Fig. B.2). This evolution of polarization fraction with wavelength can be caused by several effects: temperature distribution in the lines of sight, opacity effects, and largest aligned grains contributing more to the polarized dust emission at longer wavelengths. In addition, the plots from Figs. 6 and 7 exhibit quite a large scatter in the distribution, and the column density maps in our models tend to be more peaked compared to ALMA observations¹. Therefore, we tend to discuss the global evolution of a given 𝒮 × Ƥ_frac distribution, as the radiative transfer parameters and the wavelength of observations vary.

We find that increasing the ratio f_high−J, the maximum grain size a_max, and the central luminosity L_★ implemented in POLARIS strongly influences the grain alignment efficiency in the envelope, typically causing the efficiency to increase (see also Appendices B and C). Nevertheless, the grain-alignment efficiency values calculated from the simulations are on average lower than the values from ALMA observations of Class 0 cores, if luminosities ≤20 L_⊙ are considered. We note that our MHD model may exhibit higher inner volume gas density values (given the initial mass of the model of 30 M_⊙) compared with what is expected in typical low-mass Class 0 protostars, thus justifying the need of high luminosity to match the observed grain-alignment efficiency.

Thanks to the multiple sets of parameters we adopt to run our radiative transfer calculations (from I to VI), one can infer the optimal parameters that yield the values of grain alignment efficiency most similar to those observed in the ALMA observations. Assuming f_high−J = 1 and a_max = 10 µm, Fig. 6 shows that only the implementation of large values of L_★ (i.e., 20, 50, and 100 L_⊙) produces mean values of 𝒮 × Ƥ_frac at 0.87 mm high enough that they are potentially able to match the values from ALMA observations. At 1.3 mm, 𝒮 × Ƥ_frac values from ALMA observations are notably higher than our models.

The results presented in Fig. 7 argue that low values of f_high−J and a_max (i.e., 0.25 and 2 µm, respectively) hinder efficient grain alignment in the envelope. In other words, we cannot reproduce the average values of dust grain alignment efficiency obtained with ALMA observations if we do not implement a large fraction of perfectly aligned dust grains (i.e., f_high−J close to 1), and if we do not include a population of large dust grains in the core (≥10 µm). Figure 7 thus justifies the values of the two fixed parameters of f_high−J and a_max we have chosen in Fig. 6.

The sets of radiative transfer calculations presented here provide hints that relatively high irradiation (i.e., with central luminosities L_★ higher than 20 L_⊙ but also large grains, ≥ 10 µm, and f_high−J close to 1) needs to be implemented in our model in order to increase the estimated values of grain alignment efficiency and match ALMA observations. Overall, the ALMA values 𝒮 × Ƥ_frac remain larger than our models, suggesting that further investigations are necessary to improve our understanding of grain alignment theories and/or our knowledge of the physical conditions in the inner region (i.e., ~ 10–1000 au) of Class 0 protostellar cores.

Fig. 5 Effects of the central luminosity, L★, on the radiative transfer results at 0.87 mm, with the fixed parameters of amax = 10 µm and fhigh−J = 1. This corresponds to set III of Table 1. Each column is one radiative transfer run of POLARIS, with L★ = 1, 5, 20, 50, and 100 L⊙. Each row is a quantity provided by the radiative transfer, from the first to the sixth row: total intensity Stokes I, polarized intensity (P), polarization fraction (Ƥfrac), 2D temperature slice obtained at the center, 2D radiation field slice obtained at the center (urad/uISRF), and the 2D slice of the aalign parameter obtained at the center.

Fig. 6 Comparison of the evolution of 𝒮 × Ƥfrac as a function of normalized column density, , between ALMA observations and our models, before and after spatial filtering. The left (right) panel corresponds to the results of the radiative transfer performed at 0.87 mm (1.3 mm). The distributions of 𝒮 × Ƥfrac values from our models are in different tints of blue, while the distributions of 𝒮 × Ƥfrac values from the ALMA observations presented in Le Gouellec et al. (2020) are in gray (where we split the 0.87 and 1.3 mm observations). The solid lines are the mean of the 𝒮 × Ƥfrac values in a given bin of normalized column density, and the corresponding shaded area represents the standard deviation. These two plots correspond to set III of the models (see Table 1), i.e., we fix fhigh−J = 1 and amax = 10 µm, and we vary L★ in the range 1.0, 5.0, 20, 50, 100 L⊙. The dashed colored lines correspond to the filtered models, and the solid lines are from the non-filtered models.

Fig. 7 Comparison of the evolution of 𝒮 × Ƥfrac as a function of normalized column density, , between ALMA observations and our models, before and after spatial filtering. Same as Fig. 6 for sets IV and V (see Table 1), where we use smaller values for amax (2 µm, top row) and for fhigh−J (0.25, bottom row).

3.3 Dust grain rotational disruption

Figures 6 and 7 encourage us to consider relative high protostellar luminosity (with central luminosities L_★ higher than 20 L_⊙) so that the radiation field is sufficiently high in the core, in order to reproduce the level of grain alignment efficiency measured with ALMA. However, increasing the efficiency of radiative torques can in turn cause the rotational energy to overcome the tensile strength of dust grains, leading to the disruption of grains into fragments. Given that the large aligned grains are the ones that are efficiently aligned, dust grain disruption could deplete the large grains necessary to produce the polarized dust emission measured inside star-forming cores with ALMA.

This mechanism is called the radiative torque disruption (RATD; Hoang et al. 2019). Already detected in the dense star-forming cores of Orion KL and ρ Ophiuchus (Tram et al. 2021a,b), it can be triggered by RAT (Hoang et al. 2019; Hoang & Tram 2020; Hoang 2020) or mechanical torques (Hoang & Tram 2019; Hoang & Lee 2020). Hoang et al. (2021) derived an analytical model to investigate whether RATD is at work inside star-forming cores. Here, we attempt to use the theoretical considerations of the RATD mechanism and compute them using the results of our radiative transfer modeling, in order to check if our claims – that the presence of large grains in the inner envelope and the need of a high central luminosity, – lead to the disruption of large grains.

A dust grain of mass density ρ, of size a, rotating at the angular velocity ω experiences a tensile stress of S = ρa²ω²/4. We designate the tensile strength of dust grains as S_max. The critical angular velocity above which dust grains are rotationally disrupted is given by (1)

If the maximum rotation speed induced by RAT ω_RAT exceeds ω_disr, dust grains are rotationally disrupted. Given the two regimes of the dust grain alignment efficiency , defined by the relative values the grain size a has with respect to the mean wavelength of the radiation field spectrum received by dust grains, one can define the grain size interval [a_disr; a_disr,max], inside which the aligned dust grains are rotationally disrupted. The respective values of a_disr and a_disr,max, are obtained developing the equation ω_RΑΤ = ω_disr, in the two aforementioned regimes of the dust grain alignment efficiency. From Hoang et al. (2021), we have (2)

and (3)

where γ is the anisotropy of the radiation field, u_rad is the radiation field, is the mean wavelength of the radiation field, n_H is the gas density, T_gas is the gas temperature, m_H is the mass of a hydrogen atom, and F_IR is the relative importance of rotational damping by infrared emission (due to the emission of infrared photons emitted by the grain, which reduces the grain’s angular momentum; see Draine & Lazarian 1998; we adopt the relation of Hoang et al. 2021) with respect to the damping by gas collision². We adopt ρ = 3 g cm⁻³ (this corresponds to the grain density of the dust model we adopt in POLARIS), and use the temperature, radiation field strength, and local gas density from each of our models. Both γ and are derived in each cell of the grid during the radiative transfer calculations. However, is not present in the output files of POLARIS. Therefore, to calculate , we use the relations in the Sect. 4 of Hoang et al. (2021), who derived, thanks to an analytical model, the mean wavelength as a function of the gas column density of a typical protostellar envelope. In the 2D slice of each of our models, we derive the column density “seen” by each cell by summing the gas density on the straight line separating the cell from the center (this apparent column density is dominated by the central ~500 au region). We note that this is an approximation, as the photons received by a given cell have been reprocessed, and thus do not come “straight” from the protostar. Our values of are thus lower-limits. As both a_disr and a_disr,max are proportional to , more realistic could thus shift the size range of rotationally disrupted grains toward larger values. However, the power-law dependence of with the apparent column density is weak: ~0.65.

We calculate the values of a_disr and a_disr,max for a given value of S_max. The grains whose size falls within the size interval of aligned grains [a_align; a_max] as well as the interval of rotationally disrupted grains [a_disr; a_disr,max] can be considered as rotation-ally disrupted. Now only one parameter remains to be set for us to evaluate the Eqs. (2) and (3): the grains tensile strength S_max. The tensile strength of interstellar dust is uncertain because it depends on both the grain structure (i.e., compact versus composite) and the grain composition, both of which are not strongly constrained quantities in Class 0 envelopes. We explore different grain structures, considering range 10⁵−10⁹ erg cm⁻³ for S_max. However, given that the largest grains in protostellar envelopes that undergo collisions and coagulation during the collapse are expected to be aggregates (Jones 2016), our results will focus on the S_max = 10⁵ erg cm⁻³ (see Hoang 2019; Hoang et al. 2021). The largest value for a_max that we considering in our work is 50 µm, which implies that significant grain growth have occurred from the typical ISM a_max value of 0.5 µm. We note that the presence in protostellar envelopes of such large aggregates and underlying efficient grain growth, on which the resulting submillimeter polarization is highly dependent (Valdivia et al. 2019), is still debated (see the recent works by Guillet et al. 2020; Bate 2022; Silsbee et al. 2022 that put constraints on early dust grain growth given the initial conditions such as turbulent velocities, timescales, and gas density).

From our calculations, the intervals [a_align; a_max] and [a_disr; a_disr,max] can overlap for the largest grains, if the irradiation and tensile strength are sufficiently high and low, respectively. Figure 8 presents the results of the computation of a_disr and a_disr,max, where we vary the protostar luminosity L_★, for a tensile strength value of 10⁷ erg cm⁻³. The last row indicates whether both intervals overlap. As the radiation field strength increases, we notice that more grains get aligned (i.e., a_align decreases), but at the same time a_disr decreases and the region where a_disr is lower than a_max becomes larger. In other words, increasing the luminosity at the center causes larger and larger depletion of large grains. This depletion occurs mainly toward outflow cavities. However, as we saw earlier when implementing maximum grain size of 2 µm for example, depletion of large grains considerably reduces the amount of polarized intensity in the core. Given that only grains larger than ~10 µm can produce a sufficient submillimeter polarized emission, our models require the rotational disruption of large grains to be limited in order to reproduce the high grain alignment efficiency revealed by ALMA observations. We also note that cells with high-density material along their lines of sight to the protostar (i.e., their equatorial mid-planes; see Fig. 1) seem not to be affected by RATD. Therefore, they will be the last regions of the envelope to suffer depletion of large dust grains when the irradiation increases.

In addition, we notice that there is a slight asymmetry, from the east to the west part of the core, in the plot displaying the region affected by RATD in Fig. 8. This is due to the structure in density of the very inner part of the circumstellar disk, as can be seen in Fig. 1. When deriving the column density “seen” by a given cell, the densest cells located at the center govern the resulting column density. In the slice density map, small asymmetries toward the circumstellar disk result in differences in , from the east to the west part of the map. At a given radius from the center, this, in turn, can cause the a_disr parameter to change from one side of the core to another. Therefore, inhomogeneities in density among the densest regions of the disk can be responsible for a region to be differently affected by RATD compared to one another. This phenomenon can be a potential explanation for asymmetric polarized intensity, such as that seen in B335, where strong polarization is only seen in the northern equatorial plane. One side might be more affected by RATD, and thus depleted in large grains that are needed to produce detectable polarized dust emission. This remains, however, a hypothesis, as different dust grains characteristics between the northern and southern portion of the equatorial plane have not been probed.

Finally, we investigate the locations affected by RATD in our models for several values of tensile strength S_max, 10⁸, and 10⁹erg cm⁻³, in Fig. 9. The smaller the implemented value of S_max, the lower the luminosity required to disrupt large grains in a given region. The RATD mechanism clearly affects outflow cavities for L_★ ≥ 20 L_⊙ if S_max ≤ 10⁷ erg cm⁻³ (i.e., if grains are composite), with grains ≥ 5 µm being depleted. However, for L_★ > 20L_⊙, if S_max > 10⁸ erg cm⁻³ (i.e., if grains are extremely compact), 10–50 µm grains survive throughout all the envelope.

The RATD mechanism seems to affect one specific part of the envelope. This mechanism, if it occurs, could cause the depletion of large, moderately compact grains (i.e., a ≥ 5 µm and S_max ~ 10⁵erg cm⁻³) in outflow cavities that are directly exposed to a sufficiently strong radiation field (L_★ ≥ 20 L_⊙). As we showed earlier, this may be problematic because both large grains and significant irradiation need to be implemented in order to reproduce the dust grain alignment efficiency values observed by ALMA. Assuming that RAT is the mechanism at play in the inner regions of Class 0 protostellar cores and that our implementation of the RAT theory is accurate enough, the combination of synthetically observed MHD models (with varying grain sizes and luminosities) alongside ALMA observations could constrain the tensile strength of the dust grains responsible for the polarized dust emission we observe.

Fig. 8 Effect of the central luminosity, L★, on the RATD of dust grains, at 0.87 mm, with the fixed parameters of amax = 50 µm and fhigh−J = 1. Each plot is a 2D slice taken at the center of the core in our radiative transfer results. Each column is one radiative transfer run of POLARIS, with L★ = 1, 5, 20, 50, and 100 L⊙. The first row is the aalign parameter, i.e., the dust grain size above which dust grains are considered aligned in our radiative transfer calculations. The second and third rows are the adisr and adisr,max parameters, respectively. In RATD theory, these two values correspond to the window in grain size, inside which dust grains are rotationally disrupted. The last row indicates whether dust grains are rotationally disrupted in our models for Smax = 105 erg cm−3, i.e., if the intervals [aalign; amax] and [adisr; adisr,max] overlap. For a given pixel, if there is an overlap, the pixel is black (i.e., grains are disrupted). If it is not the case, it is white (i.e., grains are not disrupted).

Fig. 9 Effects of the central luminosity, L★, on the RATD of dust grains at 0.87 mm, with the fixed parameters of amax = 50 µm and fhigh−J = 1. Each plot is a 2D slice taken at the center of the core in our radiative transfer results. Each column is one radiative transfer run of POLARIS, with L★ = 1, 5, 20, 50, and 100 L⊙. Each row indicates whether dust grains are rotationally disrupted in our models, i.e., if the intervals [aalign; amax] and [adisr; adisr,max] overlap. For a given pixel, if there is an overlap, the pixel is black. If it is not the case, it is white. The results from each row are obtained for a given value of grain tensile strength, Smax, i.e., 105, 106, 107 (the same plots as the bottom row of Fig. 8), 108, and 109 erg cm−3.

4 Discussion

Our radiative transfer analysis highlights that in the environment of Class 0 protostellar envelopes, the radiative feedback from the accretion/ejection processes has an important role in the efficiency of grain alignment, and potentially on dust grain evolution. One needs to reconcile the constraints we obtained from these results in terms of irradiation (and grain size and paramag-neticity) with the phenomenon of grain rotational disruption and any time variation of the accretion luminosity.

4.1 The radiation field in a Class 0 core: The photons’ sources and propagation

The radiation field strength and spectrum is at the center of our study aiming at understanding the mechanism responsible for the dust polarization we observe. In the radiative transfer calculations we presented, the central protostellar embryo represented by a sink particle is the only source of high energy photons. We describe here below how this can represent a caveat.

4.1.1 Accretion luminosity

In the inner regions of a Class 0 protostellar core, the radiative energy come mainly from the accretion energy radiated away that is constantly reprocessed along its pathway. In reality, the densest regions of highest opacities in the core, whose typical size is smaller than the spatial resolution of our gridding pattern, are responsible for most of the reprocessing of the high energetic photons from the accretion shock. We note that, UV photons would be reprocessed faster than X-ray photons throughout the envelope (Stäuber et al. 2004, 2005). In our models, our implementation of the spectrum of the energy radiated away from the sink due to the mass accretion onto the protostellar embryo is incomplete. However, we notice that the high-energy photons are quickly reprocessed toward the central high density cells surrounding the sink particle, in such a way that changing the blackbody spectrum parameters (fixing the luminosity) does not affect the radiation field spectrum throughout the core.

Baraffe et al. (2009) developed an analytical model of the accretion thermal efficiency that accounts for the fraction of accretion luminosity that is radiated away, with respect to the energy absorbed by the protostar. This model also considers details of the accretion process, such as how much energy can be stored in viscous heating or rotational energy when accretion is happening via an accretion disk (Hartmann et al. 1997; see also Jensen & Haugbølle 2018; Kuffmeier et al. 2018). We attempted to use the smooth step function Baraffe et al. (2009) proposed modeling the thermal efficiency of the accretion. In our case, using the method of Baraffe et al. (2009), the accretion rate on the sink provides too high values of luminosity compared to what is tentatively measured in intermediate-mass protostars. This can be explained by the spatial resolution of our MHD simulation (~5 au) being too low to properly take into account the storage of material in a disk, which can lower the mass accretion rate onto the protostellar embryo.

In order to reproduce the averaged grain alignment efficiency measured in ALMA observations with our model, our radiative transfer results show that a black body spectrum of luminosity ≥20 L_⊙ is necessary. Indeed, the protostellar photospheric luminosity (≤1 L_⊙) alone cannot result in a significant polarized dust emission. On average, the luminosity of protostellar cores is, however, significantly below the expected values. This is the so-called luminosity problem (Kenyon et al. 1990; Kenyon & Hartmann 1995; Young & Evans 2005; Evans et al. 2009). Therefore, requiring high radiative energy in all the cores that exhibit the high grain alignment efficiency quantified in Le Gouellec et al. (2020) can be problematic because these values may be above the values observed in protostars. Nevertheless, Krumholz et al. (2012) proposed that a significant fraction of the accretion energy can escape through outflow cavities. This is also what we found when plotting the 2D slice of the radiation field in our models. The low values of measured bolometric luminosity, could thus be reconciled with the expected mass accretion rate onto Class 0 protostars, if the luminosity from the accretion is not isotropically distributed throughout the envelope.

4.1.2 Shocks along outflow cavities

The only source of photons in our radiative transfer models is from the center, where the sink particle is located. However, from the temperature maps of the RAMSES simulations, we can see that the mechanical energy deposited by the jet is also responsible for heating the gas. This potential source of UV photons can directly heat the dense walls of the cavity, in contrast to the photons from the accretion that have already been reprocessed well before reaching the cavity walls. Therefore, the presence of shocks can modify the mean wavelength received by dust grains responsible for producing the polarized dust emission. In particular, if the radiation field spectrum contains more high energy photons thanks to these additional sources, smaller grains could get aligned; these grains are less likely to be affected by RATD. These suggestions are, however, not included in our models, as the gas dynamics (and associated shocks and radiative energy) are not taken into account in POLARIS (see the theoretical work of Hoang & Tram 2019). The required subsequent investigations would require the use of the modeling of shocks in cavities in order to identify the source location and spectrum of the photons produced in shocks. Estimating their energetics and penetration path-length could then inform us how they can contribute to the radiative torques applied to dust grains. One could also implement self-consistently in the radiative transfers additional sources of photons along the interaction zones between the jet and the envelope, based on 1D modeling of self-irradiated shocks calculated with the appropriate initial conditions (Lehmann et al. 2020). We also note that such shocks in outflow cavities may affect the size distribution of dust grains via grain shattering (Jones et al. 1996; Guillet et al. 2009).

4.2 Will the dependence with radiation be canceled by RATD?

To reproduce the dust grain alignment efficiency inferred from ALMA observations of Class 0 cores, large grains (≥10 µm) need to be present in the irradiated inner envelope. However, while this irradiation increases the dust grain alignment efficiency, it can also trigger grain rotational disruption, and deplete large grains in the affected regions. This phenomenon is dependent on the density along the photons’ pathway, the local gas density, and the radiation field. Therefore, the regions that are detected in dust polarization are also the regions “resisting” RATD (see Fig. 8), which are irradiated and dense enough, such as cavity walls. Now, these considerations are made assuming a given tensile strength S_max of dust grains. As presented in Sect. 3.3, grains that are more compact (higher S_max) are better able to resist RATD. On the contrary, grains made of aggregates with low tensile strength are more easily destroyed by RATD. As significant grain growth is necessary to justify the observed polarized dust emission in Class 0 protostellar cores, one could constrain the tensile strength of dust grains studying in what conditions and via what processes they grew.

Constrains on the grain internal structures, such as the tensile strength, could thus be obtained via observations of dust emission (constraining a_max) and dust polarization (constraining a_disr with dust polarization models). In order to study the regime where RATD can disrupt dust grains within typical ISM irradiation and density conditions, Hoang (2019) explored the dynamical constraint for several interstellar dust models: a contact binary grain model made of a silicate particle stuck to a carbonaceous particle, a model of composite grains made of subfragments of silicate and carbonaceous materials, and a model of a silicate core with an amorphous carbon mantle. For a range of irradiation and density conditions, the authors obtained a variety of rotational disrupted size a_disr given the internal structure, and tensile strength derived from these dust models, highlighting the importance of understanding grain internal structure to analyze the dust polarization observations.

In dense cores, large grains (≥ 1 µm) are expected to have grown via coagulation processes such as grain sticking, and are thought to have composite structures, or compact structures with a high tensile strength (S_max = 10⁷ erg cm⁻³) with potentially a core-thick ice mantle structure (e.g., see McClure 2009; Andersen et al. 2014; Poteet et al. 2015 for observational supports of ice mantles). However, dust growth and the physical conditions of the environment where grains have grown are not completely observationally constrained, and neither is the grains’ internal structure. It is likely that grains become fluffy aggregates via low-velocity collisions, but also can go through a phase of compaction at higher velocities (Ormel et al. 2009). In the interior of a protostellar core, dust grains from the high temperature and high density regions of the inner circumstellar disk may have different internal structures compared to the grains that have grown on longer timescales within dense and cold infalling structures. Recently Garcia & Gonzalez (2020) developed a dust model that follows the grain porosity (a direct proxy for the dust grain tensile strength) during grain growth processes (in this case the formation of porous aggregates) occurring in protoplanetary disks. The high density of the disk may cause the dust aggregates to be compressed by static compression due to ram pressure of the disk gas, decreasing the grains’ porosity and thus increasing the grains’ tensile strength (Kataoka et al. 2013a,b). Therefore, if these large dust grains are efficiently lifted up to populate outflow cavity and cavity walls (Wong et al. 2016; Tsukamoto et al. 2021), their structure might allow them to resist RATD more strongly in such irradiated regions.

In other words, if RAT is the mechanism aligning dust grains in cores, and if our models, which require us to implement large dust grains and high irradiation conditions, are accurate, then in order to protect large grains from being depleted by RATD, we must explain the high tensile strength of the corresponding grains. To that end, it would be required to characterize the physical conditions experienced by grains both in the inner disk and in the cavities in order to constrain their compactness using dust grain evolution models (see Garcia & Gonzalez 2020 for example).

4.3 Effects of the radiation field strength on the grain precession timescales

4.3.1 Radiative precession timescales

Finally, we attempted several timescale comparisons to determine whether our results are coherent with one another. The main part of the irradiation in the envelope that we measure in our models directly comes from the accretion shock onto the central protostellar embryo, which has been suggested to be episodic (Vorobyov & Basu 2005, 2006, 2010; Baraffe et al. 2009, 2012, 2017; Fischer et al. 2022). During an accretion burst, the luminosity of a protostellar object can increase by a factor of a few, up to ~ 10–100. We need to determine how fast the temperature in the envelope adapts to a change in the central luminosity and how fast dust grain rotation, precession, and potential disruption react to a change in irradiation and temperature. Figure 10 presents the effects of the central luminosity value, L_★, on the different timescales concerning the precession, alignment, and disruption of dust grains. We present and discuss the results of the this timescale analysis below.

In a core, the cooling time is commonly much shorter than the free-fall time (Eq. (3) of Hennebelle et al. 2020, which consist in calculating the radiative cooling time of a cell at a given distance, with its density, opacity, and temperature). We can consider that gas and dust temperature will become equal with one another in a few gaseous collisional timescale τ_gas (see Fig. 10). In an outflow cavity wall of the fiducial radiative transfer, at ~ 1000 au from the center, the cooling time is ~ 1 yr, and the free-fall time is ~6 × 10³ yr. As a consequence, it is reasonable to consider that the temperature adjusts instantaneously to a change of the central source luminosity, and that the gas dynamics is not affected. Therefore, if there is an accretion burst, one can expect the dust grains to immediately experience the different irradiation conditions. The accretion burst duration is thought to be ≤ 100–200 yr in protostars (Vorobyov et al. 2013), and the typical interval between two bursts varies within 10² and 10⁴ yr (different observational techniques have found a large variety of typical interval between bursts; e.g., see Contreras Peña et al. 2019; Jørgensen et al. 2020; Lee et al. 2021b; Park et al. 2021; Fischer et al. 2022). However, we note that in strongly irradiated regions the gas is likely warmer than the dust, even at high density (see for example the work of Koumpia et al. 2015 who estimated a higher gas temperature by ~15 K toward the photon dominated region S 140). This may be due to shocks, and/or intense UV heating of the gas. Therefore, the temperature of the gas in outflow cavities and cavity walls may be higher than the dust temperature. This would affect our estimation of the gaseous damping timescale τ_gas (defined below), which scales with .

The timescale for a previously nonaligned grain to spin up via RAT up to a rotational velocity that would allow its alignment (at least three times the thermal velocity) is about the radiative precession timescale, τ_rad, if the grain is directly driven to a high-J attractor point (Lazarian & Hoang 2007), given as (4)

where a₋₅ = a/10⁻⁵ cm, the dust grain density, the grains’ aspect ratio, , is the anisotropy of the mean radiation field (see Eq. (20) of Tazaki et al. 2017, or Eq. (53) of Hoang et al. 2022), is the RAT efficiency (see the relation 10 in Hoang et al. 2021), and is the mean wavelength of the radiation field³. The alignment timescale for a grain directly driven to a high-J attractor point corresponds to a fast alignment, while the fraction 1 − f_high-J of grains that are driven to a low-J attractor point experiences a slow alignment, as gaseous bombardments gradually move grains from the low-J attractor toward a high-J state (Hoang & Lazarian 2008, 2016b; Lazarian & Hoang 2021). Such a slow alignment is typically achieved after 5–10 τ_gas, where τ_gas is the collisional gaseous damping time of the grain (the infrared emission damping and plasma drag are neglected): (5)

where is the thermal velocity, I₁ is the principal grain moment of inertia, and µ is the mean molecular weight per hydrogen molecule. Nevertheless, the timescale for slow alignment is reduced when considering super-paramagnetic dust grains (Hoang & Lazarian 2008, 2016b; Lazarian & Hoang 2021; Chau Giang et al. 2023), but can also increase for increasing radiation field strength because stronger radiative torques push grains to the low-J attractor point (Lazarian & Hoang 2021). Therefore, the super-paramagneticity of dust grains plays an important role in the alignment timescale because (1) it governs the fraction of grains f_high−J that can be driven to a high-J attractor point (and thus of grains that can experience a fast alignment), and (2) it diminishes the grain alignment timescale for grain experiencing slow alignment. The radiation field improves the strength of RAT and thus the maximum value of grain angular momentum (which can become high enough to trigger RATD), but forces grains to the low-J attractor point that will experience slow alignment. An increase in the irradiation field during an accretion burst will thus generate an increase in the dust polarized emission on a timescale corresponding to the one of the slow alignment for the 1 − f_high−J fraction of grains driven at low-J. As the increase in irradiation also diminishes the minimum size of grains rotational disrupted adisr, more grains could be rotationally disrupted. In the absence of high-J attractor point for grains with inefficient magnetic relaxation (low δ_m) the increase in RAT via higher irradiation may thus cause slow disruption of grains (on the order of the slow alignment, i.e., 5–10 τ_gas), because only low-J attractors would exist. This makes RATD unlikely for low-J aligned grains (Hoang & Tram 2020; Hoang et al. 2022). In our model, the observed grain alignment efficiency strongly suggests that most grains are perfectly aligned with f_high−J ≃ 1 (see Appendix C and Chau Giang et al. 2023). Therefore, in the type of environment we model, an increase in the radiation field strength can rapidly produce an increase in the dust polarized emission (i.e., on a timescale on the order of τ_rad).

The potential subsequent rotational disruption of grains that are supposed to be directly driven at high J also happens rapidly compared to the gaseous damping timescale, on the order of the timescale τ_disr (see Fig. 10), which is the characteristic timescale for fast rotational grain disruption (Hoang & Tram 2020; we note that this is much shorter than the shattering time by grain-grain collisions); it is derived as follows: (6)

where is the radiative torque average over the incident radiation spectrum (as for τ_rad, we adopt the value of from Hoang et al. 2021), and ω_disr is presented in Sect. 3.3. If grains have super-paramagnetic (iron) inclusions, such as a fast alignment, the dust polarization signal can have the time to increase during an accretion burst because τ_rad is generally much lower than the typical accretion burst duration. However, depending on the tensile strength of dust grains and luminosity reached during the burst, aligned dust grains could be rotationally disrupted if τ_disr becomes comparable or lower than an accretion burst duration. If otherwise slow alignment is occurring, slow disruption (on the order of the timescale for slow alignment, i.e., 5–10 τ_gas) can happen in the steady-state regime where a_disr ≤ a_max. However, during an accretion burst the increase in grain alignment and rotational disruption efficiency may not happen in the case of slow alignment because 5–10 τ_gas is on the order of an accretion burst duration. Even if only slow rotational disruption occurs, we note that RATD may remain an important process that can constrain early dust grain growth inside protostellar cores. Indeed, the typical timescale required to grow grains at sizes ≥10 µm inside a medium of n_H = 10⁵ cm⁻³ is on the order of 7–10 times the free-fall timescale (Ormel et al. 2009; Wong et al. 2016), namely, ~7−10 × 10⁵ years at such a density (while the density of the envelope of our model is n_H = 10⁴−10⁶ cm⁻³); this timescale is higher than the timescale of slow rotational disruption.

In addition, we note that an accretion burst may also influence the dust polarization position angle. Indeed, during such a burst the ionization increases and thus the coupling between the material and the magnetic field improves, such that a larger amount of material can dynamically affect magnetic field lines, which can be dragged more intensively inward. Such a change in dust polarization position angle given a change in magnetic field orientation would happen on a timescale on the order of the Larmor precession timescale, that is, lower than or on the order of the accretion burst duration, for grains ≥ 1 µm.

Fig. 10 Effects of the central luminosity, L★, on the alignment timescales in our radiative transfer models at 0.87 mm, with the fixed parameters of amax = 50 µm and fhigh−J = 1. Each column is associated with a central luminosity value, L★. While the plots in the three first rows represent the alignment timescales as a function of grain size, the plots in the last row represent the polarization fraction. The timescales are computed in specific locations, i.e., toward the equatorial mid-plane (first row), at an outflow cavity wall (second row), and in the middle of an outflow cavity (third row), in the central 2D slice of the model. These locations are indicated by the three colored circles in the polarization fraction maps, with the dark point for the first row (equatorial plane), the blue point for the second row (outflow cavity wall), and the green point for the third row (outflow cavity). In each plot of the first three rows, we present the Larmor precession timescale, τLarmor (with a dotted line for a grain without super-paramagnetic inclusions and with dot-dash and solid blue lines, respectively, for two values of Ncl, the number of atoms per cluster, 103 and 105), the radiative precession timescale, τrad, the collisional gaseous damping time of the grain, τgas (solid brown line), and the rotational disruption timescale, τdisr (for an Smax value of 105 erg cm−3, in solid dark yellow lines). We show the minimum grain size of aligned grains, aalign, with the vertical dotted green line. The range of grain size affected by rotationally disruption (adisr − adisr,max) for Smax = 105 erg cm−3 is shown by the horizontal dark line. At a given grain size, the shortest precession timescale between τLarmor and τrad dictates the direction of alignment for the grains aligned at low J, i.e., with the magnetic field or the radiation field.

4.3.2 Radiative alignment

Finally, if the radiation field is strong enough, grains can change their axis of alignment and become aligned with the axis of the anisotropic component of the radiation field (this is the k-RAT mechanism that we have already introduced; Lazarian & Hoang 2007). This occurs if the radiative precession timescale, τ_rad, becomes shorter than the Larmor precession timescale, which illustrates the efficiency of the interaction between the grain magnetic moment with an external magnetic field. Following Hoang & Lazarian (2014) and Tazaki et al. (2017), we have (7)

where is the magnetic field, and the grains’ paramagnetic zero-frequency susceptibility. The zero-frequency susceptibility for a super-paramagnetic grain is given by Curie’s law (Morrish 2001): (8)

where ϕ_sp is the fraction of atoms that are super-paramagnetic, and N_cl is the number of atoms per cluster. From measurement of GEMS, the typical value are ϕ_sp = 0.03 (Bradley 1994; Martin 1995; Goodman & Whittet 1995), and N_cl is expected to be N_cl = 10³−10⁵ (Kneller & Luborsky 1963; Jones & Spitzer 1967). We highlight that the need to implement f_high−J ≃ 1 in our models argue toward super-paramagnetic grains (see also Chau Giang et al. 2023). We also consider an ordinary paramagnetic grain, with (9)

where f_p is the fraction of atoms in the grain that are paramagnetic, evaluated at 10% (Draine & Weingartner 1996; Tazaki et al. 2017). For reference, Draine & Weingartner (1996) proposed; χ(0) values in the range 4 × 10⁻⁵−10⁻³ for ordinary paramagnetic grains. The evolution of τ_Larmor as a function of grain size for a super-paramagnetic grain is shown in Fig. 10. We note that, while magnetic relaxation participates to the alignment of grains’ angular momentum with magnetic field lines (Weingartner et al. 2021), grain alignment by magnetic relaxation alone is still inefficient in the absence of grain suprathermal rotation because of internal thermal fluctuations, even for super-paramagnetic grains (Hoang & Lazarian 2016b,a). For example, for typical protostellar envelope conditions Hoang et al. (2022) calculated that super-paramagnetic relaxation is negligible for grains ≳ 1 µm inside protostellar envelope where n ≳ 10⁷ cm⁻³. This makes Larmor precession the main agent in our models that controls the grain size parameter space describing which grains are aligned with magnetic field lines for grains ≳1 µm. Therefore, at a given grain size, the shortest precession timescale between τ_Larmor and τ_rad dictates the axis of alignment (i.e., the magnetic field or the anisotropic component of the radiation field).

Equation (4) is valid for grains enduring slow-alignment, aligned at the low-J attractor point. The radiative precession timescale for grains aligned in the high-J state is several orders of magnitude higher. Indeed, at the high-J attractor point the spin-up effect of RAT dominates, such as the Larmor precession is the strongest, and these grains are aligned with the magnetic field (see Sect. 5.3 of Hoang et al. 2022). As L_★ increases, the grain size at which the transition of alignment axis must occur decreases, and can reach ~1−10 µm for paramagnetic grains when L_★ ≥ 20 L_⊙, toward the most irradiated regions of the protostar (i.e., inside outflow cavities). In cavity walls and equatorial planes, where the dust polarization is detected, the transition would occur for paramagnetic grains ≥ 10 µm. Such grain sizes are viable in these environments. However, if aligned dust grains are super-paramagnetic with f_high−J ≃ 1, the grains whose size is eligible to the k-RAT mechanism would rotate rapidly enough to be aligned with the magnetic field, and be potentially affected by RATD where it is efficient (i.e., in cavities; see Figs. 9 and 10). The observed grain alignment efficiency is high, and points toward super-paramagnetic grains aligned with the magnetic field with f_high−J ≃ 1 (Le Gouellec et al. 2020; Chau Giang et al. 2023). Therefore, the k-RAT mechanism is unlikely in pro-tostellar environments because the low-J aligned grains eligible to k-RAT would produce too low polarization fraction values if they were to dominate the polarized dust emission.

With our approach, the most crucial model parameters that describe the alignment mechanism in a given protostellar envelope are thus f_high−J, L_★, and S_max. While f_high−J ≃ 1 and high irradiation (i.e., L_★ ≥ 20 L_⊙) are necessary to reproduce the observed grain alignment efficiency, one must also constrain the effect of RATD (the dependence of a_align with the radiation field strength is weaker than the dependence of a_disr). The polarized dust emission observed in Class 0 protostars is thus produced in dense enough regions, and/or from compact enough aligned grains, where RATD does not disrupt the large grains that are responsible for the dust polarization.

5 Conclusions and summary

While in Paper I an analysis of ALMA dust polarization and molecular line observations was presented, this paper explores in detail the physics of grain alignment occurring in Class 0 protostellar cores. To assess the role of the radiation field, we performed POLARIS synthetic observations of a MHD simulation of a collapsing protostellar core (which implements the radiation treatment method developed by Mignon-Risse et al. 2020 and the protostellar jet implementation of Verliat et al. 2022). We compared the grain alignment efficiency that we derived from ALMA polarization observations of Class 0 protostellar cores with the results from the radiative transfer modeling. The main results and conclusions of this paper are as follows:

1. We produced synthetic observations from MHD models of the protostellar evolution for an intermediate-mass core using POLARIS. We show that these synthetic polarized dust emission maps reproduce the enhancement of polarized dust emission along outflow cavity walls that has been observed toward several objects.

2. In the radiative transfer computations, we varied the luminosity of the central object, L_★, and studied how it favors dust grain alignment in the inner core. Comparing the dust grain alignment efficiency (as traced by 𝒮 × Ƥ_frac) derived in the radiative transfer calculations with those found in Le Gouellec et al. (2020) with ALMA observations, we find that a significant luminosity (i.e., L_★ ≥ 20 L_⊙) as well as super-paramagnetic grains (with a fraction of supra-thermally rotating grains close to 1) and large grains (with a_max ≥ 10 µm) are required for the dust grain alignment efficiency of the synthetic observations to reach observed values. The central luminosity values we implemented are dependent on the characteristics of the model we used (core initial mass and inner density structures) to reproduce the dust polarization.

3. One caveat of our modeling is the nature of the radiation field spectrum. The potential UV photons produced in shocks in jets and/or outflows and along the cavity walls and the temperature increase triggered by the shocked wind are not implemented in the radiative transfer calculations.

4. We implemented the equations of the RATD phenomenon in the radiative transfer calculations. This phenomenon triggers the rotational disruption of the largest aligned dust grains, which can cause a depletion of the large grains that are necessary to produce the observed polarized dust emission. Varying the irradiation and the tensile strength of dust grains, we find that, in our models, RATD can occur in outflow cavities for central luminosities of L_★ ≳ 20 L_⊙ and grain tensile strengths of S_max ≲ 10⁷ erg cm⁻³. The densest regions (i.e., outflow cavity walls and, especially, the equatorial mid-planes) seem more protected and are the last regions of the envelope to have their populations of large dust grains depleted by RATD.

5. Our models require both large grains and high irradiation, which raises the questions of where large grains formed in Class 0 protostellar cores, and whether their tensile strengths are high enough to resist rotational disruption. These two questions are intrinsically linked, given that the environmental conditions (such as temperature) in which large grains have grown can dictate their structure (compact versus composite) and thus their tensile strength.

6. The physical conditions that favor the polarized dust emission detected toward equatorial planes and infalling structures such as accretion streamers still represent an open question. One needs to explain how large grains (a ≥ 10 µm) form in the envelope, and how they resist disruption via RATD.

7. We studied the characteristic timescales of the grain alignment physics and rotational disruption of grains in an attempt to understand the effect that highly variable accretion rates would have on the dust grain alignment conditions. In the environment we model, during an accretion burst or a steady-state phase of high luminosity from the protostellar embryo (with L_★ ≤ 20 L_⊙), RATD could have enough time to disrupt the largest grains in irradiated regions, depending on their tensile strength. In such luminous conditions, slowly rotating grains ≥10 µm (at the low-J attractor point) could in theory preferentially align themselves with the radiation field (k-RAT grain alignment). However, given the high grain alignment efficiency observed in protostellar envelopes, having mainly rapidly rotating grains (at the high-J attractor point) aligned with the magnetic field is the most likely scenario. This suggests that the k-RAT grain alignment is not the main cause of the observed polarization dust emission in Class 0 objects.

Assuming RAT is the predominant cause of dust polarization in Class 0 sources, the next step would be to further investigate the characteristics of the dust grains populating Class 0 envelopes: in particular, their composition and helicity (see Hoang 2022). Right now we use the standard MRN distribution of grains with ISM chemical grain properties (silicates and carbonaceous grains). We may have to revisit the picture of dust grains in protostars, as their compositions and shape must favor the efficient grain alignment we observe (see for example the promising “astrodust” model; Draine & Hensley 2021; Hensley & Draine 2023). At this point, while RAT is favored by the investigations we present, one still needs to address the questions of grain growth, disruption, and structure. The radiation field strength seems to be an important agent for aligning grains in the inner core, but rotational disruption needs to be taken into account alongside grain structure to fully understand the role of irradiation in dust grain evolution. Another additional step in testing RAT theory involves the implementation of additional grain alignment theories in dust radiative transfer codes such as POLARIS: dust grains aligned with the radiation field, imperfect internal alignment (Hoang 2022; Chau Giang et al. 2023; Hoang et al. 2022), and mechanically aligned grains (Hoang et al. 2018, 2022; Reissl et al. 2023; requiring the implementation of dustgas drift in 3D protostellar collapse simulations; see Lebreuilly et al. 2019, 2020). These other grain alignment mechanisms must also be further characterized observationally (Pattle et al. 2021; Hull et al. 2022).

Acknowledgments

V.J.M.L.G. and C.L.H.H. acknowledge the ESO Studentship Program, and the guidance and support of Eric Villard. C.L.H.H. acknowledges the support of both the NAOJ Fellowship as well as JSPS KAKENHI grants 18K13586, 20K14527, and 22H01271. This work has received support from the European Research Council (ERC Starting Grant MagneticYSOs with grant agreement no. 679937). We thank B.-G. Andersson and Thiem Hoang for providing helpful explanations. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Facilities: ALMA, SMA. Software: APLpy, an open-source plotting package for Python hosted at https://aplpy.github.io (Robitaille & Bressert 2012), CASA (McMullin et al. 2007), Astropy (Astropy Collaboration 2018), POLARIS (https://portia.astrophysik.uni-kiel.de/polaris/; Reissl et al. 2016).

Appendix A Details about our synthetic observations and MHD model

This simulation belongs to a set of simulations aiming to study the evolution of protostellar cores, and especially the role of jets, initial turbulent and magnetic energy, accretion efficiency and luminosity, and protostellar outflow feedback. Constructing this model was motivated by our quest to reproduce the observable features of Class 0 cores, observed at high angular resolution with ALMA. While outflows and jets are ubiquitous in young protostellar cores (Bally 2016; Ray & Ferreira 2021), their complete modeling remains difficult. The main reason is that they are launched at sub-au scales, while the largest scale we need to implement to model an isolated collapsing core is ~0.1 pc (the tiniest resolution element of the simulation that uses the AMR is ~5 au). In addition, the physics occurring at the sub-au scales of the accretion demands very small time-steps in the calculations, which in turn, correspond to unrealistic computational times in current simulations. Magneto-centrifugal outflows can be launched (Gerrard et al. 2019); however, their energetics do not correspond to the protostellar jets observed in Class 0 cores. Therefore, a jet is implemented by hand at the creation of the sink particle. This jet ejects one third of the accreted material, with a speed of 66% of the escape speed, and with an opening angle of 30°.

The ALMA observations of Class 0 protostars used in Le Gouellec et al. (2020) to estimate an average grain alignment efficiency correspond to objects located in close-by low- and intermediate mass star-forming regions. They span a range of total envelope mass of ~1-40 M_⊙, and are likely of heterogeneous evolutionary ages. Not all of these cores exhibit pristine and symmetric envelope structures that we observe in our simulation. This suggests that the cores observed by ALMA have been formed in a variety of environments, that is, with different relative amount of large-scale turbulent, gravitational, and magnetic energies. We chose not to implement initial turbulence in order to obtain a pristine case of protostellar collapse, and understand more easily the propagation of photons in the radiative transfer calculations. If turbulence is implemented, significant precession of the outflow axis is caused by the sink particle’ spin axis changes, which are due to the turbulent accretion flow. This causes outflow broadening and additional entrained material (Rosen & Krumholz 2020). This could have limited our interpretations, especially as we aim to identify specific grain alignment environmental conditions, for example toward outflow cavities and equatorial mid-planes. The turbulence in a core is also responsible for delaying the formation of the sink in simulations. The simulation snapshot we have chosen corresponds to a relatively young star-forming core, compared to what is derived from observations for Class 0/I objects. The absence of initial turbulence (and thus the more rapid formation of the sink) can thus explain this young age. The initial mass-to-flux ratio, which varies significantly among star-forming cores (Hennebelle 2018), also plays a role in regulating the core dynamics during the collapse. We chose a reasonably magnetized core with µ = 5.

The luminosity from the accretion and the associated radiative pressure (Hennebelle et al. 2020; Rosen & Krumholz 2020), alongside the mechanical energy of the outflow/jet system, provide a radiative feedback in star-forming regions, and can be important to reproduce the star-forming activity of a given region. However, the effect of the accretion luminosity and subsequent radiative pressure on the evolution of the envelope are not implemented in this MHD simulation. We investigate the heating from the total protostellar luminosity (photospheric and accretion luminosities) in our radiative transfer calculations, and look at the resulting radiative field. From the jet and outflow launching models, the fraction of accretion energy that translates into radiative energy is still a free parameter in models (Ostriker & Shu 1995; Offner et al. 2009). Therefore, we chose to explore a broad range of protostar luminosities in Sect. 2.

The radiation field that escapes from the surface of the central protostellar embryo and the accretion shock are different in nature, compared to the radiation field escaping out of the sink particle. The radiation field escaping out of the sink particle is supposed to have already been reprocessed by high density/opacity material. In other words, at a given luminosity , the two parameters T_★ and R_★ evolve in such a way that the spectrum becomes less rich in high energy photons as they propagate within this high density material. However, the radiation field is expected to be reprocessed toward longer wavelength photons very quickly (i.e., at scales ≲ 1 au). Given the spatial scales we study, the shape of the blackbody spectrum thus does not have a significant impact on our results. Therefore, in the work presented here we fix R_★ and vary T_★ in order to vary L_★. At given luminosities, we investigated the impact of varying (R_★, T_★) and found that, indeed, no effects are seen on the grain alignment parameters or the dust polarization maps.

Appendix B Effects of the maximum grain size on the radiative transfer results

Valdivia et al. (2019) first presented the impact of the maximum dust grain size on polarized dust emission with POLARIS radiative transfers. In this Appendix, we present such effects on our model, exploring a_max = 0.5, 2, 10, 30, and 50 µm. This parameter is of great importance on the grain alignment efficiency within the dense regions of protostars. Because photons are rapidly reprocessed toward longer wavelength in the envelope of protostars, and radiative torques can only spin up dust grains whose size are larger than the wavelength of the anisotropic component of the radiation field, dust grains larger than ~ 10 µm must populate protostellar envelopes. Otherwise, the polarized dust emission would not be observable. We note that our study does not implement the maximum size of internally aligned grains. For the densest conditions of our model (i.e., n_H ≥ 10⁹ cm⁻³), grains larger than ~ 50 µm should not have efficient enough Barnett and inelastic relaxation processes to resist gaseous de-alignment. However, considering high irradiation conditions and a high level of iron inclusions (i.e., u_rad/u_ISRF ≥ 10⁶ and N_cl ~ 10⁵) can increase this upper limit (see Sect. 4 of Hoang et al. 2022).

Figure B.1 shows that changing the maximum grain size from 0.5 to 50 µm, with f_high−J = 1 and L_★ =20 L_⊙, leads to significant changes in the resulting dust polarization maps (see also Chau Giang et al. 2023). The heating of the envelope remains the same, that is to say, changing the last bins of the MRN distribution does not affect the temperature or the radiation field in our calculations. The maximum grain size has, however, a large impact on the emissivity properties of dust grains, as Stokes I derived at 0.87 mm increases with increasing a_max from 10 to 50 µm. However, the total intensity remains, on average, constant between 0.5 and 10 µm. The polarized intensity, which depends on how the fraction of aligned dust grains evolves, evolves differently than Stokes I does. While the irradiation remains the same in the cases of this figure, the evolution of the mean wavelength of the photons impinging onto dust grains could explain the evolution of the a_align. When the a_max increases, a_align also increases. However, the more the maximum grain size increases, the larger the fraction of aligned grains is, and the more numerous the grains that contribute to the polarized dust emission are. This explain why the polarized intensity increases with increasing amax. Finally, the increase in polarized intensity is different from the increase in Stokes I. At 0.87 mm, this leads to a maximum polarization fraction in the case implementing a_max = 10 µm.

Fig. B.1 Effects of the maximum dust grain size, amax, on the radiative transfer results at 0.87 mm, with the fixed parameters of fhigh−J = 1 and L★ = 20 L⊙. This corresponds to set I of Table 1. Each column is one radiative transfer run of POLARIS, with amax = 0.5, 2, 10, 30, and 50 µm. Each row is a quantity provided by the radiative transfer (from the first to the sixth row): total intensity Stokes I, polarized intensity (P), polarization fraction, 2D temperature slice obtained at the center, 2D radiation field slice obtained at the center (urad/uISRF), and the 2D slice of the aalign parameter obtained at the center.

Fig. B.2 Effects of the maximum dust grain size, amax, on the radiative transfer results at 0.87 mm, 1.3mm, and 3mm, with the fixed parameters of fhigh−J = 1 and L★ = 20 L⊙. Each column is one radiative transfer run of POLARIS, with amax = 0.5, 2, 10, 30, and 50 µm, while each row corresponds to a wavelength of observation for the radiative transfer, i.e., 0.87 mm, 1.3mm, and 3mm, as indicated on the first plot of each row. We only show the polarization fraction.

We show in Fig. B.2 the impact that the wavelength of observation (0.87, 1.3, and 3 mm) has on the polarization fraction maps, for the different values of maximum grain sizes of a_max = 0.5, 2, 10, 30, and 50 µm. Within this range of maximum grain sizes, the maximum average polarization fraction in the core is obtained at a_max = 10 µm when observing at 0.87 mm, at a_max = 30 µm when observing at 1.3 mm, and at a_max = 50 µm when observing at 3 mm. This is consistent with Valdivia et al. (2019), who mentioned that larger dust grains behave more like spherical grains, producing a smaller difference between Stokes I and the polarized intensity, thus producing a smaller polarization fraction. In other words, when observed as a function of maximum grain size a_max, the polarization fraction peak shifts to a higher value of a_max when increasing the wavelength of observations. In addition, at a given maximum grain size, the polarization fraction increases with increasing wavelength of observation, especially for a_max = 30 and 50 µm. This highlights the importance of multiwavelength observations to target the evolution of the polarization fraction. Given that the magnetic field morphology does not change, multiwavelength observations of a given core offer the possibilities to constrain the properties of the aligned dust grains. If our implementation of RAT accurately reproduce the alignment of grains in the inner core, one could fit the evolution of the polarization fraction with wavelength and constrain the maximum grain sizes in the region exhibiting polarized dust emission. This could in turn bring constrains on the minimum size of rotationally disrupted grains. However, one needs to mind that the aligned grains contributing to the polarized dust emission may come from different regions when observing at different wavelength. How vary the magnetic field morphology across these different regions of emission can also affect the multiwavelength evolution of the polarization fraction (see Valdivia et al. 2022).

Similar to the figures presented in Sect. 3.2, Fig. B.3 presents the evolution of 𝒮 × Ƥ_frac as a function of normalized column density , in ALMA observations and our models, before and after spatial filtering, for sets I and IV of radiative transfer runs. In this fig ure we vary the a_max values in the range 0.5, 2, 10, 30, 50 µm, with L_★ = 20 L_⊙. When f_high−J = 1, the implementation of large grains with a_max ≥ 10 µm allows the average level of the 𝒮 × Ƥ_frac of ALMA observations to be marginally reproduced. As noticed in Fig. B.2, we note also that the level of grain alignment efficiency traced by 𝒮 × Ƥ_frac decreases when a_max =50 µm, at 0.87 and 1.3 mm. When f_high−J = 0.25, none of the a_max values we implement produces high enough 𝒮 × Ƥ_frac values to match ALMA observations.

Fig. B.3 Comparisons of the evolution of 𝒮 × Ƥfrac as a function of normalized column density, , between ALMA observations and our models, before and after spatial filtering. The left (right) column corresponds to results of the radiative transfer performed at 0.87 mm (1.3mm). This is the same as Fig. 6 but for sets I and IV (see Table 1).

Appendix C Effects of f_high−J on the radiative transfer results

We refer to the recent study by Chau Giang et al. (2023) who thoroughly explored the polarized dust emission produced by super-paramagnetic grains. Yet we explore in this appendix the effects of changing the f_high−J parameter on the results of the radiative transfer calculations. In Fig. C.1 we vary the f_high−J parameter, from 0 to 1, with a_max = 10 µm and L_★ = 20 L_⊙. The polarized intensity and polarization fraction increase with increasing fraction of grains at the high-J attractor point, the fraction of grains (among those that are aligned) that can be considered perfectly aligned if their internal alignment is efficient (Hoang 2022; Chau Giang et al. 2023). The resulting polarization fraction is strongly affected by this parameter. Inside the central ~ 2000 au, f_high−J = 0 causes the polarization fraction to not be higher than 3%, while f_high−J = 1 causes the polarization fraction to reach ~ 12% in the outflow cavity walls.

Figure C.2 presents the evolution of 𝒮 × Ƥ_frac as a function of normalized column density , in ALMA observations and our models, before and after spatial filtering, for set II of radiative transfer runs. In this figure we vary the f_high−J values in the range 0.0, 0.25, 0.5, 0.75, 1.0., with L_★ = 20 L_⊙ and a_max = 10 µm. Only values of f_high−J close to 1 could produce high enough 𝒮 × Ƥ_frac values in order to match ALMA observations.

Fig. C.1 Effects of the fhigh−J parameter on the radiative transfer results at 0.87 mm, with the fixed parameters of amax = 10 µm and L★= 20 L⊙. This corresponds to set II of Table 1. Each column is one radiative transfer run of POLARIS, with fhigh−J = 0.0, 0.25, 0.5, 0.75, and 1.0. The top row is the polarized intensity (P), and the bottom row is the polarization fraction.

Fig. C.2 Comparisons of the evolution of 𝒮 × Ƥfrac as a function of normalized column density, , between ALMA observations and our models, before and after spatial filtering. This is the same as Fig. 6 but for set II (see Table 1).

Appendix D Correction for the optical thickness on the normalized column density of radiative transfer models

The distribution of normalized gas column density is more peaked in our model than in the distribution derived from ALMA observations in Le Gouellec et al. (2020). This difference can be explained by optically thick emission in the lines of sight located toward the total intensity peak in ALMA maps, precluding us from accessing the true value of the column density peak. This effect results in the distributions of 𝒮 × Ƥ_frac as a function of normalized column density, , for our models being shifted to the left compared to the ALMA distribution.

We applied a tentative correction to the column density map in our model to make the comparisons with the ALMA observations more reliable. This correction factor is computed from the optical depth τ obtained in each radiative transfer calculation (see Fig. 2). In each pixel, the gas column density is multiplied by , to take into account the error made in assuming an optically thin emission computing in ALMA observations. This correction is a first order approximation, and represents the error made assuming an optically thin emission while τ is no longer close to 0. Figure D.1 presents the radial distribution of τ_λ, and the distribution of original and corrected , at λ = 1.3 and 0.87 mm.

Fig. D.1 Correction for the optical thickness on the gas column density of our radiative transfer models, computed with amax = 10 µm at λ = 0.87 mm and λ = 1.3 mm. The dot-dash light blue line (solid dark blue line) shows the azimuthally averaged modified (original) gas column density as a function of radius. The solid dark red line shows the azimuthally averaged optical depth, τ, as a function of radius.

References

All Tables

All Figures

Fig. 1 Snapshot of our MHD simulation. Top panel: system seen edge-on. The color scale and white contours represent the gas column density. The yellow streamlines trace the orientation of the density-weighted average magnetic field along the line of sight. The white arrows represent the velocity field in a slice of the simulation centered on the sink particle. Bottom panel: the contours and white arrows are the same as in the top panel. The color scale represents a slice of gas mass density, centered on the sink particle. In both panels, the little white circle is the full width half maximum of a circular Gaussian kernel we use to slightly smooth the highest resolution of the AMR grid, for visibility.

In the text

Fig. 2 Dust polarization radiative transfer results of the fiducial case (amax = 10 µm, fhigh−j = 1, and L★ = 20 L⊙). The plots are from one run of POLARIS, where our simulation has been synthetically observed in dust polarization at 0.87 mm. Top left panel: the color scale and black contours represent the dust continuum total intensity (Stokes I). The white arrows represent the velocity field in a slice of the simulation centered on the sink particle. Top right panel: the color scale represents the polarized intensity (). The yellow streamlines trace the orientation of the density-weighted average magnetic field along the line of sight. The line segments represent the polarization position angle orientations. Bottom left panel: the color scale represents the fractional polarization, Ƥfrac. Bottom right panel: the color scale is the optical depth computed during the ray-tracing of the radiative transfer. In the bottom panels, the white contours trace the total intensity.

In the text

Fig. 3 Temperature maps from RAMSES and POLARIS. The left panel presents the temperature map of the RAMSES simulation cube before the POLARIS radiative transfer, from the slice centered on the sink particle (in the simulation the only central source of heating is the luminosity of the central protostellar object, with L★ = 0.58 L⊙). The right panel shows the temperature resulting from the radiative transfer performed by POLARIS for the reference case (amax = 10 µm, fhigh−J = 1, and L★ = 20 L⊙). The point is not to compare the temperature values across the two maps, but to highlight that the heating mechanisms are different - in the RAMSES simulation the dynamics of the gas induce a significant heating of a thin layer surrounding the high velocity jet.

In the text

Fig. 4 Dust polarization radiative transfer results performed at λ = 0.87 µm with fhigh−J = 1.0, L★=100 L⊙, and amax = 10 µm, after spatial filtering with the CASA simulator, mimicking ALMA interferometric observations. Top panel: synthetically observed dust continuum total intensity (Stokes I, color scale) plotted with the velocity field of the central slice in the MHD simulation (white arrows) and the density-weighted magnetic field lines (yellow lines). Stokes I is shown when I >3 σI, where σI = 0.10 mJy beam−1. Bottom panel: the color scale is the total linearly polarized intensity, which is shown where Ρ>3 σP, where σΡ = 10 µJy beam−1. The line segments represent the magnetic field orientations inferred from the dust polarization map. They are plotted where Ρ > 5 σΡ. The white circle at the bottom-left corner is the resolution of the ALMA synthetic observations and is 100 au in size.

In the text

Fig. 5 Effects of the central luminosity, L★, on the radiative transfer results at 0.87 mm, with the fixed parameters of amax = 10 µm and fhigh−J = 1. This corresponds to set III of Table 1. Each column is one radiative transfer run of POLARIS, with L★ = 1, 5, 20, 50, and 100 L⊙. Each row is a quantity provided by the radiative transfer, from the first to the sixth row: total intensity Stokes I, polarized intensity (P), polarization fraction (Ƥfrac), 2D temperature slice obtained at the center, 2D radiation field slice obtained at the center (urad/uISRF), and the 2D slice of the aalign parameter obtained at the center.

In the text

Fig. 6 Comparison of the evolution of 𝒮 × Ƥfrac as a function of normalized column density, , between ALMA observations and our models, before and after spatial filtering. The left (right) panel corresponds to the results of the radiative transfer performed at 0.87 mm (1.3 mm). The distributions of 𝒮 × Ƥfrac values from our models are in different tints of blue, while the distributions of 𝒮 × Ƥfrac values from the ALMA observations presented in Le Gouellec et al. (2020) are in gray (where we split the 0.87 and 1.3 mm observations). The solid lines are the mean of the 𝒮 × Ƥfrac values in a given bin of normalized column density, and the corresponding shaded area represents the standard deviation. These two plots correspond to set III of the models (see Table 1), i.e., we fix fhigh−J = 1 and amax = 10 µm, and we vary L★ in the range 1.0, 5.0, 20, 50, 100 L⊙. The dashed colored lines correspond to the filtered models, and the solid lines are from the non-filtered models.

In the text

	Fig. 7 Comparison of the evolution of 𝒮 × Ƥfrac as a function of normalized column density, , between ALMA observations and our models, before and after spatial filtering. Same as Fig. 6 for sets IV and V (see Table 1), where we use smaller values for amax (2 µm, top row) and for fhigh−J (0.25, bottom row).
In the text

Fig. 8 Effect of the central luminosity, L★, on the RATD of dust grains, at 0.87 mm, with the fixed parameters of amax = 50 µm and fhigh−J = 1. Each plot is a 2D slice taken at the center of the core in our radiative transfer results. Each column is one radiative transfer run of POLARIS, with L★ = 1, 5, 20, 50, and 100 L⊙. The first row is the aalign parameter, i.e., the dust grain size above which dust grains are considered aligned in our radiative transfer calculations. The second and third rows are the adisr and adisr,max parameters, respectively. In RATD theory, these two values correspond to the window in grain size, inside which dust grains are rotationally disrupted. The last row indicates whether dust grains are rotationally disrupted in our models for Smax = 105 erg cm−3, i.e., if the intervals [aalign; amax] and [adisr; adisr,max] overlap. For a given pixel, if there is an overlap, the pixel is black (i.e., grains are disrupted). If it is not the case, it is white (i.e., grains are not disrupted).

In the text

Fig. 9 Effects of the central luminosity, L★, on the RATD of dust grains at 0.87 mm, with the fixed parameters of amax = 50 µm and fhigh−J = 1. Each plot is a 2D slice taken at the center of the core in our radiative transfer results. Each column is one radiative transfer run of POLARIS, with L★ = 1, 5, 20, 50, and 100 L⊙. Each row indicates whether dust grains are rotationally disrupted in our models, i.e., if the intervals [aalign; amax] and [adisr; adisr,max] overlap. For a given pixel, if there is an overlap, the pixel is black. If it is not the case, it is white. The results from each row are obtained for a given value of grain tensile strength, Smax, i.e., 105, 106, 107 (the same plots as the bottom row of Fig. 8), 108, and 109 erg cm−3.

In the text

Fig. 10 Effects of the central luminosity, L★, on the alignment timescales in our radiative transfer models at 0.87 mm, with the fixed parameters of amax = 50 µm and fhigh−J = 1. Each column is associated with a central luminosity value, L★. While the plots in the three first rows represent the alignment timescales as a function of grain size, the plots in the last row represent the polarization fraction. The timescales are computed in specific locations, i.e., toward the equatorial mid-plane (first row), at an outflow cavity wall (second row), and in the middle of an outflow cavity (third row), in the central 2D slice of the model. These locations are indicated by the three colored circles in the polarization fraction maps, with the dark point for the first row (equatorial plane), the blue point for the second row (outflow cavity wall), and the green point for the third row (outflow cavity). In each plot of the first three rows, we present the Larmor precession timescale, τLarmor (with a dotted line for a grain without super-paramagnetic inclusions and with dot-dash and solid blue lines, respectively, for two values of Ncl, the number of atoms per cluster, 103 and 105), the radiative precession timescale, τrad, the collisional gaseous damping time of the grain, τgas (solid brown line), and the rotational disruption timescale, τdisr (for an Smax value of 105 erg cm−3, in solid dark yellow lines). We show the minimum grain size of aligned grains, aalign, with the vertical dotted green line. The range of grain size affected by rotationally disruption (adisr − adisr,max) for Smax = 105 erg cm−3 is shown by the horizontal dark line. At a given grain size, the shortest precession timescale between τLarmor and τrad dictates the direction of alignment for the grains aligned at low J, i.e., with the magnetic field or the radiation field.

In the text

Fig. B.1 Effects of the maximum dust grain size, amax, on the radiative transfer results at 0.87 mm, with the fixed parameters of fhigh−J = 1 and L★ = 20 L⊙. This corresponds to set I of Table 1. Each column is one radiative transfer run of POLARIS, with amax = 0.5, 2, 10, 30, and 50 µm. Each row is a quantity provided by the radiative transfer (from the first to the sixth row): total intensity Stokes I, polarized intensity (P), polarization fraction, 2D temperature slice obtained at the center, 2D radiation field slice obtained at the center (urad/uISRF), and the 2D slice of the aalign parameter obtained at the center.

In the text

Fig. B.2 Effects of the maximum dust grain size, amax, on the radiative transfer results at 0.87 mm, 1.3mm, and 3mm, with the fixed parameters of fhigh−J = 1 and L★ = 20 L⊙. Each column is one radiative transfer run of POLARIS, with amax = 0.5, 2, 10, 30, and 50 µm, while each row corresponds to a wavelength of observation for the radiative transfer, i.e., 0.87 mm, 1.3mm, and 3mm, as indicated on the first plot of each row. We only show the polarization fraction.

In the text

	Fig. B.3 Comparisons of the evolution of 𝒮 × Ƥfrac as a function of normalized column density, , between ALMA observations and our models, before and after spatial filtering. The left (right) column corresponds to results of the radiative transfer performed at 0.87 mm (1.3mm). This is the same as Fig. 6 but for sets I and IV (see Table 1).
In the text

	Fig. C.1 Effects of the fhigh−J parameter on the radiative transfer results at 0.87 mm, with the fixed parameters of amax = 10 µm and L★= 20 L⊙. This corresponds to set II of Table 1. Each column is one radiative transfer run of POLARIS, with fhigh−J = 0.0, 0.25, 0.5, 0.75, and 1.0. The top row is the polarized intensity (P), and the bottom row is the polarization fraction.
In the text

	Fig. C.2 Comparisons of the evolution of 𝒮 × Ƥfrac as a function of normalized column density, , between ALMA observations and our models, before and after spatial filtering. This is the same as Fig. 6 but for set II (see Table 1).
In the text

Fig. D.1 Correction for the optical thickness on the gas column density of our radiative transfer models, computed with amax = 10 µm at λ = 0.87 mm and λ = 1.3 mm. The dot-dash light blue line (solid dark blue line) shows the azimuthally averaged modified (original) gas column density as a function of radius. The solid dark red line shows the azimuthally averaged optical depth, τ, as a function of radius.

In the text