Jupiters and improved masses in 38 Kepler and K 2 small-planet systems from 3661 high-precision HARPS-N radial velocities . No excess of cold Jupiters in small-planet systems

The exoplanet population with relatively short orbital periods (P < 100 d) around solar-type stars is dominated by super-Earths and sub-Neptunes. These planets are, however, missing in our Solar System, and the reason for that is unknown. Two theoretical scenarios invoke the role of Jupiter as the possible culprit: Jupiter may have acted as a dynamical barrier to the inward migration of sub-Neptunes from beyond the water iceline or, alternatively, may have reduced considerably the inward flux of material (pebbles) required to form super-Earths inside that iceline. Both scenarios predict an anti-correlation between the presence of small planets and that of cold Jupiters in exoplanetary systems. To test that prediction, we homogeneously analyzed the radial-velocity measurements of 38 Kepler and K2 transiting small-planet systems gathered over almost 10 years with the HARPS-N spectrograph, as well as publicly available radial velocities collected with other facilities. We used Bayesian differential evolution Markov chain Monte Carlo techniques, in some cases coupled with Gaussian Process regression to model nonstationary variations due to stellar magnetic activity phenomena. We detected five cold Jupiters in three systems, two in Kepler-68, two in Kepler454, and a very eccentric one in K2-312. We also found linear trends caused by bound companions in Kepler-93, Kepler-454 and K2-12, with slopes still compatible with a planetary mass for outer bodies in the Kepler-454 and K2-12 systems. By using binomial statistics and accounting for the survey completeness, we derived an occurrence rate of 9.3+7.7 −2.9 % for cold Jupiters with 0.3 − 13 MJup and 1-10 au, which is lower but still compatible at 1.3σ with that measured from radial-velocity surveys for solar-type stars, regardless of small-planet presence. This does not allow us to draw a firm conclusion about the predicted anti-correlation between small planets and cold Jupiters, which would require a considerably larger sample. Nevertheless, we found no evidence of previous claims of an excess of cold Jupiters in small-planet systems. As an important by-product of our analyses, we homogeneously determined the masses of 64 Kepler and K2 small planets, reaching a precision better than 5, 7.5 and 10σ for 25, 13 and 8 planets, respectively. Finally, we release to the scientific community the 3661 HARPS-N radial velocities used in this work, which mostly benefit from an improved Data Reduction Software that corrects for previous subtle systematic effects.


Introduction
One of the most striking findings from the detection of almost 4000 transiting planets so far is that the most common type of exoplanets in relatively close orbits (orbital periods P < 100 d) around solar-type stars are small planets (SPs) with radii 1 < R p < 4 R ⊕ . These are hosted by about half of solar-type stars in the Milky Way (e.g., Winn & Fabrycky 2015 and references therein) and can be subdivided in two main classes: i) highdensity super-Earths with 1.0 R p 1.7 R ⊕ having a rocky composition, and ii) lower-density sub-Neptunes with 1.7 The second framework is described in Lambrechts et al. (2019) and is based on the formation of super-Earths inside the water iceline through pebble accretion (e.g., Ogihara et al. 2015). With extensive simulations, Lambrechts et al. (2019) showed that the outcome in the formation of super-Earths is strongly dependent on the amount of pebble flux drifting inwards from the outer regions of the protoplanetary disk. Low pebble fluxes would generate Mars-mass embryos, which may then grow to terrestrial planets at orbital distances a 0.4 AU through mutual collisions after the disk dissipation, as likely occurred in the Solar System (e.g., Raymond et al. 2014). Conversely, higher fluxes of pebbles would produce more massive embryos in shorter time. These could migrate towards their parent star if gas is still present in the disk, and thus form a compact system of close-in super-Earths (see their Figure 3).
If a giant planet forms, it opens a gap in the disk, considerably reducing or even halting the inward flux of pebbles from the regions outside its orbit. The formation of Jupiter might thus explain why the Solar System contains no short-period super-Earths, but terrestrial planets only: Jupiter might have reduced the flux of material required to form bigger planets within the water iceline.

Theoretical predictions and observational corroboration
Both scenarios above predict an anti-correlation between the presence of short-period small planets and that of cold Jupiters 1 (CJs) in exoplanetary systems, which could be tested observationally. Previous works by Zhu & Wu (2018) and Bryan et al. (2019) seem to contradict this anti-correlation, reporting instead an excess of CJs in small-planet systems. Specifically, by analyzing 65 transiting and non-transiting systems with low-mass (1 < M p < 10 M ⊕ ) and short-period (P < 100 d) planets, Bryan et al. (2019) found an occurrence rate f CJ|SP = 36 +7 −6 % for gaseous giant planets with masses M p = 0.5 − 13 M Jup and semi-major axis a ∼ 1 − 10 AU. This rate is higher than that found by Wittenmyer et al. (2020), i.e. f CJ = 20.2 +6.3 −3.4 % 2 , for M p = 0.3 − 13 M Jup and a = 1 − 10 AU, from the AAT, HARPS/ESO, HIRES/Keck and CORALIE radial-velocity (RV) data of solar-type stars with time spans longer than 8 years, irrespective of small-planet presence.
Other works, based on RV surveys only, have tried to estimate the frequency of CJs in low-mass planet systems and/or the frequency of low-mass planets in CJ systems, sometimes with apparently conflicting results. For instance, Barbato et al. (2018) reported no low-mass planets with M p = 10 − 30 M ⊕ and P < 150 d in 20 CJ systems around solar-type stars observed with HARPS. On the contrary, from the California Legacy Survey conducted with the HIRES/Keck and APF/Lick spectrographs Rosenthal et al. (2022) found that planets with M p = 2 − 30 M ⊕ and P 150 d may occur approximately twice as frequently around CJ-host stars. However, they observed no significant differences in the occurrence of inner low-mass planets with and without CJ siblings, when limiting their range in mass to M p = 2 − 20 M ⊕ .
As also noted by Rosenthal et al. (2022), RV surveys are sensitive to more massive inner planets than transit surveys, besides determining minimum masses only. Moreover, the adopted ranges in semi-major axes by Rosenthal et al. (2022) for both the inner low-mass planets and the outer CJs differ from those in Bryan et al. (2019), and for instance include warm Jupiters with a = 0.23 − 1 AU. This prevents a straightforward comparison of their results with those in Bryan et al. (2019).

Aim of the present work
In the present work we aim to test the theoretical predictions of Izidoro et al. (2015) and Lambrechts et al. (2019) by searching for cold Jupiters and determining their occurrence rate in 38 transiting systems, 19 observed by Kepler and 19 by K2, 14 of which are in common with the sample studied by Bryan et al. (2019). For this purpose, we used 3661 high-precision HARPS-N radial velocities, 3471 out of which were collected by the HARPS-N Guaranteed Time of Observations (GTO) consortium, and the remaining 190 RVs by other groups, mainly for the purpose of determining the masses and densities of Kepler and K2 transiting planets. Nonetheless, we monitored these systems over the years specifically to look for outer giant planets.
An important by-product of our RV analyses is the improvement in the precision and/or accuracy of planetary masses and densities, thanks to the use of a significant number of yet unpublished HARPS-N/GTO RVs as well as the first combination of HARPS-N/GTO RVs with literature measurements. studies have shown that giant planets might have masses greater than the deuterium burning limit of 13 M Jup (e.g., Chabrier et al. 2014), the formation models in Sect. 1.1 and 1.2 generally refer to cold Jupiters with M p < 13 M Jup .

Target selection and radial-velocity data
Among the Kepler and K2 systems observed by the HARPS-N/GTO program we chose those that i) host small (low-mass) planets with radius 1 < R p < 4 R ⊕ , mass 1 < M p < 20 M ⊕ and orbital period P < 100 d; and ii) have at least 15 HARPS-N RV measurements for a time span longer than ∼ 1 yr. This resulted on the selection of the vast majority of Kepler and K2 systems monitored by the HARPS-N/GTO program. Like Zhu & Wu (2018), we adopted a wider range in M p than Bryan et al. (2019), given that several SPs are known to have M p > 10 M ⊕ .
We also included Kepler-22, even though it meets the second criterion only, because the possible presence of cold Jupiters may provide valuable information on the architecture of a system with a planet in the habitable zone . However, it was not counted in the computation of f CJ|SP in shortperiod small-planet systems, because Kepler-22b has an orbital period of 289 d > 100 d.
The HARPS-N radial velocities used in this work were extracted with the original HARPS-N Data Reduction Software version 3.7 from the stellar spectra obtained before the early failure of the red side of the HARPS-N charge-coupled device (CCD) in late September 2012 (e.g., Bonomo et al. 2014), and with the updated DRS version 2.3.5 from the spectra gathered afterwards. This latter version of the pipeline, adapted from the ESPRESSO spectrograph to HARPS-N (Dumusque 2021), computes a more stable wavelength solution through a careful selection of the lines of the Thorium-Argon calibration lamp, by avoiding saturated Thorium and Argon lines, and also corrects for i) possible RV long-term variations due to changing levels in the flux of the Thorium-Argon calibration lamp with time, and ii) an offset in the DRS v3.7 data, occurring at the beginning of June 2020 for the replacement of the Thorium-Argon calibration lamp.
We observed the majority of the stars in our sample in OBJ_AB observing mode, i.e. with fiber A on the target and fiber B on the sky to monitor possible contamination by moonlight. For the brightest stars of our sample, namely /Wolf 503 and K2-312/HD 80563, we used simultaneous calibration with fiber A on the target and fiber B on the calibration Thorium-Argon or Fabry-Perot lamp, to achieve higher accuracy on the relative RVs. We extracted the RVs by cross-correlating the spectra with a stellar template close to the stellar spectral type (e.g., Pepe et al. 2002). The only exception is the early M-/late K-dwarf K2-3, for which we used the TERRA software (Anglada-Escudé & Butler 2012) to overcome the issue of distorted CCF for cooler stars (Rainer et al. 2020), thereby achieving a reduced RV scatter . For Kepler-10 we performed an additional reduction using the Yarara-v2 tool (Cretignier et al. 2021(Cretignier et al. , 2022, because it proved to slightly enhance the detectability of the planet-induced Doppler signals (Bonomo et al. in prep.).
Possible contamination of the HARPS-N spectra by moonlight was checked following Malavolta et al. (2017) and corrected by computing the cross-correlation function (CCF) after subtracting the flux of fiber B from the flux of fiber A. This procedure led to a reduced RV scatter in a few systems, the most evident cases being  Five systems, namely K2-96/HD 3167, K2-106, K2-111, K2-131 and K2-135/GJ 9827, were also observed with HARPS-N by other groups. To have homogeneous HARPS-N datasets on longer time spans and take advantages of the aforementioned improvements of the HARPS-N pipeline, we recomputed all the HARPS-N RVs with the DRS-v2.3.5 from both the spectra acquired by the GTO and the other publicly available spectra.
For each system, we also collected the published RVs gathered with spectrographs other than HARPS-N, such as HIRES/Keck, HARPS/ESO, ESPRESSO/VLT, PFS/MagellanII and APF, and analyzed them along with the HARPS-N RVs (see Sect. 3.2). This combination is needed to improve the constraints on the presence or lack of CJs, determine more precise/accurate orbital and physical parameters of the detected CJs, and achieve a better precision on the masses and densities of the inner small (low-mass) planets. A few RV datasets with a limited number of RVs and/or considerably lower precision than our HARPS-N RVs were discarded, as they do not yield any improvement in the orbital solution, while they require additional free parameters (the radial-velocity zero point and the uncorrelated jitter term, cf. Sect. 3.2). We did not use the 71 available HIRES/Keck RVs (Weiss et al. 2016) for the analysis of the Kepler-10 system, because they tend to reduce the RV semi-amplitudes of the long-period small planets Kepler-10c and Kepler-10d (Sect. 4.3) compared to those obtained with the 291 HARPS-N RV measurements only (Bonomo et al. in prep.), even though the simultaneous modeling of both datasets is mainly driven by the much more numerous HARPS-N RVs. In fact, also Weiss et al. (2016) showed that the signal of Kepler-10c is practically undetected in the HIRES data.
We searched for possible outliers in the RV datasets for each system using the Chauvenet's criterion 3 and removed them. For systems with long-period RV signals, such as long-term slopes and/or Keplerians of CJs, we applied the Chauvenet's criterion after removing those signals. Through visual inspection we checked that this criterion efficiently removes all the clear outliers. Table 1 lists the 38 systems in our sample, the multiplicity (single planet or multiple planets) of the transiting SPs, the stellar parameters, the total number of RVs and the number of HARPS-N RVs used in this work after the removal of outliers, and the total time span of the observations. For each target, Table 2 reports the epochs, values and formal uncertainties of the HARPS-N RVs, the activity indicators of the CCF, i.e. Full-Width-at-Half-Maximum (FWHM), Contrast and Bisector span of the CCF, and the spectroscopic activity indexes S-index and log R ′ HK . We warn that correlated variations of the FWHM and Contrast may have occurred at certain epochs due to changes in the focus of HARPS-N, but do not affect the RVs because the product [FWHM · Contrast] remains constant.

Updated atmospheric and physical stellar parameters
For uniformity with previous studies of the systems published by the HARPS-N/GTO consortium, we derived the atmospheric parameters, i.e. effective temperature T eff , metallicity [Fe/H] and surface gravity log g, for the host stars Kepler-22, Kepler-109, Kepler-323, Kepler-409, Kepler-1876, K2-12, K2-38, K2-106, K2-131 and K2-167 from the HARPS-N spectra. To this end, we employed two independent methods, ARES+MOOG (e.g., Mortier et al. 2014) and SPC , and adopted as final parameters and uncertainties the average values and error bars provided by the two methods (see Mortier et al. 2018 for more details).
To determine the stellar physical parameters, i.e. mass, radius and age, we used the EXOFASTv2 tool (Eastman 2017;Eastman et al. 2019), which simultaneously adjusts the stellar radius, mass and age through a Bayesian differential evolution Markov chain Monte Carlo (DE-MCMC) method (Ter Braak 2006), by simultaneously fitting the stellar Spectral Energy Distribution (SED) and employing the MESA Isochrones and Stellar Tracks (MIST) (Paxton et al. 2015). To sample the stellar SED, we utilized the WISE W1, W2, W3 and W4 infrared magnitudes (Cutri et al. 2021), the 2-MASS near-infrared J, H and K magnitudes (Cutri et al. 2003), the optical Tycho B T and V T magnitudes (Høg et al. 2000) and/or APASS Johnson B, V and Sloan g ′ , r ′ , i ′ magnitudes (Henden et al. 2016). We imposed Gaussian priors on the T eff and [Fe/H] as derived from the analyses of the HARPS-N spectra, and on the Gaia EDR3 parallax (Gaia Collaboration et al. 2021). A uniform prior was instead used for the V-band extinction, A V , with upper limits provided by reddening maps (Schlegel et al. 1998;Schlafly & Finkbeiner 2011).
For the remaining 27 systems we adopted the previously published stellar parameters, giving preference to those derived using the Gaia parallaxes and/or asteroseismic analyses of the Kepler light curves (e.g., Kepler-10, Kepler-454 and Kepler-107). In one case, i.e. Kepler-20, we redetermined the stellar parameters because the prior on the Gaia parallax yields slightly more precise and accurate stellar radius and mass.
The atmospheric and physical parameters of all the host stars in our sample are listed in Table 1.

Orbital fitting
We modeled the RV data of all the 38 systems in our sample with non-interacting Keplerian orbits and a slope, by maximizing a Gaussian likelihood function (e.g., Ford 2006) through a DE-MCMC technique. The parameters of the Keplerian model for each planet in a given system are the inferior conjunction time T c , which is equivalent to the transit midtime for transiting planets; the orbital period P; the widely adopted parameterization √ e cos(ω) and √ e sin(ω) of eccentricity e and argument of periastron ω; and the RV semi-amplitude K. We included linear slopes to check for significant long-term trends, and fitted for the RV zero point γ i and jitter term σ jit,i for the i-th RV dataset gathered with the i-th spectrograph. The jitter terms σ jit were summed in quadrature to the formal RV uncertainties σ RV (e.g., Gregory 2005) to account for additional white noise of unknown origin (stellar or instrumental).
Correlated noise in the RV time series caused by stellar magnetic activity was modeled through Gaussian Process (GP) regression (e.g., Haywood et al. 2014;Grunblatt et al. 2015) within the same DE-MCMC tools, using a covariance matrix described by either the quasi-periodic (QP) kernel in the original form of Rasmussen & Williams (2006): or the simpler squared-exponential (SE) kernel where h is the semi-amplitude of the correlated noise, λ 1 the correlation decay timescale, P rot the period of the quasi-periodic variations, and λ 2 the inverse complexity harmonic parameter. The hyper-parameters λ 1 , P rot and λ 2 can be associated respectively with the decay timescale of the active regions, the stellar rotation period and the complexity of the activity signals (with λ 2 ∼ 3 − 5 approaching simpler sinusoidal signals). We point out that the parameter λ 2 in Eq. 1 is twice the equivalent parameter w in other implementations of the quasi-periodic kernel (e.g., Grunblatt et al. 2015;Damasso et al. 2018).
The SE kernel was employed for Kepler-93 and K2-2/HIP 116454 to account for RV variations on timescales longer than the stellar rotation, which are likely due to magnetic activity cycles and/or shorter Rieger-type cycles (Rieger et al. 1984). The DE-MCMC analysis with the SE kernel proved to model such variations efficiently, producing flat and Gaussiandistributed residuals. The analysis of the same data with the QP kernel yielded very low acceptance rates, which indicates that the GP-QP model is more complex than needed to model long-term variations in both systems, leaving the QP hyperparameters P rot and λ 2 practically unconstrained.
We imposed several priors on the model parameters as well as on the GP hyper-parameters in case GP regression was used, i.e. in the presence of correlated noise. Specifically, we used Gaussian priors on the transit time T c and period P of the inner transiting SPs, as provided by the transit ephemerides derived in previous analyses of the Kepler and K2 light curves (see the second and third columns, and references in Table 10). With regard to orbital eccentricities, we adopted − circular orbits for the closest SPs, whose orbital circularization time is considerably shorter than the stellar age (e.g., Matsumura et al. 2008Matsumura et al. , 2010. Null eccentricities for these planets are also consistent with the observation of their secondary eclipses at orbital phases ∆φ ∼ 0.5 from transits (Singh et al. 2022). − half-Gaussian priors with zero mean and σ e = 0.098 for the transiting SPs in multiple systems, following the finding of Van Eylen et al. (2019). This prior prevents the fit from converging to spurious eccentricities, a well-known critical effect occurring for low signal-to-noise Doppler signals (e.g., Zakamska et al. 2011;Hara et al. 2019), being the typical RV semi-amplitudes of the transiting SPs in our sample usually comparable to the RV scatter. Spurious high eccentricities can also be unphysical as they would lead to dynamical instabilities (e.g., Giuppone et al. 2013 We adopted uninformative priors on the RV semi-amplitudes K, zero points γ, jitter terms σ jit and slopesγ. As for the GP hyper-parameters, we imposed uniform priors with bounds wide enough to encompass the expected values of λ 1 , P rot and λ 2 , and only a lower bound of 0 m s −1 for h. Table 3 lists the adopted priors on the GP hyper-parameters for the systems that required GP regression to model the non-stationary activity variations along with the Keplerian signals. For three systems, namely Kepler-102, K2-141 and K2-132, we further imposed λ 1 > P rot /2, otherwise λ 1 converged to very low values, of the order of a couple of days, making the quasi-periodic term practically irrelevant, and hence unconstrained. Nonetheless, the RV semi-amplitudes obtained with and without this prior, i.e. by also allowing λ 1 to converge towards very low values, are fully consistent. Additional signals attributed to non-transiting planets were searched for in the RV residuals with GLS periodograms, and included in the DE-MCMC RV analysis if i) their FAP < 10 −3 , ii) their periodicity does not appear in the GLS periodograms of any activity indicator, and iii) the difference in the Bayesian Information Criterion (BIC) (Burnham & Anderson 2004;Liddle 2007) between the model with the additional planet and the model without it is ∆BIC > 10 ( Kass & Raftery 1995). Similarly, possible RV long-term slopes were considered significant if they are detected at more than 3σ and the ∆BIC in favor of the model with the slope is greater than 10.
We determined the values and the 1σ uncertainties of the model and derived parameters from the medians and the 15.87%-84.14% quantiles of their posterior distributions. For distributions consistent with zero, such as those of eccentricities or RV semi-amplitudes in cases of non-detection of the Doppler signals, we provided only the 1σ upper limits defined as the 0%-68.27% quantiles. In Table 4 we report for each system the HARPS-N systemic velocities and jitter terms, the GP hyperparameters in the presence of stellar activity signals, and the linear accelerations; in Tables 6 and 9 the parameters of the nontransiting cold Jupiters and low-mass planets, respectively; in Table 10 the parameters of the 64 Kepler and K2 transiting SPs in our sample (see Sect. 4).

Survey completeness
To determine the frequency of cold Jupiters in our sample of Kepler and K2 systems, we first need to evaluate the sensitivity of our survey to the presence of such planets. Indeed, CJs might not have been detected in some systems because of limited temporal baselines, poor temporal sampling and/or relatively low precision of the RV measurements. Therefore, our measure of the occurrence rate of CJs has to take the completeness (or recovery rate) of our survey into account.
The completeness can be estimated with experiments of injection and recovery of planetary signals for each system, by considering the real times of the observations and the RV uncertainty of each measurement at time t, i.e. σ(t) = σ 2 RV (t) + σ 2 jit . Following Bryan et al. (2019), we simulated signals of CJs in a logarithmic grid of 30x30 cells of planetary mass ∆M p vs semimajor axis ∆a, covering the ranges 0.3 − 20 M Jup in M p and 0.5 − 20 AU in a. For each cell of a given system, we simulated 300 RV signals of CJs at the epochs of our RV observations by randomly varying i) M p and a within the cell bounds, ii) T c within the orbital period corresponding to a and the stellar mass M ⋆ (Table 1) from Kepler's third law; iii) cos i from 0 to 1, where i is the orbital inclination; iv) the argument of periastron ω from 0 to 2π; and drawing the orbital eccentricity e from a Beta distribution (Kipping 2013). We then shifted every RV point at time t according to a Gaussian distribution with mean equal to the RV value and standard deviation σ(t). For simplicity, we did not simulate stellar magnetic activity signals, assuming that those signals would be efficiently modeled with GP regression, as shown in Sect. 4.
To establish the recovery rate in every ∆M p -∆a cell, we fitted the injected signals with a slope, a quadratic trend and a Keplerian orbit (with input parameters close to the simulated ones), and compared these models with a constant model (i.e., no signal) through the ∆BIC criterion: if ∆BIC > 10 in favor of the model with the planet-induced signal, we recorded a detection of the simulated signal, otherwise its non-detection. Figure 1 shows the completeness of two systems, Kepler-1876 (top left panel) and Kepler-93 (top righ panel), as one of the worst and best cases in our sample, respectively, and the average completeness of the 37 systems (bottom panel).

Search for and characterization of cold Jupiters
We first present the results of the search for cold Jupiters in our survey.

Systems with no long-term trends
Thirty-one out of thirty-eight systems do not show any significant long-term trend as caused by sufficiently massive outer companions within the completeness limits (Table 4). The RVs of two of these systems, K2-110 ) and Kepler-78 Howard et al. 2013), are shown in Figures 2 and 3 as representative cases. Possible extra noise was taken into account through the white noise term σ jit only in the K2-110 system, because the host star is not magnetically active. On the contrary, GP regression with a Quasi-Periodic kernel was needed to model the strongly correlated variations of the active star Kepler-78 ( Fig. 3; see also Grunblatt et al. 2015).  . The trends in Kepler-93, Kepler-454 and K2-12 are caused by bound companions, while those in K2-96/HD 3167 and K2-262/Wolf 503 are due to stellar magnetic activity, because similar trends are also seen in the activity indicators (see Fig. 6 for the case of K2-96/HD 3167; cf. also Bourrier et al. 2022). Gaussian Processes with a SE kernel were used to model the RVs of Kepler-93 to account for longterm variations overlapping with the linear trend and the signal of the planet Kepler-93b (Fig. 4). Such variations are likely due to stellar activity, as they seem to follow those of the S-index activity indicator after 7000 BJD TDB − 2450000 with a Pearson correlation coefficient of 0.26 ± 0.03 (Fig. 5, right panel).  From the temporal baseline, ∆T , of the RV time series and the amplitude of the trend, ∆K =γ · ∆T , we could estimate the minimum orbital period P and minimum M p sin i of the companions producing the trends, by assuming circular or quasi-circular orbits, i.e. P ≥ 4 · ∆T and M p sin i ≥ (∆K/28.4 m s −1 ) · (M s /M ⊙ ) 2/3 · (P/yr) 1/3 . Table 5   erating the linear trend in Kepler-93 is a brown dwarf or a lowmass star, having M p sin i ≥ 21 M Jup , while the masses of the companions producing the slopes in the Kepler-454 and K2-12 systems are currently compatible with a planetary companion, but further monitoring is needed to unveil their nature.

Systems with Keplerian signals of cold Jupiters
Resolved Keplerian orbits of cold Jupiters are observed only in three of the thirty-eight systems, i.e. Kepler-68 (Gilliland et al. 2013), Kepler-454 ) and K2-312/HD 86053 .
The two CJs Kepler-68d and Kepler-454c with P = 633 and 524 d and M p sin i = 0.75 and 4.51 M Jup were previously discovered by Gilliland et al. (2013) and Gettel et al. (2016). Moreover, a long-term quadratic trend in Kepler-68  and a slope in Kepler-454  were also found in the RV data, revealing the presence of additional outer companions in both systems, given that no trends were seen in the activity indicators. Our analysis not only allowed us to refine the parameters of the Kepler-68d and Kepler-454c giant planets, but also unveiled two additional CJs, i.e. Kepler-68e and Kepler-454d with P ∼ 3450 and 4070 d and minimum masses of 0.27 and 2.31 M Jup , respectively (see also Margini et al. in prep. for Kepler-68). The RVs of both systems as a function of time and the phase-folded RV signals of the CJs Kepler-68d, Kepler-68e, Kepler-454c and Kepler-454d are shown in Figs. 7 and 8.
The RV monitoring of the K2-312 system revealed that the linear trend observed by Frustagli et al. (2020) is due to a very eccentric CJ, K2-312c, with P = 921 d, M p sin i = 5.41 M Jup and e = 0.85 ( Fig. 9; Poretti et al. in prep.). The parameters of the five CJs with resolved Keplerian orbits and their 1σ uncertainties are given in Table 6.

Occurrence rate of cold Jupiters
By considering the three systems with resolved Keplerian orbits (Kepler-68, Kepler-454 and K2-312) and possibly the K2-12 system, which shows a trend that is currently compatible with a planetary companion at orbital distance a < 10 AU, we can derive the occurrence rate of cold Jupiters f CJ|SP in our sample. From this sample we have to remove Kepler-22 for the reason explained in Sect. 2, which yields a total of 37 Kepler and K2 systems. We used binomial statistics, i.e.
Article number, page 7 of 21  where d is the number of systems with detected cold Jupiters, i.e. d = 3 or 4 depending on whether K2-12 is included or not, and N ⋆,eff is not just the number of systems in our sample N ⋆ = 37, but the "effective" number of stars N ⋆,eff = N ⋆ · C, where C is the average completeness obtained by computing the mean of the completeness maps of the 37 systems (Fig. 1, bottom panel).
To compare our results with those of Wittenmyer et al. (2020), we computed C and f CJ|SP for the range in planetary mass ∆M p = 0.3 − 13 M Jup and different intervals of semi-major axis ∆a = 1 − 2, 2-4, 4-10, and 1-10 AU (see Table 7). Similarly, we also derived C and f CJ|SP for ∆M p = 0.5 − 13, 0.5-20 M Jup and ∆a = 1 − 10 and 1-20 AU, for comparison with the f CJ|SP found by Bryan et al. (2019). We report our values and 1σ error bars of the occurrence rates of CJs as well as those obtained by Wittenmyer et al. (2020) and Bryan et al. (2019) in Table 8. In summary, we found f CJ|SP = 9.3 +7.7 −2.9 % for d=3 and f CJ|SP = 12.3 +8.1 −3.7 % for d=4, i.e. by considering the long-term trend in K2-12 as planetary in origin, for ∆M p = 0.3 − 13 M Jup and ∆a = 1 − 10 AU. The former value is lower than that found by Wittenmyer et al. (2020), i.e. f CJ = 20.2 +6.3 −3.4 %, by a factor two, though, given the large uncertainties, the two measures are compatible at 1.3σ. Our f CJ|SP is instead four times lower than that derived by Bryan et al. (2019).   P e = 102.09 ± 0.52 d, which is likely due to a different treatment of the activity signal: we fitted a slope to all the RVs (see Fig. 6), while Bourrier et al. (2022) included in their MCMC analysis two activity-decorrelation terms for the HARPS and HARPS-N spectrographs. In any case, as noted by Bourrier et al. (2022), the observing spectral window allows for different solutions of P e with multiple peaks in the posterior distribution (see their Fig. 12).

Non-transiting low-mass planets
Article number, page 9 of 21 A&A proofs: manuscript no. Bonomoetal_2023_astroph  We report a new planet candidate, Kepler-1876c (P = 15.8 d, M p sin i = 11.0 M ⊕ ), in the single transiting system Kepler-1876 (Fig. 10). The GLS periodogram of the HARPS-N RVs of Kepler-1876 shows a significant periodicity at P = 15.8 d with FAP of 5.6 · 10 −5 , which does not appear in any of the activity indicators. Nonetheless, the ∆BIC in favor of the 2-planet model is currently ∆BIC = 5.3 < 10, and thus more RVs are needed to confirm the planetary nature of this signal, by also checking that its phase and amplitude do not change with time.
We cannot confirm the signals of the non-transiting planet Kepler-20g ) and the non-transiting candidate K2-2c  with periods P = 35 and 45 d and RV semi-amplitudes K = 4.1 and ∼ 2 ms −1 , respectively. The former signal is no longer present after the improvement of the DRS software to extract the HARPS-N RVs (Sect. 2). The latter was likely due to a combination of the spectral window of the original RVs in Vanderburg et al. (2015) and an activity signal with periodicity of ∼ 270 d, which is also seen in the GLS periodograms of the S-index and FWHM activity indicators, and was modeled with a GP-SE approach (see Table 4). That the 45 d signal could be originated by stellar activity variations was already discussed and considered by Vanderburg et al. (2015).
As specified in Sect. 3.2, we used non-interacting Keplerians to model RV planetary signals and GLS periodograms to search for non-transiting planets in the RV residuals. This implies that

Improved physical and orbital parameters of 64 Kepler and K2 small planets
As mentioned in Sect. 1, the DE-MCMC analyses of new HARPS-N RVs, in some cases combined for the first time with the literature RVs obtained with other instruments, allowed us to update the orbital and physical parameters of the 64 transiting planets. In particular, from the stellar parameters given in Table 1, the transit parameters (P, R p and i) from the literature 4 , and our updated RV semi-amplitudes K, we re-derived the planetary masses, densities and surface gravities for all the planets as in Table 10. We also report in the same table the planet equilibrium temperature T eq , assuming a null Bond albedo and full redistribution of heat from the dayside to the nightside (e.g., López-Morales & Seager 2007), and incident flux F p . Both T eq and F p were updated in case the stellar T eff was re-determined from our HARPS-N spectra (Table 1). Table 10 also emphasizes the fundamental contribution of HARPS-N in determining precise masses and densities to infer the composition of SPs. Specifically, the masses of 25 planets are determined with a precision higher than 5σ (M p /σ M p > 5), 13 of which with M p /σ M p > 7.5, and 8 with M p /σ M p > 10. Only upper limits on M p are given for 20 planets as their induced Doppler signal is not detected. The discussion of the planet compositions from the measure of planetary masses and densities is beyond the scope of this work 5 . 4 For the systems with updated R ⋆ values in Table 1, we recomputed R p from the R p /R ⋆ transit parameter and R ⋆ . In cases where the impact parameter b only is provided in literature assuming circular orbits, we derived the orbital inclination as i = arccos (b · R ⋆ /a). 5 In particular, we refer the reader to companion papers, which discuss in detail the updated compositions of Kepler-10b and c (Bonomo et al. in prep.), Kepler-68b and c (Margini et al. in prep.) and EPIC-229004835b (Tronsgaard et al. in prep.)

Discussion and conclusions
The value of the occurrence rate of cold Jupiters from the sample of 37 Kepler and K2 systems monitored by the HARPS-N/GTO consortium over the long term, i.e. f CJ|SP = 9.3 +7.7 −2.9 %, is lower than f CJ = 20.2 +6.3 −3.4 % derived by Wittenmyer et al. (2020) for solar-type stars from RV surveys, regardless of the possible presence of inner SPs (Table 8). This might hint at the theoretical anti-correlation between the presence of SPs and CJs predicted by Izidoro et al. (2015) and Lambrechts et al. (2019). However, the large uncertainty on our f CJ|SP associated with the inevitably limited target sample does not allow us to draw a firm conclusion. Moreover, the sample considered by Wittenmyer et al. (2020) likely contains a certain fraction of stars hosting SPs, which went undetected through the Doppler method, and thus a comparison of cold-Jupiter occurrence rates for stars with and without SPs is not straightforward. Nonetheless, assuming that the aforementioned anti-correlation does hold, the possible "contamination" of the sample of Wittenmyer et al. (2020) by smallplanet systems would tend to make f CJ lower than it really is (in other words, f CJ might be higher than ∼ 20% for solar-type stars hosting no SPs).
Our results do not support the claim by Bryan et al. (2019) about an excess of CJs in small-planet systems (see Table 8): according to their f CJ|SP , we should have discovered CJs in 12 ± 2 of our systems, while we only found them in 3 systems. The reason for such a discrepancy lies at least in part in an incorrect interpretation of the planetary origin of some linear trends by Bryan et al. (2019) (see their Fig. 3): for instance, the trend in the GJ 273 system is due to secular acceleration, the slope in HD 3167 to stellar activity (Fig. 6), and those in Kepler-93b and Kepler-407b are generated by brown-dwarf or low-mass stellar companions. In addition, f CJ|SP was computed by Bryan et al. (2019) in a different way, i.e. by assuming a double power law f CJ|SP = A · M α p · a β , deriving the coefficients A, α and β with a likelihood approach, and integrating f CJ|SP over the parameter space. Nonetheless, the exclusion of approximately half of the long-term trends found by Bryan et al. (2019) as not due to CJs, would reduce considerably their f CJ|SP and make it more compatible with our estimate.
The present work will be extended to a considerably larger sample, including the TESS small planets monitored with HARPS-N since 2019 (e.g., Cloutier et al. 2021;Lacedelli et al. 2021;Naponiello et al. 2022), as well as the K2 and TESS smallplanet systems observed from the Southern hemisphere with other spectrographs such as ESPRESSO and HARPS. Performing the same analyses of RV data in Sect. 3.2 on a sample at least three times as large will allow us to (i) derive a more precise f CJ|SP in small-planet systems.
(ii) compute f CJ|SP as a function of small-planet multiplicity. Indeed, the predicted anti-correlation between CJs and SPs is expected to be more pronounced for systems with a higher level of multiplicity of SPs, because single cores have a higher probability to jump inside the orbit of the cold gas giant during their inward migration (see Sect. 1.1 and Fig. 3 in Izidoro et al. 2015 System might have prevented the Uranus and Neptune cores from migrating towards the Sun. Nevertheless, the CJ pairs we discovered in the Kepler-68 and Kepler-454 systems did not hinder the formation of the inner SPs Kepler-68b/c and Kepler-454b. (iv) determine f CJ|SP as a function of the small-planet composition. Small planets would be mainly ice-rich (sub-Neptunes) in the scenario of Sect. 1.1 (cf. also Zeng et al. 2019) or rocky (super-Earths) in the scenario of Sect. 1.2, as they are expected to form beyond or inside the water iceline, respectively. Unveiling a possible anti-correlation between CJs and predominantly ice-rich or rocky SPs may thus yield important clues on the mechanisms of formation of shortperiod SPs. An anti-correlation between rocky super-Earths and cold Jupiters would not be expected instead, if the former mainly originate inside rings of silicate-rich planetesimals at approximately 1 AU and then migrate inward (Batygin & Morbidelli 2023). Actually, a positive correlation could exist if the mass initially in the ring of silicaterich planetesimals, leading to rocky super-Earths, is correlated with the mass of the ice-rich planetesimals, leading to the formation of the cores of outer giant planets. In this regard, the improvement in the mass determination of several of the 64 Kepler and K2 transiting SPs from our RV analyses (Table 10) is useful to distinguish between rocky and non-rocky compositions.
For these purposes, we recommend continuing to follow up small-planet systems, even after achieving the desired precision on planetary masses for the investigation of small-planet compositions. This can be done at a very low cost, given that just a few RV measurements over the years would in principle be sufficient to search for outer giant planets, provided that the spectrograph is stable. For the brightest targets, additional and complementary information on the presence of cold Jupiters will also be provided by the Gaia mission (e.g., Holl et al. 2022 Article number, page 12 of 21 Table 1. Kepler and K2 systems in our sample. From left to right the columns report the name of the system, the multiplicity of transiting planets, the stellar parameters (mass, radius, effective temperature, metallicity and isochronological age), the number of radial velocities from all surveys (N RV tot), the number of HARPS-N radial velocities (N RV HN), the number of radial-velocity datasets (N Dat ), the total duration of the radial-velocity time series, and the literature references for both the stellar parameters and the radial-velocity measurements.   System  Table 3. Priors imposed on the hyper-parameters of the Gaussian Process regression for the stars showing significant magnetic activity variations according to the criteria discussed in the text. The hyper-parameters of the covariance function are the radial-velocity semi-amplitude h, the exponential decay time λ 1 , the inverse harmonic complexity term λ 2 , and the rotation period P rot . The second column indicate the use of either Squared Exponential (SE) or Quasi-Periodic (QP) kernel for the covariance function, the former having just two hyper-parameters, i.e. h and λ 1 . U stands for uniform (uninformative) prior.  Table 4. Systemic radial velocities (γ), jitter terms (σ jit ), Gaussian Process hyper-parameters -namely radial-velocity semi-amplitude h, exponential decay time λ 1 , inverse harmonic complexity term λ 2 and rotation period P rot -, and linear accelerations (γ) from the DE-MCMC radial-velocity modeling. Even though, for simplicity, γ and σ jit are reported for the HARPS-N spectrograph only, the radial-velocity modeling was carried out by including the data gathered with other spectrographs, when available. Accelerationsγ evaluated as significant, from both a confidence level higher than 3σ and an odds ratio greater than 10 from the ∆BIC in favor of the model with acceleration, are highlighted in boldface. HN-1 and HN-2 refer to the HARPS-N data collected before and after September 2012, respectively. In the last column we report comments about the origin of the detected trends, such as possible cold Jupiter (CJ), brown dwarf (BD) or low-mass star (LMS) companions or stellar magnetic activity, and about the presence of cold Jupiters in resolved orbits (ro), which were modeled with Keplerians and are thus not responsible for the trends. The radial-velocity zero point is close to zero because the radial velocities were extracted with the Yarara-v2 tool (Cretignier et al. 2021) from the HARPS-N spectra reduced with the DRS-v2.3.5 (see Bonomo et al. in prep.).    -39 ± 7 1 f CJ|SP for 1 < M p < 10 M ⊕ instead of the wider range 1 < M p < 20 M ⊕ used in both this work and Zhu & Wu (2018).
Article number, page 17 of 21 Table 10. Orbital and physical parameters of the transiting Kepler and K2 planets. From left to right the columns report the planet name, the transit mid-time, the orbital period, the planetary radius, the orbital inclination, the references for the transit parameters, the orbital eccentricity, the radial-velocity semi-amplitude, the planet mass, density and surface gravity, the semi-major axis, the equilibrium temperature by considering a null Bond albedo and full heat redistribution from the day to the night side, and the stellar incident flux.