© ESO 2021

1. Introduction

Stars are formed from the gravitational infall of material in molecular clouds, which produces many stars with masses lower than the mass of the Sun. A substantially greater amount of material is needed for stars with masses higher than the Sun, and so these stars are formed in much fewer numbers. Stars undergo nuclear reactions in their cores for masses m ≳ 0.08 M_⊙, and thus the majority of all stars have masses from this lower limit up to the solar mass.

The stellar luminosity function and the stellar mass function are used to make comparisons across all types of stars, and these functions depend sensitively on the luminosity–mass relation. This relation has long been considered a smooth and differentiable function that exhibits structure correlating to the features of the stellar luminosity function (Kroupa et al. 1990). For example, an inflection in the luminosity–mass relation at ≈0.30 M_⊙ was first discussed by Copeland et al. (1970) as the H₂ dissociation zone having a lower adiabatic gradient. The minimum in the first derivative of the luminosity–mass relation, dm/dM_bol, at this mass was found to correspond to the maximum of the observed luminosity function of low-mass stars (Kroupa et al. 1990; Kroupa & Tout 1997) due to H₂ formation as well as stars becoming fully convective at this mass. Thus, all mono-age and mono-metallicity stellar populations show a pronounced and sharp maximum in the stellar luminosity function at a visual absolute magnitude of M_V ≈ 12 (Kroupa 2002; Kroupa et al. 2013). The luminosity–mass relation has long been known to be metallicity dependant (Elson et al. 1995; Kroupa & Tout 1997).

A recently discovered and related feature, called the Gaia M-dwarf gap, or Jao gap, was found in the observational Gaia Data Release 2 (DR2) data (Jao et al. 2018) and coincides with the same stellar mass of the inflection in the luminosity–mass relation at ≈0.30 M_⊙. This feature lies at the lower end of the main sequence (MS) in the G_BP − G_RP, M_G colour-magnitude diagram (CMD) and represents a deficiency in the density of stars (and subsequently a dip in the luminosity function). This gap was recently confirmed in the Gaia Early Data Release 3 (EDR3) data (Gaia Collaboration 2021) at Gaia magnitudes 2.2 < G_BP − G_RP < 2.8 and 10.0 < M_G < 10.3 and portrays a drop in density by 17 ± 6% (Jao et al. 2018). It is postulated that this occurs due to M-dwarf stars with masses around ≈0.30−0.35 M_⊙ undergoing the transition from being partially to fully convective (Jao et al. 2018; Baraffe & Chabrier 2018).

The thermonuclear process for low-mass stars is governed by the proton-proton I (ppI) chain:

$$\begin{array}{cc}& p+p⟶d+e^{+}+{\nu }_{\mathrm{e}}\hfill \end{array}$$(1)

$$\begin{array}{cc}& p+d⟶{}^3\mathrm{He}+\gamma \hfill \end{array}$$(2)

$$\begin{array}{cc}& {}^3\mathrm{He}+{}^3\mathrm{He}⟶{}^4\mathrm{He}+2p\hfill \end{array}$$(3)

where two protons (p) combine to make deuterium (d) by releasing a positron (e⁺) and an electron neutrino (v_e), and ultimately synthesise a ⁴He nucleus with the release of a photon (γ) and via two ³He nuclei. Equation (3) becomes important once the central temperature is T ≳ 7 × 10⁶ K (Dearborn et al. 1986), that is, m ≳ 0.26 M_⊙. Without this reaction, the energy produced by the other two reactions (Eqs. (1) and (2)) is not great enough for the material in the core to become unstable to convection (Baraffe & Chabrier 2018), which occurs when the radiative temperature gradient is larger than the adiabatic gradient, ∇_rad > ∇_ad, where

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mathrm{rad}}\propto L\kappa /T^4.\end{array}$$(4)

When the reaction given by Eq. (2) dominates, the ³He abundance grows, and at a temperature T ≳ 7 × 10⁶ K, Eq. (3) produces enough energy in the core for the centre to become convective (Chabrier & Baraffe 1997). Then, for example, as a 0.30 M_⊙ star approaches the MS, it has a radiative core, but the temperature increase due to the pre-MS contraction means that the production and destruction of ³He by Eqs. (2) and (3) releases enough energy for the material in the core to become unstable to convection; as such, the star reaches the MS with a convective core and envelope, which are separated by a thin radiative layer (Ezer & Cameron 1967).

Before the discovery of the M-dwarf gap, van Saders & Pinsonneault (2012) theorised that at the mass directly above the convective boundary, stars undergo the ‘convective kissing instability’. As ³He is destroyed at a slower rate by Eq. (3) than is produced by Eq. (2) at these low-mass temperatures (Baraffe & Chabrier 2018), it cannot reach equilibrium abundance. The abundance of ³He increases, and the core grows in size until it comes into contact with the convective envelope, resulting in periods of full convection (van Saders & Pinsonneault 2012). The mixing during these periods carries ³He out towards the surface, nuclear reactions in the core subside, and the core contracts and becomes separated again from the envelope (van Saders & Pinsonneault 2012). This continues periodically, with variations in the star’s luminosity and effective temperature, until the abundance of ³He is high enough throughout the star that it remains fully convective for the remainder of its lifetime. The variations from this convective kissing instability are thought to result in the M-dwarf gap in the density of observed stars seen at these temperatures and luminosities. Baraffe & Chabrier (2018) produced models in steps of 0.01 M_⊙ and note the merging of the core and envelope in their 0.34 M_⊙ and 0.36 M_⊙ models, along with a decrease in central ³He abundance and a change in the slope of the luminosity–mass relation at these masses. MacDonald & Gizis (2018) also find this relation and a dip in the luminosity function for their 0.33−0.35 M_⊙ models. Feiden et al. (2021) reproduced the M-dwarf gap in a CMD made from population synthesis models and also find periodic pulsations in luminosity, radius, and core temperature due to the ³He instability.

The Gaia M-dwarf gap is prominent across the blueward edge of the MS, where lower-metallicity stars are found. Feiden et al. (2021) show that their models with low metallicity ([Fe/H] = −0.7) encounter this instability at lower masses (and thus cooler core temperatures) than models with higher metallicities; from this they conclude that the dependence of temperature on stellar mass is one of the critical factors for this instability. At low temperatures, the opacity dependence of the radiative temperature gradient (Eq. (4)) becomes important. For lower metallicity, the opacity is decreased and the central temperature needed for the energy transport in the core to occur by radiative transport is lower (Chabrier & Baraffe 1997) and thus requires less mass. As such, the convective kissing instability is expected to happen at lower masses for lower metallicities since the higher-mass models have the temperature necessary for radiative cores. This work aims to investigate the dependence of the M-dwarf gap on the very low metallicities found in old globular clusters. Based on this, we suggest that the M-dwarf gap should be well represented in metal-poor globular cluster CMDs, although the low stellar masses at which this effect is found will be difficult to detect in clusters at great distances.

This work uses stellar evolution models to reproduce the convection-mass transition and investigate a fine grid of masses and very low metallicities with a short-scale time step. These models are described in Sect. 2, the results are presented in Sect. 3, and we conclude with a discussion in Sect. 4.

2. Stellar models

The 1D stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA, version 15140) was used (Paxton et al. 2011, 2013, 2015, 2018, 2019). The equations of state (EOSs) used in MESA are a combination of the SCVH (Saumon et al. 1995), OPAL (Rogers & Nayfonov 2002), HELM (Timmes & Swesty 2000), FreeEOS (Irwin 2004), and PC (Potekhin & Chabrier 2010) EOSs. The opacity tables relevant to this work are given by Grevesse & Sauval (1998) and OPAL (Iglesias & Rogers 1993, 1996); additionally, the electron conduction opacities are from Cassisi et al. (2007) and the opacity tables for the low temperatures of M-dwarf atmospheres in which molecules and dust grains are able to form are from Ferguson et al. (2005). For the outer boundary conditions we adopted the prescriptions provided by the Hauschildt et al. (1999) model atmosphere tables, taken at τ = 100 (Paxton et al. 2011; Chabrier & Baraffe 1997). The nuclear reaction rates are from Cyburt et al. (2010), and we included additional weak reaction rates (Fuller et al. 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2000). Screening was included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from Itoh et al. (1996). The initial helium abundance is Y = 0.24 + 2Z by mass fraction. Material becomes unstable to convection under the Schwarzschild criterion, and convective energy transport is described using mixing length theory.

Three initial metallicities were chosen (Z = 0.01, Z = 0.001, and Z = 0.0001) in order to investigate the M-dwarf gap on the blue edge of the MS and investigate the low metallicities found in globular clusters. An initial set of models for each metallicity was computed using a mass step of 0.01 M_⊙ to determine at which mass the models begin to have radiative cores. Once found, model sets were then produced using a smaller mass step of 0.00025 M_⊙. The models were calculated from the zero-age main sequence (ZAMS) up to 10 Gyr with a maximum time step of 50 000 years. One final model with an initial mass of 0.28 M_⊙ and Z = 0.0001 was computed up to 1 Gyr using a maximum time step of 1000 years.

3. Results

The low-mass MS models slowly undergo nuclear burning via the ppI chain described by Eqs. (1)–(3). By the time they reach an age of 10 Gyr, the hydrogen abundance in the core has only just begun to deplete. The Hertzsprung-Russell (HR) diagrams in Fig. 1 (and also later in Fig. 4) show the models experiencing a period of fluctuations that appear as loops along their evolutionary tracks during the MS. The age of the model at which this occurs increases with increasing mass. The mass range for this instability occurs at a higher mass for a higher metallicity, and the range of masses is larger for the higher-metallicity models. For the masses below the convective boundary, the models are fully convective throughout their lifetimes. Between a certain mass range for each metallicity (given in Table 1), the models begin their lives with a convective core, a thin radiative layer, and a convective envelope. They then undergo a period whereby they repeatedly switch from this structure to one of full convection and then relax into a fully convective state for the remainder of their lifetimes. This is illustrated in Fig. 2, which shows the mass coordinates of the convective and radiative regions of the Z = 0.001 models over time. At 1 Gyr, the models with mass m < 0.28600 M_⊙ are fully convective, models with mass 0.28600 M_⊙ ≤ m ≤ 0.29845 M_⊙ have begun the convective kissing instability, and the models with higher masses have a stable convective core and radiative region. The 0.28550 M_⊙, 0.28900 M_⊙, and 0.29450 M_⊙ models are also fully convective at 1 Gyr, having undergone a merger of the convective core and envelope. At 3 Gyr, models with mass m < 29275 M_⊙ and the 0.29550 M_⊙ model that has merged are fully convective. At 5 Gyr, the models with mass m < 0.29500 M_⊙ along with the 0.29525 M_⊙ and 0.29550 M_⊙ models that are merged at this time are fully convective. By 9 Gyr, models with mass m < 0.29775 M_⊙ are fully convective, and models with mass m > 0.31150 M_⊙ have a radiative core.

Fig. 1. Evolutionary tracks shown in the HR diagram for three sets of stellar models with metallicities Z = 0.0001 (top), Z = 0.001 (middle), and Z = 0.01 (bottom). The models begin their evolution on the left and move upwards along the tracks to reach the final points at 10 Gyr. The stars are on the MS throughout this evolution.

Fig. 2. Mass coordinates of the convective (blue) and radiative (grey dotted) regions within the Z = 0.001 models over 1, 3, 5, 7, and 9 Gyr. The models that experience a merger of the convective core and envelope can be seen as gaps in the radiative region. A movie of the full time lapse is available online.

3.1. Convective kissing instability

The luminosity, effective temperature, and central ³He abundance over time are displayed in Fig. 3 for several models with a metallicity of Z = 0.001. For the sake of clarity, some models are not shown. At this metallicity, the models up to 0.30025 M_⊙ undergo fluctuations in luminosity and temperature, some over the course of several billion years, which correlate to the increase and decrease of ³He in the core. There is a clear correlation between the dips in central ³He abundance and both the peaks in luminosity and the troughs in temperature. The amplitude of the fluctuations generally decreases over time until the models settle into a fully convective state; however, the models also exhibit large changes in luminosity and temperature, corresponding to the larger loops seen in the HR diagrams in Fig. 1. Figure 4 shows the convective kissing instability for the 0.296 M_⊙, Z = 0.001 model. The top panel shows the HR diagram, with each point representing 50 000 years, illustrating that the loops occur over a long period of time. The next panels down in the figure give the temperature, luminosity, radius, and ³He abundance profiles over time, along with the mass coordinate positions of the top of the convective core and the bottom of the convective envelope; the region in between is radiative. Again, the correlation between these properties is clear. When the bottom of the convective envelope reaches down towards the core, the central ³He abundance drops as the material is carried out towards the surface, and correspondingly the surface ³He abundance rises. This results in a decrease in luminosity as the ³He production is reduced. The temperature rises due to the now fully convective model contracting to compensate for the loss in nuclear energy output, which is represented by the drop in radius. Once the temperature is high enough to resume ³He production, the abundance of ³He in the core begins to rise along with the luminosity, and the model can increase again in radius. As the radius increases, the surface temperature drops. This repeated process, which results in the fluctuations in temperature and luminosity, produces the loops seen in the HR diagrams in Fig. 1. After ³He reaches its equilibrium abundance, the fluctuations subside and the models remain fully convective. The time period during which the models undergo the convective kissing instability, as well as the amplitude in luminosity and temperature fluctuations, increases with increasing mass. The 1000 year time step model showed that the smallest pulsations took place on timescales of a few thousand years.

Fig. 3. Temperature (top) luminosity (middle), and central 3He abundance (by mass fraction; bottom) over time for models with a metallicity of Z = 0.001. The fluctuations occur for several billion years for the higher-mass models.

Fig. 4. HR diagram (top) for the 0.296 M⊙, Z = 0.001 model, where each point represents 50 000 years. Also displayed are: the temperature, luminosity, and radius profiles for the same model (middle); and the centre and surface 3He abundances (by mass fraction), along with the mass coordinate of the bottom of the convective envelope and the top of the convective core (bottom). A radiative region is found in between. After about 7 Gyr, this model has evolved to be fully convective.

The models additionally exhibit several events where the bottom of the convective envelope reaches down towards the core but remains disconnected from it. These near-merger events correspond to small drops in central ³He abundance and increases in surface ³He, and they occur over a period of ≈1 Myr.

3.2. Metallicity dependence

The main effect due to metallicity is the convective kissing instability occurring at larger masses for models with higher metallicity. Correspondingly, the effective temperature is also higher for models with increased metallicity. Figure 5 shows the effect of metallicity for one chosen model with mass m = 0.290 M_⊙. This figure includes the HR diagrams along with the temperature, luminosity, radius, and central ³He abundance over time for metallicities of Z = 0.001 and Z = 0.0001. A comparison to Z = 0.01 could not be made here as the lowest mass to undergo the convective kissing instability for this metallicity is 0.32650 M_⊙. The higher-Z model begins the instability at an earlier age, around 0.3 Gyr, relative to the low-Z model, which begins at around 0.8 Gyr. The fluctuation period subsides and the higher-Z model becomes fully convective by 2.5 Gyr, compared to 5.5 Gyr for the low-Z model (more than twice as long). The amplitude of the fluctuations found in the luminosity, temperature, radius, and central ³He abundance are much higher for the low-Z model.

Fig. 5. HR diagrams for the 0.290 M⊙ model with Z = 0.0001 (top left) and Z = 0.001 (top right). Each point represents 50 000 years. Also shown are the temperature, luminosity, radius, and central 3He abundance profiles (bottom grid) for the two models.

3.3. Luminosity–mass relation

Figure 6 shows the luminosity–mass relation for the three metallicities investigated, at the model ages of 1, 3, 5, 7, and 9 Gyr. A discontinuity is seen at all metallicities. The discontinuity is present over many billions of years and is more prominent at the lower age. The mass ranges in which the discontinuity occurs are given in Table 1. By 9 Gyr, the lower-mass models that underwent the convective kissing instability at 1 Gyr have settled into a fully convective state, and as such the discontinuity shifts to higher masses and higher luminosities over time. Additionally, the discontinuity occurs at higher masses and luminosities for higher metallicities. Interestingly, the discontinuity is greatly reduced at 7 Gyr for Z = 0.01 to little more than a dip in the luminosity–mass relation, but then appears again at 9 Gyr. Also for Z = 0.01, the highest-mass model to undergo the instability remains 0.34325 M_⊙ at 3, 5, and 7 Gyr, whereas for the other metallicities the highest mass increases over time.

Fig. 6. Luminosity–mass relation for metallicities Z = 0.0001 (top), Z = 0.001 (middle), and Z = 0.01 (bottom) at model ages 1, 3, 5, 7, and 9 Gyr. A movie of the full time lapse is available online for the Z = 0.001 model set.

4. Discussion

The models in this work undergo luminosity, temperature, and radial fluctuations due to the convective kissing instability that was predicted by van Saders & Pinsonneault (2012), and found in models by Baraffe & Chabrier (2018) and Feiden et al. (2021) to be responsible for the Gaia M-dwarf gap. This work finds these fluctuations in the form of loops in the evolutionary tracks in the HR diagrams for models with metallicities down to Z = 0.0001, and our results agree with the conclusion that low-mass stars undergo these periods of fluctuations due to ³He production and transport.

A dip in the luminosity–mass relation was found by MacDonald & Gizis (2018), and by using smaller mass steps, we see that there is in fact a discontinuity. We find that the discontinuity occurs at lower masses and luminosities for the lower-metallicity models, as a direct consequence of the convective kissing instability developing at higher masses for higher metallicities. The discontinuity is a result of the merging of the convective envelope and core with the production of ³He, and of the model subsequently and periodically switching between partially and fully convective states. This affects the nuclear energy output, which results in the disruption in the luminosity–mass relation as well as the Gaia M-dwarf gap. This discontinuity, which is present over many billions of years, has a great impact on stellar physics as the luminosity–mass relation is no longer differentiable at these masses. For example, in order to constrain the stellar mass function, the slope of the luminosity–mass relation is used (Kroupa et al. 1990, 2013). This stellar mass function provides information on the nature of low-mass star formation and the contribution of low-mass stars to the dark matter problem (Kroupa et al. 1990).

The extra model with a maximum time step of 1000 years that ran up to 1 Gyr was calculated in order to check if the 50 000 year time step was able to provide an acceptable resolution for the computational time taken. The 1000 year time step model did reveal that the smallest pulsations took place on timescales of a few thousand years. However, the HR diagram as well as the luminosity and temperature profiles of the 1000 year model were visually comparable to the 50 000 year time step model. As the 1000 year model took 82 hours to complete, the computationally quicker 50 000 year time step was adequate for our purposes.

Feiden et al. (2021) find that the merger events between the convective core and envelope happen four to six times before their models settle into a fully convective state; however, the small time step used in this work reveals a much higher number of merger events. Additionally, some models experience several near-merger events during which the convective envelope reaches down into the radiative layer but not all the way to the core. These mixing events transport ³He-rich material into the radiative region, into which the convective envelope reaches and pulls ³He outwards to the surface. Whilst the core and envelope do not merge during these events, there exists a million-year-long instability between the two that results in small fluctuations in the ³He abundance, as well as in the surface temperature and luminosity. Further investigation is needed to fully resolve and determine the processes occurring during these near-merger events which was beyond the scope of this work.

The convective kissing instability occurs at higher masses for higher metallicities. This was expected due to the increased opacity for higher metallicities. Convection occurs when the radiative temperature gradient is higher than the adiabatic gradient, and by increasing the opacity, a higher temperature (and thus mass) is needed for the growth of a radiative core than for a lower opacity. For lower metallicities, the central temperature required for radiative transport is lower (Chabrier & Baraffe 1997), and thus the low-Z models have a radiative core at lower masses than the models with high metallicities. For a given mass, the convective kissing instability occurs earlier in the model lifetime for a higher metallicity since the evolution of the central temperature occurs earlier (see Chabrier & Baraffe 1997, Fig. 10a). Additionally, the convective kissing instability occurs for longer portions of the model lifetime for lower metallicities (Fig. 5), in agreement with Feiden et al. (2021) and Jao et al. (2018), who note that the observational M-dwarf gap is more prominent at the blue edge of the MS in the HR diagram that represents stars of low metallicities. The instability is found in models with metallicities down to Z = 0.0001, which would be well represented in old globular clusters; this is especially true since the instability period was found to be almost twice as long for the Z = 0.0001 models relative to the Z = 0.001 models. It would be of great importance to find the M-dwarf gap in the observations of these star clusters, although the great distances of the clusters create difficulties in finding these faint low-mass stars. Infrared observations of the nearest globular clusters may be able to achieve this.

5. Summary

Stellar evolution models were computed using three different metallicities over a range of masses in order to reproduce the M-dwarf gap observed in the lower region of the MS. The models used a mass step of 0.00025 M_⊙ and a time step of 50 000 years and ran until 10 Gyr. A final model used a 1000 year time step up to 1 Gyr; however, this was computationally time consuming, and the time step of 50 000 years was adequate for the purposes of this investigation. A discontinuity was found in the luminosity–mass relation for masses 0.32650−0.34350 M_⊙ for the metallicity Z = 0.01, masses 0.28600−0.30025 M_⊙ for Z = 0.001, and masses 0.28275−0.29525 M_⊙ for Z = 0.0001. The HR diagrams for these models show loops during the MS, which results in fluctuations in luminosity and effective temperature. The lower-metallicity models undergo the convective kissing instability at lower masses than the higher-metallicity models, due to the lower opacity and thus the lower temperature (and mass) required for radiative transport in the core. The instability occurs for a longer portion of the model’s lifetime for a given mass with decreasing metallicity.

Comparing the central and surface ³He abundances, as well as the mass coordinates of the convective core and envelope, the loops and pulsations correspond to the merging of the core and envelope, resulting in the switching of the model from a partially to fully convective state and the transport of ³He from the core to the surface. These loops would lead to the deficiency of stars in the Gaia M-dwarf gap and are due to the production and mixing of ³He throughout the star from the merging of the convective core and envelope.

Movie

Movie 1 associated with Figs. 2 and 6 (0_001Z) Access here

Acknowledgments

We thank the referee for the detailed and thoughtful comments. PK acknowledges support from the Grant Agency of the Czech Republic under grant number 20-21855S.

References

All Tables

All Figures

	Fig. 1. Evolutionary tracks shown in the HR diagram for three sets of stellar models with metallicities Z = 0.0001 (top), Z = 0.001 (middle), and Z = 0.01 (bottom). The models begin their evolution on the left and move upwards along the tracks to reach the final points at 10 Gyr. The stars are on the MS throughout this evolution.
In the text

	Fig. 2. Mass coordinates of the convective (blue) and radiative (grey dotted) regions within the Z = 0.001 models over 1, 3, 5, 7, and 9 Gyr. The models that experience a merger of the convective core and envelope can be seen as gaps in the radiative region. A movie of the full time lapse is available online.
In the text

	Fig. 3. Temperature (top) luminosity (middle), and central 3He abundance (by mass fraction; bottom) over time for models with a metallicity of Z = 0.001. The fluctuations occur for several billion years for the higher-mass models.
In the text

Fig. 4. HR diagram (top) for the 0.296 M⊙, Z = 0.001 model, where each point represents 50 000 years. Also displayed are: the temperature, luminosity, and radius profiles for the same model (middle); and the centre and surface 3He abundances (by mass fraction), along with the mass coordinate of the bottom of the convective envelope and the top of the convective core (bottom). A radiative region is found in between. After about 7 Gyr, this model has evolved to be fully convective.

In the text

	Fig. 5. HR diagrams for the 0.290 M⊙ model with Z = 0.0001 (top left) and Z = 0.001 (top right). Each point represents 50 000 years. Also shown are the temperature, luminosity, radius, and central 3He abundance profiles (bottom grid) for the two models.
In the text

	Fig. 6. Luminosity–mass relation for metallicities Z = 0.0001 (top), Z = 0.001 (middle), and Z = 0.01 (bottom) at model ages 1, 3, 5, 7, and 9 Gyr. A movie of the full time lapse is available online for the Z = 0.001 model set.
In the text