# The GAPS programme at TNG. XLVIII. The unusual formation history of V1298 Tau

D. Turrini, F. Marzari, D. Polychroni, R. Claudi, S. Desidera, D. Mesa, M.

Pinamonti, A. Sozzetti, A. Suárez Mascareño, M. Damasso, et al.

► To cite this version:

D. Turrini, F. Marzari, D. Polychroni, R. Claudi, S. Desidera, et al. The GAPS programme at TNG. XLVIII. The unusual formation history of V1298 Tau. Astronomy & Astrophysics, 2023, 679, 10.1051/0004-6361/202244752. insu-04479041

HAL Id: insu-04479041

https://insu.hal.science/insu-04479041

Submitted on 10 Mar 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The GAPS programme at TNG. XLVIII. The unusual formation history of V1298 Tau

D. Turrini, F. Marzari, D. Polychroni, R. Claudi, S. Desidera, D. Mesa, M.

Pinamonti, A. Sozzetti, A. Suárez Mascareño, M. Damasso, et al.

► To cite this version:

D. Turrini, F. Marzari, D. Polychroni, R. Claudi, S. Desidera, et al. The GAPS programme at TNG. XLVIII. The unusual formation history of V1298 Tau. Astronomy & Astrophysics - A&A, 2023, 679, 10.1051/0004-6361/202244752. insu-04479041

HAL Id: insu-04479041

https://insu.hal.science/insu-04479041v1

Submitted on 10 Mar 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The GAPS programme at TNG

## XLVIII. The unusual formation history of V1298 Tau*

D. Turrini¹,²,³, F. Marzari⁴, D. Polychroni⁶, R. Claudi⁵, S. Desidera⁵, D. Mesa⁵, M. Pinamonti¹,
A. Sozzetti¹, A. Suárez Mascareño⁹,¹⁰, M. Damasso¹, S. Benatti⁷, L. Malavolta⁴, G. Micela⁷, A. Zinzi⁸,¹⁰,
V. J. S. Béjar⁹,¹⁰, K. Biazzo¹¹, A. Bignamini⁶, M. Bonavita¹²,¹³,⁵, F. Borsa¹⁴, C. del Burgo¹⁵, G. Chauvin¹⁶,
P. Delorme¹⁷, J. I. González Hernández⁹,¹⁰, R. Gratton⁵, J. Hagelberg¹⁸, M. Janson¹⁹, M. Langlois²⁰,
A. F. Lanza²¹, C. Lazzoni²², N. Lodieu⁹, A. Maggio⁷, L. Mancini²⁴,¹,²⁵, E. Molinari²³, M. Molinaro⁶,
F. Murgas⁹,¹⁰, and D. Nardiello⁵,²⁶

*(Affiliations can be found after the references)*

Received 15 August 2022 / Accepted 4 July 2023

## ABSTRACT

**Context.** Observational data from space- and ground-based campaigns have revealed that the 10-30 Ma old V1298 Tau star hosts a compact and massive system of four planets. Mass estimates are available for the two outer giant planets and point to unexpectedly high densities for their young ages.

**Aims.** We investigate the formation of these two outermost giant planets, V1298 Tau b and e, and the present dynamical state of V1298 Tau's global architecture in order to shed light on the history of this young and peculiar extrasolar system.

**Methods.** We performed detailed *N*-body simulations to explore the link between the densities of V1298 Tau b and e and their migration and accretion of planetesimals within the native circumstellar disk. We combined *N*-body simulations and the normalized angular momentum deficit (NAMD) analysis of the architecture to characterize V1298 Tau's dynamical state and connect it to the formation history of the system. We searched for outer planetary companions to constrain V1298 Tau's planetary architecture and the extension of its primordial circumstellar disk.

**Results.** The high densities of V1298 Tau b and e suggest they formed at quite a distance from their host star, likely beyond the CO₂ snowline. The higher nominal density of V1298 Tau e suggests it formed farther out than V1298 Tau b. The current architecture of V1298 Tau is not characterized by resonant chains. Planet-planet scattering with an outer giant planet is the most likely cause for the lack of a resonant chain between V1298 Tau's planets, but currently our search for outer companions using SPHERE and Gaia observations can exclude only the presence of planets more massive than $2 M_J$.

**Conclusions.** The most plausible scenario for V1298 Tau's formation is that the system formed by convergent migration and resonant trapping of planets born in a compact and plausibly massive disk. In the wake of their migration, V1298 Tau b and e would have left a dynamically excited protoplanetary disk, naturally creating the conditions for the later breaking of the resonant chain by planet-planet scattering.

**Key words.** planetary systems – stars: individual: V1298 Tau – planets and satellites: formation – planets and satellites: dynamical evolution and stability – planets and satellites: detection – chaos

## 1. Introduction

Young planets offer us the unique opportunity to study unaltered products of planet formation, before secular evolution and interactions with host stars modify or cancel the original characteristics of the planets. When they belong to multiplanet systems, young planets also represent invaluable case studies to comparatively investigate how different formation histories can shape planets born from the same stellar and disk environments.

Multiple authors (e.g., Donati et al. 2016; Yu et al. 2017; Benatti et al. 2021) have argued that the population of massive planets in close orbits around young stars can be significantly larger than that of their counterparts around older stars. The proposed roots of this difference are the mass loss in the early

evolutionary stages driven by the strong X and UV radiation from the host star, as well as the processes of orbital migration and chaotic evolution at play during the early phases of planetary systems that favor the engulfment (e.g., Donati et al. 2017; Spina et al. 2021) or the removal (e.g., Zinzi & Turrini 2017; Turrini et al. 2018, 2020, 2022) of planets.

Because young stars are usually very active, another reason for this difference could reside in the false positives induced by stellar activity. The number of known young planets has been questioned by recent works that have unveiled fake detections (e.g., Carleo et al. 2018; Donati et al. 2020; Damasso et al. 2020). In the past years, however, NASA *Kepler* and its K2 campaigns (Borucki et al. 2003; Howell et al. 2014) and recently the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) detected transits of several young and multiple planets (e.g., David et al. 2016; Mann et al. 2016; Rizzuto et al. 2020; Plavchan et al. 2020), allowing the study of the early evolution and architecture of these systems (e.g., Benatti et al. 2019, 2021).

<footnote>* Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundaciòn Galileo Galilei (FGG) of the Istituto Nazionale di Astrofisica (INAF) at the Observatorio del Roque de los Muchachos (La Palma, Canary Islands, Spain).</footnote>
# The GAPS programme at TNG

## XLVIII. The unusual formation history of V1298 Tau*

D. Turrini¹,²,³, F. Marzari⁴, D. Polychroni⁶, R. Claudi⁵, S. Desidera⁵, D. Mesa⁵, M. Pinamonti¹,
A. Sozzetti¹, A. Suárez Mascareño⁹,¹⁰, M. Damasso¹, S. Benatti⁷, L. Malavolta⁴, G. Micela⁷, A. Zinzi⁸,¹⁰,
V. J. S. Béjar⁹,¹⁰, K. Biazzo¹¹, A. Bignamini⁶, M. Bonavita¹²,¹³,⁵, F. Borsa¹⁴, C. del Burgo¹⁵, G. Chauvin¹⁶,
P. Delorme¹⁷, J. I. González Hernández⁹,¹⁰, R. Gratton⁵, J. Hagelberg¹⁸, M. Janson¹⁹, M. Langlois²⁰,
A. F. Lanza²¹, C. Lazzoni²², N. Lodieu⁹, A. Maggio⁷, L. Mancini²⁴,¹,²⁵, E. Molinari²³, M. Molinaro⁶,
F. Murgas⁹,¹⁰, and D. Nardiello⁵,²⁶

*(Affiliations can be found after the references)*

Received 15 August 2022 / Accepted 4 July 2023

## ABSTRACT

**Context.** Observational data from space- and ground-based campaigns have revealed that the 10-30 Ma old V1298 Tau star hosts a compact and massive system of four planets. Mass estimates are available for the two outer giant planets and point to unexpectedly high densities for their young ages.

**Aims.** We investigate the formation of these two outermost giant planets, V1298 Tau b and e, and the present dynamical state of V1298 Tau's global architecture in order to shed light on the history of this young and peculiar extrasolar system.

**Methods.** We performed detailed *N*-body simulations to explore the link between the densities of V1298 Tau b and e and their migration and accretion of planetesimals within the native circumstellar disk. We combined *N*-body simulations and the normalized angular momentum deficit (NAMD) analysis of the architecture to characterize V1298 Tau's dynamical state and connect it to the formation history of the system. We searched for outer planetary companions to constrain V1298 Tau's planetary architecture and the extension of its primordial circumstellar disk.

**Results.** The high densities of V1298 Tau b and e suggest they formed at quite a distance from their host star, likely beyond the CO₂ snowline. The higher nominal density of V1298 Tau e suggests it formed farther out than V1298 Tau b. The current architecture of V1298 Tau is not characterized by resonant chains. Planet-planet scattering with an outer giant planet is the most likely cause for the lack of a resonant chain between V1298 Tau's planets, but currently our search for outer companions using SPHERE and Gaia observations can exclude only the presence of planets more massive than $2 M_J$.

**Conclusions.** The most plausible scenario for V1298 Tau's formation is that the system formed by convergent migration and resonant trapping of planets born in a compact and plausibly massive disk. In the wake of their migration, V1298 Tau b and e would have left a dynamically excited protoplanetary disk, naturally creating the conditions for the later breaking of the resonant chain by planet-planet scattering.

**Key words.** planetary systems – stars: individual: V1298 Tau – planets and satellites: formation – planets and satellites: dynamical evolution and stability – planets and satellites: detection – chaos

## 1. Introduction

Young planets offer us the unique opportunity to study unaltered products of planet formation, before secular evolution and interactions with host stars modify or cancel the original characteristics of the planets. When they belong to multiplanet systems, young planets also represent invaluable case studies to comparatively investigate how different formation histories can shape planets born from the same stellar and disk environments.

Multiple authors (e.g., Donati et al. 2016; Yu et al. 2017; Benatti et al. 2021) have argued that the population of massive planets in close orbits around young stars can be significantly larger than that of their counterparts around older stars. The proposed roots of this difference are the mass loss in the early

evolutionary stages driven by the strong X and UV radiation from the host star, as well as the processes of orbital migration and chaotic evolution at play during the early phases of planetary systems that favor the engulfment (e.g., Donati et al. 2017; Spina et al. 2021) or the removal (e.g., Zinzi & Turrini 2017; Turrini et al. 2018, 2020, 2022) of planets.

Because young stars are usually very active, another reason for this difference could reside in the false positives induced by stellar activity. The number of known young planets has been questioned by recent works that have unveiled fake detections (e.g., Carleo et al. 2018; Donati et al. 2020; Damasso et al. 2020). In the past years, however, NASA *Kepler* and its K2 campaigns (Borucki et al. 2003; Howell et al. 2014) and recently the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) detected transits of several young and multiple planets (e.g., David et al. 2016; Mann et al. 2016; Rizzuto et al. 2020; Plavchan et al. 2020), allowing the study of the early evolution and architecture of these systems (e.g., Benatti et al. 2019, 2021).

<footnote>* Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundaciòn Galileo Galilei (FGG) of the Istituto Nazionale di Astrofisica (INAF) at the Observatorio del Roque de los Muchachos (La Palma, Canary Islands, Spain).</footnote>
Analyzing the *Kepler* light curve of the young star V1298 Tau (K2 Mission; Howell et al. 2014), David et al. (2019a) first discovered a Jupiter-sized planet (V1298 Tau b), and soon after David et al. (2019b) added three other planets following a forward analysis of the K2 Campaign 4 photometry. The transits of all four planets were later observed with TESS (Feinstein et al. 2022), confirming that V1298 Tau's system is host to two Neptune-sized planets (dubbed "c" and "d") and two Jovian-sized planets ("b" and "e") – names are listed in order of distance from the parent star. The four planets inhabit an orbital region comparable to that contained within Mercury's orbit in the Solar System.

In the framework of the Global Architecture of Planetary Systems project (GAPS; Covino et al. 2013; Carleo et al. 2020) and in collaboration with other groups, V1298 Tau was later observed by means of an intensive spectroscopic campaign. The GAPS campaign attained radial velocity (RV) measurements using several high-resolution spectrographs (HARPS-N, CARMENES, SES, and HERMES; Suárez Mascareño et al. 2022) in order to constrain the masses of the four planets. These observations (Suárez Mascareño et al. 2022), combined with the constraints from *Kepler* and K2, yielded mass values for both planet b (0.64 ± 0.19 $M_1$) and planet e (1.16 ± 0.30 $M_1$), while they provided only upper limits for planets c and d (<0.24 and <0.31 $M_1$, respectively). Notwithstanding the uncertainties affecting the mass estimates, the GAPS survey has revealed that V1298 Tau is one of the most massive systems among those characterized by compact orbital architectures discovered so far.

Using the Exoplanet3T online tool (Exoplanet Analysis and 3D visualization Tool<fnref n="1" />; Zinzi et al. 2021a,b) developed by the Space Science Data Center of the Italian Space Agency, we queried the NASA Exoplanet Archive<fnref n="2" /> on the 8 May 2023 to search for massive exoplanetary systems with similarly compact architectures. We searched for systems with multiplicity greater than two around solar-type stars, orbital periods lower than 100 days, and at least one planet with a mass greater than 100 $M_⊕$ being hosted. The only systems satisfying these constraints are V1298 Tau, Kepler-46 (three planets, two giants) and Kepler-256 (four planets, one giant). Adding the further requirement that the total planetary mass is greater than 1 $M_⊕$ restricted this list to V1298 Tau and Kepler-46 only.

The availability of both planetary radii and masses for planets b and e allows for the estimation of their bulk densities. The values derived by Suárez Mascareño et al. (2022) using the planetary radii from *Kepler* and K2 are 1.2 ± 0.45 and 3.6 ± 1.6 g cm⁻³ for planets b and e, respectively. While affected by large uncertainties, these values are significantly greater than those predicted by formation theories for their age. Two possible explanations are discussed by Suárez Mascareño et al. (2022): the more rapid contraction of these young planets than predicted by interior evolution models or their extreme enrichment in heavy elements.

While the first possibility implies the need to revise our current understanding of giant planet evolution (Suárez Mascareño et al. 2022), the second explanation naturally arises from the planet formation process when giant planets undergo extensive migration within their native disk (Thorngren et al. 2016; Shibata et al. 2020; Turrini et al. 2021). The conclusions drawn by Suárez Mascareño et al. (2022) are also valid for the new measurements of the radii of V1298 Tau's planets from the

photometric observations by TESS (Feinstein et al. 2022), which update the densities of planets b and e upward and downward, respectively (roughly by about 30%).

The architecture of V1298 Tau's planets initially appeared very close to being in a resonant chain. In particular, the orbit of the fourth planet was originally assessed by Suárez Mascareño et al. (2022) to be close to completing a resonant chain, either 3:2, 2:1, and 3:2 or 3:2, 2:1, and 2:1. However, the period ratio between planet b and d is too far from a 2:1 resonant ratio, effectively excluding the possibility that the system is currently in a resonant chain (Tejada Arevalo et al. 2022). Furthermore, recent data indicate an orbital period of planet e (Feinstein et al. 2022; Damasso et al. 2023) that is not in resonance with that of planet b.

The combination of the masses and densities of planets b and e and the compact architecture of V1298 Tau points toward the system having formed by convergent migration within the circumstellar disk followed by resonant trapping. The trapping in a stable resonant configuration is particularly important for planets b and e, as they would otherwise become unstable long before reaching their current compact orbits. In this scenario, however, the current non-resonant orbits of its planets require that the system underwent a phase of instability that broke the original resonant chain despite its young age (Tejada Arevalo et al. 2022).

In this work, we jointly study the formation of V1298 Tau's planets b and e and their orbital architecture with the aim of shedding light on V1298 Tau's unusual formation history. We complement our study with new observations of V1298 Tau, searching for additional outer planets to constrain the extension of the planet-forming region in V1298 Tau's circumstellar disk. As we later show, the convergent and large-scale migration required to explain the high density values and compact orbits of V1298 Tau b and e naturally creates the conditions to evolve the planetary system into its present configuration.

This paper is organized as follows: in Sect. 2 we revise the fundamental parameters of V1298 Tau's planets based on the most up-to-date observations, while in Sect. 3 we describe the numerical algorithms used in modeling their formation, their capture in resonance, and the onset of their subsequent instability. In Sect. 4, we outline our results concerning the history of the planetary system and the conditions that can lead to its present non-resonant configuration. Section 5 is devoted to observational constraints on the presence of additional planets in outer orbits obtained with SPHERE (the Spectro-Polarimetric High-contrast Exoplanet REsearch facility at the VLT telescope), while in Sect. 6 we discuss the implications of our results and combine them in a unified picture.

## 2. V1298 Tau: The star and its planets

With its estimated age ranging between 10 and 30 Ma (Suárez Mascareño et al. 2022; Maggio et al. 2022), V1298 Tau is one of the youngest solar-type planet-host stars. As its $T_{eff}$ = 5050 ± 100 K, V1298 Tau is a K1 spectral type with iron abundance [Fe/H] = 0.10 ± 0.15 dex (Suárez Mascareño et al. 2022), and it belongs to the young stellar group 29 identified by Oh et al. (2017). The star is characterized by a high and steady activity level, causing large RV activity variations, as also testified by the value of $<\log(R'_{HK})> = -4.2397 \pm 0.0179$ determined by GAPS (Suárez Mascareño et al. 2022). We report the main stellar parameters of V1298 Tau in Table 1.

<footnote n="1">¹ https://tools.ssdc.asi.it/exoplanet/</footnote>

<footnote n="2">² https://exoplanetarchive.ipac.caltech.edu/</footnote>
Analyzing the *Kepler* light curve of the young star V1298 Tau (K2 Mission; Howell et al. 2014), David et al. (2019a) first discovered a Jupiter-sized planet (V1298 Tau b), and soon after David et al. (2019b) added three other planets following a forward analysis of the K2 Campaign 4 photometry. The transits of all four planets were later observed with TESS (Feinstein et al. 2022), confirming that V1298 Tau's system is host to two Neptune-sized planets (dubbed "c" and "d") and two Jovian-sized planets ("b" and "e") – names are listed in order of distance from the parent star. The four planets inhabit an orbital region comparable to that contained within Mercury's orbit in the Solar System.

In the framework of the Global Architecture of Planetary Systems project (GAPS; Covino et al. 2013; Carleo et al. 2020) and in collaboration with other groups, V1298 Tau was later observed by means of an intensive spectroscopic campaign. The GAPS campaign attained radial velocity (RV) measurements using several high-resolution spectrographs (HARPS-N, CARMENES, SES, and HERMES; Suárez Mascareño et al. 2022) in order to constrain the masses of the four planets. These observations (Suárez Mascareño et al. 2022), combined with the constraints from *Kepler* and K2, yielded mass values for both planet b (0.64 ± 0.19 $M_1$) and planet e (1.16 ± 0.30 $M_1$), while they provided only upper limits for planets c and d (<0.24 and <0.31 $M_1$, respectively). Notwithstanding the uncertainties affecting the mass estimates, the GAPS survey has revealed that V1298 Tau is one of the most massive systems among those characterized by compact orbital architectures discovered so far.

Using the Exoplanet3T online tool (Exoplanet Analysis and 3D visualization Tool<fnref n="1" />; Zinzi et al. 2021a,b) developed by the Space Science Data Center of the Italian Space Agency, we queried the NASA Exoplanet Archive<fnref n="2" /> on the 8 May 2023 to search for massive exoplanetary systems with similarly compact architectures. We searched for systems with multiplicity greater than two around solar-type stars, orbital periods lower than 100 days, and at least one planet with a mass greater than 100 $M_⊕$ being hosted. The only systems satisfying these constraints are V1298 Tau, Kepler-46 (three planets, two giants) and Kepler-256 (four planets, one giant). Adding the further requirement that the total planetary mass is greater than 1 $M_⊕$ restricted this list to V1298 Tau and Kepler-46 only.

The availability of both planetary radii and masses for planets b and e allows for the estimation of their bulk densities. The values derived by Suárez Mascareño et al. (2022) using the planetary radii from *Kepler* and K2 are 1.2 ± 0.45 and 3.6 ± 1.6 g cm⁻³ for planets b and e, respectively. While affected by large uncertainties, these values are significantly greater than those predicted by formation theories for their age. Two possible explanations are discussed by Suárez Mascareño et al. (2022): the more rapid contraction of these young planets than predicted by interior evolution models or their extreme enrichment in heavy elements.

While the first possibility implies the need to revise our current understanding of giant planet evolution (Suárez Mascareño et al. 2022), the second explanation naturally arises from the planet formation process when giant planets undergo extensive migration within their native disk (Thorngren et al. 2016; Shibata et al. 2020; Turrini et al. 2021). The conclusions drawn by Suárez Mascareño et al. (2022) are also valid for the new measurements of the radii of V1298 Tau's planets from the

photometric observations by TESS (Feinstein et al. 2022), which update the densities of planets b and e upward and downward, respectively (roughly by about 30%).

The architecture of V1298 Tau's planets initially appeared very close to being in a resonant chain. In particular, the orbit of the fourth planet was originally assessed by Suárez Mascareño et al. (2022) to be close to completing a resonant chain, either 3:2, 2:1, and 3:2 or 3:2, 2:1, and 2:1. However, the period ratio between planet b and d is too far from a 2:1 resonant ratio, effectively excluding the possibility that the system is currently in a resonant chain (Tejada Arevalo et al. 2022). Furthermore, recent data indicate an orbital period of planet e (Feinstein et al. 2022; Damasso et al. 2023) that is not in resonance with that of planet b.

The combination of the masses and densities of planets b and e and the compact architecture of V1298 Tau points toward the system having formed by convergent migration within the circumstellar disk followed by resonant trapping. The trapping in a stable resonant configuration is particularly important for planets b and e, as they would otherwise become unstable long before reaching their current compact orbits. In this scenario, however, the current non-resonant orbits of its planets require that the system underwent a phase of instability that broke the original resonant chain despite its young age (Tejada Arevalo et al. 2022).

In this work, we jointly study the formation of V1298 Tau's planets b and e and their orbital architecture with the aim of shedding light on V1298 Tau's unusual formation history. We complement our study with new observations of V1298 Tau, searching for additional outer planets to constrain the extension of the planet-forming region in V1298 Tau's circumstellar disk. As we later show, the convergent and large-scale migration required to explain the high density values and compact orbits of V1298 Tau b and e naturally creates the conditions to evolve the planetary system into its present configuration.

This paper is organized as follows: in Sect. 2 we revise the fundamental parameters of V1298 Tau's planets based on the most up-to-date observations, while in Sect. 3 we describe the numerical algorithms used in modeling their formation, their capture in resonance, and the onset of their subsequent instability. In Sect. 4, we outline our results concerning the history of the planetary system and the conditions that can lead to its present non-resonant configuration. Section 5 is devoted to observational constraints on the presence of additional planets in outer orbits obtained with SPHERE (the Spectro-Polarimetric High-contrast Exoplanet REsearch facility at the VLT telescope), while in Sect. 6 we discuss the implications of our results and combine them in a unified picture.

## 2. V1298 Tau: The star and its planets

With its estimated age ranging between 10 and 30 Ma (Suárez Mascareño et al. 2022; Maggio et al. 2022), V1298 Tau is one of the youngest solar-type planet-host stars. As its $T_{eff}$ = 5050 ± 100 K, V1298 Tau is a K1 spectral type with iron abundance [Fe/H] = 0.10 ± 0.15 dex (Suárez Mascareño et al. 2022), and it belongs to the young stellar group 29 identified by Oh et al. (2017). The star is characterized by a high and steady activity level, causing large RV activity variations, as also testified by the value of $<\log(R'_{HK})> = -4.2397 \pm 0.0179$ determined by GAPS (Suárez Mascareño et al. 2022). We report the main stellar parameters of V1298 Tau in Table 1.

<footnote n="1">¹ https://tools.ssdc.asi.it/exoplanet/</footnote>

<footnote n="2">² https://exoplanetarchive.ipac.caltech.edu/</footnote>
From XMM-Newton observations, Maggio et al. (2022) found that this young star has a bolometric luminosity ratio of log L<sub>X</sub>/log L<sub>Bol</sub> = -3.35<sup>+0.01</sup><sub>-0.02</sub>, confirming that it is an X-ray bright young star near the saturated emission regime observed for G-K stars. Despite the close position of the four planets to their star (see Table 1), Maggio et al. (2022) found that the two outer planets (b and e) are not affected by evaporation on timescales of billions of years. The two inner planets (c and d) are also impervious to evaporation if their masses are higher than 40 M<sub>⊕</sub>.

The analysis performed by Suárez Mascareño et al. (2022) to derive the mass of planet e attributed an orbital period of 40 days to the planet. The subsequent observations by TESS led to this value being revised upward, suggesting a most probable period of 44 days and providing additional discreet solutions for larger orbital periods characterized by decreasing probabilities (Feinstein et al. 2022). The TESS observations also led to the planetary radii of planets b and e being revised downward and upward, respectively (Feinstein et al. 2022). Following these new observations, we reevaluated the mass of planet e and then updated its density as well as that of planet b.

To reassess the mass of planet e accounting for the larger possible values of its orbital period, we reanalyzed the RV dataset from Suárez Mascareño et al. (2022), imposing the additional constraint to the priors that the orbital period of planet e should not be smaller than 44 days. Feinstein et al. (2022) showed that an orbital period shorter than this value would result in an unobserved second transit within the TESS baseline. Previously, K2 data only constrained the orbital period to values higher than 36 days. Following the results of Feinstein et al. (2022), we used a half-normal distribution starting at 44 days, with a sigma of 7 days, as the prior for the period of planet e. Apart from this additional prior, our analysis is the same as described in Suárez Mascareño et al. (2022).

We fit K2 photometry, ground-based photometry, and RV time series simultaneously. We modeled the activity signals in all datasets using a Gaussian processes regression with celerite (Foreman-Mackey et al. 2017). The K2 photometry was used to obtain information on the parameters of the transits. The ground-based photometric data is contemporary to the RV data and helped constrain the parameters of the stellar activity model. To sample the posterior distribution, we relied on nested sampling (Skilling 2004) using *dynesty* (Speagle 2020). For more details and a full description of the priors used in the global model, we refer the reader to Suárez Mascareño et al. (2022).

The priors and posteriors of the analysis are shown in the top row of Fig. 1, while the bottom row shows the derived planetary and orbital parameters of planet e. The new mass value of planet e is 1.23<sup>+0.35</sup><sub>-0.48</sub> M<sub>J</sub>, which is slightly higher yet consistent within 1-σ with the previous value from Suárez Mascareño et al. (2022). Our analysis also updates the plausible periods of V1298 Tau e and shows that the most likely values are at 45 and 46 days, with two possible lower-probability solutions at 44 and 47.7 days (see Fig. 1). While multimodal posteriors, like those characterizing the orbital period, are known to be difficult to sample with traditional MCMC methods, nested sampling methods in general, and *dynesty* in particular, have been shown to be able to sample them robustly and efficiently (Speagle 2020).

None of the solutions presented above results in a resonant coupling between planets b and e. Preliminary results from ongoing observations of the ESA mission CHEOPS exclude the solutions at 44 and 46 days and are instead compatible with the one at about 45 days (Damasso et al. 2023), which we adopted as the nominal period of V1298 Tau e in this work.

<table><caption>**Table 1.** Main stellar and planetary parameters of V1298 Tau's system.</caption><thead><tr><th>Parameters</th><th>Value</th><th>Notes</th></tr></thead><tbody><tr><td colspan="3"><strong>Stellar parameters</strong></td></tr><tr><td>R<sub>⋆</sub> (R<sub>☉</sub>)</td><td>1.33<sup>+0.04</sup><sub>-0.03</sub></td><td>(1,*)</td></tr><tr><td>M<sub>⋆</sub> (M<sub>☉</sub>)</td><td>1.17 ± 0.06</td><td>(2)</td></tr><tr><td>B (mag)</td><td>11.11 ± 0.09</td><td>(3)</td></tr><tr><td>V (mag)</td><td>10.12 ± 0.05</td><td>(3)</td></tr><tr><td>G (mag)</td><td>10.0702 ± 0.0007</td><td>(4)</td></tr><tr><td>J (mag)</td><td>8.687 ± 0.023</td><td>(5)</td></tr><tr><td>Spectral type</td><td>K1</td><td>(6)</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>5050 ± 100</td><td>(2)</td></tr><tr><td>Parallax (mas)</td><td>9.258 ± 0.020</td><td>(4)</td></tr><tr><td>Distance (pc)</td><td>108.5 ± 0.7</td><td>(4)</td></tr><tr><td>⟨log(R'<sub>HK</sub>)⟩</td><td>-4.240 ± 0.018</td><td>(7)</td></tr><tr><td>Age (Ma)</td><td>10 - 30</td><td>(2,8)</td></tr><tr><td colspan="3"><strong>Planetary parameters</strong></td></tr><tr><td colspan="3"><strong>Planet b</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>⋆</sub></td><td>0.0636 ± 0.0018</td><td>(1)</td></tr><tr><td>R<sub>P</sub> (R<sub>J</sub>)</td><td>0.84 ± 0.02</td><td>(7)</td></tr><tr><td>M<sub>P</sub> (M<sub>J</sub>)</td><td>0.64 ± 0.19</td><td>(2)</td></tr><tr><td>P<sub>Orb</sub> (d)</td><td>24.1315<sup>+0.00033</sup><sub>-0.0034</sub></td><td>(1)</td></tr><tr><td>a (au)</td><td>0.1719 ± 0.0027</td><td>(2)</td></tr><tr><td>Inclination (°)</td><td>>88.7</td><td>(2)</td></tr><tr><td>Eccentricity</td><td>0.134 ± 0.075</td><td>(2)</td></tr><tr><td colspan="3"><strong>Planet c</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>⋆</sub></td><td>0.0337 ± 0.0009</td><td>(1)</td></tr><tr><td>R<sub>P</sub> (R<sub>J</sub>)</td><td>0.45 ± 0.01</td><td>(7)</td></tr><tr><td>M<sub>P</sub> (M<sub>J</sub>)</td><td>&lt; 0.24</td><td>(2)</td></tr><tr><td>P<sub>Orb</sub> (d)</td><td>8.2438<sup>+0.0024</sup><sub>-0.0020</sub></td><td>(1)</td></tr><tr><td>a (au)</td><td>0.0841 ± 0.0013</td><td>(2)</td></tr><tr><td>Inclination (°)</td><td>>87.5</td><td>(2)</td></tr><tr><td>Eccentricity</td><td>&lt;0.30</td><td>(2)</td></tr><tr><td colspan="3"><strong>Planet d</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>⋆</sub></td><td>0.0409<sup>+0.0014</sup><sub>-0.0015</sub></td><td>(1)</td></tr><tr><td>R<sub>P</sub> (R<sub>J</sub>)</td><td>0.54 ± 0.02</td><td>(7)</td></tr><tr><td>M<sub>P</sub> (M<sub>J</sub>)</td><td>&lt;0.31</td><td>(2)</td></tr><tr><td>P<sub>Orb</sub> (d)</td><td>12.3960<sup>+0.0019</sup><sub>-0.0020</sub></td><td>(1)</td></tr><tr><td>a (au)</td><td>0.1103 ± 0.0017</td><td>(2)</td></tr><tr><td>Inclination (°)</td><td>>88.3</td><td>(2)</td></tr><tr><td>Eccentricity</td><td>&lt;0.20</td><td>(2)</td></tr><tr><td colspan="3"><strong>Planet e</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>⋆</sub></td><td>0.0664<sup>+0.0025</sup><sub>-0.0021</sub></td><td>(1)</td></tr><tr><td>R<sub>P</sub> (R<sub>J</sub>)</td><td>0.88 ± 0.03</td><td>(7)</td></tr><tr><td>M<sub>P</sub> (M<sub>J</sub>)</td><td>1.23<sup>+0.35</sup><sub>-0.48</sub></td><td>(7)</td></tr><tr><td>P<sub>Orb</sub> (d)</td><td>45.46<sup>+1.10</sup><sub>-0.23</sub></td><td>(7)</td></tr><tr><td>a (au)</td><td>0.2639<sup>+0.0062</sup><sub>-0.0057</sub></td><td>(7)</td></tr><tr><td>Inclination (°)</td><td>89.295<sup>+0.205</sup><sub>-0.091</sub></td><td>(7)</td></tr><tr><td>Eccentricity</td><td>0.263<sup>+0.052</sup><sub>-0.134</sub></td><td>(7)</td></tr></tbody></table>

**Notes.** (1) Feinstein et al. (2022); (2) Suárez Mascareño et al. (2022); (3) Høg et al. (2000); (4) Gaia Collaboration (2021); (5) Skrutskie et al. (2006); (6) Nguyen et al. (2012); (7) this work; (8) Maggio et al. (2022); (*) Suárez Mascareño et al. (2022) report a slightly smaller value of R<sub>⋆</sub> = 1.29 ± 0.07. As the two values are consistent with each other within 1 σ, we adopted the value of Feinstein et al. (2022) for consistency with their fitting of TESS light curves.
From *XMM-Newton* observations, Maggio et al. (2022) found that this young star has a bolometric luminosity ratio of $log L_X / log L_{Bol} = -3.35^{+0.01}_{-0.02}$, confirming that it is an X-ray bright young star near the saturated emission regime observed for G-K stars. Despite the close position of the four planets to their star (see Table 1), Maggio et al. (2022) found that the two outer planets (b and e) are not affected by evaporation on timescales of billions of years. The two inner planets (c and d) are also impervious to evaporation if their masses are higher than $40 M_⊕$.

The analysis performed by Suárez Mascareño et al. (2022) to derive the mass of planet e attributed an orbital period of 40 days to the planet. The subsequent observations by TESS led to this value being revised upward, suggesting a most probable period of 44 days and providing additional discreet solutions for larger orbital periods characterized by decreasing probabilities (Feinstein et al. 2022). The TESS observations also led to the planetary radii of planets b and e being revised downward and upward, respectively (Feinstein et al. 2022). Following these new observations, we reevaluated the mass of planet e and then updated its density as well as that of planet b.

To reassess the mass of planet e accounting for the larger possible values of its orbital period, we reanalyzed the RV dataset from Suárez Mascareño et al. (2022), imposing the additional constraint to the priors that the orbital period of planet e should not be smaller than 44 days. Feinstein et al. (2022) showed that an orbital period shorter than this value would result in an unobserved second transit within the TESS baseline. Previously, K2 data only constrained the orbital period to values higher than 36 days. Following the results of Feinstein et al. (2022), we used a half-normal distribution starting at 44 days, with a sigma of 7 days, as the prior for the period of planet e. Apart from this additional prior, our analysis is the same as described in Suárez Mascareño et al. (2022).

We fit K2 photometry, ground-based photometry, and RV time series simultaneously. We modeled the activity signals in all datasets using a Gaussian processes regression with celerite (Foreman-Mackey et al. 2017). The K2 photometry was used to obtain information on the parameters of the transits. The ground-based photometric data is contemporary to the RV data and helped constrain the parameters of the stellar activity model. To sample the posterior distribution, we relied on nested sampling (Skilling 2004) using *dynesty* (Speagle 2020). For more details and a full description of the priors used in the global model, we refer the reader to Suárez Mascareño et al. (2022).

The priors and posteriors of the analysis are shown in the top row of Fig. 1, while the bottom row shows the derived planetary and orbital parameters of planet e. The new mass value of planet e is $1.23^{+0.35}_{-0.48} M_J$, which is slightly higher yet consistent within $1-\sigma$ with the previous value from Suárez Mascareño et al. (2022). Our analysis also updates the plausible periods of V1298 Tau e and shows that the most likely values are at 45 and 46 days, with two possible lower-probability solutions at 44 and 47.7 days (see Fig. 1). While multimodal posteriors, like those characterizing the orbital period, are known to be difficult to sample with traditional MCMC methods, nested sampling methods in general, and *dynesty* in particular, have been shown to be able to sample them robustly and efficiently (Speagle 2020).

None of the solutions presented above results in a resonant coupling between planets b and e. Preliminary results from ongoing observations of the ESA mission CHEOPS exclude the solutions at 44 and 46 days and are instead compatible with the one at about 45 days (Damasso et al. 2023), which we adopted as the nominal period of V1298 Tau e in this work.

<table><caption>**Table 1.** Main stellar and planetary parameters of V1298 Tau's system.</caption><thead><tr><th>Parameters</th><th>Value</th><th>Notes</th></tr></thead><tbody><tr><td colspan="3"><strong>Stellar parameters</strong></td></tr><tr><td>R<sub>⋆</sub> (R<sub>☉</sub>)</td><td>1.33<sup>+0.04</sup><sub>-0.03</sub></td><td>(1,*)</td></tr><tr><td>M<sub>⋆</sub> (M<sub>☉</sub>)</td><td>1.17 ± 0.06</td><td>(2)</td></tr><tr><td>B (mag)</td><td>11.11 ± 0.09</td><td>(3)</td></tr><tr><td>V (mag)</td><td>10.12 ± 0.05</td><td>(3)</td></tr><tr><td>G (mag)</td><td>10.0702 ± 0.0007</td><td>(4)</td></tr><tr><td>J (mag)</td><td>8.687 ± 0.023</td><td>(5)</td></tr><tr><td>Spectral type</td><td>K1</td><td>(6)</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>5050 ± 100</td><td>(2)</td></tr><tr><td>Parallax (mas)</td><td>9.258 ± 0.020</td><td>(4)</td></tr><tr><td>Distance (pc)</td><td>108.5 ± 0.7</td><td>(4)</td></tr><tr><td>⟨log(R'<sub>HK</sub>)⟩</td><td>-4.240 ± 0.018</td><td>(7)</td></tr><tr><td>Age (Ma)</td><td>10 - 30</td><td>(2,8)</td></tr><tr><td colspan="3"><strong>Planetary parameters</strong></td></tr><tr><td colspan="3"><strong>Planet b</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>⋆</sub></td><td>0.0636 ± 0.0018</td><td>(1)</td></tr><tr><td>R<sub>P</sub> (R<sub>J</sub>)</td><td>0.84 ± 0.02</td><td>(7)</td></tr><tr><td>M<sub>P</sub> (M<sub>J</sub>)</td><td>0.64 ± 0.19</td><td>(2)</td></tr><tr><td>P<sub>Orb</sub> (d)</td><td>24.1315<sup>+0.00033</sup><sub>-0.0034</sub></td><td>(1)</td></tr><tr><td>a (au)</td><td>0.1719 ± 0.0027</td><td>(2)</td></tr><tr><td>Inclination (°)</td><td>>88.7</td><td>(2)</td></tr><tr><td>Eccentricity</td><td>0.134 ± 0.075</td><td>(2)</td></tr><tr><td colspan="3"><strong>Planet c</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>⋆</sub></td><td>0.0337 ± 0.0009</td><td>(1)</td></tr><tr><td>R<sub>P</sub> (R<sub>J</sub>)</td><td>0.45 ± 0.01</td><td>(7)</td></tr><tr><td>M<sub>P</sub> (M<sub>J</sub>)</td><td>&lt; 0.24</td><td>(2)</td></tr><tr><td>P<sub>Orb</sub> (d)</td><td>8.2438<sup>+0.0024</sup><sub>-0.0020</sub></td><td>(1)</td></tr><tr><td>a (au)</td><td>0.0841 ± 0.0013</td><td>(2)</td></tr><tr><td>Inclination (°)</td><td>>87.5</td><td>(2)</td></tr><tr><td>Eccentricity</td><td>&lt;0.30</td><td>(2)</td></tr><tr><td colspan="3"><strong>Planet d</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>⋆</sub></td><td>0.0409<sup>+0.0014</sup><sub>-0.0015</sub></td><td>(1)</td></tr><tr><td>R<sub>P</sub> (R<sub>J</sub>)</td><td>0.54 ± 0.02</td><td>(7)</td></tr><tr><td>M<sub>P</sub> (M<sub>J</sub>)</td><td>&lt;0.31</td><td>(2)</td></tr><tr><td>P<sub>Orb</sub> (d)</td><td>12.3960<sup>+0.0019</sup><sub>-0.0020</sub></td><td>(1)</td></tr><tr><td>a (au)</td><td>0.1103 ± 0.0017</td><td>(2)</td></tr><tr><td>Inclination (°)</td><td>>88.3</td><td>(2)</td></tr><tr><td>Eccentricity</td><td>&lt;0.20</td><td>(2)</td></tr><tr><td colspan="3"><strong>Planet e</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>⋆</sub></td><td>0.0664<sup>+0.0025</sup><sub>-0.0021</sub></td><td>(1)</td></tr><tr><td>R<sub>P</sub> (R<sub>J</sub>)</td><td>0.88 ± 0.03</td><td>(7)</td></tr><tr><td>M<sub>P</sub> (M<sub>J</sub>)</td><td>1.23<sup>+0.35</sup><sub>-0.48</sub></td><td>(7)</td></tr><tr><td>P<sub>Orb</sub> (d)</td><td>45.46<sup>+1.10</sup><sub>-0.23</sub></td><td>(7)</td></tr><tr><td>a (au)</td><td>0.2639<sup>+0.0062</sup><sub>-0.0057</sub></td><td>(7)</td></tr><tr><td>Inclination (°)</td><td>89.295<sup>+0.205</sup><sub>-0.091</sub></td><td>(7)</td></tr><tr><td>Eccentricity</td><td>0.263<sup>+0.052</sup><sub>-0.134</sub></td><td>(7)</td></tr></tbody></table>

**Notes.** (1) Feinstein et al. (2022); (2) Suárez Mascareño et al. (2022); (3) Høg et al. (2000); (4) Gaia Collaboration (2021); (5) Skrutskie et al. (2006); (6) Nguyen et al. (2012); (7) this work; (8) Maggio et al. (2022); (*) Suárez Mascareño et al. (2022) report a slightly smaller value of $R_*$ = $1.29 \pm 0.07$. As the two values are consistent with each other within $1 \sigma$, we adopted the value of Feinstein et al. (2022) for consistency with their fitting of TESS light curves.
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** New analysis of the RV dataset and its results on planet e. Top: priors and posteriors of the new analysis from Suárez Mascareño et al. (2022) where we include the constraints from TESS on its orbital period (Feinstein et al. 2022). Bottom: derived distributions of the planetary and orbital parameters. We note that we used the updated planetary radius from TESS (Feinstein et al. 2022) in our analysis in place of the one appearing in this figure, which is still based on the *Kepler* and K2 data, as in Suárez Mascareño et al. (2022).</figcaption></figure>

The full updated characterization of V1298 Tau's four planets is reported in Table 1. Figure 1 and Table 1 show that, like its mass, the eccentricity of planet e is higher than that estimated by Suárez Mascareño et al. (2022), but also in this case, the two values are consistent at the 1-σ level.

## 2.1. Dynamical excitation of V1298 Tau b and e

We used the updated architecture of V1298 Tau from Table 1 to investigate the dynamical excitation of planets b and e by means of their normalized angular momentum deficit (NAMD; see Chambers 2001; Turrini et al. 2020). The NAMD provides an architecture-agnostic measure of dynamical excitation that can intuitively be interpreted as the “dynamical temperature” of planetary systems. Stated simply, the higher the value, the more excited the dynamical state of the system. We adopted the NAMD of the Solar System (1.3 × 10⁻<fnref n="3" />; Turrini et al. 2020) as the boundary between dynamically cold and hot orbits (see Carleo et al. 2021; Turrini et al. 2022 for a discussion).

Following Carleo et al. (2021), the goal of our analysis is not to pinpoint the exact value of V1298 Tau's NAMD but to verify if, notwithstanding the uncertainties on the physical and orbital parameters of its planets, its dynamical excitation is systematically higher than that of the Solar System and whether it indicates past or current phases of dynamical instability. We accounted for the uncertainties on the mass, semimajor axis, and eccentricity of the two planets following the Monte Carlo approach described in Turrini et al. (2020) and Carleo et al. (2021). To better sample the effects of the large uncertainties of these parameters on the NAMD, we used 10⁶ extractions for each parameter in place of the original 10⁴ proposed by Laskar & Petit (2017).

Based on the low mutual inclination of the two planets (<1°; see Table 1) that limitedly contributes to the NAMD, we assumed the two orbits to be coplanar. Due to the positive-defined nature of the NAMD, this choice means that we err toward lower excitations. The resulting NAMD modal value (2.2 × 10⁻²) is 17 times higher than that of the Solar System, with the 3σ range extending from five times (6.4 × 10⁻<fnref n="3" />) to 120 times

(0.16) that of the Solar System. As discussed by Turrini et al. (2020, 2022) and illustrated by the recent study by Rickman et al. (2023), such a systematically high range of NAMD values is associated with phases of dynamical instability and planet-planet scattering events (either in the past or presently).

The results of the NAMD analysis quantitatively confirm the dynamical indication supplied by the present non-resonant state of this compact system (Tejada Arevalo et al. 2022). Specifically, the dynamical state of V1298 Tau argues in favour of a formation scenario where the planetary system acquired its compact architecture by forming through convergent migration and resonant capture. The original resonant architecture was then broken during a later phase of dynamical instability. This scenario is supported by dynamical population studies of exoplanetary systems (Laskar & Petit 2017; Gajdoš & Vaňko 2023) revealing how 25-45% of known multiplanet systems, and the majority of high-multiplicity systems hosting four or more planets like V1298 Tau, show the signatures of chaos and instability in their architectures. We investigate the plausible causes of the instability in Sects. 4.4 and 4.5.

## 2.2. Density values of V1298 Tau b and e

We reevaluated the densities of planets b and e using the new mass estimate of planet e from Sect. 2 and the planetary radii provided by TESS (Feinstein et al. 2022). Specifically, we used the values of $R_p/R_*$ and $R_*$ from Feinstein et al. (2022) reported in Table 1 together with the volumetric mean radii of the Sun and Jupiter<fnref n="3" /> to estimate the planetary radii (see Table 1) and volumes. The resulting density values of the planets are 1.4 ± 0.4 g cm⁻<fnref n="3" /> and 2.4 ± 0.8 g cm⁻<fnref n="3" />, respectively. As discussed by Suárez Mascareño et al. (2022), unless the two giant planets underwent a more rapid contraction than predicted by interior evolution models, their densities can be explained by marked enrichments in heavy elements (Thorngren et al. 2016).

<footnote n="3">³ https://nssdc.gsfc.nasa.gov/planetary/planetfact.html</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** New analysis of the RV dataset and its results on planet e. Top: priors and posteriors of the new analysis from Suárez Mascareño et al. (2022) where we include the constraints from TESS on its orbital period (Feinstein et al. 2022). Bottom: derived distributions of the planetary and orbital parameters. We note that we used the updated planetary radius from TESS (Feinstein et al. 2022) in our analysis in place of the one appearing in this figure, which is still based on the *Kepler* and K2 data, as in Suárez Mascareño et al. (2022).</figcaption></figure>

The full updated characterization of V1298 Tau's four planets is reported in Table 1. Figure 1 and Table 1 show that, like its mass, the eccentricity of planet e is higher than that estimated by Suárez Mascareño et al. (2022), but also in this case, the two values are consistent at the 1-σ level.

## 2.1. Dynamical excitation of V1298 Tau b and e

We used the updated architecture of V1298 Tau from Table 1 to investigate the dynamical excitation of planets b and e by means of their normalized angular momentum deficit (NAMD; see Chambers 2001; Turrini et al. 2020). The NAMD provides an architecture-agnostic measure of dynamical excitation that can intuitively be interpreted as the “dynamical temperature” of planetary systems. Stated simply, the higher the value, the more excited the dynamical state of the system. We adopted the NAMD of the Solar System (1.3 × 10⁻<fnref n="3" />; Turrini et al. 2020) as the boundary between dynamically cold and hot orbits (see Carleo et al. 2021; Turrini et al. 2022 for a discussion).

Following Carleo et al. (2021), the goal of our analysis is not to pinpoint the exact value of V1298 Tau's NAMD but to verify if, notwithstanding the uncertainties on the physical and orbital parameters of its planets, its dynamical excitation is systematically higher than that of the Solar System and whether it indicates past or current phases of dynamical instability. We accounted for the uncertainties on the mass, semimajor axis, and eccentricity of the two planets following the Monte Carlo approach described in Turrini et al. (2020) and Carleo et al. (2021). To better sample the effects of the large uncertainties of these parameters on the NAMD, we used 10⁶ extractions for each parameter in place of the original 10⁴ proposed by Laskar & Petit (2017).

Based on the low mutual inclination of the two planets (<1°; see Table 1) that limitedly contributes to the NAMD, we assumed the two orbits to be coplanar. Due to the positive-defined nature of the NAMD, this choice means that we err toward lower excitations. The resulting NAMD modal value (2.2 × 10⁻²) is 17 times higher than that of the Solar System, with the 3σ range extending from five times (6.4 × 10⁻<fnref n="3" />) to 120 times

(0.16) that of the Solar System. As discussed by Turrini et al. (2020, 2022) and illustrated by the recent study by Rickman et al. (2023), such a systematically high range of NAMD values is associated with phases of dynamical instability and planet-planet scattering events (either in the past or presently).

The results of the NAMD analysis quantitatively confirm the dynamical indication supplied by the present non-resonant state of this compact system (Tejada Arevalo et al. 2022). Specifically, the dynamical state of V1298 Tau argues in favour of a formation scenario where the planetary system acquired its compact architecture by forming through convergent migration and resonant capture. The original resonant architecture was then broken during a later phase of dynamical instability. This scenario is supported by dynamical population studies of exoplanetary systems (Laskar & Petit 2017; Gajdoš & Vaňko 2023) revealing how 25-45% of known multiplanet systems, and the majority of high-multiplicity systems hosting four or more planets like V1298 Tau, show the signatures of chaos and instability in their architectures. We investigate the plausible causes of the instability in Sects. 4.4 and 4.5.

## 2.2. Density values of V1298 Tau b and e

We reevaluated the densities of planets b and e using the new mass estimate of planet e from Sect. 2 and the planetary radii provided by TESS (Feinstein et al. 2022). Specifically, we used the values of $R_p/R_*$ and $R_*$ from Feinstein et al. (2022) reported in Table 1 together with the volumetric mean radii of the Sun and Jupiter<fnref n="3" /> to estimate the planetary radii (see Table 1) and volumes. The resulting density values of the planets are 1.4 ± 0.4 g cm⁻<fnref n="3" /> and 2.4 ± 0.8 g cm⁻<fnref n="3" />, respectively. As discussed by Suárez Mascareño et al. (2022), unless the two giant planets underwent a more rapid contraction than predicted by interior evolution models, their densities can be explained by marked enrichments in heavy elements (Thorngren et al. 2016).

<footnote n="3">³ https://nssdc.gsfc.nasa.gov/planetary/planetfact.html</footnote>
Recent results on the evolution over time of the dust abundance in circumstellar disks (Manara et al. 2018; Mulders et al. 2021; Bernabò et al. 2022) and the radiometric ages of meteorites in the Solar System (see Scott 2007; Coradini et al. 2011; Lichtenberg et al. 2023, and references therein) argue that the bulk of the heavy elements is locked into planetesimals by the time giant planets form. As discussed by Shibata et al. (2020) and Turrini et al. (2021), the mass of planetesimals that giant planets accrete during their growth is directly proportional to the extent of their migration. The larger the migration, the greater the mass of planetesimals that can enter their feeding zone and be accreted. Conversely, giant planets undergoing little or no migration will experience limited accretion of planetesimals (Turrini et al. 2015; Shibata & Ikoma 2019).

We used the updated density values of the two giant planets to derive order-of-magnitude estimates of the masses of planetesimals that they need to accrete. We focused on the density of the solid material and not on the bulk density of the planetesimals in order to remove the issue of the unknown macroporosity of the planetesimals. We assumed the solid material accreted through the planetesimals to possess density of $3 \text{ g cm}^{-3}$, that is, to be composed half of rock and metals with an average density of $5 \text{ g cm}^{-3}$ and half of ice with a density of $1 \text{ g cm}^{-3}$ (see e.g., Turrini et al. 2021; Pacetti et al. 2022, for the mass balance between rock and ice).

For the gas composing the envelopes of V1298 Tau b and e, we assumed the same metallicity and density as that of Jupiter in the Solar System. This gas has a density of $1.33 \text{ g cm}^{-3}$, and its composition is about three times richer in heavy elements with respect to hydrogen than that of the Sun (see Atreya et al. 2018; Öberg & Wordsworth 2019, and references therein). In other words, heavy elements account for at least 4% of Jupiter's mass. As in the enrichment scenario discussed by Suárez Mascareño et al. (2022) for planets b and e, Jupiter's formation, metallicity, and enrichment in heavy elements are argued to have been shaped by large-scale migration and planetesimal accretion (Pirani et al. 2019; Öberg & Wordsworth 2019).

The density of V1298 Tau b can be obtained by adding about 11 $M_⊕$ of planetesimals to about 192 $M_⊕$ of Jovian gas. Merging the contributions in heavy elements of planetesimals and enriched Jovian gas, the resulting mixture is composed of 91% (184 $M_⊕$) H and He and 9% (19 $M_⊕$) heavy elements. The same approach applied to V1298 Tau e required combining 140 $M_⊕$ of Jovian gas and 250 $M_⊕$ of planetesimals. Grouping the heavy elements together resulted in a mixture where 66% of the mass (256 $M_⊕$) is provided by heavy elements and 34% (135 $M_⊕$) is by H and He. While V1298 Tau b is similar to Jupiter and Saturn in terms of the ratio between heavy elements and hydrogen, V1298 Tau e appears closer to a Neptunian planet, as H and He do not dominate its mass.

Since this metallicity is anomalously high for such a massive planet (Thorngren et al. 2016) and the uncertainty toward lower mass values is particularly large, we explored the case where the real mass of V1298 Tau e is 0.75 $M_⊕$ (i.e., $1-\sigma$ less than the nominal value from Table 1 and Fig. 1). This resulted in a planetary density of $1.5 \text{ g cm}^{-3}$ and required adding 29 $M_⊕$ of heavy elements to 209 $M_⊕$ of H and He. In this scenario, the metallicities of V1298 Tau b and e are similar, and both giant planets are consistent with the ratio between heavy elements and hydrogen of Jupiter and Saturn.

In the following analysis, we consider both scenarios summarized in Table 2, that is, one where V1298 Tau e is characterized by a higher mass and metallicity (HMZe in the following and in Table 2) and one where the giant planet has a lower mass

<table><caption>**Table 2.** Physical scenarios considered for the planetary masses and the enrichment in heavy elements of V1298 Tau b and V1298 Tau e.</caption><thead><tr><th rowspan="2">Scenario Id.</th><th colspan="2">V1298 Tau b</th><th colspan="2">V1298 Tau e</th></tr><tr><th>Mass (M<sub>J</sub>)</th><th>Z (M<sub>⊕</sub>)</th><th>Mass (M<sub>J</sub>)</th><th>Z (M<sub>⊕</sub>)</th></tr></thead><tbody><tr><td>LMZe</td><td>0.64</td><td>19</td><td>0.75</td><td>29</td></tr><tr><td>HMZe</td><td>0.64</td><td>19</td><td>1.23</td><td>256</td></tr></tbody></table>

and metallicity (LMZe in the following and in Table 2). Before proceeding, we note that the gas density we adopted in our back-of-the-envelope computations is the one characterizing Jupiter today. Between 10–30 Ma after its formation, the radius of Jupiter is expected to have been 1.5–1.6 times the present value (Lissauer et al. 2009; D’Angelo et al. 2021), resulting in its larger volume by a factor of three to four and correspondingly to a lower density of the gas (about 0.3–0.4 g cm$^{-3}$).

This lower gas density requires larger amounts of heavy elements to fit the estimated densities, specifically 87 $M_⊕$ for V1298 Tau b and 108 $M_⊕$ for V1298 Tau e, even in the LMZe scenario. In this case, the mass of the two giant planets would be composed of 45% heavy elements ($Z/Z_* \sim 25$ using the nominal metallicity of V1298 Tau from Sect. 2). In the HMZe scenario, V1298 Tau e would require about 305 $M_⊕$ of heavy elements, that is, H and He would supply only 20% of its mass ($Z/Z_* \sim 45$). Given the uncertainty affecting the masses of the two planets and its significant impact on their densities and metallicities (see Table 2), we did not also simulate these scenarios, but we discuss their implications in Sects. 4.1 and 6.

## 3. V1298 Tau's formation and dynamical histories: Numerical methods

### 3.1. Formation simulations of V1298 Tau b and e

We simulated the growth of V1298 Tau b and e, their migration, and their interactions with the planetesimal disk using the N-body code Mercury-Arχes (Turrini et al. 2019, 2021). The N-body simulations modeled the effects of the mass growth of the forming V1298 Tau b and e planets, their planetary radius evolution, and orbital migration as well as the dynamical evolution of their surrounding planetesimal disk under the effects of their gravitational perturbations alongside those of gas drag and the disk gravity.

We modeled the native circumstellar disk as possessing the characteristic radii $r_c = 50$ AU and gas surface density $\Sigma(r) = \Sigma_0 (r/r_c)^\gamma \exp[-(r/r_c)^{(2-\gamma)}]$ (see Sect. 4.1 for the values of $\Sigma_0$ considered in the simulations), where $\gamma = 0.8$ (Isella et al. 2016). We assumed that the disk gas is in a steady state and does not decline over time. The disk temperature profile on the midplane is $T(r) = T_0 r^{-0.6}$, where $T_0 = 200$ K (Andrews & Williams 2007; Öberg et al. 2011; Eistrup et al. 2016). The mass of V1298 Tau was set to 1.17 $M_⊕$ (Suárez Mascareño et al. 2022).

Planetesimals were included in the N-body simulations as particles possessing an inertial mass, computed assuming a common diameter of 100 km (see Klahr & Schreiber 2016; Johansen & Lambrechts 2017; Turrini et al. 2019), a bulk density of $1 \text{ g cm}^{-3}$ (see Turrini et al. 2019, 2021), and no gravitational mass. The density of the planetesimals is lower than what was used in Sect. 2 to estimate the amounts of solid material accreted by the giant planets to account for the macro-porosity of these planetary bodies. The dynamical evolution of the planetesimals
Recent results on the evolution over time of the dust abundance in circumstellar disks (Manara et al. 2018; Mulders et al. 2021; Bernabò et al. 2022) and the radiometric ages of meteorites in the Solar System (see Scott 2007; Coradini et al. 2011; Lichtenberg et al. 2023, and references therein) argue that the bulk of the heavy elements is locked into planetesimals by the time giant planets form. As discussed by Shibata et al. (2020) and Turrini et al. (2021), the mass of planetesimals that giant planets accrete during their growth is directly proportional to the extent of their migration. The larger the migration, the greater the mass of planetesimals that can enter their feeding zone and be accreted. Conversely, giant planets undergoing little or no migration will experience limited accretion of planetesimals (Turrini et al. 2015; Shibata & Ikoma 2019).

We used the updated density values of the two giant planets to derive order-of-magnitude estimates of the masses of planetesimals that they need to accrete. We focused on the density of the solid material and not on the bulk density of the planetesimals in order to remove the issue of the unknown macroporosity of the planetesimals. We assumed the solid material accreted through the planetesimals to possess density of $3 \text{ g cm}^{-3}$, that is, to be composed half of rock and metals with an average density of $5 \text{ g cm}^{-3}$ and half of ice with a density of $1 \text{ g cm}^{-3}$ (see e.g., Turrini et al. 2021; Pacetti et al. 2022, for the mass balance between rock and ice).

For the gas composing the envelopes of V1298 Tau b and e, we assumed the same metallicity and density as that of Jupiter in the Solar System. This gas has a density of $1.33 \text{ g cm}^{-3}$, and its composition is about three times richer in heavy elements with respect to hydrogen than that of the Sun (see Atreya et al. 2018; Öberg & Wordsworth 2019, and references therein). In other words, heavy elements account for at least 4% of Jupiter's mass. As in the enrichment scenario discussed by Suárez Mascareño et al. (2022) for planets b and e, Jupiter's formation, metallicity, and enrichment in heavy elements are argued to have been shaped by large-scale migration and planetesimal accretion (Pirani et al. 2019; Öberg & Wordsworth 2019).

The density of V1298 Tau b can be obtained by adding about 11 $M_⊕$ of planetesimals to about 192 $M_⊕$ of Jovian gas. Merging the contributions in heavy elements of planetesimals and enriched Jovian gas, the resulting mixture is composed of 91% (184 $M_⊕$) H and He and 9% (19 $M_⊕$) heavy elements. The same approach applied to V1298 Tau e required combining 140 $M_⊕$ of Jovian gas and 250 $M_⊕$ of planetesimals. Grouping the heavy elements together resulted in a mixture where 66% of the mass (256 $M_⊕$) is provided by heavy elements and 34% (135 $M_⊕$) is by H and He. While V1298 Tau b is similar to Jupiter and Saturn in terms of the ratio between heavy elements and hydrogen, V1298 Tau e appears closer to a Neptunian planet, as H and He do not dominate its mass.

Since this metallicity is anomalously high for such a massive planet (Thorngren et al. 2016) and the uncertainty toward lower mass values is particularly large, we explored the case where the real mass of V1298 Tau e is 0.75 $M_⊕$ (i.e., $1-\sigma$ less than the nominal value from Table 1 and Fig. 1). This resulted in a planetary density of $1.5 \text{ g cm}^{-3}$ and required adding 29 $M_⊕$ of heavy elements to 209 $M_⊕$ of H and He. In this scenario, the metallicities of V1298 Tau b and e are similar, and both giant planets are consistent with the ratio between heavy elements and hydrogen of Jupiter and Saturn.

In the following analysis, we consider both scenarios summarized in Table 2, that is, one where V1298 Tau e is characterized by a higher mass and metallicity (HMZe in the following and in Table 2) and one where the giant planet has a lower mass

<table><caption>**Table 2.** Physical scenarios considered for the planetary masses and the enrichment in heavy elements of V1298 Tau b and V1298 Tau e.</caption><thead><tr><th rowspan="2">Scenario Id.</th><th colspan="2">V1298 Tau b</th><th colspan="2">V1298 Tau e</th></tr><tr><th>Mass (M<sub>J</sub>)</th><th>Z (M<sub>⊕</sub>)</th><th>Mass (M<sub>J</sub>)</th><th>Z (M<sub>⊕</sub>)</th></tr></thead><tbody><tr><td>LMZe</td><td>0.64</td><td>19</td><td>0.75</td><td>29</td></tr><tr><td>HMZe</td><td>0.64</td><td>19</td><td>1.23</td><td>256</td></tr></tbody></table>

and metallicity (LMZe in the following and in Table 2). Before proceeding, we note that the gas density we adopted in our back-of-the-envelope computations is the one characterizing Jupiter today. Between 10–30 Ma after its formation, the radius of Jupiter is expected to have been 1.5–1.6 times the present value (Lissauer et al. 2009; D’Angelo et al. 2021), resulting in its larger volume by a factor of three to four and correspondingly to a lower density of the gas (about 0.3–0.4 g cm$^{-3}$).

This lower gas density requires larger amounts of heavy elements to fit the estimated densities, specifically 87 $M_⊕$ for V1298 Tau b and 108 $M_⊕$ for V1298 Tau e, even in the LMZe scenario. In this case, the mass of the two giant planets would be composed of 45% heavy elements ($Z/Z_* \sim 25$ using the nominal metallicity of V1298 Tau from Sect. 2). In the HMZe scenario, V1298 Tau e would require about 305 $M_⊕$ of heavy elements, that is, H and He would supply only 20% of its mass ($Z/Z_* \sim 45$). Given the uncertainty affecting the masses of the two planets and its significant impact on their densities and metallicities (see Table 2), we did not also simulate these scenarios, but we discuss their implications in Sects. 4.1 and 6.

## 3. V1298 Tau's formation and dynamical histories: Numerical methods

### 3.1. Formation simulations of V1298 Tau b and e

We simulated the growth of V1298 Tau b and e, their migration, and their interactions with the planetesimal disk using the N-body code Mercury-Arχes (Turrini et al. 2019, 2021). The N-body simulations modeled the effects of the mass growth of the forming V1298 Tau b and e planets, their planetary radius evolution, and orbital migration as well as the dynamical evolution of their surrounding planetesimal disk under the effects of their gravitational perturbations alongside those of gas drag and the disk gravity.

We modeled the native circumstellar disk as possessing the characteristic radii $r_c = 50$ AU and gas surface density $\Sigma(r) = \Sigma_0 (r/r_c)^\gamma \exp[-(r/r_c)^{(2-\gamma)}]$ (see Sect. 4.1 for the values of $\Sigma_0$ considered in the simulations), where $\gamma = 0.8$ (Isella et al. 2016). We assumed that the disk gas is in a steady state and does not decline over time. The disk temperature profile on the midplane is $T(r) = T_0 r^{-0.6}$, where $T_0 = 200$ K (Andrews & Williams 2007; Öberg et al. 2011; Eistrup et al. 2016). The mass of V1298 Tau was set to $1.17 M_⊕$ (Suárez Mascareño et al. 2022).

Planetesimals were included in the N-body simulations as particles possessing an inertial mass, computed assuming a common diameter of 100 km (see Klahr & Schreiber 2016; Johansen & Lambrechts 2017; Turrini et al. 2019), a bulk density of $1 \text{ g cm}^{-3}$ (see Turrini et al. 2019, 2021), and no gravitational mass. The density of the planetesimals is lower than what was used in Sect. 2 to estimate the amounts of solid material accreted by the giant planets to account for the macro-porosity of these planetary bodies. The dynamical evolution of the planetesimals
is affected by the gravity of the host star and of the forming planets V1298 Tau b and e, and by the disk gas through aerodynamic drag and gravity (as they possess inertial mass).

The dynamical evolution of the planetesimals is not affected by the interactions among planetesimals themselves (as they do not possess gravitational mass) nor do the planetesimals perturb the two forming planets. The damping effects of gas drag on the planetesimals were simulated following the treatment from Brasser et al. (2007) with updated drag coefficients from Nagasawa et al. (2019) accounting for both the Mach and Reynolds numbers of the planetesimals. The exciting effects of the disk self-gravity were simulated based on the analytical treatment for axisymmetric disks by Ward (1981) following Marzari (2018) and Nagasawa et al. (2019).

The formation of the two giant planets was modeled over two growth and migration phases using the parametric approach from Turrini et al. (2019, 2021). The first phase accounts for their core growth and subsequent capture of an expanded atmosphere (e.g., Bitsch et al. 2015; Johansen et al. 2019; D'Angelo et al. 2021). The planetary mass evolved as $M_p(t) = M_0 + (e/e-1)(M_1 - M_0)(1 - e^{-t/\tau_p})$, where $M_0 = 0.1 M_\oplus$ is the initial mass of the core, $M_1 = 30 M_\oplus$ is the final cumulative mass of core and expanded atmosphere at the end of the first growth phase (see Turrini et al. 2021 for further discussion), and e is the Euler number.

The constant $\tau_p$ is the duration of the first growth phase and was set to about 1 Ma for both planets based on observational and theoretical constraints from circumstellar disks (Manara et al. 2018; Mulders et al. 2021; Bernabò et al. 2022) and the Solar System (Scott 2007; Coradini et al. 2011; Lichtenberg et al. 2023). The individual values of $\tau_p$ are 1 Ma and 1.25 Ma for planets b and e, respectively. The slower growth of planet e was introduced to delay the onset of its second growth and migration phase and to allow planet b to get close to its final orbit before planet e reaches its peak rate of migration. This choice limits the chances of destabilizing close encounters between the two giant planets before they achieve their compact architecture.

The second phase of mass growth accounted for the runaway gas accretion of the two giant planets, where their mass evolves as $M_p(t) = M_1 + (M_2 - M_1)(1 - e^{-(t-\tau_p)/\tau_g})$, with $M_2$ as the final mass of the giant planets and $\tau_g$ as the e-folding time of the runaway gas accretion process. Based on the discussion in Sect. 2.2, $M_2$ was set to 184 $M_\oplus$ for planet b when aiming to reproduce both the HMZe and LMZe scenarios from Table 2. In the case of planet e, $M_2$ was set to 135 $M_\oplus$ and 209 $M_\oplus$ when the simulations focus on the HMZe and LMZe scenarios, respectively.

The value of $\tau_g$ was set to 0.1 Ma based on the results of hydrodynamic simulations (Lissauer et al. 2009; D'Angelo et al. 2021), meaning that the gas giants reach more than 99% of their final mass in about 0.5 Ma from the onset of the runaway gas accretion. During the runaway gas accretion, giant planets form a gap in the disk gas, whose width was modeled as $W_{gap} = C \cdot R_H$ (Isella et al. 2016; Marzari 2018) and where the numerical proportionality factor $C = 8$ is from Isella et al. (2016) and Marzari (2018) and $R_H$ is the planetary Hill's radius. The gas density $\Sigma_{gap}(r)$ inside the gap evolves over time with respect to the local unperturbed gas density $\Sigma(r)$ as $\Sigma_{gap}(r) = \Sigma(r) \cdot \exp[-(t-\tau_p)/\tau_g]$ (Turrini et al. 2021).

The physical radius of the growing giant planets is a critical parameter governing the accretion efficiency of planetesimals. The planetary radius, $R_p$, evolves together with the planetary

mass across the two growth phases following the approach described by Fortier et al. (2013), which is based in turn on the hydrodynamic simulations of Lissauer et al. (2009). During the first phase, the planetary core grows its extended atmosphere, and the physical radius evolves as $R_p = \frac{GM_p}{c_s^2/k_1+(GM_p)/(k_2R_H)}$, where G is the gravitational constant, $M_p$ is the instantaneous mass of the giant planet, $c_s$ is the sound speed in the protoplanetary disk at the orbital distance of the planet, $k_1 = 1$, and $k_2 = 1/4$ (Lissauer et al. 2009).

When the giant planets enter their runaway gas accretion phase (i.e., for $t > \tau_p$), the gravitational infall of the gas causes the planetary radius to shrink as $R_p = R_E - \Delta R(1 - \exp(-(t-\tau_c)/\tau_g))$, where $R_E$ is the planetary radius at the end of the extended atmosphere phase and $\Delta R = R_E - R_T$ is the decrease of the planetary radius during the gravitational collapse of the gas. We adopted as the final values of the planetary radii $R_T$ the ones recently measured from TESS light curves (see Table 1). Particles in the N-body simulations impact one of the giant planets when their relative distance from said planet is less than the planetary radius (see Turrini et al. 2021, for further discussion).

The migration of the giant planets over two growth phases was modeled after the migration tracks from Mordasini et al. (2015) following the parametric approach by Turrini et al. (2021). During the first growth phase, the planets underwent a damped Type I migration regime with a drift rate (Turrini et al. 2021) $\Delta v_1 = \frac{\Delta a_1}{a_p} \frac{\Delta t}{\tau_p} v_p$, where $\Delta t$ is the timestep of the N-body simulation, $\Delta a_1$ is the radial displacement during the first growth phase, and $v_p$ and $a_p$ are the instantaneous planetary orbital velocity and semimajor axis, respectively. During the second growth phase, encompassing the transition to full Type I regime first and Type II regime later, the drift rate became (Hahn & Malhotra 2005; Turrini et al. 2021) $\Delta v_2 = \frac{\Delta a_2}{a_p} \frac{\Delta t}{\tau_g} \exp(-(t-\tau_c)/\tau_g) v_p$, where $\Delta a_2$ is the radial displacement during this second phase (see Sect. 4.1 for the values adopted in the simulations).

We set the spatial density of planetesimals in the N-body simulations to 1000 particles/au, with the inner edge of the planetesimal disk at 1 au and the outer edge at $r_c$ (see Turrini et al. 2021 for the discussion of the choice of the disk inner edge). To compute the mass of the heavy elements accreted by V1298 Tau b and e during their growth and migration, we treated each impacting particle in the N-body simulations as a swarm of real planetesimals. The cumulative mass of each swarm was computed by integrating the disk gas density profile over a ring wide 0.1 au centered on the initial orbit of the impacting particle and multiplying the resulting gas mass by the local solid-to-gas ratio.

The solid-to-gas ratio is a function of the disk metallicity and local disk midplane temperature and is described by a simplified radial profile based on the realistic ones from Turrini et al. (2021) and Pacetti et al. (2022). The disk metallicity was set to 1.4% (Asplund et al. 2009) based on V1298 Tau's solar metallicity (Suárez Mascareño et al. 2022). The solid-to-gas ratio is 0.5 times the disk metallicity for planetesimals formed at temperatures between 1200 K and 140 K (i.e., between the condensation of silicates and that of water). The solid-to-gas ratio grows to 0.75 times the disk metallicity for planetesimals formed at temperatures comprised between 140 K and 30 K (i.e., between the snowlines of water and carbon monoxide), and it reaches 0.9 times the disk metallicity for planetesimals formed at temperatures below 30 K (i.e., beyond the carbon monoxide snowline).
is affected by the gravity of the host star and of the forming planets V1298 Tau b and e, and by the disk gas through aerodynamic drag and gravity (as they possess inertial mass).

The dynamical evolution of the planetesimals is not affected by the interactions among planetesimals themselves (as they do not possess gravitational mass) nor do the planetesimals perturb the two forming planets. The damping effects of gas drag on the planetesimals were simulated following the treatment from Brasser et al. (2007) with updated drag coefficients from Nagasawa et al. (2019) accounting for both the Mach and Reynolds numbers of the planetesimals. The exciting effects of the disk self-gravity were simulated based on the analytical treatment for axisymmetric disks by Ward (1981) following Marzari (2018) and Nagasawa et al. (2019).

The formation of the two giant planets was modeled over two growth and migration phases using the parametric approach from Turrini et al. (2019, 2021). The first phase accounts for their core growth and subsequent capture of an expanded atmosphere (e.g., Bitsch et al. 2015; Johansen et al. 2019; D'Angelo et al. 2021). The planetary mass evolved as $M_p(t) = M_0 + (e/e-1)(M_1 - M_0)(1 - e^{-t/\tau_p})$, where $M_0 = 0.1 M_\oplus$ is the initial mass of the core, $M_1 = 30 M_\oplus$ is the final cumulative mass of core and expanded atmosphere at the end of the first growth phase (see Turrini et al. 2021 for further discussion), and e is the Euler number.

The constant $\tau_p$ is the duration of the first growth phase and was set to about 1 Ma for both planets based on observational and theoretical constraints from circumstellar disks (Manara et al. 2018; Mulders et al. 2021; Bernabò et al. 2022) and the Solar System (Scott 2007; Coradini et al. 2011; Lichtenberg et al. 2023). The individual values of $\tau_p$ are 1 Ma and 1.25 Ma for planets b and e, respectively. The slower growth of planet e was introduced to delay the onset of its second growth and migration phase and to allow planet b to get close to its final orbit before planet e reaches its peak rate of migration. This choice limits the chances of destabilizing close encounters between the two giant planets before they achieve their compact architecture.

The second phase of mass growth accounted for the runaway gas accretion of the two giant planets, where their mass evolves as $M_p(t) = M_1 + (M_2 - M_1)(1 - e^{-(t-\tau_p)/\tau_g})$, with $M_2$ as the final mass of the giant planets and $\tau_g$ as the e-folding time of the runaway gas accretion process. Based on the discussion in Sect. 2.2, $M_2$ was set to 184 $M_\oplus$ for planet b when aiming to reproduce both the HMZe and LMZe scenarios from Table 2. In the case of planet e, $M_2$ was set to 135 $M_\oplus$ and 209 $M_\oplus$ when the simulations focus on the HMZe and LMZe scenarios, respectively.

The value of $\tau_g$ was set to 0.1 Ma based on the results of hydrodynamic simulations (Lissauer et al. 2009; D'Angelo et al. 2021), meaning that the gas giants reach more than 99% of their final mass in about 0.5 Ma from the onset of the runaway gas accretion. During the runaway gas accretion, giant planets form a gap in the disk gas, whose width was modeled as $W_{gap} = C \cdot R_H$ (Isella et al. 2016; Marzari 2018) and where the numerical proportionality factor $C = 8$ is from Isella et al. (2016) and Marzari (2018) and $R_H$ is the planetary Hill's radius. The gas density $\Sigma_{gap}(r)$ inside the gap evolves over time with respect to the local unperturbed gas density $\Sigma(r)$ as $\Sigma_{gap}(r) = \Sigma(r) \cdot \exp[-(t - \tau_p)/\tau_g]$ (Turrini et al. 2021).

The physical radius of the growing giant planets is a critical parameter governing the accretion efficiency of planetesimals. The planetary radius, $R_p$, evolves together with the planetary

mass across the two growth phases following the approach described by Fortier et al. (2013), which is based in turn on the hydrodynamic simulations of Lissauer et al. (2009). During the first phase, the planetary core grows its extended atmosphere, and the physical radius evolves as $R_p = \frac{GM_p}{c_s^2/k_1+(GM_p)/(k_2R_H)}$, where G is the gravitational constant, $M_p$ is the instantaneous mass of the giant planet, $c_s$ is the sound speed in the protoplanetary disk at the orbital distance of the planet, $k_1 = 1$, and $k_2 = 1/4$ (Lissauer et al. 2009).

When the giant planets enter their runaway gas accretion phase (i.e., for $t > \tau_p$), the gravitational infall of the gas causes the planetary radius to shrink as $R_p = R_E - \Delta R(1 - \exp(-(t-\tau_c)/\tau_g))$, where $R_E$ is the planetary radius at the end of the extended atmosphere phase and $\Delta R = R_E - R_T$ is the decrease of the planetary radius during the gravitational collapse of the gas. We adopted as the final values of the planetary radii $R_T$ the ones recently measured from TESS light curves (see Table 1). Particles in the N-body simulations impact one of the giant planets when their relative distance from said planet is less than the planetary radius (see Turrini et al. 2021, for further discussion).

The migration of the giant planets over two growth phases was modeled after the migration tracks from Mordasini et al. (2015) following the parametric approach by Turrini et al. (2021). During the first growth phase, the planets underwent a damped Type I migration regime with a drift rate (Turrini et al. 2021) $\Delta v_1 = \frac{\Delta a_1}{a_p} \frac{\Delta t}{\tau_p} v_p$, where $\Delta t$ is the timestep of the N-body simulation, $\Delta a_1$ is the radial displacement during the first growth phase, and $v_p$ and $a_p$ are the instantaneous planetary orbital velocity and semimajor axis, respectively. During the second growth phase, encompassing the transition to full Type I regime first and Type II regime later, the drift rate became (Hahn & Malhotra 2005; Turrini et al. 2021) $\Delta v_2 = \frac{\Delta a_2}{a_p} \frac{\Delta t}{\tau_g} \exp(-(t-\tau_c)/\tau_g) v_p$, where $\Delta a_2$ is the radial displacement during this second phase (see Sect. 4.1 for the values adopted in the simulations).

We set the spatial density of planetesimals in the N-body simulations to 1000 particles/au, with the inner edge of the planetesimal disk at 1 au and the outer edge at $r_c$ (see Turrini et al. 2021 for the discussion of the choice of the disk inner edge). To compute the mass of the heavy elements accreted by V1298 Tau b and e during their growth and migration, we treated each impacting particle in the N-body simulations as a swarm of real planetesimals. The cumulative mass of each swarm was computed by integrating the disk gas density profile over a ring wide 0.1 au centered on the initial orbit of the impacting particle and multiplying the resulting gas mass by the local solid-to-gas ratio.

The solid-to-gas ratio is a function of the disk metallicity and local disk midplane temperature and is described by a simplified radial profile based on the realistic ones from Turrini et al. (2021) and Pacetti et al. (2022). The disk metallicity was set to 1.4% (Asplund et al. 2009) based on V1298 Tau's solar metallicity (Suárez Mascareño et al. 2022). The solid-to-gas ratio is 0.5 times the disk metallicity for planetesimals formed at temperatures between 1200 K and 140 K (i.e., between the condensation of silicates and that of water). The solid-to-gas ratio grows to 0.75 times the disk metallicity for planetesimals formed at temperatures comprised between 140 K and 30 K (i.e., between the snowlines of water and carbon monoxide), and it reaches 0.9 times the disk metallicity for planetesimals formed at temperatures below 30 K (i.e., beyond the carbon monoxide snowline).
## 3.2. Resonant capture and resonance break simulations

To study the resonant capture process responsible for the present architecture of V1298 Tau, we performed *N*-body integrations on converging orbits of V1298 Tau's planets, simulating their migration due to the interaction with the disk. To perform the *N*-body simulations, we modified the RADAU algorithm (Everhart 1985) to include different damping terms in the semimajor axis, eccentricity, perihelion precession, and velocity of each planet.

The parameters of these damping terms are time dependent and can be tuned to simulate the disk dissipation timescales. To this end, we adopted damping terms that decrease exponentially on a timescale of 0.5 Ma. We assumed that the four planets are fully formed (i.e., the resonant capture simulations take place after the conclusion of the formation simulations) and that the disk is close to its final dispersion. The masses adopted for the planets b and e are the nominal ones of Table 1. For the two inner planets, V1298 Tau c and d, only the upper limits of their mass values $m_c < 0.24 M_J$ and $m_d < 0.31 M_J$ are available (Suárez Mascareño et al. 2022).

We considered two mass configurations for these planets. In the first configuration, the planets are very light, with a density of $\rho = 0.65 \text{ g cm}^{-3}$, giving mass values of $m_c = 0.045 M_J$ and $m_d = 0.077 M_J$. This choice is based on the assumption that the planets are significantly puffed up due to their young age. We also considered a case where the density for both planets is $\rho = 1.3 \text{ g cm}^{-3}$. This case results in mass values $m_c = 0.09 M_J$ and $m_d = 0.15 M_J$, which are still within the observational constraints set by Suárez Mascareño et al. (2022) as well as the upper limits set by Tejada Arevalo et al. (2022).

We used the resonant configurations obtained with these simulations to select the starting conditions to study the breakup of the resonant chain by planet-planet and planetesimal-planet scattering. The simulations of the interactions between massive planetesimal belts and V1298 Tau's planets were performed with Mercury-ArXes, using a version of its hybrid symplectic algorithm ported to GPU computing through OpenACC.

# 4. Results

## 4.1. Constraining the formation tracks of V1298 Tau b and e

The campaign of *N*-body simulations explores the characteristics of the circumstellar disk and the migration tracks that can give rise to the HMZe and LZMe scenarios in Table 2. Due to the uncertainty on the planetary masses and densities of V1298 Tau b and e, the goal of the simulations was to gather indications of how divergent the formation histories of the two giant planets are rather than to pinpoint their exact migration tracks. When comparing the planetesimals accreted in the simulations to the amounts of heavy elements reported in Table 2, we considered the two values to match when the differences are limited to 10–20% (i.e., smaller than the uncertainty on the masses).

We considered three host circumstellar disks (see Table 3). The first disk has a mass of $0.06 M_\odot$, which is comparable to that of a Minimum Mass Solar Nebula-like disk (Hayashi 1981), with the radial extension reported in Sect. 3.1 (simulations 1-4 in Table 3). The second and third disks are twice ($0.12 M_\odot$, simulations 5–12 in Table 3) and three times ($0.18 M_\odot$, simulations 13–18 in Table 3) more massive, respectively. The increasing masses of these disks impact the planetesimal accretion efficiencies of the planets both by providing more solid material and by changing the balance between gas drag and planetary

perturbations. As a consequence, the results do not simply scale linearly with the amount of solid material available.

Table 3 shows the combination of disk parameters, formation regions, and migration scenarios we explored. As illustrated by Table 3 and consistent with the results of Shibata et al. (2020), the amounts of heavy elements accreted by the giant planets can be enhanced or reduced by proportionally increasing or decreasing their orbital displacements (see Sect. 6 for further discussion). The partition of the displacements between the two migration phases was chosen to maximise the planetesimal accretion efficiency of the planets (Shibata et al. 2020; Turrini et al. 2021).

The total mass of solids contained in the first disk ($0.06 M_\odot$) is not enough to support the HMZe scenario, so in simulations 1-4, we focused exclusively on the LMZe scenario. The metallicity values of the two giant planets in the LMZe scenario were obtained when V1298 Tau b began its formation around 10 au, and V1298 Tau e was about twice as distant as showcased by simulations 3 and 4. This picture is qualitatively preserved in the two more massive disks. The LMZe scenario was reproduced by having V1298 Tau b start its formation slightly closer to the star (about 7 au in both disks) and having V1298 Tau e initially located 1.5–2 times farther away (see simulations 12 and 18 in Table 3).

The HMZe scenario could only be reproduced by the most massive disk considered in our simulations ($0.18 M_\odot$). Even in this case, the high metallicity value of V1298 Tau e required the planet to start its formation at the outer edge of the disk in order to encounter and accrete enough planetesimals (see simulations 16). In the intermediate disk ($0.12 M_\odot$), the same extreme migration track results in about half the metallicity required for V1298 Tau e in the HMZe scenario (see simulation 9 in Table 3). An additional process, such as the accretion of disk gas enriched in heavy elements (e.g., Booth & Ilee 2019; Cridland et al. 2019; Schneider & Bitsch 2021), needs to be invoked together with planetesimal accretion to explain the missing metallicity.

The larger enrichment required in the case of the lower density values of the inflated gaseous envelopes expected for the two hot and young giant planets (87 and 108 $M_\odot$ for V1298 Tau b and e in the LMZe scenario, see Sect. 2.2) point to both planets having undergone extensive migration in a disk whose mass was at least 10% that of the host star (i.e., similar or more massive than our intermediate disk). The results reported in Table 3 suggest that in this scenario, V1298 Tau b should have started its formation about midway through its host disk and V1298 Tau e about twice as far (i.e., closer to the outer edge of the disk). Lower planetary mass values than the nominal ones reported in Table 1 or the presence of massive cores (which we ignored in our computations) would lower the required planetesimal accretion and shift the initial formation regions of the planets inward.

As the previous results and discussion show, the planetary density and metallicity values allow the discussion of the formation histories of the two giant planets only from a comparative point of view. Specifically, the picture emerging from our campaign of simulations consistently indicates that V1298 Tau e should have formed further out than V1298 Tau b, likely about twice as far from their host star. The uncertainty regarding their masses and our ignorance of the characteristics of their now long-gone native disk, however, prevent us from precisely pinpointing the formation regions of the two giant planets. We discuss this issue further in Sect. 6.1 and highlight how future atmospheric characterization of the two planets will allow us to overcome it in Sect. 6.2.
## 3.2. Resonant capture and resonance break simulations

To study the resonant capture process responsible for the present architecture of V1298 Tau, we performed *N*-body integrations on converging orbits of V1298 Tau's planets, simulating their migration due to the interaction with the disk. To perform the *N*-body simulations, we modified the RADAU algorithm (Everhart 1985) to include different damping terms in the semimajor axis, eccentricity, perihelion precession, and velocity of each planet.

The parameters of these damping terms are time dependent and can be tuned to simulate the disk dissipation timescales. To this end, we adopted damping terms that decrease exponentially on a timescale of 0.5 Ma. We assumed that the four planets are fully formed (i.e., the resonant capture simulations take place after the conclusion of the formation simulations) and that the disk is close to its final dispersion. The masses adopted for the planets b and e are the nominal ones of Table 1. For the two inner planets, V1298 Tau c and d, only the upper limits of their mass values $m_c < 0.24 M_J$ and $m_d < 0.31 M_J$ are available (Suárez Mascareño et al. 2022).

We considered two mass configurations for these planets. In the first configuration, the planets are very light, with a density of $\rho = 0.65 \text{ g cm}^{-3}$, giving mass values of $m_c = 0.045 M_J$ and $m_d = 0.077 M_J$. This choice is based on the assumption that the planets are significantly puffed up due to their young age. We also considered a case where the density for both planets is $\rho = 1.3 \text{ g cm}^{-3}$. This case results in mass values $m_c = 0.09 M_J$ and $m_d = 0.15 M_J$, which are still within the observational constraints set by Suárez Mascareño et al. (2022) as well as the upper limits set by Tejada Arevalo et al. (2022).

We used the resonant configurations obtained with these simulations to select the starting conditions to study the breakup of the resonant chain by planet-planet and planetesimal-planet scattering. The simulations of the interactions between massive planetesimal belts and V1298 Tau's planets were performed with Mercury-Arχes, using a version of its hybrid symplectic algorithm ported to GPU computing through OpenACC.

# 4. Results

## 4.1. Constraining the formation tracks of V1298 Tau b and e

The campaign of *N*-body simulations explores the characteristics of the circumstellar disk and the migration tracks that can give rise to the HMZe and LZMe scenarios in Table 2. Due to the uncertainty on the planetary masses and densities of V1298 Tau b and e, the goal of the simulations was to gather indications of how divergent the formation histories of the two giant planets are rather than to pinpoint their exact migration tracks. When comparing the planetesimals accreted in the simulations to the amounts of heavy elements reported in Table 2, we considered the two values to match when the differences are limited to 10–20% (i.e., smaller than the uncertainty on the masses).

We considered three host circumstellar disks (see Table 3). The first disk has a mass of $0.06 M_\odot$, which is comparable to that of a Minimum Mass Solar Nebula-like disk (Hayashi 1981), with the radial extension reported in Sect. 3.1 (simulations 1-4 in Table 3). The second and third disks are twice ($0.12 M_\odot$, simulations 5–12 in Table 3) and three times ($0.18 M_\odot$, simulations 13–18 in Table 3) more massive, respectively. The increasing masses of these disks impact the planetesimal accretion efficiencies of the planets both by providing more solid material and by changing the balance between gas drag and planetary

perturbations. As a consequence, the results do not simply scale linearly with the amount of solid material available.

Table 3 shows the combination of disk parameters, formation regions, and migration scenarios we explored. As illustrated by Table 3 and consistent with the results of Shibata et al. (2020), the amounts of heavy elements accreted by the giant planets can be enhanced or reduced by proportionally increasing or decreasing their orbital displacements (see Sect. 6 for further discussion). The partition of the displacements between the two migration phases was chosen to maximise the planetesimal accretion efficiency of the planets (Shibata et al. 2020; Turrini et al. 2021).

The total mass of solids contained in the first disk ($0.06 M_\odot$) is not enough to support the HMZe scenario, so in simulations 1-4, we focused exclusively on the LMZe scenario. The metallicity values of the two giant planets in the LMZe scenario were obtained when V1298 Tau b began its formation around 10 au, and V1298 Tau e was about twice as distant as showcased by simulations 3 and 4. This picture is qualitatively preserved in the two more massive disks. The LMZe scenario was reproduced by having V1298 Tau b start its formation slightly closer to the star (about 7 au in both disks) and having V1298 Tau e initially located 1.5–2 times farther away (see simulations 12 and 18 in Table 3).

The HMZe scenario could only be reproduced by the most massive disk considered in our simulations ($0.18 M_\odot$). Even in this case, the high metallicity value of V1298 Tau e required the planet to start its formation at the outer edge of the disk in order to encounter and accrete enough planetesimals (see simulations 16). In the intermediate disk ($0.12 M_\odot$), the same extreme migration track results in about half the metallicity required for V1298 Tau e in the HMZe scenario (see simulation 9 in Table 3). An additional process, such as the accretion of disk gas enriched in heavy elements (e.g., Booth & Ilee 2019; Cridland et al. 2019; Schneider & Bitsch 2021), needs to be invoked together with planetesimal accretion to explain the missing metallicity.

The larger enrichment required in the case of the lower density values of the inflated gaseous envelopes expected for the two hot and young giant planets (87 and 108 $M_\odot$ for V1298 Tau b and e in the LMZe scenario, see Sect. 2.2) point to both planets having undergone extensive migration in a disk whose mass was at least 10% that of the host star (i.e., similar or more massive than our intermediate disk). The results reported in Table 3 suggest that in this scenario, V1298 Tau b should have started its formation about midway through its host disk and V1298 Tau e about twice as far (i.e., closer to the outer edge of the disk). Lower planetary mass values than the nominal ones reported in Table 1 or the presence of massive cores (which we ignored in our computations) would lower the required planetesimal accretion and shift the initial formation regions of the planets inward.

As the previous results and discussion show, the planetary density and metallicity values allow the discussion of the formation histories of the two giant planets only from a comparative point of view. Specifically, the picture emerging from our campaign of simulations consistently indicates that V1298 Tau e should have formed further out than V1298 Tau b, likely about twice as far from their host star. The uncertainty regarding their masses and our ignorance of the characteristics of their now long-gone native disk, however, prevent us from precisely pinpointing the formation regions of the two giant planets. We discuss this issue further in Sect. 6.1 and highlight how future atmospheric characterization of the two planets will allow us to overcome it in Sect. 6.2.
<table><caption>**Table 3.** Parameters of the circumstellar disks and the migration scenarios considered in the 18 formation simulations for V1298 Tau's planets b and e, and the resulting accretion of heavy elements of each scenario.</caption><thead><tr><th rowspan="2">Id.</th><th colspan="3">Disk</th><th colspan="4">V1298 Tau b</th><th colspan="4">V1298 Tau e</th></tr><tr><th>M<sub>disk</sub><br>(M<sub>☉</sub>)</th><th>Σ<sub>0</sub><br>(g cm<sup>-3</sup>)</th><th>M<sub>dust</sub><br>(M<sub>⊕</sub>)</th><th>Initial<br>seed (au)</th><th>Δa<sub>1</sub><br>(au)</th><th>Δa<sub>2</sub><br>(au)</th><th>Accreted<br>mass (M<sub>⊕</sub>)</th><th>Initial<br>seed (au)</th><th>Δa<sub>1</sub><br>(au)</th><th>Δa<sub>2</sub><br>(au)</th><th>Accreted<br>mass (M<sub>⊕</sub>)</th></tr></thead><tbody><tr><td>1</td><td>0.06</td><td>42</td><td>205</td><td>11</td><td>-6.4</td><td>-4.3</td><td>22</td><td>50</td><td>-29.6</td><td>-19.8</td><td>54</td></tr><tr><td>2</td><td>0.06</td><td>42</td><td>205</td><td>10</td><td>-5.8</td><td>-3.9</td><td>19</td><td>30</td><td>-17.6</td><td>-11.8</td><td>46</td></tr><tr><td>3</td><td>0.06</td><td>42</td><td>205</td><td>10</td><td>-5.8</td><td>-3.9</td><td>18</td><td>20</td><td>-11.6</td><td>-7.8</td><td>26</td></tr><tr><td>4</td><td>0.06</td><td>42</td><td>205</td><td>11</td><td>-6.4</td><td>-4.3</td><td>21</td><td>22</td><td>-12.8</td><td>-8.6</td><td>28</td></tr><tr><td>5</td><td>0.12</td><td>80</td><td>390</td><td>6</td><td>-3.4</td><td>-2.3</td><td>15</td><td>50</td><td>-29.6</td><td>-19.8</td><td>139</td></tr><tr><td>6</td><td>0.12</td><td>80</td><td>390</td><td>6</td><td>-3.4</td><td>-2.3</td><td>13</td><td>15</td><td>-8.6</td><td>-5.8</td><td>52</td></tr><tr><td>7</td><td>0.12</td><td>80</td><td>390</td><td>9</td><td>-5.2</td><td>-3.5</td><td>30</td><td>50</td><td>-29.6</td><td>-19.8</td><td>132</td></tr><tr><td>8</td><td>0.12</td><td>80</td><td>390</td><td>9</td><td>-5.2</td><td>-3.5</td><td>36</td><td>12.5</td><td>-7.1</td><td>-4.8</td><td>11</td></tr><tr><td>9</td><td>0.12</td><td>80</td><td>390</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>20</td><td>50</td><td>-29.6</td><td>-19.8</td><td>129</td></tr><tr><td>10</td><td>0.12</td><td>80</td><td>390</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>19</td><td>15</td><td>-8.6</td><td>-5.8</td><td>41</td></tr><tr><td>11</td><td>0.12</td><td>80</td><td>390</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>20</td><td>12.5</td><td>-7.1</td><td>-4.8</td><td>25</td></tr><tr><td>12</td><td>0.12</td><td>80</td><td>390</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>19</td><td>13</td><td>-7.4</td><td>-5.0</td><td>28</td></tr><tr><td>13</td><td>0.18</td><td>120</td><td>585</td><td>5</td><td>-2.8</td><td>-1.9</td><td>7</td><td>50</td><td>-29.6</td><td>-19.8</td><td>250</td></tr><tr><td>14</td><td>0.18</td><td>120</td><td>585</td><td>10</td><td>-5.8</td><td>-3.9</td><td>48</td><td>50</td><td>-29.6</td><td>-19.8</td><td>233</td></tr><tr><td>15</td><td>0.18</td><td>120</td><td>585</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>26</td><td>50</td><td>-29.6</td><td>-19.8</td><td>234</td></tr><tr><td>16</td><td>0.18</td><td>120</td><td>585</td><td>7</td><td>-4</td><td>-2.7</td><td>23</td><td>50</td><td>-29.6</td><td>-19.8</td><td>245</td></tr><tr><td>17</td><td>0.18</td><td>120</td><td>585</td><td>7</td><td>-4</td><td>-2.7</td><td>23</td><td>12.5</td><td>-7.1</td><td>-4.8</td><td>42</td></tr><tr><td>18</td><td>0.18</td><td>120</td><td>585</td><td>7</td><td>-4</td><td>-2.7</td><td>23</td><td>11</td><td>-6.2</td><td>-4.2</td><td>28</td></tr></tbody></table>

**Notes.** The rows highlighted in blue identify the simulations consistent with the LMZe scenario in Table 2. The row highlighted in red is the simulation consistent with the HMZe scenario, while the row highlighted in yellow identifies the simulation consistent with the HMZe scenario when assuming a second process assists in increasing the metallicity of planet e.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Dynamical excitation of the planetesimal disk after the formation and migration of V1298 Tau b and e in scenario 9 from Table 3. Left panel: inclination. The empty circles in the panels mark the initial position of the giant planets. The simulation stopped with the giant planets about three times more distant from V1298 Tau than their current orbits to prevent dynamical instabilities that would result in a system incompatible with V1298 Tau's compact architecture.</figcaption></figure>

## 4.2. Impact of the formation tracks of V1298 Tau b and e on the subsequent evolution of the system

As shown in Fig. 2, the migration histories of V1298 Tau b and e shape the characteristics of the surviving planetesimal disk in two ways. First, the two migrating giant planets leave in their wake surviving planetesimal disks that are characterized by two different dynamical regions (see Fig. 2). The innermost region,

roughly extending from the current orbits of the giant planets to the initial formation region of V1298 Tau b (see Fig. 2), is depopulated of planetesimals. Furthermore, due to the higher gas density in the inner disk regions, the dynamical excitation of the surviving planetesimals is efficiently damped by gas drag. The outermost region, roughly extending between the initial formation regions of V1298 Tau b and e (see Fig. 2), is populated by dynamically excited planetesimals.
<table><caption>**Table 3.** Parameters of the circumstellar disks and the migration scenarios considered in the 18 formation simulations for V1298 Tau's planets b and e, and the resulting accretion of heavy elements of each scenario.</caption><thead><tr><th rowspan="2">Id.</th><th colspan="3">Disk</th><th colspan="4">V1298 Tau b</th><th colspan="4">V1298 Tau e</th></tr><tr><th><i>M</i><sub>disk</sub><br>(<i>M</i><sub>☉</sub>)</th><th>Σ<sub>0</sub><br>(g cm<sup>-3</sup>)</th><th><i>M</i><sub>dust</sub><br>(<i>M</i><sub>⊕</sub>)</th><th>Initial<br>seed (au)</th><th>Δ<i>a</i><sub>1</sub><br>(au)</th><th>Δ<i>a</i><sub>2</sub><br>(au)</th><th>Accreted<br>mass (<i>M</i><sub>⊕</sub>)</th><th>Initial<br>seed (au)</th><th>Δ<i>a</i><sub>1</sub><br>(au)</th><th>Δ<i>a</i><sub>2</sub><br>(au)</th><th>Accreted<br>mass (<i>M</i><sub>⊕</sub>)</th></tr></thead><tbody><tr><td>1</td><td>0.06</td><td>42</td><td>205</td><td>11</td><td>-6.4</td><td>-4.3</td><td>22</td><td>50</td><td>-29.6</td><td>-19.8</td><td>54</td></tr><tr><td>2</td><td>0.06</td><td>42</td><td>205</td><td>10</td><td>-5.8</td><td>-3.9</td><td>19</td><td>30</td><td>-17.6</td><td>-11.8</td><td>46</td></tr><tr><td>3</td><td>0.06</td><td>42</td><td>205</td><td>10</td><td>-5.8</td><td>-3.9</td><td>18</td><td>20</td><td>-11.6</td><td>-7.8</td><td>26</td></tr><tr><td>4</td><td>0.06</td><td>42</td><td>205</td><td>11</td><td>-6.4</td><td>-4.3</td><td>21</td><td>22</td><td>-12.8</td><td>-8.6</td><td>28</td></tr><tr><td>5</td><td>0.12</td><td>80</td><td>390</td><td>6</td><td>-3.4</td><td>-2.3</td><td>15</td><td>50</td><td>-29.6</td><td>-19.8</td><td>139</td></tr><tr><td>6</td><td>0.12</td><td>80</td><td>390</td><td>6</td><td>-3.4</td><td>-2.3</td><td>13</td><td>15</td><td>-8.6</td><td>-5.8</td><td>52</td></tr><tr><td>7</td><td>0.12</td><td>80</td><td>390</td><td>9</td><td>-5.2</td><td>-3.5</td><td>30</td><td>50</td><td>-29.6</td><td>-19.8</td><td>132</td></tr><tr><td>8</td><td>0.12</td><td>80</td><td>390</td><td>9</td><td>-5.2</td><td>-3.5</td><td>36</td><td>12.5</td><td>-7.1</td><td>-4.8</td><td>11</td></tr><tr><td>9</td><td>0.12</td><td>80</td><td>390</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>20</td><td>50</td><td>-29.6</td><td>-19.8</td><td>129</td></tr><tr><td>10</td><td>0.12</td><td>80</td><td>390</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>19</td><td>15</td><td>-8.6</td><td>-5.8</td><td>41</td></tr><tr><td>11</td><td>0.12</td><td>80</td><td>390</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>20</td><td>12.5</td><td>-7.1</td><td>-4.8</td><td>25</td></tr><tr><td>12</td><td>0.12</td><td>80</td><td>390</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>19</td><td>13</td><td>-7.4</td><td>-5.0</td><td>28</td></tr><tr><td>13</td><td>0.18</td><td>120</td><td>585</td><td>5</td><td>-2.8</td><td>-1.9</td><td>7</td><td>50</td><td>-29.6</td><td>-19.8</td><td>250</td></tr><tr><td>14</td><td>0.18</td><td>120</td><td>585</td><td>10</td><td>-5.8</td><td>-3.9</td><td>48</td><td>50</td><td>-29.6</td><td>-19.8</td><td>233</td></tr><tr><td>15</td><td>0.18</td><td>120</td><td>585</td><td>7.5</td><td>-4.3</td><td>-2.9</td><td>26</td><td>50</td><td>-29.6</td><td>-19.8</td><td>234</td></tr><tr><td>16</td><td>0.18</td><td>120</td><td>585</td><td>7</td><td>-4</td><td>-2.7</td><td>23</td><td>50</td><td>-29.6</td><td>-19.8</td><td>245</td></tr><tr><td>17</td><td>0.18</td><td>120</td><td>585</td><td>7</td><td>-4</td><td>-2.7</td><td>23</td><td>12.5</td><td>-7.1</td><td>-4.8</td><td>42</td></tr><tr><td>18</td><td>0.18</td><td>120</td><td>585</td><td>7</td><td>-4</td><td>-2.7</td><td>23</td><td>11</td><td>-6.2</td><td>-4.2</td><td>28</td></tr></tbody></table>

**Notes.** The rows highlighted in blue identify the simulations consistent with the LMZe scenario in Table 2. The row highlighted in red is the simulation consistent with the HMZe scenario, while the row highlighted in yellow identifies the simulation consistent with the HMZe scenario when assuming a second process assists in increasing the metallicity of planet e.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Dynamical excitation of the planetesimal disk after the formation and migration of V1298 Tau b and e in scenario 9 from Table 3. Left panel: inclination. The empty circles in the panels mark the initial position of the giant planets. The simulation stopped with the giant planets about three times more distant from V1298 Tau than their current orbits to prevent dynamical instabilities that would result in a system incompatible with V1298 Tau's compact architecture.</figcaption></figure>

## 4.2. Impact of the formation tracks of V1298 Tau b and e on the subsequent evolution of the system

As shown in Fig. 2, the migration histories of V1298 Tau b and e shape the characteristics of the surviving planetesimal disk in two ways. First, the two migrating giant planets leave in their wake surviving planetesimal disks that are characterized by two different dynamical regions (see Fig. 2). The innermost region,

roughly extending from the current orbits of the giant planets to the initial formation region of V1298 Tau b (see Fig. 2), is depopulated of planetesimals. Furthermore, due to the higher gas density in the inner disk regions, the dynamical excitation of the surviving planetesimals is efficiently damped by gas drag. The outermost region, roughly extending between the initial formation regions of V1298 Tau b and e (see Fig. 2), is populated by dynamically excited planetesimals.
Second, the larger migration of V1298 Tau e favors its acquisition of a marked population of Trojan satellites (see the planetesimals “piled up” on V1298 Tau e in Fig. 2), planetesimals co-orbiting with the giant planet in stable orbits around its L4 and L5 Lagrangian points with the star, which is similar to what is argued in the case of the Solar System for Jupiter (Pirani et al. 2019). The formation history of V1298 Tau b appears to hinder its efficiency in capturing Trojan satellites, although this could be an artefact of the specific migration tracks we adopted. Dedicated studies of the convergent migration and resonant trapping of V1298 Tau b and e are required to assess the dynamical lifetime of possible Trojan bodies.

The population of dynamically excited planetesimals created by the forming giant planets can affect the later evolution of the planetary system in two ways. On the one hand, the resulting higher impact frequency among planetesimals promotes the growth of existing oligarchs and the formation of additional massive planetary bodies, whose mutual interactions can inject them on eccentric orbits outside the present orbits of the four known planets. Later encounters between such excited outer planets and the four resonant inner planets could result in the breakup of the resonant chain and the onset of dynamical instability. On the other hand, the high-velocity impacts between the excited planetesimals cause large-scale production of collisional dust, replenishing the dust population of the circumstellar disk, as shown by the results of Turrini et al. (2019) and Bernabò et al. (2022).

We used the DEBRIS code (Turrini et al. 2019; Bernabò et al. 2022) to estimate the possible amount of collisional dust produced in simulations 4 and 9 (see Table 3 and Fig. 2 for the excited planetesimal disk at the end of simulation 9) by the migration of the V1298 Tau b and e. We assumed the planetesimal population to be characterized by the size-frequency distribution and mechanical strength typical of primordial planetesimals formed by streaming instability (Krivov et al. 2018; Turrini et al. 2019). We refer the interested reader to Turrini et al. (2019) and Bernabò et al. (2022) for details on the collisional and dust production algorithms implemented in the DEBRIS code.

The collisional cascade among the surviving excited planetesimals is capable of converting back into dust between 50 $M_⊕$ (simulation 4) and 160 $M_⊕$ (simulation 9). Such an amount of collisional dust, while drifting inward toward the star, is large enough to impact the formation history of V1298 Tau's system. First, it promotes pebble accretion on surviving planetary embryos not accreted by the migrating planets, resulting in the formation of additional massive planets. We discuss the role of such planets in shaping the current architecture of V1298 Tau in Sect. 4.4. Second, it promotes new phases of streaming instability and can form massive planetesimal belts outside the current orbits of V1298 Tau e. The gravitational interactions between such belts and V1298 Tau's four resonant planets can also potentially break the resonance chain. We investigate this scenario in Sect. 4.5.

### 4.3. Primordial capture in four-body resonance as V1298 Tau's original architecture

As discussed in Sect. 2.1, the current architecture of V1298 Tau's planets is not characterized by the resonant chain expected for such compact systems and shows the signs of having been sculpted by dynamical instabilities. To understand the origin of these characteristics, we first explored the possible primordial architectures of the system. There are two possible Laplace resonance chains that are dynamically close to the nominal solution

and stable. They are the 3:2, 2:1, 3:2 chain and the 3:2, 2:1, 2:1 chain.

A possible scenario for the formation of either of the two resonant chains is that the planets migrated into resonance while the circumstellar disk was still present which is consistent with the migration scenarios simulated in Sect. 4.1 to explain the density values of V1298 Tau b and e. The capture in both these Laplace resonances was simulated with the parallel version of the N-body integration code RADAU15 (Everhart 1985) described in Sect. 3.2 using exponential damping terms in the semimajor axis and eccentricity with e-folding times of 0.5 Ma.

The masses adopted for the planets b and e are the nominal ones of Table 1, while for planets c and d, we considered two configurations (see Sect. 3.2): first, one configuration where the planets have a density of $\rho = 0.65 \text{ g cm}^{-3}$ and masses $m_c = 0.045 M_J$ and $m_d = 0.077 M_J$, and a second one where they have a density of $\rho = 1.3 \text{ g cm}^{-3}$ and masses $m_c = 0.09 M_J$ and $m_d = 0.15 M_J$. We found no significant differences between the two configurations.

The way in which the real capture in resonance may have occurred depends on a large number of parameters, including the initial orbits of the planets, the disk gas density profile and the timescale of its dissipation, the local effects of the planets on the gas, and the viscosity of the disk and its temperature profile. As it is impossible to perform a complete exploration of all these parameters, we first found a case in which the trapping occurs to confirm that it is indeed a realistic scenario.

Then, starting from this single case, we explored the phase space near this resonant solution in order to find all the possible resonant configurations. This approach allowed us to map the range of the planetary orbital elements in the specific Laplace resonant configuration. A random search without the knowledge of the orbital elements of at least one resonant configuration would be impossible. Depending on the above-mentioned parameters influencing the planet migration rate, different resonant configurations, among those we have found, were achieved during the evolution of the system.

The detailed exploration of the resonant phase space was performed by randomly sampling all planetary orbital elements in the proximity of the resonant solution. During the numerical integration of the planetary orbits, the critical angles of the Laplace resonance were automatically checked, and the non-librating, non-apsidal corotation cases were rejected. This last condition was dictated by the results of Beaugé et al. (2006), which suggest that as long as the migration is sufficiently slow to be approximated as an adiabatic process, all captured planets must be in apsidal corotations. The range we found in eccentricity shows the possible final orbital configurations of the four planets at the end of migration and resonance capture for different potential initial disk and planet parameters.

In Fig. 3, we illustrate the trapping in the 3:2, 2:1, 3:2 resonant chain where at the beginning of the simulation, the planets are close to the resonance and become captured during the subsequent inward migration. The eccentricity of each planet at the resonance trapping is pumped up to a value that is kept constant until the end of the migration. These high values are compatible with observations. The critical arguments of the three individual resonances librate around different values. In Fig. 4, we show the second scenario of resonance trapping where the planets are captured in a 3:2, 2:1, 2:1 resonant chain. Even in this case, the eccentricity is pumped up, and the resonant angles are all librating. We note that these simulations can be scaled in semimajor axis to get as close as possible to the observed semimajor axes of the planets.
Second, the larger migration of V1298 Tau e favors its acquisition of a marked population of Trojan satellites (see the planetesimals “piled up” on V1298 Tau e in Fig. 2), planetesimals co-orbiting with the giant planet in stable orbits around its L4 and L5 Lagrangian points with the star, which is similar to what is argued in the case of the Solar System for Jupiter (Pirani et al. 2019). The formation history of V1298 Tau b appears to hinder its efficiency in capturing Trojan satellites, although this could be an artefact of the specific migration tracks we adopted. Dedicated studies of the convergent migration and resonant trapping of V1298 Tau b and e are required to assess the dynamical lifetime of possible Trojan bodies.

The population of dynamically excited planetesimals created by the forming giant planets can affect the later evolution of the planetary system in two ways. On the one hand, the resulting higher impact frequency among planetesimals promotes the growth of existing oligarchs and the formation of additional massive planetary bodies, whose mutual interactions can inject them on eccentric orbits outside the present orbits of the four known planets. Later encounters between such excited outer planets and the four resonant inner planets could result in the breakup of the resonant chain and the onset of dynamical instability. On the other hand, the high-velocity impacts between the excited planetesimals cause large-scale production of collisional dust, replenishing the dust population of the circumstellar disk, as shown by the results of Turrini et al. (2019) and Bernabò et al. (2022).

We used the DEBRIS code (Turrini et al. 2019; Bernabò et al. 2022) to estimate the possible amount of collisional dust produced in simulations 4 and 9 (see Table 3 and Fig. 2 for the excited planetesimal disk at the end of simulation 9) by the migration of the V1298 Tau b and e. We assumed the planetesimal population to be characterized by the size-frequency distribution and mechanical strength typical of primordial planetesimals formed by streaming instability (Krivov et al. 2018; Turrini et al. 2019). We refer the interested reader to Turrini et al. (2019) and Bernabò et al. (2022) for details on the collisional and dust production algorithms implemented in the DEBRIS code.

The collisional cascade among the surviving excited planetesimals is capable of converting back into dust between 50 $M_⊕$ (simulation 4) and 160 $M_⊕$ (simulation 9). Such an amount of collisional dust, while drifting inward toward the star, is large enough to impact the formation history of V1298 Tau's system. First, it promotes pebble accretion on surviving planetary embryos not accreted by the migrating planets, resulting in the formation of additional massive planets. We discuss the role of such planets in shaping the current architecture of V1298 Tau in Sect. 4.4. Second, it promotes new phases of streaming instability and can form massive planetesimal belts outside the current orbits of V1298 Tau e. The gravitational interactions between such belts and V1298 Tau's four resonant planets can also potentially break the resonance chain. We investigate this scenario in Sect. 4.5.

### 4.3. Primordial capture in four-body resonance as V1298 Tau's original architecture

As discussed in Sect. 2.1, the current architecture of V1298 Tau's planets is not characterized by the resonant chain expected for such compact systems and shows the signs of having been sculpted by dynamical instabilities. To understand the origin of these characteristics, we first explored the possible primordial architectures of the system. There are two possible Laplace resonance chains that are dynamically close to the nominal solution

and stable. They are the 3:2, 2:1, 3:2 chain and the 3:2, 2:1, 2:1 chain.

A possible scenario for the formation of either of the two resonant chains is that the planets migrated into resonance while the circumstellar disk was still present which is consistent with the migration scenarios simulated in Sect. 4.1 to explain the density values of V1298 Tau b and e. The capture in both these Laplace resonances was simulated with the parallel version of the N-body integration code RADAU15 (Everhart 1985) described in Sect. 3.2 using exponential damping terms in the semimajor axis and eccentricity with e-folding times of 0.5 Ma.

The masses adopted for the planets b and e are the nominal ones of Table 1, while for planets c and d, we considered two configurations (see Sect. 3.2): first, one configuration where the planets have a density of $\rho = 0.65 \text{ g cm}^{-3}$ and masses $m_c = 0.045 M_J$ and $m_d = 0.077 M_J$, and a second one where they have a density of $\rho = 1.3 \text{ g cm}^{-3}$ and masses $m_c = 0.09 M_J$ and $m_d = 0.15 M_J$. We found no significant differences between the two configurations.

The way in which the real capture in resonance may have occurred depends on a large number of parameters, including the initial orbits of the planets, the disk gas density profile and the timescale of its dissipation, the local effects of the planets on the gas, and the viscosity of the disk and its temperature profile. As it is impossible to perform a complete exploration of all these parameters, we first found a case in which the trapping occurs to confirm that it is indeed a realistic scenario.

Then, starting from this single case, we explored the phase space near this resonant solution in order to find all the possible resonant configurations. This approach allowed us to map the range of the planetary orbital elements in the specific Laplace resonant configuration. A random search without the knowledge of the orbital elements of at least one resonant configuration would be impossible. Depending on the above-mentioned parameters influencing the planet migration rate, different resonant configurations, among those we have found, were achieved during the evolution of the system.

The detailed exploration of the resonant phase space was performed by randomly sampling all planetary orbital elements in the proximity of the resonant solution. During the numerical integration of the planetary orbits, the critical angles of the Laplace resonance were automatically checked, and the non-librating, non-apsidal corotation cases were rejected. This last condition was dictated by the results of Beaugé et al. (2006), which suggest that as long as the migration is sufficiently slow to be approximated as an adiabatic process, all captured planets must be in apsidal corotations. The range we found in eccentricity shows the possible final orbital configurations of the four planets at the end of migration and resonance capture for different potential initial disk and planet parameters.

In Fig. 3, we illustrate the trapping in the 3:2, 2:1, 3:2 resonant chain where at the beginning of the simulation, the planets are close to the resonance and become captured during the subsequent inward migration. The eccentricity of each planet at the resonance trapping is pumped up to a value that is kept constant until the end of the migration. These high values are compatible with observations. The critical arguments of the three individual resonances librate around different values. In Fig. 4, we show the second scenario of resonance trapping where the planets are captured in a 3:2, 2:1, 2:1 resonant chain. Even in this case, the eccentricity is pumped up, and the resonant angles are all librating. We note that these simulations can be scaled in semimajor axis to get as close as possible to the observed semimajor axes of the planets.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Migration of the four planets in resonant configuration with the planets in the 3:2, 2:1, and 3:2 resonances. The *left panel* shows the semimajor axis evolution, the *middle panel* shows the eccentricity, and the *right panel* shows the resonant critical angles during the migration and subsequent resonance capture. The gas was slowly dissipating across the simulation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Same as Fig. 3 but for the resonance sequence 3:2, 2:1, and 2:1.</figcaption></figure>

In subsequent detailed explorations of the phase space, we numerically integrated millions of initial planetary configurations to find the extension of the phase space of the Laplace resonant configurations. In Fig. 5, we show the ranges of the initial values of semimajor axis (very tiny) and eccentricity with all four planets locked in a Laplace resonance in both cases (3:2, 2:1, and 3:2 on the left and 3:2, 2:1, and 2:1 on the right). Within these ranges, there are resonant solutions that are compatible in terms of eccentricity with the observed nominal system, suggesting that the primordial system was indeed in a Laplace resonance, and it later escaped from the resonant configuration. When comparing the two different configurations, 3:2, 2:1, and 3:2 versus 3:2, 2:1, and 2:1, we find it interesting to note that the last configuration leads to higher eccentricities for the third planet.

It is worth pointing out that the resonant solutions plotted in Fig. 5 were obtained by requiring that both resonant arguments librate, a condition naturally leading to apsidal libration. This condition is in agreement with the work of Beaugé et al. (2006), who suggest that as long as the orbital migration is sufficiently slow to be approximated by an adiabatic process (which is a condition compatible with the slow migration of giant planets toward the end of their formation; see e.g., Pirani et al. 2019; Tanaka et al. 2020), all captured planets should be in apsidal co-rotations.

Even if the resonant lock provides a solution for the origin of V1298 Tau's compact architecture, the system is not in resonant chain at present, as confirmed by the extensive parameter exploration performed by Tejada Arevalo et al. (2022). The ratios between the planetary semimajor axes, according to the current nominal solutions from Table 1, are 1.312, 1.559, and 1.535,

while to be in a four-body resonance, they should be 1.310, 1.587, and 1.310 or 1.310, 1.587, and 1.587.

This mismatch can be explained if the system evolved into either of the two resonant chains in the initial stages of its formation history but the resonance lock was broken during its later evolution by some additional dynamical mechanism. According to the capture simulations with a dissipating disk (Fig. 3 and Fig. 4), it is difficult for the gas dispersal alone to break the resonance lock. Dynamically destabilizing mechanisms must be invoked, specifically the planet-planet scattering and planetesimals scattering introduced in Sect. 4.1 and explored in detail in the next sections.

### 4.4. V1298 Tau's present architecture: Resonance break by planet-planet scattering

The first possible scenario to break the resonant lock requires that additional dynamically excited planets populated, and maybe still populate, the system on orbits external to that of planet e. These excited planets would be difficult to detect through transits due to their longer orbital periods. As discussed in Sect. 4.2, the convergent migration shaping the formation path of V1298 Tau could act to promote the formation of such planets. These additional planets would have formed at later times in the wake of the passage of V1298 Tau b and e (see Sect. 4.2) and would not have migrated as close to the star as planets b and e either because of the dissipation of the gaseous disk or because they began to interact gravitationally with each other (Marzari et al. 2010). The dissipation of the disk itself can cause the onset of the dynamical instability among the resulting population of planets.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Migration of the four planets in resonant configuration with the planets in the 3:2, 2:1, and 3:2 resonances. The *left panel* shows the semimajor axis evolution, the *middle panel* shows the eccentricity, and the *right panel* shows the resonant critical angles during the migration and subsequent resonance capture. The gas was slowly dissipating across the simulation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Same as Fig. 3 but for the resonance sequence 3:2, 2:1, and 2:1.</figcaption></figure>

In subsequent detailed explorations of the phase space, we numerically integrated millions of initial planetary configurations to find the extension of the phase space of the Laplace resonant configurations. In Fig. 5, we show the ranges of the initial values of semimajor axis (very tiny) and eccentricity with all four planets locked in a Laplace resonance in both cases (3:2, 2:1, and 3:2 on the left and 3:2, 2:1, and 2:1 on the right). Within these ranges, there are resonant solutions that are compatible in terms of eccentricity with the observed nominal system, suggesting that the primordial system was indeed in a Laplace resonance, and it later escaped from the resonant configuration. When comparing the two different configurations, 3:2, 2:1, and 3:2 versus 3:2, 2:1, and 2:1, we find it interesting to note that the last configuration leads to higher eccentricities for the third planet.

It is worth pointing out that the resonant solutions plotted in Fig. 5 were obtained by requiring that both resonant arguments librate, a condition naturally leading to apsidal libration. This condition is in agreement with the work of Beaugé et al. (2006), who suggest that as long as the orbital migration is sufficiently slow to be approximated by an adiabatic process (which is a condition compatible with the slow migration of giant planets toward the end of their formation; see e.g., Pirani et al. 2019; Tanaka et al. 2020), all captured planets should be in apsidal co-rotations.

Even if the resonant lock provides a solution for the origin of V1298 Tau's compact architecture, the system is not in resonant chain at present, as confirmed by the extensive parameter exploration performed by Tejada Arevalo et al. (2022). The ratios between the planetary semimajor axes, according to the current nominal solutions from Table 1, are 1.312, 1.559, and 1.535,

while to be in a four-body resonance, they should be 1.310, 1.587, and 1.310 or 1.310, 1.587, and 1.587.

This mismatch can be explained if the system evolved into either of the two resonant chains in the initial stages of its formation history but the resonance lock was broken during its later evolution by some additional dynamical mechanism. According to the capture simulations with a dissipating disk (Fig. 3 and Fig. 4), it is difficult for the gas dispersal alone to break the resonance lock. Dynamically destabilizing mechanisms must be invoked, specifically the planet-planet scattering and planetesimals scattering introduced in Sect. 4.1 and explored in detail in the next sections.

### 4.4. V1298 Tau's present architecture: Resonance break by planet-planet scattering

The first possible scenario to break the resonant lock requires that additional dynamically excited planets populated, and maybe still populate, the system on orbits external to that of planet e. These excited planets would be difficult to detect through transits due to their longer orbital periods. As discussed in Sect. 4.2, the convergent migration shaping the formation path of V1298 Tau could act to promote the formation of such planets. These additional planets would have formed at later times in the wake of the passage of V1298 Tau b and e (see Sect. 4.2) and would not have migrated as close to the star as planets b and e either because of the dissipation of the gaseous disk or because they began to interact gravitationally with each other (Marzari et al. 2010). The dissipation of the disk itself can cause the onset of the dynamical instability among the resulting population of planets.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Random resonant solutions in the semimajor axis and eccentricity planet for the two different chains (3:2, 2:1, and 3:2 on the *left* and 3:2, 2:1, and 2:1 on the *right*). The solutions are scalable in semimajor axes, and the three inner planet semimajor axes can be shifted until they coincide with the observed ones.</figcaption></figure>

During their mutual gravitational interactions, one (or more) of such outer planets may have ended up in a highly eccentric orbit (Weidenschilling & Marzari 1996; Rasio & Ford 1996; Lin & Ida 1997 and Davies et al. 2014 for a review of the dynamical process and Limbach & Turner 2015; Zinzi & Turrini 2017; Turrini et al. 2020, 2022 for its signatures in observed multi-planet systems). This additional planet on an eccentric orbit may have become close enough to the four presently discovered planets while they were locked in a Laplace resonance to destabilize them. Encounters between this additional planet and planet e would break the resonance chain, leading the system of inner planets to instability on short timescales. The young age of the system suggests that the event triggering the instability of the outer planets, from which the fifth planet originated, could have been the dissipation of the gaseous disk, unless the planet formed directly on an eccentric orbit from the planetesimal disk after it was excited by the migration of planets b and e.

In Fig. 6, we show two examples of such a planet-planet scattering scenario. In the first case (top panel), a 200 $M_⊕$ planet is added to the system on an orbit with $a = 8$ au, $e = 0.95$, and $i = 10°$. The inner planets remain locked in the four-body resonance until a series of close encounters with the fifth planet breaks the resonance lock, the critical arguments begin to circulate, and the system is quickly destabilized. In the second case (bottom panel), a heavier planet with $M = 2M_J$ is placed on an orbit with a semimajor axis $a = 4$ au, eccentricity $e = 0.85$, and inclination $i = 2°$, and a similar evolution is observed.

These two examples are representative of the large sample of simulations we performed including outer planets with masses ranging from $m = 150 M_⊕$ to $m = 2 M_J$ and different initial eccentricities and semimajor axes, all leading to resonance break and instability. The first case has a very high eccentricity, which compensates the lower planet mass and the larger semimajor axis. It illustrates a scenario where the outer planetary system formed far away from the inner one and, becoming unstable, produced the perturbing planet. In the second case, the perturbing planet is closer but with a smaller eccentricity and a larger mass, representing a scenario where the inner and outer planets were potentially born more packed.

There is an infinite number of possible configurations of a perturbing planet on a crossing orbit with respect to the inner bodies leading to planet-planet scattering, chaotic evolution, and instability. Due to the chaotic nature of the problem, there is little reason to perform more extended explorations of the parameter space, particularly since our goal is solely to prove that this is a viable scenario. There is no way to identify the exact configuration that led to the present planetary system, as even those

closely reproducing the observed system would not be unique. The timescale for the onset of instability depends on the initial conditions and can be tuned by changing the initial semimajor axis of the fifth planet to delay the instability.

This behavior can explain why at present we observe V1298 Tau's four planets close to a resonance chain but not locked in it and is supported by the high value of the NAMD of V1298 Tau (see Sect. 2), which is highly suggestive of a period of violent planetary encounters (Carleo et al. 2021; Turrini et al. 2022). Unfortunately, it is not possible to trace back the initial conditions that led to the present unstable system with precision, even with a large sample of numerical simulations. The chaotic evolution due to planet-planet scattering can drive any putative system close to the observed one at different times during its evolution, and many different initial conditions can bring the system close to the observed one.

A possible alternative mechanism that may drive a system initially formed in a resonant chain to instability is tidal eccentricity damping, which may lead to divergence of the orbital semimajor axes of the resonant bodies (Batygin & Morbidelli 2013; Lithwick & Wu 2012). However, the eccentricity damping timescale of the innermost planet (planet c) is two orders of magnitude longer than the age of the system according to an estimate based on the tidal model of Leconte et al. (2010) and a modified tidal quality factor of $10^5$ for the planet. The mass loss of the innermost planets due to photoevaporation, caused by their proximity to the star, is an additional plausible candidate for triggering dynamical instabilities in a packed resonant chain (Goldberg et al. 2022). Unfortunately, to be effective, this mechanism also requires a timescale that is on the order of 100 Ma or more, which is significantly greater than the age of V1298 Tau (see Sect. 2 and Table 1). Finally, it is worth noting that the two inner planets are also very close to a 2:1 resonance, and it is the third that is far from the resonant ratio with respect to the inner two.

### 4.5. V1298 Tau's present architecture: Planetesimal scattering does not break the resonance

We also investigated the possibility that residual planetesimal scattering may be responsible for the resonant chain breaking and the current unstable configuration of the system. Remnant planetesimals may be present in the region where the planets migrate while in resonance. Furthermore, as discussed in Sect. 4.2, the migrating giant planets leave in their wake dynamically excited planetesimals whose high-velocity impacts convert a significant fraction of their mass into second-generation dust
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Random resonant solutions in the semimajor axis and eccentricity planet for the two different chains (3:2, 2:1, and 3:2 on the *left* and 3:2, 2:1, and 2:1 on the *right*). The solutions are scalable in semimajor axes, and the three inner planet semimajor axes can be shifted until they coincide with the observed ones.</figcaption></figure>

During their mutual gravitational interactions, one (or more) of such outer planets may have ended up in a highly eccentric orbit (Weidenschilling & Marzari 1996; Rasio & Ford 1996; Lin & Ida 1997 and Davies et al. 2014 for a review of the dynamical process and Limbach & Turner 2015; Zinzi & Turrini 2017; Turrini et al. 2020, 2022 for its signatures in observed multi-planet systems). This additional planet on an eccentric orbit may have become close enough to the four presently discovered planets while they were locked in a Laplace resonance to destabilize them. Encounters between this additional planet and planet e would break the resonance chain, leading the system of inner planets to instability on short timescales. The young age of the system suggests that the event triggering the instability of the outer planets, from which the fifth planet originated, could have been the dissipation of the gaseous disk, unless the planet formed directly on an eccentric orbit from the planetesimal disk after it was excited by the migration of planets b and e.

In Fig. 6, we show two examples of such a planet-planet scattering scenario. In the first case (top panel), a 200 $M_⊕$ planet is added to the system on an orbit with $a = 8$ au, $e = 0.95$, and $i = 10°$. The inner planets remain locked in the four-body resonance until a series of close encounters with the fifth planet breaks the resonance lock, the critical arguments begin to circulate, and the system is quickly destabilized. In the second case (bottom panel), a heavier planet with $M = 2M_J$ is placed on an orbit with a semimajor axis $a = 4$ au, eccentricity $e = 0.85$, and inclination $i = 2°$, and a similar evolution is observed.

These two examples are representative of the large sample of simulations we performed including outer planets with masses ranging from $m = 150 M_⊕$ to $m = 2 M_J$ and different initial eccentricities and semimajor axes, all leading to resonance break and instability. The first case has a very high eccentricity, which compensates the lower planet mass and the larger semimajor axis. It illustrates a scenario where the outer planetary system formed far away from the inner one and, becoming unstable, produced the perturbing planet. In the second case, the perturbing planet is closer but with a smaller eccentricity and a larger mass, representing a scenario where the inner and outer planets were potentially born more packed.

There is an infinite number of possible configurations of a perturbing planet on a crossing orbit with respect to the inner bodies leading to planet-planet scattering, chaotic evolution, and instability. Due to the chaotic nature of the problem, there is little reason to perform more extended explorations of the parameter space, particularly since our goal is solely to prove that this is a viable scenario. There is no way to identify the exact configuration that led to the present planetary system, as even those

closely reproducing the observed system would not be unique. The timescale for the onset of instability depends on the initial conditions and can be tuned by changing the initial semimajor axis of the fifth planet to delay the instability.

This behavior can explain why at present we observe V1298 Tau's four planets close to a resonance chain but not locked in it and is supported by the high value of the NAMD of V1298 Tau (see Sect. 2), which is highly suggestive of a period of violent planetary encounters (Carleo et al. 2021; Turrini et al. 2022). Unfortunately, it is not possible to trace back the initial conditions that led to the present unstable system with precision, even with a large sample of numerical simulations. The chaotic evolution due to planet-planet scattering can drive any putative system close to the observed one at different times during its evolution, and many different initial conditions can bring the system close to the observed one.

A possible alternative mechanism that may drive a system initially formed in a resonant chain to instability is tidal eccentricity damping, which may lead to divergence of the orbital semimajor axes of the resonant bodies (Batygin & Morbidelli 2013; Lithwick & Wu 2012). However, the eccentricity damping timescale of the innermost planet (planet c) is two orders of magnitude longer than the age of the system according to an estimate based on the tidal model of Leconte et al. (2010) and a modified tidal quality factor of $10^5$ for the planet. The mass loss of the innermost planets due to photoevaporation, caused by their proximity to the star, is an additional plausible candidate for triggering dynamical instabilities in a packed resonant chain (Goldberg et al. 2022). Unfortunately, to be effective, this mechanism also requires a timescale that is on the order of 100 Ma or more, which is significantly greater than the age of V1298 Tau (see Sect. 2 and Table 1). Finally, it is worth noting that the two inner planets are also very close to a 2:1 resonance, and it is the third that is far from the resonant ratio with respect to the inner two.

### 4.5. V1298 Tau's present architecture: Planetesimal scattering does not break the resonance

We also investigated the possibility that residual planetesimal scattering may be responsible for the resonant chain breaking and the current unstable configuration of the system. Remnant planetesimals may be present in the region where the planets migrate while in resonance. Furthermore, as discussed in Sect. 4.2, the migrating giant planets leave in their wake dynamically excited planetesimals whose high-velocity impacts convert a significant fraction of their mass into second-generation dust
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Evolution of the critical argument of the outer planet pair (b and e) initially locked in a 3:2 resonance. In the *upper panel*, the additional planet has a mass *m* = 200 *M*<sub>⊕</sub>, *a* = 8 au, *e* = 0.95, and *i* = 10°. In the *bottom panel*, a more massive planet on a less eccentric orbit was adopted (*m* = 2 *M*<sub>J</sub>, *a* = 4 au, *e* = 0.85, and inclination *i* = 2°). In both cases, after a period of chaotic evolution, the onset of repeated close encounters breaks the resonance lock and leads to fast instability.</figcaption></figure>

and pebbles (Turrini et al. 2012, 2019; Bernabò et al. 2022). While drifting toward the star, this second-generation dust can trigger new phases of planetesimal formation and create new planetesimal belts closer to the resonant planets.

We considered two scenarios characterized by planetesimal belts of 10 and 50 $M_\oplus$, respectively, both extending between 0.3 au and 10 au. This is the same orbital region where gas drag efficiently damps the dynamical excitation of the planetesimals (see Fig. 2), and as a result, the planetesimals populating the belts start their evolution after the dispersal of the circumstellar disk on low-eccentricity, low-inclination orbits ($e \approx i < 0.01$, where the inclination is expressed in radians).

We simulated the 10 and 50 $M_\oplus$ planetesimal belts using 1000 massive particles whose initial masses are 0.01 and 0.05 $M_\oplus$, respectively. Collisions among planetesimals were treated as inelastic mergers resulting in mass growth. The giant planets were placed in a stable resonant chain at the beginning of the simulations, and the whole system was evolved for 10 Ma (i.e., the lower bound of V1298 Tau's age in Table 1) with a timestep of 0.35 days. The simulations were performed using

the GPU-accelerated version of Mercury-Arχes N-body code (Turrini et al. 2019, 2021) discussed in Sect. 3.2.

The simulations revealed that the presence of the 10 $M_\oplus$ planetesimal belt has negligible effects on the dynamical evolution of the four planets, whose orbital architecture remain unaltered with respect to the reference scenario with no planetesimal belt. The presence of the 50 $M_\oplus$ planetesimal belt produces limited alterations to the dynamical evolution of the planets, as the encounters with the more massive planetesimals slightly damp the eccentricity of V1298 Tau e and reduce the amplitude of the eccentricity oscillations of V1298 Tau b. Aside from these limited differences, the secular variations of the orbital elements of the four planets are within the ranges of the scenarios with the 10 $M_\oplus$ belt and no planetesimal belt. Overall, even the presence of the 50 $M_\oplus$ planetesimal belt does not appear capable of breaking the resonance chain over the lifetime of the system. Unless V1298 Tau hosts a planetesimal belt more massive than those considered in this analysis, planetesimal scattering does not appear as a viable solution to break the resonance chain.

**5. Observational clues on the presence of additional planets at wide separation**

The scenario of high enrichment in heavy elements of V1298 Tau b and e we explore in Sect 4.1 is best reproduced if the two planets form and migrate across compact protoplanetary disks with a radial extension of <50–100 au. Wider disks imply lower spatial densities of the planetesimals, as the disk mass is spread over wider areas since, for the same total disk mass, they would have a lower coefficient of the power law describing the mass density profile. This, in turn, implies that migrating giant planets encounter fewer planetesimals and accrete lower masses of heavy elements over the same radial displacement.

Fitting the density values of planets b and e in such wide disks would require initial formation regions for the two planets similar to those of the giant planets revealed by ALMA surveys (e.g., Long et al. 2018; Andrews et al. 2018). However, it is currently unclear if such outer formation regions are capable of producing hot planets such as those around V1298 Tau. The presence of massive planets at large orbital separations, the signpost of an extended native disk, would therefore require revisiting the scenario discussed in Sect. 4.1.

Considering the young age of the system, the direct imaging technique allowed us to achieve sensitivities well into the planetary regime over a broad range of separations. Only moderately shallow data are available in the literature (Daemgen et al. 2015), so we observed V1298 Tau with SPHERE as part of a program on the characterization of the outer regions around young stars with transiting planets (Desidera et al., in prep.). The comparison of proper motions from different catalogs is also considered as an indicator of the possible presence of companions at moderately wide separations.

**5.1. Observations and data reduction**

We observed V1298 Tau three times with the SPHERE high-contrast imaging instrument at VLT (Beuzit et al. 2019). A first-epoch observation was acquired on 18 November 2019, exploiting the IRDIFS mode, then observing simultaneously with IFS (Claudi et al. 2008) in the *Y* and *J* bands (from 0.95 to 1.35 µm) and with IRDIS (Dohlen et al. 2008) in the *H* band using the H23 filter pair (wavelength of 1.593 µm and 1.667 µm for H2 and H3, respectively; Vigan et al. 2010). The second
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Evolution of the critical argument of the outer planet pair (b and e) initially locked in a 3:2 resonance. In the *upper panel*, the additional planet has a mass *m* = 200 *M*<sub>⊕</sub>, *a* = 8 au, *e* = 0.95, and *i* = 10°. In the *bottom panel*, a more massive planet on a less eccentric orbit was adopted (*m* = 2 *M*<sub>J</sub>, *a* = 4 au, *e* = 0.85, and inclination *i* = 2°). In both cases, after a period of chaotic evolution, the onset of repeated close encounters breaks the resonance lock and leads to fast instability.</figcaption></figure>

and pebbles (Turrini et al. 2012, 2019; Bernabò et al. 2022). While drifting toward the star, this second-generation dust can trigger new phases of planetesimal formation and create new planetesimal belts closer to the resonant planets.

We considered two scenarios characterized by planetesimal belts of 10 and 50 $M_\oplus$, respectively, both extending between 0.3 au and 10 au. This is the same orbital region where gas drag efficiently damps the dynamical excitation of the planetesimals (see Fig. 2), and as a result, the planetesimals populating the belts start their evolution after the dispersal of the circumstellar disk on low-eccentricity, low-inclination orbits ($e \approx i < 0.01$, where the inclination is expressed in radians).

We simulated the 10 and 50 $M_\oplus$ planetesimal belts using 1000 massive particles whose initial masses are 0.01 and 0.05 $M_\oplus$, respectively. Collisions among planetesimals were treated as inelastic mergers resulting in mass growth. The giant planets were placed in a stable resonant chain at the beginning of the simulations, and the whole system was evolved for 10 Ma (i.e., the lower bound of V1298 Tau's age in Table 1) with a timestep of 0.35 days. The simulations were performed using

the GPU-accelerated version of Mercury-Arχes *N*-body code (Turrini et al. 2019, 2021) discussed in Sect. 3.2.

The simulations revealed that the presence of the 10 $M_\oplus$ planetesimal belt has negligible effects on the dynamical evolution of the four planets, whose orbital architecture remain unaltered with respect to the reference scenario with no planetesimal belt. The presence of the 50 $M_\oplus$ planetesimal belt produces limited alterations to the dynamical evolution of the planets, as the encounters with the more massive planetesimals slightly damp the eccentricity of V1298 Tau e and reduce the amplitude of the eccentricity oscillations of V1298 Tau b. Aside from these limited differences, the secular variations of the orbital elements of the four planets are within the ranges of the scenarios with the 10 $M_\oplus$ belt and no planetesimal belt. Overall, even the presence of the 50 $M_\oplus$ planetesimal belt does not appear capable of breaking the resonance chain over the lifetime of the system. Unless V1298 Tau hosts a planetesimal belt more massive than those considered in this analysis, planetesimal scattering does not appear as a viable solution to break the resonance chain.

**5. Observational clues on the presence of additional planets at wide separation**

The scenario of high enrichment in heavy elements of V1298 Tau b and e we explore in Sect 4.1 is best reproduced if the two planets form and migrate across compact protoplanetary disks with a radial extension of <50–100 au. Wider disks imply lower spatial densities of the planetesimals, as the disk mass is spread over wider areas since, for the same total disk mass, they would have a lower coefficient of the power law describing the mass density profile. This, in turn, implies that migrating giant planets encounter fewer planetesimals and accrete lower masses of heavy elements over the same radial displacement.

Fitting the density values of planets b and e in such wide disks would require initial formation regions for the two planets similar to those of the giant planets revealed by ALMA surveys (e.g., Long et al. 2018; Andrews et al. 2018). However, it is currently unclear if such outer formation regions are capable of producing hot planets such as those around V1298 Tau. The presence of massive planets at large orbital separations, the signpost of an extended native disk, would therefore require revisiting the scenario discussed in Sect. 4.1.

Considering the young age of the system, the direct imaging technique allowed us to achieve sensitivities well into the planetary regime over a broad range of separations. Only moderately shallow data are available in the literature (Daemgen et al. 2015), so we observed V1298 Tau with SPHERE as part of a program on the characterization of the outer regions around young stars with transiting planets (Desidera et al., in prep.). The comparison of proper motions from different catalogs is also considered as an indicator of the possible presence of companions at moderately wide separations.

**5.1. Observations and data reduction**

We observed V1298 Tau three times with the SPHERE high-contrast imaging instrument at VLT (Beuzit et al. 2019). A first-epoch observation was acquired on 18 November 2019, exploiting the IRDIFS mode, then observing simultaneously with IFS (Claudi et al. 2008) in the *Y* and *J* bands (from 0.95 to 1.35 µm) and with IRDIS (Dohlen et al. 2008) in the *H* band using the H23 filter pair (wavelength of 1.593 µm and 1.667 µm for H2 and H3, respectively; Vigan et al. 2010). The second
<table><caption>**Table 4.** Main characteristics of the SPHERE observations of V1298 Tau used for this work.</caption><thead><tr><th>Date</th><th>Observing mode</th><th>Coronograph</th><th>DIMM seeing</th><th>τ<sub>0</sub></th><th>Wind speed</th><th>Field rotation</th><th>DIT</th><th>Total exposure</th></tr></thead><tbody><tr><td>2019-11-18</td><td>IRDIFS</td><td>N_ALC_YJH_S</td><td>0.56"</td><td>7.7 ms</td><td>7.9 m s<sup>-1</sup></td><td>16.3°</td><td>96 s</td><td>3072 s</td></tr><tr><td>2021-10-28</td><td>IRDIFS_EXT</td><td>N_ALC_YJH_S</td><td>1.10"</td><td>4.0 ms</td><td>5.1 m s<sup>-1</sup></td><td>15.0°</td><td>96 s</td><td>3072 s</td></tr><tr><td>2021-12-02</td><td>IRDIFS_EXT</td><td>N_ALC_YJH_S</td><td>0.56"</td><td>6.6 ms</td><td>6.4 m s<sup>-1</sup></td><td>16.5°</td><td>96 s</td><td>3072 s</td></tr></tbody></table>

<table><caption>**Table 5.** Astrometry and photometry of the candidate companion detected around V 1298 Tau.</caption><thead><tr><th>Date</th><th>Band</th><th>ρ (mas)</th><th>θ (deg)</th><th>Δ mag</th></tr></thead><tbody><tr><td>2019-11-18</td><td>H2</td><td>2768.32 ± 4.29</td><td>344.44 ± 0.13</td><td>11.74 ± 0.09</td></tr><tr><td>2019-11-18</td><td>H3</td><td>2770.44 ± 4.39</td><td>344.45 ± 0.14</td><td>11.71 ± 0.09</td></tr><tr><td>2021-10-28</td><td>K1</td><td>2798.73 ± 5.96</td><td>344.21 ± 0.30</td><td>11.62 ± 0.14</td></tr><tr><td>2021-10-28</td><td>K2</td><td>2804.62 ± 30.56</td><td>344.21 ± 1.36</td><td>11.70 ± 0.38</td></tr><tr><td>2021-12-02</td><td>K1</td><td>2794.16 ± 5.29</td><td>344.33 ± 0.26</td><td>11.52 ± 0.10</td></tr><tr><td>2021-12-02</td><td>K2</td><td>2793.30 ± 13.92</td><td>344.25 ± 0.79</td><td>11.20 ± 0.20</td></tr></tbody></table>

and third observations were obtained on 28 October 2021 and 2 December 2021, with the main goal of determining the status (physical companion versus background) of a faint candidate detected in the first epoch, as discussed in Sect. 5.2.

The follow-up observations were performed in the IRDIFS_EXT mode, then observing simultaneously with IFS in *Y*, *J*, and *H* bands (from 0.95 to 1.65 µm) and with IRDIS in the *K* band using the K12 filter pair (wavelength of 2.110 µm and 2.251 µm for K1 and K2 bands, respectively). We chose a different setup for the follow-up in order to obtain complementary photometric measurements for the candidate and achieve a better sensitivity for planets with very dusty atmospheres (see, e.g., Chauvin et al. 2018). The characteristics of the three datasets are summarized in Table 4.

The data were reduced through the SPHERE Data Center (Delorme et al. 2017) following the SPHERE Data Reduction and Handling pipeline (Pavlov et al. 2008) and applying the appropriate calibrations for our datasets. We then applied speckle subtraction algorithms TLOCI (Marois et al. 2014) and principal components analysis (PCA; Soummer et al. 2012) on the reduced data as implemented in the SpeCal pipeline (Galicher et al. 2018).

## 5.2. A background object projected close to V1298 Tau

A point source was detected in all the three epochs in the IRDIS field of view at a separation of about 2.77" (Fig. 7). The source has a contrast higher than 11 magnitudes with respect to the central star. The signal-to-noise ratio of the detection is 35.2, 6.8, and 9.4 for the first, second, and third epochs, respectively. This difference is reflected by the higher error bars in the astrometry and the photometry of the candidate listed in Table 5. The astrometric calibration of each epoch was obtained following the method devised by Maire et al. (2016). The photometry was calculated using the negative planet method as described in Bonnefoy et al. (2011) and in Zurlo et al. (2014), for example.

We checked the NIRI images of V1298 Tau obtained by Daemgen et al. (2015), and we confirmed that the source is not detected in this dataset, as expected from the published detection limits. We compared the astrometry in the first and the third epoch, adopting the stellar parameters listed in Table 1. We

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** SPHERE image of V1298 Tau in the H2 band with the faint background object at 2.77" highlighted with the red arrow.</figcaption></figure>

decided to exclude the second epoch from this analysis as, in that epoch, the candidate companion is just above the detection limit and the astrometric error bars are very large as a consequence.

This comparison is displayed in Fig. 8, where the green square represents the position of the candidate relative to the star in the first epoch, while the orange diamond gives the relative position of the candidate in the third epoch. The solid black line represents the course of the candidate during the considered period if it were a background object with no proper motion. The black square at the end of this line is the expected position of the candidate in this latter case at the third epoch. This plot excludes the fact that the candidate could be gravitationally bound to the central star and points to it being a background object possibly possessing a non-negligible proper motion.

## 5.3. Constraints on additional planets

No additional sources were detected either with IFS or IRDIS. In order to quantify our detection limits, we defined the contrast
<table><caption>**Table 4.** Main characteristics of the SPHERE observations of V1298 Tau used for this work.</caption><thead><tr><th>Date</th><th>Observing mode</th><th>Coronograph</th><th>DIMM seeing</th><th>τ<sub>0</sub></th><th>Wind speed</th><th>Field rotation</th><th>DIT</th><th>Total exposure</th></tr></thead><tbody><tr><td>2019-11-18</td><td>IRDIFS</td><td>N_ALC_YJH_S</td><td>0.56"</td><td>7.7 ms</td><td>7.9 m s<sup>-1</sup></td><td>16.3°</td><td>96 s</td><td>3072 s</td></tr><tr><td>2021-10-28</td><td>IRDIFS_EXT</td><td>N_ALC_YJH_S</td><td>1.10"</td><td>4.0 ms</td><td>5.1 m s<sup>-1</sup></td><td>15.0°</td><td>96 s</td><td>3072 s</td></tr><tr><td>2021-12-02</td><td>IRDIFS_EXT</td><td>N_ALC_YJH_S</td><td>0.56"</td><td>6.6 ms</td><td>6.4 m s<sup>-1</sup></td><td>16.5°</td><td>96 s</td><td>3072 s</td></tr></tbody></table>

<table><caption>**Table 5.** Astrometry and photometry of the candidate companion detected around V 1298 Tau.</caption><thead><tr><th>Date</th><th>Band</th><th>ρ<br>(mas)</th><th>θ<br>(deg)</th><th>Δ mag</th></tr></thead><tbody><tr><td>2019-11-18</td><td>H2</td><td>2768.32 ± 4.29</td><td>344.44 ± 0.13</td><td>11.74 ± 0.09</td></tr><tr><td>2019-11-18</td><td>H3</td><td>2770.44 ± 4.39</td><td>344.45 ± 0.14</td><td>11.71 ± 0.09</td></tr><tr><td>2021-10-28</td><td>K1</td><td>2798.73 ± 5.96</td><td>344.21 ± 0.30</td><td>11.62 ± 0.14</td></tr><tr><td>2021-10-28</td><td>K2</td><td>2804.62 ± 30.56</td><td>344.21 ± 1.36</td><td>11.70 ± 0.38</td></tr><tr><td>2021-12-02</td><td>K1</td><td>2794.16 ± 5.29</td><td>344.33 ± 0.26</td><td>11.52 ± 0.10</td></tr><tr><td>2021-12-02</td><td>K2</td><td>2793.30 ± 13.92</td><td>344.25 ± 0.79</td><td>11.20 ± 0.20</td></tr></tbody></table>

and third observations were obtained on 28 October 2021 and 2 December 2021, with the main goal of determining the status (physical companion versus background) of a faint candidate detected in the first epoch, as discussed in Sect. 5.2.

The follow-up observations were performed in the IRDIFS_EXT mode, then observing simultaneously with IFS in *Y*, *J*, and *H* bands (from 0.95 to 1.65 µm) and with IRDIS in the *K* band using the K12 filter pair (wavelength of 2.110 µm and 2.251 µm for K1 and K2 bands, respectively). We chose a different setup for the follow-up in order to obtain complementary photometric measurements for the candidate and achieve a better sensitivity for planets with very dusty atmospheres (see, e.g., Chauvin et al. 2018). The characteristics of the three datasets are summarized in Table 4.

The data were reduced through the SPHERE Data Center (Delorme et al. 2017) following the SPHERE Data Reduction and Handling pipeline (Pavlov et al. 2008) and applying the appropriate calibrations for our datasets. We then applied speckle subtraction algorithms TLOCI (Marois et al. 2014) and principal components analysis (PCA; Souummer et al. 2012) on the reduced data as implemented in the SpeCal pipeline (Galicher et al. 2018).

## 5.2. A background object projected close to V1298 Tau

A point source was detected in all the three epochs in the IRDIS field of view at a separation of about 2.77" (Fig. 7). The source has a contrast higher than 11 magnitudes with respect to the central star. The signal-to-noise ratio of the detection is 35.2, 6.8, and 9.4 for the first, second, and third epochs, respectively. This difference is reflected by the higher error bars in the astrometry and the photometry of the candidate listed in Table 5. The astrometric calibration of each epoch was obtained following the method devised by Maire et al. (2016). The photometry was calculated using the negative planet method as described in Bonnefoy et al. (2011) and in Zurlo et al. (2014), for example.

We checked the NIRI images of V1298 Tau obtained by Daemgen et al. (2015), and we confirmed that the source is not detected in this dataset, as expected from the published detection limits. We compared the astrometry in the first and the third epoch, adopting the stellar parameters listed in Table 1. We

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** SPHERE image of V1298 Tau in the H2 band with the faint background object at 2.77" highlighted with the red arrow.</figcaption></figure>

decided to exclude the second epoch from this analysis as, in that epoch, the candidate companion is just above the detection limit and the astrometric error bars are very large as a consequence.

This comparison is displayed in Fig. 8, where the green square represents the position of the candidate relative to the star in the first epoch, while the orange diamond gives the relative position of the candidate in the third epoch. The solid black line represents the course of the candidate during the considered period if it were a background object with no proper motion. The black square at the end of this line is the expected position of the candidate in this latter case at the third epoch. This plot excludes the fact that the candidate could be gravitationally bound to the central star and points to it being a background object possibly possessing a non-negligible proper motion.

## 5.3. Constraints on additional planets

No additional sources were detected either with IFS or IRDIS. In order to quantify our detection limits, we defined the contrast
<figure><img src="image_9.png" /></figure>

**Fig. 8.** Relative astrometric position of the proposed companion in the first and the last observing epoch. The solid black line represents the expected course of the companion if it were a stationary background object. The black square at the end of the line represents the expected position at the epoch of the last observation in this case. The minor difference between the observed and expected positions is consistent at ~1.5σ with measurement errors or could be due to a non-zero intrinsic motion of the background source.

around the central star for both instruments and for the first and third epochs, exploiting the procedure described in Mesa et al. (2015) and corrected for the small sample statistic following the method described by Mawet et al. (2014). From these contrast limits, we calculated the upper mass limits for companions around V 1298 Tau using the AMES-COND models (Allard et al. 2003) and adopting the stellar parameters by Suárez Mascareño et al. (2022). Finally, we defined the best upper limits for all three epochs and considering both IFS and IRDIS mass limits. The final results of this procedure are displayed in Fig. 9.

The present constraints rule out the presence of giant planets more massive than Jupiter beyond 50-100 au (i.e., at orbital periods greater than 300-1000 years), suggesting that the circumstellar disk from which V1298 Tau's planetary system formed was more compact and dense than those presently being observed by ALMA surveys (e.g., Andrews et al. 2018; Long et al. 2018). Such a compact disk fits the picture discussed in Sect. 4.1, as it favors the accretion of large amounts of planetesimals by the migrating planets and the production of large amounts of collisional dust from the surviving planetesimal disk.

Furthermore, the present constraints rule out only giant planets more massive than a few Jovian masses in the inner regions within 50 au. This means that between 10 and 50 au (i.e., at orbital periods between 30 and 300 yr), there may well be planets as massive as $2M_J$ on eccentric or unstable orbits that are below the detection limit of SPHERE (as discussed in Sect. 5.4, such planets would fall near or below the 50% iso-probability curve of Gaia's astrometric observations). Among these planets, the one that possibly destabilized the inner four may be still orbiting V1298 Tau on a highly eccentric orbit, unless it was ejected from the system by the planet-planet scattering process that triggered the instability.

## 5.4. Additional clues from astrometry

The star V1298 Tau was not observed by the HIPPARCOS space astrometry mission, so it is not possible to provide direct constraints on long-period orbiting companions in a regime of orbital separations complementary to the one probed by the

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

**Fig. 9.** Mass detection limits (5σ) obtained from the SPHERE data (both IRDIS and IFS) expressed in units of Jovian mass as a function of the separation expressed in astronomical units.

SPHERE direct imaging observations ($a \lesssim 10$ au) using the HIPPARCOS-Gaia proper motion anomaly technique (see, e.g., Kervella et al. 2019, 2022; Brandt 2018, 2021). However, Gaia DR3 provides the renormalized unit weight error (RUWE) statistics, which is a good proxy for the quality of the single-star astrometric solution. Typically, RUWE ≈ 1.0 indicates a good-quality single-star solution, while sources with RUWE above the threshold value ≈1.4 are typically considered to have observations inconsistent with the astrometric five-parameter model, with additional variability in the astrometry possibly due to binarity (e.g., Lindegren et al. 2018, 2021).

The reported RUWE value for V1298 Tau in Gaia DR3 is 1.035, which indicates no significant departures from a single-star model are observed, and in fact no non-single-star solution of any type has been reported. Following the approach by Belokurov et al. (2020) and Penoyre et al. (2020), for example, it is nevertheless possible to explore the regime of companion masses and orbital separations that can be excluded, as they would have produced larger RUWE values than the reported one.

We set up a numerical simulation, creating synthetic Gaia observations of V1298 Tau using the nominal astrometric parameters published in Gaia DR3. We used the Gaia observation forecast tool<fnref n="4" /> to obtain a close representation of the actual Gaia observation times, scan angle, and along-scan parallax factors encompassing the mission time span utilized in Gaia DR3. Orbital motion effects were linearly superposed considering companions with orbital periods in the range 0.5–20.0 yr and masses in the range 1.0–40.0 $M_J$.

For each period-mass pair, 100 random realizations of the remaining orbital elements were produced, all being drawn from uniform distributions within their nominal intervals. The Gaia-like observations of each star plus companion system were then perturbed with Gaussian measurement uncertainties appropriate for the case of a $G = 10$ mag star such as V1298 Tau (see, e.g., Holl et al. 2023, Fig. 3). Each time series was then fitted with a single-star model, and the corresponding RUWE value was recorded. For each period-mass pair, the fraction of systems with RUWE values exceeding the one reported for V1298 Tau in the Gaia DR3 archive was recorded. In total, the simulation was run for 4 000 000 systems.

<footnote n="4">⁴ https://gaia.esac.esa.int/gost/</footnote>
<figure><img src="image_9.png" /></figure>

**Fig. 8.** Relative astrometric position of the proposed companion in the first and the last observing epoch. The solid black line represents the expected course of the companion if it were a stationary background object. The black square at the end of the line represents the expected position at the epoch of the last observation in this case. The minor difference between the observed and expected positions is consistent at ~1.5σ with measurement errors or could be due to a non-zero intrinsic motion of the background source.

around the central star for both instruments and for the first and third epochs, exploiting the procedure described in Mesa et al. (2015) and corrected for the small sample statistic following the method described by Mawet et al. (2014). From these contrast limits, we calculated the upper mass limits for companions around V 1298 Tau using the AMES-COND models (Allard et al. 2003) and adopting the stellar parameters by Suárez Mascareño et al. (2022). Finally, we defined the best upper limits for all three epochs and considering both IFS and IRDIS mass limits. The final results of this procedure are displayed in Fig. 9.

The present constraints rule out the presence of giant planets more massive than Jupiter beyond 50-100 au (i.e., at orbital periods greater than 300-1000 years), suggesting that the circumstellar disk from which V1298 Tau's planetary system formed was more compact and dense than those presently being observed by ALMA surveys (e.g., Andrews et al. 2018; Long et al. 2018). Such a compact disk fits the picture discussed in Sect. 4.1, as it favors the accretion of large amounts of planetesimals by the migrating planets and the production of large amounts of collisional dust from the surviving planetesimal disk.

Furthermore, the present constraints rule out only giant planets more massive than a few Jovian masses in the inner regions within 50 au. This means that between 10 and 50 au (i.e., at orbital periods between 30 and 300 yr), there may well be planets as massive as $2M_J$ on eccentric or unstable orbits that are below the detection limit of SPHERE (as discussed in Sect. 5.4, such planets would fall near or below the 50% iso-probability curve of Gaia's astrometric observations). Among these planets, the one that possibly destabilized the inner four may be still orbiting V1298 Tau on a highly eccentric orbit, unless it was ejected from the system by the planet-planet scattering process that triggered the instability.

## 5.4. Additional clues from astrometry

The star V1298 Tau was not observed by the HIPPARCOS space astrometry mission, so it is not possible to provide direct constraints on long-period orbiting companions in a regime of orbital separations complementary to the one probed by the

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

**Fig. 9.** Mass detection limits (5σ) obtained from the SPHERE data (both IRDIS and IFS) expressed in units of Jovian mass as a function of the separation expressed in astronomical units.

SPHERE direct imaging observations ($a \lesssim 10$ au) using the HIPPARCOS-Gaia proper motion anomaly technique (see, e.g., Kervella et al. 2019, 2022; Brandt 2018, 2021). However, Gaia DR3 provides the renormalized unit weight error (RUWE) statistics, which is a good proxy for the quality of the single-star astrometric solution. Typically, RUWE ≈ 1.0 indicates a good-quality single-star solution, while sources with RUWE above the threshold value ≈1.4 are typically considered to have observations inconsistent with the astrometric five-parameter model, with additional variability in the astrometry possibly due to binarity (e.g., Lindegren et al. 2018, 2021).

The reported RUWE value for V1298 Tau in Gaia DR3 is 1.035, which indicates no significant departures from a single-star model are observed, and in fact no non-single-star solution of any type has been reported. Following the approach by Belokurov et al. (2020) and Penoyre et al. (2020), for example, it is nevertheless possible to explore the regime of companion masses and orbital separations that can be excluded, as they would have produced larger RUWE values than the reported one.

We set up a numerical simulation, creating synthetic Gaia observations of V1298 Tau using the nominal astrometric parameters published in Gaia DR3. We used the Gaia observation forecast tool<fnref n="4" /> to obtain a close representation of the actual Gaia observation times, scan angle, and along-scan parallax factors encompassing the mission time span utilized in Gaia DR3. Orbital motion effects were linearly superposed considering companions with orbital periods in the range 0.5–20.0 yr and masses in the range 1.0–40.0 $M_J$.

For each period-mass pair, 100 random realizations of the remaining orbital elements were produced, all being drawn from uniform distributions within their nominal intervals. The Gaia-like observations of each star plus companion system were then perturbed with Gaussian measurement uncertainties appropriate for the case of a $G = 10$ mag star such as V1298 Tau (see, e.g., Holl et al. 2023, Fig. 3). Each time series was then fitted with a single-star model, and the corresponding RUWE value was recorded. For each period-mass pair, the fraction of systems with RUWE values exceeding the one reported for V1298 Tau in the Gaia DR3 archive was recorded. In total, the simulation was run for 4 000 000 systems.

<footnote n="4">⁴ https://gaia.esac.esa.int/gost/</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** *Gaia* DR3 sensitivity to companions of a given mass (in units of Jovian mass) as a function of the orbital period (in years). Solid, dashed, and dashed-dotted lines correspond to iso-probability curves for 50%, 90%, and 99% probability of a companion of given properties to produce a RUWE value that exceeds the value reported in the *Gaia* DR3 archive (1.035).</figcaption></figure>

Figure 10 shows the iso-probability contours in companion period-mass space corresponding to different fractions of systems for which RUWE was recorded to be larger than the value reported in the *Gaia* DR3 archive for V1298 Tau. Given the time span and quality of *Gaia* DR3 observations for the star, we deduced that the *Gaia* sensitivity to super-Jupiter-type companions around V1298 Tau is presently limited to a narrow regime of orbital periods in the approximate interval of 1.5–4 yr, or roughly 1.4–2.7 au. This is the regime of orbital separations for which *Gaia* DR3-level astrometry provides mass sensitivity comparable to that achieved by SPHERE’s direct imaging data at 10–20 au.

The astrometric signals of the much shorter-period giants V1298 Tau b and e are clearly out of reach for *Gaia*, as expected. In the approximate period interval 5–10 yr, the presence of low-mass brown dwarfs can be safely ruled out. At *P* = 20 yr, a 20-*M*<sub>J</sub> companion still has a 50% chance of going undetected. As in the case of SPHERE’s observations, *Gaia’s* astrometric data leave open the possibility that planets as massive as 2*M*<sub>J</sub> linked to the destabilization of the primordial resonant chain of V1298 Tau’s four inner planets still orbit the host star on eccentric or unstable orbits below the detection limits.

5.5. Additional clues from radial velocities

To further constrain the presence of additional companions, we investigated the detection limits from the RV time series from Suárez Mascareño et al. (2022) through the Bayesian procedure adopted in Pinamonti et al. (2022). This technique consists of modeling an additional planetary signal in the RV time series by means of the publicly available *emcee* Affine Invariant MCMC Ensemble sampler (Foreman-Mackey et al. 2013) and using the posterior distribution of its orbital parameters to derive the detectability function.

The additional planet was modeled simultaneously to the
other planetary signals present in the data, both to account for
their uncertainties and to accurately estimate its detectability.
However, the strong stellar activity contamination in the RV hin-
ders a straightforward computation of the detection function of

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Detection function map of the RV time series of V1298 Tau from Suárez Mascareño et al. (2022), with the stellar activity already removed (i.e., the optimistic case). The color scale expresses the detection function, while the red dashed lines correspond to the 90% detection threshold including the stellar activity in the modeling (i.e., the pessimistic case). The blue circles correspond to the position of the measured planets in the systems, while the blue lines mark the upper limits of the planets with no mass-measurement yet. The left panel shows the Bayesian detection map computed for *P* < 400 days, while the right panel shows the inject-recovery detection map computed for 400 < *P* < 4000 days.</figcaption></figure>

the time series. This issue can be addressed by adopting the RV
residuals from Suárez Mascareño et al. (2022), as the activity
have already been removed, or by including it in the analysis via
GP regression (see Pinamonti et al. 2022). The first approach
produces an optimistic detection threshold, while the second
yields a pessimistic one, as illustrated in the left panel of Fig. 11.

The detection function can be reliably computed with the Bayesian technique only for periods up to the time span of the RV observations (≈400 days). To extend the search to longer periods, comparable with those investigated in Sects. 5.3 and 5.4, we performed an injection-recovery simulation in which synthetic planetary signals were injected in the RV residuals and the detection function was computed as the recovery rate of these signals, covering orbital periods between 400–4000 days (i.e., approximately between 1 and 5 au). The simulated signals were fit with either a circular-Keplerian orbit or a linear or quadratic trend (because a period longer than the time span might not be properly modeled by a Keplerian). We adopted the Bayesian information criterion (BIC; Schwarz 1978) in order to compare the fitted planetary model with a constant model with a threshold of ΔBIC > 10. The resulting detection function map is shown in the right panel of Fig. 11.

As can be seen from Fig. 11, in the optimistic case and con-
sidering the best-case correction of the stellar activity, we can
rule out the presence of companions of 0.5*M*<sub>J</sub> at periods shorter
than *P* = 650 d, and we can rule out the presence of additional
companions as massive as 2*M*<sub>J</sub> up to *P* = 1500 d. However,
when considering the pessimistic case and taking into account
the effect of the stellar noise in the RV time series, we can
exclude only the presence of 2*M*<sub>J</sub> companions on short-period
orbits (*P* ≲ 3 d). The real detection limit of V1298 Tau likely
falls between these two scenarios. As a result, while we can
somewhat confidently rule out the presence of massive compan-
ions at orbits close to those of the transiting planets, it is difficult
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** *Gaia* DR3 sensitivity to companions of a given mass (in units of Jovian mass) as a function of the orbital period (in years). Solid, dashed, and dashed-dotted lines correspond to iso-probability curves for 50%, 90%, and 99% probability of a companion of given properties to produce a RUWE value that exceeds the value reported in the *Gaia* DR3 archive (1.035).</figcaption></figure>

Figure 10 shows the iso-probability contours in companion period-mass space corresponding to different fractions of systems for which RUWE was recorded to be larger than the value reported in the *Gaia* DR3 archive for V1298 Tau. Given the time span and quality of *Gaia* DR3 observations for the star, we deduced that the *Gaia* sensitivity to super-Jupiter-type companions around V1298 Tau is presently limited to a narrow regime of orbital periods in the approximate interval of 1.5–4 yr, or roughly 1.4–2.7 au. This is the regime of orbital separations for which *Gaia* DR3-level astrometry provides mass sensitivity comparable to that achieved by SPHERE’s direct imaging data at 10–20 au.

The astrometric signals of the much shorter-period giants V1298 Tau b and e are clearly out of reach for *Gaia*, as expected. In the approximate period interval 5–10 yr, the presence of low-mass brown dwarfs can be safely ruled out. At *P* = 20 yr, a 20-*M*<sub>J</sub> companion still has a 50% chance of going undetected. As in the case of SPHERE’s observations, *Gaia’s* astrometric data leave open the possibility that planets as massive as 2*M*<sub>J</sub> linked to the destabilization of the primordial resonant chain of V1298 Tau’s four inner planets still orbit the host star on eccentric or unstable orbits below the detection limits.

5.5. Additional clues from radial velocities

To further constrain the presence of additional companions, we investigated the detection limits from the RV time series from Suárez Mascareño et al. (2022) through the Bayesian procedure adopted in Pinamonti et al. (2022). This technique consists of modeling an additional planetary signal in the RV time series by means of the publicly available *emcee* Affine Invariant MCMC Ensemble sampler (Foreman-Mackey et al. 2013) and using the posterior distribution of its orbital parameters to derive the detectability function.

The additional planet was modeled simultaneously to the
other planetary signals present in the data, both to account for
their uncertainties and to accurately estimate its detectability.
However, the strong stellar activity contamination in the RV hin-
ders a straightforward computation of the detection function of

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Detection function map of the RV time series of V1298 Tau from Suárez Mascareño et al. (2022), with the stellar activity already removed (i.e., the optimistic case). The color scale expresses the detection function, while the red dashed lines correspond to the 90% detection threshold including the stellar activity in the modeling (i.e., the pessimistic case). The blue circles correspond to the position of the measured planets in the systems, while the blue lines mark the upper limits of the planets with no mass-measurement yet. The left panel shows the Bayesian detection map computed for *P* < 400 days, while the right panel shows the inject-recovery detection map computed for 400 < *P* < 4000 days.</figcaption></figure>

the time series. This issue can be addressed by adopting the RV
residuals from Suárez Mascareño et al. (2022), as the activity
have already been removed, or by including it in the analysis via
GP regression (see Pinamonti et al. 2022). The first approach
produces an optimistic detection threshold, while the second
yields a pessimistic one, as illustrated in the left panel of Fig. 11.

The detection function can be reliably computed with the Bayesian technique only for periods up to the time span of the RV observations (≈400 days). To extend the search to longer periods, comparable with those investigated in Sects. 5.3 and 5.4, we performed an injection-recovery simulation in which synthetic planetary signals were injected in the RV residuals and the detection function was computed as the recovery rate of these signals, covering orbital periods between 400–4000 days (i.e., approximately between 1 and 5 au). The simulated signals were fit with either a circular-Keplerian orbit or a linear or quadratic trend (because a period longer than the time span might not be properly modeled by a Keplerian). We adopted the Bayesian information criterion (BIC; Schwarz 1978) in order to compare the fitted planetary model with a constant model with a threshold of ΔBIC > 10. The resulting detection function map is shown in the right panel of Fig. 11.

As can be seen from Fig. 11, in the optimistic case and con-
sidering the best-case correction of the stellar activity, we can
rule out the presence of companions of 0.5*M<sub>J</sub>* at periods shorter
than *P* = 650 d, and we can rule out the presence of additional
companions as massive as 2*M<sub>J</sub>* up to *P* = 1500 d. However,
when considering the pessimistic case and taking into account
the effect of the stellar noise in the RV time series, we can
exclude only the presence of 2*M<sub>J</sub>* companions on short-period
orbits (*P* ≲ 3 d). The real detection limit of V1298 Tau likely
falls between these two scenarios. As a result, while we can
somewhat confidently rule out the presence of massive compan-
ions at orbits close to those of the transiting planets, it is difficult
to accurately constrain the presence of longer-period planetary
signals underneath the stellar RV noise.

# 6. Discussion and conclusions

The recent observational campaigns aimed at characteriz-
ing the planets around V1298 Tau (David et al. 2019a,b;
Suárez Mascareño et al. 2022; Feinstein et al. 2022; Damasso
et al. 2023) reveal that this young planetary system is highly
peculiar. Its current architecture is quite compact, yet its planets
do not appear to be globally locked in resonance Tejada Arevalo
et al. (2022), and their orbital characteristics argue for the system
having crossed phases of dynamical instability (Sect. 2.1).

The dynamical state of the system, however, is not the sole
peculiar aspect of V1298 Tau. The estimated masses of its two
outer giant planets are associated with anomalously high den-
sity values for such young planetary objects. The confirmation
of its currently assessed physical characteristics would make
V1298 Tau a unique natural laboratory to study the formation
and evolution of planetary systems as well as a compelling target
for future physical and atmospheric characterization campaigns.

Specifically, unless young giant planets contract faster than
previously thought (Suárez Mascareño et al. 2022), the current
densities of V1298 Tau b and e point to their formation tracks as
being shaped by large-scale migration and enrichment in heavy
elements (Thorngren et al. 2016; Shibata et al. 2020; Turrini
et al. 2021; Pacetti et al. 2022). Consequently, reconstructing the
formation history of V1298 Tau must account for both the high
enrichment of heavy elements in its outer two planets and the
non-resonant state of the compact dynamical configuration of
the system.

In this work, we outlined the possible sequence of events that
drove the system into its present configuration. Specifically, we
argued that the convergent, large-scale migration of V1298 Tau's
planets while embedded in their native circumstellar disk can
explain both the density values of the two outermost giant planets
and the formation of V1298 Tau's compact architecture by res-
onant trapping. The dynamically excited planetesimal disk left
behind by the migrating giant planets can then support the for-
mation of additional planets that can break the original resonant
chain.

## 6.1. Assembling the mosaic of V1298 Tau's formation history

The uncertainty affecting the mass and density estimates of V1298 Tau's planets make it difficult to precisely pinpoint their formation regions. Optimistic estimates require about 20 and 30 $M_\oplus$ of heavy elements for V1298 Tau b and e, respectively. Accounting for the lower density of the hot gas composing such young giant planets brings the masses of heavy elements to 90 and 110 $M_\oplus$, respectively, making these planets some of the most metallic gas giants currently known. Adopting the nominal mass of V1298 Tau e would result in a planet whose composition is two-thirds heavy elements, and therefore more similar to Neptune than Jupiter, notwithstanding its Jovian mass.

Explaining the nominal density of V1298 Tau e in the frame-
work of the migration and enrichment scenario would require
the planet to have undergone extreme migration within a mas-
sive circumstellar disk. The intense planetesimal flux hitting the
forming planet would likely hinder the onset of its runaway gas
accretion process. This would keep the planet in the form of
a massive core surrounded by an extended atmosphere (unless

the primordial atmosphere is collisionally stripped) until the
planetesimal flux drops and the gas accretion process can restart.

A similar scenario is discussed for Jupiter by Alibert et al. (2018), but the characteristics of V1298 Tau e would invoke a significantly more extreme version than its original formulation. This scenario would make V1298 Tau e the first example of a new class of outcomes of the planet formation process. Due to the large uncertainty on its mass, however, it is more plausible that the true density of V1298 Tau e is lower than its current nominal value. Alternatively, the high density can be explained by invoking the interplay between planetesimal accretion and the accretion of disk gas enriched in heavy elements (Booth & Ilee 2019; Cridland et al. 2019; Schneider & Bitsch 2021).

Notwithstanding these uncertainties, our results describe an
overall coherent picture. Our simulations of the accretion history
of V1298 Tau b and e show that to reproduce their enrichment
in heavy elements, particularly that of planet e, their native cir-
cumstellar disk was plausibly massive, at least as massive as 0.1
times the host star. While massive, such a disk should have been
compact, with radial extension not exceeding 50–100 au. This is
consistent with the lack of massive companions at large separa-
tions that we verified by combining SPHERE’s direct imaging
and Gaia’s astrometric data.

Both planets must have crossed significant fractions of this
disk to accrete enough heavy elements to fit the current data,
as when adopting even the lowest metallicity values, their for-
mation region lies beyond Jupiter's current orbit. In all of our
simulations, V1298 Tau e should have started its formation about
twice as far as V1298 Tau b to be able to encounter and accrete
enough planetesimals. Depending on their real density values
and the interplay between planetesimal and enriched-gas accre-
tion in delivering heavy elements, their migration could have
spanned from a few astronomical units to several tens of astro-
nomical units. This uncertainty showcases the limited diagnostic
power of planetary density and metallicity (see Turrini et al. 2022
for more discussion) but also makes these planets compelling tar-
gets for future atmospheric characterization, as we discuss below
(see Sect. 6.2).

Fitting V1298 Tau's current compact architecture with the
large-scale migration of its two outer planets requires that the
system formed by convergent migration and resonant trapping
(i.e., the planets were trapped in a sequence of mean motion
resonances while migrating inwards due to disk-planet interac-
tions). This is confirmed by our simulations showing how the
planets are captured in a Laplace resonance in the presence of
dissipative forces such as those caused by their interactions with
the native circumstellar disk. The resonance trapping mechanism
naturally explains why at present V1298 Tau's four planets are
very close to a resonant configuration. Moreover, their conver-
gent migration promotes the conditions for the late formation of
massive planets in their wake, an important piece to explain the
later evolution of the system.

Resonant capture can account for the high eccentricity of
the planets, which will be excited during the dissipative reso-
nant capture. However, since the present system is only close to
but not in a Laplace resonance (Tejada Arevalo et al. 2022) and
also in a dynamically excited state, some additional mechanism
is needed to break the resonance lock after its formation. Thanks
to its well-constrained young age, V1298 Tau is an almost unique
laboratory for exploring different resonance-breaking mecha-
nisms. Specifically, V1298 Tau's age allows for the exclusion of
both the mechanism proposed by Batygin & Morbidelli (2013)
and a possible role of planetary evaporation (see Goldberg et al.
2022).