Super-Earth LHS3844b is Tidally Locked

For a given orbital configuration the incident stellar flux on the planet can be computed using Kepler's laws. Here we use the same approach as Pierrehumbert (2010). The distance between star and planet is given by

$$\begin{eqnarray}&&r(t)\,=\,a\displaystyle \frac{1-{e}^{2}}{1+e\cos ({k}_{\theta }(t))},\end{eqnarray} \tag{ 1 }$$

where k_θ (t) represents the planet–star orbital angle, a the semimajor axis, and e the planet's eccentricity. The change in k_θ (t) can be found from the conservation of angular momentum,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{k}_{\theta }(t)}{{\rm{d}}t}\,=\,\displaystyle \frac{2\pi }{P}\displaystyle \frac{{[1+e\cos ({k}_{\theta }(t))]}^{2}}{{\left(1-{e}^{2}\right)}^{3/2}},\end{eqnarray} \tag{ 2 }$$

where P denotes the planet's orbital period and k_θ (t) is numerically integrated to obtain the time-varying distance between the star and planet.

The flux absorbed by the planet at latitude x and longitude y is equal to

$$\begin{eqnarray}&&{I}_{0}(x,y,t)\,=\,\displaystyle \frac{{R}_{s}^{2}}{r{\left(t\right)}^{2}}{\mu }_{0}\times \int (1-A(\lambda ,{\mu }_{0})){I}_{s}(\lambda ){\rm{d}}\lambda ,\end{eqnarray} \tag{ 3 }$$

where λ represents wavelength, R_s the stellar radius, I_s (λ) the stellar spectrum, A(λ, μ₀) is a specific form of surface albedo called the directional-hemispherical reflectance (described in more detail below), and μ₀ the zenith angle factor. The zenith angle factor is equal to ${\mu }_{0}=\cos (x-{x}_{0}(t))\cos (y-{y}_{0}(t))$ on the dayside and μ₀ = 0 on the nightside, where x₀ and y₀ are the substellar latitude and longitude, which are generally time dependent. The above formulation assumes R_s /r ≪ 1 so all incoming stellar flux is parallel (Nguyen et al. 2020); LHS 3844b is roughly compatible with this assumption as R_s /r ∼ 0.14. The above formulation also assumes LHS 3844b's obliquity is zero, justified by previous work which showed that obliquity decays quickly around low-mass stars (Barnes et al. 2008; Heller et al. 2011).

The planetary and stellar parameters are taken from Vanderspek et al. (2019), while the stellar spectrum of LHS 3844 is identical to that used in Kreidberg et al. (2019). The stellar spectrum was fitted in Kreidberg et al. (2019) to the star's overall spectral energy distribution; however, its predicted flux of 55.9 mJy in the 4.5 μm Spitzer bandpass did not match the observed flux of 43.4 mJy. Therefore, a scaling constant was introduced to ensure that the stellar flux in the Spitzer bandpass matches the observed value. For consistency the same approach is used here—we use the unscaled stellar spectrum to compute the surface temperature distribution of LHS 3844b, but multiply the stellar spectrum by a scaling constant of 55.9/43.4 ∼ 1.3 (applied to all wavelengths) when calculating the planet–star contrast ratio. Nonetheless, the best way of dealing with the mismatch between modeled and observed fluxes for LHS 3844 remains ambiguous, and one could also address it using alternative methods (Whittaker et al. 2022).

Planets with nonzero eccentricity or in spin–orbit resonances experience tidal dissipation, which leads to internal heating. Here we model the resulting heat flux using the constant time lag model from Leconte et al. (2010). The tidal dissipation rate is given by

$$\begin{eqnarray}&&{E}_{\mathrm{tide}}=2{K}_{p}\left[{N}_{a}(e)-2N(e)\displaystyle \frac{\omega }{n}+{\rm{\Omega }}(e){\left(\displaystyle \frac{\omega }{n}\right)}^{2}\right],\end{eqnarray} \tag{ 4 }$$

where the tidal constants N(e), N_a (e), and Ω(e) describe the dependence of tidal dissipation on eccentricity, given by

$$\begin{eqnarray*}\begin{array}{rcl}N(e) & = & \displaystyle \frac{1+\tfrac{15}{2}{e}^{2}+\tfrac{45}{8}{e}^{4}+\tfrac{5}{16}{e}^{6}}{{\left(1-{e}^{2}\right)}^{6}},\\ {N}_{a}(e) & = & \displaystyle \frac{1+\tfrac{31}{2}{e}^{2}+\tfrac{255}{8}{e}^{4}+\tfrac{185}{16}{e}^{6}+\tfrac{25}{64}{e}^{8}}{{\left(1-{e}^{2}\right)}^{15/2}},\\ {\rm{\Omega }}(e) & = & \displaystyle \frac{1+3{e}^{2}+\tfrac{3}{8}{e}^{4}}{{\left(1-{e}^{2}\right)}^{9/2}}.\end{array}\end{eqnarray*}$$

Here n = 2π/P represents the planet's mean motion, ω its angular velocity, and K_p is a constant that represents the energy dissipated in each tidal cycle,

$$\begin{eqnarray}&&{K}_{p}\,=\,\displaystyle \frac{3}{2}{k}_{p}{\rm{\Delta }}{t}_{p}\left(\displaystyle \frac{{{GM}}_{{\rm{p}}}^{2}}{{R}_{p}}\right){\left(\displaystyle \frac{{M}_{{\rm{s}}}}{{M}_{{\rm{p}}}}\right)}^{2}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{6}{n}^{2}.\end{eqnarray} \tag{ 5 }$$

Above, M_s and M_p represent the star's and planet's mass, R_s and R_p represent the corresponding radii, k_p is the Love number measuring the planet's rigidity and susceptibility to tidal potential (Love 1909), and Δt represents the time lag between maximum tidal potential and the planet's tidal bulge (Leconte et al. 2010).

The parameters Δt and k_p are typically fitted to observations, with k_p = 0.3 and Δt = 630 s providing a good match for present-day Earth's ocean-dominated tidal dissipation, while Δt ∼ 60 s is more appropriate for a rocky planet with tidal dissipation in its mantle only (de Surgy & Laskar 1997). The internal structure of LHS 3844b is uncertain, including whether it might have an internal magma ocean. We therefore use Earth-like parameters as the default for LHS 3844b, but we additionally explore uncertainty in tidal heating by increasing and decreasing Δt by 1 order of magnitude. If the tidal dissipation inside the planet is spatially uniform, then the surface heat flux is ${e}_{\mathrm{tide}}={E}_{\mathrm{tide}}/4\pi {R}_{p}^{2}$.

For planets with nonzero eccentricity and outside of spin–orbit resonances, tidal dissipation is minimized in an orbital state called pseudo-synchronous rotation (Hut 1981). In this state, the planet's rotation nearly matches its orbital period during apoastron, so the planet rotates faster over the course of an orbit than it would in synchronous rotation. The rotation rate for a pseudo-synchronous state is ω_eq/n = N(e)/Ω(e), which simplifies the tidal dissipation equation to

$$\begin{eqnarray}&&{e}_{\mathrm{tide}}=2{K}_{p}\left[{N}_{a}(e)-\displaystyle \frac{N{\left(e\right)}^{2}}{{\rm{\Omega }}(e)}\right]/(4\pi {R}_{p}^{2}).\end{eqnarray} \tag{ 6 }$$

Although it has been argued that rocky planets cannot be in pseudo-synchronous rotation because this state is incompatible with realistic rheological models of solids (Makarov 2013), we still consider pseudo-synchronous rotation for LHS 3844b because it provides a lower bound on the planet's tidal dissipation. If LHS 3844b was not in pseudo-synchronous rotation it would be captured into another steady state, such as synchronous rotation or a Mercury-like spin–orbit resonance (note, a planet in synchronous rotation with nonzero eccentricity will undergo Moon-like librations). For all of these states, the planet's tidal dissipation would be even larger, which would only strengthen the constraints we derive from the planet's observed phase curve.

Heat storage and heat redistribution inside the subsurface are modeled using the thermal diffusion equation, ∂T/∂t = k/(c ρ)∇² T, where the thermal conductivity k and the volumetric heat capacity c ρ can be combined into the thermal diffusivity α = k/c ρ. On planetary scales, horizontal heat transport due to heat diffusion can be ignored (Seager & Deming 2009), so we only consider vertical heat diffusion,

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}\,=\,\alpha \displaystyle \frac{{\partial }^{2}T}{\partial {z}^{2}}.\end{eqnarray} \tag{ 7 }$$

To generate the planet's surface temperature map we numerically solve the thermal diffusion equation using the CLIMLAB climate modeling toolkit (Rose 2018) on a 360 × 180 × 10 longitude–latitude-depth grid, where the depth points are linearly distributed between the surface and the model's bottom boundary. For planets not in synchronous rotation, we set the bottom boundary to 10 times the thermal skin depth (Spencer et al. 1989). For planets in synchronous rotation the thermal skin depth is infinite, so we set the bottom boundary to 10 times the thermal skin depth of a planet in 2:1 rotation. The top and bottom boundary conditions at each latitude–longitude point are

$$\begin{eqnarray}&&k{\left(\displaystyle \frac{\partial T}{\partial z}\right)}_{z=0}\,=\,{I}_{0}-{I}_{p},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&k{\left(\displaystyle \frac{\partial T}{\partial z}\right)}_{z={z}_{\max }}\,=\,-{e}_{\mathrm{tide}}.\end{eqnarray} \tag{ 9 }$$

Here I₀ is the absorbed stellar flux and I_p is the emitted thermal flux

$$\begin{eqnarray}&&{I}_{p}\,=\,\int \varepsilon (\lambda )\pi B(\lambda ,{T}_{z=0}){\rm{d}}\lambda .\end{eqnarray} \tag{ 10 }$$

where B(λ, T) is the Planck function and ε(λ) is the hemispherical emissivity (described below). Thermal observations of the Moon and Mercury can be fitted with a surface heat capacity ${\rm{\Gamma }}=\sqrt{k\rho c}$ of about 50–80 J s^−1/2 K⁻¹ m⁻² (Winter & Saari 1969; Murdock 1974), and other data suggest Γ ∼ 100 J s^−1/2 K⁻¹ m⁻² for regolith surfaces (Pierrehumbert 2010), so here we use Γ = 100 J s^−1/2 K⁻¹ m⁻². Tests show that our results are insensitive to the exact value of Γ (see Appendix A).

Modeling the reflectance and emissivity of a bare rock surface poses a complex radiative transfer problem. Physically, an atmosphere-less body like the Moon or Mercury has a surface composed of loose unconsolidated material termed regolith (exoplanet surfaces composed of solid rock are possible but unlikely, given that space weathering via micrometeorites and stellar irradiation would convert these surfaces into regolith; see below). Stellar flux enters the regolith and then gets scattered, absorbed, and/or reemitted multiple times. This leads to a surface reflectance and emissivity that strongly depends on the surface's orientation with respect to the star and the observer, the regolith's composition, as well as interactions between regolith particles such as shadow hiding and backscattering (Hapke 2002).

Here we use a modeling approach originally developed by Hapke for solar system bodies like the Moon and Mercury (Hapke 1981, 2002, 2012), and which was subsequently applied to exoplanets (Hu et al. 2012). In this approach, one first infers the single scattering albedo spectrum of a regolith particle based on laboratory reflectance measurements. The single scattering albedo is then used in an analytical model that approximates the emergent flux of an illuminated semi-infinite regolith surface. Following Hapke (2012), the surface's bidirectional reflectance function is given by

$$\begin{eqnarray}\begin{array}{rcl}r(\lambda ,{\mu }_{0},{\mu }_{e}) & = & \displaystyle \frac{w(\lambda )}{4\pi }\displaystyle \frac{{\mu }_{0}}{{\mu }_{0}+{\mu }_{e}}[{B}_{S0}{B}_{S}({\mu }_{0},{\mu }_{e})\\ & & +\,H({\mu }_{0}/K)H({\mu }_{e}/K)].\end{array}\end{eqnarray} \tag{ 11 }$$

Here w(λ) is the wavelength-dependent single scattering albedo, μ₀ and μ_e are the cosine angles with which radiation enters and leaves the surface, B_S0 and B_S are parameters that take into account the opposition effect, H is an analytical function that arises from solving the radiative transfer equations, and K primarily depends on the size of regolith particles; for comparison, the reflectance function for a Lambert surface is r = w(λ)/(4π). We set B_S0 = 0 (no opposition effect) and K = 1 to minimize the albedo at short wavelengths and maximize the observed secondary eclipse depth in the infrared. In addition, Equation (11) assumes regolith particles scatter isotropically.

To solve the surface's energy balance one has to derive relevant albedo and emissivity quantities from the bidirectional reflectance function (Equation (11)). The relevant surface albedo in Equation (3) is the directional-hemispherical reflectance,

$$\begin{eqnarray}&&A(\lambda ,{\mu }_{0})={\int }_{0}^{1}r(\lambda ,{\mu }_{0},{\mu }_{e})/{\mu }_{0}\cdot 2\pi {\mu }_{e}d{\mu }_{e},\end{eqnarray} \tag{ 12 }$$

which describes what fraction of a collimated light beam with incident angle μ₀ is scattered back over all outgoing angles. The relevant emissivity in Equation (10) is the hemispherical emissivity ε(λ), which describes the surface's emission over all outgoing angles. Applying Kirchhoff's law to the outgoing hemisphere, this quantity is

$$\begin{eqnarray}&&\varepsilon (\lambda )\,=\,1-{\int }_{0}^{1}A(\lambda ,{\mu }_{0})\cdot 2\pi {\mu }_{0}d{\mu }_{0}.\end{eqnarray} \tag{ 13 }$$

We compute the above quantities using single scattering albedos that were derived based on laboratory measurements in Hu et al. (2012). To validate our model we compared it to a set of bidirectional reflectance simulations performed with the numerical model from Hu et al. (2012). The two models agree fairly closely, with some remaining model differences that could be either due to round-off errors in fundamental constants or different choices in model resolution (see Appendix B).

Space weathering refers to alterations that occur in the space environment over time (Pieters & Noble 2016). On an exposed surface, impacts by micrometeorites and stellar irradiation lead to comminution and melting. Besides reworking the surface and generating a regolith, these alterations can also cover the surface with molten ferric material called nanophase metallic iron. Similarly, micrometeorites can deliver exogenous material and thus modify the composition of the surface's top layer (Syal et al. 2015). In the solar system, these processes tend to lower the albedo of atmosphere-less bodies and modify their spectral signatures, so they should also be considered for atmosphere-less exoplanets. The Moon's low surface albedo is due to weathering by nanophase metallic iron (Cassidy & Hapke 1975; Pieters et al. 2000), while Mercury's low surface albedo has been attributed to mixed-in graphite (Syal et al. 2015), so we consider both possibilities here.

We first assume LHS 3844b's surface is primarily composed of basalt. Previous analyses of LHS 3844b showed that its secondary eclipse depth can be variously matched by basalt, pyrite (also previously called metal-rich), or a mixture of basalt and nanophase hematite (also called Fe-oxidized) (Kreidberg et al. 2019; Whittaker et al. 2022). Of these, basalt is the most plausible candidate (Hu et al. 2012); the lunar maria are basaltic (Pieters 1978), while the other surface types proposed for LHS 3844b are not found on atmosphere-less bodies in the solar system and so would require invoking unusual formation scenarios.

Next, to simulate how space weathering modifies the surface's spectral properties we follow the same approach as Hapke (2001) and Hapke (2012). For simplicity, we only consider uniform surfaces, so space weathering affects all latitude–longitude points equally. At a given point, the surface is assumed to consist of a mixture of host material (basalt) and absorbing particles. One then computes an effective single scattering albedo for the mixture as follows. The relationship between the absorption coefficient α and the single scattering albedo w of either host material or mixture is (Hapke 2012)

$$\begin{eqnarray}&&{\alpha }_{h}\,=\,\displaystyle \frac{1}{D}\mathrm{ln}\left[{S}_{t}+\displaystyle \frac{(1-{S}_{e})(1-{S}_{t})}{{w}_{h}-{S}_{e}}\right],\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\alpha }_{w}\,=\,\displaystyle \frac{1}{D}\mathrm{ln}\left[{S}_{t}+\displaystyle \frac{(1-{S}_{e})(1-{S}_{t})}{{w}_{w}-{S}_{e}}\right],\end{eqnarray} \tag{ 15 }$$

where subscript h refers to the host material, subscript w refers to the mixture, D is the mean photon path length through the host material, w_h and w_w are the single scattering albedos of the host material and mixture, and S_e and S_t are angle-averaged Fresnel reflection coefficients of the host material (see Hapke 2001). The single scattering albedo of the mixture w_w can then be solved for using

$$\begin{eqnarray}&&{\alpha }_{w}\,=\,{\alpha }_{h}+36\pi z\phi /\lambda ,\end{eqnarray} \tag{ 16 }$$

where ϕ is the volume fraction of absorbing particles, and z is related to the refractive indices of host and absorbing material, defined as

$$\begin{eqnarray}&&z\,=\,\displaystyle \frac{{n}_{h}^{3}{n}_{m}{k}_{m}}{{\left({n}_{m}^{2}-{k}_{m}^{2}+2{n}_{h}^{2}\right)}^{2}+{\left(2{n}_{m}{k}_{m}\right)}^{2}},\end{eqnarray} \tag{ 17 }$$

where n_h is the real part of the host material's refractive index, while n_m and k_m are the real and imaginary parts of the absorbing material's refractive index (Hapke 2001, 2012). The refractive indices of graphite and iron nanoparticles are evaluated using data from the Refractive Index database (Polyanskiy 2022).

In comparing models against data for LHS 3844b there are two options. Previous work compared various physical models to the best-fit secondary eclipse depth of 380 ± 40 ppm (all error bars in this work are 1σ uncertainty), which itself was derived from the observed Spitzer phase curve (Kreidberg et al. 2019; Whittaker et al. 2022). However, we find that the previously derived secondary eclipse depth only provides a modest constraint for discriminating between models, as all surface models we consider in detail match this depth to within 2σ (see below). Moreover, the best-fit secondary eclipse depth of 380 ± 40 ppm reported in Kreidberg et al. (2019) was based on a phase curve model with a day–night flux amplitude of 350 ± 40 ppm, and thus implicitly requires a best-fit nightside flux slightly larger than zero (although the data also matches zero nightside flux within 1σ). A second option is to directly compare models to the observed Spitzer phase curve based on χ². Doing so allows us to include all out-of-eclipse data points in the model-data comparison, and moreover, avoids any assumptions about the planet's nightside flux, which are implicit when using the derived secondary eclipse depth from Kreidberg et al. (2019). When reporting χ² and ${\chi }_{\nu }^{2}$ we use the binned phase curve from Kreidberg et al. (2019), which covers 24 orbital phase bins and is normalized with respect to the secondary eclipse, for ν = 23 degrees of freedom.

To simulate phase curves we assume an edge-on orbit. For noncircular orbits, one additionally needs to specify when the secondary eclipse occurs relative to the periastron. We assume LHS 3844b's periastron coincides with a secondary eclipse (orbital phase of 0.5) and its apoastron coincides with transit (orbital phase of 0). We do not consider alternative viewing geometries as our thermal modeling analysis constrains the planet's eccentricity to less than ∼10⁻³ (see below), for which the impact of viewing geometry is negligible. Neither do we attempt to constrain LHS 3844b's eccentricity from its observed transit duration or the timing between secondary eclipse and transit, as eccentricity constraints from these methods are typically much less precise (Van Eylen & Albrecht 2015; Sagear & Ballard 2023) than the constraint we derive from thermal modeling.

The planet's phase curve is composed of reflected and emitted light. Unlike the angle-integrated albedo and emissivity quantities that are used to compute the surface's energy balance (see above), the observed albedo and emissivity of a given latitude–longitude patch also depend on the patch's orientation with respect to the observer. We therefore compute observed reflected and emitted light components using the bidirectional reflectance function (Equation (11)), which is separately applied to each visible latitude–longitude patch on the planet's surface.

Finally, to compute the observed planet-star flux contrast in the Spitzer bandpass we use a spectral average weighted by the Spitzer response function (IRAC Instrument and Instrument Support Teams 2013). Note that in computing photometric contrasts the order of the spectral averaging procedure matters. Here we spectrally average fluxes first before taking their ratio

$$\begin{eqnarray}&&\mathrm{contrast}\,=\,\displaystyle \frac{\overline{{F}_{p}(\lambda )}}{\overline{{F}_{s}(\lambda )}}\,=\,\displaystyle \frac{\int {F}_{p}(\lambda )f(\lambda )d\lambda }{\int {F}_{s}(\lambda )f(\lambda )d\lambda },\end{eqnarray} \tag{ 18 }$$

which is not the same as the alternative of taking the ratio first before spectrally averaging, $\overline{{F}_{p}(\lambda )/{F}_{s}(\lambda )}\ne \overline{{F}_{p}(\lambda )}/\overline{{F}_{s}(\lambda )}$. The above overbars denote a weighted spectral average over the Spitzer bandpass, F_p and F_s are the observed planet and star fluxes, and f(λ) is the Spitzer response function.

We first revisit the conclusion of Kreidberg et al. (2019) that the secondary eclipse of LHS 3844b is explained well by a planet in synchronous rotation, with zero eccentricity, and a basalt surface. Figure 1 shows that although this scenario is consistent with the planet's secondary eclipse, it does not provide the best fit to the overall phase curve. Focusing on the secondary eclipse, a pure basalt surface would have an eclipse depth of 339 ppm, which compares well against the previously derived secondary eclipse depth of 380 ± 40 ppm (see Table 1). For comparison, a hypothetical blackbody would produce a secondary eclipse depth of 434 ppm, which also matches the previously derived secondary eclipse depth to better than 2σ. Nevertheless, the fit to the overall phase curve is notably different between these two models. For the basalt surface, we find a fit of only χ² = 75.78 (${\chi }_{\nu }^{2}=3.3$, p-value = 1.5 × 10⁻⁷), whereas a blackbody would have a much better fit of χ² = 44.0 (${\chi }_{\nu }^{2}=1.9$, p-value = 5 × 10⁻³). The difference is visible in Figure 1, as a basalt surface matches the observed fluxes clearly worse than a blackbody, particularly around quadrature (near orbital phases 0.25 and 0.75).

Zoom In Zoom Out Reset image size

Based on Figure 1, one might further expect that an asymmetric surface model could fit the data even better than a symmetric one; the observed flux immediately before the secondary eclipse is ∼50–100 ppm lower than the flux right after the eclipse. However, the apparent dayside asymmetry is likely an artifact. First, previous work found that the data bin was down as expected based on photon noise statistics, whereas a robust asymmetry should have led to a noticeable deviation from photon noise at large bin sizes (see Extended Data Figure 3 in Kreidberg et al. 2019). Second, Kreidberg et al. (2019) reported that the best phase curve fit was consistent with no hot spot offset, again supporting symmetric models. In this work, we therefore focus on models with uniform surface properties, such that any phase curve asymmetry only enters through the surface's heat capacity (in practice this effect is negligible; see below).

Overall, while a synchronously rotating planet with a basalt surface is thus able to match LHS 3844b's reported secondary eclipse depth, its poor fit to the planet's phase curve prompts us to explore alternative hypotheses that could also match the data.

Next, we revisit LHS 3844b's rotation and eccentricity. To establish an upper bound for both parameters, we compare the observed dayside brightness temperature of 1040 ± 40 K (Kreidberg et al. 2019) with the planet's theoretical dayside equilibrium temperature T_day. Tidal dissipation resulting from nonzero eccentricity or nonsynchronous rotation will cause the equilibrium temperature to rise and possibly exceed the observed dayside brightness temperature, which we use here to limit the plausible parameter space.

The dayside equilibrium temperature, incorporating tidal dissipation, is equal to

$$\begin{eqnarray}&&\sigma {T}_{\mathrm{day}}^{4}={e}_{\mathrm{tide}}+\displaystyle \frac{{R}_{s}^{2}}{2{\overline{r}}^{2}}\int \int (1-A){\mu }_{0}{I}_{s}d\lambda d{\mu }_{0}.\end{eqnarray} \tag{ 19 }$$

This form assumes a uniform dayside and no heat redistribution to the nightside. The surface albedo is computed using the single scattering albedo of basalt (for simplicity here we replace A's angular dependence with Lambert scattering), $\overline{r}$ is the time-averaged planet–star distance, and e_tide is the tidal dissipation, which is a function of planetary eccentricity and rotation rate. As e_tide also depends on the planet's internal structure, which is uncertain, we vary the tidal time lag parameter Δt (see Equation (5)) by 1 order of magnitude upward and downward.

Figure 2 shows the planet's dayside equilibrium temperature T_day for a planet that is either pseudo-synchronized (blue curve) or in a Mercury-like 3:2 spin–orbit resonance (green curve). The colored regions indicate the impact of varying Δt by 1 order of magnitude. We find that, regardless of the exact value of Δt, if LHS 3844b were in a 3:2 resonance, tidal heating would raise the dayside temperature to over 7000 K, which is significantly higher than its host star temperature. This large level of heating is inconsistent with both the Spitzer data and the lack of stellar-temperature companions in the original TESS phase curve (Vanderspek et al. 2019).

Zoom In Zoom Out Reset image size

Additionally, Figure 2 constrains the planet's eccentricity. Assuming LHS 3844b is in pseudo-synchronous rotation and for the default value of Δt (solid blue curve), the eccentricity must be smaller than 0.002 for T_day to be within 3σ of the dayside brightness temperature. Even if Δt is decreased by 1 order of magnitude to allow for slightly larger eccentricities, only eccentricities below 0.005 are permitted. For comparison, the Moon and Io have eccentricities of 0.05 and 0.004, respectively (Peale et al. 1979; Meeus 1997). This estimate is conservative as it assumes pseudo-synchronous rotation. As explained above, if LHS 3844b was in synchronous rotation the resulting tidal dissipation would be higher, and the permitted eccentricity lower. Even without any further analysis, the absence of large observed tidal dissipation on LHS 3844b thus strongly rules out Mercury-like rotation states and constrains its orbit to be about as circular as Io's.

For the nonzero eccentricities consistent with LHS 3844b's dayside brightness temperature in Figure 2, which one provides the best fit to the Spitzer phase curve? To address this question we use our thermal model to simulate detailed latitude–longitude temperature maps and phase curves. The calculations assume the planet is pseudo-synchronized, its surface composition is basalt, and Δt = 630 s. Figure 3 shows that as eccentricity increases, the flux curves all rise in response to increasing tidal dissipation. If we broadly consider phase curves with ${\chi }_{\nu }^{2}\lt 2$ acceptable, so their fit is at least as good as that of a zero-eccentricity blackbody (see Table 1), then the acceptable eccentricity range for LHS 3844b is 4 × 10⁻⁴ < e < 10⁻³, or less than 25% of Io's eccentricity. Within this eccentricity range, LHS 3844b's rotation has to be slower than ω/n < 1.000006, or about one planet rotation every 200 yr. This means even if LHS 3844b was in pseudo-synchronous rotation, its rotation would be so slow as to be effectively tidally locked.

Zoom In Zoom Out Reset image size

Figure 3 additionally shows that, for a surface made of basalt, LHS 3844b's phase curve also requires tidal dissipation. The best fit occurs at e = 0.0008, that is, at nonzero eccentricity and tidal heating. The tidal heat flux associated with this eccentricity is about 10,700 W m⁻², significantly smaller than LHS 3844b's stellar constant of about 95,600 W m⁻². The resulting fit, χ² = 32.74 (${\chi }_{\nu }^{2}=1.4$, p-value = 0.09), is significantly better than the above-described fit for a planet with a basalt surface and zero eccentricity, χ² = 75.78 (see Table 1). Based on LHS 3844b's extremely short tidal circularization timescale of ∼500 yr (see Section 1), a nonzero eccentricity for LHS 3844b almost certainly requires the presence of a planetary companion. We discuss this possibility in Section 4. However, a surface made of pure basalt might not provide the best fit to LHS 3844b's phase curve either, so we revisit the planet's surface composition next.

Motivated by the Moon and Mercury, we consider how space weathering via nanophase metallic iron particles (npFe) and graphite would affect LHS 3844b's surface. What is the maximum impact these processes could have? If space weathering primarily reduces surface iron into npFe (Domingue et al. 2014), then an upper bound for LHS 3844b is the limit in which all surface iron has been reduced. Earth's crust contains ∼5% iron, which is on the same order of magnitude as the surface iron abundances of the Moon, Mercury, and Mars (Fleischer 1954; Lucey et al. 1995; Robinson & Taylor 2001; Taylor et al. 2006). A highly space-weathered surface with a composition similar to solar system bodies can therefore reasonably contain up to ∼5% npFe. Correspondingly, we compute the surface albedo of LHS 3844b for mixtures ranging from an unweathered surface (pure basalt) up to highly weathered surfaces with a 5% volume fraction of either npFe or graphite.

The left column of Figure 4 shows the single scattering albedo spectrum for different amounts of space weathering. Regardless of whether space weathering occurs via npFe particles or graphite, albedo is always dramatically reduced at wavelengths shorter than 5 μm. The change depends nonlinearly on the amount of darkening material—for example, at 1 μm only 0.5% of added npFe is necessary to reduce the surface's single scattering albedo by half, but the effect saturates above 5% npFe. The nonlinearity arises because small amounts of a dark absorber inside a reflective medium have a nonlinear effect on the mixture's overall absorption (Hapke 2001). Per volume fraction of absorbing material, graphite is more effective at reducing albedo than npFe particles. In the near-infrared the surface's single scattering albedo still exceeds 0.5 with 5% of npFe added, whereas the single scattering albedo is lowered to less than 0.3 at all wavelengths with 5% of graphite.

Zoom In Zoom Out Reset image size

The right column of Figure 4 shows the corresponding thermal phase curves, assuming LHS 3844b has zero eccentricity and is synchronously rotating. Because space weathering reduces the surface albedo, it always increases the phase curve amplitude. In all cases, the largest flux originates at the substellar point, consistent with fits to the Spitzer phase curve, which constrain the hot spot to lie within −6° ± 6° of the substellar point (Kreidberg et al. 2019).

Compared to a pure basalt surface, we find that models with space weathering are essentially indistinguishable in terms of their fit to the previously derived secondary eclipse depth, but they significantly improve the fit to the overall phase curve. Table 1 shows that all models with space weathering match the secondary eclipse depth of 380 ± 40 ppm to better than 2σ. In contrast, the fit to the phase curve improves monotonically with the addition of space weathering. For reference, for a pure basalt surface χ² = 75.78 (see Table 1). Strong space weathering via the addition of 5% npFe improves the phase curve fit to χ² = 53.07 (${\chi }_{\nu }^{2}=2.4$, p-value = 4 × 10⁻⁴), while strong space weathering with 5% graphite produces an even better fit, χ² = 46.48 (${\chi }_{\nu }^{2}=2.0$, p-value = 3 × 10⁻³). Note that the latter model with 5% graphite produces a phase curve that is almost identical to that of a blackbody (χ² = 46.48 versus χ² = 44.0; see Figure 4). We conclude that, similar to tidal heating, the addition of space weathering significantly improves the fit to LHS 3844b's phase curve.