© ESO 2019

1 Introduction

Theoretical studies of planet formation suggest that during their early phase of evolution, rocky planets undergo multiple magma ocean stages as a result of impacts (e.g., Benz & Cameron 1990; Melosh 1990; Tonks & Melosh 1993), internal heating from the decay of short-lived radioisotopes (Elkins-Tanton 2012), and the release of potential energy during core formation (Flasar & Birch 1973; Sasaki & Nakazawa 1986). As an example, the giant impact that formed the Moon (Hartmann & Davis 1975) melted a large portion of the Earth’s mantle to produce a global magma ocean (e.g., Canup & Asphaug 2001; Nakajima & Stevenson 2015) that subsequently cooled and crystallized. Processes taking place during magma ocean solidification are of fundamental importance for the subsequent evolution of a planet, as they determine its early thermal and chemical structure, atmospheric composition, tectonic behavior, and ultimately its habitability (Massol et al. 2016; Ikoma et al. 2018).

Planets covered by a magma ocean have not yet been observed directly. Nevertheless, hot molten bodies are likely common in the Universe; one such example is the tidally locked super-Earth 55 Cnc e. Heat maps of 55 Cnc e indicate a temperature on the star-facing side that is sufficient to melt silicates (Demory et al. 2016; Hammond & Pierrehumbert 2017). Employing direct imaging techniques to newly forming planets could lead to the detection of another type of molten surface: surface magma oceans that are the product of a recent giant impact. Such a detection would provide direct observational constraints for theoretical models of planet formation, interior and atmospheric dynamics, as well as insights into the origin and diversification of planets in the solar system and extrasolar systems.

The composition of the atmosphere plays an important role in the direct detection of hot molten protoplanets. In general, a primary atmosphere grows by degassing of over-saturated volatiles (e.g., H2O, CO2, CH4, N2) during the crystallization of a magma ocean. The accumulation of greenhouse gases within the pre-existing (primordial) atmosphere exerts a thermal blanketing effect which inhibits heat radiation to space and hence slows down the cooling of the interior relative to no atmosphere being present (Matsui & Abe 1986; Elkins-Tanton 2008; Lebrun et al. 2013; Hamano et al. 2013). The composition and evolution of the outgassed atmosphere thus strongly controls the longevity and observability of molten planets.

On the one hand, prolonged cooling timescales (up to 100 Myr, see Hamano et al. 2013) increase the probability of detection because a planet remains hotter (and hence brighter) for longer. On the other hand, dense and optically thick atmospheres can hinder the direct observation of planetary surfaces. This is because the observed (brightness) temperatureof a planet, as dictated by the convective-radiative equilibrium of its atmosphere, may be lower than its actual surface temperature. Hence magma ocean bodies may appear less bright and thus be more challenging to detect.

Following a giant impact, the molten surface of a hot protoplanet is expected to be between 1000 and 2000 K and thus emits primarily at infrared (IR) wavelengths (Stern 1994; Zhang & Sigurdsson 2003; Mamajek & Meyer 2007). As a magma ocean cools it reaches the so-called “rheological transition” around a critical melt fraction of 40%; at this transition the viscosity of the solid–melt mixture abruptly increases by orders of magnitude (Solomatov 2000; Costa et al. 2009) and the magma ocean behaves rheologically as a solid. Past studies have addressed the detectability of magma ocean planets by predicting the wavelengths that are able to penetrate thick outgassed atmospheres. Miller-Ricci et al. (2009) explored thermal emission spectra of magma ocean bodies for a range of atmospheric compositions, showing that thermal radiation can pass through a series of atmospheric near-IR windows and potentially be observed by future telescopes. Lupu et al. (2014) further investigated high-resolution spectra for systems composed of atmospheres in equilibrium with magmas of different compositions, and concluded that wavelengths of 1–4 μm are the most favorable for the direct detection of such bodies.

Spectral features are expected to evolve during magma ocean cooling due to the evolution of the surface and atmosphere of the hot protoplanets. In particular, the atmosphere is replenished by volatiles outgassing from the interior, and eroded by hydrodynamic escape caused by extreme ultraviolet (XUV) radiation emitted by the young host star. Hamano et al. (2013, 2015) linked the thermal evolution of magma ocean planets to their spectral variations, and found that bodies located beyond a critical orbital distance a_cr (~ 0.8 AU for a solar-type star), so-called type I planets, solidify fast and display thermal emissions that decay within several million years. Most of the water acquired during planetary formation is preserved and forms primordial water oceans. Conversely, on type II planets, located inside the critical orbital distance a_cr, magma oceans can be maintained for much longer–perhaps more than 100 Myr. The timescale depends on the initial water inventory and the hydrodynamic escape process, which progressively leads to the desiccation of the planet. Bodies with extended solidification times emit significant thermal radiation. Based on planet-to-star contrast estimates, Hamano et al. (2015) found that the K_s (2.16 μm) and L bands (3.55 μm) are the most favorable for detection.

With the advent of a number of imminent direct imaging facilities and potential future missions, a comprehensive and quantitative assessment of the detectability of bodies in their early formation stage is desired. Past investigations have focused on magma ocean evolution and its associated spectra without providing clear prospects for future detectability based on telescope specifications. These studies also typically include an advanced description of energy transfer in Earth-like atmospheres at the expense of a more complete model of interior dynamics.

Here, we aim to quantify the likelihood of observing magma ocean planets by convolving the expected occurrence rate of giant impacts with the frequency and ages of close-by stellar associations in order to find the best targets and telescope specifications for future observations. We use N-body simulations (Grimm & Stadel 2014) to constrain the occurrence rate and timing of energetic collisions inducing global magma oceans during terrestrial planet formation around different stellar types. Subsequently, we chart the evolution of the surface temperature of these molten protoplanets using an interior model devised to describe the cooling and crystallization history of rocky bodies (Bower et al. 2018). The detectability of these bodies is then assessed by employing performance estimates for instruments that are to be installed at ESO’s 39-m Extremely Large Telescope (ELT) and a hypothetical space-based mid-infrared interferometer (Kammerer & Quanz 2018; Quanz et al. 2018; Defrère et al. 2018). Finally, combining the results of the aforementioned modeling steps enables us to compute the likelihood of detecting at least one magma ocean planet within nearby young stellar associations.

2 Methods

2.1 Detectability assessment

The ability of future telescopes to detect protoplanetary collisional afterglows is assessed using technical specifications or performanceestimates for ELT instruments and a potential space telescope, called Large Interferometer for Exoplanets (LIFE¹). LIFE is inspired by the Darwin mission concept, which was originally proposed to ESA in 2007 (Leger & Herbst 2007) with the goal of directly detecting and characterizing extrasolar planets at mid-infrared wavelengths using a free-flying interferometer in space (Cockell et al. 2009). However, Darwin, and similarly NASA’s TPF-I concept, was never realized due to technical challenges and uncertainties related to the expected scientific yield. The growing interest in the search for extrasolar life and significantly improved yield estimates based on results from NASA’s Kepler spacecraft (e.g., Winn & Fabrycky 2015) motivate a reconsideration of similar space missions for future explorations (Kammerer & Quanz 2018; Defrère et al. 2018; Quanz et al. 2018).

Young stellar associations represent promising candidates for harboring forming planetary systems subject to energetic encounters between planetary bodies. We consider ten young stellar groups in the solar vicinity (Table 1) that could be searched for magma ocean planets in the future.

We assume that all stars belonging to a given stellar association have the same age as the association and have planets forming around them. The number of G-type members in Table 1 also comprises F- and early K-type stars (up to K4) for simplicity.

We calculate the probability P_MO to detect at least one magma ocean planet in a given set of observed stars in a stellar association using $$P_{\text{MO}}({\lambda }_{\text{cen}},\text{d},{\tau }_{\text{*}},ϵ)=1-{\displaystyle \prod _{i=1}^{i=n_{*}}\left(1-\frac{{\overline{n}}_{\text{GI},\text{i}}\cdot \Delta t_{{\text{MO}}_{\overline{\text{R}}},\text{i}}}{\Delta t_{\text{bin},\text{i}}}\right)},$$(1)

which implicitly depends on λ_cen, the central wavelength of the considered filter, the distance d to the stellar association, the age of the stellar sample τ_*, and the assumed planetary atmospheric emissivity ϵ. n_* is the number of stars of a given spectral type (A, G or M) located inside a given stellar association, according to Table 1. ${\overline{n}}_{\text{GI},\text{i}}$ is the number of detectable giant impacts within a specified time interval of planet formation (Δt_bin,i = 20 Myr), which accounts for the age of the considered stellar association. $\Delta t_{{\text{MO}}_{\overline{\text{R}}},\text{i}}$ indicates the time interval within which a magma ocean planet of radius $\overline{R}$ is bright enough to be directly observed by future instruments. P_MO depends on τ_* through ${\overline{n}}_{\text{GI},\text{i}}$ (depending on the age of a stellar association, and hence the stage of planet formation, planetary objects will experience a different number of magma ocean inducing giant impacts), and on ϵ through $\Delta t_{{\text{MO}}_{\overline{\text{R}}},\text{i}}$ (depending on the thickness of its atmosphere, a magma ocean planet will display a more or less strong radiative signature). ${\overline{n}}_{\text{GI},\text{i}}$ and $\Delta t_{{\text{MO}}_{\overline{\text{R}}},\text{i}}$ are determined from N-body simulations (Sect. 2.2) and a magma ocean evolution model (Sect. 2.3), respectively. Both parameters also depend on the imaging filter wavelength (λ_cen) and corresponding telescope performance estimates (Table 2), as well as on the distance of the targeted stellar association (d).

The telescope parameters most relevant for the detectability assessment are the angular resolution and the sensitivity (Table 2). A successful imaging program requires both high angular resolution and high sensitivity. In this study we do not explicitly consider contrast-limited observations, but comment on the relevance of contrast performance in Sect. 3.4. The limiting angular resolution of a single-aperture telescope is defined by its inner working angle (IWA), that is, the smallest angular distance still allowing a clear spatial separation of the planet from its host star. For the ELT, we assume an IWA of 2λ_cen∕D, where λ_cen is the central wavelength of the filter and D is the aperture size of the telescope. For LIFE, a free-flying nulling interferometer, we assume that a nulling baseline of up to B = 168 m can be accommodated (Leger & Herbst 2007; Cockell et al. 2009), resulting theoretically in a best possible IWA of λ_cen ∕(4B) = 3 mas at a wavelength of 10 μm. To be conservative, the IWAs assumed in our analysis for LIFE are almost a factor of three larger than this estimate. The sensitivity quantifies the faintest signal that a telescope can detect at a certain wavelength within a specified time interval. Increasing the observation time reduces the noise of the thermal background emission and hence enables the detection of progressively fainter bodies. We assume that background-noise limited performance can be achieved at separations as small as the assumed IWAs.

For the ELT we consider two filters from the Mid-infrared ELT Imager and Spectrograph (METIS, Brandl et al. 2016, 2018), one in the L band and one in the N band, with central wavelengths λ_cen of 3.8 μm and 11.6 μm, respectively.In addition, we perform our analysis for a generic instrument operating in the K band filter (λ_cen = 2.2 μm). For the LIFE telescope, we consider two filters, F560W and F1000W, having central wavelengths λ_cen of 5.6 and 10 μm, respectively.The assumed sensitivities for LIFE are based on the faint source detection limits of the Mid-Infrared Instrument (MIRI) installed on the James Webb Space Telescope (JWST), but take into account the reduced throughput of a nulling interferometer (cf. Glasse et al. 2015; Kammerer & Quanz 2018).

In our analysis, a planet is considered as detectable if its angular separation θ = arctana∕d (with a being the planet’s semi-major axis, and d the distanceto the stellar system) exceeds the instrument’s IWA, and if its flux observed in a specific band is higher than the instrument’s sensitivity. Inclination effects are ignored by implicitly assuming face-on orbits. The comparison of angular separation and observed flux with the assumptions for IWA and sensitivity of the considered instruments (Table 2) enables the determination of parameters related to the detectable bodies, such as their size and brightness temperature.

2.2 Timing and energy of giant impacts

We model the giant impact stage of planet formation using the N-body code GENGA (Gravitational ENcounters with Gpu Acceleration, Grimm & Stadel 2014). GENGA is a hybrid symplectic N-body integrator running fully on GPUs and is specifically designed for simulations of planet formation and orbital evolution. We investigate the occurrence rate of magma ocean-inducing impacts by simulating planet formation around A- (2 M_⊙), G- (1 M_⊙), and M-type (0.5 M_⊙) stars during 200 Myr of evolution, focusing on the giant impact phase. Due to the stochastic nature of N-body models (Kokubo et al. 2006; Hoffmann et al. 2017) we perform five simulations for each stellar type. The initial conditions (Table 3) differ between each stellar mass, but are identical for all simulations of a given star type. The initial disk masses and surface density profiles are based on terrestrial planet formation around low-mass stars (Raymond et al. 2007).

The initial planetesimal disk around each star contains a number of planetesimals N_{planetesimals} with total mass M_disk distributed within a circumstellar annulus R_disk. The planetesimals are distributed according to a surface density profile Σ ∝ r⁻¹. During the subsequent orbital evolution, any body passing within 0.1 AU or more than 50 AU from the central star is removed from the N-body integration. In the former case the body is considered to have collided with the central star, whereas in the latter case it is considered ejected. During the subsequent integration, planetesimals and protoplanets grow larger through pairwise collisions between bodies. In our simulations, colliding bodies undergo perfect accretion (i.e., disruptive and partial accretion events are not considered), with the mass of the resulting body corresponding to the sum of the masses of the colliding planetesimals.

To quantify which of these collisions are sufficiently energetic to form magma oceans, we calculate the specific impact energy Q in the center of mass frame as defined by Leinhardt & Stewart (2012). The specific impact energy can be interpreted as the kinetic energy of the impact distributed throughout the mass of the resulting body, $$Q=\frac{\mu v_{\text{imp}}^2}{2\left(m_{\text{imp}}+m_{\text{tar}}\right)},$$(2)

where m_imp is the mass of the impactor, m_tar is mass of the target, v_imp is impact speed, and μ is reduced mass, $$\mu =\frac{m_{\text{imp}} m_{\text{tar}}}{\left(m_{\text{imp}}+m_{\text{tar}}\right)}.$$(3)

We set the critical specific impact energy Q_crit = 0.3 × 10⁶ J kg⁻¹, above which a magma ocean is formed. This value is consistent with energy estimates of the canonical Moon-forming giant impact (Canup & Asphaug 2001). Regarding the frequency of giant impacts, Quintana et al. (2016) consider a giant impact energy threshold Q_crit = 2 × 10⁶ J kg⁻¹, corresponding to the energy needed to strip half of a planet’s atmosphere. Here we consider a smaller value, since atmospheric loss is not needed to generate a magma ocean. We calculate Q for every collision in each simulation to determine the occurrence rate of giant impacts and the radii of the post-impact bodies around each stellar mass.The radii of the post-impact bodies are determined using a mass-radius relationship for rocky planets (Zeng et al. 2016) $$\frac{R}{R_{\oplus }}=1.0154{\left(\frac{M}{M_{\oplus }}\right)}^{1/3.7},$$(4)

where R is the planetary radius, R_⊕ is Earth’s radius (6371km), M is planetary mass, and M_⊕ = 6 × 10²⁴ kg is Earth’s mass. The magma ocean lifetimes of the post-impact bodies are then calculated using an interior evolution model (Sect. 2.3).

2.3 Lifetimes of magma oceans

We investigate the thermal evolution of magma ocean bodies using an energy balance model with self-consistent thermodynamics of melt and solid silicate phases (Bower et al. 2018; Wolf & Bower 2018). This model provides an improved description of the interior dynamics of magma oceans, which are often simplified in previous studies that focus mainly on atmospheric evolution (Elkins-Tanton et al. 2011; Hamano et al. 2013, 2015; Lebrun et al. 2013; Lupu et al. 2014). It includes a newly developed equation of state that provides the thermo-physical properties of silicate melt throughout the pressure-temperature range of Earth’s mantle (Wolf & Bower 2018). The interior evolution model assumes 1D spherically symmetric geometry and accounts for heat transfer by conduction, convection, mixing (latent heat transport), and gravitational separation. We consider bottom-up crystallization where solids first form at the base of the mantle and crystallization proceeds towards the surface.

The rate at which a magma ocean can transfer energy through its outgassed atmosphere constitutes the rate-limiting factor controlling the cooling timescale of the molten body. We consider two different models for the efficiency of energy transport in the atmosphere: (1) a greybody, which parameterizes the insulating effect of the atmosphere using an effective emissivity ϵ (e.g., Elkins-Tanton 2008), and (2) a steam atmosphere parameterization (Zahnle et al. 1988; Solomatov & Stevenson 1993). Moreover, our models include a parameterization of the temperature drop that occurs in the ultra-thin (≲2 cm thickness) thermal boundary layer at the surface of the magma ocean (Abe 1993). The temperature drop scales as the surface temperature cubed (e.g., Reese & Solomatov 2006).

We define the lifetime of a magma ocean as the time required for the surface of an initially molten planet to reach the rheological transition (around 40% melt fraction; Solomatov 2000). Beyond this time it is unclear if a cold lid forms at the surface or if the surface remains partially molten due to the delivery of heat by the convecting magma below. Furthermore, a cool surface leads to a cool upper atmosphere that is expected to harbor thick water clouds for H2 O–CO2 atmospheres, thus efficiently masking the surface (Marcq et al. 2017). Hence, our definition of magma ocean lifetime is designed as a minimum duration to give a conservative estimate of the potential for detectability. Other studies typically define the lifetime as the time taken for the planet to become dominantly (~98%) solid (e.g., Elkins-Tanton 2008).

We compute surface temperature evolutions and magma ocean cooling timescales for planetary radii R ∈ [0.1, 0.3, 0.5, 0.7, 1] × R_⊕ and atmospheric emissivities ϵ ∈ [0.01, $\sqrt{{10}^{-3}}$, 0.1, $\sqrt{{10}^{-1}}$, 0.5, 1] for the gray atmosphere model. For an Earth-like planet with a 1 bar atmosphere and absorption coefficients of H2 O of 0.01 m² kg⁻¹ (Yamamoto & Onishi 1952) and CO2 of 0.05 m² kg⁻¹ (Pujol & North 2003), an emissivity ϵ = 0.01 corresponds to an H2O atmosphere of 1.18 × 10²⁰ kg or a CO2 atmosphere of 5.30 × 10¹⁹ kg. By contrast, ϵ = 1 refers to planetary bodies with no atmosphere (i.e., a blackbody). The obtained cooling timescales are then interpolated for radii 0.1 ≤ R ≤ 1 R_⊕ and emissivities 0.01 ≤ ϵ ≤ 1 using a 2D surface fit.

Magma ocean lifetimes are also calculated for planets with a steam atmosphere (Solomatov & Stevenson 1993). For these models, we back-compute the effective emissivity as a function of time using the surface heat flux and surface temperature calculated by the interior model. The emissivity decreases as a result of volatiles outgassing from the magma ocean into the atmosphere. For planetary radii 0.1 R_⊕, 0.5 R_⊕ and 1 R_⊕, the mean emissivity is $\overline{ϵ}\approx 0.001$, with a standard deviation of σ ≈ 10^-4. Since the average emissivities are comparable, we choose a fixed emissivity of ϵ = 0.001 for the subsequent analyses of the steam atmosphere cases.

For each parameter combination, we estimate planetary brightness temperature for a given surface temperature and observation wavelength. Brightness temperature is a measure of the thermal radiation emitted by a planet and recorded by imaging instruments, and thus directly determines the potential for detection. Employing the definition of emissivity, i.e., the ratio of the spectral radiance at a given wavelength to that of a blackbody (e.g., Catling & Kasting 2017), we map surface to brightness temperature using $$T_{\text{b}}^{-1}=\frac{k\lambda }{hc}\cdot \text{ln}\left(\frac{e^{\frac{hc}{\lambda k{T}_{\text{surf}}}}-1}{ϵ}+1\right),$$(5)

where T_b and T_surf are brightness and surface temperature, respectively, k is Boltzmann’s constant, λ is wavelength, h is Planck’s constant, c is the speed of light, and ϵ is the emissivity of the atmosphere. The brightness temperature is then translated into an estimate of flux density for different time intervals after a giant impact. We determine the time interval $\Delta t_{{\text{MO}}_{\overline{\text{R}}}}$ within which a magma ocean planet is bright enough to be detected (Eq. (1)) by comparing the flux density evolution of a planet of radius $\overline{R}$ with a giventelescope’s sensitivity (Table 2).

3 Results and discussion

3.1 Frequency of giant impacts during planet formation

For each spectral type (A, G, and M), we performed five N-body simulations using the same initial conditions (Table 3). This allows us to constrain the statistical occurrence rate of magma ocean-inducing giant impacts as a function of time for a given star type. The specific impact energies of the collisions occurring in all simulations are shown in Fig. 1A, together with the number of giant impacts for time intervals of 20 Myr shown in Fig. 1B. The latter are median values computed over all simulations for a given spectral type. The dashed horizontal line in Fig. 1A marks the critical specific impact energy Q_crit = 0.3 × 10⁶ J kg⁻¹. Most giant impacts take place within the first ~20 Myr (Fig. 1B), with similar trends for all stellar types. Simulations of systems forming around lower mass stars feature higher initial numbers of planetesimals and smaller protoplanetary disks, thus exhibiting an increased likelihood of close encounters and hence a greater number of giant impacts. We would like to point out that the time binning in 20 Myr intervals introduces a bias to stellar associations (like β Pictoris with an age of 23 Myr) that lie close to edges of these bins. In addition, towards higher ages, the low number statistics for cases like AB Doradus (with an age of 150 Myr) affect the resulting number of giant impacts. Future work is needed to more realistically constrain the giant impact rate for the individual members of each association in order to overcome these limitations.

Previous work investigating the late stages of planet formation (e.g., Leinhardt & Stewart 2012; Quintana et al. 2016) has shown that Earth-like bodies on average experience ≈10–20 giant impacts during accretion. The giant impact rates in our simulations are in good agreement with the aforementioned studies (especially for A- and G-stars), with the noted exception of the first 20 Myr bin. The early impact rates in N-body simulations are strongly dependent on the initial configuration of their constituent planetesimals. As it is notoriously difficult to construct dynamically self-consistent initial architectures for N-body systems, each simulation will experience an initial period of dynamical instability. The high rate of impacts in the first 20 Myr bin reflects this initial period of instability and should therefore be taken as an overestimate. However, while the first 20 Myr bin is biased towards higher impact rates, in general our study tends to underestimate the giant impact rate due to the way collisions are handled. By assuming perfect merging, we are likely underestimating the giant impact rate (Quintana et al. 2016).

Our simulations and those performed in Leinhardt & Stewart (2012) and Quintana et al. (2016) begin from different sets of initial conditions, including different numbers and masses of planetesimals. Thus, the dynamical configurations of the systems are not identical and the initial stages of the simulations are expected to differ, which in particular influences the early impact rate. This can also be interpreted as starting from different stages in the planet formation process. As computing power increases, the early giant impact rate following disk dispersal will become better constrained by studies starting from an increasing number of smaller bodies, allowing the simulated systems to reach more realistic dynamical configurations prior to the onset of the giant impact phase. Future studies will also benefit from the ability to run an increasing number of simulations for each initial condition, allowing researchers to better understand the effect of stochasticity on the giant impact rate.

The radii of the post-impact bodies are shown in Fig. 2. In general, perfect accretion of the colliding bodies results in planets with masses of less than 1.5 M_⊕ and radii less than 1.1 R_⊕. Protoplanetary disk mass scales with stellar mass and, therefore, the amount of material available for planet formation is greater around higher mass stars (Table 3). The post-impact bodies around A-type stars are thus expected to be larger and more massive on average. In terms of detectability, larger planets are favored as they are likely to emit an increased amount of thermal radiation and have extended magma ocean lifetimes. The specific impact energies are normalized by the total mass of the colliding bodies. Hence, even though post-impact bodies are on average smaller around lower-mass stars, they still experience similar specific impact energies to planets around A-type stars (Fig. 1A).

Fig. 1 Panel A: specific impact energies for all simulations (A: A-stars, G: G-stars, M: M-stars) as a function of time. The dashed line denotes the energy threshold (Qcrit = 0.3 × 106 J kg−1), above which we assume that a global magma ocean is generated from the impact. Panel B: median number of giant impacts (GI) occurring within intervals of 20 Myr decays over time due to the decreasing number of planetesimals.

Fig. 2 Median radii of bodies generated from the giant impacts in Fig. 1B for different time intervals of planet formation around different star types (A: A-stars, G: G-stars, M: M-stars). Planetary sizes are computed using a mass-radius relationship for rocky planets (Eq. (4)).

3.2 Magma ocean cooling timescales

Magma ocean lifetimes are important to assess the direct detectability of planets in their early evolution stage, and determine which young stellar associations can be considered for future direct observations. As an example, if the magma ocean lifetime is on the order of millions of years or shorter, young stellar associations with ages < 100 Myr are the prime target for observations. The evolution of surface temperature and resulting magma ocean lifetimes are shown in Figs. 3A and B, respectively. We track the lifetime of the magma ocean from an initially fully molten state until the time when the rheological transition (40% melt fraction) reaches the surface. Larger planets experience longer magma ocean lifetimes due to their reduced surface area to volume ratio (i.e., a larger energy budget relative to radiative heat loss). Moreover, the emissivity strongly influences magma ocean cooling timescales. A low emissivity is indicative of a thick outgassed atmosphere, whereas a high value is indicative of a thin atmosphere. For a given planet size, a decreasing emissivity from ϵ = 1 to ϵ = 0.01 causes the cooling timescale to increase by ~2 orders of magnitude. The longest magma ocean lifetimes are observed for cases with a steam atmosphere (≈ 2–3 × 10⁶ yr).

Previous studies of coupled interior-atmosphere evolution reveal typical magma ocean lifetimes of 10^-2−10⁶ yr (Elkins-Tanton 2008), and possibly up to hundreds of millions of years (Hamano et al. 2013, 2015). Our magma ocean lifetimes (10² −3 × 10⁶ yr) compare well to these previous studies and are typically slightly shorter due to the different definition of the cooling timescale. Elkins-Tanton (2008) and Hamano et al. (2013) define the magma ocean lifetime based on the time required for a fully molten magma ocean to reach the solidus. In our modeling, we calculate the time taken for the rheological transition to reach the surface. Beyond this time, the luminosity of the planet is uncertain since it is unclear whether the planetary surface remains hot and molten or if a cold rigid lid forms. Hence, by design, we compute conservative estimates of magma ocean detectability based on minimum estimates of magma ocean lifetimes. This is further compatible with our exclusion of internal heat sources, which would extend magma ocean lifetimes.

Fig. 3 Panel A: surface temperature evolution of magma ocean planets with different radii and atmospheric emissivities (or steam atmosphere). The horizontal dashed black line denotes the rheological transition (RT), below which a magma ocean may not be detectable due to a cold rigid lid forming at the surface. Hence, the RT criterion gives the minimum lifetime of a magma ocean that is likely to be detectable. Panel B: magma ocean lifetimes ΔtMO obtained for the gray atmosphere model (black circles) and the 2D surface fit for 0.1 R⊕ ≤ R ≤ 1 R⊕ and 0.01 ≤ ϵ ≤ 1.

Fig. 4 Fraction of giant impacts (GI) occurring around A-, G-, and M-stars that can be detected by LIFE and ELT.Here, we only consider the IWA of the instruments without taking into account the effect of the planetary atmosphere and the exposure time. The positions of the colored data points on the x-axis correspond to the distances of the young stellar associations in Table 1. The orbital range that can be resolved by a telescope decreases with increasing distance from the Sun. Hence, the fraction of giant impacts that can be detected is smaller for stellar associations located farther away. Due to its smaller IWA, and thus higher angular resolution, LIFE can image higher fractions of impacts.

3.3 Telescope performance

The IWA and sensitivity of the ELT filters and LIFE (Table 2) are taken as detectability thresholds for constraining the observability of giant collision afterglows. Using the results of the N -body simulations (Sect. 3.1), we determine the fraction of detectable giant impacts within the young stellar groups (Table 1). These have ages ranging from 10 to 150 Myr and are thus associated with different numbers of collisions and post-impact body characteristics according to the time bins in Figs. 1B and 2. Figure 4 illustrates the fraction of detectable giant impacts around A-, G-, and M-type stars as a function of distance from the Sun. We only consider the IWA of the instruments, thus not accounting for atmospheric effects and exposure times, as well as the ages of the stellar associations. For each telescope filter, a data point at a given distance corresponds to one of the ten stellar associations in Table 1.

Due to its small IWA and hence high angular resolution, LIFE is able to capture giant impacts that are close (≈ 0.1 AU) to the central star. Therefore, compared to the ELT, it displays a shallow decrease in the amount of detectable impacts with increasing distance from the Sun. Furthermore, the baseline of the individual free-flying LIFE telescopes could be increased, thus achieving higher spatial resolution. This would enable the imaging of the innermost region of planet forming regions. A- and G-stars both exhibit a similar decay in the number of observable giant impacts. Conversely, for M-type stars the fraction of detectable collisions decreases more rapidly due to the smaller radial extent of the disk (Table 3). This causes giant collisions to occur close to the central star, rendering the task of imaging them challenging for LIFE and almost impossible for the ELT, especially at longer wavelengths.

During magma ocean cooling, the outgoing radiation flux decays due to the decrease of surface temperature. Figure 5 shows flux evolutions following a giant impact for magma ocean bodies located inside the stellar associations β Pictoris (37 pc), TW Hydrae (53 pc), and η Chamaeleontis (94 pc), as they reasonably sample the distance range (up to 100 pc) considered here. Planetary sizes are extracted from Fig. 2, according to the spectral type of the host star and the age of the stellar association. The energy radiated by a planet is quantified by its brightness temperature (Eq. (5)), which depends on the atmospheric emissivity; surface and brightness temperature are equivalent for a planet without an atmosphere. For observation times of five hours and emissivities ϵ = 1, planetary bodies within the three associations are detectable over the entirety or majority of the lifetime of a magma ocean, except for ELT 11.6 μm. However, these planets tend to cool faster, thus displaying shorter time windows during which they are observable ($\Delta t_{{\text{MO}}_{\overline{\text{R}}}}$, see Eq. (1)). By contrast, bodies with lower atmospheric emissivities or a steam atmosphere never or barely exceed the telescopes’ sensitivity thresholds, but could be observable for an extended duration.

If the observation time is increased by a factor of ten (50 h instead of 5 h), the sensitivity thresholds decrease by a factor $\sqrt{10}$, leading to an improved detection of magma oceans. For observations of 50 h, most planets located within TW Hydrae, and η Chamaeleontis would still not be detectable due to their distance from the Sun and hence apparent faintness. However, most of the telescope filters would be able to observe magma oceans within β Pictoris at least once during their lifetimes. Magma ocean bodies with or without an atmosphere are too faint to be detected with the ELT at 11.6 μm due to high background emission at these wavelengths.

Fig. 5 Flux density evolution for planets located within the β Pictoris (37 pc), TW Hydrae (53 pc), and η Chamaeleontis (94 pc) associations during the cooling of a magma ocean. The subplots and subplot panels indicate different telescope filter wavelengths and atmospheric emissivities, respectively. Each set of bars consists of flux evolutions for bodies located inside a given stellar association, orbiting central stars of different spectral type (A, G or M). The colors indicate different time intervals after the occurrence of a giant impact (GI). Each vertical bar terminates according to the cooling timescale of the magma ocean ΔtMO (see Figs. 3A and B). Planetary sizes used for the flux calculations are taken from the data in Fig. 2, based on the spectral type of the host star and the age of the considered stellar association. The black and green horizontal dashed lines mark the sensitivity of the telescope filters for integration times of 5 and 50 h (Table 2), respectively. Magma ocean bodies displaying fluxes greater than these thresholds are detectable.

3.4 Likelihood of detecting magma ocean planets

The probabilities of detecting at least one magma ocean planet in young nearby stellar associations (Table 1) are shown in Fig. 6 for integration times of 5 h and 50 h, for various combinations of atmospheric emissivities and instruments. These were calculated according to Eq. (1). A young stellar age translates to a high number of expected giant impacts (Fig. 1B). Therefore, the younger a stellar association and the more stars of a given type it contains, the higher the probability of detecting at least one magma ocean planet. Bodies with lower emissivities are the most likely to be directly observed due to their long-lived magma oceans, but the atmosphere could be sufficiently dense to prevent the near surface from being imaged. In contrast, a less dense atmosphere will pose a lower barrier to observing the surface, but it will enable the magma ocean to cool faster.

Generally, the likelihoods of detecting hot molten planets within observation times of five hours (Fig. 6A) are highest for the β Pictoris association, due to its proximity to the Sun, its relatively high number of stellar members and its young age (Table 1). Beyond this stellar group, mostly bodies with high atmosphericemissivities (i.e., higher brightness temperatures) can be detected. If planetary fluxes for a given atmospheric emissivity are higher than a given instrument’s sensitivity over the whole magma ocean cooling timescale (e.g., for LIFE 5.6 μm, ϵ = 1), the highest probabilities are obtained for the association containing more stars (Table 1) and/or experiencing the highest number of detectable giant impacts (Fig. 4). The preferred instruments for the detection of magma oceans are ELT 2.2 μm, followed by LIFE 5.6 μm (see Table 4). This assumes that the full number of young planetary systems in star-forming stellar regions can be utilized. If longer integration times of 50 h (Fig. 6B) are allocated, the prospects of detection are improved, especially for planets with low atmospheric emissivities. In both integration scenarios, the most promising stellar targets to be explored are β Pictoris (37 pc, 23 Myr), Columba (50 pc, 42 Myr), TW Hydrae (53 pc, 10 Myr), and Tucana-Horologium (48 pc, 45 Myr).

ELT 2.2 μm is the preferred filter for future magma ocean observations. We remind the reader, however, that in our assumptions for ELT K band observations we assumed that the instrument is equipped with a laser tomography adaptive optics (LTAO) system (see Table 2). Changing this assumption to seeing limited performance or ground-layer adaptive optics (AO) significantly reduces the achieved sensitivity. In addition, at such short wavelengths, it is important to compare the contrast between the stars in a given association and the magma ocean planets. We therefore convert the K band magnitudes of three stars (one A-, one G-, and one M-type star) from the β Pictoris moving group, the most promising stellar group for magma ocean detections, into flux densities (Table 5). We find that stars are around 5–6 orders of magnitude brighter than magma ocean bodies shortly after the occurrence of a giant impact in the best case scenario (i.e., for ϵ = 1) and the contrast further increases by factors of ten going to ϵ = 0.1, ϵ = 0.01 and a steam atmosphere. While it is not trivial to achieve this performance from the ground at the ELT at 2.2 μm, dedicated high-contrast instruments equipped with extreme AO systems aim at achieving near-infrared contrast levels of 10⁻⁸–10⁻⁹ (e.g., Kasper et al. 2010).

Finally, it is worth pointing out that while we were interested in the detectability of magma oceans in a statistical sense, there are individual stars in eachstellar association that have an intrinsically higher probability than others given their smaller-than-average distance. This is particularly true for β Pictoris moving group members that show a large dispersion in their distances.

Fig. 6 Probability of detecting at least one magma ocean planet in nearby young stellar associations (Table 1) for observation times of 5 h (panel A) and 50 h (panel B). The colors and shapes refer to telescope filters and atmospheric emissivities, respectively. The prime stellar target for the detection of magma oceans is the β Pictoris association (37 pc, 23 Myr), followed by Columba (50 pc, 42 Myr), TW Hydrae (53 pc, 10 Myr), and Tucana-Horologium (48 pc, 45 Myr). Due to the low probability of giant impacts at their respective ages (Fig. 1B), no data points are shown for stars located inside the AB Doradus association (20 pc, 150 Myr). ELT 11.6 μm is not able to detect any magma ocean planets.

3.5 Comparison with atmospheric calculations

We compare flux density evolutions and detection probabilities presented in Sects. 3.3 and 3.4 with those obtained using data based on radiative and convective equilibrium (atmospheric) calculations. Figure 7 in Lupu et al. (2014) shows brightness temperature as a function of wavelength, for surface temperatures of 1200, 1600, and 2200 K and surface pressures of p_surf = 10 bar and p_surf = 100 bar. These relations are derived by considering an atmosphere in equilibrium with a magma ocean that has a compositionsimilar to bulk silicate Earth (BSE). Using Eq. (9) in Elkins-Tanton (2008), we calculate the mass of an H2 O atmosphere and its associated surface pressure for atmospheric emissivities of ϵ = 0.001, ϵ = 0.01 and ϵ = 0.1. We consider an H2O atmosphere since this is the dominant molecule in the Lupu et al. (2014) model.

For ϵ = 0.001 and ϵ = 0.1, we obtain surface pressures outside the range considered by Lupu et al. (2014; i.e., p_surf = 232 bar and p_surf = 2 bar, respectively). However, for ϵ = 0.01 we compute p_surf = 23 bar, and by linearly interpolating the results of Lupu et al. (2014) we map surface to brightness temperature. Flux density evolutions for this alternative model are presented in Fig. 7 for ELT and LIFE filters, and the detection probabilities are summarized in Table 6 for the most promising stellar targets. While ELT is the preferred instrument in the previous scenario (using Eq. (5) for the calculation of brightness temperature), in this alternative model the best observation strategy is represented by both LIFE 10 μm and LIFE 5.6 μm. The ELT filters yield substantially lower chances for detection. This dependence on the chosen relation between surface temperatureand brightness temperature provides additional motivation to develop fully coupled models of the interior, surface, and atmosphere to track the evolution of a rocky body during the magma ocean stage.

Fig. 7 Flux density evolutions during the cooling of a magma ocean with an atmospheric emissivity ϵ = 0.01, computed using brightness temperatures from Lupu et al. (2014, Fig. 7). As for Fig. 5, planets in the β Pictoris (37 pc), TW Hydrae (53 pc), and η Chamaeleontis (94 pc) associations are considered. The subplots indicate different telescope filter wavelengths. Each set of bars consists of flux evolutions for bodies in a given stellar association, orbiting central stars of spectral type (A, G or M). The colors indicate different time intervals after the occurrence of a giant impact (GI). Each vertical bar terminates according to the cooling timescale of the magma ocean ΔtMO (see Fig. 3B). Planetary sizes used for the flux calculations are taken from the data in Fig. 2, based on thespectral type of the host star and the age of the considered stellar association. The black and green horizontal dashed linesmark the sensitivity of the telescope filters for integration times of 5 and 50 h (Table 2), respectively.Magma ocean bodies displaying fluxes greater than these thresholds are detectable.

4 Summary and conclusion

During the final stages of terrestrial planet formation, growing planets experience violent pair-wise collisions, known as giant impacts. These collisions release a substantial amount of energy, which can lead to extensive melting of planetary bodies. The post-impact surface magma ocean may be detectable as an “afterglow”. We present a modeling framework that provides first-order predictions of the detectability of such events using performance estimates of next-generation direct imaging instruments, simulations of giant impact occurrence rates during planet formation around A-, G-, and M-stars, and models of magma ocean cooling for planets of various sizes and atmospheric emissivities. With sufficient observing time to fully exploit the number of planetary systems in nearby stellar associations, the observability of giant impacts dominantly depends on two parameters: (1) the expected number of giant impacts for a given stellar age, and (2) the magma ocean lifetime. The latter crucially depends on the atmospheric emissivity and hence the mass and composition of the planetary atmosphere, which controls the rate at which a magma ocean can lose its heat to space. Planets with high emissivities (i.e., thinner atmospheres) are easier to detect, but at the expense of a shorter magma ocean lifetime. Conversely, dense atmospheres delay the cooling of magma oceans but will also hinder the imaging of planetary surfaces.

We find that target selection favoring young and nearby stellar associations, which contain a large number of stars, significantly increases the likelihood of detecting a magma ocean planet. Furthermore, for sufficiently long integration times, the β Pictoris association is best suited for potential future observations of molten protoplanets, followed by Columba, TW Hydrae, and Tucana-Horologium. Within the modeling framework employed in this work, K band filters of the ELT and mid-infrared wavelength filters of a space-based interferometer similar to LIFE are best suited for the exploration of protoplanetary collision afterglows.

Our work highlights the need to explore observing strategies specifically devised for the stars (and their associations) that are most likely to host molten planetary bodies. For integration times of 50 h, we predict that detection probabilities can reach values as high as ≳ 60% for young stellar groups, such as β Pictoris, depending on the statistical distribution of atmospheric properties and volatile constituents of molten protoplanets. Our estimates can be further enhanced by constraining the expected distribution of planetary (and thus atmospheric) compositions within the suggested stellar associations, and by providing more precise stellar ages and estimates of their total mass. Finally, an optimal observation strategy can be devised by focusing on the most prospective stellar members of the considered associations, and further constraining the anticipated giant impact rate per individual star in a stellar group.

Our work motivates the use of ELT instruments to search for molten rocky protoplanets, and highlights the advantages of a space-based mid-infrared nulling interferometer similar to LIFE for future explorations devoted to the study of the formation and early evolution of planetary bodies. Such a telescope would complement upcoming missions that use both direct and indirect exploration approaches for exoplanet science.

Acknowledgements

The authors thank the referee Markus Kasper for a timely and constructive report that lead to considerableimprovements of the manuscript. We would like to thank Taras Gerya for helpful discussions and for reading an earlier version of the manuscript. I.B. thanks Matthieu Laneuville for valuable suggestions, and acknowledges financial support from a scholarship by the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT). I.B. and T.L. thank Patrick Sanan for technical guidance on using the SPIDER code. T.L. thanks Michael Meyer for his persistent enthusiasm about protoplanet collision afterglows and acknowledges financial support from ETH Zurich Research Grant ETH-17 13-1, a MERAC travel grant of the Swiss Society for Astrophysics and Astronomy, and a postdoctoral fellowship of the Swiss National Science Foundation (SNSF). D.J.B. acknowledges SNSF Ambizione Grant 173992. Parts of this work have been carried out within the framework of the National Center for Competence in Research PlanetS supported by the SNSF. This work made use of infrastructure provided by S3IT (www.s3it.uzh.ch), the Service and Support for Science IT team at the University of Zurich. The numerical simulations were analyzed using the open source software environment MATPLOTLIB (Hunter 2007, matplotlib.org) and the ASTROPY package (Astropy Collaboration 2018).

References

All Tables

All Figures

Fig. 1 Panel A: specific impact energies for all simulations (A: A-stars, G: G-stars, M: M-stars) as a function of time. The dashed line denotes the energy threshold (Qcrit = 0.3 × 106 J kg−1), above which we assume that a global magma ocean is generated from the impact. Panel B: median number of giant impacts (GI) occurring within intervals of 20 Myr decays over time due to the decreasing number of planetesimals.

In the text

	Fig. 2 Median radii of bodies generated from the giant impacts in Fig. 1B for different time intervals of planet formation around different star types (A: A-stars, G: G-stars, M: M-stars). Planetary sizes are computed using a mass-radius relationship for rocky planets (Eq. (4)).
In the text

Fig. 3 Panel A: surface temperature evolution of magma ocean planets with different radii and atmospheric emissivities (or steam atmosphere). The horizontal dashed black line denotes the rheological transition (RT), below which a magma ocean may not be detectable due to a cold rigid lid forming at the surface. Hence, the RT criterion gives the minimum lifetime of a magma ocean that is likely to be detectable. Panel B: magma ocean lifetimes ΔtMO obtained for the gray atmosphere model (black circles) and the 2D surface fit for 0.1 R⊕ ≤ R ≤ 1 R⊕ and 0.01 ≤ ϵ ≤ 1.

In the text

Fig. 4 Fraction of giant impacts (GI) occurring around A-, G-, and M-stars that can be detected by LIFE and ELT.Here, we only consider the IWA of the instruments without taking into account the effect of the planetary atmosphere and the exposure time. The positions of the colored data points on the x-axis correspond to the distances of the young stellar associations in Table 1. The orbital range that can be resolved by a telescope decreases with increasing distance from the Sun. Hence, the fraction of giant impacts that can be detected is smaller for stellar associations located farther away. Due to its smaller IWA, and thus higher angular resolution, LIFE can image higher fractions of impacts.

In the text

Fig. 5 Flux density evolution for planets located within the β Pictoris (37 pc), TW Hydrae (53 pc), and η Chamaeleontis (94 pc) associations during the cooling of a magma ocean. The subplots and subplot panels indicate different telescope filter wavelengths and atmospheric emissivities, respectively. Each set of bars consists of flux evolutions for bodies located inside a given stellar association, orbiting central stars of different spectral type (A, G or M). The colors indicate different time intervals after the occurrence of a giant impact (GI). Each vertical bar terminates according to the cooling timescale of the magma ocean ΔtMO (see Figs. 3A and B). Planetary sizes used for the flux calculations are taken from the data in Fig. 2, based on the spectral type of the host star and the age of the considered stellar association. The black and green horizontal dashed lines mark the sensitivity of the telescope filters for integration times of 5 and 50 h (Table 2), respectively. Magma ocean bodies displaying fluxes greater than these thresholds are detectable.

In the text

Fig. 6 Probability of detecting at least one magma ocean planet in nearby young stellar associations (Table 1) for observation times of 5 h (panel A) and 50 h (panel B). The colors and shapes refer to telescope filters and atmospheric emissivities, respectively. The prime stellar target for the detection of magma oceans is the β Pictoris association (37 pc, 23 Myr), followed by Columba (50 pc, 42 Myr), TW Hydrae (53 pc, 10 Myr), and Tucana-Horologium (48 pc, 45 Myr). Due to the low probability of giant impacts at their respective ages (Fig. 1B), no data points are shown for stars located inside the AB Doradus association (20 pc, 150 Myr). ELT 11.6 μm is not able to detect any magma ocean planets.

In the text

Fig. 7 Flux density evolutions during the cooling of a magma ocean with an atmospheric emissivity ϵ = 0.01, computed using brightness temperatures from Lupu et al. (2014, Fig. 7). As for Fig. 5, planets in the β Pictoris (37 pc), TW Hydrae (53 pc), and η Chamaeleontis (94 pc) associations are considered. The subplots indicate different telescope filter wavelengths. Each set of bars consists of flux evolutions for bodies in a given stellar association, orbiting central stars of spectral type (A, G or M). The colors indicate different time intervals after the occurrence of a giant impact (GI). Each vertical bar terminates according to the cooling timescale of the magma ocean ΔtMO (see Fig. 3B). Planetary sizes used for the flux calculations are taken from the data in Fig. 2, based on thespectral type of the host star and the age of the considered stellar association. The black and green horizontal dashed linesmark the sensitivity of the telescope filters for integration times of 5 and 50 h (Table 2), respectively.Magma ocean bodies displaying fluxes greater than these thresholds are detectable.

In the text