© ESO 2021

1 Introduction

In 2016, a terrestrial-type planet was discovered orbiting the nearest star to the Sun, Proxima Centauri (hereafter Proxima Cen), in its habitable zone (Anglada-Escudé et al. 2016). The detection of such a near target (the distance between Proxima Cen and the Sun is only ~1.3 pc or ~4.2 ly) opens new and exciting prospects and unique opportunities for follow-up observations, spectroscopy, and planetary characterisation. Proxima Cen b has been discovered by means of radial velocity and has a minimum mass of M sin i = 1.27 M_⊕. Suárez Mascareño et al. (2020) suggest a potentially even lower minimal mass of M sin i = 1.173 ± 0.086 M_⊕. The actual planetary mass may be well above that limit, since the inclination i of the system is unknown. Anglada et al. (2017) proposed an inclination of the system of 45° due to the possible detection of a dust belt, which is, however, not yet confirmed. Two companion planets were recently suggested but are still unconfirmed. A less massive planet with minimal mass of M sin i = 0.29 M_⊕ was suggested to exist at an even tighter orbit than Proxima Cen b, namely with an orbital period of 5.15 Earth days compared to Proxima cen b’s orbital period of 11.2 days (Suárez Mascareño et al. 2020). Damasso et al. (2020) instead suggest a second companion further out of the system with an orbital period of 5.21 yr for a planet with a minimum mass of M sin i = 5.8 ± 1.9 M_⊕, and they suggest an eccentricity of Proxima Cen b of 0.17 (Anglada-Escudé et al. 2016 already indicated a possible eccentricity below 0.3). The age of Proxima Cen is often assumed to be similar to α Cen A with an approximate age of ~5 Gyr (Bazot et al. 2016). However, Beech et al. (2017) suggested a model age of 7–8 Gyr to better represent Proxima’s observations, which would be possible if Proxima Cen had been captured by α Cen A and B (Feng & Jones 2018). In our models, we ran our simulations for the first 5 Gyr of Proxima Cen b’s evolution, but we do not claim that the planet could not in principle be older than that. Table 1 lists available observations and constraints for Proxima Cen and Proxima Cen b.

Although the information about the properties of Proxima Cen b is rather sparse at the moment, the available data led to a vivid discussion in the literature about possible evolution and properties of the planet (Barnes et al. 2016; Ribas et al. 2016; Zuluaga & Bustamante 2018), including possible climates and observability (Turbet et al. 2016; Meadows et al. 2018; Del Genio et al. 2019), taking into account the tremendous differences between a low-mass star such as Proxima and the Sun. Studies such as Barnes et al. (2016) and Ribas et al. (2016) presented a comprehensive analysis of Proxima Cen b and its possible evolution. Barnes et al. (2016) considered such processes as the evolution of the galactic environment of Proxima Cen, stellar radiation, volatile inventory of the planet, radiogenic evolution of planetary interiors, etc. Ribas et al. (2016) studied possible evolution scenarios of irradiation, rotation and volatile inventory from formation of the planet to the present. Del Genio et al. (2019) investigated possible climate scenarios including a surface ocean. The strong flares emitted from Proxima Cen not only challenge the potential surface habitability of the planet (Vida et al. 2019), but would – in the case that an atmosphere can prevail or be re-established by volcanic outgassing (Godolt et al. 2019) – also influence potential biosignatures in the atmosphere (Scheucher et al. 2020). The question of whether an atmosphere can exist ona planet in close orbit around an active star depends, among other things, on the surface gravity (influencing the escape velocity) and the atmosphere composition, both of which are unknown for Proxima Cen b. A CO₂ -dominated atmosphere,for example, might be able to survive on Proxima Cen b rather than in an Earth-like atmosphere (Johnstone et al. 2019), under which circumstances life might be able to exist on its surface and even cope with the high UV irradiation (Abrevaya et al. 2020). The question of whether Proxima Cen b can be habitable is therefore strongly linked to the evolution of its atmosphere, especially its composition and the outgassing fluxes from the interior, which are feeding the atmosphere over time, compensating for potential atmosphere losses into space (Noack et al. 2014; Godolt et al. 2019).

As part of this work, we investigated the possible influence that different heating sources (radioactive heating, induction heating, andtidal heating) may have had on the evolution of the planet. To calculate electromagnetic induction heating in the planet, we also account for the evolution of the stellar magnetic field. Due to the close proximity of the planet to its host star, the planet is embedded in a strong magnetic field in the course of its orbital motion, which varies depending on the stellar rotation and orbital phase. Kislyakova et al. (2017) showed that an inclined stellar magnetic dipole can lead to strong induction heating in the uppermost part of the mantle for planets orbiting late-orbit M dwarfs. Similarly, planets on an inclined orbit can experience strong induction heating, which in some cases leads to interior heating exceeding Earth’s present-day radioactive heating by several orders of magnitude, hence either melting the planetary mantle until a shallow magma ocean forms, or leading to extreme volcanic events on the planetary surface (Kislyakova et al. 2018), and to observable effects (Guenther & Kislyakova 2020; Kislyakova et al. 2019). Later, Kislyakova & Noack (2020) showed that induction heating can be an important driver of volcanic activity on massive rocky planets, which potentially makes it relevant to Proxima Cen b. Induction heating depends on the strength of the magnetic dipole of the star, orbital distance, and stellar rotation rate, and it is therefore strongest for close-in planets around fast-rotating stars with a high magnetic field strength. Although at present Proxima Cen exhibits a slow rotation period of ~ 83 days, it likely has rotated much faster in the past (Yadav et al. 2016). We therefore studied both the present conditions in the system as well as possible evolutionary tracks during Proxima Cen’s lifetime. The different methods used to model the planetary interior, induction heating, the rotational evolution of Proxima Cen, and the long-term evolution models are introduced in Sect. 2. The results are presented in Sect. 3 and discussed in Sect. 4, followed by a concluding section.

2 Methods

2.1 Predicted planet composition and interior structure

The composition of planets in the inner part of the accretion disc is expected to reflect the abundances of refractory elements with similar condensation temperatures as measured from stellar spectra (Dorn et al. 2015; Hinkel & Unterborn 2018), as is the case for the main planet-building elements: magnesium, iron, and silicon. Güdel et al. (2004) and Fuhrmeister et al. (2011) used optical and X-ray spectroscopy with VLT/UVES and XMM-Newton to obtain element abundances in the stellar spectrum during different stages of quiescence and large flares. Fuhrmeister et al. (2011) find slightly elevated Mg/Fe and Si/Fe values during flares, and element ratios are solar-like during quiescence. In contrast, Güdel et al. (2004) reported elevated Mg and Si values compared to solar abundances during quiescence, but solar-like values and even Mg- and Si-reduced abundances during flares. Within the errors of the observations, on average, Mg/Fe/Si ratios seem to be similar to solar abundances. Therefore, due to the high number of uncertainties in stellar abundances, we adopted an Earth-like Mg–Si–e–O-based mantle mineralogy and a core consisting of pure iron. However, to assess the influence of different iron abundances on planet composition and evolution, we varied the iron weight fraction from Earth-like (40 wt% iron) to an iron-depleted (20 wt%) and an iron-enriched (60 wt%) case. We note that values for an Earth-like iron content given in the literature (e.g. Wagner et al. 2011) typically refer only to the core-mass fraction and are therefore below 40 wt%, whereas we consider that some of the iron can also be stored in mantle minerals (such as Fayalite and Wüstite to give only two examples). We assume here an Earth-like magnesium number (Mg# = 0.9) for the ratio of magnesium-to-iron end-member minerals.

The abundance of volatile elements, on the other hand, cannot be predicted by observations of the stellar composition since the various volatile species condense at different locations within the disc, far away from each other and from where rocky and metal-rich material condenses. A planet forming in (or migrating into) the inner rockier part of the accretion disc may either contain minimal to no amounts of volatile material (depending on where the accreted material came from, i.e. inward scattering from the outer disc), or may be composed of large amounts of water and other volatile material if formed outside the snow line. With an orbital period of 11.2 days, Proxima Cen b lies in the rock-forming, inner part of the system. In our study, we therefore assumed that the amount of volatile material delivered to the planet is in the Earth range and that Proxima Cen b is a rocky planet.

Radial velocity measurements (Anglada-Escudé et al. 2016) can only give us a minimum planetary mass, as long as no additional information is available to constrain the inclination of the star–planet system with respect to our observational plane. No transit of the planet has been observed so far, leading to the conclusion that the inclination is likely below 90°, and therefore the actual mass of the planet should be larger than 1.27 M_⊕. Anglada et al. (2017) proposed the existence of a dust belt at an orbit of about 30 AU from the star, which, if confirmed, would suggest an inclination of the system of 45° with respectto the plane of the sky. Since confirmation has not been possible to date, we assumed four different inclinations spanning the wide range of likely inclinations and hence actual planet masses: 90, 60, 45, and 30°. The corresponding planet masses would then be 1.27, 1.47, 1.8, and 2.54 M_⊕ (see also Table 2). The table lists in addition the predicted planet and core radii determined with our interior structure model (Noack et al. 2017) as well as the assumed initial temperature at the core–mantle boundary after planet accretion and solidification based on Stixrude (2014) and Noack & Lasbleis (2020). Both the iron content and the planet mass have a strong effect on the planetary structure and temperature profile. The resulting initial temperature profile in mantle and core, as well as the density, gravitational acceleration, pressure, and electrical conductivity profiles are shown in Fig. 1 for all investigated planet cases. The thermodynamic parameters (density as well as heat capacity and thermal expansivity of the material) were calculated following Noack et al. (2017) and using equations ofstate from Stixrude & Lithgow-Bertelloni (2011) for the mantle and Bouchet et al. (2013) for the metal core (here composed only of pure iron). The electrical conductivity profile is taken to be Earth-like from Xu et al. (2000). Thermal conductivity was calculated with the parametrisations derived in Tosi et al. (2013). For simplicity, the profiles for the different mantle properties were fixed during the evolution and did not take into account the influence of local melting or evolving chemical heterogeneity. The initial temperature profile, which is used in our interior evolution studies (described in Sect. 2.4), resembles the temperature profile depicted in Fig. 1. The temperature at the core–mantle boundary, however, and therefore the adiabatic temperature in the core, was scaled with the core–mantle boundary (CMB) temperature jump, which was varied in our study to investigate its effect on outgassing efficiency. Furthermore, we did not assume the existence of a primary atmosphere (see discussion in Sect. 4.1). During the evolution of Proxima Cen b, changes in atmospheric composition and mass can affect the surface temperature. However, since we did not couple our outgassing model to an atmospheric evolution model, surface temperatures remain constant throughout the evolution.

Fig. 1 Interior profiles for density, pressure, gravitational acceleration, post-magma-ocean temperature and electrical conductivity in the mantle for Proxima Cen b. Colours indicate planet masses between 1.27 (orange) and 2.54 (dark purple) Earth masses (related to different observation inclinations), and line styles indicate different iron mass fractions (from 20 to 60 wt%) for each of the investigated masses.

2.2 Evolution of the rotation and the magnetic field of Proxima Cen

On the main sequence, the activity level of a star is determined primarily by its mass and rotation rate. The general trend is that more rapidly rotating stars have stronger magnetic fields and emit more X-rays, up until a certain rotation rate called the saturation threshold, where activity no longer depends on rotation (Pizzolato et al. 2003; Vidotto et al. 2014; Johnstone et al. 2021). This dependence was recently confirmed for late M dwarfs (Wright et al. 2018). As stars age, contraction causes their rotation rates to increase on the pre-main sequence, and angular momentum removal by magnetised stellar winds causes their rotation rates to decrease once the stars reach the main sequence (Bouvier et al. 2014). At young ages, stars can reach rotation rates of ~ 100 Ω_⊙, where Ω_⊙ is the solar rotation rate.

Since the rotation rate of the star determines the magnetic field strength and the frequency of the magnetic field variation as seen by the planet and thus influences the induction heating, we needed to reconstruct the rotational evolution of Proxima Cen. Young stars exhibit a large spread in their rotation rates, and although they spin down due to angular momentum carried away by the stellar wind and finally converge into a single rotational track, they follow very different evolutionary tracks (Johnstone et al. 2015; Tu et al. 2015). Although it is possible that a subset of M dwarfs remain rapidly rotating for their entire lifetimes (Irwin et al. 2011), this has not been the case of Proxima Cen given its current slow rotation rate. Stellar rotation is also directly connected to stellar activity and a short wavelength radiation level (Tu et al. 2015; Johnstone et al. 2021), but in this study we are primarily interested in the magnetic field.

Stellar rotation influences induction heating in two ways. First, it influences the strength of the stellar magnetic field, and second, it determines the penetration depth of the variable magnetic field into the conducting planet’s mantle (see Kislyakova et al. 2017 for the discussion). In general, multiple pathways for the evolution of the stellar rotation are possible. They mostly depend on initial conditions, with some stars born as fast, slow, or intermediate rotators (Henderson & Stassun 2012; Johnstone et al. 2015, 2021). In the case of M dwarfs, the difference between slow and fast rotators in terms of XUV radiation is less significant than for more massive G stars (Johnstone et al. 2021). However, as we show below, the difference in the stellar rotation still matters for the magnitude of induction heating.

To estimate possible pathways of the rotational evolution of Proxima Cen, we used the model presented in Johnstone et al. (2015). The rotation tracks were fit to observational constraints on the rotational evolution of slow and fast rotating M dwarfs (Irwin et al. 2011). Figure 2 illustrates two possible evolutionary tracks of Proxima Cen, with the purple line showingthe fast rotator case, and the green line the slow rotator. Rotation tracks in between the fast and slow rotators are also possible, but for simplicity we focused on these two possibilities. As one can see, for such a low-mass star as Proxima Cen, the rotation tracks converge first at the age of a few Gyr. The black star marks the observed present-day rotation rate of the star. As one can see, Proxima Cen might have been a very fast rotator in the past, with a maximum rotation rate ~ 10 (slow rotator) and ~100 (fast rotator) times faster than the current Sun. Due to contraction, stellar rotation accelerates within the first ~ 150 Myr. After that, it slowly decelerates due to angular momentum loss from the stellar wind.

At present, Proxima Cen has a very long rotation period (83.5 days according to Kiraga & Stepien 2007 and 89.8 days according toKlein et al. 2021). It is highly likely that this star had a much shorter rotation period, and therefore a much stronger magnetic field in the past, which has important implications on the amplitude of induction heating. Previously, Reiners & Basri (2008) observed the average magnetic field of Proxima Cen of 600 G. Proxima Cen also exhibits a photometric cycle spanning approximately seven years (Yadav et al. 2016); therefore, one can expect the stellar magnetic field to also exhibit cyclic behaviour. Recently, Yadav et al. (2015) reproduced observed magnetic fields of M dwarfs using their dynamo model. They showed that a typical magnetic field of a young fast rotating M dwarf is about 2100 G, which is the field we adopted for the earlier stages of the star’s and the planet’s evolutions. Furthermore, Yadav et al. (2016) used Proxima Cen as an example to show that the stellar magnetic field starts declining once the rotational period of the star reaches the critical length of about 20 days. According to our rotation model shown in Fig. 2, Proxima Cen reaches this rotation period at the ages of 3100 Myr or 1800 Myr if it was born as a fast or slow rotator, respectively. To account for the evolution of the stellar magnetic field beyond this point, we used the relation presented by Vidotto et al. (2013), their Eq. (11), for the decay of the stellar magnetic field. In this model, the magnetic field decays as a power law. We have determined the steepness of this power law knowing the strength of the magnetic field we adopted before the magnetic field starts to decline (2100 G) and the observed strength of the magnetic field today (600 G). Recent observationsin radio frequencies also confirm the magnetic field strength of 600 G (Pérez-Torres et al. 2021). The resulting evolution of the magnetic field for the fast and slow rotator cases is presented in Fig. 3.

Recent observations using Zeeman-Doppler Imaging by Klein et al. (2021) indicate that the dipolar magnetic field strength of the star observed near the maximum of its magnetic cycle equals 135 G, and the dipole is tilted by 51° with respectto the star’s rotational axis. We considered a stronger field of 600 G at present following Reiners & Basri (2008) and Pérez-Torres et al. (2021), but various inclination angles of the magnetic field which results in similar magnetic fields at the planet’s orbit. We do not consider the influence of the stellar cycle on the magnetic field even after the rotation period reaches a value of 20 days. However, as we show below, induction heating mostly influences the interior evolution of the planet early on, when the stellar rotation is fast and the magnetic field of the star is strong; therefore, our results for the present-day state of the planet can be treated as an upper limit.

Fig. 2 Possible rotational evolution tracks of Proxima Cen (1 Ω⊙ = 2.67 × 10−6 rad s−1). The black star shows the modern observed rotation rate of Proxima Cen, which corresponds to a rotation period of ~ 82.5 days. The triangles indicate the observational constraints on the 10th and 90th percentiles of the rotational distributions from Irwin et al. (2011). The purple line corresponds to the evolution track of a fast rotator, the green line shows the evolution track of a slow rotator. All values in between are also possible.

Fig. 3 Possible evolution tracks of the magnetic field of Proxima Cen for the fast and slow rotators.

2.3 Induction heating of planetary interiors

We calculate the strength of induction heating in the interior based on the model by Kislyakova et al. (2017) and Kislyakova et al. (2018). In our work, we assumed that Proxima Cen b is a mostly dry, rocky planet, and we investigated induction heating for each of the planet masses and compositions given in Sect. 2.1. Moreover, we considered different inclinations between the stellar magnetic dipole and stellar rotation axis (we remark that this is not the same inclination as that of the planet’s orbit as seen from the Earth). The effect of the inclination of the stellar dipole on the strength of induction heating has been discussed in detail in Kislyakova et al. (2017). For Proxima Cen b, the measured inclination of the stellar magnetic dipole with respect to the stellar rotational axis is 51° (Klein et al. 2021). However, this value might change both in the short term due to a stellar cycle, and in the long term due to the evolution of the stellar magnetic field. For this reason, we compare the effect of different magnetic field inclinations ranging from 10° to 80°.

In our model, we used an approach for a sphere divided into shells with varying conductivity, presented by Parkinson (1983). One can find energy release in a spherical layer knowing the induction current flowing inside of it, which can be found based on the vector potential of the magnetic field. We assumed that the star has a dipole-dominated magnetic field, which becomes uniform at distances greater than the source surface due to super-radial expansion (Johnstone 2012). In this case,the induced magnetic field inside the planet can be described by a spherical harmonics of the first order, that is, the induced field is a dipole. If one adopts a spherical coordinate system with standard coordinates r, ϕ, and θ for radius, azimuthal, and polar angles, then the vector potential inside the planet only has a ϕ-component and is given by (1)

where the function F(r) depends only on r and is a surface spherical harmonics of the first order, where n and m are the spherical harmonic’s degree and order, respectively. Here and below, the equations are given in SI units. We assume that the angle θ is measured from the direction of the magnetic field outside the planet.

On the one hand, the current flowing within each layer can be found as ∇ × B = μ₀j. In our approximation of a skin effect, we assume that electromagnetic waves do not propagate into the medium, but can only penetrate up to a skin depth δ. We therefore do not consider wave propagation effects, and the electric field is not included in this equation (Laine et al. 2008). On the other hand, ∇×B = ∇²A = k²A, where A is given by Eq. (1), μ₀ is the magnetic constant, j is the current density, and B is the magnetic field. Therefore, we obtain (2)

as the onlycomponent of the current flowing within the lth layer. Here, ω is the frequency of the magnetic field variation, k² = −iωμσ is the wave number, μ is the magnetic permeability, σ_l is the electrical conductivityof rocks in the lth layer, and i is the imaginary unit. The energy release in the layer can be obtained as (3)

where we integrate over the volume V of the layer. H^I is then defined as the total energy released in the entire mantle. One can see that the energy release will also depend on the angle θ in the same way as the current does (see Eq. (2)). In Eq. (3), the factor of 1/2 originates from the averaging of the energy release over the period. For the detailed derivation of the equations, we invite the reader to consult Kislyakova et al. (2017). In our study, we neglected the influence of melt or variable iron content on the electrical conductivity profile, which is considered to be Earth-like (Xu et al. 2000). The potential influence of melt on the conductivity profile and hence induction heating is discussed in detail in Kislyakova et al. (2017).

In our model approach, we go beyond our previous studies (Kislyakova et al. 2017, 2018; Kislyakova & Noack 2020) by (a) modelling the stellar rotation and predicted magnetic field strength over time (described in Sect. 2.2), and (b) studying the local, latitude-dependent heating effect instead of assuming global, latitude-averaged values. Specifically, a planet will receive stronger induction heating at the equator than at its poles and show stronger volcanic activity. In this work, we investigated how strong the effect of latitude-dependent heating is on the local thermal-magmatic evolution (see Sect. 3.1 and Fig. 6). The frequency of the magnetic field variation ω and the strength of the external magnetic field are determined by the evolution of the star’s rotation and magnetic field described in Sect. 2.2.

As the star ages, its rotation decelerates, and at some point the planetary orbital period and the stellar rotation become synchronised. In our model, this leads to the energy release vanishing, because the frequency of the magnetic field variation is given as a residual of the frequency of the planet’s orbital motion and the stellar rotation frequency due to the planet’s prograde motion. After the synchronisation, the star continues to decelerate, and ω is then mostly determined by the planetary orbital motion. This resumes the induction energy release in planetary interiors, which, however, never reaches the highest values again (see Fig. 4 in the results section for energy release rates inside Proxima Cen b for different magnetic fields).

2.4 Interior evolution and outgassing model

Startingfrom the initial profiles derived in Sect. 2.1, we modelled the long-term evolution of Proxima Cen b for different possible observational inclinations, leading to planetary masses between 1.27 and 2.54 Earth masses. We assumed different iron weight fractions (20, 40, and 60 wt%; see Table 2), leading to planet radii between 1 and 1.35 Earth radii.

Several parameters that strongly influence the long-term evolution of Proxima Cen b are currently impossible to constrain, such as its exact surface temperature, temperatures in the interior, mantle rheology, volatile content and concentration of heat-producing, radioactive elements. We therefore defined a range of possible values for each of these unknown parameters and used a random selection routine to provide four main example parameter cases that we investigated in more detail for all planet masses and compositions. All varied parameters are listed in Table 3. To determine the mantle thermodynamic properties, the thermal and electrical conductivity, rheology parameters, as well as the melting temperatures, we assume as reference a rather dry mantle comparable to Earth’s bulk mantle as in Noack et al. (2017) and Kislyakova et al. (2017). Volatiles are only present as trace amounts (varied as listed in Table 3). They are incompatible elements in mantle minerals and preferentially partition into melt upon melting, leading to volcanic outgassing at the surface. The redox factor determines how strongly the oxidation state of the mantle deviates from the Iron–Wustite buffer, which affects the outgassed species entering the atmosphere (Ortenzi et al. 2020). Here, we consider the outgassing of H₂, H₂O, CO, and CO₂ depending on the volatile content in the melt and the mantle redox state. Carbon partitions into melt also depending on the redox state (Holloway et al. 1992; Grott et al. 2011; Ortenzi et al. 2020), whereas partitioning of water into the melt depends mainly on the melt fraction (Katz et al. 2003). The melt and outgassing models are described in detail in Noack et al. (2017) and Guimond et al. (2021).

We used a 2D convection code (Noack et al. 2016) to model the heat transfer and chemical mixing in the mantle assuming a compressible fluid under the truncated anelastic liquid approximation (TALA, benchmarked with King et al. 2010)in a half-spherical annulus (Hernlund & Tackley 2008). We modelled a viscous mantle without plastic deformation, meaning we do not consider that plate tectonics may occur on the planet. We use an Arrhenius viscosity law following Karato & Wu (1993) for upper mantle minerals and Tackley et al. (2013) for the lower mantle (using the original viscosity parameters and not applying an additional pre-factor of 100 in contrast with Tackley et al. 2013).

The radial resolution of the cells varies with depth, with the highest resolution of 10 km applied close to the surface of the planet and the coarsest resolution of 100 km applied at the bottom of the mantle. This non-uniform grid allows for a high-resolution treatment of the melt-producing region. Additionally, it is numerically inexpensive, since the average resolution of all cells is 50 km in lateral and radial directions.

We assumed that the convecting mantle is separated from the surface by an isolating lithosphere (here called a thermal boundary layer). The initial temperature increases linearly from surface temperature to the upper mantle temperature set at the bottom of the initial thermal boundary layer. Throughout the mantle, the temperature increases along an adiabat. The core is assumed to be super-heated after magma ocean solidification, and a steep temperature increase is applied at the bottom of the mantle over the thermal boundary layer at the CMB. Table 3 lists the initial temperature values and initial thermal boundary layer thickness. The core is allowed to cool over time depending on the heat flowing into the mantle. We note that no freezing of the core was considered here, which could lead to an additional heating of the mantle due to the release of latent heat upon crystallisation (Stevenson et al. 1983).

For Earth-like mantle rocks, melt is chemically buoyant only up to a specific pressure, which is called the density cross-over pressure and depends on planet composition. The assumed density cross-over pressure varied between 8.4 and 13.6 GPa in our parameter cases, which is in the range of values suggested for Earth and Mars (Ohtani et al. 1995). Melt can form deeper in the mantle but may not be able to reach the surface. We therefore only account for melting in the uppermost layer of the planet where melt can be buoyant, since we are interested in the outgassing efficiency of Proxima Cen b over time. Outgassing is also linked to the fraction of extrusive volcanism in comparison to intrusive magma, that does not reach the surface. We varied the proportion of extrusive melt according to the range of values observed for Earth (Crisp 1984).

Trace elements in Earth’s mantle are very important for its long-term evolution: these include the radioactive elements (Th, U, and K) and volatile elements (water, carbonates). The abundances of radioactive elements have not been determined yet for Proxima Cen (and may not directly correspond to abundances in the planet depending on the condensation temperature and planet formation history). In the absence of precise measurements, we assumed values between 0.5 and 1.5 times present-day Earth-like abundances (listed in Table 3) and calculated them back for 4.5 Gyr to obtain initial abundances of radioactive elements and an initial heating rate H^R (t = 0), which then decays over time in our simulations (Schubert et al. 2001; Breuer 2009). Heat sources are distributed homogeneously in mantle and crust. In reality, radioactive heat sources partition into the crust over time, therefore especially influencing the later stages of thermal evolution of the mantle (Laneuville et al. 2013). The mechanism of heat source redistribution is not considered to simplify the comparison of the effect of induction heating and planetary parameters on outgassing.

Similar to radiogenic heating, tidal heating is added as a mantle-averaged value, calculated for each planet individually.We determined the strength of tidal heating using a planet-averaged formula for tidal energy release (H^T, often denoted in the literature as Ė) derived for a homogeneous planet: (4)

where obliquity I is here set to zero. The term k₂∕Q can be derived using a complex formulation of the shear modulus (see e.g. Henning et al. 2009): (5)

where μ is the material rigidity and is written as a complex number representing energy storage and energy loss (Henning et al. 2009): (6)

for mantle-averaged shear modulus M, viscosity η, and angular frequency ω. The term k₂ is the Love number and is calculated as (7)

such that k₂ is also a complex number and (8)

Q is the quality factor and depends on the shear modulus, the viscosity, and the orbital frequency.

Heat produced by induction heating H^I (calculated for the initial mantle property profiles and according to time-dependent stellar variations as described in Kislyakova et al. (2017) and in Sect. 2.3) is added as an additional heat source in the thermal evolution models. We remark that there is no feedback from the mantle convection model to the a priori determined induction heating evolution.

3 Results

3.1 Evolution of the induction heating in the mantle of Proxima Cen b

In this section, we present our calculations of energy release by induction heating H^I in the planet’s mantle. Figure 4 shows the evolution of the induction heating in the mantle and illustrates the influence of different stellar and planetary parameters on the heating magnitude. Upper left panel of Fig. 4 compares the heating rates assuming the same mass and iron content of the planet (1.8 M_⊕, 20 wt%), but different rotation evolution tracks of the host star. For the fast and slow rotator cases, the magnetic field evolves as shown in Fig. 3. Heating due to electromagnetic induction is quite different in these two cases, with the maximum heating in the case of the fast rotator exceeding the heating in the slow rotator case by approximately a factor of seven. In this case, the time period when the heating is strong is longer, and the corotation point when the period of the stellar rotation equals the planet’s orbital period is reached later (2760 Myr and 1471 Myr for the fast and slow rotator cases, respectively). This result emphasises the importance of the stellar rotation for the planet’s evolution and shows that faster rotating stars not only produce higher levels of XUV (X-ray and extreme ultraviolet) radiation that can drive atmospheric escape, but also have stronger influence on the planet’s interiors through induction heating.

The upper right panel of Fig. 4 shows how the planet’s mass influences the heating. In this figure, we show the computed induction heating assuming different inclinations of the planet’s orbit as seen from Earth, which result in different estimates for the planet’s mass. In all cases, the star was assumed to be a fast rotator, and the iron content in the planet was equal to 20 wt%. This result shows that the heating is the strongest for the largest mass of the planet, if all other parameters are kept equal.

In the lower left panel of Fig. 4, we have investigated the influence of the iron fraction on total energy release by induction heating. Even though the iron content of the mantle is not varied in our study and therefore has no direct influence on the electrical conductivity (Fig. 1) or energy release by induction heating, an increase in bulk planet iron content leads to a larger iron core and decreased planet radius (see Table 2). The largest energy release is therefore obtained for the iron mass fraction of 20 wt% and hence the largest planet radius and surface area.

Finally, the lower right panel of Fig. 4 illustrates the influence of the inclination of the star’s dipolar axis with respect to the star’s rotation axis under the assumption that the planet orbits in the equatorial plane of thestar. Since a dipolar field that has a higher inclination produces a larger variation of the magnetic field at the planet’s orbit, energy release is consequently the strongest for the highest considered inclination of 80°. The inclination angle of the stellar dipole can change both on long timescales during the star’s evolution and on short timescales due to stellar cycles (e.g. Fares et al. 2017; Boro Saikia et al. 2018).

In all tracks presented in Fig. 4, one can see that the heating first increases and reaches the maximum, which coincides with the fastest rotation of the host star. Then, the heating decreases and even vanishes at the corotation point. The two local maxima on both sides of the corotation point arise due to an interplay between the magnetic field’s penetration depth into the mantle and the heating rate per volume. In general, the heating rate per volume is lower for lower frequencies of the magnetic field variation. However, as this frequency decreases, the penetration depth of the variable magnetic field in the conducting medium increases, which leads to energy release taking place in a larger part of the mantle and therefore to a second increase in energy release. However, at the corotation point, induction heating is negligible, leading to a sudden decrease before and a third increase of energy release after the corotation phase. This effect produces the two local maxima on both sides of the corotation point (at approximately 2600 Myr and 3100 Myr for the fast rotator and 1200 Myr and 1900 Myr for the slow rotator). In our model, at the corotation point energy release completely vanishes because the frequency of magnetic field variation becomes increasingly long. In all cases shown in Fig. 4, energy release is insignificant at ages later than 3 Gyr due to a relatively weak stellar magnetic field. As stated above, we assume that the star’s magnetic field starts to decrease after a certain time (approximately 3500 Myr for the fast and 1900 Myr for the slow rotator).

Figure 5 illustrates the distribution of the energy release rates inside Proxima Cen b’s mantle at different ages. Here, we assumed an identical interior composition, but we compared the interior heating rates for the slow and fast rotators. Energy release close to the surface strongly depends on the local conditions, specifically the rock density and electrical conductivity (see Kislyakova et al. 2017). Mineral phase changes lead to sudden jumps in both density (Stixrude & Lithgow-Bertelloni 2011) and electrical conductivity (Xu et al. 2000), and hence to a more step-wise function for the energy release rate. The main difference between the heating rates presented in the left and right panels are the heights of the maxima at the young ages (black and red lines). A higher heating rate is maintained longer if the star evolves as a fast rotator. It should be noted that here we assumed a solidified mantle to calculate the expected energy release rates, and that we do not model the magma ocean phase of Proxima Cen b. The second difference at the young age is the different depth at which the main heating takes place. In the fast rotator case, the frequency of the magnetic field variation is higher, so that the field doesn’t penetrate deeply into the planet’s mantle. This leads to most of the heat being released closer to the planet’s surface than in the case of slow rotator. The energy release rate is very low near the corotation point (2760 Myr for the fast rotator, 1471 Myr for the slow rotator). When the star’s rotation is exactly synchronised with the planet’s orbital motion, energy release completely vanishes. Here, we chose a time point close to it to illustrate that magnetic field penetrates very deeply into the planetary mantle, but produces a very low heating rate. Finally, at later ages both cases show very similar low heating rates. After several Gyr, both rotational tracks converge, which leads to very close rotation rates and magnetic field strengths. For this reason, energy release distributions at the age of 5 Gyr are practically identical in both cases.

In Fig. 6, we show exactly where induction heating reaches its maximum inside a rocky body (left panel) for a 1.8 M_⊕ planet (i.e. an inclination of 45°), an iron content of 20 wt%, and we applied parameter case 3 for the fast rotator scenario. Since the magnetic field of the star does not penetrate deeply into the interior, energy is released in the uppermost part of the mantle, where we would normally expect a lithosphere to form (for stagnant lid planets or underneath continental crust for plate tectonics planets). Induction heating is weaker at the poles than at the equator (see Sect. 2.3). As shown in the right panel of the figure, local temperatures can increase by several hundreds of K compared to a planet heated only by radioactive heating and secular cooling. Induction heating and temperature variations are shown after 1 Gyr of thermal evolution, i.e. after the highest induction heating occurred.

Close tothe equator, where the (latitude-dependent) induction heating is strongest, temperatures increase in this case by up to 400 K. Atthe poles, heating is weaker but still considerably exceeds the reference case without any induction heating. The temperaturevariation scale is normalised to 0 K. In some regions in the mantle, temperatures are locally colder for the simulation with induction heating compared to the simulation without considering this heating source, indicated by shades of blue. Such variations are expected in the mantle, which convects vigorously and shows different convective patterns for different simulations. The shades of blue and red in the convecting mantle also indicate that the additional heat created in the uppermost mantle for the induction heating case does not easily conduct downwards into the mantle.

Fig. 4 Energy release rates inside Proxima Cen b for different model setups compared to our reference case (1.8 M⊕, fast rotator, 20 wt% iron and inclination between the star’s magnetic field axis and the rotational axis of 60°). Top left: two different stellar rotation histories. Arrows indicate corotation time, when the stellar rotation period equals the planetary orbital period of 11.2 days. Top right: different observational inclinations of the Proxima Cen system, leading to different actual planet masses (M = Mmin/sin(i)). Bottom left: influence of planet iron content. Bottom right: influence of inclination between the star’s magnetic field axis and rotational axis.

Fig. 5 Induction heating inside Proxima Cen b averaged over latitude at different evolutionary stages. The heating has been calculated for a 1.8 M⊕ and 1.2 R⊕ planet with an iron content of 20 wt%. Left panel: distribution of induction heating inside Proxima Cen b if Proxima Cen has evolved as a fast rotator. Right panel: the same if Proxima Cen has evolved as a slow rotator. The given times are as follows: i) during the early pre-MS phase (1.8 Myr); ii) at the time of the maximum heating, which coincides with the fastest stellar rotation (219 Myr for the fast rotator, 190 Myr for slow rotator); iii) near the corotation time, when the stellar rotation period nearly equals the planetary orbital period of 11.2 days and the energy release is near its minimum (2760 Myr for the fast rotator, 1471 Myr for the slow rotator); iv) at an old age of 5004 Myr, when the energy release is determined by the planetary orbital motion. We note that the energy release at time step iv) is the same for the fast rotator and for the slow rotator. For higher rotation rates, the energy release declines faster inside the planetary interiors due to shallower skin depths.

Fig. 6 Local induction heating and resulting temperature variations compared to a simulation without induction heating after 1 Gyr of thermal evolution. We show an example for 1.8 M⊕, iron content of 20 wt% and parameter case 3 for a fast rotator. Here, we only show the silicate mantle from the core radius of 2922 km to the planet radius of 7840 km.

3.2 Comparison of tidal heating, induction heating, and radioactive heating

In this section, we investigate how the strength of induction heating compares to other heat sources such as radioactive heating in the mantle, secular cooling from the core, or possibly tidal heating due to a non-circular orbit or non-zero obliquity. It is not yet entirely clear whether Proxima Cen b does have an eccentric orbit (e.g. due to possible companion planets, should they be confirmed, Damasso et al. 2020).

If we assume an eccentricity of 0.17 (Damasso et al. 2020), we can determine how strong tidal heating would be in the interior of Proxima Cen b depending on the average mantle viscosity and shear modulus. Assuming that Proxima Cen b is rather dry in its interior, we can assume an average mantle viscosity of the order of 10²¹ Pa s (Karato & Wu 1993). For the shear modulus, we assumed a value of M = 100 GPa, similar to Earth’s rock values. The resulting values for the Love number k₂, the quality factor Q and the resulting tidal energy released in the interior of Proxima Cen b are listed in Table 4. We note that the produced tidal energy critically depends on the mantle rheology. A more reduced viscosity in the mantle would lead to stronger tidal heating than assumed here.

Comparing the total energy release in TW for tidal heating to radiogenic heating and induction heating (Table 4), one can see that (a) the total energy release for all three heating contributions increases with planet mass, (b) radiogenic heating is (at least when averaged over time) the strongest heat source, and (c) the total energy release is of the same order of magnitude for tidal heating and induction heating. It should be noted, however, that in the case of radiogenic and tidal heating, we assumed an equal distribution of energy release everywhere inside the planet, whereas induction heating actslocally and only close to the surface. Furthermore, radiogenic heat sources decrease over time and induction heating is strongly time dependent, whereas for tidal heating we used constant values over time. For induction heating, local values are therefore much higher than for radiogenic or tidal heating. They range from 61 to 141 pW kg⁻¹ at the point of fastest stellar rotation, leading to more partial melting of the mantle close to the surface. Radioactive heating, in comparison, contributes between 3.5 and 30 pW kg⁻¹ in our set ofsimulations, depending on the time (with higher values earlier on) and assumed concentration of radioactive heat sources (see Table 3). The effect of radioactive heating is therefore much stronger than tidal heating over the entire evolution (locally and globally), while induction heating is the dominant heat source in the upper mantle during fastest stellar rotation. Below we compare the contribution of the different heating sources to volcanic outgassing.

It should be noted that the purpose of this study is not to carry out a detailed investigation of tidal heating or interactions between different heating sources. We therefore chose to compare the different energy sources in the mantle (radioactive heating with secular cooling of the core; tidal heating; induction heating) only for one simple reference case, where we assumed an inclination of observation of 45° (i.e. 1.8 Earth masses) and an iron content of 20 wt%, parameter case 3 and a magnetic field inclination of 60°.

Figure 7 shows how strong depletion of the mantle by partial melting (i.e. dehydration and efficiency of crust formation) depends on mantle heating sources for this example case. Blue colours indicate primordial mantle material. We allowed for a maximal depletion of the mantle of 30%, which would lead to a depleted harzburgite layer as suggested for the uppermost mantle of Earth due to the extraction of basaltic material forming a crust (which is not specifically modelled here). The four panels show the state of mantle depletion by melting at 2 Gyr for a case where we only assumed radioactive heating and secular cooling of the core (top), added tidal heating, added induction heating, and a combination of all heating sources. The strength of depletion increases from the first to the last panel. The corresponding total outgassed pressure within the first 2 Gyr increases from 10.8 bar (top panel) to 27.3 bar (bottom panel), since the two added heat sources (tidal heating and induction heating) are most important close to the surface, leading to more melting in the shallow mantle. In the investigated case, induction heating contributed to an additional outgassing of 8.3 bar, while tidal heating led to additional outgassing of 7 bar in the same time frame, leading to more remnant primordial material in the upper mantle than with the added influence of induction heating.

In general, we see that radioactive heating and secular cooling dominate the global energy release, whereas tidal heating and induction heating lead to more melting close to the surface (see Fig. 7). In addition, since the induction heating varies with latitude, we can see an increase of mantle depletion (and hence a stronger influence of induction heating) in the centre of our half-sphere, which relates to the equator of the planet.

Another influence on the depletion of the mantle comes from mineral physics; as can be seen in the melt depletion plots in Fig. 7, depleted material is transported away from the melt region, which lies close to the surface at the bottom of the lithosphere, into the mantle. However, we can see that at least initially a two-layer convection pattern evolves due to the phase transition of perovskite to post-perovskite when assuming that diffusion creep dominates the lower mantle (Tackley et al. 2013). The depleted material resides in the upper half of the mantle during the first two billion years, even though we did not account for additional buoyancy of depleted material (which may be slightly lighter than primordial surrounding rocks). During the later evolution, however, the depleted material is mixed into the lower mantle as well, allowing for an almost complete depletion and dehydration of the mantle (see Fig. 10).

3.3 Volcanic outgassing scenarios

In all simulations discussed in this section, we considered radioactive heating homogeneously in the mantle, as well as cooling from an initially super-heated core (see Fig. 1). Tidal heating was not considered here. We describe the outgassing evolution scenarios obtained for models with and without induction heating as an additional energy source in the (uppermost) mantle. We began our simulations from the solid-state stage and did not take into account any influences of our parameter cases and heating mechanisms on the magma ocean stage (see Sect. 4.1). Instead, the initial state of the mantle – in terms of both temperature and composition – is varied following the four parameter cases listed in Table 3.

In Figs. 8–10, we compare different volatile outgassing scenarios depending on time. Figure 8 compares the total accumulated amount of volatiles outgassed from the mantle for parameter case 4 for a planet mass of 1.8 M_⊕ and a low iron content of 20 wt%. The total outgassed pressure is calculated from the mass of volatiles (considering here H₂, H₂O, CO, and CO₂, following Guimond et al. 2021), assuming that no volatiles condense out of the atmosphere, react with the surface, or are lost to space. No influence of the atmosphere on the surface temperature is taken into account, which instead is set to the fixed value listed in Table 3. Such outgassing simulations therefore only serve as upper estimates of the actual atmospheric pressure butcan be considered as input for atmospheric loss studies.

The ratio of outgassed species depends on the assumed redox state of the mantle, which is varied here with respect to the IW buffer between reduced mantles (IW+ 0.6, case 2) to oxidised mantles (IW+ 3.52, case 3). In parameter case 1, the molar fractions of H₂, H₂O, CO, and CO₂ are 0.062/0.755/0.053/0.13. For case 2, they are 0.26/0.64/0.07/0.03, for case 3 they are 0.01/0.53/0.03/0.43, and for case 4 they are 0.03/0.52/0.11/0.34. In the most reduced case (parameter case 3), partitioning of carbonates into the melt is inefficient, leading mostly to the outgassing of hydrogen species (see also Ortenzi et al. 2020). The mantle redox state therefore has a first-order impact on the composition of the outgassed volatiles and influences the atmospheric pressure (since CO₂, for example, is much heavier than H₂). However, the redox state does not influence the amount of melt that is transported to the surface. This is instead varied via the density cross-over pressure, from which we assume melt to be denser than solid rock. The lowest density cross-over pressure that we considered was 8.4 GPa (case 4) compared to the largest value of 13.6 GPa (case 3), leading to more melt being transported towards the surface in the latter case. However, we also considered different ratios of extrusive to intrusive melt, where the largest amount of melt that actually reaches the surface is in case 4 with a percentage of extrusive volcanism of 27%. Both the density cross-over pressure and the amount of extrusive volcanism do not seem to have a first-order effect on the total outgassed atmospheric pressure (see Fig. 11). The strongest outgassing is observed in most simulations for the first parameter case, which can be explained by the large temperature jump at the core–mantle boundary, increased radiogenic heat sources, as well as a low mantle viscosity. All three factors lead to efficient convection in the mantle and hence more volcanic outgassing (Dorn et al. 2018).

In parameter case 4, the mantle is initially moderately warm and heats up further during the first few hundreds of Myr. The applied viscosity pre-factor of 1.32 (see Table 3), comparable to a rather dry mantle, furthermore leads to slower heat transport through the mantle than would be the case for any of the other parameter cases. In cases 1 and 3, the viscosity is much lower, leading to faster convection, and in case 2 the mantle is initially much hotter, also leading to faster convection. Heating only initiates here after about 1–1.5 Gyr of thermal evolution. In the left panel of Fig. 8, we therefore compare the outgassing pressure starting at 1 Gyr and until a 50-bar-thick atmosphere is reached (again not considering any sinks or losses; e.g. escape to space for the atmosphere for simplification). The fast rotator scenario leads to strong local heating in the uppermost mantle (see Fig. 6), and first volcanic outgassing about 550 Myr earlier (at 0.87 Gyr) than in either the slow rotator scenario (1.42 Gyr) or the reference case without any induction heating (1.43 Gyr). The effect of induction heating for the slow rotator case is less strong than for the fast rotator case. To better visualise the outgassing evolution over time, in the right panel of Fig. 8 we show the difference between total outgassed pressure calculated with induction heating (for both the fast and the slow rotator cases) compared to the simulation without considering induction heating over the entire evolution period of 5 Gyr. Fluctuations in the pressure difference over time can be attributed to the irregularly occurring volcanic events, which add more volatiles to the atmosphere at different times for the different simulations. Outgassing increases significantly when accounting for induction heating. Between about 2 Gyr and the end of the simulation, the fast rotator case shows an almost constant increase in outgassing of 15 bar compared to the case without induction heating. For the slow rotator case, the difference in outgassing also reaches about 5 bar at the end of the simulation. In the following, we concentrate only on the fast rotator case.

Figure 9 compares different magnetic field inclinations for the fast rotator case, for the same scenario as in Fig. 8. Volcanic outgassing starts between 0.6 Gyr and 1.41 Gyr depending on the inclination of the stellar magnetic dipole. This is directly related to the induction heating strength depicted in the bottom right panel in Fig. 4, showing that the planet-wide heating by induction changes from 8.67 TW or 0.93 pW kg⁻¹ for a magnetic field inclination of 80° to 0.27 TW for 10° at the time of the highest energy release at around 220 Myr. The difference in outgassed pressure compared to the simulation without induction heating reaches 20–25 bar for the largest magnetic field inclination and less than 5 bar for the smallest inclination of 10°. We subsequently limited our investigation to a magnetic field inclination of 60°, which is close to the measured present-day inclination (Klein et al. 2021).

To better understand the strength of the influence of induction heating on outgassing at different times, in Fig. 10 we showthe efficiency of water outgassing over time together with the evolving pressure. The efficiency factor is the percentage of volatiles that can be directly outgassed from the interior. It takes into account that not all melt produced in the interior actually reaches the surface and hence contributes to outgassing (depending on the fraction of extrusive volcanism, Table 3). A value of 100% would therefore indicate maximal possible depletion and dehydration of the mantle since here no replenishment of volatiles is considered. This is due to the fact that we only investigated the scenario of a stagnant-lid planet, without subduction of water or carbonates back into the mantle and without delamination of the crust.

This comparison shows that the outgassing efficiency is mostly defined by other parameters such as the mass of the planet (hence observational inclination). However, induction heating typically leads to (a) faster outgassing, (b) stronger outgassing at earlier times, and (c) typically more outgassing over the entire thermal evolution.

In general, the effect of induction heating on outgassing is strongest for the cases where melting is in general more restricted (e.g. for higher iron content leading to higher pressures in the mantle), and where additional heat can increase the amount of melt in the mantle. In some cases, outgassing is increased by induction heating by a factor of ~ 10. For cases where outgassing is already moderately efficient in the absence of an additional heating source, induction heating leads to a weaker increase, with up to 30% additional melt. In some cases (e.g. for the lowest iron content and smallest investigated planet masses), melting is already very efficient without induction heating (as also shown in Noack et al. 2017 and Dorn et al. 2018), and a further increase in mantle temperature (due to induction heating or other heating processes) has a less prominent effect. On the other hand, for a few cases (especially for more massive planets with high iron content), even with induction heating almost no outgassing takes place. In these cases, additional heating sources such as tidal heating (Sect. 4) might become very important.

Furthermore, we see from the comparison of the outgassing efficiency and the total outgassed pressure that even though outgassing efficiency increases with decreasing planet mass (increasing observational inclination of the system), the total outgassed pressure is not necessarily the highest compared to the other planetmass cases. This is because the total pressure is calculated as the mass of the atmosphere M_atm times the surface gravitational acceleration g divided by the surface area of the planet: (9)

where M_atm, g, and R_p are given in SI units. The influence of planet radius and mass (included in the gravitational acceleration) on atmospheric pressure is therefore non-linear, leading to a smaller total outgassed pressure for the smallest planet mass (observational inclination of 90°) compared to the next largest mass (observational inclination of 60°), even though outgassing efficiency is higher.

We find that for the planet with 1.27 Earth masses, water outgassing efficiency (i.e. the total mass of water outgassed divided by the initial mass of water in the mantle) is in general high. The highest observed values reach 96% of water outgassing. Averaged over all parameter cases, 73% of all volatiles initially stored in the mantle are released to the atmosphere. For the planet with 1.47 Earth masses, outgassing efficiency decreases to about 61% in our investigated parameter space. For 1.8 Earth masses, in some cases no outgassing occurs at all, and the average outgassing efficiency is at 50%. For the most massive investigated planet with 2.55 Earth masses, in the investigated parameter space we see no outgassing in several cases, with on average only 25% of the volatiles released to the atmosphere.

A similar trend can be observed for the influence of iron. For the smallest iron content considered (20 wt%), average water outgassing efficiency is at 69%, decreasing to 57% and 31% when increasing the iron content to 40 and 60 wt%, respectively. This shows that the planet mass and the iron content remain the most important factors for volcanic outgassing, similarly to what has been found in previous mantle convection studies investigating the few factors that can influence volcanic outgassing (Noack et al. 2017; Dorn et al. 2018; Ortenzi et al. 2020).

In Fig. 11, we summarise the outgassing scenarios for all four parameter cases (Table 3). In the first three panels on the left, we plot the total accumulated outgassed atmospheric pressure for 20, 40, and 60 wt% iron (with considering induction heating). The three panels on the right show, for the same respective iron fractions, the difference in outgassing between simulations considering induction heating and cases without induction heating. Depending on the parameter case and planet mass, outgassing starts shortly after the beginning of the simulation or after about two billion years, if at all. During the early evolution, due to the additional heat provided by induction heating, induction heating produces stronger outgassing than for simulations where induction heating is not considered – leading to differences in the outgassing over time of up to 25 bar (and even 32 bar in one case). Induction heating is most effective in triggering increased volcanic activity and outgassing for the low-mass planets (1.27 and 1.47 M_⊕) for cases 2 and 4. For cases with efficient outgassing without additional heat sources (see left panels in Fig. 10), simulations without induction efficiently heat outgas volatiles from the interior at a later stage of the evolution, leading to smaller differences in the final total outgassed pressure after 5 Gyr than during the earlier stages of evolution. Induction heating therefore mostly drives early outgassing.

Fig. 7 Effect of radioactive heating (HR), tidal heating(HT), and induction heating (HI) on mantle depletion (and hence volcanic outgassing) for the example case from Fig. 6 after 2 Gyr. Blue colours show volatile-rich material and yellow to pink colours indicate depleted, volatile-deprived material. The half-sphere extends from the pole over the equator to the other pole. From top to bottom: only radioactive heating and secular cooling of the core (corresponding to 10.8 bar volcanic outgassed pressure), effect of adding tidal heating (17.7 bar), effect of adding induction heating (19.1 bar), and effect of adding both tidal and induction heating (27.3 bar).

Fig. 8 Influence of fast rotator (purple line) or slow rotator (green line) stellar evolution in comparison to a simulation without considering induction heating (dotted black line) on the outgassing evolution for a planet mass of 1.8 M⊕, a magnetic field inclination of 60°, 20 wt% iron content, and parameter case 4. Left: total outgassed pressure over time (limited between 1 and 3 Gyr for a better comparison). Right: difference in outgassing over time compared to a simulation where no induction heating was considered.

Fig. 9 Influence of the magnetic field inclination and hence induction heating strength. Shown is the case with an observation inclination of 45° (1.8 M⊕), 20 wt% iron content, and parameter case 4. The dotted line indicates the simulation without induction heating and almost coincides with the weakest induction heating case (magnetic field inclination of 10°). Left: total outgassed pressure over time (limited between 1 and 3 Gyr for a better comparison). Right: difference in outgassing over time compared to a simulation where no induction heating was considered.

Fig. 10 Influence of planet mass (i.e. inclination of observation) on outgassing efficiency (left, as the example shown for H2O) and resulting accumulated outgassed pressure (right) from the iron-poor (top) to the iron-rich scenario (bottom). Dotted lines indicate the corresponding simulations without induction heating.

Fig. 11 Accumulated outgassed pressure (left) when considering induction heating as an additional heat source, and difference inoutgassing to the same simulations without induction heating (right). The iron content decreases from top to bottom from 60 wt% to 20 wt%. Colours indicate different planet masses, whereas the line style changes by parameter case.

4 Discussion

4.1 Early magma ocean stage and atmosphere

All simulations presented in this study start from an already solidified mantle, that is, after the magma ocean phase. For a more realistic scenario, the early stage of planet accretion and the magma ocean (potentially multiple times) would need to be coupled with a sophisticated atmospheric model including climate calculations and atmospheric escape. In the absence of such a model, we can only speculate how the solidification time of the magma ocean would change when considering the different internal heating sources discussed here. The solubility of volatiles in the magma ocean at the surface of the planet depends largely on the atmospheric pressure and composition. Atmospheric losses to space lead to the cooling of the surface temperature, and hence the cooling and solidification of the mantle. Solidification of the magma ocean leads to an enrichment of volatiles in the residual magma, and volatile concentrations rise above the solubility threshold resulting in a catastrophic release of volatiles from the interior to the surface (Elkins-Tanton 2011). Prolonging the length of the fully molten magma ocean stage (e.g. due to strong additional internal heating sources) should not therefore, in principle, affect the volatile concentrations remaining in the solidified mantle (though this also depends on the evolution of the atmospheric pressure). The speed of magma ocean cooling, however, is more important as a fast-solidifying magma ocean may trap volatile-rich melts in the solid mantle (Hier-Majumder & Hirschmann 2017). However, the strength of both tidal heating and induction heating change dramatically from a solid-state mantle to a magma ocean. It remains to be investigated whether either of these heat sources would have a larger impact on the length of the magma ocean phase, as both tidal friction and induction heating (due to a strongly increased electrical conductivity in the magma) would be strongly reduced in the interior (Beuthe 2013; Kislyakova et al. 2017). In addition, the depth of the mantle (scaled here via the planet iron content), the magma ocean rheology, as well as the composition of the atmosphere, would have a first-order influence on the surface temperature and hence the length of the magma ocean phase (Nikolaou et al. 2019; Lichtenberg et al. 2021). We used variable initial volatile concentrations in our four parameter cases but kept the volatile budget in the mantle after solidification identical for our models with and without induction heating and/or tidal heating.

We point out that since we only consider a solid-state mantle here, the energy release calculated for induction heating at 1.8 Myr in Fig. 5 should only be taken as an example calculation indicating the change in induction heating strength due to the stellar rotational evolution. In reality, at this time the planet was likely still in a magma ocean phase (and may well have remained so for a few tens of Myr, depending on atmosphere escape efficiency), which would lead to quite different energy release profiles, as discussed above.

One additional assumption remains to be discussed, and this is the assumption that Proxima Cen b is in a stagnant-lid phase. There is currently no clear observational method available that could give us information on the tectonic state of any exoplanet – though some indirect measurements have been suggested, such as white dwarf pollution in continental crust components (Jura et al. 2014), or the observation of strong volcanic activity for cool, massive super-Earths, which would be difficult to explain in the absence of plate tectonics (Noack et al. 2017). However, looking at the rocky planets and moons in the Solar System, Earth seems to be the exception with ongoing crustal recycling via subduction. It is therefore not unreasonable to assume that plate tectonics on rocky planets is the exception rather than the rule. However, if Proxima Cen b had initiated plate tectonics in its past, this would have severe consequences for both the mantle volatile concentrations (due to replenishment from the surface) as well as the mantle temperature (due to the efficient cooling of a subducting cold crust and lithosphere into the mantle). The outgassing simulations shown here are therefore not comparable to a planet exhibiting plate tectonics (Noack et al. 2014), which remains subject of further studies, then including also the return flux of volatiles back into the mantle (Höning et al. 2014; Nakagawa & Spiegelman 2017).

4.2 Strength of induction heating compared to other heat sources

Internal heating mechanisms are very important for planetary evolution, and a thorough comparison of different heating mechanisms (specifically global vs. local heat sources) is needed. In this study, we compared Earth-like radioactive heating of the mantle, secular cooling of the core depending on various initial core temperatures (see Table 3), tidal heating, and induction heating. For the induction heating of planetary interiors, stellar rotational evolution is very important. In this paper, we consider the magnetic field to be saturated for stellar rotation periods shorter than 20 days, and that they decline after that, so the magnetic field drops from an initial 2100 G to 600 G at present (following Reiners & Basri 2008). This is higher than the most recent observation of the dipolar component of the magnetic field of Proxima Cen (135 G; Klein et al. 2021). On the other hand, we also considered different inclination angles for the stellar dipole around the observed one (51° with respectto the star’s rotation axis). We also compared two different stellar evolution tracks: a fast rotator and a slow rotator scenario. Asexpected, the interiors of planets orbiting rapidly rotating stars experience significantly stronger heating (seeFig. 2). It is, however, unclear, which scenario best describes the actual evolution of Proxima Cen.

In general, we find the influence of induction heating to be less strong at present than suggested for other planets (Kislyakova et al. 2017; Kislyakova & Noack 2020). The strongest heating effect is observed in the early evolution of the star, accumulating over the first few hundreds of Myr up to about 1 Gyr, and leading to local strongly enhanced temperatures by several hundred K (Fig. 6). At present, for both the fast rotator and the slow rotator scenarios, and for all investigated parameter cases, planet masses, and magnetic field inclinations, energy release inside Proxima Cen b should be negligible compared to other internal heating sources. We also note that we didn’t take into account stellar cycles that could further decrease the heating. Therefore, we conclude that induction heating has the strongest effect on the planet early in its evolution. However, for several simulations, the final amount of outgassed volatiles is much higher when induction heating was considered compared to simulations where only radioactive heating and secular cooling was included. Also, due the initially intense heating phase by induction heating, melting scenarios are in general shifted by several hundred Myr to earlier (and stronger) outgassing events in the case where induction heating is considered, even cases where the final amount of outgassed volatiles after 5 Gyr did not vary strongly according to the internal heat source.

Recently, Kislyakova & Noack (2020) showed that induction heating can be substantial even for planets orbiting solar-like stars with weak global magnetic fields of the order of a few Gauss. Moreover, they showed that induction heating might be one of the dominant heating mechanisms that can trigger outgassing on massive planets where volcanism is otherwise suppressed due to high pressures in the planets’ mantles. However, they considered a hot super-Earth HD 3167b orbiting very close to its host star. Here, we show that for planets in the habitable zone such as Proxima Cen b, induction heating can be relevant mostly early on in its history. One can safely assume that the effect would be stronger for planets orbiting at the inner edge of the habitable zone.

5 Conclusion

In this work, we studied the influence of induction heating and randomly selected initial planetary conditions on the evolution of volcanic outgassing of planets orbiting in the habitable zone of an M dwarf with the example of Proxima Cen b. We show that even though induction heating inside this planet is now moderate (mostly due to a relatively large orbital distance of this planet to its host star and the present-day stellar magnetic field strength), it was much stronger during early evolution, leading at that time to a local mantle temperature increase of several hundreds of Kelvin. Our models imply that induction heating is as significant as tidal heating and radioactive heating for the long-term evolution of Proxima Cen b, especially the depletion of volatiles by melting from the interior. In general, we observe that melting in the mantle and therefore volcanic outgassing at the surface occur earlier when accounting for induction heating (or tidal heating, if Proxima Cen b indeed had an eccentric orbit). In addition, for some of our investigated parameter cases we observed rather inefficient volcanic outgassing when not considering induction or tidal heating. For these cases, the outgassed pressure could change dramatically (by more than one order of magnitude) when considering the expected strength of induction heating on the upper mantle of Proxima Cen b.

In addition to the influence of different heat sources on outgassing efficiency, we observe that if the planet had a higher iron content than suggested by the Mg–Fe ratio observed in the stellar spectrum (Güdel et al. 2004), we would expect less outgassing than for an iron-poor and magnesium-rich composition. Also, for an increasing planet mass (correlating here to a decreasing inclination of the star–planet system with respect to our observational plane), our simulations would suggest a less efficient outgassing and hence lower accumulated outgassed pressure of a secondary atmosphere when not considering any atmosphere loss processes, which we did not study in this paper.

Acknowledgements

We would like to thank the anonymous reviewer for their suggestions that helped to improve the manuscript. LN acknowledges support by the German Research Foundation (DFG) for project no. 1324/6-1. The authors would like to thank the HPC Service of ZEDAT, Freie Universität Berlin, for computing time. KK and MG acknowledge the support by the Austrian Research Promotion Agency (FFG) project 873671 “SmileEarth”. This study was further supported by the Austrian Science Fund (FWF) project S11601-N16 “Pathways to Habitability: From Disk to Active Stars, Planets and Life” and the related subprojects S11604-N16 and S11607-N16. The authors thank the Erwin Schrödinger Institute (ESI) of the University of Vienna for hosting the meetings of the Thematic Program “Astrophysical Origins: Pathways from Star Formation to Habitable Planets” and Europlanet for providing additional support for this program. The authors also acknowledge the International Space Science Institute for the support of ISSI team 370.

References

All Tables

All Figures

Fig. 1 Interior profiles for density, pressure, gravitational acceleration, post-magma-ocean temperature and electrical conductivity in the mantle for Proxima Cen b. Colours indicate planet masses between 1.27 (orange) and 2.54 (dark purple) Earth masses (related to different observation inclinations), and line styles indicate different iron mass fractions (from 20 to 60 wt%) for each of the investigated masses.

In the text

Fig. 2 Possible rotational evolution tracks of Proxima Cen (1 Ω⊙ = 2.67 × 10−6 rad s−1). The black star shows the modern observed rotation rate of Proxima Cen, which corresponds to a rotation period of ~ 82.5 days. The triangles indicate the observational constraints on the 10th and 90th percentiles of the rotational distributions from Irwin et al. (2011). The purple line corresponds to the evolution track of a fast rotator, the green line shows the evolution track of a slow rotator. All values in between are also possible.

In the text

	Fig. 3 Possible evolution tracks of the magnetic field of Proxima Cen for the fast and slow rotators.
In the text

Fig. 4 Energy release rates inside Proxima Cen b for different model setups compared to our reference case (1.8 M⊕, fast rotator, 20 wt% iron and inclination between the star’s magnetic field axis and the rotational axis of 60°). Top left: two different stellar rotation histories. Arrows indicate corotation time, when the stellar rotation period equals the planetary orbital period of 11.2 days. Top right: different observational inclinations of the Proxima Cen system, leading to different actual planet masses (M = Mmin/sin(i)). Bottom left: influence of planet iron content. Bottom right: influence of inclination between the star’s magnetic field axis and rotational axis.

In the text

Fig. 5 Induction heating inside Proxima Cen b averaged over latitude at different evolutionary stages. The heating has been calculated for a 1.8 M⊕ and 1.2 R⊕ planet with an iron content of 20 wt%. Left panel: distribution of induction heating inside Proxima Cen b if Proxima Cen has evolved as a fast rotator. Right panel: the same if Proxima Cen has evolved as a slow rotator. The given times are as follows: i) during the early pre-MS phase (1.8 Myr); ii) at the time of the maximum heating, which coincides with the fastest stellar rotation (219 Myr for the fast rotator, 190 Myr for slow rotator); iii) near the corotation time, when the stellar rotation period nearly equals the planetary orbital period of 11.2 days and the energy release is near its minimum (2760 Myr for the fast rotator, 1471 Myr for the slow rotator); iv) at an old age of 5004 Myr, when the energy release is determined by the planetary orbital motion. We note that the energy release at time step iv) is the same for the fast rotator and for the slow rotator. For higher rotation rates, the energy release declines faster inside the planetary interiors due to shallower skin depths.

In the text

	Fig. 6 Local induction heating and resulting temperature variations compared to a simulation without induction heating after 1 Gyr of thermal evolution. We show an example for 1.8 M⊕, iron content of 20 wt% and parameter case 3 for a fast rotator. Here, we only show the silicate mantle from the core radius of 2922 km to the planet radius of 7840 km.
In the text

Fig. 7 Effect of radioactive heating (HR), tidal heating(HT), and induction heating (HI) on mantle depletion (and hence volcanic outgassing) for the example case from Fig. 6 after 2 Gyr. Blue colours show volatile-rich material and yellow to pink colours indicate depleted, volatile-deprived material. The half-sphere extends from the pole over the equator to the other pole. From top to bottom: only radioactive heating and secular cooling of the core (corresponding to 10.8 bar volcanic outgassed pressure), effect of adding tidal heating (17.7 bar), effect of adding induction heating (19.1 bar), and effect of adding both tidal and induction heating (27.3 bar).

In the text

Fig. 8 Influence of fast rotator (purple line) or slow rotator (green line) stellar evolution in comparison to a simulation without considering induction heating (dotted black line) on the outgassing evolution for a planet mass of 1.8 M⊕, a magnetic field inclination of 60°, 20 wt% iron content, and parameter case 4. Left: total outgassed pressure over time (limited between 1 and 3 Gyr for a better comparison). Right: difference in outgassing over time compared to a simulation where no induction heating was considered.

In the text

Fig. 9 Influence of the magnetic field inclination and hence induction heating strength. Shown is the case with an observation inclination of 45° (1.8 M⊕), 20 wt% iron content, and parameter case 4. The dotted line indicates the simulation without induction heating and almost coincides with the weakest induction heating case (magnetic field inclination of 10°). Left: total outgassed pressure over time (limited between 1 and 3 Gyr for a better comparison). Right: difference in outgassing over time compared to a simulation where no induction heating was considered.

In the text

	Fig. 10 Influence of planet mass (i.e. inclination of observation) on outgassing efficiency (left, as the example shown for H2O) and resulting accumulated outgassed pressure (right) from the iron-poor (top) to the iron-rich scenario (bottom). Dotted lines indicate the corresponding simulations without induction heating.
In the text

	Fig. 11 Accumulated outgassed pressure (left) when considering induction heating as an additional heat source, and difference inoutgassing to the same simulations without induction heating (right). The iron content decreases from top to bottom from 60 wt% to 20 wt%. Colours indicate different planet masses, whereas the line style changes by parameter case.
In the text