Prospects from TESS and Gaia to Constrain the Flatness of Planetary Systems

The mutual inclination between planets orbiting the same star provides key information to understand the formation and evolution of multiplanet systems. In this work, we investigate the potential of Gaia astrometry in detecting and characterizing cold Jupiters in orbits exterior to the currently known Transiting Exoplanet Survey Satellite (TESS) planet candidates. According to our simulations, out of the ∼3350 systems expected to have cold Jupiter companions, Gaia, by its nominal 5 yr mission, should be able to detect ∼200 cold Jupiters and measure the orbital inclinations with a precision of σcosi<0.2 in ∼120 of them. These numbers are estimated under the assumption that the orbital orientations of the CJs follow an isotropic distribution, but these only vary slightly for less broad distributions. We also discuss the prospects from radial velocity follow-ups to better constrain the derived properties and provide a package to do quick forecasts using our Fisher matrix analysis. Overall, our simulations show that Gaia astrometry of cold Jupiters orbiting stars with TESS planets can distinguish dynamically cold (mean mutual inclination ≲5°) from dynamically hot systems (mean mutual inclination ≳20°), placing a new set of constraints on their formation and evolution.


INTRODUCTION
Over 5,100 exoplanets have been confirmed (Akeson et al. 2013), the majority of which were discovered through transits or radial velocities.Around 380 of the known exoplanets were discovered by the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015), and it has ∼ 5,900 candidates yet to be confirmed.These known exoplanets have revealed rich information about the occurrence rate, architecture, and theoretical implications of the planetary systems in general (see a recent review by Zhu & Dong 2021).
Compared to other detection methods such as transits and radial velocities, astrometry has a controversial past as nearly all claimed planet detections have been dis-Corresponding author: Juan I. Espinoza-Retamal jiespinozar@uc.clcarded by subsequent measurements 1 (e.g., Bean et al. 2010).However, this panorama is expected to be changed in the next few years as the Gaia astrometry mission is going to release about 20,000 giant planet detections with its upcoming data release (DR) 4 (Perryman et al. 2014).In fact, with the release of the Gaia DR3 (Gaia Collaboration et al. 2022a), we already have dozens of astrometric candidates in the substellar regime (see, e.g., Gaia Collaboration et al. 2022b), and a few systems have been identified using Hipparcos and Gaia astrometry, and confirmed by direct imaging (e.g., Currie et al. 2023;Mesa et al. 2023;De Rosa et al. 2023).Astrometry will be especially useful as it can provide us measurements of the orbital inclinations and true masses of planets (see, e.g., Sozzetti et al. 2001;Casertano et al. 2008;Sozzetti et al. 2014;Perryman et al. 2014).For example, with data from the Gaia EDR3 (Gaia Collab-1 According to NASA Exoplanet Archive, the only exoplanet discovered using astrometry is DENIS-P J082303.1-491201b (Sahlmann et al. 2013) but due to its high mass of ∼ 28 M J , it is debatable if this object can be regarded as a planet or a brown dwarf.
oration et al. 2021), Brandt et al. (2021a) measured the true mass of the planet HR 8799 e, and thanks to this, they estimated its age at ∼ 42 Myr.Combining radial velocities with Hipparcos and Gaia astrometry, Venner et al. (2021) measured the orbital inclination and the true mass of the companion to the star HD 92987.They found that, in fact, that object was not a planet but rather a star of ∼ 0.2 M ⊙ in a nearly-polar orbit.
There are multiple pieces of evidence suggesting that planetary systems are not always as flat.Some protoplanetary disks exhibit significant internal misalignments, either warps or disks broken in pieces with different orientations, as evidenced by multiple observations, including scattered light observations (shadows; e.g., Casassus et al. 2018), gas kinematics (e.g., Marino et al. 2015), dust emission from ALMA images (e.g., Francis & van der Marel 2020), and periodic light extinction caused by dusty disks (e.g., Ansdell et al. 2016).Also, we have found planets orbiting stars with large obliquities (angles between the host star's equator and the planetary orbit).This includes planets from nearly polar to fully retrograde orbits as measured for transiting exoplanets from spectroscopy (see the review by Albrecht et al. 2022), with the Rossiter-McLaughlin effect (e.g., Lendl et al. 2014), spot-crossing events (e.g., Sanchis-Ojeda et al. 2013), stellar rotation (e.g., Winn et al. 2017), and stellar variability (e.g., Mazeh et al. 2015;Li & Winn 2016).On the population level, statistical studies of the planetary systems found by the Kepler transit survey have suggested that a large fraction of the mature planetary systems probably have substantial mutual inclinations, as revealed from the observed planet multiplicity distributions and timing of the transits (Zhu et al. 2018;He et al. 2020;Millholland et al. 2021).
More recently, using radial velocities and astrometry, both Xuan &Wyatt (2020) andDe Rosa et al. (2020) measured the orbital inclination of the cold Jupiter (CJ) in the π Men system (Jones et al. 2002;Huang et al. 2018), and combining that with TESS data of this system they found a large mutual inclination between the transiting super-Earth and its outer giant companion.From this type of measurement, a set of question that motivate our work arise: how many more π Men-like systems will we find?More concretely, for how many planetary systems that have been discovered with TESS will we be able to measure the mutual inclination between the transiting planet and its possible outer companion using astrometry?How can we best exploit these upcoming datasets to understand the evolution of planetary systems?How important are radial velocity followups to better constrain the parameters?

Mutual Inclinations and formation histories
Mutual inclination measurements can give us indications of past interactions that happened to form the architectures of planetary systems that we see today.These interactions range from violent giant impacts or gravitational scattering (e.g., Huang et al. 2017;Gratia & Fabrycky 2017;Mustill et al. 2017;Pu & Lai 2021) to long-term chaotic diffusion (e.g., Wu & Lithwick 2011;Hamers et al. 2017;Petrovich et al. 2019).
By measuring mutual inclinations in systems with a transiting planet and its outer companion, we may constrain their formation pathway.For instance, in systems composed of two gas giants, including a transiting hot or warm Jupiter (HJ/WJ) and a CJ, their mutual inclinations can constrain the migration mechanism.If the migration was produced by angular momentum exchanges with the protoplanetary disk (e.g., Goldreich & Tremaine 1980;Ward 1997;Baruteau et al. 2014), we should expect low mutual inclinations.In turn, if the migration was produced by high-eccentricity migration (e.g., Rasio & Ford 1996;Wu & Murray 2003;Petrovich & Tremaine 2016), we generally expect high mutual inclinations2 .
Other systems of interest are the short-period transiting super-Earth/mini-Neptune (sub-Jovians, SJs) and outer cold Jupiters.As eccentricities in these systems are generally small due to tidal circularization (or stability considerations) and/or hard to constrain by radial velocities due to their low masses (e.g., MacDougall et al. 2021), we may gauge the level of dynamical upheaval using mutual inclinations.

Structure
In this paper, we estimate the number of TESS Objects of Interest (TOIs) for which Gaia astrometric observations should detect an outer companion and the number of those that will have a well-constrained orbital inclination.In Section 2, we describe the methodology used for the simulations.In Section 3, we present the results.In Section 4 we discuss how much more we can improve the results if we add information from radial velocity (RV) measurements.In Section 5, we discuss how our results will change with model assumptions, especially the underlying mutual inclination distribution.We conclude in Section 6.

METHODS
We use the TOI catalog that was obtained from the Exoplanet Follow-up Observing Program for TESS (ExoFOP-TESS) 3 on August 23, 2023.Although there will be more detections from the ongoing TESS extended mission and dedicated searches from transit signal, a significant fraction of the identified TOIs are, or will be, false positives and thus not transiting planets (e.g., Guerrero et al. 2021).Thus, our catalog suits the purpose of the present work, namely, to estimate the number of planetary systems with detections from both TESS transit and Gaia astrometry.We did not consider in the analysis TOIs without reported stellar mass or planetary radius.Also, we did not consider stars with masses greather than 2 M ⊙ to avoid unreliable measurements.By applying these filters, we end up with 5,864 TOIs from 5,625 unique stars.
The probability of having an exterior cold Jupiter depends on the properties of the inner planet.In this work, we adopt the following conditional probabilities (Zhu & Wu 2018) 0.75 if in = HJ (Bryan et al. 2016) 0.49 if in = WJ (Bryan et al. 2016) (1) Here "SJ", "HJ", and "WJ" stand for sub-Jovian, hot Jupiter, and warm Jupiter, respectively.We classify the TOIs into these categories based on the measured planet size and semi-major axis: a SJ if R p,in < 6 R ⊕ , a HJ if R p,in > 6 R ⊕ and a p,in < 0.1 au, and a WJ if R p,in > 6 R ⊕ and a p,in > 0.1 au.
While the conditional probabilities given above are derived from observations, those studies also report nonnegligible uncertainties around these benchmark values.Furthermore, different studies also reported different values for these conditional probabilities.For example, the conditional rate of CJs on inner SJs is reported to be lower in Bonomo et al. (2023) (but see Zhu 2023).These uncertainties on the conditional probabilities will affect the expected numbers of the CJ detections, so the exact number of detections will be useful to further refine the conditional probabilities.For constraining the flatness of the planetary systems, we expect the results of different mutual inclination distributions to be affected in the same way, so our result on the mutual inclination distribution may remain largely unaffected.
3 https://exofop.ipac.caltech.edu/tess/We injected the signal of a CJ into each TOI and attempted to recover it using Gaia astrometry in order to assess whether could detect the CJ and the precision with which we could measure the inclination of its orbit.We assumed that each Gaia measurement would have 1-D astrometric precision σ fov , which only depended on the magnitude G of the star (Perryman et al. 2014).To obtain realistic estimates of the times in which Gaia will observe each star, we used the HTOF tool (Brandt et al. 2021b).Epochs taken before January 25, 2020, are considered in order to have a close match with the upcoming Gaia DR44 .We randomly reject 20% of the Gaia epochs because this fraction of Gaia observations is shown to be problematic due to satellite dead times, unusable observations, or observations rejected as astrometric outliers (see, e.g., Lindegren et al. 2018;Boubert et al. 2020;Brandt et al. 2021b).After applying these rejections, we obtained realistic epochs for 3,350 unique stars.According to Perryman et al. (2014), the number of measurements is primarily dependent on the ecliptic latitude of the target, so we divided stars in bins of 5 • based on this value and selected the TOI with the median value of observations in each bin.We use the epochs of this median TOI in each bin as the epochs for the remaining stars in the same bin without epochs.The HTOF tool can also give the scanning law of the Gaia satellite, but for simplicity, we do not use this information.Instead, we include exactly half of the two-dimensional information of the astrometric measurements in the Fisher matrix analysis.See Appendix A for details.
For the injected CJs, their physical and orbital properties were randomly sampled from the following distributions: • The mass-ratio q ≡ M p /M ⋆ follows a broken power-law distribution with a break at q break = 1.7 × 10 −4 .The power-law indexes above and below the break were -0.93 and 0.6, respectively (Suzuki et al. 2016).We worked with planetary masses between 0.3 and 15 M J .The lowest mass ratio used in our simulations was ∼ 1.4 × 10 −4 when M ⋆ ∼ 2 M ⊙ .
• The orbital period P follows a broken power-law distribution with a break at P break = 1717 days.The power-law indices above and below the break were -1.22 and 0.53, respectively (Fernandes et al. 2019).We worked with periods between 100 and 10000 days (∼ 0.27 − 27.4 yrs).
• The orbital eccentricity e follows a Beta distribution with parameters a = 1.12 and b = 3.09 (Kipping 2013).
• The orbital inclination i is uniform in cos i between 0 and 1.
• The argument of periapsis ω and the longitude of ascending node Ω both follow a uniform distribution between 0 and 2π.
Once the properties of the injected CJs were known, we modeled their astrometric signals in the standard way.Specifically, the astrometric motion of the host star along two perpendicular directions is given by (2) Here A, B, F, G are the so-called Thiele-Innes elements: where ρ is the semi-amplitude of the astrometric motion that can be written in terms of the mass-ratio q, semimajor axis a, and stellar distance d as: The eccentric anomaly E is related with the mean anomaly M by: The mean anomaly is defined as: For a chosen reference time t 0 , the astrometric motion can be modeled by a set of 9 parameters: the systemic velocities µ x and µ y , the semi-amplitude of the astrometric motion ρ, the orbital period P and eccentricity e, the reference position of the planet M 0 , and the three angles of orientation of the orbit ω, cos i and Ω.Note that we choose cos i instead of i just for simplicity.We choose not to perform a joint modeling of the stellar parallactic motion because parallax is much better determined and not correlated with the binary astrometric motion in the frequency domain.For many of the stars studied here, other means of distance determination may be available to further improve the parallax determination.
We use the Fisher matrix analysis to evaluate the detectability of the astrometric signal and the uncertainties on individual model parameters.This approach is more computationally efficient than a Markov chain Monte Carlo (MCMC) approach by a factor of ∼ 3,000.The details of the fisher matrix analysis are given in Appendix A. For each TOI, we carried out 10 4 simulations and considered that the outer giant was detected if ρ/σ ρ > 3 and that the orbital inclination was well constrained if σ cos i < 0.2.This implies an uncertainty of ∼ 11 • if the orbit is edge-on and ∼ 34 • if the orbital inclination is 20 • .

RESULTS
For each TOI, we obtained a distribution for the uncertainty in ρ and calculated the probability of detecting the CJ if it exists (i.e., ρ/σ ρ > 3).From Equation 1, the probability of the existence of the CJ is related to the type of planet that exists in the inner part of the system.The total number of CJs that should exist around TOI hosts is estimated to be: The number of these CJs that could be detected using Gaia astrometry is then As shown in Figure 1, the probability of detecting the CJ is a strong function of the stellar distance, and the probability is higher for nearby (≲ 100 pc) M-dwarfs.About half of the CJs will be detected in systems with SJs, whereas the remaining half in systems with giant planets (HJ or WJ).
From the distribution obtained for σ cos i , we also calculated the probability of having the inclination well constrained (i.e., σ cos i < 0.2) for each TOI system.With this information, we then estimated the number of CJs that would have the inclination well constrained The distribution of the size of the inner transiting planets is shown in the left panel of Figure 2. According to our definitions of small and large planets, 72 and 46 of the CJs with inclination measurements are from systems with SJs and HJs/WJs, respectively.These numbers

COMPLEMENTARY RADIAL VELOCITIES
In the astrometry method, orbital parameters describing the sky-projected motion of an elliptical orbit can be correlated.Specifically, the orbital inclination is correlated with several of the other parameters, of which the most important one is the astrometric amplitude ρ mainly due to the planet mass.Therefore, additional constraints on planet properties can help improve the constraints on the orbital inclination.Here, we assess to what level our results improve by adding information from complementary RV observations.
In appendix B, we show how our Fisher matrix analysis is modified to obtain the uncertainties in a model combining astrometry and RV measurements.This model parameters increases to 10: the previous nine from Section 2 and the systemic velocity in the z-axis (the line-of-sight direction), µ z .
In Figure 3 we show three examples of how the uncertainty in the inclination improves as the number of radial velocity measurements increases and for different representative precisions.We assume that the radial velocity measurements are taken uniformly over the 5 years after the last epoch of the astrometric observations.Also, we assume that the signal of the transiting planet was removed from the radial velocities, and the only signal present is the one of the CJ.Because RV observations alone provide no information on the orbital inclination, the best constraint one can achieve on the orbital inclination is limited by the information available from Gaia astrometry.As a result, there is a theoretical limit on the statistical uncertainty of the cos i parameter.This limit is given by (see Appendix C): There are a few things to notice from Figure 3. First, supplementary RVs will always be useful, even in systems that can be well-constrained by astrometric observations.Second, for systems that cannot be wellconstrained by astrometry, supplementary RVs can be crucial in confirming the planet signal and refining the system configurations.In fact, as the middle and right panels indicate, the orbital inclination can be much better constrained with only a few RV observations.Last but not least, RV observations with higher precision are always better.Since all the analysis will depend on the campaign and instruments chosen to carry out the follow-up we decided to make public a python script called Fisher for astrometry and RV5 with which it is possible to estimate the uncertainties that would be obtained for a system using the methodology described in Appendices A and B. We expect that the code will help observers know the precision level they will achieve in the parameters of a given system if they try different observing strategies.

DISCUSSION
Our simulations show that if the orbital inclination of the CJs is isotropic, Gaia should detect CJ companions in ∼ 206 TOI systems out of the over 5,600 TOI targets.A CJ is considered detectable if its astrometric amplitude is three times the per-measurement uncertainty, namely ρ/σ ρ > 3.Among these CJ detections, we expect that 118 will have well-constrained orbital inclinations (i.e., σ cos i < 0.2).The majority of CJs with well-constrained inclinations are found in systems with inner sub-Jovian planets, and nearby M-dwarfs are preferred for CJ detections and inclination measurements.
Additionally, we find that complementary RVs will always be useful, even in systems that can be well constrained by astrometric observations.For systems that cannot be well constrained by astrometry, complementary RVs can be crucial in confirming the planet signal and refining the system configuration.RV observations with higher precision require fewer measurements to improve the precision in parameters of the planet.

Comparison with previous works
Several studies have investigated the potential of Gaia astrometry in exoplanet study, including a few that looked into its capability of constraining the mutual inclination.Sozzetti et al. (2001) evaluated the capability of Gaia to detect planets around solar-type stars in the Solar neighborhood.Using the ν And system as the case of their study, they conclude that Gaia should be able to detect the outer two planets in the system and provide estimates of the full set of orbital elements accurate to better than 1 − 10%.Casertano et al. (2008) studied in more detail the detectability of planets around FGK dwarfs, finding that under favorable orbital configurations (both planets with P ≤ 4 yr and ρ/σ fov ≥ 10) Gaia could measure their orbital elements to better than 10% accuracy in more than 90% of the time.Using a Galaxy model (Besançon, e.g., Robin et al. 2003) their estimated yield is ∼ 8, 000 Gaia-detected planets and ∼ 4, 000 of them with accurately measured orbital parameters, including inclinations.Sozzetti et al. (2014) extended that study to close M-dwarfs concluding that in a sample of ∼ 3,150 M-dwarfs within 33 pc, Gaia should detect ∼ 100 CJs and almost all of them with good quality orbital solutions.Also, as mentioned in the introduction, Perryman et al. (2014) estimated that ∼ 20,000 giant exoplanets should be detected using Gaia astrometry.
Similar to these previous works, we also studied the capability of Gaia in detecting planets and measuring orbital inclinations, but now for a sample of stars in which we know, thanks to TESS, that there are transiting planets at close-in orbits.The advantage of trying to measure orbital inclinations in those systems is that we can put constraints on the mutual inclination between the transiting planet and its outer companion, allowing us to explore the parameter space.With Gaia alone, one Figure 3. Uncertainty in inclination as a function of the number of RV data taken for 3 fixed systems.Left: A system detectable with only astrometry and with the inclination well constrained (σ astro cos i ≈ 0.003).Center: A system detectable with astrometry but with the inclination not well constrained (σ astro cos i ≈ 0.45).Right: A system not detectable with only astrometry (σ astro cos i ≈ 1.5).Different colors represent different precisions for the instrument used to measure the RV.The black dashed line corresponds to the analytic limit in Equation 10 reachable in each case.
can only detect and measure orbital inclinations of the relatively long-period planets, whereas, with Gaia and TESS combined, one can constrain the mutual inclinations between planets in the inner and the outer parts of the system, which are likely related (e.g., Masuda et al. 2020;Zhu & Dong 2021).

Constraining the flatness of planetary systems
The astrometry method is more sensitive to more massive planets at relatively large orbital distances.If there is a second planet in the system detected with transits, we can constrain the mutual inclination between planets, i mut , defined as: cos i mut = cos i in cos i CJ + sin i in sin i CJ cos (∆Ω), (11) where i in and i CJ are the orbital inclinations of the inner planet and the Gaia CJ, respectively.In deriving the mutual inclination, we assume that the difference in longitudes of ascending nodes, ∆Ω, follows a uniform distribution between 0 and 2π.
Until now, the inclination of the CJ has been assumed to follow an isotropic distribution (see Section 2), and thus the mutual inclination also follows an isotropic distribution.To see if we could distinguish between isotropic and, for example, Rayleigh distributions for the mutual inclination, we repeated the same simulations but considering that the mutual inclination followed a Rayleigh distribution with σ = 5 and 20 • (hereafter R5 and R20).Using equation ( 11) and setting i in = 90 • (transiting), we obtained a new distribution for the inclination of the CJs.With these new distributions, we re-run the simulations and obtained that we should detect 191 and 202 CJs companions to TOIs for R5 and R20, respectively, compared to 206 in the isotropic case.Out of these detections, we expect to have the inclination well constrained for 149 and 121 of them for R5 and R20, respectively, compared to 118 in the isotropic case.In other words, because of the correlation in orbital inclinations between inner and outer planets that forces the CJ to have more inclined orbits (more edgeon), it becomes slightly more difficult to detect the CJ, but once detected it is easier to measure its inclination.
We generated random samples following those distributions (Uniform, R5, and R20) with their respective number of inclinations well constrained (118, 149, and 121) to compare them (see an example in Figure 4) via Kolmogorov-Smirnov (KS) tests.In a single KS test, the null hypothesis was that the two samples were drawn from the same underlying distribution.We set the threshold to be p > 0.05 if the hypothesis is to be accepted.Based on 100 simulations, we find the null hypothesis can always be rejected for KS tests between the R5 model and any of the other two models, whereas the null hypothesis is rejected 90% of the time for KS tests between the R20 and the Uniform models.We conclude that, with the expected numbers estimated in this paper, we will always be able to distinguish between R5 and the other 2 and between R20 and the Uniform models most of the time.The conclusion remains unchanged even if the numbers of inclinations well measured are all cut half.
If we restrict our sample to two gas giants-namely a transiting HJ/WJ and a Gaia CJ-we expect to have 62, 48, and 46 systems with well-constrained inclinations for the R5, R20, and uniform inclination distributions, respectively.These numbers allow us to always distinguish between R5 and the other two distributions, as well as between R20 and the uniform distribution most of the time.If the numbers of well-measured inclinations are cut half, R20 and Uniform will be distinguishable only ∼ 30% of the time.In turn, if we restrict our samples to an inner SJ and a CJ, we should have the inclination of the CJ wellconstrained for 87, 73, and 72 systems if the mutual inclination follows R5, R20, or uniform, respectively.Similar to the whole sample and to the case of HJs/WJs, with those numbers, we will always be able to distinguish between R5 and the other 2 and between R20 and the uniform models most of the time.If the numbers of well-measured inclinations are cut half, R20 and Uniform will be distinguishable only ∼ 30% of the time.

Caveats
In this work, we studied the capability of Gaia to detect CJs in the current population of TOIs with the idea of constraining the mutual inclination between the transiting planet and its outer giant companion.A strong correlation between the inner planets and the outer giant ones has been adopted in our work.Although a correlation is supported by several pieces of observational evidence, there is an ongoing debate regarding the strength of this correlation and whether it should apply to all types of stellar host types (e.g., Bryan et al. 2016;Zhu & Wu 2018;Bryan et al. 2019;Herman et al. 2019;Masuda et al. 2020;Rosenthal et al. 2022).This leads to an additional source of uncertainty in the derived numbers of systems with mutual inclination measurements.We will not explore this uncertainty further in the current work, as our primary goal is to investigate the power of Gaia in constraining the flatness of the planetary system.Nevertheless, it is worth noting that the number of actual detections should provide useful constraints on the strength and generality of the inner-outer correlation as well.
Also, we have not considered the possibility that the same systems contain additional planets and their impact on our results so far.In principle, there could be planets that are either undetectable or marginally detectable, such as in the case of π Men, where recently Hatzes et al. (2022) revealed the presence of a third planet on a 125-day orbit.Because only CJs are detectable with Gaia observations, only the presence of a second CJ in the system can affect the measurements of parameters for the detected planet.But, as we argue next the signal contamination from these potential second CJs is expected to be low.
Ground-based RV observations have enabled studies of the fraction of systems with multiple CJs.Recent work by Zhu (2022) analyzed the California Legacy Survey data (Rosenthal et al. 2021) and derived the intrinsic multiplicity distributions of different planet classes.According to that study, about 27% of CJ systems have at least two CJs.This serves as a theoretical upper limit if one is to estimate the fraction of Gaia CJ systems with multiple planet detections.Furthermore, considering that the ground-based RV surveys have better coverage in the planet mass-semi-major axis plane than Gaia astrometry, the above upper limit can be further reduced.According to Zhu (2022), there are only eight two-CJ systems out of the 49 systems with CJs in the CLS sample.This puts an upper limit of ∼ 16% on the fraction of CJ systems with multiple CJ detections in the Gaia sample.
From a theoretical point of view, the presence of two giant planets in the same system may be unstable.The star with the median probability of detecting the CJ in this study was TOI-5612 and the typical planet detected here was a ∼ 9 M J at 3.3 AU with an orbital eccentricity of 0.2.Using this planet as the one detected, we studied the stability of the system if there was another CJ drawn from the same population.From a population of 100,000 CJs, and using the stability criterion from Petrovich (2015b), we found that only ∼ 20% of the simulated two-planet systems are stable.Furthermore, only in ∼ 10% of the stable systems does the second planet produce a comparable astrometric signal compared to the first planet.Therefore, the fraction of systems that will be affected by planet multiplicity is small.We leave a detailed study of these multi-planet systems to some future study.

CONCLUSIONS
We have performed injection-recovery simulations of the Gaia astrometric observations for the current sample of TOIs (5,625) in order to estimate the detection yields of CJs in these systems as well as their sky-projected inclinations, thereby constraining the mutual inclination between the transiting planet and its outer companion.We find the following results: • Under the assumption that the mutual inclination distribution is isotropic, out of the estimated 3,340 TOIs with CJ companions, Gaia should detect 206 and have the inclination well constrained for 118 of them.Nearly ∼ 60% (72/118) of these correspond to TOIs with sub-Jovian size candidates.
• If the mutual inclination follows a Rayleigh distribution with σ = 5 • and 20 • (R5 and R20), Gaia should detect 191 and 202 CJs and have the inclination well-constrained for 149 and 121 of them, respectively.With those numbers, we can confidently distinguish between the R5 model and the models with broader distributions (R20 and Uniform), while R20 and the Uniform models can be distinguished most of the time.These conclusions remain unchanged even if the numbers of wellmeasured inclinations are all cut by half.
• The uncertainties in the CJ inclinations can be reduced significantly if complementary RV observations are taken on the Gaia targets.This is especially true for systems in which astrometry alone provides a poor constraint.The RV followup strategy should be assessed on a case-by-case basis.We provide a Python script to compute the expected uncertainties using our Fischer matrix formalism quickly.
Overall, our simulations show that Gaia's astrometric measurements of planet-hosting stars from TESS will constrain the flatness of systems hosting inner transiting planets and outer cold Jupiters at levels that can distinguish dynamically cold (mean mutual inclination ≲ 5 • ) from dynamically hot system (mean mutual inclination ≳ 20 • ), thus placing a new set of constraints on their formation and evolution.
Given a series of measurements at (t 1 , t 2 , ..., t N ), each with a precision σ fov for the astrometric motion, and a series of measurements at (t ′ 1 , t ′ 2 , ..., t ′ L ), each with a precision σ RV for the RV, the individual element of the Fisher Matrix is given by Here ⃗ θ = (µ x , µ y , µ z , ρ, P, e, M 0 , ω, cos i, Ω).

C. ANALYTICAL APPROACH
The key information from RV is its constraint on the minimum mass ρ sin i through the RV semi-amplitude where K 0 is a quantity that captures the dependence of K on the other parameters such as orbital period, eccentricity, etc.The Fisher information matrix for the selected variables, (ρ, cos i), is then Here we have made use of the following relations If the RV semi-amplitude is measured to a fraction precision of ϵ K the matrix is further simplified to The above Fisher information matrix should be added to the corresponding rows and columns of the full Fisher matrix from astrometry, to determine the final uncertainties on the Keplerian parameters.Below we provide a simplified version to gain some insight.
From astrometry alone, one can determine the covariance matrix of (ρ, cos i), which we denote as Here σ ρ and σ cos i are the uncertainties on ρ and cos i in the astrometry-alone case, respectively, and r is the correlation coefficient between ρ and cos i.The corresponding Fisher information matrix is Combining the Fisher information of RV and astrometry, we have where we have rewritten σ ρ = ρϵ ρ .The covariance matrix in the RV+astrometry joint case is , (C19) where det The ratio between the new uncertainties on ρ and cos i is After some rearrangements, this yields If RV provides good enough constraints on K, the above equation reduces to This is rewritten in the form of Equation ( 10) with the introduction of σ new ρ ≡ ρϵ new ρ .

Figure 1 .
Figure 1.Left: Histogram of the number of cold Jupiters that should be detected as a function of the size of the inner planet.Right: Scatter plot of stellar distance vs. stellar effective temperature of all TOIs.Color represents the probability of detecting the cold Jupiter.

Figure 2 .
Figure 2. Left: Histogram of the number of cold Jupiters that should have the inclination well constrained as a function of the size of the inner planet.Right: Scatter plot of stellar distance vs. stellar effective temperature of all TOIs.The color represents the probability of having the inclination of the cold Jupiter well constrained.suggest that systems like π Men, which has mutual inclination contrained between the inner transiting super-Earth and the outer CJ (Xuan & Wyatt 2020; De Rosa et al. 2020), will not be uncommon.In terms of the stellar properties, the probability is higher for nearby M-dwarfs to have the CJ inclination well constrained, as shown in the right panel of Figure 2.

Figure 4 .
Figure 4.A random cumulative distribution for the absolute value of the cosine of the inclination of the cold Jupiter generated assuming that the mutual inclination follows a Rayleigh distribution with σ = 5 and 20 • , and Uniform.