© The Authors 2025

1 Introduction

The GHZ was defined by Gonzalez et al. (2001) as the region within the Galaxy where the abundance of metals is sufficiently high to favour the formation and evolution of Earth-like planets. In this context, the chemical evolution of the Galaxy plays a crucial and prominent role, as demonstrated in the work of Lineweaver et al. (2004). Their study presented, for the first time, a map of the habitable zone (HZ) as a function of both time and Galactocentric distance. It revealed that the most likely regions for the presence of life in our Galaxy are concentrated in an annular region, centred at 8 kpc from the Galactic Centre, 2 kpc wide, which gradually expands over time. They emphasised that in addition to the presence of a host star and a sufficient quantity of heavy elements to form terrestrial planets, the absence of the disruptive effects of nearby supernova (SN) explosions is a key factor to allow for habitability. Indeed, it is well established that a SN emits intense radiation that can reach the atmosphere of Earth-like planets, causing stratospheric ozone depletion. This allows for ultraviolet (UV) flux from the host star to penetrate the surface and oceans, potentially damaging genetic material (DNA), which can lead to mutations, cell death, and ultimately planetary sterilisation (Gehrels et al. 2003; Ellis & Schramm 1995). Moreover, effects are also expected at the global climate level and, consequently, across the entire biosphere - and not necessarily limited to the ozone layer (e.g. Kirkby et al. 2011; Svensmark 2022, 2023). However, the precise effects of SN explosions are not yet fully understood and the optimal method for incorporating them into GHZ models remains an open issue. Even if planetary life were to be destroyed by a nearby SN explosion, it could reappear after several hundred million years. Life shows remarkable resilience, and a cosmic catastrophe may even accelerate the evolution of life forms that are presently unknown (i.e. Krug & Jablonski 2012; Sloan et al. 2017). Moreover, life in deep water could evolve without suffering the destructive effects of a SN explosion.

In recent years, several chemical evolution models have built upon the pioneering work of Lineweaver et al. (2004) presenting new maps of Galactic habitability (e.g. Prantzos 2008; Carigi et al. 2013; Spitoni et al. 2014, 2017; Spinelli et al. 2021; Scherf et al. 2024). In particular, Spitoni et al. (2014, 2017) applied detailed chemical evolution models, including also dust evolution, to identify habitable zones in the Milky Way for M and FGK stars and Andromeda. Galactic habitability has also been explored in relation to the bulge and the impact of a central AGN (e.g. Balbi et al. 2020; Ambrifi et al. 2022). On the other side, cosmological hydrodynamical simulations within the ΛCDM framework have studied the concept of habitability on larger scales by investigating the types of halos that can produce galactic structures conducive to the formation of habitable planets (Forgan et al. 2017; Gobat & Hong 2016; Zackrisson et al. 2016; Vukotić et al. 2016).

In most of the models discussed above, whether focused on Galactic chemical evolution or based on cosmological frameworks, the possibility of forming planetary systems with terrestrial planets, but without gas giants is considered. The presence of gas giants or hot Jupiters could, in principle, disrupt or even destroy Earth-like planets during their formation and evolution. However, there is no clear consensus in the literature regarding the precise effects of gas giants on terrestrial planet formation. In fact, it is well known that gas giant planets could exert a gravitational influence that promotes the accretion of planetesimals into larger bodies in the inner disk, while simultaneously mitigating excessive dynamical excitation that could otherwise impede their growth. Indeed, He & Weiss (2023) confirmed that the presence of outer gas giants increases the ‘gap complexity’ in planetary systems, thereby facilitating more efficient formation of terrestrial planets. Hence, in this study, we also aim to explore the possibility that giant planets may enhance the probability of terrestrial planet formation.

Moreover, for the first time, this paper considers the effect of stellar migration to develop a habitability map of our Galaxy within the framework of chemical evolution models, using the parametric approach proposed by Frankel et al. (2018) and later used in Palla et al. (2022). This is complementary to the work of Baba et al. (2024), who have recently focused on the orbital migration of the solar system. Stellar migration in the Galactic disc has been the subject of several investigations in recent years. Lynden-Bell & Kalnajs (1972) demonstrated that stellar bars and spiral arms cause stars to migrate away from their birth radii through the process known as radial migration. More recently, various explanations of the radial migration mechanism have been proposed, such as bar-induced migration (e.g. Solway et al. 2012; Sellwood & Carlberg 2014), and multiple bar and spiral modes (e.g. Minchev & Quillen 2006; Roškar et al. 2008). Minchev et al. (2013, 2015) and Kubryk et al. (2013, 2015) studied the effects of stellar migration on the chemical properties of galaxies, highlighting how the redistribution of stars influences the metallicity gradients and overall chemical evolution. Moreover, chemical abundance trends exhibit significant variations when considering stars’ guiding or current radii instead of their birth locations (Ratcliffe et al. 2023, 2024). In turn, the stellar migration process is also expected to play a crucial role in shaping the conditions necessary for planet formation and the emergence of habitable environments, thereby critically influencing the Galactic habitable zone (GHZ). In a related context, Baba et al. (2024) highlighted the dynamic nature of Galactic habitability, revealing how a star’s journey through the Milky Way profoundly influences its surrounding environment and, consequently, the potential for life. They introduced the concept of ‘galactic habitable orbits’, which considers the evolving structures of the Galaxy and their impact on stellar and planetary systems.

This paper is organised as follows: in Section 2, we outline the chemical evolution models used in this study and explain how stellar migration has been incorporated. Section 3 presents our GHZ model, while Section 4 details our results. Finally, in Section 5, we draw our concluding remarks. In Appendices A and B more details and results on the reference chemical evolution model and GHZ are reported. In addition, in Appendix C, we present the observational data used for comparison with the metallicity distribution predicted by our models.

2 Chemical evolution model for the Galactic disc

In Section 2.1, we describe the reference two-infall model adopted in this study. In Section 2.2, we indicate how stellar migration was implemented.

2.1 The reference multi-zone chemical evolution model

The adopted chemical evolution model is based on the two-infall model originally developed by Chiappini et al. (1997). We use the revised version by Molero et al. (2023, 2025), focusing exclusively on the thick and thin discs, without considering the evolution of the Galactic halo. The two-infall model assumes that the Galaxy forms as a result of two distinct gas accretion episodes. The first episode forms the chemical thick disc¹, while the second, delayed with respect to the first, leads to the formation of the thin disc. The infalling gas is assumed to be of primordial composition and the Galactic disc is modeled as independent concentric rings of width 2 kpc. The star formation rate (SFR) is the one parametrised by the Schmidt-Kennicutt law (Kennicutt 1998): $\psi (R,t)\propto \nu (R){\sigma }_{\mathrm{g}\mathrm{a}\mathrm{s}}(R,t{)}^k$, where σ_gas(R, t) is the surface gas density, k = 1.5 is the power-law index, and ν is the star formation efficiency, expressed in Gyr^-1. The parameter ν(R) varies with Galactocentric radius according to the prescription by Palla et al. (2020, as seen in their Table 3). The initial mass function (IMF) considered here is taken from Kroupa et al. (1993), which was specifically derived from star counts in the solar neighbourhood field. The gas accretion rate in the form of element i, in the two-infall model is: $$I(R,t{)}_{i,\mathrm{i}\mathrm{n}\mathrm{f}}=A(R)X_{i,\mathrm{i}\mathrm{n}\mathrm{f}}{e}^{-t/{\tau }_1}+\theta (t-t_{\mathrm{m}\mathrm{a}\mathrm{x}})B(R)X_{i,\mathrm{i}\mathrm{n}\mathrm{f}}{e}^{-(t-t_{\mathrm{m}\mathrm{a}\mathrm{x}})/{\tau }_2},$$(1)

where X_i,inf is the composition of the infalling gas, assumed to be primordial. The timescales, τ₁ and τ₂(R), describe the accretion duration for the thick and thin discs, respectively. We adopted τ₁ = 1 Gyr and allow τ₂(R) to vary with radius according to the inside-out formation scenario (e.g. Matteucci & Francois 1989; Chiappini et al. 2001): τ₂(R) = (1.033 R [kpc] - 1.267) [Gyr].

We recall here that the θ quantity in Eq. (1) is the Heaviside step function. The parameter t_max denotes the time of maximum accretion onto the thin disc, corresponding to the end of the thick disc formation and the start of the second infall episode. While earlier studies (e.g. Chiappini et al. 2001; Spitoni et al. 2009; Romano et al. 2010; Grisoni et al. 2018) have typically assumed t_max ~ 1 Gyr, recent works suggest a longer delay to reproduce stellar abundances and asteroseismic ages (Spitoni et al. 2019). Here, we adopted t_max ≃ 3.25 Gyr, following the prescriptions of Palla et al. (2020, 2022, 2024) and Molero et al. (2023); this value is also consistent with what was claimed in Nissen et al. (2020) and Spitoni et al. (2024) studies. In Eq. (1), the coefficients A(R) and B(R) are calibrated to reproduce the present-day surface mass densities of the thick and thin discs as functions of radius. Assuming exponential profiles, we adopt: $${\sigma }_{\mathrm{t}\mathrm{h}\mathrm{i}\mathrm{n}}(R)={\sigma }_{8,\mathrm{t}}{e}^{-(R-R_8)[\text{kpc}]/3.5}\text{ and }{\sigma }_{\mathrm{t}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{k}}(R)={\sigma }_{8,\mathrm{T}}{e}^{-(R-R_8)[\text{kpc}]/2.3},$$(2)

where σ_8,t and σ_8,t are the total surface mass densities at 8 kpc. We imposed the condition that σ_8,t + σ_8,τ = 47.1 M_⊙ pc^-2, as suggested by McKee et al. (2015) and assumed by Molero et al. (2025). Following Spitoni et al. (2020, 2022), Mackereth et al. (2017), and Molero et al. (2023), we supposed that the ratio between the two disc components at 8 kpc is ${\sigma }_{8,\mathrm{t}}/{\sigma }_{8,\mathrm{T}}\sim 4$. An important constraint for the chemical evolution model is the present-time stellar surface mass density in the solar vicinity. The proposed model predicts a value of 34.3 M pc^-2, in excellent agreement with the value of 33.4 ± 3 M pc^-2 suggested by McKee et al. (2015). Appendix A.1 presents the adopted nucleosynthesis prescriptions.

2.2 The implementation of stellar radial migration

As mentioned in the introduction, our analysis includes stellar radial migration in the multi-zone chemical evolution model of the Galactic disc components. We implemented it by including the approach used in Palla et al. (2022). This framework follows the parametric prescriptions proposed by Frankel et al. (2018), where migration is seen as a result of a diffusion process and treated in a parametric way (see also Frankel et al. 2020). In particular, the probability that is star would currently be located at a Galactocentric radius, R_f, given that it was born at R_i and has a specific age, τ, can be written as: $$\ln p(R_f|R_i,\tau )=\ln (c_3)-\frac{(R_f-R_i{)}^2}{2{\sigma }_{RM}\tau /10\mathrm{G}\mathrm{y}\mathrm{r}},$$(3)

where σ_RM is the radial migration strength expressed in kpc and c₃ a normalisation constant ensuring that stars do not migrate to negative radii (see Frankel et al. 2018 for details). For σ_RM we adopt as reference a value of 3.5 kpc, as used by Palla et al. (2022) and very similar to that found in Frankel et al. (2018). However, we also show results for σ_RM = 6 kpc to test the effect of an extreme case for stellar migration on the GHZ. As for the prescriptions from Palla et al. (2022), we assumed stellar migration for the two-infall model only after the onset of the second infall (~10 Gyr ago). This is also similar to the approach in Frankel et al. (2018), where migration was not considered for an ‘old disc’ (age > 8 Gyr in their case) component.

We recall that Frankel et al. (2018, 2020) applied their analysis to fit the observed age-metallicity distribution of low-α red clump stars within Galactocentric radii of 5-14 kpc, as derived from the Apache Point Observatory Galactic Evolution Experiment project (APOGEE DR12, Alam et al. 2015, DR14 Abolfathi et al. 2018). In our study, we extended this relation to a broader range, from 3 to 19 kpc. We acknowledge that this methodology may lead to an underestimation of stellar migration from the innermost regions of the Galaxy; however, in our parametric approach, this represents the best possible approximation. Moreover, the above-mentioned extreme migration scenario (σ_RM = 6) is also meant to overcome such a limitation, showing the effect of a larger number of migrators from the Galaxy inner regions. In this way, we would also expect to encompass the effects of a σ_RM varying in time. We also assumed that radial migration is only a “passive” tracer of chemical evolution, namely, the metal enrichment of the ISM in the solar vicinity is practically unaffected by the stars born in other regions.

The choice of the above-mentioned prescriptions to model the radial migration effects is driven by the analogue chemical evolution model framework adopted in Palla et al. (2022), which is already tested in the light of abundance data in the solar vicinity. In particular, Palla et al. (2022) used a chemical evolution model characterised by a two-infall scenario with a significant delay (3.25 Gyr) between the accretion episodes, reproducing the abundance diagrams as traced by the AMBRE:HARPS observational data sample (Santos-Peral et al. 2020, 2021) within the AMBRE project (de Laverny et al. 2013), assuming a radial migration strength of σ_RM = 3.5. This assures us to probe the effects of stellar migration on the GHZ in a reasonable and solid physical framework, already tested in the context of Galactic observations.

3 The GHZ model

The main goal of this work is to provide updated GHZ maps compared to those presented in Spitoni et al. (2014, 2017) by taking into account the effects of stellar migration and different prescriptions for gas giants. In the following, we present the main ingredients required to obtain a map of Galactic habitability as a function of the Galactic time and Galactocentric radius using the formalism introduced by Spitoni et al. (2014).

3.1 The GHZ model without stellar migration

We define P_E as the probability of forming a terrestrial planet around a star. Buchhave et al. (2012), analysing data from the Kepler mission, found that the occurrence rate of Earth-sized planets around single stars remains nearly independent of the host star’s metallicity for [M/H] values below 0.5 dex (see also Buchhave et al. 2018). This finding has been confirmed by other analyses based on Kepler and K2 (e.g. Petigura et al. 2018; Zink et al. 2023). A minimum metallicity threshold for planet formation was identified by Johnson & Li (2012) by comparing the timescales of dust grain growth and settling, with that of disk photoevaporation. Their theoretical study, which considered a wide range of circumstellar disc models and dust grain properties, indicated the existence of a critical metallicity above which planets can form as a function of the distance at which the planet orbits its host star. In particular, their calculations support the hypothesis that the first Earth-like planets likely formed from circumstellar discs with [Fe/H] ≃ -1.0 dex, but planets are unlikely to form at lower metallicities. To maintain consistency with these findings, Spitoni et al. (2017) assumed a probability of 0.4 for the formation of Earth-like planets around stars with [Fe/H] > -1 dex, while setting P_E = 0 for stars with lower metallicity. The specific choice of P_E = 0.4 aligns with the metallicity-integrated probability proposed by Lineweaver et al. (2004).

In Spitoni et al. (2014, 2017), the authors explored a scenario in which gaseous giant planets orbiting the same host star are able to disrupt terrestrial planets, particularly during their migration around the host star. Armitage (2003) highlighted the potential dangers posed by gas giant migration to the formation processes of Earth-like planets, suggesting that such planets are more likely to form in systems where massive giants have undergone a minimal migration. Similarly, Matsumura et al. (2013) investigated the orbital evolution of terrestrial planets in dynamically unstable systems with gas giants, showing that even Earth-like planets located far from the giants can be ejected.

We adopt the probability of formation of a gaseous giant planet as the function of the host star mass and iron abundance from Ghezzi et al. (2018): $$P_{GGP}(M_{\star },\text{[Fe/H}])={0.085}_{-0.010}^{+0.008}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{{1.05}_{-0.24}^{+0.28}}{10}^{{1.05}_{-0.17}^{+0.21}\text{[Fe/H]}}.$$(4)

In particular, this relation provides a good fit for the occurrence frequency of gas giant planets across a wide range of stellar host masses, including M, FGK, and retired A stars. To take into account of the quantity P_GGP in our GHZ model, we weighted the term (M_★/M_⊙)^1.05 on the Kroupa et al. (1993) IMF for the different mass ranges spanned by the considered stellar types². Assuming that the range of masses spanned by M-type stars is 0.1 ≤ (M_★/M_⊙) ≤ 0.45, for FGK stars is 0.45 ≤ (M_★/M_⊙) ≤ 1.40 and for retired A stars is 1.40 ≤ (M_★/M_⊙) ≤ 2.40, we find that the IMF-weighted P_GGP becomes $$<P_{GGP}(M_{\star },\text{[Fe/H]}){>}_{IMF}=\{\begin{array}{ll}M,& 0.085\times 0.205\times {10}^{1.05\text{[Fe/H],}}\\ FGK,& 0.085\times 0.723\times {10}^{1.05\text{[Fe/H],}}\\ retiredA,& 0.085\times 1.834\times {10}^{1.05\text{[Fe/H]}}.\end{array}$$(5)

In the upper panel of Fig. 1, we show the <P_GGP (M_★, [Fe/H])>_IMF quantity for M, FGK, and retired A stars. As a reference, in the same plot, we draw the Fischer & Valenti (2005) P_GG probability where no dependence on the stellar mass was considered. Hence, the probability of forming an Earth-like planet, P_E, as a function of the host star [Fe/H] abundance, assuming no formation of gas giants, is $$P_E(\text{[Fe/H}])=0.4\times (1-<P_{GGP}(M_{\star },\text{[Fe/H}]){>}_{IMF}),$$(6)

as shown in the lower panel of Fig. 1 by the solid lines labelled as ‘GG BAD’ for M, FGK, and retired A stars. It is worth commenting on the baseline frequency of Earth-like planets (excluding giant planets). We acknowledge that several studies in the literature have focussed on deriving this value in different cases. For instance, Bergsten et al. (2022) estimated a lower occurrence rate for Sun-like stars of around 0.15-0.2, while for M dwarfs, this value appears to be higher, namely, in the range of 0.5-0.8. However, to ensure consistency with previous works on the modeling of the GHZ (Lineweaver et al. 2004; Prantzos 2008; Spitoni et al. 2014, 2017), we adopted a uniform value of 0.4 for all stellar types. Having expressed the probability, P_E, we can then define the total number of stars formed at a certain time t and Galac-tocentric distance, R, hosting Earth-like planets with minimum habitability conditions, N_{★ mHC} (R, t), as $$N_{\star mHC}(R,t)=P_{GHZ}(R,t)\times N_{\star tot}(R,t),$$(7)

where N_★tot (R, t) is the total number of stars formed up to time t (and still alive at that evolutionary time) at the Galactocentric

distance, R, while P_GHZ (R, t) is the fraction of all stars having Earths (but no gas giants), that had survived SN explosions as a function of the Galactocentric distance and time. This latter quantity can be expressed as: $$P_{GHZ}(R,t)=\frac{{\int }_0^t\psi (R,t^{\prime })P_E(R,t^{\prime })P_{SN}(R,t^{\prime })d{t}^{\prime }}{{\int }_0^t\psi (R,t^{\prime })d{t}^{\prime }},$$(8)

and must be interpreted as the relative probability of having minimum habitability conditions around one star at a given position (see Prantzos 2008). In Eq. (8), we have ψ(R, t′) as the star formation rate (SFR) at a time, t′, and Galactocentric distance, R; P_SN(R, t′) is the probability of surviving a SN explosion (other catastrophic events, such as nearby nova outbursts or compact object mergers, are not taken into account because of the rarity of their occurrence). Due to uncertainties on the actual impact of SNe on life destruction, we consider two possible scenarios, referred to as case 1 and case 2, as presented in Spitoni et al. (2014). In these scenarios, SN destruction is assumed to be effective if the total SN rate (i.e. CC + Type Ia) at any given time and location exceeds the average SN rate in the solar vicinity over the past 4.5 Gyr of the Milky Way’s life (<R_{S N,⊙}>) for case 1, and if it is larger than 2 × < R_{S N,⊙}> for case 2. Therefore, the associated probability, P_SN(R, t), is $$P_{SN}(R,t)=\{\begin{array}{ll}0,& \text{if } R_{SN}(R,t) > <R_{SN,\odot }> \text{ (case 1),}\\ & \vee R_{SN}(R,t) >2\times <R_{SN,\odot }>\text{ (case 2),}\\ 1,& \text{otherwise.}\end{array}$$(9)

For <R_SN,⊙>, we adopted the value of 0.023 Gyr^-1 pc^-2 as reported in the results of the reference chemical evolution model (presented in Section 2.1). To understand how the average <R_SN,⊙> compares to the SN rates at different Galactocentric radii, the left panel of Fig. 2 presents the case 1 and case 2 thresholds alongside the predicted time evolution of total SN rates. It is important to stress that in this work we are exploring the possibility of Earth-like planets located in regions safe from SN damage. However, the number of life-hosting planets among them may be significantly reduced when considering additional key requirements for life development; specifically, N₂-O₂-dominated atmospheres with constrained CO₂ levels, as recently discussed in Scherf et al. (2024). In addition, it is worth noting that a detailed analysis of the effects of SN explosions on the GHZ is beyond the scope of this paper. Our primary aim is to highlight the impact of stellar migration on previously presented models. A more precise and comprehensive treatment of SN-induced disruptions to habitability will be addressed in future work.

Fig. 1 Probability of forming planets as a function of metallicity and stellar type. Upper panel: probabilities of forming gas giant planets as a function of [Fe/H] for M, FGK and retired A proposed by Ghezzi et al. (2018) are shown with the black, blue and red lines, respectively. The probabilities have been IMF weighted according to the procedure described in Section 3.1. Fischer & Valenti (2005) relation is also shown with the dashed grey line. Lower panel: probability, PE, of forming Earth-like planets and not gas giants (GG BAD case) as a function of [Fe/H], and indicated in Eq. (6) for different stellar spectral types, shown as solid lines. In contrast, the dashed lines illustrate the new scenario discussed in Section 3.3 and called the GG GOOD case, which assumes that the presence of gas giants enhances the likelihood of terrestrial planet formation. All the cases are computed using the gas giant probabilities of Ghezzi et al. (2018).

Fig. 2 Predictions of our multi-zone chemical evolution model presented in Section 2.1 as a function of the evolutionary time t and the Galacto-centric distance. Left panel: evolution of the total (CC+Type Ia) SN rates. The red and green dashed horizontal lines indicate <RS N, ⊙> (case 1) and 2 × <RSN,⊙> (case 2) values, respectively, representing the minimum SN rate thresholds tested for SN explosion-induced destruction effects in this study. Right panel : predicted age-metallicity relation. The red-shaded region marks [Fe/H] < -1 dex where, as suggested by Johnson & Li (2012), conditions are considered unfavourable for the formation of terrestrial planets.

3.2 GHZ formalism in the presence of the stellar migration

In the presence of the stellar migration, the definition of P_GHZ and the formalism to compute N_★,mHC must be revisited to account for the increased complexity added to the system. Here, we adopt the following expression: $$N_{\star mHC}(R_f,t)=\sum _{R_i}(P_{GHZ}(R_i\to R_f,t)\times N_{\star }(R_i\to R_f,t)),$$(10)

where P_GHZ(R_i → R_f, t) is defined as: $$P_{GHZ}(R_i\to R_f,t)=\frac{{\int }_0^t\psi (R_i,t^{\prime })P_E(R_i,t^{\prime })P_{SN}(R_f,t^{\prime }) d{t}^{\prime }}{{\int }_0^t\psi (R_i,t^{\prime }) d{t}^{\prime }}.$$(11)

The quantity P_GHZ(R_i → R_f, t) expresses the fraction of all stars born at the Galactocentric distance Ri that have migrated to R_f by the time t and could potentially host a terrestrial planet. In Eq. (11), the abundance ratio [Fe/H] used to compute the probability P_E (see Eq. (6)), is determined by the value predicted by the chemical evolution model at the star birth radius, R_i. This is because the stellar metallicity is inherited by the ISM abundance at stellar birthplace and birthtime (e.g. Trimble 1996; Pagel 2006; Matteucci 2021). However, the SN explosion damage effect, represented by P_SN, must be evaluated at the star final location, R_f, as SN explosions happening at the time, t, are impacting planetary systems in their surroundings, without any relations with the stellar planetary host system birthplace. We emphasise that for stars formed in situ (R_i = R_f), the expressions given in Eqs. (10) and (11) are equivalent to those in Eqs. (7) and (8) presented in Section 3.1 for the reference model without stellar migration. In summary, the total number of stars hosting habitable terrestrial planets that have migrated from R_i to R_f by the evolutionary time, t, is obtained by multiplying P_GHZ(R_i → R_f, t) by the number of migrating stars, N_★(R_i → R_f, t). To account for all stars that have reached R_f with the potential to support life, we sum over the different birth radii, R_i, as expressed in Eq. (10).

3.3 Gas giants as catalysts for terrestrial planet formation

In previous GHZ models gas giants have been considered an important hazard for the creation and evolution of terrestrial planets. For instance, Walsh et al. (2011) discussed how Jupiter’s migration sculpted the inner solar system, limiting material available for terrestrial planets. On the other hand, gas giants clear material in their orbital vicinity, creating gaps in the protoplanetary disk. This can enhance planet formation by concentrating planetesimals and embryos into specific regions, as described in the ‘Nice Model’ and similar frameworks (e.g. Tsiganis et al. 2005). For instance, simulations have shown that this process can help stabilise inner regions for terrestrial planet growth (Morbidelli et al. 2012). The gravitational influence of gas giants promotes the accretion of planetesimals into larger bodies in the inner disk while also mitigating excessive dynamical excitation that might hinder their formation. Zhu & Wu (2018) found a positive correlation between the presence of an Earth-like planet and the existence of a cold giant within the same system. He & Weiss (2023) confirmed that the presence of outer giants increases the ‘gap complexity’, favouring more efficient terrestrial planet formation. Gas giants can also scatter icy bodies from the outer disk inwards, delivering essential volatiles such as water to form terrestrial planets. This is supported by isotopic analysis of Earth’s water, which is aligned with an outer Solar System origin, corroborated by O’Brien et al. (2018).

All these findings highlight the dual role of gas giants as either facilitators or potential disruptors of Earth-like planet formation, depending on their location, mass, and dynamical evolution. For this reason, we also explored the opposite scenario to the one proposed in Spitoni et al. (2014, 2017); namely, considering the probability of forming an Earth-like planet, P_E, in the extreme case where the formation is pumped up by the presence of gas giant planets. The dashed lines in Fig. 1 represent this extreme case for different stellar spectral types, which can be expressed as $$P_E(\text{[Fe/H}])=0.4\times (1+<P_{GGP}(M_{\star },\text{[Fe/H}]){>}_{IMF}).$$(12)

Throughout this article, we refer to this case as GG GOOD. In Fig. 1, we note that differences between GG BAD and GG GOOD scenarios are more pronounced for retired A stars, particularly at super-solar values. For the first time in GHZ map modelling, we tested this assumption using chemical evolution models, showing the results in Section 4.3. In Table 1, we summarise all the GHZ models we have considered in this study, obtained by varying: the probability of finding terrestrial planets as a function of the metallicity in the presence of gas giants (GG GOOD or GG BAD), the SN damage scenario (case 1 or case 2), and the stellar migration parameters.

4 Results

In this section, we report our main results based on the model prescriptions introduced in Sections 2 and 3. We present and discuss the results for the ‘classical’ chemical evolution model, which does not account for stellar migration, in Section 4.1. Then, in Section 4.2, we report the main effects of stellar migration on the computation of the GHZ. Finally, GHZ maps obtained under the assumption that gas giant planets favour the formation of Earth-like planets are shown in Section 4.3.

4.1 Reference model results without stellar migration

In Fig. 2, we show the temporal evolution of total SN rates and the age-metallicity relation of the multizone model presented in Section 2.1. Concerning the SN rates, we also indicate the <R_SN,⊙> and 2 × < R_SN,⊙> quantities defined in Section 3. For most of the thin disc evolution in the solar vicinity (the annular region centred at 8 kpc), the calculated SN rates exceed the average threshold, <R_SN,⊙>. However, when the higher limit of 2 × < R_SN,⊙> is imposed, the evolution remains entirely below this threshold.

The predicted age-metallicity relation³ shown in the right panel of Fig. 2 highlights the fact that during the first Gyr of Galactic chemical evolution in the thick disc phase, the ISM undergoes a rapid metal enrichment. Notably, at an evolutionary time of 0.25 Gyr, the [Fe/H] = -1 dex threshold, identified by Johnson & Li (2012) as the minimum metallicity required for the formation of terrestrial planets, is already reached. It is worth noting that the dilution phase characterising the early evolution of the thin disc (at an evolutionary time of ~3.25 Gyr) leads to a decrease in metallicity across all Galactic annular regions, with the effect being more pronounced at 18 kpc. However, this dilution phase does not reduce the [Fe/H] abundance below the -1 dex threshold, as shown in the lower panel of Fig. 2. In conclusion, throughout the evolution of the Galactic disc at all radii (except for the first 0.25 Gyr) the metallicity remains above the threshold identified by Johnson & Li (2012).

In Fig. 3, we show for the different Galactocentric distances the evolution of the probability P_E for stars with different spectral types introduced by Eq. (6), where gas giants are assumed to be hazards for life. It is worth noting that in the 3D space formed by P_E, [Fe/H] and evolutionary time, the transition between the thick and thin disc phases emerges clearly. In fact, the stronger metal dilution in outer Galactic regions (see also Fig. 2), together with the inside-out formation growth of the disc, implies a small effect of gas giants on the probability of finding Earth-like planets P_E. Conversely, the inner regions achieve lower values of P_E due to the higher metallicities reached during the thin disc evolution phase. This effect is more evident for retired A stars.

Keeping in mind the evolution of the probability, P_E, presented in Fig. 3, we compare in Fig. 4 the number of stars hosting habitable Earth-like planets but not gas giant planets (N_★,mHC) as a function of Galactic time and the Galac-tocentric distance, as predicted by model 1 for M, FGK and retired A stars, respectively. Consistently with the findings of Spitoni et al. (2014, 2017), the distribution of N_★mHC over Galactic time and across Galactocentric distances reveals lower values in the innermost regions, primarily due to the high rates of SN events. Similarly, in the outskirts, the number remains low because of the reduced number of stars predicted by the model under the inside-out formation scenario. Model 1, which assumes the case 2 SN damage scenario, shows the peak of Galactic habitability at 8 kpc for all spectral types. The numbers of M and FGK stars suitable for habitability (N_★,mHC) are compatible with Spitoni et al. (2017) predictions. The main difference lies in the ratio between M and FGK stars, which is higher in this study due to the different adopted IMF (Kroupa et al. 1993 versus Scalo 1986). Regarding retired A-type stars, we emphasise that the lifetime of the average stellar mass (weighted by the IMF (i.e. 1.78 M_⊙) is 2.4 Gyr. This is the reason why in Fig. 4 and in Fig. B.1 of Appendix B, the distribution of stars hosting habitable terrestrial planets closely follows the star formation history, revealing two distinct clumps associated with the peak gas infall rates during the thick and thin disk phases. We emphasise that due to the short lifetimes of retired A-type stars, the development of complex life forms is unlikely on these objects. In Fig. 5, we show the evolution of the number of FGK stars hosting habitable Earth-like planets but not gas giant planets (N_★,mHC) as a function of Galactic time and the Galactocentric distance, predicted by model 2. In contrast to the results for model 1, model 2, which incorporates the case 1 SN effect, predicts peak habitability at 10 kpc. This difference arises because (as noted earlier in this paper) the average SN rate in the solar vicinity over the past 4.5 Gyr exceeds the predicted rates during the thin-disc phase for most of the Galaxy’s evolution at R = 8 kpc. Importantly, at 4 kpc, no stars can host terrestrial planets capable of supporting life because of the extremely high SN rates over the course of Galactic history. Different thresholds in SN rate also impact the absolute values in the number of stars hosting terrestrial planets with conditions favourable for the development of life (N_★,mHC). In model 1, the present-day value of N_★,mHC in the solar vicinity (annular region centred at 8 kpc from the Galactic centre) is 1.12 × 10⁹⁴. At the same Galactocentric distance and time, model 2 predicts 3.91 × 10⁸ (N_★,mHC), as highlighted in Fig. 5.

In Fig. B.2 of Appendix B, we present the results for a ‘modified model 1,’ which accounts for the time necessary for a sustained increase in atmospheric O₂ to significant levels. On Earth, this process needed approximately 2.5 Gyr (Lyons et al. 2014). Specifically, we considered delays of 1 Gyr and 3 Gyr in our analysis. At a distance of 8 kpc, we observe that the present-day total number of habitable stars, N_★,mHC, decreases as the delay increases compared to the reference model 1. For a 1 Gyr delay, N_★,mHC is 1.08 × 10⁹, representing a 3.57% decrease from the reference model 1. For a 3 Gyr delay, N_★,mHC decreases further to 9.48 × 10⁸, corresponding to a 21.00% reduction.

Fig. 3 Evolution of the probability, PE, of forming Earth-like planets and not gas giants (the ‘GG BAD’ case illustrated in Eq. (6) and in Fig. 1), as predicted by our reference multi-zone chemical evolution model presented in Section 2.1 in the 3D space formed by the PE, [Fe/H] and evolutionary time t for different stellar types.

Fig. 4 Total number of different spectral stellar types hosting Earthlike planets with minimum habitability conditions (N★,mHC in Eq. (8)) as a function of Galactic distance and time predicted by the chemical evolution model without stellar migration (model 1; see Section 3.1 and Table 1). The values of (N★,mHC) are computed within concentric rings of 2 kpc width. Results for M stars are drawn in the upper panel. FGK and retired A are shown in the middle and the lower ones, respectively.

Fig. 5 Total number of FGK stars with minimum habitability conditions (N★,mHC in Eq. (8)) as a function of Galactocentric distance and Galactic time predicted by the model 2; i.e. case 1 SN damage scenario and no stellar migration (see Section 3.1 and Table 1).

4.2 Model results with stellar migration

In this section, we present our findings in the presence of stellar migration. In Fig. 6, we show the temporal evolution of the total number of FGK stars with habitable terrestrial planets (N_★,mHC) in the solar vicinity, highlighting the contributions from migrators originating in different Galactic regions. In the presence of stellar migration, we assume case 2 for the SN destruction scenario, adopting σ_RM = 3.5 and σ_RM = 6 (model 3 and model 4, see Table 1). We stress that in our models, stellar migration is considered to begin with the onset of the thin disc phase, starting 3.25 Gyr after the beginning of Galactic evolution. Aside from the stars formed in situ, namely, those born at 8 kpc, the majority of FGK stars with habitable planets originate from 4, 6, and 10 kpc. Looking in more detail at model 3, at the present time, the total number of predicted stars with habitable terrestrial planets (currently residing in the solar vicinity) is N_★,mHC = 1.14 × 10⁹. Of these, 43.78% were formed in situ, 17.16% originated from 4 kpc, 24.53% from 6 kpc, and 11.57% from 10 kpc.

Similarly, in model 4, the total number of stars hosting habitable planets is 1.01 × 10⁹ N_★,mHC. Of these, 38.82% were formed in situ, 23.21% originated from 4 kpc, 22.51% from 6 kpc, and 10.27% from 10 kpc. Concerning the total number of stars hosting life, we notice that the computed present-day value at 8 kpc is similar to that predicted by the model without stellar migration examined in the previous section. This result can be explained by the fact that in the solar vicinity, we are able to reach a ‘balance’ between stars leaving the solar vicinity and stars migrating towards this Galactic region (see e.g. Palla et al. 2022). However, as we show later in the text, such a balance is not present ubiquitously across the Galactic disc.

In the left panel of Fig. 7, we present the cumulative metal-licity distribution function (CMDF) predicted by model 3 for M and FGK stars possibly hosting life (N_★,mHC) in the solar vicinity. Our results are compared with the observed distribution of a sample of planet-hosting stars in the solar vicinity, comprising 124 stars (106 FGK-type and 18 M-type) that host planets with masses of M_P < 30 M_⊕. This sample is divided into local (born within 1 kpc of the Sun’s location), inner (originating 1-2 kpc closer than the Sun to the Galactic centre), and outer (located 1-2 kpc beyond the Sun) populations. This separation is based on the orbital properties of each star, from which we can extrapolate their origin in different regions of the Galactic disk. Further details on the stellar sample are provided in the Appendix C. To ensure consistency with the observed local, inner, and outer populations (see Fig. C.1 in Appendix C), we consider only predicted stars born at 6, 8, and 10 kpc in this analysis. The overall distribution of these populations is well reproduced by model 3 predictions, although the individual ones corresponding to stars originating from 6 and 10 kpc show that the model predicts a higher number of migrators. The data exhibit a higher contribution from the local sample (76.5%) compared to model 3 (54.2%), indicating that a lower migration rate is needed to reproduce the data. However, we must keep in mind that the sample used for comparison with the model accounts only for part of stellar migration. Indeed, data can account only for the ‘blurring’ component (increase in the orbital eccentricity of a star, while maintaining the mean radius), while no direct observables are available for ‘churning’ (change in mean radius and angular momentum of the orbit). Moreover, it is worth stressing that the sample is affected by observational biases, not representing all disc stars with minimal habitability conditions. In the right panel of Fig. 7, we explore a scenario with a weaker radial migration, adopting σ_RM = 1 while keeping all other model parameters identical to those of model 3 (‘model weak’). In this case, the agreement with the observed sample across the different populations is nearly perfect. The reduction of the migration strength of a factor larger than 3 is actually in line with what was found in previous works (e.g. Frankel et al. 2020), where blurring accounts for a minor contribution in the radial migration process. However, we emphasise once again that a direct comparison with the observational sample is not entirely straightforward, due to internal biases and unknown selection effects in the latter.

In Fig. 8, we provide the full map of habitability in the 3D plot as functions of the Galactocentric distance and evolutionary times comparing models with and without migrations (model 1, model 3, and model 4) for FGK stars. As discussed above, at 8 kpc the total number of stars is not dramatically affected by the migration. However, it is possible to note that the distribution of stars with habitability conditions in their hosting planets is wider towards external Galactocentric distances. For instance, at 12 kpc, the present-day number of stars hosting habitable terrestrial planets increases from N_★,mHC = 4.18 × 10⁸ (model 1, without migration) to 4.92 × 10⁸ (model 3, migration with σ_RM = 3.5 kpc) and 5.78 × 10⁸ (model 4, migration with σ_RM = 6 kpc).

For this reason, in Fig. 9, we present the ratio of N_★,mHC values predicted by models with stellar migration to those without, as a function of evolutionary time and Galactocentric distance. In this way, it is possible to better visualise the impact of stellar migration on the redistribution of stars hosting habitable planets. This is particularly important in the outer Galactic regions, where the absolute value of the predicted number of stars is significantly smaller compared to the inner regions. In particular, we see that model 4, which is characterised by σ_RM = 6, shows a substantial increase in the number of stars hosting habitable planets in the outer Galactic regions. This occurs because stellar migration enables stars to move towards more external parts, characterised by lower star formation activity, but reduced SN rates (see Fig. 2). This creates more favourable conditions for life to develop. For example, in model 4, at a Galactocentric distance of 18 kpc, the number of habitable planets increases by a factor of approximately 4.9 at around t = 6.1 Gyr of Galactic evolution for FGK stars⁵. It is important to note that at 4 kpc, ratio values are available only for the last 2 Gyr of Galactic evolution, during which the model without stellar migration exhibits a SN rate below the threshold defined by the case 2 SN destruction prescriptions (see Fig. 2). The increased number of stars hosting habitable planets in the outer Galaxy due to stellar migration widens the outer boundaries of the GHZ. Such a result is particularly relevant considering that the abundances of organic molecules found in outer Galaxy star-forming regions are comparable to those measured near the Sun (e.g. Bernal et al. 2021; Fontani et al. 2022, 2024), despite the sub-Solar metal-licity. Some of these molecules may have been the seeds of organic compounds that favoured the emergence of life in the Solar System and elsewhere in the Universe. Therefore, both findings go in the same direction that the conditions in the outer Galaxy are likely more favourable to host habitable planets than we previously thought.

As discussed in Section 2.1, our analysis adopts a delay of t_max = 3.25 Gyr between the two gas infall episodes to match the APOGEE observational constraints (Palla et al. 2020). To assess the sensitivity of our results to this parameter, we also explored models with t_max = 1 and 5 Gyr. Relative to model 4 (see the lower panel of Fig. 8), we find that at the Galactocentric distance of 12 kpc the number of FGK stars with minimum habitability conditions (N_★,mHC) increases by approximately 6.4% for t_max = 1 Gyr, and decreases by about 8.0% for t_max = 5 Gyr. This behaviour reflects our assumption, following the prescriptions of Frankel et al. (2018), that stellar migration occurs only during the thin-disc phase, which is longer in the t_max = 1 Gyr model and shorter in the t_max = 5 Gyr case.

Fig. 6 Temporal evolution of the total number of FGK stars hosting minimum conditions for life (N★,mHC) in the solar vicinity (black dashed lines). Both panels show the contribution from stars born in situ within the solar annulus (Ri = Rf = 8 kpc, dark blue lines) as well as those migrating from different Galactic regions (Ri ≠ Rf). We consider the case 2 SN destruction scenario in the presence of stellar migration, with radial migration strength of σRM = 3.5 kpc (model 3, upper panel) and σRM = 6 kpc (model 4, lower panel), respectively.

Fig. 7 CMDFs for FGK and M stars hosting habitable terrestrial planets (N★,mHC) in the solar vicinity predicted by models with stellar migration (solid lines). Model results are compared with the A&S data sample described in Appendix C for stars hosting only low mass planets (Mp < 30M⊕, dashed lines). We display the contribution from predicted stars that were born in situ within the solar annulus (Ri = Rf = 8 kpc, solid blue lines), as well as those that migrated from 6 and 10 kpc (solid red and grey line, respectively) to ensure consistency with the observed local, inner, and outer populations. Model 3 (σRM=3.5) results are reported in the left panel. In the right one, we show the ‘model weak’ approach, adopting a weaker radial migration (σRM=1), while keeping all other model parameters identical to those of model 3.

Fig. 8 Total number of FGK stars having Earths (N★,mHC) as a function of the Galactocentric distance and time considering the case 2 SN destruction scenario. Upper panel: results from the reference model without including stellar migration (model 1). Middle panel: results from the model incorporating stellar migration as described in Section 2.2, with the radial migration strength fixed at σRM = 3.5 kpc (model 3). Lower panel: same as the middle panel but with a radial migration strength of σRM = 6 kpc (model 4). In each panel, we highlight with a green cross the maximum number of FGK stars hosting habitable terrestrial planets predicted at 12 kpc (grey point in the 2D projection).

Fig. 9 Ratios of the number of FGK stars hosting Earth-like planets (N★,mHC) in models with stellar migration to the reference model without migration, as a function of Galactocentric distance and Galactic time. All considered models adopt case 2 for SN damage. Upper panel: stellar migration model with σRM = 3.5 (ratio of model 3 to model 1 in Table 1). Lower panel : stellar migration model with σRM = 6 (ratio of model 4 to model 1 in Table 1).

4.3 Results with gas giants as catalysts for terrestrial planets

As anticipated in Section 3.3, another main purpose of this work is to test the scenario where gas giants act as catalysts for terrestrial planet formation using the new probability, P_E GG GOOD, introduced in Eq. (12) and also given in the lower panel of Fig. 1. Figure 10 shows the ratios of the predicted number of M, FGK and retired A stars hosting habitable Earth-like planets (N_★,mHC) under two distinct P_E scenarios: GG GOOD and GG BAD, as defined in Fig. 1. Specifically, the following cases are analysed: i) no migration (ratio between model 7 and model 1 from Table 1); ii) stellar migration with σ_RM = 3.5 (the ratio between model 5 and model 3 from Table 1). As noted in Section 3.3, the two P_E probabilities, GG GOOD and GG BAD, exhibit significant differences only at super-solar metallicities. Consequently, we expect disparities primarily in the innermost regions of the Galactic disc, where higher metallicities are predicted. For all the stellar spectral types considered, the differences between the GHZ maps in the GG GOOD and GG BAD scenarios, as shown in Fig. 10, are primarily concentrated in the annular region centred at 4 kpc. As expected from Fig. 1, the highest GG GOOD and GG BAD ratio is found for retired A stars, reaching 3.3 in the model without stellar migration. For FGK stars, the maximum ratio is 1.53, while for M stars - where P_E remains nearly constant across all metallicity ranges in both scenarios - it is 1.29. In Fig. B.1 of Appendix B, we show the total number of retired A stars having Earths (N_★,mHC) as a function of the Galac-tocentric distance considering different gas giant effects. It is possible to appreciate that at 8 kpc the number of stars within the GG GOOD scenario increased by a factor of 1.23. However, these factors decrease when stellar migration is included, as illustrated in the lower panels of Fig. 10. This reduction is a natural consequence of the net effect of stellar migration, which (on average) redistributes stars from the innermost regions towards outer Galactic areas. In general, the GG GOOD and GG BAD scenarios yield minimal differences in numbers of the predicted stars hosting habitable planets (N_★,mHC). In fact, regardless of the presence of stellar migration, the enhancement factor is practically negligible in the regions where the peak of habitability occurs (i.e. in the solar vicinity for the case 2 SN damage model).

Fig. 10 Ratios between the predicted number of stars having Earths (N★,mHC) of the two different prescriptions for PE as in Fig. 1: GG GOOD and GG BAD cases, respectively. Results for M, FGK, and retired A are reported in the left, middle, and right panels, respectively. Ratios are computed as functions of the Galactocentric distance and Galactic time and all the displayed models consider the case 2 for the SN damage. Upper panels: reference models without stellar migration (ratio between model 7 and model 1 in Table 1). Lower panels: models with stellar migration and σRM = 3.5 (ratio between model 5 and model 3 in Table 1). In each panel the maximum value of the ratio is highlighted with an orange cross.

5 Conclusions

In this paper, we present new GHZ maps for the Milky Way disc constructed by means of detailed and well-tested chemical evolution models (see Molero et al. 2023, 2025). We have taken into account, for the first time, the effects of stellar migration, following the parametric approach suggested by Frankel et al. (2018) and Palla et al. (2022). We explored different scenarios for the effect of the presence of gas giant planets on the formation of terrestrials, considering both positive and negative effects. Our main results are summarised below:

• The Galactic habitability maps are highly sensitive to the adopted threshold for SN rates required for destruction to be effective. Model 1, which assumes the highest SN threshold (case 2), shows the peak of Galactic habitability at 8 kpc. In contrast, model 2, which adopts the case 1 SN effect, predicts the peak shifted at 10 kpc;

• At the present time, in the solar vicinity, the total number of predicted stars with habitable terrestrial planets, N_★,mHC , is similar in models with or without stellar migration. However, stellar migration leads to a redistribution of stars and a large fraction of them originated in different Galactic regions. For instance, in model 3 (where the migration strength is fixed at the value of σ_RM = 3.5 kpc, as suggested by chemical evolution studies) 17.16% of FGK stars hosting habitable planets are born at 4 kpc, 24.53% at 6 kpc, 11.57% at 10 kpc, and 43.78% are formed in situ;

• Stellar migration has a larger impact in the outskirts of the Galactic disc. In the annular region centred at 18 kpc, the number N_★,mHC of stars hosting habitable planets is increased by a factor of ~5 compared to the reference model without stellar migration when the extreme migration case is considered (model 4);

• The hypothesis that gas giant planets facilitate the formation of terrestrial planets has the most pronounced effect in the inner Galactic disc (in the annular region centred at 4 kpc from the Galactic centre). However, we find that this increased probability is also mitigated by stellar migration. In particular, at the present time, when stellar migration is taken into account, the number of FGK stars hosting habitable terrestrial planets in the inner ring centred at 4 kpc is approximately 1.4 times higher than in scenarios where gas giants are assumed to hinder the formation and evolution of Earth-like planets. Without stellar migration, this factor increases to 1.5. Even larger ratios are predicted for terrestrial planets orbiting retired A stars, reaching 2.8 in models with stellar migration and 3.3 in models without it. In general, the two scenarios (i.e. GG BAD and GG GOOD cases in Table 1) yield to minimal differences in terms of absolute numbers of the predicted stars hosting habitable planets (N_★,mHC). In fact, regardless of the presence of stellar migration, the enhancement factor is practically negligible in the regions where the peak of habitability occurs forM and FGK stars (i.e. in the solar vicinity).

In conclusion, here we have significantly expanded the exploration of the parameter space defining the GHZ, compared to previous analyses present in literature (e.g. Lineweaver et al. 2004; Prantzos 2008; Spitoni et al. 2017). Our findings are particularly relevant in the context of upcoming space missions, such as the ESA PLAnetary Transits and Oscillations of Stars (PLATO; Rauer et al. 2024), the ESA Ariel space mission (Tinetti et al. 2018, 2022), and Large Interferometer For Exoplanets (LIFE, Quanz et al. 2022). These missions will deliver unprecedented data on planetary properties, orbital architectures, and atmospheric compositions. To fully interpret this information, it will be essential to place exoplanets within a broader Galactic framework that accounts for stellar chemical composition, formation environment, and radial migration. This work contributes to the foundational understanding of planetary formation and habitability from circumstellar to Galactic scales, forming part of the necessary groundwork to interpret the influx of data expected from upcoming discoveries.

Acknowledgements

The authors thank the referee, R. Heller, for helpful suggestions that improved the quality of the manuscript. The analysis and work presented in this paper were developed during and after discussions at the “Molecules and planets in the outer Galaxy: is there a boundary of the Galactic Habitable Zone?” workshop, funded by INAF. E.S., F.M. and G.C. thank I.N.A.F. for the 1.05.24.07.02 Mini Grant - LEGARE “Linking the chemical Evolution of Galactic discs AcRoss diversE scales: from the thin disc to the nuclear stellar disc” (PI E. Spitoni). E.S., L.M., M.T. and G.C. thank I.N.A.F. for the 1.05.23.01.09 Large Grant - Beyond metallicity: Exploiting the full POtential of CHemical elements (EPOCH) (ref. Laura Magrini). M.P. and D.R. acknowledge financial support from the project “LEGO - Reconstructing the building blocks of the Galaxy by chemical tagging” granted by the Italian MUR through contract PRIN2022LLP8TK_001 (PI: A. Mucciarelli). D.R. acknowledges I.N.A.F. for financial support through the 1.05.12.06.08 Theory Grant - An in-depth theoretical study of CNO element evolution in galaxies (ref. Donatella Romano). F.M. thanks I.N.A.F. for the 1.05.12.06.05 Theory Grant - Galactic archaeology with radioactive and stable nuclei. F.M. thanks also for support from Project PRIN MUR 2022 (code 2022ARWP9C) “Early Formation and Evolution of Bulge and HalO (EFEBHO)” (PI: M. Marconi), M.M. thanks for the support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -- Project-ID 279384907 - SFB 1245, the State of Hessen within the Research Cluster ELEMENTS (Project ID 500/10.006). L.M. and G.C. thank INAF for the support with the Mini-Grants Checs (1.05.23.04.02), and the financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 104 published on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union - NextGenerationEU - Project ‘Cosmic POT’ Grant Assignment Decree No. 2022X4TM3H (PI: L. Magrini) by the Italian Ministry of the University and Research (MUR). C.D. acknowledges financial support from the INAF initiative “IAF Astronomy Fellowships in Italy”, grant name GExoLife. This work used the Ariel Stellar Catalogue developed by the Stellar Characterisation WG in preparation of the ESA Ariel space mission. The catalogue can be found at bit.ly/ArielStellarCatalogue. Supported by Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union - NextGenerationEU - and National Recovery and Resilience Plan (NRRP) - Mission 4 Component 2 within the activities of Spoke 3 (Astrophysics and Cosmos Observations).

Appendix A Chemical evolution model

A.1 Nucleosynthesis prescriptions

The nucleosynthesis prescriptions for the reference chemical evolution model adopted in this work and described in Section 2.1 follow those of the MWG11 model from Romano et al. (2019). For low- and intermediate-mass stars, we adopt the stellar yields from Ventura et al. (2013). Our models also account for the nucleosynthetic outcomes of binary systems that give rise to Type Ia SNe, assuming the single-degenerate scenario for progenitors as described by Matteucci & Recchi (2001) and references therein, with nucleosynthetic yields from Iwamoto et al. (1999). The yield for massive stars in the presence of stellar rotation are the ones of Limongi & Chieffi (2018). A variable initial rotational velocity for massive stars is considered, favouring high rotational speeds at low metallicities while assuming no rotation at solar metallicity (see also Prantzos et al. 2018; Rizzuti et al. 2019; Molero et al. 2024). Specifically, we adopt v_rot = 300 km s^-1 for [Fe/H]< -1, v_rot = 0 for [Fe/H] ≥ -1 dex.

Fig. A.1 Predicted time evolution of the SFR of our multi-zone chemical evolution model presented in Section 4.1 as a function of the evolutionary time, t, and the Galactocentric distance.

A.2 Other model predictions

In Fig. A.1, we show the predicted time evolution of the SFR of our multi-zone chemical evolution model presented in Section 2.1. In Fig. A.2, we compare the total number of stars predicted by the reference model without stellar migration to those where radial migration is taken into account. For the latter, we fix the migration strength σ_RM to two values: 3.5 kpc, as used in Palla et al. (2022) and Frankel et al. (2018), and 6 kpc, representing an extreme case of stellar radial migration. As expected, the stellar migration has the effect of increasing the number of stars in the outer parts of the Galactic disc. For instance, at a Galactocentric distance of 12 kpc, the present day total number of stars predicted by the reference model without migration is 3.71 × 10⁹. When migration is considered with σ_RM = 3.5, the number increases to 4.41 × 10⁹, and in the extreme case with σ_RM = 6, the total reaches 5.25 × 10⁹.

Fig. A.2 Total number of predicted stars (N★) as a function of the Galactocentric distance and the Galactic time. Upper panel: Results from the reference model without including stellar migration. Middle panel: Results from the model incorporating stellar migration, as described in Section 2.2, with the radial migration strength fixed at σRM = 3.5. Lower panel : Same as the middle panel but with a radial migration strength of σRM = 6. In each panel, we highlight with a green cross the maximum number of stars computed at 12 kpc.

Appendix B Other GHZ results

In this appendix, we provide more results on the GHZ without stellar migration. In Fig. B.1, we show the total number of retired A stars having Earths (N_★,mHC) as a function of the Galactocen-tric distance considering different gas giant effects. In the upper panel, the model 1 (GG BAD) results are reported and model 7 (GG GOOD) in the lower one. In Fig. B.2, we show habitability maps for a modified model 1, which accounts for the time necessary for a sustained increase in atmospheric O₂ to significant levels. On Earth, this process needed approximately 2.5 Gyr (Lyons et al. 2014). Specifically, we evaluated delays of 1 Gyr and 3 Gyr.

Fig. B.1 Total number of retired stars having Earths (N★,mHC) as a function of the Galactocentric distance considering different gas giant effects. In the upper panel, the model 1 (GG BAD) results are reported and model 7 (GG GOOD) results in the lower one. Neither model includes stellar migration.

Fig. B.2 Same as the upper panel of Fig. 5 (model 1), but accounting for the time required to form planetary atmospheres. In the upper panel, it is assumed that 1 Gyr is needed for an atmosphere to form after the stellar host’s formation, while in the lower panel, this timescale is extended to 3 Gyr.

Appendix C Observational data

In this appendix, we provide detailed information on the A&S data sample used for the CMDF comparison with model predictions given in Fig. 7 and in Section 4.2. From the Ariel Stellar Catalogue⁶ (Magrini et al. 2022; Tsantaki et al. 2025), we selected stars whose planet(s) have masses M_P < 30 M_⊕, for a total of 57 stars. Then, we complemented such a sample with the SWEET-Cat Catalogue⁷ (Sousa et al. 2021; Santos et al. 2013), selecting stars that are not in the Ariel sample, for a total of 67 stars hosting low mass planets. Both of these catalogues provide precise and, most importantly, homogeneous stellar properties (effective temperature, surface gravity and [Fe/H]), obtained from high-signal to-noise and high resolution spectra, which are essential characteristics for conducting robust population studies (e.g. Adibekyan 2019; Danielski et al. 2022). In particular, we note that the [Fe/H] determination within both catalogues has been tested to be consistent (Brucalassi et al. 2022). Our final sample consists of a total of 124 stars: 106 FGK and 18 M type stars.

Fig. C.1 Orbital radial variation (Δ(Rmean - RGC)) versus galactocen-tric distance (RGC) of the sample of homogeneous determined stars. The marker size shows the orbital eccentricity of the stars (i.e. smaller size means smaller eccentricity), while the color refers to the metallicity [Fe/H] as reported in the colourbar. The blue line marks a null difference between the Rmean of the orbit and the current RGC, and marks the local region i.e. 8 ± 1 kpc. Grey line marks those stars that likely migrated from the outer disc (Δ > 1 kpc), while red dashed line marks those stars that likely migrated from the inner disc (Δ < -1 kpc).

To set such stars within the Galactic context we computed the velocity components and the orbital parameters with the GALPY package (Bovy 2015)⁸, using the Gaia DR3 data. We followed the same procedure as described in Magrini et al. (2022) and Tsantaki et al. (2025). Consequently, we compared their current Galactocentric position with the average position of their orbits (calculated as the average of apogalactic and perigalactic radii). We show in Fig. C.1 the orbital radial variation (R_mean - R_GC) versus their Galactocentric distance (R_GC). While the local population (i.e. those with variation of less than 1 kpc) is the only one showing very low eccentricities, other more eccentric stars show larger variations, which indicates that they have probably originated in the outer or inner part of the Galactic disc. We used such radial variations to divide our stars within local and inner and outer populations.

References

All Tables

All Figures

Fig. 1 Probability of forming planets as a function of metallicity and stellar type. Upper panel: probabilities of forming gas giant planets as a function of [Fe/H] for M, FGK and retired A proposed by Ghezzi et al. (2018) are shown with the black, blue and red lines, respectively. The probabilities have been IMF weighted according to the procedure described in Section 3.1. Fischer & Valenti (2005) relation is also shown with the dashed grey line. Lower panel: probability, PE, of forming Earth-like planets and not gas giants (GG BAD case) as a function of [Fe/H], and indicated in Eq. (6) for different stellar spectral types, shown as solid lines. In contrast, the dashed lines illustrate the new scenario discussed in Section 3.3 and called the GG GOOD case, which assumes that the presence of gas giants enhances the likelihood of terrestrial planet formation. All the cases are computed using the gas giant probabilities of Ghezzi et al. (2018).

In the text

Fig. 2 Predictions of our multi-zone chemical evolution model presented in Section 2.1 as a function of the evolutionary time t and the Galacto-centric distance. Left panel: evolution of the total (CC+Type Ia) SN rates. The red and green dashed horizontal lines indicate <RS N, ⊙> (case 1) and 2 × <RSN,⊙> (case 2) values, respectively, representing the minimum SN rate thresholds tested for SN explosion-induced destruction effects in this study. Right panel : predicted age-metallicity relation. The red-shaded region marks [Fe/H] < -1 dex where, as suggested by Johnson & Li (2012), conditions are considered unfavourable for the formation of terrestrial planets.

In the text

	Fig. 3 Evolution of the probability, PE, of forming Earth-like planets and not gas giants (the ‘GG BAD’ case illustrated in Eq. (6) and in Fig. 1), as predicted by our reference multi-zone chemical evolution model presented in Section 2.1 in the 3D space formed by the PE, [Fe/H] and evolutionary time t for different stellar types.
In the text

Fig. 4 Total number of different spectral stellar types hosting Earthlike planets with minimum habitability conditions (N★,mHC in Eq. (8)) as a function of Galactic distance and time predicted by the chemical evolution model without stellar migration (model 1; see Section 3.1 and Table 1). The values of (N★,mHC) are computed within concentric rings of 2 kpc width. Results for M stars are drawn in the upper panel. FGK and retired A are shown in the middle and the lower ones, respectively.

In the text

	Fig. 5 Total number of FGK stars with minimum habitability conditions (N★,mHC in Eq. (8)) as a function of Galactocentric distance and Galactic time predicted by the model 2; i.e. case 1 SN damage scenario and no stellar migration (see Section 3.1 and Table 1).
In the text

Fig. 6 Temporal evolution of the total number of FGK stars hosting minimum conditions for life (N★,mHC) in the solar vicinity (black dashed lines). Both panels show the contribution from stars born in situ within the solar annulus (Ri = Rf = 8 kpc, dark blue lines) as well as those migrating from different Galactic regions (Ri ≠ Rf). We consider the case 2 SN destruction scenario in the presence of stellar migration, with radial migration strength of σRM = 3.5 kpc (model 3, upper panel) and σRM = 6 kpc (model 4, lower panel), respectively.

In the text

Fig. 7 CMDFs for FGK and M stars hosting habitable terrestrial planets (N★,mHC) in the solar vicinity predicted by models with stellar migration (solid lines). Model results are compared with the A&S data sample described in Appendix C for stars hosting only low mass planets (Mp < 30M⊕, dashed lines). We display the contribution from predicted stars that were born in situ within the solar annulus (Ri = Rf = 8 kpc, solid blue lines), as well as those that migrated from 6 and 10 kpc (solid red and grey line, respectively) to ensure consistency with the observed local, inner, and outer populations. Model 3 (σRM=3.5) results are reported in the left panel. In the right one, we show the ‘model weak’ approach, adopting a weaker radial migration (σRM=1), while keeping all other model parameters identical to those of model 3.

In the text

Fig. 8 Total number of FGK stars having Earths (N★,mHC) as a function of the Galactocentric distance and time considering the case 2 SN destruction scenario. Upper panel: results from the reference model without including stellar migration (model 1). Middle panel: results from the model incorporating stellar migration as described in Section 2.2, with the radial migration strength fixed at σRM = 3.5 kpc (model 3). Lower panel: same as the middle panel but with a radial migration strength of σRM = 6 kpc (model 4). In each panel, we highlight with a green cross the maximum number of FGK stars hosting habitable terrestrial planets predicted at 12 kpc (grey point in the 2D projection).

In the text

Fig. 9 Ratios of the number of FGK stars hosting Earth-like planets (N★,mHC) in models with stellar migration to the reference model without migration, as a function of Galactocentric distance and Galactic time. All considered models adopt case 2 for SN damage. Upper panel: stellar migration model with σRM = 3.5 (ratio of model 3 to model 1 in Table 1). Lower panel : stellar migration model with σRM = 6 (ratio of model 4 to model 1 in Table 1).

In the text

Fig. 10 Ratios between the predicted number of stars having Earths (N★,mHC) of the two different prescriptions for PE as in Fig. 1: GG GOOD and GG BAD cases, respectively. Results for M, FGK, and retired A are reported in the left, middle, and right panels, respectively. Ratios are computed as functions of the Galactocentric distance and Galactic time and all the displayed models consider the case 2 for the SN damage. Upper panels: reference models without stellar migration (ratio between model 7 and model 1 in Table 1). Lower panels: models with stellar migration and σRM = 3.5 (ratio between model 5 and model 3 in Table 1). In each panel the maximum value of the ratio is highlighted with an orange cross.

In the text

	Fig. A.1 Predicted time evolution of the SFR of our multi-zone chemical evolution model presented in Section 4.1 as a function of the evolutionary time, t, and the Galactocentric distance.
In the text

Fig. A.2 Total number of predicted stars (N★) as a function of the Galactocentric distance and the Galactic time. Upper panel: Results from the reference model without including stellar migration. Middle panel: Results from the model incorporating stellar migration, as described in Section 2.2, with the radial migration strength fixed at σRM = 3.5. Lower panel : Same as the middle panel but with a radial migration strength of σRM = 6. In each panel, we highlight with a green cross the maximum number of stars computed at 12 kpc.

In the text

	Fig. B.1 Total number of retired stars having Earths (N★,mHC) as a function of the Galactocentric distance considering different gas giant effects. In the upper panel, the model 1 (GG BAD) results are reported and model 7 (GG GOOD) results in the lower one. Neither model includes stellar migration.
In the text

	Fig. B.2 Same as the upper panel of Fig. 5 (model 1), but accounting for the time required to form planetary atmospheres. In the upper panel, it is assumed that 1 Gyr is needed for an atmosphere to form after the stellar host’s formation, while in the lower panel, this timescale is extended to 3 Gyr.
In the text

Fig. C.1 Orbital radial variation (Δ(Rmean - RGC)) versus galactocen-tric distance (RGC) of the sample of homogeneous determined stars. The marker size shows the orbital eccentricity of the stars (i.e. smaller size means smaller eccentricity), while the color refers to the metallicity [Fe/H] as reported in the colourbar. The blue line marks a null difference between the Rmean of the orbit and the current RGC, and marks the local region i.e. 8 ± 1 kpc. Grey line marks those stars that likely migrated from the outer disc (Δ > 1 kpc), while red dashed line marks those stars that likely migrated from the inner disc (Δ < -1 kpc).

In the text