Constraints on the early luminosity history of the Sun: applications to the Faint Young Sun problem

ABSTRACT

Stellar evolution theory predicts that the Sun was fainter in the past, which can pose difficulties for understanding Earth’s climate history. One proposed solution to this Faint Young Sun (FYS) problem is a more luminous Sun in the past. In this paper, we address the robustness of the solar luminosity history using the yrec code to compute solar models including rotation, magnetized winds, and the associated mass-loss. We present detailed solar models, including their evolutionary history, which are in excellent agreement with solar observables. Consistent with prior standard models, we infer a high solar metal content. We provide predicted X-ray luminosities and rotation histories for usage in climate reconstructions and activity studies. We find that the Sun’s luminosity deviates from the standard solar model trajectory by at most 0.5 per cent during the Archean (corresponding to a radiative forcing of 0.849 W m|$^{-2}$|⁠). The total mass-loss experienced by solar models is modest because of strong feedback between mass and angular momentum loss. We find a maximum mass-loss of |$1.35 \times 10^{-3} \,{\rm M}_\odot$| since birth, at or below the level predicted by empirical estimates. The associated maximum luminosity increase falls well short of the level necessary to solve the FYS problem. We present compilations of paleotemperature and CO|$_2$| reconstructions. One-dimensional ‘inverse’ climate models demonstrate a mismatch between the solar constant needed to reach high temperatures (e.g. 60–80 |$^{\circ }$|C) and the narrow range of plausible solar luminosities determined in this study. Maintaining a temperate Earth, however, is plausible given these conditions.

1 INTRODUCTION

The Faint Young Sun (FYS) problem is a long-standing problem in solar and climate modelling first introduced by Sagan & Mullen (1972). According to classical solar models, the Sun was 20–25 per cent less luminous 3.8–2.5 billion years ago (Ga) during the Archean Eon (Bahcall, Pinsonneault & Basu 2001). Assuming otherwise present-day conditions on Earth, this would result in a surface temperature for the Earth that was too low for liquid oceans in the distant past. This assumption is in disagreement with geological evidence both for liquid water on the Earth’s surface and for warm surface temperatures soon after Earth’s formation 4.568 Ga. The solution to the FYS problem is of interest for understanding the early habitability of the planets in our own Solar system and has implications for that of exoplanets.

The earliest indirect evidence for liquid water on Earth’s surface comes through Hadean zircons (4.4–4.2 Ga; Mojzsis, Harrison & Pidgeon 2001; Wilde et al. 2001) and the earliest direct evidence in pillow basalts and layered banded iron formations (3.8 Ga; Holland 1984; Appel, Rollinson & Touret 2001) as well as fluid inclusions in hydrothermal quartz (3.5 Ga; Foriel et al. 2004). Regarding temperature, isotopic ratios of |$^{18}$|O/|$^{16}$|O imply global mean surface temperatures of 60–80 |$^{\circ }$|C (Knauth & Epstein 1976; Knauth & Lowe 2003; Tartese et al. 2017). These studies and their methods remain controversial, relying on assumptions such as the Archean oceans having a similar |$\delta ^{18}$|O to a modern, ice-free ocean (Hren, Tice & Chamberlain 2009; Blake, Chang & Lepland 2010; Sengupta et al. 2020). In contrast to the high temperatures implied by isotopic ratios, various glacial rocks imply global mean temperatures that fell below 20 |$^{\circ }$|C at times during the Archean, with the Huronian glaciations marking the end of the Archean at 2.5 Ga (Kasting & Catling 2003).

Several solutions to the FYS problem have been proposed (see e.g. Sagan & Mullen 1972; Kasting & Catling 2003; Feulner 2012; Charnay et al. 2020). Increasing the atmospheric abundances of greenhouse gases such as CO|$_2$| or CH|$_4$| is most often invoked as a solution (Feulner 2012). Other proposed ideas include increasing atmospheric pressure to broaden absorption lines (Goldblatt et al. 2009), though there are claims that air pressure was lower during this time (e.g. Som et al. 2016); lowering the Earth’s surface albedo, such as by constraining the emerged land fraction (Jenkins, Marshall & Kuhn 1993; Rosing et al. 2010; Goldblatt & Zahnle 2011b); investigating changes in Archean cloud properties and distribution which could result in stronger warming in climate models (Goldblatt & Zahnle 2011a; Goldblatt, McDonald & McCusker 2021; Yan et al. 2022); and making the early Sun more luminous (Graedel, Sackmann & Boothroyd 1991; Whitmire et al. 1995). In this paper, we focus on the idea that the early Sun could have been more luminous than expected based on standard solar models. We aim to quantify the uncertainties in the Sun’s history due to the input physics in the models, the adopted solar abundances, and the solar mass-loss history. This should be considered more broadly in the context of geological records that constrain atmospheric composition, surface albedo, and surface temperature.

There are sound theoretical reasons to expect the Sun to have been fainter in the past. As the Sun fuses hydrogen into helium in its core, the number of particles decreases and the mean molecular weight increases. In response, the core contracts and reaches higher temperatures in order for the pressure to maintain its balance against gravity. Energy generation is a strong function of temperature, so the Sun’s luminosity increases. For the same reason, more massive stars, which have higher core temperatures, are more luminous. Some authors have argued that the Sun could have been more luminous in the past than expected based on standard solar models if it were more massive and subsequently lost that mass in a magnetized wind. Even in the absence of the solar wind, the mass of the sun will still decrease because the product of hydrogen fusion is less massive than the reactants. The difference in masses is converted to energy via |$E=mc^2$|⁠, and the Sun will eventually radiate this energy via its luminosity.

Mass-loss affects the flux received by the Earth from the Sun in two ways. First, the luminosity of a star is a strong function of the mass (Parker 1965). These are classically related on the main sequence via a power law |$L \propto M^{\alpha }$|⁠, where a typical value |$\alpha \sim 3.5$| is determined through an empirical fit to observed data. The exact relationship is determined via a combination of the energy generation, which increases with mass and temperature, and the opacity, which decreases with temperature allowing for more efficient energy transport. Values of |$\alpha$| from 3 to 5 are often cited in FYS papers. Near the mass of the Sun (⁠|$1.0 \lt M/{\rm M}_\odot \le 1.1$|⁠) the power is locally higher than across the entire main sequence, with |$\alpha \sim 4.5$| from empirically fitting to MIST isochrones (Choi et al. 2016). Secondly, mass-loss can alter the Earth–Sun distance such that a more massive Sun was closer to Earth. Assuming slow (adiabatic) mass-loss, the orbit of the Earth will increase over time such that |$a_{\oplus } \propto M^{-1}$|⁠, and thus the flux received at the Earth will increase as |$M^2$|⁠. The extra power of 2 for the flux results in |$F \propto M^{\alpha +2}$|⁠. Even a slightly more massive Sun could therefore dramatically increase luminosity of the Sun and the flux at Earth. Only a |$\sim$|6 per cent increase in mass would be enough to give the Sun the same luminosity 4.568 Ga, as it has today (Charnay et al. 2020). However, this would require an average mass-loss rate over the main-sequence lifetime of the Sun more than 600 times higher than the current rate.

There are also sound reasons to expect that the Sun experienced a higher mass-loss rate in the past. The activity, mass-loss, and rotation of a star are coupled due to the magnetized solar wind (e.g. Hartmann & Noyes 1987). The stellar dynamo generates the Sun’s magnetic field and hot stellar coronae, which are linked to X-ray emission and mass-loss. Charged particles are carried away from the Sun along magnetic field lines, conserving angular momentum and applying a torque that slows the Sun’s rotation. Skumanich (1972) found that there was a relationship between stars’ rotational velocities v and their ages t such that the velocity decreases with age as |$v \propto t^{-1/2}$|⁠. The Sun’s faster rotation when it was younger corresponds to higher levels of activity and higher X-ray fluxes. Higher X-ray fluxes are correlated to higher mass-loss rates in solar analogues (Wood et al. 2002). The majority of the Sun’s mass-loss occurs early in its lifetime due to its more rapid rotation during this time.

The present-day solar mass-loss rate is small, |$\dot{\rm M}_\odot = 2 \times 10^{-14}\, {\rm M}_\odot \,yr^{-1}$|⁠, so mass-loss has traditionally been neglected in solar models. However, both observations and models predict higher mass-loss rates for young stars. Most papers in the literature find that these mass-loss rates are none the less insufficient to explain the FYS problem. Wood et al. (2002) empirically derived mass-loss rates from solar analogues and found that the Sun’s total mass-loss was less than 0.03 M|$_\odot$| (with a best estimate of 6–7 |$\times 10^{-3}$| M|$_\odot$|⁠) with the majority of the mass-loss occurring within the first few 100 Myr. Johnstone et al. (2015) calculated an expected mass-loss < |$3 \times 10^{-4}\,{\rm M}_\odot$| in a magnetized wind.

Others have argued in favour of greater mass-loss. For example, Martens (2017) defended the plausibility of greater mass-loss for a more luminous Sun. They argued for sustained mass-loss in a stronger wind based on an order of magnitude estimate relating the torque due to mass-loss to the angular momentum of the Sun. However, they did not perform a fully self-consistent calculation including the rotation history of the Sun. In their model, the current solar mass-loss rate also would induce a torque much smaller than that which is observed. Papers that advocate for arbitrary mass-loss to satisfy the FYS problem (e.g. by saying that a few per cent more massive Sun could compensate for the lower luminosity of the young Sun) fail to account for the fact that the Sun’s rotation history will set bounds on the allowed amount of mass-loss with feedback based on the Sun’s angular momentum. This feedback − that high mass-loss rates imply high angular momentum loss rates, which in turn reduce mass-loss − has not had a full, self-consistent treatment in models previously used to address the FYS problem.

There is also uncertainty in the Sun’s mass-loss history because the birth conditions of the Sun − whether it was born as a rapid, moderate, or slow rotator − are still under debate. Solar mass stars arrive on the main sequence with a range of rotation rates, and the majority of them are already slow rotators, with correspondingly modest winds, by the age of the Pleiades cluster (125 Myr; Stauffer, Schultz & Kirkpatrick 1998; Rebull et al. 2016). The origin of the range of rotation rates is a matter of active research, but the standard assumption adopted in angular momentum evolution modelling involves magnetic coupling between protostars and their accretion discs (Koenigl 1991; Shu et al. 1994), enforcing a nearly constant rotation rate during the massive gaseous disc phase. A star with a short disc lifetime will be released from its disc when it is large, and it will become a rapid rotator; a star with a long-lived or massive disc will be released when it is smaller, producing a slow rotator. Although the detailed mechanism is controversial (Matt & Pudritz 2005), something like this is needed to understand the existence of young slow rotators in star clusters. While enhanced lithium depletion in the Sun suggests that it was born as a rapid rotator (Pinsonneault, DePoy & Coffee 2001), another recent paper argued that a rapid rotator was disfavoured to avoid a more active Sun leading to rapid escape of the atmosphere on Earth (Johnstone et al. 2021). A long Jupiter formation time-scale (Kruijer et al. 2017) may suggest a long disc locking time and thus slow or moderate early rotation. Because a rapid rotating Sun has more activity and higher mass-loss than a slow-rotating Sun, the birth rotation period of the Sun is one of the largest uncertainties in the Sun’s early luminosity history.

The luminosity history of the Sun could also be altered by incorporating other new physics such as including star-spots, changing nuclear cross-sections, or changing the composition or opacities used in the models. Their exact impacts may be quantified through perturbing standard solar models. The accuracy with which the luminosity evolution is known changes based on the uncertainties of the input parameters. Past work in the literature has shown that deviant models (e.g. no diffusion) only change the Sun’s luminosity by up to |$\sim$|1 per cent (Bahcall et al. 2001). Bahcall & Serenelli (2005) demonstrated that there were uncertainties in the calculated neutrino fluxes of the Sun due to changing the abundances of certain heavy elements. In addition, other papers have noted that changing the chemical composition of the Sun can lead to disagreements with helioseismic constraints on the Sun’s structure such as the depth of the convective zone or the sound speed (e.g. Asplund, Amarsi & Grevesse 2021).

Conversely, if the solar history is well constrained, we can explore the consequences for models of the ancient Earth. Such reconstructions will depend on both the time history of the Sun and that of the climate on Earth, but the constraints that emerge are none the less powerful.

As a by-product of our efforts to update classical solar models to include effects not traditionally included in solar models, we have developed an improved standard model of the evolution and interior of the sun. Improved standard solar models are of general interest to investigate the effects of refining parameter estimates or expanding the scope of the models to include additional physical effects. For example, models used for neutrino and helioseismology studies have not typically included the effects of rotation, mass-loss, or magnetic fields.

In this work, we first present a non-rotating, classical solar model as a reference model in Section 2. Any subsequent changes to the luminosity are measured relative to this model. In Section 3, we then investigate the effects of changing the solar composition and perturbing the input physics on the solar luminosity history and to what extent this can contribute to solving the FYS problem. We construct new standard solar models including updated compositions, rotation, and for the first time, a fully self-consistent treatment of mass-loss and the angular momentum evolution of the Sun. In Section 4, we use reconstructions of the temperature and atmospheric composition of the Earth to investigate the surface temperatures that can be maintained based on our luminosities and the corresponding solar constant. Finally, we present the characteristics of our new standard solar model including rotation and mass-loss in Section 5, which we believe is the best model of the Sun to date.

2 THE CLASSICAL SOLAR MODEL

In this section, we present our classical standard solar model (e.g. Bahcall & Pinsonneault 1992), which we adopt and henceforth refer to as our reference model. Our reference model follows the traditional structure of solar models computed for studies of helioseismology and neutrinos in order to facilitate comparison with the literature. We begin in the pre-main sequence, at the end of the deuterium burning birthline. This anchors the start of our calculations at the formation epoch of the oldest meteorites. We then evolve models to the age, luminosity, and radius of the present-day Sun. We focus on the time evolution of the bolometric luminosity, the most important component for the FYS problem.

Our reference model and all subsequent variations on this model are calibrated to the current solar age, luminosity, and radius (see Table 1; age of 4.568 Gyr, L|$_\odot = 3.828 \times 10^{33}$| erg s⁻¹, R|$_\odot = 6.957 \times 10^{10}$| cm; Prša et al. 2016). We adopt a solar mass of |$1.988 \times 10^{33}$| g (Prša et al. 2016). We explore the uncertainties in our reference model due to uncertainties in the prior assumptions. In Section 2.2, we will perturb the input physics of the models for comparison with our reference model. In Section 3, we will add physical effects that may affect the time evolution of the Sun, but are not traditionally included in standard solar models. Examples of such effects include rotation, magnetic fields, and mass-loss, and their inclusion allows us to compute our best case standard solar model.

2.1 Reference model

We construct our reference model using the Yale Rotating Evolution Code (yrec; van Saders & Pinsonneault 2012). Our overall approach is similar to that employed by Bahcall & Pinsonneault (1992) and subsequent papers. We adopt our reference case with state-of-the-art physics. We adopt opacities from the Opacity Project (Badnell et al. 2005; Delahaye et al. 2016) with meteoritic metal abundances from Magg et al. (2022). We include gravitational settling and thermal diffusion of helium and heavy elements using the Thoul, Bahcall & Loeb (1994) formalism, Kurucz atmosphere and boundary conditions (Castelli & Kurucz 2003), the OPAL equation of state (Rogers, Swenson & Iglesias 1996; Rogers & Nayfonov 2002), and nuclear reaction cross-sections from Adelberger et al. (2011). We calibrate our models by adjusting the mixing length and surface helium abundance in order to match the present-day luminosity and radius of the Sun (see Tables 7 for adopted values). Changing the mixing length primarily calibrates the radius, whereas the surface helium abundance primarily calibrates the luminosity. We also require our models to match the observed surface metal-to-hydrogen ratio Z/X, which differs from the birth ratio because of gravitational settling and thermal diffusion (Bahcall & Pinsonneault 1992; Bahcall, Pinsonneault & Wasserburg 1995).

Fig. 1 shows the evolution of the solar luminosity scaled to the current solar luminosity as a function of time predicted by our reference models and thus quantifies the magnitude of the FYS problem. The increase in luminosity during the main sequence is the origin of the FYS problem, with the predicted luminosity being 18–25 per cent fainter at the bounds of the Archean (3.8–2.5 Ga). The horizontal dashed line shows a luminosity of 1 L|$_\odot$|⁠, corresponding to the luminosity required for the present-day surface temperature assuming a constant albedo and atmosphere for Earth. Fractional changes in luminosity of |$\delta$|L/L = 0.33 and 0.21 at the beginning and end of the Archean are required to reach this value. In order to be sufficient to solely explain the FYS problem, the luminosity of the Sun would have to be underestimated by as much as 21–33 per cent in our reference model, while still reaching the same present-day value.

Another common set of units used to describe the FYS problem instead uses the mean annual global solar insolation as shown in Fig. 1. This is distinct from the solar constant S, which measures the flux of the Sun at the distance of Earth. The mean annual global solar insolation averages the flux the Earth receives from the Sun over the Earth’s surface area, accounting for the albedo of the Earth. It is expressed as

where L is the luminosity of the Sun, |$a=1$| au is the semimajor axis for the Earth orbiting the Sun, |$R_{e}$| is the radius of the Earth, and |$A_B\sim 0.294$| is the Bond albedo for Earth (Wild et al. 2012). We assume a and |$A_B$| to be constant for Fig. 1. The radiative forcing is the difference in mean annual global insolation for a perturbed model and our reference model. The present-day value of the mean annual global solar insolation is 240 W m|$^{-2}$|⁠. Values of 181 W m|$^{-2}$| (3.8 Ga) and 198 W m|$^{-2}$| (2.5 Ga) corresponding to the bounds of the Archean show that there is a deficit of 42–59 W m|$^{-2}$| for our reference model. Cumulatively, solutions to the FYS problem must contribute at this level to have maintained present-day conditions on Earth. Moving forward, we will report the radiative forcing alongside the fractional change in luminosity for the reader’s convenience.

The evolution of the solar radius and |$T_{\rm eff}$| as a function of time follow a similar trend as the luminosity with a monotonic increase to the present age of the Sun. The present-day effective temperature is |$T_{\rm eff}$| = 5772 K (Prša et al. 2016). |$T_{\rm eff}$| is defined in terms of the radius and luminosity of the Sun assuming it is a blackbody, ensuring consistency between these values.

2.2 Variations on reference model inputs

Intrinsic sources of uncertainty in the predicted luminosity evolution of the sun arise from the uncertainties in the input physics, including the equation of state, opacities of matter to light, energy generation, and energy transport. There are also measured uncertainties in values such as the nuclear reaction cross-sections or the surface metal to hydrogen ratio (often referred to as Z/X, where Z is the mass fraction of all elements heavier than helium and X is the mass fraction of hydrogen). The work in this section investigates the influence that perturbing different inputs to our reference model in Section 2 has on the output luminosity. Each set of models varies an individual input parameter. These models are then compared to our reference model to investigate potential uncertainties in the luminosity of the reference solar model due uncertainties in the input physics. Based on extensive solar model studies, non-linear interactions induced by simultaneous modest changes in input physics are small, justifying our perturbation approach. We assume that changes are uncorrelated, add them in quadrature, and treat systematic errors in input physics as effective 2|$\sigma$| or 3|$\sigma$| uncertainties. See Bahcall & Pinsonneault (1992) for a detailed justification of these assumptions. Impacts on the other aspects of the models (e.g. sound speed) are described in Section 5.

2.2.1 Changing solar compositions

Our reference model in Section 2.1 uses meteoritic metal abundances from Magg et al. (2022, hereafter M22M). In this subsection, we investigate models using the classic metal abundances presented in Grevesse & Sauval (1998, GS98 hereafter) and other recent abundances from Asplund et al. (2021, AAG21 hereafter) and photospheric abundances from Magg et al. (2022, M22P hereafter) (see Table 2). Each set of metal abundances uses its corresponding opacity tables in our models with surface Z/X identified in Table 2. GS98 and AAG21 use a mix of meteoritic and photospheric abundances, while M22M adopts meteoritic abundances where available and M22P adopts photospheric abundances where available. Even for M22M, e.g. photospheric abundances must be adopted for volatile elements like C, N, and O, which are depleted in meteorites.

For almost 20 yr, there has been a persistent disagreement in the literature on the heavy element content of the Sun. GS98 summarized the standard solar composition based on (at the time) updated modelling of new solar spectra. The authors used the empirical Holweger & Mueller (1974) atmospheric structure, and assumed LTE for important species. A series of changes to model inputs led to a different family of solutions, first popularized in Asplund (2005). Hydrodynamic simulations, rather than a fixed microturbulence parameter and a 1D empirical temperature structure, were used; in addition, departures from LTE and improvements in the atomic physics were included. The net result was a substantial reduction in the overall predicted heavy element content of the Sun, as seen in the highly influential Asplund et al. (2009) summary.

AAG21 revised the composition of the solar photosphere using 3D non-LTE calculations for additional elements compared to their previous work (Asplund et al. 2009; see Table 1 of AAG21). AAG21 adopted meteoritic abundances for 15 elements. They argued that the older stellar abundances using 1D LTE calculations like GS98 should be discarded, however, a new generation of independent model atmosphere calculations (Caffau et al. 2008; Caffau et al. 2011) including non-LTE effects and 3D model atmospheres has emerged, and the solar abundances derived from these calculations are systematically higher than to the Asplund et al. (2009) scale. This is in part tied to judgement calls concerning which features to include (see Pinsonneault & Delahaye 2009). M22 used up-to-date 3D non-LTE calculations and simulations, derived from an independent code. M22 argued that many of the differences in elemental abundances between their method and that of AAG21 can be attributed to artefacts or data reduction, particularly for weak features.

Modern stellar atmosphere models have therefore been used to advocate for both a higher (M22M) and lower (AAG21) abundance scale. We thus treat the solar heavy element mixture as a source of systematic error. We adopt M22M as our base case because solar models computed with this heavy element mixture are in excellent agreement with helioseismic data (see our discussion of the solar sound speed profile in Section 5; Fig. 13). By contrast, the lower metallicities in AAG21 lead to a much worse agreement with these data.

Fig. 2 shows the fractional difference in the luminosities as a function of age for the solar models with the adopted metal abundances labelled in the legend. Each luminosity is measured relative to our reference model. The fractional difference in the luminosity of the model adopting the AAG21 abundances increases approximately linearly with time before the present up to |$\sim$|0.4  per cent around |$\sim$|100 Myr. This fractional difference increases rapidly during the pre-main sequence because the luminosity is changing rapidly with time, and thus even small changes in the composition can lead to large fractional differences in the luminosity. Nevertheless, we find that the fractional difference in the luminosity is still always below |$\sim$|1  per cent. In contrast, the fractional changes for the models adopting the GS98 and M22P abundances are |$\lt 0.05$| per cent for essentially all of the Sun’s history, with the GS98 model exhibiting the similar rapid increase prior to |$\sim$|100 Myr. The choice of chemical composition can slightly alter the Sun’s luminosity history. In particular, adopting lower metal abundances can lead to a higher luminosity on the order of a few tenths of a per cent. None the less, with such a small luminosity change compared to that which is required, we conclude that the solar heavy element mixture which we choose to adopt in our models does not affect the Sun’s luminosity at a level necessary to solve the FYS problem.

2.2.2 Perturbing input physics

The work that follows builds upon our reference model, aiming to investigate deviations caused by the uncertainties in the input physics. We vary the individual parameters of our reference model to determine the influence of each parameter on the luminosity history.

Table 2 summarizes the parameter variations and includes references for their values and uncertainties. Historically, there were modest zero-point uncertainties in the solar luminosity (0.4 per cent in Bahcall & Pinsonneault 1992) and age (30 Myr). Over the past decade, the uncertainties in both have declined dramatically. With current data, neither is a major contributor to the error budget.

The effects of diffusion and gravitational settling could potentially have a greater influence on the solar luminosity by, like the metal abundances, changing the interior densities, temperatures, and opacities. The p–p and |$^{14}$|N–p nuclear reaction cross-sections directly influence the luminosity by changing the rate of energy generation. Feedback can cause certain parameters to counterintuitively impact the luminosity. At fixed temperature, a higher p–p cross-section increases the luminosity of a solar model. However, a calibrated solar model has the solar luminosity fixed, so increasing the p–p cross-section has the effect of lowering the model central temperature. On the other hand, the |$^{14}$|N-p cross-section does not have the same counterintuitive impact because the CNO cycle is a small part of the Sun’s energy generation. A higher diffusion rate increases the central metal content relative to the surface value, which is fixed in the calibration process. This leads to a higher core temperature and a greater energy generation rate and luminosity. Helium diffusion has a similar effect.

We compute the partial derivatives of the luminosity with respect to the various input parameters. These are calculated through comparing the result of a variant model with an individual parameter varied by 3|$\sigma$| to our reference model. Note that diffusion is an exception which we do not vary by 3|$\sigma$|⁠, but instead by the 15  per cent 1|$\sigma$| uncertainty. These partial derivatives are a function of time with the partial derivatives equal to 0 at the present day due to the models’ construction with a solar calibration. Table 3 presents these partial derivatives multiplied by the uncertainties, then added in quadrature and converted to the appropriate radiative forcings at 3.8 and 2.5 Ga, the bounds of the Archean Eon when liquid water existed on Earth’s surface.

Fig. 3 illustrates the difference in solar luminosity compared to the standard solar model as a function of age considering variations in the input physics. The top line (blue) is the upper limit on solar luminosity. To obtain these uncertainties, three times the partial derivatives times the uncertainties in the parameters (illustrated by the other lines) are added in quadrature. The individual parameters are identified in the legend. Adding the luminosities of the different models in quadrature assumes that all of these uncertainties are independent and Gaussian, which is undoubtedly an incorrect assumption. We do not know how these uncertainties are distributed. One could make an argument that adding the absolute value of each of the partial derivatives would provide an estimate of the maximum change in luminosity that is possible, but the central limit theorem suggests that our estimate should be approximately correct.

The full distribution of solar luminosities is roughly symmetric about our reference model, so it is possible for uncertainties in the input parameters to lead to a decrease in the inferred luminosity instead. This figure shows a maximum increase (or decrease) in solar luminosity of about 0.25 per cent at early ages, which lies between roughly the same order of magnitude to several times less than the potential change induced by the solar composition. As mentioned above, the central limit theorem suggests that the parameter errors may combine in such a way that the luminosity history of the Sun is not substantially different than our reference model. Nonetheless, the case that the luminosity corresponding to each of these parameters is at a maximum is a useful upper limit on the Sun’s plausible luminosity history, and thus an upper limit on the contribution of solar model inputs to the FYS problem. With a radiative forcing of 0.382–0.299 W m|$^{-2}$| as an upper limit, this is insufficient to solely explain the FYS problem.

3 ANGULAR MOMENTUM AND MASS-LOSS VIA MAGNETIZED WIND

Stars like the sun with thick outer convective zones are thought to generate strong ordered magnetic fields via the dynamo mechanism (e.g. Parker 1955; Kraft 1967). These magnetic fields, coupled with mass-loss due to an ionized stellar wind, cause these stars to spin down over time (e.g. Weber & Davis 1967). Thus, the mass-loss rates of such stars are coupled to their spin down rates. Since a change in the mass of star leads to a change in its luminosity, the luminosity evolution can be tied to the spin down rate. Furthermore, the magnetic fields in sun-like stars leads to X-ray emission, with the strength of this emission coupled to the strength of the magnetic field, which itself depends on the rotation rate of the star. As we have measurements of the rotation rates and X-ray emission of sun-like stars at a wide variety of ages, we can calibrate physical models of mass-loss due to a magnetized wind to empirical data.

Here we model angular momentum loss from birth to the present-day age of the Sun using the Matt et al. (2012) wind model, using the van Saders & Pinsonneault (2013) model (hereafter VP13) to relate the mass-loss rate and the magnetic field strength to stellar observables. We calibrate this model using observations of the rotation rate of sunlike stars in clusters of a variety of ages.

Following VP13, the magnetic field strength scales as

where |$P_{phot}$| is the photospheric pressure, and the Rossby number |$Ro$| is defined as

where |$\omega$| is the angular velocity and |$\tau _{cz}$| is the convective overturn time-scale.

Through the generation of magnetic fields, the stellar dynamo is linked to the X-ray emission in the Sun. There is an empirical relationship between the X-ray luminosity, the bolometric luminosity, and Rossby number (Wright et al. 2013):

where the value of |$\beta$|⁠, which describes the strength of the dynamo, has typically been assumed to be about −2, although some papers have argued for values closer to −3 (e.g. |$\beta = -2.70\pm 0.13$| from Wright et al. 2013). Like X-ray emission, mass-loss occurs due to the magnetic field, which generates a hot, ionized stellar coronae. It is therefore not surprising that empirical relations between mass-loss and X-ray flux have been found in the form of

with values of |$\gamma$| ranging from 0.66 to |$1.34 \pm 0.18$| (Wood et al. 2005; See et al. 2017; Vidotto 2021).

In general, these quantities are close to linearly proportional. We assume mass-loss is linearly proportional to the X-ray luminosity in our models, i.e. the exponent on the right-hand size of equation (5) is equal to 1.

Combining the above equations and assuming the mass-loss to be directly proportional to the X-ray luminosity, the mass-loss rate is given by

Because |$\beta$| describes the strength of the dynamo, it is also the parameter in our models which describes the slope of spin-down versus time. We calibrate |$\beta$| to the median rotation rates of 1 |${\rm M}_\odot$| stars in the Pleiades, M37, Praesepe, and NGC 6811 reported in Godoy-Rivera et al. (2021; see Fig. 4). We expect our answer to lie in the range of |$\beta$| = –2 to –3 (see e.g. the discussion in Wright et al. 2013).

Torque occurs on the Sun due to mass-loss from its magnetized stellar wind. From Matt et al. (2012), the torque is related to the mass-loss rate and dipole magnetic field strength by

where |$B_*$| is the magnetic field strength at the stellar equator, |$\dot{M}_W$| is the total mass-loss rate in the wind, |$R_*$| is the stellar radius, |$\Omega _*$| is the angular rotation rate of the stellar surface, |$v_{esc}$| is the gravitational escape velocity, and |$K_1 \sim 1.30$|⁠, |$K_2 \sim 0.05$|⁠, and m |$\sim$| 0.22 are dimensionless constants which were numerically determined through their simulations. While we could have chosen from different forms to describe the torque, all modern prescriptions are tied to the same empirical constraints on the time history of rotation of solar analogues. As a result, they reduce to similar predicted time dependencies for mean field strength and mass-loss rate.

Substituting the prescriptions for mass-loss rate and the magnetic field, this results in the VP13 parametrization for angular momentum loss in a magnetized wind, which is expressed as

The parameter |$f_K$| is the solar calibration constant and is calibrated in our individual models to ensure that rotation rates are equal to the solar value at present day. The parameter |$K_M$| is given by

where |$c(\omega)$| is the centrifugal correction term which is roughly 1.

As seen in equation (8), the angular momentum loss is saturated above some critical angular velocity (⁠|$\omega _{crit}$|⁠, the saturation threshold) scaled by the Rossby number, given by the ratio of the rotation period and the convective overturn time-scale. This form is consistent with the behaviour of a variety of diagnostics of stellar activity, such as X-ray emission, chromospheric activity indicators, and mean surface magnetic field strengths inferred directly or through star spot filling factors. All saturate in young, rapidly rotating stars. We adopt a saturation threshold of log |$Ro = -0.677$| from Cao & Pinsonneault (2022), which was obtained by fitting to star-spot data. Above the saturation threshold |${\rm d}J/{\rm d}t \propto \omega _{\rm crit}^2 \omega$|⁠, while below this threshold |${\rm d}J/{\rm d}t \propto \omega ^3$|⁠. The effect of this is that initially fast rotating stars can remain fast rotating for longer, which further increases the mass-loss observed in the highest mass-loss models.

3.1 Rotation

The exact birth rotation of our Sun is unknown, so we consider a range of slow, moderate, and rapid rotators in our analysis. For clarity, slow, moderate, or rapid rotators roughly correspond to the 25th, 50th, or 90th percentiles of observed rotation rates, respectively. The youngest observed protostars have massive accretion discs that can exchange angular momentum with the star. Modern angular momentum calculations bypass the complexities of star–disc coupling by using empirical data in systems older than 10 Myr, where massive accretion discs are relatively rare. In this study, we explored a set of birth conditions for the Sun with initial rotation periods corresponding to observed rotation periods of 0.9–1.1 M|$_\odot$| stars in Upper Sco. Models were initialized with a slower rotation rate at an earlier age and allowed to spin up to the desired rotation rate at 10 Myr. While it seems likely that the Sun was a slow or moderate rotator (e.g. only 15 per cent of solar mass stars in the Pleiades are rapidly rotating), we cannot rule out the possibility that the Sun could have been a rapid rotator, which would maximize its early mass and luminosity. We aim to address the question of whether any plausible solar history would be sufficient to explain the FYS problem.

We explore both rotation as a solid body and rotation with partial core envelope coupling, also allowing the time-scale for coupling to vary (see Somers & Pinsonneault 2016 for a detailed description of the models). The mean angular velocity in the surface convection zone is only weakly dependent on depth. We therefore assume solid body rotation in convective regions at all times.

The radiative interior is more complex. From helioseismic data, the solar core is rotating nearly uniformly below a shear layer at the base of the solar convection zone. One limiting case model for internal transport of angular momentum is to assume solid body rotation in radiative regions; this would occur, for example, with strong magnetic coupling. However, there are other models where the spin-down time-scale is longer – for example, time-scales of 10s of Myr are expected from wave-driven angular momentum transport. The spin-down of stars is in better agreement with models including this ‘core-envelope decoupling’, and the required time-scales for the Sun are consistent with the observed flat rotation curve in the interior (Denissenkov et al. 2010; Somers & Pinsonneault 2016; Eggenberger, Buldgen & Salmon 2019). The results shown in our figures use rotational evolution with core envelope coupling. Our solid body models fail to properly reproduce the observed rotation periods of solar mass stars in young open clusters and thus we did not consider them further.

The parameters in our rotational models have been chosen to reproduce the spin-down of the open clusters (Fig. 4). The open cluster data is represented by the black lines. The |$\beta$| value is calibrated to the median rotators (given by the green line) and |$f_K$| is calibrated for each individual model. The starting point of these models is based on rotation data for 0.9–1.1 M|$_\odot$| stars from Upper Sco (Rebull et al. 2018) with a median rotation period of 4.627 d (and a range of 14.33–1.52d for 10th–90th percentiles) at an age of 10 Myr, and all models have the solar rotation period of 25.4 d at 4.568 Gyr. The best-fitting value for |$\beta = -2.29$|⁠. As specified in the previous section, log Ro = −0.677 corresponding to |$\omega _{crit} = 9.79 \omega _\odot$| is adopted from Cao & Pinsonneault (2022).

The Sun is highly depleted in lithium compared to meteorites (Lodders 2021). We calibrate the surface lithium abundance to the solar log relative abundance of A(Li) = 0.96 (Asplund et al. 2021) on the scale where the abundance of hydrogen is 12 via a mixing efficiency parameter. Holding this parameter fixed instead would have resulted in the lowest surface abundance of lithium in the rapid rotating case, and the highest surface lithium abundance in the slowest rotating case, so a range of solar rotation rates in practice translates to a range of efficiencies (see the discussion in Somers & Pinsonneault 2016). The mixing efficiencies |$f_C$| and solar calibration constants |$f_K$| for each model are shown in Table 4. The lithium abundances are shown in Fig. 5. Some lithium is depleted in the early stages of solar evolution, prior to the main sequence. However, solar analogues in the Pleiades have lithium of about 2.6–2.9 at 125 Myr (see fig. 3 of Bouvier et al. 2018), implying 1.6 dex, or a factor of 40, depletion on the main sequence. Standard solar models predict little or no main sequence depletion, and the most promising explanation is rotationally induced mixing (Pinsonneault 1997). There is a drastic difference between the non-rotating and rotating models that highlights the need for additional mixing processes to reproduce the present solar lithium abundance. Our reference and slowest rotating models show lithium depletion close to those seen in the Pleiades at the corresponding age, while our fastest rotating models are deficient in lithium.

With these calibrations in mind, we move towards computing standard solar models which adopt parameters that are chosen to reproduce the behaviour of solar analogues. Our new standard solar models adopt much of the same physics as our reference model, with a few key exceptions. Our standard solar models once again use yrec and this time include a full, self-consistent treatment of rotation, magnetism, and mass-loss based on the physics presented in Section 3. Now that we consider rotation histories in our yrec models, we must compute a suite of models from slow to rapid rotators to capture all plausible solar histories.

The total angular momentum of the Sun is set by a combination of an assumed birth rotation rate and a disc coupling time-scale; we can therefore achieve any given target rotation rate with a combination of the two. We adopt a fixed disc lifetime (0.1 Myr) and a range of birth rotation rates to reproduce the observed range in Upper Sco, our chosen system for initial conditions. We do so to avoid numerical problems in our most rapidly rotating cases. Our results for a range of birth rates and for a range of disc lifetimes are similar at later ages, so this choice does not affect our conclusions.

Fig. 6 shows the fractional difference luminosities for the 10th–90th percentile rotators with a constant mass of 1|${\rm M}_\odot$|⁠. Through much of the Sun’s history, the luminosity of the Sun is consistent with non-rotating models. The slowest rotation periods (10th percentile) have the greatest change in luminosity at |$\sim$|100 Myr at roughly a tenth of a per cent. There are two key competing effects which dominate at different times. First, rotation reduces the efficacy of diffusion on the main sequence, effectively making the star behave as if it had a lower metallicity and resulting in a hotter, more luminous star. Secondly, the stars are rotating relatively rapidly on the pre-main sequence and the centrifugal force from rotation partially counteracts the star’s gravity. The star behaves like a slightly lower mass star, maintaining hydrostatic equilibrium at lower core temperatures, and is less luminous.

3.2 Mass-loss

We investigate the effect of mass-loss on the Sun’s luminosity using a tracer calculation. Determining the mass of the Sun as a function of time from our rotational evolution calculation allows us to interpolate the luminosity as a function of time from stellar models. We construct two non-rotating stellar models, the first with a mass of 1 |${\rm M}_\odot$| and the second with a mass of |$1\, {\rm M}_\odot + {\delta }M$|⁠, where |${\delta }M = 1.35 \times 10^{-3} \,{\rm M}_\odot$| is the maximum total mass-loss in any of the rotational models. Using the mass as a function of time to compute the linear interpolation of the two models’ luminosity tracks, the resulting luminosity as a function of age is shown in Fig. 7. The colors in this plot each show a luminosity history corresponding to a different rotational model. Models that were initially faster rotating lost the most mass and show a greater increase in luminosity. Additional mass-loss occurs as fusion converts some mass to energy, |$m = E/c^2$|⁠, which is lost via the Sun’s luminosity. Due to a technical detail of yrec, it is difficult to incorporate mass-loss occurring in the core of the Sun into the full stellar models, but it is included in our tracer calculation in Fig. 7.

We observe an increase in luminosity of less than 0.5 per cent for the highest mass-loss (fastest rotating) models. This is comparable to the difference between the AAG21 and M22M luminosities in Fig. 2 at very early ages and less than that during the Archean. Moreover, even in the lowest mass-loss (slowest rotating) models, the fractional change in luminosity is comparable to what is observed in Fig. 3, demonstrating that realistic mass-loss has a similar impact on the Sun’s early luminosity to uncertainties in the input physics and a similar or smaller impact than uncertainties in the solar mixture. There are opposite trends in Figs 6 and 7. Because the slowest rotator has the highest luminosity in the rotation-only case, it counteracts that the slowest rotator has the lowest luminosity in the mass-loss-only case. Effects of rotation dominate and the slowest rotator is the most luminous as a result. All cases including rotation and mass-loss are more luminous than the reference model during the Archean.

Johnstone et al. (2015) also used rotational evolution models to estimate the mass-loss in the stellar wind of the Sun. They used solid body models and observations of stellar rotation rates in young clusters to fit the mass-loss rate as a function of radius, rotation rate, and mass. They found a dependence of mass-loss |$\dot{M}\propto \omega ^{1.33}$| compared to our effective |$\dot{M}\propto \omega ^{2.29}$|⁠. The early mass-loss rate from our calculation is roughly 3–5× higher at 300 Myr (see fig. 10 of Johnstone et al. 2015) and even greater at earlier times. Their integrated mass-losses are less than |$3 \times 10^{-4}\,{\rm M}_\odot$| even in their fastest rotating models, 4× smaller than we find in this paper. While they integrated from 100 Myr to 5 Gyr, which is a different time-frame than our work, assuming a reasonable saturation rate of mass-loss for the first 100 Myr still leaves their total mass-loss well short of our result. It is interesting to note that because of the feedback between mass-loss and angular momentum loss, we arrived at estimates of total mass-loss within a factor of 3, despite our having adopted much higher mass-loss rates for rapid rotators.

Another notable difference is that they adopted solid body models, while we adopted core-envelope decoupling. Core-envelope decoupling matters less after the first |$\sim$|100 Myr because the evolution of the star’s rotation becomes long compared to the coupling time-scale (⁠|$\sim$|20 Myr). This explains why they found good agreement between their observations and solid body models, while our solid body models did not reproduce the cluster rotation data from Godoy-Rivera et al. (2021). However, both our calculation and theirs set strict limits on mass-loss that are far less than those required to substantively impact the FYS problem.

3.3 X-rays

The X-ray luminosities of solar-type stars as a function of age have been measured (e.g. Pizzolato et al. 2003). The model that we adopt is consistent with the trends seen in the data, with higher X-ray luminosities and the saturation of X-rays among the most rapidly rotating stars. In addition to reporting other properties, we can therefore use our models to trace the energetic photon environment for the Sun in the past, and we use |$L_x/L_{bol}$| as a convenient proxy. The X-ray luminosity of the Sun (and in particular the high-energy photons produced) could have a substantial impact on the atmospheres of the early planets through atmospheric escape or photodissociation, changing the early atmospheric chemistry of the planets.

In our rotational evolution models, the X-ray luminosity is calculated using equation (4). The chosen solar parameters are |$L_{X,\odot } = 10^{27.6}$| erg s|$^{-1}$|⁠, |$\omega _\odot = 25.4$| d, and |$\tau _{cz,\odot }=12.3$|d. There is uncertainty in the choice of |$L_{X,\odot }$| and |$\omega _\odot$|⁠, while |$\tau _{cz,\odot }$| is a parameter computed via models. The X-ray luminosity of the Sun varies over the solar cycle between |$L_{X,\odot } = 10^{26.8}$| and |$L_{X,\odot } = 10^{27.9}$| (Judge, Solomon & Ayres 2003). While several choices for the solar X-ray luminosity could be adopted, we chose the arithmetic mean of these two values. The rotation period of the Sun varies based on latitude from 24.5 d at the equator to 33.4 d near the poles. We chose the Carrington rotation rate of 25.4 d, which is the rotation period observed for low-latitude sunspots. Fig. 8 shows the ratio of the X-ray to bolometric luminosities in our models. The proportion of X-rays is highest in the fastest rotating models with about an order of magnitude in spread, indicating substantial uncertainty for the young Sun based on its rotation history.

Tu et al. (2015) used a rotational evolution model to calculate the X-ray emission as a function of time for solar-like stars. They adopt a range of possible rotation rates based on observations of h Per and NGC 6530, a slightly lower Ro = 0.13 (log R_o = −0.89) as a saturation threshold, and a present day |$L_X = 10^{27.2}$| erg s|$^{-1}$|⁠. Their X-ray luminosity histories are computed using methods that are similar to ours and are consistent with our estimates.

4 IMPLICATIONS FOR CLIMATE MODELLING AND PALEOTEMPERATURE RECONSTRUCTIONS

Earth’s temperature evolution throughout the Archean and Proterozoic remains a topic of debate in the literature. A compilation of paleotemperature reconstructions is given in Fig. 9.

Geochemical temperature reconstructions based on |$\delta ^{18}$|O (Knauth & Epstein 1976; Knauth & Lowe 1978) and |$\delta ^{30}$|Si (Robert & Chaussidon 2006) isotopic data from marine charts have been interpreted as indicating temperatures > 30 |$^{\circ }$|C throughout most of the Proterozoic, even exceeding 60 |$^{\circ }$|C in the Archean and early Proterozoic. Various studies have questioned these result, citing e.g. the possibility of hydrothermal alterations, diagenetic processes, or challenging the central assumption that |$\delta ^{18}$|O in seawater did not change over time (van den Boorn et al. 2007; Johnson & Wing 2020; Sengupta et al. 2020). A newer analysis of charts by Hren et al. (2009) uses |$\delta ^{18}$|O in conjunction with hydrogen isotopes (⁠|$\delta$|D) and estimates ocean temperatures <40|$^{\circ }$|C. Likewise, |$\delta ^{18}$|O data derived from phosphates in sediments point to a temperate Archean with temperatures of 26–35 |$^{\circ }$|C at 3.2–3.5 Ga Blake et al. (2010).

Despite the above criticisms and alternative analyses, interpretations of oxygen isotope data (⁠|$\delta ^{18}$|O, |$\delta ^{17}$|O) as indicative of a very hot Archean |$\sim$| 66–100 |$^{\circ }$|C have persisted in the literature (Tartese et al. 2017; Lowe et al. 2020; McGunnigle et al. 2022). These claims of moderate to high paleotemperatures based on isotope ratios are in conflict with evidence for glaciations occurring throughout the Archean and extending to the Paleoproterozoic (e.g. Young et al. 1998; Williams 2005; Kuipers et al. 2013; Ojakangas et al. 2014; de Wit & Furnes 2016). Later in Earth’s history, a dearth of evidence for glacial deposits during the ‘boring billion’ (ca. 1.8–0.8 Ga) suggests a prolonged period of elevated sea surface temperatures exceeding >20 |$^\circ$|C. To explain why the modern length of day (lod) is 24 h, Wu et al. (2023) propose surface temperatures averaging 40–45 |$^{\circ }$|C during the boring billion, creating a resonance between the solar driven thermal tide and an atmospheric oscillation (Zahnle & Walker 1987). The resulting increased thermal torque would have effectively balanced the lunar tidal torque that would have otherwise lengthened the day to >65 h, instead leading to a nearly constant lod for an extended time (lod lock). However, geological evidence for a lod lock during the Precambrian is highly uncertain (Laskar et al. 2024) and the tidal evolution of the Earth–Moon system may be explained without a resonant atmospheric tidal lock (Farhat et al. 2022). Farhat et al. (2024) find that the atmospheric resonance should have occurred later – in the Phanerozoic – and that its amplitude is not high enough to cause a lod lock.

In recent years, a biochemical line of evidence has been invoked to bolster claims of gradual cooling of Earth’s surface through time. Based on thermostability measurement of experimentally reconstructed ancestral enzymes, Garcia et al. (2017) report a decline in environmental temperature of photic zone organisms from |$\sim$| 65 to 80 |$^{\circ }$|C in the Archean to |$\sim$| 20–60 |$^{\circ }$|C in the Neoproterozoic. Similarly, Harada et al. (2019) report environmental temperatures |$\sim$|40 |$^{\circ }$|C during the Proterozoic for marine phototrophic and planktonic organisms.

Some paleotemperature reconstructions from Fig. 9 may be hard to reconcile with other constraints to Archean or Proterozoic Earth’s properties proposed in the literature. According to nitrogen and argon isotope analysis of fluid inclusions in Archean hydrothermal quartz, the partial pressure of N|$_{2}$| was <1.1 bar and may have been as low as 0.5 bar at 3.0–3.5 Ga (Marty et al. 2013). With the same method, pN|$_2$| was constrained to similar or lower than modern at 3.3 Ga (Avice et al. 2018). Moreover, Som et al. (2012) limit air pressure at 2.7 Ga to less than 2 bar by studying fossilized raindrop imprints. The size distribution of gas bubbles in basaltic lava flows may even indicate an upper limit of half of modern air pressure at 2.7 Ga (Som et al. 2016).

Past CO|$_2$| concentrations may have reached much higher values than the preindustrial atmospheric CO|$_2$| level (PAL) of 280 ppm as shown in Fig. 10. Atmospheric escape modelling implies that the Archean atmosphere at 3.8 Ga may have required a minimum of 40 per cent CO|$_2$| to prevent rapid atmospheric loss (Johnstone et al. 2021). Hessler et al. (2004) derive a minimum of |$2.51\times 10^{-3}$| bar CO|$_2$| at 3.2 Ga by evaluating weathering rinds of river gravels. 2.7 Ga old fossil micrometeorites exhibiting oxidation (Tomkins et al. 2016) have spurred a number of studies trying to constrain pCO|$_2$| with inferred lower abundance limits ranging from 23 to 70 per cent (Lehmer et al. 2020; Payne, Brownlee & Kasting 2020; Huang et al. 2021), and some authors additionally considering low pN|$_2$| (Rimmer, Shorttle & Rugheimer 2019; Payne et al. 2020) as explanation. Sheldon (2006) developed a weathering mass-balance approach to estimate pCO|$_2$| from paleosols that was subsequently used in other studies (Driese et al. 2011; Sheldon 2013; Medaris et al. 2017). Estimates from this method reach up to |$\sim$| 40–45 PAL in the Archean and Paleoproterozoic. A different method by Kanzaki & Murakami (2015) relies on the cation concentrations in porewaters at the time of weathering of paleosols to estimate pCO|$_2$| between 2.77 and 1.85 Ga. Their estimates imply a gradual decline in pCO|$_2$| during the Neoarchean to Paleoproterozoic period and exceed the values calculated by Sheldon (2006)’s method roughly up to one order of magnitude. Consistent with the low CO|$_2$| values found by Sheldon (2013) for the late Mesoproterozoic, cyanobacteria microfossils finds exhibiting calcification – a process thought to only occur below the threshold of |$\sim$| 13 PAL – further only permit moderately elevated CO|$_2$| concentrations for 1.2 Ga (Kah & Riding 2007) and 0.4 Ga (Liu et al. 2020).

In Fig. 11, we present a set of simple 1D climate model calculations (see Appendix  A for a model description). The ‘inverse’ climate model (originally by Kasting et al., used in e.g. Kopparapu et al. 2013) receives a parametrized temperature–pressure profile and a set of greenhouse gas concentrations as input and calculates the solar constant S needed to sustain the given surface temperature. We calculate required solar constant values to reach various target surface temperatures for the CO|$_2$| constraints sets from Fig. 10 using a fixed pN|$_2$| value of 1 bar. We additionally apply a conservative CH|$_4$| amount of 20 mbar for all calculations. This value is consistent with the lower limit for Archean sulphur mass-independent fractionation (Zahnle, Claire & Catling 2006) and roughly corresponds to predictions by biogeochemical box modelling coupled to photochemistry (Claire, Catling & Zahnle 2006).

Our results show that solar luminosity is high enough to maintain temperate conditions above freezing or similar to modern Earth’s 288 K, even with our conservative assumptions on greenhouse warming by CH|$_4$|⁠. However, we struggle to reconcile CO|$_2$| constraints and conservative CH|$_4$| assumptions with very high Archean temperatures or high Proterozoic temperatures. This is not unexpected. For example, Charnay et al. (2017)’s GCM calculations need 1 bar CO₂ and 2 |$\times 10^{-3}$| bar CH|$_4$| with 1 bar N|$_2$| to reach 68 |$^{\circ }$|C at 3.8 Ga. To obtain a 40–45 |$^{\circ }$|C temperature average over the boring billion, Wu et al. (2023) used a surface pressure of 2 bar in conjunction with 100–200 mbar of CO|$_2$|⁠, vastly exceeding Sheldon (2006)’s range (a factor of |$\sim$|500) and the upper limit set by Kah & Riding (2007) (a factor of |$\sim$| 40).

Our 1D climate calculations demonstrate a clear mismatch between the solar constant values needed to reach high temperatures >|$\sim$|30 |$^{\circ }$|C and the narrow range of possible solar luminosity determined in this study. In contrast, the difference in required luminosities implied by the errors given in CO|$_2$| reconstructions alone, e.g. by Kanzaki & Murakami (2015) may reach 10 per cent S. Any attempt of reconciling very high temperature reconstructions, a low-pressure atmosphere and literature CO|$_2$| reconstructions via the minimal uncertainties in solar luminosity our modelling predicts seems unviable. That said, maintaining cooler-but-still clement conditions that would avoid global glaciation (⁠|$\sim$|0 |$^{\circ }$|C < T < |$\sim$|30 |$^{\circ }$|C) seems plausible given the solar constants in this work, and our simple climate model simulations.

5 THE NEW STANDARD SOLAR MODEL

In Section 3, we demonstrated that the impact of non-standard physics on classical solar models like our reference model is generally a small perturbation over the time evolution of the models. None the less, there are certain properties of the Sun that can only be reproduced via models which include a more complete understanding of physics. For example, we showed that our work agrees with prior results that rotational mixing is necessary to reproducing the lithium depletion in the Sun. In this section, we describe the results from our rotating yrec models, which includes rotation and mixing in a set of standard solar models for the first time. The inclusion of rotation and magnetism in our models also allows us to also include calculations of the X-ray luminosity as a proxy for solar activity, making the outputs of this set of standard solar models useful to a broader set of studies. Our results are expanded to a standard solar evolutionary history rather than simply a present-day standard solar model. Classical solar models match today’s properties of the Sun, but our models are also constrained by relevant observations – rotation periods in particular. The slope of spin-down has been calibrated to data from open clusters, the initial rotation rates of our models match stars from Upper Sco at 10 Myr, the saturation threshold is obtained via star-spot data, and the mixing efficiency is calibrated to match surface lithium.

We will highlight a median rotator as the best singular result for the Sun’s evolutionary history, citing the range from the 10th to 90th percentiles of rotation rates as our uncertainty. Because it was used for the calibration, the median rotation in Fig. 4 most closely matches medians in cluster data (compared to slow or rapid rotators and their matching percentiles). It also naturally lies in the intermediate region between the slow and rapid rotators. We again caution that the Sun’s past rotation history is uncertain and not well constrained, so the reader may wish to consider whether the full range of allowed rotation is relevant to their studies. Regardless of which initial rotation rate we choose, the overall picture of the Sun’s evolution will be similar because stellar rotation quickly converges on the main sequence within a few 100 Myr.

Tables 5, 6, and 7 summarize the relevant results from all of our solar models and selected results from the literature. Results from yrec models adopting the M22M abundances and including rotation are grouped together, separated by a horizontal line. The values are the percentiles of rotation rates, with the 90th percentile corresponding to rapid rotators. Results from non-rotating models adopting various metal abundances are grouped together, separated by another horizontal line from results from the literature.

Table 5 describes the neutrino fluxes in our models. Neutrino fluxes in this work are consistent with those found in the literature, with the biggest discrepancy generally being the non-rotating model adopting AAG21 metal abundances and mostly showing slightly lower (but still consistent) neutrino fluxes. Additionally, any rotation slightly decreases the neutrino fluxes for heavier elements since they do not settle as much.

Table 6 describes the properties at the centre of the Sun. Our central helium abundance is lower than in the literature with a fractional change of |$\sim 0.01{\!-\!}0.02$|⁠. Stellar interiors models use a Taylor series expansion in the deep core, and the precise central abundance depends on the inner fitting point. The central temperatures and density, by contrast, are integrated to the centre and insensitive to the fitting point location. As a consequence, there is slight disagreement for the central abundances since they are reported from the inner fitting point.

Table 7 describes the properties at the surface of the Sun and the base of the convective zone. The rotation and no rotation cases have generally good agreement with the surface helium (⁠|$Y_s \sim 0.248$|⁠) and the radius of the base of the convective zone (0.7132 R|$_\odot$|⁠). The fractional changes in sound speed, |$\lt \delta c/c\gt $|⁠, are similar in each case. Only the cases with rotation reproduce the solar lithium abundance. When higher mixing efficiencies are needed to match the lithium (see Table 4), it also increases the surface helium abundance. Helium will otherwise undergo settling in the Sun and a greater mixing efficiency returns more helium to the surface.

Fig. 12 shows the present day sound speed profile of the Sun in our reference model. It peaks above 500 km s⁻¹ in the core and decreases to less than 100 km s⁻¹ in the outer envelope. The sound speed profiles can be used as one test to help identify which solar compositions best match solar observations, so the sound speed presented here serves as a baseline for reference in later plots. Sound speed difference plots are measured relative to the observed sound speed profile of the Sun. Fig. 12 differs from the observed sound speed profile at roughly the tenth of a per cent level.

As in the literature, we observe a difference in the sound speed profiles of our models. Figs 13 (no rotation) and 14 (rotation) compare the sound speeds in the various models to the sound speed as presented in Basu, Pinsonneault & Bahcall (2000). The observed sound speed profile of the Sun is determined via helioseismic inversions using solar frequencies from the Michelson Doppler Imager (MDI) data set. A change of 0 is shown by the dotted line. In Fig. 13, the agreement between our models and observations is best for GS98 and worst for AAG21, with the difference emphasized near the base of the convective zone at a radius R = 0.7132 R|$_\odot$| (also see |$\lt \delta c/c\gt $|⁠, the root-mean-square of the fractional difference in sound speed, in Table 7). This result is in good agreement with results presented in other works (e.g. Magg et al. 2022, others). The sound speed agreement favours higher metallicity compositions and factors into our selection of the M22M abundances for our standard solar models.

Fig. 14 shows the effects of rotation on the sound speed profile. The rotating models show less severe peaks and valleys. Near the base of the convection zone, this occurs because rotational mixing smooths the effects of diffusion. If you change the helium content there, it must change elsewhere, and we see an increase in |$\delta c /c$| in the deeper layers of the Sun. Changes in the sound speed profile are correlated due to mass conservation and the dependence of the sound speed on density (⁠|$c_s^2 = {\rm d}P/{\rm d}\rho$|⁠). In both the rotating and the non-rotating cases, the disagreement between the observed and modelled sound speeds increases near the surface of the Sun. This is due to the choice of radius to which our models are calibrated, specifically that we adopted the radius from Prša et al. (2016) (a smaller radius shows better agreement with the observed sound speed profile).

6 DISCUSSION

The FYS problem could be viewed as a problem for our understanding or the Sun, or a problem for our understanding of the Earth. When viewed through the lens of stellar evolution theory, the evolution of the Sun is driven by fundamental physical properties, and there are powerful reasons to expect solar evolution to be insensitive to changes in standard input physics. The solar calibration itself is an important consideration: solar models must reproduce the current solar properties, dramatically reducing the dynamic range of possible model trajectories. In addition, important current solar global properties – mass, age, radius and luminosity – are now exceptionally well measured. We confirm that the level that the uncertainties in solar modelling contributes to the solution is small. At most, perturbing the input physics to our models show a maximum variation in luminosity of 0.2 per cent. Recent controversies about the solar heavy element mixture modestly impact the luminosity at a given time by up to 0.3 per cent. Neither effect is substantive enough to impact the FYS problem.

However, there have been models proposed where physical processes neglected in the standard framework, in particular mass-loss, could lead to very different histories. We have therefore investigated mass-loss in the context of a physically self-consistent model that also incorporates rotation. We find that mass-loss and rotation can only induce small changes, at the 0.3 per cent level in L. The underlying reason is that the current solar mass-loss rate is negligible, so the only avenue for a high solar mass-loss rate is to have the Sun experience much stronger winds in the past. Wind strength is, indeed, higher in rapid rotators; however, most solar analogues rotate slowly for almost their entire main sequence lifetime, and virtually all rotate only |$\sim 3$| times the solar rate at roughly fifteen percent of the solar age. Any self-consistent solar wind model is thus constrained to minimal main sequence mass-loss, severely limiting the impact on the history of the Sun. Independent calculations by Johnstone et al. (2015) found even smaller effects. Mass-loss rates derived from full evolutionary models are therefore below empirical estimates (a best estimate of |$6{\!-\!}7 \times 10^{-3}\, {\rm M}_\odot$|⁠; Wood et al. 2002). Combining all of the effect that we identified these leads to an increase in luminosity at any given time at the percent level or less. We conclude that the Sun does not cause the FYS problem, and the solution should be found closer to home.

Fortunately, we can now build terrestrial climate models with a robust and well predicted solar luminosity evolution. Incorporating paleotemperature reconstructions and predictions for greenhouse gas concentrations into climate models demonstrates that is reasonable for our solar luminosity history to reproduce temperate conditions on Earth. However, given the small uncertainties in luminosity, such models fall well short of reproducing the higher temperature reconstructions during the Archean.

We find that the most substantial increase in luminosity in our models occurs before the Archean. The mass-loss that occurs according to a prescription for mass-loss via a magnetized wind primarily occurs early in the life of the Sun (see Fig. 7). Additionally, the compositional uncertainties and input physics uncertainties show the potential for their greatest changes early in the Sun’s history (Figs 2 and 3). The potential impact of this result is greater for Mars than it is for Earth. As Mars had conditions for liquid water at an even earlier age than Earth, the FYS problem is more severe for the red planet and contributions from any source may be significant.

One avenue through which the Sun could lose additional mass that we did not explore in this paper is coronal mass ejections (CMEs). Mass-loss due to CMEs may be greater than mass-loss in the steady solar wind (see Cranmer 2017). We caution that any mass-loss due to CMEs will, however, result in angular momentum loss and spin-down, lowering the mass-loss rate due to the solar wind. There is a finite angular momentum for young stars and as such, a finite amount of angular momentum that can be lost to arrive at the present-day Sun. The interplay between CMEs and angular momentum loss may make the total mass-loss comparable to our models which only consider angular momentum loss due to a wind.

Given that a changing mass for the Sun changes its orbital separation with the planets, there may be observable signatures on Earth or Mars that could reveal the Sun’s early history if its mass-loss is large enough. Spalding, Fischer & Laughlin (2018) found that measuring Milankovitch frequencies could test the Sun’s mass to a precision of |$\sim$|1  per cent using a link between the banding of sediments and the orbital properties of Earth or Mars. We note that the |$\delta$|M we calculated is much less than the precision of such a measurement. None the less, it would be useful to learn whether a more massive young Sun is observationally disfavoured.

The stellar spectrum is an important ingredient for the habitability of planets. In particular, the high-energy UV and X-rays dictates much about the atmospheres. The high-energy flux can heat and inflate the atmosphere, leading to increased rates of evaporation (e.g. Wordsworth & Pierrehumbert 2013). In addition, it impacts atmospheric chemistry, photodissociating molecules or driving chemical reactions, and groups have derived estimates of the early solar spectrum (e.g. Claire et al. 2012). Knowledge of the UV is crucial for the interpretation of possible biosignatures. The concentrations to which gases such as O|$_2$| and CH|$_4$| build up in an atmosphere are tied to their destruction rates through photochemistry as well as their production rates (Meadows et al. 2018).

Star spots could further complicate the stellar spectrum. The fraction of a star’s surface that is covered by star-spots is higher for rapidly rotating stars and the radius of these stars is also inflated. Cao & Pinsonneault (2022) used a two-temperature model to describe the appearance of these stars. Studying the effect this has on stellar spectra could be an interesting direction for future study. We note that star-spot filling factors are already low for solar mass stars in young open clusters, however, so the effect will likely be modest. Finally, a direct model including mass-loss could capture non-linear effects that could be of interest in precise solar models.

7 CONCLUSIONS

The FYS problem is a long-standing problem in astronomy. Several solutions have been proposed. While the most likely solution is a combination of different effects that must balance surface temperatures, atmospheric compositions, etc. with the solar luminosity, this paper aims to address the extent to which uncertainties in solar modelling can contribute with a particular interest in the mass-loss history of the Sun.

We use yrec stellar models and perturb the input physics to investigate the effect on the evolutionary tracks. The result is that the uncertainty is small – less than 0.2  per cent (0.382–0.299 W m|$^{-2}$| from 3.8–2.5 Ga) even at the beginning of the Archean when the uncertainty is greatest – which is far less than the 21–33  per cent increase in luminosity (59–42 W m|$^{-2}$|⁠) that would be required for solar modelling to solely solve the FYS problem. We also investigated the effects of composition, finding a slightly larger increase in luminosity of 0.3  per cent (0.596–0.421 W m|$^{-2}$|⁠) for models adopting the AAG21 metal abundances.

We then consider rotation and mass-loss. Stars are born with a range of rotation rates, which are mostly slow even at young ages because of interactions between the star and its disc. The feedback between mass-loss and angular momentum leads to spin down on the main sequence, with rotation rates converging even for relatively young open clusters (e.g. Praesepe, 700 Myr). For slow rotators, Skumanich implies |$\dot{J} \propto \omega ^3$|⁠, with wind loss rates scaling with rotation at comparable rates (see Section 3). This saturates at a level around 10x solar, and the saturated domain corresponds to high mass-loss. However, the duration of the saturated phase is brief for solar analogues, and mass-loss rates rapidly approach the low solar value once it has ended.

Previous standard solar models assumed a constant mass for the star over the course of its evolution. We find that including mass-loss can increase the luminosity of the Sun up to 0.3 per cent (0.468–0.335 W m|$^{-2}$|⁠). The magnitude of this increase is dependent on the birth conditions of the Sun’s rotation, with rapid early rotation corresponding to greater mass-loss and thus a higher luminosity. Combining this and the previous effects yields a total increase in luminosity less than 0.5 per cent (0.849–0.616 W m|$^{-2}$|⁠).

We present compilations of paleotemperature reconstructions and atmospheric compositions for the Earth in Section 4. We use inverse climate models to compute the solar constant S required to produce temperate surface temperatures. We find that our solar luminosity and its uncertainty is consistent with a temperate Earth, but it would not easily reproduce higher temperatures. There is controversy over temperature reconstructions |$\gt 60\,^\circ$|C, but separate from this controversy it is worth noting that temperature estimates are |$\gt 20\,^\circ$|C through much of Earth’s history. This indicates that the magnitude of the FYS problem is actually worse than simply reproducing a temperature comparable to the present-day surface temperature of 288 K.

We present the luminosity, solar constant, X-ray luminosity, and the |$T_{\rm eff}$| as a function of age for our solar models. The increase in these properties (other than |$T_{\rm eff}$|⁠) at early ages could be interesting with regard to the atmospheres and habitability of young exoplanets as the increase in high-energy photons could lead to more atmospheric evaporation or change the atmospheric chemistry. Future work could aim to include the effects of star-spots on the stellar spectrum and whether that alters the impact of young host stars on their exoplanets.

ACKNOWLEDGEMENTS

B. Scott Gaudi was supported by the Thomas Jefferson Chair for Space Exploration endowment from the Ohio State University. Shawn D. Domagal-Goldman and Sandra T. Bastelberger acknowledge support from the Goddard Space Flight Center (GSFC) Sellers Exoplanet Environments Collaboration (SEEC), which is supported by NASA’s Planetary, Astrophysics and Heliophysics Science Divisions’ Research Program; and from the Virtual Planetary Laboratory Team, a member of the NASA Nexus for Exoplanet System Science, funded via NASA Astrobiology Program Grant No. 80NSSC18K0829.

DATA AVAILABILITY

The outputs of our best yrec standard solar model and evolutionary history are available online. The properties in the online trackfiles are described in Table A1. Other data are available upon request to the authors.

REFERENCES

APPENDIX A: CLIMATE MODELS

We use Hayworth et al. (2020)’s version of the climate model originally developed by Kasting et al. In this version, the k-coefficients for H|$_2$|O and CO|$_2$| have been updated with HITRAN2016 line list (Gordon et al. 2017, see Vidaurri et al. 2022 for details). The climate code does not include explicit clouds. Their effects on climate is accounted for implicitly by increasing the surface albedo to a higher value of 0.24 (compared to the modern value of |$\sim$| 0.125), which we kept constant. This value was chosen as it reproduces modern Earth’s mean temperature of 288 K with modern atmospheric composition. Tropospheric water vapour was calculated according to the adapted version of the Manabe & Wetherald (1967) parametrization – which is widely used in 1D models – Geology from Kasting & Ackerman (1986). Their adapted version is based on Cess (1976) and allows the humidity in the upper troposphere to increase in response to higher temperatures.

For our inverse runs, a moist adiabat was used in the troposphere. In the stratosphere, the temperature profile was assumed to be isothermal at 171 K. This value was determined by running a set of 23 forward runs for solar constant values spanning 0.7–1.0 and CO|$_2$| concentrations between 0.1 and 10 mbar as well as 20 mbar. The inverse model was then run with the surface temperatures determined in the forward runs. The stratospheric temperatures were adjusted to reproduce the input solar constant to the fourth decimal with values ranging from 163 to 180 K, yielding an average of 171 K.

1D climate models without explicit cloud representation tend to overestimate forcing by greenhouse gases and insolation changes (e.g. Fauchez et al. 2018), and tend to overestimate greenhouse gas concentrations needed to maintain a temperate young Earth compared to more sophisticated GCMs (e.g. Charnay et al. 2020). However, we note that our forward climate model reproduces the GCM results of Charnay et al. (2017) well enough for the questions addressed in this work (see Fig. A1).