A Models \mathcal{L} and \mathcal{H}

Cosmic-ray ionisation in collapsing clouds

Key Words.:
ISM: cosmic rays – ISM: clouds, magnetic fields


Context:Cosmic rays play an important role in dense molecular cores, affecting their thermal and dynamical evolution and initiating the chemistry. Several studies have shown that the formation of protostellar discs in collapsing clouds is severely hampered by the braking torque exerted by the entrained magnetic field on the infalling gas, as long as the field remains frozen to the gas.

Aims:In this paper we examine the possibility that the concentration and twisting of the field lines in the inner region of collapse can produce a significant reduction of the ionisation fraction.

Methods:To check whether the cosmic-ray ionisation rate can fall below the critical value required to maintain good coupling, we first study the propagation of cosmic rays in a model of a static magnetised cloud varying the relative strength of the toroidal/poloidal components and the mass-to-flux ratio. We then follow the path of cosmic rays using realistic magnetic field configurations generated by numerical simulations of a rotating collapsing core with different initial conditions.

Results:We find that an increment of the toroidal component of the magnetic field, or, in general, a more twisted configuration of the field lines, results in a decrease in the cosmic-ray flux. This is mainly due to the magnetic mirroring effect that is stronger where larger variations in the field direction are present. In particular, we find a decrease of the cosmic-ray ionisation rate below  s in the central  AU, where density is higher than about  cm. This very low value of the ionisation rate is attained in the cases of intermediate and low magnetisation (mass-to-flux ratio and 17, respectively) and for toroidal fields larger than about 40% of the total field.

Conclusions:Magnetic field effects can significantly reduce the ionisation fraction in collapsing clouds. We provide a handy fitting formula to compute approximately the attenuation of the cosmic-ray ionisation rate in a molecular cloud as a function of the density and the magnetic configuration.

1 Introduction

The study of the interaction of cosmic rays (CRs) with the interstellar matter is a multi-disciplinary task that involves the analysis of several physical and chemical processes: ionisation of atomic and molecular hydrogen, energy loss by elastic and inelastic collisions, energy deposition by primary and secondary electrons, -ray production by pion decay, generation of small-scale turbulence by streaming instabilities, and the production of light elements by spallation reactions. CR ionisation activates the rich chemistry of dense molecular clouds and determines the degree of coupling of the gas with the local magnetic field, which in turn controls the ambipolar diffusion timescale and the star-formation efficiency of a molecular cloud.

In recent years a wealth of observations from the ground and from space has provided information and constraints on the flux and the ionisation rate of cosmic rays. Detections of large abundances of H in diffuse clouds (e.g. Indriolo et al. (2012)), observations of OH and HO in low H fraction regions (Neufeld et al. (2010), Gerin et al. (2010)), estimates of enhanced values of the CR ionisation rate in a molecular cloud close to a supernova remnant (Ceccarelli et al. (2011)) as well as the measurement of the luminosity of molecular clouds (e.g. Montmerle (2010)) raised the questions about the origin of the CR flux that generates such a high ionisation rate and how to reconcile these values with those ones measured in denser clouds that are more than one order of magnitude lower. A number of studies approached this problem using different strategies, analysing the effects of Alfvén waves on CR streaming (Skilling & Strong (1976), Hartquist et al. (1978), Padoan & Scalo (2005), Rimmer et al. (2012)), magnetic mirroring and focusing (Cesarsky & Völk (1978), Chandran (2000), Padovani & Galli (2011)), or the possible existence of a low-energy flux of CR particles able to ionise diffuse but not dense clouds (Takayanagi (1973), Umebayashi & Nakano (1981), McCall et al. (2003), Padovani et al. (2009)).

Disc formation is another integral aspect of star formation. One of the main concerns is the so-called “catastrophic magnetic braking problem” that suppresses the formation of a rotationally supported disc in the ideal MHD limit during the protostellar accretion phase of a low-mass forming star (Allen et al. (2003), Galli et al. (2006), Mellon & Li (2008), Hennebelle & Fromang (2008)). Given the observational evidence of discs on 100 AU or even larger scales at least around Class I–II protostars, a number of possible solutions to the problem of catastrophic magnetic braking have been invoked, including: (i) non-ideal MHD effects such as Ohmic dissipation (Shu et al. (2006), Dapp & Basu (2010)) and Hall effect (Krasnopolsky et al. (2011), Braiding & Wardle (2012a); (2012b)); (ii) the possible misalignment between the rotation axis and the magnetic field direction that acts reducing the braking torque (Hennebelle & Ciardi (2009)); (iii) depletion of the infalling envelope that anchors the magnetic field braking (Mellon & Li (2009), Machida et al. (2011)); (iv) turbulent diffusion of the magnetic field (Seifried et al. (2012), Santos-Lima et al. (2013)).

Mellon & Li (2009) advanced the interesting possibility that a reduction of the CR ionisation rate, , corresponding to a decrease of the ionisation fraction by a factor , could result in sufficient ambipolar diffusion to allow the formation of a rotationally supported disc. They concluded that both a suppression of the CR flux and a low level of magnetisation (measured by the non-dimensional mass-to-flux ratio ) were needed in order to circumvent the magnetic braking problem. Although they did not perform a detailed exploration of the parameter space of their models, Mellon & Li (2009) found that a value of  s was needed to spin up the gas significantly during collapse if (but no disc larger than  AU was formed in this case), or to form a fully rotationally supported disc of radius  AU if .

Following the suggestion by Mellon & Li (2009), in this paper we focus on the influence of different magnetic field and density configurations on the CR propagation following the conclusions achieved in Padovani et al. ((2009)(2011)). Our aim is to show that in the inner regions of a cloud, where the disc is formed, magnetic and column-density effects can indeed cause a significant decrease of the interstellar CR ionisation rate and consequently of the ionisation degree, helping to decouple the gas from the magnetic field.

The paper is organised as follows. In Section 2 we provide a detailed description of the method used to calculate the CR ionisation rate. In Section 3 we analyse a semi-analytical model of a singular isothermal toroid threaded by a toroidal magnetic field with the purpose of understanding the role of column density versus magnetic effects. In Section 4 we explore the evolution of the CR ionisation rate on the initial conditions (mass-to-flux ratio and alignment between rotation axis and magnetic field direction) for a number of numerical simulations. In Section 5 we discuss the variations of the CR ionisation rate in discs and in Section 6 we give a fitting formula to compute the CR ionisation rate accounting for the magnetic field configuration. In Section 7 we summarise our conclusions. Comments on other possible models of the CR ionisation rate are provided in Appendix A.

2 Method

Cosmic rays, being charged particles, perform an helicoidal motion around the magnetic field lines and we follow their path starting from the outer boundary of the computational domain throughout the core. As it is well known, a charged particle of mass and velocity spiraling along a magnetic field of increasing strength must increase its pitch angle (the angle between the particle’s velocity and the field direction) as a consequence of the conservation of kinetic energy and magnetic moment . In particular, for a particle starting from the intercloud medium (ICM) with a pitch angle and a magnetic field strength , the pitch angle is given by


where is the focusing factor (see e.g. Desch et al. (2004) and Padovani & Galli (2011), hereafter PG11). The assumption of energy conservation along the particle’s trajectory is clearly violated in the presence of collisional losses. In principle, Eq. 1 should be replaced by an equation for the time evolution of the pitch angle , including the effect of the magnetic field as well as the diffusion induced by collisional ionisation of H molecules. In the present study we neglect these aspects, and while we assume that kinetic energy is conserved along each individual trajectory, we take into account energy losses globally by propagating the CR spectrum inside the cloud. We have verified that this approximation is valid for pitch angles not too small (that evolve to before substantial energy losses occur) and for proton energies larger than about 10 MeV. For CR protons of lower energies, our treatment overestimates somewhat the efficiency of magnetic mirroring. Since the bulk of ionisation from CR protons at the typical column densities of molecular clouds is produced by particles in the energy range 1-300 MeV (Padovani et al. (2009), hereafter PGG09), our approximation is satisfactory. For CR electrons our approximation is satisfied for all energies of interest. We also assume that the number of particles is conserved, ignoring electron capture reactions of CR protons with H and He as well as the fusion reactions that form Li and Li, because of the small cross sections (Meneguzzi et al. (1971)).

The column density of H passed through by the particle is given by


where is the maximum depth reached inside the core and is the H volume density. If the pitch angle remains smaller than along the entire particle’s trajectory, CRs of sufficient energy can cross the whole core. Vice versa, if reaches at some point inside the cloud, the particle is mirrored and it will follow the same field line backwards.

Once the column density is known for each value of the initial pitch angle and for each field line, we compute the CR ionisation rate using the CR propagation theory developed in PGG09. Since two of their models, those where the ionisation degree is dominated by CR electrons, are very similar, we assume three possible trends for the ionisation rate as a function of the column density, indicated by with , , and , where =low, =medium, and =high correspond to the models W98+C00, M02+C00, and the average between the two models W98+E00 and M02+E00 detailed in PGG09, respectively. PGG09 shows that the decrease of the CR ionisation rate with increasing penetration in the cloud can be described by a power-law at column densities in the range  cm,


and by an exponential attenuation for  cm,


Padovani et al. ((2013), hereafter PGG13) give a simple fitting formula which combines in a single expression the low- and high-column density approximations above,


where is the effective surface density seen by a CR propagating with pitch angle , the proton mass and the molecular weight for the assumed fractional abundances of H and He. We have used this fitting formula for the three models adopted in this work (see Fig. 1). The fitting coefficients of Eq. (5) are given in PGG09 and PGG13.

Figure 1: CR ionisation rate as a function of the molecular hydrogen column density for the three models described in the text.

In order to evaluate the ionisation rate along a field line, we average the contribution of all the particles with different initial pitch angles over the solid angle


Finally, to account for magnetic focusing, we multiply by . In the following, we choose as a “fiducial” spectrum the one corresponding to the case and we neglect the subscript . Appendix A shows the results for the other two models .

3 Semi-analytic model

In PG11, we modelled molecular cloud cores as singular isothermal toroids, namely scale-free, axisymmetric equilibrium configurations of an isothermal gas cloud under the influence of self-gravity, gas pressure, and magnetic forces (Li & Shu (1996)) These models are uniquely characterised by the mass-to-flux ratio defined by


where is the gravitational constant, the magnetic flux, and the mass contained in the flux tube . In this model the purely poloidal magnetic field threading the core takes the form


where is the magnetic flux function. Separation of variables requires the equilibrium solution to take the self-similar form


where is the sound speed and is a dimensionless function.

If the core rotates with respect to the ambient medium, a toroidal component of the magnetic field is expected to arise. We introduce a time-independent toroidal magnetic field component by adding a term


with constant to the magnetic field. This term is curl-free, divergence-free, and scales as the poloidal component of the magnetic field. Additionally, as is unchanged by this term, the equilibrium equations of the core are independent of the toroidal component and, therefore, the solutions for the density and flux functions obtained by Li & Shu ((1996)) remain formally valid. We are aware that, in general, one does not expect the toroidal component of the field generated by the rotation of the core to be curl-free. Nevertheless, it is very instructive to look at these idealised configurations with the aim of disentangling column-density from magnetic-field effects in the propagation of CR particles.

Substitution of Eq. (9) into Eq. (8), and the addition of the toroidal component (Eq. 10) results in




is the ratio between the strength of the toroidal, , and the poloidal, , field in the cloud’s midplane. Notice that the magnetic field diverges on the polar axis. However, this is not a problem for our calculations, as the polar axis contains no mass.

The field lines that cross the midplane at and , are then be given by




As shown by Fig. 2, the magnetic field lines lie on the surfaces of nested flux tubes defined by the condition , and their twisting increases as the distance from the axis of symmetry decreases.

Figure 2: Magnetic field lines of the toroid with (purely poloidal, upper plot), (middle plot), and (lower plot). Iso-density contours are shown in grey scales in unit of cm in logarithmic scale, while the magnetic field module is shown in colour scale in unit of G. Field line plots are generated using the visualisation software package MAYAVI (Ramachandran & Varoquaux (2011)).

3.1 Column-density versus magnetic-field effects on cosmic-ray ionisation rate

During the propagation of a CR, the ionisation rate decreases due to two main factors: CRs lose energy, thus their capability of ionising hydrogen molecules, because of the increasing column density “seen” by the particles themselves (PGG09); magnetic mirroring reduces the ionisation rate more than magnetic focusing amplifies it (PG11). In order to distinguish the origin of the variation in the ionisation rate during the propagation, we assume several values for and resulting in different density profiles and magnetic configurations.

We assume three different values for the mass-to-flux ratio () that acts modifying the density distribution as well as the pinching of the magnetic field lines. corresponds to the case of a “roundish” core, with a density distribution almost spherically symmetric; when the density is distributed in a disc-like configuration; finally, represents the intermediate case. The pinching of magnetic field lines increases for decreasing .

We start assuming (purely poloidal field) whose effects on CR propagation are described in PG11, then we increase the toroidal field from (poloidal and toroidal fields have comparable strengths) to . We follow the propagation of CRs entering the cloud from any direction with all possible initial pitch angles applying the method described in Sect. 2. With respect to PG11 we also consider the fact that the mirrored CRs can still ionise while propagating backwards. As expected, we find that this contribution to the ionisation is stronger in the outer part of the cloud, namely in the region also crossed by CRs which propagate backwards. This further ionisation is not crucial since these CRs have already passed through a large column density, losing their ionising power, but now it is included in the model. Figure 3 shows the order of magnitude of this contribution for the model : we plot the trend of the CR ionisation rate for a flux tube enclosing a mass of including all the contributions deriving from mirroring, focusing, and backward ionisation.

Figure 3: Contributions to the ionisation rate as a function of the position along the symmetry axis for the model and a flux tube enclosing 1 .

3.2 Dependence of on the magnetic field configuration

Figure 4: CR ionisation rate profiles for the model in the plane crossing the symmetry axis and perpendicular to the midplane . The mass-to-flux ratio is and the strength of the toroidal field increases from left to right.

In order to investigate on the effects of a given magnetic field configuration on the variations of , we assume a mass-to-flux ratio so that the density shape is fixed and we increase the magnitude of the toroidal component (Fig. 2). We find that increasing values of correspond to a decrease of , in fact the area with low ionisation rate becomes larger and larger reflecting the signatures of the magnetic field configuration (Fig. 4). Since , the intensity of the magnetic field increases towards the symmetry axis. Therefore, the focusing factor becomes larger closer to the cloud’s axis, and this explains the increase of in this region.

As increases, the pitch angle gets closer to quickly and the mirroring becomes more and more effective. However, as we know from PG11, increasing the toroidal component of does not change significantly the relative importance of mirroring versus focusing effects, but it increases the effective column density as seen by a CR gyrating around the field. If this increase in reaches the regime of exponential attenuation ( cm, see Fig. 1), the effect can be dramatic. However, this hardly happens in our semi-analytical model, which is taken as representative of a pc-scale clump of modest column density. Using profiles of density and magnetic field strength valid for smaller spatial scales taken from simulations of collapsing clouds, we will show that a reduction of of orders of magnitude can be achieved (see Sect. 4).

Figure 5 shows a zoom in the central region at the core-size scale of 0.1 pc radius for the model in order to appreciate more explicitly the reduction in the ionisation rate: when the toroidal field is 50 times stronger than the poloidal field component, a value of lower than the canonical value of  s ( s) is readily reached in the region close to the cloud’s midplane.

Figure 5: CR ionisation rate maps for the model in the central 0.1 pc region in the plane for and increasing . Black and white solid lines show the iso-ionisation rate contours.

3.3 Dependence of on the density profile

To examine the effects of the density configuration on the ionisation rate, we fix a value of the toroidal-to-poloidal ratio and we vary the mass-to-flux-ratio obtaining density profiles that span from a “roundish” core to a disc-like configuration (see Sect. 3.1). The effects of the column density can be promptly recognised by noticing that the ionisation rate profile follows the shape of the density profile. In Fig. 6 we superpose the iso-density contours to the ionisation rate maps. The regions where iso-density contours follow the ionisation rate profile reveal column-density effects. Conversely, the areas where the iso-density contours depart from the spatial distribution of are representative of magnetic-field effects. In fact, for high values of the column density crossed is larger and decreases rapidly (see PGG09). Said in a different way, in the model a particle “sees” a larger amount of and it loses more energy than in the the flatter model where the density along the field direction increases only close to the cloud’s midplane.

Figure 6: CR ionisation rate maps for the model in the plane for a fixed toroidal-to-poloidal ratio and different values of the mass-to-flux-ratio. White and black contours represent the iso-density contours and the labels show .

4 Numerical models

In this Section we describe the ionisation rate maps obtained from a number of numerical simulations related to a collapsing rotating core performed with the AMR code RAMSES1 (Teyssier (2002), Fromang et al. (2006)) and detailed in Joos et al. ((2012)). They considered a spherical cloud with an initial density profile


where  cm and  pc according to observations (André et al. (2000), Belloche et al. (2002)). The thermal-to-gravitational energy ratio is about 0.25 and the rotational-to-gravitational energy ratio is about 0.03. We select a series of simulations from Joos et al. ((2012)) varying the mass-to-flux ratio and the angle between the initial magnetic field direction and the initial rotation axis . Table 1 lists the parameters.

Case Disc ?
[rad] [kyr] (Y / N / K)
A 5 0 0.824 N
A 5 0 11.025 0.26 0.05 N
B 5 /4 7.949 0.23 0.15 Y
C 5 /2 10.756 0.46 0.28 K
D 2 0 5.702 0.24 N
E 17 0 6.620 0.43 0.15 K

A disc with flat rotation curve is formed (Fig. 15 in Joos et al. (2012)).
No significant disc is formed ().
A keplerian disc is formed (Fig. 14 in Joos et al. (2012)).

Table 1: Parameters of the simulations described in the text (from Joos et al. (2012)): mass-to-flux ratio, initial angle between the magnetic field direction and the rotation axis, time after the formation of the first Larson’s core (core formed in the centre of the pseudo-disc with  cm and  AU), maximum mass of the protostellar core and of the disc. Last column gives information about the disc formation.

According to the method described in Sect. 2 we compute the ionisation rate maps making use of the model . For each case in Table 1, we show the maps of in a plane parallel and perpendicular to the main direction of the magnetic field (always along the axis) and containing the density peak, namely the and the plane, respectively. Besides, in Figures 712 we plot the results for the whole computational domain (upper and middle left plots) and zooming into the inner 1000 AU (upper and middle right plots) plus a graph showing the magnetic field line morphology in the central 600 AU (lower plot).

We are aware of the fact that the ionisation rate should be computed simultaneously to the MHD model. This will be the subject of a future work, but the results described in this paper can be considered an important proof of concept that very low ionisation rate can be achieved in the inner regions of a collapsing cloud.

4.1 Intermediate magnetisation ()

We consider a couple of outputs (cases A and A in Table 1) for an aligned rotator () and super-critical cloud () in agreement with observations (Crutcher (1999)). We choose two different times after the formation of the first Larson’s core in order to understand the effect of the tangling of magnetic field lines on CR propagation when a disc is not formed. Unlike the semi-analytical case (Sect. 3) where the density does not even achieve  cm and the symmetry of the magnetic field configuration is conserved for any toroidal-to-poloidal ratio, in all the numerical models presented here the central density reaches the value of about  cm and in general it is higher than  cm in the inner  AU radius. Besides, the symmetry of magnetic field lines is broken very soon with time. This is likely due to the development of the interchange instability (Krasnopolsky et al. (2012)).

We know from the semi-analytical model (Sect. 3) that it is not possible to disentangle column-density from magnetic effects, but both intervene on the decrease of the ionisation rate. In Section 6 we give an estimate of the relative incidence of these two effects. As explained in Sect. 3.3, we can interpret the deviations between the iso-density contours and maps as due to magnetic imprints. As an instance, in the upper right plot of Fig. 7 there is a clear departure between the contours at  cm and the shape of the region where reaches values of  s (in yellow in the Figure). Another example is the upper left plot of the same Figure where the region in red with  s extends for densities spanning from to more than  cm in a horizontal “strip” of about 5000 AU along the axis. This calls to mind the magnetic field configuration (see lower panel of Fig. 7), in fact this is the region where field lines start to be twisted due to rotation. The middle right panel of Fig. 7 shows that less than  yr after the formation of the first Larson’s core, a central region with  AU is characterised by  s.

Figure 7: CR ionisation maps and iso-density contours (black solid lines) for the case A in Table 1. Upper and middle left panels show the entire computational domain while upper and middle right panels show a zoom in the inner region. Upper panels show a slice of a plane parallel to the magnetic field direction, while middle panels refer to a slice of a perpendicular plane. Both planes contain the density peak and labels show . The lower plot shows the magnetic field line morphology in the central 600 AU of the RAMSES data cube. The colour bar shows the magnetic field module in Gauss units.

We compute for a later-time configuration (case A) with the same initial conditions (, ). As previously mentioned, we lose the symmetry of field lines (see lower panel of Fig. 8) as well as of the density profile (see upper and middle right panels of Fig. 8). At large scales (upper and middle left plots), we notice that the region with  s is less extended along the axis and it is elongated parallel to the axis (i.e., parallel to the magnetic field). At small scales (upper and middle right plots), a flattened structure formed almost perpendicularly to the rotation axis can be noticed. In fact, even if the disc is not formed, the plane perpendicular to the rotation axis shows the presence of a “ring” with densities between 10 and 10 cm and average ionisation rate of about 10 s circumscribing the density peak up to a radius of about 300 AU.

Figure 8: CR ionisation maps and iso-density contours (black solid lines) for the case A in Table 1. See Fig. 7 for further information.

We account for a configuration with the initial rotation axis twisted by towards the axis, while the initial magnetic field direction is still along the axis (case B in Table 1). This case predicts the formation of a disc with a flat rotation curve in the plane perpendicular to the rotation axis. In the central  AU region of both the and the plane, we find a CR ionisation rate lower than about  s. In the inner 500 AU region (upper and middle right panels of Fig. 9) we can still identify a relationship between and , even if the density profile becomes more and more irregular and field lines almost lose memory of their initial configuration (see lower panel of Fig. 9).

Figure 9: CR ionisation maps and iso-density contours (black solid lines) for the case B in Table 1. See Fig. 7 for further information.

Finally, we consider the case of a perpendicular rotator () so that initially the rotation axis is along the axis. In this circumstance (case C in Table 1), a keplerian disc perpendicular to the rotation axis is predicted. Figure 10 shows this configuration, namely a face-on view in the plane and an edge-on view of the disc in the plane. It is worth noting in the face-on view (upper right panel of Fig. 10) a large region of  AU and  cm with  s. Here the ratio between the toroidal component and the total magnetic field, , is larger than about 0.4. We reach even lower values, with a minimum of  s in the inner area that has an extent of a few tens of AU. The CR ionisation rate is so low that we can assume the gas to be effectively decoupled with the magnetic field. A similar behaviour can be also appreciated in the edge-on view (lower right panel of Fig. 10). The very low value of found in this collapse region corresponds to the condition required by Mellon & Li ((2009)) for the formation of a 10 AU disc when . However, we stress that in this case the formation of the disc is made possible by the misalignment of the angular momentum and the magnetic field of the cloud, whereas in the situation considered by Mellon & Li ((2009)) the two vectors are aligned and disc formation is enabled by the enhanced ambipolar diffusion resulting from the lower value of . Finally, the lower panel of Fig. 10 shows how the rotation perpendicular to the axis forces the field lines (initially with a poloidal configuration along the axis) to wrap around the rotation axis.

Figure 10: CR ionisation maps and iso-density contours (black solid lines) for the case C in Table 1. See Fig. 7 for further information.

4.2 Strong magnetisation ()

We analyse the case of an aligned rotator with strong magnetic field for a late-time configuration (case D of Table 1). A stronger field is more resistant to line twisting caused by rotation. In fact, in this case, the poloidal configuration can be still identified after about  kyr from the formation of the first Larson’s core (see lower panel of Fig. 11). As expected, the effect of magnetic braking is remarkable and the disc is not formed in the plane perpendicular to the rotation axis (see middle right panel of Fig. 11). In this case we find the ionisation rate to be more uniform and close to  s, while in the inner  AU is of the order of  s in the plane and  s in the plane. Density contours in the plane reveal high density fluctuations not allowing the formation of large high-density structures.

Since most molecular cloud cores appear to have significant levels of magnetisation (Crutcher 1999), this case is probably the most relevant for modelling cloud collapse. Our findings that the CR ionisation rate is one-two orders of magnitudes below the interstellar value in a region around the accreting protostar of radius  AU, could indicate that a centrifugally supported disc of this size might form even in the strongly magnetised case. Of course, this question can only be answered by performing a fully self-consistent dynamical calculation of CR propagation during cloud collapse.

Figure 11: CR ionisation maps and iso-density contours (black solid lines) for the case D in Table 1. See Fig. 7 for further information.

4.3 Weak magnetisation ()

We model the case of weak magnetisation for an aligned rotator and we consider a late-time configuration (case E in Table 1). As shown in the lower panel of Fig.12, since the magnetic field braking is very weak, rotation is able to strongly twist the field lines. It is interesting to notice how the contribution to the decrease of is substantial perpendicularly to the disc plane. In fact, the upper right panel of Fig. 12 shows a large region along the axis where  s. This region is not limited to the high-density domain ( cm), but it broadens out along the rotation axis where the magnetic mirroring due to field line tangling up is very marked. As for case C (Sect. 4.1), the region with  cm and  s is characterised by . In the plane one can see the presence of the face-on disc and the rapid decrease of CR ionisation rate that reaches about  s in the inner  cm iso-density contour. This is compatible with the results of Mellon & Li ((2009)) who find that a centrifugally supported disc of radius  AU is formed by the collapse of a cloud with , close to our value of 17, if  s, which is exactly the value that we find in this region. Again, the agreement found in this case may not be particularly significant, because in clouds characterised by an initially weak field, the magnetic braking may result inefficient regardless of the degree of ionisation of the gas. We stress again that, in order to draw firm conclusions on the role of CR ionisation in the resolution of the so-called “magnetic braking problem”, a fully self-consistent calculation including CR propagation and magnetic diffusion is required. In Sect. 6 we outline a procedure to include approximate treatment of CR transport in a MHD simulation.

Figure 12: CR ionisation maps and iso-density contours (black solid lines) for the case E in Table 1. See Fig. 7 for further information.

5 Cosmic-ray ionisation rate at high column densities

The regime of high column densities,  cm, corresponding to surface mass densities larger than a few g cm, is interesting for applications to protoplanetary discs, where the ionisation rate plays a major role in determining the “dead zones” where the generation of turbulence by the magnetorotational instability is relatively inefficient (see e.g. Sano et al. (2000), Okuzumi & Hirose (2011)). In these applications, the CR ionisation rate is usually assumed to be exponentially attenuated within the disc, with a decay length  g cm (Umebayashi & Nakano (1981))2. Umebayashi & Nakano ((2009)) propose an empirical formula for the CR ionisation rate as a function of the depth from the disc surface, assuming a geometrically thin disc and taking into account the fact that CRs penetrate the disc almost isotropically. Their formula in cylindrical coordinates reads


where  s and and are the vertical gas surface densities measured from the upper () and the lower boundary () of the computational domain, respectively, namely




where is the proton mass and is the molecular weight. We consider the numerical models described in Sect. 4 that allow the formation of a keplerian disc (case C and E in Table 1) comparing our CR ionisation rate maps with those obtained by using Eq. (16).

In order to make a consistent comparison, we let CRs propagate along straight lines without including magnetic effects, computing using Eq. (5). This equation gives similar results to Eq. (16), the small variations being due to the difference in the assumed attenuation length (see left panels of Fig. 13). Nevertheless, the picture changes dramatically when considering the true path of CRs along the field lines with the inclusion of magnetic and focusing effects. In fact, the region with  s grows considerably (see right panels of Fig. 13).

Figure 13: Comparison between the CR ionisation rate distribution obtained by using the fitting formula (Eq. 16) from Umebayashi & Nakano ((2009)) considering rectilinear propagation, black contours, and Eq. (5), white contours, for rectilinear propagation (left panels) and following the path of CRs along field lines with the inclusion of mirroring and focusing effects (right panels). The three concentric black contours in all the panels as well as the white contours in the left panels refer to going inwards. Labels are on a logarithmic scale and the colour coding shows the density distribution. Upper panels: model C; lower panels: model E of Table 1.

This comparison represents the ultimate proof of the role of magnetic fields on CR propagation. Magnetic effects cannot be neglected since they set the extent to which CRs can determine the coupling between magnetic fields and the gas. In the next Section, we give a new fitting formula for that accounts for the magnetic configuration.

6 Magnetic effects on the reduction of the CR ionisation rate: a useful fitting formula

In order to compare the relative importance of density and magnetic effects on the decrease of in the central region of a core ( cm), we let CRs propagate along field lines without accounting for mirroring and focusing, we refer to this setting as the “non-magnetic case”. This is easily done by following the path of CRs with initial pitch angle . According to Eq. (1), the pitch angle during the propagation of the particle remains equal to zero. We generate CR ionisation maps as those presented in Figures 7 – 12 and then we compute the quantity defined as the ratio between the rates obtained including and neglecting magnetic effects as in PG11.

We focus on cases C and E (see Table 1), corresponding to the cases where a keplerian disc is formed. In these models, the density reaches values larger than  cm inside a region of a few 100 AU in size. Thus, one can be led to conclude that magnetic effects may be negligible since CRs have crossed a large column density so that the regime of exponential attenuation ( cm) is attained. On the contrary, we find that even at very high densities the role of magnetic fields in removing CRs and attenuating is substantial. The results are shown in Fig. 14, where we plot the ratio of the CR ionisation rate computed keeping into account the evolution of the pitch angle to the same quantity computed in the non-magnetic case. We see that magnetic effects cause a net reduction of of a factor of at large scales up to a factor of 3 to 10 in the inner regions.

We give a convenient fitting formula in order to reproduce the ionisation rate maps without running the full code described in the paper. When in presence of a magnetic field, the effective column density, , seen by a CR can be much larger than that obtained through a rectilinear propagation (see also Sect. 5): the more complex the morphology of the magnetic field, the larger is the path covered by a charged particle. If is the average column density seen by an isotropic flux of CRs, has the form


The factor depends on the ratio between the toroidal and the poloidal components of the magnetic field, , as well as on its module. It reads




and being the minimum and the maximum value of in the whole data cube, respectively. When the magnetic field strength is negligible, CRs propagate along straight lines. In this case and , otherwise and . Comparing the column density seen by a CR evaluated with our code and with Eq. (19), we find that it is safe to approximate the poloidal configuration to a rectilinear path. This explains why we introduce the dependence of on . In fact, varies between 0 and 1, being equal to 0 for a purely poloidal configuration. We also find that the higher the density, the stronger is the role of the magnetic field in increasing . This justifies the presence in Eq. (19) of the power that reads


and being the minimum and the maximum value of the density in the whole data cube, respectively. The factor 0.7 in Eq. (22) has been introduced to reproduce the CR ionisation maps in the upper and middle panels of Figures 712. Notice also that the factor depends on the local value of the magnetic field, since and are computed in a given point, but it also accounts for the large scale configuration by means of and .

Figure 14: Maps of the ratio between the CR ionisation rates in the magnetic and non-magnetic case for the case E in Table 1. Left panels show the entire computational domain while right panels show a zoom in the inner region. Upper panels show a slice of a plane parallel to the magnetic field direction, while lower panels refer to a slice of a perpendicular plane. Iso-contours of are shown in cyan (), yellow (), orange (), and purple ().

Once evaluated the effective column density, the corresponding (effective) CR ionisation rate, , is obtained by


where is computed using Eq. (5) after replacing and with and , respectively. The factor is given by


and it represents the correction for magnetic effects. Equation (23) gives a correct result within a factor of 3 and it holds for magnetic field strengths smaller than about 1 G. It can be very helpful for non-ideal MHD simulations so as to compute diffusion coefficients and it can be also implemented in chemical models in order to have a more precise description of the observational results. Figure 15 shows the goodness of the fit for the inner region of case C in Table 1.

Figure 15: Comparison between the CR ionisation rate evaluated with the code described in the paper, solid contours, and with Eq. (23), dashed contours for the case C in Table 1. The grey shaded area shows the region where the difference among contours is larger than a factor of 3 and the black solid contour refers to the region where  G. Green, orange, blue, cyan, and red contours are related to and , respectively.

7 Conclusions

In this study we explored the distribution of the CR ionisation rate in molecular clouds, examining and extending the results obtained in PGG09 and PG11. We employed a semi-analytical model in order to understand when the variations of can be attributed to energy losses due to the increasing column density passed through or to magnetic mirroring and focusing. The main conclusion is that an increment of the toroidal component, and in general a more tangled magnetic field, corresponds to a decrease of because of the growing preponderance of mirroring over focusing. That is to say, large variations of the direction of field lines cause a rapid increase of the cosmic-ray pitch angle towards the mirroring angle. Conversely, fixing a magnetic field line configuration while varying the density profile allowed us to identify the weakening of the ionisation power caused by energy losses. As expected, moving from a “roundish” core towards a disc-like configuration, the distribution of follows the density profile shape.

In the second part of the paper we analysed a number of numerical simulation outputs related to a rotating collapsing core following the propagation of cosmic rays at different time steps, varying the degree of magnetisation and the initial orientation of the main magnetic field direction with respect to the rotation axis. Being aware of the fact that the correct manner of dealing with CR propagation should be computing their distribution simultaneously with the MHD simulation, we believe that our conclusions represent an important proof of concept. In particular, in the central 100 AU region, the number density is higher than  cm and a H column density larger than 10 cm is promptly reached. This actually means entering the exponential attenuation regime where the CR ionisation rate is independent of the CR spectrum assumed and drops below  s. However, we prove that even in the high-density region, the presence of a magnetic field can reduce up to a factor larger than 10. As for the semi-analytical model, we also conclude that the morphology of the maps depends both on the density profile and on the magnetic field line configuration.

Finally we focused on the morphology of in the inner region where a keplerian disc is formed. We found that the inclusion of magnetic field effects in the calculation of brings to the formation of a large central region of 100 to 200 AU where the CR ionisation rate is well below the ordinarily used value of  s. This provides support to the hypothesis of Mellon & Li ((2009)) according to which the magnetic braking efficiency can be reduced if  s allowing the formation of the disc. In order to test this hypothesis a self-consistent MHD collapse calculation including CR propagation is needed.

We compared our results with those obtained by Umebayashi & Nakano ((2009)) in the case of an unmagnetised disc, showing that the exclusion of CRs resulting from magnetic mirroring deeply affects the CR ionisation rate pattern in the collapse region. To account for the effects of the magnetic configuration, we formulated a general fitting expression to approximately compute as a function of the column density, magnetic field strength, and toroidal-to-poloidal magnetic field ratio. This empirical expression reproduces quite accurately our numerical results.

Non-ideal MHD models predict that magnetic braking becomes inefficient at densities  cm, when magnetic field diffusion becomes faster than the dynamical evolution (Dapp et al. (2012)). In our models we observe that the drop in takes place in some cases even at lower densities ( cm), resulting in very low ionisation fractions. The consequences of the reduced CR ionisation rate on the magnetic diffusion coefficients (ambipolar, Hall, and Ohm) will be the subject in a forthcoming paper.

MP thanks Marc Joos and Andrea Ciardi for their help in accessing RAMSES simulations and Natalia Dzyurkevich for illuminating discussions about ionisation in discs. MP and PH acknowledge the financial support of the Agence National pour la Recherche (ANR) through the COSMIS project.

Appendix A Models and

For the sake of completeness, we show the CR ionisation rate obtained adopting the spectra corresponding to the cases and in Fig. 1. These represent two extreme behaviours of CR ionisation as a function of column density. For both cases, we compute the dependence of on the density profile and the magnetic field configuration for the semi-analytical cloud models described in Sect. 3. In particular we show results for , and 1.63 with (Fig. 16) and for with , and 50 (Fig. 17). As expected, the CR ionisation rate for the case is independent of the value of the mass-to-flux ratio and the toroidal-to-poloidal ratio since the exponential regime ( cm) is not attained and is substantially independent of the column density.

On the other hand, the CR ionisation rate for the model is comparable with that of the model described in the main text. However, as shown in Fig. 17, even when the toroidal field is strong (), in the inner region of the toroid below 0.1 pc, is of the order of  s, about one order of magnitude higher than the values estimated from observations. For this reason, the model has been chosen as the fiducial model.

Figure 16: CR ionisation rate maps for the model (left column) and (right column) in the plane for a fixed toroidal-to-poloidal ratio and different values of the mass-to-flux-ratio. White and black contours represent the iso-density contours and the labels show .
Figure 17: CR ionisation rate profiles for the model (top row) and (bottom row) in the plane crossing the symmetry axis and perpendicular to the midplane . The mass-to-flux ratio is and the strength of the toroidal field increases from left to right.


  1. RAMSES simulations were analysed using PyMSES (Labadens et al. (2011)).
  2. The calculations of PGG09 and PGG13 suggest however that the attenuation length for protons may be a factor of larger than that derived by Umebayashi & Nakano ((1981)).


  1. Allen, F. C., Li, Z.-Y. & Shu, F. H. 2003, ApJ, 599, 363
  2. André, P., Ward-Thompson, D. & Barsony, M. 2000, Protostars and Planets IV, 59
  3. Belloche, A., André, P., Despois, D. et al. 2002, A&A, 393, 927
  4. Braiding, C. R. & Wardle, M. 2012a, MNRAS, 422, 261
  5. Braiding, C. R. & Wardle, M. 2012b, MNRAS, 427, 3188
  6. Ceccarelli, C., Hily-Blant, P., Montmerle, T. et al. 2011, ApJ, 740, L4
  7. Cesarsky, C. J. & Völk, H. J. 1978, A&A, 70, 367
  8. Chandran, B. D. G. 2000, ApJ, 529, 513
  9. Crapsi, A., Caselli, P., Walmsley, C. M. et al. 2007, A&A, 470, 221
  10. Crutcher, R. M. 1999, ApJ, 520, 706
  11. Dapp, W. B. & Basu, S. 2010, A&A, 521 256
  12. Dapp, W. B., Basu, S. & Kunz, M. W. 2012, A&A, 541, A35
  13. Desch, S. J., Connolly, H. C., Jr., & Srinivasan, G. 2004, ApJ, 602, 528
  14. Draine, B. T. & Sutin, B. 1987, ApJ, 320, 803
  15. Fromang, S., Hennebelle, P. & Teyssier, R. 2006, A&A, 457, 371
  16. Galli, D., Lizano, S., Shu, F. H. & Allen, A. 2006, ApJ, 647, 374
  17. Gerin, M., De Luca, M., Black, J. et al. 2010, A&A, 518, L110
  18. Hartquist, T. W., Doyle, H. T. & Dalgarno, A. 1978, A&A, 68, 65
  19. Hennebelle, P. & Fromang, S. 2008, A&A, 477, 9
  20. Hennebelle, P. & Ciardi, A. 2009, A&A, 506, L29
  21. Indriolo, N. & McCall, B. J. 2012, ApJ, 745, 91
  22. Joos, M., Hennebelle, P. & Ciardi, A. 2012 A&A, 543, 128
  23. Krasnopolsky, R., Li, Z.-Y. & Shang, H. 2011 ApJ, 733, 54
  24. Krasnopolsky, R., Li, Z.-Y. & Shang, H. et al. 2012 ApJ, 757, 77
  25. Labadens, M., Chapon, D., Pomaréde, D. et al. 2011, ADASS XXI Proceedings, 461, 837
  26. Li, Z.-Y. & Shu, F. H. 1996, ApJ, 472, 211
  27. Machida, M. N., Inutsuka, S.-I. & Matsumoto, T. 2011, PASJ, 63, 555
  28. McCall, B. J., Geballe, T. R., Hinke, K. H., et al. 1998, Science, 279, 1910
  29. McCall, B. J., Huneycutt, A. J., Saykally, R. J. et al. 2003, Nature, 422, 500
  30. Mellon, R. R. & Li, Z.-H. 2008, ApJ, 681, 1356
  31. Mellon, R. R. & Li, Z.-H. 2009, ApJ, 698, 922
  32. Meneguzzi, M., Adouze, J. & Reeves, H. 1971, A&A, 15, 377
  33. Montmerle, T. 2010, in ASP Conf. Ser. 422, High Energy Phenomena in Massive Stars, ed. J. Martí, P. L. Luque-Escamilla, & J. A. Combi (San Francisco, CA: ASP), 85
  34. Nakano, T. 1984, Fund. Cosmic Phys., 9, 139
  35. Neufeld, D. A. and 52 co-authors 2010, A&A, 521, 10
  36. Okuzumi, S. & Hirose, S., ApJ, 742, 65
  37. Padoan, P. & Scalo, 2005, ApJ, 624, L97
  38. Padovani, M., Galli, D. & Glassgold, A. E. 2009, A&A, 501, 619 (PGG09)
  39. Padovani, M. & Galli, D. 2011, A&A, 530, A109 (PG11)
  40. Padovani, M., Galli, D. & Glassgold, A. E. 2013, A&A, 549, C3 (PGG13)
  41. Ramachandran, P. & Varoquaux, G. 2011, IEEE Computing in Science & Engineering, 13 (2), 40
  42. Rimmer, P. B., Herbst, E., Morata, O. et al. 2012, A&A, 537, 7
  43. Sano, T., Miyama, S. M., Umebayashi, T. et al. 2000, ApJ, 543, 486
  44. Santos-Lima, R., de Gouveia Dal Pino, E. M. & Lazarian, A. 2013, MNRAS, 429, 3371
  45. Seifried, D., Banerjee, R., Pudritz, R. E. et al. 2012, MNRAS, 423, 40
  46. Shu, F. H., Galli, D., Lizano, S. & Cai, M. 2006, ApJ, 647, 382
  47. Skilling, J. & Strong, A. W. 1976, A&A, 53, 253
  48. Takayanagi, M. 1973, PASJ, 25, 327
  49. Teyssier, R. 2002, A&A, 385, 337
  50. Umebayashi, T. & Nakano, T. 1981, PASJ, 33, 617
  51. Umebayashi, T. & Nakano, T. 2009, ApJ, 690, 69
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Comments 0
The feedback must be of minumum 40 characters
Add comment

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question