Numerical Simulation and Completeness Survey of Bubbles in the Taurus and Perseus Molecular Clouds

Stellar feedback plays a crucial role in the dynamics and energy balance of the interstellar medium (ISM; Zinnecker & Yorke 2007). Feedback associated with protostars injects momentum and energy to the parent molecular cloud, altering the velocity field and density distribution of the cloud (Arce & Goodman 2001, 2002; Arce et al. 2010), contributing to the mass loss of the surrounding dense gas (Fuller & Ladd 2002; Arce & Sargent 2006; Arce et al. 2010; Offner & Arce 2015), affecting star formation efficiency (Frank et al. 2014; Federrath 2015; Kroupa et al. 2018), inducing changes in the chemical composition of the impacted media (Bally et al. 2007), sustaining or generating turbulence (e.g., Fukui et al. 1986; Arce & Goodman 2001; Arce et al. 2011; Nakamura et al. 2011, 2011; Plunkett et al. 2013; Federrath 2015; Feddersen et al. 2018—henceforth Ar11, Fe18, see also Li et al. 2015—henceforth L15), resisting gravitational collapse and even disrupting the surrounding gas to limit the lifetime of their parent molecular cloud (Solomon et al. 1981; Arce & Goodman 2001; Hartmann et al. 2001; Arce & Goodman 2002; Duarte-Cabral et al. 2012; Plunkett et al. 2013).

The primary manifestations of stellar feedback are molecular outflows and bubbles. On the Galactic scale, superbubbles are shown to be globally important, with structures spanning hundreds of parsecs, such as chimneys and large cavities (Heiles 1979; Norman & Ikeuchi 1989). Abundant studies of outflows have been carried out (Kwan & Scoville 1976; Snell et al. 1980; Lada & Harvey 1981; Arce et al. 2010; Nakamura et al. 2011, 2011; Narayanan et al. 2012; Mottram et al. 2017). In local to main star-forming molecular cloud complexes, relatively less attention has been paid to "local" molecular bubbles that are hollowed out by low- to intermediate-mass stars. Despite being expected to be less influential than superbubbles, they are substantial in number and could inject much more energy back into their natal clouds than that from outflows (Ar11, L15, Fe18).

Molecular bubbles are partially or fully enclosed three-dimensional structures whose projections on the sky resemble partial or full rings (Churchwell et al. 2006). Young stellar objects have sufficient stellar winds to entrain and accelerate ambient gas to sculpt spherical or ring-like cavities in their surrounding molecular clouds. Compared to collimated outflows, bubbles affect a larger volume of ambient molecular gas. Compared to supernova remnants, bubbles occur around more stars and persist for a longer period of time (Matzner 2002; Arce et al. 2011). In the Taurus molecular cloud complex (TMC), the kinetic energy and energy injection rate of identified bubbles are larger than those of outflows, implying a substantial input of mechanical energy into the TMC (L15). Due to their complex morphology, bubbles are more difficult to identify (Beaumont et al. 2014) in comparison to the relatively clear characteristics of outflows. In this work, we examine the completeness of bubble identification from previous studies.

The Five College Radio Astronomical Observatory (FCRAO) CO survey (Goldsmith et al. 2008) remains one of the largest molecular spectral-line maps of a continuous cloud complex. Following an empirical and iterative procedure, L15 identified 37 bubbles in the TMC. Meanwhile, Ar11 detected 12 bubbles in the Perseus molecular cloud complex (PMC). The identified Taurus bubbles inject ∼9.2 × 10⁴⁶ erg into the surrounding ISM. The energy injection rate was measured to be ∼6.4 × 10³³ erg s⁻¹, surprisingly comparable to the turbulence dissipation rate. In the PMC, the kinetic energy of identified bubbles is around ∼7.6 × 10⁴⁶ erg with an energy input rate of around ∼8.9 × 10³³ erg s⁻¹, which is similar to the turbulence dissipation rate. Small or slowly expanding bubbles may be missed. However, large or irregular bubbles may also be missed due to confusion. We approach this problem with semiempirical numerical simulations of artificial molecular bubbles based on the bubble model from Cazzolato & Pineault (2005). We generated artificial spherical and partially spherical ¹³CO bubbles and embedded them into the real data cubes of Taurus and Perseus. We then carried out parameter studies including average antenna temperature, number of pixels, and expansion velocity to examine the completeness of the empirical identification procedures.

In Section 2, we provide a basic description of the Taurus and Perseus data sets and the bubble identification routines of previous surveys. The bubble model is presented in Section 3, including a detailed description of our method for simulating bubbles and artificial injection into the data sets. In Section 4, we describe the bubble detection procedures and quantify detectability. The properties of missing bubbles and the implication for stellar feedback in Taurus and Perseus are discussed in Section 5. In Section 6, we discuss and summarize the conclusions which may be drawn from this research.

¹³CO(1−0) is a tracer frequently used to trace molecular hydrogen and study the dynamic morphology of the ISM due to its much smaller opacities than those of ¹²CO(1−0). The Taurus ¹³CO(1−0) and the Perseus ¹³CO(1−0) data sets we used are from the FCRAO CO surveys conducted between 2003 and 2005 (Ridge et al. 2006; Goldsmith et al. 2008; Narayanan et al. 2008; COMPLETE team 2011).

The common empirical procedures used to identify bubble structures in molecular clouds are as follows:

• 1.  
  To search the regions with circular (or arc) structures brighter than the surrounding molecular gas (Ar11, L15, Fe18).
• 2.  
  To visualize the high-velocity features of these regions by plotting a position–velocity (P–V) diagram. If there is an expanding bubble with a redshifted part or a blueshifted part, we can identify the ∪ or ∩-shaped feature on the P–V diagram of each region (Ar11, L15, Fe18).
• 3.  
  To search for candidate sources such as pre-main-sequence stars or protostars associated with the cloud (Ar11, L15, Fe18).
• 4.  
  To compare the bubble structures in ¹³CO(1−0) with an infrared map to determine whether they have similar morphologies (Ar11, L15, Fe18).

The more conditions identified bubbles satisfy, the more likely they are to be real dynamic feedback driven by protostars rather than superpositions of unrelated patterns in the cloud. The radius and thickness of each bubble can be estimated by fitting a Gaussian to the azimuthally averaged profile of the CO integrated intensity map (Ar11, L15).

We simulated spherical and partially spherical ¹³CO(1−0) bubbles based on the model proposed by Cazzolato & Pineault (2005) for the TMC and the PMC. This model relies on three fundamental assumptions:

• 1.  
  The density of the surrounding medium and within the bubble is different but homogeneous.
• 2.  
  The bubble expansion is isotropic.
• 3.  
  All pixels on the bubble are at the same distance from Earth.

There are eight significant parameters in this model, listed in Table 1. Since identified bubbles in Taurus and Perseus often have partial circular or arc-like morphology, we employed completeness of bubble β as another parameter. The detailed bubble simulation calculation is shown in the Appendix A.

Table 1.  Artificial Bubble Parameters^a

Parameters	Definition	Taurus	Perseus
R	Bubble radius	[0.28 pc–1.90 pc]	[0.14 pc–2.79 pc]
ΔR	Bubble thickness	[0.03 pc, 0.44 pc]	[0.1 pc, 0.68 pc]
Vexp	Bubble expansion velocity	[1.0 km s−1–3.3 km s−1]	[1.2 km s−1–6.0 km s−1]
σ	Velocity dispersion	∼1.4 km s−1	∼2.0 km s−1
β	The completeness of bubble	[90°, 360°]	[90°, 360°]
${n}_{0}^{{\prime} }$	Bubble H2 number density	35 cm−3	280 cm−3
	with assumed 13CO abundance of $1.43\times {10}^{-6}$.
	This is estimated based on intensity of
	the high-velocity line wings and used to
	calculate the amount of gas being moved by
	bubbles (see more explanation in
	the Discussion section).
ΔV	Velocity width of a spectrometer	0.266 km s−1	0.066 km s−1
	channel
apc	The bubble physical depth along	Calculations are presented in Appendix A of Cazzolato & Pineault (2005)
	the LOS in parsecs.

Note. ^aParameter name, definition, values for Taurus, and values for Perseus. First five parameters are measured values from L15 and Ar11, while the last four parameters are adopted values from calculation and observation property.

Download table as:  ASCIITypeset image

These artificial bubbles were embedded individually into the Taurus and Perseus ¹³CO(1−0) data cubes at random positions and channels. We compared inserted artificial bubbles with real identified bubbles by morphology and high-velocity features in Figures 1–2. Many of the identified bubbles in Taurus are ambiguous and difficult to identify, as in the right column of Figure 1. The identified bubbles in Perseus tend to have more prominent and identifiable apparent characteristics.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We simulated artificial bubbles with a similar radius, thickness, β, and expansion velocity for each cloud, as shown in the left columns of Figures 1–2. Circular structures brighter than the surrounding molecular gas may be identified in the intensity and channel maps. Meanwhile, the P–V diagram illustrates the ∪- or ∩-shaped features. We found that in general artificial bubbles are somewhat easier to detect than real ones because they are not really embedded within the surrounding environment but are added on to the original images, and their morphology is necessarily independent of the surrounding environment. This means that our detectability function actually presents an upper limit for the bubble detectability function with the same parameters. We found it challenging to simulate realistic artificial bubbles in part to the ambiguous morphology of identified bubbles.

Identifying the circular structure in the channel maps is the first step in bubble identification. Therefore, we searched for the artificial bubbles by going through the channel maps. When we detected a circular structure, we compared the position of detected structure with the position of the inserted artificial bubble. If they are consistent with each other, we consider the artificial bubble detectable. When we cannot detect any circular structure through channel maps, we consider the artificial bubble undetectable.

We generated 500 artificial bubbles with different parameters for each cloud. Each bubble was embedded into Taurus and Perseus separately in random positions and channels one at a time. We undertook a blind search for each inserted simulated bubble to determine the detectability function. The artificial bubble identification routine largely follows the previous bubble manual identification procedures in Section 2.

4.1. Description of Variable Parameters

4.2. Analysis Method

The previous bubble identification procedures are subjective and depending on the apparent detection threshold resulting from the bubble identification procedures in L15, Ar11, and Fe18. We quantified the bubble detection results of previous bubble identification procedures with the detectability function. Since the bubble detection results correspond to a Bernoulli distribution, which is one of the exponential family distributions, and the experimental bubble parameters can be considered as independent to each other, we constructed the bubble detectability function using generalized linear models (GLMs; Nelder & Wedderburn 1972). GLMs are widely used regression models for dependent variables which follow an exponential family distribution, such as Gaussian distribution, Poisson distribution, and chi-squared distribution. Derivation of the detectability function is discussed in Appendix B. The detectability function is in the form of

$$\begin{eqnarray}&&\begin{array}{l}h({\mu }_{{T}_{a}},{\mu }_{N},{\mu }_{{V}_{e}})\,=\,\displaystyle \frac{1}{1+{e}^{-({\alpha }_{0}+{\alpha }_{1}^{* }{\mu }_{{T}_{a}}+{\alpha }_{2}^{* }{\mu }_{N}+{\alpha }_{3}^{* }{\mu }_{{V}_{e}})}},\end{array}\end{eqnarray} \tag{ 1 }$$

where α₀ to α₃ are constants that need to be fitted, ${\mu }_{{T}_{a}}$, μ_N, and ${\mu }_{{V}_{e}}$ are the bubble average antenna temperature, number of pixels, and expansion velocity divided by their maximum values (see Appendix B). The maximum value of T_a, N, and V_e are obtained from real observed bubbles.

4.3. Parametric Detectability

We performed 500 trials in both the TMC and the PMC to obtain the data set of N, T_a, and V_e, and the detection result (0 or 1). N, T_a, and V_e are then divided by their maximum values. The maximum values are obtained from real observed bubbles. We performed tenfold cross-validation on the experimental data sets to fit the constants, while testing how well the detectability function performs on the bubble identification and detectability prediction in each cloud. In 10-fold cross-validation, the data sets are randomly partitioned into 10 equal size subsets. One single subset is used as the validation data for testing the model, while the remaining nine subsets are used as training data to fit the model. We used the maximum likelihood estimator (MLE; see Appendix B) to get an estimate for each parameter from training sets and use the remaining subset for validation. We repeated this procedure 10 times. Our estimate for the parameters is the average over the 10 training runs, while the error is the average over the 10 validation runs. Tenfold cross-validation is frequently used when evaluating performance of models on multiclass data (James et al. 2013).

4.3.1. Results in Taurus

We performed tenfold cross-validation on the Taurus experimental data set with random segment selection on all the data for each fold to fit the detectability function. We adopt the average value of 10 times the fitting results for α₀ to α₃, where α₀ = −6.0, α₁ = 112.0, α₂ = 42.9, and α₃ = 3.3. The uncertainties are derived from the standard deviation of 10 times the fitting results for α₀ to α₃, which are 0.3, 8.3, 3.9, and 0.5, respectively. The detectability function in Taurus is described as

$$\begin{eqnarray}&&{h}^{T}({\mu }_{{T}_{a}},{\mu }_{N},{\mu }_{{V}_{e}})=\displaystyle \frac{1}{(1+{e}^{6.0-112.0* {\mu }_{{T}_{a}}-42.9* {\mu }_{N}-3.3* {\mu }_{{V}_{e}}})},\end{eqnarray} \tag{ 2 }$$

where ${\mu }_{{T}_{a}}$, μ_N, and ${\mu }_{{V}_{e}}$ are scaled experimental bubble parameters—average antenna temperature, number of pixels, and expansion velocity—in Taurus.

Meanwhile, we estimated the training error from the training set and generalization error from the testing set to analyze how well h^T performs on bubble identification and prediction in Taurus. The average training error is about 0.14 which is the probability that h^T misclassifies samples in training sets. The average generalization error for h^T is about 0.12 which is the probability if we draw a new set of bubble parameters and bubble detection results, h^T misclassifies it. In Taurus bubble experiments, there are 330 identified experimental bubbles, which are 66% of total experiments. As long as the detectability function can fit and predict the bubble detection result with correctness larger than 66%, we can consider that the detectability function can be used for bubble identification. According to the average training error and generalization error, the correctness of h^T to fit and predict the bubble detection result is 86% and 88%, respectively, which indicates that h^T can well fit and predict the bubble detection result in Taurus.

4.3.2. Results in Perseus

We performed tenfold cross-validation on the Perseus experimental data set with random segment selection on all the data for each fold to fit the detectability function. We adopt the average value of 10 times fitting results for α₀ to α₃, where α₀ = −10.2, α₁ = 121.3, α₂ = 6.3, and α₃ = 2.9. The uncertainty is derived from the standard deviation of 10 times the fitting results for α₀ to α₃, which are 1.5, 23.6, 1.1, and 0.2, respectively. The detectability function in Perseus is described as

$$\begin{eqnarray}&&{h}^{P}({\mu }_{{T}_{a}},{\mu }_{N},{\mu }_{{V}_{e}})=\displaystyle \frac{1}{(1+{e}^{10.2-121.3* {\mu }_{{T}_{a}}-6.3* {\mu }_{N}-2.9* {\mu }_{{V}_{e}}})},\end{eqnarray} \tag{ 3 }$$

where ${\mu }_{{T}_{a}}$, μ_N, and ${\mu }_{{V}_{e}}$ are scaled experimental bubble parameters—average antenna temperature, number of pixels, and expansion velocity—in Perseus.

Meanwhile, we estimated the training error from the training set and generalization error from the testing set to analyze how well h^P performs on bubble identification and prediction in Perseus. The average training error is about 0.06 which is the probability that h^P misclassifies samples in training sets. The average generalization error for h^P is about 0.08 which is the probability if we draw a new set of bubble parameters and bubble detection result, h^P would misclassify it. In Perseus bubble experiments, there are 321 identified experimental bubbles, which are 64.2% of total experiments. As long as the detectability function can fit and predict the bubble detection result with correctness larger than 64.2%, we can consider that the detectability function can be used for bubble identification. According to the average training error and generalization error, the correctness of h^P to fit and predict bubble detection result is 94% and 92%, respectively, which indicate that h^P can well fit and predict the bubble detection result in Perseus.

4.4. Comparison

The detectability functions for bubble average antenna temperature, number of pixels, and expansion velocity for each cloud are illustrated in Figure 3. In Taurus, the change from yellow to deep blue is more gradual than in Perseus, which is caused by the larger training error in Taurus. We found that bubble detectability in Taurus and Perseus both strongly depend on average antenna temperature and number of pixels. The brightness of the bubble understandably dominates the detectability of the bubbles. In the case of a large, relatively regularly shaped bubble, the detectability roughly scaled with the number of pixels. Weak dependence is seen to occur for expansion velocity. However, we still find that bubbles with larger expansion velocity are easier to detect. In general, bubbles in Perseus are easier to detect than Taurus with the same N, T_a, and V_e. According to the behaviors of the detectability function, bubbles with low average antenna temperature, small spatial size, and slow expansion velocity are more likely to be missed during manual identification. We quantify the kinetic energy of the missing bubble in the following section.

Zoom In Zoom Out Reset image size

Bubble kinetic energy is associated with the input of mechanical energy into the parent molecular cloud, which is the crucial parameter to evaluate the impact of bubbles on the surrounding ISM. We estimated the completeness of previous bubble surveys by comparing the kinetic energy of missing bubbles with identified bubbles in Taurus and Perseus. The identified bubbles here and afterward represent the real identified bubbles in previous bubble surveys.

In order to estimate the kinetic energy distribution of all bubbles, including identified and missing bubbles in each molecular cloud, we assumed that the kinetic energy distribution of identified bubbles results from the kinetic energy distribution of all bubbles modified by the bubbles' detectability. Therefore, the number of the missing bubble can be estimated by the difference between the kinetic energy distribution of identified bubbles and the total number of bubbles predicted by our model.

The relation between the kinetic energy distribution of identified bubbles and all bubbles in the molecular cloud can be expressed as

$$\begin{eqnarray}&&\int {N}_{D}\cdot p({E}_{k}){{dE}}_{k}=\int {N}_{A}\cdot {p}^{\mathrm{tot}}({E}_{k})\cdot h({E}_{k}){{dE}}_{k},\end{eqnarray} \tag{ 4 }$$

where N_D is the number of real identified bubbles, N_A is the number of all bubbles in a given cloud, $p({E}_{k})$ is the probability density function (PDF) of identified bubbles with respect to kinetic energy, ${p}^{\mathrm{tot}}({E}_{k})$ is the PDF of all bubbles with respect to kinetic energy within the cloud, and $h({E}_{k})$ denotes the detectability function with respect to kinetic energy.

The number of missing bubbles ${N}^{\mathrm{miss}}$ for kinetic energy is described as

$$\begin{eqnarray}\begin{array}{rcl}{N}^{\mathrm{miss}} & = & {\displaystyle \int }_{{E}_{k}^{\min }}^{{E}_{k}^{\max }}[{N}_{A}\cdot {p}^{\mathrm{tot}}({E}_{k})-{N}_{D}\cdot p({E}_{k})]{{dE}}_{k},\\ & = & {N}_{D}\cdot {\displaystyle \int }_{{E}_{k}^{\min }}^{{E}_{k}^{\max }}[p({E}_{k})/h({E}_{k})-p({E}_{k})]{{dE}}_{k}.\end{array}\end{eqnarray} \tag{ 5 }$$

The kinetic energy of missing bubbles ${E}_{k}^{\mathrm{miss}}$ is written by

$$\begin{eqnarray}\begin{array}{rcl}{E}_{k}^{\mathrm{miss}} & = & {\displaystyle \int }_{{E}_{k}^{\min }}^{{E}_{k}^{\max }}{E}_{k}[{N}_{A}\cdot {p}^{\mathrm{tot}}({E}_{k})-{N}_{D}\cdot p({E}_{k})]{{dE}}_{k},\\ & = & {N}_{D}\cdot {\displaystyle \int }_{{E}_{k}^{\min }}^{{E}_{k}^{\max }}{E}_{k}[p({E}_{k})/h({E}_{k})-p({E}_{k})]{{dE}}_{k}.\end{array}\end{eqnarray} \tag{ 6 }$$

We constructed the kinetic energy of each experimental bubble based on their number of pixels, average antenna temperature, and expansion velocity. The detectability of each experimental bubble was evaluated by the detectability function of Equations (2) and (3). By fitting E_k with detectability using MLE, we estimate parameters of the kinetic energy detectability function. The parameters of PDF of the kinetic energy of identified bubbles is estimated from the MLE. A detailed description of the fitting method is illustrated in Appendix C.

5.1. The Kinetic Energy of Missing Bubbles in Taurus

There are 37 identified bubbles N_D in Taurus. Their total kinetic energy is about 9.2 × 10⁴⁶ erg. The maximum and minimum kinetic energies of identified bubbles are ${E}_{k}^{\max T}\,=4.18\times {10}^{46}$ erg and ${E}_{k}^{\min T}=0.2\times {10}^{44}$ erg, respectively. The detectability of each experimental bubble in Taurus is evaluated from Equation (2). We plotted E_k versus detectability for both experimental bubbles and identified bubbles in the right panel of Figure 4. Our detectability function apparently overestimates the detectability of the bubble with low kinetic energy. However, such overestimation does not affect energy calculation related to feedback, as these bubbles add up to negligible energies compared to bright ones. One possibility is that there is a large number of small, low-energy, undetectable bubbles that might contribute significant energy to the clouds. In order for this to be true, there would need to be at least 4600 such bubbles in the TMC and 379 in the PMC. We find this to be unlikely. By fitting E_k with detectability using MLE, we obtain parameters of the kinetic energy detectability function by

$$\begin{eqnarray}&&h({E}_{k})={\tanh }(925.5\times {E}_{k}^{u}),\end{eqnarray} \tag{ 7 }$$

where ${E}_{k}^{u}=\tfrac{{E}_{k}}{{E}_{k}^{\max T}}$ refers to uniformed kinetic energy. We chose the tanh function because it changes from 0 to 1 gradually when kinetic energy is larger than 0, which is consistent with the detectability of all experimental bubbles from small kinetic energy to large kinetic energy. The distribution of detectability of all experimental bubbles is not a Bernoulli distribution. Therefore, we cannot use the logistic function to fit the kinetic energy detectability directly.

Zoom In Zoom Out Reset image size

The bubble detectability increases as the kinetic energy increases, which makes the bubble with low kinetic energy hard to identify. The number of identified bubbles with low kinetic energy should be small. Although the bubble with high kinetic energy is easier to identify, it is unlikely that low-mass star formation in Taurus and Perseus would generate a lot of high kinetic energy bubbles. Therefore, the number of identified bubbles with high kinetic energy should be small as well. Lognormal distribution satisfies the above conditions and can fit the real identified bubbles well. The probability distribution of kinetic energy for identified bubbles $p({E}_{k})$ can be well fitted in truncated lognormal distribution using MLE (see Appendix C) by

$$\begin{eqnarray}&&p({E}_{k})=\displaystyle \frac{-\sqrt{2}{e}^{-\tfrac{1}{2}\tfrac{1}{{\sigma }^{2}}{\left(\mathrm{ln}\left(\tfrac{{E}_{k}}{m}\right)\right)}^{2}}}{\sqrt{\pi \sigma }\left(\mathrm{erf}\left({a}_{3}\right)-\mathrm{erf}\left({a}_{8}\right)\right){E}_{k}},\end{eqnarray} \tag{ 8 }$$

where σ is the fitting shape parameter which is 1.7 ± 0.80 in logarithm, erf is the error function, we defined $\mu =\mathrm{ln}(m)$ to be the logarithmic fitting scale parameter which is 102.8 ± 0.27, a₃ and a₈ are compact notations which are given by

$$\begin{eqnarray*}&&{a}_{3}=\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}\left(\mathrm{ln}\left({E}_{k}^{\min T}\right)-\mathrm{ln}\left(m\right)\right)}{\sigma },\end{eqnarray*}$$

$$\begin{eqnarray*}&&{a}_{8}=\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}\left(\mathrm{ln}\left({E}_{k}^{\max T}\right)-\mathrm{ln}\left(m\right)\right)}{\sigma }.\end{eqnarray*}$$

The fitting uncertainty was derived from 1000 sets of Monte Carlo experiments, each of size 37 (the number of identified bubbles in Taurus) with the kinetic energies ranging from ${E}_{k}^{\min T}$ to ${E}_{k}^{\min T}$. For each Monte Carlo experiment, the kinetic energy is randomly generated from the truncated lognormal distribution with μ = 102.8, σ = 1.7. We applied the same fitting algorithm to those 37 random samples to estimate μ and σ and repeated the Monte Carlo experiment 1000 times to get sets of the estimated value of μ and σ. We took the standard deviation of the estimated 1000 sets of μ and σ to be the fitting uncertainty. Since there are only 37 identified bubbles in Taurus, we did not plot the PDF of E_k of identified bubbles, instead we plotted the cumulative distribution function (CDF) of E_k. The CDF of E_k for identified bubbles and CDF of E_k of fitting distribution for Taurus are illustrated in the left panel of Figure 4.

By combining Equations (7) and (8) with Equation (5), the number of missing bubbles in Taurus can be derived, which ranges from 2 to 3. This is the probable number of bubbles not identified in surveys of L15 in Taurus. The total kinetic energy of the missing bubbles is estimated by combining Equations (7) and (8) with Equation (6). However, the kinetic energy of most missing bubbles aggregates at less than 1 × 10⁴⁴ erg as illustrated in the green line of the right panel in Figure 6, which leads to the total kinetic energy of missing bubbles is ranging from 7.2 × 10⁴³ to 8.6 × 10⁴³ erg. This corresponds to about 0.01% of the kinetic energy of identified bubbles. Therefore, the Taurus bubble survey can be considered to be energetically complete.

5.2. Kinetic Energy of Missing Bubbles in Perseus

There are 12 identified bubbles N_D in Perseus. Their total kinetic energy is about 7.58 × 10⁴⁶ erg. The kinetic energy of identified bubbles ranges from 2 × 10⁴⁴ to 1.88 × 10⁴⁶ erg. We adopt the minimum and maximum kinetic energy of identified bubbles in Taurus (which has a wider range than in Perseus) to be the range used in estimating the missing kinetic energy in Perseus.

The detectability of each experimental bubble in Perseus is evaluated from Equation (3). We plotted E_k versus detectability for both experimental bubbles and identified bubbles in the right panel of Figure 5. Similar to Taurus, we obtain parameters of the kinetic energy detectability function by fitting E_k with detectability using MLE, which can be written as

$$\begin{eqnarray}&&h({E}_{k})={\tanh }(79.34\times {E}_{k}^{u}),\end{eqnarray} \tag{ 9 }$$

where ${E}_{k}^{u}=\tfrac{{E}_{k}}{{E}_{k}^{\max T}}$ refers to uniformed kinetic energy.

Zoom In Zoom Out Reset image size

Similarly, the probability distribution of kinetic energy for identified bubbles in Perseus can also be well fitted in truncated lognormal distribution using MLE shown in Equation (8). The fitting shape parameter is 2.10 ± 0.43 in logarithm. The logarithmic fitting scale parameter is 119.92 ± 0.22. The fitting uncertainty was derived from 1000 sets of Monte Carlo experiments, each of size 12 (the number of identified bubbles in Perseus) with the kinetic energies ranging from ${E}_{k}^{\min P}$ to ${E}_{k}^{\min P}$. For each Monte Carlo experiment, the kinetic energy is randomly generated from the truncated lognormal distribution with μ = 119.92 and σ = 2.10. We applied the same fitting algorithm to those 12 random samples to estimate μ and σ then repeated the Monte Carlo experiment 1000 times to get sets of the estimated value of μ and σ. The fitting uncertainty is the standard deviation of the estimated 1000 sets of μ and σ. There are only 12 identified bubbles in Perseus. Similar to Taurus, we plotted up the CDF of E_k for identified bubbles and CDF of E_k of the fitting distribution for Perseus, which are illustrated in the left panel of Figure 5.

By combining Equations (8) and (9) with Equation (5), the number of missing bubbles in Perseus can be derived which is ranging from 3 to 14. This is the probable number of bubbles not identified in the Perseus survey by Ar11. The total kinetic energy of the missing bubbles is estimated by combining Equations (9) and (8) with Equation (6). However, the kinetic energy of most missing bubbles aggregates at less than 2 × 10⁴⁴ erg as illustrated in the green line of the left panel of Figure 6, which leads to the total kinetic energy of missing bubbles is ranging from 1.4 × 10⁴⁴ to 8.0 × 10⁴⁴ erg. This corresponds to about 1% of the kinetic energy of identified bubbles. Therefore, the Perseus bubble survey can be considered to be energetically complete.

Zoom In Zoom Out Reset image size

Stellar feedback has a significant impact on the surrounding gas. The total kinetic energy of bubbles detected in Taurus and Perseus are ∼3.9 × 10⁴⁶ erg and ∼7.6 × 10⁴⁶ erg, respectively. The gravitational potential energy and turbulence energy (E_G, E_tur) of Taurus and Perseus are about 1.5 × 10⁴⁸ erg, 3.2 × 10⁴⁷ erg, and 6 × 10⁴⁷, 1.6 × 10⁴⁷ erg, respectively. The energy contained in bubbles are orders of magnitude smaller than those of either gravity or turbulence, which are expected to be equal in a virialized supersonic cloud. The maximum energy injection (into the natal cloud) rate of molecular bubbles can be estimated from their momentum. They are ∼6.4 × 10³³ erg s⁻¹ and ∼2 × 10³³ erg s⁻¹ for bubbles in Taurus and Perseus, respectively. Following the methods laid out in Ar11 and L15, we estimated the turbulence dissipation rate of Taurus and Perseus to be ∼3.1 × 10³³ erg s⁻¹ and ∼1 × 10³³ erg s⁻¹, respectively. The observed injection rates are thus similar to or even slightly larger than the dissipation rates.

Typically, when discussing molecular clouds we refer to total proton volume densities of 10³ cm⁻³ for regions with sufficient volume density to produce ¹³CO emission. This is typical for regions with A_v ≥ 1. However, the volumes of gas that are affected by bubbles are large and not necessarily located in the dense regions. In fact, they almost necessarily include large volumes of very diffuse gas. The gas that is dynamically affected ("moved") by bubble expansion and can be identified as such has to exhibit a velocity offset (Figure 7). Such "high" velocity gas has much lower density. Due to this, we chose to use our estimates of the average H₂ volume densities within Taurus and Perseus (35 cm⁻³ and 280 cm⁻³ respectively) when estimating the energy injected into the cloud by each bubble. As a comparison, if the main cloud density was assumed to be 3.5 × 10³ cm⁻³, the energy injected into the cloud by star formation feedback in Taurus and Perseus would be ∼3.9 × 10⁴⁸ erg and ∼9.1 × 10⁴⁷ erg, respectively, which could disperse the whole cloud. Similarly, if we use such a high volume density for the undisturbed gas in our simulations, then the resulting bubbles are much brighter than any that are actually observed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We simulated bubbles for the Taurus and Perseus molecular cloud. We randomly inserted artificial bubbles into the Taurus and Perseus ${}^{13}\mathrm{CO}$ data cubes from the 13.7 m FCRAO telescope. With changing average antenna temperature, the number of pixels, and expansion velocity of artificial bubbles, we parameterized bubble detections according to the bubble detectability functions to evaluate how the detectability is affected by the parameters for two molecular clouds. According to the properties of identified bubbles in Taurus and Perseus and the detectability functions, we estimated the energetic completeness of previously stellar feedback studies. Our conclusions regarding the completeness of bubble identification in Taurus and Perseus and their properties are as follows:

• 1.  
  The detectability of bubbles can be described as a logistic function, which is derived from GLMs. The distribution of bubble detection results is a Bernoulli distribution. Therefore we adopted GLMs to fit the bubble detectability functions. In the TMC, the detectability function can be described as ${h}^{T}({\mu }_{{T}_{a}},{\mu }_{N},{\mu }_{{V}_{e}})\,=\tfrac{1}{(1+{e}^{6.0-112.0* {\mu }_{{T}_{a}}-42.9* {\mu }_{N}-3.3* {\mu }_{{V}_{e}}})}$, while in the PMC as ${h}^{P}({\mu }_{{T}_{a}},{\mu }_{N},{\mu }_{{V}_{e}})=\tfrac{1}{(1+{e}^{10.2-121.3* {\mu }_{{T}_{a}}-6.3* {\mu }_{N}-2.9* {\mu }_{{V}_{e}}})}$. We then fitted the bubble kinetic energy distributions in Taurus and Perseus with truncated lognormal distribution.
• 2.  
  The number of missing bubbles in Taurus is less than 8% of the number of identified bubbles. The number of missing bubbles in Perseus is about 25%–125% of the number of identified bubbles.
• 3.  
  We used bubble kinetic energy distributions and the bubble detectability functions to estimate the total kinetic energy of missing bubbles, which suggests that although the numbers of missing bubbles are large, their kinetic energies are relatively small (usually less than 2 × 10⁴⁴ erg). The total kinetic energies of missing bubbles in Taurus and Perseus during manual identification range from 7.2 × 10⁴³ erg to 8.6 × 10⁴³ erg and 1.4 × 10⁴⁴ erg to 8.0 × 10⁴⁴ erg, respectively. Such potential incompleteness only accounts for ∼0.2% and ∼1% of the total kinetic energy of identified bubbles in Taurus and Perseus, respectively.
• 4.  
  The empirical surveys in L15 and Ar11 for identifying bubble structures in Taurus and Perseus can be considered as energetically complete.
• 5.  
  The total energy of bubbles in a cloud is orders of magnitude smaller than those of either turbulence or gravity. The bubbles cannot generate the observed turbulence or disperse the cloud. The observed energy injection rate from bubbles, now considered complete, is similar to the turbulence dissipation rate. We conclude that, even in low-mass star-forming regions, the feedback from star formation is sufficient to sustain turbulence at ranges from ∼0.1 pc to ∼2.8 pc scales.

This work is supported by the National Natural Science Foundation of China grant No. 11988101, No. 11725313, No. 11721303, the International Partnership Program of Chinese Academy of Sciences grant No. 114A11KYSB20160008, the National Science Foundation of China No. 11721303, and the National Key R&D Program of China No. 2016YFA0400702.

The ¹³CO column densities may be described by

$$\begin{eqnarray}&&\ {N}_{{}^{13}\mathrm{CO}}=1.18\times {10}^{15}{T}_{b}{\rm{\Delta }}V,\end{eqnarray} \tag{ 10 }$$

where T_b is the brightness temperature of each pixel, ${N}_{{}^{13}\mathrm{CO}}$ is the column density, and ΔV represents the velocity width of a spectrometer channel which is listed in Table 1. This relation is derived under several assumptions. We assumed that the excitation temperature T_ex of ¹³CO is 25 K, the background temperature is 2.7 K, the ¹³CO emission from the bubble is generally optically thin ($\tau {(}^{13}\mathrm{CO})\ll 1$) and in local thermal equilibrium. ${N}_{{}^{13}\mathrm{CO}}$ can also be expressed as

$$\begin{eqnarray}&&\ {N}_{{}^{13}\mathrm{CO}}=3.0857\times {10}^{18}{n}_{s}{a}_{\mathrm{pc}}\,{\mathrm{cm}}^{-2},\end{eqnarray} \tag{ 11 }$$

where a_pc describes the bubble physical depth along the line of sight (LOS) in parsecs. By combining Equation (10) with Equation (11), we obtain the T_b for each pixel in the artificial bubbles with known a_pc, ΔV, and n_s, which are given by

$$\begin{eqnarray}&&{T}_{b}=2.615* {10}^{3}\displaystyle \frac{{n}_{s}{a}_{\mathrm{pc}}}{{\rm{\Delta }}V}.\end{eqnarray} \tag{ 12 }$$

We assume that the ¹³CO abundance relative to H₂ is 1.43 × 10⁻⁶. The densities of ¹³CO affected by bubble expansion n₀ are 5 × 10⁻⁵ cm⁻³ and 4 × 10⁻⁴ cm⁻³ for Taurus and Perseus, respectively. If all of the ¹³CO within the radius R + ΔR is distributed within the bubble, we can get n_s/n₀ = 1-ξ⁻³, where ξ = (R + ΔR)/R.

Figure 8 illustrates how a_pc is calculated. It shows a ¹³CO bubble sampled by five velocity channels in the left hemisphere as seen by an observer located to the far left. The starting and ending angles of the second velocity channel are θ_i and θ_i + 1, respectively. The angular dimension of each velocity channel is given by $\cos ({\theta }_{i+1})=\cos ({\theta }_{i})+{\rm{\Delta }}V/{V}_{\exp }$. Thereby, the angle of any LOS within the velocity channel can be determined. The physical depth a_pc of each LOS can be calculated from the angle of this LOS, bubble radius, and thickness. The detailed calculations are presented in Appendix A of Cazzolato & Pineault (2005). By estimating the physical depth a_pc of each LOS for a given expansion velocity, and knowing n_s of the bubble, we can obtain antenna temperature for each LOS using Equation (12).

According to the definition, the experimental bubble detection result, y, is either 0 or 1 which corresponds to a Bernoulli distribution. The probability mass function p of the detection results y is given by

$$\begin{eqnarray*}\begin{array}{rcl}p(y=1;\phi ) & = & \phi ,\\ p(y=0;\phi ) & = & 1-\phi ,\end{array}\end{eqnarray*}$$

where $p(y=1;\phi )$ is the probability of detecting a bubble and $p(y=0;\phi )$ is the probability of not detecting a bubble. p can be also expressed as (Nelder & Wedderburn 1972)

$$\begin{eqnarray*}\begin{array}{rcl}p(y;\phi ) & = & {\phi }^{y}{\left(1-\phi \right)}^{1-y}\\ & = & \exp (\mathrm{log}({\phi }^{y}{\left(1-\phi \right)}^{1-y}))\\ & = & \exp (y\mathrm{log}(\phi )+(1-y)\mathrm{log}(1-\phi ))\\ & = & \exp \left(\mathrm{log}\left(\displaystyle \frac{\phi }{1-\phi }\right)y+\mathrm{log}(1-\phi )\right).\end{array}\end{eqnarray*}$$

The exponential family is written as

$$\begin{eqnarray*}&&f(y;\eta )=b(y)\exp ({\eta }^{T}T(y)-a(\eta )),\end{eqnarray*}$$

where η is the natural parameter of the distribution; T(y) is the sufficient statistic; and a(η) is the log partition function. When $\eta =\mathrm{log}\left(\tfrac{\phi }{1-\phi }\right)$, $a(\eta )=-\mathrm{log}(1-\phi )$, T(y) = y, and b(y) = 1, we can derive the p(y; ϕ) from f(y; η). This indicates that the probability mass function of the experimental bubble detection result is one of exponential family. The experimental bubble parameters: N, T_a and V_e can be considered as independent to each other, so that we can combine them linearly after scaling. The scaling rule is as follows:

$$\begin{eqnarray}&&{\mu }_{\mathrm{scale}}=\displaystyle \frac{\mu }{{\mu }_{\max }}\end{eqnarray} \tag{ 13 }$$

where μ represents the bubble parameter, μ_min and μ_max refer to the minimum and the maximum value of parameter.

Thus, η can be constructed as η = α^T X. Here, X denotes a three-dimensional vector containing N, T_a, and V_e. α denotes a four-dimensional coefficient vector containing an intercept α₀. Our goal is to predict the expected value of the probability of detecting the bubble y with given X, which means we would like the detectability function h(X) output to satisfy $h(X)=E[y| X]$. In our formulation of the Bernoulli distribution as an exponential family distribution, we had $\eta ={log}?\left(\tfrac{\phi }{1-\phi }\right)$ which can be written as $\phi =\tfrac{1}{1+{e}^{-\eta }}$ and $E[y| X;\alpha ]=\phi$. So that the bubble detectability function h(X) can be expressed as

$$\begin{eqnarray}&&h(X)=E[y| X;\alpha ]\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&=\,\phi \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{1}{1+{e}^{-{\alpha }^{T}X}}.\end{eqnarray} \tag{ 16 }$$

h(X) is the logistic function. With given experimental bubble parameter and bubble detection result data sets, we can estimate the α^T for h(X) based on the maximizing the likelihood function $l(\theta )={\sum }_{i=1}^{\mathrm{rn}}{y}^{(i)}\mathrm{log}(h({X}^{(i)}))+(1-{y}^{(i)})\mathrm{log}(1-h({X}^{(i)})$, where rn is the row number of the data set.

We adopt the MLE to estimate the parameters of PDF of kinetic energy for both Taurus and Perseus using the truncated lognormal distribution. The derivation for the MLE of the truncated lognormal distribution largely follows Zaninetti (2017). Consider a sample ${ \mathcal X }={x}_{1},{x}_{2},...,{x}_{n}$. The maximum and minimum value of sample can be expressed as x_l and x_u, respectively, which are given by

$$\begin{eqnarray}&&{x}_{l}=\max ({x}_{1},{x}_{2},...,{x}_{n}),\qquad {x}_{u}=\min ({x}_{1},{x}_{2},...,{x}_{n}).\end{eqnarray} \tag{ 17 }$$

The PDF can be expressed as

$$\begin{eqnarray}&&\mathrm{PDF}(x;m,\sigma ,{x}_{l},{x}_{u})=\displaystyle \frac{-\sqrt{2}{e}^{-\tfrac{1}{2}\tfrac{1}{{\sigma }^{2}}{\left(\mathrm{ln}\left(\tfrac{x}{m}\right)\right)}^{2}}}{\sqrt{\pi }\sigma \left(\mathrm{erf}\left({a}_{3}\right)-\mathrm{erf}\left({a}_{8}\right)\right)x},\end{eqnarray} \tag{ 18 }$$

where m is the scale parameter, σ is the shape parameter, a₃ and a₈ are compact notations which are given by

$$\begin{eqnarray*}&&{a}_{3}=\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}\left(\mathrm{ln}\left({x}_{l}\right)-\mathrm{ln}\left(m\right)\right)}{\sigma },\end{eqnarray*}$$

$$\begin{eqnarray*}&&{a}_{8}=\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}\left(\mathrm{ln}\left({x}_{u}\right)-\mathrm{ln}\left(m\right)\right)}{\sigma }.\end{eqnarray*}$$

The CDF can be expressed as

$$\begin{eqnarray}&&\mathrm{CDF}(x;m,\sigma ,{x}_{l},{x}_{u})=\displaystyle \frac{-\mathrm{erf}\left(\tfrac{1}{2}\,\tfrac{\sqrt{2}}{\sigma }\mathrm{ln}\left(\tfrac{x}{m}\right)\right)+\mathrm{erf}\left({a}_{3}\right)}{\mathrm{erf}\left({a}_{3}\right)-\mathrm{erf}\left({a}_{8}\right)}.\end{eqnarray} \tag{ 19 }$$

The MLE is obtained by maximizing

$$\begin{eqnarray}&&{\rm{\Lambda }}=\displaystyle \sum _{i}^{n}\mathrm{ln}(\mathrm{PDF}({x}_{i};m,\sigma ,{x}_{l},{x}_{u})).\end{eqnarray} \tag{ 20 }$$

The two derivatives $\tfrac{\partial {\rm{\Lambda }}}{\partial m}=0$ and $\tfrac{\partial {\rm{\Lambda }}}{\partial \sigma }=0$ generate two nonlinear equations in m and σ which can be solved numerically,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\rm{\Lambda }}}{\partial m}=\left(\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{l})-\mathrm{ln}(m))}{\sigma }\right)\right.\\ -\left.\,\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{u})-\mathrm{ln}(m))}{\sigma }\right)\right)\\ \times \,\left(n\sqrt{2}\sigma \,{e}^{-\tfrac{1}{2}\tfrac{{\left(\mathrm{ln}({x}_{l})-\mathrm{ln}(m)\right)}^{2}}{{\sigma }^{2}}}-n\sqrt{2}\sigma \,{e}^{-\tfrac{1}{2}\tfrac{{\left(\mathrm{ln}({x}_{u})-\mathrm{ln}(m)\right)}^{2}}{{\sigma }^{2}}}\right.\\ -\,\sqrt{\pi }\left(\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{l})-\mathrm{ln}(m))}{\sigma }\right)\right.\\ -\,\left.\left.\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{u})-\mathrm{ln}(m))}{\sigma }\right)\right)(n\mathrm{ln}(m)-\displaystyle \sum _{i=1}^{n}\mathrm{ln}({x}_{i}))\right)=0,\end{array}\end{eqnarray} \tag{ 21 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Lambda }}}{\partial \sigma }=\displaystyle \frac{N}{D}=0,\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}N=\mathrm{ln}\left({x}_{u}\right)\sqrt{2}{e}^{-\tfrac{1}{2}\tfrac{{\left(\mathrm{ln}({x}_{u})-\mathrm{ln}(m)\right)}^{2}}{{\sigma }^{2}}}n\sigma -\mathrm{ln}({x}_{l})\sqrt{2}{e}^{-\tfrac{1}{2}\tfrac{{\left(\mathrm{ln}({x}_{l})-\mathrm{ln}(m)\right)}^{2}}{{\sigma }^{2}}}n\sigma \\ +\,\sqrt{2}{e}^{-1/2\tfrac{{\left(\mathrm{ln}({x}_{l})-\mathrm{ln}(m)\right)}^{2}}{{\sigma }^{2}}}\mathrm{ln}(m)n\sigma -\sqrt{2}{e}^{-\tfrac{1}{2}\tfrac{{\left(\mathrm{ln}({x}_{u})-\mathrm{ln}(m)\right)}^{2}}{{\sigma }^{2}}}\mathrm{ln}(m)n\sigma \\ +\,n{\left(\mathrm{ln}(m)\right)}^{2}\sqrt{\pi }\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{u})-\mathrm{ln}(m))}{\sigma }\right)\\ -\,n{\sigma }^{2}\sqrt{\pi }\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{u})-\mathrm{ln}(m))}{\sigma }\right)\\ -\,n{\left(\mathrm{ln}(m)\right)}^{2}\sqrt{\pi }\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{l})-\mathrm{ln}(m))}{\sigma }\right)\\ +\,n{\sigma }^{2}\sqrt{\pi }\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{l})-\mathrm{ln}(m))}{\sigma }\right)\\ +\,\displaystyle \sum _{i=1}^{n}\mathrm{ln}({x}_{i})(\mathrm{ln}({x}_{i})-2\,\mathrm{ln}(m))\sqrt{\pi }\\ \,\times \mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{u})-\mathrm{ln}(m))}{\sigma }\right)\\ -\,\displaystyle \sum _{i=1}^{n}\mathrm{ln}({x}_{i})(\mathrm{ln}({x}_{i})-2\,\mathrm{ln}(m))\sqrt{\pi }\\ \,\times \mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{l})-\mathrm{ln}(m))}{\sigma }\right),\end{array}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}\begin{array}{rcl}D & = & \sqrt{\pi }\left(-\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{l})-\mathrm{ln}(m))}{\sigma }\right)\right.\\ & & +\,\left.\mathrm{erf}\left(\displaystyle \frac{1}{2}\,\displaystyle \frac{\sqrt{2}(\mathrm{ln}({x}_{u})-\mathrm{ln}(m))}{\sigma }\right)\right){\sigma }^{3}.\end{array}\end{eqnarray} \tag{ 24 }$$