Super-Earth LHS3844b is Tidally Locked

Short-period exoplanets on circular orbits are thought to be tidally locked into synchronous rotation. If tidally locked, these planets must possess permanent day- and night-sides, with extreme irradiation on the dayside and none on the nightside. However, so far the tidal locking hypothesis for exoplanets is supported by little to no empirical evidence. Previous work showed that the super-Earth LHS 3844b likely has no atmosphere, which makes it ideal for constraining the planet’s rotation. Here we revisit the Spitzer phase curve of LHS 3844b with a thermal model of an atmosphere-less planet and analyze the impact of nonsynchronous rotation, eccentricity, tidal dissipation, and surface composition. Based on the lack of observed strong tidal heating we rule out rapid nonsynchronous rotation (including a Mercury-like 3:2 spin–orbit resonance) and constrain the planet's eccentricity to less than ∼0.001 (more circular than Io's orbit). In addition, LHS 3844b’s phase curve implies that the planet either still experiences weak tidal heating via a small-but-nonzero eccentricity (requiring an undetected orbital companion), or that its surface has been darkened by space weathering; of these two scenarios we consider space weathering more likely. Our results thus support the hypothesis that short-period rocky exoplanets are tidally locked, and further show that space weathering can significantly modify the surfaces of atmosphere-less exoplanets.


INTRODUCTION
Exoplanets on short period orbits are widely believed to be tidally locked to their host stars.Tidal interactions between planet and star tend to align the planet's rota-tion period about its own axis with its orbital period.The final outcome of this process is that the planet's orbital and rotation periods become equal, a state called 1:1 tidal locking or synchronous rotation (Dole 1964;Kasting et al. 1993;Barnes 2017).In such a state the planet always faces the star with the same hemisphere, creating a permanent day-and nightside.The resulting instellation asymmetry is expected to significantly impact the climates of short period exoplanets, including their thermal emission maps, their atmospheric dynamics, and their potential habitability (e.g., Joshi et al. 1997;Showman & Guillot 2002;Yang & Abbot 2014).
Not all short period exoplanets need to be synchronously rotating, however.A planet's spin evolution is complicated by many factors including the planet's internal structure, the presence of an atmosphere or ocean, and the potential influence of companion planets (Barnes 2017).As a historical example, Mercury was long believed to be synchronously rotating based on early tidal theory and observations of its surface features (Lowell 1898;Sheehan et al. 2011).Despite these expectations, radar observations decades later revealed that Mercury is in a 3:2 spin-orbit resonance instead (Colombo & Shapiro 1965).
Similarly, tidal models predict a wide range of possible spin states for short period exoplanets.These states include Mercury-like spin-orbit resonances (Makarov 2012;Makarov et al. 2012); chaotic rotation (Correia & Delisle 2019); and complex non-synchronous states for planets with atmospheres (Leconte et al. 2015) or oceans (Auclair-Desrotour et al. 2019), with similar nonsynchronous states plausible for planets that have internal magma oceans or other fluid layers.Whether any particular exoplanet is in such a state is uncertain.For example, Mercury-like spin-orbit resonances are not possible if a planet's eccentricity is zero.However, tidal models permit stable spin-orbit resonances at eccentricities as low as e = 0.001−0.01(Correia et al. 2014;Correia & Delisle 2019), which is smaller than the precision to which one can typically constrain exoplanet eccentricities from radial velocity measurements or transit light curves (Shen & Turner 2008;Kane et al. 2012;Van Eylen & Albrecht 2015;Sagear & Ballard 2023).As a result it is difficult to assert that short period exoplanets are synchronously rotating solely based on tidal models.
Ideally, one would thus constrain exoplanet rotation rates empirically.For hot Jupiters at least the bulk atmospheric rotation can be inferred via Doppler broadening of absorption lines or via phase curve measurements (Snellen et al. 2014;Kempton et al. 2014;Rauscher & Kempton 2014;Brogi et al. 2016).In practice, however, it is difficult to obtain precise constraints from such observations.For example, the Doppler broadening of HD 189733b is compatible with synchronous rotation, but to within 1σ the planet could also be in 3:2 or 2:1 rotation (Brogi et al. 2016).Similarly, the unusual westward hot spots found in the phase curves of hot Jupiters Kepler-7b, Kepler-12b, Kepler-41b, and CoRoT-2b might be due to non-synchronous rotation, but they could alternatively also be caused by clouds or magnetic drag (Rauscher & Kempton 2014;Hu et al. 2015;Shporer & Hu 2015;Dang et al. 2018).
In contrast to hot Jupiters, it should be relatively straightforward to infer the rotation of atmosphere-less exoplanets (bare rocks for short).The thermal phase curve of a bare rock primarily depends on the planet's rotation, its orbit, the surface's thermal inertia, and any potential tidal heating (Selsis et al. 2013).Many of these parameters are easier to constrain than the atmospheric circulation and cloud coverage of hot Jupiters, which means thermal phase curves of bare rocks are a promising way to test the tidal locking hypothesis.
Here we pursue this idea using observations of the super-Earth LHS 3844b.At a radius of 1.32 R and an orbital period of 11.1 hours, LHS 3844b is one of the smallest short period exoplanets with a measured thermal phase curve.The planet's Spitzer observations at 4.5 µm rule out many likely atmospheric scenarios and constrain the atmosphere's overall thickness to less than 1-10 bars (Kreidberg et al. 2019;Whittaker et al. 2022), compatible with the planet's transmission spectrum which rules out H 2 -rich and H 2 O-rich atmospheres in excess of 0.1 bars (Diamond-Lowe et al. 2020b).These observations are consistent with atmospheric escape models, which predict the planet lost its atmosphere due to the star's strong stellar wind and high XUV flux (Kreidberg et al. 2019;Diamond-Lowe et al. 2020a;Kane et al. 2020).The most plausible interpretation is thus that LHS 3844b is a bare rock, which removes any potential degeneracy between non-synchronous rotation and atmospheric heat redistribution, and makes LHS 3844b ideal for testing the tidal locking hypothesis.
Previous analyses of LHS 3844b found that, besides ruling out a thick atmosphere, its phase curve is consistent with a circular orbit, synchronous rotation, and a relatively low-albedo surface composed of basalt (Kreidberg et al. 2019;Whittaker et al. 2022).Nevertheless, it is unclear whether such a state is necessarily the best fit to the data.A priori one can expect that LHS 3844b should be in a tidally evolved state, such as synchronous rotation, as the planet should no longer retain any initial rotation or eccentricity.We estimate it would only take ∼ 1.7 × 10 10 sec = 530 years for LHS 3844b's orbit to circularize1 , with even less time required to spin down LHS 3844b's rotation.However, previous radial velocity measurements of LHS 3844 did not detect any periodic variations and only placed an upper bound of ∼ 1 M Jup on LHS 3844b's mass (Vanderspek et al. 2019), which means the system could host one or more companion planets that could keep LHS 3844b's eccentricity excited and allow the planet to still be in non-synchronous rotation.Similarly, previous work showed that LHS 3844b's secondary eclipse requires the surface to be quite dark, with basalt a likely candidate for the surface's composition (Kreidberg et al. 2019;Whittaker et al. 2022).However, previous conclusions about LHS 3844b's surface composition were based on geologically fresh materials and did not address the potential impact of space weathering, even though space weathering is known to strongly modify the properties of atmosphere-less bodies in the Solar System (Pieters & Noble 2016).The goal of this paper is thus to revisit the previous conclusions about LHS 3844b, and to explore what alternative scenarios might also fit the data.
The rest of this paper proceeds as follows.We develop a global thermal model for an atmosphere-less planet and compare its output to the Spitzer observations of LHS 3844b.The model is described in Section 2 and combines features of two previous exoplanet bare rock models (Hu et al. 2012;Selsis et al. 2013), including sub-surface heat diffusion, tidal heating, and bidirectional surface reflectance and emissivity.In contrast to previous studies of LHS 3844b, when comparing models against data we use the planet's full phase curve instead of just the derived secondary eclipse depth (see Section 2).Section 3 finds that the Spitzer data rule out Mercury-like spin-orbit resonances as well as eccentricities larger than e ≳ 0.001, largely because such spin and orbital states would cause intense tidal heating which is incompatible with the observations.In addition, we identify two scenarios which fit LHS 3844b's phase curve significantly better than a basalt planet on a circular orbit.The first requires that LHS 3844b still experiences weak tidal heating through a small-but-nonzero eccentricity (requiring a companion planet); the second requires that the surface of LHS 3844b has been substantially darkened by space weathering.Section 4 discusses why we cannot rule out tidal heating but consider space weathering more plausible, before Section 5 concludes the paper.

METHODS
Our goal is to simulate the thermal phase curve of LHS 3844b for arbitrary rotational and orbital states.As long as the planet has no atmosphere, or the atmosphere is optically thin and does not significantly redistribute heat, the phase curve of a bare rock planet will be dominated by only a handful of processes, which we model as follows.

Stellar flux
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 where k θ (t) represents the planet-star orbital angle, a the semi-major axis, and e the planet's eccentricity.The change in k θ (t) can be found from conservation of angular momentum, 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 where λ represents wavelength, R s the stellar radius, I s (λ) the stellar spectrum, A(λ, µ 0 ) is a specific form of surface albedo called the directional-hemispherical reflectance (described in more detail below), and µ 0 the zenith angle factor.The zenith angle factor is equal to µ 0 = cos(x − x 0 (t)) cos(y − y 0 (t)) on the dayside and µ 0 = 0 on the nightside, where x 0 and y 0 are the substellar latitude and longitude which are generally timedependent.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 lowmass 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).

Tidal dissipation
Planets with non-zero 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 where the tidal constants N (e), N a (e), and Ω(e) describe the dependence of tidal dissipation on eccentricity, given by (1−e 2 ) 9/2 .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, 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 default for LHS 3844b, but we additionally explore uncertainty in tidal heating by increasing and decreasing ∆t by one order of magnitude.If the tidal dissipation inside the planet is spatially uniform, then the surface heat flux is e tide = E tide /4πR 2 p .For planets with non-zero 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 pseudosynchronous state is ω eq /n = N (e)/Ω(e), which simplifies the tidal dissipation equation to 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 non-zero 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.

Thermal diffusion
Heat storage and heat redistribution inside the subsurface are modeled using the thermal diffusion equation, ∂T /∂t = k/(cρ)∇ 2 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, 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 latitudelongitude point are Here I 0 is the absorbed stellar flux and I p is the emitted thermal flux 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 Γ = √ kρc of about 50 to 80 J Winter & Saari 1969;Murdock 1974), and other data suggest Γ ∼ 100 J s −1/2 K −1 m −2 for regolith surfaces (Pierrehumbert 2010), so here we use Γ = 100 J s −1/2 K −1 m −2 .Tests show that our results are insensitive to the exact value of Γ (see Appendix A).

Bidirectional surface reflectance and emissivity
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 shadowhiding and back-scattering (Hapke 2002).
Here we use a modeling approach originally developed by Hapke for Solar System bodies like the Moon and Mercury (Hapke 1981(Hapke , 2002(Hapke , 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 which approximates the emergent flux of an illuminated semi-infinite regolith surface.Following Hapke (2012), the surface's bidirectional reflectance function is given by r(λ, µ 0 , µ e ) = w(λ) 4π Here w(λ) is the wavelength-dependent single-scattering albedo, µ 0 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 (Eqn.11).The relevant surface albedo in Equation 3 is the directional-hemispherical reflectance, which describes what fraction of a collimated light beam with incident angle µ 0 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 We compute the above quantities using singlescattering 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
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 exogeneous 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 mare 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) 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 where ϕ is the volume fraction of absorbing particles, and z is related to the refractive indices of host and absorbing material, defined as 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(Hapke , 2012)).The refractive indeces of graphite and iron nanoparticles are evaluated using data from the Refractive Index database (Polyanskiy 2022).

Comparison against observations
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 χ 2 .Doing so allows us to include all out-of-eclipse data points in the modeldata 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 χ 2 and χ 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 non-circular orbits one additionally needs to specify when the secondary eclipse occurs relative to periastron.We assume LHS 3844b's periastron coincides with 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 −3 (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 (Eqn.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 (Instrument 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 contrast = F p (λ) which is not the same as the alternative of taking the ratio first before spectrally averaging, F p (λ)/F s (λ) ̸ = F p (λ)/F s (λ).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.

Previous conclusions
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 χ 2 = 75.78(χ 2 ν = 3.3, p-value = 1.5 × 10 −7 ), whereas a blackbody would have a much better fit of χ 2 = 44.0(χ 2 ν = 1.9, p-value = 5 × 10 −3 ).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).
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 secondary eclipse is ∼50-100 ppm lower than the flux right after eclipse.However, the apparent dayside asymmetry is likely an artifact.First, previous work found that the data bin 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 Fig 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 Fig. 3).
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.

Non-zero eccentricity and non-synchronous rotation
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 non-zero 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 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), 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 Eqn. 5) by one order of magnitude upwards and downwards.
Figure 2 shows the planet's dayside equilibrium temperature T day for a planet that is either pseudosynchronized (blue curve) or in a Mercury-like 3:2 spinorbit resonance (green curve).The colored regions indicate the impact of varying ∆t by one 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 7,000 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).
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 one 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 (Meeus 1997;Peale et al. 1979).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 non-zero 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 χ 2 ν < 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 LHS3844b is 4 × 10 −4 < e < 10 −3 , 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 years.This means even if LHS 3844b was in pseudo-synchronous rotation, its rotation would be so slow as to be effectively tidally locked.
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 10700 W m −2 , significantly smaller than LHS 3844b's stellar constant of about 95600 W m −2 .The resulting fit, χ 2 = 32.74(χ 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, χ 2 = 75.78(see Table 1).Based on LHS 3844b's extremely short tidal circularization timescale of ∼ 500 years (see Section 1), a non-zero 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.

Space Weathering
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 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 non-linearly 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 non-linearity 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.
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 degrees 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 χ 2 = 75.78(see Table 1).Strong space weathering via the addition of 5% npFe improves the phase curve fit to χ 2 = 53.07(χ 2 ν = 2.4, p-value = 4 × 10 −4 ), while strong space weathering with 5% graphite produces an even better fit, χ 2 = 46.48(χ 2 ν = 2.0, pvalue = 3 × 10 −3 ).Note that the latter model with 5% graphite produces a phase curve that is almost identical to that of a blackbody (χ 2 = 46.48versus χ 2 = 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.

DISCUSSION
Our results show a pure basalt surface with zero tidal heating provides a poor fit to LHS 3844b's phase curve (χ 2 ν = 3.3), while the fit is significantly improved if we include either small tidal heating associated with a residual eccentricity (best model, χ 2 ν = 1.4), or space weathering (best model, χ 2 ν = 2.0).Of these two hypotheses, which one is more likely?
If we take the observational errors from Kreidberg et al. (2019) at face value then tidal heating provides a better fit to the data.One way to evaluate the different models in Table 1 is to impose a p-value cutoff.Considering only models with p-value above 0.05 as acceptable (χ 2 < 35.17), the only model that provides a reasonable fit to the data is the model with tidal heating and e = 0.0008.As an alternative, one can also use the minimum χ 2 method (Avni 1976).Treating the degree of tidal dissipation and space weathering in our thermal model as two independent parameters, one can reject with 0.99 significance level (roughly 3σ) all models that have ∆χ 2 > 9.21 relative to the model which minimizes χ 2 .Again, we find that in this case only models with tidal heating fall inside the 0.99 significance level relative to the best-fitting model with e = 0.0008.Based on LHS 3844b's short tidal circularization timescale, any nonzero tidal heating for the planet would require a third body that keeps its orbit excited, such as a planetary companion.Given the large uncertainties in the system's radial velocity data, a companion to LHS 3844b cannot be ruled out (Vanderspek et al. 2019).More precise ra-  dial velocity measurements of LHS 3844 would thus be highly desirable.
Nevertheless, we do not consider tidal heating necessary to explain LHS 3844b's phase curve.Tidal heating requires one to invoke a still-undiscovered planetary companion, whereas by analogy with the Solar System space weathering is highly likely to have occurred on LHS 3844b.As a pessimistic estimate we inflate the observational errors from Kreidberg et al. (2019) by 38.3%, so a blackbody would result in a phase curve fit of χ 2 ν = 1 (see Table 1, right columns).In this case a wide range of models with tidal heating and space weathering are able to match the data well, both based on their absolute pvalues and ∆χ 2 .Note that even with inflated error bars a pure basalt surface without tidal heating provides a poor fit to the data (p-value=0.02,∆χ 2 = 22.5, so > 3σ worse compared to the best-fitting model).We infer that LHS 3844b's phase curve indicates either tidal heating or space weathering.The first of these hypotheses requires strong assumptions about the existence of an additional planet, whereas the second hypothesis is physically plausible and can also match the data well if previous studies slightly underestimated the true observational error (also see Hansen et al. 2014).
Finally, to what extent are our conclusions affected by modeling uncertainty?The main conclusion, namely that LHS 3844b has to be in synchronous rotation (or nearly so) and has to have a small eccentricity, is sensitive to the assumed tidal model.In the constant time lag model uncertainty about the planet's internal structure is parameterized by the parameter ∆t (Leconte et al. 2010), but alternative modeling approaches exist.We tested the tidal parameterization from the VPLanet model (Barnes et al. 2020), and found tidal heat fluxes consistent with the constant time lag model to within a factor of two.Nevertheless, the planet's tidal heating could vary more strongly than anticipated.For example, tidal heating on Io is thought to generate a shallow subsurface magma layer (Khurana et al. 2011;Tyler et al. 2015).Depending on internal heat sources, LHS 3844b might similarly sustain an internal magma ocean which could strongly affect the planet's tidal response.Future work with more detailed tidal models is thus desirable to support our conclusions.

CONCLUSION
We have re-analyzed the Spitzer thermal phase curve of LHS 3844b using a global bare rock thermal model.We Table 1.Secondary eclipse depth and χ 2 phase curve fit for different surface compositions.All models match the previously-derived secondary eclipse depth of 380 ± 40 ppm to within 2σ (Kreidberg et al. 2019), but differ in their fits to the overall phase curve.Bold marks models that minimize χ 2 .Center columns use the original error bars from Kreidberg et al. (2019), while right columns inflate error bars by 1.383 to achieve χ 2 ν = 1 for a black body.P-values show the probability of obtaining a χ 2 value at least as large as the one observed, assuming the respective model is correct.find that the Spitzer data rule out non-synchronous rotation, such as a Mercury-like 3:2 spin-orbit resonance, because such rotation would generate tidal heating far in excess of the planet's observed flux.The lack of observed large tidal heating also constrains the planet's eccentricity, albeit with some uncertainty due to the planet's uncertain tidal response.As long as LHS 3844b's tidal dissipation efficiency lies within one order of magnitude of Earth's, its eccentricity has to be smaller than ∼ 0.001.This means that even if LHS 3844b was in pseudo-synchronous rotation, it has to rotate so slowly as to be effectively 1:1 tidally locked; more likely, given that rocky planets probably cannot be in pseudo-synchronous rotation (Makarov 2013), LHS 3844b has to be 1:1 tidally locked.
In addition the Spitzer data provide evidence for either tidal heating or space weathering on LHS 3844b.The secondary eclipse depth only provides a weak constraint, as all models we consider in detail match the secondary eclipse to within 2σ (Table 1).Instead we compare models directly to the observed phase curve and find that a geologically fresh basalt surface on a circular orbit provides a poor fit to the data.We therefore propose two hypotheses: either LHS 3844b still retains a small-butnonzero eccentricity which drives weak tidal heating, requiring a third body such as a companion planet, or LHS 3844b's surface has been darkened by space weathering.
Current data are insufficient to rule out tidal heating, but we consider space weathering more plausible.
Our results thus suggest that LHS 3844b is a potential exoplanet analog to the Moon and Mercury in our own Solar System, with a similarly darkened and spaceweathered surface.Future observations will be able to test and refine this interpretation in a number of ways.Radial velocity measurements can be used to constrain the planet's eccentricity and rule out the presence of planetary companions.Similarly, JWST observations can identify space weathering through its impact on the planet's secondary eclipse spectrum, while tidal heating can be inferred through more precise thermal phase curves.Such measurements should be possible not only for LHS 3844b, but also for other exoplanet bare rock candidates including GJ 1252b (Crossfield et al. 2022), TRAPPIST-1b (Greene et al. 2023), and TRAPPIST-1c (Zieba et al. 2023).
DK was supported by the National Natural Science Foundation of China (NSFC) under grant no.42250410318.RYH was supported by a grant from the STScI (JWST Program 1846) under NASA contract NAS5-03127.Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).Figure A5:.Phase curve of LHS3844b for different rotational states and for two different surface heat capacities.Left panel assumes a heat capacity of rock, Γ = 2420 J s −1/2 K −1 m −2 , right panel assumes a heat capacity of regolith Γ = 100 J s −1/2 K −1 m −2 .Tidal heating is disabled.All other parameters are the same as in Figure 1.

B. MODEL VALIDATION
To validate our bidirectional reflectance calculations, we compare them against a set of calculations performed with the model from Hu et al. (2012) for a hypothetical bare rock exoplanet.The assumed planet and star parameters correspond to a roughly Earth-sized around a cool M-type star, T eff = 2800K, M s = 0.1009M , R s = 0.1230R , R p = 0.977R .
Dashed black lines in Figure 6 show emission spectra computed with the numerical model of Hu et al. (2012), solid colored lines show our corresponding results.Both sets of calculations use the same stellar spectrum and single scattering albedo for a given surface material.We find that both sets of calculations agree closely on the variation of contrast ratio with respect to wavelength, and moreover agree on the relative flux difference between different surface types.In terms of model differences, we find that our calculations systematically produce a slightly lower flux than the calculations with Hu et al. (2012)'s model.The maximum flux difference at a single wavelength in the infrared occurs in the wiggles around 5 micron and amounts up to 10% flux difference, presumably due to the models' different spectral resolutions or spectral averaging.The maximum broadband-averaged flux difference is much smaller, about 2%.This level of model agreement is sufficient, given that it is much smaller than the typical precision of rocky planet emission spectra with JWST (Greene et al. 2023;Zieba et al. 2023).
To explain the remaining model differences, we hypothesize two plausible sources of error.First, we find that the relative offset between the two models is also robust at short wavelengths (less than 2µm), where the contrast ratio is dominated by reflected light.A model bias that is independent of wavelength suggests that the two models are using slightly different parameters or constants; for example, rounding off a fundamental constant like the Astronomical Unit, the Solar Radius, or Earth's Radius at the second digit would induce an O(1)% error, which in turn is squared and thus could result in up to 10% error in the calculated contrast ratio.Second, Figure 6 also shows some wavelengthdependent model differences.For example, for a Granitoid surface the two models essentially agree at 9 µm but slightly disagree at 14 µm, whereas for an Ultramafic surface the offset is roughly the same at both wavelengths.Unfortunately, we are not able to explain this source of error, as it could be caused by a large number of possible issues such as differences in horizontal (latitude-longitude) resolution, spectral resolution, or (in our model's case) biases in the surface temperature due to the model's finite vertical resolution.

Figure 1 .
Figure 1.Phase curve of LHS3844b.Dots show Spitzer data, lines are bare rock models which assume the planet is synchronously rotating and on a circular orbit.

Figure 2 .
Figure2.Observed dayside brightness temperature of LHS3844b (red; 1σ range), versus the planet's theoretical dayside equilibrium temperature for two different rotation states (solid curves).The blue curve assumes pseudo-synchronous rotation, the green curve assumes Mercury-like 3:2 rotation.Shaded blue and green regions indicate the effect of increasing or decreasing LHS 3844b's tidal dissipation efficiency by one order of magnitude.

Figure 3 .
Figure 3. Simulated phase curves as a function of eccentricity, assuming LHS 3844b is in pseudo-synchronous rotation and has a basalt surface.Blue curve shows zero eccentricity and is identical to the curve in Fig. 1.

Figure 4 .
Figure4.Effect of space weathering on the single scattering albedo of basalt (left column), and resulting thermal phase curves (right column).Top row shows the effect of space weathering via the formation of nanophase metallic iron particles (npFe), bottom row shows the effect of space weathering via graphite.Percentages show the added volume fraction of darkening material.
Figure B6:.Comparison of emission spectra for a hypothetical rocky exoplanet according to our model (colored lines) versus the model of Hu et al. (2012) (black lines).