Analytic Models of Brown Dwarfs and The Substellar Mass Limit

We present the current status of the analytic theory of brown dwarf evolution and the lower mass limit of the hydrogen burning main sequence stars. In the spirit of a simplified analytic theory we also introduce some modifications to the existing models. We give an exact expression for the pressure of an ideal non-relativistic Fermi gas at a finite temperature, therefore allowing for non-zero values of the degeneracy parameter ($\psi = \frac{kT}{\mu_{F}}$, where $\mu_{F}$ is the Fermi energy). We review the derivation of surface luminosity using an entropy matching condition and the first-order phase transition between the molecular hydrogen in the outer envelope and the partially-ionized hydrogen in the inner region. We also discuss the results of modern simulations of the plasma phase transition, which illustrate the uncertainties in determining its critical temperature. Based on the existing models and with some simple modification we find the maximum mass for a brown dwarf to be in the range $0.064M_\odot-0.087M_\odot$. An analytic formula for the luminosity evolution allows us to estimate the time period of the non-steady state (i.e., non-main sequence) nuclear burning for substellar objects. Standard models also predict that stars that are just above the substellar mass limit can reach an extremely low luminosity main sequence after at least a few million years of evolution, and sometimes much longer. We estimate that $\simeq 11 \%$ of stars take longer than $10^7$ yr to reach the main-sequence, and $\simeq 5 \%$ of stars take longer than $10^8$ yr.


Introduction
One of the most interesting avenues in the study of stellar models lies in understanding the physics of objects at the bottom of and below the hydrogen burning main sequence stars. The main obstacle in the study of very low mass (VLM) stars and substellar objects is their low luminosity, typically of order 10 −4 L ⊙ , which makes them difficult to detect. There is also a degeneracy between mass and age for these objects, which have a luminosity that decreases with time. This makes a determination of the initial mass function (IMF) difficult in this mass regime. However in the last two decades there has been substantial observational evidence that supports the existence of faint substellar objects. Since the first discovery of a brown dwarf (Oppenheimer et al. 1995;Rebolo et al. 1995) several similar objects were identified in young clusters  and Galactic fields (Ruiz et al. 1997) and have generated great interest among theorists and observational astronomers. The field has matured remarkably in recent years and recent summaries of the observational situation can be found in Luhman et al. (2007) and Chabrier et al. (2014).
In two consecutive papers, Kumar (1963a) and Kumar (1963b) revolutionized the understanding of low mass objects by studying the Kelvin-Helmholtz time scale and structure of very low mass stars. He successfully estimated that stars below about 0.1M ⊙ contract to a radius of about 0.1R ⊙ in about 10 9 years, which was a correction to the earlier estimate of 10 11 years. The earlier calculation was based on the understanding that low mass stars evolve horizontally in the H-R diagram and thus evolve with a low luminosity for a long period of time. However, Hayashi and Nakano (1963) showed that such low mass stars remain fully convective during the pre-main sequence evolution and are much more luminous than the previously accepted model based on radiative equilibrium. Kumar's analysis showed that for a critical mass of 0.09M ⊙ , the time scale has a maximum value that decreases on either side. Although this crude model neglected any nuclear reactions, it did give a very close estimate of the time scale. The second paper (Kumar 1963b) gave more detailed insight on the structure of the interior of low mass stars. This model was based on the non-relativistic degeneracy of electrons in the stellar interior. Kumar's extensive numerical analysis for a particular abundance of hydrogen, helium and other chemical compositions yielded a limiting mass below which the central temperature and density are never high enough to fuse hydrogen. A more exact analysis required a detailed understanding of the atmosphere and surface luminosity of such contracting stars.
The next major breakthrough in theoretical understanding came from the work of Hayashi and Nakano (1963), who studied the pre-main sequence evolution of low mass stars in the degenerate limit. Although it was predicted that there exist low mass objects that cannot fuse hydrogen, the internal structure of these objects remained a mystery. A com-plete theory demanded a better understanding of the physical mechanisms which govern the evolution of these objects. It became essential to develop a complete equation of state (EOS). D' Antona & Mazzitelli (1985) used numerical simulations to study the evolution of VLM stars and brown dwarfs for Population I chemical composition (Y = 0.25, Z = 0.02) and different opacities. Their model showed that for the same central condition (nuclear output) an increasing opacity reduces the surface luminosity. Thus, a lower opacity causes a greater surface luminosity and subsequent cooling of the object. Their results implied that the hydrogen burning minimum mass is M = 0.08M ⊙ for opacities considered in their model. Furthermore, they showed that objects with mass close to M = 0.08M ⊙ spend more than a billion years at a luminosity of ∼ 10 −5 L ⊙ . Burrows et al. (1989) modeled the structure of stars in the mass range 0.03M ⊙ −0.2M ⊙ . They used a detailed numerical model to study the effects of varying opacity, helium fraction, and the mixing length parameter, and compared their results with the existing data. Their important modification was that they considered thermonuclear burning at temperatures and densities relevant for low masses. A detailed analysis of the equation of state was performed in order to study the thermodynamics of the deep interior, which contained a combination of pressure ionized hydrogen, helium nuclei, and degenerate electrons. This analysis clearly expressed the transition from brown dwarfs to very low mass stars. These two families are connected by a steep luminosity jump of two orders of magnitude for masses in the range of 0.07M ⊙ − 0.09M ⊙ . There was a clear indication that masses in that intermediate regime do ignite hydrogen but that it eventually subsides to yield a brown dwarf. Saumon & Chabrier (1989) proposed a new EOS for fluid hydrogen that, in particular, connects the low density limit of molecular and atomic hydrogen to the high density fully pressure-ionized plasma. They used the consistent free energy model but with the added prediction of a first order "plasma phase transition" (PPT) (Saumon & Chabrier 1989) in the intermediate regime of the molecular and the metallic hydrogen. As an application of this EOS, they modelled the evolution of a hydrogen and helium mixture in the interior of Jupiter, Saturn, and a brown dwarf (Chabrier & van Horn 1992;). They adopted a compositional interpolation between the pure hydrogen EOS and a pure Helium EOS to obtain a H/He mixed EOS. This was based on the additive volume rule for an extensive variable (Fontaine et al. 1977) and allowed calculations of the H/He EOS for any mixing ratio of hydrogen and helium. Their analysis suggested that the cooling of a brown dwarf with a PPT proceeds much more slowly than in previous models (Burrows et al. 1989). Stevenson (1991) presented a detailed theoretical review of brown dwarfs. His simplified EOS related pressure and density for degenerate electrons and for ions in the ideal gas approximation. Although corrections due to Coulomb pressure and exchange pressure are of physical relevance, they together contribute less than 10% in comparison to the other dominant term in the pressure-density relationship for massive brown dwarfs (M ≥ 0.04 M ⊙ ). The theoretical analysis gave a very good understanding of the behavior of the central temperature T c as a function of radius and degeneracy parameter ψ. Stevenson (1991) discussed the thermal properties of the interior of brown dwarfs and provided an approximate expression for the entropy in the interior and in the atmosphere of a brown dwarf. He also derived an expression for the effective temperature as a function of mass.
A method to use the surface lithium abundance as a test for brown dwarf candidates was proposed by Rebolo et al. (1992). Lithium fusion occurs at a temperature of about 2.5 × 10 6 K, which is easily attainable in the interior of the low mass stars. However brown dwarfs below the mass of 0.065M ⊙ never develop this core temperature. They will then have the same lithium abundance as the interstellar medium independent of their age. However, for objects slightly more massive than 0.065M ⊙ , the core temperature can eventually reach 3 × 10 6 K. They deplete lithium in the core and the entire lithium content gets exhausted rapidly due to the convection. This causes significant change in the observable photospheric spectra. Thus lithium can act as a brown dwarf diagnostic (Basri et al. 1996) as well as a good age detector (Stauffer et al. 1998).
Following this, an extensive review on the analytic model of brown dwarfs was presented by Burrows and Liebert (1993). They presented an elaborate discussion on the atmosphere and the interior of brown dwarfs and the lower edge of the hydrogen burning main sequence. Based on the convective nature of these low mass objects, they modelled them as polytropes of order n = 1.5. Once again the atmospheric model was approximated based on a matching entropy condition of the plasma phase transition between molecular hydrogen at low density and ionized hydrogen at high density. The polytropic approximation enabled the calculation of the nuclear burning luminosity within the core adiabatic density profile (Fowler and Hoyle 1964). While the luminosity did diminish with time in the substellar limit, the model did show that brown dwarfs can undergo hydrogen burning for a substantial period of time before it eventually ceases. The critical mass deduced from this model did indeed match that obtained from more sophisticated numerical calculations (Burrows et al. 1989).
In this work we give a general outline of the analytic model of the structure and the evolution of brown dwarfs. We advance some aspects of the existing analytic model by introducing a modification to the equation of state. We also discuss some of the unresolved problems like estimation of the surface temperature and the existence of PPT in the brown dwarf environment. Our paper is organized as follows. In Section 2 we discuss the derivation of a more accurate equation of state for a partially degenerate Fermi gas. We incorporates a finite temperature correction to the expression for the Fermi pressure to give a more general solution to the Fermi integral. In Section 3 we discuss the scaling laws for various thermodynamic quantities for an analytic polytrope model of index n = 1.5. In Section 4 we discuss the derivation of the equations (Burrows and Liebert 1993) connecting the photospheric (surface) temperature with density, where the entropies at the interior and the exterior are matched using the first order phase transition. In the spirit of an analytic model we derive simplified analytic expressions for the specific entropies above and below the PPT. We also highlight the need to seek alternate methods given current concerns about the relevance of the PPT in BD interiors. We discuss the nuclear burning rates for low mass objects in Section 5 and determine the nuclear luminosity L N (Fowler and Hoyle 1964). In Section 6 we estimate the range of minimum mass required for stable sustainable nuclear burning. In Section 7, we discuss a cooling model and examine the evolution of photospheric properties over time. In Section 8, we estimate the number fraction of stars that enters the main sequence after more than a million years. In the concluding section, we discuss further possibilities for an improved and generalized theoretical model of brown dwarfs.

Equation of state
In main sequence stars, the thermal pressure due to nuclear burning balances the gravitational pressure and the star can sustain a large radius and nondegenerate interior for a long period of time. However substellar objects like brown dwarfs fail to have a stable hydrogen burning sequence and instead derive their stability from electron degeneracy pressure. A simple but accurate model needs to have a good equation of state that incorporates the degeneracy effect and the ideal gas behaviour at a relative higher temperature. Burrows and Liebert (1993) give a pressure law that applies to both the extremes but has a poor connection in the intermediate zone. Stevenson (1991) also gives an empirical relation for the pressure that does include an approximate correction term to connect the two extremes. Here, in order to obtain a more accurate analytic expression for the pressure, we integrate the Fermi-Dirac integral exactly using the polylogarithm functions Li s (x). The most general expression for the pressure is (1) (Padmanabhan 1999), where b = 1 for the Fermi gas and ǫ(p) = p 2 c 2 + m 2 c 4 −mc 2 and the other variables are the standard constants. For substellar objects, the electrons are mainly non-relativistic due to relatively low temperature and density. In the non-relativistic limit, i.e. m 2 c 4 ≫ p 2 c 2 , the energy density reduces to ǫ(p) ≃ p 2 2m . Now rewriting Eq. (1) in terms of the energy density gives where a = 2 3 4π(2m) 3 2 (2π ) 3 , β = (k B T ) −1 , and we have taken g s = 2. In the limit T → 0, for all ǫ < µ the argument of the exponential is negative and hence the exponential goes to zero as β → ∞. Thus the integral reduces to the Fermi pressure at zero temperature. However in a physical situation at finite temperature, the integral can be solved analytically using the polylogs. The details of the exact derivation for a general Fermi integral are shown in Appendix A. The expression for the pressure of a degenerate Fermi gas at finite temperature is The above expression for pressure is the most general analytic relation for the pressure of a degenerate Fermi gas at a finite temperature. The first term is the zero temperature pressure and the subsequent terms are the corrections due to the finite temperature of the gas, and include Li s , the polylogarithm functions of different orders s. The expression is terminated after the fourth term as the polylogs fall off exponentially as the gas becomes more and more degenerate. Eq. (3) is a natural extension of the first-order Sommerfeld correction (Sommerfeld 1928).
The central temperature of VLM stars and brown dwarfs are of the same order as that of the electron Fermi temperature and thus the degeneracy parameter ψ is defined as where µ F is the electron Fermi energy in the degenerate limit and 1 µe = X + Y 2 is the number of baryons per electron and X and Y are the mass fractions of hydrogen and helium respectively. Other constants have their standard meaning. Rewriting Eq. (3) in terms of the degeneracy parameter ψ and retaining terms only up to second order, we arrive (for µ = µ F ) at 2me is a constant. However, the interior of a brown dwarf is also composed of ionized hydrogen and helium. The total pressure is a combined effect of both electrons and ions, i.e., P = P F + P ion , where P F is the Fermi pressure for an ideal non-relativistic gas at a finite temperature. The pressure due to ions for an ionized hydrogen gas can be approximated as Therefore the final equation of state for the combined pressure is where α = 5µe 2µ 1 and µ 1 is the mean molecular weight for helium and ionized hydrogen mixture and is expressed as where x H + is the ionization fraction of hydrogen. It should be noted that x H + changes as one moves from the core (completely ionized) to the surface which is mainly composed of molecular hydrogen and helium.
There are several corrections to the EOS that can be considered. The Coulomb pressure and the exchange pressure (see Eq. (13) in Stevenson (1991)) are two important corrections to Eq. (7). However, as stated earlier they are less important for more massive brown dwarfs. Hubbard (1984) presents the contribution due to the electron correlation pressure, which depends on the logarithm of r e , the mean distance between electrons. Stolzmann et al. (1996) present an analytic formulation of the EOS for fully ionized matter to study the thermodynamic properties of stellar interiors. They show that the inclusion of both electron and the ionic correlation pressure results in a ∼ 10% correction to the EOS. Furthermore, Gericke et al. (2010) state that the main volume of the brown dwarfs and the interior of giant gas planets are in a warm dense matter state, where correlation energy, effective ionization energy and the electron Fermi energy are of the same order of magnitude, making it effectively a strongly correlated quantum system. Becker et al. (2014) give an EOS for hydrogen and helium covering a wide range of density and temperature. They extend their ab intio EOS to the strongly correlated quantum regime and connect it with the data derived using other methods for the neighboring regions of the ρ − T plane. These simulations are within the framework of density functional theory molecular dynamics (DFT-MD) and give a detailed description of the internal structure of brown dwarfs and giant planets. This leads to a 2.5% − 5% correction in the mass-radius relation. The study of the EOS of brown dwarfs will help in understanding degenerate bodies in the thermodynamic regime that is not so close to the high pressure limit of a fully degenerate Fermi gas. In this context, the Mie-Grueneisen equation of state is of relevance to test the validity of the assumption that the Grueneisen parameter γ = ∂ log T ∂ log ρ s is independent of the temperature T (Anderson 2000) at a constant volume V . The brown dwarf regime is in a way more interesting than the white dwarf regime since it is not so close to the limit of a fully degenerate Fermi gas. In Appendix C we have provided analytic expressions for two parameters that are of particular relevance for the brown dwarfs; the specific heat (C v or C p ) and the Grueneisen parameter.

An analytic model for brown dwarfs
In this section, we derive some of the essential thermodynamic properties of a polytropic gas sphere based on the discussion in Chandrasekhar (1939). As is evident from Eq. (7), the P − ρ relation for a brown dwarf is a polytrope where the index n = 3/2. K is a constant depending on the composition and degeneracy and can be expressed (from Eq. (7)) as where for a simplified presentation we represent the correction terms as and C 1 = 2 5 aA 5 2 is a constant. The solution to the Lane-Emden equation subject to the zero pressure outer boundary condition can be used to arrive at useful results for R, ρ c and P c for the polytropic equation of state, Eq. (7) The radius can be expressed as 1 On using the values of natural constants we get C = 10 13 cm 4 g −2/3 s −2 . (Chandrasekhar 1939). On substituting Eq. (10) for K, the radius for a brown dwarf can be expressed as the function of degeneracy and mass, Similarly, the expressions for the central density ρ c and central pressure are given by the relations ρ c = δ n ( 3M 4πR 3 ) and P c = W n GM 2 R 4 , where the constant δ n = 5.991 and W n = 0.77 for the polytrope of n = 1.5 (Chandrasekhar 1939). On substituting the expression for R (Eq. 13) in these relations we get and These are the scaling laws of the density and pressure in the interior core of a brown dwarf. Interestingly, these vary with the degeneracy parameter ψ that is a function of time. Thus a very simple polytropic model can yield the time evolution of the internal thermodynamical conditions of a brown dwarf. From the definition of the degeneracy parameter in Eq. (4) and using Eq. (14), the central temperature can be expressed as a function of ψ: The central temperature has a maximum for a certain value of ψ, and it increases for greater values of M. Further, using Eq. (13) we have shown the variation of central temperature T c as a function of radius R. T c increases as the object contracts under the influence of gravity. It peaks at a certain R and then cools over time. The maximum peak temperature increases for heavier objects and also depends on the extent of ionization of hydrogen and helium. Figure 1 shows the variation of the central temperature as a function of radius for different mass ranges. If the critical temperature for thermonuclear reactions is around 3 × 10 6 K, we can roughly estimate the critical mass for the main sequence as ∼ 0.085 M ⊙ . This is similar to the estimated critical mass (∼ 0.084 M ⊙ ) for the main sequence (see Figure 1 in Stevenson (1991)). However it should be noted that the estimate of minimum mass is very sensitive to the mean molecular mass µ 1 . In the Fig. 1 we have used µ 1 = 1.23 for fully neutral gas (similar toĀ ∼ 1.24 for cosmic mixture as used in Stevenson (1991)). Depending upon the value of µ 1 the minimum mass may vary significantly. For example, if we consider a fully ionised gas, i.e., µ 1 = 0.59, it yields a minimum mass of 0.12 M ⊙ .

Surface properties
In this section we discuss a very simple but crude model which is broadly based on the phase transition proposed by  and the isentropic nature of the interior of brown dwarfs. The development of a theoretical model for studying the variation of surface luminosity over time for low mass stars (LMS) and brown dwarfs is a great challenge. There is no stable phase of nuclear burning for brown dwarfs and the luminosity gradually decreases with time. This leads to an age-mass degeneracy in observational determinations. Our lack of knowledge in understanding the physics of the interior of brown dwarfs restricts the development of a comprehensive model. However, extensive simulations were done on the molecular-metallic transition of hydrogen for LMS and planets Chabrier & van Horn 1992). The  model predicts a first order transition for the metallization of hydrogen at a pressure of ∼ 1 Mbar and critical temperature of ∼ 15300 K. Such pressure and temperature values are appropriate for giant planets and brown dwarfs. Modern numerical simulations (Yang et al. 2015;Morales et al. 2010) do confirm the existence of such phase transitions at the same pressure range but predict a much different range of temperature ∼ 2000 K − 3000 K. This new temperature regime is certainly too low for brown dwarfs. Although the pressure estimate is relatively well established in these numerical simulations, the phase transition temperature is still a matter of continuing investigation (Becker et al. 2014). Having noted these caveats, we present the existing model for the surface temperature, based on the  and Burrows and Liebert (1993). We also introduce a simpler treatment of the specific entropy.
The PPT occurs over a narrow range of densities near 1.0 g/cm 3 from a partially ionized phase (x H + ∼ 0.5) to a neutral molecular phase (x H + < 10 −3 ). For massive brown dwarfs the phase transition occurs nearer to the surface. Burrows and Liebert (1993) used the following approximate analytic expressions for the specific entropy (Eqs. (2.48) and (2.49) in Burrows and Liebert (1993)) for the two phases of the PPT: Similar expressions for the entropy at the interior and the atmosphere are given in Eqs. (21) and (22) in Stevenson (1991). We use a simplified approach to make the origin of the above equation clear in the spirit of an analytic model. We derive analytic expressions for the entropy of the ionized and the molecular hydrogen separately for the two phases (similar to equations (2.48) and (2.49) in Burrows and Liebert (1993)) and match them via the phase transition. It is assumed that the presence of helium does not affect the hydrogen PPT ).
The region between the strongly correlated quantum regime and the ideal gas limit can be modelled with corrections to the ideal gas equation. Such correction terms can be expressed by virial coefficients (see Eq. (1) in Becker et al. (2014)). For simplicity, we ignore such corrections in our EOS (Eq. (7)) and consider only the contribution of electron pressure, Eq.(5), and the ion pressure, Eq. (6), of the partially ionized hydrogen (of ionization fraction x H + ) and helium mixture.
The total entropy for our EOS (Eq. (7)) is the sum of the entropies of the atomic/molecular gas and the degenerate electrons. The internal energy per gram for the monoatomic gas particles is U = 3 2 Ideally, we can consider the total energy as a combination of kinetic energy, radiation energy and the ionization energy (Eq. B3). But as the electron gas is degenerate, the radiation pressure is relatively unimportant. Furthermore, as shown in Appendix B, we may even neglect the contribution of the ionization energy as the gas is only partially ionized. Taking the partial derivatives of the expression for the internal energy Eq. (19) we can express the change in heat dQ using the first law of thermodynamics, Eq. (B2). However, as the ionization fraction x H + (ρ, T ) is a function of density and temperature we further use Saha's ionization equations (Eqs. B5, B6) to get . A more generalized version including the radiation and the ionization terms are shown in Eq. (B7) in Appendix B. For T dS = dQ we integrate Eq.
(20) to express the entropy in the interior of the brown dwarf as where and µ 1 is different for each model in this region and is calculated using x H + from Table 1. However, the contributions to entropy due to radiation and the degenerate electrons (see Eq 2-145 in Clayton (1968)) are negligible in the range of temperature and density applicable for brown dwarfs. Based on a similar argument, the analytic expression for the entropy of non-ionized molecular hydrogen and helium mixture at the photosphere is expressed as where 1 µ 2 = X 2 + Y 4 is the mean molecular weight for the hydrogen and helium mixture. The expression for entropy S 2 is derived using the first law of thermodynamics (Appendix B) and the relation of the internal energy of diatomic molecules (U = 5 2 k B N 0 T µ 2 ). Here we have considered only five degrees of freedom as the temperatures are just sufficient to excite the rotational degrees of H 2 but the vibrational degrees remain dormant. It should be noted that Eqs. (21) and (23) are just simplified forms of the entropy expressions presented in Burrows and Liebert (1993) and Stevenson (1991).
Thus the entropy in the two phases is dominated by contributions from the ionic and molecular gas, respectively. Using the same argument as Burrows and Liebert (1993), that the two regions of different temperature and density are separated by a phase transition of order one we can estimate the surface temperature. Using the expression for the degeneracy parameter ψ from Eq. (4) we can simplify Eq. (21) to be Furthermore, the jump of entropy (see Table 1) for the phase transition at each point of the coexistence curve of PPT (Saumon & Chabrier 1989), is used to estimate the relation | C 1 − C 2 | between the two constants in Eqs. (21) and (23). For T = T eff and ρ 2 = ρ e in Eq (23), we can use Eqs. (24) and (23) in Eq. (25) and the value of | C 1 − C 2 | to obtain a wide range of possible values of surface temperature T eff in terms of the degeneracy parameter and photospheric density ρ e : The values of the parameters b 1 and v for different models are shown in Table 1. According to the Chabrier model, the critical temperature ∼ 1.53 × 10 4 K and critical density ∼ 0.35 g cm −3 marks the end of the phase transition with ∆σ = ∆S k B N = 0. In the following discussion we briefly summarise the steps from Burrows and Liebert (1993). We replace Eq. (2.50) in Burrows and Liebert (1993) by Eq. (26) to estimate the surface luminosity. As an example we select a particular phase transition point (model D) and show the derivation of surface luminosity using hydrostatic equilibrium and the ideal gas approximation. The photosphere of a brown dwarf is located at approximately the τ = 2 3 surface, where is the optical depth. Using the general equation for hydrostatic equilibrium, dP = −(GM/r 2 )ρ dr, and Eq. (27), the photospheric pressure can be expressed as where κ R is the Rosseland mean opacity and the other variables have their standard meanings. Furthermore, our EOS (Eq. (7)) in the approximation of negligible degeneracy pressure near the photosphere gives the photoshperic pressure as Now using the expression for radius R (Eq. 13) in Eq. (28), we can calculate the external pressure P e as a function of M and ψ: Similarly, the surface temperature can be evaluated for all the other models. Since the procedure is same for all the models in Table 1, we just show one calculation. For this range of surface temperatures, the Stefan-Boltzmann law, L = 4πR 2 σT 4 eff , yields a set of possible values of the surface luminosity L as a function of the degeneracy parameter ψ. The luminosity for model D using Eq. (13) and Eq. (32) is where σ is the Stefan-Boltzmann constant. Substellar objects below the main sequence mass gradually evolve towards complete degeneracy and a state of stable equilibrium as their luminosity decreases over time. In the following sections we show that the degeneracy parameter ψ is a function of time and it evolves toward ψ = 0 over the lifetime of brown dwarfs. This gives us an estimate of the luminosity at different epochs of time.

Validity of PPT in brown dwarfs
There is a distinction between the temperature-driven PPT with a critical point at ∼ 0.5 Mbar and between 10000 K and 20000 K as predicted by the chemical models (Saumon et al. 1995), and the pressure-driven transition from an insulating molecular liquid to a metallic liquid with a critical point below 2000 K at pressures between 1 and 3 Mbars. The latter is predicted e.g., by Lorenzen Knudson et al. (2015). Figure 1 in Knudson et al. (2015) shows the melting line (black) as well as the different predictions for the coexistence lines for the first order transition (green curves). Brown dwarf interior temperatures are far above these estimates for a first order transition from the insulating to the metallic system. The same is true for Jupiter. Of course a continuous transition may be possible in Jupiter and brown dwarfs, but a first order transition may not be possible. Thus the determination of the range of temperature of this transition provides a much needed benchmark for the theory of the standard models for the internal structures of the gas-giant planets and low mass stars.   . This gives the possible range of surface temperature depending on the phase transition points. For different values of temperature and density at which the phase transition takes place, the effective surface temperature is calculated using Eq. (26).

Nuclear processes
VLM stars and brown dwarfs contract during their evolution due to gravitational collapse. The core temperature increases and the contraction is halted by either the degeneracy pressure of the electrons or the onset of the nuclear burning, whichever comes first. In the first case, the brown dwarf continues to lose energy through radiation and cools down with time without any further compression. However, massive brown dwarfs or stars at the edge of the main sequence can burn hydrogen for a very long time before they either cease nuclear burning or settle into a steady state main sequence. The thermonuclear reactions suitable for the brown dwarfs and VLM stars are p + d → 3 He + γ.
As the central temperature is not high enough to overcome the Coulomb barrier of the 3 He − 3 He reaction and the p-p chain is truncated, 4 He is not produced. Most of the thermonuclear energy is produced from the burning of the primordial deuterium, Eq. (35). The energy generation rates of the above processes are given aṡ ǫ pp = 2.5 × 10 6 ρX 2 T 2/3 6 e −33.8 T 1/3 6 erg/g · s ,  (Fowler et al. 1975). However one can fit the thermonuclear rates to a power law in T and ρ in terms of the central temperature (T c ) and density (ρ c ) as in Fowler and Hoyle (1964): where u ≃ 2.28 and s = 6.31 are constants that depend on the core conditions (Burrows and Liebert 1993). To obtain the luminosity due to the nuclear burning L N = ǫ n dm, we use the power law form for the nuclear burning rateǫ n , (Eq. 38), and making the polytropic approximation ρ = ρ c θ n and setting T Tc = ρ ρc 2/3 , we obtain where r = aζ (Chandrasekhar 1939). Inserting Eqs. (14) and (16)  (40)

Estimate of the minimum mass
In this section we estimate the minimum main sequence mass by comparing the surface luminosity (Eq. 33), with the luminosity (L N ) due to nuclear burning at the core of LMS and brown dwarfs. Instead of just quoting one value as the critical mass, we have presented a range of values depending on the various phase transition points listed in Table 1. This will give us a range of values for the minimum critical mass that is sufficient to ignite hydrogen burning. Model H marks the end of the phase transition and gives a lower limit of the critical mass. However, we calculate the mass limit for model B and model D only. Equating L N of Eq. (40) with L of Eq. (33) gives us where α = 2.32 for µ e = 1.143 and µ 1 = 1.24, which are the number of baryons per electron and the mean molecular weight of neutral (x H + = 0) hydrogen and helium, respectively. These mass densities are evaluated for hydrogen and helium mass fractions of X = 0.75 and Y = 0.25, respectively. The right hand side of Eq. (41) has a minimum at a certain value of ψ. This gives the lowest mass at which Eq. (41) has a solution and this corresponds to the boundary of brown dwarfs and VLM stars. The minimum of F (ψ) is at ψ min = 0.042. Substituting this in Eq. (41) and for κ R = 0.01 cm 2 /g, the minimum mass (model D) is A similar analysis for model B gives the value of minimum mass of M = 0.085M ⊙ for ψ min = 0.042. For the other models the minimum main sequence mass is in the range of 0.064 − 0.087M ⊙ . The solution is relatively independent of mean molecular weight µ 1 . For example, using partially ionised gas i.e. µ 1 = 0.84 in α, it increases the minimum stellar mass by only ∼ 5%.

A cooling model
A simple cooling model for a brown dwarf is presented in both Burrows and Liebert (1993) and Stevenson (1991). In this section we review some of these steps using our more exact EOS Eq. (7) and represent the evolution of the brown dwarfs over time. Using the first and the second law of thermodynamics, the time varying energy equation for a contracting star is expressed as where S is the entropy per unit mass and the other symbols have their standard meaning. The energy generation termǫ is ignored. On integrating over mass we get, where L is the surface luminosity and σ = S k B N A . Now replacing T in terms of the degeneracy parameter ψ in Eq. (4), and using the polytropic relation P = Kρ The variation of the entropy with time (Eq. 44) can be expressed as the rate of change of degeneracy over time. As the star collapses, the gas in the interior becomes more and more degenerate and finally the degeneracy pressure halts further contraction. A completely degenerate star (ψ = 0) becomes static and cools with time. Thus by substituting the time variation of entropy, using Eq. (24), i.e (1 + γ + αψ) 1.715 ψ 4.5797 .
This is a nonlinear differential equation of ψ for the model D, and an exact solution can only be obtained numerically. However we use some very simple and physical approximations to solve this differential equation to yield a simple relation of ψ as a function of time and mass M. As we are trying to estimate the critical mass for the hydrogen burning, it is safe to ignore the early evolution of VLM stars and brown dwarfs. Thus we will solve the differential equation (48) with the assumption ψ ≪ 1. Thus we can drop the term (1 + γ + αψ) 1.715 as it is almost unity in the range 0 < ψ < 0.1 and integrate Eq. (48) in the above limit to obtain Similarly, we can solve for ψ for all other models and obtain the evolution of degeneracy over time. We use this expression of ψ, for model D, and can express luminosity as a function of time t and mass M. The time evolution of luminosity (model D) is represented in Figure 2. It is evident that such low mass objects continue to have low luminosity for millions of years before they gradually start to cool. For t > 10 7 yr, the luminosity declines as a function of time, L ≃ t −1.2 , as shown in Figure 2 (dashed black line). A simplified expression of the variation of luminosity after 10 7 yr for model D is Our luminosity model is consistent with the simulation results of the present day stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA) (see figure 17 in Paxton et al. (2011)). Paxton et al. (2011) use a one-dimensional stellar evolution module, MESA star, to study evolutionary phases of brown dwarfs, pre-main-sequence stars, and LMS.
In Figs. 3 and 4 the ratio L N /L is plotted against time for different masses in the substellar regime for models B and D, respectively. As evident for both the models there is a non-steady state of substantial nuclear burning for millions of years for substellar objects. For a critical mass of 0.085M ⊙ (model B) and 0.078M ⊙ (model D), the ratio L N /L approaches 1 in about a few billion years and marks the beginning of main sequence nuclear burning 2 . Stars with greater mass reach a steady state where the thermal energy balances the gravitational collapse. However, as our model does not consider any feedback from hydrogen fusion, the curves do not stabilize to a steady state main sequence regime, in which L N /L remains 1 until nuclear burning stops. Interestingly, the ratio L N /L is close to unity for many objects below the main sequence transition mass. This suggests that they burn nuclear fuel for a part of their evolutionary cycle but do not have enough mass to sustain a steady state. Note that the results of Becker et al. (2014) discussed in Section 2 would affect the luminosity L by less than a few percent. For example if the radius R increases (or decreases) by 2.5 % for a constant value of mass M, K in Eq. (12) increases by 2.5 % and T eff (Eq. 32) decreases by 1.5 %, therefore L decreases by ∼ 1 %.

Brown Dwarfs as clocks
Interestingly the cooling properties of brown dwarfs (Fig. 2) can be calibrated to serve as an astronomical clock. As the electron degeneracy pressure puts a lower limit to the size of the dwarf, it cools slowly and radiates its internal energy. The luminosity of a brown dwarf is the most directly accessible observable quantity. As luminosity is a time variable, one can get important information on the age of a brown dwarf depending on its mass and the cooling rate. As evident from Figure 2, given mass and the luminosity one can roughly identify the age of the dwarfs. However it is still a challenge to estimate the mass of a brown dwarf. An essential part of the solution is to find brown dwarfs in a binary system where one can get an accurate estimate of the mass and then compare its luminosity against available models. Newly discovered brown dwarfs in eclipsing binaries (Stassun et al. 2006;David et al. 2015) can provide a data set of directly measured mass and radii. This can yield an empirical mass-radius relation that also tests the prediction of the theoretical models. Furthermore, lithium in brown dwarfs has been used as a clock to obtain the ages of young open clusters as originally suggested by Martin et al. (1994) and Basri et al. (1996), and most recently applied to the Pleiades by Dahm (2015). Massive brown dwarfs (M > 0.065M ⊙ ) deplete their lithium on a longer time scale, but VLM stars and objects above the hydrogen burning limit fuse lithium on a much shorter time scale (Magazzu et al. 1993). A limitation of our model is that it does not include rotation. But the mechanical equilibrium in our models may not be significantly affected by this. For example,in model D,using Eq. (49) in Eq. (13) we find the radius of a 10 7 yr old brown dwarf of mass 0.075M ⊙ to be ∼ 8 × 10 9 cm. This implies that for a median observed rotational period (2π/ω) of one day (Scholz et al. 2015) the ratio of magnitude of the rotational energy (∼ MR 2 ω 2 ) to gravitational energy(∼ GM 2 /R) for a 0.075M ⊙ brown dwarf is ∼ 10 −4 . However, convective mixing and the consequent lithium abundance has a strong connection to the rotation rate of brown dwarfs and pre-main-sequence stars (Martin & Claret 1996). Fast rotators are lithium-rich compared to their slow rotating counterparts, indicating a connection between the lithium content and the spin rate of young pre-main-sequence stars and brown dwarfs (Bouvier et al. 2016). Rapid rotation reduces the convective mixing, resulting in a higher lithium abundance in fast rotating pre-main-sequence stars. Thus, the rotational evolution of a brown dwarf can potentially be used as a clock as discussed in Scholz et al. (2015).

Low mass stars that reach the main sequence
In Fig. 5 we plot the time required by low mass stars to reach the main sequence. The curves represent the steady state limit where L N /L = 1 for four different models. It is interesting to note that objects of masses at the critical mass boundary between brown dwarfs and main sequence stars, e.g., 0.078M ⊙ for model D, reach the steady state in about 2.5 Gyr. This suggests that objects just below the critical mass undergo nuclear burning for an extended period of time but fail to enter the main sequence. Furthermore, stars in the mass range 0.078M ⊙ − 0.086M ⊙ for model D take more than 10 8 yr to reach the main sequence. Depending on the phase transition points for different models, these numbers vary but the fact that stars close to the minimum mass limit can take an extended amount of time to reach the main sequence means that young stellar clusters may contain a significant fraction of objects that are still in a phase of decreasing luminosity, and behave like brown dwarfs. These objects will ultimately settle into an extremely low luminosity main sequence. Here, we estimate the fraction of stars that take more than a specified time to reach the main sequence by using the modified lognormal power-law (MLP) probability distribution function of Basu et al. (2015). Their cumulative distribution function is We use the best fit MLP parameters corresponding to the  initial mass function (IMF) as obtained by Basu et al. (2015), where µ 0 = −2.404, σ 0 = 1.044, and α = 1.396. We then use the cumulative function to calculate the fraction of stars taking more than either 10 7 yr, 10 8 yr or 10 9 yr to reach the main sequence. Table 2 contains our results for models A to H. It turns out that about 0.2% of stars take more than 10 9 yr to reach the main sequence and about 4% take longer than 10 8 yr and about 12% take longer than 10 7 yr, for model D. Some of these objects will end up on an extremely low luminosity main sequence, and a sample of luminosity values when an object just above the substellar limit achieves L N /L = 1, i.e., reaches the main sequence, is given in Table 3.  Table 1.  Table 1. b M 9 , M 8 and M 9 are the masses up to which the stars take at least 10 9 yr, 10 8 yr and 10 7 yr, respectively, to reach the main sequence. c N 9 , N 8 and N 7 are the number fraction of stars reaching the main sequence in more than 10 9 yr, 10 8 yr and 10 7 yr, respectively. d Note that for Model A, low mass stars reach main sequence in < 10 9 yr.

Discussion and conclusions
This paper presents a simple analytic model of substellar objects. A focus of the paper was to revisit both the development and the shortcomings of the theoretical understanding of the physics governing the evolution of low mass stars and substellar objects over the last 50 years. We have also made some modifications to the existing models to better explain the physics using analytic forms. Although observational constraints hinder our understanding, a simple analytic model can answer many questions. We have summarized the method of determining the minimum mass for sustained hydrogen burning. Objects in the mass range 0.064M ⊙ − 0.087M ⊙ mark this critical boundary between brown dwarfs and the main sequence stars. We have derived a general equation of state using polylogarithm functions (Tanguay et al. 2010) to obtain the P − ρ relation in the interior of brown dwarfs. The inclusion of the finite temperature correction gives us a much more complete and sophisticated analytic expression of the Fermi pressure (Eq. 3). The application of this relation can extend to other branches of physics, especially for semiconductor and thermoelectric materials (Molli et al. 2011). The estimate of the surface luminosity is a challenge given our limitations in understanding the physics inside such low mass objects. Also, it is still an open question if a phase transition actually occurs in a brown dwarf. The results of modern day simulations (Yang et al. 2015;Morales et al. 2010) do raise doubts about the relevance of phase transitions in the brown dwarf scenario. We are not aware of well defined analytic models that have a unique way of estimating the surface luminosity apart from using the PPT technique as given in Burrows and Liebert (1993). In this work, rather than considering a single value for the phase transition point we have used the entire range of temperatures from the phase transition coexistence table . These are within the uncertainty range of the critical temperature of the PPT as proposed by the recent simulations (Yang et al. 2015). Thus, considering the large uncertainties involved in such models, this range of values of the minimum mass is much more acceptable than a single distinct transition mass. However, the next step forward is to develop an analytic model for surface temperature that is independent of the PPT.
We estimate that ≃ 5% of stars take more than 10 8 yr to reach the main sequence, and ≃ 11% of stars take more than 10 7 yr to reach the main sequence ( Table 2). The stars in these categories have mass very close to the minimum hydrogen burning limit, and will eventually settle into an extremely low luminosity main sequence with L/L ⊙ in the range ≈ 10 −5 − 10 −4 . The very low luminosity non-main-sequence hydrogen burning in substellar objects and the pre-main-sequence nuclear burning in very low mass stars are very interesting to study further, and our simplified model can certainly be improved in its ability to estimate the time evolution.
We thank Dr. Andrea Becker and Dr. Gilles Chabrier for their valuable comments and input during the preparation of this article. We also thank Anushrut Sharma for contributing in the initial stages of this project. Our special thanks go to the anonymous reviewer for a thorough critique and constructive suggestions on the manuscript. SB and SRV were supported by a Discovery Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada. SRV also thanks King's University College for its continued support for his research endeavours.

A. Appendix A: Pressure Integral
The general Fermi Integral can be written as: On substituting x = −β(ǫ − µ) in the second term and x = β(ǫ − µ) in the third term, we arrive at Let us expand the numerator of the second and the third term of the above integral and retain up to the first three terms: µ ± x β n ≃ µ n ± nµ n−1 x β + n(n − 1) 2 µ n−2 x β 2 + . . . .
Substituting this in the integral for the pressure, we can proceed as follows: x 2 dx e x + 1 .
Substituting n = 3 2 for the Fermi pressure (Eq. 2) and evaluating these integrals using the polylogs (Tanguay et al. 2010) we arrive at a simplified form P F ≃ a 2 5 µ terms of linear order in T . On integrating and simplifying the above expression we get the entropy for the partially ionized hydrogen and helium gas, At low temperature (T < 4000K) and low pressure, hydrogen is predominantly molecular and fluid. Repeating the above derivation for the non-ionized diatomic hydrogen gas with energy U = 5 2 k B T N 0 µ 2 and molecular helium we can arrive at an expression for the entropy, In the above expressions for entropy, ρ and T are the density and the temperature, respectively. Other variables are as described in the text.

C. Appendix C: Thermal Properties
We discuss some of the more accurate expressions for the important thermal properties of a degenerate system.

C.1. Fermi energy
Using Eq. (A6) for the number density we can write the general expression for the chemical potential µ in terms of the Fermi energy µ F at T = 0 (Feynman 1972). Considering only the first three terms Eq. (A6) and for ρ = 2 3 µ 3 2 F , we find The second and the third terms are the correction factor C to the zero temperature Fermi energy. For 0.03 < 1 µ F β < 0.20 the correction factor C will be in the range ∼ 8 × 10 −4 < C < 4 × 10 −2 .

C.2. Specific Heat
In the nondegenerate completely ionized limit, the specific heat C v ∼ 3k B N/2. At finite temperatures, the value of the specific heat is less than the limiting value, i.e., C v < 3N k B 2 . The specific heat of the ideal Fermi gas decreases monotonically. At low but finite temperatures, A detailed analysis in the calculation of the specific heat shows that the numerical coefficients in the expansion approached a limiting value of 2.

C.3. Grueneisen Parameter
Applying the condition of constant entropy to Eq. (21) leads to the condition where C is a constant. The Grueneisen parameter γ is given by the expression γ = ∂ log T ∂ log ρ s .
Using Eq. (C3) in the above expression we estimate the value of the Grueneisen parameter γ to be 2 3 . This value is in approximate accord with Stevenson (1991), who indicated that γ ≃ 0.6 in a dense Coulomb plasma when obtained from computer simulations.

C.4. Ionic correlation
Ionic correlation is an important contribution as considered by Stolzmann et al. (1996), Becker et al. (2014), Hubbard et al (1985) and Gericke et al. (2010) to name a few. Stolzmann et al. (1996) use the method of Pade's approximations to provide explicit expressions for the fully ionized plasma of the Helmholtz free energy and the pressure. They have considered the nonideal effects of different correlations such as the electron-electron, ion-ion, electron-ion, as well as exchange contribution for a wide range of values of the Coulomb coupling parameter Γ, which is the ratio of the Coulomb to thermal energy: for ions of different species i. The greatest effect is in the relative pressure contribution P ii /P ideal of the ionic correlation term for hydrogen at T = 10 5 K, estimated to be ∼ − 0.1 at ρ ∼ 10 3 g/cm 3 and a minimum of ∼ − 0.2 at ρ ∼ 10 g/cm 3 .