© The Authors 2026

1 Introduction

The HD 332231 system was first flagged as a promising candidate when the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) recorded a single transit with a depth of approximately 0.5% (∼ 5000 ppm) in stellar flux during its two-year Prime Mission (Guerrero et al. 2021). Subsequent high-precision radial-velocity (RV) observations obtained through the TESS–Keck Survey (Chontos et al. 2022) definitively confirmed the planetary nature of the signal (Dalba et al. 2020). These measurements revealed a companion with a mass of 80 M_⊕ and radius of 9.5 R_⊕ on a nearly circular orbit with a period of 18.7 days. Designated HD 332231 b, the planet has a mean density of 0.5 g cm⁻³, placing it firmly in the category of warm, sub-Saturnian worlds.

The host star is a 3 Gyr-old dwarf with a mass, M_*, of 1.2 M_⊙; radius, R_*, of 1.3 R_⊙; and effective temperature of 6130 K (MacDougall et al. 2023). The study by Knudstrup & Albrecht (2022) shows that the orbital angular momentum of HD 332231 b is aligned with the stellar spin. The orbital eccentricity was found to be consistent with zero within 2 σ (Dalba et al. 2020; Knudstrup & Albrecht 2022; Polanski et al. 2024), and it must be primordial because the tidal circularisation timescale is of the order of several hundreds of gigayears, which is much longer than the age of the system.

This orbital configuration of HD 332231 b supports quiescent formation scenarios, in which the planet was formed in situ or ex situ – beyond the water ice line – and smoothly migrated inward while embedded in a protoplanetary disc (see Dawson & Johnson 2018 for a detailed review). In such scenarios, any lower mass planet, if formed together with HD 332231 b, might still be preserved. This finding motivated us to select the HD 332231 system for a project aimed at searching for additional planets accompanying gas giants in wide and circular orbits.

Furthermore, Knudstrup & Albrecht (2022) reported possible variations in observed transit times, and Polanski et al. (2024) flagged HD 332231 b as exhibiting significant transit timing variations (TTVs). Although these irregularities in orbital motion remained un-characterised due to a limited number of observations, gravitational perturbations from an undetected planetary companion were invoked as a possible explanation (Knudstrup & Albrecht 2022). These findings further boosted our motivation to revisit the system properties using TESS observations obtained over the last five years.

2 TESS observations and data reduction

The apparent magnitude of HD 332231 in the TESS passband is 8.0417 ± 0.0060 (Stassun et al. 2019). This brightness ensures sub-millimagnitude precision on one-minute timescales in the photometric time series. The system was observed in Sectors 14, 15, 41, 55, 75, 81, and 82 in short-cadence mode (120 second exposures). For Sectors 75, 81, and 82, ultra-short-cadence data with 20-second resolution are also available. However, these were not incorporated into the final analysis as their inclusion does not influence the results, and a homogeneous data-reduction procedure was preferred in all sectors. Details on individual observing runs are summarised in Table A.1.

Fluxes were extracted from target pixel files using custom scripts based on the Lightkurve package (version 2.5.0; Lightkurve Collaboration 2018). Raw light curves, shown in Fig. A.1, were produced via aperture photometry, employing the optimal masks provided within the target pixel files. Scattered light was mitigated using the RegressionCorrector tool, which outperformed the standard sky-background subtraction based on median filtering. This approach reduced the out-of-transit (OOT) scatter by 2−3% in Sectors 41−82 and by 13% and 16% in Sectors 14 and 15, respectively. To ensure that the final analysis relied solely on high-quality data, we discarded measurements with cadence quality flags other than zero.

The HD 332231 b transits were initially masked using the ephemeris of Knudstrup & Albrecht (2022) and subsequently manually refined to account for variations in the transit timing. Long-term trends and low-frequency systematics were removed using the flatten function, which applies a Savitzky–Golay filter (Savitzky & Golay 1964) with a window width of 12 hours (361 data points). A 5 σ clipping was then applied to the OOT data to eliminate outliers. Finally, the light curves were corrected for flux contamination from nearby stars using the CROWDSAP values provided by the pipeline for each sector.

3 Modelling of TESS transit data

Given the 18.7-day orbital period of HD 332231 b, typically only a single transit falls within each TESS sector. Exceptions include Sectors 41 and 81, where two transits were recorded. Although the second event in Sector 81 was only partially observed, our subsequent analysis confirmed its usefulness. In contrast, the transit in Sector 14 was affected by scattered light, distorting its shape, and it was excluded from the final analysis. Further discussion is provided in Appendix B.

We used the Transit Analysis Package (TAP; Gazak et al. 2012) to derive the best-fitting transit model. TAP was applied to data segments centred on the mid-transit times. Each segment was five times the transit duration in length, ensuring a sufficient OOT baseline to accurately characterise photometric noise. This process was carried out in two iterations: the first used the initial ephemeris, and the second employed an updated ephemeris that includes a periodic term (see Sect. 4), enabling improved centring of the data chunks.

TAP generates analytic transit models following the formalism by Mandel & Agol (2002). The main model parameters describing transit geometry include the orbital inclination (i_b), the scaled semi-major axis (a_b/R_*), and the planet-to-star radius ratio (R_b/R_*). Stellar limb darkening was modelled using a quadratic law (Kopal 1950), with linear (u₁) and non-linear (u₂) coefficients. These five parameters were jointly fitted using all available light curves. For each dataset, TAP also independently fitted the mid-transit time (T_mid), uncorrelated and correlated photometric noise, and any remnant quadratic trends in the flux baseline. The orbital eccentricity (e_b=0.029 ± 0.019) and argument of periastron (ω_b=39^° ± 31^°) were incorporated as Gaussian priors taken from Knudstrup & Albrecht (2022, the shadow model in their Table 2).

We employed ten Markov chain Monte Carlo chains, each with 10⁶ steps, to explore the posterior distributions of all free parameters. Because the posteriors of some parameters were significantly skewed, we adopted the mode (rather than the median) of the histograms (constructed following the Freedman–Diaconis rule) as the best-fit value. Uncertainties were estimated from the 68% highest density intervals surrounding the mode.

We also analysed posterior distributions for the transit impact parameter, $$b_{\mathrm{b}}=\frac{a_{\mathrm{b}}}{R_{\star }}\frac{1-e_{\mathrm{b}}^2}{1+e_{\mathrm{b}}sin{\omega }_{\mathrm{b}}}cos{i}_{\mathrm{b}},$$(1) and the total transit duration, $${\tau }_{14,\mathrm{b}}=\frac{P_{\mathrm{b}}}{\pi }\frac{\sqrt{1-e_{\mathrm{b}}^2}}{1+e_{\mathrm{b}}sin{\omega }_{\mathrm{b}}}\sqrt{{(1+\frac{R_{\mathrm{b}}}{R_{\star }})}^2-b_{\mathrm{b}}^2},$$(2) where P_b is the orbital period. Both quantities depend on i_b and a_b/R_*, which are correlated (see Fig. C.1). As discussed in Appendix D, a test run in which transit parameters were allowed to vary independently across light curves revealed no statistically significant variations in these parameters.

Individual light curves and their best-fitting models are shown in Fig. 1. The inferred transit parameters are reported in Table 1.

Fig. 1 Transit light curves of HD 332231 b observed with TESS in Sectors 15–82. Left: Individual photometric time series sorted by the epoch of observation, with numbering consistent with the refined ephemeris given in Sect. 4. Best-fitting models are overlaid in red. Right: Corresponding photometric residuals from the transit models.

4 Refining the transit ephemeris for HD 332231 b

The TESS observations yielded mid-transit times for eight epochs, spanning 97 orbital cycles, that is, equivalent to 1814 days, or nearly five years. A trial fit using a constant-period ephemeris in the classic form of $$T_{\mathrm{mid},\mathrm{b}}\left(E\right)=T_{0,\mathrm{b}}+P_{\mathrm{b}}·E,$$(3) where E is the transit (cycle) number counted from the reference epoch T_{0, b} (set to the mid-point of the transit observed in Sector 15), resulted in a reduced chi-squared (${\chi }_{\text{red}}^2$) value of approximately 440, revealing the residuals of the order of tens of minutes (see Fig. E.1a). A periodogram analysis using the analysis of variance (AoV) algorithm (Schwarzenberg-Czerny 1996) revealed a signal with a period of approximately 100 orbital cycles (Fig. E.1b).

Knudstrup & Albrecht (2022) previously used four midtransit times from Sectors 14, 15, and 41 to refine the linear ephemeris. They reported that their spectroscopic transit observation, taken between Sectors 15 and 41, aligned with the model only when the mid-point was allowed to shift by around +20 minutes relative to the linear prediction. However, their transit timing dataset was too sparse to permit further interpretation. We incorporated this mid-transit time into our dataset, adopting the value derived from the planetary shadow model (see Table 2 of Knudstrup & Albrecht 2022). As described in Appendix F, we also attempted to determine a mid-transit time from the RV time series provided by Sedaghati et al. (2022). Our reanalysis yielded an inconclusive result, and we therefore excluded it from further consideration.

To account for the timing variations, we expanded the ephemeris model by including a periodic term: $$T_{\mathrm{mid},\mathrm{b}}\left(E\right)=T_{0,\mathrm{b}}+P_{\mathrm{b}}·E+A_{\mathrm{TTV}}·sin\left[2\pi (\frac{E}{P_{\mathrm{TTV}}}-{\phi }_{\mathrm{TTV}})\right],$$(4) where A_TTV, P_TTV, and ϕ_TTV denote the amplitude, period, and phase of the TTV signal, respectively. These parameters were extracted from posterior distributions built using 100 Markov chains, each with 10⁴ steps (excluding the first 1000 as burn-in). The median values and uncertainties (15.9, 50, and 84.1 percentiles) are presented in Table 1. The best-fitting model yields a ${\chi }_{\text{red}}^2$ of 0.39, suggesting either slight over-fitting or that our timing uncertainties may be overestimated. In either case, the residuals leave no room for a more complex model at this stage.

Figure 2 presents the measured TTV signal, with the linear term of the ephemeris removed in upper panel, and the residuals after fitting the full model in lower panel. The best-fitting periodic model is plotted, along with 100 randomly drawn solutions from the posterior distribution to illustrate model’s uncertainty. HD 332231 b exhibits TTVs with an amplitude of $46.3_{-5.7}^{+7.4}$ minutes – over an order of magnitude larger than the uncertainties in the individual T_mid,b measurements – and with a period of $131.0_{-7.6}^{+8.9}$ orbital cycles $(={2450}_{-140}^{+170}$ days $=6.7_{-0.4}^{+0.5}$ years).

Further follow-up observations are expected to improve the precision of the TTV parameters, particularly the period P_TTV. The grey bands in Fig. 2 mark a transit window in 2026, observable from the ground. As the system is not scheduled for additional TESS coverage before August 2027¹, groundbased transit monitoring could be valuable. While scheduling such observations is challenging due to the long P_b and τ_{14, b} of approximately 6 hours, even partial transits – centred near ingress or egress – can help constrain the ephemeris more tightly.

Fig. 2 TTVs of HD 332231 b. Top: residuals after subtracting the linear ephemeris. The sinusoidal TTV model and 100 posterior models are shown. The grey shaded region marks the predicted 2026 observing window available from the ground. Bottom: final residuals after subtracting the periodic TTV model, provided by Eq. (4).

5 Search for additional transiting planets

The well-aligned, nearly edge-on, and circular orbit of HD 332231 b offers favourable conditions for detecting additional, smaller transiting planets – provided they exist and share a similar coplanar alignment. To explore this possibility, we analysed the joint OOT photometric time series using the Analysis of Variance for planetary Transits method (AoVtr; Schwarzenberg-Czerny & Beaulieu 2006). This algorithm is a specialized adaptation of the classical AoV periodogram, tailored to identifying box-shaped periodic signals characteristic of planetary transits.

To construct the OOT dataset, we removed the transits of HD 332231 b from the original light curves, applying a 15-minute margin around the predicted transit windows (based on the ephemeris in Eq. (4)). The light curves from all individual sectors were then combined, yielding a final dataset of 114 061 measurements.

We searched for periodic signals over trial periods ranging from 0.4 to 100 days, corresponding to orbital separations from roughly 2 R_* outward. The frequency resolution was set to 2 × 10⁻⁵ d⁻¹. The AoVtr algorithm folds and bins the light curve for each trial period, fitting a negative-pulse model to identify potential transit-like events. Its sensitivity depends on both the number of phase bins and the transit duration. To optimise performance across the period range, we varied the bin count from 10 to 200 (in steps of 5) for periods shorter than that of HD 332231 b and from 80 to 400 (in steps of 10) for longer periods.

In each periodogram, we identified the strongest peak and evaluated its significance using the signal-detection efficiency metric (Alcock et al. 2000; Kovács et al. 2002), in the generalized form of Ofir (2014). This metric quantifies the signal-to-noise ratio (S/N) in the power spectrum as $$\mathrm{S}/\mathrm{R}=\frac{\mathrm{PWR}-{\mathrm{PWR}}_{\mathrm{MED}}}{{\sigma }_{\mathrm{MAD}}},$$(5) where PWR is the peak power, PWR_MED is the median power in a frequency window around the peak, and σ_MAD is the noise level estimated via the median absolute deviation. As noted by Gondhalekar et al. (2023), there is no strict prescription for defining the width of the median filter. We set its width to one-fiftieth of the frequency span of the periodogram, providing a balance between smoothing the noise and preserving long-term trends. Injection-recovery tests (see below) indicated that signals with S/N>30 are likely to pass robust visual confirmation.

The analysis returned no statistically significant signals indicative of additional transiting planets, either interior or exterior to HD 332231 b. To estimate the detection sensitivity of the current photometric dataset, we performed an injection-recovery analysis. Artificial box-shaped transits were inserted into the combined OOT light curve over periods from 0.4 to 100 days, which were incremented logarithmically (step size of 0.01 dex). We began with transit depths of 10 ppm, increasing by 10 ppm until the S/N exceeded the detection threshold. The injected transit duration was fixed at two hours. For each trial, the number of bins used by the AoVtr algorithm was calculated as the ratio of the trial period to the duration of the artificial transits, constrained between 10 and 400 bins. A periodogram was computed for each synthetic light curve, and the transit depth required to exceed an S/N of 30 was recorded and converted into a planetary radius limit using R_*.

For interior orbits, our analysis was sensitive to planets as small as 0.8 R_⊕ on the shortest period orbits, increasing to 1.9 R_⊕ for ten-day orbits. For exterior configurations, planets up to 2.5−3.0 R_⊕ could be detected out to 35 days. Beyond this, data gaps between TESS sectors reduced detection efficiency, with limiting radii varying between 3 and 13 R_⊕ depending on the extent of transit overlap with observed cadences.

This approach assumes that hypothetical planets maintain stable orbital periods, allowing for coherent phase-folding. However, TTVs comparable to or larger than the transit duration can smear the phase-folded signal and hinder detection. As discussed in Sect. 6, the TTVs observed in HD 332231 b could be caused by an additional planetary companion. For near-circular orbits, the TTV amplitudes of two interacting planets are inversely proportional to their mass ratio and exhibit anti-correlated phases (Deck & Agol 2016; Nesvorný & Vokrouhlický 2016). If such a companion is smaller and less massive than HD 332231 b, its transits could exhibit TTVs with amplitudes significantly exceeding A_TTV, thus complicating detection with our method.

To mitigate this, we reanalysed four shorter data subsets where the orbital period could be considered approximately constant. These subsets combined data from Sectors 14+15, 41+55, 81+82, and 75+81+82. This alternative approach also yielded null results.

6 Interpretation of the TTV signal

Transit timing variations offer a powerful tool for detecting additional planets that gravitationally perturb the orbits of known transiting planets – particularly when those companions are otherwise undetectable through other observational methods (Holman & Murray 2005; Agol et al. 2005). While early attempts to detect TTVs focused on massive, short-period planets using ground-based observations, it was the advent of space-based photometry – most notably from the Kepler mission – that enabled the first practical applications of this technique (Holman et al. 2010; Ballard et al. 2011). However, interpreting TTVs remains a complex inverse problem: different orbital architectures can produce similar timing signals, making it difficult to uniquely constrain the properties of the perturbing body.

Moreover, not all TTVs are necessarily of planetary origin. Other astrophysical mechanisms may also cause variations in transit timing. In the following subsections, we explore both scenarios: non-planetary explanations are considered in Sect. 6.1, while possible planetary perturbations are discussed in Sect. 6.2.

6.1 Non-planetary scenarios

The observed variations in the orbital period of HD 332231 b could be apparent rather than physical, potentially arising from non-planetary phenomena such as the light travel time effect. In the presence of a third body, the barycentre of the star–planet pair would oscillate around the system’s common centre of mass. This geometric effect would manifest as periodic shifts in transit timing, with the orbital period of the third body matching the observed P_TTV. To produce a signal with an amplitude of A_TTV, the radial displacement of the system would need to be approximately A_TTV · c ≈ 5.6 au, corresponding to an RV amplitude of around 25 km s⁻¹.

The longest RV time series was reported by Dalba et al. (2020) and was acquired with the Levy Spectrograph on the 2.4 m Automated Planet Finder (APF) telescope (Radovan et al. 2014; Vogt et al. 2014) at the Lick Observatory. It covers only about 6% of the TTV signal’s phase. On one hand, these observations may have missed such significant RV variations, particularly near the extrema of the signal. This observational limitation is consistent with the absence of long-term RV trends noted by Dalba et al. (2020). On the other hand, to generate the required light travel time effect signal, the mass of the hypothetical companion would need to be at least 5.6 M_⊙ – which is significantly more massive than HD 332231 itself. This makes the scenario implausible because spectroscopic analyses show no evidence of a secondary stellar component (Dalba et al. 2020; MacDougall et al. 2023).

We also examined the possibility of a gravitationally bound brown dwarf or stellar companion in a wide orbit, constituting a hierarchical triple system and capable of producing the observed P_TTV. Using the formalism of Borkovits et al. (2011) and assuming circular orbits for both HD 332231 b and the hypothetical companion, we found no viable solutions that could reproduce the observed A_TTV. Even under the most favourable configuration – mutually perpendicular orbital planes and a companion mass of 50 M_Jup – the predicted TTV amplitude was only half the observed value.

Another potential explanation involves apsidal precession of an eccentric orbit of HD 332231 b. In this scenario, the variations in P_b could be an apparent effect arising from the precession of the planetary orbit. The required eccentricity, estimated as 2 π A_TTV/P_b ≈ 0.01, lies within the 2 σ confidence interval reported by Polanski et al. (2024). We calculated precession timescales considering contributions from general relativity, as well as tidal and rotational deformations of the host star and planet, using the equations of Ragozzine & Wolf (2009). For these calculations, we adopted a stellar tidal Love number of k_{2, *} ≈ 0.0087 (Claret 1995) and a planetary value of k_{2, b}=0.39 for Saturn (Lainey et al. 2017). Assuming a planetary-rotation period of ten hours – consistent with Solar System gas giants – we found that the shortest precession periods, driven by relativistic and rotational effects, are of the order of 10⁵ years. Tidal effects led to even longer timescales, exceeding 10⁷ years. These are vastly inconsistent with the observed TTV period, ruling out precession as a viable explanation.

Lastly, we considered magnetic-activity cycles as a possible cause of timing variations. In such cases, cyclic changes in the host star’s gravitational potential – induced by shape deformations linked to magnetic activity – can modulate the orbital period of a companion (Applegate 1992). Using Eq. (13) from Watson & Marsh (2010), together with the stellar parameters from MacDougall et al. (2023) and a rotation period inferred from the projected equatorial velocity (Knudstrup & Albrecht 2022), we found the expected timing variations to be ∼ 10⁻² seconds only. This is five orders of magnitude smaller than the observed A_TTV. Even when adopting less conservative assumptions regarding the stellar quadrupole moment (Lanza et al. 1998), the effect remains negligible.

Given the implausibility of these non-planetary explanations – whether geometric, dynamical, or magnetic in nature – we were led to consider alternative causes. Planetary perturbations thus emerge as the most likely origin of the observed TTV signal.

6.2 Planetary perturbations

To explore whether gravitational perturbations from an undetected planetary companion could explain the observed TTVs, we performed a suite of reconnaissance simulations using the TTVFast code (Deck et al. 2014). It employs a numerically efficient symplectic integrator to compute transit times in simulated planetary systems. Synthetic mid-transit times were used to construct observed minus calculated (O–C) diagrams, with transit ephemerides refined via linear regression. We then analysed the residuals for sinusoidal patterns using the AoV algorithm, searching over periods from 2 to 7000 epochs at a frequency resolution of 1 × 10⁻⁴ epoch ⁻¹.

We considered perturbers with an orbital period, P_c, ranging from 3 to 350 days (in steps of 2 × 10⁻⁴ days); eccentricity, e_c, from 0 to 0.85 (step size 0.05); argument of periastron from 0^° to 330^° (in steps of 30^°); and a mass, M_c, of 5,10,25, and 50 M_⊕. The mass of HD 332231 b was calculated by adopting the RV amplitude K_b=17.5 ± 1.1 m s⁻¹ from Knudstrup & Albrecht (2022) and the stellar mass M_*=1.16 ± 0.03 M_⊙ from MacDougall et al. (2023).

Each simulation spanned 10 000 days – approximately four full TTV cycles – with integration steps of 0.025 days for inner configurations and 0.1 days for outer ones. Due to the system’s alignment, only prograde and coplanar orbits were considered. Configurations resulting in orbit crossings were excluded. For each model, we assessed whether the simulated TTV period matched P_TTV within uncertainties and whether the required perturber mass to reproduce A_TTV remained below the RV detection threshold defined as $$K_{lim}=\sqrt{2}{\sigma }_{\mathrm{jitter},\text{HIRES}}\approx 4.3\mathrm{m}{\mathrm{s}}^{-1},$$(6) where σ_jitter,HIRES represents the RV jitter in the high-precision measurements obtained with the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on the Keck I telescope, as reported by Dalba et al. (2020). For the purposes of this exercise, we approximated the RV signature of the perturber as a sinusoid for simplicity.

Many configurations successfully reproduced P_TTV, typically clustering near low- and high-order orbital resonances. These are visualised in Fig. 3 (upper panels).

To assess whether the existing RV data could constrain the perturber’s properties more precisely, we attempted to identify the best-fitting two-planet solutions near inner and outer x:1, x:2, and x:3 resonances, where x denotes relevant integer values across the explored period range. Here, the notation p:q refers to the period ratio of inner (p) to outer (q) planets. We combined the TTVFast algorithm with the dynamic nested sampling routine from the dynesty package (Speagle 2020) to explore the likelihood function in its general form: $$ln\mathcal{L}=-\frac{1}{2}{\sum }_i[\frac{{(O_i-{\mathcal{M}}_i)}^2}{{\sigma }_i^2}+ln\left(2\pi {\sigma }_i^2\right)],$$(7) where O_i and σ_i are the observed values and uncertainties, ℳ_i are the corresponding model predictions, and the sum runs over all data points. In addition to the transit timing dataset for HD 332231 b, the observational data were supplemented with the OOT RV data taken from Dalba et al. (2020), comprising 13 high-precision HIRES observations and 68 moderate-precision APF measurements. Data from the SONG telescope, also reported by those authors, were excluded due to their large uncertainties, which were comparable to the RV amplitude induced by HD 332231 b. All timestamps were converted to BJD_TDB to ensure consistent timing for TTV analysis.

The model included 14 free parameters: the orbital period, mass, eccentricity, argument of periastron, and mean anomaly for both HD 332231 b and the trial planet, as well as the offsets and jitter terms for each of the two RV datasets (added in quadrature to the formal uncertainties). Priors were assigned uniformly within limits individually optimised based on our reconnaissance simulations. For each considered resonance, we explored solutions slightly below and above the exact orbital period commensurability. We employed 7000 live points, allowing the parameter space to be sampled until the change in log-evidence fell below Δ ln 𝒵 of 0.01. The solutions with the highest ln ℒ were subsequently fine-tuned using the Nelder–Mead minimisation routine implemented in SciPy (Virtanen et al. 2020), driven by the following negative log-likelihood-based cost function: $$C=-2ln\mathcal{L}.$$(8) The variation of −ln ℒ with the perturber’s orbital period is shown in the middle panels of Fig. 3. In our approach, the RVs primarily constrain the planetary masses and orbital parameters, whereas the transit data tightly define the TTV signature. Relative to the single-planet solution, in which the TTV signal was approximated by a pure sinusoid as in Eq. (4), two-planet configurations with orbital periods of P_c ≲ 60 days do not improve the fit. As shown in the lower panels of Fig. 3, this is primarily because the plausible perturbers are of low mass (1−10 M_⊕, with a median of 5 M_⊕), producing RV amplitudes of only 0.2−3.3 m s⁻¹ (median 1 m s⁻¹). Near the 2:1 resonances, even Earth-mass perturbers are capable of reproducing the observed TTV signal, but their induced Doppler signals fall below current detection thresholds. Similarly, as demonstrated in Sect. 5, the TESS photometry lacks the sensitivity to detect such perturbers unless they are low-density, gaseous planets with radii larger than 1.5 R_⊕ (for inner orbits) or 3 R_⊕ (for outer orbits). Rocky planets in these regimes remain undetectable with current instrumentation.

For longer orbital periods, that is P_c ≳ 60 days, we observed a modest improvement in ln ℒ, forming a shallow plateau for P_c ≳ 120 days. Although no long-term trend is evident in the full RV dataset (Dalba et al. 2020), a linear trend of 0.19 ± 0.04 m s⁻¹ day ⁻¹ appears in the HIRES data from BJD 2 458 815 to 2 458 858. Three additional HIRES measurements near BJD 2 458 917.15 deviate from this trend, suggesting possible curvature that could be indicative of an underlying periodic signal. In Fig. 4, we present three illustrative solutions for a perturber orbiting near the 13:3 resonance (a), 21:2 resonance (b), and 33:2 resonance (c), also indicated in Fig. 3. In each case, the model RV curves reproduce the local HIRES trend, while the TTV curves match the observed mid-transit times. The TTVs become increasingly non-sinusoidal with longer P_c, featuring characteristic flat segments associated with enhanced perturbations near periastron. However, the statistical preference for these long-period configurations over more compact solutions remains limited.

Although our systematic transit search returned no detections, we additionally inspected both de-trended and raw TESS light curves, searching for monotransits, as expected for solutions such as those presented in Fig. 4. For example, model (a) predicts two transit windows (in Sectors 14 and 41), while models (b) and (c) each predict a single window (in Sectors 14 and 55, respectively). Although transits in these cases are effectively ruled out, long-period transits could still evade detection due to incomplete photometric coverage.

Fig. 3 Exploration of potential planetary perturbers responsible for the observed TTV signal in the HD 332231 system. Upper panels: map of simulated configurations in the parameter space of perturber’s orbital period and eccentricity. Red pixels denote solutions that reproduce both the observed TTV period and amplitude while remaining below the RV detection threshold. Violet pixels indicate configurations that match the TTV signal and may be confirmed with the current RV data. Yellow regions correspond to dynamically unstable (orbit-crossing) configurations. Middle panels: Negative log-likelihood metric, −ln ℒ, as a function of the perturber’s orbital period for various resonance zones. The horizontal dashed line indicates the −ln ℒ value for a single-planet model. Bottom panels: perturber mass required to reproduce the observed ATTV across the tested orbital periods. Orange points and dotted vertical lines indicate the three illustrative solutions presented in Fig. 4.

Fig. 4 Comparison of three illustrative two-planet dynamical models fitted to HD 332231. Each row corresponds to a different perturber configuration with orbital periods of about 81 days (a), 196 days (b), and 307 days (c); the corresponding masses and orbital eccentricities are also provided. Left panels: modelled RV signals overlaid on HIRES and APF data points (red and orange points, respectively). Right panels: corresponding TTV model compared to observed mid-transit times. In analogy to Fig. 2, blue points come from TESS light curves, and the green square is from spectral observations by Knudstrup & Albrecht (2022).

7 Conclusions

Our analysis reveals that HD 332231 b exhibits strong periodic TTVs with a cycle of approximately 6.7 years and an amplitude of about 45 minutes. We provide an updated transit ephemeris that includes a periodic term, which will be important for planning future transit observations. Dynamical simulations suggest that the observed TTVs can be explained by gravitational perturbations from an additional, as-yet undetected planetary companion. Although many orbital configurations can reproduce the TTV signal, current photometric and RV data are insufficient to uniquely constrain their properties. However, the RV data show a modest improvement in likelihood for solutions with an outer companion on an orbit longer than about 60 days, but the statistical preference over more compact configurations remains limited.

The system is not expected to be revisited by TESS until August 2027, limiting access to space-based transit data in the short term. In the meantime, ground-based transit observations could significantly improve constraints on the TTV period. Likewise, acquiring new high-precision RV measurements would help test the plausibility of outer perturber scenarios, particularly those near high-order mean-motion resonances. If such a companion exists, our modelling indicates that it would likely occupy a relatively eccentric orbit – strikingly different from the inner planet’s dynamically calm configuration.

The system HD 332231 thus remains a compelling target for continued monitoring. Coordinated efforts involving ground-based transit observations, enhanced RV coverage, and refined dynamical modelling could yield valuable insights into the formation and evolution of multi-planet systems featuring TTVs and eccentric outer perturbers.

Acknowledgements

I would like to thank the anonymous referee for constructive and insightful comments, which helped improve the robustness of the analysis. I acknowledge the financial support from the National Science Centre, Poland, through grant no. 2023/49/B/ST9/00285. This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This research made use of Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration 2018). This research has made use of the SIMBAD database and the VizieR catalogue access tool, operated at CDS, Strasbourg, France, and NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A TESS observations

Fig. A.1 TESS instrumental light curves of HD 332231 in individual observing sectors. Black points represent photometric measurements with quality flag QF=0, which were included in the final analysis. Red points correspond to measurements with QF>0 that were excluded. The segments centred on transits of planet b used in transit modelling (Sect. 3), are shown in blue, with individual transits marked by their epoch numbers.

Appendix B Modelling of the transit light curve in Sector 14

Dalba et al. (2020) identified the transit of HD 332231 b in Sector 14 data that had been flagged as affected by scattered light and were therefore excluded from the original Pre-search Data Conditioning Simple Aperture Photometry (PDC-SAP) analysis. Post-transit observations in this sector also display signs of an Argabrightening event, with erratic flux variations near BTJD 1693.2 approximately seven hours after fourth contact. To recover these compromised data, Dalba et al. (2020) employed a standard aperture photometry approach using Lightkurve, augmented with an advanced systematic correction based on spacecraft quaternion time series within each exposure. The resulting Sector 14 light curve was then included in their analysis and contributed to the refinement of the transit ephemeris.

However, visual inspection of Fig. 4 in Dalba et al. (2020) suggests that the Sector 14 transit appears anomalously long. To further investigate this anomaly, we reprocessed the Sector 14 data using our pipeline, deliberately bypassing the quality flag criteria and retaining all available exposures. A preliminary fit performed with TAP indicated that the Sector 14 transit was ${66}_{-48}^{+11}$ minutes longer and ${1300}_{-900}^{+300}\mathrm{ppm}$ deeper compared to transits observed in other sectors. Although these differences are not statistically significant (at the ∼ 1.5 σ level), they raise concerns regarding the reliability of this dataset. Figures B.1a and B.1b illustrate the Sector 14 transit overlaid on the phase-folded light curves from other sectors, along with their respective best-fitting models, clearly showing the discrepancy. As shown in panel c, the raw instrumental flux exhibits a distinct offset in the OOT baseline, accompanied by mild undulations and clusters of negative outliers.

We conducted a suite of tests, including masking groups of outlying data points and applying relative flux offsets to correct for potential discontinuities between segments. These attempts did not yield satisfactory improvements, suggesting that the underlying systematic effects are more complex than assumed. This exercise ultimately demonstrated that the quality of the transit light curve in Sector 14 renders it unsuitable for detailed analysis. Nevertheless, the initial detection of the transit in this sector played a crucial role in constraining P_b and establishing the preliminary transit ephemeris by Dalba et al. (2020).

Fig. B.1 Transit of HD 332231 b in Sector 14. Panels a and b compare the Sector 14 transit light curve (red points) with the phase-folded data from other sectors (grey points), zoomed in on ingress and egress to highlight discrepancies. The corresponding best-fitting models are shown as red and black/white lines. Panel c shows the raw instrumental flux around the Sector 14 transit, illustrating systematics affecting the measurements.

Appendix C TAP modelling posterior distributions

Fig. C.1 Corner plot showing the posterior distributions of the fitted parameters for the transit model of HD 332231 b, based on the TESS data. The diagonal panels show marginalized 1D distributions, while the off-diagonal panels display joint 2D projections to illustrate parameter covariances. Contours correspond to 1 σ, 2 σ, and 3 σ confidence levels. The orbital parameters eb and ωb were sampled under Gaussian priors taken from Knudstrup & Albrecht (2022). The fit was performed using the TAP model described in Sect. 3, and informed the dynamical simulations discussed in later sections.

Appendix D Search for the evolution of the transit parameters

In a test run, we allowed i_b, a_b/R_*, and R_b/R_* to be independently fitted to the individual transit light curves. We considered only the complete ones to base the subsequent analysis on the most reliable data that limited our sample to 7 light curves. The remaining free parameters were set as in the model described in Sect. 3. The results are collected in Table D.1 together with derived individual values for b_b and τ_{14, b}. As illustrated in Fig. D.1, these transit parameters do not reveal any variation correlated with time, represented by transit epochs E.

Fig. D.1 Parameters for individual complete transits observed with TESS. The best-fitting values for the joint model and their 1 σ uncertainties, reported in Sect. 3, are marked with horizontal continuous and dashed lines, respectively.

Appendix E Preliminary transit timing analysis

Fig. E.1 Results of the preliminary analysis of the transit-timing dataset. Panel a: Transit-timing residuals relative to a trial linear ephemeris. The symbol coding follows that in Fig. 2 and is explained in the legend. Timing uncertainties for the TESS data are smaller than the symbol size. The uncertainty of the trial ephemeris is illustrated by 100 posterior realisations, shown in orange. Panel b: AoV periodogram of the residuals obtained after subtracting the refined linear ephemeris shown in panel a. Horizontal lines indicate the false-alarm probability (FAP) levels, estimated using a bootstrap procedure applied to 100 000 resampled datasets. The strongest peak, near 96 epochs, has a FAP of 1.9% and is marked by a vertical line. Several harmonics of this signal are also visible.

Appendix F Revisiting the Rossiter-McLaughlin effect of HD 332231 b

In their study, Knudstrup & Albrecht (2022) observed a transit of HD 332231 b on UT 2020 August 5 using the High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPS-N; Cosentino et al. 2012) on the 3.6 m Telescopio Nazionale Galileo. The 6.4-hour observing run began near first contact and concluded shortly after last contact. The Doppler data showed a characteristic distortion due to the Rossiter–McLaughlin (RM) effect. Three modelling techniques were applied to extract the sky-projected obliquity λ_b – the angle between the projected stellar spin axis and the planet’s orbital axis. The most precise method, which analysed spectral line distortions, yielded λ_b=−2 ± 6 degrees, indicating that the planet’s orbital plane and the stellar equator are aligned within 1 σ.

In contrast, Sedaghati et al. (2022) reported a significantly different result, based on a transit observed on UT 2020 October 18 using the CARMENES high-resolution spectrograph (Quirrenbach et al. 2014) on the 3.5 m telescope at Calar Alto Observatory. Their observations spanned from shortly before second contact to just after fourth contact. They also collected RV measurements on the nights before and after the transit, each lasting 1.75 hours, to better constrain the RV baseline. Their analysis produced λ_b=−42 ± 11 degrees, suggesting a moderate misalignment.

Knudstrup & Albrecht (2022) attributed this discrepancy to Sedaghati et al.’s use of a linear ephemeris based solely on TESS photometry from Sectors 14 and 15, along with contemporaneous RV data. Indeed, for that epoch, E=22, their ephemeris differed by 14 minutes from that given in Eq. (4) (Sect. 4). To improve on this, Knudstrup & Albrecht incorporated two additional mid-transit times from TESS Sector 41 and introduced a free parameter, Δ T_{0, b}, to account for a potential mid-transit time offset. Their best-fitting model yielded Δ T_{0, b} of ∼ 20 minutes, detected with >5 σ significance. When they fixed T_{0, b} to their updated ephemeris, that is, adopting Δ T_{0, b}=0, the derived obliquity shifted to λ_b=−31 ± 6 degrees, consistent with Sedaghati et al.’s value. This finding underscored the risk of obtaining spurious λ_b estimates when TTVs are neglected. Knudstrup & Albrecht (2022) also questioned the higher projected rotation velocity, $vsin{i}_{\star }=16.3_{-4.4}^{+6.9}\mathrm{km}{\mathrm{s}}^{-1}$, reported by Sedaghati et al. (2022), noting its strong correlation with λ_b (see Albrecht et al. 2011).

To revisit these results, we reanalysed both in-transit RV datasets using the allesfitter package (Günther & Daylan 2019, 2021), employing the ephemeris with the periodic term from Eq. (4). We additionally incorporated the OOT RV data reported by Dalba et al. (2020). The composition of this RV dataset and the criteria used for selecting the measurements are described in Sect. 6.2.

In a trial run, we fitted the transits from Sectors 15–82 to independently determine transit parameters. The results agreed with those from TAP (Sect. 3) within 1 σ, and we subsequently adopted the TAP-derived values as Gaussian priors on R_b/R_*, the normalised sum of radii (R_b+R_*)/a_b, and the inclination through cos i_b. Free parameters in our model included λ_b, the RV semi-amplitude K_b, the orbital eccentricity e_b and argument of periastron ω_b represented by $\sqrt{e_{\mathrm{b}}}cos{\omega }_{\mathrm{b}}$ and $\sqrt{e_{\mathrm{b}}}sin{\omega }_{\mathrm{b}}$, as well as instrumental offsets and jitter terms. Following Knudstrup & Albrecht (2022), we modelled stellar limb darkening using a quadratic law with coefficients u₁ and u₂, adopting nominal values interpolated from Claret & Bloemen (2011) for the V band. These were allowed to vary with Gaussian priors (σ=0.1 for u₁ and 0.2 for u₂; Patel & Espinoza 2022). The alternative reparametrisation, based on the q-space (Kipping 2013), yielded consistent results, so we retained the physical u-form to match our TAP modelling. The space of free parameters was explored with dynamic nested sampling with 500 live points and default settings (Günther & Daylan 2021).

We attempted to jointly model the RM signals from both HARPS-N and CARMENES. However, the resulting fits were unsatisfactory. Figures F.1a and F.1b show clearly that the RM amplitudes differ markedly. Since the RM amplitude scales approximately as ${(R_{\mathrm{b}}/R_{\star })}^2\sqrt{1-b_{\mathrm{b}}^2}(vsin{i}_{\star })$ (Triaud 2018), and the planetary/stellar radii and impact parameter must be identical, this discrepancy implies a difference in v sin i_*. As such, we proceeded to model the transits separately.

Fig. F.1 Rossiter-McLaughlin effect observed during two transits of HD 332231 b. Panel a: HARPS-N observations from Knudstrup & Albrecht (2022), showing the in-transit RV anomaly consistent with a well-aligned orbit. Blue points denote the measured RVs with 1 σ uncertainties. The solid red line shows the best-fitting model, while the orange-shaded curves represent 20 posterior samples illustrating the range of plausible solutions. The dashed gray line indicates the higher-amplitude RM model based on the CARMENES data. Panel b: CARMENES observations from Sedaghati et al. (2022), shown as green points with the corresponding best-fitting model. The dashed gray line shows the HARPS-N-based solution for comparison. Panel c: Same CARMENES data, re-fitted with an added linear RV trend (dashed gray line), yielding a reduced RM amplitude. Panel d: OOT CARMENES measurements taken the night after the transit, showing a non-zero RV slope, also indicated with a dashed gray line. The RV changes induced by the orbital motion of HD 332231 b are shown in red and remain negligible in this scale.

Although v sin i_* can, in principle, be derived from RM data alone, central transits (low b_b or cos i_b near zero) introduce a degeneracy between v sin i_* and λ_b, degrading parameter precision (Albrecht et al. 2011). Our tests confirmed this effect for HD 332231 b, where the posterior distribution of v sin i_* showed a long high-velocity tail, inconsistent with spectral line measurements. We therefore adopted a Gaussian prior of v sin i_*=5.3 ± 1.0 km s⁻¹, based on the HIRES spectra (Dalba et al. 2020).

Key results are shown in Table F.1. All other parameters agreed with literature values within uncertainties and are omitted for brevity. The best-fitting models (solid lines) and 20 posterior draws (shaded lines) are shown in Fig. F.1. To highlight differences in the RM amplitudes, dashed lines show the alternative model in each panel.

In both transits, λ_b is consistent with 0^°, supporting alignment and confirming findings of Knudstrup & Albrecht (2022). Our HARPS-N value of v sin i_* agrees within 0.8−1.1 σ with the values of 5.3 ± 1.0,5.4 ± 1.0, and 7.0 ± 0.5 km s⁻¹, derived from spectral analysis by Dalba et al. (2020). It also agrees within 0.9−1.5 σ with the values of 5.64 ± 0.14,5.63 ± 0.11, and $5.{89}_{-0.13}^{+0.12}\mathrm{km}{\mathrm{s}}^{-1}$, obtained using three methods by Knudstrup & Albrecht (2022).

However, our CARMENES-derived v sin i_*=10.0 ± 0.6 km s⁻¹ is higher by 4−7 σ, suggesting potential systematics. The residuals in Fig. F.1b reveal outliers at the end of the CARMENES series. Sedaghati et al. (2022) masked the final four OOT points, citing increased airmass but without identifying a cause. We explored this further by fitting a model with a linear baseline trend (“CARMENES+slope”; Fig. F.1c). This model yielded ${\lambda }_{\mathrm{b}}=0_{-20}^{+21}$ degrees and a lower v sin i_*=6.3 ± 0.7 km s⁻¹, consistent with other datasets. The OOT RVs taken one day later with CARMENES (Fig. F.1d) also show a slope of −27.4 ± 9.2 ms⁻¹ day⁻¹, differing from zero at nearly 3 σ, which supports using a baseline slope in the RM model.

Mid-transit times from our RM analysis are also listed in Table F.1. For HARPS-N, our result matches Knudstrup & Albrecht (2022) within 1 σ and is more precise. Nevertheless, we use their value in the transit timing analysis, as it was obtained via the planetary shadow method, which provides robust results.

For CARMENES, the situation is more complex. Although including a baseline trend reduced the RM amplitude, it doubled the timing uncertainty and shifted T_{mid, b} by 26 minutes (≈ 2.2 σ effect). Given these uncertainties, we excluded the CARMENES data from our timing analysis. Actual in-night trends may have a more complex nature, and a refined reduction of the original observations may help clarify this issue.

All Tables

All Figures

	Fig. 1 Transit light curves of HD 332231 b observed with TESS in Sectors 15–82. Left: Individual photometric time series sorted by the epoch of observation, with numbering consistent with the refined ephemeris given in Sect. 4. Best-fitting models are overlaid in red. Right: Corresponding photometric residuals from the transit models.
In the text

	Fig. 2 TTVs of HD 332231 b. Top: residuals after subtracting the linear ephemeris. The sinusoidal TTV model and 100 posterior models are shown. The grey shaded region marks the predicted 2026 observing window available from the ground. Bottom: final residuals after subtracting the periodic TTV model, provided by Eq. (4).
In the text

Fig. 3 Exploration of potential planetary perturbers responsible for the observed TTV signal in the HD 332231 system. Upper panels: map of simulated configurations in the parameter space of perturber’s orbital period and eccentricity. Red pixels denote solutions that reproduce both the observed TTV period and amplitude while remaining below the RV detection threshold. Violet pixels indicate configurations that match the TTV signal and may be confirmed with the current RV data. Yellow regions correspond to dynamically unstable (orbit-crossing) configurations. Middle panels: Negative log-likelihood metric, −ln ℒ, as a function of the perturber’s orbital period for various resonance zones. The horizontal dashed line indicates the −ln ℒ value for a single-planet model. Bottom panels: perturber mass required to reproduce the observed ATTV across the tested orbital periods. Orange points and dotted vertical lines indicate the three illustrative solutions presented in Fig. 4.

In the text

Fig. 4 Comparison of three illustrative two-planet dynamical models fitted to HD 332231. Each row corresponds to a different perturber configuration with orbital periods of about 81 days (a), 196 days (b), and 307 days (c); the corresponding masses and orbital eccentricities are also provided. Left panels: modelled RV signals overlaid on HIRES and APF data points (red and orange points, respectively). Right panels: corresponding TTV model compared to observed mid-transit times. In analogy to Fig. 2, blue points come from TESS light curves, and the green square is from spectral observations by Knudstrup & Albrecht (2022).

In the text

Fig. A.1 TESS instrumental light curves of HD 332231 in individual observing sectors. Black points represent photometric measurements with quality flag QF=0, which were included in the final analysis. Red points correspond to measurements with QF>0 that were excluded. The segments centred on transits of planet b used in transit modelling (Sect. 3), are shown in blue, with individual transits marked by their epoch numbers.

In the text

Fig. B.1 Transit of HD 332231 b in Sector 14. Panels a and b compare the Sector 14 transit light curve (red points) with the phase-folded data from other sectors (grey points), zoomed in on ingress and egress to highlight discrepancies. The corresponding best-fitting models are shown as red and black/white lines. Panel c shows the raw instrumental flux around the Sector 14 transit, illustrating systematics affecting the measurements.

In the text

Fig. C.1 Corner plot showing the posterior distributions of the fitted parameters for the transit model of HD 332231 b, based on the TESS data. The diagonal panels show marginalized 1D distributions, while the off-diagonal panels display joint 2D projections to illustrate parameter covariances. Contours correspond to 1 σ, 2 σ, and 3 σ confidence levels. The orbital parameters eb and ωb were sampled under Gaussian priors taken from Knudstrup & Albrecht (2022). The fit was performed using the TAP model described in Sect. 3, and informed the dynamical simulations discussed in later sections.

In the text

	Fig. D.1 Parameters for individual complete transits observed with TESS. The best-fitting values for the joint model and their 1 σ uncertainties, reported in Sect. 3, are marked with horizontal continuous and dashed lines, respectively.
In the text

Fig. E.1 Results of the preliminary analysis of the transit-timing dataset. Panel a: Transit-timing residuals relative to a trial linear ephemeris. The symbol coding follows that in Fig. 2 and is explained in the legend. Timing uncertainties for the TESS data are smaller than the symbol size. The uncertainty of the trial ephemeris is illustrated by 100 posterior realisations, shown in orange. Panel b: AoV periodogram of the residuals obtained after subtracting the refined linear ephemeris shown in panel a. Horizontal lines indicate the false-alarm probability (FAP) levels, estimated using a bootstrap procedure applied to 100 000 resampled datasets. The strongest peak, near 96 epochs, has a FAP of 1.9% and is marked by a vertical line. Several harmonics of this signal are also visible.

In the text

Fig. F.1 Rossiter-McLaughlin effect observed during two transits of HD 332231 b. Panel a: HARPS-N observations from Knudstrup & Albrecht (2022), showing the in-transit RV anomaly consistent with a well-aligned orbit. Blue points denote the measured RVs with 1 σ uncertainties. The solid red line shows the best-fitting model, while the orange-shaded curves represent 20 posterior samples illustrating the range of plausible solutions. The dashed gray line indicates the higher-amplitude RM model based on the CARMENES data. Panel b: CARMENES observations from Sedaghati et al. (2022), shown as green points with the corresponding best-fitting model. The dashed gray line shows the HARPS-N-based solution for comparison. Panel c: Same CARMENES data, re-fitted with an added linear RV trend (dashed gray line), yielding a reduced RM amplitude. Panel d: OOT CARMENES measurements taken the night after the transit, showing a non-zero RV slope, also indicated with a dashed gray line. The RV changes induced by the orbital motion of HD 332231 b are shown in red and remain negligible in this scale.

In the text