© ESO 2021

1. Introduction

Young stellar objects grow mainly by accreting mass from a circumstellar or protoplanetary disk. Planet formation theories predict that the composition of the disk should vary with time owing to the growth of dust grains into centimeter-sized grains, or “pebbles”, that drift toward the central star, the formation of planetesimals and planets, and the presence of disk winds. In our previous studies, we provided constraints on the evolutionary models of protostellar and pre-main-sequence (pre-MS) stars, including accretion (Kunitomo et al. 2017), and also investigated the impact of the evolving composition of the accreting material on the composition of the stellar surface (Kunitomo et al. 2018). We showed that some of the anomalies observed in the chemical composition of stellar clusters, λ Boo stars, and binary stars can be explained by planet formation processes. Moreover, our studies revealed that the large deficit in the Sun’s refractory elements compared to that in solar twins, as found by Meléndez et al. (2009), may be explained by the formation of giant planets in the Solar System with a very high rock-to-ice ratio (Kunitomo et al. 2018) (see Booth & Owen 2020, for an alternative explanation). Given the detailed information and constraints available regarding the composition and structure of the Sun, it would be interesting to investigate the imprint of planet formation processes in the Sun that may be detected by observations.

The main characteristics of the Sun, namely, its radius, luminosity, and age, are precisely known. However, the theoretical modeling of the Sun’s atmospheric composition, which is determined using high-resolution spectroscopy, has significantly evolved in the past decade. This is owing to the replacement of old one-dimensional atmospheric models by more accurate three-dimensional models with new atomic data that account for convection and nonlocal thermodynamic equilibrium effects, which substantially improved the fit to the observed spectral lines (Asplund et al. 2005, 2009). As a consequence, the inferred present-day solar surface metallicity decreases significantly from Z_surf = 0.018 (Grevesse & Sauval 1998, hereafter GS98) to Z_surf = 0.013 (Asplund et al. 2009, hereafter AGSS09).

The best constraints on the internal structure of the Sun are obtained from helioseismology, which helps determine not only the location of the base of the convective zone (CZ) but also the sound speed from the radius r ≈ 0.1 R_⊙ to the solar surface, where R_⊙ is the solar radius (e.g., Basu 2016). However, it has been shown that the old high-Z solar models result in significantly better fits to the observed sound speed compared to the low-Z solar models. This is the so-called solar abundance problem (see reviews by Asplund et al. 2009; Serenelli 2016), which has been the focus of many studies. Potential solutions to the solar abundance problem include the consideration of increased opacities (Bahcall et al. 2005; Christensen-Dalsgaard et al. 2009; Villante 2010; Ayukov et al. 2013; Bailey et al. 2015; Vinyoles et al. 2017), increased efficiencies of diffusion (so-called extra mixing; Christensen-Dalsgaard et al. 2018; Buldgen et al. 2019), diffusion due to the solar rotation (so-called rotational mixing; Yang 2019), helium-poor accretion (Zhang et al. 2019), an updated solar composition (Young 2018), and revised nuclear reaction rates (Ayukov & Baturin 2017). Accretion models with a time-dependent composition, in particular, low-Z accretion in the late pre-MS phase, have also been investigated but were found to be unsuccessful in solving the solar abundance problem (Castro et al. 2007; Guzik & Mussack 2010; Serenelli et al. 2011).

To investigate the effect of planet formation on the structure and composition of the Sun requires a quantitative analysis of: (1) various planet formation processes and (2) solar evolution models that consider the accretion of matter with a time-dependent composition to reproduce all observational constraints. In the present work, which is an extension of our previous works (Kunitomo et al. 2017, 2018), we conducted a thorough search for solutions to the solar abundance problem by applying our accretion models with a time-dependent composition to optimized solar evolution models.

The remainder of this paper is organized as follows. In Sect. 2, we discuss how the growth of dust grains, the radial drift of pebbles, and planet formation affect the composition of the gas accreted by the proto-Sun during its protostellar and pre-MS evolution. In Sect. 3, we describe the method used to compute our stellar evolution models and the procedure followed to optimize these models to reproduce the observational constraints. In Sect. 4, we present the results of different solar evolution models and compare them with observations; in addition, we show that an increased opacity appears to be the most likely solution for the solar abundance problem. Using these models, in Sect. 5, we examine the effect of planet formation on the composition and structure of the present-day Sun. Our results are summarized in Sect. 6.

2. Context

2.1. Evolution of the composition of protoplanetary disks

During the collapse of a molecular cloud core, a large fraction of the material first forms a circumstellar disk around the central star (or stars) (Inutsuka 2012). Next, the transport of angular momentum from the inner region of the disk to the outer region causes the matter to be accreted by the central star(s) (Lynden-Bell & Pringle 1974) on a timescale of order 1–10 Myr (e.g., Hartmann et al. 2016). During this time, initially micron-sized dust grains grow larger in size, as indicated by the dust emission at millimeter wavelengths (e.g., Beckwith et al. 1990). A comparison of the mass of heavy elements present in exoplanetary systems and that inferred to be present in protoplanetary disks indicates that these grains must grow further in size to form planetesimals and even planets during the early stages of disk evolution (Manara et al. 2018; Tychoniec et al. 2020). Dust growth in the Class 0/I phases of young stellar objects has also been investigated in several theoretical studies (e.g., Tsukamoto et al. 2017; Tanaka & Tsukamoto 2019).

The formation of planetesimals and planets is likely to lead to a decrease in the metallicity of the gas accreted by the central star(s). Kunitomo et al. (2018) estimated that approximately 97 to 168 M_⊕ of heavy elements in the Solar System were either used to form planets or were ejected from the system. This is a small value compared to the total mass of heavy elements present in the Sun, which is estimated to be approximately 5000 M_⊕; however, it can still cause the metallicity of the accreted gas to change, particularly during the later stages of stellar accretion. Moreover, it has been proposed that the peculiar composition of λ Boo stars can be explained by the presence of giant planets in the surrounding protoplanetary disks (Kama et al. 2015).

It is also worth noting that as the grains grow in size, they drift rapidly inward in the protoplanetary disk (Adachi et al. 1976; Weidenschilling 1977). This can create a wave of solids that may be lost to the central star in the absence of any other processes (Stepinski & Valageas 1996; Garaud 2007). Circumstellar disks have been proposed to be large (100 au or more in size) and expanding (e.g., Hueso & Guillot 2005), which implies that the outer part of the disk acts as a reservoir in which the dust grains grow and from which they are gradually released. By assuming the standard theory of grain growth and adopting the α-disk model (Shakura & Sunyaev 1973), Garaud (2007) demonstrated that the timescale in which the solids are released from the outer disk is determined by a balance between the outward propagation of the disk and the growth of the dust grains.

The metallicity Z_acc of the gas accreted by the proto-Sun can be calculated from the ratio between the mass fluxes of the dust Ṁ_d and gas Ṁ_acc, such that

(1)

The above equation assumes that the metallicity is equal to the amount of condensing material. We note that because most elements heavier than hydrogen and helium condense in the protoplanetary disk, the above assumption is justified in our case.

Figure 1 shows the evolution of Z_acc with time and the mass of the proto-Sun, as obtained from various studies on the evolution of dust in protoplanetary disks. In the beginning, Z_acc = Z_proto, where Z_proto is the primordial metallicity of the molecular cloud core. Eventually, Z_acc increases with time due to the incoming pebble wave before decreasing sharply owing to the removal of all or most of the solids from the disk.

Fig. 1. Evolution of the metallicity Zacc of the accreted gas in units of the initial metallicity Zproto as a function of time (left panel) and mass of the proto-Sun (right panel). Several evolutionary models of circumstellar disks that include both gas and dust are shown (see text for details). The protostar mass was calculated assuming some mass loss (e.g., due to partial photoevaporation of the circumstellar disk) except in two cases, which are represented by dashed lines.

The first set of disk models shown in Fig. 1 corresponds to the analytical solutions obtained by Garaud (2007) for α = 10⁻² and 10⁻³ and R₀ = 30 au and 1000 au, where α is the dimensionless turbulent viscosity parameter and R₀ is the initial disk radius. The second model corresponds to the numerical solution obtained by Appelgren et al. (2020) for a standard one-dimensional disk model. In this case, the authors assumed a variable α such that its value increased in the presence of gravitational instabilities but had a fixed minimum value of α_min = 10⁻². Finally, the third set of models corresponds to the solutions obtained by Elbakyan et al. (2020) for a two-dimensional disk model. The accretion is highly variable in this case, and hence, we plot only the mean value over a timescale of 2.5 × 10⁴ years.

The left panel of Fig. 1 shows that the pebble wave appears at less than 1 Myr and lasts until approximately 10 Myr. This timescale is extremely difficult to estimate because it depends on several factors, including the turbulent viscosity in the disk, initial angular momentum of the disk material, overall structure of the disk, and the model used to calculate the grain growth. The peak value of Z_acc/Z_proto varies from 1.5, for models L and M of Elbakyan et al. (2020), to approximately 5, for the model proposed by Appelgren et al. (2020). It should be noted that Mousis et al. (2019) found an even higher metallicity peak in the range 10–20. The peak in metallicity is followed by a sudden drop to a value close to the initial metallicity in the models proposed by Elbakyan et al. (2020) and to zero in the remaining models; however, in reality, the final metallicity value should lie somewhere in between. This is because we know from observations that protoplanetary disks are never completely cleared of dust (e.g., Dullemond & Dominik 2005), indicating that the perfect depletion obtained by Garaud (2007) and Appelgren et al. (2020) is due to the simplifications adopted in their models to account for the grain size distribution. Conversely, these models do not consider planetesimal and planet formation processes, which can remove solids from the system and thus lead to an additional decrease in the metallicity of the accreted gas. In addition, we note that fully formed planets can efficiently filter out pebbles and prevent them from reaching the protostar (Guillot et al. 2014).

For the purpose of this study, it is also important to estimate how the metallicity of the accreted gas varies as a function of the mass of the protostar. The disk models above neglect the change in the mass of the protostar owing to disk accretion. Therefore, we estimated the mass of the protostar as a function of time, as follows:

(2)

where t is time, M_g(t) is the gas mass in the disk, M_final is the final mass of the star (note that M_final = 1 M_⊙ for the Sun), and M_XY, lost is the mass ejected from the disk (e.g., via photoevaporation or magnetohydrodynamic (MHD) winds; see Hollenbach et al. 2000; Suzuki & Inutsuka 2009; Alexander et al. 2014). We stop the calculation of M_⋆ at a time t_final when M_g(t = t_final)=M_XY, lost. Following Guillot & Hueso (2006), we assumed that predominantly hydrogen and helium were lost from the disk (also, see Gorti et al. 2015; Miyake et al. 2016). We note that photoevaporation is more effective than MHD disk winds in the selective removal of hydrogen and helium during the later stages of accretion (Kunitomo et al. 2020). Although the total mass of the disk lost via winds is not well determined, we expect M_XY, lost ≲ 0.05 M_⊙.

The right panel of Fig. 1 shows that the peak in the metallicity of the accreted gas occurs after approximately 90% of the mass has been accreted onto the protostar. For most of the models, the peak occurs at a high value of the protostellar mass, namely, between 95% and 98%. For the disk model proposed by Appelgren et al. (2020) with no mass loss, a pebble wave occurs when the mass of the protostar is approximately 85% of the final mass of the star. This is because in their model the disk extends beyond 1000 au. However, such a large disk is likely to get photoevaporated rapidly due to external radiation (e.g., Guillot & Hueso 2006; Anderson et al. 2013). Thus, for our simulations of the Sun, we considered a situation in which the metallicity starts to increase after 0.85 M_⊙ of the gas has been accreted, peaks between 0.9 M_⊙ and 0.98 M_⊙, and finally decreases to a very low value.

2.2. Internal structure of the pre-main-sequence Sun

One may expect the metallicity of the accreted gas to have a significant impact on the composition of the Sun’s atmosphere (see, e.g., Chambers 2010). However, this is not the case because of the internal structure of the proto-Sun. Studies on the evolution of the Sun, which is in agreement with the evolution of stars in young clusters of similar mass as the Sun (Kunitomo et al. 2017), indicate that the Sun must have been almost fully convective for the first 1–2 Myr. Then, the CZ starts to recede slowly as the Sun evolves toward the main sequence (MS; which takes approximately 40 Myr), finally reaching a stage where it is only 2.5% of the total mass of the present Sun.

Figure 2 shows the first 30 Myr of the evolution of the solar CZ for three models adopted from Kunitomo et al. (2018): a standard accretion model with ξ = 0.5, which corresponds to the maximum accretion efficiency; a limiting model with ξ = 0.1, which is in agreement with the observational data from young clusters but corresponds to an accretion efficiency of only 10%; and a cold accretion model with ξ = 0 for reference (see Kunitomo et al. 2017, for a complete discussion). The ξ value corresponds to the ratio of the accretion heat injected into the protostar to the liberated gravitational energy of accreting materials (Kunitomo et al. 2017). In the present work, we adopted models equivalent to the standard model with ξ = 0.5 (see Sect. 3.1.2 for details).

Fig. 2. Early evolution of the solar interior. The total mass of the accreting proto-Sun is indicated by the dashed line. The surface CZ is depicted as a cloudy region (see Kippenhahn & Weigert 1990) delimited by a radiative zone that grows outward from the center of the Sun. Three evolution models are shown corresponding to different values of the accretion efficiency parameter ξ: a classical evolution model with ξ = 0.5 (red line), a cold accretion model with ξ = 0.1 (green line) that is a lower limit of the ξ value in the observations of young clusters (Kunitomo et al. 2017), and a model with ξ = 0 (gray line) corresponding to a theoretical proto-Sun formed in the absence of any accretion heat. The right-hand y-axis provides the radius corresponding to the mass on the left-hand y-axis of the present-day Sun (which is 4.567 Gyr old). Bottom panel: highlights the key physical processes that occur during the growth of the proto-Sun: the presence of a circumstellar gas disk (during the first million years); an increase in the metallicity Zacc of the accreted gas due to a pebble wave; the concomitant formation of planetesimals and planets; and a sudden decrease in Zacc.

The lifetime of protoplanetary disks is still uncertain but is thought to be less than 10 Myr (see, e.g., Haisch et al. 2001; Kennedy & Kenyon 2009; Fedele et al. 2010; Hartmann et al. 2016). In fact, the magnetism observed in meteorites indicates that the lifetime of the protoplanetary disk in the Solar System was ≃4 Myr (see, e.g., Wang et al. 2017; Weiss et al. 2021). Thus, during the protoplanetary disk phase, the interior of the Sun was largely convective, with the CZ encompassing at least 50% of the Sun’s mass, whereas the typical disk mass is ∼0.01 M_⊙ (see, e.g., Williams & Cieza 2011). This implies that any anomalies in the solar composition due to the accretion of the protoplanetary disk gas with the pebble wave or with the depletion of heavy elements by the formation of planetesimals and planets were likely to be suppressed by the large CZ.

Thus, we can expect two consequences from the above discussion: (1) the signature left by both the high-Z and low-Z accretion is suppressed by roughly one to two orders of magnitude and (2) convective mixing leads to a uniform composition in a large fraction of the solar interior. As seen in Fig. 2, a mass coordinate of 0.5 M_⊙ corresponds to approximately 27% of the present-day solar radius, which indicates that the signatures of the planet formation processes are buried deep in the solar interior (mostly in its nuclear burning core).

3. Computation method

In this section, we describe how we (i) simulate stellar evolution, (ii) compare our results with observations, and (iii) minimize the total χ² value by changing the input parameters.

3.1. Stellar evolution with accretion

We used the one-dimensional stellar evolution code MESA version 12115 (Paxton et al. 2011, 2013, 2015, 2018, 2019). For details of the computational method used in this work, we refer the readers to Kunitomo et al. (2017, 2018), and the series of papers by Paxton et al. Below we briefly summarize the method and the various parameters used.

3.1.1. Initial conditions

Stars are formed via the gravitational collapse of a molecular cloud core. A protostar (or second hydrostatic core) forms after the formation of a transient hydrostatic object (the so-called first core; see Larson 1969; Inutsuka 2012). Radiation hydrodynamic simulations have suggested that the initial mass of a protostar is typically ∼0.003 M_⊙ (Masunaga & Inutsuka 2000; Vaytet & Haugbølle 2017). In this study, we used a stellar seed of mass 0.1 M_⊙, radius 4 R_⊙, and metallicity Z = 0.02, to avoid numerical convergence issues at very low protostellar masses in the new version of the MESA code. Recent studies have shown that the thermal evolution of the protostar depends on the entropy s_acc of the accreted material (e.g., Hartmann et al. 1997; Baraffe et al. 2009; Hosokawa et al. 2011; Tognelli et al. 2015; Kunitomo et al. 2017; Kuffmeier et al. 2018). The initial seed corresponds to a protostar of mass 0.1 M_⊙ that has evolved via hot accretion from its birth mass of ∼0.003 M_⊙¹.

For the non-accreting cases (see Sect. 4.1), we used an initial stellar seed of mass 1 M_⊙ and a central temperature of 3 × 10⁵ K (i.e., corresponding to the top of the Hayashi track). The seed has a uniform composition. The mass fractions of hydrogen, helium, and metals of the seed are X_proto, Y_proto, and Z_proto, respectively. We note that the Kelvin-Helmholtz timescale at the top of the Hayashi track is short (see, e.g., Stahler & Palla 2004), and hence, the pre-MS evolution without accretion is not sensitive to the choice of the initial central temperature.

3.1.2. Mass accretion

We adopted the following mass accretion rate Ṁ_acc in the protostellar and pre-MS phases²:

(3)

(see Hartmann et al. 1998, for details about the exponent −1.5). In our fiducial model, we set t₁ = 31, 160 yr so that M_⋆ could reach 1 M_⊙ at the end of the accretion phase, that is, at t_acc  =  10 Myr. The evolution of M_⋆ and Ṁ_acc in our fiducial model is shown in Fig. 3. In some of our simulations, we changed the accretion timescale t_acc and t₁, but we always used the same exponent (i.e., −1.5) for t > t₁. As discussed in Sect. 2.2, accretion generally ends several million years; therefore, t_acc  =  10 Myr corresponds to the upper limit of the observational constraint.

Fig. 3. Evolution of stellar mass M⋆ (red solid line) and accretion rate Ṁacc (green dashed line) with time in the fiducial case. The gray vertical lines indicate t1 and tacc in the fiducial case. The gray dotted lines indicate M⋆ and Ṁacc for tacc = 3 and 20 Myr.

We note that previous studies, such as Serenelli et al. (2011) and Zhang et al. (2019), considered accretion only in the pre-MS phase, whereas, in this study, we considered accretion in the protostellar phase, which led to a larger mass of the accreted material. Moreover, Serenelli et al. (2011, see their Sect. 3.2.3) and Zhang et al. (2019, see their Sect. 2.4) introduced accretion in their models after a certain time³ to allow the pre-MS Sun to develop a radiative core. These studies needed a non-accreting phase because they adopted arbitrary initial conditions, whereas our model is based on the current understanding of star formation. We also note that the non-accreting timescale in some of the models of Serenelli et al. (2011) exceeded the typical disk lifetime (i.e., several million years).

In the present work, we did not consider mass loss. Vigorous stellar winds have been known to flow out from young stars (Wood et al. 2005; Suzuki et al. 2013); moreover, Zhang et al. (2019) suggested that mass loss has the potential to alter the surface composition of the star. Although in this study we focused only on the effects of planet formation and opacity enhancement on the stellar composition, the effect of mass loss should be considered in future studies.

We used the same accretion model characterized by the parameter s_acc as in Kunitomo et al. (2017, 2018), who assumed that the accretion heat injected into the protostar L_add is a fraction of the gravitational energy liberated by accretion. Consequently, we assumed that L_add = ξGM_⋆Ṁ_acc/R_⋆, where ξ is a dimensionless parameter, G is the gravitational constant, and R_⋆ is the stellar radius. The accretion heat is assumed to be uniformly distributed throughout the star (see related discussion in Kunitomo et al. 2017). Although we set ξ = 0.1 (see Kunitomo et al. 2017) in our accretion models, we chose a high-entropy initial seed (see Sect. 3.1.1). This implies that our thermal evolution calculations are equivalent to those of standard accretion models (Stahler et al. 1980; ξ = 0.5 models in Kunitomo et al. 2017). Thus, the evolution of the CZ in our models corresponds to the ξ = 0.5 line shown in Fig. 2 (also, see Appendix A). In this context, we can expect the calculations with ξ = 0.1 and a low-entropy initial seed to capture a slightly larger effect due to planet formation processes. We leave this for future work.

3.1.3. Convection

We adopted the mixing-length theory for convection proposed by Cox & Giuli (1968). We assumed the ratio of the mixing length to the local pressure scale height (H_p) to be α_MLT. We also considered the effect of the composition gradient on the convective stability (i.e., the Ledoux criterion).

Previous studies have suggested that additional mixing (e.g., convective overshooting and tachocline circulation), especially at the base of the CZ, plays an important role in solar sound-speed anomalies (see, e.g., Christensen-Dalsgaard et al. 2018; Buldgen et al. 2019; Zhang et al. 2019). However, these processes are poorly understood and are not well constrained. Christensen-Dalsgaard et al. (2018) showed that extending the Sun’s CZ or adding extra mixing in the tachocline region leads to a better agreement between theoretical models and helioseismic constraints. Zhang et al. (2019) included convective overshoot in their models by adopting an exponentially decreasing diffusion coefficient and considering a change in the luminosity due to kinetic energy transfer in the overshooting region. They showed that the radial location R_CZ of the base of the CZ is sensitive to the underlying overshooting model.

In this study, we considered the conventional convective overshooting model proposed by Herwig (2000), in which the diffusion coefficient decreases exponentially from the convective–radiative boundary. We assumed the e-folding length of the diffusion coefficient to be f_overshootH_p. In addition, we adopted the same overshooting parameters below and above the CZ, irrespective of the presence of nuclear burning. Although we also considered semiconvection in our models, we confirmed that it has little effect on the results. In this study, α_MLT and f_overshoot are the two free parameters. Based on the results of Zhang et al. (2019) and given the fact that a detailed exploration of convective overshooting models is beyond the scope of the present work, we slightly relaxed the constraint on R_CZ when seeking the best solutions (see Sect. 3.2).

3.1.4. Abundance tables

We used the abundance table presented in AGSS09 unless otherwise mentioned. We note that recent studies have suggested some modifications to the tables in AGSS09⁴. However, these modifications are relatively minor considering the goals of the present study and are therefore not considered here (see Buldgen et al. 2019, for the effect of modified solar abundances).

We also performed simulations using the abundance table presented in GS98 for comparison with previous studies. We adopted the table in which the abundances of refractory elements were modified using meteorites.

3.1.5. Opacity

In stellar evolution calculations, the Rosseland-mean opacity κ is determined by interpolating the opacities listed in standard opacity tables using local thermodynamic and compositional quantities, such that κ = κ(ρ, T, X, Z), where ρ is the density, T is the temperature, X is the hydrogen mass fraction, and Z is the metallicity. We adopted the OPAL opacity table (Iglesias & Rogers 1996) for the cases using the AGSS09 composition and the OP opacity table (Seaton 2005) for the cases using the GS98 composition. The opacity table presented in Ferguson et al. (2005) was used to determine the opacities in the low-temperature regions for both cases. We refer the readers to Paxton et al. (2011) for more details. We confirmed that the simulations with the OP table yielded very similar results to those with the OPAL table for the AGSS09 composition (Buldgen et al. 2019, also, see Appendix A).

As mentioned in Sect. 1, we aim to investigate the effect of opacity enhancement. We modeled the opacity enhancement using a dimensionless factor δ_κ, such that

(4)

Below, we describe our δ_κ model. Laboratory experiments have not been able to reproduce the real conditions at the base of the solar CZ, which has led to uncertainties in the opacity. Bailey et al. (2015) performed experiments with conditions close to those in the real Sun and showed that the wavelength-dependent opacity of iron was 30%–400% higher than previously thought, and also, the Rosseland-mean opacity was increased by 7 ± 3% (see also Nagayama et al. 2019). Le Pennec et al. (2015, see their Fig. 2) obtained the contribution of each element to the opacity. Iron was shown to have three peaks at log(T/K) = 5.66, 6.45, and 7.18⁵. We fitted the contribution of iron to the opacity by the sum of three Gaussian functions of T and we used the same function to model δ_κ, considering the possibility that further uncertainties may remain in the iron opacity and assuming that the iron abundance is inhomogeneous in the solar interior. Thus, we modeled the opacity increase as

(5)

where A_i is a free parameter, and b_i and c_i were derived by fitting Fig. 2 of Le Pennec et al. (2015), leading to (b₁, b₂, b₃)=(5.66, 6.45, 7.18) and (c₁, c₂, c₃)=(0.22, 0.18, 0.25).

We note that Villante (2010) suggested another δ_κ model that increases linearly with T. We note that our opacity increase model is based on actual data of the contribution to the opacity, whereas the model in Villante (2010) is more ad hoc in nature. Although in this study we only present the results using Eq. (5), we confirmed that both opacity-increase formalisms improve the solar evolution model in a similar manner.

3.1.6. Composition of accreted material

For all the models with protostellar accretion, we fixed the mass fraction of deuterium to X_D = 28 ppm (see Kunitomo et al. 2017, and references therein) and that of ³He also to 28 ppm (Mahaffy et al. 1998). In the non-accreting cases, we set X_D = 0, assuming that the deuterium was completely depleted in the protostellar phase.

For the models with accretion, we considered three models to account for the composition of the accreted material (i.e., ¹H, ⁴He, and metals): homogeneous accretion, metal-poor accretion, and helium-poor accretion (see Table 1). In the first model, the helium mass fraction Y_acc and metallicity of the accreted materials Z_acc are constant with time; therefore, Y_acc = Y_proto and Z_acc = Z_proto, where the subscript “proto” indicates the initial value.

In the metal-poor accretion model, Z_acc is time-dependent because of the underlying planet formation processes. We adopted a simple model for Z_acc to capture the following two processes: pebble drift and planetesimal formation (see Sect. 2.1). Figure 4 shows the evolution model of Z_acc used in this work. In the early phase (Phase I; M_⋆ ≤ M₁), Z_acc = Z_proto is constant. Once the Stokes number of dust grains increases and the grains start to migrate inward, Z_acc increases monotonically with time. We refer to this phase as Phase II or the pebble-accretion phase. Planetesimals begin to form when M_⋆ = M₂. Subsequently, in Phase III (i.e., the metal-poor accretion phase), drifting dust grains are captured by the planetesimals and are not allowed to reach the proto-Sun, causing Z_acc to abruptly decrease to 0. When M_⋆ reaches 1 M_⊙, the accretion stops abruptly (i.e., Ṁ_acc = 0). After the accretion phase, the proto-Sun evolves to become the Sun (Phase IV).

Fig. 4. Sketch of the evolution of the metallicity Zacc of the accreted material with stellar mass M⋆. In this work, we used the model depicted by the blue solid line and assumed that the difference between this model and the one depicted by the green dotted line has little impact on the internal structure of the Sun (see text for details).

When dZ ≡ Z_acc − Z_proto ≠ 0 (i.e., M_⋆ > M₁), we assumed that X_acc = X_proto − 0.7dZ and Y_acc = Y_proto − 0.3dZ, where X_acc is the hydrogen mass fraction of accreted materials and X_proto is the initial value of X_acc. We note that because X_D and the ³He abundance were kept constant, only the ¹H and ⁴He abundances changed with time.

We point out that our assumption of Z_acc = 0 in Phase III is a major simplification. In reality, grains are not completely filtered by planets and planetesimals. Some of the grains collide with the planetesimals leading to a small but nonzero Z_acc, which corresponds to the green dotted line in Fig. 4 (see also Sect. 2.1). However, because the mass of the solar CZ is large during the accretion phase (see Fig. 2), the difference between these two models has little impact on the results.

In the He-poor accretion model, the helium abundance Y_acc of the accreted material was assumed to vary with time, such that

(6)

The He-poor accretion was proposed by Zhang et al. (2019) as a solution to the solar sound-speed anomaly. Zhang et al. (2019) considered a constant Y_acc but had a larger initial stellar mass in their simulations than in the present study (see Sect. 3.1.2); thus, their simulations are equivalent to ours in the late phase (i.e., M_⋆ > M₁). In the He-poor accretion models, we assumed that Z_acc = Z_proto and X_acc = X_proto − dY, where dY ≡ Y_acc − Y_proto.

3.1.7. Other input physics

We adopted the OPAL equation-of-state tables updated in 2005 (Rogers & Nayfonov 2002). We used the nuclear reaction rates in Caughlan & Fowler (1988) and Angulo et al. (1999) with certain modifications (see Paxton et al. 2011, for details). The prescription for element diffusion presented in Thoul et al. (1994) was also used. We note that the JINA reaction table is also available in MESA; however, we confirmed that the results are insensitive to the choice of the reaction table (see Appendix A).

The outer boundary condition used in this study was given by the atmospheric tables (see Sect. 2.7 of Kunitomo et al. 2018). Although previous studies adopted different boundary conditions⁶, we confirmed that the outer boundary conditions did not have a significant impact on the results.

3.1.8. Summary of the input parameters

Table 2 summarizes the input parameters used in this work. In the runs with opacity enhancement and Z_acc evolution, the number of input parameters can be up to 10.

3.2. Comparison with observations

When the elapsed time of the stellar evolution simulations described in Sect. 3.1 reached the solar age, we compared the simulation results with spectroscopic and helioseismic observations of: the ratio of the surface metallicity to the surface hydrogen abundance (Z/X)_surf; the surface helium abundance Y_surf; the location of the convective–radiative boundary R_CZ; the root mean square (rms) of δc_s (see Eq. (7)); the bolometric luminosity L_⋆; and the effective temperature T_eff. We define the observed minus calculated sound speed as

(7)

To obtain the rms(δc_s), we compared our simulated results with the observed data provided in Basu et al. (2009, see their Table 3). We interpolated the c_s profile obtained at the solar age (typically ∼3000 grids) to the locations of 37 points given in Basu et al. (2009) and thus calculated the rms value of δc_s (see discussions in Sect. 4.6).

Table 3 summarizes the observed parameters and their 1σ uncertainties with three exceptions. First, although the uncertainty of R_CZ in the helioseismic observations is 0.001 R_⊙ in Bahcall et al. (2005), we relaxed this constraint by a factor of 10. This is because the R_CZ value is likely to be affected by the uncertainties in the convective overshooting model, which are beyond the scope of this study (see Sect. 3.1 of Zhang et al. 2019, and Sect. 3.1.3 in this work). Second, we chose the rms(δc_s) uncertainty (10⁻³) based on the results obtained with the GS98 composition presented in Serenelli et al. (2009) (see their Table 2), rather than from the helioseismic observations. Finally, the uncertainties in L_⋆ and T_eff are arbitrary; however, we confirmed that our results are not sensitive to these values because all the models reproduce these quantities well (< 0.001 dex and 3 K, respectively).

We assumed the solar age to be 4.567 Gyr (Amelin et al. 2002) after Ca-Al-rich inclusions (CAIs) were condensed. We note that there could be a time difference between the formation of CAIs and time taken for M_⋆ to reach 0.1 M_⊙. Assuming that CAIs were probably formed in the protosolar nebula, this time difference is at most several million years, and therefore, we neglect it in the present study.

3.3. Chi-squared tests

By comparing our simulation results with observations, we derived the χ² value, which is given by

(8)

where q_i is the simulated value of a particular quantity, q_i,target is the corresponding observed (or target) value, N = 6 is the total number of target quantities, and σ is the uncertainty of each quantity listed in Table 3. We aimed to search for a set of input parameters that minimized the χ² value. To do so, we used the downhill simplex method (Nelder & Mead 1965). Typically, ∼200 simulations were needed for the minimization. We performed such simulations under a variety of settings.

In the simplex algorithm, we imposed constraints on the parameter space composed of Z_acc,max, M₁, and M₂. We did not consider any cases where the total metal mass of the accreted material was outside the range [0.5, 2] × Z_proto(1 M_⊙ − 0.1 M_⊙).

It should be noted that in general standard solar models are optimized in the non-accreting case with three input parameters, namely α_MLT, Y_proto, and Z_proto, and three target quantities, namely L_⋆, either T_eff or R_⋆, and (Z/X)_surf (see, e.g., Vinyoles et al. 2017). Farag et al. (2020, see their Tables 1 and 2) showed that, if the above three input and target quantities are used, then the optimized solutions using the MESA code are similar to those in Vinyoles et al. (2017). We confirmed that, using these three target quantities, we could reproduce exactly the results of Farag et al. (2020). However, we found that in some cases this can result in values of Y_surf outside the range of observed values (see, e.g., Table 4 of Vinyoles et al. 2017). To search for the solutions that explain all the observable quantities, we used as constraints the six quantities in Table 3.

We also performed simulations in which the 37 points used to calculate the δc_s values (see Sect. 3.2) were considered as independent and were added to the five other target values in Table 3 (thus N = 42; see, e.g., Villante et al. 2014). We found that the models that poorly fit the available constraints were not improved and that our best models were not changed by this new approach.

We note that the solutions with the simplex method sometimes fall into local minima. Although we carefully chose the initial values of the input parameters to avoid this problem, future simulations using the Markov chain Monte Carlo method are encouraged.

4. Solar models fitting the observational constraints

In this section, we present the results of the simulations optimized using the simplex method for different conditions. First, in Sects. 4.1 and 4.2, we show that our results are in agreement with those of previous studies. Next, we present our results for the metal-poor accretion (Sect. 4.3) and opacity increase (Sect. 4.4) models. The parameter settings are summarized in Table 1. The detailed results are provided in Tables A.1 and A.2.

4.1. Non-accreting models

To date, most studies have performed solar evolution simulations without including accretion. In this work, we simulated such non-accreting models to ensure consistency with the previous studies.

We performed four sets of simulations (Runs 1–4 in Table 1) adopting either the GS98 or AGSS09 abundance tables, and with and without convective overshooting. We started the simulations with a 1 M_⊙ star having a fully convective structure (see Sect. 3.1.1). The initial stellar composition, α_MLT, and f_overshoot were optimized using the simplex method.

We found that the converged c_s profiles are in good agreement with those of previous studies (see Fig. A.1), which validates our simulations with the MESA code. The optimized input parameters and corresponding best results are shown in Fig. 5. The maximum values of δc_s for the cases with the GS98 and AGSS09 composition were approximately 0.004 and 0.009, respectively, irrespective of the inclusion of overshooting (see Fig. A.1). The rms(δc_s) value for the AGSS09 composition was found to be in the range (3.2–3.4)×10⁻³, which is approximately twice as large as that for the GS98 composition, that is, (1.6–1.8)×10⁻³ (see Table A.2). The other results obtained with the AGSS09 composition (i.e., surface composition and R_CZ) were also found to be worse than those obtained with the GS98 composition, as indicated by the four times larger total χ² value in the AGSS09 case.

Fig. 5. Optimized input parameters (left) and results at the solar age (right). The left panels show αMLT, fovershoot, Yproto, Zproto, and stellar mass (M1 and M2) for nine models from top to bottom. The right panels show the χ2 value (see Eq. (8)), rms(δcs), (Z/X)surf, Ysurf and RCZ from top to bottom. The five models shown in each panel are: non-accreting models with GS98 composition (Runs 2 and 4); non-accreting models with AGSS09 composition (Runs 1 and 3); accreting models with time-dependent Zacc (Runs 7 and 8); models with He-poor accretion (Runs 5 and 6); and K2 model (Run 9) with A2 = 0.12, from left to right (see Table 1 for details). Squares and crosses represent the cases with and without overshooting, respectively. Diamonds and asterisks represent the cases with He-poor accretion for tacc = 12 and 20 Myr, respectively. The shaded regions in the right panels indicate χ2 < 1 or the 1σ uncertainties (see Table 3).

4.2. Helium-poor accretion

Next, we performed simulations with the He-poor accretion following Zhang et al. (2019). In this case, we started simulations with a 0.1 M_⊙ seed (see Sect. 3.1.1) and simulated its evolution from the protostellar to the solar age. The input parameters optimized using the simplex method were α_MLT, f_overshoot, initial composition of the accreted material, M₁, and Y_acc,min (see Eq. (6)). Considering that Zhang et al. (2019) required a longer accretion timescale of t_acc ≥ 12 Myr⁷, we set t_acc = 12 and 20 Myr in this work. We note that our fiducial model corresponds to t_acc = 10 Myr (see Table 1).

Figure 5 shows that He-poor accretion improves the rms(δc_s) value, as claimed by Zhang et al. (2019). The rms(δc_s) value was found to be 2.8 × 10⁻³ and 1.4 × 10⁻³ for t_acc = 12 and 20 Myr, respectively. Zhang et al. (2019) obtained even better (but qualitatively similar) rms(δc_s) values by considering mass loss by stellar winds and a more complex overshooting model.

Despite the above results, we regard the possibility of He-poor accretion to be unlikely. First, the giant planets in our Solar System, which also captured gas in the protoplanetary disk, do not show a large depletion in their helium abundance (see, e.g., Guillot & Gautier 2015). Jupiter and Saturn’s atmospheres are characterized by slightly lower-than-protosolar helium abundances, which is consistent with the helium settling in their interiors from an initial protosolar value (Mankovich & Fortney 2020). In contrast, Uranus and Neptune appear to have a helium abundance compatible with the protosolar value (Guillot & Gautier 2015).

Second, Zhang et al. (2019) surmised that the high first ionization potential of helium might lead to the accumulation of helium at the inner edge of the disk followed by He-poor accretion onto the proto-Sun. However, this would lead to ≃0.015 M_⊙ of the helium remaining in the disk (see their Table 4). It is difficult to understand why such a large amount of helium would not have eventually accreted onto the proto-Sun.

4.3. Metal-poor accretion

The possibility that metal-poor accretion may affect the structure of the Sun was first suggested by Guzik et al. (2005) and then tested by Castro et al. (2007), Guzik & Mussack (2010), Serenelli et al. (2011), and Hoppe et al. (2020). Although these studies observed some improvement due to the inclusion of metal-poor accretion, they failed to find a solution as good as the models with high-Z abundances (i.e., Grevesse & Noels 1993, GS98) with respect to δc_s, R_CZ, and Y_surf.

In this work, we revisited metal-poor accretion in a larger parameter space and in the framework of recent planet formation theories (see Sect. 2.1). While the aforementioned studies considered metal-poor accretion onto a pre-MS or MS Sun, we considered metal-poor accretion during the protostellar phase. We assumed the composition of the accreted material to vary with time, as shown in Fig. 4. We ran simulations with model MZvar that had seven input parameters: α_MLT, f_overshoot, Y_proto, Z_proto, Z_acc,max, M₁, and M₂. We also performed a series of calculations without including overshooting (model MZvar-noov).

Figure 5 shows the results obtained with models MZvar and MZvar-noov. We found that the results are almost identical to those of the non-accreting models, with no significant improvement in the χ² value. This indicates that metal-poor accretion cannot be a solution to the solar abundance problem. This is consistent with the findings of previous studies (e.g., Serenelli et al. 2011) and can be easily explained. As shown in Fig. 2, the proto-Sun has a nearly fully convective interior, implying that the effect of metal-poor accretion is heavily suppressed therein. This, in turn, limits the signatures of metal-poor accretion on the internal structure of the present-day Sun.

Although metal-poor accretion does not have an impact on the χ² test, it affects the metallicity profile of the solar interior. We will revisit this issue in Sect. 5.

4.4. Opacity increase

Before the experiments conducted by Bailey et al. (2015) showed that iron opacities in stellar interiors were underestimated, solar models with a ∼10–30% ad hoc increase in the opacities had already been examined and have even been shown to partly solve the solar abundance problem (e.g., Christensen-Dalsgaard et al. 2009; Serenelli et al. 2009; Christensen-Dalsgaard & Houdek 2010; Villante 2010). Guided by the results of Bailey et al. (2015), we adopted a physically motivated opacity-increase factor δ_κ that depends on three adjustable parameters A₁, A₂, and A₃, as defined by Eqs. (4) and (5). Plausible values of these parameters were obtained by minimizing the χ² value using the simplex method.

First, we obtained results that were insensitive to the value of the low-temperature parameter A₁. This is because this parameter modifies the opacities in the CZ, where temperature changes are insensitive to the opacity increase. Therefore, we set A₁ = 0 in our models.

Next, we constructed models with different values of A₂ in the range 0 to 0.22 and five variable input parameters for the simplex method: α_MLT, f_overshoot, Y_proto, Z_proto, and A₃ (corresponding to model K23 in Table 1). Another series of models was constructed by setting A₃ = 0 (model K2).

Figure 6 shows the results as a function of A₂. When A₂ = A₃ = 0, the results obtained with the AGSS09 abundances were found to be clearly much worse than those obtained with the GS98 abundances. However, when A₂ was increased (in both models K23 and K2), the quality of the fit to the observational constraints improved significantly, with the best results being obtained for A₂ ≃ 0.12–0.18. Modifying A₃ in addition to A₂ (as in model K23) was found to be helpful but not necessary. The χ² and rms(δc_s) values decreased in this case because of the additional degree of freedom. However, additional simulations with A₃ = ±0.05 and A₂ = 0.10 (model K23′) confirmed that the effect of A₃ was marginal. Therefore, we did not explore the effect of A₃ further and simulated instead the K2 models (with A₃ = 0).

Fig. 6. Similar to Fig. 5 but showing the dependence on A2 in Runs 9, 10, 11, 13, and 14; in addition, the left bottom panel shows A3. The orange and blue colors indicate the cases with varying A3 and A3 = 0, respectively. The circles and plus signs denote the cases with and without planet formation, respectively (i.e., Zacc either evolves with time or is a constant). The gray triangles at A2 = 0.1 denote the cases with A3 = ±0.05. The two gray vertical lines demarcate the opacity enhancement (i.e., 7 ± 3%) suggested by Bailey et al. (2015).

The iron opacities measured in laboratory experiments at a temperature of ∼2 × 10⁶ K should correspond to an increase of 7 ± 3% in the Rosseland-mean opacities over standard opacities (Bailey et al. 2015). This temperature is consistent with an increase in A₂ of similar magnitude (see Sect. 3.1.5 and Le Pennec et al. 2015). Based on this and the best-fit results from Fig. 6, we adopted A₂ = 0.12 for our fiducial model. We note that this value of A₂ is also consistent with that used in models with increased opacities in previous studies (Serenelli et al. 2009; Buldgen et al. 2019).

We observed a clear correlation between Z_proto and A₂, such that Z_proto decreased monotonically with A₂. This is because both parameters have the same effect on the opacity in the solar interior. The correlations between A₂ and the other parameters (i.e., α_MLT, f_overshoot, and Y_proto) are more complicated. This is because we used the simplex method to fit six observational constraints on which the input parameters affect differently (see, e.g., Henyey et al. 1965; Kippenhahn et al. 2012).

4.5. Opacity increase and metal-poor accretion

Finally, we performed simulations by including both opacity increase and metal-poor accretion (i.e., the Z_acc evolution as shown in Fig. 4). The blue and orange circles in Fig. 6 show the results of models K2-MZvar (with A₃ = 0) and K23-MZvar (with varying A₃), respectively, for A₂= 0.12, 0.15, and 0.18. We found that the minimized χ² and rms(δc_s) values were almost the same as those obtained using models K2 and K23 (i.e., with homogeneous Z_acc), as expected from the results presented in Sect. 4.3. This again points to the conclusion that planet formation processes do not affect the solar sound speed profile.

4.6. The solar abundance problem

It is interesting to note that the K2 model with A₂ = 0.12 fits the observational constraints, and in particular, rms(δc_s), even better than the models using GS98 abundances (see Fig. 6). This is better demonstrated in Fig. 7, where we plot δc_s as a function of radius. For our fiducial model (i.e., model K2 with A₂ = 0.12), the peak value at the base of the CZ is δc_s = 0.003, which is much smaller than the peak value obtained for standard models with the AGSS09 composition (δc_s = 0.009), and even better than that obtained for models with the GS98 composition.

Fig. 7. Radial profiles of observed minus calculated sound speed δcs. The blue solid line with circles shows the result obtained in this study using model K2 with A2 = 0.12. The gray solid and dashed lines show the results of non-accreting models with GS98 and AGSS09 compositions, respectively (see Fig. A.1). The other solid lines show the results in the literature: the orange, red, purple, and green lines show the results of: Christensen-Dalsgaard et al. (2018, see their Fig. 7), Buldgen et al. (2019, see the magenta line in their Fig. 9), Model AGSSr2a of Yang (2019), and Model TWA of Zhang et al. (2019).

Figure 7 also compares our model with other models in the literature that include various physical processes. These models are qualitatively equivalent fits to the sound speed constraint. Zhang et al. (2019) claimed the importance of helium-poor accretion (however, see discussion in Sect. 4.2) and the improved overshooting models (see Sect. 3.1.3). In the models by Christensen-Dalsgaard et al. (2018) and Buldgen et al. (2019), both an ad hoc opacity increase and extra mixing around the base of the CZ were considered. Yang (2019) showed that rotational mixing is also promising. We note that in the present work, we did not investigate these possibilities, but Fig. 7 shows that all the aforementioned models also fit the helioseismic constraints better than the models with the high-Z GS98 abundances.

We note that the smoothness of our δc_s profile differs from that of other models in the literature. We derived δc_s simply by comparing our simulated c_s profile with that given in Basu et al. (2009) (see Sect. 3.2). Conversely, in some of the studies, δc_s was derived by comparing the oscillation modes from the calculated solar structure using an inversion method with the observed modes (see, e.g., Buldgen et al. 2019). We already confirmed that our δc_s profiles for the non-accreting models are in good agreement with those of previous studies, and therefore, the difference in the derivation of δc_s does not change the conclusions of this study.

5. Consequences of planet formation

Using our fiducial K2 model, shown in Sect. 4 to be our best fit to the observational constraints, we now examine the impact of an inhomogeneous Z_acc on the following: present-day solar central metallicity Z_center (see Sect. 5.2), primordial Solar-System metallicity Z_proto (see Sect. 5.3), and primordial Solar-System helium abundance Y_proto (see Sect. 5.4). Before discussing Z_center, Z_proto, and Y_proto, in Sect. 5.1, we first discuss the consequences of planet formation on the solar interior, as inferred from the results presented in Sect. 4.

5.1. Planet formation and the solar interior

As discussed in Sect. 2, planet formation is associated with both an increase in Z_acc (i.e., the “pebble wave” phase) and a late metal-poor accretion phase during which planetesimals or planets are formed. The largely convective structure of the young Sun partly erases the signature of planet formation; however, a small trace remains and can be highlighted by comparing solar interior models calculated with and without planet formation mechanisms.

Figure 8a shows the metallicity profile in the interior of the present-day Sun for the models K2-MZvar and K2 (i.e., with and without a time-varying Z_acc, respectively) with A₂ = 0.12. We observe that planet formation increases the central metallicity Z_center by 9 × 10⁻⁴ (≃5%); however, in the outer part of the Sun (r ≳ 0.2 R_⊙), the metallicity is almost the same. This is because, as shown in Fig. 2, the proto-Sun at ≤10 Myr has a small radiative core (i.e., a non-mixing region) and the signature of an initial high metallicity remains there until the present day. Therefore, realistic protostellar and pre-MS evolution models are crucial for obtaining realistic Z_center values.

Fig. 8. Comparison between models K2 (gray dashed line) and K2-MZvar (blue solid line) with A2 = 0.12. Left panel: metallicity profile in the interior of present-day Sun. Right panel: time evolution of the surface metallicity Zsurf. The gray vertical lines denote M1(= 0.904 M⊙), M2(= 0.962 M⊙), and the end of the accretion phase for model K2-MZvar. The target value of Zsurf at 4.567 Gyr is 0.0134 ± 0.0007 (see Table 3).

Figure 8b shows the time evolution of Z_surf. For model K2-MZvar with A₂ = 0.12, Z_surf varies considerably from 0.0147 to 0.0162 in the accretion phase because of planet formation; however, it decreases monotonically in the MS phase due to gravitational settling. We note that Z_surf decreases slightly with time in the very early phase (e.g., Z_surf decreases from 0.0149 to 0.0147 during the time 0.1 to 0.7 Myr). This is because the 0.1 M_⊙ initial stellar seed has Z = 0.02, which is larger than Z_proto = 0.0140.

To evaluate the effects of planet formation on the properties of the solar interior, we examined models having different values of the three parameters, M₁, M₂, and ΔZ_acc ≡ Z_acc,max − Z_proto (see Fig. 4), as listed in Table 4. We refer to these new simulation models as K2-MZ. We set the fiducial t_acc for these models to be 5 Myr because the typical protoplanetary disk lifetime (i.e., half-life period) is several million years (see Sect. 2.2). We set t_acc = 10 Myr in Sect. 4 to investigate the maximum impact of planet formation on the c_s profile. In contrast, in this section, our aim is to explore the realistic extent of the influence of planet formation. In Sects. 5.2–5.4 we evaluate the effects of planet formation on the Z_center, Z_proto, and Y_proto values.

It is important to link the parameters M₁, M₂, and ΔZ_acc to M_Z, planet, which is the mass retained by planets in the form of “metals” (i.e., excluding hydrogen and helium but including all the other elements) and to M_XY, lost, which is the mass in metal-poor gas (i.e., containing hydrogen and helium only) that is lost due to selective photoevaporation and/or MHD disk winds (see Sect. 2.1). The conservation of the metal mass implies that

(9)

We note that in the above equation, M_XY, lost and M_Z, planet are degenerate, and hence, we assumed M_Z, planet = 150 M_⊕ following Kunitomo et al. (2018), unless otherwise noted. Table 4 also lists the values obtained for M_XY, lost using the optimized Z_proto and Eq. (9). If we obtained M_XY, lost < 0, we set M_XY, lost = 0 and consequently increased M_Z, planet. We note that the value of M_Z, planet in the Solar System is rather uncertain, but we estimate the upper limit to be ≃168 M_⊕ (Kunitomo et al. 2018).

5.2. Impact on the central solar metallicity

To what extent can planet formation modify Z_center? To address this question, we focused on the results of the K2-MZ models (see Sect. 5.1 and Table 4) in addition to those of the K2 and K2-MZvar models presented in Sect. 4. We explore in Fig. 9 the dependence of Z_center on the parameters A₂, t_acc, f_overshoot, and planet formation scenarios.

Fig. 9. Dependence of Zcenter on (a) A2, (b) tacc, (c) fovershoot, and (d) planet formation models (corresponding to different values of M1, M2, and Zacc,max; see Table 4). In panel a, the results of models K2 and K2-MZvar are shown, in addition to those of models K2′, K2-MZ1′, and K2-MZ8′, with fovershoot = 0.01 and tacc = 10 Myr (see Table 1). In panels b and c, the green and orange curves indicate models K2′ and K2-MZ1′, respectively. The orange dots in panel d are the K2-MZ simulation results. The two horizontal gray dashed lines in all the panels denote the results of the non-accreting models without overshooting for the GS98 (top) and AGSS09 (bottom) compositions. The horizontal green dashed line in panel d denotes the result of model K2′ with A2 = 0.12, tacc = 5 Myr, and fovershoot = 0.01.

Figure 9a shows that Z_center is negatively correlated with A₂; moreover, Z_center increases by ≃0.008 because of planet formation processes, irrespective of the value of A₂. This negative correlation arises owing to the negative correlation between Z_proto and A₂ (see Sect. 4.4). We note that, in the K2 and K2-MZvar models, f_overshoot, M₁, M₂, and ΔZ_acc were allowed to vary in order to minimize χ². We performed additional simulations (K2′, K2-MZ1′, and K2-MZ8′ with t_acc = 10 Myr; see Table 1) by fixing these parameters. The M₁, M₂, and ΔZ_acc values of K2-MZ1′ and K2-MZ8′ are the same as those of K2-MZ1 and K2-MZ8, respectively (see Table 4). The results of the different models were quite similar, confirming that the final Z_center value is mostly determined by the opacity (A₂) of the solar interior and planet formation processes (in particular the metal-poor accretion).

Next, we investigated the dependence of Z_center on t_acc (see Sect. 2.2). Figure 9b shows the results of models K2-MZ1′ and K2′ for f_overshoot = 0.01 and A₂ = 0.12 with different t_acc values. We find that for model K2-MZ1′, Z_center increases with t_acc. This is because the size of the radiative core in the accretion phase increases with time (Sect. 2.2); thus, in the longer t_acc case, the metal-poor accretion (see Fig. 4) has a larger effect and Z_proto (and therefore Z_center) can be higher. However, for model K2′ with a homogeneous Z_acc, Z_center is independent of t_acc. Although Serenelli et al. (2011, see their Figs. 4 and 7) and Zhang et al. (2019, see their Figures 15 and 18) have already shown that Z_acc has the potential to modify Z_center, the accretion history in these studies (and t_acc of Serenelli et al. 2011) is not based on the standard model of star formation (see Sect. 3.1.2). Considering the strong dependence of Z_center on t_acc, a realistic accretion history is crucial for accurately evaluating Z_center.

Figure 9c shows that Z_center decreases with f_overshoot for the model with planet formation. This is because if a star has a more vigorous overshooting, then the radiative core becomes smaller and is developed in a slightly later phase. However, for the cases with homogeneous Z_acc, Z_center is insensitive to f_overshoot because the internal metallicity profile is homogeneous until gravitational settling sets in (i.e., in ∼1 Gyr). Although recent studies have attempted to constrain the efficiency of overshooting by hydrodynamic simulations (Freytag et al. 1996; Korre et al. 2019; Higl et al. 2021) and observations (e.g., Deheuvels et al. 2016), it remains uncertain. Further studies to constrain the mixing in stellar interiors are necessary.

Finally, Fig. 9d shows Z_center for the K2-MZ models. We find that Z_center is not sensitive to planet formation scenarios. This is related to the process used to satisfy the constraints in this work. We fixed M₁, M₂, and ΔZ_acc, and chose Z_proto such that it satisfied the constraints (in particular, (Z/X)_surf; see Fig. A.4b) imposed by the simplex method. Consequently, Z_center does not have much variation with respect to different planet formation scenarios.

We note that in models with a higher Z_center, the central opacity is higher, leading to an increased radiative temperature gradient (e.g., Kippenhahn & Weigert 1990), and thus, to a higher central temperature. For models with A₂ = [0.12, 0.18] and A₃ = 0, the central temperature ranges from 1.559 × 10⁷ to 1.567 × 10⁷ K for values of Z_center in the range 0.0155–0.0171.

In summary, the value of Z_center is larger for a lower A₂, longer t_acc, and lower f_overshoot. For A₂ = 0.12 and t_acc = 5 Myr, planet formation enhances Z_center up to 0.01686. The Z_center values for models with A₂ = [0.12, 0.18] (i.e., the models that reproduce the sound speed profile) and planet formation are ≃5% higher than those of homogeneous Z_acc models (see Fig. 9a). Interestingly, recent observations of solar neutrino fluxes reaching the Earth have also pointed to high values of the metallicity at the solar center (Agostini et al. 2018; Borexino Collaboration 2020); however, their interpretation remains tentative at present. More accurate observations are required in the future.

5.3. Impact on the constraints for the primordial metallicity

Vinyoles et al. (2017) obtain constraints on the primordial metallicity Z_proto that mostly depend on the assumed abundance table from their non-accreting models, that is, Z_proto = 0.0149 for AGSS09 and Z_proto = 0.0187 for GS98. For non-accreting models, we obtain Z_proto = 0.0163 and Z_proto = 0.0187, respectively. We find that the difference in the AGSS09 case arises from the different number of target values, N. As described in Sect. 3.3, N = 3 in Vinyoles et al. (2017, and most previous studies), whereas N = 6 in this study, leading to a much better fit of Y_surf. We stress that these solutions, however, do not represent the best fits to the observational constraints.

Figure 10 shows the variation in the primordial metallicity of the protosolar molecular cloud core Z_proto (see Sect. 2.1) for the same models and parameters as in Fig. 9, including our best models. Figure 10a shows that Z_proto is negatively correlated with A₂, as described in Sect. 4.4. Although this behavior is the same as that of Z_center, Fig. 10a shows some differences compared to Fig. 9a. Most importantly, while planet formation processes always increase Z_center, they can either increase or decrease Z_proto. As shown by Figs. 10b and c, this cannot be explained by varying t_acc or f_overshoot, which have little impact on Z_proto.

Fig. 10. Same as Fig. 9 except the vertical axis is the initial metallicity Zproto, which corresponds to the primordial metallicity of the protosolar molecular cloud core; and the horizontal axis of (d) is MXY, lost, which was calculated using Eq. (9) and MZ, planet = 150 M⊕. We note that the unphysical MXY, lost < 0 values are due to our assumption that MZ, planet = 150 M⊕. These solutions are equivalent to solutions with MXY, lost = 0 and MZ, planet > 150 M⊕ (see Table 4).

Figure 10d shows that Z_proto decreases with M_XY, lost. When the mass of the hydrogen and helium that is selectively removed from the disk is high, the primordial metallicity Z_proto must be low to compensate for the accretion of the disk gas that becomes relatively metal-rich and to account for the present-day observations. Conversely, in models with a low value of M_XY, lost, Z_proto must be high to compensate for the metal-poor accretion. We note that the retention of heavy elements by planet formation has the opposite effect (see Eq. (9)) so that a low value of M_XY, lost is equivalent to a high value of M_Z, planet (we set M_Z, planet = 150 M_⊕ in Fig. 10d; see Sect. 5.1). If M_XY, lost ≃ 0.03 M_⊙, then Z_protoM_XY, lost compensates for M_Z, planet = 150 M_⊕ and the value of Z_proto is the same as that obtained for model K2′ (i.e., with homogeneous Z_acc).

The above behavior can also explain the non-monotonic relation between Z_proto and A₂ for the K2-MZvar models in Fig. 10a. Indeed, the models with A₂ = 0.12, 0.15, and 0.18 correspond to M_XY, lost = 0.13, 0.18, and 0.13 M_⊙, which explains the much lower Z_proto value for the K2-MZvar model with A₂ = 0.15.

For our preferred models with A₂ = [0.12, 0.18], which successfully reproduce the observational constraints (see Sect. 4.4), Z_proto ranges from 0.0127 to 0.0157. A slightly narrower range may be obtained by imposing tighter constraints on M_XY, lost and M_Z, planet, namely, M_XY, lost ≲ 0.05 M_⊙ and M_Z, planet = 97–168 M_⊕ (see Sect. 2.1). In this case, Z_proto lies approximately in the range 0.0140–0.0153.

5.4. Impact on the constraints for the primordial helium abundance

Estimating the primordial helium abundance (i.e., the helium abundance in the protosolar molecular cloud core Y_proto) accurately is of particular interest, because it provides reliable constraints on the internal structure and composition of Jupiter and Saturn (Guillot et al. 2018). Serenelli & Basu (2010) obtained a primordial helium abundance of 0.273 ± 0.006 by using solar evolution models with turbulent mixing below the CZ. Vinyoles et al. (2017) obtained a lower value of Y_proto = 0.2613 for the AGSS09 composition, which resulted in a poor fit of Y_surf (see Sect. 3.3).

Figure 11 shows Y_proto for the same models as in Figs. 9 and 10. Again, the most significant correlation is obtained between Y_proto and the opacity-enhancement parameter A₂. We observe that Y_proto decreases from 0.274 for the standard model (K2 model with A₂ = 0) to 0.267 for the highest value of A₂ = 0.22. Figures 11b–d show that Y_proto has a weak dependence on t_acc, f_overshoot, and M_XY, lost (for the models with planet formation).

Fig. 11. Same as Fig. 10 except the vertical axis is the initial helium mass fraction Yproto, which corresponds to the primordial helium abundance of the protosolar molecular cloud core. The gray open circles in panel d show the results of models K2, K2′, K2-MZvar, K2-MZ1′, and K2-MZ8′ with A2 = [0.12, 0.18].

However, Fig. 11a shows that planet formation processes result in an increase in Y_proto, independent of M_XY, lost. This is similar to what is observed for Z_center but different from what is observed for Z_proto. The reason for such behavior is threefold. First, hydrogen and helium are believed to have a common evolutionary history that differs from that of metals (see Sect. 3.1.6). This implies that planet formation processes (i.e., pebble wave, planetesimal formation, and disk winds) affect Z_acc while conserving the mass ratio of helium to hydrogen (Y/X) in the accreted material. Second, when planet formation results in a high Z_center, the central temperature T_center of the star increases because of the changes in the opacity, which in turn increases the radiative temperature gradient (see Sect. 5.2).

Third, this increase in temperature affects the nuclear burning rate r_pp because (Kippenhahn & Weigert 1990), where X_center ∼ 1 − Y_center is the hydrogen abundance in the Sun’s nuclear burning core. Given the global constraints on the structure of the present-day Sun, to compensate for the higher T_center, X_center needs to be lower on average on the MS, which in turn implies a lower X_proto and therefore a higher Y_proto. For these reasons, planet formation processes lead to higher values of Z_center, T_center, and hence, Y_proto.

In addition, Y_proto is positively correlated with A₃. Therefore, if the high-temperature opacities are changed, then Y_proto would be modified by an absolute factor, which is estimated to be δY_proto ∼ A₃/20.

Overall, our results for the models that fit the observational constraints (i.e., A₂ = [0.12, 0.18]) imply that Y_proto ranges from 0.268 to 0.274. This range compatible with, but somewhat more tightly constrained than the range obtained by Serenelli & Basu (2010), which is 0.267–0.278. The Y_proto/(1 − Z_proto) value was found to range from 0.272 to 0.278.

6. Conclusion

In this work, we studied how the formation of the Solar System affected the composition and internal structure of the Sun. The Sun was formed owing to the collapse of a molecular cloud core, and it grew via accretion of the gas in the circumsolar disk over several million years while planet formation processes were at play. According to protoplanetary disk evolution models, when the proto-Sun was approximately between 90% and 98% of its final mass (see Fig. 1), a pebble wave led to a phase of accretion of high-metallicity gas. This was followed by a phase of metal-poor accretion due to the formation of planetesimals and planets, while the hydrogen and helium may have been selectively lost from the disk atmosphere via photoevaporation and/or MHD disk winds. Therefore, the Sun grew via the accretion of gas with an evolving composition.

To study the evolution of the Sun and reproduce the present-day constraints on its structure and composition obtained from spectroscopy and helioseismology, we performed an extensive ensemble of simulations. Our simulations included accretion in the protostellar and pre-MS phases. The input parameters were adjusted using χ² minimization to best reproduce the present-day constraints (i.e., after 4.567 Gyr of evolution) on the luminosity, effective temperature, surface composition (Y_surf and (Z/X)_surf), CZ radius (R_CZ), and sound speed profile of the Sun. The input parameters used were the mixing length parameter α_MLT, overshooting parameter f_overshoot, initial helium abundance Y_proto, and initial metallicity Z_proto. A second set of adjustable parameters were introduced to modify the opacity, namely, A₁, A₂, and A₃. A third set of parameters were introduced to modify the composition of the accreted material, M₁, M₂, Z_acc,max, and Y_acc,min.

Several scenarios were tested. Classical non-accreting models with old abundances (GS98) are known to provide better fits to the helioseismic constraints than those using more recent abundances (AGSS09), leading to the so-called “solar abundance problem.” Models that consider the accretion of gas with an evolving composition (i.e., include planet formation processes) do not improve the χ² fit significantly. This is because the proto-Sun has an almost fully convective interior in the accretion phase, implying that the accreted gas is heavily diluted and therefore the changes in the structure of the present-day Sun are quite limited. Models involving helium-poor accretion, as proposed by Zhang et al. (2019), indeed improve the fit. However, such models are probably unlikely because the giant planets in our Solar System contain significant amounts of helium in their interiors. We note that other possibilities, such as extra mixing, were not tested in this work (see Christensen-Dalsgaard et al. 2018; Buldgen et al. 2019; Yang 2019). Our best models were found to be those with the AGSS09 composition and a 12%–18% opacity increase centered at T = 10^6.4 K. This is slightly higher but qualitatively in good agreement with the high iron opacities measured by Bailey et al. (2015) at this temperature range. The models with an opacity increase in the range 12%–18% represent better fits to the observations than those using old abundances, and are therefore a promising solution to the solar abundance problem.

The impact of planet formation on the solar structure can be examined by using the aforementioned models to globally fit the observational constraints and by modifying the parameters M₁, M₂, and Z_acc,max. We find that despite the negligibly small effect on the sound speed profile (and therefore on the helioseismic constraints), planet formation processes lead to a limited but real (up to 5%) increase in the metallicity Z_center of the deep solar interior (r ≲ 0.2 R_⊙). The increase in Z_center is smaller for a larger overshooting parameter f_overshoot but larger for a longer disk-accretion timescale t_acc. Qualitatively, this is in accordance with recent solar neutrino measurements that appear to favor high-Z values in the deep solar interior (Agostini et al. 2018; Borexino Collaboration 2020). We will examine this issue in future work.

We also constrained the primordial composition of the molecular cloud core that gave birth to the Sun by using models that best reproduced all the present-day constraints. We found that the protosolar metallicity Z_proto ranged from 0.0127 to 0.0157; moreover, the retention of heavy elements due to planet formation yielded higher values of Z_proto, while significant disk winds led to lower values. Similarly, the protosolar helium mass fraction Y_proto was found to slightly increase because of planet formation processes and its value ranges from 0.268 to 0.274.

In conclusion, we expect that a combined investigation of the solar interior, solar evolution models (in particular, those that include improved opacities and a more sophisticated treatment for mixing), and observational constraints (from surface abundances, helioseismic observations, and neutrino fluxes) will help constrain the processes that led to the formation of the Solar System.

Acknowledgments

This paper is dedicated to the memories of Rudolf Kippenhahn, whose diagrams were a source of inspiration for this work, and Jean-Paul Zahn, who taught one of us how to read them (upside-down!). We are grateful to Gaël Buldgen, Jørgen Christensen-Dalsgaard, Johan Appelgren, and Vardan Elbakyan for kindly providing their simulation results. We also thank Lionel Bigot for fruitful discussions and comments. This work was supported by the JSPS KAKENHI Grant (number 20K14542), the Astrobiology Center Program of National Institutes of Natural Sciences (NINS; grant number AB301023), and a long-term visitor JSPS fellowship to TG. Numerical computations were carried out on the PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan, and in part on the “Mesocentre SIGAMM” machine hosted by the Observatoire de la Cote d’Azur. Software: MESA (version 12115; Paxton et al. 2011, 2013, 2015, 2018, 2019); Numpy (van der Walt et al. 2011); and Scipy (Virtanen et al. 2020).

References

Appendix A: Details of simulation results

In this Appendix, we provide the details of our simulation results. Tables A.1 and A.2 show the input parameters and results at 4.567 Gyr, respectively, which were optimized using the simplex method.

Figure A.1 presents a comparison of the δc_s profiles obtained by previous studies and the present study for the non-accreting models with the GS98 and AGSS09 compositions. As mentioned in Sect. 4.1, our results agree well with those of previous studies.

Fig. A.1. Radial profiles of observed minus calculated sound speed δcs for the non-accreting models with the GS98 (left panel) and AGSS09 (right) compositions. The blue solid and dashed lines with circles indicate the results with and without overshooting, respectively. The other solid lines show the results in the literature: the red, green, and brown lines indicate the results of Zhang et al. (2019, Model SSM98 for the left panel and SSM09 for the right panel), Buldgen et al. (2019, Models “GS98-Free-OPAL” and “AGSS09-OPAL-OPAL”; see their Figure 1), and Serenelli (2013, see their Figure 1).

Figure A.2 shows the δc_s profiles for non-accreting models with the AGSS09 composition for different settings of the opacity and nuclear reaction tables. Only marginal differences exist between the different models (see also Serenelli et al. 2011; Buldgen et al. 2019; Salmon et al. 2021).

Fig. A.2. Same as Fig. A.1(b) (i.e., non-accreting models with the AGSS09 composition). The pink and gold lines indicate the cases with the OP opacity table (with overshooting) and JINA nuclear reaction table (without overshooting), respectively. The thin blue lines are the same as the thick blue lines in Fig. A.1(b).

Figure A.3 shows the metallicity profiles at 4.567 Gyr. Since the observed (Z/X)_surf value for the GS98 composition is higher than that for the AGSS09 composition, the corresponding Z_center is also higher. As shown in Sect. 5, a metal-rich structure is obtained for the case with metal-poor accretion (see model MZvar).

Fig. A.3. Radial metallicity profiles in the solar interior at 4.567 Gyr for non-accreting models (left panel) and accreting models (right panel). The red and gray lines denote the cases with the GS98 and AGSS09 compositions, respectively (models noacc-GS98-noov, noacc-GS98, noacc-noov, and noacc, from top to bottom). The target values of Zsurf at 4.567 Gyr are 0.0168 ± 0.0007 and 0.0134 ± 0.0007, respectively (see Table 3). The brown and purple lines in the right panel show the results of models MZvar and He20Myr, respectively (see Table 1 and Fig. 5).

Figure A.4 shows the metallicity profiles at 4.567 Gyr and the time evolution of Z_surf for the K2-MZ simulations with different Z_acc models (see Table 4). Although there exists a substantial difference in the Z_surf evolution in the accretion phase, it converges in the MS phase to meet the observed constraints. The δc_s profiles are also shown in Fig. A.5. As expected, the different Z_acc models do not affect the δc_s profile (see Sect. 4.3).

Fig. A.4. Results of the K2-MZ model with different Zacc (models MZ1–MZ9; see Tables 1 and 4). The left and right panels show the radial metallicity profile in the solar interior at 4.567 Gyr and the time evolution of the surface metallicity Zsurf, respectively.

Fig. A.5. Radial sound speed profiles of the K2-MZ model with different Zacc; the lines follow the same color-coding as in Fig. A.4. The gray solid and dashed lines show the results of the non-accreting models with overshooting and with AGSS09 and GS98 abundances (i.e., noacc and noacc-GS98 in Table 1), respectively.

Figures A.6(a) and (b) show various evolutionary tracks on the Hertzsprung-Russell diagram and time evolution of the mass coordinate of the base of the CZ M_BCZ, respectively. We found that the tracks are insensitive to A₂, t_acc, Z_acc, and Y_acc. In Fig. A.6, we also show the tracks obtained by Kunitomo et al. (2018, calculated with a 0.01 M_⊙ seed using the MESA code version 6596; see their Table C.1)⁸. As described in Sect. 3.1.1, our initial condition corresponds to a high-entropy 0.1 M_⊙ stellar seed. Therefore, the initial location of the evolutionary track (i.e., log T_eff = 3.41 and log L_⋆ = −0.2) coincides with the track for ξ = 0.5 obtained by Kunitomo et al. (2018). Because we adopted ξ = 0.1 in this work, our tracks deviate from the one with ξ = 0.5, and the difference in luminosity can be up to ≃1 dex. However, after the stellar mass reaches ≃1 M_⊙ and they reach their Hayashi track, the tracks agree well. We note that one finds a small offset in T_eff, which arises due to the switching of the outer boundary conditions (see Section 2.7 of Kunitomo et al. 2018). Figure A.6(b) shows that the evolution of M_BCZ is also insensitive to A₂, t_acc, Z_acc, and Y_acc, and is almost same as that for the case with ξ = 0.5 obtained by Kunitomo et al. (2018).

Fig. A.6. Evolutionary tracks on the Hertzsprung-Russell diagram (left panel) and time evolution of the mass coordinate of the base of the CZ MBCZ (right), for models K2 with A2 = 0 (blue line), K2-MZvar with A2 = 0.12 (orange line), K2-MZ1′ with fovershoot = 0.01 and tacc = 3 Myr (green line), and He12Myr (red line) (see Table 1). The filled circles indicate the times at which accretion ends. The open circle in the left panel indicates the position of the Sun. The gray dashed lines show the evolutionary tracks obtained by Kunitomo et al. (2018, for ξ = 0.5, 0.1, and 0).

Additional supplemental materials are available at the CDS (see the links provided in the footnote on the first page). These include a csv file, summarizing the optimized input parameters and results of all the simulation models, and structure and evolution data of the optimized cases of models K2, K2-MZ5, and K23-MZvar with A₂ = 0.12. In addition, plots showing how the input parameters and results are changed with runs by the simplex method, animations of the evolutions of the optimized case of each model, and the input files of the MESA code are also available at the Zenodo. In the csv file, we list the values of M_XY, lost [M_⊙], assuming M_Z, planet = 150 M_⊕, and M_Z, planet [M_⊕], assuming M_XY, lost = 0, for the models with planet formation. Neutrino fluxes, calculated using the subroutine provided in Farag et al. (2020), are also included in the csv file.

All Tables

All Figures

Fig. 1. Evolution of the metallicity Zacc of the accreted gas in units of the initial metallicity Zproto as a function of time (left panel) and mass of the proto-Sun (right panel). Several evolutionary models of circumstellar disks that include both gas and dust are shown (see text for details). The protostar mass was calculated assuming some mass loss (e.g., due to partial photoevaporation of the circumstellar disk) except in two cases, which are represented by dashed lines.

In the text

Fig. 2. Early evolution of the solar interior. The total mass of the accreting proto-Sun is indicated by the dashed line. The surface CZ is depicted as a cloudy region (see Kippenhahn & Weigert 1990) delimited by a radiative zone that grows outward from the center of the Sun. Three evolution models are shown corresponding to different values of the accretion efficiency parameter ξ: a classical evolution model with ξ = 0.5 (red line), a cold accretion model with ξ = 0.1 (green line) that is a lower limit of the ξ value in the observations of young clusters (Kunitomo et al. 2017), and a model with ξ = 0 (gray line) corresponding to a theoretical proto-Sun formed in the absence of any accretion heat. The right-hand y-axis provides the radius corresponding to the mass on the left-hand y-axis of the present-day Sun (which is 4.567 Gyr old). Bottom panel: highlights the key physical processes that occur during the growth of the proto-Sun: the presence of a circumstellar gas disk (during the first million years); an increase in the metallicity Zacc of the accreted gas due to a pebble wave; the concomitant formation of planetesimals and planets; and a sudden decrease in Zacc.

In the text

	Fig. 3. Evolution of stellar mass M⋆ (red solid line) and accretion rate Ṁacc (green dashed line) with time in the fiducial case. The gray vertical lines indicate t1 and tacc in the fiducial case. The gray dotted lines indicate M⋆ and Ṁacc for tacc = 3 and 20 Myr.
In the text

	Fig. 4. Sketch of the evolution of the metallicity Zacc of the accreted material with stellar mass M⋆. In this work, we used the model depicted by the blue solid line and assumed that the difference between this model and the one depicted by the green dotted line has little impact on the internal structure of the Sun (see text for details).
In the text

Fig. 5. Optimized input parameters (left) and results at the solar age (right). The left panels show αMLT, fovershoot, Yproto, Zproto, and stellar mass (M1 and M2) for nine models from top to bottom. The right panels show the χ2 value (see Eq. (8)), rms(δcs), (Z/X)surf, Ysurf and RCZ from top to bottom. The five models shown in each panel are: non-accreting models with GS98 composition (Runs 2 and 4); non-accreting models with AGSS09 composition (Runs 1 and 3); accreting models with time-dependent Zacc (Runs 7 and 8); models with He-poor accretion (Runs 5 and 6); and K2 model (Run 9) with A2 = 0.12, from left to right (see Table 1 for details). Squares and crosses represent the cases with and without overshooting, respectively. Diamonds and asterisks represent the cases with He-poor accretion for tacc = 12 and 20 Myr, respectively. The shaded regions in the right panels indicate χ2 < 1 or the 1σ uncertainties (see Table 3).

In the text

Fig. 6. Similar to Fig. 5 but showing the dependence on A2 in Runs 9, 10, 11, 13, and 14; in addition, the left bottom panel shows A3. The orange and blue colors indicate the cases with varying A3 and A3 = 0, respectively. The circles and plus signs denote the cases with and without planet formation, respectively (i.e., Zacc either evolves with time or is a constant). The gray triangles at A2 = 0.1 denote the cases with A3 = ±0.05. The two gray vertical lines demarcate the opacity enhancement (i.e., 7 ± 3%) suggested by Bailey et al. (2015).

In the text

Fig. 7. Radial profiles of observed minus calculated sound speed δcs. The blue solid line with circles shows the result obtained in this study using model K2 with A2 = 0.12. The gray solid and dashed lines show the results of non-accreting models with GS98 and AGSS09 compositions, respectively (see Fig. A.1). The other solid lines show the results in the literature: the orange, red, purple, and green lines show the results of: Christensen-Dalsgaard et al. (2018, see their Fig. 7), Buldgen et al. (2019, see the magenta line in their Fig. 9), Model AGSSr2a of Yang (2019), and Model TWA of Zhang et al. (2019).

In the text

Fig. 8. Comparison between models K2 (gray dashed line) and K2-MZvar (blue solid line) with A2 = 0.12. Left panel: metallicity profile in the interior of present-day Sun. Right panel: time evolution of the surface metallicity Zsurf. The gray vertical lines denote M1(= 0.904 M⊙), M2(= 0.962 M⊙), and the end of the accretion phase for model K2-MZvar. The target value of Zsurf at 4.567 Gyr is 0.0134 ± 0.0007 (see Table 3).

In the text

Fig. 9. Dependence of Zcenter on (a) A2, (b) tacc, (c) fovershoot, and (d) planet formation models (corresponding to different values of M1, M2, and Zacc,max; see Table 4). In panel a, the results of models K2 and K2-MZvar are shown, in addition to those of models K2′, K2-MZ1′, and K2-MZ8′, with fovershoot = 0.01 and tacc = 10 Myr (see Table 1). In panels b and c, the green and orange curves indicate models K2′ and K2-MZ1′, respectively. The orange dots in panel d are the K2-MZ simulation results. The two horizontal gray dashed lines in all the panels denote the results of the non-accreting models without overshooting for the GS98 (top) and AGSS09 (bottom) compositions. The horizontal green dashed line in panel d denotes the result of model K2′ with A2 = 0.12, tacc = 5 Myr, and fovershoot = 0.01.

In the text

Fig. 10. Same as Fig. 9 except the vertical axis is the initial metallicity Zproto, which corresponds to the primordial metallicity of the protosolar molecular cloud core; and the horizontal axis of (d) is MXY, lost, which was calculated using Eq. (9) and MZ, planet = 150 M⊕. We note that the unphysical MXY, lost < 0 values are due to our assumption that MZ, planet = 150 M⊕. These solutions are equivalent to solutions with MXY, lost = 0 and MZ, planet > 150 M⊕ (see Table 4).

In the text

	Fig. 11. Same as Fig. 10 except the vertical axis is the initial helium mass fraction Yproto, which corresponds to the primordial helium abundance of the protosolar molecular cloud core. The gray open circles in panel d show the results of models K2, K2′, K2-MZvar, K2-MZ1′, and K2-MZ8′ with A2 = [0.12, 0.18].
In the text

Fig. A.1. Radial profiles of observed minus calculated sound speed δcs for the non-accreting models with the GS98 (left panel) and AGSS09 (right) compositions. The blue solid and dashed lines with circles indicate the results with and without overshooting, respectively. The other solid lines show the results in the literature: the red, green, and brown lines indicate the results of Zhang et al. (2019, Model SSM98 for the left panel and SSM09 for the right panel), Buldgen et al. (2019, Models “GS98-Free-OPAL” and “AGSS09-OPAL-OPAL”; see their Figure 1), and Serenelli (2013, see their Figure 1).

In the text

	Fig. A.2. Same as Fig. A.1(b) (i.e., non-accreting models with the AGSS09 composition). The pink and gold lines indicate the cases with the OP opacity table (with overshooting) and JINA nuclear reaction table (without overshooting), respectively. The thin blue lines are the same as the thick blue lines in Fig. A.1(b).
In the text

Fig. A.3. Radial metallicity profiles in the solar interior at 4.567 Gyr for non-accreting models (left panel) and accreting models (right panel). The red and gray lines denote the cases with the GS98 and AGSS09 compositions, respectively (models noacc-GS98-noov, noacc-GS98, noacc-noov, and noacc, from top to bottom). The target values of Zsurf at 4.567 Gyr are 0.0168 ± 0.0007 and 0.0134 ± 0.0007, respectively (see Table 3). The brown and purple lines in the right panel show the results of models MZvar and He20Myr, respectively (see Table 1 and Fig. 5).

In the text

	Fig. A.4. Results of the K2-MZ model with different Zacc (models MZ1–MZ9; see Tables 1 and 4). The left and right panels show the radial metallicity profile in the solar interior at 4.567 Gyr and the time evolution of the surface metallicity Zsurf, respectively.
In the text

	Fig. A.5. Radial sound speed profiles of the K2-MZ model with different Zacc; the lines follow the same color-coding as in Fig. A.4. The gray solid and dashed lines show the results of the non-accreting models with overshooting and with AGSS09 and GS98 abundances (i.e., noacc and noacc-GS98 in Table 1), respectively.
In the text

Fig. A.6. Evolutionary tracks on the Hertzsprung-Russell diagram (left panel) and time evolution of the mass coordinate of the base of the CZ MBCZ (right), for models K2 with A2 = 0 (blue line), K2-MZvar with A2 = 0.12 (orange line), K2-MZ1′ with fovershoot = 0.01 and tacc = 3 Myr (green line), and He12Myr (red line) (see Table 1). The filled circles indicate the times at which accretion ends. The open circle in the left panel indicates the position of the Sun. The gray dashed lines show the evolutionary tracks obtained by Kunitomo et al. (2018, for ξ = 0.5, 0.1, and 0).

In the text