© The Authors 2026

1 Introduction

Detecting and characterising Earth-like planets is a key step towards understanding planetary formation and evolution. Photometric surveys, such as the Transiting Exoplanet Survey Satellite (TESS) mission (Ricker et al. 2014), have identified hundreds of Earth-sized planets and candidates, many of which are subsequently characterised with ground-based radial velocities measurements (e.g. Hacker et al. 2024; Serrano Bell et al. 2024; Armstrong et al. 2025). Ultra-short-period planets (USPs) provide a unique opportunity for precise mass and radius measurements, as their RV semi-amplitudes are typically an order of magnitude larger than those of longer-period terrestrial planets (Sanchis-Ojeda et al. 2014). Additionally, USPs allow for the development of specific observational techniques to efficiently mitigate the effects of stellar rotational modulation (Hatzes et al. 2011). Accurately determining their fundamental parameters enables the measurement of bulk densities, which is crucial for constraining planetary compositions and internal structure. Understanding the compositions of these planets, in turn, informs models of planetary formation and evolution.

USPs are terrestrial worlds with orbital periods shorter than one day (Sahu et al. 2006; Sanchis-Ojeda et al. 2014; Winn et al. 2018). They are relatively rare, occurring around ~1% of stars (Sanchis-Ojeda et al. 2014; Petrovich et al. 2019), and are generally smaller than 2 R_⊕ (Rappaport et al. 2013). The USPs are typically found in multi-planet systems, favouring hosts with solar or sub-solar metallicity (Winn et al. 2017). Their architectures often show enhanced mutual inclinations compared to compact systems of longer periods, suggesting dynamical excitation accompanying orbital shrinkage (Dai et al. 2018). Ultra-short-period planets are also believed to have undergone significant migration, as their present-day orbits lie within the dust sublimation radius for solar-type hosts (Millholland & Spalding 2020), and recent statistical studies indicate that their occurrence increases with stellar age and that the architecture of USP systems evolves on gigayear timescales, consistent with inward tidal migration (Tu et al. 2025). Their bulk densities suggest compositions consistent with terrestrial planets, often resembling Earth-like or iron-rich interiors (Rappaport et al. 2013; Dai et al. 2019; Uzsoy et al. 2021). Given their extreme irradiation (>650 S_⊕ for Sun-like hosts), USPs are believed to be stripped rocky cores (Sanchis-Ojeda et al. 2014; Lundkvist et al. 2016), a hypothesis supported by transmission and emission spectroscopy (Crossfield et al. 2022; Kreidberg et al. 2019; Zhang et al. 2024). Precise measurements of mass and radius are particularly important in this context, as they allow for direct tests of atmospheric loss mechanisms and internal differentiation processes.

While most USPs lack significant volatile envelopes, some larger ones (e.g. TOI-561 b and 55 Cnc e) exhibit anomalously low densities that suggest the possible presence of residual atmospheres (Brinkman et al. 2023; Demory et al. 2016; Tsiaras et al. 2016; Dai et al. 2019), and recent observations indeed indicate a secondary atmosphere in 55 Cnc e (Hu et al. 2024). Precise composition constraints on these planets provide valuable insight into the mass dependence of planetary differentiation, particularly during the final stages of the giant impact phase (Scora et al. 2022). Furthermore, their short orbital periods make them ideal candidates for detailed surface characterisation with the James Webb Space Telescope (JWST), through the analysis of the phase light curve.

One of the earliest planetary discoveries from TESS was LHS 3844b, a highly irradiated terrestrial exoplanet orbiting the nearby mid-M dwarf LHS 3844 (Vanderspek et al. 2019). Located 15 parsecs away, this USP planet completes an orbit in approximately 11 hours and has a radius of 1.303 ± 0.022 R_⊕, placing it within the terrestrial planet regime. LHS 3844b is a key target for atmospheric and compositional studies because it orbits a small (0.189 ± 0.006 R_⊙) star nearby, making it a spectroscopically accessible terrestrial exoplanet.

The significance of LHS 3844b as a target for exoplanetary studies was significantly raised by multiple observational and theoretical efforts. Infrared phase curve measurements from Spitzer indicated that the planet does not have a thick atmosphere, down to a limit of 10 bars, consistent with efficient atmospheric erosion over its lifetime (Kreidberg et al. 2019). This conclusion was reinforced by Diamond-Lowe et al. (2020), who conducted ground-based optical transmission spectroscopy with Magellan/LDSS3C and found a featureless spectrum, disfavouring clear, low mean molecular weight atmospheres at surface pressures of 0.1 bars or greater. The role of high-energy stellar irradiation in driving atmospheric loss has been investigated, with estimates suggesting that the high XUV irradiation from LHS 3844 would have stripped a primordial hydrogen-helium atmosphere within 54 Myr (Diamond-Lowe et al. 2021). In contrast, Kane et al. (2020) explored an alternative scenario based on stellar age estimates (7.8 ± 1.6 Gyr) and geodynamical models, concluding that the volatile-poor nature of LHS 3844b may be the result of in situ formation within the system's snow line rather than post-formation atmospheric loss, although a giant impact stripping the planet's atmosphere is also possible. More recently, Kokori et al. (2023) refined the transit ephemeris of LHS 3844b, providing updated estimates of the orbital period and mid-transit time, P_t = 0.46292964 ± 0.00000033 d and t_0t = 2 458 828.93037 ± 0.00025 d (BJD), which we adopt in this work.

The planet has also been a prime target of the JWST. More than 40 hours were granted during Cycle 1 and 2 to study the planet mid-infrared spectrum with MIRI¹ and the planetary surface roughness with NIRSpec².

Despite the extensive observational efforts dedicated to LHS 3844b, a key piece of information remains unknown: its mass. A precise mass measurement is essential to constrain the planet's bulk composition and to determine whether it aligns with the terrestrial planets of the Solar System or represents a different class of rocky world. Moreover, independent determination of surface gravity is critical for interpreting atmospheric data gathered in transmission, as it directly influences atmospheric scale height and composition retrievals (e.g. Heng & Kitzmann 2017). Measuring the dynamical mass of LHS 3844b would therefore provide a fundamental constraint for both interior and atmospheric studies, refining our understanding of terrestrial planets around M dwarfs. In this work, we present the first mass determination of LHS 3844b, placing it in context with the growing population of well-characterised terrestrial exoplanets.

In Sect. 2 we describe the observational data and its reduction process. Section 3 deals with the characterisation of the central star of the system. In Sect. 4 we perform a preliminary exploration of the radial velocity data, including a GLS periodogram analysis and tests of additional linear covariates and mild orbital eccentricity. Section 5 describes our radial velocity modelling strategy, consisting of a Bayesian analysis of 15 alternative models combining different Gaussian process covariance structures and prescriptions for long-term velocity drifts. The results of the model comparison and the inferred planetary and orbital parameters are presented in Sect. 6, including an assessment of additional periodic signals and an interior characterisation of the planet. Finally, Sect. 7 summarises our main findings and conclusions.

2 Observations and data

2.1 Observations

LHS3844 was observed with the Échelle Spectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO), a fibre-fed echelle spectrograph installed at the European Southern Observatory Very Large Telescope (VLT) of the Paranal Observatory. A technical description of the instrument is provided by Mégevand et al. (2014); Pepe et al. (2014, and references therein), and its on-sky performance is evaluated by Pepe et al. (2021). The spectrograph works between 378.2 and 788.7 nm, and has a resolving power between 130 000 and above 190 000, depending on the observing mode, when the instrument is used in single-telescope mode³. Pepe et al. (2021) demonstrated that ESPRESSO provides precise radial velocity measurements at the level of 10 cm s⁻¹. Since its commission, ESPRESSO has been used in a large number of studies (e.g. Faria et al. 2020; Hobson et al. 2024; Suárez Mascareno et al. 2020, 2024; Ramírez et al. 2025; Cristiani et al. 2024).

LHS 3844 was targeted over two seasons: from November 2 2018 to December 3 2018 (14 visits; program 0102.C-0496, PI: Astudillo-Defru) and from June 9 2019 to September 6 2019 (11 visits; program 0103.C-0849, PI: Astudillo-Defru). In total, 25 spectra were acquired, with an exposure time of 2400 seconds for each in slow readout mode, and using 2x1 binning (mode HR21) of the instrument. In this mode, the spectrograph resolving power is around 138 000. The second fibre was used for sky subtraction.

In June 2019, the ESPRESSO fibre link was upgraded, probably introducing a radial velocity offset. We therefore treat the data before and after the upgrade separately, allowing for a velocity offset between them, as suggested by Pepe et al. (2021). There are 15 and ten measurements taken before and after the upgrade, respectively. The 15 ESPRESSO pre-upgrade spectra reached at 610 nm a median S/N of 16.1, with a minimum S/N of 11.8 and a maximum S/N of 18.7. The ten ESPRESSO post-upgrade spectra reached a median S/N of 19.6, with a minimum S/N of 15.4 and a maximum S/N of 24.3 at the same wavelength.

2.2 Data reduction and exploratory analysis

The ESPRESSO spectra were reduced using the ESO data pipeline version 3.0.0, described in Pepe et al. (2021). The extracted wavelength-calibrated two-dimensional spectra resulting from extraction were used to measure the radial velocities using a template-matching technique (Astudillo-Defru et al. 2015, 2017b). Due to the ESPRESSO intervention, we processed data from pre- and post-intervention as independent datasets. The S/N ratio reached translates into a median photon noise uncertainty of 60 and 52 cm s⁻¹, for the velocities before and after the upgrade operation, respectively, in line with the ESPRESSO performances described in Sect. 2.

In addition to the radial velocities measured with the template-matching technique, we computed the cross-correlation function (CCF) of the spectra with a numerical mask (Baranne et al. 1996; Pepe et al. 2002), from which we derived the full width at half maximum (FWHM) and contrast. We also measured activity indexes from the ESPRESSO spectra. The Ha index was computed following the prescriptions of Cincunegui et al. (2007) and Gomes da Silva et al. (2011). The H_β, H_γ, and Na D indices were computed using fluxes integrated within narrow wavelength bands centred on the respective spectral features, which are summarised in Table 1.

The resulting dataset is presented in Table C.2. None of the observations is expected to be contaminated by moonlight. In all observations, the Moon was further than 50 degrees from the target star. Also, no measurements are expected to be contaminated by sunlight either, because no observation was made during twilight. However, an exploratory data analysis revealed an unusually small FWHM for the observation on night BJD 2 458 643.91, the last visit before the fibre upgrade. For this observation, the contrast was not computed by the pipeline. This observation was performed closest to the morning twilight, but the Sun was still more than 21 degrees below the horizon. No contamination is expected and no trend in velocities is seen with Sun altitude, but the coincidence is noteworthy. A deeper study of the effect of sunlight on ESPRESSO observations is beyond the scope of this paper. The secular acceleration was computed from LHS 3844 proper motion and parallax (Gaia Collaboration 2021), yielding 0.219 m/s/yr, equivalent to 0.185 m/s over the time span of observations. The effect of secular acceleration was not removed from the data. The velocity dispersion is 5.4 and 5.8 m s⁻¹, for the pre- and post-upgrade sets, respectively. Prior to the modelling stage, the mean radial velocity of each instrumental sub-set (pre- and post-upgrade) was subtracted from the corresponding measurements to improve numerical stability. Since independent systemic velocity offsets are included in all models, this normalisation does not affect the inferred physical parameters.

3 Stellar characterisation

LHS 3844 is an M 4.5 or M 5 star (Vanderspek et al. 2019), located at a distance d = 14.89 ± 0.01 pc from the Sun (Gaia Collaboration 2018). Vanderspek et al. (2019) determined that the mass and radius of LHS 3844 are 0.151 ± 0.014 M_⊙ and 0.189 ± 0.006 R_⊙, based on empirical relations.

The rotational period of the star was measured to be P_rot = 128 ± 24 days by Vanderspek et al. (2019) using 1935 photometric measurements from the MEarth-South observatory. These data were obtained over more than two and a half years, but are mostly concentrated around the observing season between March and September 2018. This value agrees with the upper limit on the rotational rate given by the rotation-activity relation of Kiraga & Stepien (2007) using the X-ray upper limits of LHS 3844 from Diamond-Lowe et al. (2021).

The combined ESPRESSO spectra also allow for a determination of the stellar atmospheric parameters, following the procedure by Astudillo-Defru et al. (2017a). From the Ca II H&K analysis we obtain log($\left({R^{\prime }}_{\text{HK}}\right)$) = −5.534 ± 0.217. Although the low flux level at that spectral range translates into a large uncertainty, $\left({R^{\prime }}_{\text{HK}}\right)$ provides an independent estimation of 95 ± 32 days for the stellar rotation (Astudillo-Defru et al. 2017a), consistent within ~1σ with the photometric determination. Given the higher precision and long baseline of the MEarth data, we adopt P_rot = 128 ± 24 days as the reference value in this work.

4 RV preliminary analysis

We carried out a preliminary exploration of the ESPRESSO RV data independently of the transit information. This stage relied on a simplified linear modelling framework and aimed at identifying the dominant periodicities and assessing potential additional contributions to the signal. It consisted of a GLS periodogram analysis, statistical tests of additional linear covariates (stellar activity indicators and a temporal trend), and an evaluation of mild orbital eccentricity.

4.1 Linear modelling framework

Given an observed RV time series, {t_i, υ_i, σ_i}, for i = {1,…, N}, where σ_i is the velocity uncertainty based on photon noise statistics and calibration error, the velocity measurement at time t_i is modelled as $${\upsilon }_i=m_i={ϵ}_i,$$

where m_i is the prediction of the deterministic model at time t_i and ϵ_i is the error term. The variations in the RV measurements are caused by a combination of instrumental systematics, stellar activity, an additional physical acceleration, and/or the gravitational influence of orbiting companions.

As a basis for the exploratory analysis described in this section, we adopt a simplified linear model for the radial velocity time series. In this framework, we assume that the noise terms ϵ_i are not correlated and are normally distributed (white noise), the planet follows a circular orbit, and the orbital period is fixed. While not intended for the final parameter estimation, this approach provides a convenient and computationally efficient way to test for additional signals and assess the significance of potential model extensions.

The RV signal is modelled as $$m_i={\gamma }_0+{{\gamma }^{\prime }}_0{\delta }_{\text{post},i}+{\upsilon }_{\text{kep}}\left(t_i\right),$$(1)

where γ₀ is the global systemic velocity offset common to all observations, and ${{\gamma }^{\prime }}_0$ is the differential offset applied only to the post-upgrade ESPRESSO dataset. The indicator function δ_post,i equals 1 if the i-th observation belongs to the post-upgrade dataset and 0 otherwise. Under the assumption of a circular orbit with fixed period P, v_kep is given by $${\upsilon }_{\text{kep}}\left(t\right)={\beta }_0\text{\hspace{0.17em}sin}\left(\frac{2\pi t}{P}\right)+{\beta }_1\text{\hspace{0.17em}cos}\left(\frac{2\pi t}{P}\right),$$

from which the semi-amplitude K of the radial velocity variation can be expressed as $K^2={\beta }_0^2+{\beta }_1^2$.

The initial assumption of circularity is justified by the ultrashort orbital period of the planet, which leads to strong tidal interactions that damp orbital eccentricity and drive the system towards a nearly circular configuration (e.g. Hut 1980), but we lift this constraint below.

4.2 Periodogram analysis

In Fig. 1 we present the generalised Lomb Scargle (GLS) periodogram (Zechmeister & Kürster 2009) of velocities and stellar activity indicators. We have allowed for a frequency-dependent offset between the measurements acquired before and after the instrumental upgrade of June 2019 for each time series, as shown in Eq. (1).

In the RV periodogram (topmost panel) we detect a peak at P_max = 0.4628 ± 0.0007 d above the 0.1 per cent false alarm probability (FAP), which lies within 0.1σ of the transit-determined period. A second peak above the 0.1% FAP level appears at 0.86 d, corresponding to a first harmonic alias of the 0.46 d signal. After removing the best-fit sinusoidal model for the stellar RV variation at that period, the periodogram of the RV residuals (second panel) shows no further significant peaks. Within the scope of this simple model, we therefore find no evidence of an additional planetary signal in this analysis.

With the exception of the Na I D index (penultimate panel), no other activity index exhibits significant periodicity at the rotational rate of the star. In addition, no clearly significant peaks are seen around the period of the planetary companion, but the periodogram of the Hy index shows a peak with an amplitude at the level of 0.1%.

Fig. 1 GLS periodograms (Zechmeister & Kürster 2009) for the ESPRESSO radial velocities, their residuals after removing the best fit model for a circular orbit, and several stellar activity indicators measured for LHS 3844. Vertical dashed lines mark the transit period P = 0.46 days (green) and the rotational period of the star Prot = 128 days (red; Vanderspek et al. 2019). The 10, 1, and 0.1 percent false alarm probabilities (FAP) are shown as horizontal dashed lines.

4.3 Testing additional covariates

To further examine the effects of stellar activity on the RV data, we tested if the RV curves of our adopted models described in Sect. 5 should include additional linear terms associated with stellar activity indicators or a velocity trend over time. For this, we built simple linear models in which the radial velocity was modelled as a function of each activity indicator separately, in addition to the base components (instrumental offsets). Each extended model was compared with the corresponding base model using F-tests. In all cases, the resulting p-values exceeded 0.01. In addition, neither the Akaike nor the Bayesian information criterion (AIC, BIC, respectively) indicated any significant improvement.

We then incorporated a circular orbit signal with a fixed period corresponding to the highest peak in the GLS periodogram into each model and repeated the comparisons against the circular orbit-only model. Once again, none of the activity indicators yielded p-values below 0.01, and the AIC and BIC values similarly showed no evidence of improved fit. This shows that, at least within the domain of the tested linear model, it is not necessary to incorporate linear terms with these covariates to explain the RV data.

4.4 Testing mild orbital eccentricity

We also explored the possibility that the planetary orbit exhibits a small but non-zero eccentricity. To keep the model linear, we used the first order Taylor expansion in eccentricity of the Keplerian model, which led to the inclusion of two additional periodic terms at the first harmonic of the orbital period (see Appendix A.2). This approach allowed us to test for eccentricity effects with all the machinery of linear models.

We compared the extended mildly eccentric model (with the additional harmonic terms) to the simpler circular orbit model using F-tests and information criteria (AIC and BIC). The results did not show statistically significant improvement: all p-values were greater than 0.01, and both the AIC and the BIC values favoured the circular model. We therefore conclude that a small orbital eccentricity is not supported by the data and adopt the simpler circular model in the subsequent analysis.

5 Radial velocity modelling

We now move beyond the simplified linear framework adopted in the preliminary analysis and develop a comprehensive Bayesian model for the ESPRESSO radial velocities. This framework simultaneously accounts for correlated noise through Gaussian process regression, instrumental systematics, possible long-term velocity drifts, and informative priors on the orbital parameters. All model parameters and hyperparameters are treated as free and inferred jointly from the data.

5.1 Deterministic model

The deterministic component of the radial velocity signal consists of a Keplerian contribution and instrumental offsets, optionally supplemented by linear or quadratic velocity drifts. The component of the model not associated with the planetary signal is expressed as $${\gamma }_0+{{\gamma }^{\prime }}_0{\delta }_{\text{POST},i}+{\gamma }_1\left(t_i-t_{\text{ref}}\right)+{\gamma }_2{\left(t_i-t_{\text{ref}}\right)}^2,$$

where γ₁ and γ₂ correspond to the linear and quadratic drift terms, respectively. We adopted t_ref = 2458 578.66 d (BJD), chosen close to the midpoint of the observations to reduce the parameter correlations.

We define three nested versions of this component: a constant-only model (d₀), a model that includes a linear drift term (d₁), and a model that incorporates both linear and quadratic drift terms (d₂).

To incorporate prior constraints from transit photometry, we express the circular Keplerian signal in terms of the semiamplitude K, orbital period P, and phase φ $${\upsilon }_{\text{kep}}\left(t_i\right)=K\text{\hspace{0.17em}sin}\left(\frac{2\pi }{P}{t}_i+\varphi \right),$$

where, for a circular orbit, the phase relates to the transit mid-time as $\varphi =\pi -\frac{2\pi }{P}{t}_0$. This choice allows external information to be directly imposed through φ and P.

5.2 Error model and Gaussian process regression

The error terms ϵ = (ϵ₁,…, ϵ_N)^T account for both uncorrelated (white noise) and correlated noise sources. The uncorrelated noise component includes the measurement uncertainties σ_i and an additional jitter term σ_Jt for each instrument. Correlated noise is modelled using a GP with a covariance computed from a kernel function k(t_i, t_j), which here is stationary and depends exclusively on the time between observations |t_i – t_j|.

We explore the dependence of the inferred planetary mass with the following types of kernels: white noise (WN), squared-exponential (SE), periodic (Per), and quasi-periodic (QP) kernels. The kernel functions depend on the amplitude A, a characteristic timescale τ, a periodic timescale 𝒫, and a smoothness parameter λ, with functional forms given by $$\begin{array}{lll}{k}_{\text{WN}}\left(t_i,t_j\right)\hfill & =\hfill & A^2{\delta }_{ij},\hfill \\ k_{\text{SE}}\left(t_i,t_j\right)\hfill & =\hfill & A^2\text{\hspace{0.17em}exp}\left(-\frac{{\left(t_i-t_j\right)}^2}{2{\tau }^2}\right),\hfill \\ k_{\text{Per}}\left(t_i,t_j\right)\hfill & =\hfill & A^2\text{\hspace{0.17em}exp}\left(-\frac{2{\text{\hspace{0.17em}sin}}^2\left(\frac{\pi \left|t_i-t_j\right|}{𝒫}\right)}{{\lambda }^2}\right),\hfill \\ k_{\text{QP}}\left(t_i,t_j\right)\hfill & =\hfill & A^2\text{\hspace{0.17em}exp}\left(-\frac{{\left(t_i-t_j\right)}^2}{2{\tau }^2}\right)\text{exp}\left(-\frac{2{\text{\hspace{0.17em}sin}}^2\left(\frac{\pi \left|t_i-t_j\right|}{𝒫}\right)}{{\lambda }^2}\right).\hfill \end{array}$$

To characterise a broad range of temporal correlations present in RV time series, we test the following five kernel compositions, using the fact that the sum of valid kernel functions remains a valid kernel function.

1. WN

2. WN + SE

3. WN + QP

4. WN + SE + QP

5. WN + SE + Per

In all models we include both an instrumental jitter term for each version of ESPRESSO and a common white-noise term. The latter is intended to capture stellar variability that is equally present in both instruments but not adequately sampled by our data (e.g. granulation effects).

Finally, the covariance matrix for the error terms is given by $${\text{Σ}}_{ij}=k\left(t_i,t_j\right)+{\delta }_{ij}\left({\sigma }_i^2+{\sigma }_{Jt}^2\right),$$(2)

where k(t_i, t_j) is the covariance function constructed using one of the kernel combinations above.

5.3 Model priors and likelihood

We explored five alternative GP kernel prescriptions combined with three nested versions of the deterministic drift component (d₀, d₁, d₂), resulting in a total of 15 models. All parameters and hyperparameters were assigned the prior distributions summarised in Table 2, which were consistently applied across all models.

Informative priors were adopted for the orbital period P and transit mid-time t₀, based on the results of the transit analysis (Kokori et al. 2023). The reported transit epoch was propagated to the temporal midpoint of the radial velocity observations by shifting it by an integer number of orbital periods. The corresponding uncertainty was updated by propagating the period error accordingly. This choice ensures that the phase prior is referenced to the time span of the RV data, reducing parameter correlations and improving numerical stability. Through the phase relation defined in the previous subsection, the prior in t₀ is directly propagated into the radial velocity model. For kernels including a quasi-periodic component, the periodic timescale was assigned a prior informed by the stellar rotation period measured in Sect. 3.

Linear parameters $\beta =\left\{{\gamma }_0,{{\gamma }^{\prime }}_0,{\gamma }_1,{\gamma }_2,K\right\}$ were assigned multivariate normal priors, β ~ 𝒩(μ_β, B). Under this assumption, the marginal likelihood obtained after analytical integration over the linear parameters (see Rasmussen & Williams (2005), Sects. 2.7.1 and 5.4.1) can be written as $$\text{ln\hspace{0.17em}}p\left(v|\theta \right)=-\frac{1}{2}\left[{\left(m-v\right)}^{⊺}{K}_{\upsilon }^{-1}\left(m-v\right)+\text{ln}\left|K_{\upsilon }\right|+n\text{\hspace{0.17em}ln}\left(2\pi \right)\right],$$

where v is the vector of observed radial velocities, m = Χμ_β corresponds to the model evaluated at the prior mean of the linear parameters, n denotes the number of observations and K_υ = Σ + XBX^T, with X being the design matrix and Σ the covariance matrix given by Eq. (2). The vector of non-linear parameters θ includes the orbital period P, transit mid-time t₀, instrumental jitter terms, and the hyperparameters of the GP kernels. These parameters enter the likelihood through both the deterministic and stochastic components of the model.

5.4 Posterior distribution sampling

The posterior distributions of each model's parameters were explored using a Markov chain Monte Carlo (MCMC) algorithm, sampling only the non-linear parameter vector θ through the marginal likelihood defined in Sect. 5.3, where the linear parameters have been analytically integrated out. We used exoplanet to first obtain the maximum a posteriori (MAP) solution for each model by maximising the log-probability. The resulting parameter values were then used to initialise the PyMC3 sampler, which draws samples from the posterior using the No U-Turn Sampler (NUTS), a variant of Hamiltonian Monte Carlo. For each of the 15 models, we ran ten independent chains. We allowed 15 000 tuning steps, which were discarded, followed by 30 000 posterior samples per chain. Convergence was assessed using the rank-normalised split-$\widehat{R}$ statistic (Vehtari et al. 2021). For all models, the vast majority of parameters exhibit $\widehat{R}<1.01$, indicating satisfactory mixing and convergence.

For models that include a quadratic drift component (d₂), a small number (1–3) of chains were discarded due to persistent divergences during sampling. In these cases, the remaining chains showed stable behaviour. Additionally, in the d₂ models, the white-noise hyperparameters (instrumental jitter terms and the common white-noise amplitude) exhibited slightly elevated $\widehat{R}$ values, reaching up to 1.02. This behaviour is likely due to increased posterior correlations between curvature terms and noise amplitudes, which can introduce mild sampling inefficiencies without significantly affecting the inferred planetary parameters. All other parameters in these models satisfy $\widehat{R}<1.01$.

Subsequently, posterior samples of the linear parameters β were obtained by drawing from their conditional multivariate normal distribution, p(β | v, θ), evaluated in each posterior sample of θ. In this way, full posterior samples for all model parameters were recovered, while retaining the computational advantages of marginalising over the linear components during the MCMC exploration.

To compute the marginal likelihood of each model, we employed the estimator proposed by Perrakis et al. (2014), first implemented in an astronomical context by Diaz et al. (2016). The estimator relies on importance sampling using draws from the marginal posterior distribution of the non-linear parameters. The estimator relies on importance sampling using draws from the product of the marginal posterior distributions of the nonlinear parameters. To approximate independent samples from this distribution, the MCMC chains were randomly permuted independently for each parameter prior to evaluation, thereby breaking the correlations present in the joint posterior samples. For each model, the log-marginal likelihood was computed from sub-samples of increasing size (from 500 up to 60 000 draws), with each configuration repeated ten times using independent re-samplings. We adopt as the final estimate the mean value obtained with 60 000 samples. In all cases, the estimator showed clear convergence as the number of samples increased, with the dispersion decreasing monotonically, indicating the numerical stability of the evidence estimates.

Fig. 2 Posterior probability matrix for the 15 evaluated models, spanning all combinations of three secular acceleration models and five kernels. The numbers inside the cells indicate the normalised posterior probability of each model. The rightmost column and bottom row show the marginalised probabilities over drift and kernel types, respectively.

6 Results and discussion

6.1 Model comparison

Using the marginal likelihood estimates described in Sect. 5.4, we compare the full set of models explored in this work. Figure 2 presents the posterior probability of each model within the explored set of 15 kernel-drift combinations. These probabilities were obtained by normalising the marginal likelihoods under the assumption of equal prior probability for all models in the model space considered. The last row and column show the marginalised probabilities over the drift and kernel dimensions, respectively. Models without drift (d0) consistently dominate, regardless of the chosen kernel. This outcome is consistent with the preliminary linear modelling analysis (Sect. 4.3), where adding temporal drift terms did not produce statistically significant improvements based on F-tests or changes in AIC and BIC. Among kernel configurations, those including periodic or quasi-periodic components are favoured over simpler white-noise or squared-exponential kernels.

6.2 Inferred semi-amplitude and dominant models

Figure 3 shows the distribution of the RV semi-amplitude (K) as inferred from the full set of models explored in this work. Each violin plot represents the posterior distribution of K for a specific combination of the GP kernel and drift term, while the rightmost violin shows the marginal distribution across all models, weighted by their marginal likelihood. Despite substantial differences in the adopted noise prescriptions (from simple white noise to multi-component quasi-periodic kernels) and in the inclusion or exclusion of temporal drift terms, all models yield very similar posterior distributions for K. In practice, the inferred semi-amplitudes consistently fall in the range 6–7 m s⁻¹, indicating that the estimate is largely insensitive to the adopted noise model.

As an additional robustness check, we repeated the analysis using a maximum likelihood type-II (ML-II) approach, described in Appendix B. Although this method fixes the orbital period and epoch and optimises the marginal likelihood with respect to the hyperparameters rather than sampling them, the resulting estimates of K are fully consistent with those obtained from the full Bayesian MCMC framework. This agreement across fundamentally different inference schemes further supports the stability of the inferred semi-amplitude.

The posterior predictive distributions for the highest-probability model (WN+SE+QP: d₀) are shown in Fig. 4, which presents the phase-folded orbital component, and in Fig. 6 in Appendix C, which shows the predictive fits of GP to the residuals.

6.3 Inferred planetary properties of LHS 3844 b

Using the posterior samples described in Sect. 5.4, we derived the physical parameters of LHS 3844 b by marginalising over the full set of models through ancestral sampling, adopting the marginal likelihood as the weighting factor in the model averaging. The reported values correspond to the median and the 16th and 84th percentiles of the resulting marginalised posterior distributions. A summary of the planetary parameters is given in Table 3, including the estimates obtained in this work, values derived by combining our results with the transit analysis of Vanderspek et al. (2019) and Kokori et al. (2023), and those reported directly from the transit fit.

The orbital period is determined to be P = 0.4629296 ± 0.0000003 days, consistent with the prior based on transit observations, as expected given that the RV data alone carry less information on the planet's periodicity than the transit measurements. From the period we also derive the orbital semi-major axis a and the equilibrium temperatures for different albedo assumptions, which are necessarily compatible with the transit-based results. The short semi-major axis may suggest the possibility of star-planet interactions (Shkolnik et al. 2003, 2008). Although we tentatively detect variations in the Hy index in the orbital period of the planet (Sect. 4.2), such behaviour is not observed in other indicators, including those tracing regions higher in the stellar atmosphere. The available evidence therefore remains too limited to draw firm conclusions at this stage.

The RV semi-amplitude obtained in this work is $K={6.65}_{-0.53}^{+0.49}{\text{m\hspace{0.17em}s}}^{-1}$, yielding a minimum mass of M_p sin i; = 2.27 ± 0.23 M_⊕. The quoted uncertainty includes contributions from the errors in K and P (this work) and the stellar mass (from Vanderspek et al. 2019). This corresponds to a relative uncertainty of roughly 10%. Using the inclination measured from the transit geometry (i = 88.5°), the absolute mass is effectively indistinguishable from the minimum mass, as the difference is much smaller than the associated uncertainty (0.02 σ). Finally, the bulk density is found to be $\rho ={5.67}_{-0.61}^{+0.65}{\text{g\hspace{0.17em}cm}}^{-3}$ by combining our RV-based mass with the planetary radius derived from transit observations. This value indicates a composition consistent with a rocky super-Earth.

In Fig. 5 we place LHS 3844 b in the context of the known population of small exoplanets with well-constrained masses (σ_Μ/M < 0.3). As of February 2026, among the 6107 confirmed exoplanets, only 1072 (17.5%) have mass determinations with a relative uncertainty better than 10%, the precision we obtain here for LHS 3844 b. The planet lies close to the mass–radius relation expected for terrestrial iron–silicate compositions, but its peak posterior probability favours a slightly lower bulk density than an Earth-like composition, consistent with a rocky interior between the Earth's composition and a pure silicate (100% mantle) model. Taken together, our results confirm LHS 3844 b as an ultra-short-period, rocky exoplanet with well-constrained mass and density, likely of terrestrial composition.

Fig. 3 Comparison of the inferred RV semi-amplitude (K) across different models, ordered by increasing marginal likelihood. Each violin plot represents the posterior distribution of K for a given model, with the mean (solid line), median (dotted line), and 1st and 3rd quartiles (green). The rightmost violin corresponds to the weighted average model. The overlaid bar plot (grey) shows the relative log-likelihood difference for each model. The numerical annotations indicate the mean, standard deviation, and relative uncertainty of K for each model.

6.4 Linear and non-linear parameters

Summaries of the inferred posterior distributions for both nonlinear and linear model parameters are provided in Appendix C. Table C.3 reports the marginalised means and standard deviations of the non-linear hyperparameters grouped by kernel type, while Table C.4 presents the corresponding results for linear parameters, grouped by acceleration model. These tables offer a detailed overview of model dependencies and provide additional context for interpreting model selection based on marginal likelihood computation, as discussed above.

We find that the model requires distinct jitter terms for the pre- and post-upgrade datasets. The posterior distributions of ln σ_J (pre) cluster around −2.0 to −2.3 (i.e. ~0.1 m s⁻¹), whereas ln σ_J (post) spans ~ −0.2 to 0.7 (i.e. ~1–2 m s⁻¹), indicating a substantially larger excess variance in the post-upgrade data. This effect is most pronounced for kernels without a periodic component (WN and WN+SE), where the post-upgrade jitter dominates the variance budget. In contrast, the amplitudes of the time-correlated GP components (SE or QP/Per) are small compared to the post-upgrade jitter and typically comparable to or below the white-noise amplitude. As a result, the correlation hyperparameters (τ, λ, 𝒫_GP) remain weakly constrained and largely reflect their priors, consistent with the limited dataset.

Although periodic kernels introduce a mild preference for a characteristic timescale (see Fig. 6), the correlated component is weakly required. Overall, the variance is dominated by instrument-dependent jitter, and the RV signal is driven primarily by the planetary component, with no strong evidence for coherent stellar activity.

Regarding drift terms, the coefficients of the finear (γ₁) and quadratic (γ₂) trends show at most marginal deviations from zero (at the ≲2σ level), and their inclusion leads to a significant drop in marginal likelihood (see Fig. 2). This reinforces the conclusion from the linear models and the Bayesian model comparison that long-term trends are not supported by the data.

Fig. 4 Phase-folded RV curve of LHS 3844, based on the best-fit model with the highest marginal likelihood: WN + SE + QP Kernel with no drift (d0). The data points correspond to ESPRESSO observations obtained before (green circles) and after (red squares) the instrumental upgrade. The data were corrected by subtracting the per-instrument velocity offsets but not the noise component of the GP model. Phase folding was performed using the median of the posterior distribution for the orbital period. Error bars represent the quadrature sum of the reported uncertainties and the empirically inferred per-instrument jitter. The black curve shows the median orbital model (circular orbit), while the grey shaded region represents the 68% credible interval derived from posterior samples.

Fig. 5 Mass-radius diagram for small exoplanets (Mp < 5 Μ⊕) with precisely measured masses (σΜ/Μ < 0.3). Grey points represent the observed exoplanets with error bars, while the red point marks LHS 3844 b. The shaded background corresponds to a kernel density estimation (KDE) of the joint posterior distribution of mass and radius for LHS 3844 b. Solid lines indicate theoretical mass-radius relations for different compositions from Zeng et al. (2016). Known planets were sourced from the NASA Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu/) on February 26 2026.

6.5 Second periodic signal

To investigate the presence of additional periodic variability in the data, we computed a generalised Lomb–Scargle periodogram of the RV residuals obtained after subtracting the median posterior predictive deterministic component of the model with the highest marginal likelihood (WN+SE+QP kernel without drift terms; Sect. 5). The resulting periodogram (Fig. 7) reveals a peak at Ρ = 6.89 ± 0.15 days above the 0.1% FAP. This signal is not related to the stellar rotation period and no clear feature is seen at this period in the ESPRESSO window function (Fig. C.1). However, it could be related to the second harmonic of the alias produced by the orbital period of LHS 3844 b and the sidereal-day sampling, which would be expected at P_alias = 6.47 ± 0.30 days. This value lies within ~1.5σ of the detected peak. We also note that the same peak is present in the residuals of all other models considered in Sect. 5. The RV residuals after subtracting the dominant periodic component, folded in this period, are shown in Fig. 8. The periodogram also displays a secondary peak at Ρ = 14.85 ± 0.72 days, which corresponds to the first harmonic of the main signal and drops below the FAP thresholds once the 6.89-day periodic component is removed from the data.

Motivated by this result, we tested an extended model that included an additional periodic signal. We adopted the same Bayesian framework and modelling choices described in Sect. 5, starting from the configuration with highest marginal likelihood (WN+SE+QP - d₀) and adding a second Keplerian component with circular orbit. The orbital phase parameter was assigned a uniform prior, while the period was given a moderately informative prior, modelled as a normal distribution in ln Ρ centred on the period indicated by the residual periodogram (~6.886 days) with a dispersion of σ = 1 in log-period space.

The Markov chains converge well ($\widehat{R}<1.01$ for all parameters), indicating satisfactory mixing and convergence. The model recovers a second periodic signal with P₂ = 6.694 ± 0.007 days and semi-amplitude K₂ = 3.0 ± 0.4 m s⁻¹, corresponding to a minimum mass of M₂ sin i = 2.5 ± 0.4 M_⊕. In this configuration, the inferred semi-amplitude of LHS 3844 b increases to Κ = 7.85 ± 0.49 m s⁻¹. This would imply a larger planetary mass and would place LHS 3844 b more firmly along the Earth-like composition curve in the mass-radius diagram of Fig. 5. In contrast to the models discussed in Sect. 5, the noise budget becomes strongly dominated by the periodic component of the GP. The amplitudes of the jitter terms and of the white-noise and squared-exponential components are all consistent with zero within the uncertainties, while the quasi-periodic component remains significantly non-zero at the level of ~1−4 m s⁻¹. Its periodicity hyperparameter converges to ${𝒫}_{\text{QP}}={133}_{-19}^{+14}$ days, consistent with the prior expectation based on the stellar rotation period. The corresponding GP predictive model is shown in Fig. 9, where the correlated variability at the stellar rotation timescale becomes apparent once both periodic components are removed from the data.

The persistence of the ~6.9-day signal across the residuals of all tested models, together with its absence in the ESPRESSO window function and in the periodograms of the stellar activity indicators (Fig. 1), raises the possibility that it may correspond to an additional planetary companion in the system. If interpreted as such, the signal would be consistent with a non-transiting planet LHS 3844 c on a short-period orbit. This scenario is not implausible in the context of USP systems, whose architectures often show enhanced mutual inclinations relative to longer-period compact systems, potentially leading to non-transiting companions (Dai et al. 2018).

However, statistical support for this additional signal remains inconclusive. The marginal likelihood of the two-planet model is slightly higher than that of the single-planet models, but this comparison is affected by the informative prior adopted for the period of the second signal. Because this prior is centred on the detected peak, it likely increases the resulting model evidence. In addition, the detected period lies within ~1.5σ of the value expected for the second harmonic of an alias involving the orbital period of LHS 3844b and the sidereal-day sampling, further complicating its interpretation. For this reason, we do not include this model in the evidence-weighted averaging used to derive our final estimates of Κ and the other planetary parameters⁴. In addition, the current dataset comprises only 25 radial-velocity measurements, limiting our ability to robustly disentangle additional periodic signals from stellar variability and correlated noise. Additional RV observations will therefore be required to confirm whether the ~6.9-day signal corresponds to a second planet in the LHS 3844 system or arises from residual stellar or instrumental effects.

Fig. 6 RV residuals of LHS 3844 after subtracting the per-instrument offsets and the orbital component, sampled from the posterior distribution of the model with highest marginal likelihood: WN+SE+QP kernel without drift terms (d0). The blue curve represents the median of the GP predictive distribution, computed over the full posterior of the GP hyperparameters and marginalised over the linear parameters and orbital period. The shaded region denotes the 68% confidence interval. Data points correspond to ESPRESSO observations taken before (green circles) and after (red squares) the instrumental upgrade. Error bars include both the reported measurement uncertainties and the empirically derived per-instrument jitter, added in quadrature.

Fig. 7 Generalised Lomb–Scargle periodogram of the radial-velocity residuals after subtracting the median posterior predictive deterministic component (instrumental offsets and orbital signal) of the model with highest marginal likelihood (WN+SE+QP kernel without drift terms, d0). Horizontal dashed lines indicate the FAP levels of 10, 1, and 0.1%. The vertical red line marks the strongest peak in the periodogram at Ρ = 6.886 days.

6.6 Interior characterisation

We quantify the confidence regions of possible interior properties using the approach in Dorn et al. (2015, 2017) with recent updates from Dorn & Lichtenberg (2021); Luo et al. (2024). We use a surrogate-accelerated Bayesian inference framework (De Wringer et al. 2025), which replaces the computationally expensive physics-based forward model with a fast polynomial chaos-Kriging (PCK) surrogate directly within a Markov chain Monte Carlo (MCMC) sampling loop (Schobi et al. 2015; Marelli & Sudret 2014). For our inference, the surrogate model provides high quality fits with R-squared values (coefficient of determination) of 0.98 and 0.99 for the planetary mass and radius, respectively. Also, root mean square errors are well below the observational uncertainties 0.03 and 0.004 for the planetary mass and radius, respectively. Those errors of the model uncertainty are accounted for in the likelihood function.

The interior model assumes a rocky interior with the addition of water that is distributed between the different parts of a planet. The model consists of three layers: an iron-dominated core, a silicate mantle, and a H₂O steam atmosphere. We assume an adiabatic temperature profile for the core and mantle with possible temperature jumps at the core-mantle-boundary depending on melt temperatures. Further, we allow for both liquid and solid phases in those two layers. For liquid iron and iron alloys, we use the equation of state (EOS) by Luo et al. (2024). For solid iron, we use the EOS for hexagonal close- packed iron (Hakim et al. 2018; Miozzi et al. 2020). For pressures below -125 GPa, the solid mantle mineralogy is modelled using the thermodynamical model PERPLE_X (Connolly 2009) considering the system of MgO, SiO₂, and FeO. At higher pressures we define the stable minerals a priori and use their respective EOS from various sources (Hemley et al. 1992; Fischer et al. 2011; Faik et al. 2018; Musella et al. 2019). The liquid mantle is modelled as a mixture of Mg₂SiO₄, SiO₂, and FeO (Melosh 2007; Faik et al. 2018; Ichikawa & Tsuchiya 2020; Stewart et al. 2020), and mixed using the additive volume law.

Water can be added to the mantle and the core melts, depending on its solubility and partitioning behaviour, for which we follow (Dorn & Lichtenberg 2021; Luo et al. 2024). The addition of water reduces the density, for which we follow Bajgain et al. (2015) and decrease the melt density per wt% water by 0.036 g cm⁻³. For small water mass fractions, this reduction is nearly independent of pressure and temperature. The addition of water in core melts lowers the density as described in Luo et al. (2024). The effect of dissolved water on the melting temperature is accounted for. Beyond dissolved water in the deep interior, water can be in the solid, supercritical and gas phase, for which we employed the EOS compilation in Haldemann et al. (2020).

The transit radius of a planet is assumed to be at a pressure of P_Transit = 1 mbar. This is a simplification as the transit radius depends on temperature; however, the effect on the planet of interest is small. The thermal profile is assumed to be fully adiabatic, except for pressures less than the pressure at the tropopause (here fixed at 0.1 bar), where we keep an isothermal profile equal to the equilibrium temperature. We use an equilibrium temperature of 737 ± 25K, assuming an albedo of 0.3.

We first consider a purely rocky interior, where we fix an iron-free mantle and use a simple uniform prior (Table 4). Figure 10 shows the ID and 2D posterior distributions. In that case, the core mass fraction is well constrained to 0.29, which is similar to Earth's 0.325 estimate. There is a high posterior probability for core mass fractions lower than Earth's value (see Fig. 10), which is reflected in the lower mean density compared to an Earth-like composition. The possibility of an iron-rich super-Mercury interior is not supported by our new mass estimates (Kane et al. 2020). A purely rocky interior, as assumed here, with the absence of any thick atmosphere is disfavoured given atmospheric observations (Diamond-Lowe et al. 2020).

We also tested for a more general interior composition that allows for water distributed between core, mantle, and surface, as well as a variable mantle composition (see Appendix C.2). For this, we list the prior in Table C.1. This more general composition comes with more degeneracy, and we add the constraint that the bulk Fe/Si ratio is similar to Earth with 1.69 ± 0.35. We find that possible global water masses reach a few percents $\left({0.01}_{-0.01}^{+0.20}\%\right)$, although the interior is also consistent in this case with a dry interior. The core mass fraction has a slightly higher estimate of 0.3 in that case, compensating for the addition of water.

In summary, the mass-radius data support a purely rocky interior with a core mass fraction that is similar to or tendentially smaller than that of Earth. Mass-radius data alone cannot distinguish the existence or absence of a volatile layer on top of the rocky interior; however, such a volatile layer is disfavoured by phase curve measurements (Kreidberg et al. 2019), transmission spectroscopy (Diamond-Lowe et al. 2020), and atmospheric erosion considerations (Diamond-Lowe et al. 2021).

Fig. 8 RV residuals phase-folded at the period corresponding to the strongest peak in the residual periodogram. The blue curve shows the best-fitting sinusoidal model obtained from a weighted linear fit. Green circles and red squares correspond to ESPRESSO observations obtained before and after the instrumental upgrade, respectively. Error bars include the quadrature sum of the measurement uncertainties and the adopted instrumental jitter.

7 Conclusion

We have presented the first dynamical mass determination for the ultra-short-period planet LHS 3844 b, based on 25 high-resolution ESPRESSO spectra. Because the scope of the paper is to determine the mass of the planet, we focused our efforts on the RV data. The information of the transits is included as priors on model parameters.

To ensure a robust inference we adopted a fully Bayesian modelling framework and explored 15 competing radial-velocity models combining different GP covariance kernels and polynomial drift prescriptions. Marginal likelihoods were computed for all models and used for Bayesian model averaging and evidence-weighted parameter estimation. As a consistency check, we also derived analytical solutions using the ML-II approach, obtaining compatible mass estimates.

The planetary signal is robustly detected throughout all models. From our marginalised, evidence-weighted posterior samples we derive a minimum mass of M_p sin i = 2.27 ± 0.23 M_⊕ and, combining with the transit inclination and the adopted stellar mass, an absolute mass consistent with this minimum mass. The resulting bulk density, ρ = 5.67 ± 0.65 g cm⁻³, is consistent with a predominantly rocky composition.

Model comparison (marginal likelihoods) systematically favours kernels that include periodic or quasi-periodic components and disfavours models with additional long-term linear and/or quadratic drifts. Models without drift (d₀) dominate the posterior model probability across kernel choices, and three kernel configurations (WN+QP, WN+SE+QP, WN+SE+Per) contribute the bulk of the evidence. Importantly, the inferred Κ is stable across every model, demonstrating that the mass estimate is not sensitive to the precise choice of kernel or drift model.

The analysis of the RV residuals also reveals aperiodic signal at ~6.9 days that appears consistently across the residuals of all tested models. An extended two-planet model recovers this signal with a semi-amplitude of ~3 m s⁻¹, but its statistical support remains inconclusive given the limited number of observations and the possible impact of sampling and prior assumptions on inference. While this signal may hint at the presence of an additional non-transiting companion in the system, further RV measurements will be necessary to confirm its planetary nature.

Interior-structure inference indicates a predominantly rocky composition for LHS 3844 b with a core mass fraction similar to or slightly lower than Earth. Allowing for water-bearing interiors yields only trace water mass fractions within the posterior; a thick volatile envelope is disfavoured by the available constraints and by previous atmospheric observations (e.g. Kreidberg et al. 2019; Diamond-Lowe et al. 2020).

Taken together, the results robustly establish LHS 3844 b as a terrestrial, ultra-short-period planet with a well-constrained mass and density. Although the current RV dataset suffices to measure the planetary mass at the ~10% level, further improvements in precision and/or cadence (and additional epochs) would help to reduce the remaining uncertainty on K, better characterise the GP hyperparameters associated with stellar variability, and test for lower-amplitude companions. Finally, the precise mass and density reported here strengthen the scientific case for the continued atmospheric and surface characterisation of LHS 3844 b with JWST. Our dynamical constraints provide a necessary anchor for interpreting emission and phase-curve observations and for placing the planet in the broader context of terrestrial exoplanet composition and evolution.

During the preparation of this manuscript, we became aware of an independent and contemporaneous analysis by Nagel et al. (University of Göttingen), who also report a mass measurement of LHS 3844 b using ESPRESSO and CRIRES+ data.

Fig. 9 RV residuals of LHS 3844 after subtracting the per-instrument offsets and the two orbital components, sampled from the posterior distribution of the model described in Sect. 6.5. The blue curve represents the median of the GP predictive distribution, computed over the full posterior of the GP hyperparameters and marginalised over the linear parameters and orbital period. The shaded region denotes the 68% confidence interval. Data points correspond to ESPRESSO observations taken before (green circles) and after (red squares) the instrumental upgrade. Error bars include both the reported measurement uncertainties and the empirically derived per-instrument jitter, added in quadrature.

Fig. 10 Marginalised posterior distribution of interior parameters in ID and 2D. The shown posterior is calculated assuming a purely rocky interior with a fixed iron-free mantle composition and is constrained by planet mass, radius, and equilibrium temperature.

Acknowledgements

N.A-D. acknowledges the support of FONDECYT project 1240916. C.D acknowledges support from the Swiss National Science Foundation under grant TMSGI2_211313. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grant 51NF40_205606. X.B., X.D. & T.F. acknowledge funding from the French ANR under contract number ANR18-CE310019 (SPlaSH), and the French National Research Agency in the framework of Weighing the mass of LHS 3844 b the Investissements d'Avenir program (ANR-15-IDEX-02), through the funding of the "Origin of Life" project of the Grenoble-Alpes University.

References

Appendix A Orbital model

The effect of the planet's gravitational influence is modelled using a single Keplerian function, assuming one orbiting companion. The Keplerian signal is parameterised by the orbital period P, the semi-amplitude K, the eccentricity e, the argument of periastron ω, and the mean anomaly M₀. The RV contribution from the planet at time ti is $${\upsilon }_{\text{kep}}=K\left[\text{cos}\left(f_i+\omega \right)+e\text{\hspace{0.17em}cos}\left(\omega \right)\right],$$(A.1)

where f_i = f(t_i; P, e, ω, M₀) is the true anomaly at time t_i.

Appendix A.1 Circular orbit

Assuming a circular orbit, which corresponds to e = 0, the dependence on the argument of periastron ω is eliminated, and f_i becomes equal to M₀. Then, the radial velocity contribution simplifies to $${\upsilon }_{\text{kep}}=K\text{\hspace{0.17em}cos}\left(\omega \right)\text{cos}\left(M_0\right)+K\text{\hspace{0.17em}sin}\left(\omega \right)\text{sin}\left(M_0\right).$$

Since we can write the mean anomaly as $M_0=\frac{2\pi }{P}\left(t-\tau \right)$, with τ being the periastron passage time, we can rewrite the contribution from the circular orbit as $${\upsilon }_{\text{kep}}={\beta }_0\text{\hspace{0.17em}sin}\left(\frac{2\pi t}{P}\right)+{\beta }_1\text{\hspace{0.17em}cos}\left(\frac{2\pi t}{P}\right),$$

where ${\beta }_0\text{\hspace{0.17em}sin}\left(\omega -\frac{2\pi }{P}\tau \right)$ and ${\beta }_1=K\text{\hspace{0.17em}cos}\left(\omega -\frac{2\pi }{P}\tau \right)$. This way we can recover the amplitude K of the RV from the coefficients of the model, through $K=\sqrt{{\beta }_0^2+{\beta }_1^2}$.

Appendix A.2 Mildly eccentric orbit

We can lift the assumption of circular orbit while keeping the linear nature of the model and all of its advantages by expanding the Keplerian expression of the radial velocity (Eq. A.1) as a series of the orbital eccentricity. To first order, we obtain (Lucy & Sweeney 1971) $${\upsilon }_{\text{kep}}=K\text{\hspace{0.17em}cos}\left(M_0\right)+eK\text{\hspace{0.17em}cos}\left(2M_0-\omega \right)+O\left(e^2\right).$$(A.2)

Expanding the second term and neglecting terms of order O(e²) or higher, this simplifies to $$V=V_0+K^{\prime }\text{\hspace{0.17em}cos}\left(M_0\right)+e{K}^{\prime }\text{\hspace{0.17em}cos}\left(\omega \right)\text{sin}\left(2M_0\right)+e{K}^{\prime }\text{\hspace{0.17em}sin}\left(\omega \right)\text{cos}\left(2M_0\right),$$(A.3)

where K′ is related to K by $$K^{\prime }=\frac{\text{cos}\left(\omega -\frac{2\pi }{P}\tau \right)}{\text{cos}\left(\frac{2\pi }{P}\tau \right)}K.$$(A.4)

Since 2M₀(P) = M₀(P/2), the eccentricity-dependent terms introduce periodic components at the first harmonic P/2 instead of the orbital period P. Thus, by fitting these periodic terms to the data, one can estimate their coefficients and infer the orbital eccentricity.

Appendix B Maximum likelihood type-II robustness analysis

As an additional robustness check of the planetary signal characterized in Sect. 5, we implemented an alternative inference scheme based on the Maximum Likelihood Type-II approach (ML-II; Rasmussen & Williams 2005). In this framework, the likelihood marginalised over the linear parameters is maximised with respect to the hyperparameters, which are fixed at the values that maximise the marginal likelihood of the model.

We adopted the same deterministic and noise models described in Sect. 5, assuming Gaussian priors for the linear parameters. However, in contrast to the full Bayesian MCMC analysis, the planetary period P and transit epoch t₀ were fixed to their transit-derived values. Under these assumptions, and given that the measurement errors are modelled as multivariate normal, the marginalisation over the linear parameters can be performed analytically.

For each of the five Gaussian process kernels considered in the main analysis, we applied an ML-II optimisation procedure. First, the marginal likelihood was maximised with respect to the GP hyperparameters of the selected kernel. Once the optimal hyperparameters were obtained, the model was conditioned on these values and a Bayesian linear regression was performed to derive the posterior distribution of the linear parameters, including the planetary semi-amplitude K and instrumental offsets. The semi-amplitude K values and the corresponding log-marginal-likelihoods for each kernel model are reported in Table B.1.

To facilitate comparison with the fully Bayesian results, sampling distributions for K were generated for each ML-II model by drawing from normal distributions centred on the analytical estimates and with standard deviations given by their corresponding uncertainties. The models were ordered according to increasing marginal likelihood, and a model-averaged estimate was constructed by weighting each model by its marginal likelihood. The resulting distributions are shown in Fig. B.1. The inferred values of K are fully consistent with those obtained from the MCMC analysis presented in the main text, confirming the robustness of the detected planetary signal under this alternative inference framework.

Appendix C Additional plots and tables

The ESPRESSO data described in Sect. 2 are listed in Table C.2, and their window function is shown in Fig. C.1. Further details on the radial velocity modelling are provided in Tables C.3 and C.4, which summarise the posterior distributions of the non-linear and linear model parameters, respectively. Finally, Table C.1 lists the prior parameter distributions adopted in the interior modelling (described in Sect. 6.6), while Fig. C.2 shows the resulting marginalised posterior distributions.

Fig. B.1 Comparison of the inferred RV semi-amplitude (K) across the ML-II analytical models, ordered by increasing marginal likelihood. For each model, the violin represents a sampling distribution obtained by drawing from a normal distribution with mean K and standard deviation σK derived from the analytical ML-II results. The mean (solid line), median (dotted line), and 1st/3rd quartiles are indicated. The second violin on the right corresponds to the ML-II model-averaged estimate, constructed by weighting each model according to its marginal likelihood, while the rightmost violin shows the posterior model-averaged estimate derived from the MCMC analysis. Numerical annotations report the mean, standard deviation, and relative uncertainty of K for each case.

Fig. C.1 Window function of the ESPRESSO data, dominated by a strong peak at the sidereal day.

Fig. C.2 Marginalised posterior distribution of interior parameters in 1D and 2D. The shown posterior is calculated assuming a rocky interior with the possibility of water distributed between core, mantle, and surface. The mantle composition can vary in terms of Fe/Si and Mg/Si ratios. The posterior is constrained by planet mass, radius and equilibrium temperature, and an Fe/Si bulk ratios similar to Earth (see main text).

All Tables

All Figures

Fig. 1 GLS periodograms (Zechmeister & Kürster 2009) for the ESPRESSO radial velocities, their residuals after removing the best fit model for a circular orbit, and several stellar activity indicators measured for LHS 3844. Vertical dashed lines mark the transit period P = 0.46 days (green) and the rotational period of the star Prot = 128 days (red; Vanderspek et al. 2019). The 10, 1, and 0.1 percent false alarm probabilities (FAP) are shown as horizontal dashed lines.

In the text

	Fig. 2 Posterior probability matrix for the 15 evaluated models, spanning all combinations of three secular acceleration models and five kernels. The numbers inside the cells indicate the normalised posterior probability of each model. The rightmost column and bottom row show the marginalised probabilities over drift and kernel types, respectively.
In the text

Fig. 3 Comparison of the inferred RV semi-amplitude (K) across different models, ordered by increasing marginal likelihood. Each violin plot represents the posterior distribution of K for a given model, with the mean (solid line), median (dotted line), and 1st and 3rd quartiles (green). The rightmost violin corresponds to the weighted average model. The overlaid bar plot (grey) shows the relative log-likelihood difference for each model. The numerical annotations indicate the mean, standard deviation, and relative uncertainty of K for each model.

In the text

Fig. 4 Phase-folded RV curve of LHS 3844, based on the best-fit model with the highest marginal likelihood: WN + SE + QP Kernel with no drift (d0). The data points correspond to ESPRESSO observations obtained before (green circles) and after (red squares) the instrumental upgrade. The data were corrected by subtracting the per-instrument velocity offsets but not the noise component of the GP model. Phase folding was performed using the median of the posterior distribution for the orbital period. Error bars represent the quadrature sum of the reported uncertainties and the empirically inferred per-instrument jitter. The black curve shows the median orbital model (circular orbit), while the grey shaded region represents the 68% credible interval derived from posterior samples.

In the text

Fig. 5 Mass-radius diagram for small exoplanets (Mp < 5 Μ⊕) with precisely measured masses (σΜ/Μ < 0.3). Grey points represent the observed exoplanets with error bars, while the red point marks LHS 3844 b. The shaded background corresponds to a kernel density estimation (KDE) of the joint posterior distribution of mass and radius for LHS 3844 b. Solid lines indicate theoretical mass-radius relations for different compositions from Zeng et al. (2016). Known planets were sourced from the NASA Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu/) on February 26 2026.

In the text

Fig. 6 RV residuals of LHS 3844 after subtracting the per-instrument offsets and the orbital component, sampled from the posterior distribution of the model with highest marginal likelihood: WN+SE+QP kernel without drift terms (d0). The blue curve represents the median of the GP predictive distribution, computed over the full posterior of the GP hyperparameters and marginalised over the linear parameters and orbital period. The shaded region denotes the 68% confidence interval. Data points correspond to ESPRESSO observations taken before (green circles) and after (red squares) the instrumental upgrade. Error bars include both the reported measurement uncertainties and the empirically derived per-instrument jitter, added in quadrature.

In the text

Fig. 7 Generalised Lomb–Scargle periodogram of the radial-velocity residuals after subtracting the median posterior predictive deterministic component (instrumental offsets and orbital signal) of the model with highest marginal likelihood (WN+SE+QP kernel without drift terms, d0). Horizontal dashed lines indicate the FAP levels of 10, 1, and 0.1%. The vertical red line marks the strongest peak in the periodogram at Ρ = 6.886 days.

In the text

Fig. 8 RV residuals phase-folded at the period corresponding to the strongest peak in the residual periodogram. The blue curve shows the best-fitting sinusoidal model obtained from a weighted linear fit. Green circles and red squares correspond to ESPRESSO observations obtained before and after the instrumental upgrade, respectively. Error bars include the quadrature sum of the measurement uncertainties and the adopted instrumental jitter.

In the text

Fig. 9 RV residuals of LHS 3844 after subtracting the per-instrument offsets and the two orbital components, sampled from the posterior distribution of the model described in Sect. 6.5. The blue curve represents the median of the GP predictive distribution, computed over the full posterior of the GP hyperparameters and marginalised over the linear parameters and orbital period. The shaded region denotes the 68% confidence interval. Data points correspond to ESPRESSO observations taken before (green circles) and after (red squares) the instrumental upgrade. Error bars include both the reported measurement uncertainties and the empirically derived per-instrument jitter, added in quadrature.

In the text

	Fig. 10 Marginalised posterior distribution of interior parameters in ID and 2D. The shown posterior is calculated assuming a purely rocky interior with a fixed iron-free mantle composition and is constrained by planet mass, radius, and equilibrium temperature.
In the text

Fig. B.1 Comparison of the inferred RV semi-amplitude (K) across the ML-II analytical models, ordered by increasing marginal likelihood. For each model, the violin represents a sampling distribution obtained by drawing from a normal distribution with mean K and standard deviation σK derived from the analytical ML-II results. The mean (solid line), median (dotted line), and 1st/3rd quartiles are indicated. The second violin on the right corresponds to the ML-II model-averaged estimate, constructed by weighting each model according to its marginal likelihood, while the rightmost violin shows the posterior model-averaged estimate derived from the MCMC analysis. Numerical annotations report the mean, standard deviation, and relative uncertainty of K for each case.

In the text

	Fig. C.1 Window function of the ESPRESSO data, dominated by a strong peak at the sidereal day.
In the text

Fig. C.2 Marginalised posterior distribution of interior parameters in 1D and 2D. The shown posterior is calculated assuming a rocky interior with the possibility of water distributed between core, mantle, and surface. The mantle composition can vary in terms of Fe/Si and Mg/Si ratios. The posterior is constrained by planet mass, radius and equilibrium temperature, and an Fe/Si bulk ratios similar to Earth (see main text).

In the text