TOI-791 b recovery and a publication-informed c-period ranking¶

Question: Can a compact, transparent Box Least Squares (BLS) pipeline broadly recover TOI-791 b and show how the two c-like TESS event windows rank under a publication-informed period search?

This notebook deliberately stops short of claiming an independent planet confirmation. The peer-reviewed study used detailed transit-timing models, ground-based follow-up, high-resolution imaging, and spectroscopy. Here, the narrower goal is to test whether the periods themselves are visible in the public space-based photometry.

Primary references:

  • Peer-reviewed study: https://doi.org/10.1093/mnras/stag864
  • NASA announcement: https://science.nasa.gov/missions/tess/nasas-tess-mission-reveals-the-puffiest-planets-ever-found/
  • TESS archive/data products: https://archive.stsci.edu/missions-and-data/tess/data-products
  • ExoFOP target page: https://exofop.ipac.caltech.edu/tess/target.php?id=306472057
In [1]:
from __future__ import annotations

import csv
import hashlib
import json
import warnings
from datetime import datetime, timezone
from importlib.metadata import version
from pathlib import Path

warnings.filterwarnings("ignore", message="Warning: the tpfmodel")

import astropy.units as u
import lightkurve as lk
import matplotlib.pyplot as plt
import numpy as np
from astropy.timeseries import BoxLeastSquares

PROJECT_ROOT = Path(__file__).resolve().parent if "__file__" in globals() else Path.cwd()
CACHE_DIR = PROJECT_ROOT / ".cache" / "lightkurve"
OUTPUT_DIR = PROJECT_ROOT / "outputs"
CACHE_DIR.mkdir(parents=True, exist_ok=True)
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

TARGET = "TIC 306472057"  # TOI-791
QUERY_TIMESTAMP_UTC = datetime.now(timezone.utc).isoformat(timespec="seconds")
CADENCE_SECONDS = 120
BIN_SIZE = 30 * u.min
FLATTEN_WINDOW = 4001  # about 5.6 days at two-minute cadence
SEARCH_DURATIONS = np.linspace(0.30, 0.90, 9) * u.day
MASK_HALF_WIDTH_DAYS = 0.80

PUBLISHED = {
    "b": {
        "period_days": 139.29931,
        "duration_hours": 12.366,
        "depth_ppt": 5.428,
        "radius_jupiter": 0.993,
        "mass_earth": 9.5,
        "density_g_cm3": 0.037,
    },
    "c": {
        "period_days": 232.01570,
        "duration_hours": 11.926,
        "depth_ppt": 6.78,
        "radius_jupiter": 1.155,
        "mass_earth": 18.6,
        "density_g_cm3": 0.047,
    },
}

print(f"Lightkurve {lk.__version__}")
print(f"Target: {TARGET}")
print("Archive query details are recorded in the machine-readable manifest.")

Lightkurve 2.6.0
Target: TIC 306472057
Archive query details are recorded in the machine-readable manifest.

1. Discover and download the public TESS products¶

The paper reports 44 observed sectors, including long-cadence products. For a uniform analysis, this notebook keeps only official SPOC light curves at two-minute cadence. quality_bitmask="default" removes cadences carrying the standard set of serious TESS quality flags. Lightkurve selects PDCSAP_FLUX for these products, which has already been corrected for common instrumental trends and crowding effects by the mission pipeline.

In [2]:
search_result = lk.search_lightcurve(TARGET, mission="TESS", author="SPOC")
two_minute_mask = np.isclose(np.asarray(search_result.table["exptime"], dtype=float), CADENCE_SECONDS)
two_minute_result = search_result[two_minute_mask]

if len(two_minute_result) == 0:
    raise RuntimeError("MAST returned no two-minute SPOC light curves for TOI-791.")

lightcurves = two_minute_result.download_all(
    download_dir=str(CACHE_DIR),
    quality_bitmask="default",
)
if lightcurves is None or len(lightcurves) != len(two_minute_result):
    raise RuntimeError("One or more TESS products did not download successfully.")


def sector_number(mission_name: str) -> int:
    """Extract the sector number from a string such as 'TESS Sector 04'."""

    return int(str(mission_name).split()[-1])


sectors = [sector_number(value) for value in two_minute_result.table["mission"]]
print(f"Downloaded {len(lightcurves)} two-minute SPOC light curves.")
print("Sectors:", sectors)

# Save a compact data manifest so the exact archive response is auditable.
manifest_path = OUTPUT_DIR / "toi791_tess_product_manifest.csv"
manifest_columns = [
    name
    for name in ("mission", "author", "exptime", "year", "productFilename", "dataURI")
    if name in two_minute_result.table.colnames
]
with manifest_path.open("w", newline="", encoding="utf-8") as handle:
    writer = csv.writer(handle)
    writer.writerow([*manifest_columns, "sha256", "query_timestamp_utc"])
    for row, lightcurve in zip(two_minute_result.table, lightcurves, strict=True):
        local_path = Path(lightcurve.filename)
        digest = hashlib.sha256(local_path.read_bytes()).hexdigest()
        writer.writerow(
            [*[row[name] for name in manifest_columns], digest, QUERY_TIMESTAMP_UTC]
        )

Downloaded 40 two-minute SPOC light curves.
Sectors: [4, 5, 6, 8, 9, 10, 12, 13, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 61, 62, 63, 64, 65, 66, 67, 68, 69, 87, 88, 89, 90, 93, 94, 95, 96, 97, 98]

2. Clean each sector without erasing the long transit¶

Each sector is processed separately:

  1. remove missing values;
  2. normalize its median flux to one;
  3. remove slow trends with a roughly 5.6-day Savitzky-Golay window;
  4. clip strong upward outliers while using a loose lower threshold so real transits are retained;
  5. stitch the sectors and bin to 30 minutes.

The planets' published transits last about 12 hours, so the detrending window is intentionally much longer than a transit. Even so, this processing can attenuate transit depths. The analysis therefore treats period recovery, not depth measurement, as its primary result.

In [3]:
def valid_odd_window(length: int, requested: int = FLATTEN_WINDOW) -> int:
    """Return a valid odd Savitzky-Golay window shorter than the light curve."""

    candidate = min(requested, length - 1 if length % 2 == 0 else length - 2)
    candidate = candidate if candidate % 2 == 1 else candidate - 1
    if candidate < 101:
        raise ValueError(f"Light curve is too short to flatten safely ({length} cadences).")
    return candidate


cleaned_sectors = []
for lightcurve in lightcurves:
    if lightcurve.meta.get("FLUX_ORIGIN") != "pdcsap_flux":
        raise RuntimeError("Expected PDCSAP flux, but Lightkurve selected another flux column.")

    cleaned = lightcurve.remove_nans().normalize()
    cleaned = cleaned.flatten(
        window_length=valid_odd_window(len(cleaned)),
        polyorder=2,
        break_tolerance=5,
        niters=3,
        sigma=5,
    )
    cleaned = cleaned.remove_outliers(sigma_lower=20, sigma_upper=5)
    cleaned_sectors.append(cleaned)

stitched = lk.LightCurveCollection(cleaned_sectors).stitch().remove_nans()
binned = stitched.bin(time_bin_size=BIN_SIZE).remove_nans()

time_days = np.asarray(binned.time.value, dtype=float)
relative_flux = np.asarray(binned.flux.value, dtype=float)
relative_flux_error = np.asarray(binned.flux_err.value, dtype=float)
finite = (
    np.isfinite(time_days)
    & np.isfinite(relative_flux)
    & np.isfinite(relative_flux_error)
    & (relative_flux_error > 0)
)
time_days = time_days[finite]
relative_flux = relative_flux[finite]
relative_flux_error = relative_flux_error[finite]

# Prevent a small number of unusually tiny pipeline uncertainties from
# dominating the BLS statistic.
relative_flux_error = np.maximum(relative_flux_error, np.percentile(relative_flux_error, 5))

print(f"Cleaned cadences: {len(stitched):,}")
print(f"Thirty-minute bins: {len(time_days):,}")
print(f"Time baseline: {np.ptp(time_days):,.1f} days")

Cleaned cadences: 656,247
Thirty-minute bins: 44,149
Time baseline: 2,635.0 days

3. Search for the strongest long-period box-shaped signal¶

Box Least Squares tests many candidate periods and transit durations for a repeating, box-shaped dip. The first search spans 80–280 days. Its power is a useful ranking statistic, but it is not a formal discovery significance: it does not include all look-elsewhere effects, instrumental systematics, or astrophysical false-positive tests. A dense local grid then refines the numerical maximum, but no formal period uncertainty is estimated.

In [4]:
def scalar(value) -> float:
    """Convert an Astropy scalar or NumPy scalar to a Python float."""

    return float(value.value if hasattr(value, "value") else value)


def run_bls(
    time: np.ndarray,
    flux: np.ndarray,
    flux_error: np.ndarray,
    minimum_period: float,
    maximum_period: float,
    grid_size: int,
):
    periods = np.linspace(minimum_period, maximum_period, grid_size) * u.day
    model = BoxLeastSquares(time * u.day, flux, flux_error)
    result = model.power(periods, SEARCH_DURATIONS, objective="snr")
    best_index = int(np.nanargmax(np.asarray(result.power)))
    return model, result, best_index


def refine_bls_peak(
    model: BoxLeastSquares,
    coarse_result,
    coarse_best_index: int,
    half_width_days: float = 0.10,
    grid_size: int = 4001,
):
    """Refine a coarse BLS maximum while keeping its preferred duration fixed.

    The dense grid locates the numerical peak more accurately, but it does not
    provide a formal period uncertainty. Results are therefore reported to only
    two decimal places in human-facing text.
    """

    center = scalar(coarse_result.period[coarse_best_index])
    duration = scalar(coarse_result.duration[coarse_best_index]) * u.day
    periods = np.linspace(center - half_width_days, center + half_width_days, grid_size) * u.day
    result = model.power(periods, duration, objective="snr")
    best_index = int(np.nanargmax(np.asarray(result.power)))
    return result, best_index


def recovered_signal(name: str, result, best_index: int, coarse_spacing_days: float) -> dict:
    """Create a labeled record for a refined BLS peak."""

    return {
        "planet": name,
        "period_days": scalar(result.period[best_index]),
        "duration_days": scalar(result.duration[best_index]),
        "transit_time_btjd": scalar(result.transit_time[best_index]),
        "depth_relative_flux": scalar(result.depth[best_index]),
        "bls_depth_snr": scalar(result.power[best_index]),
        "coarse_grid_spacing_days": coarse_spacing_days,
        "reporting_precision_days": 0.01,
        "formal_period_uncertainty_estimated": False,
    }


bls_b, periodogram_b, best_b_index = run_bls(
    time_days,
    relative_flux,
    relative_flux_error,
    minimum_period=80,
    maximum_period=280,
    grid_size=8000,
)
periodogram_b_refined, best_b_refined_index = refine_bls_peak(
    bls_b, periodogram_b, best_b_index
)

recovered_b = recovered_signal(
    "TOI-791 b",
    periodogram_b_refined,
    best_b_refined_index,
    coarse_spacing_days=200 / (8000 - 1),
)

print("First-pass strongest signal:", recovered_b)

First-pass strongest signal: {'planet': 'TOI-791 b', 'period_days': 139.30531342667834, 'duration_days': 0.44999999999999996, 'transit_time_btjd': 1427.6069603162136, 'depth_relative_flux': 0.0044422361635901705, 'bls_depth_snr': 31.32059829125975, 'coarse_grid_spacing_days': 0.025003125390673835, 'reporting_precision_days': 0.01, 'formal_period_uncertainty_estimated': False}

4. Mask TOI-791 b, then search for a second signal¶

The first planet creates harmonics and aliases in an irregularly sampled survey. In a deliberately publication-informed step, remove a ±0.8-day window around each predicted TOI-791 b transit and repeat the BLS search from 170–280 days. This ranks candidate periods; it does not independently resolve every integer-cycle alias for the sparsely observed c-like events.

In [5]:
def signed_phase_days(time: np.ndarray, period: float, transit_time: float) -> np.ndarray:
    """Return signed phase in days, centered on a transit at phase zero."""

    return ((time - transit_time + 0.5 * period) % period) - 0.5 * period


phase_b_days = signed_phase_days(
    time_days,
    recovered_b["period_days"],
    recovered_b["transit_time_btjd"],
)
keep_after_b_mask = np.abs(phase_b_days) > MASK_HALF_WIDTH_DAYS

bls_c, periodogram_c, best_c_index = run_bls(
    time_days[keep_after_b_mask],
    relative_flux[keep_after_b_mask],
    relative_flux_error[keep_after_b_mask],
    minimum_period=170,
    maximum_period=280,
    grid_size=7000,
)
periodogram_c_refined, best_c_refined_index = refine_bls_peak(
    bls_c, periodogram_c, best_c_index
)

recovered_c = recovered_signal(
    "TOI-791 c",
    periodogram_c_refined,
    best_c_refined_index,
    coarse_spacing_days=110 / (7000 - 1),
)

print("Second-pass strongest signal:", recovered_c)


def observed_event_windows(
    time: np.ndarray,
    flux: np.ndarray,
    signal: dict,
    minimum_depth_ppt: float = 2.0,
) -> list[dict]:
    """Identify observed transit-like windows aligned by a candidate period.

    This is an event inventory, not a planet-validation test. It makes the
    sparse support for the long-period c signal explicit.
    """

    phase = signed_phase_days(time, signal["period_days"], signal["transit_time_btjd"])
    cycles = np.rint(
        (time - signal["transit_time_btjd"]) / signal["period_days"]
    ).astype(int)
    events = []
    for cycle in np.unique(cycles[np.abs(phase) < 0.75]):
        in_transit = (cycles == cycle) & (np.abs(phase) < 0.25)
        local_baseline = (
            (cycles == cycle)
            & (np.abs(phase) >= 0.60)
            & (np.abs(phase) < 0.75)
        )
        if in_transit.sum() < 5 or local_baseline.sum() < 5:
            continue
        depth_ppt = (
            np.median(flux[local_baseline]) - np.median(flux[in_transit])
        ) * 1000
        if depth_ppt >= minimum_depth_ppt:
            events.append(
                {
                    "cycle_index": int(cycle),
                    "predicted_center_btjd": signal["transit_time_btjd"]
                    + cycle * signal["period_days"],
                    "in_transit_bins": int(in_transit.sum()),
                    "local_baseline_bins": int(local_baseline.sum()),
                    "local_median_depth_ppt": float(depth_ppt),
                }
            )
    return events


c_event_windows = observed_event_windows(time_days, relative_flux, recovered_c)
c_events_path = OUTPUT_DIR / "toi791_c_event_windows.csv"
c_event_fieldnames = [
    "cycle_index",
    "predicted_center_btjd",
    "in_transit_bins",
    "local_baseline_bins",
    "local_median_depth_ppt",
]
with c_events_path.open("w", newline="", encoding="utf-8") as handle:
    writer = csv.DictWriter(handle, fieldnames=c_event_fieldnames)
    writer.writeheader()
    writer.writerows(c_event_windows)

print(f"TOI-791 c event windows in selected data: {len(c_event_windows)}")
for event in c_event_windows:
    print(event)

Second-pass strongest signal: {'planet': 'TOI-791 c', 'period_days': 232.00771453064723, 'duration_days': 0.39, 'transit_time_btjd': 1505.9549603162136, 'depth_relative_flux': 0.004206795943946173, 'bls_depth_snr': 17.423797867038008, 'coarse_grid_spacing_days': 0.015716530932990427, 'reporting_precision_days': 0.01, 'formal_period_uncertainty_estimated': False}
TOI-791 c event windows in selected data: 2
{'cycle_index': 3, 'predicted_center_btjd': np.float64(2201.9781039081554), 'in_transit_bins': 24, 'local_baseline_bins': 8, 'local_median_depth_ppt': 4.4568415874635114}
{'cycle_index': 7, 'predicted_center_btjd': np.float64(3130.0089620307444), 'in_transit_bins': 24, 'local_baseline_bins': 14, 'local_median_depth_ppt': 5.940103055872337}

5. Compare the recovered periods with the peer-reviewed values¶

A dense local grid refines each numerical peak after the broad search, but this notebook does not estimate formal period uncertainties. Human-facing results are therefore rounded to 0.01 day. The study also reports transit-timing variations of up to 50 minutes, and its physical model is much more sophisticated than the box model used here.

In [6]:
recovered = {"b": recovered_b, "c": recovered_c}
comparison_rows = []
for label in ("b", "c"):
    measured = recovered[label]
    reference = PUBLISHED[label]
    comparison_rows.append(
        {
            "planet": measured["planet"],
            "refined_bls_peak_period_days": measured["period_days"],
            "reported_recovered_period_days": round(measured["period_days"], 2),
            "published_period_days": reference["period_days"],
            "absolute_peak_difference_days": abs(
                measured["period_days"] - reference["period_days"]
            ),
            "coarse_grid_spacing_days": measured["coarse_grid_spacing_days"],
            "reporting_precision_days": measured["reporting_precision_days"],
            "formal_period_uncertainty_estimated": measured[
                "formal_period_uncertainty_estimated"
            ],
            "bls_box_depth_ppt": measured["depth_relative_flux"] * 1000,
            "published_transit_depth_ppt": reference["depth_ppt"],
            "bls_box_duration_hours": measured["duration_days"] * 24,
            "published_duration_hours": reference["duration_hours"],
            "bls_depth_snr_search_statistic": measured["bls_depth_snr"],
        }
    )

results_path = OUTPUT_DIR / "toi791_recovery_results.csv"
with results_path.open("w", newline="", encoding="utf-8") as handle:
    writer = csv.DictWriter(handle, fieldnames=list(comparison_rows[0].keys()))
    writer.writeheader()
    writer.writerows(comparison_rows)

summary_path = OUTPUT_DIR / "toi791_recovery_summary.json"
package_versions = {
    package: version(package)
    for package in ("astropy", "lightkurve", "matplotlib", "numpy")
}
summary = {
    "target": TARGET,
    "query_timestamp_utc": QUERY_TIMESTAMP_UTC,
    "package_versions": package_versions,
    "flux_column": "pdcsap_flux",
    "number_of_two_minute_spoc_products": len(lightcurves),
    "sectors": sectors,
    "cleaned_cadences": len(stitched),
    "thirty_minute_bins": len(time_days),
    "baseline_days": float(np.ptp(time_days)),
    "results": comparison_rows,
    "toi791_c_event_windows": c_event_windows,
    "scientific_scope": (
        "Approximate period-peak recovery from public TESS photometry. The c "
        "result is a publication-informed ranking supported by two TESS event "
        "windows, not independent alias resolution or planet confirmation."
    ),
}
summary_path.write_text(json.dumps(summary, indent=2), encoding="utf-8")

for row in comparison_rows:
    print(
        f"{row['planet']}: refined numerical peak "
        f"{row['refined_bls_peak_period_days']:.6f} d; report as "
        f"{row['reported_recovered_period_days']:.2f} d; published "
        f"{row['published_period_days']:.5f} d; no formal uncertainty estimated"
    )

TOI-791 b: refined numerical peak 139.305313 d; report as 139.31 d; published 139.29931 d; no formal uncertainty estimated
TOI-791 c: refined numerical peak 232.007715 d; report as 232.01 d; published 232.01570 d; no formal uncertainty estimated

6. Publication-quality diagnostic chart¶

The top panel shows the actual observation coverage. The middle panels show the two BLS searches. The bottom panels fold all observations onto each recovered period so repeated transit events line up at zero hours.

In [7]:
def robust_phase_bins(
    phase_hours: np.ndarray,
    flux_ppt: np.ndarray,
    limit_hours: float = 36,
    width_hours: float = 0.75,
):
    """Median-bin a folded light curve without assuming independent errors."""

    edges = np.arange(-limit_hours, limit_hours + width_hours, width_hours)
    centers = 0.5 * (edges[:-1] + edges[1:])
    bin_index = np.digitize(phase_hours, edges) - 1
    medians = np.full_like(centers, np.nan, dtype=float)
    counts = np.zeros_like(centers, dtype=int)
    for index in range(len(centers)):
        values = flux_ppt[bin_index == index]
        values = values[np.isfinite(values)]
        counts[index] = len(values)
        if len(values) >= 3:
            medians[index] = np.median(values)
    valid = np.isfinite(medians)
    return centers[valid], medians[valid], counts[valid]


NAVY = "#15253D"
TEAL = "#148A8A"
CORAL = "#DB5A42"
GOLD = "#D6A21E"
LIGHT = "#DDE6EC"

plt.rcParams.update(
    {
        "font.family": "DejaVu Sans",
        "font.size": 10,
        "axes.titlesize": 12,
        "axes.labelsize": 10,
        "axes.edgecolor": NAVY,
        "axes.labelcolor": NAVY,
        "xtick.color": NAVY,
        "ytick.color": NAVY,
        "text.color": NAVY,
    }
)

figure = plt.figure(figsize=(14, 13), constrained_layout=True)
grid = figure.add_gridspec(3, 2, height_ratios=[0.8, 1, 1])
ax_coverage = figure.add_subplot(grid[0, :])
ax_bls_b = figure.add_subplot(grid[1, 0])
ax_bls_c = figure.add_subplot(grid[1, 1])
ax_fold_b = figure.add_subplot(grid[2, 0])
ax_fold_c = figure.add_subplot(grid[2, 1])

flux_ppt = (relative_flux - 1.0) * 1000
coverage_low, coverage_high = np.percentile(flux_ppt, [0.5, 99.5])
ax_coverage.scatter(time_days, flux_ppt, s=2, color=NAVY, alpha=0.22, rasterized=True)
ax_coverage.set(
    title=f"Public TESS coverage: {len(lightcurves)} two-minute SPOC products",
    xlabel="Time (BTJD = BJD - 2,457,000)",
    ylabel="Detrended flux (ppt)",
    ylim=(coverage_low, coverage_high),
)
ax_coverage.grid(color=LIGHT, linewidth=0.7, alpha=0.8)

period_b_values = np.asarray(periodogram_b.period.value, dtype=float)
power_b_values = np.asarray(periodogram_b.power, dtype=float)
ax_bls_b.plot(period_b_values, power_b_values, color=TEAL, linewidth=1.2)
ax_bls_b.axvline(PUBLISHED["b"]["period_days"], color=NAVY, linestyle=":", linewidth=1.6)
ax_bls_b.scatter(
    [recovered_b["period_days"]],
    [recovered_b["bls_depth_snr"]],
    color=CORAL,
    s=50,
    zorder=5,
    label=f"Approx. peak: {recovered_b['period_days']:.2f} d",
)
ax_bls_b.set(
    title="First BLS search: strongest signal is TOI-791 b",
    xlabel="Candidate period (days)",
    ylabel="BLS depth S/N (search statistic)",
)
ax_bls_b.legend(frameon=False, loc="upper right")
ax_bls_b.grid(color=LIGHT, linewidth=0.7, alpha=0.8)

period_c_values = np.asarray(periodogram_c.period.value, dtype=float)
power_c_values = np.asarray(periodogram_c.power, dtype=float)
ax_bls_c.plot(period_c_values, power_c_values, color=GOLD, linewidth=1.2)
ax_bls_c.axvline(PUBLISHED["c"]["period_days"], color=NAVY, linestyle=":", linewidth=1.6)
ax_bls_c.scatter(
    [recovered_c["period_days"]],
    [recovered_c["bls_depth_snr"]],
    color=CORAL,
    s=50,
    zorder=5,
    label=f"Approx. peak: {recovered_c['period_days']:.2f} d",
)
ax_bls_c.set(
    title="Targeted search after masking b: c ranks highest",
    xlabel="Candidate period (days)",
    ylabel="BLS depth S/N (search statistic)",
)
ax_bls_c.legend(frameon=False, loc="upper right")
ax_bls_c.grid(color=LIGHT, linewidth=0.7, alpha=0.8)


def draw_folded_panel(
    axis,
    label: str,
    measured: dict,
    color: str,
    bls_model: BoxLeastSquares,
):
    phase_hours = (
        signed_phase_days(time_days, measured["period_days"], measured["transit_time_btjd"])
        * 24
    )
    near_transit = np.abs(phase_hours) <= 36
    axis.scatter(
        phase_hours[near_transit],
        flux_ppt[near_transit],
        s=8,
        color=NAVY,
        alpha=0.16,
        rasterized=True,
        label="30-min observations",
    )
    centers, medians, _ = robust_phase_bins(
        phase_hours[near_transit], flux_ppt[near_transit]
    )
    axis.plot(
        centers,
        medians,
        "o",
        markersize=4,
        color=color,
        label="45-min robust median",
    )
    model_phase_hours = np.linspace(-36, 36, 1000)
    model_times = measured["transit_time_btjd"] + model_phase_hours / 24
    model_flux = bls_model.model(
        model_times * u.day,
        measured["period_days"] * u.day,
        measured["duration_days"] * u.day,
        measured["transit_time_btjd"] * u.day,
    )
    axis.plot(
        model_phase_hours,
        (np.asarray(model_flux, dtype=float) - 1.0) * 1000,
        color=CORAL,
        linestyle="--",
        linewidth=1.7,
        label="Fitted BLS model",
    )
    axis.axvline(0, color=NAVY, linestyle=":", linewidth=1)
    axis.set(
        title=f"Folded signal: TOI-791 {label}",
        xlabel="Hours from recovered transit center",
        ylabel="Detrended flux (ppt)",
        xlim=(-36, 36),
    )
    axis.grid(color=LIGHT, linewidth=0.7, alpha=0.8)
    axis.legend(frameon=False, loc="lower left", fontsize=8)
    if label == "c":
        axis.text(
            0.98,
            0.96,
            f"{len(c_event_windows)} TESS event windows\npublication-informed period ranking",
            transform=axis.transAxes,
            ha="right",
            va="top",
            fontsize=8,
            color=NAVY,
        )


draw_folded_panel(ax_fold_b, "b", recovered_b, TEAL, bls_b)
draw_folded_panel(ax_fold_c, "c", recovered_c, GOLD, bls_c)

figure.suptitle(
    "TOI-791 b recovery and publication-informed c period ranking",
    fontsize=18,
    fontweight="bold",
    color=NAVY,
)
figure.text(
    0.5,
    -0.01,
    "Data: NASA TESS / MAST public 2-minute SPOC PDCSAP products. "
    "Dotted lines mark peer-reviewed periods; BLS models are diagnostic, not physical transit fits.",
    ha="center",
    va="bottom",
    fontsize=9,
    color=NAVY,
)

chart_path = OUTPUT_DIR / "toi791_transit_recovery.png"
figure.savefig(chart_path, dpi=220, bbox_inches="tight", facecolor="white")
if "agg" in plt.get_backend().lower():
    plt.close(figure)
else:
    plt.show()

print(f"Saved chart: {chart_path.relative_to(PROJECT_ROOT)}")
print(f"Saved results: {results_path.relative_to(PROJECT_ROOT)}")
print(f"Saved manifest: {manifest_path.relative_to(PROJECT_ROOT)}")
print(
    "Saved TOI-791 c event inventory: "
    f"{c_events_path.relative_to(PROJECT_ROOT)}"
)
No description has been provided for this image 
Saved chart: outputs/toi791_transit_recovery.png
Saved results: outputs/toi791_recovery_results.csv
Saved manifest: outputs/toi791_tess_product_manifest.csv
Saved TOI-791 c event inventory: outputs/toi791_c_event_windows.csv

Interpretation and limitations¶

  • TOI-791 b: the strongest broad-search BLS peak is approximately 139.31 days, consistent with the paper's 139.29931-day period at the reporting precision used here. This is a numerical peak location, not a formal period measurement with an uncertainty.
  • TOI-791 c: after b is masked, a publication-informed 170–280-day search ranks approximately 232.01 days highest. The selected TESS data contain two c-like event windows separated by about 928 days, or four 232-day cycles. TESS-only folding does not independently eliminate every integer-cycle alias; the paper's additional observations resolve that ambiguity.
  • Folded panels: folding stacks observations separated by one recovered period. For b, multiple events align. For c, the two event windows align at zero phase. The dashed curve is only the fitted BLS approximation; real transits have curved ingress, egress, and limb-darkened bottoms.
  • Do not overread the depths: sector detrending and 30-minute binning attenuate them. The peer-reviewed values, not this notebook's BLS boxes, should be used for physical interpretation.
  • Not a confirmation analysis: BLS cannot rule out eclipsing binaries, nearby contaminants, or every instrumental effect. The study's ground-based and dynamical follow-up is what establishes the planetary interpretation.
  • Not a blinded discovery search: the period ranges and the second-stage masking strategy were chosen with the published system in mind. This is a transparent recovery exercise, not independent alias resolution.

Possible next steps¶

  1. Repeat the recovery across several defensible detrending windows and bin sizes to measure how sensitive the periods, depths, and rankings are to preprocessing choices.
  2. Run injection/recovery tests and odd/even transit checks, then inspect TESS pixel-level centroids to quantify signal completeness and contamination risk.
  3. Fit individual transit times from later sectors and compare them with a linear ephemeris before attempting a full transit-timing-variation model.