EFFECTS OF STELLAR FLUX ON TIDALLY LOCKED TERRESTRIAL PLANETS: DEGREE-1 MANTLE CONVECTION AND LOCAL MAGMA PONDS

Throughout the last decade there has been escalating work to model the internal structure, composition, evolution, and habitability of super-Earths (e.g., Kuchner 2003; Léger et al. 2004; Ehrenreich et al. 2006; Valencia et al. 2006, 2007a, 2007b). The extremely small semimajor axes of many exoplanet orbits suggest several new processes in planetary evolution that do not appear to have occurred in our solar system. Several studies, for example, characterize the effects of tidal heating on rocky planets (Jackson et al. 2008a, 2008b, 2008c; Barnes et al. 2008) and work to constrain surface temperatures from observables such as star size, planet size, semimajor axis, eccentricity, and obliquity. Recent geodynamical studies have studied the possibility and detectability of volcanism on rocky extrasolar planets (e.g., Kite et al. 2009; Kaltenegger et al. 2010). Others have addressed problems specific to tidally locked planets, such as magnetosphere (Grieβmeier et al. 2005) and climate (Joshi 2003) dynamics. The aim of this work is to further add to this library with particular attention to the effect of a tidally locked orbit on the mantle convection structure and long-term geodynamical evolution of terrestrial exoplanets.

Synchronous rotation refers to an orbit in which one rotation of the body takes precisely the same amount of time as one revolution about the orbital center. This also means that, in the case of a planet, one day is equal to one year, and one hemisphere will always be star-facing, and the other will never be. Tidal-lock occurs when a star's rotation and mass raise tides on a planet, decaying its orbit, and despinning it.

One of the main hypotheses tested here is whether low-degree mantle convection may occur due to fixed stellar heating on a synchronously rotating planet. Spherical harmonic degree-1 convection refers to a single pair of convection cells in a spherical body, and implies the presence of a single upwelling and a single downwelling in opposite poles. Previous work has suggested and explored low-degree convection, and is briefly reviewed in the discussion.

In the following section we formulate a numerical model to test tidal-lock boundary conditions on a convecting mantle. This simple model does not treat internal frictional tidal heating or atmospheric heat transport or insulation. We analyze mantle convection patterns through deep geologic time and calculate regions of possible stable liquid water near the planetary surface. From this we discuss the possibility of extreme temperatures that would permit localized magma oceans around the substellar point. Finally, we explore implications for the habitability of tidally locked Earth-like planets.

Numerical models have been formulated using SSAXC, a spherical axisymmetric version of the two-dimensional finite element code ConMan (King et al. 1990). ConMan models incompressible, infinite-Prandtl number thermochemical convection by solving the equations for conservation of mass,

$$\begin{equation} \nabla \cdot {\bf u} = 0, \end{equation} \tag{ 1 }$$

momentum,

$$\begin{equation} \nabla ^2 {\bf u} = \nabla p + Ra \, T \, {\bf \hat{z}}, \end{equation} \tag{ 2 }$$

and energy,

$$\begin{equation} \frac{\partial T}{\partial t} = - {\bf u} \cdot \nabla T + \nabla ^2 T, \end{equation} \tag{ 3 }$$

where u is the dimensionless velocity, T is the dimensionless temperature, p is the dimensionless pressure, ${\bf \hat{z}}$ is the unit vector in the vertical direction, and t is the dimensionless time. The Prandtl is the ratio of momentum diffusivity to thermal diffusivity: Pr = $\frac{\nu }{\alpha }$, where ν is kinematic viscosity and α is thermal diffusivity. In mantle silicates, viscosity is much greater than thermal diffusivity, and so SSAXC is appropriate for this application. All material properties used in these convection calculations are combined into the thermal Rayleigh number, Ra, which indicates the presence and vigor of convection, given by

$$\begin{equation} Ra = \frac{\rho g \alpha \Delta T d^3}{\kappa \eta _0}, \end{equation} \tag{ 4 }$$

where g is the acceleration due to gravity, α is the coefficient of thermal expansion, ΔT is the temperature drop across the mantle, d is the depth of the mantle, κ is the thermal diffusivity, and η₀ is the reference viscosity. The reference viscosity is the dimensional value of the viscosity low when the mantle is at the reference temperature. At temperatures below the reference temperature, the mantle viscosity rises, while at higher temperatures, viscosity falls. The reference viscosity is therefore a rough, average viscosity for the convecting mantle, neglecting pressure or viscosity law dependencies. Here we use a value of 10²² Pa s, which is reasonable for the Earth's mantle (neglecting the asthenosphere and hydrous material in subduction zones, which would have a few orders of magnitude lower viscosity than the rest of the mantle; Turcotte & Schubert 2002). SSAXC calculates a temperature-dependent viscosity based on the Arrhenius equation:

$$\begin{equation} \eta (T) = \eta _0 e^{E^*\left(\frac{1+T_0}{T-T_0} - 1\right)}, \end{equation} \tag{ 5 }$$

where E* is activation energy and T₀ is the temperature offset, set to 273 K.

Numerical experiments were benchmarked against Earth before being used to model exoplanetary systems. Intense stellar flux and attendant steep temperature gradients near the planetary surface requires high resolution grids, consisting here of 300 tangential nodes by 150 radial nodes. In most models presented here, initial internal mantle temperatures were set to a potential temperature of 1300 K. We assumed an initial lithosphere over the top 5% of the planetary radius, where lithosphere is defined simply by the region where heat is transferred conductively due to high viscosity produced by low surface temperatures. Our initial temperature conditions for this lithosphere consisted of a uniform temperature gradient between the mantle temperature at depth and 273 K at the surface. The evolutionary pathways of super-Earths prior to tidal-lock are not known, and this forms a simple and plausible initial condition. Two additional experiments were conducted with an initial internal mantle potential temperature at 1700 K, which explores the sensitivity of our results on initial internal thermal conditions. No significant difference was observed in the resulting planetary evolution between these two initial mantle potential temperature assumptions.

The stellar flux arriving at the planet is strongest at the substellar point (θ = 0) and decays to zero at the terminator ($\theta = \frac{\pi }{2}$); the far-side ($\frac{\pi }{2} < \theta < \pi$) will receive no stellar flux. Every surface element radiates as a blackbody at each time step. Radiative fluxes are based on the Stefan–Boltzmann equation: F = σT⁴, where F is the radiative energy flux, σ is the Stefan–Boltzmann constant, and T is the temperature of the radiative body or element. From this flux we calculate the temperature of the surface of the solid silicate planet using ΔF = mC_pΔT, where m is the mass (or density times volume: ρ dx² dz), C_p is the specific heat capacity, and ΔT is the change in temperature; finally, dx and dz are the length and depth (respectively) of the surface elements in the numerical model grid.

For the stellar temperature addition, T_in, on the star-side of the planet at each time step, this flux must be multiplied by a term describing the ratio of energy radiated at the surface of the star and that at the planet's distance: $\frac{4 \pi R_{*}^2}{4 \pi a^2}$, where R_* is the star's radius and a is the planet's semimajor axis. This gives

$$\begin{equation} T_{\rm in} = \frac{R_{*}^2 \, \sigma \, T_{*}^4 \, {dt}}{a^2 \, \rho \, {dx}^2 \, {dz} \, C_p}, \end{equation} \tag{ 6 }$$

where T_* is the star's effective temperature, dt is the Courant time step, ρ is the density of rock in the planet, and dx and dz reference the width and depth of an element, respectively. An analogous derivation for the planet's radiative flux out yields

$$\begin{equation} T_{{\rm out}} = \frac{\sigma \, T_{p}^4 \, {dt}}{\rho \, {dx}^2 \, {dz} \, C_p}, \end{equation} \tag{ 7 }$$

where T_p is the temperature of each element at the beginning of the calculation of that time step. Equations (6) and (7) are added and subtracted (respectively) to each temperature calculation for the surface elements on the star-side; only Equation (7) is subtracted from the temperature calculations for the far-side surface elements.

Since this study is focused in understanding the isolated long-term effects of stellar flux and blackbody radiation during synchronous rotation on Earth-like planets, we make several simplifying assumptions. First, internal compositional variations are assumed to be negligible, and therefore a uniform approximate dry peridotite mantle is used. No interior phase changes were included in the mantle. Material input parameters and other constants throughout these numerical experiments are listed in Table 1.

We also do not include atmospheres in these simulations. Any atmospheres on these relatively small planets may be removed by high surface temperatures and stellar bombardment. The presence of an atmosphere on a super-Earth has not yet been observed, though several theoretical approaches for atmosphere detection have been proposed (e.g., Miller-Ricci et al. 2009; Seager & Deming 2009). Importantly, for cooler planets, particularly those orbiting M-stars, the atmosphere-less assumption may be poor, albeit necessary in our models. However, our results can be extrapolated to hotter super-Earths, for which this assumption may be justified. Future work should consider the validity of this assumption, since a significant atmospheric presence would redistribute heat and dampen the effect of a strong, fixed stellar heating.

Finally, previous work has shown that the effect of orbital tides on a planet without a perfectly circular orbit may be a potentially large source for internal heating (Jackson et al. 2008a, 2008b, 2008c; Barnes et al. 2008). Here we assume negligible deviations in both obliquity and eccentricity. Future work should compare the effect of tidal heating with fixed stellar heating to determine the total effect on planets with noncircular orbits at various semimajor axes around massive stars.

Our numerical experiments explore the parameter spaces of both Rayleigh number and stellar flux. A starting Ra was calculated for an Earth-like planet using parameters defined in Table 1. Variations in Ra spanned three orders of magnitude, from 3 × 10⁶ to 1 × 10⁹; this wide range can be produced by changes in several parameters, including mantle viscosity, interior planetary temperature, and planetary radius. Here we assume a temperature-dependent mantle viscosity with a reference value of 10²² Pa s, and an approximately adiabatic mantle at either 1300 K or 1700 K, which allows us to use variations in Ra as a way to investigate varying planetary radius. Stellar flux variations spanned 550–5400 Wm⁻². Note that Earth, and the newly discovered Gliese 581 g, correspond to approximately 1300 W m^–2, and from Table 2, Gliese 581 d and c to 550 and 5400 W m^–2, respectively. A summary of experimental conditions and results can be found in Table 3.

In all experiments the large-scale degree-1 mantle convective structure developed, and we will generally use the timescale associated with development of this structure for comparison and evaluation of our results. However, in some cases, this timescale far exceeded the age of the universe, and hence not all planets in tidal-lock may have developed degree-1 convection yet. For a planet to develop degree-1 convection in under approximately 5 billion years, both the stellar flux the planet receives and the Ra must be at the high end of the ranges explored in this study. This result is examined conceptually in Section 4.1.

In all simulations, we have identified two main evolutionary stages in mantle convection. The first comprises a period of high-degree convection dictated by the Ra. During this stage, upwellings and downwellings are pervasive, while the substellar surface radiatively warms and the antistellar surface radiatively cools while slowly growing thicker. In high Ra regimes, narrow upwellings and downwellings occur throughout the mantle with no clear periodic spacing. At 10 evenly spaced time step panels throughout this stage in the Ra = 1 × 10⁹ simulation, we counted downwellings and characterized them as "star-side" or "dark-side." We find that generally more lithospheric drips occurred on the star-side of the planet than on the dark-side (95 to 77 drips, respectively). Some of these convective cells span the entire depth of the mantle, while others are confined largely to the upper third of the mantle. There did not appear to be any correlation of the depth or width of a lithospheric drip with the hemispheric surface dichotomy. In the low Ra end member (Ra = 3 × 10⁶), as expected, thick, periodic, alternating upwellings and downwellings form.

The second stage of evolution is the formation of degree-1 mantle convection. Figure 1 shows the time evolution of an Earth-like planet receiving a flux corresponding to Gliese 581 c (5400 W m^–2). Figure 2 is the time evolution of a planet with three orders of magnitude higher Ra, but again receiving the same stellar flux as Gliese 581 c. Comparison between these figures shows the important control Ra has on development of degree-1 convection and overall mantle cooling patterns. Since both material properties and mantle geometry are collapsed in the Ra, any combination of parameters which raise this value (such as larger mantle depth, a greater temperature gradient, or lower mantle viscosity) will promote a faster development of degree-1 convection. Based on these results, while stellar flux influences the long-term structure of mantle convection with a hot substellar region and a cool antistellar region, the Ra plays a primary role in determining the timescale of onset of this structure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows near-surface temperatures (one node from the surface boundary) for both Earth-like planets (Ra = 3 × 10⁶) receiving a range of stellar fluxes (colored curves), and planets with larger Ra (black curves). In the case of an Earth-like set of mantle parameters (Ra = 3 × 10⁶), generally, near-surface temperatures remain constant throughout the evolution of the planet. While lithospheric thicknesses change, and hence the geothermal gradient changes, due to the overwhelming control of both stellar heating and radiative heat loss to space, surface temperatures can remain largely stable, even with a striking dichotomy between substellar to antistellar regions. Substellar temperatures range from approximately 600 K in the Gliese 581 d case (flux is 550 W m^–2) to 830 K in the Gliese 581 c case (flux is 5400 W m^–2). Temperatures at the terminator vary far less, from 310 to 350 K, and temperatures at the substellar point are independent of stellar flux. Figure 3 also highlights the range of temperatures which permit stable liquid water at one bar pressure. For all cases at the low Ra end member, temperatures from the substellar point to intermediate star-side regions fall within this tentative range allowing the possibility of stable liquid water at or near the surface. It should be repeated, however, that the presence of an atmosphere, not modeled here, would more evenly distribute heat across the surface of the planet and, depending on the potential atmospheric composition, may also insulate and warm the planet. Our geophysical results would therefore be most applicable for hotter planets whose atmospheres may not be sustainable due to stellar bombardment.

Zoom In Zoom Out Reset image size

The surface temperatures of the high Ra end member receiving a flux equivalent to Gliese 581 c (5400 W m^–2) are plotted in black in Figure 3. Substellar temperatures are approximately 100 K cooler than in the lower Ra case, and have a significant temporary plateau corresponding to the initiation of degree-1 convection. Temperatures at the terminator are constant, but approximately 100 K warmer than the lower Ra case, while antistellar temperatures, on the other hand, are somewhat cooler. Similar to the low Ra cases, there is a region which would tentatively allow for stable liquid water near the surface in the intermediate regions of the "dark" side of the planet.

4.1. Degree-1 Mantle Convection

4.2. Surficial Magma Ponds

While in the models presented here temperatures do not rise above approximately 850 K, several detected exoplanets receive higher stellar fluxes than we are able to model (Table 2) and would be tentatively expected to have significantly higher substellar temperatures. Depending on mineralogy and volatile content of the upper mantle and lithosphere, near-solidus partial melting may begin at temperatures as low as 870°C for felsic compositions, or 1370 K for mafic compositions. If substellar temperatures are high enough, and stable for long periods of time, they may partially melt the surface in the hottest regions of the planet, producing stable, long-term "magma ponds," or localized magma oceans.

Localized magma oceans have been studied previously in order to understand the effects of large impactors early in terrestrial planet evolution (e.g., Tonks & Melosh 1993; Reese & Solomatov 2006; Watters et al. 2009). A significant difference between the dynamics of impact-induced transient molten regions and stellar flux-induced magma ponds is long-term evolution and stability. In the former, crystallization timescales depend on the depth, extent, composition, and, in general, the cooling efficiency of the region, while in the latter, the same is true but, given enough stellar heating, large-scale magma ponds may persist indefinitely. While estimating the dimensions and dynamics of such a pond is complex (for example, considering various compositions and their solidi, extent of melting, melt extraction efficiency, isostatic rebound of the solid floor, and internal convection), if the arriving stellar flux delivers more energy than is lost by the convecting magma pond, we can assume that the feature will be stable over long timescales.

According to Kraichnan (1962), the heat lost from the top of a cooling magma pond, F_m, in a high Prandtl number convective regime is

$$\begin{equation} F_m \sim 0.089 Ra^{1/3}. \end{equation} \tag{ 8 }$$

After Reese & Solomatov (2006), we combine Equations (4) and (8) and redimensionalize to obtain

$$\begin{equation} F_m = 0.089k(T_p - T_s)^{4/3}\left(\frac{\rho _l \alpha g}{\eta _l \kappa }\right)^{1/3}, \end{equation} \tag{ 9 }$$

where k is thermal conductivity, ρ_l is the silicate liquid density, α is the coefficient of thermal expansion, η_l is the silicate liquid dynamic viscosity, κ is the thermal diffusivity, T_p is the potential temperature, and T_s is the surface temperature. To obtain values for these simple calculations, we assume a sufficiently high surface temperature, T_s ≈ 1470 K, such that peridotite would be molten. We assume that the liquid magma pond will have a thin surface boundary layer, with a change in viscosity by a factor of 10 between the boundary layer and the lower convecting material (e.g., Zaranek & Parmentier 2004). A mafic silicate melt at 1470 K at the surface corresponds to a viscosity of ∼10^−0.2; once we decrease this by an order of magnitude, we expect the viscosity of the convecting fluid to be ∼10^−1.2, which correlates to a potential temperature, T_p, of ∼1670 K (Bottinga & Weill 1972, Figure 7, and references therein). A summary of the parameters and values used in these calculations is given in Table 4.

Table 4. Magma Pond Calculation Parameters

Symbol	Parameter	Value
k	Thermal conductivity	4W (m · K)−1
ρl	Silicate liquid density	3300 kg m−3
α	Coefficient of thermal expansion	2 × 10−5 K−1
ηl	Silicate liquid dynamic viscosity	10−1.2 Pa · s
κ	Thermal diffusivity	1 × 10−6 m2 s−1
Tp	Potential temperature	1670 K
Ts	Surface temperature	1470 K

Note. T_p and η_l for magma with peridotite composition obtained from Bottinga & Weill (1972 and references therein).

Download table as:  ASCIITypeset image

To a first order, if the flux released by the partially molten surface remains less than or equal to the flux received from the star, F_s = σT⁴_*, then the magma ocean should be stable over long timescales. Though still neglecting the presence of an insulating atmosphere, a simple calculation of maximum semimajor axes for the stability of an exposed molten surface around various stars which have detected super-Earths is shown in Figure 4. For our solar system, in order for a localized magma ocean to be stable on an Earth-like planet, the semimajor axis must be less than ∼0.12 AU. The known planets which, based on this calculation, are able to sustain a magma ocean for long timescales from fixed stellar heating are CoRoT 7 b, CoRoT-7 c, HD 7924 b, 55 Cnc e, 61 Vir b, HD 1461 b, and HD 40307 b. Though several other exoplanets actually have significantly smaller semimajor axes than the above, they orbit significantly weaker stars (for example, Gliese 581 e is around an M-star) which do not deliver enough stellar flux, even for close-in planets, to match the efficiency of heat loss from a magma ocean.

Zoom In Zoom Out Reset image size

A planet with a large localized magma ocean will become asymmetric based on the differences in density between molten and solid silicates. Analogous to glacial rebound on Earth, the solid mantle underlying the magma pond will isostatically rebound under the lighter load, causing magma to spill out of the pond and spread beyond the region of stellar heat sufficient to melt silicates. The overflowing magma will crystallize, forming lava levées around the magma pond. At the same time, the magma pond will be melting back into the planetary interior to its original depth, determined by stellar flux.

During melting the magma will lose its more volatile components, beginning with water, carbon, and oxygen and sulfur species, and at higher temperatures continuing to sodium, fluorine, chlorine, silicon, and potassium, possibly even to iron (for more details, see Schaefer & Fegley 2009). These species, released into a transient or permanent atmosphere, leave behind an increasingly refractory liquid and solid silicate planet. The magma pond will therefore require increasing stellar flux to produce the same depth of melt. If the volatile components cannot be recycled back into the planetary interior through a process like plate tectonics, then the magma pond will eventually solidify as the planetary interior composition's melting temperature rises above the input from the star. Any degassing of nonrecycled species may be detectable indications of this process. An additional source of mantle degassing may be volcanism, which we would qualitatively expect to be focused near the substellar point. However, the presence or extent of any plate tectonic activity could greatly influence the effect of this potential devolatilization source.

We have explored the effect of synchronous rotation on terrestrial exoplanetary mantle and surface evolution. While planets around M-stars such as Gliese 581 may be expected to have atmospheres in general (which were not modeled here), our results provide a baseline geodynamical framework upon which such complexities may be added. Exploring a range of orbital conditions, from fluxes corresponding to Gliese 581 d to Gliese 581 c, we found two main phases of planetary evolution: a transient stage of pervasive high-degree convection and a culminating period of degree-1 convection. The timescale dictating the transition between these two evolutionary stages is controlled largely by the Rayleigh number. We find a gradient of surface temperatures ranging from 850 K to below 270 K, and hence regions either at the surface, or shallow subsurface, which may support stable liquid water. Lastly, we extrapolate our results to planets receiving higher stellar fluxes. With sufficiently high temperatures, a silicate lithosphere may melt and even support persisting localized magma oceans.

Two aspects of planetary evolution explored here may be directly detectible. A large temperature dichotomy on an atmosphere-less tidally locked planet may be detectible by measuring a thermal phase curve as described by Seager & Deming (2009). Second, with advances in detection techniques, potentially degassing magma ponds may provide a transient or permanent atmosphere which can be detected in the future, as has been done for transiting giant gas planets (e.g., Charbonneau et al. 2002; Vidal-Madjar et al. 2010).

The search for exoplanets has produced an array of exotic planets whose orbital conditions and surface environments are not found anywhere in our solar system. At the core of the exoplanet search is the hunt for habitable environments beyond Earth. Many workers propose that liquid water is a pre-condition for life due to its efficient transport of biochemical material and link to carbon-based biology (e.g., Hart 1978; Kasting & Catling 2003). This has led many to define a habitable zone: a range of semimajor axes in which stellar heating allows an Earth-like planet to have surface temperatures permitting liquid water. Using this definition, weak stars such as Gliese 581 may support tidally locked planets just outside the habitable zone which would normally be dismissed as uninhabitable. However, based on these results, a planet which does not strictly lie within the habitable zone may still support liquid water on its surface, or shallow subsurface, in certain regions of the planet. Thus, instead of assuming or precluding liquid water on planets like Gliese 581 c and d, tidally locked mantle and climate patterns (Joshi 2003) must be combined and assessed to determine the surface environment, keeping in mind this may vary greatly from the substellar to antistellar regions.

We gratefully acknowledge funding from NSF Astronomy CAREER grant 0747154, the MIT Undergraduate Opportunities Program, and from NSF REU. We also acknowledge preliminary work done by Anita Ganesan. This manuscript benefited from the thoughtful review of an anonymous reviewer. Lastly, we thank Steven Kawaler for his editorial efforts.