# Planet Hunters TESS. V. A Planetary System Around a Binary Star, Including a Mini-Neptune in the Habitable Zone

Nora Eisner, Samuel Grunblatt, Oscar Barragán, Thea Faridani, Chris Lintott, Suzanne Aigrain, Cole Johnston, Ian Mason, Keivan Stassun, Megan Bedell, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Nora Eisner, Samuel Grunblatt, Oscar Barragán, Thea Faridani, Chris Lintott, et al.

HAL Id: hal-04599976

https://u-paris.hal.science/hal-04599976

Submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# Planet Hunters TESS. V. A Planetary System Around a Binary Star, Including a Mini-Neptune in the Habitable Zone

Nora L. Eisner<sup>1,2,33,34</sup>,<sup>5</sup> Samuel K. Grunblatt<sup>3</sup>,<sup>5</sup> Oscar Barragán<sup>4</sup>,<sup>5</sup> Thea H. Faridani<sup>5</sup>,<sup>6</sup> Chris Lintott<sup>4</sup>,<sup>6</sup> Suzanne Aigrain<sup>4</sup>,<sup>6</sup> Cole Johnston<sup>6</sup>,<sup>7</sup> Ian R. Mason<sup>7</sup>, Keivan G. Stassun<sup>8</sup>,<sup>7</sup> Megan Bedell<sup>1</sup>,<sup>7</sup> Andrew W. Boyle<sup>9</sup>,<sup>7</sup> David R. Ciardi<sup>10</sup>,<sup>7</sup> Catherine A. Clark<sup>11,12</sup>,<sup>13</sup>, Guillaume Hebrard<sup>13,14</sup>, David W. Hogg<sup>1</sup>,<sup>13</sup> Steve B. Howell<sup>15</sup>,<sup>16</sup> Baptiste Klein<sup>4</sup>, Joe Llama<sup>16</sup>,<sup>17</sup> Joshua N. Winn<sup>2</sup>,<sup>17</sup> Lily L. Zhao<sup>1</sup>,<sup>17</sup> Joseph M. Akana Murphy<sup>17,35</sup>,<sup>18</sup> Corey Beard<sup>18,36</sup>,<sup>19</sup> Casey L. Brinkman<sup>19</sup>,<sup>19</sup> Ashley Chontos<sup>2,37</sup>,<sup>19</sup> Pia Cortes-Zuleta<sup>20</sup>,<sup>19</sup> Xavier Delfosse<sup>21</sup>, Steven Giacalone<sup>22</sup>,<sup>19</sup> Emily A. Gilbert<sup>23</sup>,<sup>19</sup> Neda Heidari<sup>20</sup>, Rae Holcomb<sup>18</sup>,<sup>19</sup> Jon M. Jenkins<sup>15</sup>,<sup>19</sup> Flavien Kiefer<sup>24,25</sup>,<sup>19</sup> Jack Lubin<sup>26</sup>,<sup>19</sup> Eder Martioli<sup>27,28</sup>,<sup>19</sup> Alex S. Polanski<sup>29</sup>,<sup>19</sup> Nicholas Saunders<sup>19,38</sup>,<sup>19</sup> Sara Seager<sup>30,31,32</sup>,<sup>19</sup> Avi Shporer<sup>30</sup>,<sup>19</sup> Dakotah Tyler<sup>5</sup>,<sup>19</sup> Judah Van Zandt<sup>5</sup>,<sup>19</sup> Safaa Alhassan<sup>7</sup>, Daval J. Amratlal<sup>7</sup>,<sup>19</sup> Lais I. Antonel<sup>7</sup>, Simon L. S. Bentzen<sup>7</sup>, Milton K. D Bosch<sup>7</sup>,<sup>19</sup> David Bundy<sup>7</sup>, Itayi Chitsiga<sup>7</sup>, Jérôme F. Delaunay<sup>7</sup>, Xavier Doisy<sup>7</sup>, Richard Ferstenou<sup>7</sup>, Mark Fynø<sup>7</sup>, James M. Geary<sup>7</sup>, Gerry Haynaly<sup>7</sup>, Pete Hermes<sup>7</sup>, Marc Huten<sup>7</sup>, Sam Lee<sup>7</sup>, Paul Metcalfe<sup>7</sup>, Garry J. Pennell<sup>7</sup>, Joanna Puszkarska<sup>7</sup>, Thomas Schäfer<sup>7</sup>, Lisa Stiller<sup>7</sup>,<sup>19</sup> Christopher Tanner<sup>7</sup>, Allan Tarr<sup>7</sup>, and Andrew Wilkinson<sup>7</sup>

<sup>1</sup>Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA

<sup>2</sup>Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

<sup>3</sup>Department of Physics and Astronomy, Johns Hopkins University, 3400 N Charles Street, Baltimore, MD 21218, USA

<sup>4</sup>Oxford Astrophysics, Denys Wilkinson Building, University of Oxford, OX1 3RH, Oxford, UK

<sup>5</sup>Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA

<sup>6</sup>Radboud University Nijmegen, Department of Astrophysics, IMAPP, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands

<sup>7</sup>Citizen Scientist, Zooniverse c/o University of Oxford, Keble Road, Oxford OX1 3RH, UK

<sup>8</sup>Vanderbilt University, Department of Physics & Astronomy, 6301 Stevenson Center Lane, Nashville, TN 37235, USA

<sup>9</sup>Department of Astronomy, California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA 91125, USA

<sup>10</sup>NASA Exoplanet Science Institute-Caltech/IPAC, Pasadena, CA 91350 USA

<sup>11</sup>Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109 USA

<sup>12</sup>NASA Exoplanet Science Institute, IPAC, California Institute of Technology, Pasadena, CA 91125 USA

<sup>13</sup>Institut d'Astrophysique de Paris, UMR 7095 CNRS, Université Pierre & Marie Curie, 98 bis boulevard Arago, F-75014 Paris, France

<sup>14</sup>Observatoire de Haute-Provence, Université d'Aix-Marseille & CNRS, F-04870 Saint Michel l'Observatoire, France

<sup>15</sup>NASA Ames Research Center, Moffett Field, CA 94035, USA

<sup>16</sup>Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA

<sup>17</sup>Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA

<sup>18</sup>Department of Physics & Astronomy, The University of California, Irvine, Irvine, CA 92697, USA

<sup>19</sup>Institute for Astronomy, University of Hawaii at Mānoa, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

<sup>20</sup>Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

<sup>21</sup>University Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

<sup>22</sup>Department of Astronomy, University of California Berkeley, Berkeley, CA 94720, USA

<sup>23</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

<sup>24</sup>LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France

<sup>25</sup>American University of Paris, 5, boulevard de La Tour-Maubourg F-75007 Paris, France

<sup>26</sup>Department of Physics & Astronomy, University of California Irvine, Irvine, CA 92697, USA

<sup>27</sup>Laboratório Nacional de Astrofísica, Rua Estados Unidos 154, 37504-364, Itajubá - MG, Brazil

<sup>28</sup>Institut d'Astrophysique de Paris, CNRS, UMR 7095, Sorbonne Université, 98 bis bd Arago, F-75014 Paris, France

<sup>29</sup>Department of Physics and Astronomy, University of Kansas, Lawrence, KS, USA

<sup>30</sup>Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

<sup>31</sup>Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

<sup>32</sup>Department of Aeronautics and Astronautics, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA

*Received 2023 October 20; revised 2024 January 4; accepted 2024 January 9; published 2024 April 30*

## Abstract

We report on the discovery and validation of a transiting long-period mini-Neptune orbiting a bright ($V=9.0$ mag) G dwarf (TOI 4633; $R=1.05 R_⊕$, $M=1.10 M_⊕$). The planet was identified in data from the Transiting Exoplanet Survey Satellite by citizen scientists taking part in the Planet Hunters TESS project. Modelling of the transit events yields an orbital period of 271.9445 ± 0.0040 days and radius of $3.2 \pm 0.20 R_⊕$. The Earth-like orbital period and

<figure><img src="image_2.png" /></figure>

<footnote n="33"><sup>33</sup>Flatiron Research Fellow.</footnote>

<footnote n="34"><sup>34</sup>Henry Norris Russel Fellow.</footnote>

<footnote n="35"><sup>35</sup>NSF Graduate Research Fellow.</footnote>

<footnote n="36"><sup>36</sup>NASA FINESST Fellow.</footnote>

<footnote n="37"><sup>37</sup>Henry Norris Russel Fellow.</footnote>

<footnote n="38"><sup>38</sup>NSF Graduate Research Fellow.</footnote>
# Planet Hunters TESS. V. A Planetary System Around a Binary Star, Including a Mini-Neptune in the Habitable Zone

Nora L. Eisner<sup>1,2,33,34</sup>,<sup>5</sup> Samuel K. Grunblatt<sup>3</sup>,<sup>5</sup> Oscar Barragán<sup>4</sup>,<sup>5</sup> Thea H. Faridani<sup>5</sup>,<sup>6</sup> Chris Lintott<sup>4</sup>,<sup>6</sup> Suzanne Aigrain<sup>4</sup>,<sup>6</sup> Cole Johnston<sup>6</sup>,<sup>7</sup> Ian R. Mason<sup>7</sup>, Keivan G. Stassun<sup>8</sup>,<sup>7</sup> Megan Bedell<sup>1</sup>,<sup>8</sup> Andrew W. Boyle<sup>9</sup>,<sup>8</sup> David R. Ciardi<sup>10</sup>,<sup>9</sup> Catherine A. Clark<sup>11,12</sup>,<sup>13</sup>, Guillaume Hebrard<sup>13,14</sup>, David W. Hogg<sup>1</sup>,<sup>15</sup> Steve B. Howell<sup>15</sup>,<sup>16</sup> Baptiste Klein<sup>4</sup>, Joe Llama<sup>16</sup>, Joshua N. Winn<sup>2</sup>,<sup>16</sup> Lily L. Zhao<sup>1</sup>,<sup>16</sup> Joseph M. Akana Murphy<sup>17,35</sup>,<sup>16</sup> Corey Beard<sup>18,36</sup>,<sup>16</sup> Casey L. Brinkman<sup>19</sup>,<sup>16</sup> Ashley Chontos<sup>2,37</sup>,<sup>16</sup> Pia Cortes-Zuleta<sup>20</sup>,<sup>16</sup> Xavier Delfosse<sup>21</sup>, Steven Giacalone<sup>22</sup>,<sup>16</sup> Emily A. Gilbert<sup>23</sup>,<sup>16</sup> Neda Heidari<sup>20</sup>, Rae Holcomb<sup>18</sup>,<sup>16</sup> Jon M. Jenkins<sup>15</sup>,<sup>16</sup> Flavien Kiefer<sup>24,25</sup>,<sup>16</sup> Jack Lubin<sup>26</sup>,<sup>16</sup> Eder Martioli<sup>27,28</sup>,<sup>16</sup> Alex S. Polanski<sup>29</sup>,<sup>16</sup> Nicholas Saunders<sup>19,38</sup>,<sup>16</sup> Sara Seager<sup>30,31,32</sup>,<sup>16</sup> Avi Shporer<sup>30</sup>,<sup>16</sup> Dakotah Tyler<sup>5</sup>,<sup>16</sup> Judah Van Zandt<sup>5</sup>,<sup>16</sup> Safaa Alhassan<sup>7</sup>, Daval J. Amratlal<sup>7</sup>,<sup>16</sup> Lais I. Antonel<sup>7</sup>, Simon L. S. Bentzen<sup>7</sup>, Milton K. D Bosch<sup>7</sup>,<sup>16</sup> David Bundy<sup>7</sup>, Itayi Chitsiga<sup>7</sup>, Jérôme F. Delaunay<sup>7</sup>, Xavier Doisy<sup>7</sup>, Richard Ferstenou<sup>7</sup>, Mark Fynø<sup>7</sup>, James M. Geary<sup>7</sup>, Gerry Haynaly<sup>7</sup>, Pete Hermes<sup>7</sup>, Marc Huten<sup>7</sup>, Sam Lee<sup>7</sup>, Paul Metcalfe<sup>7</sup>, Garry J. Pennell<sup>7</sup>, Joanna Puszkarska<sup>7</sup>, Thomas Schäfer<sup>7</sup>, Lisa Stiller<sup>7</sup>,<sup>16</sup> Christopher Tanner<sup>7</sup>, Allan Tarr<sup>7</sup>, and Andrew Wilkinson<sup>7</sup>

<sup>1</sup>Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA

<sup>2</sup>Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

<sup>3</sup>Department of Physics and Astronomy, Johns Hopkins University, 3400 N Charles Street, Baltimore, MD 21218, USA

<sup>4</sup>Oxford Astrophysics, Denys Wilkinson Building, University of Oxford, OX1 3RH, Oxford, UK

<sup>5</sup>Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA

<sup>6</sup>Radboud University Nijmegen, Department of Astrophysics, IMAPP, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands

<sup>7</sup>Citizen Scientist, Zooniverse c/o University of Oxford, Keble Road, Oxford OX1 3RH, UK

<sup>8</sup>Vanderbilt University, Department of Physics & Astronomy, 6301 Stevenson Center Lane, Nashville, TN 37235, USA

<sup>9</sup>Department of Astronomy, California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA 91125, USA

<sup>10</sup>NASA Exoplanet Science Institute-Caltech/IPAC, Pasadena, CA 91350 USA

<sup>11</sup>Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109 USA

<sup>12</sup>NASA Exoplanet Science Institute, IPAC, California Institute of Technology, Pasadena, CA 91125 USA

<sup>13</sup>Institut d'Astrophysique de Paris, UMR 7095 CNRS, Université Pierre & Marie Curie, 98 bis boulevard Arago, F-75014 Paris, France

<sup>14</sup>Observatoire de Haute-Provence, Université d'Aix-Marseille & CNRS, F-04870 Saint Michel l'Observatoire, France

<sup>15</sup>NASA Ames Research Center, Moffett Field, CA 94035, USA

<sup>16</sup>Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA

<sup>17</sup>Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA

<sup>18</sup>Department of Physics & Astronomy, The University of California, Irvine, Irvine, CA 92697, USA

<sup>19</sup>Institute for Astronomy, University of Hawaii at Mānoa, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

<sup>20</sup>Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

<sup>21</sup>University Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

<sup>22</sup>Department of Astronomy, University of California Berkeley, Berkeley, CA 94720, USA

<sup>23</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

<sup>24</sup>LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France

<sup>25</sup>American University of Paris, 5, boulevard de La Tour-Maubourg F-75007 Paris, France

<sup>26</sup>Department of Physics & Astronomy, University of California Irvine, Irvine, CA 92697, USA

<sup>27</sup>Laboratório Nacional de Astrofísica, Rua Estados Unidos 154, 37504-364, Itajubá - MG, Brazil

<sup>28</sup>Institut d'Astrophysique de Paris, CNRS, UMR 7095, Sorbonne Université, 98 bis bd Arago, F-75014 Paris, France

<sup>29</sup>Department of Physics and Astronomy, University of Kansas, Lawrence, KS, USA

<sup>30</sup>Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

<sup>31</sup>Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

<sup>32</sup>Department of Aeronautics and Astronautics, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA

*Received 2023 October 20; revised 2024 January 4; accepted 2024 January 9; published 2024 April 30*

## Abstract

We report on the discovery and validation of a transiting long-period mini-Neptune orbiting a bright ($V = 9.0$ mag) G dwarf (TOI 4633; $R = 1.05 R_⊕$, $M = 1.10 M_⊕$). The planet was identified in data from the Transiting Exoplanet Survey Satellite by citizen scientists taking part in the Planet Hunters TESS project. Modelling of the transit events yields an orbital period of 271.9445 ± 0.0040 days and radius of $3.2 ± 0.20 R_⊕$. The Earth-like orbital period and

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<footnote n="33"><sup>33</sup>Flatiron Research Fellow.</footnote>

<footnote n="34"><sup>34</sup>Henry Norris Russel Fellow.</footnote>

<footnote n="35"><sup>35</sup>NSF Graduate Research Fellow.</footnote>

<footnote n="36"><sup>36</sup>NASA FINESST Fellow.</footnote>

<footnote n="37"><sup>37</sup>Henry Norris Russel Fellow.</footnote>

<footnote n="38"><sup>38</sup>NSF Graduate Research Fellow.</footnote>
an incident flux of $1.56_{-0.16}^{+0.20} F_⊕$ places it in the optimistic habitable zone around the star. Doppler spectroscopy of the system allowed us to place an upper mass limit on the transiting planet and revealed a non-transiting planet candidate in the system with a period of $34.15 \pm 0.15$ days. Furthermore, the combination of archival data dating back to 1905 with new high angular resolution imaging revealed a stellar companion orbiting the primary star with an orbital period of around 230 yr and an eccentricity of about 0.9. The long period of the transiting planet, combined with the high eccentricity and close approach of the companion star makes this a valuable system for testing the formation and stability of planets in binary systems.

*Unified Astronomy Thesaurus concepts: Exoplanet systems (484); Exoplanet dynamics (490); Transit photometry (1709); Binary stars (154); Habitable zone (696); Mini Neptunes (1063); Radial velocity (1332)*

# 1. Introduction

The advancement of space based photometric exoplanet missions, such as CoRoT (Auvergne et al. 2009), Kepler (Borucki et al. 2010) and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) has significantly improved our understanding of extrasolar planetary systems, including our understanding of planet occurrence rates and system architectures. However, the detection of planets using the transit method is inherently biased toward short-period planets. This is in part due to the fact that the transit probability of a planet decreases with increased orbital distance from the host star, and in part due to the fact that automated transit detection pipelines typically require two or more transit events in order to reach the signal-to-noise level required for detection and to achieve confidence that the signal is periodic.

As a result, only 9.7% of all confirmed transiting planets have orbital periods longer than 50 days, and 1.6% have orbital periods longer than 200 days. Similarly, only around 2.5% of known planets with a semimajor axis greater than 1 au were detected using the transit method, with over 70% of them having been detected using Radial Velocity (RV) observations (NASA Exoplanet Archive).<fnref n="39" /> While the RV method can yield planet properties such as the orbital period and minimum mass measurements, without the detection of a transit event the planet radius, and therefore the bulk density, cannot be constrained. Furthermore, the detection of a transit event helps constrain the system inclination, thus enabling an absolute (instead of minimum) mass measurement. Similarly, atmospheric characterization via transmission spectroscopy is only possible for transiting planets.

Transiting planets on long orbital periods, in particular, allow for new investigations into the formation, migration, and long-term stability of planetary systems. The comparison between planets with short and long orbital periods, for example, allows us to probe how equilibrium temperatures affect planet formation (e.g., Lopez & Rice 2018; Fernandes et al. 2019).

The long-term stability and evolution of planetary systems can also be affected by stellar binarity (Veras 2016; Hamer & Schlaufman 2019). As shown by Raghavan et al. (2010), $54\% \pm 2\%$ of solar-type stars are single, with the rest existing in pairs or higher order multiple systems. These companion stars can perturb planet orbits resulting in high eccentricity tidal migration which can produce hot Jupiters (Naoz et al. 2012; Vick et al. 2023), truncate protoplanetary disks and shorten disk lifespans (Kraus et al. 2012; Manara et al. 2019; Winter et al. 2020; Zagaria et al. 2022), and limit the formation of terrestrial planets (when the binary separation is less than around 10 au; Quintana et al. 2007). Due to the high fraction of stars that are part of binaries, a thorough understanding of how

binary interactions affect planet formation, migration and long-term stability is important to constrain the underlying planet population in our Galaxy (Moe & Kratter 2021a).

In this paper we present the detection and validation of a transiting mini-Neptune (hereafter TOI 4633 c) orbiting a bright (Vmag = 9.0), nearby ($d = 95$ pc) G Dwarf, that was detected by citizen scientists. The period of 272 days makes this planet the second longest-period confirmed planet identified in the TESS data to date (with the longest being TOI 4600c, Mireles et al. 2023), and only one of five confirmed TESS planets with orbital periods longer than 100 days (Dalba et al. 2022; Heitzmann et al. 2023). Furthermore, the long orbital period and incident flux of $\sim 1.6 F_⊕$ places it in the habitable zone of its host star, making it only the fourth habitable zone planet identified in the TESS data to date, following TOI 700 d (Gilbert et al. 2020; Rodriguez et al. 2020), TOI 700 e (Gilbert et al. 2023), and TOI 715 b (Dransfield et al. 2024).

RV monitoring revealed an additional, non-transiting, planet candidate with a 34 days period (hereafter TOI 4633 b). Furthermore, the combination of newly obtained speckle imaging and archival high-contrast imaging data dating back to 1905 revealed a bound stellar companion with a period of $\sim 230$ yr. As such, we present a bright multi-planet, multi-star system.

The paper is structured as follows. In Section 2 we describe the discovery of TOI 4633 c in the TESS data. In Section 3 we present the spectroscopic and imaging observations and in Section 4 we discuss activity indicators and statistically validate the transit signal. In Section 5 we present the derivation of the parameters of both stars and the orbital properties of the binary system. In Section 6 we discuss the planet parameters of TOI 4633 c and planet candidate TOI 4633 b and their long-term stability within the binary system. Finally, in Section 7, we place TOI 4633 into the context of other long-period multi-planet systems as well as into the context of the population of confirmed planets in binaries.

# 2. Photometry and Discovery of TOI 4633 c

TOI 4633 (TIC 307958020; Stassun et al. 2019) is located at high ecliptic latitude (near the ecliptic pole) and was observed by TESS nearly continually during years 2, 4 and 5 of the mission (Sectors 14–26, 40, 41, 47–53, 55, 56, 58, 59). Following the identification of the first two transit events, we proposed for the target to be observed at the shortest TESS cadence. Therefore, from Sector 49 and onward, the data have been obtained at a 20 s cadence (proposal ID: DDT054; PI: Eisner). Prior to this sector all observations were obtained at a 2 minutes cadence.

The light curve of TOI 4633 exhibits three transit events located in Sectors 20, 40 and 50. The former two transit events were identified by citizen scientists taking part in the Planet Hunters TESS citizen science project (Eisner et al. 2021). The project, which is hosted by the Zooniverse platform (Lintott et al. 2008, 2011), engages over 40,000 citizen scientists in the task of

<footnote n="39">³⁹ https://exoplanetarchive.ipac.caltech.edu</footnote>
an incident flux of $1.56_{-0.16}^{+0.20} F_⊕$ places it in the optimistic habitable zone around the star. Doppler spectroscopy of the system allowed us to place an upper mass limit on the transiting planet and revealed a non-transiting planet candidate in the system with a period of $34.15 \pm 0.15$ days. Furthermore, the combination of archival data dating back to 1905 with new high angular resolution imaging revealed a stellar companion orbiting the primary star with an orbital period of around 230 yr and an eccentricity of about 0.9. The long period of the transiting planet, combined with the high eccentricity and close approach of the companion star makes this a valuable system for testing the formation and stability of planets in binary systems.

*Unified Astronomy Thesaurus concepts: Exoplanet systems (484); Exoplanet dynamics (490); Transit photometry (1709); Binary stars (154); Habitable zone (696); Mini Neptunes (1063); Radial velocity (1332)*

# 1. Introduction

The advancement of space based photometric exoplanet missions, such as CoRoT (Auvergne et al. 2009), Kepler (Borucki et al. 2010) and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) has significantly improved our understanding of extrasolar planetary systems, including our understanding of planet occurrence rates and system architectures. However, the detection of planets using the transit method is inherently biased toward short-period planets. This is in part due to the fact that the transit probability of a planet decreases with increased orbital distance from the host star, and in part due to the fact that automated transit detection pipelines typically require two or more transit events in order to reach the signal-to-noise level required for detection and to achieve confidence that the signal is periodic.

As a result, only 9.7% of all confirmed transiting planets have orbital periods longer than 50 days, and 1.6% have orbital periods longer than 200 days. Similarly, only around 2.5% of known planets with a semimajor axis greater than 1 au were detected using the transit method, with over 70% of them having been detected using Radial Velocity (RV) observations (NASA Exoplanet Archive).<fnref n="39" /> While the RV method can yield planet properties such as the orbital period and minimum mass measurements, without the detection of a transit event the planet radius, and therefore the bulk density, cannot be constrained. Furthermore, the detection of a transit event helps constrain the system inclination, thus enabling an absolute (instead of minimum) mass measurement. Similarly, atmospheric characterization via transmission spectroscopy is only possible for transiting planets.

Transiting planets on long orbital periods, in particular, allow for new investigations into the formation, migration, and long-term stability of planetary systems. The comparison between planets with short and long orbital periods, for example, allows us to probe how equilibrium temperatures affect planet formation (e.g., Lopez & Rice 2018; Fernandes et al. 2019).

The long-term stability and evolution of planetary systems can also be affected by stellar binarity (Veras 2016; Hamer & Schlaufman 2019). As shown by Raghavan et al. (2010), 54% ± 2% of solar-type stars are single, with the rest existing in pairs or higher order multiple systems. These companion stars can perturb planet orbits resulting in high eccentricity tidal migration which can produce hot Jupiters (Naoz et al. 2012; Vick et al. 2023), truncate protoplanetary disks and shorten disk lifespans (Kraus et al. 2012; Manara et al. 2019; Winter et al. 2020; Zagaria et al. 2022), and limit the formation of terrestrial planets (when the binary separation is less than around 10 au; Quintana et al. 2007). Due to the high fraction of stars that are part of binaries, a thorough understanding of how

binary interactions affect planet formation, migration and long-term stability is important to constrain the underlying planet population in our Galaxy (Moe & Kratter 2021a).

In this paper we present the detection and validation of a transiting mini-Neptune (hereafter TOI 4633 c) orbiting a bright (Vmag = 9.0), nearby ($d = 95$ pc) G Dwarf, that was detected by citizen scientists. The period of 272 days makes this planet the second longest-period confirmed planet identified in the TESS data to date (with the longest being TOI 4600c, Mireles et al. 2023), and only one of five confirmed TESS planets with orbital periods longer than 100 days (Dalba et al. 2022; Heitzmann et al. 2023). Furthermore, the long orbital period and incident flux of ~1.6 $F_⊕$ places it in the habitable zone of its host star, making it only the fourth habitable zone planet identified in the TESS data to date, following TOI 700 d (Gilbert et al. 2020; Rodriguez et al. 2020), TOI 700 e (Gilbert et al. 2023), and TOI 715 b (Dransfield et al. 2024).

RV monitoring revealed an additional, non-transiting, planet candidate with a 34 days period (hereafter TOI 4633 b). Furthermore, the combination of newly obtained speckle imaging and archival high-contrast imaging data dating back to 1905 revealed a bound stellar companion with a period of ~230 yr. As such, we present a bright multi-planet, multi-star system.

The paper is structured as follows. In Section 2 we describe the discovery of TOI 4633 c in the TESS data. In Section 3 we present the spectroscopic and imaging observations and in Section 4 we discuss activity indicators and statistically validate the transit signal. In Section 5 we present the derivation of the parameters of both stars and the orbital properties of the binary system. In Section 6 we discuss the planet parameters of TOI 4633 c and planet candidate TOI 4633 b and their long-term stability within the binary system. Finally, in Section 7, we place TOI 4633 into the context of other long-period multi-planet systems as well as into the context of the population of confirmed planets in binaries.

# 2. Photometry and Discovery of TOI 4633 c

TOI 4633 (TIC 307958020; Stassun et al. 2019) is located at high ecliptic latitude (near the ecliptic pole) and was observed by TESS nearly continually during years 2, 4 and 5 of the mission (Sectors 14–26, 40, 41, 47–53, 55, 56, 58, 59). Following the identification of the first two transit events, we proposed for the target to be observed at the shortest TESS cadence. Therefore, from Sector 49 and onward, the data have been obtained at a 20 s cadence (proposal ID: DDT054; PI: Eisner). Prior to this sector all observations were obtained at a 2 minutes cadence.

The light curve of TOI 4633 exhibits three transit events located in Sectors 20, 40 and 50. The former two transit events were identified by citizen scientists taking part in the Planet Hunters TESS citizen science project (Eisner et al. 2021). The project, which is hosted by the Zooniverse platform (Lintott et al. 2008, 2011), engages over 40,000 citizen scientists in the task of

<footnote n="39">³⁹ https://exoplanetarchive.ipac.caltech.edu</footnote>
visually inspecting TESS data in the search for transit events. At any given time, Planet Hunters TESS only ever displays the data from a single TESS sector. As such, the former two transit events were independently identified by 15 citizen scientists who were randomly presented the light curve of TOI 4633. The data consist of the Presearch Data Conditioning (PDC) 2 minutes cadence observations, which are produced by the Science Processing Operations Center (SPOC; Smith et al. 2012; Stumpe et al. 2012, 2012, 2014; Jenkins et al. 2016) pipeline. Individual measurements that are flagged by the SPOC pipeline as being affected by various instrumental anomalies are not shown on Planet Hunters TESS. Unfortunately, the data around the time of the third transit event (in Sector 50) was identified as being affected by scattered light and as such this transit was not seen or identified by citizen scientists. For more details regarding the PHT pipeline and the identification of transit-like signals we refer the reader to Eisner et al. (2021).

Once the transit events were identified, we analyzed both the 2 minutes cadence and 20 s cadence PDC SPOC light curves. As mentioned above, the 2 minutes and 20 s cadence data around the time of the third transit event (BTJD~2680.6 days) were affected by scattered light. In order to recover the third transit, we performed a tailored correction of the Sector 50 light curve using cotrending basis vectors (CBVs). We followed a similar approach as the one presented by Barragán et al. (2022b), which makes use of the *lightkurve* package (Lightkurve Collaboration et al. 2018). In brief, we first created a light curve from the target pixel file (TPF) using the nominal TESS aperture and did not remove data that had been flagged as “bad” by the SPOC pipeline (flag: bit 13, value 4096). The CBVs provided with the TPF were then used to correct the light curve using the built-in correction function (CBVCorrec-tor<fnref n="40" />) in *lightkurve* and allowing for interpolation. This generates a light curve where the large scale trends are removed, including the trends around the times of the transit events. Finally, we performed a crowding correction to account for extra flux from nearby stars that may be present in the Simple Aperture Photometry (SAP) mask. To do this we used the nominal crowding values given in the TPF to account and correct for the light curve contamination.<fnref n="41" /> We note that this method did not correct for the light contribution from the close, bound companion star (see Section 5). As a test, we also used this method to extract and detrend the light curves around the time of the other two transit events. We found there to be no significant difference between these and the SPOC light curves. As such, for these two former transits we use the SPOC data.

The absence of further transit events in the TESS data allows us to confirm that the transiting object (TOI 4633 c) has an orbital period of ~271.94 days. We confirm that all shorter aliases of this orbital period would result in at least one observable transit in the available TESS data. A subset of the 2 minutes cadence light curve is shown in Figure 1 and the three individual transits are shown in Figure 2.

In addition to the TESS data, there is archival data from the All Sky Automated Survey for SuperNovae (ASAS-SN), consisting of 3482 observations obtained between 2012 March 28 and 2023 November 10, with 2481 observations in the *g* filter, and 1001 in the *V* filter (Shappee et al. 2014; Kochanek

et al. 2017). However, we note that the *g* filter data show a large change in the root mean square scatter before and after 2459330 HJD of around 105 and 170 mJy, respectively. We, therefore, do not use the *g* filter data for any further analysis. We find no evidence of any further transit signals in the ASAS-SN data

## 2.1. Excluding False Positive Scenarios

In order to help rule out instrumental and astrophysical false positive scenarios including background eclipsing binaries, systematic effects, and background events such as asteroids passing through the field of view we performed a number of standard diagnostic tests using the TESS data. We used the publicly available Lightcurve Analysis Tool for Transiting Exoplanets (LATTE; Eisner et al. 2020a) for this analysis.<fnref n="42" /> For a full description of the diagnostic tests we refer the reader to Eisner et al. (2020a); however, in brief, the tests allowed us to ensure that:

1. the transit events do not coincide with the times of the periodic momentum dumps that introduce spurious signals into the data.

2. the *x* and *y* centroid positions are smoothly varying with time in the vicinity of the transit events and thus the transit events are unlikely to be caused by systematic effects or by background eclipsing binaries.

3. the light curves of the five nearest two-minute cadence TESS stars do not show similar signals at the same time (projected distances between these stars and TOI 4633 range from 2.7 to 28.7).

4. the signal is on target by investigating the light curve extracted for each pixel surrounding the target in order to ensure that the transit events are not caused by spurious background signals on other pixels.

5. there are no spurious signals, such as sudden jumps or strong variations, in the background flux at the same time as the event.

6. the transit shapes and depths when extracted with different aperture sizes are consistent.

7. the signal is on target by comparing the average in-transit and average out-of-transit flux, as well as the difference between them, which indicates the location of the change in flux (i.e., the location of the transit event).

We used these tests to show that the transit signals are unlikely to be systematic or astrophysical false positives. We uploaded TOI 4633 c to the Exoplanet Follow-up Observing Program for TESS (ExoFOP-TESS) site on 2020 May 27 as a community TESS Object of Interest (cTOI). The planet was later promoted to a priority 1 (1 = highest priority, 5 = lowest priority) TOI candidate (TOI 4633.01) on 2022 December 14.

While the SPOC did not recover the true period of TOI 4633 c, the first two transits were detected as a duo-transit in the two multi-sector searches conducted over sectors 14–23 and 14–26. Transit searches conducted with a noise-compensating matched filter over sectors 14–41, 14–50, 14–55, and 14–59 (Jenkins 2002; Jenkins et al. 2010, 2020) identified the first two transits at twice the actual orbital period. A limb-darkened transit model was fit to the transit signal in each case (Li et al. 2019) and a suite of diagnostic tests were performed to assess the nature of the signal (Twicken et al. 2018). The signal

<footnote n="40">⁴⁰ See *lightkurve* documentation https://docs.lightkurve.org/tutorials/2-creating-light-curves/2-3-how-to-use-cbvcorrector.html</footnote>

<footnote n="41">⁴¹ See https://heasarc.gsfc.nasa.gov/docs/tess/UnderstandingCrowding.html for more details on TESS crowding correction.</footnote>

<footnote n="42">⁴² http://latte-online flatironinstitute.org/</footnote>
visually inspecting TESS data in the search for transit events. At any given time, Planet Hunters TESS only ever displays the data from a single TESS sector. As such, the former two transit events were independently identified by 15 citizen scientists who were randomly presented the light curve of TOI 4633. The data consist of the Presearch Data Conditioning (PDC) 2 minutes cadence observations, which are produced by the Science Processing Operations Center (SPOC; Smith et al. 2012; Stumpe et al. 2012, 2012, 2014; Jenkins et al. 2016) pipeline. Individual measurements that are flagged by the SPOC pipeline as being affected by various instrumental anomalies are not shown on Planet Hunters TESS. Unfortunately, the data around the time of the third transit event (in Sector 50) was identified as being affected by scattered light and as such this transit was not seen or identified by citizen scientists. For more details regarding the PHT pipeline and the identification of transit-like signals we refer the reader to Eisner et al. (2021).

Once the transit events were identified, we analyzed both the 2 minutes cadence and 20 s cadence PDC SPOC light curves. As mentioned above, the 2 minutes and 20 s cadence data around the time of the third transit event (BTJD~2680.6 days) were affected by scattered light. In order to recover the third transit, we performed a tailored correction of the Sector 50 light curve using cotrending basis vectors (CBVs). We followed a similar approach as the one presented by Barragán et al. (2022b), which makes use of the *lightkurve* package (Lightkurve Collaboration et al. 2018). In brief, we first created a light curve from the target pixel file (TPF) using the nominal TESS aperture and did not remove data that had been flagged as “bad” by the SPOC pipeline (flag: bit 13, value 4096). The CBVs provided with the TPF were then used to correct the light curve using the built-in correction function (CBVCorrec-tor<fnref n="40" />) in *lightkurve* and allowing for interpolation. This generates a light curve where the large scale trends are removed, including the trends around the times of the transit events. Finally, we performed a crowding correction to account for extra flux from nearby stars that may be present in the Simple Aperture Photometry (SAP) mask. To do this we used the nominal crowding values given in the TPF to account and correct for the light curve contamination.<fnref n="41" /> We note that this method did not correct for the light contribution from the close, bound companion star (see Section 5). As a test, we also used this method to extract and detrend the light curves around the time of the other two transit events. We found there to be no significant difference between these and the SPOC light curves. As such, for these two former transits we use the SPOC data.

The absence of further transit events in the TESS data allows us to confirm that the transiting object (TOI 4633 c) has an orbital period of ~271.94 days. We confirm that all shorter aliases of this orbital period would result in at least one observable transit in the available TESS data. A subset of the 2 minutes cadence light curve is shown in Figure 1 and the three individual transits are shown in Figure 2.

In addition to the TESS data, there is archival data from the All Sky Automated Survey for SuperNovae (ASAS-SN), consisting of 3482 observations obtained between 2012 March 28 and 2023 November 10, with 2481 observations in the g filter, and 1001 in the V filter (Shappee et al. 2014; Kochanek

et al. 2017). However, we note that the g filter data show a large change in the root mean square scatter before and after 2459330 HJD of around 105 and 170 mJy, respectively. We, therefore, do not use the g filter data for any further analysis. We find no evidence of any further transit signals in the ASAS-SN data

## 2.1. Excluding False Positive Scenarios

In order to help rule out instrumental and astrophysical false positive scenarios including background eclipsing binaries, systematic effects, and background events such as asteroids passing through the field of view we performed a number of standard diagnostic tests using the TESS data. We used the publicly available Lightcurve Analysis Tool for Transiting Exoplanets (LATTE; Eisner et al. 2020a) for this analysis.<fnref n="42" /> For a full description of the diagnostic tests we refer the reader to Eisner et al. (2020a); however, in brief, the tests allowed us to ensure that:

1. the transit events do not coincide with the times of the periodic momentum dumps that introduce spurious signals into the data.

2. the x and y centroid positions are smoothly varying with time in the vicinity of the transit events and thus the transit events are unlikely to be caused by systematic effects or by background eclipsing binaries.

3. the light curves of the five nearest two-minute cadence TESS stars do not show similar signals at the same time (projected distances between these stars and TOI 4633 range from 2.7 to 28.7).

4. the signal is on target by investigating the light curve extracted for each pixel surrounding the target in order to ensure that the transit events are not caused by spurious background signals on other pixels.

5. there are no spurious signals, such as sudden jumps or strong variations, in the background flux at the same time as the event.

6. the transit shapes and depths when extracted with different aperture sizes are consistent.

7. the signal is on target by comparing the average in-transit and average out-of-transit flux, as well as the difference between them, which indicates the location of the change in flux (i.e., the location of the transit event).

We used these tests to show that the transit signals are unlikely to be systematic or astrophysical false positives. We uploaded TOI 4633 c to the Exoplanet Follow-up Observing Program for TESS (ExoFOP-TESS) site on 2020 May 27 as a community TESS Object of Interest (cTOI). The planet was later promoted to a priority 1 (1 = highest priority, 5 = lowest priority) TOI candidate (TOI 4633.01) on 2022 December 14.

While the SPOC did not recover the true period of TOI 4633 c, the first two transits were detected as a duo-transit in the two multi-sector searches conducted over sectors 14–23 and 14–26. Transit searches conducted with a noise-compensating matched filter over sectors 14–41, 14–50, 14–55, and 14–59 (Jenkins 2002; Jenkins et al. 2010, 2020) identified the first two transits at twice the actual orbital period. A limb-darkened transit model was fit to the transit signal in each case (Li et al. 2019) and a suite of diagnostic tests were performed to assess the nature of the signal (Twicken et al. 2018). The signal

<footnote n="40">⁴⁰ See *lightkurve* documentation https://docs.lightkurve.org/tutorials/2-creating-light-curves/2-3-how-to-use-cbvcorrector.html</footnote>

<footnote n="41">⁴¹ See https://heasarc.gsfc.nasa.gov/docs/tess/UnderstandingCrowding.html for more details on TESS crowding correction.</footnote>

<footnote n="42">⁴² http://latte-online flatironinstitute.org/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Normalized photometric data from TESS sectors 20–22 to illustrate the stellar variability of the light curve. The first of the three transit events observed by TESS to date is highlighted by the vertical gray column. The normalized 2 minutes cadence flux binned to 10 minutes and unbinned are shown in maroon and gray, respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The three transits of TOI 4633 c. Nominal TESS observations are shown in light gray while solid circles represent 60 minutes binned data.</figcaption></figure>

passed all the diagnostic tests, including the difference image centroiding test, which localized the source of the transit signal to within 2''3 ± 4''1 of the presumed host star. The TESS Science Office reviewed the Data Validation reports and issued an alert for TOI 4633c on 2021 November 19 (Guerrero et al. 2021).

### 3. Follow-up Observations

In this section we outline the spectroscopic data, as well as the imaging data that revealed a close, bound stellar

companion. The analysis of these data is discussed in Sections 5 and 6.

#### 3.1. Spectroscopic RV Monitoring

We acquired high-resolution ($R \sim 120,000$) spectra with the High Resolution Echelle Spectrometer (HIRES) mounted on the 10 m Keck telescope on Maunakea, Hawaii (Vogt & Donald Penrod 1988, Program numbers: 2022A_N93, 2022B_N025, and 2023A_N085). The instrument has a wavelength coverage between ~350 and ~620 nm. The instrument passes the star's light through a heated iodine cell, allowing for precise wavelength calibration when determining relative RVs (Howard et al. 2010).

We obtained 54 spectra using the iodine cell between 2022 February 21 and 2023 July 22 (mean per pixel S/N ~167 at 550 nm). Radial velocity measurements were determined following the approach of Vogt et al. (1994) and then fit using the publicly available software package RadVel via comparison to a S/N = 260 template spectrum of G5 standard star HD 162232 taken on 2008 June 20 (Fulton et al. 2018). The RV analysis was also carried out using a template spectrum of TOI 4633 taken on 2021 November 24 (per pixel S/N = 46), which yielded equivalent results within the uncertainties. The analysis of the HIRES spectra showed no evidence of a companion star, allowing us to place an upper limit of $\Delta RV \le 10 \text{ km s}^{-1}$ on the relative motion of the two stars known to be in this system (Kolbl et al. 2015). Furthermore, we note that we do not see any evidence for multiple stars in our radial velocity template or in the analysis of the measured radial velocity values alone. Due to the lack of the obvious signature of both stars in the spectra, we are unable to use a two-star template spectrum to extract the RVs.

Additional data were obtained with the SOPHIE high-resolution fiber-fed echelle spectrograph mounted on the 1.93 m Observatoire de Haute-Provence (OHP) telescope (Perruchot et al. 2008). We used the high-resolution mode, which delivers a spectral resolution of $R \sim 75,000$ across the wavelength range of 387–694 nm. We obtained 20 observations between 2021 November 15 and 2023 March 5 (mean per pixel S/N ~33 at 555 nm). Standard stars that were observed at the same epochs using the same SOPHIE mode did not show significant instrumental drifts. The spectra were reduced using the standard SOPHIE RV reduction pipeline (Bouchy et al. 2009), including
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Normalized photometric data from TESS sectors 20–22 to illustrate the stellar variability of the light curve. The first of the three transit events observed by TESS to date is highlighted by the vertical gray column. The normalized 2 minutes cadence flux binned to 10 minutes and unbinned are shown in maroon and gray, respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** The three transits of TOI 4633 c. Nominal TESS observations are shown in light gray while solid circles represent 60 minutes binned data.</figcaption></figure>

passed all the diagnostic tests, including the difference image centroiding test, which localized the source of the transit signal to within $2''3 \pm 4''1$ of the presumed host star. The TESS Science Office reviewed the Data Validation reports and issued an alert for TOI 4633c on 2021 November 19 (Guerrero et al. 2021).

### 3. Follow-up Observations

In this section we outline the spectroscopic data, as well as the imaging data that revealed a close, bound stellar

companion. The analysis of these data is discussed in Sections 5 and 6.

#### 3.1. Spectroscopic RV Monitoring

We acquired high-resolution ($R \sim 120,000$) spectra with the High Resolution Echelle Spectrometer (HIRES) mounted on the 10 m Keck telescope on Maunakea, Hawaii (Vogt & Donald Penrod 1988, Program numbers: 2022A_N93, 2022B_N025, and 2023A_N085). The instrument has a wavelength coverage between $\sim 350$ and $\sim 620$ nm. The instrument passes the star’s light through a heated iodine cell, allowing for precise wavelength calibration when determining relative RVs (Howard et al. 2010).

We obtained 54 spectra using the iodine cell between 2022 February 21 and 2023 July 22 (mean per pixel S/N $\sim 167$ at 550 nm). Radial velocity measurements were determined following the approach of Vogt et al. (1994) and then fit using the publicly available software package RadVel via comparison to a S/N = 260 template spectrum of G5 standard star HD 162232 taken on 2008 June 20 (Fulton et al. 2018). The RV analysis was also carried out using a template spectrum of TOI 4633 taken on 2021 November 24 (per pixel S/N = 46), which yielded equivalent results within the uncertainties. The analysis of the HIRES spectra showed no evidence of a companion star, allowing us to place an upper limit of $\Delta RV \le 10 \text{ km s}^{-1}$ on the relative motion of the two stars known to be in this system (Kolbl et al. 2015). Furthermore, we note that we do not see any evidence for multiple stars in our radial velocity template or in the analysis of the measured radial velocity values alone. Due to the lack of the obvious signature of both stars in the spectra, we are unable to use a two-star template spectrum to extract the RVs.

Additional data were obtained with the SOPHIE high-resolution fiber-fed echelle spectrograph mounted on the 1.93 m Observatoire de Haute-Provence (OHP) telescope (Perruchot et al. 2008). We used the high-resolution mode, which delivers a spectral resolution of $R \sim 75,000$ across the wavelength range of 387–694 nm. We obtained 20 observations between 2021 November 15 and 2023 March 5 (mean per pixel S/N $\sim 33$ at 555 nm). Standard stars that were observed at the same epochs using the same SOPHIE mode did not show significant instrumental drifts. The spectra were reduced using the standard SOPHIE RV reduction pipeline (Bouchy et al. 2009), including
<table><caption>**Table 1**

Spectroscopic Data Obtained with the OHP/SOPHIE (Mean per Pixel S/N ~33 at 555 nm) and Keck/HIRES (mean per Pixel S/N ~167 at 550 nm)</caption><thead><tr><th>Time<br>(BJD-2457000)</th><th>RV<br>(m s<sup>-1</sup>)</th><th>RV error<br>(m s<sup>-1</sup>)</th><th>SNR</th><th>Instrument</th><th>Time<br>(BJD-2457000)</th><th>RV<br>(m s<sup>-1</sup>)</th><th>RV error<br>(m s<sup>-1</sup>)</th><th>SNR</th><th>Instrument</th></tr></thead><tbody><tr><td>2534.3461</td><td>-47.40</td><td>7</td><td>34</td><td>OHP</td><td>2780.9910</td><td>14.50</td><td>1.60</td><td>216</td><td>HIRES</td></tr><tr><td>2601.7097</td><td>-6.40</td><td>5.00</td><td>32</td><td>OHP</td><td>2782.4014</td><td>9.60</td><td>5.00</td><td>31</td><td>OHP</td></tr><tr><td>2606.7163</td><td>9.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2785.9500</td><td>29.50</td><td>1.81</td><td>215</td><td>HIRES</td></tr><tr><td>2622.7200</td><td>33.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2786.9510</td><td>28.40</td><td>1.76</td><td>215</td><td>HIRES</td></tr><tr><td>2623.6940</td><td>10.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2789.8280</td><td>26.75</td><td>1.68</td><td>216</td><td>HIRES</td></tr><tr><td>2631.6864</td><td>-9.40</td><td>5.00</td><td>31</td><td>OHP</td><td>2792.8320</td><td>13.20</td><td>1.50</td><td>214</td><td>HIRES</td></tr><tr><td>2632.0840</td><td>-18.11</td><td>2.67</td><td>215</td><td>HIRES</td><td>2799.9740</td><td>-0.47</td><td>3.26</td><td>75</td><td>HIRES</td></tr><tr><td>2648.6661</td><td>21.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2806.8260</td><td>-15.09</td><td>3.04</td><td>75</td><td>HIRES</td></tr><tr><td>2657.1140</td><td>22.26</td><td>2.23</td><td>214</td><td>HIRES</td><td>2809.8690</td><td>-12.29</td><td>2.97</td><td>75</td><td>HIRES</td></tr><tr><td>2661.0570</td><td>1.01</td><td>2.49</td><td>214</td><td>HIRES</td><td>2812.7900</td><td>-13.26</td><td>2.76</td><td>76</td><td>HIRES</td></tr><tr><td>2662.6285</td><td>-8.40</td><td>3.00</td><td>55</td><td>OHP</td><td>2822.8680</td><td>38.68</td><td>2.66</td><td>75</td><td>HIRES</td></tr><tr><td>2672.0130</td><td>-16.19</td><td>2.20</td><td>216</td><td>HIRES</td><td>2824.4013</td><td>12.60</td><td>5.00</td><td>31</td><td>OHP</td></tr><tr><td>2681.0310</td><td>48.24</td><td>2.47</td><td>174</td><td>HIRES</td><td>2826.7810</td><td>14.72</td><td>2.53</td><td>76</td><td>HIRES</td></tr><tr><td>2683.6198</td><td>14.60</td><td>5.00</td><td>30</td><td>OHP</td><td>2828.7990</td><td>16.66</td><td>3.02</td><td>76</td><td>HIRES</td></tr><tr><td>2687.6359</td><td>8.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2829.7400</td><td>4.81</td><td>3.26</td><td>76</td><td>HIRES</td></tr><tr><td>2690.0110</td><td>9.61</td><td>2.11</td><td>216</td><td>HIRES</td><td>2831.8460</td><td>5.52</td><td>2.79</td><td>73</td><td>HIRES</td></tr><tr><td>2695.0180</td><td>-11.66</td><td>2.05</td><td>214</td><td>HIRES</td><td>2833.7860</td><td>-4.98</td><td>2.74</td><td>75</td><td>HIRES</td></tr><tr><td>2700.9730</td><td>-44.11</td><td>2.09</td><td>213</td><td>HIRES</td><td>2834.7620</td><td>-10.65</td><td>2.87</td><td>76</td><td>HIRES</td></tr><tr><td>2710.0270</td><td>-11.04</td><td>1.76</td><td>215</td><td>HIRES</td><td>2835.7440</td><td>-8.97</td><td>2.69</td><td>77</td><td>HIRES</td></tr><tr><td>2712.0200</td><td>3.27</td><td>2.27</td><td>163</td><td>HIRES</td><td>2838.8030</td><td>-8.72</td><td>2.87</td><td>75</td><td>HIRES</td></tr><tr><td>2712.9410</td><td>18.14</td><td>1.87</td><td>215</td><td>HIRES</td><td>2840.7390</td><td>-10.70</td><td>3.22</td><td>76</td><td>HIRES</td></tr><tr><td>2713.5387</td><td>6.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2858.7850</td><td>37.97</td><td>3.17</td><td>74</td><td>HIRES</td></tr><tr><td>2715.9930</td><td>-2.48</td><td>1.79</td><td>214</td><td>HIRES</td><td>2879.2694</td><td>-27.40</td><td>5.00</td><td>31</td><td>OHP</td></tr><tr><td>2725.5476</td><td>17.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2892.3184</td><td>-9.40</td><td>5.00</td><td>33</td><td>OHP</td></tr><tr><td>2730.0720</td><td>5.38</td><td>1.82</td><td>214</td><td>HIRES</td><td>3009.6665</td><td>-5.40</td><td>6.00</td><td>31</td><td>OHP</td></tr><tr><td>2736.5730</td><td>-32.40</td><td>5.00</td><td>31</td><td>OHP</td><td>3043.0680</td><td>-20.97</td><td>2.65</td><td>143</td><td>HIRES</td></tr><tr><td>2738.9790</td><td>-17.85</td><td>1.54</td><td>215</td><td>HIRES</td><td>3045.0940</td><td>-38.01</td><td>2.18</td><td>216</td><td>HIRES</td></tr><tr><td>2741.9610</td><td>-19.38</td><td>1.51</td><td>215</td><td>HIRES</td><td>3047.1150</td><td>-36.49</td><td>2.07</td><td>202</td><td>HIRES</td></tr><tr><td>2748.8440</td><td>-2.73</td><td>1.63</td><td>204</td><td>HIRES</td><td>3068.0640</td><td>-7.83</td><td>2.08</td><td>153</td><td>HIRES</td></tr><tr><td>2749.8980</td><td>-0.83</td><td>1.48</td><td>216</td><td>HIRES</td><td>3070.9230</td><td>-31.26</td><td>2.01</td><td>214</td><td>HIRES</td></tr><tr><td>2750.4574</td><td>12.60</td><td>4.00</td><td>42</td><td>OHP</td><td>3089.8800</td><td>-9.88</td><td>1.75</td><td>216</td><td>HIRES</td></tr><tr><td>2768.9970</td><td>1.38</td><td>1.57</td><td>215</td><td>HIRES</td><td>3101.0680</td><td>15.11</td><td>1.51</td><td>210</td><td>HIRES</td></tr><tr><td>2769.9960</td><td>-2.72</td><td>1.89</td><td>214</td><td>HIRES</td><td>3108.0060</td><td>-9.65</td><td>1.53</td><td>214</td><td>HIRES</td></tr><tr><td>2773.4748</td><td>-11.40</td><td>6.00</td><td>31</td><td>OHP</td><td>3122.0480</td><td>5.74</td><td>1.63</td><td>171</td><td>HIRES</td></tr><tr><td>2775.9770</td><td>-6.20</td><td>1.58</td><td>216</td><td>HIRES</td><td>3132.8510</td><td>6.65</td><td>1.51</td><td>216</td><td>HIRES</td></tr><tr><td>2776.8030</td><td>-11.86</td><td>1.58</td><td>210</td><td>HIRES</td><td>3138.9630</td><td>5.20</td><td>1.86</td><td>196</td><td>HIRES</td></tr><tr><td>2779.9810</td><td>3.22</td><td>1.69</td><td>214</td><td>HIRES</td><td>3148.0280</td><td>1.62</td><td>1.68</td><td>206</td><td>HIRES</td></tr></tbody></table>

Note. The reported RV values were systematically shifted to be centered around 0, as determined using a Keplerian fit. The low resolution, iodine-free reconnaissance spectra obtained with Keck/HIRES are not included.

CCD charge transfer inefficiency correction (Bouchy et al. 2013). Following the method described, e.g., in Pollacco et al. (2008) and Hebrard et al. (2008), we estimated and corrected for the sky background contamination (mainly due to the Moon) using the second SOPHIE fiber aperture, which is targeted 2' away from the first one pointing toward the star. We estimated that only one of the 20 exposures was significantly polluted by sky background. We note that the analysis of the SOPHIE spectra showed no evidence of a companion star. The extracted HIRES and SOPHIE RV observations are listed in Table 1 and their implications for the system's architecture discussed in Section 6.

The RV data allowed us to place an upper mass limit on planet c and revealed an additional non-transiting planet candidate (with a period of ~34 days; planet b) in the system, as discussed in Section 6.

## 3.2. Imaging Observations and Companion Star.

Archival data listed in the Washington Double Star Catalog (WDS) show that TOI 4633 has a bound companion star.

Hereafter, the stars will be referred to as star A and star B for the more and less massive star, respectively. The 10 archival observations were obtained between 1905 February and 2011 June and their position angles (PA) and angular separations (ρ) are listed in Table 2.

In addition to archival data, we obtained high-contrast imaging observations of TOI 4633 to further constrain the orbit of the two stars. High-contrast imaging was performed using the NESSI high-resolution speckle imaging instrument mounted on the 3.5 m WIYN telescope located at Kitt Peak National Observatory (Scott et al. 2018) on 2022 April 21; the PHARO adaptive optics instrument mounted on 5.1 m Hale telescope at Palomar Observatory (Hayward et al. 2001) on 2021 November 11; the 'Alopeke speckle instrument mounted on the 8.1 m Gemini North telescope on Maunakea (Howell et al. 2011; Matson et al. 2019) on 2022 May 9; and the NIRC2 adaptive optics instrument mounted on the 10 m Keck II (Wizinowich et al. 2000) located on Maunakea, Hawaii, on 2023 April 26. While the former two observations were unable to resolve the two stars, the 'Alopeke and NIRC2 observations
<table><caption>**Table 1**

Spectroscopic Data Obtained with the OHP/SOPHIE (Mean per Pixel S/N ~33 at 555 nm) and Keck/HIRES (mean per Pixel S/N ~167 at 550 nm)</caption><thead><tr><th>Time<br>(BJD-2457000)</th><th>RV<br>(m s<sup>-1</sup>)</th><th>RV error<br>(m s<sup>-1</sup>)</th><th>SNR</th><th>Instrument</th><th>Time<br>(BJD-2457000)</th><th>RV<br>(m s<sup>-1</sup>)</th><th>RV error<br>(m s<sup>-1</sup>)</th><th>SNR</th><th>Instrument</th></tr></thead><tbody><tr><td>2534.3461</td><td>-47.40</td><td>7</td><td>34</td><td>OHP</td><td>2780.9910</td><td>14.50</td><td>1.60</td><td>216</td><td>HIRES</td></tr><tr><td>2601.7097</td><td>-6.40</td><td>5.00</td><td>32</td><td>OHP</td><td>2782.4014</td><td>9.60</td><td>5.00</td><td>31</td><td>OHP</td></tr><tr><td>2606.7163</td><td>9.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2785.9500</td><td>29.50</td><td>1.81</td><td>215</td><td>HIRES</td></tr><tr><td>2622.7200</td><td>33.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2786.9510</td><td>28.40</td><td>1.76</td><td>215</td><td>HIRES</td></tr><tr><td>2623.6940</td><td>10.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2789.8280</td><td>26.75</td><td>1.68</td><td>216</td><td>HIRES</td></tr><tr><td>2631.6864</td><td>-9.40</td><td>5.00</td><td>31</td><td>OHP</td><td>2792.8320</td><td>13.20</td><td>1.50</td><td>214</td><td>HIRES</td></tr><tr><td>2632.0840</td><td>-18.11</td><td>2.67</td><td>215</td><td>HIRES</td><td>2799.9740</td><td>-0.47</td><td>3.26</td><td>75</td><td>HIRES</td></tr><tr><td>2648.6661</td><td>21.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2806.8260</td><td>-15.09</td><td>3.04</td><td>75</td><td>HIRES</td></tr><tr><td>2657.1140</td><td>22.26</td><td>2.23</td><td>214</td><td>HIRES</td><td>2809.8690</td><td>-12.29</td><td>2.97</td><td>75</td><td>HIRES</td></tr><tr><td>2661.0570</td><td>1.01</td><td>2.49</td><td>214</td><td>HIRES</td><td>2812.7900</td><td>-13.26</td><td>2.76</td><td>76</td><td>HIRES</td></tr><tr><td>2662.6285</td><td>-8.40</td><td>3.00</td><td>55</td><td>OHP</td><td>2822.8680</td><td>38.68</td><td>2.66</td><td>75</td><td>HIRES</td></tr><tr><td>2672.0130</td><td>-16.19</td><td>2.20</td><td>216</td><td>HIRES</td><td>2824.4013</td><td>12.60</td><td>5.00</td><td>31</td><td>OHP</td></tr><tr><td>2681.0310</td><td>48.24</td><td>2.47</td><td>174</td><td>HIRES</td><td>2826.7810</td><td>14.72</td><td>2.53</td><td>76</td><td>HIRES</td></tr><tr><td>2683.6198</td><td>14.60</td><td>5.00</td><td>30</td><td>OHP</td><td>2828.7990</td><td>16.66</td><td>3.02</td><td>76</td><td>HIRES</td></tr><tr><td>2687.6359</td><td>8.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2829.7400</td><td>4.81</td><td>3.26</td><td>76</td><td>HIRES</td></tr><tr><td>2690.0110</td><td>9.61</td><td>2.11</td><td>216</td><td>HIRES</td><td>2831.8460</td><td>5.52</td><td>2.79</td><td>73</td><td>HIRES</td></tr><tr><td>2695.0180</td><td>-11.66</td><td>2.05</td><td>214</td><td>HIRES</td><td>2833.7860</td><td>-4.98</td><td>2.74</td><td>75</td><td>HIRES</td></tr><tr><td>2700.9730</td><td>-44.11</td><td>2.09</td><td>213</td><td>HIRES</td><td>2834.7620</td><td>-10.65</td><td>2.87</td><td>76</td><td>HIRES</td></tr><tr><td>2710.0270</td><td>-11.04</td><td>1.76</td><td>215</td><td>HIRES</td><td>2835.7440</td><td>-8.97</td><td>2.69</td><td>77</td><td>HIRES</td></tr><tr><td>2712.0200</td><td>3.27</td><td>2.27</td><td>163</td><td>HIRES</td><td>2838.8030</td><td>-8.72</td><td>2.87</td><td>75</td><td>HIRES</td></tr><tr><td>2712.9410</td><td>18.14</td><td>1.87</td><td>215</td><td>HIRES</td><td>2840.7390</td><td>-10.70</td><td>3.22</td><td>76</td><td>HIRES</td></tr><tr><td>2713.5387</td><td>6.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2858.7850</td><td>37.97</td><td>3.17</td><td>74</td><td>HIRES</td></tr><tr><td>2715.9930</td><td>-2.48</td><td>1.79</td><td>214</td><td>HIRES</td><td>2879.2694</td><td>-27.40</td><td>5.00</td><td>31</td><td>OHP</td></tr><tr><td>2725.5476</td><td>17.60</td><td>5.00</td><td>31</td><td>OHP</td><td>2892.3184</td><td>-9.40</td><td>5.00</td><td>33</td><td>OHP</td></tr><tr><td>2730.0720</td><td>5.38</td><td>1.82</td><td>214</td><td>HIRES</td><td>3009.6665</td><td>-5.40</td><td>6.00</td><td>31</td><td>OHP</td></tr><tr><td>2736.5730</td><td>-32.40</td><td>5.00</td><td>31</td><td>OHP</td><td>3043.0680</td><td>-20.97</td><td>2.65</td><td>143</td><td>HIRES</td></tr><tr><td>2738.9790</td><td>-17.85</td><td>1.54</td><td>215</td><td>HIRES</td><td>3045.0940</td><td>-38.01</td><td>2.18</td><td>216</td><td>HIRES</td></tr><tr><td>2741.9610</td><td>-19.38</td><td>1.51</td><td>215</td><td>HIRES</td><td>3047.1150</td><td>-36.49</td><td>2.07</td><td>202</td><td>HIRES</td></tr><tr><td>2748.8440</td><td>-2.73</td><td>1.63</td><td>204</td><td>HIRES</td><td>3068.0640</td><td>-7.83</td><td>2.08</td><td>153</td><td>HIRES</td></tr><tr><td>2749.8980</td><td>-0.83</td><td>1.48</td><td>216</td><td>HIRES</td><td>3070.9230</td><td>-31.26</td><td>2.01</td><td>214</td><td>HIRES</td></tr><tr><td>2750.4574</td><td>12.60</td><td>4.00</td><td>42</td><td>OHP</td><td>3089.8800</td><td>-9.88</td><td>1.75</td><td>216</td><td>HIRES</td></tr><tr><td>2768.9970</td><td>1.38</td><td>1.57</td><td>215</td><td>HIRES</td><td>3101.0680</td><td>15.11</td><td>1.51</td><td>210</td><td>HIRES</td></tr><tr><td>2769.9960</td><td>-2.72</td><td>1.89</td><td>214</td><td>HIRES</td><td>3108.0060</td><td>-9.65</td><td>1.53</td><td>214</td><td>HIRES</td></tr><tr><td>2773.4748</td><td>-11.40</td><td>6.00</td><td>31</td><td>OHP</td><td>3122.0480</td><td>5.74</td><td>1.63</td><td>171</td><td>HIRES</td></tr><tr><td>2775.9770</td><td>-6.20</td><td>1.58</td><td>216</td><td>HIRES</td><td>3132.8510</td><td>6.65</td><td>1.51</td><td>216</td><td>HIRES</td></tr><tr><td>2776.8030</td><td>-11.86</td><td>1.58</td><td>210</td><td>HIRES</td><td>3138.9630</td><td>5.20</td><td>1.86</td><td>196</td><td>HIRES</td></tr><tr><td>2779.9810</td><td>3.22</td><td>1.69</td><td>214</td><td>HIRES</td><td>3148.0280</td><td>1.62</td><td>1.68</td><td>206</td><td>HIRES</td></tr></tbody></table>

Note. The reported RV values were systematically shifted to be centered around 0, as determined using a Keplerian fit. The low resolution, iodine-free reconnaissance spectra obtained with Keck/HIRES are not included.

CCD charge transfer inefficiency correction (Bouchy et al. 2013). Following the method described, e.g., in Pollacco et al. (2008) and Hebrard et al. (2008), we estimated and corrected for the sky background contamination (mainly due to the Moon) using the second SOPHIE fiber aperture, which is targeted 2' away from the first one pointing toward the star. We estimated that only one of the 20 exposures was significantly polluted by sky background. We note that the analysis of the SOPHIE spectra showed no evidence of a companion star. The extracted HIRES and SOPHIE RV observations are listed in Table 1 and their implications for the system's architecture discussed in Section 6.

The RV data allowed us to place an upper mass limit on planet c and revealed an additional non-transiting planet candidate (with a period of ~34 days; planet b) in the system, as discussed in Section 6.

## 3.2. Imaging Observations and Companion Star.

Archival data listed in the Washington Double Star Catalog (WDS) show that TOI 4633 has a bound companion star.

Hereafter, the stars will be referred to as star A and star B for the more and less massive star, respectively. The 10 archival observations were obtained between 1905 February and 2011 June and their position angles (PA) and angular separations ($\rho$) are listed in Table 2.

In addition to archival data, we obtained high-contrast imaging observations of TOI 4633 to further constrain the orbit of the two stars. High-contrast imaging was performed using the NESSI high-resolution speckle imaging instrument mounted on the 3.5 m WIYN telescope located at Kitt Peak National Observatory (Scott et al. 2018) on 2022 April 21; the PHARO adaptive optics instrument mounted on 5.1 m Hale telescope at Palomar Observatory (Hayward et al. 2001) on 2021 November 11; the 'Alopeke speckle instrument mounted on the 8.1 m Gemini North telescope on Maunakea (Howell et al. 2011; Matson et al. 2019) on 2022 May 9; and the NIRC2 adaptive optics instrument mounted on the 10 m Keck II (Wizinowich et al. 2000) located on Maunakea, Hawaii, on 2023 April 26. While the former two observations were unable to resolve the two stars, the 'Alopeke and NIRC2 observations
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Contrast curves showing the 5σ detection sensitivity and speckle auto-correlation functions obtained using Gemini North/‘Alopeke (left) and Keck II/ NIRC2 (right).</figcaption></figure>

<table><caption>**Table 2**
Archival and new Imaging of TOI 4633</caption><thead><tr><th>Date</th><th>PA (deg)</th><th>ρ (arcsec)</th><th>Δ mag</th><th>Aperture (m)</th><th>λ (nm)</th><th>Method<sup>a</sup></th><th>Reference</th><th>Notes</th></tr></thead><tbody><tr><td>1905.170</td><td>124.7</td><td>0.42</td><td>0.5</td><td>0.9</td><td>...</td><td>M</td><td>Hussey (1905)</td><td>WDS</td></tr><tr><td>1922.876</td><td>124.1</td><td>0.38</td><td>0.4</td><td>1.0</td><td>...</td><td>M</td><td>van Biesbroeck (1927)</td><td>WDS</td></tr><tr><td>1959.440</td><td>124.9</td><td>0.4</td><td>0.3</td><td>0.4</td><td>...</td><td>M</td><td>Couteau (1960)</td><td>WDS</td></tr><tr><td>1963.790</td><td>124.8</td><td>0.47</td><td>0.2</td><td>0.4</td><td>...</td><td>M</td><td>Baize (1967)</td><td>WDS</td></tr><tr><td>1966.672</td><td>122.8</td><td>0.42</td><td>0.2</td><td>0.7</td><td>...</td><td>M</td><td>Worley (1972)</td><td>WDS</td></tr><tr><td>1976.460</td><td>122.7</td><td>0.42</td><td>...</td><td>0.5</td><td>...</td><td>M</td><td>Muller (1978)</td><td>WDS</td></tr><tr><td>1982.078</td><td>124</td><td>0.3</td><td>0.3</td><td>0.7</td><td>...</td><td>M</td><td>Worley (1989)</td><td>WDS</td></tr><tr><td>1984.520</td><td>127.8</td><td>0.28</td><td>...</td><td>0.6</td><td>...</td><td>M</td><td>Heintz (1985)</td><td>WDS</td></tr><tr><td>1993.100</td><td>119</td><td>0.26</td><td>...</td><td>0.5</td><td>...</td><td>M</td><td>Muller (1997)</td><td>WDS</td></tr><tr><td>2021.8616</td><td>...</td><td>&lt;0.1</td><td>...</td><td>5.1</td><td>1000–2500</td><td>AO</td><td>...</td><td>Palomar/PHARO</td></tr><tr><td>2022.303</td><td>...</td><td>&lt;0.06</td><td>...</td><td>3.5</td><td>812–852</td><td>S</td><td>...</td><td>WIYN/NESSI</td></tr><tr><td>2022.352</td><td>300.5 ± 0.5<sup>b</sup></td><td>0.05 ± 0.01</td><td>0.4</td><td>8.1</td><td>535–589</td><td>S</td><td>...</td><td>Gemini/‘Alopeke</td></tr><tr><td>2023.316</td><td>303.18 ± 1.29</td><td>0.062 ± 0.01</td><td>0.23</td><td>10</td><td>2256–2285</td><td>AO</td><td>...</td><td>Keck II/NIRC2</td></tr></tbody></table>

**Notes.**

revealed the close companion at angular separations of 51 and 62 mas, respectively. The observations showed that the two stars have a brightness ratio $F_B/F_A = 0.7$, as shown in Figure 3. The former two instruments were unable to resolve the two stars due to a lack of angular resolution to resolve a separation of <100 mas. Similarly, neither Gaia DR2 nor Gaia DR3, which obtained observations of the target in 2015 and 2016, respectively, were able to resolve the two stars due to the small angular separation of the stars at the time of the observations. Furthermore, due to the proximity of the two stars, the Gaia DR3 noise metric RUWE~1.15 is consistent with a single-star solution, which highlights the need for high-resolution imaging. Finally, we note that in addition to the archival data listed in Table 2, WDS records one additional observation obtained by Gili et al. (2021) in 2011. Their recorded values of position angle (37.6 deg) and angular separation (0.′′233) are

inconsistent with both the remainder of the archival and our newly obtained data. More specifically, there is no stable binary orbit that would be able to explain this measurement in combination with the archival and newly obtained measurements. Furthermore, the analysis method that was used to reduce these observations follows a unique procedure that is ambiguously defined and distinct from the processes used for any of the other imaging data, and thus we do not include it in any further analysis.

Given the available Gaia data, the proper motions of the two stars cannot be disentangled. However, given the proper motion of the system (R.A.: −23.056, decl.: −67.0188 mas yr⁻¹) and the range of angular separations of the two stars reported between 1905 and 2023 (470–50 mas), we confirm that these two stars are bound and not a chance alignment along our line of sight.

<footnote><sup>a</sup> M indicates that the observation were obtained using a micrometer on a refractor telescope; S indicates that the observations were obtained with the Speckle technique; and AO indicated the use of adaptive optics.</footnote>

<footnote><sup>b</sup> 180° ambiguity in the PA.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Contrast curves showing the 5σ detection sensitivity and speckle auto-correlation functions obtained using Gemini North/‘Alopeke (left) and Keck II/ NIRC2 (right).</figcaption></figure>

<table><caption>**Table 2**
Archival and new Imaging of TOI 4633</caption><thead><tr><th>Date</th><th>PA (deg)</th><th>ρ (arcsec)</th><th>Δ mag</th><th>Aperture (m)</th><th>λ (nm)</th><th>Method<sup>a</sup></th><th>Reference</th><th>Notes</th></tr></thead><tbody><tr><td>1905.170</td><td>124.7</td><td>0.42</td><td>0.5</td><td>0.9</td><td>...</td><td>M</td><td>Hussey (1905)</td><td>WDS</td></tr><tr><td>1922.876</td><td>124.1</td><td>0.38</td><td>0.4</td><td>1.0</td><td>...</td><td>M</td><td>van Biesbroeck (1927)</td><td>WDS</td></tr><tr><td>1959.440</td><td>124.9</td><td>0.4</td><td>0.3</td><td>0.4</td><td>...</td><td>M</td><td>Couteau (1960)</td><td>WDS</td></tr><tr><td>1963.790</td><td>124.8</td><td>0.47</td><td>0.2</td><td>0.4</td><td>...</td><td>M</td><td>Baize (1967)</td><td>WDS</td></tr><tr><td>1966.672</td><td>122.8</td><td>0.42</td><td>0.2</td><td>0.7</td><td>...</td><td>M</td><td>Worley (1972)</td><td>WDS</td></tr><tr><td>1976.460</td><td>122.7</td><td>0.42</td><td>...</td><td>0.5</td><td>...</td><td>M</td><td>Muller (1978)</td><td>WDS</td></tr><tr><td>1982.078</td><td>124</td><td>0.3</td><td>0.3</td><td>0.7</td><td>...</td><td>M</td><td>Worley (1989)</td><td>WDS</td></tr><tr><td>1984.520</td><td>127.8</td><td>0.28</td><td>...</td><td>0.6</td><td>...</td><td>M</td><td>Heintz (1985)</td><td>WDS</td></tr><tr><td>1993.100</td><td>119</td><td>0.26</td><td>...</td><td>0.5</td><td>...</td><td>M</td><td>Muller (1997)</td><td>WDS</td></tr><tr><td>2021.8616</td><td>...</td><td>&lt;0.1</td><td>...</td><td>5.1</td><td>1000–2500</td><td>AO</td><td>...</td><td>Palomar/PHARO</td></tr><tr><td>2022.303</td><td>...</td><td>&lt;0.06</td><td>...</td><td>3.5</td><td>812–852</td><td>S</td><td>...</td><td>WIYN/NESSI</td></tr><tr><td>2022.352</td><td>300.5 ± 0.5<sup>b</sup></td><td>0.05 ± 0.01</td><td>0.4</td><td>8.1</td><td>535–589</td><td>S</td><td>...</td><td>Gemini/‘Alopeke</td></tr><tr><td>2023.316</td><td>303.18 ± 1.29</td><td>0.062 ± 0.01</td><td>0.23</td><td>10</td><td>2256–2285</td><td>AO</td><td>...</td><td>Keck II/NIRC2</td></tr></tbody></table>

**Notes.**

<sup>a</sup> M indicates that the observation were obtained using a micrometer on a refractor telescope; S indicates that the observations were obtained with the Speckle technique; and AO indicated the use of adaptive optics.

<sup>b</sup> 180° ambiguity in the PA.

revealed the close companion at angular separations of 51 and 62 mas, respectively. The observations showed that the two stars have a brightness ratio $F_B/F_A = 0.7$, as shown in Figure 3. The former two instruments were unable to resolve the two stars due to a lack of angular resolution to resolve a separation of <100 mas. Similarly, neither Gaia DR2 nor Gaia DR3, which obtained observations of the target in 2015 and 2016, respectively, were able to resolve the two stars due to the small angular separation of the stars at the time of the observations. Furthermore, due to the proximity of the two stars, the Gaia DR3 noise metric RUWE~1.15 is consistent with a single-star solution, which highlights the need for high-resolution imaging. Finally, we note that in addition to the archival data listed in Table 2, WDS records one additional observation obtained by Gili et al. (2021) in 2011. Their recorded values of position angle (37.6 deg) and angular separation (0.′′233) are

inconsistent with both the remainder of the archival and our newly obtained data. More specifically, there is no stable binary orbit that would be able to explain this measurement in combination with the archival and newly obtained measurements. Furthermore, the analysis method that was used to reduce these observations follows a unique procedure that is ambiguously defined and distinct from the processes used for any of the other imaging data, and thus we do not include it in any further analysis.

Given the available Gaia data, the proper motions of the two stars cannot be disentangled. However, given the proper motion of the system (R.A.: −23.056, decl.: −67.0188 mas yr⁻¹) and the range of angular separations of the two stars reported between 1905 and 2023 (470–50 mas), we confirm that these two stars are bound and not a chance alignment along our line of sight.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Generalized Lomb-Scargle Periodogram of the HIRES and SOPHIE RV data (navy and maroon bold panels, respectively). The non bold panels show periodograms of various activity indicators for each instrument, while the vertical gray line indicates a period of 34 days. The horizontal dotted orange and dashed gray lines represent the 10% and 30% FAPs, respectively.</figcaption></figure>

**4. Activity Indicators and Statistical Validation of the Transit Signals**

In this section we discuss how the stellar activity affects our confidence in the planetary nature of planet candidate b. Furthermore, we present statistical validation of planet c.

**4.1. *Stellar Activity Indicators and Planet Candidate TOI-4633 b***

The RV data exhibit a periodic signal with a period of ~34 days. This can be seen in Figure 4, which shows periodograms of the RV observations obtained with HIRES (navy) and SOPHIE (maroon), calculated using the

PyAstronomy Generalized Lomb-Scargle Periodogram (Czesla et al. 2019). We carried out a number of tests in order to determine whether this periodic signal seen in the RV data is caused by a planetary body or by stellar activity. This is important as magnetic stellar activity manifests itself by producing brighter or darker regions on the surface of the star, such as faculae, starspots, and plages. This, in turn, affects the observed stellar spectra and can induce RV signals that mimic those from planets (e.g., Queloz et al. 2001; Figueira et al. 2010; Boisse et al. 2011; Díaz et al. 2016).

For the HIRES data we investigate the chromospheric stellar activity using the S-index and its corresponding derivative $log R'_{HK}$ value, a parameter that provides a proxy for the level
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Generalized Lomb-Scargle Periodogram of the HIRES and SOPHIE RV data (navy and maroon bold panels, respectively). The non bold panels show periodograms of various activity indicators for each instrument, while the vertical gray line indicates a period of 34 days. The horizontal dotted orange and dashed gray lines represent the 10% and 30% FAPs, respectively.</figcaption></figure>

**4. Activity Indicators and Statistical Validation of the Transit Signals**

In this section we discuss how the stellar activity affects our confidence in the planetary nature of planet candidate b. Furthermore, we present statistical validation of planet c.

**4.1. *Stellar Activity Indicators and Planet Candidate TOI 4633 b***

The RV data exhibit a periodic signal with a period of ~34 days. This can be seen in Figure 4, which shows periodograms of the RV observations obtained with HIRES (navy) and SOPHIE (maroon), calculated using the

PyAstronomy Generalized Lomb-Scargle Periodogram (Czesla et al. 2019). We carried out a number of tests in order to determine whether this periodic signal seen in the RV data is caused by a planetary body or by stellar activity. This is important as magnetic stellar activity manifests itself by producing brighter or darker regions on the surface of the star, such as faculae, starspots, and plages. This, in turn, affects the observed stellar spectra and can induce RV signals that mimic those from planets (e.g., Queloz et al. 2001; Figueira et al. 2010; Boisse et al. 2011; Díaz et al. 2016).

For the HIRES data we investigate the chromospheric stellar activity using the S-index and its corresponding derivative $\log R'_{HK}$ value, a parameter that provides a proxy for the level
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Normalized light curve (left), Generalized Lomb-Scargle periodogram (middle), and light curve phase folded at the period of planet candidate b ($P_b$ = 34.15 ± 0.15 days) for the 2 minutes cadence TESS data (top) and the V-band ASAS-SN data. The vertical blue line in the middle panel indicates the period of planet candidate b, while the horizontal dotted orange and dashed gray lines represent the 10% and 30% FAPs, respectively.</figcaption></figure>

of magnetic activity on the surface of a star (e.g., Vaughan et al. 1978; Noyes et al. 1984; Saar et al. 1998; Santos et al. 2000, 2010; Cincunegui et al. 2007). In brief, $R'_{HK}$ is the ratio of the emission in the cores of the Ca II H and K lines (at 3933 Å and 3968 Å) to the total bolometric flux of the star, where the Ca II H and K lines probe the temperature of the chromosphere of the star. As magnetic activity primarily heats up the chromosphere (resulting in stronger emission features), the ratio of the intensity of the emission feature in these lines to the bolometric flux probes the magnetic activity of the star.

We used the $H\alpha$ Balmer line of hydrogen at 6563 Å as an additional indicator of the chromospheric activity of the star. Using the HIRES data and following the methodology described by Gomes da Silva et al. (2011) we measured the ratio of the flux within ±0.8 Å of the $H\alpha$ line at 6562.808 Å to the flux of the two wavelength regions of 6550.87 ± 5.375 Å and 6580.31 ± 4.375 Å. Prior to measuring the flux ratios, we derived a wavelength solution for the $H\alpha$ order of the spectrum and shifted each spectrum to the star's rest frame by cross correlating each spectra with the spectrum from the NSO solar atlas. Both the NSO solar atlas spectrum and the HIRES spectra were re-sampled to have wavelength steps of 0.003 Å. The spectra were normalized with a third order polynomial (where the $H\alpha$ line ±5 Å was masked out).

Similarly, we computed activity indicators using the SOPHIE data by measuring the ratio of the fluxes in the cores of the $H\alpha$ (6562.808 Å), $H\beta$ (4861.363 Å) and Na I D1 & D2 (5895.92 & 5889.95 Å, respectively) lines to the flux in continuum regions around each of the lines. The same wavelength regions as described above were used for the $H\alpha$ index. For the $H\beta$ index we measured the flux within a 1.2 Å window centered on the absorption line and divided by the two reference regions defined as 4855.0 ± 2.5 Å and 4870.0 ± 2.5 Å (Klein et al. 2022). The Na I index was defined as the flux within a 0.5 Å window centered on each line divided by the two reference regions of 5805.0 ± 5 Å and 6090.0 ± 10 Å (Gomes da Silva et al. 2011). In addition to these stellar activity indicators, we made use of the fact that periodic changes in the shapes of absorption lines (e.g., FWHM, BIS) can be indicative of whether RV signals are caused by companions or by stellar activity. We note that different

activity indicators are used for the HIRES and SOPHIE data due to the different spectral ranges that both of these instruments cover.

We used the PyAstronomy Generalized Lomb-Scargle Periodogram (Czesla et al. 2019) to search for periodicities in the chromospheric activity indicators, FWHM and BIS. As shown in Figure 4, where the dotted orange and dashed gray lines show the 10% and 30% False Alarm Probability (FAP), respectively, there are no significant trends that coincide with the periodic signals seen in the RV data (indicated by the gray vertical line). We do note that the S-index shows a local power maximum at 37 days, which is near the radial velocity power maximum at 34 days. This power maximum is, however, not the largest peak in the S-index power spectrum and has an FAP of 99.96%.

We also determine the degree of correlation between the RV time series and each of the activity indicators (independent for the HIRES and SOPHIE data) by computing Pearson's correlation coefficients *r* and their *p*-values. The test showed that there are no correlations between the RVs and any of the activity indicators with *r* ranging from -0.28 to 0.14 and *p*-values ranging from 0.23 to 0.29. As all of the *p*-values are significantly greater than 0.05, we conclude that there is no significant evidence for a correlation between the RV time series and the activity indicators for both the HIRES and SOPHIE data.

Finally, we ensured that the PyAstronomy Generalized Lomb-Scargle periodograms of the TESS light curve and the ASAS-SN light curve (where the transit events in the TESS data are removed) show no significant periodic signals. The Lomb-Scargle periodograms of the TESS and ASAS-SN data are shown in Figure 5, alongside the light curves phase folded at the period of TOI 4633 b ($P$ = 34.15 ± 0.15 days). This shows that there are no significant trends at this period with an FAP better than 30% that could result from additional transit events, stellar rotation, or asteroseismic pulsations in either the TESS or the ASAS-SN data. This is further supported by the lack of a coherent signal in the phase folded data at the period of planet candidate b in the rightmost column of Figure 5.

In summary, based on our analysis of the stellar activity indicators we have been able to extract from the HIRES and
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Normalized light curve (left), Generalized Lomb-Scargle periodogram (middle), and light curve phase folded at the period of planet candidate b ($P_b$ = 34.15 ± 0.15 days) for the 2 minutes cadence TESS data (top) and the V-band ASAS-SN data. The vertical blue line in the middle panel indicates the period of planet candidate b, while the horizontal dotted orange and dashed gray lines represent the 10% and 30% FAPs, respectively.</figcaption></figure>

of magnetic activity on the surface of a star (e.g., Vaughan et al. 1978; Noyes et al. 1984; Saar et al. 1998; Santos et al. 2000, 2010; Cincunegui et al. 2007). In brief, $R'_{HK}$ is the ratio of the emission in the cores of the Ca II H and K lines (at 3933 Å and 3968 Å) to the total bolometric flux of the star, where the Ca II H and K lines probe the temperature of the chromosphere of the star. As magnetic activity primarily heats up the chromosphere (resulting in stronger emission features), the ratio of the intensity of the emission feature in these lines to the bolometric flux probes the magnetic activity of the star.

We used the $H\alpha$ Balmer line of hydrogen at 6563 Å as an additional indicator of the chromospheric activity of the star. Using the HIRES data and following the methodology described by Gomes da Silva et al. (2011) we measured the ratio of the flux within ±0.8 Å of the $H\alpha$ line at 6562.808 Å to the flux of the two wavelength regions of 6550.87 ± 5.375 Å and 6580.31 ± 4.375 Å. Prior to measuring the flux ratios, we derived a wavelength solution for the $H\alpha$ order of the spectrum and shifted each spectrum to the star's rest frame by cross correlating each spectra with the spectrum from the NSO solar atlas. Both the NSO solar atlas spectrum and the HIRES spectra were re-sampled to have wavelength steps of 0.003 Å. The spectra were normalized with a third order polynomial (where the $H\alpha$ line ±5 Å was masked out).

Similarly, we computed activity indicators using the SOPHIE data by measuring the ratio of the fluxes in the cores of the $H\alpha$ (6562.808 Å), $H\beta$ (4861.363 Å) and Na I D1 & D2 (5895.92 & 5889.95 Å, respectively) lines to the flux in continuum regions around each of the lines. The same wavelength regions as described above were used for the $H\alpha$ index. For the $H\beta$ index we measured the flux within a 1.2 Å window centered on the absorption line and divided by the two reference regions defined as 4855.0 ± 2.5 Å and 4870.0 ± 2.5 Å (Klein et al. 2022). The Na I index was defined as the flux within a 0.5 Å window centered on each line divided by the two reference regions of 5805.0 ± 5 Å and 6090.0 ± 10 Å (Gomes da Silva et al. 2011). In addition to these stellar activity indicators, we made use of the fact that periodic changes in the shapes of absorption lines (e.g., FWHM, BIS) can be indicative of whether RV signals are caused by companions or by stellar activity. We note that different

activity indicators are used for the HIRES and SOPHIE data due to the different spectral ranges that both of these instruments cover.

We used the PyAstronomy Generalized Lomb-Scargle Periodogram (Czesla et al. 2019) to search for periodicities in the chromospheric activity indicators, FWHM and BIS. As shown in Figure 4, where the dotted orange and dashed gray lines show the 10% and 30% False Alarm Probability (FAP), respectively, there are no significant trends that coincide with the periodic signals seen in the RV data (indicated by the gray vertical line). We do note that the S-index shows a local power maximum at 37 days, which is near the radial velocity power maximum at 34 days. This power maximum is, however, not the largest peak in the S-index power spectrum and has an FAP of 99.96%.

We also determine the degree of correlation between the RV time series and each of the activity indicators (independent for the HIRES and SOPHIE data) by computing Pearson's correlation coefficients *r* and their *p*-values. The test showed that there are no correlations between the RVs and any of the activity indicators with *r* ranging from -0.28 to 0.14 and *p*-values ranging from 0.23 to 0.29. As all of the *p*-values are significantly greater than 0.05, we conclude that there is no significant evidence for a correlation between the RV time series and the activity indicators for both the HIRES and SOPHIE data.

Finally, we ensured that the PyAstronomy Generalized Lomb-Scargle periodograms of the TESS light curve and the ASAS-SN light curve (where the transit events in the TESS data are removed) show no significant periodic signals. The Lomb-Scargle periodograms of the TESS and ASAS-SN data are shown in Figure 5, alongside the light curves phase folded at the period of TOI 4633 b ($P$ = 34.15 ± 0.15 days). This shows that there are no significant trends at this period with an FAP better than 30% that could result from additional transit events, stellar rotation, or asteroseismic pulsations in either the TESS or the ASAS-SN data. This is further supported by the lack of a coherent signal in the phase folded data at the period of planet candidate b in the rightmost column of Figure 5.

In summary, based on our analysis of the stellar activity indicators we have been able to extract from the HIRES and
SOPHIE spectra, there is no indication that the 34 days signal is caused by stellar activity, and we therefore consider the hypothesis that this signal is of planetary origin to be more likely. However, the limited precision of the HIRES activity indicators, and the limited number of SOPHIE observations, preclude a joint analysis of the RVs and activity indicators to disentangle the contributions of activity and planet(s) to the RVs. Consequently, we report the 34 days signal as a planet candidate, rather than a confirmed planet detection.

## 4.2. Statistical Validation of Planet c

The standard diagnostic tests presented in Section 2.1 allowed us to rule out a number of instrumental and astrophysical false positive scenarios that could mimic the transit-like signals seen in the TESS data. Furthermore, the spectroscopic follow-up observations presented in Sections 3.1 allowed us to place an upper mass limit of 123 $M_{\oplus}$ (99 percentile of the credible interval) on TOI 4633 c, as discussed in Section 6.1. As we do not have a $3\sigma$ mass measurement of planet c, we carry out a statistical analysis of the likelihood that the transit signals are caused by a planet as opposed to the range of alternative, astrophysical false positive scenarios.

We used the open-source package TRICERATOPS (Giacalone et al. 2021), which was specifically developed to aid in the vetting and validation of transit-like signals identified in the TESS data, to calculate the false positive probability (FPP) of the observed transit signals of TOI 4633 c. In brief, TRICERATOPS uses a Bayesian framework that incorporates prior knowledge of the target star, planet occurrence rates, and stellar multiplicity to calculate the probability that the transit signal is due to a transiting planet. It also makes use of the TRILEGAL (Girardi et al. 2005) galactic model to simulate a population of stars around the line of sight of the target.

The resulting FPP quantifies the probability that the observed transit signal can be attributed to something other than a transiting planet. As inputs to the code we used the 2 minutes cadence and 20 s cadence (where available) data (Section 2), combined with the contrast curves obtained using ‘Alopeke and NIRC2 (Section 3.2). We found the FPP to be $0.0003 \pm 0.0005$ and the nearby FPP (NFFP) to be $6.6 \times 10^{-21} \pm 4.6 \times 10^{-21}$. Both of these are better than the commonly accepted validation threshold of FPP < 0.015 and NFPP < $10^{-3}$, as defined by Giacalone et al. (2021), allowing us to conclude that TOI 4633 c is a non-self-luminous object transiting one of the two stars in the binary. Furthermore, combined with the upper mass limit provided by the RV monitoring, we consider TOI 4633 c to be a confirmed planet.

For the remainder of this paper, we consider the 34 days RV-detected signal to be a likely planet candidate (TOI 4633 b), and the 272 days transit and RV-detected signal to be a confirmed planet (TOI 4633 c).

# 5. Analysis of Stellar System

The combination of archival and new high-resolution spectral imaging data revealed a bound companion star. Given the brightness ratio of the two stars of $F_R/F_A = 0.7$, we expect both stars to contribute to the spectra. However, there are no observable RV shifts from the companion star, detectable changes in the shape of the absorption lines, nor evidence for double lined spectra. As discussed in Section 3.1, HIRES is able to detect companion stars where $\Delta \text{RV} \ge 10 \text{ km s}^{-1}$. We

assume for the remainder of this paper that the obtained HIRES and SOPHIE spectra are a composite of the light from both stars, where the relative RV shift between the two stars is less than $10 \text{ km s}^{-1}$. Similarly, we assume that the TESS light curve is a composite of the light from both stars. In this section we discuss the properties of the two stars and the stellar configuration.

## 5.1. Stellar Parameter Determination

In order to determine stellar parameters for this system, we use precise multi-wavelength photometric measurements as well as high-resolution spectra. We note that while the multi-wavelength photometric fits account for the multiple stars in the system, the spectroscopic solutions only consider one star and are subject to error.

### 5.1.1. Spectroscopic Parameter Determination

To extract stellar parameter values from the spectra, we used a moderate signal-to-noise (per pixel S/N = 46) iodine-free observation obtained as a reconnaissance observation using the HIRES spectrograph on the Keck I telescope (Vogt et al. 1994). We measured the effective temperature ($T_{\text{eff}}$), surface gravity (log g), iron abundance ([Fe/H]), and projected rotational velocity of the star using the tools available in the SpecMatch software package (Petigura 2015). We first corrected the observed wavelengths to be in the observer's rest frame by cross correlating a solar model with the observed spectrum. Then, we fit for $T_{\text{eff}}$, log g, [Fe/H], vsini, and the instrumental point-spread function (PSF) using the underlying Bayesian differential-evolution Markov Chain Monte Carlo (MCMC) machinery of ExoPy (Fulton et al. 2013). At each step in the MCMC chains, a synthetic spectrum is created by interpolating the Coelho (2014) grid of stellar models for a set of $T_{\text{eff}}$, log g, and [Fe/H] values and solar alpha abundance. We convolved this synthetic spectrum with a rotational plus macroturbulence broadening kernel using the prescriptions of Valenti & Fischer (2005) and Hirano et al. (2011). Finally, we performed another convolution with a Gaussian kernel to account for the instrumental PSF, and compared the synthetic spectrum with the observed spectrum to assess the goodness of fit. The priors are uniform in $T_{\text{eff}}$, log g, and [Fe/H], but we assign a Gaussian prior to the instrumental PSF that encompasses the typical variability in the PSF width caused by seeing changes and guiding errors. Five echelle orders of the spectrum were fit separately and the resulting posterior distributions were combined before taking the median values for each parameter. Parameter uncertainties were estimated as the scatter in spectroscopic parameters given by SpecMatch relative to the values for 352 stars in the Valenti & Fischer (2005) sample and 76 stars in the Huber et al. (2013) asteroseismic sample. Systematic trends in SpecMatch values as a function of $T_{\text{eff}}$, log g, and [Fe/H] relative to these benchmark samples were fit for and removed in the final quoted parameter values.

Although the approach described above performs well for main sequence, single-star parameter estimation, calibrations using an empirical spectral library can often result in more robust parameters. As such, we independently determined the stellar parameters of TOI 4633 using SpecMatch-Emp, which follows a similar procedure as that described above but using an empirical library of stellar spectra taken with Keck/HIRES (Yee et al. 2017). While we find that our determinations of $T_{\text{eff}}$, log g, and
SOPHIE spectra, there is no indication that the 34 days signal is caused by stellar activity, and we therefore consider the hypothesis that this signal is of planetary origin to be more likely. However, the limited precision of the HIRES activity indicators, and the limited number of SOPHIE observations, preclude a joint analysis of the RVs and activity indicators to disentangle the contributions of activity and planet(s) to the RVs. Consequently, we report the 34 days signal as a planet candidate, rather than a confirmed planet detection.

## 4.2. Statistical Validation of Planet c

The standard diagnostic tests presented in Section 2.1 allowed us to rule out a number of instrumental and astrophysical false positive scenarios that could mimic the transit-like signals seen in the TESS data. Furthermore, the spectroscopic follow-up observations presented in Sections 3.1 allowed us to place an upper mass limit of 123 $M_{\oplus}$ (99 percentile of the credible interval) on TOI 4633 c, as discussed in Section 6.1. As we do not have a $3\sigma$ mass measurement of planet c, we carry out a statistical analysis of the likelihood that the transit signals are caused by a planet as opposed to the range of alternative, astrophysical false positive scenarios.

We used the open-source package TRICERATOPS (Giacalone et al. 2021), which was specifically developed to aid in the vetting and validation of transit-like signals identified in the TESS data, to calculate the false positive probability (FPP) of the observed transit signals of TOI 4633 c. In brief, TRICERATOPS uses a Bayesian framework that incorporates prior knowledge of the target star, planet occurrence rates, and stellar multiplicity to calculate the probability that the transit signal is due to a transiting planet. It also makes use of the TRILEGAL (Girardi et al. 2005) galactic model to simulate a population of stars around the line of sight of the target.

The resulting FPP quantifies the probability that the observed transit signal can be attributed to something other than a transiting planet. As inputs to the code we used the 2 minutes cadence and 20 s cadence (where available) data (Section 2), combined with the contrast curves obtained using ‘Alopeke and NIRC2 (Section 3.2). We found the FPP to be $0.0003 \pm 0.0005$ and the nearby FPP (NFFP) to be $6.6 \times 10^{-21} \pm 4.6 \times 10^{-21}$. Both of these are better than the commonly accepted validation threshold of FPP < 0.015 and NFPP < $10^{-3}$, as defined by Giacalone et al. (2021), allowing us to conclude that TOI 4633 c is a non-self-luminous object transiting one of the two stars in the binary. Furthermore, combined with the upper mass limit provided by the RV monitoring, we consider TOI 4633 c to be a confirmed planet.

For the remainder of this paper, we consider the 34 days RV-detected signal to be a likely planet candidate (TOI 4633 b), and the 272 days transit and RV-detected signal to be a confirmed planet (TOI 4633 c).

# 5. Analysis of Stellar System

The combination of archival and new high-resolution spectral imaging data revealed a bound companion star. Given the brightness ratio of the two stars of $F_R/F_A = 0.7$, we expect both stars to contribute to the spectra. However, there are no observable RV shifts from the companion star, detectable changes in the shape of the absorption lines, nor evidence for double lined spectra. As discussed in Section 3.1, HIRES is able to detect companion stars where $\Delta \text{RV} \ge 10 \text{ km s}^{-1}$. We

assume for the remainder of this paper that the obtained HIRES and SOPHIE spectra are a composite of the light from both stars, where the relative RV shift between the two stars is less than $10 \text{ km s}^{-1}$. Similarly, we assume that the TESS light curve is a composite of the light from both stars. In this section we discuss the properties of the two stars and the stellar configuration.

## 5.1. Stellar Parameter Determination

In order to determine stellar parameters for this system, we use precise multi-wavelength photometric measurements as well as high-resolution spectra. We note that while the multi-wavelength photometric fits account for the multiple stars in the system, the spectroscopic solutions only consider one star and are subject to error.

### 5.1.1. Spectroscopic Parameter Determination

To extract stellar parameter values from the spectra, we used a moderate signal-to-noise (per pixel $S/N = 46$) iodine-free observation obtained as a reconnaissance observation using the HIRES spectrograph on the Keck I telescope (Vogt et al. 1994). We measured the effective temperature ($T_{\text{eff}}$), surface gravity (log g), iron abundance ([Fe/H]), and projected rotational velocity of the star using the tools available in the SpecMatch software package (Petigura 2015). We first corrected the observed wavelengths to be in the observer's rest frame by cross correlating a solar model with the observed spectrum. Then, we fit for $T_{\text{eff}}$, log g, [Fe/H], vsini, and the instrumental point-spread function (PSF) using the underlying Bayesian differential-evolution Markov Chain Monte Carlo (MCMC) machinery of ExoPy (Fulton et al. 2013). At each step in the MCMC chains, a synthetic spectrum is created by interpolating the Coelho (2014) grid of stellar models for a set of $T_{\text{eff}}$, log g, and [Fe/H] values and solar alpha abundance. We convolved this synthetic spectrum with a rotational plus macroturbulence broadening kernel using the prescriptions of Valenti & Fischer (2005) and Hirano et al. (2011). Finally, we performed another convolution with a Gaussian kernel to account for the instrumental PSF, and compared the synthetic spectrum with the observed spectrum to assess the goodness of fit. The priors are uniform in $T_{\text{eff}}$, log g, and [Fe/H], but we assign a Gaussian prior to the instrumental PSF that encompasses the typical variability in the PSF width caused by seeing changes and guiding errors. Five echelle orders of the spectrum were fit separately and the resulting posterior distributions were combined before taking the median values for each parameter. Parameter uncertainties were estimated as the scatter in spectroscopic parameters given by SpecMatch relative to the values for 352 stars in the Valenti & Fischer (2005) sample and 76 stars in the Huber et al. (2013) asteroseismic sample. Systematic trends in SpecMatch values as a function of $T_{\text{eff}}$, log g, and [Fe/H] relative to these benchmark samples were fit for and removed in the final quoted parameter values.

Although the approach described above performs well for main sequence, single-star parameter estimation, calibrations using an empirical spectral library can often result in more robust parameters. As such, we independently determined the stellar parameters of TOI 4633 using SpecMatch-Emp, which follows a similar procedure as that described above but using an empirical library of stellar spectra taken with Keck/HIRES (Yee et al. 2017). While we find that our determinations of $T_{\text{eff}}$, log g, and
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Spectral energy distribution of TIC 307958020. Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the passband. Blue symbols are the model fluxes from the best-fit Kurucz atmosphere model (gray) for the combined light of the binary star system. Cyan and red curves represent the Kurucz atmosphere models corresponding to the warmer (larger) and cooler (smaller) stellar components of the system, respectively.</figcaption></figure>

[Fe/H] are in excellent agreement between the two methods within their respective uncertainties, the determined values of the stellar radius are discrepant ($R_{\text{SpecMatch-Syn}} = 1.50 \pm 0.04 R_⊙$ and $R_{\text{SpecMatch-Emp}} = 1.11 \pm 0.18 R_⊙$). In addition to the discrepancy between these two methods, both are affected by the assumption that the light originates from a single star. We, therefore, independently derived the stellar parameters using the spectral energy distribution (SED), as this is able to account for the light contributions of both stars in the system.

### 5.1.2. SED Fitting

As an independent determination of the stellar parameters, we performed an analysis of the broadband SED of the star together with the Gaia DR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016); Stassun et al. (2017, 2018). We pulled the $B_T V_T$ magnitudes from Tycho-2, the JHK$_S$ magnitudes from 2MASS, the W1–W4 magnitudes from WISE, the $G$, $G_{BP}$, and $G_{RP}$ magnitudes from Gaia, and the FUV and NUV magnitudes from GALEX. Together, the available photometry spans the full stellar SED over the wavelength range 0.2–22 µm (see Figure 6).

For an initial fit, we treated the SED as arising from a single star, with the fitted parameters being the effective temperature ($T_{\text{eff}}$), surface gravity (log g), and metallicity ([Fe/H]), as well as the extinction $A_V$, which we limited to the maximum line-of-sight value from the Galactic dust maps of Schlegel et al. (1998). For the initial single-star model fit, we adopted the $T_{\text{eff}}$, log g, and [Fe/H] from the spectroscopically determined values. The goodness of fit was, not surprisingly, only marginally good ($\chi^2 = 4.5$) and, given the stringent Gaia distance, implied an oversized star with $R_* \approx 1.5 R_⊙$.

Next, we performed a two-component fit using the methodology of Stassun & Torres (2016), solving for the component $T_{\text{eff}}$ and $R_*$ by requiring the flux-weighted average $T_{\text{eff}}$ to agree with the spectroscopic value used above while also requiring the flux ratio $F_B/F_A$ in the visible to match that determined from the high-contrast imaging, and of course requiring agreement with the strict Gaia distance.

The result, shown in Figure 6, has $\chi^2 = 1.1$ and best-fit parameters of $A_V = 0.07 \pm 0.02$, $T_{\text{eff,A}} = 5800 \pm 50$ K, $R_A = 1.05 \pm 0.05 R_⊙$, $T_{\text{eff,B}} = 5600 \pm 50$ K, and $R_B = 0.98 \pm 0.05 R_⊙$. Based on the empirical relations of Torres et al. (2010), the hotter/larger star appears to have a mass $M_A = 1.10 \pm 0.06 M_⊙$, and the cooler/smaller star has $M_B = 1.05 \pm 0.06 M_⊙$, consistent with an interpretation of main-sequence stars.

Finally, the vsini measurement reported by APOGEE DR16 ($\text{vsini} = 4.39 \text{ km s}^{-1}$) together with the brighter star's radius implies a maximum rotation period of $12.1 \pm 1.3$ days. With the gyrochronology relations of Mamajek & Hillenbrand (2008), this implies a minimum age of $1.3 \pm 0.3$ Gyr, again consistent with unevolved main-sequence stars. As a comparison, the vsini derived using the SpecMatch software package from the HIRES data is $4.25 \text{ km s}^{-1}$, which is also consistent with the system containing main-sequence stars. We note that the vsini estimates are likely strongly affected by the contamination of the companion star. Finally, as mentioned in Section 4.1, the TESS and ASAS-SN data show no signs of measurable rotation and thus a more precise age of this system cannot be determined using gyrochronology.

As the two-component SED analysis is able to account for the light contribution of both stars, we adopt these stellar parameters for the remainder of the paper. All parameters are listed in Table 3.

## 5.2. Binary Orbit Modeling

The orbital parameters of the binary system (TOI 4633 AB) were determined by Bayesian parameter estimation using the open-source software “ORBITIZE!” (Blunt et al. 2020). The position angles and angular separations extracted from all of the available imaging data (listed in Table 2) were used as input data. Observations obtained prior to 2011 did not report uncertainties on their measurements. We, therefore, adopted large uncertainties for these archival position angles and angular separations of $\pm 10$ deg and $\pm 50$ mas, respectively. As discussed in Section 3.1, the HIRES data allows us to place an upper limit on the relative motion of the two stars of $10 \text{ km s}^{-1}$. As such, RV values of $0 \pm 5 \text{ km s}^{-1}$ at the times of the HIRES observations are used as input RV measurements for the ORBITIZE! model.

We used the parallel-tempered Affine-invariant sampler PTEMCEE (Foreman-Mackey et al. 2013; Vousden et al. 2016) and adopted priors of $10.55 \pm 0.013$ mas (Gaia Collaboration et al. 2021) for the parallax and $1.10 \pm 0.2 M_⊙$ and $1.05 \pm 0.2 M_⊙$ for the masses of the primary and secondary stars, respectively (see Section 5.1.2). Due to the lack of more constraining RV measurements, values of $\omega$ and $\Omega$ that are separated by 180 deg are degenerate, as discussed in Blunt et al. (2020). In order to account for this, we used uniform priors between 0 and 180 deg for both of these parameters. All priors are listed in Table 3.

The sampler was run using 40 temperatures (using an exponential ladder, with each temperature increasing by a factor of $\sqrt{2}$, so the highest temperature is $\sqrt{2}^{40}$), 1000 walkers per temperature and 50 million steps per walker. Convergence was assessed by visual inspection of the chains. Due to the ambiguity in the two position angles ($\pm 180$ deg ambiguity) derived from the speckle observations obtained in 2011 and 2022, we ran the ORBITIZE! model four times to account for each possible combination of position angles. Convergence was not reached after 50 million runs when either of the two
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Spectral energy distribution of TIC 307958020. Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the passband. Blue symbols are the model fluxes from the best-fit Kurucz atmosphere model (gray) for the combined light of the binary star system. Cyan and red curves represent the Kurucz atmosphere models corresponding to the warmer (larger) and cooler (smaller) stellar components of the system, respectively.</figcaption></figure>

[Fe/H] are in excellent agreement between the two methods within their respective uncertainties, the determined values of the stellar radius are discrepant ($R_{\text{SpecMatch-Syn}} = 1.50 \pm 0.04 R_⊙$ and $R_{\text{SpecMatch-Emp}} = 1.11 \pm 0.18 R_⊙$). In addition to the discrepancy between these two methods, both are affected by the assumption that the light originates from a single star. We, therefore, independently derived the stellar parameters using the spectral energy distribution (SED), as this is able to account for the light contributions of both stars in the system.

### 5.1.2. SED Fitting

As an independent determination of the stellar parameters, we performed an analysis of the broadband SED of the star together with the Gaia DR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016); Stassun et al. (2017, 2018). We pulled the $B_T V_T$ magnitudes from Tycho-2, the JHK$_S$ magnitudes from 2MASS, the W1–W4 magnitudes from WISE, the G, $G_{BP}$, and $G_{RP}$ magnitudes from Gaia, and the FUV and NUV magnitudes from GALEX. Together, the available photometry spans the full stellar SED over the wavelength range 0.2–22 µm (see Figure 6).

For an initial fit, we treated the SED as arising from a single star, with the fitted parameters being the effective temperature ($T_{\text{eff}}$), surface gravity (log g), and metallicity ([Fe/H]), as well as the extinction $A_V$, which we limited to the maximum line-of-sight value from the Galactic dust maps of Schlegel et al. (1998). For the initial single-star model fit, we adopted the $T_{\text{eff}}$, log g, and [Fe/H] from the spectroscopically determined values. The goodness of fit was, not surprisingly, only marginally good ($\chi^2 = 4.5$) and, given the stringent Gaia distance, implied an oversized star with $R_* \approx 1.5 R_⊙$.

Next, we performed a two-component fit using the methodology of Stassun & Torres (2016), solving for the component $T_{\text{eff}}$ and $R_*$ by requiring the flux-weighted average $T_{\text{eff}}$ to agree with the spectroscopic value used above while also requiring the flux ratio $F_B/F_A$ in the visible to match that determined from the high-contrast imaging, and of course requiring agreement with the strict Gaia distance.

The result, shown in Figure 6, has $\chi^2 = 1.1$ and best-fit parameters of $A_V = 0.07 \pm 0.02$, $T_{\text{eff,A}} = 5800 \pm 50$ K, $R_A = 1.05 \pm 0.05 R_⊙$, $T_{\text{eff,B}} = 5600 \pm 50$ K, and $R_B = 0.98 \pm 0.05 R_⊙$. Based on the empirical relations of Torres et al. (2010), the hotter/larger star appears to have a mass $M_A = 1.10 \pm 0.06 M_⊙$, and the cooler/smaller star has $M_B = 1.05 \pm 0.06 M_⊙$, consistent with an interpretation of main-sequence stars.

Finally, the vsini measurement reported by APOGEE DR16 ($\text{vsini} = 4.39 \text{ km s}^{-1}$) together with the brighter star's radius implies a maximum rotation period of $12.1 \pm 1.3$ days. With the gyrochronology relations of Mamajek & Hillenbrand (2008), this implies a minimum age of $1.3 \pm 0.3$ Gyr, again consistent with unevolved main-sequence stars. As a comparison, the vsini derived using the SpecMatch software package from the HIRES data is $4.25 \text{ km s}^{-1}$, which is also consistent with the system containing main-sequence stars. We note that the vsini estimates are likely strongly affected by the contamination of the companion star. Finally, as mentioned in Section 4.1, the TESS and ASAS-SN data show no signs of measurable rotation and thus a more precise age of this system cannot be determined using gyrochronology.

As the two-component SED analysis is able to account for the light contribution of both stars, we adopt these stellar parameters for the remainder of the paper. All parameters are listed in Table 3.

## 5.2. Binary Orbit Modeling

The orbital parameters of the binary system (TOI 4633 AB) were determined by Bayesian parameter estimation using the open-source software "ORBITIZE!" (Blunt et al. 2020). The position angles and angular separations extracted from all of the available imaging data (listed in Table 2) were used as input data. Observations obtained prior to 2011 did not report uncertainties on their measurements. We, therefore, adopted large uncertainties for these archival position angles and angular separations of $\pm 10$ deg and $\pm 50$ mas, respectively. As discussed in Section 3.1, the HIRES data allows us to place an upper limit on the relative motion of the two stars of $10 \text{ km s}^{-1}$. As such, RV values of $0 \pm 5 \text{ km s}^{-1}$ at the times of the HIRES observations are used as input RV measurements for the ORBITIZE! model.

We used the parallel-tempered Affine-invariant sampler PTEMCEE (Foreman-Mackey et al. 2013; Vousden et al. 2016) and adopted priors of $10.55 \pm 0.013$ mas (Gaia Collaboration et al. 2021) for the parallax and $1.10 \pm 0.2 M_⊙$ and $1.05 \pm 0.2 M_⊙$ for the masses of the primary and secondary stars, respectively (see Section 5.1.2). Due to the lack of more constraining RV measurements, values of $\omega$ and $\Omega$ that are separated by 180 deg are degenerate, as discussed in Blunt et al. (2020). In order to account for this, we used uniform priors between 0 and 180 deg for both of these parameters. All priors are listed in Table 3.

The sampler was run using 40 temperatures (using an exponential ladder, with each temperature increasing by a factor of $\sqrt{2}$, so the highest temperature is $\sqrt{2}^{40}$), 1000 walkers per temperature and 50 million steps per walker. Convergence was assessed by visual inspection of the chains. Due to the ambiguity in the two position angles ($\pm 180$ deg ambiguity) derived from the speckle observations obtained in 2011 and 2022, we ran the ORBITIZE! model four times to account for each possible combination of position angles. Convergence was not reached after 50 million runs when either of the two
<table><caption>**Table 3**
Stellar System Parameters</caption><thead><tr><th>Identifiers</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td>HU</td><td>918</td><td>Hussey (1905)</td></tr><tr><td>TOI</td><td>4633</td><td>...</td></tr><tr><td>TIC</td><td>307958020</td><td>Stassun et al. (2019)</td></tr><tr><td>Gaia DR3</td><td>1630906044157332224</td><td>Gaia DR3<sup>a</sup></td></tr><tr><td>2MASS</td><td>J17072238 + 6228330</td><td>2MASS<sup>b</sup></td></tr><tr><td>Astrometry</td><td>Value</td><td>Source</td></tr><tr><td>α<sub>J2000</sub></td><td>17:07:22.396</td><td>Gaia eDR3<sup>a</sup></td></tr><tr><td>δ<sub>J2000</sub></td><td>62:28:33.011</td><td>Gaia eDR3<sup>a</sup></td></tr><tr><td>Distance (pc)</td><td>95.20 ± 0.24</td><td>Bailer-Jones et al. (2018)</td></tr><tr><td>π (mas)</td><td>10.552 ± 0.014</td><td>Gaia eDR3<sup>a</sup></td></tr><tr><td>Spectral Type</td><td>early G Dwarf</td><td>...</td></tr><tr><td>Photometry</td><td>Magnitude</td><td>Source</td></tr><tr><td>B</td><td>9.767 ± 0.033</td><td>Tycho-2<sup>c</sup></td></tr><tr><td>V</td><td>9.017 ± 0.002</td><td>Tycho-2<sup>c</sup></td></tr><tr><td>J</td><td>7.723 ± 0.030</td><td>2MASS<sup>c</sup></td></tr><tr><td>H</td><td>7.448 ± 0.044</td><td>2MASS<sup>c</sup></td></tr><tr><td>K</td><td>7.349 ± 0.024</td><td>2MASS<sup>c</sup></td></tr><tr><td>W1</td><td>7.230 ± 0.039</td><td>WISE<sup>c</sup></td></tr><tr><td>W2</td><td>7.326 ± 0.020</td><td>WISE<sup>c</sup></td></tr><tr><td>W3</td><td>7.305 ± 0.016</td><td>WISE<sup>d</sup></td></tr><tr><td>SED derived properties</td><td>Star A</td><td>Star B</td></tr><tr><td>Effective temperature T<sub>eff</sub><br>(K)</td><td>5800 ± 50</td><td>5600 ± 50</td></tr><tr><td>Stellar mass M<sub>☉</sub> (M<sub>☉</sub>)</td><td>1.10 ± 0.06</td><td>1.05 ± 0.06</td></tr><tr><td>Stellar radius R<sub>☉</sub> (R<sub>☉</sub>)</td><td>1.05 ± 0.05</td><td>0.98 ± 0.05</td></tr><tr><td>ORBITIZE! binary<br>parameters</td><td>Prior<sup>e</sup></td><td>Derived value</td></tr><tr><td>Semimajor axis a<sub>bin</sub> (au)</td><td>U[0, 1000]</td><td>48.6<sup>+4.4</sup><sub>-3.5</sub></td></tr><tr><td>Eccentricity e<sub>bin</sub></td><td>U[0, 1]</td><td>0.91<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>Inclination i<sub>bin</sub> (deg)</td><td>U[0, 180]</td><td>90.1<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>ω (deg)</td><td>U[0, 180]</td><td>110.5<sup>+2.1</sup><sub>-2.1</sub></td></tr><tr><td>Ω (deg)</td><td>U[0, 180]</td><td>123.5<sup>+3.3</sup><sub>-2.9</sub></td></tr><tr><td>Phase</td><td>U[0, 1]</td><td>0.42<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>Period (years)</td><td>derived</td><td>231<sup>+32</sup><sub>-24</sub></td></tr><tr><td>Periastron (au)</td><td>derived</td><td>4.5<sup>+2.1</sup><sub>-1.5</sub></td></tr></tbody></table>

**Notes.**

position angles obtained from speckle observations were rotated by 180 deg from the values listed in Table 2. As such, we will assume for the remainder of this paper that the position angles listed in this table are the most likely to be correct.

We found the binary semimajor axis, eccentricity and inclination to be constrained to $a_{\text{bin}} = 48.6_{-3.5}^{+4.4}$ au, $e_{\text{bin}} = 0.91_{-0.03}^{+0.03}$, and $i_{\text{bin}} = 90.1_{-0.4}^{+0.4}$ deg (see Figure 7 for the posterior distributions of these three parameters). This corresponds to a stellar orbital period of $231_{-24}^{+32}$ yr. Figure 8 shows 100 model fits to the position angles and angular separation that were randomly sampled from the posteriors. The derived RV model indicates a semi-amplitude of $7.2 \text{ km s}^{-1}$, with a predicted RV shift of less

than $0.2 \text{ km s}^{-1}$ during the ground-based RV observational time base.

## 6. Analysis of Planets

In this Section we discuss the properties of the two observed planetary signals, where we make the assumption that both signals originate from planets that are orbiting the same star. Given the current small angular separation of the two stars, we are unable to determine around which of the two stars either of the planets are orbiting. For the remainder of the paper we report all planet properties under the assumption that star A is the host star for both planets. We note that due to the similarities and large uncertainties in the derived properties of the two stars ($M_A = 1.10 \pm 0.06 \text{ M}_\odot$, $R_A = 1.05 \pm 0.05 \text{ R}_\odot$; $M_B = 1.05 \pm 0.06 \text{ M}_\odot$, $R_B = 0.98 \pm 0.05 \text{ R}_\odot$), the derived planet properties agree to within their uncertainties when derived with either of the two stars as the host.

### 6.1. Joint Transit and RV Modeling

The transit and RV data were jointly modeled using *pyaneti* (Barragán et al. 2019, 2022a). This open-source software creates marginalized posterior distributions for different orbital parameters by sampling the parameter space using an MCMC approach. The transits are modeled using the limb darkened quadratic models by Mandel & Agol (2002) while the RV data are fit with two Keplerian RV models. In order to account for an RV offset between the HIRES and SOPHIE data, we allow for a systematic velocity for each instrument and include a jitter term per instrument to account for imperfections in our transit and RV model. For each sector, we used either the 2 minutes cadence or, where available, the 20 s cadence data for the transit model (see Section 2) and both the HIRES and SOPHIE data for the RV fit. All fitted parameters and priors used for the joint modeling are presented in Table 4.

The parameter space was sampled using an MCMC approach with 250 individual chains and posterior distributions were generated using 5000 iterations of converged chains with a thin factor of 10. The inferred parameters extracted from the posteriors are shown in Table 4 while the inferred models are shown in Figures 9 and 10. We note that we recover relatively large jitter terms in the RV modeling of $12.36_{-1.20}^{+1.36} \text{ m s}^{-1}$ and $13.40_{-2.39}^{+3.10} \text{ m s}^{-1}$ for HIRES and SOPHIE, respectively. This can be caused by a sub-optimal RV model or unknown systematics affecting the RV data. We suspect that the contamination from the second star also contributes to the observed scatter. We see no evidence of a long-term RV trend.

We checked for additional periodic signals in the data using the PyAstronomy Generalized Lomb-Scargle Periodogram on the residual RV data. We find a significant peak at a period of 182.8 days with a false alarm probability of ~1.5%. However, as this is half of the orbital period of the Earth, we consider this to most likely be a systematic effect and do not discuss this signal any further.

### 6.1.1. Effects of Companion Star on Planet Radius and Mass

Planet radii are calculated based on the observed transit depth and the measured radius of the host star. However, for multi-stellar systems, the relationship of the observed transit depth and the true planet radius depends on the brightness ratio of the brightness of the star which is being transited to the total brightness of all the stars in the system (Furlan et al. 2017).

<footnote><sup>a</sup> Gaia early Data Release 3 (eDR3; Gaia Collaboration et al. 2021).</footnote>

<footnote><sup>b</sup> Two-micron All Sky Survey (2MASS; Cutri et al. 2003).</footnote>

<footnote><sup>c</sup> Tycho-2 catalog (Hög et al. 2000).</footnote>

<footnote><sup>d</sup> Wide-field Infrared Survey Explorer catalog (WISE; Cutri et al. 2013).</footnote>

<footnote><sup>e</sup> U[a, b] refers to uniform priors between a and b. The uncertainties in the derived values represent the 68% confidence interval.</footnote>
<table><caption>**Table 3**
Stellar System Parameters</caption><thead><tr><th>Identifiers</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td>HU</td><td>918</td><td>Hussey (1905)</td></tr><tr><td>TOI</td><td>4633</td><td>...</td></tr><tr><td>TIC</td><td>307958020</td><td>Stassun et al. (2019)</td></tr><tr><td>Gaia DR3</td><td>1630906044157332224</td><td>Gaia DR3<sup>a</sup></td></tr><tr><td>2MASS</td><td>J17072238 + 6228330</td><td>2MASS<sup>b</sup></td></tr><tr><td>Astrometry</td><td>Value</td><td>Source</td></tr><tr><td>α<sub>J2000</sub></td><td>17:07:22.396</td><td>Gaia eDR3<sup>a</sup></td></tr><tr><td>δ<sub>J2000</sub></td><td>62:28:33.011</td><td>Gaia eDR3<sup>a</sup></td></tr><tr><td>Distance (pc)</td><td>95.20 ± 0.24</td><td>Bailer-Jones et al. (2018)</td></tr><tr><td>π (mas)</td><td>10.552 ± 0.014</td><td>Gaia eDR3<sup>a</sup></td></tr><tr><td>Spectral Type</td><td>early G Dwarf</td><td>...</td></tr><tr><td>Photometry</td><td>Magnitude</td><td>Source</td></tr><tr><td>B</td><td>9.767 ± 0.033</td><td>Tycho-2<sup>c</sup></td></tr><tr><td>V</td><td>9.017 ± 0.002</td><td>Tycho-2<sup>c</sup></td></tr><tr><td>J</td><td>7.723 ± 0.030</td><td>2MASS<sup>c</sup></td></tr><tr><td>H</td><td>7.448 ± 0.044</td><td>2MASS<sup>c</sup></td></tr><tr><td>K</td><td>7.349 ± 0.024</td><td>2MASS<sup>c</sup></td></tr><tr><td>W1</td><td>7.230 ± 0.039</td><td>WISE<sup>c</sup></td></tr><tr><td>W2</td><td>7.326 ± 0.020</td><td>WISE<sup>c</sup></td></tr><tr><td>W3</td><td>7.305 ± 0.016</td><td>WISE<sup>d</sup></td></tr><tr><td>SED derived properties</td><td>Star A</td><td>Star B</td></tr><tr><td>Effective temperature T<sub>eff</sub><br>(K)</td><td>5800 ± 50</td><td>5600 ± 50</td></tr><tr><td>Stellar mass M<sub>☉</sub> (M<sub>☉</sub>)</td><td>1.10 ± 0.06</td><td>1.05 ± 0.06</td></tr><tr><td>Stellar radius R<sub>☉</sub> (R<sub>☉</sub>)</td><td>1.05 ± 0.05</td><td>0.98 ± 0.05</td></tr><tr><td>ORBITIZE! binary<br>parameters</td><td>Prior<sup>e</sup></td><td>Derived value</td></tr><tr><td>Semimajor axis a<sub>bin</sub> (au)</td><td>U[0, 1000]</td><td>48.6<sup>+4.4</sup><sub>-3.5</sub></td></tr><tr><td>Eccentricity e<sub>bin</sub></td><td>U[0, 1]</td><td>0.91<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>Inclination i<sub>bin</sub> (deg)</td><td>U[0, 180]</td><td>90.1<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>ω (deg)</td><td>U[0, 180]</td><td>110.5<sup>+2.1</sup><sub>-2.1</sub></td></tr><tr><td>Ω (deg)</td><td>U[0, 180]</td><td>123.5<sup>+3.3</sup><sub>-2.9</sub></td></tr><tr><td>Phase</td><td>U[0, 1]</td><td>0.42<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>Period (years)</td><td>derived</td><td>231<sup>+32</sup><sub>-24</sub></td></tr><tr><td>Periastron (au)</td><td>derived</td><td>4.5<sup>+2.1</sup><sub>-1.5</sub></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Gaia early Data Release 3 (eDR3; Gaia Collaboration et al. 2021).

<sup>b</sup> Two-micron All Sky Survey (2MASS; Cutri et al. 2003).

<sup>c</sup> Tycho-2 catalog (Hög et al. 2000).

<sup>d</sup> Wide-field Infrared Survey Explorer catalog (WISE; Cutri et al. 2013).

<sup>e</sup> U[a, b] refers to uniform priors between a and b. The uncertainties in the derived values represent the 68% confidence interval.

position angles obtained from speckle observations were rotated by 180 deg from the values listed in Table 2. As such, we will assume for the remainder of this paper that the position angles listed in this table are the most likely to be correct.

We found the binary semimajor axis, eccentricity and inclination to be constrained to $a_{\text{bin}} = 48.6_{-3.5}^{+4.4}$ au, $e_{\text{bin}} = 0.91_{-0.03}^{+0.03}$, and $i_{\text{bin}} = 90.1_{-0.4}^{+0.4}$ deg (see Figure 7 for the posterior distributions of these three parameters). This corresponds to a stellar orbital period of $231_{-24}^{+32}$ yr. Figure 8 shows 100 model fits to the position angles and angular separation that were randomly sampled from the posteriors. The derived RV model indicates a semi-amplitude of $7.2 \text{ km s}^{-1}$, with a predicted RV shift of less

than $0.2 \text{ km s}^{-1}$ during the ground-based RV observational time base.

## 6. Analysis of Planets

In this Section we discuss the properties of the two observed planetary signals, where we make the assumption that both signals originate from planets that are orbiting the same star. Given the current small angular separation of the two stars, we are unable to determine around which of the two stars either of the planets are orbiting. For the remainder of the paper we report all planet properties under the assumption that star A is the host star for both planets. We note that due to the similarities and large uncertainties in the derived properties of the two stars ($M_A = 1.10 \pm 0.06 \text{ M}_\odot$, $R_A = 1.05 \pm 0.05 \text{ R}_\odot$; $M_B = 1.05 \pm 0.06 \text{ M}_\odot$, $R_B = 0.98 \pm 0.05 \text{ R}_\odot$), the derived planet properties agree to within their uncertainties when derived with either of the two stars as the host.

### 6.1. Joint Transit and RV Modeling

The transit and RV data were jointly modeled using *pyaneti* (Barragán et al. 2019, 2022a). This open-source software creates marginalized posterior distributions for different orbital parameters by sampling the parameter space using an MCMC approach. The transits are modeled using the limb darkened quadratic models by Mandel & Agol (2002) while the RV data are fit with two Keplerian RV models. In order to account for an RV offset between the HIRES and SOPHIE data, we allow for a systematic velocity for each instrument and include a jitter term per instrument to account for imperfections in our transit and RV model. For each sector, we used either the 2 minutes cadence or, where available, the 20 s cadence data for the transit model (see Section 2) and both the HIRES and SOPHIE data for the RV fit. All fitted parameters and priors used for the joint modeling are presented in Table 4.

The parameter space was sampled using an MCMC approach with 250 individual chains and posterior distributions were generated using 5000 iterations of converged chains with a thin factor of 10. The inferred parameters extracted from the posteriors are shown in Table 4 while the inferred models are shown in Figures 9 and 10. We note that we recover relatively large jitter terms in the RV modeling of $12.36_{-1.20}^{+1.36} \text{ m s}^{-1}$ and $13.40_{-2.39}^{+3.10} \text{ m s}^{-1}$ for HIRES and SOPHIE, respectively. This can be caused by a sub-optimal RV model or unknown systematics affecting the RV data. We suspect that the contamination from the second star also contributes to the observed scatter. We see no evidence of a long-term RV trend.

We checked for additional periodic signals in the data using the PyAstronomy Generalized Lomb-Scargle Periodogram on the residual RV data. We find a significant peak at a period of 182.8 days with a false alarm probability of ~1.5%. However, as this is half of the orbital period of the Earth, we consider this to most likely be a systematic effect and do not discuss this signal any further.

### 6.1.1. Effects of Companion Star on Planet Radius and Mass

Planet radii are calculated based on the observed transit depth and the measured radius of the host star. However, for multi-stellar systems, the relationship of the observed transit depth and the true planet radius depends on the brightness ratio of the brightness of the star which is being transited to the total brightness of all the stars in the system (Furlan et al. 2017).
<figure><img src="image_9.png" /><figcaption>Figure 7. Corner plot of the posteriors of the ORBITIZE! results.</figcaption></figure>

Assuming that the planet transits star A and given $F_B/F_A = 0.7$, we correct the observed planet radius of TOI 4633 c ($R_{c,obs} = 2.42_{-0.14}^{+0.15} R_⊕$) by a factor of $\sqrt{\frac{F_A + F_B}{F_A}} = \sqrt{1.7}$, resulting in a true planet radius of $3.2_{-0.19}^{+0.20} R_⊕$.

Similarly, the presence of the second star acts to reduce the observed amplitude of the RV shifts. The amount by which the observed signal is reduced depends on the relative strengths of the absorption lines (line depths) of both stars, as well as on their rotational velocities (line widths). As we are not able to disentangle the two spectra we cannot quantify the effect of the second star on the RV amplitude. As such, the observed RV amplitudes, and therefore the derived planet masses, are lower limits. However, given the radius of planet c and the known distribution of planet masses for a given planetary radius (as shown in Figure 11), we consider it unlikely for planet c to have a mass significantly greater than $47 M_⊕$.

## 6.2. Search and Recovery of Transit Signals

In order to search for transits of planet candidate b using the TESS data, and to search for additional transit signals, we searched the full TESS light curve using the Box Least Squares (BLS; Kovács et al. 2002) algorithm. Before running the BLS search, we masked the transit signals of TOI 4633 c and used an iterative nonlinear filter to subtract residual systematics on timescales >1.7 days (Aigrain & Irwin 2004). We carried out the BLS search on an evenly sampled frequency grid ranging from 0.00125 to $1\,\text{d}^{-1}$ (1–800 days). The signal detection efficiency (SDE), defined as the ratio of the highest peak in the SNR periodogram relative to its standard deviation, was used to

determine the significance of the recovery of the signal. The algorithms found no additional signals above an SDE of 7.6.

Furthermore, we used an injection and recovery test to quantify the detectability of additional planets in the TESS data, following the methodology outlined in Eisner et al. (2020b). In brief we injected transit signals generated using the BATMAN package (Kreidberg 2015) into the PDC TESS light curve. The injected transit signals corresponded to planets with radii ranging from 1 to 12.5 $R_⊕$ and periods ranging from 1 to 300 days, both sampled randomly from a log-uniform distribution. The impact parameter and eccentricity were assumed to be zero for simplicity. We used the stellar parameters given in Table 3 and adopted a quadratic limb-darkening law with $q_1$ and $q_2$ of 0.16 and 0.59, respectively, as taken from Table 15 in Claret (2017).

We simulated and injected transits for 500,000 planets and used the BLS methodology described above to try to recover each injected signal. For each simulation we identified the highest peak in the BLS periodogram. The signal was considered to be correctly identified when the corresponding period and orbital phase were within 1% of the injected values. The fraction of recovered signals over a grid of period and radius bins was used to evaluate the completeness of the injection and recovery search. The radius and period bins have widths of $0.2 R_⊕$ and 10 days respectively, as shown in Figure 12. Using the planet mass–radius scaling relations for volatile-rich planets by Otegi et al. (2020) and the measured minimum planet mass of $109 M_⊕$, planet candidate b has an estimated radius of ~13.7 $R_⊕$ corresponding to a minimum density of $0.11\,\text{g cm}^{-3}$ (assuming $b=0$). Given the light
<figure><img src="image_11.png" /><figcaption>Figure 7. Corner plot of the posteriors of the ORBITIZE! results.</figcaption></figure>

Assuming that the planet transits star A and given $F_B/F_A = 0.7$, we correct the observed planet radius of TOI 4633 c ($R_{c,obs} = 2.42_{-0.14}^{+0.15} R_⊕$) by a factor of $\sqrt{\frac{F_A + F_B}{F_A}} = \sqrt{1.7}$, resulting in a true planet radius of $3.2_{-0.19}^{+0.20} R_⊕$.

Similarly, the presence of the second star acts to reduce the observed amplitude of the RV shifts. The amount by which the observed signal is reduced depends on the relative strengths of the absorption lines (line depths) of both stars, as well as on their rotational velocities (line widths). As we are not able to disentangle the two spectra we cannot quantify the effect of the second star on the RV amplitude. As such, the observed RV amplitudes, and therefore the derived planet masses, are lower limits. However, given the radius of planet c and the known distribution of planet masses for a given planetary radius (as shown in Figure 11), we consider it unlikely for planet c to have a mass significantly greater than $47 M_⊕$.

## 6.2. Search and Recovery of Transit Signals

In order to search for transits of planet candidate b using the TESS data, and to search for additional transit signals, we searched the full TESS light curve using the Box Least Squares (BLS; Kovács et al. 2002) algorithm. Before running the BLS search, we masked the transit signals of TOI 4633 c and used an iterative nonlinear filter to subtract residual systematics on timescales >1.7 days (Aigrain & Irwin 2004). We carried out the BLS search on an evenly sampled frequency grid ranging from 0.00125 to $1\,\text{d}^{-1}$ (1–800 days). The signal detection efficiency (SDE), defined as the ratio of the highest peak in the SNR periodogram relative to its standard deviation, was used to

determine the significance of the recovery of the signal. The algorithms found no additional signals above an SDE of 7.6.

Furthermore, we used an injection and recovery test to quantify the detectability of additional planets in the TESS data, following the methodology outlined in Eisner et al. (2020b). In brief we injected transit signals generated using the BATMAN package (Kreidberg 2015) into the PDC TESS light curve. The injected transit signals corresponded to planets with radii ranging from 1 to 12.5 $R_⊕$ and periods ranging from 1 to 300 days, both sampled randomly from a log-uniform distribution. The impact parameter and eccentricity were assumed to be zero for simplicity. We used the stellar parameters given in Table 3 and adopted a quadratic limb-darkening law with $q_1$ and $q_2$ of 0.16 and 0.59, respectively, as taken from Table 15 in Claret (2017).

We simulated and injected transits for 500,000 planets and used the BLS methodology described above to try to recover each injected signal. For each simulation we identified the highest peak in the BLS periodogram. The signal was considered to be correctly identified when the corresponding period and orbital phase were within 1% of the injected values. The fraction of recovered signals over a grid of period and radius bins was used to evaluate the completeness of the injection and recovery search. The radius and period bins have widths of $0.2 R_⊕$ and 10 days respectively, as shown in Figure 12. Using the planet mass–radius scaling relations for volatile-rich planets by Otegi et al. (2020) and the measured minimum planet mass of $109 M_⊕$, planet candidate b has an estimated radius of ~13.7 $R_⊕$ corresponding to a minimum density of $0.11\,\text{g cm}^{-3}$ (assuming $b=0$). Given the light
<figure><img src="image_10.png" /><figcaption>**Figure 8.** 100 orbits randomly drawn from the posterior distribution of the ORBITIZE! model (gray lines) for the angular separations (top panel) and position angle (bottom panel) of TOI 4633 AB. The best-fit models constructed from the median of the posteriors of each parameter are shown in black, and the residuals to these models are shown below each panel. Archival astrometric data, which have no uncertainties reported in the literature, are shown by the black points. We adopted large uncertainties for these archival position angles and angular separations of ± 10 deg and ± 50 mas, respectively. All newly obtained data with reported uncertainties are shown in orange.</figcaption></figure>

contribution of the companion star which would act to dilute
the transit, this would result in an observed radius of ~10.7 R<sub>⊕</sub>
(indicated by the black cross). Figure 12 shows that the
properties of planet candidate *b* are such that if the planet were
transiting we would be able to identify the transit events in the
TESS data. This, in turn, allows us to constrain the inclination
of planet candidate *b* to *i* < 88.5 deg.

6.3. Dynamical Stability

The orbital periods of planet candidate b and planet c are
much shorter than the period of the binary orbit, indicating that
this is a circumstellar, or “s-type” system. Furthermore, the
large period ratio between the planets’ orbits and the binary
orbits indicate that this system is in a hierarchical configuration.
Additionally, the significant eccentricity of the binary orbit and
that the mutual inclination between the planets and the binary
orbit is not known indicates that the stability of this system can
be analyzed in light of the Eccentric Kozai–Lidov mechanism
(see e.g., Kozai 1962; Lidov 1962; Naoz 2016), where the
system is unstable if the mutual inclination between the binary
and the planets lies between the “Kozai Angles,” that is
39°2 ≤ *i*mut ≤ 140°77. Within these angles, the planets will

undergo large oscillations in their eccentricities and inclina-
tions. We explore the stability of orbital configurations with
low mutual inclinations using the framework of Quarles et al.
(2020), who numerically constrain the stability regime of
Earth-mass planets in circumstellar orbits in binary systems
with approximately 700 million *N*-body simulations. The grid
map of Quarles et al. (2020) indicates that, in the prograde orbit
case, without TOI 4633 *b*, the transiting planet TOI 4633 *c*
remains stable against perturbation by the binary orbit if its
mutual inclination with the binary orbit (*i*<sub>m,bin-c</sub>) is between 0°
and 45°.

However, as this is likely a two-planet system, planet–planet
interactions should be considered. It has been shown that
planet–planet interactions (and general-relativistic precession
of the periapsis) can suppress eccentricity oscillations and
destabilization from the influence of an outer perturber (see
e.g., Naoz 2016; Denham et al. 2019; Wei et al. 2021; Faridani
et al. 2022). In this system, however, the period ratio between
the two planets is too large for planet–planet interactions to
meaningfully stabilize TOI 4633 c against the perturbations
from the binary if the binary and the planets have meaningful
mutual inclination (as calculated using Equation (10) from
<figure><img src="image_12.png" /><figcaption>**Figure 8.** 100 orbits randomly drawn from the posterior distribution of the ORBITIZE! model (gray lines) for the angular separations (top panel) and position angle (bottom panel) of TOI 4633 AB. The best-fit models constructed from the median of the posteriors of each parameter are shown in black, and the residuals to these models are shown below each panel. Archival astrometric data, which have no uncertainties reported in the literature, are shown by the black points. We adopted large uncertainties for these archival position angles and angular separations of ± 10 deg and ± 50 mas, respectively. All newly obtained data with reported uncertainties are shown in orange.</figcaption></figure>

contribution of the companion star which would act to dilute
the transit, this would result in an observed radius of ~10.7 R<sub>⊕</sub>
(indicated by the black cross). Figure 12 shows that the
properties of planet candidate *b* are such that if the planet were
transiting we would be able to identify the transit events in the
TESS data. This, in turn, allows us to constrain the inclination
of planet candidate *b* to *i* < 88.5 deg.

6.3. Dynamical Stability

The orbital periods of planet candidate b and planet c are
much shorter than the period of the binary orbit, indicating that
this is a circumstellar, or “s-type” system. Furthermore, the
large period ratio between the planets’ orbits and the binary
orbits indicate that this system is in a hierarchical configuration.
Additionally, the significant eccentricity of the binary orbit and
that the mutual inclination between the planets and the binary
orbit is not known indicates that the stability of this system can
be analyzed in light of the Eccentric Kozai–Lidov mechanism
(see e.g., Kozai 1962; Lidov 1962; Naoz 2016), where the
system is unstable if the mutual inclination between the binary
and the planets lies between the “Kozai Angles,” that is
39°2 ≤ *i*mut ≤ 140°77. Within these angles, the planets will

undergo large oscillations in their eccentricities and inclina-
tions. We explore the stability of orbital configurations with
low mutual inclinations using the framework of Quarles et al.
(2020), who numerically constrain the stability regime of
Earth-mass planets in circumstellar orbits in binary systems
with approximately 700 million *N*-body simulations. The grid
map of Quarles et al. (2020) indicates that, in the prograde orbit
case, without TOI 4633 *b*, the transiting planet TOI 4633 *c*
remains stable against perturbation by the binary orbit if its
mutual inclination with the binary orbit (*i*<sub>m,bin-c</sub>) is between 0°
and 45°.

However, as this is likely a two-planet system, planet–planet
interactions should be considered. It has been shown that
planet–planet interactions (and general-relativistic precession
of the periapsis) can suppress eccentricity oscillations and
destabilization from the influence of an outer perturber (see
e.g., Naoz 2016; Denham et al. 2019; Wei et al. 2021; Faridani
et al. 2022). In this system, however, the period ratio between
the two planets is too large for planet–planet interactions to
meaningfully stabilize TOI 4633 c against the perturbations
from the binary if the binary and the planets have meaningful
mutual inclination (as calculated using Equation (10) from
<table><caption>**Table 4**
System Parameters</caption><thead><tr><th>Parameter</th><th>Prior<sup>a</sup></th><th>Value<sup>b</sup></th><th>Comments</th></tr></thead><tbody><tr><td colspan="4"><strong>Model Parameters for TOI 4633 b</strong></td></tr><tr><td>Orbital period P<sub>orb</sub> (days)</td><td>U[33, 36]</td><td>34.15 ± 0.15</td><td></td></tr><tr><td>Time of min. conjunction T<sub>0</sub> (BJD—2457000)</td><td>U[2787.32, 2807.32]</td><td>2796.64<sup>+1.10</sup><sub>-1.19</sub></td><td></td></tr><tr><td>Parameterization √e sin ω</td><td>U[-1, 1]</td><td>-0.12<sup>+0.26</sup><sub>-0.22</sub></td><td>The code ensures e &lt; 1</td></tr><tr><td>Parameterization √e cos ω</td><td>U[-1, 1]</td><td>0.07<sup>+0.22</sup><sub>-0.24</sub></td><td>The code ensures e &lt; 1</td></tr><tr><td>Doppler semi-amplitude K (m s<sup>-1</sup>)<sup>c</sup></td><td>U[0, 30]</td><td>19.97<sup>+2.29</sup><sub>-2.30</sub></td><td></td></tr><tr><td colspan="4"><strong>Model Parameters for TOI 4633 c</strong></td></tr><tr><td>Orbital period P<sub>orb</sub> (days)</td><td>U[271.9, 272.1]</td><td>271.9445<sup>+0.0039</sup><sub>-0.0040</sub></td><td></td></tr><tr><td>Transit epoch T<sub>0</sub> (BJD—2457000)</td><td>U[1863.81, 1865.81]</td><td>1864.8265<sup>+0.0088</sup><sub>-0.0083</sub></td><td></td></tr><tr><td>Parameterization √e sin ω</td><td>U[-1, 1]</td><td>-0.06<sup>+0.23</sup><sub>-0.27</sub></td><td>The code ensures e &lt; 1</td></tr><tr><td>Parameterization √e cos ω</td><td>U[-1, 1]</td><td>0.04<sup>+0.31</sup><sub>-0.35</sub></td><td>The code ensures e &lt; 1</td></tr><tr><td>Observed scaled planet radius R<sub>p</sub>/R<sub>⋆</sub></td><td>U[0, 0.1]</td><td>0.0211<sup>+0.00076</sup><sub>-0.00069</sub></td><td></td></tr><tr><td>Impact parameter b</td><td>U[0, 1.1]</td><td>0.33<sup>+0.24</sup><sub>-0.21</sub></td><td></td></tr><tr><td>Doppler semi-amplitude K (m s<sup>-1</sup>)<sup>c</sup></td><td>U[0, 30]</td><td>4.58<sup>+2.56</sup><sub>-2.29</sub></td><td>11.5 m s<sup>-1</sup>, 99% percent upper limit</td></tr><tr><td colspan="4"><strong>Other Parameters</strong></td></tr><tr><td>Stellar density ρ<sub>⋆</sub> (g cm<sup>-3</sup>)</td><td>N[1.34, 0.23]</td><td>1.35 ± 0.21</td><td></td></tr><tr><td>Parameterized limb-darkening coefficient q<sub>1</sub></td><td>U[0, 1]</td><td>0.39<sup>+0.37</sup><sub>-0.23</sub></td><td>q<sub>1</sub> parameter as in Kipping (2013)</td></tr><tr><td>Parameterized limb-darkening coefficient q<sub>2</sub></td><td>U[0, 1]</td><td>0.34<sup>+0.31</sup><sub>-0.22</sub></td><td>q<sub>2</sub> parameter as in Kipping (2013)</td></tr><tr><td>Offset velocity HIRES (km s<sup>-1</sup>)</td><td>U[-0.50, 0.50]</td><td>0.0023<sup>+0.0019</sup><sub>-0.0020</sub></td><td></td></tr><tr><td>Offset velocity SOPHIE (km s<sup>-1</sup>)</td><td>U[-0.50, 0.50]</td><td>-0.0007 ± 0.0034</td><td></td></tr><tr><td>Jitter HIRES (m s<sup>-1</sup>)</td><td>J[1, 100]</td><td>12.36<sup>+1.36</sup><sub>-1.20</sub></td><td></td></tr><tr><td>Jitter SOPHIE (m s<sup>-1</sup>)</td><td>J[1, 100]</td><td>13.40<sup>+3.10</sup><sub>-2.39</sub></td><td></td></tr><tr><td>Jitter TESS (ppm)</td><td>J[1, 100]</td><td>204<sup>+14</sup><sub>-15</sub></td><td></td></tr><tr><td colspan="4"><strong>Derived Parameters TOI 4633 b</strong></td></tr><tr><td>Planet minimum mass M<sub>p</sub> sin i (M<sub>⊕</sub>)<sup>c</sup></td><td>...</td><td>106.8<sup>+13.0</sup><sub>-12.8</sub></td><td></td></tr><tr><td>Eccentricity e</td><td>...</td><td>0.096<sup>+0.102</sup><sub>-0.065</sub></td><td></td></tr><tr><td>Argument of periastron w (deg)</td><td>...</td><td>-43.9<sup>+104.8</sup><sub>-72.8</sub></td><td></td></tr><tr><td colspan="4"><strong>Derived Parameters TOI 4633 c</strong></td></tr><tr><td>Planet mass (M<sub>⊕</sub>)<sup>c</sup></td><td>...</td><td>47.8<sup>+27.6</sup><sub>-23.8</sub></td><td>123 M<sub>⊕</sub>, 99% percent upper limit</td></tr><tr><td>Observed planet radius (R<sub>⊕</sub>)</td><td>...</td><td>2.42<sup>+0.15</sup><sub>-0.14</sub></td><td></td></tr><tr><td>Corrected planet radius (R<sub>⊕</sub>)</td><td>...</td><td>3.2<sup>+0.20</sup><sub>-0.19</sub></td><td>Light contribution corrected</td></tr><tr><td>Semimajor axis a (au)</td><td>...</td><td>0.847 ± 0.061</td><td></td></tr><tr><td>Eccentricity e</td><td>...</td><td>0.117<sup>+0.186</sup><sub>-0.085</sub></td><td></td></tr><tr><td>Argument of periastron w (deg)</td><td>...</td><td>-21<sup>+131</sup><sub>-108</sub></td><td></td></tr><tr><td>Transit duration τ (hours)</td><td>...</td><td>11.45<sup>+0.46</sup><sub>-0.28</sub></td><td></td></tr><tr><td>Orbit inclination i (deg)</td><td>...</td><td>89.888<sup>+0.069</sup><sub>-0.064</sub></td><td></td></tr><tr><td>Insolation F<sub>p</sub> (F<sub>⊕</sub>)</td><td>...</td><td>1.56<sup>+0.20</sup><sub>-0.16</sub></td><td></td></tr></tbody></table>

**Notes.** All parameters are calculated based on the assumption that both planets orbit star A.

Denham et al. 2019). Therefore, destabilizing planet–planet interactions, such as scattering, are considered.

Figure 13 shows the results of 600 N-body simulations of the TOI 4633 system (with TOI 4633 A hosting the planets) under three scenarios. The scenarios are (from left-to-right) (1) where TOI 4633 b is not included, (2) where TOI 4633 b has an initial inclination that narrowly avoids transit (i ~ 88°) but initially has a mutual inclination of 2° or less with TOI 4633 c, and (3) where TOI 4633 b is initially coplanar with the binary orbit. For each scenario, the average eccentricity of TOI 4633 c during the last 10% of each simulation is plotted against TOI 4633 c's initial mutual inclination with the binary orbit. In each run, the longitude of ascending node and argument of periapsis of TOI 4633 c are randomized, its eccentricity is

initialized to 0.118, and TOI 4633 b (if present) has the same longitude of ascending node and argument of periapsis as TOI 4633 c in the second scenario or the binary orbit in the third scenario. The runs were performed using the N-body code MERCURY (Chambers 1999). This version of MERCURY includes the first post-Newtonian term accounting for general-relativistic precession, (M. Payne, private communication).

The possible orbits of the planets in the system can be constrained by eliminating orbits that rapidly become unstable. The runs presented in Figure 13 only ran for 10⁶ yr, a short time compared to the likely age of the system (estimated to be around 1.3 Gyr), meaning that unstable initial conditions can be excluded from consideration, and the orbits of the planets can be constrained to only orbits that remain stable. We find that

<footnote><sup>a</sup> U[a, b] refers to an uniform prior between a and b, N[a, b] to a Gaussian prior with mean a and standard deviation b, and J[a, b] to the modified Jeffrey's prior as defined by (Gregory 2005 Equation (16)).</footnote>

<footnote><sup>b</sup> Inferred parameters and errors are defined as the median and 68.3 % credible interval of the posterior distribution.</footnote>

<footnote><sup>c</sup> Due to the unknown effect of the second star on the observed semi-amplitude, these are lower limits.</footnote>
<table><caption>**Table 4**
System Parameters</caption><thead><tr><th>Parameter</th><th>Prior<sup>a</sup></th><th>Value<sup>b</sup></th><th>Comments</th></tr></thead><tbody><tr><td colspan="4"><strong>Model Parameters for TOI 4633 b</strong></td></tr><tr><td>Orbital period P<sub>orb</sub> (days)</td><td>U[33, 36]</td><td>34.15 ± 0.15</td><td></td></tr><tr><td>Time of min. conjunction T<sub>0</sub> (BJD—2457000)</td><td>U[2787.32, 2807.32]</td><td>2796.64<sup>+1.10</sup><sub>-1.19</sub></td><td></td></tr><tr><td>Parameterization √e sin ω</td><td>U[-1, 1]</td><td>-0.12<sup>+0.26</sup><sub>-0.22</sub></td><td>The code ensures e &lt; 1</td></tr><tr><td>Parameterization √e cos ω</td><td>U[-1, 1]</td><td>0.07<sup>+0.22</sup><sub>-0.24</sub></td><td>The code ensures e &lt; 1</td></tr><tr><td>Doppler semi-amplitude K (m s<sup>-1</sup>)<sup>c</sup></td><td>U[0, 30]</td><td>19.97<sup>+2.29</sup><sub>-2.30</sub></td><td></td></tr><tr><td colspan="4"><strong>Model Parameters for TOI 4633 c</strong></td></tr><tr><td>Orbital period P<sub>orb</sub> (days)</td><td>U[271.9, 272.1]</td><td>271.9445<sup>+0.0039</sup><sub>-0.0040</sub></td><td></td></tr><tr><td>Transit epoch T<sub>0</sub> (BJD—2457000)</td><td>U[1863.81, 1865.81]</td><td>1864.8265<sup>+0.0088</sup><sub>-0.0083</sub></td><td></td></tr><tr><td>Parameterization √e sin ω</td><td>U[-1, 1]</td><td>-0.06<sup>+0.23</sup><sub>-0.27</sub></td><td>The code ensures e &lt; 1</td></tr><tr><td>Parameterization √e cos ω</td><td>U[-1, 1]</td><td>0.04<sup>+0.31</sup><sub>-0.35</sub></td><td>The code ensures e &lt; 1</td></tr><tr><td>Observed scaled planet radius R<sub>p</sub>/R<sub>⋆</sub></td><td>U[0, 0.1]</td><td>0.0211<sup>+0.00076</sup><sub>-0.00069</sub></td><td></td></tr><tr><td>Impact parameter b</td><td>U[0, 1.1]</td><td>0.33<sup>+0.24</sup><sub>-0.21</sub></td><td></td></tr><tr><td>Doppler semi-amplitude K (m s<sup>-1</sup>)<sup>c</sup></td><td>U[0, 30]</td><td>4.58<sup>+2.56</sup><sub>-2.29</sub></td><td>11.5 m s<sup>-1</sup>, 99% percent upper limit</td></tr><tr><td colspan="4"><strong>Other Parameters</strong></td></tr><tr><td>Stellar density ρ<sub>⋆</sub> (g cm<sup>-3</sup>)</td><td>N[1.34, 0.23]</td><td>1.35 ± 0.21</td><td></td></tr><tr><td>Parameterized limb-darkening coefficient q<sub>1</sub></td><td>U[0, 1]</td><td>0.39<sup>+0.37</sup><sub>-0.23</sub></td><td>q<sub>1</sub> parameter as in Kipping (2013)</td></tr><tr><td>Parameterized limb-darkening coefficient q<sub>2</sub></td><td>U[0, 1]</td><td>0.34<sup>+0.31</sup><sub>-0.22</sub></td><td>q<sub>2</sub> parameter as in Kipping (2013)</td></tr><tr><td>Offset velocity HIRES (km s<sup>-1</sup>)</td><td>U[-0.50, 0.50]</td><td>0.0023<sup>+0.0019</sup><sub>-0.0020</sub></td><td></td></tr><tr><td>Offset velocity SOPHIE (km s<sup>-1</sup>)</td><td>U[-0.50, 0.50]</td><td>-0.0007 ± 0.0034</td><td></td></tr><tr><td>Jitter HIRES (m s<sup>-1</sup>)</td><td>J[1, 100]</td><td>12.36<sup>+1.36</sup><sub>-1.20</sub></td><td></td></tr><tr><td>Jitter SOPHIE (m s<sup>-1</sup>)</td><td>J[1, 100]</td><td>13.40<sup>+3.10</sup><sub>-2.39</sub></td><td></td></tr><tr><td>Jitter TESS (ppm)</td><td>J[1, 100]</td><td>204<sup>+14</sup><sub>-15</sub></td><td></td></tr><tr><td colspan="4"><strong>Derived Parameters TOI 4633 b</strong></td></tr><tr><td>Planet minimum mass M<sub>p</sub> sin i (M<sub>⊕</sub>)<sup>c</sup></td><td>...</td><td>106.8<sup>+13.0</sup><sub>-12.8</sub></td><td></td></tr><tr><td>Eccentricity e</td><td>...</td><td>0.096<sup>+0.102</sup><sub>-0.065</sub></td><td></td></tr><tr><td>Argument of periastron w (deg)</td><td>...</td><td>-43.9<sup>+104.8</sup><sub>-72.8</sub></td><td></td></tr><tr><td colspan="4"><strong>Derived Parameters TOI 4633 c</strong></td></tr><tr><td>Planet mass (M<sub>⊕</sub>)<sup>c</sup></td><td>...</td><td>47.8<sup>+27.6</sup><sub>-23.8</sub></td><td>123 M<sub>⊕</sub>, 99% percent upper limit</td></tr><tr><td>Observed planet radius (R<sub>⊕</sub>)</td><td>...</td><td>2.42<sup>+0.15</sup><sub>-0.14</sub></td><td></td></tr><tr><td>Corrected planet radius (R<sub>⊕</sub>)</td><td>...</td><td>3.2<sup>+0.20</sup><sub>-0.19</sub></td><td>Light contribution corrected</td></tr><tr><td>Semimajor axis a (au)</td><td>...</td><td>0.847 ± 0.061</td><td></td></tr><tr><td>Eccentricity e</td><td>...</td><td>0.117<sup>+0.186</sup><sub>-0.085</sub></td><td></td></tr><tr><td>Argument of periastron w (deg)</td><td>...</td><td>-21<sup>+131</sup><sub>-108</sub></td><td></td></tr><tr><td>Transit duration τ (hours)</td><td>...</td><td>11.45<sup>+0.46</sup><sub>-0.28</sub></td><td></td></tr><tr><td>Orbit inclination i (deg)</td><td>...</td><td>89.888<sup>+0.069</sup><sub>-0.064</sub></td><td></td></tr><tr><td>Insolation F<sub>p</sub> (F<sub>⊕</sub>)</td><td>...</td><td>1.56<sup>+0.20</sup><sub>-0.16</sub></td><td></td></tr></tbody></table>

**Notes.** All parameters are calculated based on the assumption that both planets orbit star A.

Denham et al. 2019). Therefore, destabilizing planet–planet interactions, such as scattering, are considered.

Figure 13 shows the results of 600 N-body simulations of the TOI 4633 system (with TOI 4633 A hosting the planets) under three scenarios. The scenarios are (from left-to-right) (1) where TOI 4633 b is not included, (2) where TOI 4633 b has an initial inclination that narrowly avoids transit (i ~ 88°) but initially has a mutual inclination of 2° or less with TOI 4633 c, and (3) where TOI 4633 b is initially coplanar with the binary orbit. For each scenario, the average eccentricity of TOI 4633 c during the last 10% of each simulation is plotted against TOI 4633 c's initial mutual inclination with the binary orbit. In each run, the longitude of ascending node and argument of periapsis of TOI 4633 c are randomized, its eccentricity is

initialized to 0.118, and TOI 4633 b (if present) has the same longitude of ascending node and argument of periapsis as TOI 4633 c in the second scenario or the binary orbit in the third scenario. The runs were performed using the N-body code MERCURY (Chambers 1999). This version of MERCURY includes the first post-Newtonian term accounting for general-relativistic precession, (M. Payne, private communication).

The possible orbits of the planets in the system can be constrained by eliminating orbits that rapidly become unstable. The runs presented in Figure 13 only ran for 10⁶ yr, a short time compared to the likely age of the system (estimated to be around 1.3 Gyr), meaning that unstable initial conditions can be excluded from consideration, and the orbits of the planets can be constrained to only orbits that remain stable. We find that

<footnote><sup>a</sup> U[a, b] refers to an uniform prior between a and b, N[a, b] to a Gaussian prior with mean a and standard deviation b, and J[a, b] to the modified Jeffrey's prior as defined by (Gregory 2005 Equation (16)).</footnote>

<footnote><sup>b</sup> Inferred parameters and errors are defined as the median and 68.3 % credible interval of the posterior distribution.</footnote>

<footnote><sup>c</sup> Due to the unknown effect of the second star on the observed semi-amplitude, these are lower limits.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** RV time-series (top) and phase-folded plots for TOI 4633 b (lower left) and TOI 4633 c (lower right). SOPHIE (blue circles) and HIRES (orange diamonds) RV measurements are shown following the subtraction of the systemic velocities for each instrument. The light colored error bars show the uncertainties accounting for the jitter. Solid black lines show the inferred median model. Gray lines show models made with 100 random samples from the posterior distributions.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Phase-folded light curve for TOI 4633 c. 2 minutes and 20 s cadence TESS observations are shown in light gray. Solid color circles represent 60 minutes binned data. The inferred transit model is shown with a solid black line.</figcaption></figure>

without TOI 4633 b, TOI 4633 c is stable for a wide range of mutual inclinations with the binary—including many retrograde orbits, where mutual inclinations with the binary below 45° or above 125° (retrograde) remain stable. If TOI 4633 b is included coplanar with TOI 4633 c, TOI 4633 c again remains stable if its mutual inclination with the binary is below 45°, but its range of stable retrograde orbits is lessened to above 145°. When TOI 4633 b is initially coplanar with the binary, the stable regime of TOI 4633 c expands slightly for both prograde orbits—stable up to ~55°, but not for retrograde orbits where it is still only stable if its mutual inclination with the binary is above 145°.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Planet radius vs. planet mass for all confirmed planets listed in the NASA exoplanet archive with a mass measurement with 30% uncertainty or better. Bright systems ($m_V < 10$) are highlighted with a black outline. TOI 4633 c is shown by the blue triangle, where uncertainties correspond to our radial velocity fit. Although the existence of a second star in this system means that the radial velocity amplitude is diluted and thus the true planet mass may be larger, the lack of planets with the radius of TOI 4633 c at masses larger than 47 $M_⊕$ suggests that the planet is unlikely to be significantly more massive than what is reported here.</figcaption></figure>

Between the three scenarios, we find that the presence of TOI 4633 b is a lightly destabilizing presence in the system, reducing the number of stable retrograde orbital configurations. Planet–planet interactions are not sufficient to suppress eccentricity excitation of TOI 4633 c against perturbations from the binary if there is significant mutual inclination between the planets and the binary. Moreover, because planet–planet interactions are not stabilizing TOI 4633 c against
<figure><img src="image_13.png" /><figcaption>**Figure 9.** RV time-series (top) and phase-folded plots for TOI 4633 b (lower left) and TOI 4633 c (lower right). SOPHIE (blue circles) and HIRES (orange diamonds) RV measurements are shown following the subtraction of the systemic velocities for each instrument. The light colored error bars show the uncertainties accounting for the jitter. Solid black lines show the inferred median model. Gray lines show models made with 100 random samples from the posterior distributions.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Phase-folded light curve for TOI 4633 c. 2 minutes and 20 s cadence TESS observations are shown in light gray. Solid color circles represent 60 minutes binned data. The inferred transit model is shown with a solid black line.</figcaption></figure>

without TOI 4633 b, TOI 4633 c is stable for a wide range of mutual inclinations with the binary—including many retrograde orbits, where mutual inclinations with the binary below 45° or above 125° (retrograde) remain stable. If TOI 4633 b is included coplanar with TOI 4633 c, TOI 4633 c again remains stable if its mutual inclination with the binary is below 45°, but its range of stable retrograde orbits is lessened to above 145°. When TOI 4633 b is initially coplanar with the binary, the stable regime of TOI 4633 c expands slightly for both prograde orbits—stable up to ~55°, but not for retrograde orbits where it is still only stable if its mutual inclination with the binary is above 145°.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Planet radius vs. planet mass for all confirmed planets listed in the NASA exoplanet archive with a mass measurement with 30% uncertainty or better. Bright systems ($m_V < 10$) are highlighted with a black outline. TOI 4633 c is shown by the blue triangle, where uncertainties correspond to our radial velocity fit. Although the existence of a second star in this system means that the radial velocity amplitude is diluted and thus the true planet mass may be larger, the lack of planets with the radius of TOI 4633 c at masses larger than 47 $M_⊕$ suggests that the planet is unlikely to be significantly more massive than what is reported here.</figcaption></figure>

Between the three scenarios, we find that the presence of TOI 4633 b is a lightly destabilizing presence in the system, reducing the number of stable retrograde orbital configurations. Planet–planet interactions are not sufficient to suppress eccentricity excitation of TOI 4633 c against perturbations from the binary if there is significant mutual inclination between the planets and the binary. Moreover, because planet–planet interactions are not stabilizing TOI 4633 c against
<figure><img src="image_14.png" /><figcaption>**Figure 12.** The recovery completeness of injected transit signals into the light curve of TOI 4633 as a function of the radius and orbital period. The radius of planet candidate b was estimated using the planet mass–radius scaling relations for volatile-rich planets by Otegi et al. (2020). The signals were recovered using a BLS search. The figure clearly highlights that if planet candidate b were aligned such that it would transit star A, the transit events would have been detected by the BLS search algorithm.</figcaption></figure>

perturbation, the influence of the masses of TOI 4633 b and TOI 4633 c within their allowed ranges are minor. This is because scattering is caused by perturbations by the binary rather than arising from the planet–planet interactions, so whether scattering occurs is not affected by varying the planet masses within their allowed uncertainties.

We note that the choice of host star (whether the planets orbit TOI 4633 A or TOI 4633 B) is not expected to significantly affect these results. Using the grids by Quarles et al. (2020), we show that for mutual inclination ($i_{m,bin-c}$) between 0° and 45°, planet c would be stable orbiting around either of the two stars in the binary. As such, with the current available data, we cannot use stability arguments to determine around which of the two stars the planets are orbiting. Furthermore, as previously discussed, due to the similar masses and radii of the two stars in the binary, the derived planet properties are the same to within their uncertainties with either of the two stars as the planet host.

# 7. Discussion and Conclusion

The TOI 4633 system consists of two G-type stars, a 272 days period transiting planet and a 34 days period non-transiting planet candidate. The system stands out from other transiting planet systems due to (i) the long orbital period of planet c which places it in the habitable zone; (ii) the brightness of the system; and (iii) the similar-mass, bound stellar companion.

## 7.1. Long-period Planet around a Bright Star

To date, there is a distinct lack of confirmed transiting planets with long orbital periods, with only 175 planets with periods longer than 100 days, and 40 with periods longer than 250 days (only 5 of which are brighter than Vmag = 12).<sup>43</sup> As such, with a period of ~272 days and a host star brightness of Vmag = 9.0, planet c lies in an under-explored region of parameter space, as shown in Figure 14. The figure shows the planet mass versus orbital period for all confirmed planets listed in the NASA exoplanet archive with a mass measurement with 30% uncertainty or better as gray points, and highlights bright systems ($m_V < 10$) in color. Bright systems with only

one planet detected are shown as orange circles, while bright systems with more than one planet are shown as purple squares. Planet candidate b and c are shown by the black cross and triangle, respectively. The figure highlights a lack of confirmed planets with mass measurements to better that 30% residing in bright ($V<10$) multi-planet systems. Furthermore, it shows that planet c is the second longest-period planet known around a bright star; and that if both planets are orbiting around the same star, TOI 4633 would be one of only a few bright, multi-planet systems. The difficulty of detecting long-period planets is illustrated by the fact that the transit probability of a planet at the semimajor axis of TOI 4633 c is 0.6%. Overall, this highlights the importance of further characterization of TOI 4633 in order to help further our understanding of planetary system demographics. Particularly, the brightness of the host star makes TOI 4633 c a prime candidate for further ground- and space-based characterization.

## 7.2. Two Planets, Two Stars

The combination of new and archival high-contrast imaging data dating back to 1905 showed that the system is comprised of two stars ($a_{bin} = 48.6_{-4.4}^{+3.5}$ au, $e_{bin} = 0.91_{-0.03}^{+0.03}$, $P_{bin} = 231_{-32}^{+24}$ yr). Due to the proximity of the two stars, we are unable to determine around which star the planets orbit, or whether the two planets orbit the same star.

Despite the large fraction of stars that reside in binaries, the known sample of confirmed planets in binaries remains limited. The catalog of exoplanets in binary star systems (Schwarz et al. 2016) lists 154 systems containing a total of 217 planets. Out of these, 27 are P-type (“circumbinary planets”) and 190 are S-type (“circumstellar planets”). The properties of the S-type planets are shown in Figure 15 on the binary semimajor axis versus planet orbital period plane. Planets listed as having been detected using the transit method (74 planets) and the radial velocity method (107 planets) are depicted by black triangles and blue crosses, respectively. The properties of TOI 4633 b and c are shown in red. As highlighted by this figure, the stellar semimajor axis of TOI 4633 AB, of ~48.6 au, places TOI 4633 in an under-sampled region of parameter space, with only 18 confirmed planets around a star with a binary semimajor axis less than that of TOI 4633 AB, only two of which are transiting. Furthermore, it highlights that planet c has the longest orbital period of any confirmed transiting planet in an S-type binary. In addition, there are currently only a handful of circumstellar systems where both stars are known to host a planet (e.g., Teske et al. 2016a, 2016b). Thus, distinguishing around which star these planets orbit will, in the future, allow for a study of how differences in the host star properties (e.g., chemical abundances) can result in different planet properties.

Stellar multiplicity is widely believed to affect planet formation and evolution. However, the extent and details of the effect of a companion star remain topics of debate. For example, studies have suggested that giant planets on short orbital periods are preferentially found in systems with wide stellar companions (e.g., Wang et al. 2015a; Ngo et al. 2016; Ziegler et al. 2018; Moe & Kratter 2021a; Ziegler et al. 2021) due to the effect of the companion star triggering planet inward migration. Similarly, Moe & Kratter (2021a) combined a variety of RV and high-resolution imaging surveys to show that planet occurrence rates are suppressed as a function of binary separation (i.e., smaller binary semimajor axes result in a lower occurrence of giant planets). Moe & Kratter (2021a) also show

<footnote n="43">43 NASA Exoplanet Archive.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** The recovery completeness of injected transit signals into the light curve of TOI 4633 as a function of the radius and orbital period. The radius of planet candidate b was estimated using the planet mass–radius scaling relations for volatile-rich planets by Otegi et al. (2020). The signals were recovered using a BLS search. The figure clearly highlights that if planet candidate b were aligned such that it would transit star A, the transit events would have been detected by the BLS search algorithm.</figcaption></figure>

perturbation, the influence of the masses of TOI 4633 b and TOI 4633 c within their allowed ranges are minor. This is because scattering is caused by perturbations by the binary rather than arising from the planet–planet interactions, so whether scattering occurs is not affected by varying the planet masses within their allowed uncertainties.

We note that the choice of host star (whether the planets orbit TOI 4633 A or TOI 4633 B) is not expected to significantly affect these results. Using the grids by Quarles et al. (2020), we show that for mutual inclination ($i_{m,bin-c}$) between 0° and 45°, planet c would be stable orbiting around either of the two stars in the binary. As such, with the current available data, we cannot use stability arguments to determine around which of the two stars the planets are orbiting. Furthermore, as previously discussed, due to the similar masses and radii of the two stars in the binary, the derived planet properties are the same to within their uncertainties with either of the two stars as the planet host.

# 7. Discussion and Conclusion

The TOI 4633 system consists of two G-type stars, a 272 days period transiting planet and a 34 days period non-transiting planet candidate. The system stands out from other transiting planet systems due to (i) the long orbital period of planet c which places it in the habitable zone; (ii) the brightness of the system; and (iii) the similar-mass, bound stellar companion.

## 7.1. Long-period Planet around a Bright Star

To date, there is a distinct lack of confirmed transiting planets with long orbital periods, with only 175 planets with periods longer than 100 days, and 40 with periods longer than 250 days (only 5 of which are brighter than Vmag = 12).<sup>43</sup> As such, with a period of ~272 days and a host star brightness of Vmag = 9.0, planet c lies in an under-explored region of parameter space, as shown in Figure 14. The figure shows the planet mass versus orbital period for all confirmed planets listed in the NASA exoplanet archive with a mass measurement with 30% uncertainty or better as gray points, and highlights bright systems ($m_V < 10$) in color. Bright systems with only

one planet detected are shown as orange circles, while bright systems with more than one planet are shown as purple squares. Planet candidate b and c are shown by the black cross and triangle, respectively. The figure highlights a lack of confirmed planets with mass measurements to better that 30% residing in bright ($V<10$) multi-planet systems. Furthermore, it shows that planet c is the second longest-period planet known around a bright star; and that if both planets are orbiting around the same star, TOI 4633 would be one of only a few bright, multi-planet systems. The difficulty of detecting long-period planets is illustrated by the fact that the transit probability of a planet at the semimajor axis of TOI 4633 c is 0.6%. Overall, this highlights the importance of further characterization of TOI 4633 in order to help further our understanding of planetary system demographics. Particularly, the brightness of the host star makes TOI 4633 c a prime candidate for further ground- and space-based characterization.

## 7.2. Two Planets, Two Stars

The combination of new and archival high-contrast imaging data dating back to 1905 showed that the system is comprised of two stars ($a_{bin} = 48.6_{-4.4}^{+3.5}$ au, $e_{bin} = 0.91_{-0.03}^{+0.03}$, $P_{bin} = 231_{-32}^{+24}$ yr). Due to the proximity of the two stars, we are unable to determine around which star the planets orbit, or whether the two planets orbit the same star.

Despite the large fraction of stars that reside in binaries, the known sample of confirmed planets in binaries remains limited. The catalog of exoplanets in binary star systems (Schwarz et al. 2016) lists 154 systems containing a total of 217 planets. Out of these, 27 are P-type (“circumbinary planets”) and 190 are S-type (“circumstellar planets”). The properties of the S-type planets are shown in Figure 15 on the binary semimajor axis versus planet orbital period plane. Planets listed as having been detected using the transit method (74 planets) and the radial velocity method (107 planets) are depicted by black triangles and blue crosses, respectively. The properties of TOI 4633 b and c are shown in red. As highlighted by this figure, the stellar semimajor axis of TOI 4633 AB, of ~48.6 au, places TOI 4633 in an under-sampled region of parameter space, with only 18 confirmed planets around a star with a binary semimajor axis less than that of TOI 4633 AB, only two of which are transiting. Furthermore, it highlights that planet c has the longest orbital period of any confirmed transiting planet in an S-type binary. In addition, there are currently only a handful of circumstellar systems where both stars are known to host a planet (e.g., Teske et al. 2016a, 2016b). Thus, distinguishing around which star these planets orbit will, in the future, allow for a study of how differences in the host star properties (e.g., chemical abundances) can result in different planet properties.

Stellar multiplicity is widely believed to affect planet formation and evolution. However, the extent and details of the effect of a companion star remain topics of debate. For example, studies have suggested that giant planets on short orbital periods are preferentially found in systems with wide stellar companions (e.g., Wang et al. 2015a; Ngo et al. 2016; Ziegler et al. 2018; Moe & Kratter 2021a; Ziegler et al. 2021) due to the effect of the companion star triggering planet inward migration. Similarly, Moe & Kratter (2021a) combined a variety of RV and high-resolution imaging surveys to show that planet occurrence rates are suppressed as a function of binary separation (i.e., smaller binary semimajor axes result in a lower occurrence of giant planets). Moe & Kratter (2021a) also show

<footnote n="43"><sup>43</sup> NASA Exoplanet Archive.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Outcome of simulations organized by mass and initial mutual inclination between TOI 4633c and the binary orbit. We show 600 MERCURY integrations run for $10^6$ yr each across three scenarios indicated by the labels on the figures: (1) without the 34 days planet, (2) where TOI 4633b has initial inclination of ~88° narrowly avoiding transit and is coplanar with TOI 4633c, and (3) where TOI 4633b is coplanar with the binary orbit. The x-axis is the initial mutual inclination set between TOI 4633c and the binary orbit, and the y-axis is the average eccentricity of TOI 4633c during the final 10% of the integration time. Each simulation is represented by a point-blue for runs that remain stable and red for unstable.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Planet mass vs. planet orbital period for all systems with masses measured to 30% or better uncertainty. We have highlighted planets around bright stars in color. Within the bright star sample, we distinguish between single-planet systems (orange circles) and systems with more than one planet (purple squares). TOI 4633 c is the second longest-period planet in a bright host star system (however, we do note that planet c does not have a mass measured to 30% or better uncertainty). Furthermore, both planets in the system are fairly high mass and likely not rocky, which is unusual compared to the confirmed population of multi-planet systems.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Binary semimajor axis vs. planet period for planetary systems with more than one star. Transiting planets are shown as triangles, while radial velocity-detected planets are shown as crosses. TOI 4633 planets are shown in red. TOI 4633 has a smaller binary separation than most planet-hosting binary systems, and TOI 4633 c is the longest-period transiting planet currently known in a binary system.</figcaption></figure>

that at a binary separation of TOI 4633 AB ~48.6 au, the observed occurrence rate of planets is around 50% less compared to field stars. Other studies suggest that there is a suppression of transiting planets in S-type systems, due to the stellar companion disrupting the orbital coplanarity (e.g., Wang et al. 2015a).

Recent studies have also been investigating whether binary orbital parameters (e.g., binary inclination and eccentricity) affect planet formation and observed planet properties. For example, Behmard et al. (2022), Christian et al. (2022), Lester et al. (2023), and Dupuy et al. (2022) provide observational evidence that suggests that the orbital planes of the binary stars are preferentially aligned with the orbital planes of the planet(s).

In addition to allowing for population studies that can inform theories of planet formation and migration, planets in binary

systems are interesting due to the fact that the two stars formed simultaneously and thus are expected to have the same chemical abundance at the time of formation. As such, any observed differences in the chemical composition of the stars could be related to differences in the outcomes of planet formation and therefore may help provide constraints on how planet formation affects stellar properties (e.g., Teske et al. 2016a). In turn, the identification of a correlation between certain chemical abundances and planet formation could help improve target selection of future space missions searching for planets, based on a single spectrum of the star. Alternatively, observed differences in chemical compositions between the two stars could indicate recent planet engulfment (Oh et al. 2018; Behmard et al. 2023). Even though the two stars in this binary are currently too close to one another to be able to separate their spectra and therefore to measure their individual chemical abundances, our model of the orbit indicates that in ~30 yr time the two stars will be separated by ~150 mas, which will allow us to observe the individual stars using
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Outcome of simulations organized by mass and initial mutual inclination between TOI 4633c and the binary orbit. We show 600 MERCURY integrations run for $10^6$ yr each across three scenarios indicated by the labels on the figures: (1) without the 34 days planet, (2) where TOI 4633b has initial inclination of ~88° narrowly avoiding transit and is coplanar with TOI 4633c, and (3) where TOI 4633b is coplanar with the binary orbit. The x-axis is the initial mutual inclination set between TOI 4633c and the binary orbit, and the y-axis is the average eccentricity of TOI 4633c during the final 10% of the integration time. Each simulation is represented by a point-blue for runs that remain stable and red for unstable.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Planet mass vs. planet orbital period for all systems with masses measured to 30% or better uncertainty. We have highlighted planets around bright stars in color. Within the bright star sample, we distinguish between single-planet systems (orange circles) and systems with more than one planet (purple squares). TOI 4633 c is the second longest-period planet in a bright host star system (however, we do note that planet c does not have a mass measured to 30% or better uncertainty). Furthermore, both planets in the system are fairly high mass and likely not rocky, which is unusual compared to the confirmed population of multi-planet systems.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Binary semimajor axis vs. planet period for planetary systems with more than one star. Transiting planets are shown as triangles, while radial velocity-detected planets are shown as crosses. TOI 4633 planets are shown in red. TOI 4633 has a smaller binary separation than most planet-hosting binary systems, and TOI 4633 c is the longest-period transiting planet currently known in a binary system.</figcaption></figure>

that at a binary separation of TOI 4633 AB ~48.6 au, the observed occurrence rate of planets is around 50% less compared to field stars. Other studies suggest that there is a suppression of transiting planets in S-type systems, due to the stellar companion disrupting the orbital coplanarity (e.g., Wang et al. 2015a).

Recent studies have also been investigating whether binary orbital parameters (e.g., binary inclination and eccentricity) affect planet formation and observed planet properties. For example, Behmard et al. (2022), Christian et al. (2022), Lester et al. (2023), and Dupuy et al. (2022) provide observational evidence that suggests that the orbital planes of the binary stars are preferentially aligned with the orbital planes of the planet(s).

In addition to allowing for population studies that can inform theories of planet formation and migration, planets in binary

systems are interesting due to the fact that the two stars formed simultaneously and thus are expected to have the same chemical abundance at the time of formation. As such, any observed differences in the chemical composition of the stars could be related to differences in the outcomes of planet formation and therefore may help provide constraints on how planet formation affects stellar properties (e.g., Teske et al. 2016a). In turn, the identification of a correlation between certain chemical abundances and planet formation could help improve target selection of future space missions searching for planets, based on a single spectrum of the star. Alternatively, observed differences in chemical compositions between the two stars could indicate recent planet engulfment (Oh et al. 2018; Behmard et al. 2023). Even though the two stars in this binary are currently too close to one another to be able to separate their spectra and therefore to measure their individual chemical abundances, our model of the orbit indicates that in ~30 yr time the two stars will be separated by ~150 mas, which will allow us to observe the individual stars using
<figure><img src="image_18.png" /></figure>

**Figure 16.** Stellar effective temperature vs. effective incident flux on a planet. The boundaries of the habitable zone defined by Kopparapu et al. (2014) are illustrated as shaded regions. The black circles show all transiting exoplanets listed in the NASA exoplanet archive, where the filled circles show only planets that have a mass measurement to better than 30% accuracy. The orange triangle shows the properties of TOI 4633 c. There are no other transiting planets within the habitable zone that are brighter than 11th magnitude, and only four that are brighter than 13th magnitude.

instruments such as the Keck Planet Imager and Characterizer
(Delorme et al. 2020).

7.3. TOI4633 c: A Mini-Neptune in the Habitable Zone

The search for exoplanets in the habitable zones of their host
stars has been a focal point in the field of exoplanet research,
where the habitable zone is defined as the region around a star
where the surface temperature is favorable for liquid water to
exist. We estimate the location of the habitable zone around
TOI 4633 following the prescriptions defined by Kopparapu
et al. (2014). The habitable zone boundaries of “Recent
Venus,” “Maximum Greenhouse,” and “Runaway Green-
house” conditions are shown in Figure 16 (for more details
on how these boundaries are defined see Kopparapu et al.
2014). Using this definition of the habitable zone we find that
TOI 4633 c, with an insolation of 1.56<sup>+0.20</sup><sub>-0.16</sub> F<sub>⊕</sub>, lies inside the
inner edge of the potentially habitable region of the star
between the Recent Venus and Maximum Greenhouse
boundaries. Using the Stefan–Boltzman law and assuming an
albedo equivalent to that of Neptune (0.29), we show that the
effective surface temperature of planet c is ~290 K when only
the energy contribution of star A is considered. When
accounting for the energy contribution of stars A and B at
the time of periastron (stellar separation ~4.5 au), the effective
surface temperature increases to ~420 K.

According to the Catalog of Habitable Zone Exoplanets by Hill et al. (2023), TOI 4633 is the brightest star known to host a transiting planet in the habitable zone, with no other transiting planets within the habitable zone that are brighter than 11th magnitude, five that are brighter than 13th magnitude, and 29 that are brighter than 15th magnitude. Even though planet c is at a distance from its host star where liquid water could potentially reside, the density of the planet (>0.11 g cm<sup>-3</sup>) suggests that it has a large and dense atmosphere that likely makes surface level liquid water impossible. However, as shown independently by Sucerquia et al. (2019) and Dobos et al. (2021), the probability of a planet hosting a satellite increases with increased orbital period. As such, the brightness of the system and the long orbital period of planet c makes TOI 4633 c a valuable target for further characterization in the future and, in particular, for the search of planet satellites.

7.4. Conclusion

We present the discovery and validation of a transiting mini-
Neptune ($P_c = 271.9445_{-0.0040}^{+0.0039}$ days, $R_c = 3.2_{-0.19}^{+0.20} R_⊕$, $M_c =
47.8_{-23.8}^{+27.6} M_⊕$) that was discovered by citizen scientists taking
part in the Planet Hunters TESS citizen science project. The
planet's long orbital period places it in the optimistic habitable
zone of its host star, with an insolation of $1.56_{-0.16}^{+0.20} F_⊕$. RV
monitoring of the system with Keck/HIRES and OHP/
SOPHIE revealed an additional 34 days periodic signal that is
not seen in the photometric data. Due to a lack of strong
evidence that this periodic signal is caused by stellar activity,
we tentatively consider this signal to be a planet candidate
($M_b \text{sini} < 109 M_⊕$). Furthermore, we used high-contrast
imaging observations spanning over 117 yr (1905–2023) to
constrain the orbit of the bound stellar companion star ($a_{\text{bin}} =
48.6_{-3.5}^{+4.4}$ au, $e_{\text{bin}} = 0.91_{-0.03}^{+0.03}$).

Dynamical simulations were used to constrain the mutual
inclination between the binary orbit and that of TOI 4633 c.
The *N*-body simulations showed that if the 34 days planet
candidate is inclined such that it only narrowly avoids
transiting, then the mutual inclination between planet *c* and
the binary orbit has to be less than ~45° to ensure that the
system is stable for longer than $10^6$ yr.

Overall, TOI 4633 stands out due to its brightness, the long
orbital period of the transiting planet TOI 4633 c that places it
in the optimistic habitable zone, and the presence of the bound,
near equal mass, stellar companion. TOI 4633 c is currently
only the fourth habitable zone planet identified in the TESS
data, and the only one transiting a G dwarf star. The system is
valuable for studying planet formation and evolution in
binaries, as well as for being an excellent target for detailed
investigation of planets in habitable zones, particularly in
multiple star systems. While the two stars in the system cannot
currently be resolved from the ground (due to their proximity),
the binary orbit modeling showed that in around 30 yr time the
stars will be separated by more than 150 mas on the sky, a
distance large enough to be able to resolve with modern
spectrographs such as the Keck Planet Imager and Character-
izer. This will enable further characterization of the two stars
and the planets in the system.

8. Software and Data Availability

The TESS data used within this article are hosted and made
publicly available by the Mikulski Archive for Space
Telescopes (MAST, http://archive.stsci.edu/tess/). All the
TESS data used in this paper can be found in
MAST:10.17909/pepe-f853. This work also used data from
the NASA Exoplanet Archive (NASA Exoplanet
Archive 2019), and ExoFOP (ExoFOP 2019).

This work made use of Astropy, a community-developed
core Python package for Astronomy (Astropy Collaboration
et al. 2013), matplotlib (Hunter 2007), pandas (McKinney et al.
2010), NumPy (Walt 2011), astroquery (Ginsburg et al. 2019),
sklearn (Pedregosa et al. 2011), emcee (Foreman-Mackey et al.
2013), PTEMCEE (Vousden et al. 2016), TRICERATOPS
(Giacalone et al. 2021), LATTE (Eisner et al. 2020a), lightkurve
(Lightkurve Collaboration et al. 2018), pyaneti (Barragán et al.
2022a), ORBITIZE! (Blunt et al. 2020), Mercury (Chambers
1999), corner (Foreman-Mackey 2016).