Simulating Brown Dwarf Observations for Various Mass Functions, Birthrates, and Low-mass Cutoffs

A crucial parameter of brown dwarf formation is the value of the low-mass cutoff, which has been shown to be no higher than $\sim 10M_{Jup}$ (Kirkpatrick et al. 2021). Objects such as WISE J085510.83$-$071442.5, which is estimated to a have a mass between $1.5M_{Jup}$ and $8M_{Jup}$ depending on its age (Leggett et al. 2017), along with objects identified in young moving groups (see below), almost certainly push this limit lower, as seen by derived low-mass cutoffs of $\sim 4M_{Jup}$ in Bate & Bonnell (2005).

The lowest mass at which brown dwarfs form is a fundamental property of star formation at the very edge of our theoretical understanding of brown dwarfs. Notably, a lower mass cutoff not only extends the mass range in which brown dwarfs may form, but also shifts the mode of the distribution to lower masses. These faint objects that would populate the low mass end of the mass function are predominantly late-T and Y dwarfs – as seen in Figure 6 of Burgasser (2004). The observed space density in temperature bins below $750K$ has the greatest deciding influence on the value of the low-mass cutoff, as lower-mass cutoffs will more heavily populate objects at the lowest temperatures. Thus, data at these coolest temperatures will be most influential in determining the low-mass cutoff.

Due to the faint nature of late-T and Y dwarfs it is difficult to complete a volume-limited sample with sufficient statistics to provide a robust space density measurement. Since the lowest temperature bins are of paramount importance for the evaluation of the low-mass cutoff, we therefore need additional discoveries of faint, cold Y dwarfs in order to further constrain the value of the low-mass cutoff.

For this study, we choose $0.01\text{ M}_{\odot}$, $0.005\text{ M}_{\odot}$, and $0.001\text{ M}_{\odot}$ as the low-mass cutoffs within our framework, as was done in Kirkpatrick et al. (2021). These values produce populations that include low-mass brown dwarfs that either straddle or are below the deuterium-burning limit ($\sim 13M_{Jup}$; Spiegel et al. 2011). There are precedents for such brown dwarfs. Take, for example, the low-mass brown dwarfs SIMP J013656.5+093347.3 and 2MASSW J2244316+204343. SIMP J013656.5+093347.3 is a young early-T dwarf with an estimated mass of $12.7\pm 1.0M_{Jup}$, derived using its moving group association and evolutionary models (Saumon & Marley 2008), and a trigonometric distance of $6.139\pm 0.037\text{ pc}$ (Gagné et al. 2017). 2MASSW J2244316+204343 is a mid-L dwarf with a mass of $10.46\pm 1.49M_{Jup}$ (Faherty et al. 2016), also derived from evolutionary models, and a kinematic distance of $18.5\pm 1.2\text{ pc}$ (Liu et al. 2016). However, both objects are close to the Solar System and therefore less of a challenge for the current instrumentation to observe, unlike further, fainter brown dwarfs.

Integrating Equation 2, supposing $\alpha\neq 1$, to find our Cumulative Distribution Function (CDF), we get the following expression, where $M$ is the mass parameter. Our CDF, once normalized and inverted, will serve as a key component of the inverse transform method which we utilize. Here, $m_{2}$ is the high mass cutoff, defined to be $0.1M_{\odot}$, and $m_{1}$ we vary between $0.01M_{\odot}$, $0.005M_{\odot}$, $0.001M_{\odot}$.

In order to derive our constant of normalization, $C_{N}$, we need our CDF to evaluate to $1$ given the higher mass limit. Solving for the constant, we find the following:

Similarly, when $\alpha=1$ in Equation 2, the constant of normalization is the following.

Once inverted and with the value of $C_{N}$ inserted, the equation for the $\text{CDF}^{-1}$ becomes the following (Kirkpatrick et al. 2019).

Here, $x\in\mathcal{U}_{[0,1]}$, meaning $x$ is randomly sampled from the uniform distribution between 0 and 1. Histograms of the $\text{CDF}^{-1}$ sampled $10^{6}$ times per low-mass cutoff threshold with various $\alpha$ values are shown in Figure 1.

The constant birthrate function is a common starting point by virtue of its inherent simplicity. A constant distribution implies that the Galaxy’s star formation processes have been consistently efficient and have had sufficient star forming material from its nascence to present-day. We adopt the following functional form for the constant distribution, in which $C_{\mathcal{T}}$ is the eponymous constant of star formation.

The value of $C_{\mathcal{T}}$ is not of much significance in our study as we ultimately normalize our CDF.

Given that the constant age distribution is, as its name implies, constant, it does not depend on any time parameter $\mathcal{T}$. Therefore in order to convert it from an age distribution we simply change $C_{\mathcal{T}}$ to $C_{\mathcal{A}}$ to indicate the change from a constant of time to a constant of age, instead of executing the coordinate transformation outlined in Equation 7.

The $\text{CDF}^{-1}$ for our constant distribution can be attained in a manner similar to §2.2.

We choose to leave $CDF^{-1}_{C}$ in terms of $a_{1}$ to $a_{2}$, the minimum and maximum ages in Gyr respectively, as our study explores more than one age range, see Appendix A.

The Inside-Out age distribution represents a sample population where star formation initiates within the central regions of the Galaxy and propagates outward with time (Bird et al. 2013). The functional form of the distribution is the following (Johnson et al. 2021).

We choose to linearly approximate this time distribution using Equation 12, since inverting the CDF of Equation 11 requires the use of the computationally expensive error function $\left(\text{erf}(x)=\frac{2}{\sqrt{\pi}}\int^{x}_{0}e^{-t^{2}}dt\right)$. Moreover, the original functional form is already sufficiently linear between our age bounds of 0 to 10 Gyr such that our approximation retains much of the original shape of the function.

By using Equation 7 to convert this time distribution to an age distribution, we arrive at the following functional forms:

Integrating and normalizing this PDF yields the following CDF for the Inside-Out birthrate.

We solve for $A$ as we have done previously to derive the inverse form of the CDF. The negative branch of the solution is disregarded as a negative age is physically inconceivable.

Galactic star formation need not have followed a constant rate, or even one that varies linearly like the Inside-Out. In the past 10 Gyr it is possible that periods of the star formation history of our Galaxy have been more intense than others, with otherwise linear reduction of the birthrate, as seen in the Inside-Out. This manifests itself as bursts of increased star formation, possibly due to gravitational perturbances from the Sagittarius dwarf galaxy (Ruiz-Lara et al. 2020) or from an earlier galactic merger incident inducing a starburst on our Galaxy (Helmi 2020).

The Late-Burst model accounts for such a period of starburst with a spike in the disk’s total birthrate between ages of 2.65 and 5.10 Gyr. Equation 16 displays the mathematical expression for the Late-Burst model. For ease of inversion, we approximate the Late-Burst as Equation 17, defined in terms of the two previous birthrates, $PDF_{C}$ and $PDF_{IO}$.

In order to extract a meaningful $CDF^{-1}$ from this, we integrate each piece of the Late-Burst and allocate samples to preserve the ratio between the two pieces. Thus, the $CDF^{-1}_{LB}$ is a mixture of $CDF^{-1}_{IO}$ and $CDF^{-1}_{C}$ whose individual sample sizes depend on the area under each individual PDF.

We choose $\alpha$ values ranging from 0.3 to 0.8 in increments of 0.1 as they envelop the previous best $\alpha$ value of 0.6 (Kirkpatrick et al. 2019) on either side. Using the $\text{CDF}^{-1}$ of our assumed mass distribution, we pull an object at random from the distribution to assign it a mass. This is done via a Monte Carlo draw from 0 to 1, and we perform $10^{6}$ of these draws to build a population with statistical robustness. We repeat this procedure for each value of $\alpha$, and for each value of $\alpha$ we repeat the procedure for each of our three assumed low-mass cutoffs. In total, we build eighteen simulated populations, each having masses for $n=10^{6}$ objects. To be specific, the code samples the mass function again for each different combination of birthrate and evolutionary model, so in total there are 162 simulated populations (Six values of $\alpha\times$Three mass cutoffs$\times$Three birthrates$\times$Three evolutionary models).

We similarly use the inverse transform method to pull random ages from each of our three assumed age distributions. This methodology allows for the creation of brown dwarf populations of arbitrary size whose mass distribution will converge to the shape of the transformed function as the sample size approaches a statistically significant value.

For each population, the $i\textsuperscript{th}$ element in each mass list and the $i\textsuperscript{th}$ element in the age list become the mass and age of the $i\textsuperscript{th}$ object in the simulated population.

For each object, we wish to find its current-day temperature using our assumed evolutionary model. However, given the discrete sampling of our evolutionary model grids, the simulated values of age and mass for our object are unlikely to have been included directly in the models. We therefore use bilinear interpolation to fill in the sample space between the model’s reference points. Not all points in the mass-age domain can be mapped using bilinear interpolation. At the edges of the space sampled by each evolutionary model there exist mass-age regions with points that cannot be enclosed within a rectangle of reference points. As shown in the left column of Figure 3, each evolutionary model has loci in which we cannot interpolate stellar temperatures. In such cases we simply disregard the star and assign to the star’s temperature value the number $-1$ to indicate that a temperature could not be interpolated. The extent to which we lose objects during the interpolation depends on the (non-)rectilinearity of the provided evolutionary sample set in mass-age space. Both the Sonora and Phillips models are fairly well sampled and do not drop many brown dwarf samples along the whole range of possible mass and age values. In contrast, the Saumon & Marley (2008) model drops the most objects of our three evolutionary models, especially those sample points which are either young and massive, or old and light, as seen in the top left and bottom right of Figure 3. However, it should be noted that for the temperature range we consider for our fitting, namely 450K to 2100K, there are extremely few samples dropped due to a lack of rectilinear bounding reference points (see top row of Figure 3)

Although we state in § 5.1, that we simulated $n=10^{6}$ objects, in practice we simulated $n>10^{6}$ objects, and kept only the first $10^{6}$ objects for which temperatures could be obtained via bilinear interpolation.

For our Late-Burst birthrate, we simulate brown dwarf subpopulations as explained at the end of §3.3. We shuffle these proportionally sampled constant and Inside-Out birthrate subpopulations and select the first $10^{6}$ objects to create the Late-Burst birthrate.

The fact that the evolutionary model grids are not sampled over the entire mass and age space needed means that our final, simulated populations contain small biases. See Figure 4 for an example of how our original mass distribution changes after interpolation.

We now shift focus to comparing our simulated populations to the empirical temperature distribution. We analyze our results in two steps. First, in §6.1.1 we describe our methods for comparing the simulated and observed temperature distributions to determine which low-mass function $\alpha$ value fits best. Then, in §6.1.2, we evaluate which mass cutoff leads to the best fit.

A key feature of the Saumon & Marley (2008) evolutionary model is its incorporation of cloud formation at the L-T transition, in which L dwarfs cool into T dwarfs. This process is mainly limited to the temperature range from 1200-1350K, and in this temperature bin there is a noted increase of objects, forming a bump in the temperature distribution (Figure 6 and Figure 13 of Kirkpatrick et al. 2019). The exact physical conditions and processes that lead to this surplus of objects at 1200-1350K are not yet fully understood, but one theory suggests that the dispersion of the cloud layers could be driven by a radiative cloud-induced variability (Tan & Showman 2019).

Figures 7 and 8 show the temperature distribution of our best-fit model colored by object age, showing an excess of young brown dwarfs at the L/T transition bin (1200-1350K). This trend of younger brown dwarfs around the L/T transition is also visible in Figure 9, where we display the median age and its standard deviation per bin. These predictions suggest a pile-up of young objects just prior to the L/T transition, indicating that the cooling time of a brown dwarf is significantly slowed in this region. Observational confirmation of this effect may be possible once we are able to collect a large field sample of brown dwarf age estimates or, perhaps more easily, measuring the temperature distribution of L and T dwarfs belonging to young clusters and associations of known age.

The composition of each of our simulated populations depends heavily on our choice of mass function, birthrate, low-mass cutoff, and evolutionary model. In this section, we show the variation in the resulting temperature distribution when we hold all but one of these parameters constant.

The greatest change in the temperature distribution results from a change in the evolutionary model. Notably, simulations from the Saumon & Marley (2008) models possess a bump at the L/T transition, whereas those from both the Sonora and Phillips models do not.

As we vary the mass function $\alpha$ parameter, the shape of the temperature distribution also predictably varies (Figure 10). Flatter mass functions ($\alpha\sim 0.4$) lead to temperature distributions with relatively hotter objects whereas steeper mass functions ($\alpha\sim 0.7$) lead to a greater abundance of cooler objects. Also, flatter mass functions imply a larger concentration of objects at the L/T transition with a lesser low-temperature peak (300K-600K) and vice versa for steeper mass functions. It should be noted that other differences are marginal everywhere except the low-temperature peak, where the difference between the $\alpha=0.3$ and $\alpha=0.8$ distributions is pronounced. Fundamentally, increasing the mass function’s steepness serves to skew the resulting temperature distribution towards the cooler end of the temperature regime.

Differences in birthrate functions have already been shown to affect the resulting temperature distribution only marginally (Burgasser 2004). Our findings corroborate this (Figure 11).

Nonetheless, the Inside-Out age function, for example, allows for a flatter mass function to fit the empirical temperature distribution. The $\alpha=0.6$ value taken as the ideal mass function steepness in Kirkpatrick et al. (2019) assumed a constant birthrate, and our findings in Table 1 replicate that while also showing that a combination of an $\alpha=0.4$ or $\alpha=0.5$ paired with an Inside-Out birthrate also fit the empirical data quite well. This is because the declining birthrate represented by the Inside-Out function paired with a less steep (lower $\alpha$) mass function can create as many present-day late-T and Y dwarfs as a constant birthrate paired with a steeper mass function.

The mass cutoff also does not in any significant way affect the shape of the simulated temperature distribution except at the coldest temperatures (Figure 12).