Massive core/star formation triggered by cloud-cloud collision: Effect of magnetic field

We study effect of magnetic field on massive dense core formation in colliding unequal molecular clouds by performing magnetohydrodynamic simulations with sub-parsec resolution (0.015 pc) that can resolve the molecular cores. Initial clouds with the typical gas density of the molecular clouds are immersed in various uniform magnetic fields. The turbulent magnetic fields in the clouds consistent with the observation by Crutcher et al. (2010) are generated by the internal turbulent gas motion before the collision, if the uniform magnetic field strength is 4.0 $\mu$G. The collision speed of 10 km s$^{-1}$ is adopted, which is much larger than the sound speeds and the Alfv\'{e}n speeds of the clouds. We identify gas clumps with gas densities greater than 5$\times$10$^{-20}$ g cm$^{-3}$ as the dense cores and trace them throughout the simulations to investigate their mass evolution and gravitational boundness. We show that a greater number of massive, gravitationally bound cores are formed in the strong magnetic field (4.0 $\mu$G) models than the weak magnetic field (0.1 $\mu$G) models. This is partly because the strong magnetic field suppresses the spatial shifts of the shocked layer that should be caused by the nonlinear thin shell instability. The spatial shifts promote formation of low-mass dense cores in the weak magnetic field models. The strong magnetic fields also support low-mass dense cores against gravitational collapse. We show that the numbers of massive, gravitationally bound cores formed in the strong magnetic field models are much larger than the isolated, non-colliding cloud models, which are simulated for comparison. We discuss the implications of our numerical results on massive star formation.


Introduction
Massive stars have fundamental influence over the interstellar medium and galactic evolution. They ionize the surrounding gas and deposit energies by strong stellar winds and supernovae. Massive stars supply new material to surrounding gas, which is available for nextgeneration stars. Massive star formation still remains poorly understood despite decades of work (Zinnecker & Yorke 2007;Tan et al. 2014). The accretion rate,Ṁ , onto a protostar is where c T = (k B T /m) 1/2 is the isothermal sound speed of gas of mean molecular mass, m, at gas temperature, T , k B is Boltzmann constant, and G is the gravitational constant (Shu 1977;Stahler et al. 1980). Stahler et al. (1980) proposed that a low-mass star is formed in a molecular core with T = 10 K. McKee & Tan (2002) proposed that a higher rate of gas accretion onto a protostar is needed to form a massive star in a molecular core, since ram pressure associated with the accretion can exceed the radiation pressure from the protostar.
Such a high accretion rate can be realized in a massive dense core with large turbulence or high Alfvén wave speed. The study of the formation and evolution of such massive cores remains an integral part of understanding the massive star formation process.
Cloud-cloud collisions (CCCs) are strong candidates for massive star formation. If two molecular clouds collide at supersonic speed, a shock wave is produced at the interface of the colliding clouds. Since the density in the shocked region increases further due to radiative cooling, a dense clump can be formed due to the enhanced self-gravity (Gilden 1984). If the sizes of the colliding clouds are different, a converging flow appears in the shocked layer, and it can increase the mass of dense clump (Habe & Ohta 1992). If the dense cores formed in the massive clump are massive enough, we can expect massive star formation in them.
Observational evidence of massive star formation by CCCs are reported in various star formation regions (e.g., Hasegawa et al. (1994); Furukawa et al. (2009); Ohama et al. (2010); Dewangan et al. (2016); Torii et al. (2017); Fukui et al. (2018a); Fukui et al. (2018b)). These are super star clusters, the H II regions, the Spitzer bubbles in the Milky Way disk and star formation regions in the Large Magellanic Cloud. Observational evidence consists of two distinct velocity components of molecular gas with a large velocity difference (more than 10 km s −1 ) and bridge features in the position-velocity diagrams of observed molecular gas in these regions. These features were found in the recent numerical simulations of CCCs (Takahira et al. 2014;Haworth et al. 2015a;Haworth et al. 2015b;Takahira et al. 2018). The typical collision speeds observed in these regions are in the range of 10 km s −1 to 20 km s −1 , which are much higher than the gas speeds due to internal motion in molecular clouds.
Recently, hydrodynamic (HD) simulations of CCCs have been carried out to study the formation of dense cores, assuming turbulence in the clouds (Takahira et al. 2014;Takahira et al. 2018;Shima et al. 2018). They discuss that dense cores formed in the shocked layer can increase their mass by accretion of dense gas if turbulence and magnetic fields support the dense cores. Simulations with greater spatial resolution are needed to study the effect of turbulence and magnetic fields on dense core formation.
Magnetohydrodynamic (MHD) simulations of colliding molecular clouds have been performed by several authors (Inoue & Fukui 2013;Chen & Ostriker 2014;Chen & Ostriker 2015;Wu et al. 2017a;Wu et al. 2017b;Inoue et al. 2018). Chen & Ostriker (2014) and Chen & Ostriker (2015) have shown that the magnetic fields do not affect the typical mass of pre-stellar cores formed in the colliding clumps in a large molecular cloud. They found the formation of low-mass stars in their simulations. They selected collision speed as large as the typical turbulent velocity observed in molecular clouds. Since the observed CCC speeds are much larger than the typical turbulent velocity, there is an obvious need to study the collision of magnetized clouds with the collision speeds more typical of those observed. Wu et al. (2017a) and Wu et al. (2017b) have shown that the magnetic fields do not affect star formation in the CCC in their numerical simulations with more typical collision speed of CCCs, and they investigated observational signatures and physical properties of dense regions in the colliding magnetized clouds.
In their simulations, they used a spatial resolution of ∼ 0.1 pc that may not be high enough to resolve the dense cores of which typical scale is ∼ 0.1 pc (Bergin & Tafalla 2007), and a higher spatial resolution is required to properly reproduce their physical properties. Federrath et al. (2011) proposed that a minimum resolution criterion of 30 cells per Jeans length in hydrodynamical simulations of self-gravitating gas in order to resolve turbulent motion on the Jeans scale. Similar number of cells should be used to resolve the dense cores. We use a minimum cell size of 0.015 pc after some numerical tests by different choice of the minimum cell size (see section 2.1) and define dense cores using a selection criteria in line with the typical density in observed dense cores, described in detail in section 3.1.2. Simulations of the collision of magnetized clouds with typical speed and sufficient spatial resolution were performed by Inoue & Fukui (2013) and Inoue et al. (2018). Their results are in favor of massive star formation due to the role of magnetic fields. Inoue et al. (2018) simulated a collision of dense regions in a giant molecular cloud (GMC), assuming an uniform initial magnetic field 20 µG perpendicular to the collision velocity. Though the studies mentioned above have reached sufficient resolutions and have used realistic collision speeds, there is a dearth of studies that have investigated this while considering alternative magnetic field configurations (both in strength and orientation with respect to the CCC axis) and the typical GMC scale CCCs.
We present a study of the effects of magnetic fields on the formation of massive dense cores in magnetized, turbulent, and colliding clouds using high spatial resolution simulations.
We perform simulations assuming magnetic fields of varying strengths and directions to understand the role of magnetic fields on massive dense core formation. In section 2 we describe the numerical method and models, in section 3 we present our numerical results, in section 4 we discuss our results, and in section 5 we summarize our study.

Numerical Method
We use simulation code Enzo, a three-dimensional MHD adaptive mesh refinement (AMR) code (Bryan et al. 2014). We assume ideal MHD in our simulations. The code solves the MHD equations using the MUSCL 2nd-order Runge-Kutta temporal update of the conserved variables with the Harten-Lax-van Leer (HLL) method and a piecewise linear reconstruction method (PLM). The hyperbolic divergence cleaning method of Dedner et al. (2002) is adopted to ensure the solenoidal constraint on the magnetic field.
We describe numerical methods of the cooling, the pressure floor, and the Alfvén speed limiter used in our simulations. Radiative cooling of gas is calculated down to 10 K by using the cooling table made by the CLOUDY cooling code (Ferland et al. 1998) with the solar metallicity and a density n H = 100 cm −3 . For example, the cooling time of gas with a density n H = 100 cm −3 from 100 K to 10 K is estimated to be less than 0.1 Myr by using the cooling table. Due to this rapid cooling of molecular gas, the dense gas reaches the typical cloud temperature of GMC, 10 K. Photoelectric heating with the rate of 1.2 × 10 −25 (n H /1 cm −3 ) erg s −1 cm −3 is applied to gas (Tasker & Bryan 2008). Self-gravity in the gas is included in our simulations.
The pressure floor is applied for the cell in the finest grid level in which gravitational energy of gas is greater than its internal energy (Machacek et al. 2001). The pressure floor kicks in at gas density of ∼ 10 −15 g cm −3 for a cell at gas temperature of 10 K in the finest grid level. Tests with a higher (lower) value of the pressure floor parameter by a factor of 10 (1/3) resulted in a very similar core population and probability density functions, and so our simulation results are unlikely to be influenced by numerical effects of our default pressure floor parameter. Alfvén where B is magnetic field strength, ρ is density, and CGS system of units are used in all equations related to magnetic field in this paper. The maximum Alfvén speed is set to be 20 km s −1 to avoid a very short time steps created due to high Alfvén speeds in low gas density regions, which is effectively limited by increasing the gas density in such regions. We find that the Alfvén speed limiter works in very small regions with much lower gas density than the initial density of our model clouds. The minimum density in our simulations is selected as the initial density of the ambient medium of 1.69×10 −23 g cm −3 as given in section 2.2.1. Tests with a higher value of the maximum Alfvén speed by a factor of 5 and with a lower value of the minimum density by a factor of 1/10 resulted in a very similar core population and probability density functions to our simulation results with our default values of the maximum Alfvén speed and the minimum density.

Initial cloud structure and collision setup
We adopt initial conditions for clouds based on properties of observed GMCs (Heyer et al. 2009;Murray 2011). Two uniform molecular clouds, a small cloud and a large cloud, are initialized with density 3.67 ×10 −22 g cm −3 of which free-fall time is 3.5 Myr and with masses 972 M and 7774 M , respectively. We stop our simulations at t = 3.0 Myr, which is earlier than the free-fall time. These cloud masses are rather small in comparison to observed GMCs in the Milky Way of which the mass range is 10 3 M to 10 6 M (Murray 2011). We select small clouds to achieve high spatial resolution needed to study the effect of the magnetic field on the formation of massive dense cores in the colliding clouds with a rather small simulation box that we describe in the next paragraph. We adopt initial temperatures of clouds as 68 K and 273 K for the small cloud and the large cloud, respectively. While such temperatures are sufficient to provide initial pressure support, the dense gas in the clouds rapidly cools down to 10 K due to the radiative cooling during evolution. We adopt turbulence in clouds (see section 2.2.2). Parameters for each clouds are summarized in table 1. Typical collision speed of 10 km s −1 is used. The ambient medium have a density of 1.69 × 10 −23 g cm −3 and a temperature of 800 K. This high density of the ambient medium is used to avoid high Alfvén speeds in the ambient medium.
Six different colliding clouds models are considered, each with differing initial magnetic field strength and direction (see section 2.2.2 and table 2). Our simulation domain encompasses (32 pc) 3 with root grids 128 3 , and we use four refinement levels based on the condition of minimum baryon mass of 0.05 M for refinement. This gives the minimum cell size of 0.015 pc at the maximum refinement level. We have tested our simulations with more higher refinement levels, five and six. The minimum cell size is 0.0075 pc for the refinement level five and 0.0037 pc for the refinement level six. The minimum cell size of 0.0037 pc satisfies the minimum resolution criterion in HD simulations of self-gravitating gas proposed by Federrath et al. (2011) for the typical core size of 0.1 pc. We find that core mass functions and core properties in simulations with these higher refinement levels are very similar to the simulation results with our default refinement levels.
Additionally, two different isolated, non-colliding cloud models are used for comparison with the results of colliding clouds models (e.g., the population of dense cores). This isolated cloud has a sum of the masses of both the small and the large clouds, and it also includes turbulent motion (see table 1). Two different initial magnetic field strengths are selected similarly to the colliding clouds (see section 2.2.2). We summarize the simulation results of the isolated cloud models in the Appendix.

Magnetic field and Turbulence in clouds
The clouds are immersed in an initial uniform magnetic field, B 0 , and turbulent motions develop inside them from t = 0 Myr to t = 0.5 Myr. Turbulent magnetic field is generated inside clouds before the small cloud begins to move towards the large cloud at t = 0.5 Myr along collision axis (positive x-axis of the simulation box). We select two strengths of B 0 , B 0 = 0.1 µG (weak) and B 0 = 4.0 µG (strong) and three directions of B 0 , which are parallel to the collision axis, perpendicular (along positive y-axis) to the collision axis, and oblique to the collision axis. The angle, θ, between B 0 direction and collision axis for oblique B 0 model is 45 • .
Additional isolated cloud models have B 0 with magnetic field strengths B 0 = 0.1 µG (weak) and 4.0 µG (strong) similar to those selected for colliding cloud models. The direction of B 0 is taken along positive y-axis of the simulation box. We name the simulation models, as shown in table 2.  Turbulent velocities are generated as to be consistent with the Larson relation (Larson 1981;Heyer et al. 2009) at t = 0 Myr, by imposing a velocity field with power spectrum v k 2 ∝ k −4 . We define k = (2π)(n xî + n yĵ + n zk )/(aR), where R is the radius of each cloud and a ∼ 4.
We use integers for n x , n y , and n z as n 2 min ≤ n 2 x + n 2 y + n 2 z ≤ n 2 max , where n min and n max is 6 and 12 for the small cloud and 8 and 15 for both the large cloud and isolated cloud, respectively. If the initial kinetic energy of turbulence is in virial equilibrium with the self-gravitational energy, the velocity dispersion, σ v , of clouds due to the turbulent motion is σ v = 0.86 km s −1 for the small cloud and σ v = 1.72 km s −1 for the large cloud. We slightly increase these values to σ v = 1.0 km s −1 for the small cloud and σ v = 2.0 km s −1 for the large cloud to be consistent with the observational results by Heyer et al. (2009). 12 8 4 0 4 8 12 16 x (pc) 16 + By 2 ) 1/2 , By/(Bx 2 + By 2 ) 1/2 ). Color bar of the gas density is shown on the right edge.
We show the effect of turbulence on the initial magnetic field before the collision starts. to the effect of turbulence generated, and the extended tail is due to the effect of self-gravity (Kritsuk et al. 2011;Takahira et al. 2014). The log-normal part is narrower in the strong B 0 than the weak B 0 . This result is qualitatively consistent with Padoan & Nordlund (2011) who simulated supersonic, self-gravitating, MHD turbulence. In the Ystrong model, the extended tail is much narrower than that in the Yweak model. This is due to the higher magnetic field pressure that suppresses density enhancement of gas in the Ystrong model compared to that in the Yweak model. spheres of radii equal to 85 % of the initial cloud radii centered at the initial cloud centers for PDFs.

Numerical Results
We show numerical simulation results of colliding cloud models. Isolated cloud numerical simulation results are used for comparison. Numerical simulation results of isolated cloud models are given in Appendix. to the collision axis located along the collision interface with scales of less than 1 pc at t = 2.0 Myr. Figure 5 shows the close-up slice images of gas (left panels) and magnetic field strength (right panels) near the shocked layer at t = 1.7 Myr (upper panels) and 2.2 Myr (lower panels). The spatial shifts should be due to the nonlinear thin shell instability (NTSI) (Vishniac 1994;Anathpindika 2010). Figure 5 shows that the spatial shifts develop with time between t = 1.7 Myr and t = 2.2 Myr. Dense gas concentrations are formed at the extremes of the spatial shifts of the shocked layer. The distribution of magnetic field strength is similar to the gas distribution as shown in figure 5. This implies that the magnetic field plays a minor role in the evolution of the shocked layer. This can be explained as follows. If a magnetic field is strong enough to control the gas flow, the gas flow along the magnetic field is easier than the transverse flow. In this case, the gas distribution can be very different from the magnetic field distribution. If magnetic fields are too weak to affect the gas flow, the gas flow will change the magnetic field structure, and the distribution of magnetic field strength can be similar to the gas distribution for a highly turbulent magnetic field as in the post-shock gas in the clouds.
At t = 2.5 Myr, the shocked layer shape is more concave than t = 2.0 Myr with additional substructures developing in the shocked layer. This is due to the shrinking of the central part of the shocked layer along the collision axis caused by converging flow in the shocked layer.
From t = 2.5 Myr to 3.0 Myr, the concave shape develops further with additional substructures developing in the shocked layer and dense mass concentrations are formed near x ∼ 5 pc and y ∼ 0 pc in figure 4, which is the bottom of the concave shape. A similar evolution is found in the numerical results of the other models with the weak B 0 , implying that the magnetic field plays only a minor role in the evolution of the shocked layer in the weak B 0 models.

Dense core formation and evolution
In order to study dense core formation and evolution, we define a dense core by the threshold density, ρ th = 5 × 10 −20 g cm −3 , which is in the range of the typical density of molecular cores (Bergin & Tafalla 2007). We define dense cores by following steps, 1) selection of cells with ρ ≥ ρ th as dense cells, 2) grouping together with neighboring dense cells and 3) rejection of those groups with cell number less than 27. This minimum cell number condition is used for a good resolution of dense cores. We show a cumulative dense core mass distribution, N (≥ M core ), which is number of cores with mass more than M core , at t = 1.5 Myr and 2.0 Myr to 3.0 Myr in figure 6.
We also show the cumulative dense core mass distribution of the isolated cloud model with the weak magnetic field by filled triangles in figure 6. Figure 6 clearly shows that more massive dense cores are formed in the colliding clouds model than the isolated model. The dense cores are already formed at t = 1.5 Myr and the core number rapidly increases from t = 2.0 Myr to t = 2.5 Myr in the colliding clouds. In the isolated cloud, the first dense core forms at later epoch than in the colliding clouds, and the total number of dense cores is also much smaller compared to the colliding clouds at each epoch shown. At t = 3.0 Myr, four massive dense cores with their mass more than 10 M are formed in the colliding clouds.
We check the gravitational stability of each dense core by comparing its turbulent energy, E turb , its magnetic field energy, E mag , and an absolute value of its self-gravitational energy, |E grav |. Here we ignore its thermal energy, since temperature of the dense core is nearly 10 K and the sound speed of gas at this temperature is much smaller than the turbulent velocity and the Alfvén speed. E turb is given for a dense core by where i is an index of a dense cell in the dense core, the sum is made over all cells belonging to the dense core, m i is the mass of the dense cell i and v mean is the mean velocity of the dense core given by E mag is given by where B i and V i are the magnetic field flux density vector and volume, respectively, of the dense core i. We estimate |E grav | by where G is the gravitational constant, M core is mass of the dense core, and R is given by where V core is total volume of the dense core. In figure 7, we show these energies. If (E turb +E mag ) ≤ |E grav |, we can expect that the dense core is gravitationally bound and we call such a core a gravitationally bound core. Since a dense core is defined by the condition of ρ ≥ ρ th , its  figure   3) as in the MHD bow shock by the solar wind past a planet (Slavin et al. 1984). The shocked layer produced by the CCC is much thicker than that seen in the Yweak model. This may be due to smaller M A and the larger magnetic pressure in the shocked layer in the Ystrong model than the Yweak model.  The density and magnetic field enhancements are not coincident, as shown in figure 9.
Distribution of the magnetic field strength is much smoother than the density enhancements.
This indicates that the turbulent magnetic fields on a large scale are enhanced by the CCC flow, and gas moves along the smaller scale magnetic fields to create further density enhancements in the shocked layer. In this way, a difference in density and magnetic field enhancements is produced.

Dense core formation and evolution
We show the time evolution of the cumulative core mass distribution in the Ystrong model in figure 10. The total number of dense cores formed in the Ystrong model is less than that in the Yweak model (see figure 6) during the early phases of collision (t = 1.5 Myr). This  is due to the suppression of the NTSI in the smaller scales in the Ystrong model. More massive cores are formed in the Ystrong model than the Yweak model at t > ∼ 2.0 Myr, and the number of the massive cores with masses more than 10 M is also greater than the Yweak model at t = 3.0 Myr. The massive core formation is due to gas accumulation to massive cores in the dense gas regions in the thick shocked layer, as shown in figure 9.
The time evolution of the magnetic field energy, the turbulent energy, and the absolute value of self-gravitational energy of each dense core with its mass is shown in figure 11 from t = 1.5 Myr to t = 3.0 Myr. At t = 1.5 Myr and 2.0 Myr, the magnetic field energies are larger than the absolute self-gravitational energies in all cores. However, at t = 2.5 Myr and 3.0 Myr, the absolute self-gravitational energies are dominant in massive cores. This change is mainly due to their mass increase by gas accumulation.
The time evolution of the top ten massive cores and the epoch for them to become gravitationally bound is shown in figure 12. In this figure, we also show merger trees that indicate dense core mergers. The mass growth of any given dense core is a combination of the accretion of surrounding gas and mergers with other dense cores. If the mass contribution by mergers is not enough to explain the mass growth of the core, the gas accretion to the core should be a dominant process for the mass growth. Figure 12 shows that merger events play more of a secondary effect in increasing the mass of the top ten most massive cores, with mass growth primarily a smoother function of time implying accretion dominated, evolution. Since the free-fall time of these massive cores is less than 0.3 Myr, we expect that protostars form in these cores.
We highlight the gravitationally bound cores using larger open circles in figure 10. As shown in figure 10, three massive cores become gravitationally bound at t = 2.5 Myr, although these are not gravitationally bound at t = 2.0 Myr. Their mass is larger than 10 M . At t = 3.0 Myr, most cores (9 out of the 10) with more than 10 M become gravitationally bound. The number of bound cores in the Ystrong model is much larger than the Yweak model in which only two bound cores with a mass more than 10 M are found at t = 3.0 Myr as shown in figure 6. This may be due to the thick shocked layer caused by the strong magnetic field in the Ystrong model. The NTSI grows faster for a small scale shift of the shocked layer. In the Yweak model, the NTSI develops in the small scale shift and results in the accumulation of low-mass gas at the extremes of the quasi-periodic shifts, as shown in figure 5. These gas concentrations move with large irregular velocities, which may suppress further gravitational gas accumulation to the gas concentrations. As a result, this mass growth by gas accretion will be suppressed. In the Ystrong model, since the shocked layer is thickened by the strong magnetic fields as shown in figure 9, and the dense cores have small irregular velocities by suppression of the NTSI, the dense cores can acquire mass by accretion in the thick shocked layer. The cumulative core mass distribution in the isolated cloud model with the strong magnetic field B 0 = 4 µG is shown for comparison in figure 10 using filled triangle symbols. The bound cores form earlier than the Ystrong model, and the mass of bound cores just formed is less than 3 M . We can expect intermediate-mass star formation in these cores in the free-fall time ∼ 0.3 Myr assuming that the gas mass in these cores accretes to form a single star. This is very different from colliding clouds models, which hosted bound cores formed with masses greater than 10 M .
The resulting formation of protostars can be studied using sink particle models, though this is beyond the scope of the work presented here (Federrath et al. 2010;Shima et al. 2018

Core mass distribution and gravitationally bound cores
We show the core mass functions of all models at t = 3.0 Myr in figure 13. We find many gravitationally bound cores in the strong B 0 models at t = 3.0 Myr. The core mass distributions of the weak B 0 models are very similar to each other. In each weak B 0 model, there is one exceptionally massive core. In the strong B 0 models, the number of dense cores with more than 10 M is larger than the weak B 0 models. This indicates that the strong B 0 contributes to the formation of a greater number of massive, dense cores. The number of massive, dense cores in the Xstrong model is slightly smaller than the Ystrong and XYstrong models, yet it is still greater than those in the weak B 0 models (Xweak, Yweak, and XYweak models). These results indicate that the strong magnetic field parallel to the collision axis (as in the Xstrong) is less effective in suppressing the NTSI and in keeping the shocked layer thick, compared with the strong magnetic field with orientations oblique or perpendicular to the collision axis (as in the Ystrong and XYstrong models).
The mass of the most massive cores attains more than 100 M in the last 1 Myr in all models, since we cannot find such massive cores at t = 2.0 Myr. Since the free-fall time of these cores is less than 0.3 Myr, we can expect rapid protostar formation before t = 3.0 Myr.
If massive stars are formed in those colliding clouds, we can expect very strong feedback from these massive stars, as shown by Shima et al. (2018).
In figure 14, we show the position of cores with masses greater than 10 M in the strong B 0 models, the Xstrong, Ystrong, and XYstrong models, and the weak B 0 model, the Yweak model, overlaid on the plot of column density looking from the collision axis (x-axis). Large crosses show cores with more than 100 M , and small crosses show cores with 10 M < M core < 100 M . These massive cores are distributed along the filaments with column densities greater than 10 −1 g cm −2 in the strong B 0 models. In the Ystrong and XYstrong models, the filaments hosting the massive cores are roughly perpendicular to the normalized mass-weighted magnetic field lines, of which directions are shown as unit arrows in figure 14.
If the core mass function, φ, is defined as where M core is the core mass, dN is the number of cores with masses between M core and M core +dM core , and γ is power index of core mass function, the cumulative core mass distribution, N (≥ M core ), is given by, where α = -(γ-1). The least square fits α using equation 9 with standard deviation for cumulative mass distributions of cores with masses greater than 10 M for all models at t = 3.0 Myr are shown in figure 13 . The power indexes of core mass functions, γ, are γ ∼ 1.3 -1.4 in the weak B 0 models and γ ∼ 1.5 -1.9 in the strong B 0 models. The strong B 0 models have slightly larger γ than that in the weak B 0 models. These γ values are similar to those of HD simulations performed by Takahira et al. (2014) and Takahira et al. (2018). These values are also closer to the observed power indexes of core mass functions (Ikeda et al. 2007;Uehara et al. 2019), indicating the possible integral role of magnetic fields in determining the observed core mass function.

Discussion
Our simulation results have shown that a greater number of massive cores (> 10 M ) form in the strong B 0 models than the weak B 0 models. In the weak B 0 models, NTSI develops in the shocked layer produced by the CCC and induces the quasi-periodic shifts of the shocked layer in the early stage of CCC. Gas concentrations develop at the extremes of the shifts. In these gas concentrations, dense cores of small mass form earlier than the strong B 0 models.
In the strong B 0 models, the turbulent magnetic fields suppress such small scale NTSI. The turbulent magnetic fields increase the typical scale of NTSI and the thickness of the shocked layer. The turbulent magnetic fields also contribute to the increase of mass of dense cores in the thick shocked layer. Both effects contribute to the formation of massive dense cores in the strong B 0 models. The suppression effect of the magnetic field on NTSI in a shocked region is studied by Heitsch et al. (2007). They have reported that magnetic fields parallel to the shock are more effective in suppression of the NTSI than magnetic fields normal to the shock.
This effect can be the reason why the direction of B 0 to the collision axis affects the core mass functions. (vx/(vx 2 + vy 2 ) 1/2 , vy/(vx 2 + vy 2 ) 1/2 ). Color bar on the right edge shows the gas density.
In figure 15, we show the concave structures of the shocked layer created by the small cloud penetration into the large cloud in the strong B 0 models at t = 2.0 Myr, as well as for the Yweak model for comparison (all shock fronts in the weak B 0 models look essentially identical).
The  We discuss a possible role of magnetic field on self-gravitational instability in the shocked layer. The gravitational instability condition for a disk with a magnetic field B n is given by where B n is the perpendicular magnetic field to the disk and Σ is surface density of the disk (Nakano & Nakamura 1978;Tomisaka & Ikeuchi 1983). A disk with a parallel magnetic field is gravitationally unstable for a perturbation with a wave number vector parallel to the magnetic field (Tomisaka & Ikeuchi 1983;Nagai et al. 1998). We estimate the gravitational unstable scale of this case according to Tomisaka (2014). In Tomisaka (2014), the maximum mass per unit length, λ max , of a gravitational equilibrium filament perpendicular to the magnetic field B is given as for the limiting case of magnetic field energy much larger than the internal gas energy in the filament, where R 0 is the filament radius. The gravitational instability condition of a filament with a width h and with surface density, Σ, is since λ of this filament is given as λ = hΣ. We estimate the averaged magnetic field, B , and the averaged surface density, Σ , of the shocked layer formed in our numerical simulations by using a thick disk region with a radius of 4 pc and with a thickness of 2.5 pc that can contain the shocked layer. In the Xstrong, Ystrong and XYstrong models at t = 1.7 Myr (nearly epoch at which the small cloud completely penetrates the large cloud), the averaged B x = 11.8 µG, 16.4 µG, and 15.8 µG, the averaged B y = 10.1 µG, 25.9 µG, and 21.4 µG, and the averaged surface density is Σ = 0.018 g cm −2 , 0.015 g cm −2 , and 0.016 g cm −2 , respectively. We note that B x and B y in the weak B 0 models are much smaller than in the strong B 0 models. Using these values, from equation (10), we find that B x cannot suppress gravitational instability of the shocked layers in the strong B 0 models. Using B y and Σ , we get h min as 0.65 pc, 2.0 pc, and 1.6 pc for the Xstrong, Ystrong, and XYstrong models, respectively. These values of h min are less than the typical, lateral size (∼ 8 pc) of the shocked layers. This suggests that the shocked layers are gravitaionally unstable and that a typical mass of a fragment formed by the gravitational instability is Σh 2 min , which is much larger than masses of dense cores formed in early stage of the collisions of clouds. From the discussion, main effects of magnetic fields on the CCCs are the suppression of the NTSI and the increase of the thickness of the shocked layers formed in strong B 0 models as shown in section 3.
We estimate magnetic field strength B 0 that can suppress NTSI in the shocked layer from our simulation results. We can expect the magnetic field suppression effect on the NTSI, where B and ρ are magnetic field strength and gas density, respectively, in the shocked layer, λ is the typical scale of the NTSI, and ∆v is the perturbed velocity induced by the NTSI. If the magnetic pressure enhanced by the shock compression is dominant in the shocked layer, we can expect and from the MHD shock wave condition in the rest frame of the shock front, where v sh and v are pre-shock and post-shock gas velocities, respectively, and αB 0 and ρ 0 are pre-shock parallel component of the magnetic field to the shock front and gas density, respectively. Here, α is an enhancement factor of the parallel component of the magnetic field to the shock front induced by the turbulent motion in the clouds before the collision.
The suppression condition of NTSI can be given as, from equation (13), using equation (14) and equation (15). If ∆v/v sh ∼ 0.2 and α ∼ 1, we have Here, we use that v sh is roughly half of the collision speed. This estimated value is consistent with our numerical results, since the strong B 0 is much larger than this value and the weak B 0 is much less than this value.
If the magnetic field pressure is dominant in the shocked layer than the effects of the turbulent motions and the thermal gas, we estimate magnetic field strength, B, and gas density, ρ, in the shocked layer as, and ρ = 4.39 × 10 −21 ρ 0 3.67 × 10 −22 g cm −3 using equation (14) and equation (15).
After the whole small cloud penetrates the shocked layer, the shocked layer will change its structure in the free-fall timescale of the shocked layer ∼ 1 Myr, since the ram pressure by the small cloud gas does not push the shocked layer after the whole small cloud penetrates the large cloud. In this stage, dense core formation and core mass evolution proceed to form massive bound cores in the timescale of t ff , as shown in section 3.2. We thus propose that if the shocked layer moves out of the larger cloud in less than t ff after the whole small cloud penetrates, then such massive core formation will not proceed. We have adopted a collision speed of 10 km s −1 in this study. We briefly discuss the consequences of a higher collision speed, since the observed collision speeds are in the range of 10 km s −1 to 20 km s −1 (Fukui et al. 2018a). If we assume the collision speed of 20 km s −1 , v sh = 10 km s −1 , then the condition of suppression of NTSI will be from equation (18). This result indicates that stronger B 0 than used in this study is required to suppress the NTSI. The free-fall time of the shocked layer is t ff ∼ 0.7 Myr for the collision speed of 20 km s −1 . If the penetration time of the small cloud is and the crossing time of the small cloud to the large cloud is t cross is comparable to sum of t penetration and t ff . This means that larger cloud sizes of the large cloud are needed to induce massive dense core formation. In our galaxy, various sizes and masses of GMCs are observed. Magnetic field strength depends on location in our galaxy (Beck 2015). We will extend our study to a higher collision speed case with larger cloud sizes and stronger magnetic fields in our future works. We will study protostar formation and stellar feedback effects on massive star formation by CCCs, using sink particles in our future works.

Summary
We have performed magnetohydrodynamic simulations of the cloud-cloud collision to study the role of the magnetic field on massive core formation in the colliding clouds. We selected two clouds with masses of 972 M and 7774 M with the typical density of giant molecular clouds and with internal turbulence such that the clouds are in the virial equilibrium.
Two cases of uniform magnetic field strengths, B 0 = 4.0 µG (strong) and 0.1 µG (weak), and three cases of uniform magnetic field directions, parallel, perpendicular, and oblique to collision axis were studied. Magnetic fields were modified by internal turbulent motion in the clouds.
The distribution of magnetic field strength and gas density in the clouds in the strong B 0 model is consistent with the observed relation by Crutcher et al. (2010). The small cloud is given a collision speed of 10 km s −1 after the turbulent magnetic field generation in the clouds. We have also simulated the evolution of the isolated clouds with the same uniform magnetic fields as in colliding clouds for comparison. Our main conclusions are as follows: 1. In the weak B 0 models, quasi-periodic shifts with small size appear in the shocked layer formed by cloud-cloud collision and develop with time. The quasi-periodic spatial shifts may be due to nonlinear thin shell instability. Dense cores are formed at the extremes of the shifts. In the strong B 0 models, such shifts are suppressed by the stronger magnetic field, and a greater number of massive dense cores are formed than the weak B 0 models. The number of massive bound cores in which self-gravitational energy dominates over turbulent energy and magnetic field energy is also larger than the weak B 0 models. In the massive bound cores with more than 10 M , we can expect massive star formation, since the free-fall time of these cores is less than 0.3 Myr. In isolated cloud models, the bound cores form earlier and are less massive than the colliding clouds models. Since their masses are less than 3 M and their free-fall times are less than 0.3 Myr, we can expect only intermediate-mass star formation in these cores.
2. The cumulative mass distributions of dense cores formed in our simulation models clearly show that a greater number of massive cores are formed in the strong B 0 models than in the weak B 0 models.
3. In the strong B 0 models, massive cores distribute in dense gas filaments of which directions are roughly normal to the direction of B 0 except for the Xstrong (strong B 0 parallel to collision axis) model.

4.
We give a simple analytic model for the magnetic field strength needed to suppress the instability of the shocked layer formed by colliding clouds that leads to small mass core formation. The magnetic field strength is related to collision speed and cloud size. Testing this model further will be the subject of future works.  show higher gas density contrast in the ISweak model than the ISstrong model. The PDFs at later epochs (t = 2.5 Myr and 3.0 Myr) show power-law tail due to the effect of self-gravity and the formation of dense gas with a density greater than ρ th in the ISweak and ISstrong models.