GJ 273: on the formation, dynamical evolution, and habitability of a planetary system hosted by an M dwarf at 3.75 parsec (original) (raw)

A&A 641, A23 (2020)

1,2, Juan C. Suárez3,4, Gonzalo C. de Elía5,6, Zaira M. Berdiñas7, Andrea Bonfanti1,9, Agustín Dugaro5,6, Michaël Gillon2, Emmanuël Jehin1, Maximilian N. Günther8,★, Valérie Van Grootel1, Lionel J. Garcia2, Antoine Thuillier1,10, Laetitia Delrez1,2 and Jose R. Rodón4

1Space sciences, Technologies & Astrophysics Research (STAR) Institute, Université de Liège, 4000 Liège, Belgium
e-mail: fjpozuelos@uliege.be
2EXOTIC Lab, UR Astrobiology, AGO Department, University of Liège, 4000 Liège, Belgium
3Dpt. Física Teórica y del Cosmos, Universidad de Granada, Campus de Fuentenueva s/n, 18071 Granada, Spain
4 Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain
5 Instituto de Astrofísica de La Plata, La Plata, Argentina
6Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, LaPlata, Argentina
7Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Casilla 36-D, Santiago, Chile
8Department of Physics, and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
9Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042 Graz, Austria
10Observatoire de Paris, 61 Avenue de l’Observatoire 75014 Paris, France

Received: 29 March 2020
Accepted: 12 June 2020

Abstract

Context. Planets orbiting low-mass stars such as M dwarfs are now considered a cornerstone in the search for planets with the potential to harbour life. GJ 273 is a planetary system orbiting an M dwarf only 3.75 pc away, which is composed of two confirmed planets, GJ 273b and GJ 273c, and two promising candidates, GJ 273d and GJ 273e. Planet GJ 273b resides in the habitable zone. Currently, due to a lack of observed planetary transits, only the minimum masses of the planets are known: _M_b sin _i_b = 2.89 _M_⊕, _M_c sin _i_c = 1.18 _M_⊕, _M_d sin _i_d = 10.80 _M_⊕, and _M_e sin _i_e = 9.30 _M_⊕. Despite its interesting character, the GJ 273 planetary system has been poorly studied thus far.

Aims. We aim to precisely determine the physical parameters of the individual planets, in particular, to break the mass–inclination degeneracy to accurately determine the mass of the planets. Moreover, we present a thorough characterisation of planet GJ 273b in terms of its potential habitability.

Methods. First, we explored the planetary formation and hydration phases of GJ 273 during the first 100 Myr. Secondly, we analysed the stability of the system by considering both the two- and four-planet configurations. We then performed a comparative analysis between GJ 273 and the Solar System and we searched for regions in GJ 273 which may harbour minor bodies in stable orbits, that is, the main asteroid belt and Kuiper belt analogues.

Results. From our set of dynamical studies, we find that the four-planet configuration of the system allows us to break the mass–inclination degeneracy. From our modelling results, the masses of the planets are unveiled as: 2.89 ≤ _M_b ≤ 3.03 _M_⊕, 1.18 ≤ _M_c ≤ 1.24 _M_⊕, 10.80 ≤ _M_d ≤ 11.35 _M_⊕, and 9.30 ≤ _M_e ≤ 9.70 _M_⊕. These results point to a system that is likely to be composed of an Earth-mass planet, a super-Earth and two mini-Neptunes. Based on planetary formation models, we determine that GJ 273b is likely an efficient water captor while GJ 273c is probably a dry planet. We find that the system may have several stable regions where minor bodies might reside. Collectively, these results are used to offer a comprehensive discussion about the habitability of GJ 273b.

Key words: planets and satellites: dynamical evolution and stability / planets and satellites: formation

★

Juan Carlos Torres Fellow.

© ESO 2020

1 Introduction

In the last decade, interest in planetary systems around M-type stars has risen significantly. The source of this scientific interest is twofold. On the one hand, M-type stars are the most abundant star type in the Solar neighbourhood and, indeed, in the Universe (see e.g. Henry et al. 1994; Chabrier & Baraffe 2000; Winters et al. 2015). However, their physical properties are poorly known since they are generally quite faint, as compared with other stellar types, and they have only recently become a cornerstone for the detection of low-mass planets. Indeed, recent studies based on radial–velocity (RV) and transit photometry surveys have hinted at the high occurrence of low-mass planets orbiting M dwarfs, with a rate of at least one planet per star (e.g. Bonfils et al. 2013; Dressing & Charbonneau 2013, 2015; Tuomi et al. 2014). This fact is also supported by the number of dedicated surveys targeting M dwarfs, notably in the search of rocky planets in their habitable zones (HZs). Some examples are the SPECULOOS (Gillon 2018; Burdanov et al. 2018; Delrez et al. 2018; Jehin et al. 2018), CARMENES (Quirrenbach et al. 2010; Reiners et al. 2018), and MEarth (Nutzman & Charbonneau 2008) projects.

More interestingly, the great number of M-type stars has significantly increased the chances of finding Earth-mass rocky planets in temperate orbits, that is, within the host star’s HZ. Nevertheless, planets in the HZs of M dwarfs are close to the star and, thus, they are generally exposed to intense radiation and magnetic fulgurations. Whether a planet can remain habitable in the presenceof UV and X-ray flares is still a matter of debate: while a single flare probably would not compromise a planet’s habitability, planets orbiting M dwarfs may suffer frequent and strong flares at younger ages (Armstrong et al. 2016).

Among all the planets discovered so far orbiting M dwarfs, some of the most interesting are also the ones closest to Earth. One such example is the temperate Earth-mass planet orbiting Proxima Centauri-b, which is located only 1.2 pc away (Anglada-Escudé et al. 2016), while the seven Earth-size planets orbiting the TRAPPIST-1 system at a distance of 12.1 pc, where at least three of them are in the HZ (Gillon et al. 2016, 2017), are equally alluring. In this context, GJ 273 is a highly interesting M dwarf (also known as Luyten’s star) hosting a multi-planet system, which is one of the closest stars to the Sun located at only 3.75 pc (Gatewood 2008).

Compared with other M-type stars, GJ 273 is a bright star (9.87 Vmag, Koen et al. 2010), with a main-sequence spectral type of M3.5V (Hawley et al. 1996), which shows a small projected rotational velocity (v sin i < 2.5 km s−1, Browning et al. 2010). In the work by Astudillo-Defru et al. (2017b), the authors made use of a data set of 280 RVs measured with HARPS (mounted on ESO’s 3.5-m La Silla telescope in Chile) which covered 12 yr of intense monitoring of GJ 273. Astudillo-Defru et al. (2017b) reported the presence of two short-period planets, GJ 273b and GJ 273c, with minimum masses of 2.89 and 1.18 _M_⊕ and orbital periods of 18.6 and 4.7 d, respectively. Planet GJ 273b was shown to be located in the HZ of GJ 273 which, along with its close proximity to Earth, sparked interest in this particular system. In addition, in the recent study presented by Tuomi et al. (2019), the authors made use of a data set that combined HARPS, HIRES (Keck), PFS (Subaru), and APF (Lick), for a total of 466 RVs measurements, and fully confirmed the existence of the two planets previously detected by Astudillo-Defru et al. (2017b) in addition to suggesting the presence of two extra planets at larger heliocentric distances: GJ 273d, with a minimum mass of 10.8 _M_⊕ and an orbital period of 413.9 d, and GJ 273e with a minimum mass of 9.3 _M_⊕ and an orbital period of 542 d1. The authorsclaimed these new signals as “planetary candidates” which need more data to verify their planet nature. In addition, Astudillo-Defru et al. (2017b) reported the presence of long-period signals with P ~ 420 and 700 d. However, both signals were not considered by the authors as planetary candidates but were, instead, attributed to stellar activity. Prior to that report, Bonfils et al. (2013) also found a controversial signal of a possible planet at P ~ 420 d, however, due to their poor phase coverage, the authors considered it an unclear detection.

Collectively, the present data have prevented us from making any definite statements regarding the precise number of planets in the GJ 273 system. Therefore, since the existence of the two outermost planets, GJ 273d and GJ 273e, is still under debate, in the present study we consider two options: two– and four–planet configurations. In Table 1, we present the known physical characteristics of the planets and the star.

The paper is organised as follows: in Sect. 2, we attempt to reproduce the system architecture using accretion simulations where we track the content of water of the planets. In Sect. 3, we review the current observational data set that consists of RVs and photometric measurements. In Sect. 4, we analyse the stability of the system and provide extra constraints on the planetary parameters and the general system architecture, while in Sect. 5 we explore the similarities of GJ 273 with the Solar System while searching for minor-body reservoir analogues. In Sect. 6, we discuss the habitability of planet GJ 273b and, finally, in Sect. 7, we present our conclusions.

2 Formation and evolution of the planetary system GJ 273: _N_-body experiments

In the present section, we analyse the processes of formation and evolution among the planetary system GJ 273 using _N_-body simulations. To do this, we first describe the properties of the protoplanetary disk and then determine the physical and orbital initial conditions that are to be used in our model. Then we define the different scenarios to be modelled and, finally, we discuss the results obtained from our _N_-body experiments.

2.1 Protoplanetary disk: properties

An important parameter that determines the distribution of the material in a protoplanetary disk is the surface density. The gas-surface density profile Σg(R) adopted in the present research is given by:Σg(R) = Σ0g(RRc)−γe−(R/Rc)2−γ,(1)

where R is the radial coordinate in the disk mid-plane, Σ0g a normalisation constant, _R_c the characteristic radius, and γ an exponent that determines the density gradient. In the same way, the solid-surface density profile Σs (R) is written as:Σs(R) = Σ0sηice(RRc)−γe−(R/Rc)2−γ,(2)

where Σ0s is a normalisation constant and parameter _η_ice represents the change in the amount of solid material due to the condensation of water beyond the snow line, which is assumed to be initially located at _R_ice.

To determine the normalisation constant Σ0g, we integrated the expression concerning the mass of the protoplanetary disk assuming axial symmetry. From this, we can express the total mass of the disk _M_d as:Md = ∫0∞2 πRΣg(R) dR,(3)

from which we obtainΣ0g = (2−γ)Md2πRc2.(4)

The normalisation constant Σ0s can be obtained making use of the relationΣ0s = z0Σ0g10−[Fe/H],(5)

where [Fe/H] is at stellar metallicity and _z_0 is the primordial abundance of heavy elements in the sun, which has a value of _z_0 = 0.0153 (Lodders et al. 2009).

To quantify the gas- and solid-surface density profiles involved in our model, we must assign values to several parameters such as the stellar metallicity [Fe/H], the mass of the disk _M_d, the characteristic radius _R_c, the exponent γ, the location of the snow line _R_ice and the _η_ice parameter. On the one hand, the GJ 273 star has a mass of _M_⋆ = 0.29 _M_⊙ and a metallicity of [Fe/H] = 0.09 (Astudillo-Defru et al. 2017b). On the other hand, our model uses a disk mass of _M_d = 5% _M_⋆, a characteristic radius of _R_c = 20 au, and an exponent of γ = 1, all of which is consistent with observational studies of protoplanetary disks carried out as part of several studies (e.g. Andrews et al. 2010; Testi et al. 2016; Cieza et al. 2019) and with population synthesis models of protostellar discs (Bate 2018). The position of the snow line _R_ice must be carefully chosen. To do this, we followed the procedure proposed by Ciesla et al. (2015), who studied the delivery of volatiles to planets from water-rich planetesimals around low-mass stars. That work assumed that _R_ice’s location evolves in a disk as a result of the variations over time of the stellar radius and effective temperature. Based on such an assumption, the snow line is about at 0.65 au (disk temperature of 140 K, age of 1 Myr) when the stellar parameters given by Siess et al. (2000) are considered. Although the location of the sow line evolves over time, we followed Ciesla et al. (2015) by assuming a fixed value for _R_ice for the whole duration of our _N_-body experiments (see Apendix A).

It is also important to consider that the location of the snow line generally represents the transition from dry to water-rich bodies in the disk. However, any increase in the amount of solid material due to the condensation of water beyond the snow line is strongly dependent on the existence of short-lived radionuclides, notably Al26. According to the work carried out by Lodders (2003) on Solar System abundances, we consider that the _η_ice parameter has a value of 1 and 2 inside and beyond _R_ice, respectively. From this, we assume that bodies initially located inside the snow line are primarily made out of rocky material with only 0.01% water by mass, while the initial content of water by mass of those bodies associated with the outer region of the system beyond the snow line is of the order of 50%.

Our investigation is focused on the formation, evolution, and survival of the confirmed planets in GJ 273, especially the one close to the HZ. Thus, it is necessary to specify the limits of such a region in the system under study. We derive conservative and optimistic estimates for the width of the HZ around GJ 273 from the models developed by Kopparapu et al. (2013a) of 0.1–0.19 au and 0.08–0.2 au, respectively. We assume that a planet is in the HZ of the system if its whole orbit is contained inside the optimistic edges. According to this, a planet is considered to be in the HZ of the system if its pericentric distance was q ≥ 0.08 au and its apocentric distance was Q ≤ 0.2 au.

Based on the data considered in this section, we assume that the extension of the region of study is between 0.01 au and 1.5 au since it includes the HZ, the snow line, and an outer water-rich embryo population.

2.2 _N_-body experiments: characterisation, parameters, and initial conditions

In this sub-section, we analyse the process of formation and evolution of terrestrial-like planets in the GJ 273 system from _N_-body experiments. The MERCURY _N_-body code was used to develop our simulations (Chambers 1999). In particular, we made use of a hybrid integrator, which uses a second-order mixed variable symplectic algorithm to treat the interactions between objects with separations greater than 3 Hill radii, and the Bulirsch–Stoer method for resolving closer encounters.

In general terms, the accretion models were aimed at studying the formation of terrestrial-like planets via _N_-body experiments to track the evolution of planetary embryos and planetesimals in a given system. However, _N_-body simulations that include a population of embryos and thousands of planetesimals are very costly numerically. In the present research, we consider a simplified model, which assumes a disk only composed of planetary embryos. The MERCURY code computes the temporal evolution of the orbits of the embryos by simulating the gravitational interactions among them. Such interactions produce orbits that cross each other, which results in successively close encounters between them which, in turn, lead to collisions among them, ejections from the system, and collisions with the central star. Several works have shown that collisions between gravity-domained bodies yield different outcomes depending on the target size, projectile size, impact velocity, and impact angle (Leinhardt & Stewart 2012; Chambers 2013; Quintana et al. 2016; Dugaro et al. 2019). In the present study, all collisions were treated as perfect mergers, conserving the total mass of the interacting bodies in each impact event. This approach has been the main hypothesis of multiple works concerning accretion simulations aimed at studying the process of planetary formation in systems with different dynamical scenarios (e.g. Raymond et al. 2004, 2009; O’Brien et al. 2006; Lissauer 2007; De Elía et al. 2013; Dugaro et al. 2016; Darriba et al. 2017; Zain et al. 2018).

Our _N_-body simulations assumed 53 planetary embryos between 0.01 au and 1.5 au, 28 (25) of which were initially located before (beyond) the snow line at 0.65 au, and possessed individual masses and physical densities of _m_p = 0.102(0.28) _M_⊕ and _ρ_p = 3 (1.5) g cm−3, respectively. The embryos were distributed across the region of study using the acceptance–rejection method developed by Jonh von Neumann from the functiondN(R) = 2πmpRΣs(R)dR,(6)

where Σs(R) is represented by Eq. (2). Given that the planetary embryos were assumed to start the simulation in quasi-circular and coplanar orbits, the values of the variable, R, generated from Eq. (6) were adopted as their initial semi-major axes a. In fact, the initial eccentricities, e, and inclinations, i, of the embryos were calculated randomly from uniform distributions assuming values lower than 0.02 and 0.5°, respectively. In the same way, the argument of the pericenter, ω, the longitude of the ascending node, Ω, and the mean anomaly, M, of the planetary embryos were chosen randomly from uniform distributions between 0° and 360°.

The _N_-body experiments developed using the hybrid integrator of the MERCURY code required a good selection of the time step to carry out a correct integration of the system studied. Based on this, we used a time step of 0.03 d, which is shorter than 1/20th of the orbital period of the innermost body in our scenario of work. Moreover, in order to avoid the integration of orbits with very small pericenters, our model assumed a non-realistic value for the radius of the central star of 0.005 au. Finally, we considered an ejection distance of 1000 au in all our _N_-body experiments.

In order to analyse the formation and evolution history of the planets that compose system GJ 273, we carried out two sets of _N_-body accretion simulations:

• Set 1, which analyses the post-gas phase;
• Set 2, which includes the effects of the gas disk.

“Set 2” _N_-bodysimulations included modelling the effects of the gaseous component of the disk in the early evolution of the planetary embryos. To do this, we developed an improved version of the MERCURY code, which included analytical prescriptions for the orbital migration rates and eccentricity and inclination damping rates derived by Cresswell & Nelson (2008) from the formulae for damping rates obtained by Tanaka & Ward (2004) and the migration rates obtained by Tanaka et al. (2002). We give a detailed description of this in Sect. 2.4.

Due to the stochastic nature of the accretion process, we carried out between 10 and 15 runs for each set of numerical simulations. We must remark that, in both of them, each _N_-body experiment was integrated for 100 Myr.

2.3 “Set 1”: collisional accretion of protoplanets after gas dissipation

In this set of _N_-body experiments, we considered that the initial physical and orbital parameters assigned to the planetary embryos in Sect. 2.2 were associated with the post-gas phase. Thus, here, our main goal is to study the formation and evolution of terrestrial-like planets via the collisional accretion of embryos once the gas has dissipated from the system.

The left panel of Fig. 1 shows the mass distribution as a function of the semi-major axis of the planets formed in 9 _N_-body simulations of “Set 1” after 100 Myr of integration. Moreover, the light blue shaded area illustrates the optimistic limits of the HZ of the system. The distribution of planets offers two important results. On the one hand, the most massive planets in this scenario formed at locations close to the snow line, which was initially at around 0.65 au. In fact, our _N_-body experiments produced massive planets of about 5.4 _M_⊕ and semi-major axes ranging between 0.50 au and 0.55 au. On the other hand, our _N_-body simulations did not efficiently produce planets in the HZ nor in the inner region of the system. In fact, only three planets formed in the HZ with masses between 1.6 _M_⊕ and 2.1 _M_⊕, while no planet survived in the inner region of the system at 100 Myr. According to this, this scenario was not able to produce planets analogous to those observed in the real system around GJ 273.

2.4 “Set 2”: effect of the gas disk

In this set of _N_-body experiments, we tested the effects of the presence of a gas disk in the early evolution of the planetary embryos. To do this, we included additional forces to model the effects of a dissipating gaseous disk on the orbital parameters of the planetary embryos. In particular, our improved version of the MERCURY code included effects that led to changes of the semi-major axis (hereafter, “type-1 migration”) as well as decays in eccentricity and inclination (hereafter, “type-1 damping”) of the embryos.

The type-1 migration and type-1 damping were modelled based on hydrodynamic simulations focused on planetary embryos embedded within an idealised isothermal disk (Tanaka et al. 2002; Tanaka & Ward 2004). However, those works are only valid for embryos that are assumed to have no disk gap and with low values associated with their eccentricities and inclinations. In order to develop a more general treatment, we modified the model to include the additional terms derived by Cresswell & Nelson (2008), which allowed us to describe the type-1 migration and type-1 damping for embryos evolving on orbits with large eccentricities and inclinations.

In idealised vertically isothermal disks, the type-1 migration derived by Tanaka et al. (2002) can lead to the very fast inward orbital migration of the planetary embryos. In general terms, the migration times calculated with this model were comparable to or smaller than the known lifetimes of gas disks, which may be at odds with current planetary formation models. Type-1 migration tends to be more complex when a more detailed treatment concerning the physics of the disk is considered (e.g. Paardekooper et al. 2010, 2011; Guilet et al. 2013; Benítez-Llambay et al. 2015). Beyond this, several authors have carried out population synthesis studies using the expressions of Tanaka et al. (2002) concerning type-1 migration rates, with an added reduction factor in order to reproduce observations (e.g. Alibert et al. 2005; Mordasini et al. 2009; Miguel et al. 2011; Ronco et al. 2017).

Here, we considered the type-1 migration rate for isothermal disks of Tanaka et al. (2002) and Cresswell & Nelson (2008), and we also incorporated a factor, _f_1, that ranged between 0 and 1. According to this, the type-1 migration time is inversely proportional to such a factor, for which the larger the value of _f_1, the faster the orbital migration. From this, we carried out _N_-body experiments aimed at studying the gas effects on the formation and evolution of the planetary system under consideration for values of _f_1 between 0.01 and 1. Naturally, our results showed that the larger the value of _f_1, the smaller the semi-major axes of the resulting planets.

The gas disk was modelled by adopting the gas-surface density profile Σg (R) given by Eq. (1) and a heightscale of H = 0.0361(R/1ua)5/4 au. In this scenario, we considered that the gaseous component dissipates exponentially over 2.5 Myr, which is consistent with studies concerning the lifetimes of primordial disks associated with young stellar clusters developed by Mamajek et al. (2009). According to this, it is very important to remark that the zero-time of our _N_-body experiments corresponds to 1 Myr for which the gas effects were modelled for 1.5 Myr in all the numerical simulations.

The right panel of Fig. 1 illustrates the mass of the planets formed in 14 _N_-body simulations of “Set 2” after 100 Myr of integration as a function of the semi-major axis, assuming a reduction factor, _f_1, for a type I migration of 0.04. According to this, the most massive planets produced in this scenario were located in regions of the system that were more internal than those formed in “Set 1”, which of course did not include any gaseous effects. In fact, the most massive planets resulting from “Set 2” were produced around the HZ of the system, with masses between 2.2 _M_⊕ and 3.22 _M_⊕. These planets showed physical and orbital parameters similar to those associated with GJ 273b, which is a super-Earth with a minimum mass of 2.89 ± 0.26 _M_⊕ orbiting around its parent star with a semi-major axis of 0.09110 ± 0.00002 au, and an eccentricity of 0.10−0.07+0.09 (see Table 1).

Moreover, the _N_-body experiments of “Set 2” efficiently produced a close-in planet population within the inner edge of the HZ with masses ranging between 0.1 _M_⊕ and 1.3 _M_⊕. The physical and orbital properties of the most massive planets of this close-in population are consistent with those of GJ 273c, which has a minimum mass of 1.18 ± 0.16 _M_⊕, a semi-major axis of 0.036467 ± 0.000002 au, and an orbital eccentricity of 0.17−0.12+0.13 (see Table 1).

According to “Set 2”, which included gas dynamics, simulated planets with properties that were not dissimilar to the two known planets in the GJ 273 system were formed. It is worth remarking that gas effects play a key role in the final planetary architecture of simulated systems. In fact, the location of the peak of the planet mass distribution at more internal regions of the system in comparison with that associated with “Set 1”, and the formation of a close-in planet population within the inner edge of the HZ, are a clear consequence of incorporating type-1 migration and type-1 damping in the _N_-body experiments of “Set 2”.

The right panel of Fig. 1 also allows us to derive important conclusions concerning the final water content of the resulting planets. On the one hand, all planets with masses greater than 1 _M_⊕ that survived in the HZ of the system became true water worlds. In fact, in general terms, the accretion seeds2 started the simulation beyond the snow line with 50% of water by mass. Throughout the 100 Myr of evolution, they were impacted by planetary embryos associated with different regions of the disk, reaching final masses between 1.16 _M_⊕ and 3.22 _M_⊕, and final contents of 18–50% of water by mass. On the other hand, all planets with masses less than 1 _M_⊕ that survived in the HZ had accretion seeds that started the simulation inside the snow line with 0.01% water by mass. These planets had final masses between 0.3 _M_⊕ and 1.02 _M_⊕ and kept their very low initial fraction of water throughout their entire evolution since they only accreted planetary embryos from the inner region to the snow line. The right panel of Fig. 1 shows that GJ 273b is a potentially habitable planet, which is located in a region of the mass–semi-major axis plane associated with the presence of water worlds. According to this model, GJ 273b should have a very high water content by mass.

In this line of research, we suggest that planet GJ 273c should have physical properties different to those shown by GJ 273b concerning their final water contents. In fact, the right panel of Fig. 1 shows that GJ 273c is located in a region of the system where our _N_-body experiments produced planets with very low final water content.

Figure 2 shows four snapshots in time of the semi-major axis–eccentricity plane of the evolution of a given _N_-body simulationof “Set 2”, which exposes the best scenario that we found. The solid circles represent the planetary embryos (and planets) and the colour code refers to their fraction of water by mass. The solid blue lines illustrate curves of constant pericentric distance (q = 0.08 au) and apocentric distance (Q = 0.2 au), while the light blue shaded area represents the optimistic HZ of the system. Our study shows two relevant results. On the one hand, the “Set 2” simulations were able to produce planets with physical and orbital properties consistent with those associated with the planets around GJ 273 (Astudillo-Defru et al. 2017b). In fact, this _N_-body experiment formed a planet analogous to GJ 273b(c) with a mass of 2.25 (1.13) _M_⊕, a semi-major axis of 0.1095 (0.0482) au, and an orbital eccentricity that oscillates between 0.0007 (0.0081) and 0.1 (0.176) during the last 50 Myr of evolution. On the other hand, the right and bottom panel of Fig. 2 shows that additional planets have formed in different regions of the system after100 Myr of evolution.This result suggests that more planets with different physical and orbital properties might be orbiting around GJ 273.

It is important to remark that the final water contents of the planets formed in our simulations should be interpreted as upper limits. On the one hand, our _N_-body experiments treat collisions as perfect mergers, which conserve the mass and the water content of the interacting bodies. Thus, we did not account for water loss during impacts, which could be relevant for small impact angles and high velocities Dvorak et al. (2015). On the other hand, our model did not include processes such as photolysis and hydrodynamic escape, which could reduce the fraction of water of the planets (Luger & Barnes 2015; Johnstone 2020). From this last consideration, GJ 273c could be certainly a dry planet since it is interior to the HZ, while GJ 273b could contain less water than that obtained in our simulations since the initial boundaries of the HZ are located in more external regions of the disk, and then, they evolve inward significantly in the first 100 Myr of evolution for a 0.3 _M_⊙ M-dwarf (Luger & Barnes 2015). A detailed analysis about how the water loss processes mentioned in this paragraph affect the planets of our simulations is beyond the scope of this work.

Fig. 1Mass distribution of the planets formed in the _N_-body simulations after 100 Myr of Set 1 (left panel), which only model the post-gas stage, and Set 2 (right panel), which include the effects of the gas disk. In every panel, the black circles illustrate the two planets confirmed that orbit GJ 273 (Astudillo-Defru et al. 2017b), while the color code represents the final fraction of water by mass of the planets resulting from our numerical experiments. Finally, the light blue shaded region refers to the optimistic HZ derived from the model developed by Kopparapu et al. (2013a,b).

Fig. 2Temporal evolution of semi-major axis as a function of the eccentricity in the “Set 2” _N_-body simulations, which illustrates the most likely planetary-formation scenario. The planetary embryos (and the resulting planets) are represented by circles, which show a colour code that indicates their fraction of water by mass. The solid blue lines illustrate curves of constant pericentric distance (q = 0.08 au) and apocentric distance (Q = 0.2 au), while the light blue shaded region represents the optimistic HZ of the system.

3 Observational constraints: radial velocities and photometric measurements

3.1 HARPS observations

The simulations performed in Sect. 2 work well in predicting the formation of two planets at the orbits of GJ 273b and GJ 273c. However, our results for the best model indicated that more planets might be orbiting GJ 273 (see Fig. 2, panel corresponding to 100 Myr). The aim of this subsection is to demonstrate that even if these planets exist, they would not be detected given the precision and cadence of the published data used to confirm GJ 273b and GJ 273c. The fact that Tuomi et al. (2019) is still undergoing the referee process at the time of writing prevented us from using their data set in this experiment. Therefore, to test the hypothesis of the existence of extra planets in the system, we obtained the HARPS RVs provided by Astudillo-Defru et al. (2017b) in their Table A.1. First, we subtracted the secular acceleration following Eq. (2) from Zechmeister et al. (2009). Then we obtained the residuals to a model containing all signals present in the HARPS data. This included, in addition to the 18.6 and 4.7 d signals of GJ 273b and GJ 273c, two large periodicities at ~420 and ~ 700 d that the authors interpreted as stellar activity cycles. After having obtained the residuals, we injected the RV time series representing the planets predicted by the aforementioned dynamical simulations. The goal was to check whether such signals could be detected. In particular, we defined sinusoidal models at the time epoch of HARPS for the planet candidates represented in the last panel of Fig. 2. The periods of the sinusoids were calculated from the semi-major axes of the planets predicted by the dynamic simulations using Kepler’s third law and assuming circular orbits. Additionally, we included a white noise component calculated by randomly swapping the RV data among different observing epochs and adding the resulting RVs to the sinusoidal model (via bootstrapping method). This process was repeated 1000 times for each period (P i) and amplitude (K i) thus obtaining 1000 simulated signals with random phases. Therefore, our goal was to progressively increase the amplitude K i and calculate the rate of simulated signals that could be detected. Such a rate is related to the minimum mass that the hypothetical planet must have in order to be detected with the cadence and precision of the HARPS observations, assuming circular orbits. To this aim, we calculated the completeness (C) as the fraction of positive experiments (_N_pos) among a total of 1000experiments (N) for which we recovered a signal equivalent to the sinusoid injected:C = 100×(Npos/N).(7)

To consider a recovered signal equivalent to the injected signal, meaning a positive detection, we established that three criteria had to be satisfied: (1) the period recovered (P_r) had to match the injected one within the range given by the resolution frequency (P_r = [1∕(P i + 1∕_T), 1∕(P i − 1∕_T)]; where T refers to thetotal time span of the data set), (2) the highest peak of the differential likelihood periodogram had to be above the 1%-False-Alarm-Probability3, and (3) the difference between the initial and recovered amplitudes had to be below the RMS of the initial residuals (K _r_− K i < RMSres).

Once we obtained the completeness, we considered that signals with C ≥ 90% could be confidently recovered. This approach allowed us to obtain the lowest amplitude at which an RV signal could be detected, _K_min. Figure 3 shows how the completeness increases with the initial amplitude for one of the cases that we tested, that is, the set of simulated sinusoids with P = 652.3 d. Finally, we obtained the lower limit that the minimum mass of a planet had to reach to be detected as:Msin i = Kmina M⋆G,(8)

where G is the gravitational constant, _M_⋆ is the mass of the star (0.29 _M_⊙), and a is the semi-major axis obtained by making use of the Kepler’s third law.

The results are summarised in Fig. 4. There we can see how the masses of the hypothetical planets (blue crosses) are below the detection region (grey area). On the contrary, the GJ 273b and GJ 273c analogues predicted by our dynamical simulations (red crosses) are on or above the detection limit (black line) and, as consequence, they could be detected using the RV data set. The large peaks and white gaps at ~ 0.013 au and at large semi-major axis indicate misleading regions where the temporal aliases prevent a good monitoring. This is a consequence of the spectral window shown in the upper part of the figure.

	Fig. 3Percentage of well-detected sinusoids from a total of N = 1000 experiments (i.e. completeness, C) versus the initial amplitude. The case of simulated sinusoids with P = 652.3 d is illustrated here as an example of how the C is calculate. The black line indicates the polynomial model fitted to get the initial amplitude (K i) that corresponds to a 90% completeness (red cross).

3.2 TESS observations

TESS (the Transiting Exoplanet Survey Satellite; Ricker et al. 2014) is currently performing a two-year survey of nearly the entire sky, with the main goal of detecting exoplanets smaller than Neptune around bright and nearby stars. It has four 24° × 24° field-of-view cameras, each containing four 2 k × 2 k CCDs, with a pixel scale of 21′′. At the time of writing, TESS has discovered several planetary systems such as L 98-56 (Kostov et al. 2019), TOI-125 (Quinn et al. 2019), TOI-270 (Günther et al. 2019), LP 791-18 (Crossfield et al. 2019), and LP 729-54 (Nowak et al. 2020) among others. GJ 273 was observed by TESS in sector 7 with CCD camera 1 from 2019-Jan.-07 to 2019-Feb.-02 (TIC 318686860). An alert was not issued by TESS for this object, due to the lack of clearly identified transits by both the QLP (Quick Look Pipeline) and the SPOC (Science Processing Operations Center). With this in mind, we performed our own search for threshold-crossing events that might not have been detected by the aforementioned automated pipelines.

Our custom-made pipeline SHERLOCK4 (Searching for Hints of Exoplanets fRom Lightcurves Of spaCe-based seeKers) made use of the LIGHTKURVE package (Lightkurve Collaboration 2018) to obtain the TESS data of GJ 273 from the NASA Mikulski Archive for Space Telescope (MAST). We used thePre-search Data Conditioning Simple APerture (PDC-SAP) flux given by the SPOC. First, we removed the outliers, defined as data points >3_σ_ above the running mean. Then, we made use of WOTAN (Hippke et al. 2019) to detrend the flux using the bi-weight method. This choice was based on its good speed–performance ratio compared with other methods such as Gaussian processes. The bi-weight method is a time-windowed slider, where shorter windows can efficiently remove stellar variability, but there is an associate risk of removing any actual transit signal. In order to prevent this issue, our pipeline explored twelve cases with different window sizes (see Fig. B.1). The minimum value was obtained by computing the transit duration (_T_14) of a similar planet to the outermost one confirmed in the system, that is, a hypothetical Earth-size planet with a period of 19 d orbiting the given star. To protect a transit of this duration, we chose a minimum window size of 3 × _T_14. The maximum value explored was fixed to 20 × _T_14, which seemed enough to remove the variability of GJ 273. For all the detrended light curves, the transit search was performed using the TRANSIT LEAST SQUARE (TLS) package, which uses an analytical transit model based on stellar parameters, and is optimised for the detection of shallow periodic transits (Hippke & Heller 2019). We searched for signals with periods ranging from 0.5 to 25 d, with signal-to-noise ratios (S_∕_N) of ≥ 5 and signal detection efficiencies (SDEs) of ≥5, following the vetting process described by Heller et al. (2019). Bearing in mind that planet GJ 273b has an orbital period of 18.6 d, since the TESS data span 24.45 d, including a gap of 1.67 d (TJD = 1503.04–1504.71) caused by the data download during telescope apogee, this means that only a single transit may have been observed for this planet. The same situation may also apply for the outermost planets GJ 273d and GJ 273e whose orbital periods exceed 400 d, which means that in the best-case scenario only a single transit might be detected for each planet in the current set of data. We did not find any signal which fulfilled the criteria of the vetting process, which agrees with the negative results yield so far by the automatic pipelines of TESS. This is a hint that GJ 273 is not a transiting system.

However, in order to quantify the detectability of the planets in the existing set of data, we carried out injection-and-recovery tests. These consisted of injecting synthetic planetary signals into the PDC-SAP flux light curve before the detrending and the TLS search. Since there are two confirmed planets and two candidates, with significant differences in their physical properties (i.e. the innermost planets are likely Earth- and super Earth-mass planets with orbital periods < 20 d, and the outermost planets are likely mini-Neptunes with orbital periods > 400 d), we performed three independent injection-and-recovery tests: (1) the first one focused on the innermost planets GJ 273b and GJ 273c, where we studied the _R_planet–_P_planet parameter space in the ranges of 0.8–3.0 _R_⊕ with steps of 0.05 _R_⊕ and 1–25 d with steps of 1 d, that is, we explored 1080 different scenarios; (2) second, we explored the parameter space around the planet GJ 273d, over the radius range 1.0–4.0 _R_⊕, with steps of 0.05 _R_⊕ and periods of 400–425 d with steps of 1 d; therefore, we explored a total of 1625 scenarios. (3) For the planet GJ 273e, we studied radii from 1.0–4.0 _R_⊕ with steps of 0.05 _R_⊕ and periods of 530–555 d with steps of 1 d, which yielded a total of 1625 scenarios. The range of the planetary radii explored here was chosen from the dynamical arguments presented in Sect. 45. For simplicity, we assumed the impact parameters and eccentricities were zero. In these inject-and-recovery tests, we considered that an injected signal was recovered when a detected epoch matched the injected epoch within one hour and if a detected period matched any half-multiple of the injected period to better than 5%. As we injected the synthetic signals in the PDC-SAP light curve, the signals were not affected by the PDC-MAP systematic corrections; therefore, the detection limits should be considered as the most optimistic scenario (Eisner et al. 2020). For test (1), we found that the innermost planet GJ 273c had a recovery rate of 100%, and the planet GJ 273b had a recovery rate of about 40%, where these results are shown in Fig. 5. On the other hand, for tests (2) and (3), the planet signals were recovered at very low rates < 2%. This means that a transit of planet GJ 273c can be ruled out by the TESS data, while transits for the other planets in the system might have been missed due to low-S/N data or because of the limited time-coverage of the data.

Fig. 4Minimum (lower-limit) masses of the planetary signals depending on their semi-major axes (black dots). Circular orbits were assumed. Only planets within the grey areas could be detected using the HARPS data set that confirmed GJ 273b and GJ 273c. Blue crosses, which correspond to the masses of the hypothetical planets predicted by the dynamical models (see last panel of Fig. 2) are belowthis detection threshold. The inset shows a close-up of the GJ 273b and GJ 273c analogues predicted by the dynamical models. The minimum GJ 273b and GJ 273c masses reported by Astudillo-Defru et al. (2017b) are highlighted with red crosses and appear on or above the detection threshold. Open circles indicate experiments for which 90% completeness was not reached even when running the experiments with K i up to 6.75 ms−1. They correspond to high power regions of the spectral window shown in the top panel.

4 Dynamical evolution

4.1 Stellar age

Stellar age is a key value used to evaluate the dynamical stability of the whole planetary system. In our case, GJ 273 is an M3.5V star, whose mass has been estimated to be 0.29 _M_⊙ thanks to the empirical mass–luminosity relation presented by Delfosse et al. (2000). The evolution of such a low-mass star is very slow and, hence, stellar evolutionary models are rather useless with regard to establishing its accurate age. According to the PARSEC evolutionary tracks and isochrones (Bressan et al. 2012; Chen et al. 2014), the stellar effective temperature, _T_eff, and luminosity, _L_⋆, of a 0.3 _M_⊙ main sequence star vary by ~0.1% and ~10%, respectively,over the entire age of the Universe. On the other hand, it is well known that stellar age correlates with the stellar rotation rate and stellar activity, as demonstrated by several empirical relations provided in the literature. In fact, considering that GJ 273 has a rotational period of _P_⋆ = 99 d (Astudillo-Defru et al. 2017b), the gyrochronological relation presented by Barnes (2010) yields an age of _t_gyro = 8.4 Gyr. If, instead, we consider logR′HK = −5.56 (Astudillo-Defru et al. 2017b) as the index of chromospheric activity, the relation found by Mamajek & Hillenbrand (2008) yields an age _t_HK = 8.7 Gyr.

It is important to stress that both the empirical age-activity and gyrochronological relations are calibrated using stars (usually belonging toclusters or stellar associations) whose ages are already known thanks to isochrone fitting (i.e. evolutionary models) and whose ages are generally young. In fact, as stated by Soderblom (2010), for old stars, the portion of the chromospheric signal that can be attributed to age is uncertain because it is difficult to detect and calibrate the continuing decline of activity in stars older than ~ 2 Gyr. In addition, there is a lack of open clusters for calibrating the gyrochronological ages of old stars. Meibom et al. (2015) confirmed a well-defined period–age relation up to 2.5 Gyr, which can be considered as the upper limit at which gyrochronological ages are conservatively reliable. Keeping these conservative constraints in mind, the low stellar activity level and rotational rate, and the consistency among _t_gyro and _t_HK, all suggest an age higher than 8 Gyr.

Fig. 5Injection-and-recovery test performed to check the detectability of the innermost planets GJ 273b and GJ 273c in the current set of TESS data. We injected synthetic transits into the available photometric data corresponding to each planet by taking into account the uncertainties in their physical properties. We explored a total of1080 different scenarios. We found that planet GJ 273c showed a recovery rate of 100%, which clearly suggests that this planet is not transiting. On the other hand, planet GJ 273b showed a recovery rate of ~ 40%, implying its possible detection in future observations. We carried out two other experiments to check the detectability of the two outermost planets GJ 273d and GJ 273e, finding recovery rates < 2%; see the main text for more details.

4.2 Stability of the system

In order to explore the stability of the system, we used the mean exponential growth factor of nearby orbits, Y (t) (MEGNO, Cincotta & Simó 1999, 2000; Cincotta et al. 2003). MEGNO is a chaos index which evaluates the stability of a body’s trajectory after a small perturbation and predefined initial conditions. MEGNO is closely related to the maximum Lyapunov exponent and the Lyapunov characteristic number, which are well-known, powerful tools for studying the long-term evolution of Hamiltonian systems (Benettin et al. 1980; Goździewski et al. 2001). Therefore, MEGNO provides a clear picture of resonant structures and establishes the locations of stable and unstable orbits. To calculate MEGNO, Y (t), each body’s six-dimensional displacement vector, δ i, is considered as a dynamical variable in both position and velocity from its shadow particle, that is, a particle with slightly perturbed initial conditions. Then, for each δ i differential equation, the variational principle is applied to the trajectories of the original bodies. Then, MEGNO is straightforwardly computed from the variations as:Y(t) = 2t∫t0t‖δ˙(s)‖‖δ(s)‖ s ds,(9)

along with its time-average mean value:〈Y(t)〉 = 1t∫t0tY (s) ds(10)

Thanks to the time-weighted factor, the amplified stochastic behavior allows for the detection of hyperbolic regions in the time interval (t_0, t). To distinguish between chaotic and quasi-periodic trajectories in phase space, we evaluated ⟨_Y (t)⟩. Here, if ⟨Y (t)⟩→∞ for t → ∞ the system is chaotic. On the other hand, if ⟨Y (t)⟩→ 2 for t →∞ then the motion is quasi-periodic. These criteria have been extensively used within dynamical astronomy, for both the Solar System and extrasolar planetary systems (e.g. Jenkins et al. 2009; Hinse et al. 2015; Wood et al. 2017; Horner et al. 2019). We applied the MEGNO implementation within the N_-body integrator REBOUND (Rein & Liu 2012), which made use of the Wisdom-Holman WHfast code (Rein & Tamayo 2015). This allowed us to explore the whole parameter space at a small computational cost. As a first approach we only considered the two planets already confirmed in the system, that is, GJ 273b and GJ 273c. We then explored the stability of the system for different planetary eccentricities and inclinations following Jenkins et al. (2019). The natural trend of planetary systems is to reside in a nearly co-planar configuration, hence, we assumed planets b and c as being co-planar, and the initial longitude of the ascending node Ω for both planets was taken as zero. This choice was motivated by the results yielded by Kepler and the HARPS survey, which suggest that the mutual inclinations of multi-planets systems are of the order ≲10° (see e.g. Fabrycky et al. 2014; Tremaine & Dong 2012; Figueira et al. 2012). We then conducted dynamical simulations to explore the minimum values of the inclinations for planets GJ 273b and GJ 273c. Furthermore, their masses are dependent on their orbital inclination angles, the so called mass–inclination degeneracy of radial velocity data, which means that any uncertainties in their inclination angles will seriously affect their determined masses. Hence, exploring the minimum values of their inclination angles is equivalent to exploring the maximum values of their masses. We evaluated their inclinations 1000 times, ranging from 10° and 90° for three different set of eccentricities corresponding with the nominal values and ± 1_σ uncertainties given in Table 1: (1) _e_b = 0.03 and _e_c = 0.05; (2) _e_b = 0.10 and _e_c = 0.17; and (3) _e_b = 0.19 and _e_c = 0.30. Our choice to limit the inclination to 10° was motivated by geometrical arguments, which for randomly oriented systems, disfavour small values of i (see e.g. Dreizler et al. 2020). Furthermore, orbital inclinations close to zero (face on) disfavour the detection of planetary systems via RVs. The integration time was set to 106 times the orbital period of the outermost planet, and the integration time-step was 5% the orbital period of the innermost one. We found that the system was unstable when i ≲ 50° for set (3), that is, for the maximum eccentricity values. On the other hand, the system tolerated the full range of inclinations for sets (1) and (2) (see Fig. 6). This dynamical robustness impedes strong constraints on the system’s inclination and, consequently, on the planetary masses.

Next, we explored the stability of the system in a four-planet configuration. In the previous analysis of the two-planet configuration, we found that the system became unstable when we considered the +1_σ_ values of the eccentricities, that is, set (3) for i < 50°. In order to explore the effect of the inclinations on the stability of the two outermost planets, the two innermost planets were fixed to their nominal values, which yielded stable solutions. Then we ensure that if any instability is reached is due to the new planets added to the system. With this in mind we tested three different scenarios: (1) _e_d = 0.0 and _e_e = 0.0; (2) _e_d = 0.17 and _e_e = 0.03; (3) _e_d = 0.35 and _e_e = 0.23. We found that scenarios (2) and (3) are highly unstable for the whole range of inclinations explored. Only for the case (1), for angles ranging from 90° to ~ 72°, was the system stable. These provide strong constraints on the four-planet configuration, which allowed us to determine the planetary masses in the following ranges: 1.18 ≤ _M_c ≤ 1.24 _M_⊕; 2.89≤ _M_b ≤ 3.03 _M_⊕; 10.80 ≤ _M_d ≤ 11.35 _M_⊕; and 9.30≤ _M_e ≤ 9.70 _M_⊕. These results hint that for the GJ 273 system, if the two outermost planets are eventually confirmed, the entire system is likely composedofan Earth-mass, a super-Earth, and two mini-Neptunes in the outskirts. The stability of the system is mainly controlled by the two mini-Neptunes, which should reside in nearly circular orbits. Since the planets have not been found transiting so far, their expected radii are still highly uncertain. Based on the M–R relationships for ice-rock-iron planets given by Fortney et al. (2007), we estimate the radius of GJ 273c to be in the range of 1.04≤ _R_c ≤ 1.18 _R_⊕, considering two possible compositions: Earth-like (rock-iron as 67–33%) and 100% rocky. For GJ 273b we obtained a size range of 1.32≤ _R_b ≤ 1.83 _R_⊕, considering compositions of either 100% rocky, or a 50–50% mix of ice-rock. These choices for the compositions of GJ 273c and GJ 273b were motivated by the results obtained in Sect. 2, where we found that GJ 273c is likely a very dry object, while GJ 273b was an efficient water captor at early time scale, so consequently, it might be highly hydrated. For the outermost planets we followed the recently revisited relation given by Otegi et al. (2020), who established an M–R diagram with two populations: small-terrestrial planets (_ρ_ > 3.3 g cm−3) and volatile-rich planets (ρ < 3.3 g cm−3). This revisited M–R diagram showed that the two populations overlap in both mass (5–25 _M_⊕) and radius (2–3 _R_⊕), where planetary mass or radius alone cannot be used to distinguish between the two populations. This fact hints that for the case of the outermost planets in GJ 273, it is not possible to distinguish the two populations, consequently their radii are more degenerated than in the case of the inner planets. Considering both options, that is, small-terrestrial planets and volatile-rich planets, their corresponding radius ranges are 2.00 ≤ _R_d ≤ 3.68 _R_⊕ and 1.91 ≤ _R_e ≤ 3.35 _R_⊕, respectively.

Fig. 6Stability analysis based on the MEGNO parameter for different scenarios for the two-planetconfiguration. We explored 1000 different values of the orbital inclination ranging from 90° to 10°, for three different sets of eccentricities corresponding to their nominal values ±1_σ_ uncertainties: (1) the solid-green line corresponding to _e_b = 0.03 and _e_c = 0.05; (2) the solid-yellow line, which corresponds to _e_b = 0.10 and _e_c = 0.17; and (3) the solid-red line for _e_b = 0.19 and e_c = 0.30. For clarity, scenarios (1) and (2) have been shifted by 0.5 units from ⟨_Y (t)⟩ = 2. On the _x_-axis, for each inclination the corresponding planetary mass is shown due to the mass-inclination degeneracy.

4.3 Tidal evolution

In closely-packed systems like GJ 273 tidal interactions between the planets and the star may influence the evolution of the planetary parameters and their orbits. However, the time scales differ from one parameter to other. For example, the obliquity and planetary rotational period evolve fast, while the eccentricity and semi-major axes require time scales which may last up to the age of the host star. As a planetary orbit evolves through tidal coupling, the orbital energy is converted into tidal energy, and substantial internal heating of the planet can occur. Here, we studied the evolution of the tides in the context of the two-planet system. This choice was motivated by the fact that the outermost planets, if they exist, must not be affected by tides due their large heliocentric distances. We refer the reader to Barnes (2017) and references therein for a detailed review of the existing models used to study the evolution of tides in exoplanets.

In order to explore the evolution of GJ 273b and GJ 273c under tidal interactions, we employed the constant time-lag model (CTL), where bodies were modelled as a weakly viscous fluid (see e.g. Mignard 1979; Hut 1981; Eggleton et al. 1998; Leconte et al. 2010). We made use of the implementation of this model in MERCURY-T (Bolmont et al. 2015) and Posidonius (Blanco-Cuaresma & Bolmont 2017). Our planetary-formation model favored small planets in the range of their uncertainties, i.e., GJ 273c in the range of 1.18–2.50 _M_⊕ and GJ 273b in 2.89–6.00 _M_⊕. Hence, GJ 273c may be considered an Earth-mass/rocky planet, while GJ 273c might be a super-Earth or a mini-Neptune composed of a mix of rock and ice. With these in mind, for GJ 273c we assumed the product of the potential Love number of degree 2 and a time-lag corresponding to Earth’s value of _k_2,⊕ Δ_τ_⊕ = 213 s (Néron De Surgy & Laskar 1997). Furthermore, we assumed that the potential Love number of degree 2 and the fluid number were equal. For rocky-icy planets, the dissipation factor strongly depends on the ratio of rock to ice, where a larger amount of ice results in a larger dissipation factor. Unfortunately, these values are unknown. For example, for planets with 100% ice, the dissipation factor is considered to be larger than Earth’s by up to 10 times. For 100% rocky planets, this factor is assumed to be roughly 0.1–0.5 Earth’s value. Given that based on our formation models, the planet GJ 273b seems to be an efficient water captor at early times, we considered that the ratio of rock and ice was skewed towards ice-rich. Hence, for GJ 273b we assumed 3 × _k_2,⊕Δ_τ_⊕ (Bolmont et al. 2015, 2014; McCarthy & Castillo-Rogez 2013). This choice, despite having been made without taking into account the results found in the four-planet configuration, it is still in agreement with the results of the latter configuration, where it was foundthat GJ 273c is likely an Earth-mass planet and GJ 273b is a super-Earth. A detailed study of the planetary compositions would be desirable in order to better understand the real nature of the planets and the validity of the values assumed here.

Fig. 7Obliquity and rotational period evolution due to the effect of tides for the innermost planets GJ 273b and GJ 273c. Each set of line-style curves represents a different initial spin-rate and each set of coloured curves represents a different initial obliquity. It is shown that for any possible combination of initial conditions, in a short-time scale compared with the age of the star, both planets reached a pseudo-rotational state.

4.3.1 Obliquity and planetary rotational periods

In order to explore the evolution of the obliquity (ɛ) and the rotational periods of the planets (_P_rot), we performed a set of simulations with different initial conditions: planetary rotational periods of 10 h, 100 h, and 1000 h, combined with obliquities of 15°, 45°, and 75°. The rest of the values were fixed to the nominal values in Table 1. We found that planet GJ 273c evolved over a short time-scale ranging between 104 to 105 yr to a state with ɛ → 0° and _P_rot = 113 h. For planet GJ 273b, the evolution was slower and it needed a time-scale of about 106 −107 yr to evolve towards a state with a ɛ → 0° and _P_rot = 447 h. The results of this set of simulations are displayed in Fig. 7. Since GJ 273 is much older than these time scales, we concludethat the planets of the system are in pseudo-rotational state, that is, the planet’s rotational spin is oriented with that of the star (the obliquity is zero) and it reaches a constant rotational period. That is, the planets are tidally-locked and fairly well aligned with the host star. However, it is also important to note that the presence of a relatively large moon might perturb the pseudo-rotational state, thereby locking the obliquity into another value or generating chaotic fluctuations (see e.g. Laskar & Robutel 1993; Lissauer et al. 2012).

	Fig. 8Temporal evolution of the semi-major axis (top panel), eccentricity (mid-panel) and mean tidal heating (bottom-panel) for the innermost planets GJ 273b and GJ 273c in the system due to the effects of tides.

4.3.2 Eccentricity, semi-major axes, and tidal heating

To study the evolution of the orbits, we ran another suite of simulations, with integrations up to ~ 108 yr, assuming the nominal values given in Table 1, and both planets had starting obliquities of 15°, and rotational periods of 10 hr. From the previous analysis, we note that these choices do not make a substantial differences since the time scales for the evolution of these parameters are much shorter. The results of these simulations are displayed in Fig. 8. The semi-major axes of the planets do not show any remarkable evolution, which agrees with results presented by other authors where it is shown that the evolution of the semi-major axis due to the effect of tides might last as long as the life of the star (see e.g. Bolmont et al. 2015, 2014).

Regarding the eccentricities, for the inner most planet, GJ 273c shrinks from 0.17 to 0.12 at the end of the integration. For planet GJ 273b, the eccentricity maintains its initial value of 0.10. This hints that the tides are not strong enough to provoke the circularisation of the orbits in the integration time lapse considered in this study. Fromthis slow evolution, we can infer that the circularisation spans from hundreds to thousands millions of years.

Concerning the tidal heating, once a pseudo-rotational state is reached by a planet, which we know from the previous sub-section happens over a short time scale of <107 yr, tidal heating remains at play while the orbits are eccentric. We have shown that the period of circularisation lasts at least several 108 yr. A similar result was obtained by Barnes (2017) for Proxima Centauri-b, who also used the CTL model. The planet was found to reach the tidal-locking configuration in a time-scale of 104 –105 yr, but the orbit was likely non-circularised. Taken together, the results found here suggest that: (1) tidal heating acts over long-time scales of at least several ~108 yr or more, and (2) it is likely that the planets are tidally locked but the same side may not always face the star. For the innermost planet, the tides add significant extra heating which oscillates in the range of 100–1000 Wm−2. On the other hand, GJ 273b received a tidal heating of ~0.1 Wm−2.

5 Minor-body reservoir analogues

In the Solar System, minor bodies are the building blocks of planets. It is extensively accepted that they are the leftovers of planetary formation, which, through impacts, played an important role in the hydration process and the delivery of a variety of organic volatiles to the early Earth (see e.g. Owen & Bar-Nun 1995; Morbidelli et al. 2000; Horner et al. 2009; Hartogh et al. 2011; Snodgrass et al. 2017). This hints that minor bodies might have played a key role in the emergence of life on Earth (see e.g Oberbeck & Fogleman 1989a,b; Nisbet et al. 2001). On the other hand, it is also important to consider that a too-high impact flux of minor bodies might erode planetary atmospheres, thereby provoking oceanic evaporation, sterilising planetary surfaces, and even causing mass extinctions (see e.g. Alvarez et al. 1980; Chyba 1993; Becker et al. 2001; Ryder 2002; Hull et al. 2020). Regarding exoplanetary systems, two populations of comets orbiting _β_-Pictoris have been identified (Kiefer et al. 2014a). In addition, there is strong evidence of falling evaporating bodies, that is, exocomets, in other young systems such as Fomalhaut and HD 172555 (see e.g. Kiefer et al. 2014b; Matrà et al. 2017). There are also high-contrast, spatially resolved observations which showed an excess in the near infrared that has been attributed to the presence of exozodiacal light in the main-sequence stars, whose origin may have been the same as that which produced the zodiacal cloud in the Solar System: minor bodies’ dust particles coming from asteroid and cometary activity or disruptions (see e.g. Nesvorný et al. 2010; Absil et al. 2013; Sezestre et al. 2019). Furthermore, the recent discoveries of the interstellar objects 1I/2017 U1 (see e.g. Meech et al. 2017; Jewitt et al. 2017) and 2I/Borisov (see e.g. Fitzsimmons et al. 2019; de León et al. 2019) clearly hint at the existence of minor bodies in other planetary systems, which can be ejected presumably by some kind of dynamical perturbation, similarly to what happens in the Solar System. Moreover, Ballesteros et al. (2018) proposed that a swarm of Trojan-like asteroids could explain the mysterious behaviour of KIC 8462852 (see e.g. Boyajian et al. 2016, 2018). Despite being a very new and emerging field6, the formation of minor bodies together with planets around host stars appears to serve as evidence and, therefore, almost every planetary system may harbour such remnants, which may reside in stable orbits after the planetary-formation process. In this context, we explored the potential existence and location of regions in GJ 273 where minor bodies could reside in stable orbits for long periods of time. In our Solar System, these stable regions are the minor-body reservoirs, such as the main asteroid belt (MAB; DeMeo et al. (2015); Morbidelli et al. (2015)), the Kuiper belt (KB; Jewitt & Luu 1993; Brown et al. 2001; Morbidelli & Nesvorný 2020) and the Oort cloud (OC; Dones et al. 2004; Fernández 2005). For a full review of the characteristics and properties of these structures, we refer the reader to de León et al. (2018) and references therein. In this current analysis, we investigated the similarities between GJ 273 and our Solar System.

To perform this study, we made use of the MEGNO criterion, previously described in Sect. 2. We carried out a stability analysis by introducing a massless particle, such as a comet or an asteroid, into the GJ 273 system. Initially, we only considered the two-planet configuration. The particle was then evaluated in the a_–_e parameter space using initial values of _a_0 ∈ [0.0, 0.2] au and _e_0 ∈ [0.0, 0.99]. For simplicity, we considered the particle as coplanar, that is, i = 0°, meaning we only investigated moderately flat structures, such as MAB and KB analogues. This choice was motivated by the complex formation process of the OC, which implies the migrations of giant planets in early ages of the Solar System. Due to the many uncertainties in the planetary parameters for the two-planet configuration, we repeated this analysis for different scenarios by considering the maximum and minimum values of their mass and eccentricity, e, as constrained by the dynamic considerations presented in Sect. 4, but we kept their semi-major axis, a, and mean longitude, λ, constant.

We tested four extreme configurations (see Table 2), which covered all the possible parameter variations. The integration time was setto 104× the orbital period of the outermost body in the simulation. The time-step for the integration was set to 5% of the orbital period of the innermost body. The size of the parameter space explored was 500 × 500 pixels, that is, we studied the a_–_e space up to 250 000 times. In order to save computational time, we prematurely ended the simulation when ⟨Y (t)⟩ > 5, or when a collision of two bodies or the ejection of a particle took place. For all these chaotic scenarios we assigned a value of ⟨Y (t)⟩ = 5. We found that there are three main regions (Fig. 9) where a particle remained stable. These are labelled as Regions A, B, and C. Region A is placed between the host star and the inner planet GJ 273b, from ~ 0.01 au to ~ 0.03 au and e ≲ 0.4. Region B is situated between the two planets, from ~0.040 au to ~0.07 au, and e < 0.3. The last stable region, C, is the largest one, which is located beyond the outermost planet from ~ 0.1 au to ~ 0.2 au, where e increases progressively from 0 up to 0.4. In the Solar System, the MAB is a very stable region located between Mars and Jupiter, which contains roughly 3 × 1024 g (~ 4% mass of the Moon) of material spread among countless bodies with a size distribution given by a power law with an incremental size index between −2.5 and −3, most of whom have eccentricities less than 0.4 andinclinations below 20° (see e.g. de León et al. 2018). The objects in this region are composed mainly of silicates, metal, and carbon. Since Region B is inner to the snow line of the system (~0.65 au), we might consider this as a hot-MAB analogue, where minor bodies may reside, trapped in between the two planets. This region, while being sculpted slightly differently in each scenario, robustly remained in the system throughout all the possible configurations explored, which hints that it is a very suitable place where minor bodies may exist. Furthermore, our planetary formation models presented in Sect. 2 revealed that bodies in this region are very dry, similarly to what is found in our MAB. On the other hand, Region A suffered dramatic changes depending on the configuration considered, even disappearing when the maximum eccentricities were considered, namely, in scenarios (b) and (d). Therefore, the existence of objects in region A seems less possible, but perhaps mildly plausible. However, in a more realistic scenario, the tidal interactions with the host star might prevent the accumulation of minor bodies in this region. In the Solar System there is no population of objects inhabiting stable orbits between the innermost planet, Mercury, and the Sun, which hints that this structure does not have any known analogue. Finally, we identified a region which expanded from the outermost planet to the edge of the system considered at 0.2 au, namely, Region C. This region changed slightly in each scenario, but it remained present in all of them. Its shape was sculpted by planet GJ 273b, which was also responsible for several spines of stability, corresponding to mean motion mesonances (MMRs). The most evident MMRs have been labelled in Fig. 9, but many others with lower intensities were also present. Similar structures were also found in Kepler-47 by Hinse et al. (2015). In order to identify these MMRs, we made use of the ATLAS2BGENERAL code7, which allowed us to place resonances between a planet of mass _m_p and a mass-less particle (Gallardo 2006). Region C seems to be similar to the reservoir present in the outskirts of the Solar System, that is, the KB. In the Solar System, this region is differentiated into two areas: the trans-Neptunian belt and the scattered disk (SD). The first contains objects with eccentricities up to 0.2, while in the SD the bodies may have eccentricities up to 0.5; in both cases the inclinations are mostly below 40° (see e.g. de León et al. 2018). The objects contained in this region are largely composed of ices, such as water, nitrogen, or methane, mixed with silicates and refractory non-volatiles species. From our planetary formation models in Sect. 2, we found that Region C is composed mainly by high-hydrated objects, what could be considered also as volatile-rich. Since this region is also inner to the snow line, it might be a hot-KB analogue. Among all the resonances found in Region C, it is interesting to remark the 2:3 MMR with GJ 273b, which is analogous to the well-known 2:3 resonance of minor bodies with Neptune, which is one of the most populated resonances in the Solar System, where the dwarf planet Pluto and the Plutinos family orbit. It has been suggested that trans-Neptunian objects larger than 500 km may have a rocky core, an icy mantle, and a volatile-rich crust (McKinnon et al. 2008). In fact, the NASA’s New Horizon spacecraft visited the Pluto–Charon system, it found strong indications of the presence of a subsurface ocean on Pluto; such oceans may also be present in others dwarf planets (Moore et al. 2016). In this context, it is interesting to note that for GJ 273, this hot-KB analogue resides in the HZ of the system.

In the second step of our analysis, we considered the four-planet configuration. In this case, as shown in Sect. 4, we found that the system remained stable only for values of inclinations ranging from 90° to ~72°; hence the masses of the planets were highly constrained. We also found that the two outermost planets likely resided in near-circular orbits. Therefore, since the planetary parameters were well constrained, for the four-planet configuration we only studied one scenario (Table 2). We followed the same strategy described for the two-planet configuration but extended the initial semi-major axis range to 1.4 au. We assumed an intermediate value of the inclination of 80°. A stability map of this configuration is displayed in Fig. 10. We found that the innermost Region A was almost removed in this scenario, similar to that found for scenarios (b) and (d) shown in Fig. 9. Region B remained well established, which lends credence to the conclusion that is a very suitable region where minor bodies may reside. On the other hand, Region C, which we considered to be a hot-KB analogue in the two-planet configuration, consisted of a large stable region spanning between 0.1 to 0.6 au, with a maximum eccentricity of ~0.55 at 0.3 au. This structure may be considered as the second hot-MBA analogue in the system. From ~ 0.55 to ~ 1.0 au, there was an empty region where minor bodies would be expelled easily, and from ~ 1.0 au to the edge of the simulation, there was a stable region that may be a KB analogue, which we have labelled as Region D. Therefore, in the four-planet configuration, GJ 273b resided in between the two hot-MBA analogues. Regions C and D were sculpted by the pair of mini-Neptunes in the outskirts, which are likely to have provoked long-terms perturbations that altered the state of the minor bodies hosted therein.

This technique may be used not only to identify the location(s) of minor-body reservoirs such as MAB or KB analogues in exoplanetary systems, but also to predict the locations and ascertain the characteristics of debris disks. An example is the analysis of the dynamic environment of HR 8799 made by Contro et al. (2016), where the authors modelled the innermost debris disk and concluded that it exhibited similar features as our MAB, including the presence of Kirkwood gaps. Indeed, in the case of GJ 273, the system might harbour debris disks in these stable locations. However, the GJ 273 system is old and it is expected that the primordial dust and gas were removed a long time ago, as is the case in our own Solar System. If in the future infrared or millimetre observations of this system reveal the presence of debris disks, we may consider that these regions are collisionally active, meaning that collisions between minor bodies happen very often and produce vast quantities of dust, commonly enough to prevent its dissipation via a variety of mechanisms, such as Poynting-Robertson drag or radiation pressure, which act on short time scales (see e.g. Wyatt & Whipple 1950; Wyatt et al. 2011).

Fig. 9MEGNO maps in the a_–_e parameter space for the two-planet configuration for four different scenarios (as shown in Table 2). The dark-purple areas correspond to stable orbits where minor bodies may reside. On the other hand, yellow areas correspond to unstable regions from which a minor body would be expelled. Three main stable regions are found: A, B, and C. Vertical dashed-lines represent the MMRs found in the system between the minor bodies and the planets GJ 273b and GJ 273c.

	Fig. 10MEGNO map in the a_–_e parameter space for the four-planet configuration, as shown in Table 2. The colour scheme is the same as given in Fig. 9.

6 Habitability of GJ 273

The importance of M dwarfs in the search for life-bearing planets has been long debated. It was suggested that for many aspects, planets orbiting in the liquid-water HZs of these stars might provide stable environment for life, especially just after their first 0.5–1.0 Gyr, which consists of a period of intense radiation which could erode and damage planetary atmospheres. This, along with the feasibility of transit detection compared with more massive stars, has prioritised an M-dwarf’s planets as targets in the search of other biospheres beyond the Solar System (see e.g. Scalo et al. 2007; Yang et al. 2013; Shields et al. 2016a; O’Malley-James & Kaltenegger 2019).

In this context, the GJ 273 system is particularly interesting because it is one of the closest planetary systems to Earth hosted by an M dwarf and it contains a confirmed exoplanet in the inner region of the HZ. Indeed, the first factor to consider whether a planet is potentially habitable is its distance to the host star, in which a terrestrial-mass planet with an Earth-like atmospheric composition (CO2, H2 O, and N2) can sustain liquid water on its surface. Following Kopparapu et al. (2013a, 2014), the optimistic HZ in the GJ 273 system spans between 0.08 and 0.2 au, which establishes that GJ 273b is close to the inner border, but still remaining within it.

However, many different factors play important roles when describing the potential habitability of an exoplanet and for a detailed discussion, we refer the reader to Meadows & Barnes (2018) and references therein. Unfortunately, it is not possible to address each of these factors at the same time.

In this particular study, from the planetary formation modelling performed in Sect. 2, we found that GJ 273b could have been an efficient water captor during its first 100 Myr. Since we are considering perfect merger collisions, the water mass fraction found for this planet is an upper limit and it is likely that the real content of water is lower, but still significant, which implies a certain level of hydration.

From the dynamic analysis performed in Sect. 4, we found that GJ 273b, if the four-planet configuration is confirmed, would be compatible with a super-Earth with a mass range of 2.89 ≤ _M_b ≤ 3.03 _M_⊕, a size of 1.32 ≤ _R_b ≤ 1.83 _R_⊕, and it resides in a stable orbit. The orbital parameters of a given exoplanet, such as its semi-major axis, eccentricity, and obliquity, govern the radiation received by the planet along its orbit, which impacts its planetary climate (see e.g. Spiegel et al. 2009; Armstrong et al. 2014; Shields et al. 2016b). In the work by Kane & Torres (2017), the authors described the maximum flux received by a planet as a function of the latitude, β, the orbital eccentricity, the orbital phase, ϕ, and the obliquity, ɛ, as follows:F = L⋆4πr2cos| β−δ |,(11)

where _L_⋆ is the stellar luminosity, and δ is the stellardeclination given by:δ = ɛcos[2π(ϕ−Δϕ)],(12)

for which Δ_ϕ_ is the offset in phase between periastron and the highest solar declination in the Northern hemisphere. As we showed in Sect. 4, due of the action of tides, GJ 273b is likely in a pseudo-rotational state, that is, the obliquity is zero and it reached a constant rotational period. However, the planet’s orbit might not yet be completely circularised. Consequently, Eqs. (11) and (12) can be simplified as:δ = 0→F = L⋆4πr2cos| β |.(13)

This means that the orbital phase is no longer relevant and the distribution of flux along the latitude is the same for the whole orbit, which depends only on the planet’s eccentricity. Therefore, the maximum change in flux for a given latitude would take place between the periastron and the apastron:ΔF = L⋆4πa2e1−e2cos| β |.(14)

Considering the nominal value of the eccentricity and semi-major axis given in Table 1, we found that GJ 273b receives, at the equator 2088–1051 Wm−2 (1.53–0.77 _F_⊕), 1476–743 Wm−2 (1.08–0.54 _F_⊕) at β = ± 45°, and ~ 0 Wm−2 (0.0 _F_⊕) at the poles. These values suggest that the stellar flux received by the planet may vary by as much as 50% along its orbit, which might have strong implications in terms of the global climate of the planet.

Furthermore, in Sect. 4, we obtained the tidal heating experienced by the innermost planets in the system. The effect of tidal heating in rocky or terrestrial planets and exo-satellites may have significant implications for their habitability, notably for M dwarfs (Shoji & Kurita 2014) where the HZ and tide-regions overlap. Thus, tidal heating could either favor or work against the habitability of GJ 273b. As an example and in making an analogy with the Jupiter-moons system, Europa is a rocky body covered by a 150-km thick water-ice crust, tidal heating may maintain a subsurface water ocean. In the case of the Jovian satellite Io, the extreme violence of the tides provokes intense global volcanism and rapid resurfacing, ruling out any possibility of habitability. Indeed, in the case of Io, the tidal heating was found to be h = 2 Wm−2 and it is considered that larger heating rates would result in uninhabitable planets (Blaney et al. 1995; Bagenal et al. 2004). On the other hand, internal heating is used to explain Earth’s plate tectonics, which may actually enable planetary habitability since it drives the carbon–silicate cycle, which stabilises the temperature of the atmosphere and the levels of CO2. Indeed, from Martian geophysics studies it has been suggested that tectonic activity ceased when the internal heat dropped below 0.04 Wm−2 (Williams et al. 1997). Based on these arguments, Barnes et al. (2009) established the limits of tidal heating which may be compatible with a biosphere in the range of 0.04–2.0 Wm−2. In this context, we found that GJ 273b, with tidal heating of h = 0.1 Wm−2, is compatible with the development of life.

In Sect. 5, we highlighted the importance of minor bodies on the emergence and maintenance of life on a planet, and we identified stable locations in the system where minor bodies may reside after the planetary formation. It is beyond the scope of this paper to explore more deeply how the potential hosted bodies of the stable regions evolve and behave due to their gravitational interactions with other planets in the system, such as those that occur in the Solar System. But this study would be highly interesting due to the location of GJ 273b in between two minor-body reservoirs. An example of this are the theoretical studies of three systems HD 10180, 47 UMa, and HD 141399 by Dvorak et al. (2020), where the authors studied the dynamical evolution of hypothetical exocomet populations, and the study presented by Schwarz et al. (2017), where the authors explored the water delivery to the planet Proxima Centauri-b in consideration of a hypothetical population of exocomets. A similar study was performed by Dencs et al. (2019), where the authors investigated a synthetic late-heavy-bombardment-like event in the TRAPPIST-1 system, which could have transported a significant amount of water to the habitable planets discovered in the system. All these studies could be combined with the recent prescriptions given by Wyatt et al. (2020) to perform predictions about the evolution of GJ 273b’s atmosphere due to planetesimal impacts.

Another important factor to consider in terms of habitability are stellar flares and coronal mass ejections (CMEs), where the former deliver strong and isotropic ultraviolet (UV) radiation, and the latter form directional streams of charged particles. While a potential atmosphere on GJ 273b could protect the surface from harmful UV light, flares pose an imminent danger to surface life if these absorbers are depleted. For example, Tilley et al. (2019) found that the CMEs associated with strong flares can deplete the ozone in a potential atmosphere, leaving the flares’ UV radiation to hit the surface nearly unabsorbed. On the other hand, prebiotic chemistry research has shown that a certain amount of flaring might, in fact, be necessary to form possible precursors for life in the first place (Rimmer et al. 2018). This suggests that surface life might rely on a host stars with the right balance of flaring, not too much and not too little (see e.g. Günther et al. 2020).

7 Conclusions

In this study, we set out to better understand the real nature of the planetary system GJ 273. This system is of special interest because it is the fourth closest planetary system to the Sun orbiting an M dwarf (3.75 pc) known to host a planet in the HZ, just after Proxima Centaury (1.30 pc, Anglada-Escudé et al. 2016), Ross-128 (3.37 pc, Bonfils et al. 2018), and GJ 1061 (3.67 pc, Dreizler et al. 2020). The planets in GJ 273 have been detected so far only via RVs: two confirmed (GJ 273b and GJ 273c), and two candidates (GJ 273d and GJ 273e). We performed a comprehensive study of GJ 273’s planetary formation and water delivery history in its early times, a dynamical stability of the system and the effects of tides on the planets, as well as the dynamical environment, where we identified regions where the system may harbour minor bodies in stable orbits, namely, MAB and KB analogues. Each of these parameters has an substantial influence on the potential habitability of the planet located in the HZ, GJ 273b. We also reviewed the published observational data obtained via HARPS and TESS, where we explored the detectability of the known candidates and extra planets in the system. More specifically:

• 1.
  We analysed the process of formation and evolution of the planetary system around GJ 273, making use of _N_-body simulations. Our research allowed us to derive three important results: (1) the gaseous component of the disk played a key role in determining the physical and orbital properties of the planets of the system; (2) GJ 273c is a planet with very low water content, while GJ 273b seems to be a water world evolving within the HZ; and (3) the planetary system around GJ 273 might host more planets (in addition to the two confirmed) with different physical and orbital properties.
• 2.
  From our review of the observational data sets, we found that these other possible planets predicted by the _N_-body simulations could not have been detected with the precision and cadence of the HARPS RV data-set that confirmed GJ 273b and GJ 273c. From the TESS data, however, we found that a transit of the innermost planet GJ 273c can be ruled out. On the other hand, for the rest of the planets in the system, the TESS photometric data were insufficient to confidently rule out the presence of any transit signals.
• 3.
  We tested the stability of the system using the MEGNO criterion. We found that in the two-planet configuration we could not constrain the inclination of the system and consequently the masses of the planets remain unknown. However, in the four-planet configuration, the system was stable only for near-circular orbits of the two outermost planets, GJ 273d and GJ 273e, and minimum inclinations of the whole system down to ~72°. Hence, this scenario allowed us to unveil the system as likely composed by an Earth-mass planet, a super-Earth, and two mini-Neptunes with masses of: 1.18 ≤ _M_c ≤ 1.24 _M_⊕, 2.89 ≤ _M_b ≤ 3.03 _M_⊕, 10.80 ≤ _M_d ≤ 11.35 _M_⊕, and 9.30 ≤ _M_e ≤ 9.70 _M_⊕, respectively. In this scenario, using M–R relationships for different bulk compositions we found that the radii of the planets would be in the following ranges: 1.04 ≤ _R_c ≤ 1.18 _R_⊕, 1.32 ≤ _R_b ≤ 1.83 _R_⊕, 2.00 ≤ _R_d ≤ 3.68 _R_⊕, and 1.91 ≤ _R_e ≤ 3.35 _R_⊕, respectively.
• 4.
  Since the two innermost planets occupy close-in orbits, we performed a tidal-evolution analysis where we found that both GJ 273b and GJ 273c are likely in a pseudo-rotational state, meaning that their spins are aligned with that of the hoststar and they have reached constant rotational periods. This state was reached over a short time scale of <107 yr, which is much shorter than the age of the star, estimated to be older than a few Gyr. However, we found that the time required to circularise the orbits was much longer, which could last several 108or 109 yr, which hints that the orbits might not be fully circular and, hence, a certain amount of tidal heating may be induced.
• 5.
  Using the MEGNO criterion we studied the similarities of GJ 273 with the Solar System in terms of minor-body reservoirs in both the two- and the four-planet configurations. We found several regions where minor bodies might reside in stable orbits, forming structures like the MAB or KB. Each structure was sculpted differently in each scenario and, in particular, there were two structures present in all simulations: (1) region B, which was trapped in between planets GJ 273b and GJ 273c and may be considered as an hot-MBA analogue; and (2) region C, which may be a hot-KB analogue in the two-planet configuration or a hot-MBA analogue trapped between GJ 273b and GJ 273d in the four-planet configuration. All these structures identified in GJ 273 seem to have some similarities with those of the Solar System, but at the same time, they point out the large diversity of minor-body reservoirs that may exist.

Finally, taking together all the analyses performed in this study and considering whether the four-planet configuration is eventually confirmed, we conclude that GJ 273b is likely a super-Earth residing in a stable orbit close to the inner border of the HZ of its host star. The planet might be hydrated and it is likely in a pseudo-rotational state with in a non-circular orbit. Such a small eccentricity could provoke remarkable variations of the stellar flux received that could result in strong implications in terms of the global climate of the planet. Nevertheless, we also found that GJ 273b undergoes tidal heating compatible with the development and existence of a biosphere. The reservoirs of minor bodies could also play a role in the emergence and maintenance of life in GJ 273b. In fact, it would be desirable to further investigate how minor bodies such as asteroids and comets evolve and behave in relation to GJ 273b, in particular the impact-flux ratio, where aspects such as the watering, the delivery of organics, and the evolution of planetary atmosphere take place.

Acknowledgements

We thank the anonymous referee for the helpful comments. J.C.S. acknowledges funding support from Spanish public funds for research under projects ESP2017-87 676-2-2 and RYC-2012-09913 (Ramón y Cajal programme) of the Spanish Ministry of Science and Education. Z.M.B. acknowledges CONICYT- FONDECYT/Chile Postdoctorado 3 180 405. M.G. and E.J. are FNRS Senior Research Associates. V.V.G. is a FNRS Research Associate. M.N.G. acknowledges support from MIT’s Kavli Institute as a Torres postdoctoral fellow. We acknowledge the use of public data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center.

Appendix A Locationof the snow line

In our modelling, we followed the procedure proposed by Ciesla et al. (2015) to compute the location of the snow line, which sets the boundary between dry embryos and water-rich embryos in each of our _N_-body simulations. This procedure is described in detail below.

To do this, we assumed that water condenses in a protoplanetary disk, where its partial pressure exceeds the saturation pressure, which is associated with temperatures between 140 K and 180 K for typical disk conditions. In our work, we adopt a value of T = 140 K to define the position of the snow line in the protoplanetary disk. To determine the distance R from the central star in the disk mid-plane, where T = 140 K, we adopt the temperature profile suggested by Chiang & Goldreich (1997), which is given by:T(R) = (α4)1/4(R⋆R)1/2T⋆,(A.1)

where R is the radial coordinate in the disk mid-plane, _R_⋆ and _T_⋆ are the radius and effective temperature of the star, respectively, and α is the grazing angle at which the starlight strikes the disk. According to this, the location of the snow line (_R_ice) will be given by:Rice = (α4)1/2R⋆(T⋆140 K)2.(A.2)

As the reader can see, _R_ice strongly depends on the radius and the effective temperature of the central star, which evolve over time. Following the stellar evolutionary models developed by Siess et al. (2000), we made use of Eq. (A.2) to determine the temporal evolution of the snow line in a protoplanetary disk around a 0.3 _M_⊙ star, which is illustrated in Fig. A.1. The five curves presented in the figure correspond to different values of the star’s initial metallicity, Z, (in terms of mass fraction) adopted in the stellar models of Siess et al. (2000). It is important to note that we assume a value of α = 0.045 in the construction of these five curves, which is in agreement with the expressions given by Chiang & Goldreich (1997), and with the value used by Ciesla et al. (2015).

Like Ciesla et al. (2015), we set the location of the snow line at 1 Myr. We consider that this choice may be somewhat arbitrary since the snow line will continually migrate inward as the system evolves. However, our selection simply defines a boundary between the dry embryos and water-rich embryos in our simulated protoplanetary disk. According to the results illustrated in Fig. A.1, the snow line at 1 Myr is located at 0.65 au for all values of the star’s initial metallicity, Z, adopted by Siess et al. (2000). Table A.1 shows the values of the luminosity, temperature, and radius of a 0.3 _M_⊙ star at 1 Myr used in our work to compute the location of the snow line, which were extracted from the different stellar models of Siess et al. (2000).

	Fig. A.1Location of the snow line as a function of time for a 0.3 _M_⊙ star based on the stellar models developed by Siess et al. (2000). The five curves correspond to different values of the star initial metallicity, Z, (in terms of mass fraction). The gray dashed line indicates a value of 0.65 au.

Table A.1

Luminosity, temperature, and radius at 1 Myr of a 0.3 _M_⊙ star for different values of the star initial metallicity (in terms of mass fraction) (Siess et al. 2000).

From Fig. A.1, it is evident that the snow line changes over time. However, we follow Ciesla et al. (2015) and other previous studies based on _N_-body simulations, and assign a single location for the snow line (_R_ice = 0.65 au) for the whole duration of the numerical experiments. As we previously mentioned, the snow line is defined to determine a boundary between dry embryos and water-rich embryos in the protoplanetary disk.

Appendix B Detrended models

Fig. B.1Detrended models applied to the TESS data in our search for threshold-crossing events. We used the bi-weight method with different window-sizes. Each panel indicates the window-size (in units of d) in the top. The black points show the PDC-SAP fluxes, the solid-orange line are the identified trends, and the teal points are the detrended data used as input into the TLS. In all cases, the _y_-axis shows the normalised flux, and the _x_-axis is units of Barycentric TESS Julian Date.

References

1. Absil, O., Defrère, D., Coudé Du Foresto, V., et al. 2013, A&A, 555, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
2. Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
3. Alvarez, L. W., Alvarez, W., Asaro, F., & Michel, H. V. 1980, Science, 208, 1095 [NASA ADS] [CrossRef] [Google Scholar]
4. Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2010, ApJ, 723, 1241 [NASA ADS] [CrossRef] [Google Scholar]
5. Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
6. Armstrong, J., Barnes, R., Domagal-Goldman, S., et al. 2014, Astrobiology, 14, 277[Google Scholar]
7. Armstrong, D. J., Pugh, C. E., Broomhall, A. M., et al. 2016, MNRAS, 455, 3110 [NASA ADS] [CrossRef] [Google Scholar]
8. Astudillo-Defru, N., Delfosse, X., Bonfils, X., et al. 2017a, A&A, 600, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
9. Astudillo-Defru, N., Forveille, T., Bonfils, X., et al. 2017b, A&A, 602, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
10. Bagenal, F., Dowling, T. E. T. E., & McKinnon, W. B. 2004, Jupiter : the Planet, Satellites, and Magnetosphere (Cambridge: Cambridge University Press), 719[Google Scholar]
11. Ballesteros, F. J., Arnalte-Mur, P., Fernandez-Soto, A., & Martínez, V. J. 2018, MNRAS, 473, L21 [NASA ADS] [CrossRef] [Google Scholar]
12. Barnes, R. 2017, Celest. Mech. Dyn. Astron., 129, 509[Google Scholar]
13. Barnes, R., Jackson, B., Greenberg, R., & Raymond, S. N. 2009, ApJ, 700, L30 [NASA ADS] [CrossRef] [Google Scholar]
14. Barnes, S. A. 2010, ApJ, 722, 222 [NASA ADS] [CrossRef] [Google Scholar]
15. Bate, M. R. 2018, MNRAS, 475, 5618[Google Scholar]
16. Becker, L., Poreda, R. J., Hunt, A. G., Bunch, T. E., & Rampino, M. 2001, Science, 291, 1530 [CrossRef] [Google Scholar]
17. Benettin, G., Galgani, L., Giorgilli, A., & Strelcyn, J. M. 1980, Meccanica, 15, 9 [NASA ADS] [CrossRef] [Google Scholar]
18. Benítez-Llambay, P., Masset, F., Koenigsberger, G., & Szulágyi, J. 2015, Nature, 520, 63 [NASA ADS] [CrossRef] [Google Scholar]
19. Blanco-Cuaresma, S., & Bolmont, E. 2017, EWASS Special Session 4 (2017): Star-planet interactions (EWASS-SS4-2017)[Google Scholar]
20. Blaney, D., Johnson, T. V., Matson, D. L., & Veeder, G. J. 1995, Icarus, 113, 220 [NASA ADS] [CrossRef] [Google Scholar]
21. Bolmont, E., Raymond, S. N., von Paris, P., et al. 2014, ApJ, 793, 3 [CrossRef] [Google Scholar]
22. Bolmont, E., Raymond, S. N., Leconte, J., Hersant, F., & Correia, A. C. M. 2015, A&A, 583, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
23. Bonfils, X., Delfosse, X., Udry, S., et al. 2013, A&A, 549, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
24. Bonfils, X., Astudillo-Defru, N., Díaz, R., et al. 2018, A&A, 613, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
25. Boyajian, T. S., Von Braun, K., Van Belle, G., et al. 2012, ApJ, 757, 112 [NASA ADS] [CrossRef] [Google Scholar]
26. Boyajian, T. S., LaCourse, D. M., Rappaport, S. A., et al. 2016, MNRAS, 457, 3988 [NASA ADS] [CrossRef] [Google Scholar]
27. Boyajian, T. S., Alonso, R., Ammerman, A., et al. 2018, ApJ, 853, L8 [NASA ADS] [CrossRef] [Google Scholar]
28. Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127[Google Scholar]
29. Brown, M. E., & Brown, E. M. 2001, AJ, 121, 2804[Google Scholar]
30. Browning, M. K., Basri, G., Marcy, G. W., West, A. A., & Zhang, J. 2010, AJ, 139, 504 [NASA ADS] [CrossRef] [Google Scholar]
31. Burdanov, A., Delrez, L., Gillon, M., & Jehin, E. 2018, in Handbook of Exoplanets (Berlin: Springer International Publishing), 1007 [CrossRef] [Google Scholar]
32. Chabrier, G., & Baraffe, I. 2000, ARA&A, 38, 337 [NASA ADS] [CrossRef] [Google Scholar]
33. Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
34. Chambers, J. E. 2013, Icarus, 224, 43 [NASA ADS] [CrossRef] [Google Scholar]
35. Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525 [NASA ADS] [CrossRef] [Google Scholar]
36. Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [NASA ADS] [CrossRef] [Google Scholar]
37. Chyba, C. F. 1993, Geochim. Cosmochim. Acta, 57, 3351 [CrossRef] [Google Scholar]
38. Ciesla, F. J., Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 804, 9 [NASA ADS] [CrossRef] [Google Scholar]
39. Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698 [NASA ADS] [CrossRef] [Google Scholar]
40. Cincotta, P., & Simó, C. 1999, Celest. Mech. Dyn. Astron., 73, 195 [NASA ADS] [CrossRef] [Google Scholar]
41. Cincotta, P. M., & Simó, C. 2000, A&AS, 147, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
42. Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Phys. D Nonlinear Phenomena, 182, 151 [NASA ADS] [CrossRef] [Google Scholar]
43. Contro, B., Horner, J., Wittenmyer, R. A., Marshall, J. P., & Hinse, T. C. 2016, MNRAS, 463, 191 [NASA ADS] [CrossRef] [Google Scholar]
44. Cresswell, P., & Nelson, P. R. 2008, A&A, 482, 677 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
45. Crossfield, I. J. M., Waalkes, W., Newton, E. R., et al. 2019, ApJ, 883, L16 [CrossRef] [Google Scholar]
46. Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, Vizier: II/246[Google Scholar]
47. Darriba, L. A., De Elía, G. C., Guilera, O. M., & Brunini, A. 2017, A&A, 607, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
48. De Elía, G. C., Guilera, O. M., & Brunini, A. 2013, A&A, 557, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
49. de León, J., Licandro, J., & Pinilla-Alonso, N. 2018, in Handbook of Exoplanets (Berlin: Springer International Publishing), 395 [CrossRef] [Google Scholar]
50. de León, J., Licandro, J., Serra-Ricart, M., et al. 2019, Res. Notes AAS, 3, 131[Google Scholar]
51. Delfosse, X., Forveille, T., Ségransan, D., et al. 2000, A&A, 364, 217 [NASA ADS] [Google Scholar]
52. Delrez, L., Gillon, M., Queloz, D., et al. 2018, SPIE Conf. Ser., 10700, 107001I[Google Scholar]
53. DeMeo, F. E., Alexander, C. M., Walsh, K. J., Chapman, C. R., & Binzel, R. P. 2015, Asteroids IV (Tucson: University of Arizona Press), 13[Google Scholar]
54. Dencs, Z., Regály, Z., Dencs, Z., & Regály, Z. 2019, MNRAS, 487, 2191[Google Scholar]
55. Dones, L., Weissman, P. R., Levison, H. F., & Duncan, M. J. 2004, Comets II (Tucson: University of Arizona Press), 153[Google Scholar]
56. Dreizler, S., Jeffers, S. V., Rodríguez, E., et al. 2020, MNRAS, 493, 536 [NASA ADS] [CrossRef] [Google Scholar]
57. Dressing, C. D., & Charbonneau, D. 2013, ApJ, 767, 95 [NASA ADS] [CrossRef] [Google Scholar]
58. Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45 [NASA ADS] [CrossRef] [Google Scholar]
59. Dugaro, A., De Elía, G. C., Brunini, A., & Guilera, O. M. 2016, A&A, 596, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
60. Dugaro, A., de Elía, G. C., & Darriba, L. A. 2019, A&A, 632, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
61. Dvorak, R., Maindl, T. I., Burger, C., Schäfer, C., & Speith, R. 2015, Nonlinear Phenomena in Complex Systems, 18, 310[Google Scholar]
62. Dvorak, R., Loibnegger, B., & Cuntz, M. 2020, The Trans-Neptunian Solar System (Amsterdam: Elsevier), 331 [CrossRef] [Google Scholar]
63. Eggleton, P. P., Kiseleva, L. G., & Hut, P. 1998, ApJ, 499, 853 [NASA ADS] [CrossRef] [Google Scholar]
64. Eisner, N. L., Barragán, O., Aigrain, S., et al. 2020, MNRAS, 494, 750 [NASA ADS] [CrossRef] [Google Scholar]
65. Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146[Google Scholar]
66. Fernández, J. A. 2005, Astrophys. Space Sci. Lib., 328, 103[Google Scholar]
67. Figueira, P., Marmier, M., Boué, G., et al. 2012, A&A, 541, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
68. Fitzsimmons, A., Hainaut, O., Meech, K., et al. 2019, ApJ, 885, L9 [NASA ADS] [CrossRef] [Google Scholar]
69. Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661 [NASA ADS] [CrossRef] [Google Scholar]
70. Gaidos, E., Mann, A. W., Lépine, S., et al. 2014, MNRAS, 443, 2561 [NASA ADS] [CrossRef] [Google Scholar]
71. Gallardo, T. 2006, Icarus, 184, 29[Google Scholar]
72. Gatewood, G. 2008, AJ, 136, 452 [NASA ADS] [CrossRef] [Google Scholar]
73. Gillon, M. 2018, Nat. Astron., 2, 344 [NASA ADS] [CrossRef] [Google Scholar]
74. Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
75. Gillon, M., Triaud, A. H., Demory, B. O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
76. Goździewski,K., Bois, E., Maciejewski, A. J., & Kiseleva-Eggleton, L. 2001, A&A, 378, 569 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
77. Guilet, J., Baruteau, C., & Papaloizou, J. C. 2013, MNRAS, 430, 1764 [NASA ADS] [CrossRef] [Google Scholar]
78. Günther, M. N., Pozuelos, F. J., Dittmann, J. A., et al. 2019, Nat. Astron., 3, 1099[Google Scholar]
79. Günther, M. N., Zhan, Z., Seager, S., et al. 2020, AJ, 159, 60 [CrossRef] [Google Scholar]
80. Hartogh, P., Lis, D. C., Bockelée-Morvan, D., et al. 2011, Nature, 478, 218 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
81. Hawley, S. L., Gizis, J. E., & Reid, I. N. 1996, AJ, 112, 2799 [NASA ADS] [CrossRef] [Google Scholar]
82. Heller, R., Hippke, M., & Rodenbeck, K. 2019, A&A, 627, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
83. Henry, T. J., Kirkpatrick, J. D., & Simons, D. A. 1994, AJ, 108, 1437 [NASA ADS] [CrossRef] [Google Scholar]
84. Hinse, T. C., Haghighipour, N., Kostov, V. B., & Gózdziewski, K. 2015, AJ, 799. 88 [NASA ADS] [CrossRef] [Google Scholar]
85. Hippke, M., & Heller, R. 2019, A&A, 623, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
86. Hippke, M., David, T. J., Mulders, G. D., & Heller, R. 2019, AJ, 158, 143 [NASA ADS] [CrossRef] [Google Scholar]
87. Horner, J., Mousis, O., Petit, J. M., & Jones, B. W. 2009, Planet. Space Sci., 57, 1338 [NASA ADS] [CrossRef] [Google Scholar]
88. Horner, J., Wittenmyer, R. A., Wright, D. J., et al. 2019, AJ, 158, 100 [NASA ADS] [CrossRef] [Google Scholar]
89. Hull, P. M., Bornemann, A., Penman, D. E., et al. 2020, Science, 367, 266 [NASA ADS] [CrossRef] [Google Scholar]
90. Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
91. Jehin, E., Gillon, M., Queloz, D., et al. 2018, Messenger, 174, 2 [NASA ADS] [Google Scholar]
92. Jenkins, J. S., Jones, H. R., Goździewski, K., et al. 2009, MNRAS, 398, 911 [NASA ADS] [CrossRef] [Google Scholar]
93. Jenkins, J. S., Pozuelos, F. J., Tuomi, M., et al. 2019, MNRAS, 490, 5585 [CrossRef] [Google Scholar]
94. Jewitt, D., & Luu, J. 1993, Nature, 362, 730 [NASA ADS] [CrossRef] [Google Scholar]
95. Jewitt, D., Luu, J., Rajagopal, J., et al. 2017, ApJ, 850, L36 [NASA ADS] [CrossRef] [Google Scholar]
96. Johnstone, C. P. 2020, ApJ, 890, 79[Google Scholar]
97. Kane, S. R., & Torres, S. M. 2017, AJ, 154, 204 [NASA ADS] [CrossRef] [Google Scholar]
98. Kiefer, F., Des Etangs, A. L., Boissier, J., et al. 2014a, Nature, 514, 462 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
99. Kiefer, F., Lecavelier Des Etangs, A., C Augereau, J., et al. 2014b, A&A, 561, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
100. Koen, C., Kilkenny, D., van Wyk, F., & Marang, F. 2010, MNRAS, 403, 1949[Google Scholar]
101. Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013a, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
102. Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013b, ApJ, 770, 11 [NASA ADS] [CrossRef] [Google Scholar]
103. Kopparapu, R. K., Ramirez, R. M., Schottelkotte, J., et al. 2014, ApJ, 787, L29 [NASA ADS] [CrossRef] [Google Scholar]
104. Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, AJ, 158, 32 [NASA ADS] [CrossRef] [Google Scholar]
105. Laskar, J., & Robutel, P. 1993, Nature, 361, 608 [NASA ADS] [CrossRef] [Google Scholar]
106. Leconte, J., Chabrier, G., Baraffe, I., & Levrard, B. 2010, A&A, 516, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
107. Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, 745, 79 [NASA ADS] [CrossRef] [Google Scholar]
108. Lightkurve Collaboration (L., Cardoso, et al.) 2018, Astrophysics Source Code Library [record: 1812.013][Google Scholar]
109. Lissauer, J. J. 2007, ApJ, 660, L149 [NASA ADS] [CrossRef] [Google Scholar]
110. Lissauer, J. J., Barnes, J. W., & Chambers, J. E. 2012, Icarus, 217, 77 [NASA ADS] [CrossRef] [Google Scholar]
111. Lodders, K. 2003, ApJ, 591, 1220 [NASA ADS] [CrossRef] [Google Scholar]
112. Lodders, K., Palme, H., & Gail, H.-P. 2009, in Landolt Börnstein (Berlin: Springer), 4B, 712[Google Scholar]
113. Luger, R., & Barnes, R. 2015, Astrobiology, 15, 119 [NASA ADS] [CrossRef] [Google Scholar]
114. Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [NASA ADS] [CrossRef] [Google Scholar]
115. Mamajek, E. E., Usuda, T., Tamura, M., & Ishii, M. 2009, AIP Conf., 1158, 3 [NASA ADS] [CrossRef] [Google Scholar]
116. Matrà, L., MacGregor, M. A., Kalas, P., et al. 2017, ApJ, 842, 9 [NASA ADS] [CrossRef] [Google Scholar]
117. McCarthy, C.,& Castillo-Rogez, J. C. 2013 The Science of Solar System Ices (New York, NY: Springer), 183[Google Scholar]
118. McKinnon, W. B., Prialnik, D., Stern, S. A., & Coradini, A. 2008, The Solar System Beyond Neptune (Tucson: University of Arizona Press), 213[Google Scholar]
119. Meadows, V. S., & Barnes, R. K. 2018, in Handbook of Exoplanets (Cham: Springer International Publishing), 2771 [CrossRef] [Google Scholar]
120. Meech, K. J., Weryk, R., Micheli, M., et al. 2017, Nature, 552, 378 [NASA ADS] [CrossRef] [Google Scholar]
121. Meibom, S., Barnes, S. A., Platais, I., et al. 2015, Nature, 517, 589 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
122. Mignard, F. 1979, Moon Planets, 20, 301 [NASA ADS] [CrossRef] [Google Scholar]
123. Miguel, Y., Guilera, O. M., & Brunini, A. 2011, MNRAS, 417, 314 [CrossRef] [Google Scholar]
124. Moore, J. M., McKinnon, W. B., Spencer, J. R., et al. 2016, Science, 351, 1284 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
125. Morbidelli, A.,& Nesvorný, D. 2020, The Trans-Neptunian Solar System (Amsterdam: Elsevier), 25[Google Scholar]
126. Morbidelli, A., Chambers, J., Lunine, J. I., et al. 2000, Meteorit. Planet. Sci., 35, 1309 [NASA ADS] [CrossRef] [Google Scholar]
127. Morbidelli, A., Walsh, K. J., O’Brien, D. P., Minton, D. A., & Bottke, W. F. 2015, Asteroids IV (Tucson: University of Arizona Press), 493[Google Scholar]
128. Mordasini, C., Alibert, Y., Benz, W., & Naef, D. 2009, A&A, 501, 1161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
129. Néron De Surgy, O., & Laskar, J. 1997, A&A, 318, 975[Google Scholar]
130. Nesvorný, D., Jenniskens, P., Levison, H. F., et al. 2010, ApJ, 713, 816 [NASA ADS] [CrossRef] [Google Scholar]
131. Neves, V., Bonfils, X., Santos, N. C., et al. 2013, A&A, 551, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
132. Nisbet, E. G., Sleep, N. H., Nisbet, E. G., & Sleep, N. H. 2001, Nature, 409, 1083 [NASA ADS] [CrossRef] [Google Scholar]
133. Nowak, G., Luque, R., Parviainen, H., et al. 2020, A&A, submitted [arXiv:2003.01140][Google Scholar]
134. Nutzman, P., & Charbonneau, D. 2008, PASP, 120, 317 [NASA ADS] [CrossRef] [Google Scholar]
135. Oberbeck, V. R., & Fogleman, G. 1989a, Nature, 339, 434 [NASA ADS] [CrossRef] [Google Scholar]
136. Oberbeck, V. R., & Fogleman, G. 1989b, Orig. Life Evol. Biosph., 19, 549 [NASA ADS] [CrossRef] [Google Scholar]
137. O’Brien, D. P., Morbidelli, A., & Levison, H. F. 2006, Icarus, 184, 39 [NASA ADS] [CrossRef] [Google Scholar]
138. O’Malley-James, J. T., & Kaltenegger, L. 2019, MNRAS, 485, 5598 [CrossRef] [Google Scholar]
139. Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
140. Owen, T., & Bar-Nun, A. 1995, Icarus, 116, 215 [NASA ADS] [CrossRef] [Google Scholar]
141. Paardekooper, S. J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950 [NASA ADS] [CrossRef] [Google Scholar]
142. Paardekooper, S. J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293 [NASA ADS] [CrossRef] [Google Scholar]
143. Quinn, S. N., Becker, J. C., Rodriguez, J. E., et al. 2019, AJ, 158, 177 [CrossRef] [Google Scholar]
144. Quintana, E. V., Barclay, T., Borucki, W. J., Rowe, J. F., & Chambers, J. E. 2016, ApJ, 821, 126 [NASA ADS] [CrossRef] [Google Scholar]
145. Quirrenbach, A., Amado, P. J., Mandel, H., et al. 2010, SPIE, 7735, 773513 [NASA ADS] [Google Scholar]
146. Raymond, S. N., Quinn, T., & Lunine, J. I. 2004, Icarus, 168, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
147. Raymond, S. N., O’Brien, D. P., Morbidelli, A., & Kaib, N. A. 2009, Icarus, 203, 644 [NASA ADS] [CrossRef] [Google Scholar]
148. Rein, H., & Liu, S. F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
149. Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 [NASA ADS] [CrossRef] [Google Scholar]
150. Reiners, A., Zechmeister, M., Caballero, J. A., et al. 2018, A&A, 612, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
151. Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, J. Astron. Teles. Instrum. Syst., 1, 014003 [NASA ADS] [CrossRef] [Google Scholar]
152. Rimmer, P. B., Xu, J., Thompson, S. J., et al. 2018, Sci. Adv., 4, 3302 [NASA ADS] [CrossRef] [Google Scholar]
153. Ronco, M. P., Guilera, O. M., & de Elía, G. C. 2017, MNRAS, 471, 2753 [NASA ADS] [CrossRef] [Google Scholar]
154. Ryder, G. 2002, J. Geophys. Res. E Planets, 107, 6 [NASA ADS] [CrossRef] [Google Scholar]
155. Scalo, J., Kaltenegger, L., Segura, A., et al. 2007, Astrobiology, 7, 85 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
156. Schwarz, R., Bazso, A., Georgakarakos, N., et al. 2017, MNRAS, 480, 3595[Google Scholar]
157. Sezestre, É., Augereau, J.-C., & Thébault, P. 2019, A&A, 626, A2 [CrossRef] [EDP Sciences] [Google Scholar]
158. Shields, A. L., Ballard, S., & Johnson, J. A. 2016a, Phys. Rep., 663, 1 [NASA ADS] [CrossRef] [Google Scholar]
159. Shields, A. L., Barnes, R., Agol, E., et al. 2016b, Astrobiology, 16, 443 [NASA ADS] [CrossRef] [Google Scholar]
160. Shoji, D., & Kurita, K. 2014, ApJ, 789, 3 [NASA ADS] [CrossRef] [Google Scholar]
161. Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593[Google Scholar]
162. Snodgrass, C., Agarwal, J., Combi, M., et al. 2017, A&ARv, 25, 1 [NASA ADS] [CrossRef] [Google Scholar]
163. Soderblom, D. R. 2010, ARA&A, 48, 581 [NASA ADS] [CrossRef] [Google Scholar]
164. Spiegel, D. S., Menou, K., & Scharf, C. A. 2009, ApJ, 691, 596[Google Scholar]
165. Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388 [NASA ADS] [CrossRef] [Google Scholar]
166. Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef] [Google Scholar]
167. Testi, L., Natta, A., Scholz, A., et al. 2016, A&A, 593, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
168. Tilley, M. A., Segura, A., Meadows, V., Hawley, S., & Davenport, J. 2019, Astrobiology, 19, 64[Google Scholar]
169. Tremaine, S., & Dong, S. 2012, AJ, 143, 94 [NASA ADS] [CrossRef] [Google Scholar]
170. Tuomi, M., Jones, H. R. A., Barnes, J. R., Anglada-Escudé, G., & Jenkins, J. S. 2014, MNRAS, 441, 1545 [NASA ADS] [CrossRef] [Google Scholar]
171. Tuomi, M., Jones, H. R. A., Butler, R. P., et al. 2019, AAS J., submitted [arXiv:1906.04644][Google Scholar]
172. Williams, D. M., Kasting, J. F., & Wade, R. A. 1997, Nature, 385, 234 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
173. Winters, J. G., Hambly, N. C., Jao, W. C., et al. 2015, AJ, 149, 1 [NASA ADS] [CrossRef] [Google Scholar]
174. Wood, J., Horner, J., Hinse, T. C., & Marsden, S. C. 2017, AJ, 153, 245 [NASA ADS] [CrossRef] [Google Scholar]
175. Wyatt, S. P., & Whipple, F. L. 1950, ApJ, 111, 134 [NASA ADS] [CrossRef] [Google Scholar]
176. Wyatt, M. C., Clarke, C. J., & Booth, M. 2011, Celest. Mech. Dyn. Astron., 111, 1 [NASA ADS] [CrossRef] [Google Scholar]
177. Wyatt, M. C., Kral, Q., & Sinclair, C. A. 2020, MNRAS, 491, 782 [NASA ADS] [Google Scholar]
178. Yang, J., Cowan, N. B., & Abbot, D. S. 2013, ApJ, 771, L45 [NASA ADS] [CrossRef] [Google Scholar]
179. Zain, P. S., De Elía, G. C., Ronco, M. P., & Guilera, O. M. 2018, A&A, 609, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
180. Zechmeister, M., Kürster, M., & Endl, M. 2009, A&A, 505, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

2

Following Raymond et al. (2009), a planet’s accretion seed is defined as the larger body in each of its collisions.

3

For our periodograms we used a statistical estimator to the maximise the differential likelihood function (D ln L = ln _L_null − ln _L_mod). In particular, it evaluated how well represented the data were when fitting the model (+1 planet, ln _L_mod) compared to the initial hypothesis (i.e. +0 planets, ln _L_mod). See Anglada-Escudé et al. (2016) for more details.

4

SHERLOCK code is available upon request.

5

The dynamic simulations in the context of the four-planet configuration yielded as the most likely scenario a planetary system with minimum inclinations down to 72°, i.e., composed of an Earth-mass planet, a super-Earth and two mini-Neptunes. To obtain the appropriate range of the radii, we applied the corresponding mass–radius (M–R) relationships; see Sect. 4 for more details.

6

From 13 May through May 17, 2019 at Leiden in the Netherlands, an exocomet conference was held for first time: a Lorentz workshop entitle ExoComets–Understanding the composition of Planetary Building Blocks which brought together the Solar System comet community, exoplanet researchers, and astrobiologists.

All Tables

Table A.1

Luminosity, temperature, and radius at 1 Myr of a 0.3 _M_⊙ star for different values of the star initial metallicity (in terms of mass fraction) (Siess et al. 2000).

All Figures

Fig. 1Mass distribution of the planets formed in the _N_-body simulations after 100 Myr of Set 1 (left panel), which only model the post-gas stage, and Set 2 (right panel), which include the effects of the gas disk. In every panel, the black circles illustrate the two planets confirmed that orbit GJ 273 (Astudillo-Defru et al. 2017b), while the color code represents the final fraction of water by mass of the planets resulting from our numerical experiments. Finally, the light blue shaded region refers to the optimistic HZ derived from the model developed by Kopparapu et al. (2013a,b).

In the text

Fig. 2Temporal evolution of semi-major axis as a function of the eccentricity in the “Set 2” _N_-body simulations, which illustrates the most likely planetary-formation scenario. The planetary embryos (and the resulting planets) are represented by circles, which show a colour code that indicates their fraction of water by mass. The solid blue lines illustrate curves of constant pericentric distance (q = 0.08 au) and apocentric distance (Q = 0.2 au), while the light blue shaded region represents the optimistic HZ of the system.

In the text

	Fig. 3Percentage of well-detected sinusoids from a total of N = 1000 experiments (i.e. completeness, C) versus the initial amplitude. The case of simulated sinusoids with P = 652.3 d is illustrated here as an example of how the C is calculate. The black line indicates the polynomial model fitted to get the initial amplitude (K i) that corresponds to a 90% completeness (red cross).
In the text

Fig. 4Minimum (lower-limit) masses of the planetary signals depending on their semi-major axes (black dots). Circular orbits were assumed. Only planets within the grey areas could be detected using the HARPS data set that confirmed GJ 273b and GJ 273c. Blue crosses, which correspond to the masses of the hypothetical planets predicted by the dynamical models (see last panel of Fig. 2) are belowthis detection threshold. The inset shows a close-up of the GJ 273b and GJ 273c analogues predicted by the dynamical models. The minimum GJ 273b and GJ 273c masses reported by Astudillo-Defru et al. (2017b) are highlighted with red crosses and appear on or above the detection threshold. Open circles indicate experiments for which 90% completeness was not reached even when running the experiments with K i up to 6.75 ms−1. They correspond to high power regions of the spectral window shown in the top panel.

In the text

Fig. 5Injection-and-recovery test performed to check the detectability of the innermost planets GJ 273b and GJ 273c in the current set of TESS data. We injected synthetic transits into the available photometric data corresponding to each planet by taking into account the uncertainties in their physical properties. We explored a total of1080 different scenarios. We found that planet GJ 273c showed a recovery rate of 100%, which clearly suggests that this planet is not transiting. On the other hand, planet GJ 273b showed a recovery rate of ~ 40%, implying its possible detection in future observations. We carried out two other experiments to check the detectability of the two outermost planets GJ 273d and GJ 273e, finding recovery rates < 2%; see the main text for more details.

In the text

Fig. 6Stability analysis based on the MEGNO parameter for different scenarios for the two-planetconfiguration. We explored 1000 different values of the orbital inclination ranging from 90° to 10°, for three different sets of eccentricities corresponding to their nominal values ±1_σ_ uncertainties: (1) the solid-green line corresponding to _e_b = 0.03 and _e_c = 0.05; (2) the solid-yellow line, which corresponds to _e_b = 0.10 and _e_c = 0.17; and (3) the solid-red line for _e_b = 0.19 and e_c = 0.30. For clarity, scenarios (1) and (2) have been shifted by 0.5 units from ⟨_Y (t)⟩ = 2. On the _x_-axis, for each inclination the corresponding planetary mass is shown due to the mass-inclination degeneracy.

In the text

Fig. 7Obliquity and rotational period evolution due to the effect of tides for the innermost planets GJ 273b and GJ 273c. Each set of line-style curves represents a different initial spin-rate and each set of coloured curves represents a different initial obliquity. It is shown that for any possible combination of initial conditions, in a short-time scale compared with the age of the star, both planets reached a pseudo-rotational state.

In the text

	Fig. 8Temporal evolution of the semi-major axis (top panel), eccentricity (mid-panel) and mean tidal heating (bottom-panel) for the innermost planets GJ 273b and GJ 273c in the system due to the effects of tides.
In the text

Fig. 9MEGNO maps in the a_–_e parameter space for the two-planet configuration for four different scenarios (as shown in Table 2). The dark-purple areas correspond to stable orbits where minor bodies may reside. On the other hand, yellow areas correspond to unstable regions from which a minor body would be expelled. Three main stable regions are found: A, B, and C. Vertical dashed-lines represent the MMRs found in the system between the minor bodies and the planets GJ 273b and GJ 273c.

In the text

	Fig. 10MEGNO map in the a_–_e parameter space for the four-planet configuration, as shown in Table 2. The colour scheme is the same as given in Fig. 9.
In the text

	Fig. A.1Location of the snow line as a function of time for a 0.3 _M_⊙ star based on the stellar models developed by Siess et al. (2000). The five curves correspond to different values of the star initial metallicity, Z, (in terms of mass fraction). The gray dashed line indicates a value of 0.65 au.
In the text

Fig. B.1Detrended models applied to the TESS data in our search for threshold-crossing events. We used the bi-weight method with different window-sizes. Each panel indicates the window-size (in units of d) in the top. The black points show the PDC-SAP fluxes, the solid-orange line are the identified trends, and the teal points are the detrended data used as input into the TLS. In all cases, the _y_-axis shows the normalised flux, and the _x_-axis is units of Barycentric TESS Julian Date.

In the text