ON THE COAGULATION AND SIZE DISTRIBUTION OF PRESSURE CONFINED CORES

Dense molecular cores are long known to be embedded in molecular clouds. These cores are particularly interesting in the context of star formation because they are thought to be the progenitors of individual stars or small stellar systems. Recent infrared measurements of dust extinction, CO and NH₃ maps provide means for observers to systematically and quantitatively analyze the physical properties of these cores. Among the first studies of this kind is that carried out by Alves et al. (2007) on the Pipe Nebula. These observations revealed a population of dense, starless cloud cores. The similarity between the density, temperature, and thermal pressure of cores at different locations in the nebula suggests that they are confined by a global external pressure, although the most massive cores have sufficient self-gravity to hold themselves together. Prestellar cores are also found in the Aquila Nebula, Polaris Nebula, and Ophiuchus main cloud (Könyves et al. 2010; André et al. 2010, 2007), though there exists some diversities in the environment and the observed mass distribution.

One scenario for the confinement of these cores in Pipe is the possible existence of a three-phase medium which includes (1) a hot component ionized by the nearby B2 IV β Cephei star, θ Oph, (2) a cold atomic gas (∼100 K, we refer to it as warm medium to distinguish with the molecular phase), and (3) the cold molecular cores (∼10 K) (Gritschneder & Lin 2012). We associate this multi-phase medium in the Pipe Nebula with the byproduct of a thermal instability which is excited when thermal energy in a slightly cooler region is radiated away faster than its surrounding (Field 1965). A quasi-hydrostatic pressure balance is maintained over regions smaller than the field length scale for sound waves to travel within the instantaneous cooling time. These regions undergo isobaric cooling such that their density contrast to the background grows in time. Meanwhile the field length declines with decreases in both the sound speed and cooling timescale. Eventually, density difference between the cold regions and the medium which engulf them either become nonlinear or reach an asymptotic limit over regions larger than the shrinking field length (Burkert & Lin 2000).

Thermal instability is determined by local conditions and its growth rate does not depend on the wavelength. Nevertheless, the initial perturbation spectrum does determine the range of length scale over which the density contrast may become nonlinear during the onset of thermal instability. The interstellar medium is turbulent and attains a power-law perturbation spectrum, similar to the Kolmogorov law, over an extended range of length scales (Lazarian 2009). However, the cores' mass function (CMF) in the Pipe Nebula appears to resemble more closely with the stellar initial mass function (IMF; Alves et al. 2007; Salpeter 1955; Kroupa 2002) than with the single power-law spectrum of the interstellar medium. The similarity between the CMF and IMF has also been observed in the Aquila Nebula and Ophiuchus main cloud (Könyves et al. 2010; André et al. 2007). Palaris is an exception of this resemblance (André et al. 2010). The relationship between the CMF and IMF has been investigated by many authors using different approaches (Smith et al. 2009; Dib et al. 2010; Veltchev et al. 2011; Anathpindika 2011). While there remains diverse opinions on which physical processes affect the evolution of the CMF and the transition to the IMF, it is generally recognized that such a process critically determines star formation rate and efficiency.

Many past discussions have been focused on the origin of the shape of CMF (see Padoan & Nordlund 2002; Hennebelle & Chabrier 2008; Smith et al. 2009; Offner et al. 2008; Dib et al. 2008, etc.), which may be related to the turbulent and self-gravitating flows (McKee & Ostriker 2007). Such flows are complex because nonlinear structure of multi-phase medium extends over a large dynamical range, from the inertial to dissipative scales, throughout the star-forming region. Several numerical methods, including different versions of the adaptive mesh refinement scheme, have been constructed to simulate the emergence of bound and unbound clusters of cores in a turbulent environment. However, none of these attempts can fully resolve cloud complexes in which many core coexist with their confining medium. In addition, these simulations are time consuming and have been carried out for limited sets of model parameters. It is not clear whether they depend on the assumed history of the system or are applicable in general.

In this paper, we aim to explain the present-day CMF observed in the Pipe Nebula, which is just prior to the onset of star formation. We suggest that the characteristic CMF spectrum is established through the coalescence and ablation of cores rather than precipitated directly from a turbulent gas cloud through thermal instability. In Section 2, we construct a coagulation equation which takes into account various relevant physical processes. We design a numerical scheme (presented in the Appendix) to efficiently solve this integral differential equation. This method allows us to follow the cores and the surrounding medium as two separate, but coexisting, fluids. We describe our detailed modeling based on the properties relevant to the Pipe Nebula in Section 3. In Section 4, we show that we can reproduce the observed CMF, including the observed peak of mass distribution which occurs in a mass range three times heavier than that for IMF. Finally, we summarize our findings and discuss their implication in Section 5.

We assume all but the largest cloud cores are confined by the pressure of an external hot and warm medium. For simplicity, we do not distinguish between contributions from the hot and warm medium. In this paper, we are mainly concerned with the interaction between the cloud cores and the inter-core gas. In this context, we take into account two main physics processes: (1) close-encounters and collisions between the cold cores, due to the velocity dispersion between the cores and (2) the dynamic interaction between the cold cores and the background gas due to the drag on the cores by the inter-core gas. Process (1) leads to the coagulation between the cores, which induces a shift in the mass distribution toward higher masses. Process (2) leads to ablation of the cores, due to the Kelvin–Helmholtz (K-H) instability which causes the mass distribution to evolve toward the lower masses. In this section we will provide the basic equations for the two processes and order of magnitude estimation of the timescale. In addition, we briefly discuss the contribution of conduction at the end of this section.

2.1. Coagulation Equation

Since we are primarily interested in the bulk properties rather than the detailed phase space distribution of the cloud population, we adopt a conventional approach to examine the evolution of the distribution function

$$\begin{equation} f(m, t) = \frac{dN (m, t)}{dm} \end{equation} \tag{ 1 }$$

where N(m, t) is the number of cores in the mass range m and m+dm at a time t. We define that the minimum and maximum masses to be m_μ and m_max, respectively. The rate of change d(dN/dm)/dt can be expressed in term of a general coagulation equation (see Murray & Lin 1996; Shadmehri 2004; Dib et al. 2007) such that

$$\begin{eqnarray} \left(\!\frac{df(m)}{dt}\!\right)_{\rm coag}\! &=& \frac{1}{2}\!\int _{m_{{\rm min}}}^{m-m_{{\rm min}}}f(m^{\prime })\!f(m-m^{\prime })G(m,m-m^{\prime })\,dm^{\prime } \nonumber\\ &&-\int _{m_{\rm min}}^{m_{\rm max}}f(m)f(m^{\prime })G(m,m^{\prime })\,dm^{\prime }. \end{eqnarray} \tag{ 2 }$$

For computational simplicity, we assume there is no mass loss during each collision. We also assume the physical condition of the inter-core gas, cores' mass and kinematic distribution are homogenous and isotropic. These assumptions can be relaxed in a series global simulations to be presented in the future. We adopt an average velocity dispersion σ for the cores, which is observed to be σ ∼ 1 km s⁻¹ in Pipe Nebula (Onishi et al. 1999). Hydrodynamic simulations study the kinematics of molecular cloud cores in the presence of turbulence also suggest similar values for the velocity dispersion between cores (Offner et al. 2008). Considering gravitational focusing for massive cores, we adopt a collisional kernel

$$\begin{eqnarray*} &&G(m,m^{\prime })=(1+\Theta) \pi (R_c(m)+R_c(m^{\prime }))^2\sigma \end{eqnarray*}$$

where $\Theta =({c_s^2}/{2\sigma ^2})(({m^{2/3}+m^{\prime 2/3}})/({m_{\rm BE}^{2/3}}))$, following Equation (7.195) in Binney (2008). M_BE is the Bonnor–Ebert mass for the cores:

$$\begin{equation} M_{\rm BE}=1.15\, M_{\odot }\left(\frac{c_s}{0.2\,{\rm km\, s^{-1}}}\right) \left(\frac{P_{\rm ext}/\kappa }{10^5\, {\rm K\, cm^{-3}}}\right)^{-0.5} \end{equation} \tag{ 3 }$$

where P_ext is the external pressure of the inter-core gas and κ is the Boltzmann constant (Bonnor 1956; Ebert 1955).

Cores gain a significant fraction of their own mass through cohesive collisions mostly with cores of comparable sizes with a top heavy distribution of CMF. The characteristic growth timescale for cores to double their masses is τ_coag = λ/σ where the cores' mean free path λ = 1/(πR_c(m)²n_c). n_c is number density of cores with sizes R_c(m) and mass m. The number of cores decreases with merger events.

Under these conditions, for a typical 1 m_☉ core, we estimate

$$\begin{equation} \tau _{\rm coag}\sim \frac{1}{{\pi }R_c^2n_{c}\sigma } = 1.3\ {\rm Myr}\left(\frac{R_c}{0.1\ {\rm pc}}\right) \left(\frac{f_f}{0.1}\right)^{-1}\! \left(\frac{\sigma }{\ {\rm km\,s^{-1}}}\right)^{-1}. \end{equation} \tag{ 4 }$$

In the above estimate, we replace n_c with $3 f_f / 4 \pi R_c^3$, where the magnitude of the volume filling factor f_f can be obtained directly from the observed area filling factor.

2.2. Core Ablation as a Result of the Kelvin–Helmholtz Instability

The intrinsic ablation of the cores also contributes to the coagulation equation. When cores attain a relative velocity with respect to the inter-core gas, they encounter a drag force. The shear at the gas-core interface leads to the excitation of a K-H instability (Murray et al. 2004). Within the growth timescale τ_KH for the longest unstable wavelength (a significant fraction of R_c(m)), about half of the core mass break up into smaller cores (Murray et al. 1993). Cores' mass loss rate is ∂m/∂t ≃ −m/τ_KH, and the mass at time t due to ablation alone can be approximated by m ∼ m(t₀)exp (t₀ − t)/τ_KH. For a density contrast between the cores and the diffuse background, D_ρ,

$$\begin{eqnarray} \tau _{\rm KH}&=&\left(\frac{\sigma }{R_c(m)D_{\rho }^\frac{1}{2}}\right)^{-1}\nonumber\\ &=&1\ {\rm Myr}\left(\frac{\sigma }{\rm km\, s^{-1}}\right)^{-1}\left(\frac{R_c(m)}{0.1\ {\rm pc}}\right)\left(\frac{D_{\rho }}{100}\right)^{1/2} \end{eqnarray} \tag{ 5 }$$

In the above expression, we assume the relative velocity between the cores and the gas v_r to be comparable to the cores' velocity dispersion σ. It is possible that the turbulence speed of the gas v_g is directly determined by some other stirring mechanism such as external magnetic fields (Lazarian & Vishniac 1999). On the turnover timescale of the largest eddies (comparable to the size of the Pipe Nebula R_p), gas drag would lead to a terminal relative velocity v_r ∼ (D_ρR_c(m)/R_p)^1/2v_g (Lin & Murray 2000) and

$$\begin{eqnarray} &&\tau _{\rm KH}\sim \frac{R_c(m)^{\frac{1}{2}}R_p^{\frac{1}{2}}}{v_g}\nonumber\\ &&\quad =1\ {\rm Myr}\left(\frac{R_c(m)}{0.1\ {\rm pc}}\right)^{\frac{1}{2}}\left(\frac{R_p}{10\ {\rm pc}}\right)^{\frac{1}{2}} \left(\frac{v_g}{1\ {\rm km\,s^{-1}}}\right)^{-1} . \end{eqnarray} \tag{ 6 }$$

We include this ablation effect into the analysis of CMF's evolution as a source (due to the increase in the number of cores in the mass range M − M_BE) and sink terms (due to the decrease in the number of cores of mass M by ablation) such that

$$\begin{equation} \left(\frac{df(m)}{dt}\right)_{\rm abla}=\int ^{M_{\rm BE}}_{m}\frac{f(m^{\prime })}{\tau _{{\rm KH},m^{\prime }}}\frac{dm^{\prime }}{m} -\left(\frac{m-m_\mu }{m}\right)\frac{f(m)}{\tau _{{\rm KH},m}}. \end{equation} \tag{ 7 }$$

K-H instability would be suppressed by the cores' self-gravity if their mass M > M_BE (Murray et al. 1993). They undergo gravitational contraction and collapse rather than ablation. The magnitude of M_BE is approximately 2 M_☉ in the Pipe Nebula case (Lada et al. 2008). This critical mass for ablation partly causes a peak in our simulated core mass distribution function.

2.3. Other Relevant Processes

3.1. Observed Properties of Pipe Nebula

3.2. Numerical Setup

Instead of investigating the evolution of individual cores, we evolve the mass distribution function with the processes described in Section  2 with the following set up.

In our standard model, we assume the warm medium is uniformly distributed with a total mass of 10⁴ M_☉ within a sphere of radius ∼5 pc. The cores have an internal density of 7 × 10³ cm⁻³ so that their M_BE ∼ 2 M_☉, which is consistent with that derived from the observation of starless clouds in the Pipe Nebula (Lada et al. 2008).

The initial mass distribution of the cores is set up with a single power-law $dN/d\log M \propto M_c^{-\alpha }$. Analytical calculations predict the spectrum of non-self-gravitating structures to be $dN/dM \propto M^{2-n^{\prime }/3}$, where n' is the three-dimensional power spectrum index of the log density field Hennebelle & Chabrier (2008). When taking n' to be the Kolmogorov index 11/3, we found α to be 0.67. We only allow cores in the mass range between M_min (set to be 0.05 M_☉) and 1 M_☉ at the onset of the calculation. Their total mass is set to be 250 M_☉. For the standard model, we adopt D_ρ = 15 and a constant (everywhere and throughout the calculation) transonic velocity dispersion (∼1 km s⁻¹) for all cores, regardless of their size (see Table 1). For a robustness test, we carried calculations with other values of model parameters such as (1) the cores' velocity dispersion, ranging from 0.5 km s⁻¹ to 2 km s⁻¹ and (2) the slope of cores' initial mass distribution. The final result of the core distribution does not strongly depend on the above parameters in the range we tested.

Table 1. Parameters in Standard Model

mtotal/mtotalc	ρc	dρ = ρc/ρb	dt	mμ	mmax	σ	MBE
104 M☉/250 M☉	3 × 10−20 g cm−3	15	10 yr	0.05 M☉	100 M☉	1 km s−1	2 M☉

Download table as:  ASCIITypeset image

The high-mass (M_c > M_BE) cores are gravitationally unstable and they undergo gravitation contraction and eventually will evolve into stars. Based on an empirical model (Krumholz & Tan 2007), we adopt a simplified prescription for the rate of star formation ${\dot{N}} (M_c) = N (M_c)/\tau _\ast$. We set τ_* to be a M_c-independent characteristic timescale which is two orders of magnitude longer than their dynamical freefall timescale. We do not account for mass loss and possible feedback to background mass during the star formation process. The star formation process contributes equally to all the cores exceed M_BE. During the entire evolution the total mass of the stellar population is negligible compared with that of the cores.

We solve the combination of coagulation (Equation (3)) and ablation (Equation (7)) equation as described in the Appendix.

$$\begin{equation} \frac{df(m)}{dt}=\left(\frac{df(m)}{dt}\right)_{\rm coag}+\left(\frac{df(m)}{dt}\right)_{\rm abla}. \end{equation} \tag{ 8 }$$

The solution is then corrected with the effect of other processes (condensation, evaporation and star formation) in every time step. During their coagulation and ablation, the cores' mass is conservatively redistributed. The evaporation and condensation process do not significantly modify the slope of CMF in the most relevant mass ranges (0.3–10 M_☉). It mainly changes the mass ratio between total core mass and the warm gas mass. While the mass fraction between warm medium and cold cores evolves with time, the total mass in the two components is essentially constant. The determining factors of CMF are still dominated by coagulation and ablation over condensation, evaporation, and star formation.

The general shape of the observed CMF in the Pipe Nebula matches that of the IMF of young stellar clusters obtained by Kroupa & Boily (2002). But the cores' mass at the peak of the CMF is around the Bonnor–Ebert mass, a factor of three larger than that of the stellar IMF (at least in the Pipe Nebula). We simulate models to reproduce the observed CMF.

In our model, we adopt a set of idealized initial conditions with a population of low-mass cores as described in Section 3.2. We found that the asymptotic mass function of the cores does not strongly depend on this choice of initial condition, albeit in the standard model, it takes ∼1 Myr for the cores with M_c > M_BE to emerge and another ∼2–3 Myr for the CMF to establish a smooth distribution. Due to the low star formation efficiency in our model, stellar mass in the system is almost zero until around 5–6 Myr, at which we terminated the evolution for not many propostellar cores have been found in Pipe yet.

We first compare the evolution of the simulated CMF from the standard run with the observed CMF in the Pipe Nebula. The darkest black line in Figure 1 represents the stage after the emergence of gravitationally unstable cores but before the onset of star formation, as is observed in the Pipe Nebula today. For the mass range 0.3–10 M_☉, this black line matches closely with the observed core distribution (gray line). We reproduced the broad peak of distribution around the Bonnor–Ebert Mass, a flat distribution around 1 M_☉, and power-law cut offs at both low mass end and high mass end. The mismatch at the low mass end M_c < 0.3 M_☉ is partly due to observation bias as well as the simplified evaporation model we use. There are big uncertainties at the high mass end M_c > 10 M_☉, for in this region, both observation and our model need to deal with the problem of small number statistics.

Zoom In Zoom Out Reset image size

We note that the cores at the high mass end of the CMF are gravitationally bound. Their free-fall timescale is typically less than 1 Myr. In order to match the shape of the CMF and accounts for the absents of stars in them, we adopt the assumption that they evolve on a timescale which is ∼30 Myr. In the case of Pipe, magnetic field ∼17–65 μ G is observed inside the nebula, which might play an important role in the pressure budget (Alves et al. 2008). The existence of magnetic field will provide additional pressure support against gravity so that the cores might collapse on a much longer timescale compare to free-fall time.

We show the total mass evolution of both components in Figure 2. During the early stage, when the total core mass is concentrated in small mass end, the effect of evaporation win over condensation, so that the total mass of core decrease in the first 1 Myr. After the CMF approaches from an initial arbitrary distribution to a broad self similar form, most of cores are not affected by evaporation anymore. The total mass of cores increase due to condensation. Since we terminated the evolution at the stage when stellar mass can be negligible, the background gas density show inverse evolution path compare to the total core mass.

Zoom In Zoom Out Reset image size

By comparing the evolution of rate of coagulation and ablation at different mass bins in Figure 3, we can see the contribution of each processes to the final distribution. In the low mass end, coagulation contributes as sink term and ablation contributes as source term. Although ablation always dominates coagulation in mass range smaller than 0.2 M_☉, the cores of mass are small enough so that they are subjected to the influence of evaporation. Thus the overall population in this low mass range is reduced during the evolution. From 0.2 M_☉ to 1 M_☉, the number density decrease at almost constant rate since both coagulation and ablation act as sink terms. From 1 to 2 M_☉, coagulation starts to contribute positively to the mass distribution. For massive cores, they are determined by a self similar solution from a coagulation equation only. As the system evolve, the rate of ablation increase due to decrease of background density, the peak of gain from coagulation moves toward higher mass along with the total distribution.

Zoom In Zoom Out Reset image size

There are some freedom in the detailed prescription of the ablation process. In reality, the outer shell of a gravitationally bound Bonnor–Ebert sphere can still be stripped by K-H instability. However, as shown in Murray et al. (1993), the mass loss rate of core is really sensitive to the strength of gravitational field. Only 2% of total mass got stripped by K-H instability over 3.2 τ_KH in their simulation with a marginally gravitational bound cloud. In our model, for the gravitational bound clouds, the mass gain from coagulation will dominate the mass loss due to the stripping of surface layer, so that a strict terminate of K-H process at M_BE could be a good enough assumption.

Note that the timescale for transition from the initial to self similar mass distribution depends mostly on the initial filling factor and the velocity dispersion. The poorly known initial condition introduces additional uncertainties in the entire evolution timescale for the system.

We explore a range of initial conditions in Figure 4. We vary one parameter at a time, with all the other set ups the same as the standard model. We terminate all the runs at 5 Myr to make the comparison easier. The left panel of Figure 4(a) presents the result (at 5 Myr) from runs initialed with different core distribution at zero age. The power-law index of distribution dN/dlog M∝ = M^−α is varied from 0 to 1.33. As expected, a flatter initial distribution (lighter lines) assists coagulation so that the result distribution at 5 Myr is shifted toward higher mass end. The initial slope does not strongly influence the shape of final distribution. In contrast, as shown in the right panel, Figure 4(b), a different velocity dispersion (σ varied from 0.5 to 2 km s⁻¹) changes both the overall shape and the amount of mass in the massive cores. The slope of final distribution between 0.1 and 2 M_☉ is mainly due to the balance of coagulation and ablation. As we increase the velocity dispersion, the slope varies in this mass range. A higher velocity dispersion will also help with the emerge of the massive cores, lead to broader distribution with more mass in the high mass end. Nevertheless, in the allowed parameter space, the final CMF is weakly dependent on initial conditions.

Zoom In Zoom Out Reset image size

In other star formation regions, the environment is somewhat different from the case of Pipe. Ophiuchus has a much smaller global velocity dispersion σ < 0.4 km s⁻¹. If only consider coagulation and ablation, the evolution timescale will be longer, making them unlikely to be the only cause of the distribution. Other processes such as global gravitational contraction and accretion onto the core complex may also play important roles in the cores' velocity dispersion. But as we shown in Figure 4(b), within 5 Myr, a system of velocity dispersion of σ = 0.5 km s⁻¹ can still evolve into a flat and broad peak around 1–2 M_☉, demonstrating the two process could still contribute to the final distribution in such a system if the stellar less cores have longer lifetime due to the pressure support of either turbulence or magnetic field. Furthermore, the cores' size-dependent velocity dispersion is expected to have some influence on the power law index of CMF. This size dependence may originate from the cascade of turbulence in the warm gas and its ram pressure acceleration of the cores. Besides velocity dispersion, the diversity in Bonnor–Ebert mass in different regions is also likely to cause the diversity in the CMF.

In several nebula, including the Pipe, Aquila, Polaris, and Ophiuchus, a population of starless cores have been observed, at the low mass end are confined by external pressure where as the most massive cores are bound by their own gravity. The origin of these CMFs is a crucial issue in understanding of star formation in these regions.

In this paper, we present our calculation based on a two-phase-media model. We show that the dominant physical processes in this environment are coagulation, ablation and thermal conduction. We demonstrate an IMF-like CMF can be generated with a set of appropriate values. Through the exploration of a range of parameters reasonable for the physical condition of Pipe, we find this result is quite general. This robust form of the CMF is an indication that it is determined by a dynamical equilibrium, i.e., a balance between two competing processes: ablation and coagulation.

The simulated CMF closely matches the observed CMF in Pipe, and coincidentally matches the Kroupa IMF of the young stars. This result suggests that the ablation clouds occur prior to their collapse and the onset of star formation. The peak of the CMF is established near the cores' Bonnor–Ebert mass mainly because ablation is suppressed and collisional cross section is enlarged from their physical size for cores with M_c > M_BE. A similar peak is observed in the stellar IMF, albeit with a corresponding stellar mass which is one-third that of the most populous cores. If the a few of these cores capture a large fraction of their masses, they would form massive stars with copious sources of ionizing photons. Their feedback on the background gas and the nearby cores may modify the Bonnor–Ebert mass and lead to a reconciliation between the stellar IMF and the cores' CMF. We will explore these issues in our next contribution.

We thank Drs. Matthias Gritschneder and Herbert Lau for useful conversations. We also thank the anonymous referee for valuable comments. D.N.C.L. acknowledges support from NASA grant NNX08AL41G.

In general, the evolution of f(m) can be determined by the equation:

$$\begin{eqnarray} \frac{\frac{dN}{dm}}{dt}&=&\frac{df_j^i}{dt} =\frac{1}{2}\int _{m_\mu }^{m_j-m_\mu }f_k^{i+1}f_{j-k}^{i+1}G_{k,j-k}\,dm_k\nonumber\\ && -\int _{m_\mu }^{m_{{\rm max}}}f_j^{i+1}f_k^{i+1}G_{j,k}\,dm_k \nonumber \\ &&+\int ^{m_{{\rm BE}}}_{m_j}\frac{f_k^{i+1}}{\tau _{{\rm KH},k}}\frac{dm_k}{m_j} -\left(\frac{m_j-m_\mu }{m_j}\right)\frac{f_j^{i+1}}{\tau _{{\rm KH},j}} \quad \end{eqnarray} \tag{ A1 }$$

In the function, the superscript denotes the time levels and the subscript stands for different mass grids. The first two terms are due to coagulation, and the last two terms describe ablation. Change m into non-dimensional x (x_j = m_j/M_☉), and collect terms,

$$\begin{eqnarray} \frac{df_j^i}{dt} &=&A\left[\int _{x_\mu }^{x_j-x_\mu }f_k^{i+1}f_{j-k}^{i+1}\left(x_k^\frac{1}{3}+x_{j-k}^\frac{1}{3}\right)^2 \,dx_k\right. \left.\right.\nonumber\\ &&\left. -\, 2\int _{x_\mu }^{x_{\rm max}}f_j^{i+1}f_k^{i+1}\left(x_k^\frac{1}{3}+x_j^\frac{1}{3}\right)^2\,dx_k\right] \nonumber \\ &&+C\left[\int ^{x_{\rm BE}}_{x_j}\frac{f_k^{i+1}}{x_k^{\frac{1}{6}}}\frac{dx_k}{x_j}-\left(\frac{x_j-x_\mu }{x_j}\right)\frac{f_j^{i+1}}{x_j^{\frac{1}{6}}}\right]. \quad \end{eqnarray} \tag{ A2 }$$

In Equation (A2), we have

$$\begin{equation} A=\frac{\pi }{2}{f_{f}}\frac{\sigma }{R_{c0}},\,\quad \quad C=\tau _{\rm KH}^{-1}=\frac{\nu _f}{R_{c0}D_{\rho }^\frac{1}{2}}=\left(\frac{\nu _c}{R_{c0}R_b}\right)^{\frac{1}{2}} \end{equation} \tag{ A3 }$$

while the volume filling factor (f_f) equals

$$\begin{equation} f_f=\left(\frac{M_b\rho _c}{NM_c\rho _b}\right)^{-1}. \end{equation} \tag{ A4 }$$

R_c0 is the radius of a 1 M_☉ core. If we assume R_c0 to be invariant, and adopt ν_c and D_ρ to be constants, C is irrelevant to any other parameter. A is only related to the volume filling factor f_f under this assumption, which is related to the thermal interaction between two phase media, requiring further investigation in future work.

We solve Equation (A2) numerically implicitly with first order accuracy in time by discretizing the solution in both space and time. We use an N = 100 discrete mass bin ranging from minimum mass m_μ to maximum mass m_max in logarithmic scale. (See Table 1, showing all chosen parameters.)

$$\begin{equation} x_j=0.05*(1+0.09)^{j-1} \quad \quad (1<j<100). \end{equation} \tag{ A5 }$$

Introducing $\delta _j^i=f_j^{i+1}-f_j^i$, and ignoring the second order term of δ, we obtain a linear function of δ.

$$\begin{eqnarray} \frac{\delta _j^i}{\Delta t}&=&A\left[\sum _{x_k=x_\mu }^{x_k=x_j-x_\mu }\!\! \left(f_k^{i}f_{j-k}^{i}+\delta _k^{i}f_{j-k}^{i}+\delta _{j-k}^{i}f_k^{i}\right) \left(x_k^\frac{1}{3}+x_{j-k}^\frac{1}{3}\right)^2\,\Delta x_k\right. \nonumber \\ &&\left.-\, 2\sum _{x_k=x_\mu }^{x_k=x_{{\rm max}}}\left[f_k^{i}f_{j}^{i}+\delta _k^if_{j}^{i}+\delta _{j}^if_k^{i}\right] \left(x_k^\frac{1}{3}+x_j^\frac{1}{3}\right)^2\,\Delta x_k\right]\nonumber\\ && +C\!\left[\!\sum ^{x_k=x_{\rm BE}}_{x_k=x_j} \left(\!\frac{\delta _k^i+f_k^{i}}{x_k^{\frac{1}{6}}}\!\right)\!\frac{\Delta x_k}{x_j}\,{-}\,\left(\!\frac{x_j-x_\mu }{m_j}\!\right) \!\!\left(\!\frac{\delta _j^i+f_j^{i}}{x_j^{\frac{1}{6}}}\!\right)\!\right]\nonumber \\ \end{eqnarray} \tag{ A6 }$$

which can be reduced to the form of

$$\begin{equation} \sum [\Psi I_{j,k}+\Phi _{j,k}]\,\delta _{k}=B_j. \end{equation} \tag{ A7 }$$

Here, we take I_{j, k} as an unit matrix, where Ψ and Φ are some constant only related to mass distribution in the old time level.

Instead of implementing a thoroughly numerical analysis for the stability properties of the explicit method in Equation (A2), we use the result from this implicit method to justify the robustness of the solution. For reasonable time steps and initial conditions, we are also able to obtain a consistent result from the explicit numerical method. To reduce the numerical time cost, we will switch to the explicit method for further parameter space investigation in the verified region.