differences. We also investigated if there might a relation between the the
RVs and the SNR in the spectra (e.g. due to sky-brightness from the Moon) but
there is no correlation apparent (for the RVs taken with the most frequent
exposure time of 1200 sec, a correlation coefficient of -0.05 was found).
Hence, an _identifiable_ effect of the Moon-reflected solar spectrum onto the
measured RVs is hence very minor.

However, we consider that the 29.5 d peak in the _HARPS-N_ spectrograms
remains of questionable nature, and now turn our attention to the neighbouring
peak at $\approx$ 27.4 days. Fig. 24 shows BGLS spectrograms of all
contributing instruments, and it is of note that data from the APF -- which
contributed with the the second largest set of RVs -- have their strongest
peak at 26.8 d. Also, the HIRES RVs show a peak in the same period-range,
which is very broad due to the small sample of only 12 RVs. However, peaks in
that range are absent in periodograms from CARMENES and iSHELL data. We note
that these are also the instruments whose spectral coverage is the most red-
wards (see Table 3), and a Moon-reflected reflected solar spectrum would
generate a weaker signal in them. Lastly, in a combination of all available
RVs, the peak at 27.4 d is the highest overall, and is significantly stronger
than the one at 29.5 d.

In Fig. 25 we provide plots of the development of the SNR (signal to noise
ratio) and the best-fitting amplitude $K$ of RV signals at the 27.4 d and
29.53 d periods, versus the number of RV points (counting from the first
measurement), following the precepts of Mortier & Collier Cameron (2017).
These plots are shown for the two largest sets of RVs, those from _HARPS-N_
and from the APF. In the plots for _HARPS-N_ , the SNR degrades at either
period near the 40th point, which is likely due to a lesser consistency of
these signals across RV coverages spanning more than one observing season. On
the other hand, in the plots for the APF (which cover only one observing
season), the 27.4 d signal shows a steady increase in SNR and a rather
constant amplitude $K$, whereas the 29.5 d signal shows a less consistent
picture, more similar to the one from _HARPS-N_. As is stated in Sect. 5.2,
fits of Keplerian orbits to the 27.4 d signal where however significantly
worse than those to the 29.5 d one.

In summary, we cannot decide on a clear preference that either of these
signals represent a true signal from TOI-1416, nor about their actual nature,
and conclude that the RV signals with a 27.4 or 29.5 day period are at most
tentative of a further planet at either of these periods.

Figure 24: BGLS periodograms of the RVs of all contributing instruments. In
the left column are those that show a peak near 27 or 29 days (vertical red
lines); in the center are those that don’t, and the right panel shows a
periodogram with the RVs from all instruments. The blue numbers indicate the
number of RVs used. Figure 25: Development of the SNR (signal to noise ratio)
and the best-fitting amplitude ’bestK’ of RV signals of the 27.4d (top panels)
and the 29.5 d period (bottom) , versus the number of RV points since the
first measurement. The left panels are based on RVs data from _HARPS-N_ and
the right ones on data from the APF.

##  Appendix C Modelling of the lightcurve with UTM/UFIT

The normalised and gradient-corrected lightcurve around transits of planet
$b$, whose preparation is described in Sect. 2, was used for transit fits
using the Universal Transit Modeller / Fitter (UTM/UFIT, Deeg,
2014)222222Available at https://github.com/hdeeg/utm_ufit/. In brief, UTM is a
lightcurve modeller for all kinds of eclipsing or transiting configurations
between any number and kind of objects, such as stars, planets, moons and
rings. UFIT was developed as a wrapper to UTM to perform fits, albeit it has
been extended to accept several further modelling modules (such as the one
used for the FCO fit described in Appendix D). As the core modelling engine,
UTM may use either pixelised object representations suitable for arbitrary
configurations of multiple occulters, or an analytical ’fast mode’ suitable
for basic transit configurations, which employs the exofast_occultquad.pro
routine from the EXOFAST library (Eastman et al., 2013, 2019); the latter one
was used in this work. UFIT (Universal Fitter) permits the fitting of any of
UTM’s input parameters, either with the Amoeba algorithm or through the
generation of MCMC chains using the Differential Evolution Markov Chain method
of Ter Braak (2006). Its implementation is based on the EXOFAST_DEMC routine
from the same library, but with an extension that permits the constraining of
free parameters by several types of symmetric and asymmetric priors.

UTM permits the modelling of transit curves from any set of input parameters
that is fully able to describe an orbiting system. For this work, we modelled
the light-curves curves against the following set of parameters (these were
also free parameters in the fits): Orbital period $P_{\mathrm{orb}}$, transit
epoch $T_{0}$, scaled planetary radius $R_{\mathrm{p}}/R_{\star}$, the stellar
density232323The stellar density was mainly chosen for compatibility with the
set of input parameters used by pyaneti, described in Sect. 5.2. A pre-
processor routine to UTM converts the stellar density into the usually used
ratio of the semi-major axis versus the stellar radius,
$a_{\mathrm{p}}/R_{\star}$, using e.g. Eq. (31) of Barragán et al. (2019).
$\rho_{\star}$ and the transit impact parameter $b$. Initial values for these
fits were taken from the SPOC’s data validation summary for TOI-1416. For the
stellar limb-darkening (LD), a quadratic LD law was used, albeit for the
fitting we used the $q_{1}$ and $q_{2}$ coefficients for an optimised sampling
introduced by Kipping (2013). An absolute offset in flux values was a further
free parameter in our fits, in order to account for potential errors in the
normalisation of the flux described in Sect. 2. The orbital eccentricity was
kept to zero.

Initial fits were performed with the Amoeba algorithm, leading to an
intermediate transit-model that was used for the initial parameters for an
MCMC sequence. First efforts without constraints on the input parameters
showed significant correlation between the impact parameter, the stellar
density, and the planet radius. Also, the limb-darkening parameters could only
be poorly constrained from the fits. We therefore choose to impose Gaussian
priors on the stellar density, taken from the adopted value in Table 5 and on
the limb-darkening. For the later, we used the tabulation of LD coefficients
for the TESS satellite bandpass by Claret (2017, Table 25 for the ATLAS model
with plane-parallel geometry) and an interpolation to the adopted stellar
parameters from Tables 4 and 5. The obtained values for a square LD law as
defined in Eq. (2) of Claret (2017) were $u_{1}=0.4545\pm 0.05$ and
$u_{2}=0.1880\pm 0.05$, which were converted into priors for the $q_{1}$ and
$q_{2}$ coefficients, given in Table 7.

The final MCMC sequence consisted of 16 parallel chains that were iterated
until a sufficient mixing of parameters was achieved, based on the Gelman-
Rubin statistics following the precepts of Ford (2006) and Eastman et al.
(2013). The resultant values (included in Table 7) were then derived from the
posterior distributions of the parameters, based on 4280 steps, after a burn-
in period of $\approx$ 1000 steps. These distributors were in all cases close
to Gaussian shapes; see Fig. 26 for this and several further diagnostic plots
from the MCMC sequence. The best-fit transit model against the phase-folded
input lightcurve is also included in Fig. 1.

Figure 26: Graphical output from MCMC sequence performed by UFIT that led to
the results reported in Table 7. In all panels, the parameters are indicated
by the keywords used in UFIT: 1period: planet period, 1trepoch: epoch of
transit, 1radi: relative planet radius, 0densCGS: stellar density in CGS
units, 1impact: impact parameter, 0limbd and 0limbd2: Limdarkening
coefficients $q_{1}$ and $q2$, ooff: off-transit flux-offset against zero. The
sub-figures are: a) Values of the parameters against link-number of the MCMC
sequence, excluding burn-in. Each MCMC chain is shown by a different color.
The lowest panel shows the evaluation of the $\chi^{2}$ values. b) Scatter
plot of parameters versus the $\chi^{2}$ value. c) Histograms of parameter
distributions. The median value is shown by the vertical black line; the
dashed lines delimit the 68.3% credible interval and the red line gives the
value of the best fit. d) Cornerplot of the correlations among parameters. The
red crosses give the values of the best fit.

##  Appendix D Detection of the transiting planet in RV data by FCO analysis

The FCO (Floating Chunk Offset) method, developed by Hatzes et al. (2010) and
Hatzes (2014), is best suited for the determination of RV amplitudes of short-
periodic planets whose nightly RV variations are expected to be larger than
the individual RV measures’ uncertainties. In short, sets of nightly RV data
-- with at least two well separated data-points per night -- are treated as
independent chunks of data with unknown (free) RV offsets. Systematics (both
instrumental and effects from other planets or stellar activity) on time-
scales larger than a single night are therefore suppressed by the FCO method.
RV offsets for each nightly set of RVs are then fitted against an RV model and
the parameters of the best fit are obtained.

Only the RVs from _HARPS-N_ were used in this analysis. The FCO method could
not be applied to the data from the other telescopes, because all their RVs
are single data-points in a given night (with the exception of two nights from
APF, where RVs spaced about 20 min apart were obtained, which is too short a
separation to be suitable for the FCO analysis). In the _HARPS-N_ data, there
are 28 nights in which two or more RVs were obtained, which enabled the use of
77 out of the 93 RVs from _HARPS-N_. For the fitting of the RVs, we used the
same UFITfitter as described in Appendix C, but with a modelling module
ufit_rvcurve for the generation of a Keplerian RV model from any number of RV
data-sets, each with its own RV-offset $\gamma_{i}$. For the FCO method, each
’chunk’ with a nightly set of two to six RVs is considered an independent set
of data.

A fit using the FCO method for a candidate with an ephemeris known from
transits and assuming a circular orbit is in principle very simple, as it
contains as free parameters only the RV amplitude $K$ and the nightly RV
offsets $\gamma_{i}$, with $i=1,...,28$ indexing the individual nights.
However, when using UFIT with both the AMOEBA or the MCMC fitter, resultant RV
amplitudes tended to be stuck close to the amplitude’s initial value. This
behaviour was caused by the large number of 28 nights, where each one
corresponds to a free parameter $\gamma_{i}$. Due to this, either fitter found
it difficult to vary the RV amplitude, since any improvement in the fit
requires that most of the nightly RV offsets are changed _simultaneously_ by
the correct amounts. This is difficult to achieve for any fitting algorithm,
and is an expected behaviour when free parameters are strongly correlated. We
therefore kept the RV amplitude -- and hence the entire RV model -- fixed and
fitted only for the RV offsets $\gamma_{i}$, which reduces the fitting task to
a set of 28 simple linear fits. The input RV amplitude was then stepped
through a series of suitable values and the $\chi^{2}$ of each corresponding
fit was logged, which led to the curve shown in Fig 27. The minimum of this
curve, and the range where $\chi^{2}$ increases by 1, indicate an amplitude of
$K_{b}=2.14\pm 0.35$ m s-1. For the best-fit model with $K_{b}=2.14$ m s-1,
the reduced $\chi^{2}$ is 0.87 and the residuals of the RVs against the model
have an rms of 0.82 m s-1. Fig. 28 shows the RVs of the nightly chunks against
the best RV model, with a zoomed-out section across three nights of the RV
time-series showing the excellent fit in that range. We note that a further
FCO fit of the _HARPS-N_ RVs with the open source code pyaneti242424Available
at https://github.com/oscaribv/pyaneti (Barragán et al., 2019, 2022, see also
Sect. 5.2) gave a nearly identical result, of $2.12\pm 0.36$ ${\rm
m\,s^{-1}}$.

Figure 27: Best-fit $\chi^{2}$ from FCO fits of _HARPS-N_ data against RV
models of TOI-1416b with fixed RV-amplitudes $K$.

Figure 28: Upper panel: The phase-folded RV model (green line) of TOI-1416 b
which corresponds to the best fit from the FCO method, obtained by vertically
offsetting nightly chunks of RV data against the model. RV points from the
same nights have identical colours. Lower panel: A small section of the RV
model plotted against time, with RVs from three different nights.

##  Appendix E The Neptune desert in radius or mass versus insolation or
effective temperature

Figure 29: Similar to Fig. 16, but with planet radii (top panel) and masses
(bottom panel) plotted against the planets’ insolation (lower X-axis) and
their effective temperature (upper X-axis). The solid black lines show the
delineation of the Neptune Desert against insolation resp. $T_{\mathrm{eff}}$,
whereas the horizontal dotted black lines show the same lower limits to the
Neptune Desert as those proposed for periods of $P\lesssim 2d$. The dashed
orange line in the upper panel indicates the radius valley against insolation
from Petigura et al. (2022). TOI-1416 $b$ is indicated by the filled red
circle and $c$ by the unfilled one.

In Fig. 29 we show plots similar to Fig. 16, but plotting the planets’ radii
and masses against the incident bolometric flux or insolation, instead of
orbital period (Fig. 16). The insolation was calculated from first principles
from values obtained from the NASA Exoplanet Archive. It is of note that the
upper boundary of the Neptune desert is significantly sharper than in the
plots against orbital period, whereas the lower boundary remains diffuse; in
particular for the plot of planet masses. Given the sharper upper boundary, we
propose an upper radius-limit of the Neptune desert against insolation as:

$\log R_{hi}/R_{\oplus}=0.248\log S+0.33\ ,\hskip 14.22636ptS\gtrsim 150,$ (4)

where $S$ is the insolation relative to the Earth’s insolation. The
corresponding lower limit is then

$\log R_{lo}/R_{\oplus}=\left\\{\begin{array}[]{rl}-0.51\log S+1.74\
,&150\lesssim S\lesssim 1000\\\ 0.20\ ,&\text{$S\gtrsim
1000$},\end{array}\right.$ (5)

where the same limit of $R_{lo}=1.60R_{\oplus}$ as given in Sect. 6 for very
short orbital periods applies also to the strongest insolations. For the limit
of the desert against mass, there are much fewer planets with mass
measurements, and we only derive an upper limit of

$\log M_{hi}/M_{Jup}=0.74\log S-2.35\ ,\hskip 14.22636ptS\gtrsim 150,$ (6)

whereas a lower limit cannot be discerned with reliability, given the small
sample of known short-period low-mass planets, that is furthermore suffering a
strong selection effect against detectability towards smaller masses. In
Fig.29 we hence indicate only the same lower mass limit that is given in Sect.
6 for very short orbital periods, namely $M_{lo}=8.9M_{\oplus}$, resp. $\log
M_{lo}/M_{\mathrm{jup}}=-1.55$. In above limits for both radius and mass, we
maintain the gradients of the limits against period by Mazeh et al. (2016),
after multiplication of $\log P$ with a factor of -4/3, which arises from the
dependency of $S$ on the period.

For convenience, we also provide the same limits against the planets’
effective temperature $T_{\mathrm{eff}}$, using the conversion
$T_{\mathrm{eff}}=(S)^{1/4}\ 255$ K, with the 255 K corresponding to the
effective temperature of the Earth. We then obtain for the limits of radius
against $T_{\mathrm{eff}}$:

$\log R_{hi}/R_{\oplus}=0.99\log T_{\mathrm{eff}}-2.72\ ,\hskip
19.91684ptT_{\mathrm{eff}}\gtrsim 900K\\\ $ (7) $\log
R_{lo}/R_{\oplus}=\left\\{\begin{array}[]{rl}-2.04\log T_{\mathrm{eff}}-3.17\
,&900K\lesssim T_{\mathrm{eff}}\lesssim 1450K\\\ 0.20\
,&\text{$T_{\mathrm{eff}}\gtrsim 1450K$}\end{array}\right.$ (8)

and for the upper mass limit of the desert:

$\log M_{hi}/M_{Jup}=3.0\log T_{\mathrm{eff}}-9.5\ ,\hskip
19.91684ptT_{\mathrm{eff}}\gtrsim 900K,$ (9)

with the same lower limit as indicated previously against insolation. With
$logM/M_{\oplus}=2.50\log M/M_{Jup}$, we may easily convert the desert mass
limits to the units of Earth masses used in Fig 29.

##  Appendix F Further figures mentioned in main text

  * •

Fig. 30 gives a comparison of the RVs measured with _HARPS-N_ using the DRS
and the serval pipelines.

  * •

Fig. 31 shows the output of the pyaneti joint-fit for Model 3, with two
planets $b$ and $c$.

Figure 30: Comparison between _HARPS-N_ RVs measured with the _HARPS-N_ DRS
from CCFs (open red circles) and with serval (open green circles). The red
crosses show the difference between the two data-sets. The five points in
which this difference is significantly negative correspond to exposures that
were prematurely terminated, and which are not correctly processed by the DRS.
Both data-sets have been averaged to zero without considering these five
points. The difference between the two data-sets has a standard deviation of
0.92 m s-1 (excluding again these five points).

Figure 31: Like Fig. 13, but for Model 3 with the additional fit for a
Keplerian signal with the lunar synodic period (29.53d). The lower left panel
is again for the transiting planet $b$, while the lower right panel shows the
RVs folded over the period of the additional signal. A similar plot for Model
2 does not show any relevant differences to the shown one.