11institutetext: Astronomical Institute, Czech Academy of Sciences, Fričova
298, 251 65, Ondřejov, Czech Republic 22institutetext: Center for Astrophysics
${\rm\mid}$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA
33institutetext: Facultad de Ingeniera y Ciencias, Universidad Adolfo Ibáñez,
Av. Diagonal las Torres 2640, Peñalolén, Santiago, Chile 44institutetext:
Millennium Institute for Astrophysics, Chile 55institutetext: Data Observatory
Foundation, Chile 66institutetext: Thühringer Landessternwarte Tautenburg,
Sternwarte 5, 07778 Tautenburg, Germany 77institutetext: Department of
Theoretical Physics and Astrophysics, Faculty of Science, Masaryk University,
Kotlarska 2, CZ-611 37, Brno, Czech Republic 88institutetext: Department of
Electrical Engineering and Centre of Astro-Engineering, Pontificia Universidad
Catolica de Chile, Av. Vicuña Mackenna 4860, Santiago, Chile 99institutetext:
Centro de Astrobiología (CSIC-INTA), ESAC, Camino bajo del castillo s/n, 28692
Villanueva de la Cañada, Madrid, Spain 1010institutetext: NASA Ames Research
Center, Moffett Field, CA 94035, USA 1111institutetext: Department of Physics
and Astronomy, University of New Mexico, 210 Yale Blvd NE, Albuquerque, NM
87106, USA 1212institutetext: Department of Physics and Kavli Institute for
Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA 1313institutetext: Department of Earth, Atmospheric
and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA
02139, USA 1414institutetext: Department of Aeronautics and Astronautics, MIT,
77 Massachusetts Avenue, Cambridge, MA 02139, USA 1515institutetext: SETI
Institute, Mountain View, CA 94043 USA/NASA Ames Research Center, Moffett
Field, CA 94035 USA 1616institutetext: Department of Astrophysical Sciences,
Princeton University, Princeton, NJ 08544, USA

# BD-14 3065b (TOI-4987b): from giant planet to brown dwarf: evidence for
deuterium burning in old age?

Ján Šubjak 1122 0000-0002-5313-9722 David W. Latham 22 0000-0001-9911-7388
Samuel N. Quinn 22 0000-0002-8964-8377 Perry Berlind 22 0009-0005-7108-9502
Michael L. Calkins 22 0000-0002-2830-5661 Gilbert A. Esquerdo 22
0000-0002-9789-5474 Rafael Brahm 334455 0000-0002-9158-7315 Eike Guenther 66
0000-0002-9130-6747 Jan Janík 77 0000-0002-6384-0184 Petr Kabáth 11
0000-0002-1623-5352 Leonardo Vanzi 88 José A. Caballero 99 0000-0002-7349-1387
Jon M. Jenkins 1010 0000-0002-4715-9460 Ismael Mireles 1111
0000-0002-4510-2268 Sara Seager 121213131414 0000-0002-6892-6948 Avi Shporer
1212 0000-0002-1836-3120 Stephanie Striegel 1515 0009-0008-5145-0446 Joshua N.
Winn 1616 0000-0002-4265-047X

(; )

The present study reports the confirmation of BD-14 3065b, a transiting
planet/brown dwarf in a triple-star system, with a mass near the deuterium
burning boundary. BD-14 3065b has the largest radius observed within the
sample of giant planets and brown dwarfs around post-main-sequence stars. Its
orbital period is 4.3 days, and it transits a subgiant F-type star with a mass
of $M_{\star}=1.41\pm 0.05$ M⊙, a radius of $R_{\star}=2.35\pm 0.08$ R⊙, an
effective temperature of $T_{\rm eff}=6935\pm 90$ K, and a metallicity of
$-0.34\pm 0.05$ dex. By combining TESS photometry with high-resolution spectra
acquired with the TRES and Pucheros+ spectrographs, we measured a mass of
$M_{p}=12.37\pm 0.92\,\,{\rm M_{Jup}}$ and a radius of $R_{p}=1.926\pm
0.094\,\,{\rm R_{Jup}}$. Our discussion of potential processes that could be
responsible for the inflated radius led us to conclude that deuterium burning
is a plausible explanation resulting from the heating of BD-14 3065b’s
interior. Detection of the secondary eclipse with TESS photometry enables a
precise determination of the eccentricity $e_{p}=0.066\pm 0.011$ and reveals
BD-14 3065b has a brightness temperature of $3520\pm 130$ K. With its unique
characteristics, BD-14 3065b presents an excellent opportunity to study its
atmosphere through thermal emission spectroscopy.

###### Key Words.:

giant planets – brown dwarfs – techniques: photometric – techniques:
spectroscopic – techniques: radial velocities

##  1 Introduction

The question of distinguishing between planets and brown dwarfs has been
debated at length, with a proposed solution being the use of formation history
as the most fundamental differentiator (Burrows et al., 2001). According to
this approach, substellar objects formed through disk instability similar to
stars are classified as brown dwarfs, while those formed via core accretion
are classified as planets. However, determining the formation history of
individual systems observationally is a difficult task. Some authors used the
lack of known transiting brown dwarfs in the mass range between $35\leq\,{\rm
M_{Jup}}\sin{i}\leq 55$ to motivate the idea that this gap separates two
different BD populations that result from different BD formation processes (Ma
& Ge, 2014). Recent discoveries of new transiting brown dwarfs in this region
have revealed that these two formation mechanisms do not form such a well-
defined mass boundary. Furthermore, it was theorized that disk instability can
produce objects down to a few masses of Jupiter, while core accretion can
produce objects as massive as tens of Jupiter (Mordasini et al., 2012).

Another proposed way to differentiate between planets and brown dwarfs is
predicated on the deuterium reactions occurring within their interior, as
outlined by Saumon et al. (1996), Chabrier & Baraffe (2000), or Burrows et al.
(2001). The mass of an object serves as the primary determinant in this
definition, allowing for straightforward classification in most cases.
However, objects that fall near the dividing line pose difficulties,
particularly given the uncertainties surrounding their mass and chemical
composition. Many authors have adopted the dividing line of $13\,\,{\rm
M_{Jup}}$; nevertheless, every attempt to classify objects close to this
boundary is akin to shooting in the dark. As demonstrated by Spiegel et al.
(2011), or Mollière & Mordasini (2012), deuterium burning is reliant on the
deuterium abundance, helium mass fraction, and metallicity, leading to a
blurring of the dividing line from $11-16\,\,{\rm M_{Jup}}$ depending on the
specific chemical composition.

In their paper, Hatzes & Rauer (2015) presented an argument for defining
planets and brown dwarfs based on observables such as mass and radius. While
objects in the mass range of 0.3 to 80 Jupiter masses follow a consistent
trend in the mass-radius diagram, the authors contend that there is no
compelling reason to classify brown dwarfs differently. An observer might
argue that we should investigate whether objects above and below $13\,\,{\rm
M_{Jup}}$ are distinct. Therefore, objects that straddle this threshold, such
as BD-14 3065b, present an intriguing opportunity for further study. If these
objects can undergo radius re-inflation during a post-main-sequence phase of
the host star due to deuterium burning, they may indeed differ from their
lower-mass counterparts.

A critical inquiry shared by transiting giant planets and brown dwarfs
pertains to their inflated radii, given that many of these objects appear to
have an unexpectedly large radius. Notably, the typical radii of known brown
dwarf companions range from 0.7 to 1.4 Jupiter radii (Šubjak et al., 2020;
Carmichael et al., 2021; Šubjak et al., 2023), while some hot Jupiters are
approaching a radius of 2 Jupiter radii (e.g., Anderson et al., 2011; Fulton
et al., 2015; Almenara et al., 2015; Yee et al., 2023). Numerous studies have
been conducted to investigate hot Jupiters, revealing a correlation between
the incident stellar irradiation and the measured radii of such planets (e.g.,
Enoch et al., 2012; Sestovic et al., 2018). Theoretical studies have explored
the influence of stellar irradiation on the radii of companions, which has
been incorporated into various evolutionary models of giant planets and brown
dwarfs (Baraffe et al., 2003; Fortney et al., 2007; Baraffe et al., 2008). In
addition to this, several other factors may also play important roles in the
inflation of planet radii. These include tidal dissipation (Bodenheimer et
al., 2001, 2003; Arras & Socrates, 2010; Jermyn et al., 2017), kinetic heating
(a fraction of the incoming stellar flux is assumed to be converted into
kinetic energy and then dissipated at the center of the planet) (Guillot &
Showman, 2002), enhanced atmospheric opacities (Burrows et al., 2007), double-
diffusive convection (Chabrier & Baraffe, 2007; Kurokawa & Inutsuka, 2015),
vertical advection of potential temperature (Youdin & Mitchell, 2010; Tremblin
et al., 2017), Ohmic heating through magnetohydrodynamic effects (Batygin &
Stevenson, 2010; Perna et al., 2010; Wu & Lithwick, 2013; Ginzburg & Sari,
2016), or deuterium burning (Mollière & Mordasini, 2012). The mass of giant
planets and brown dwarfs is a key factor in determining the extent of
inflation, with the response decreasing as mass increases.

We have confirmed BD-14 3065b as an inflated giant planet/brown dwarf orbiting
a sub-giant F-type star. Our paper includes a description of the observations
in Section 2, data analysis in Section 3, a discussion in Section 4, and a
summary of our results in Section 5.

##  2 Observations

###  2.1 TESS light curves

TESS monitored BD-14 3065 in three different sectors. Its first observation
was in sector 9, from February 28, 2019, to March 26, 2019, with 1800 s time
sampling. The second observation was in sector 35, from February 09, 2021, to
March 07, 2021, with 600 s time sampling. Lastly, it was observed in sector
62, from February 12, 2023, to March 10, 2023, with 120 s time sampling. The
data for BD-14 3065 are publicly available on the Mikulski Archive for Space
Telescopes
(MAST)111https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html and are
provided by the TESS Science Processing Operations Center (SPOC) located at
NASA Ames Research Center. The transit signature of BD-14 3065b was first
detected by the MIT Quick Look Pipeline (QLP; Huang et al., 2020a, b) and
alerted by the TESS Science Office on January 06, 2022 (Guerrero et al.,
2021).

The light curves (LCs) from sectors 9 and 35 were processed by the QLP
pipeline, using the difference imaging technique. This involves subtracting
each frame from a reference frame generated by combining many comparison
frames. Once the difference image is obtained, the QLP performs aperture
photometry and scales the difference relative to the star’s average flux. This
ensures that any contamination from other stars within the same aperture is
automatically accounted for. The LC from TESS sector 62 processed by the SPOC
pipeline (Jenkins et al., 2016) was downloaded directly from the MAST archive,
using the lightkurve package (Lightkurve Collaboration et al., 2018). The SPOC
pipeline aims to remove systematic errors related to the instrument. The Pixel
Response Functions (PRFs) were used to analyze the crowding, and a crowding
correction was included in the PDCSAP flux time series. Neglecting such a
correction can result in underestimating the planet’s radius. The pipeline
shows that 0.997 of the light in the optimal aperture comes from the target
rather than other stellar sources, indicating insignificant dilution due to
faint background stars. However, the analysis cannot rule out other stellar
sources unresolved by Gaia (see Section 2.2). Moreover, according to Data
Validation’s (Twicken et al., 2018) difference image centroiding analysis, the
source of the transit signature was within $1.0\pm 2.5$ arcsec of the target
star.

We used the python package called citlalicue (Barragán et al., 2022) to
correct for systematics and eliminate stellar variability in the PDCSAP LCs
(Smith et al., 2012; Stumpe et al., 2012, 2014) obtained through lightkurve.
citlalicue combines a Gaussian Process regression with transit models computed
via the pytransit code (Parviainen, 2015) to generate a model comprising both
the light curve variability and the transits. The resulting model is then used
to eliminate the variability and produce a flattened LC containing only the
transit photometric variations. Fig. 1 displays the LCs before and after this
process. We detected 16 transits, with six each in Sectors 9 and 62 and four
in Sector 35. Additionally, the SPOC data validation report has reported the
secondary eclipse signature in sector 62 with a depth of $553\pm 183$ ppm. We
validate this signal and discuss the transit and secondary eclipse modeling in
Section 3.4.

We used a tool called tpfplotter (Aller et al., 2020) to compare the Gaia DR2
catalog with the TESS target pixel file (tpf). It allowed us to detect any
potential stars that could dilute the TESS photometry, with a limit of up to 6
magnitudes difference. The resulting tpf image can be viewed in Fig. 2. We
found no additional sources in the TESS aperture that could dilute TESS LCs.
It is worth noting that the Gaia DR3 catalog doesn’t contain any additional
bright sources within the TESS aperture. Therefore, the use of its older
release in our analysis is justified.

The SPOC LCs revealed additional stellar variability in the form of
pulsations. To analyze these pulsations, we utilized maximum-likelihood (MLP)
and generalized Lomb-Scargle (GLS) periodograms, which detected two signals
representing pulsations with periods of 0.35 and 0.39 days, consistent across
all three sectors. The GLS periodogram is shown in Figure 3. We subtracted
these signals by fitting two simple sinusoidal functions to the data. As a
result, the final LCs appeared flat, without any signs of stellar variability.

Table 1: Basic parameters for BD-14 3065.

Parameter Description Value Source $\alpha_{\rm J2000}$ Right Ascension (RA)
10 14 40.73 1 $\delta_{\rm J2000}$ Declination (Dec) -15 38 34.33 1 $V_{T}$
Tycho-2 $V_{T}$ mag $11.146\pm 0.098$ 2 $G$ Gaia $G$ mag $10.921\pm 0.002$ 1
$T$ TESS $T$ mag $10.631\pm 0.008$ 3 $J$ 2MASS $J$ mag $10.143\pm 0.026$ 4
$\mu_{\alpha}\cos{\delta}$ PM in RA (mas/yr) $-4.990\pm 0.074$ 1
$\mu_{\rm\delta}$ PM in DEC (mas/yr) $7.629\pm 0.062$ 1 $\pi$ Parallax (mas)
$1.765\pm 0.055$ 1 Other identifiers: TOI-4987 TIC 293607057 TYC 6059-01368-1
2MASS J10144071-1538344 Gaia DR2 3751877374435102720

References: 1 - $Gaia$ DR3, Gaia Collaboration et al. (2021);
2 - Tycho-2, Høg et al. (2000); 3 - TESS, Stassun et al. (2018);
4 - 2MASS, Skrutskie et al. (2006)

Figure 1: SPOC PDCSAP LCs for BD-14 3065 from TESS sector 9 (top panel),
sector 35 (middle panel), and sector 62 (bottom panel). Grey points represent
the TESS observations, while red lines correspond to the out-of-transit GP
models created with citlalicue to capture the variability in the LCs. Datasets
were divided by these models leading to flattened TESS LCs (blue points) with
transit models (orange lines).

Green triangles indicate the positions of transits. The LCs still show
pulsations, but these were later removed.

###  2.2 Contamination from nearby sources

As part of follow-up observations coordinated by the TESS Follow-up Observing
Program (TFOP) High-Resolution Imaging Sub-Group 3 (SG3), we obtained high
angular-resolution imaging of BD-14 3065. Observations were made on April 15,
2022, using the High-Resolution Camera (HRCam; Tokovinin & Cantarutti, 2008)
speckle interferometry instrument on the Southern Astrophysical Research
(SOAR) 4.1m telescope. The observation strategy and data reduction procedures
are described in Tokovinin (2018), or Ziegler et al. (2021). Observations were
made in the Cousins $I$ filter with a resolution of 36 mas. The final
reconstructed image, shown in Fig. 4, reaches a contrast of $\Delta$$mag=5.1$
at a separation of 0.5 arcsec and has an estimated PSF that is 0.064 arcsec
wide. The speckle imaging of BD-14 3065 detected a companion at an angular
separation of 0.92 arcsec from the primary (translated into a projected
physical separation of 520 au), which has $\Delta$$I=2.3$ mag, and position
angle of 210.5 or 30.5 degrees. We note that Gaia did not report this
companion. However, the Gaia Renormalised Unit Weight Error (RUWE) value of
3.5 for BD-14 3065 suggests the presence of a close multiplicity. We discuss
the treatment of the companion in Section 3.1.

###  2.3 Ground-based light curves

As part of the TFOP, ground-based photometric data were collected for BD-14
3065. The observations were scheduled using Transit Finder, a customized
version of the Tapir software (Jensen, 2013), and the photometric data was
extracted using AstroImageJ (Collins et al., 2017). However, the light curve
precision is not sufficient to improve the parameters from the TESS LCs.
Moreover, with some precision, different dilution parameters would need to be
adopted as these observations were made in different filters. Therefore, it
was decided not to consider ground-based photometry in this paper further.

Figure 2: Gaia DR2 catalog overplotted to the TESS TPF image. Figure 3: The GLS periodogram of SPOC LCs after substruction of the 4.3-d rotational signal using Gaussian processes. Two dominant peaks represent stellar pulsations. Figure 4: High-resolution imaging data for BD-14 3065 in the $I$ band. SOAR HRCam speckle sensitivity curve (solid line) and autocorrelation function (ACF, an inside image). The detected companion appears as two dots in the ACF image. Table 2: Multi-order relative radial velocities of BD-14 3065 from TRES and Pucheros+. $\rm BJD_{\rm TDB}$ | RV (m s-1) |  $\sigma_{\rm RV}$ (m s-1) | Instrument
---|---|---|---
2459601.934527 | 5017 | 144 | TRES
2459605.955459 | 4746 | 111 | TRES
2459616.890452 | 7113 | 144 | TRES
2459617.956610 | 5910 | 162 | TRES
2459618.882581 | 4777 | 155 | TRES
2459619.825369 | 4944 | 101 | TRES
2459620.899334 | 6977 | 143 | TRES
2459621.881599 | 6721 | 116 | TRES
2459622.928761 | 4436 | 87 | TRES
2459623.911069 | 4917 | 112 | TRES
2459624.901062 | 6643 | 105 | TRES
2459625.944000 | 6830 | 125 | TRES
2459626.800265 | 5017 | 107 | TRES
2459635.850469 | 4864 | 168 | TRES
2459922.954247 | 1358 | 172 | TRES
2459923.986662 | 257 | 173 | TRES
2459924.995943 | 2389 | 156 | TRES
2459930.979474 | 2002 | 113 | TRES
2459932.988929 | 252 | 93 | TRES
2459933.984636 | 2123 | 141 | TRES
2459934.972431 | 2169 | 156 | TRES
2459936.000038 | 32 | 102 | TRES
2459936.965439 | -380 | 122 | TRES
2459944.924846 | 37 | 160 | TRES
2459948.961852 | 67 | 107 | TRES
2459957.965266 | -70 | 152 | TRES
2459963.949222 | 2463 | 205 | TRES
2460013.776225 | 57 | 153 | TRES
2460014.767013 | 1006 | 119 | TRES
2460016.749179 | 1564 | 113 | TRES
2460017.788237 | -82 | 94 | TRES
2460028.772600 | 1471 | 129 | TRES
2460029.856313 | 1056 | 127 | TRES
2460030.758805 | 32 | 145 | TRES
2460071.660621 | 1984 | 196 | TRES
2460072.653619 | 1507 | 133 | TRES
2460073.691053 | 457 | 129 | TRES
2460074.681585 | 962 | 140 | TRES
2460075.670057 | 2148 | 118 | TRES
2460076.660340 | 2017 | 142 | TRES
2460077.697400 | 721 | 198 | TRES
2460077.486053 | 2324 | 1122 | Pucheros+
2460078.557373 | 715 | 1114 | Pucheros+
2460078.580683 | 2164 | 1455 | Pucheros+
2460080.569803 | 3522 | 973 | Pucheros+
2460080.590706 | 3201 | 1857 | Pucheros+
2460085.523472 | 3111 | 867 | Pucheros+
2460085.544838 | 4610 | 901 | Pucheros+
2460086.569549 | 1530 | 1006 | Pucheros+
2460086.590926 | 2332 | 1117 | Pucheros+
2460088.507060 | 5548 | 1197 | Pucheros+
2460088.528426 | 3275 | 1204 | Pucheros+
2460089.519780 | 4380 | 1107 | Pucheros+
2460094.491713 | 3257 | 794 | Pucheros+
2460096.531829 | 3333 | 374 | Pucheros+
2460102.509028 | 4203 | 551 | Pucheros+

###  2.4 TRES spectra

We obtained spectra of BD-14 3065 using the Tillinghast Reflector Echelle
Spectrograph (TRES) on the 1.5-m telescope at the Whipple Observatory on Mount
Hopkins, Arizona, from January 22, 2022 to May 13, 2023. The spectrograph has
a resolving power of R $\approx$ 44 000 and covers wavelengths ranging from
390 nm to 910 nm. A total of 41 spectra of BD-14 3065 were captured with TRES,
with exposure times varying from 350–3600 s depending on the weather
conditions. The S/N ratio per resolution element of these spectra ranged from
24 to 49. We derived the relative RVs from the TRES spectra by considering
multiple echelle orders from each spectrum, which were cross-correlated with a
high SNR median spectrum produced using the complete set of observed spectra.
We have also adjusted the spectra to account for the small drifts in the TRES
zero point over the period of the observations. The TRES fiber diameter
measures 2.4 arcseconds and contamination from a stellar companion around
BD-14 3065, as resolved by speckle imaging, may be expected. The exact level
of contamination also depends on the seeing. The final list of radial
velocities (RVs) can be found in Table 2.

###  2.5 Spectroscopy from the PLATOSpec project

We acquired additional spectra from the PLATOSpec project using the fully
operational telescope E152 at La Silla, Chile. An interim spectrograph,
PUCHEROS+, has been installed for scientific use. PUCHEROS+ is an improved
version of a PUCHEROS spectrograph (Vanzi et al., 2012), with a resolving
power of R $\approx 18,000$. The Andor iKon M detector provides a wavelength
range of over 400 to 700 nm. Typically, the spectrograph’s radial velocity
stability is about 100 m s-1 over a month for cool main sequence stars with
V$=8$ mag. PUCHEROS+ uses the ceres+ pipeline (Brahm et al., 2017), which
produces wavelength calibrated and optimally extracted order-by-order spectra.
The radial velocities are obtained using the cross-correlation method from the
$\mathrm{ceres+}$ pipeline. The PUCHEROS+ fiber diameter has 2.0 arcseconds,
and the expected contamination from a stellar companion around BD-14 3065 is
smaller than for the TRES spectra. While the data will be available in the ESO
archive later, they are currently accessible upon request from the PLATOSpec
team. The target was observed over the course of 16 nights between May and
June 2023, with each exposure lasting 1800 seconds.

##  3 Analysis

###  3.1 Modeling Stellar Parameters

We analyzed the spectra of BD-14 3065 to determine the spectral parameters of
the star. Then, we used the catalog photometry from Gaia, 2MASS, Tycho-2,
Galex, and WISE, along with the $\Delta$$mag$ measured by SOAR HRCam, to model
blended two-component and three-component models of the SED using the
isochrones package (Morton, 2015; Dotter, 2016; Choi et al., 2016). The two-
component model accounts for the stellar companion resolved in the speckle
imaging, while the three-component model additionally accounts for the long-
term signal detected in the radial velocity observations. Further elaboration
on this matter is presented below. We assumed that the stars are physically
associated and have the same parallax as measured for the primary by Gaia, as
is the case for most stars separated by $\leq 1\arcsec$ (Horch et al., 2014;
Matson et al., 2018).

Table 3: Properties of the primary star and the stellar companions. The upper
panel describes the results from the two-component model, with the final
column listing the observed photometry. The lower panel describes the results
from the three-component model.

Parameter BD-14 3065 A BD-14 3065 B BD-14 3065 AB $M_{\star}$ ($\rm\,{\rm
M_{\odot}}$) $1.41\pm 0.05$ $0.99\pm 0.03$ - $R_{\star}$ ($\rm\,{\rm
R_{\odot}}$) $2.35\pm 0.08$ $0.96\pm 0.03$ - $\log{g}$ $3.85\pm 0.03$ $4.47\pm
0.02$ - $T_{\rm eff}$ (K) $6935\pm 90$ $6190\pm 150$ - $[{\rm Fe/H}]$
$-0.34\pm 0.05$ $-0.37\pm 0.05$ - $v_{\rm rot}\sin{i_{\star}}$ (km s-1) $22\pm
3$ - - $P_{\rm rot}$ (days) $4.3\pm 0.2$ - - log L ($\,{\rm L_{\odot}}$)
$1.06\pm 0.03$ $0.09\pm 0.05$ - log Age (Gyr) $9.35\pm 0.05$ $9.35\pm 0.05$ -
$A_{V}$ (mag) $0.19\pm 0.01$ $0.19\pm 0.01$ - Galex $FUV$ mag $18.62\pm 0.06$
$23.67\pm 0.53$ $18.69\pm 0.09$ Galex $NUV$ mag $14.49\pm 0.03$ $17.70\pm
0.23$ $14.49\pm 0.01$ Tycho-2 $B_{T}$ mag $11.51\pm 0.02$ $14.19\pm 0.12$
$11.56\pm 0.09$ Gaia $BP$ mag $11.22\pm 0.02$ $13.73\pm 0.10$ $11.046\pm
0.003$ Tycho-2 $V_{T}$ mag $11.13\pm 0.02$ $13.62\pm 0.10$ $11.15\pm 0.10$
Gaia $G$ mag $11.02\pm 0.02$ $13.44\pm 0.10$ $10.921\pm 0.003$ Gaia $RP$ mag
$10.67\pm 0.02$ $12.98\pm 0.09$ $10.510\pm 0.004$ TESS $T$ mag $10.66\pm 0.02$
$12.97\pm 0.09$ $10.631\pm 0.008$ 2MASS $J$ mag $10.30\pm 0.02$ $12.47\pm
0.08$ $10.14\pm 0.03$ 2MASS $H$ mag $10.10\pm 0.02$ $12.17\pm 0.07$ $9.98\pm
0.03$ 2MASS $K_{S}$ mag $10.07\pm 0.02$ $12.13\pm 0.07$ $9.93\pm 0.03$ WISE1
mag $10.05\pm 0.02$ $12.11\pm 0.07$ $9.84\pm 0.02$ WISE2 mag $10.05\pm 0.02$
$12.11\pm 0.07$ $9.87\pm 0.02$ WISE3 mag $10.03\pm 0.02$ $12.09\pm 0.07$
$9.85\pm 0.05$ Parameter BD-14 3065 A BD-14 3065c BD-14 3065 B $M_{\star}$
($\rm\,{\rm M_{\odot}}$) $1.40\pm 0.06$ $1.02\pm 0.03$ $0.98\pm 0.02$
$R_{\star}$ ($\rm\,{\rm R_{\odot}}$) $2.18\pm 0.08$ $1.02\pm 0.05$ $0.96\pm
0.03$ $\log{g}$ $3.91\pm 0.03$ $4.43\pm 0.03$ $4.47\pm 0.02$ $T_{\rm eff}$ (K)
$6985\pm 90$ $6300\pm 150$ $6175\pm 120$ $[{\rm Fe/H}]$ $-0.34\pm 0.04$
$-0.39\pm 0.05$ $-0.37\pm 0.04$ $v_{\rm rot}\sin{i_{\star}}$ (km s-1) $22\pm
3$ - - $P_{\rm rot}$ (days) $4.3\pm 0.2$ - - log L ($\,{\rm L_{\odot}}$)
$1.00\pm 0.03$ $0.16\pm 0.07$ $0.08\pm 0.05$ log Age (Gyr) $9.36\pm 0.06$
$9.36\pm 0.06$ $9.36\pm 0.06$ $A_{V}$ (mag) $0.19\pm 0.01$ $0.19\pm 0.01$
$0.19\pm 0.01$ Galex $FUV$ mag $18.63\pm 0.07$ $22.92\pm 0.67$ $23.77\pm 0.49$
Galex $NUV$ mag $14.57\pm 0.04$ $17.33\pm 0.33$ $17.75\pm 0.22$ Tycho-2
$B_{T}$ mag $11.66\pm 0.04$ $13.97\pm 0.20$ $14.22\pm 0.12$ Gaia $BP$ mag
$11.37\pm 0.04$ $13.55\pm 0.18$ $13.76\pm 0.10$ Tycho-2 $V_{T}$ mag $11.28\pm
0.04$ $13.44\pm 0.17$ $13.65\pm 0.10$ Gaia $G$ mag $11.17\pm 0.04$ $13.27\pm
0.16$ $13.47\pm 0.10$ Gaia $RP$ mag $10.83\pm 0.04$ $12.83\pm 0.15$ $13.01\pm
0.09$ TESS $T$ mag $10.82\pm 0.04$ $12.82\pm 0.15$ $13.01\pm 0.09$ 2MASS $J$
mag $10.47\pm 0.04$ $12.34\pm 0.14$ $12.50\pm 0.08$ 2MASS $H$ mag $10.28\pm
0.04$ $12.05\pm 0.13$ $12.20\pm 0.07$ 2MASS $K_{S}$ mag $10.25\pm 0.04$
$12.02\pm 0.12$ $12.16\pm 0.07$ WISE1 mag $10.23\pm 0.04$ $12.00\pm 0.12$
$12.14\pm 0.07$ WISE2 mag $10.23\pm 0.04$ $11.99\pm 0.12$ $12.14\pm 0.07$
WISE3 mag $10.21\pm 0.04$ $11.97\pm 0.12$ $12.12\pm 0.07$

Figure 5: Comparison of the observed spectrum (blue line) with the G-type star
synthetic spectrum (green line) shifted by offset. The orange spectrum
represents the scaled green spectrum with respect to the brightness of the
blended stellar companion. The red spectrum represents the spectrum of the
primary star.

####  3.1.1 iSpec and isochrones Stellar Parameters

To obtain the spectral parameters of the host stars, we used the iSpec
framework (Blanco-Cuaresma et al., 2014; Blanco-Cuaresma, 2019), which
incorporates the SME radiative transfer code and the MARCS atmosphere models,
together with version 5 of the GES atomic line list. We co-added all TRES
spectra to obtain the high S/N final spectrum. This spectrum provided a higher
S/N ratio than the PUCHEROS+ spectrum and was therefore used to determine
stellar parameters, including the effective temperature $T_{\rm eff}$,
metallicity [Fe/H], surface gravity $\log{g}$, and the projected stellar
equatorial velocity $v\sin{i}$. We used a spectral fitting technique that
minimized the chi-squared value between the calculated synthetic and observed
spectrum in the interval from 490 to 560 nm.

The broad-band catalog and time-series photometry of BD-14 3065 need to be
corrected for the contamination of the second star found in the speckle
imaging. Hence, we used the python package called isochrones to create a two-
component model of the SED, taking into account the priors set on the derived
spectral parameters, available catalog photometry, $\Delta$$I$ mag difference
from speckle imaging, and Gaia parallax. We also considered the line-of-sight
extinction by placing an upper limit derived from the Bayestar dust maps of
Green et al. (2019) incorporated into the dustmaps code (Green, 2018). We have
obtained the stellar parameters of both stars reported in Table 3.

We performed a test to determine if the G-type companion star, discovered
through speckle imaging, impacts our spectral parameters. To do this, we
synthesized the spectrum of a G-type star with a $\Delta$$I$ of 2.3 mag. In
the wavelength range of $490-560$ nm, the $\Delta$mag is 2.5, indicating that
the companion star is ten times fainter. However, since a G-type star has
deeper spectral lines compared to the F-type primary star, the companion star
can still affect our findings. We also considered that a G-type star has a
slower rotation rate. Figure 5 shows the spectra of both stars. Our findings
revealed that the contamination could cause an underestimation of the
effective temperature by 200 K and an overestimation of the value of $\log{g}$
by 0.1 dex. It is consistent with the work of El-Badry et al. (2018), who
studied the implications of unresolved binaries in stellar spectra for
spectral fitting. They presented a figure that describes the effects of
unresolved binarity on the accuracy of stellar labels recovered from spectral
fitting as a function of the mass ratio between primary and secondary.
Therefore, we adjusted our priors placed on the derived spectral parameters in
SED modeling to take this effect into account. We plot the best-fit SED model
in the Fig. 6, and present the stellar properties of the primary and secondary
components in Table 3. We also overplot in Fig. 7 luminosities and effective
temperatures of both stars with the MIST stellar evolutionary tracks (Choi et
al., 2016). According to the analysis of the stellar isochrones, the system
appears to have an age of approximately $2.25$ Gyr. In order to corroborate
this finding, we searched for the lithium $670.8$ nm feature in the co-added
TRES spectrum, yet did not detect any indication of it, further supporting the
system’s old age. Additionally, the synchronized rotation period of the star,
found in Section 3.1.2, with the planet’s orbital period suggests tidal
evolution, which also aligns with the old age of the system. Finally, the
PyAstronomy python package was utilized to determine the $U$,$V$,$W$
velocities for BD-14 3065, using the
gal_uvw222https://pyastronomy.readthedocs.io/en/latest/pyaslDoc/aslDoc/gal_uvw.html
function. Based on our analysis, we found that BD-14 3065 is associated with
the thin disk (e.g., Leggett, 1992; Bensby et al., 2003). This outcome is
predictable since the thick disk is mainly composed of stars that are at least
8 Gyr old (e.g., Fuhrmann, 1998).

Figure 6: Two-component SED fit for the BD-14 3065 system. The blue points
with vertical error bars show the observed catalog fluxes and uncertainties
from Galex, Gaia, TESS, Tycho-2, 2MASS, and WISE, while the horizontal
errorbars illustrate the width of the photometric band. The red points show
the model fluxes of the primary and secondary, while the yellow points show
the combined flux of both stars. For illustrative purposes, we overplot
extinction-corrected Phoenix atmospheric models (Husser et al., 2013) for the
two stellar components, although these models were not used directly in the
fit, which was based on MIST bolometric correction tables. Figure 7:
Luminosity vs. effective temperature plot for the two-component model. Curves
represent MIST isochrones for ages: 2.0 Gyr (cyan), 2.8 Gyr (purple), 3.2 Gyr
(yellow), and for [Fe/H] = -0.5. Red points represent the positions of
TOI-1268 A and TOI-1268 B with their errorbars.

In Section 3.4 we show an observed long-period trend in radial velocities that
could be attributed to the presence of an additional star even closer than 0.9
arcsec. We can rule out that the trend is linear; however, we cannot determine
its period as it remains unconstrained. There are two possible scenarios: the
star may be too faint to cause significant dilution, or it could be more
massive and blended in catalog photometry. To account for the second option,
we modeled a three-component spectral energy distribution with the primary
star being an unresolved binary. Similarly to the two-component SED model, the
input parameters for this model included the derived spectral parameters of
the primary, available catalog photometry, $\Delta$$I$ mag difference from
speckle imaging, and Gaia parallax. However, in this case, the $\Delta$ mag
difference is between the star located at 0.9 arcsec and the blended binary
star at a smaller separation. The best-fit SED model is shown in Figure 8,
along with the stellar properties of all three components in Table 3. To
provide further context, we also overplot stars’ luminosities and effective
temperatures with the MIST stellar evolutionary tracks in Figure 9.

####  3.1.2 Stellar rotation

The rotation period of the star was established using the SPOC TESS PDCSAP
LCs. It is important to note that the long trends (see Fig. 1) and transits
have been removed prior to analysis. To achieve this, we made use of the
Gaussian Process (GP) Regression library, Celerite. In Foreman-Mackey et al.
(2017), the authors provide a detailed explanation of the library, including
the interpretation of various kernels. To account for the variations in LCs
due to heterogeneous surface features such as spots and plages, we utilized a
rotational kernel function defined as:

Figure 8: Three-component SED fit for the BD-14 3065 system. The blue points
with vertical error bars show the observed catalog fluxes and uncertainties
from Galex, Gaia, TESS, Tycho-2, 2MASS, and WISE, while the horizontal
errorbars illustrate the width of the photometric band. The red points show
the model fluxes of the three components, while the yellow points show the
combined flux of three stars. For illustrative purposes, we overplot
extinction-corrected Phoenix atmospheric models for the three stellar
components, although these models were not used directly in the fit, which was
based on MIST bolometric correction tables. Figure 9: Luminosity vs. effective
temperature plot for the three-component model. Curves represent MIST
isochrones for ages: 2.0 Gyr (cyan), 2.8 Gyr (purple), 3.2 Gyr (yellow), and
for [Fe/H] = -0.5. Red points represent the positions of all three stars in
the system with their errorbars. Figure 10: Top: TESS SPOC PDCSAP LCs (black
points) with the MAP model prediction. It is important to note that the long
trends and transits have been removed prior to analysis. The red line shows
the predictive mean, and the red contours show the predictive standard
deviation. Bottom left: probability density of the rotation period. The period
is the parameter $P_{rot}$ in Equation 1. The 1$\sigma$ error bar is indicated
with the dashed black lines. Bottom right: GLS periodogram for all three
sectors combined.

$k(\tau)={\frac{A}{2+B}e^{-\tau/L}\left[cos\left(\frac{2\pi\tau}{P_{rot}}\right)+(1+B)\right]}\,,$
(1)

where variables $A$ and $B$ define the amplitude of the GP, while $\tau$
represents the time-lag and $L$ denotes the amplitude-modulation timescale of
the GP. $P_{rot}$ indicates the rotational period. The maximum a posteriori
(MAP) parameters are estimated using the L-BFGS-B nonlinear optimization
routine (Byrd et al., 1995; Zhu et al., 1997). To obtain the marginalized
posterior distributions of free parameters, we run MCMC using emcee (Goodman &
Weare, 2010; Foreman-Mackey et al., 2013). The analysis results in a derived
rotational period of $P_{rot}=4.3\pm 0.2$ days, which is shown in Fig. 10 as
the probability density of $P_{rot}$ alongside the MAP model prediction.
Additionally, we apply the generalized Lomb-Scargle (GLS) periodograms
(Zechmeister & Kürster, 2009) to the TESS LCs as an independent check of the
derived $P_{rot}$. We observe strong peaks around 4.3 days in individual
sectors and a forest of peaks around 4.3 days when using all sectors together.
We consider this to be the rotation period, and the periodogram is shown in
Fig. 10. Table 3 reports the derived rotation period. By using combination of
$P_{rot}$, $v\sin{i}$, and $R_{\star}$ we can constrain stellar inclination
angle, which is $i_{\star}\,=\,59_{-12}^{+17}$ for the two-star model and
$i_{\star}\,=\,65_{-14}^{+16}$ for the three-star model. Our approach follows
the procedure described in Masuda & Winn (2020) and used in Bowler et al.
(2023). This method properly accounts for the correlation between stellar
equatorial velocity and projected rotational velocity. Posterior distributions
of stellar inclination for both models can be seen in Fig. 11. It should be
noted that the used value of $v\sin{i}$ represents an upper limit, as the
spectral contamination of BD-14 1065 can lead to an overestimation of the
$v\sin{i}$ value by a few km s-1. Correction for this overestimation would
lead to an even higher probability of spin-orbit misalignment.

Figure 11: Posterior distributions for line-of-sight stellar inclination in
the form of cos i⋆ for the two-star (left) and three-star (right) models. Blue
shadow regions show 1-$\sigma$ and 2-$\sigma$ credible intervals.

###  3.2 System architecture

In order to model the parameters of the transiting object, we need to
distinguish if it orbits the primary F-type star or, rather, the wide G-type
companion. A critical study that shows the significance of such analyses is
presented in the work of Mandushev et al. (2005). The GSC01944-02289 system
bears a strong resemblance to the architecture of the BD-14 3065 system, with
an F-type primary star and a wide G-type companion. The research revealed that
the potential brown dwarf around the F star is, in reality, an M-dwarf star
transiting the G star.

The high brightness ratio of two stars posed significant challenges to
revealing the true nature of the BD-14 3065 system. To do so, we adopted the
approach detailed in Mandushev et al. (2005). Initially, we computed the
bisector spans (BIS) of cross-correlation functions (CCFs) to examine the
asymmetries of spectral lines. We then determined the periodogram of the BIS
time series (see Section 3.3), which revealed two peaks at 4.3 and 3.3 days.
These peaks indicated that two CCFs are moving with respect to each other with
the period of orbital motion, causing asymmetries in the blended CCF. However,
here, we observed a key difference between the GSC01944-02289 and BD-14 3065
systems. While BIS showed a semiamplitude of more than one
$\mathrm{km\,s^{-1}}$ in the case of GSC01944-02289, BD-14 3065 BIS only
showed a semiamplitude of less than 300 $\mathrm{m\,s^{-1}}$, and the peak was
barely significant, albeit visible.

Subsequently, we performed the two-dimensional cross-correlation algorithm
TODCOR (Zucker & Mazeh, 1994). TODCOR cross-correlates the observed spectra
against a composite template obtained by adding together two distinct
templates chosen to match each star. We used stellar parameters from the two-
star SED modeling and assumed that the stars in the potential eclipsing binary
would likely have their rotation synchronized with the orbital motion due to
tidal forces (translating to the vsini of 11 $\mathrm{km\,s^{-1}}$). We were
able to measure RVs of both stars. An illustration of this detection is given
in Fig. 12, which shows the two-dimensional correlation function as well as
cross-sections of the two-dimensional correlation function at a fixed velocity
for the G star and F star, respectively. The periodogram of the F-type star
(see Section 3.3) confirmed the substellar nature of the transiting object.

Figure 12: Upper panels: Two-dimensional cross-correlation surface computed
with TODCOR as a function of the F and G star velocities for one spectrum. The
white lines represent the adopted velocities for both stars. Bottom panels:
Cross-sections of the two-dimensional correlation functions at a fixed
velocity for the G star and F star, respectively.

To further validate our findings, we conducted additional tests. Firstly, we
allowed the light fraction factor to be a free parameter, and the results
indicated that the best value was close to the expected value of 0.1 from
speckle imaging. Secondly, we recomputed the periodogram of BIS from the
cross-section of 2D CCF of the F-type star. This revealed a significant
reduction in the scatter of data, which was now on the same level as observed
in the Sun. It also eliminated the peak at the period of orbital motion,
suggesting that we had successfully extracted the contamination of the G-type
star. Finally, we compared the photometry of GSC01944-02289 and BD-14 3065. In
the TESS LCs of GSC01944-02289, clear variations with the period of half the
orbital period of the eclipsing star are visible, with maxima in phases 0.25
and 0.75. The stellar deformation of close binary stars causes these
variations. We would, therefore, expect to observe similar variations for
BD-14 3065. However, the LCs instead show variations with a period similar to
the transiting companion’s period and a clear evolution of this variability
through TESS sectors, which we interpret as the rotation period of the F-type
star (see Fig. 1). Such a rotation period would be an improbable coincidence
if the companion orbited a different star and no tidal interactions were in
play. In conclusion, the transiting object transits an F-type primary star
rather than a wide stellar companion.

###  3.3 Frequency Analysis

We performed a frequency analysis to detect all signals in RVs. Our
observations had a baseline of around 500 days, providing a frequency
resolution of approximately $0.002d^{-1}$. We observed a non-linear long trend
in the dataset, with an unconstrained period, as seen in the periodogram.
After fitting the trend, we identified a significant period of 4.3 days and
its one-day alias, consistent with the transits from TESS photometry (see Fig.
13).

To test the effect of the unconstrained period/mass of the close stellar
companion BD-14 3065c on the parameters of the transiting planet BD-14 3065b,
we gradually increased the mass of the companion from a low-mass M dwarf star
to a G-type star (the result from a blended three-component SED model). We
initialized the MCMC run for each mass using exostriker (Trifonov, 2019) to
derive posterior distributions of the transiting planet parameters. We
conclude that these parameters are well-constrained and consistent within
their error bars, regardless of the mass of the close stellar companion.

Figure 13: Generalized Lomb-Scargle periodograms of TRES RVs of BD-14 3065:
(a) TRES RVs, (b) TRES RVs minus long trend, (c) TRES BIS CCF, (d) TRES BIS
CCF from TODCOR, (e) TRES RVs of the F star from TODCOR, (f) TRES RVs of the G
star from TODCOR. The vertical green lines highlight the orbital period of the
transiting companion and its 1-day alias. Horizontal dashed lines show the
theoretical FAP levels of 10%, 1%, and 0.1% for each panel.

###  3.4 Joint RVs and transits modelling

The process of modeling transits and radial velocities was carried out
simultaneously using the allesfitter package. The fitting routine, as
described in Günther & Daylan (2021), is based on several packages, including
ellc for light-curve and RV models (Maxted, 2016), aflare for flare modeling
(Davenport et al., 2014), dynesty for static and dynamic Nested Sampling
(Speagle, 2020), emcee for MCMC sampling (Foreman-Mackey et al., 2013), and
celerite for GP models (Foreman-Mackey et al., 2017). The allesfitter package
offers a range of orbital and transit/eclipse models, accommodating multiple
exoplanets, multi-star systems, transit-timing variations, phase curves,
stellar variability, starspots, stellar flares, and various systematic noise
models, including Gaussian processes.

To obtain the optimal conjunction time, period, inclination, eccentricity, and
argument of periastron, which were later used as priors for joint modeling, we
performed a fit of photometry alone, as the secondary eclipse is visible in
TESS photometry (see Fig. 14). The LCs were already corrected for variability
using the procedure outlined in Section 2.1, and we employed the quadratic
limb darkening law with coefficients q1 and q2. To ensure accuracy and avoid
degeneracy between these coefficients and transit depth, we adopted Gaussian
priors within $0.1\,\sigma$ around the values from Claret (2017), which
presents limb-darkening coefficients specifically computed for the TESS
photometry for a wide range of effective temperatures, gravities,
metallicities, and microturbulent velocities. Although modeling the pulsations
simultaneously with the transits may be preferred, it will not impact the
results since pulsations are much slower than the transit duration. We
conducted two separate fits based on the stellar parameters and dilution
factors obtained from two-star and three-star SED modeling. Table 4 shows the
priors and posteriors of fitted parameters, while Table 5 presents the transit
parameters and physical parameters derived from Tables 3 and 4, respectively.
The results are plotted in Fig. 14.

The stellar radius is an important parameter that leads to the companion
radius computed through the transit depth. We need to double-check that the
radius of BD-14 3065b is correct, as the discussion of this work is built upon
its value. As a first check, we modeled the 2-star SED, assuming that the star
found in the high-resolution imaging is an unassociated companion with its own
age and distance. However, such a model led to similar parameters of the
primary star ($M_{\star}=1.38\pm 0.03$ M⊙, a radius of $R_{\star}=2.40\pm
0.08$ R⊙). As a second check, to ensure the reliability of reported stellar
spectral parameters, which are crucial in our SED modeling, we conducted a
search for values in the relevant literature. The Sixth Data Release of the
Radial Velocity Experiment (RAVE) reported stellar parameters of $T_{\rm
eff}=6809\pm 140$ K and $\log{g}=3.78\pm 0.19$. Gaia reported stellar
parameters of $T_{\rm eff}=6791\pm 32$ K and $\log{g}=3.90\pm 0.01$ using the
photometry333https://gea.esac.esa.int/archive/documentation/GDR3/Data_analysis/chap_cu8par/sec_cu8par_apsis/ssec_cu8par_apsis_gspphot.html.
These values are consistent with our own findings. Additionally, it is
important to note that unresolved binarity is not expected to impact these
values significantly (El-Badry et al., 2018). We are, therefore, confident
that the reported stellar parameters are reliable. Finally, the impact of an
unresolved binary star on the radius of a transiting planet has been the
subject of much discussion in the literature. The ratio of the planet’s true
radius $R_{P,\,true}$, to the observed radius $R_{P,\,observed}$ when assuming
a single star with no companions, can be described as:

$X_{R}\equiv{\frac{R_{P,\,true}}{R_{P,\,observed}}}\equiv\left(\frac{R_{\star,\,t}}{R_{\star,\,1}}\right)\,\sqrt{\frac{F_{total}}{F_{t}}}$
(2)

where $R_{\star,\,1}$ represents the radius of the primary star (assumed to be
a single star), and $F_{t}$ and $R_{\star,\,t}$ represent the brightness and
radius, respectively, of the star that is being transited by the planet. When
the true stellar radius is determined, unaffected by stellar companions, for
example, through asteroseismic analysis or spectroscopic analysis together
with stellar isochrones, we have $R_{\star,\,1}\equiv R_{\star,\,t}$ for a
planet transiting the primary star (Hirsch et al., 2017). In our study, we
derived stellar parameters from SED modeling, and
${R_{\star,\,1}}>{R_{\star,\,t}}$. However, assuming that our effective
temperature is not largely underestimated,
$\frac{R_{\star,\,1}}{R_{\star,\,t}}\sim\sqrt{\frac{F_{total}}{F_{t}}}$, and
we do not anticipate an overestimation of the radius of BD-14 3065b.

Figure 14: The phased LC of the BD-14 3065, fitted with the allesfitter as
described in Section 3.4. The grey points represent TESS data, and the blue
points are TESS binned data in the phase curve with a bin width of 15 minutes
in the phase. The red line represents the best transit and eclipse model.
Figure 15: The radial velocity curve of BD-14 3065, fitted with the
allesfitter as described in Section 3.4. The blue points represent TRES and
Pucheros+ data, and the red line represents the best model that results from
the same global fit of both datasets together. For illustration purposes, we
have plotted the shortest feasible circular orbit for BD-14 3065c. Hence, we
have plotted a two-star model that assumes that the stellar companion BD-14
3065c does not cause significant dilution. We have subtracted a value of
2,457,000 from BJD time. Table 4: Fitted parameters from the allesfitter
analysis.

Parameter | Value | Unit | fit/fixed
---|---|---|---
Fitted parameters: 2-star model
$R_{b}/R_{\star}$ | $0.0843_{-0.0022}^{+0.0023}$ |  | fit
$(R_{\star}+R_{b})/a_{b}$ | $0.1804_{-0.0028}^{+0.0032}$ |  | fit
$\cos{i_{b}}$ | $0.1602_{-0.0033}^{+0.0036}$ |  | fit
$T_{0;b}$ | $2313.48578\pm 0.00076$ | $\mathrm{BJD}$ | fit
$P_{b}$ | $4.2889731\pm 0.0000052$ | $\mathrm{d}$ | fit
$K_{b}$ | $1.210\pm 0.071$ | $\mathrm{km\,s^{-1}}$ | fit
$\sqrt{e_{b}}\cos{\omega_{b}}$ | $0.152\pm 0.013$ |  | fit
$\sqrt{e_{b}}\sin{\omega_{b}}$ | $0.208_{-0.033}^{+0.024}$ |  | fit
$D_{\mathrm{0;TESS}}$ | $0.106\pm 0.010$ |  | fixed
$q_{1;\mathrm{TESS}}$ | $0.111_{-0.110}^{+0.073}$ |  | fit
$q_{2;\mathrm{TESS}}$ | $0.209\pm 0.093$ |  | fit
$J_{b;\mathrm{TESS}}$ | $0.176_{-0.054}^{+0.070}$ |  | fit
${\sigma_{\mathrm{TESS}}}$ | $1332\pm 11$ | ${\mathrm{ppm}}$ | fit
${\sigma_{\mathrm{jitter;TRES}}}$ | $0.32\pm 0.05$ | ${\mathrm{km\,s^{-1}}}$ | fit
Fitted parameters: 3-star model
$R_{b}/R_{\star}$ | $0.0901\pm 0.0022$ |  | fit
$(R_{\star}+R_{b})/a_{b}$ | $0.1711_{-0.0029}^{+0.0031}$ |  | fit
$\cos{i_{b}}$ | $0.1495\pm 0.0038$ |  | fit
$T_{0;b}$ | $2313.48571\pm 0.00078$ | $\mathrm{BJD}$ | fit
$P_{b}$ | $4.2889734\pm 0.0000051$ | $\mathrm{d}$ | fit
$K_{b}$ | $1.200\pm 0.071$ | $\mathrm{km\,s^{-1}}$ | fit
$\sqrt{e_{b}}\cos{\omega_{b}}$ | $0.148_{-0.013}^{+0.014}$ |  | fit
$\sqrt{e_{b}}\sin{\omega_{b}}$ | $0.220_{-0.036}^{+0.025}$ |  | fit
$D_{\mathrm{0;TESS}}$ | $0.226\pm 0.10$ |  | fixed
$q_{1;\mathrm{TESS}}$ | $0.143_{-0.088}^{+0.098}$ |  | fit
$q_{2;\mathrm{TESS}}$ | $0.210\pm 0.099$ |  | fit
$J_{b;\mathrm{TESS}}$ | $0.160_{-0.049}^{+0.067}$ |  | fit
${\sigma_{\mathrm{TESS}}}$ | $1330\pm 11$ | ${\mathrm{ppm}}$ | fit
${\sigma_{\mathrm{jitter;TRES}}}$ | $0.32\pm 0.05$ | ${\mathrm{km\,s^{-1}}}$ | fit

Comments: a – BJD - 2,457,000, b – dilution factor, c – surface brightness
ratio between the companion and star

Table 5: Derived parameters from the allesfitter analysis.

Parameter Value Derived parameters: 2-star model Host radius over semi-major
axis b; $R_{\star}/a_{\mathrm{b}}$ $0.1664_{-0.0032}^{+0.0037}$ Semi-major
axis b over host radius; $a_{\mathrm{b}}/R_{\star}$ $6.01_{-0.13}^{+0.12}$
Companion radius b over semi-major axis b; $R_{\mathrm{b}}/a_{\mathrm{b}}$
$0.01403_{-0.00047}^{+0.00060}$ Companion radius b; $R_{\mathrm{b}}$
($\mathrm{R_{\oplus}}$) $21.59\pm 1.05$ Companion radius b; $R_{\mathrm{b}}$
($\mathrm{R_{jup}}$) $1.926\pm 0.094$ Semi-major axis b; $a_{\mathrm{b}}$
($\mathrm{R_{\odot}}$) $14.11\pm 0.57$ Semi-major axis b; $a_{\mathrm{b}}$
(AU) $0.0656\pm 0.0026$ Inclination b; $i_{\mathrm{b}}$ (deg)
$80.78_{-0.30}^{+0.24}$ Eccentricity b; $e_{\mathrm{b}}$
$0.0660_{-0.011}^{+0.0098}$ Argument of periastron b; $w_{\mathrm{b}}$ (deg)
$54.0_{-9.1}^{+6.0}$ Mass ratio b; $q_{\mathrm{b}}$
$0.00838_{-0.00052}^{+0.00055}$ Companion mass b; $M_{\mathrm{b}}$
($\mathrm{M_{\oplus}}$) $3932_{-280}^{+290}$ Companion mass b;
$M_{\mathrm{b}}$ ($\mathrm{M_{jup}}$) $12.37_{-0.87}^{+0.92}$ Companion mass
b; $M_{\mathrm{b}}$ ($\mathrm{M_{\odot}}$) $0.01181_{-0.00083}^{+0.00088}$
Impact parameter b; $b_{\mathrm{tra;b}}$ $0.9099_{-0.0079}^{+0.0092}$ Total
transit duration b; $T_{\mathrm{tot;b}}$ (h) $3.086\pm 0.041$ Full-transit
duration b; $T_{\mathrm{full;b}}$ (h) $0.61_{-0.22}^{+0.19}$ Epoch occultation
b; $T_{\mathrm{0;occ;b}}$ $2315.7361_{-0.0074}^{+0.0084}$ Impact parameter
occultation b; $b_{\mathrm{occ;b}}$ $1.013_{-0.024}^{+0.020}$ Host density
from orbit b; $\rho_{\mathrm{\star;b}}$ (cgs) $0.223_{-0.014}^{+0.013}$
Companion density b; $\rho_{\mathrm{b}}$ (cgs) $2.15_{-0.48}^{+0.53}$
Companion surface gravity b; $g_{\mathrm{b}}$ (cgs) $10536\pm 950$ Equilibrium
temperature b; $T_{\mathrm{eq;b}}$ (K) $2001\pm 31$ Transit depth (undil.) b;
$\delta_{\mathrm{tr;undil;b;TESS}}$ (ppt) $6.49_{-0.20}^{+0.23}$ Transit depth
(dil.) b; $\delta_{\mathrm{tr;dil;b;TESS}}$ (ppt) $5.80_{-0.11}^{+0.14}$
Occultation depth (undil.) b; $\delta_{\mathrm{occ;undil;b;TESS}}$ (ppt)
$0.482_{-0.11}^{+0.091}$ Occultation depth (dil.) b;
$\delta_{\mathrm{occ;dil;b;TESS}}$ (ppt) $0.431_{-0.096}^{+0.081}$ Derived
parameters: 3-star model Host radius over semi-major axis b;
$R_{\star}/a_{\mathrm{b}}$ $0.1570_{-0.0027}^{+0.0028}$ Semi-major axis b over
host radius; $a_{\mathrm{b}}/R_{\star}$ $6.37\pm 0.11$ Companion radius b over
semi-major axis b; $R_{\mathrm{b}}/a_{\mathrm{b}}$
$0.01414_{-0.00045}^{+0.00047}$ Companion radius b; $R_{\mathrm{b}}$
($\mathrm{R_{\oplus}}$) $21.41\pm 1.05$ Companion radius b; $R_{\mathrm{b}}$
($\mathrm{R_{jup}}$) $1.910\pm 0.094$ Semi-major axis b; $a_{\mathrm{b}}$
($\mathrm{R_{\odot}}$) $13.88\pm 0.57$ Semi-major axis b; $a_{\mathrm{b}}$
(AU) $0.0646\pm 0.0026$ Inclination b; $i_{\mathrm{b}}$ (deg) $81.40\pm 0.22$
Eccentricity b; $e_{\mathrm{b}}$ $0.0702_{-0.011}^{+0.0096}$ Argument of
periastron b; $w_{\mathrm{b}}$ (deg) $56.2_{-7.2}^{+4.7}$ Mass ratio b;
$q_{\mathrm{b}}$ $0.00826\pm 0.00056$ Companion mass b; $M_{\mathrm{b}}$
($\mathrm{M_{\oplus}}$) $3850_{-290}^{+300}$ Companion mass b;
$M_{\mathrm{b}}$ ($\mathrm{M_{jup}}$) $12.11_{-0.91}^{+0.95}$ Companion mass
b; $M_{\mathrm{b}}$ ($\mathrm{M_{\odot}}$) $0.01156_{-0.00087}^{+0.00091}$
Impact parameter b; $b_{\mathrm{tra;b}}$ $0.8959\pm 0.0067$ Total transit
duration b; $T_{\mathrm{tot;b}}$ (h) $3.052\pm 0.042$ Full-transit duration b;
$T_{\mathrm{full;b}}$ (h) $0.78_{-0.17}^{+0.14}$ Epoch occultation b;
$T_{\mathrm{0;occ;b}}$ $2315.7365\pm 0.0070$ Impact parameter occultation b;
$b_{\mathrm{occ;b}}$ $1.007_{-0.024}^{+0.019}$ Host density from orbit b;
$\rho_{\mathrm{\star;b}}$ (cgs) $0.266\pm 0.014$ Companion density b;
$\rho_{\mathrm{b}}$ (cgs) $2.16_{-0.31}^{+0.37}$ Companion surface gravity b;
$g_{\mathrm{b}}$ (cgs) $10270\pm 850$ Equilibrium temperature b;
$T_{\mathrm{eq;b}}$ (K) $1957\pm 28$ Transit depth (undil.) b;
$\delta_{\mathrm{tr;undil;b;TESS}}$ (ppt) $7.42_{-0.12}^{+0.13}$ Transit depth
(dil.) b; $\delta_{\mathrm{tr;dil;b;TESS}}$ (ppt) $5.741_{-0.092}^{+0.10}$
Occultation depth (undil.) b; $\delta_{\mathrm{occ;undil;b;TESS}}$ (ppt)
$0.57\pm 0.12$ Occultation depth (dil.) b; $\delta_{\mathrm{occ;dil;b;TESS}}$
(ppt) $0.442\pm 0.089$

Comments: a – computed assuming zero albedo

##  4 Discussion

###  4.1 Mass-radius diagram

Combining data from the TESS mission and ground-based spectroscopy confirmed
the substellar nature of a $P=4.29$ days candidate around the F-type subgiant
BD-14 3065. We found that the physical parameters of BD-14 3065b
($M_{p}=12.37\pm 0.92\,\,{\rm M_{Jup}}$, $R_{p}=1.926\pm 0.094\,\,{\rm
R_{Jup}}$) are consistent with an object at the boundary between giant planets
and brown dwarfs based on the deuterium burning criterium. A comparison of the
mass and radius of BD-14 3065b with the population of well-characterized brown
dwarfs and giant planets with a mass greater than $5\,\,{\rm M_{Jup}}$ is
presented in Fig. 16. Three subplots are created, with colors representing
different parameters: age, orbital period, and equilibrium temperature. The
position of BD-14 3065b on the mass-radius diagram is exceptional. It is one
of two objects with a significantly larger radius compared to the rest of the
population. The second object is the brown dwarf RIK 72b (David et al., 2019)
in the Upper Scorpius OB association, with an age of approximately 5 Myrs. RIK
72b transits the pre-main sequence M-dwarf star with a period of $\sim 98$
days, and thus, most of its internal thermal energy originates from its
contraction. The authors demonstrate that current evolutionary models that
predict the size and luminosity of brown dwarfs as a function of age are
consistent with the derived parameters. Conversely, BD-14 3065 is an old star
without lithium in the spectrum. Old age of about 2.3 Gyr was derived using
stellar isochrones for the subgiant primary star. Explaining the observed
radius for such an old age is challenging using substellar evolutionary
models, even when considering stellar irradiation, which slows down
contraction.

The research conducted by Baraffe et al. (2003) delved into the influence of
stellar irradiation on the radius of giant planets with a maximum mass of
$10\,\,{\rm M_{Jup}}$. It is noteworthy that BD-14 3065b is slightly more
massive. However, the effect of irradiation decreases as the mass of the
planet increases, given a certain incident flux. According to their models,
for an object with a mass of about $10\,\,{\rm M_{Jup}}$ and an age similar to
that of BD-14 3065, irradiation would inflate its radius by approximately
$0.1\,\,{\rm R_{Jup}}$. It is important to note that they hypothesized a lower
stellar temperature of 6000 K but also assumed a smaller irradiated planet
separation of 0.046 au. Even when considering the effects of stellar
irradiation, the theoretical models can only match the radius of BD-14 3065b
for a few million years. To demonstrate this conclusion, we utilized the
thermal evolution model outlined by Fortney et al. (2007), which also
considers the impact of stellar irradiation. We assumed a core mass of
$10\,M_{\oplus}$ and adjusted the model to BD-14 3065’s luminosity and BD-14
3065b’s orbital distance. We found only a small inflation of the radius;
hence, the irradiation cannot replicate the observed radius of the companion.

In Fig. 16, in the plot where colors represent object equilibrium temperature
for a zero albedo, we can see that objects with larger equilibrium
temperatures typically tend to have a higher radius. This suggests that the
radius in many systems is affected by stellar irradiation. Equilibrium
temperature refers to the theoretical temperature an object reaches when it is
in thermal equilibrium, assuming it is only being heated by its parent star
and following black body re-radiation.

We can compare BD-14 3065b with two systems hosting the same mass companions:
HATS-70 (Zhou et al., 2019) and XO-3 (Johns-Krull et al., 2008). HATS-70b is
an object with a mass of $12.9\,\,{\rm M_{Jup}}$ that transits a metallic-line
A-type star with a period of 1.9 days. Due to its proximity to the star, it
receives a higher amount of stellar radiation than BD-14 3065b, resulting in
an equilibrium temperature of 2750 K. While the authors suggest that the
stellar irradiation may account for the observed inflated radius of
$1.38\,\,{\rm R_{Jup}}$, they note that other well-documented effects may also
play a role, such as enhanced atmospheric metallicity and opacity (e.g.,
Burrows et al., 2007, 2011) or ohmic dissipation (e.g., Batygin & Stevenson,
2010). XO-3b is an object with a mass of $11.8\,\,{\rm M_{Jup}}$ that transits
an F-type star with a period of 3.2 days. Dang et al. (2022) employed Bayesian
analysis to determine the amount of additional heating required to explain its
radius of $1.3\,\,{\rm R_{Jup}}$ as a fraction of the incident flux. The
results indicate that approximately 20% of additional heating from an external
source is necessary, and the authors discussed potential sources such as tidal
heating or deuterium burning. Additionally, Yang & Wei (2022) identified
evidence of XO-3b’s tidal evolution through transit timing variations,
suggesting that tidal heating could be a source of additional heating.

The inflated radii of many hot Jupiters have been attributed to several
mechanisms, the intensity of which is often scaled to explain the observed
radii. However, even among the population of inflated gas giants and brown
dwarfs, BD-14 3065b stands out with its unusually large radius for such a
massive and old object. The potential reasons behind the unusually strong
intensity of these mechanisms for BD-14 3065b, or whether there is an unusual
combination of several mechanisms at play, are beyond the scope of this work.
Instead, we focus on the two most obvious differences between BD-14 3065b and
other systems in the literature. Firstly, BD-14 3065b is situated at the
boundary of deuterium burning, and secondly, its primary star is an evolved
subgiant.

Figure 16: The population of known companions within the mass interval
$5-70\,\,{\rm M_{Jup}}$. The position of BD-14 3065b is highlighted together
with the three objects mentioned in the discussion. The left subplot is
colored with respect to the age of the systems, the middle one with respect to
the orbital period of the planets, and the right one with respect to the
equilibrium temperature of the planets.

###  4.2 Deuterium burning

It is generally accepted that objects with masses below $13\,\,{\rm M_{Jup}}$
are incapable of burning deuterium or can only do so to a small extent and are
therefore classified as planets. Mollière & Mordasini (2012) developed a model
to determine the initial deuterium burned based on object mass, helium
fraction, metallicity, gas accretion rate, and deuterium abundance. According
to their modeling, BD-14 3065b may still retain a significant portion of its
deuterium despite the system’s age. Hence, deuterium represents an energetic
reservoir inside the object. Additionally, their model suggests that low-mass
objects, such as BD-14 3065b, continue to burn a small fraction of deuterium
over several Gyr, following a strong burning phase that lasts approximately
100 Myr. The rate of deuterium reaction is highly dependent on the interior
temperature, as expressed by the equation
$\epsilon\,\propto\,[D/H]\,{\rho}\,T^{11.8}$. Therefore, increasing the
object’s interior temperature through tidal interactions with an evolving
post-main-sequence star could lead to an increase in the rate of deuterium
burning.

Vissapragada et al. (2022) observed evidence of tidally-driven inspiral of the
Kepler-1658b planet around an evolved subgiant star, similar to BD-14 3065,
with an orbital period of 3.85 days, which is not far from BD-14 3065b’s
period. Furthermore, using a secondary eclipse from Kepler data, they suggest
that the planet is tidally superheated, with approximately ten percent of the
energy from the shrinking orbit being dissipated in the planet. WASP-12b is
another planet that has been observed to experience orbital decay (Bailey &
Goodman, 2019; Yee et al., 2020). The planet is orbiting an F-type star, which
could be a subgiant. The planet has a highly inflated radius, and its
brightness temperature has been determined to be $3640\pm 230$ K through
secondary eclipse observations (Crossfield et al., 2012). This temperature is
notably higher than the planet’s equilibrium temperature of $2990\pm 110$ K
(Crossfield et al., 2012). This suggests that the planet may be experiencing
tidal superheating, similar to Kepler-1658b.

Tidal evolution should eventually lead to planet engulfment if the orbital
angular momentum is too low to allow the star to achieve spin-orbit
synchronization (Hut, 1980). However, using TESS photometry, we have
determined that BD-14 3065b is synchronized with the rotation of its host
star. Hence, the current tidally-driven inspiral of the planet is not
expected. To confirm this, we studied the transit timings of BD-14 3065b with
juliet (Espinoza et al., 2019) and found no evidence of significant evolution.
We note that the time sampling of 1800 s in TESS sector 9 and 600 s in sector
35 results in large uncertainties in the derived transit middle times. TTVs
are shown in Fig. 17. It is also important to note that the radius of
Kepler-1658b is not comparable to that of BD-14 3065b.

Figure 17: Transit timing variations of BD-14 3065b.

In order to investigate the possibility of overheating in BD-14 3065b, we
derived its thermal emission through secondary eclipse depth analysis. The
depth is proportional to the square of the radius ratio and the temperature
ratio of the host and transiting companion. We analyzed the TESS LCs and found
a secondary eclipse depth of $\delta_{S}\,=\,0.5\pm 0.1$ parts per thousand
(ppt). The secondary eclipse fit can be observed in Fig. 14. It is important
to note that the TESS bandpass
window444https://archive.stsci.edu/files/live/sites/mast/files/home/missions-
and-data/active-missions/tess/_documents/TESS_Instrument_Handbook_v0.1.pdf is
centered at roughly 800 nm and spans $600-1000$ nm, which means that the
relative brightness of the planet to the star is only measured at these
wavelengths. We used this relative brightness to determine the effective
temperature of BD-14 3065b, employing the equation used for the secondary
eclipse depth from Charbonneau et al. (2005):

$\delta_{S}\approx\bigg{(}\frac{R_{b}}{R_{\star}}\bigg{)}^{2}\frac{\int
F_{b}(\lambda)\,S(\lambda)\,\lambda\,d\lambda}{\int
F_{\star}(\lambda)\,S(\lambda)\,\lambda\,d\lambda}$ (3)

where $S(\lambda)$ is the TESS spectral response function and $F(\lambda)$ is
the flux of companion or star. In order to derive BD-14 3065b’s effective
temperature, we simplify the equation by assuming blackbody profiles to
describe fluxes in the TESS bandpass window and zero geometric albedo. Our
findings indicate that BD-14 3065b has an effective blackbody temperature of
$T_{eff}\,=\,3520\pm 130\,K$. This temperature is notably higher than the
equilibrium temperature that would result from star heating alone.
Consequently, the secondary eclipse has revealed that the planet is
experiencing overheating, which would increase the deuterium burning rate. The
validity of this conclusion is contingent upon the assumption of BD-14 3065b’s
blackbody spectrum. However, it remains uncertain whether the eclipse spectra
of the hottest atmospheres align with blackbody curves. Studies of individual
planets with HST and Spitzer have provided observational evidence in support
of proximity to blackbody (e.g., Kreidberg et al., 2018; Mansfield et al.,
2018; Arcangeli et al., 2019). Nevertheless, recent research has highlighted a
rise in brightness for wavelengths shorter than 1.3 $\mu$m due to the combined
opacities of H-, TiO, VO, and FeH (e.g., Changeat et al., 2022; Coulombe et
al., 2023). On the other hand, according to observations from Coulombe et al.
(2023), TESS observation of planet thermal emission aligns well with the black
body spectrum. Consequently, we are of the opinion that our assumption is
reasonable.

Most models of deuterium burning consider the formation of objects through
gravitational instability. Since such objects are both large and hot,
deuterium burning can halt their contraction, which makes it challenging to
evaluate the impact of deuterium burning on compact objects like BD-14 3065b.
Mollière & Mordasini (2012) modeled the effect of deuterium burning on objects
formed through core accretion with a cold start, where gas is accreted with
low entropy. Although this formation model is theoretical, it allows us to
create an analogy with BD-14 3065b. This model shows that deuterium burning in
dense, compact objects would result in an increase in their radius,
luminosity, and effective temperature. Definitive proof would require the
modeling of the interior of BD-14 3065b, taking into account all of the
effects mentioned above. It is a difficult task that is beyond the scope of
this work.

The process of heating a planet’s atmosphere can result in its expansion, but
such a phenomenon does not fully account for the radii of many of the larger
hot Jupiters. It is highly probable that planets with the largest observed
radii have inflated interiors. The subsequent section will expound on a
potential solution to the inflated radius of BD-14 3065b by examining the
effects of an overheated interior followed by increased deuterium burning.

###  4.3 Overheating of BD-14 3065b

It is a well-established fact that the increase in entropy due to the interior
heating of a planet leads to an increase in its radius (Arras & Bildsten,
2006; Marleau & Cumming, 2014). In this regard, we aim to investigate the
possibility of BD-14 3065b being heated due to the interaction and dissipation
between the star and the planet. Analysis of the TESS photometry data reveals
that the stellar spin frequency is in sync with the planet’s orbital mean
motion, indicating that the synchronization timescale is comparable to, or
smaller than, the star’s age. The synchronization timescale for the planet’s
spin angular frequency is smaller than the orbital synchronization timescale
(Bodenheimer et al., 2001). Therefore, the synchronization process of planet
rotation cannot contribute significantly to the heating of the planetary
interior.

Another potential heating source is the dissipation of a residual
eccentricity, as indicated by a slightly eccentric orbit observed from a
secondary eclipse. The rate of eccentric orbital energy dissipation inside the
planet is proportional to the value of eccentricity, and this rate would
decrease over time as the orbit becomes more circularized. The rate also
depends on the circularization timescale or, more specifically, on the
planet’s and star’s tidal quality factors, and while the circularization
timescale is several orders smaller as stars become sub-giant, most of the
eccentric orbital energy is dissipated inside the star rather than the planet.
Therefore, while the planet may indeed be heated from a residual eccentricity,
this source of heating weakens over time as the orbit becomes less eccentric.
As this heating source was stronger in the past and affects many eccentric
systems in literature, there is no obvious connection to BD-14 3065’s large
radius. To calculate the rate of dissipation of the planet’s eccentric orbital
energy within the planet’s interior, we used an equation from Bodenheimer et
al. (2001). Assuming the planet’s tidal quality factor $log\,Q^{\prime}\approx
5.0$, we computed the value of $\dot{E}_{\rm D}\approx 3.6.10^{27}$ erg s-1.

Other possible sources could be dissipative processes powered by the stellar
irradiation flux. These processes have received significant attention, yet the
mechanism to transport a portion of the incident flux into the planet’s
interior remains an open question, and several processes have been proposed.
One proposed mechanism is atmospheric circulation, which leads to the thermal
dissipation of kinetic energy into the planet’s interior (Guillot & Showman,
2002; Showman & Guillot, 2002). Another proposed mechanism is ohmic
dissipation, which is based on the idea that fast winds generated by the
irradiation drive through the planet’s magnetic fields to create currents that
dissipate ohmically in the interior (e.g., Batygin & Stevenson, 2010; Huang &
Cumming, 2012; Wu & Lithwick, 2013; Ginzburg & Sari, 2016). The connection
between dissipative processes and the evolutionary stage of the subgiant BD-14
3065, which recently increased its radiation, makes these processes appear
very realistic.

To compare the rate of energy transfer to the planet’s interior to the one
coming from an eccentric orbit and determine if the observed radius of BD-14
3065b could be matched for a realistic fraction of stellar irradiation
transported into the planet’s interior, we used the hierarchical Bayesian code
developed by Sarkis et al. (2021). The code utilizes the planetary evolution
model completo21 (Mordasini et al., 2012) and fully non-gray atmospheric
models from the petitCODE (Mollière et al., 2015) to connect the interior
structure models to the observed physical properties such as the planet’s
radius and mass, and the stellar luminosity and semimajor axis. The code
assumes the planet is in a steady state and thus does not calculate the
planetary thermal evolution. The planet models assume no central core, and the
heavy elements (made up entirely of water) are equally distributed throughout
the interior. A study by Sarkis et al. (2021) found that the heating
efficiency increases with increasing equilibrium temperature and reaches a
maximum of 2.5% at around 1860 K, beyond which the efficiency decreases. When
assuming the highest heating efficiency of 2.5%, we computed the value of
$\dot{E}_{\rm T}\approx 2.3\times 10^{28}$ erg s-1, which is an order of
magnitude higher than the energy dissipated from an eccentric orbit. However,
this amount of energy can produce a radius only of about $1.2\,\,{\rm
R_{Jup}}$, assuming the stellar fraction of heavy element mass. Even if we
assume a heating efficiency of 20%, we can barely reach $1.6\,\,{\rm
R_{Jup}}$. There is no known process that can transfer such a high fraction of
stellar irradiation into the interior of the planet. Therefore, we conclude
that while several anticipated processes may increase the internal temperature
of the planet, they would need to be unreasonably efficient to account for the
observed radius. However, such processes would also significantly increase the
rate of deuterium reactions. This appears to be the most plausible scenario.
Alternatively, we can speculate that the planet may have experienced a recent
collision that increased its internal temperature.

Figure 18: Values of the tidal circularisation timescales as a function of
planet tidal quality factor for BD-14 3065b. Different colors represent the
dependence on stellar tidal quality factors. Figure 19: Simulated emission
spectrum of BD-14 3065b as a function of wavelength. The model was generated
by picaso. The strongest atmospheric features are marked.

###  4.4 Tidal interaction

Assuming that the observed eccentricity is residual from formation, an
analysis of the circularization timescale can provide valuable insight into
the star and planet tidal quality factors. Specifically, one may compare the
circularization timescale with the age of the system. It is important to note
that the circularization timescale represents the exponential decay of
eccentricity rather than the time required to reach zero eccentricity.

As noted by Jackson et al. (2008), the timescale required for orbital
circularization of a close-in companion can be calculated using the following
equation:

$\frac{1}{\tau_{e}}=\bigg{[}\frac{63}{4}\sqrt{GM_{\star}^{3}}\frac{R_{\mathrm{p}}^{5}}{Q^{\prime}_{\mathrm{p}}M_{\mathrm{p}}}+\frac{171}{16}\sqrt{G/M_{\star}}\frac{R_{\star}^{5}M_{\mathrm{p}}}{Q^{\prime}_{\star}}\bigg{]}a^{-\frac{13}{2}},$
(4)

where $Q^{\prime}_{\star}$ and $Q^{\prime}_{\mathrm{p}}$ are the tidal quality
factors of the star and planet, respectively. The tidal quality factor
represents the ratio of the elastic energy stored to the energy dissipated per
cycle of the tide, and larger values of the factor lead to longer
circularization timescales. Equation 4 further incorporates the Love number,
the correction factor for the tidal-effective rigidity of the body and its
radial density distribution, resulting in the expression
$Q^{\prime}=3Q/2k_{2}$ (Goldreich & Soter, 1966; Jackson et al., 2008). For a
homogeneous fluid body, $k_{2}=3/2$.

Measuring tidal quality factors is a challenging task. The number of brown
dwarfs with published constraints on Q’ is limited, and the available data
indicate $log\,Q^{\prime}_{\rm BD}\approx 4.15-4.5$ (Heller et al., 2010;
Beatty et al., 2018). Numerous studies have been conducted to derive tidal
quality factors for hot Jupiters, revealing typical values of
$log\,Q^{\prime}_{\rm J}\approx 5.0-6.5$ (e.g., Jackson et al., 2008; Lainey
et al., 2009; Quinn et al., 2014). There is a lack of consensus in the
literature concerning the values of $Q^{\prime}_{\star}$, which may range from
$10^{4}$ to $10^{8}$. Furthermore, a single system’s Q value can change over
time due to the tidal forcing changes caused by both the system’s orbital
evolution and the evolution of the star (e.g., Jackson et al., 2008; Penev et
al., 2012).

In Figure 18, we display the tidal circularisation timescale as a function of
the tidal quality factors of the planet and star. It can be observed that
within the typical range of $log\,Q^{\prime}_{\rm J}\approx 5.0-6.5$, as
reported in the literature, the circularisation timescale is only larger than
system’s age for the highest values of $Q^{\prime}_{\star}$. This outcome
aligns with the equilibrium tide values reported in the literature (Collier
Cameron & Jardine, 2018), which should provide a reasonable estimate for a
massive star like BD-14 3065 (Šubjak et al., 2020). The tidal quality factor
of subgiant stars, as measured by Vissapragada et al. (2022), is approximately
$log\,Q^{\prime}_{\star}\approx 4$, and for such a value, the circularisation
timescale would decrease to a few million years.

###  4.5 Follow-up prospects

BD-14 3065b presents itself as a highly promising target for further
atmospheric characterization through secondary eclipse observations. The
planet’s equilibrium temperature classifies it as an ultra-hot object that is
comparable to WASP-18b, whose thermal emission spectrum was previously
observed using the James Webb Space Telescope (Coulombe et al., 2023). The
most prominent features that can be observed in such an atmosphere include
H2O, CO, FeH, TiO, and VO absorption. A simple atmospheric model of the
anticipated thermal emission spectrum in the NIRSpec and NIRISS wavelength
ranges was created (see Fig. 19) using the picaso code (Batalha et al., 2019)
alongside the Sonora brown dwarf evolutionary models (Marley et al., 2021). In
a similar fashion to WASP-18b, the James Webb Space Telescope would be able to
determine the metallicity and carbon-to-oxygen (C/O) ratio of BD-14 3065b.

C/O ratios have been measured for several transiting planets (e.g.,
Madhusudhan et al., 2011; Moses et al., 2013; Line et al., 2014; Arcangeli et
al., 2018; Ahrer et al., 2023). Beatty et al. (2018) conducted a study on the
brown dwarf CWW 89Ab transiting a 3 Gyr old Sun-like star. The team used the
Spitzer/IRAC to observe the secondary eclipses, and they discovered a
significant $9.3$-$\sigma$ overluminosity of CWW 89Ab compared to the
predicted values from the evolutionary models. The authors discussed the
potential causes of this overluminosity, ruling out both stellar irradiation
and tidal heating as likely explanations. They suggested that a dayside
temperature inversion, resulting from a super-stellar C/O ratio, could be
responsible for this phenomenon. If this is indeed the case, it would imply
that CWW 89Ab is a $36.5\,\,{\rm M_{Jup}}$ object that formed via the core
accretion process. The measurement of the C/O ratio would be essential to gain
a better understanding of this system. It is also worth noting that, in
contrast to BD-14 3065b, CWW 89Ab does not exhibit an inflated radius. We may
also mention a recent study by Yan et al. (2023) in which authors used ground-
based high-resolution spectrograph CRIRES+ to observe dayside thermal emission
spectra of two ultra-hot Jupiters and detected strong CO emission lines in
both planets. BD-14 3065b is an ideal target for such a study.

The C/O ratio is an important characteristic proposed to constrain the
formation of giant planets (Öberg et al., 2011). Generally, a stellar C/O
ratio is expected for companions that form through gravitational instabilities
in a cloud, where all materials mix, and planets that form within the snow
lines. A superstellar C/O ratio could indicate that the atmosphere formed from
disk gas outside the water snowline, possibly via core accretion. It was later
discussed that such simple division is not realistic and that several complex
and interconnected processes affect the planetary formation, such as disc
elemental composition, chemical evolution of the protoplanetary disc,
planetary formation and migration, atmospheric evolution, and others (see
Mollière et al. (2022) for more details). It is still a long way off, if even
possible at all, to connect atmospheric abundances with formation processes
for individual systems. However, we can use the observed atmospheric
composition to study the importance of various formation aspects in an attempt
to reproduce observations. Furthermore, the atmospheric abundance constraints
obtained with new and upcoming instruments will be crucial in developing
planet formation models. Instruments such as GRAVITY, CRIRES+, and recent and
upcoming facilities such as JWST, ARIEL, and the ELT will obtain atmosphere
composition for many planets.

BD-14 3065b would be the most massive transiting companion with a measured C/O
ratio to date, contributing to our understanding of giant planets and brown
dwarfs. Furthermore, while WASP-18b, which orbits much closer to its star, has
an atmosphere that is heated by stellar irradiation, resulting in an inverted
temperature-pressure profile, BD-14 3065b is expected to have an additional
internal source. Hence, studying the temperature-pressure profile of its
atmosphere can assist in understanding the processes involved in atmosphere
heating.

##  5 Summary

We have presented an analysis of the substellar object BD-14 3065b transiting
a subgiant F-type star. BD-14 3065b has an orbital period of P = 4.2889731 ±
0.0000052 days, and the estimated age of the system is about 2 Gyr, making
BD-14 3065b an extremely inflated object for its mass and age. This study
focused on the two most notable differences between BD-14 3065b and other
systems in the literature, attempting to link them to the observed radius.
Firstly, BD-14 3065b is located at the deuterium burning boundary, and
secondly, its primary star is an evolved subgiant. The conclusion drawn is
that the radius could be explained by the enhanced rate of deuterium burning
caused by sufficient heating of BD-14 3065b’s interior. Finally, detection of
the secondary eclipse with TESS photometry enables a precise determination of
the eccentricity and reveals BD-14 3065b has a brightness temperature of
$3520\pm 130$ K. Hence, BD-14 3065b is considered an excellent target for
future atmospheric characterization through secondary eclipse observations.
Such an investigation would provide valuable insights into the processes
responsible for the planet’s atmospheric heating, as well as its formation
history, by comparing its C/O ratio to that of its host star. This would be of
significant value, given that BD-14 3065b sits at the giant planet/brown dwarf
boundary.

###### Acknowledgements.

J.Š. would like to acknowledge the support from GACR grant 23-06384O. P.K. is
grateful for the support from grant LTT-20015. L.V. acknowledges ANID projects
Fondecyt n. 1211162, QUIMAL ASTRO20-0025, and BASAL FB210003. R.B.
acknowledges support from FONDECYT Project 11200751 and additional support
from ANID – Millennium Science Initiative – ICN12_009. J.A.C. acknowledges
financial support from the Spanish Agencia Estatal de Investigación
(AEI/10.13039/501100011033) of the Ministerio de Ciencia e Innovación and the
European Regional Development Fund “A way of making Europe” through project
PID2022-137241NB-C42. Funding for the TESS mission is provided by NASA’s
Science Mission Directorate. We acknowledge the use of public TESS data from
pipelines at the TESS Science Office and at the TESS Science Processing
Operations Center. Resources supporting this work were provided by the NASA
High-End Computing (HEC) Program through the NASA Advanced Supercomputing
(NAS) Division at Ames Research Center for the production of the SPOC data
products. This paper includes data collected by the TESS mission that are
publicly available from the Mikulski Archive for Space Telescopes (MAST). The
PLATOSpec team would like to thank observers Marek Skarka, Jiří Srba, Zuzana
Balkóvá, Petr Škoda and Luděk Řezba. This publication was produced within the
framework of institutional support for the development of the research
organization of Masaryk University. This research has made use of the Simbad
and Vizier databases, operated at the centre de Données Astronomiques de
Strasbourg (CDS), and of NASA’s Astrophysics Data System Bibliographic
Services (ADS). This work has made use of data from the European Space Agency
(ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia
Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular, the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Ahrer et al. (2023) Ahrer, E.-M., Stevenson, K. B., Mansfield, M., et al. 2023, Nature, 614, 653
  * Aller et al. (2020) Aller, A., Lillo-Box, J., Jones, D., Miranda, L. F., & Barceló Forteza, S. 2020, A&A, 635, A128
  * Almenara et al. (2015) Almenara, J. M., Damiani, C., Bouchy, F., et al. 2015, A&A, 575, A71
  * Anderson et al. (2011) Anderson, D. R., Smith, A. M. S., Lanotte, A. A., et al. 2011, MNRAS, 416, 2108
  * Arcangeli et al. (2018) Arcangeli, J., Désert, J.-M., Line, M. R., et al. 2018, ApJ, 855, L30
  * Arcangeli et al. (2019) Arcangeli, J., Désert, J.-M., Parmentier, V., et al. 2019, A&A, 625, A136
  * Arras & Bildsten (2006) Arras, P. & Bildsten, L. 2006, ApJ, 650, 394
  * Arras & Socrates (2010) Arras, P. & Socrates, A. 2010, ApJ, 714, 1
  * Bailey & Goodman (2019) Bailey, A. & Goodman, J. 2019, MNRAS, 482, 1872
  * Baraffe et al. (2008) Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315
  * Baraffe et al. (2003) Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701
  * Barragán et al. (2022) Barragán, O., Aigrain, S., Rajpaul, V. M., & Zicher, N. 2022, MNRAS, 509, 866
  * Batalha et al. (2019) Batalha, N. E., Marley, M. S., Lewis, N. K., & Fortney, J. J. 2019, ApJ, 878, 70
  * Batygin & Stevenson (2010) Batygin, K. & Stevenson, D. J. 2010, ApJ, 714, L238
  * Beatty et al. (2018) Beatty, T. G., Morley, C. V., Curtis, J. L., et al. 2018, AJ, 156, 168
  * Bensby et al. (2003) Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527
  * Blanco-Cuaresma (2019) Blanco-Cuaresma, S. 2019, MNRAS, 486, 2075
  * Blanco-Cuaresma et al. (2014) Blanco-Cuaresma, S., Soubiran, C., Heiter, U., & Jofré, P. 2014, A&A, 569, A111
  * Bodenheimer et al. (2003) Bodenheimer, P., Laughlin, G., & Lin, D. N. C. 2003, ApJ, 592, 555
  * Bodenheimer et al. (2001) Bodenheimer, P., Lin, D. N. C., & Mardling, R. A. 2001, ApJ, 548, 466
  * Bowler et al. (2023) Bowler, B. P., Tran, Q. H., Zhang, Z., et al. 2023, AJ, 165, 164
  * Brahm et al. (2017) Brahm, R., Jordán, A., & Espinoza, N. 2017, PASP, 129, 034002
  * Burrows et al. (2011) Burrows, A., Heng, K., & Nampaisarn, T. 2011, ApJ, 736, 47
  * Burrows et al. (2001) Burrows, A., Hubbard, W. B., Lunine, J. I., & Liebert, J. 2001, Reviews of Modern Physics, 73, 719
  * Burrows et al. (2007) Burrows, A., Hubeny, I., Budaj, J., & Hubbard, W. B. 2007, ApJ, 661, 502
  * Byrd et al. (1995) Byrd, R., Lu, P., Nocedal, J., & Zhu, C. 1995, SIAM Journal of Scientific Computing, 16, 1190
  * Carmichael et al. (2021) Carmichael, T. W., Quinn, S. N., Zhou, G., et al. 2021, AJ, 161, 97
  * Chabrier & Baraffe (2000) Chabrier, G. & Baraffe, I. 2000, ARA&A, 38, 337
  * Chabrier & Baraffe (2007) Chabrier, G. & Baraffe, I. 2007, ApJ, 661, L81
  * Changeat et al. (2022) Changeat, Q., Edwards, B., Al-Refaie, A. F., et al. 2022, ApJS, 260, 3
  * Charbonneau et al. (2005) Charbonneau, D., Allen, L. E., Megeath, S. T., et al. 2005, ApJ, 626, 523
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Claret (2017) Claret, A. 2017, A&A, 600, A30
  * Collier Cameron & Jardine (2018) Collier Cameron, A. & Jardine, M. 2018, MNRAS, 476, 2542
  * Collins et al. (2017) Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman, F. V. 2017, AJ, 153, 77
  * Coulombe et al. (2023) Coulombe, L.-P., Benneke, B., Challener, R., et al. 2023, Nature, 620, 292
  * Crossfield et al. (2012) Crossfield, I. J. M., Barman, T., Hansen, B. M. S., Tanaka, I., & Kodama, T. 2012, ApJ, 760, 140
  * Dang et al. (2022) Dang, L., Bell, T. J., Cowan, N. B., et al. 2022, AJ, 163, 32
  * Davenport et al. (2014) Davenport, J. R. A., Hawley, S. L., Hebb, L., et al. 2014, ApJ, 797, 122
  * David et al. (2019) David, T. J., Hillenbrand, L. A., Gillen, E., et al. 2019, ApJ, 872, 161
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8
  * El-Badry et al. (2018) El-Badry, K., Rix, H.-W., Ting, Y.-S., et al. 2018, MNRAS, 473, 5043
  * Enoch et al. (2012) Enoch, B., Collier Cameron, A., & Horne, K. 2012, A&A, 540, A99
  * Espinoza et al. (2019) Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, celerite: Scalable 1D Gaussian Processes in C++, Python, and Julia, Astrophysics Source Code Library, record ascl:1709.008
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Fortney et al. (2007) Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661
  * Fuhrmann (1998) Fuhrmann, K. 1998, A&A, 338, 161
  * Fulton et al. (2015) Fulton, B. J., Collins, K. A., Gaudi, B. S., et al. 2015, ApJ, 810, 30
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Ginzburg & Sari (2016) Ginzburg, S. & Sari, R. 2016, ApJ, 819, 116
  * Goldreich & Soter (1966) Goldreich, P. & Soter, S. 1966, Icarus, 5, 375
  * Goodman & Weare (2010) Goodman, J. & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Green (2018) Green, G. M. 2018, The Journal of Open Source Software, 3, 695
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * Guerrero et al. (2021) Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39
  * Guillot & Showman (2002) Guillot, T. & Showman, A. P. 2002, A&A, 385, 156
  * Günther & Daylan (2021) Günther, M. N. & Daylan, T. 2021, ApJS, 254, 13
  * Hatzes & Rauer (2015) Hatzes, A. P. & Rauer, H. 2015, ApJ, 810, L25
  * Heller et al. (2010) Heller, R., Jackson, B., Barnes, R., Greenberg, R., & Homeier, D. 2010, A&A, 514, A22
  * Hirsch et al. (2017) Hirsch, L. A., Ciardi, D. R., Howard, A. W., et al. 2017, AJ, 153, 117
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 357, 367
  * Horch et al. (2014) Horch, E. P., Howell, S. B., Everett, M. E., & Ciardi, D. R. 2014, ApJ, 795, 60
  * Huang et al. (2020a) Huang, C. X., Vanderburg, A., Pál, A., et al. 2020a, Research Notes of the American Astronomical Society, 4, 204
  * Huang et al. (2020b) Huang, C. X., Vanderburg, A., Pál, A., et al. 2020b, Research Notes of the American Astronomical Society, 4, 206
  * Huang & Cumming (2012) Huang, X. & Cumming, A. 2012, ApJ, 757, 47
  * Husser et al. (2013) Husser, T.-O., von Berg, S. W., Dreizler, S., et al. 2013, Astronomy & Astrophysics, 553, A6
  * Hut (1980) Hut, P. 1980, A&A, 92, 167
  * Jackson et al. (2008) Jackson, B., Greenberg, R., & Barnes, R. 2008, ApJ, 678, 1396
  * Jenkins et al. (2016) Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9913, Software and Cyberinfrastructure for Astronomy IV, ed. G. Chiozzi & J. C. Guzman, 99133E
  * Jensen (2013) Jensen, E. 2013, Tapir: A web interface for transit/eclipse observability
  * Jermyn et al. (2017) Jermyn, A. S., Tout, C. A., & Ogilvie, G. I. 2017, MNRAS, 469, 1768
  * Johns-Krull et al. (2008) Johns-Krull, C. M., McCullough, P. R., Burke, C. J., et al. 2008, ApJ, 677, 657
  * Kreidberg et al. (2018) Kreidberg, L., Line, M. R., Parmentier, V., et al. 2018, AJ, 156, 17
  * Kurokawa & Inutsuka (2015) Kurokawa, H. & Inutsuka, S.-i. 2015, ApJ, 815, 78
  * Lainey et al. (2009) Lainey, V., Arlot, J.-E., Karatekin, Ö., & van Hoolst, T. 2009, Nature, 459, 957
  * Leggett (1992) Leggett, S. K. 1992, ApJS, 82, 351
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library
  * Line et al. (2014) Line, M. R., Knutson, H., Wolf, A. S., & Yung, Y. L. 2014, ApJ, 783, 70
  * Ma & Ge (2014) Ma, B. & Ge, J. 2014, MNRAS, 439, 2781
  * Madhusudhan et al. (2011) Madhusudhan, N., Harrington, J., Stevenson, K. B., et al. 2011, Nature, 469, 64
  * Mandushev et al. (2005) Mandushev, G., Torres, G., Latham, D. W., et al. 2005, ApJ, 621, 1061
  * Mansfield et al. (2018) Mansfield, M., Bean, J. L., Line, M. R., et al. 2018, AJ, 156, 10
  * Marleau & Cumming (2014) Marleau, G. D. & Cumming, A. 2014, MNRAS, 437, 1378
  * Marley et al. (2021) Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85
  * Masuda & Winn (2020) Masuda, K. & Winn, J. N. 2020, AJ, 159, 81
  * Matson et al. (2018) Matson, R. A., Howell, S. B., Horch, E. P., & Everett, M. E. 2018, AJ, 156, 31
  * Maxted (2016) Maxted, P. F. L. 2016, A&A, 591, A111
  * Mollière et al. (2022) Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74
  * Mollière & Mordasini (2012) Mollière, P. & Mordasini, C. 2012, A&A, 547, A105
  * Mollière et al. (2015) Mollière, P., van Boekel, R., Dullemond, C., Henning, T., & Mordasini, C. 2015, ApJ, 813, 47
  * Mordasini et al. (2012) Mordasini, C., Alibert, Y., Georgy, C., et al. 2012, A&A, 547, A112
  * Morton (2015) Morton, T. D. 2015, isochrones: Stellar model grid package, Astrophysics Source Code Library, record ascl:1503.010
  * Moses et al. (2013) Moses, J. I., Madhusudhan, N., Visscher, C., & Freedman, R. S. 2013, ApJ, 763, 25
  * Öberg et al. (2011) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16
  * Parviainen (2015) Parviainen, H. 2015, MNRAS, 450, 3233
  * Penev et al. (2012) Penev, K., Jackson, B., Spada, F., & Thom, N. 2012, ApJ, 751, 96
  * Perna et al. (2010) Perna, R., Menou, K., & Rauscher, E. 2010, ApJ, 724, 313
  * Quinn et al. (2014) Quinn, S. N., White, R. J., Latham, D. W., et al. 2014, ApJ, 787, 27
  * Sarkis et al. (2021) Sarkis, P., Mordasini, C., Henning, T., Marleau, G. D., & Mollière, P. 2021, A&A, 645, A79
  * Saumon et al. (1996) Saumon, D., Hubbard, W. B., Burrows, A., et al. 1996, ApJ, 460, 993
  * Sestovic et al. (2018) Sestovic, M., Demory, B.-O., & Queloz, D. 2018, A&A, 616, A76
  * Showman & Guillot (2002) Showman, A. P. & Guillot, T. 2002, A&A, 385, 166
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Smith et al. (2012) Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132
  * Spiegel et al. (2011) Spiegel, D. S., Burrows, A., & Milsom, J. A. 2011, ApJ, 727, 57
  * Stassun et al. (2018) Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018, AJ, 156, 102
  * Stumpe et al. (2014) Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100
  * Stumpe et al. (2012) Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985
  * Tokovinin (2018) Tokovinin, A. 2018, PASP, 130, 035002
  * Tokovinin & Cantarutti (2008) Tokovinin, A. & Cantarutti, R. 2008, PASP, 120, 170
  * Tremblin et al. (2017) Tremblin, P., Chabrier, G., Mayne, N. J., et al. 2017, ApJ, 841, 30
  * Trifonov (2019) Trifonov, T. 2019, The Exo-Striker: Transit and radial velocity interactive fitting tool for orbital analysis and N-body simulations, Astrophysics Source Code Library, record ascl:1906.004
  * Twicken et al. (2018) Twicken, J. D., Catanzarite, J. H., Clarke, B. D., et al. 2018, PASP, 130, 064502
  * Vanzi et al. (2012) Vanzi, L., Chacon, J., Helminiak, K. G., et al. 2012, Monthly Notices of the Royal Astronomical Society, 424, 2770
  * Vissapragada et al. (2022) Vissapragada, S., Chontos, A., Greklek-McKeon, M., et al. 2022, ApJ, 941, L31
  * Šubjak et al. (2023) Šubjak, J., Lodieu, N., Kabáth, P., et al. 2023, A&A, 671, A10
  * Šubjak et al. (2020) Šubjak, J., Sharma, R., Carmichael, T. W., et al. 2020, AJ, 159, 151
  * Wu & Lithwick (2013) Wu, Y. & Lithwick, Y. 2013, ApJ, 763, 13
  * Yan et al. (2023) Yan, F., Nortmann, L., Reiners, A., et al. 2023, A&A, 672, A107
  * Yang & Wei (2022) Yang, F. & Wei, X. 2022, PASP, 134, 024401
  * Yee et al. (2023) Yee, S. W., Winn, J. N., Hartman, J. D., et al. 2023, ApJS, 265, 1
  * Yee et al. (2020) Yee, S. W., Winn, J. N., Knutson, H. A., et al. 2020, ApJ, 888, L5
  * Youdin & Mitchell (2010) Youdin, A. N. & Mitchell, J. L. 2010, ApJ, 721, 1113
  * Zechmeister & Kürster (2009) Zechmeister, M. & Kürster, M. 2009, A&A, 496, 577
  * Zhou et al. (2019) Zhou, G., Bakos, G. Á., Bayliss, D., et al. 2019, AJ, 157, 31
  * Zhu et al. (1997) Zhu, C., Byrd, R. H., Lu, P., & Nocedal, J. 1997, ACM Trans. Math. Softw., 23, 550–560
  * Ziegler et al. (2021) Ziegler, C., Tokovinin, A., Latiolais, M., et al. 2021, AJ, 162, 192
  * Zucker & Mazeh (1994) Zucker, S. & Mazeh, T. 1994, ApJ, 420, 806

##  Appendix A Additional material

Figure A.20: The correlations between the free parameters of the 2-star model
from the Nested Sampling. At the end of each row is shown the derived
posterior probability distribution. Figure A.21: The correlations between the
free parameters of the 3-star model from the Nested Sampling. At the end of
each row is shown the derived posterior probability distribution.