© The Authors 2024

1 Introduction

The architectures of planetary systems are tightly connected with the properties of their host stars, including their effective temperature, age, and metallicity. This relationship is most pronounced in the established correlation between stellar metallicity [Fe/H] and the occurrence rates of giant planets. Multiple observational studies have shown that close-in giant planets are more commonly found around metal-rich stars (Gonzalez 1997; Fischer & Valenti 2005; Johnson et al. 2010). Additionally, observations suggest a decline in the fraction of stars hosting planets and a decrease in planet multiplicity in line with the increasing mass of the host star (Yang et al. 2020; He et al. 2021). Across the galaxy, stellar populations exhibit significant variations in both age and chemistry (Gilmore & Reid 1983; Prochaska et al. 2000; Reddy et al. 2006; Schuster et al. 2006; Bensby et al. 2014; Masseron & Gilmore 2015), implying that planet demographics could strongly depend on the local galactic environment.

As galaxies grow and evolve, their local chemical profiles undergo changes due to stellar evolution (McWilliam 1997). Throughout their life cycles, stars enrich the interstellar medium (ISM) with heavy elements. This process, along with the hierarchical growth of galaxies, results in a metallicity gradient across the galaxy (Lemasle et al. 2007; Luck & Lambert 2011; Genovali et al. 2014), as well as the formation of distinct stellar populations (Gilmore & Reid 1983; Bland-Hawthorn & Gerhard 2016). The Milky Way (MW), for instance, is typically divided into four components: the galactic bulge, thin disc, thick disc, and halo, each characterized by unique stellar populations with varied chemical compositions. The halo is particularly known for containing some of the oldest and most metal-poor stars in the galaxy, with metallicities ranging from −7 < [Fe/H] < −0.5 (Beers & Christlieb 2005; Frebel & Norris 2015). In contrast, the thin disc hosts a younger stellar population with metallicities near solar levels, [Fe/H] ≈ 0 (Kilic et al. 2017). The thick disc, though not as metal-rich as the thin disc, is distinguished by its α-enhanced stars (Wallerstein 1962). The bulge, representing the most heterogeneous of these populations, spans a broad range of metallicities, including both extremely metal-poor and metal-rich stars (McWilliam 2016).

The majority of exoplanet studies to date have concentrated on stars in the nearby thin disc (Bashi & Zucker 2019). This focus is primarily due to observational constraints, as planets orbiting nearby stars are more readily detected via transit and radial velocity (RV) methods. Despite these limitations, a comprehensive picture of planet formation is emerging (e.g. Pollack et al. 1996; Bitsch et al. 2019; Johansen et al. 2019; Drążkowska et al. 2023). Planetary systems are understood to originate from the protoplanetary discs that encircle young host stars. The prevailing theory is that most planets form through the core-accretion model, where large planet embryos grow by accumulating planetesimals and pebbles. Considering that the chemical composition of these protoplanetary discs reflects that of the host star, it is reasonable to expect a relationship between the architecture of planetary systems and the metallicity of their host stars, and such correlations have indeed been observed (Gonzalez 1997; Santos et al. 2001; Mulders et al. 2015; Petigura et al. 2018).

Planet formation is an intricate process that involves a wide range of physical phenomena, and considerable effort has been dedicated to modeling and simulation work. Key aspects of this process include the formation of planetesimals and planetary embryos, the evolution and eventual dissipation of the protoplanetary disc, and the growth of planets through mechanisms such as accretion, migration, and the development of planetary interiors and atmospheres (e.g. Pollack et al. 1996; Ida & Lin 2004a,b; Bitsch et al. 2019; Johansen et al. 2019; Alessi & Pudritz 2018; Alessi et al. 2020). To align these models with observational data, planet population synthesis has become a crucial methodology (Ida & Lin 2004a; Mordasini et al. 2009a,b). This approach involves generating representative samples of synthetic exoplanets for specific stellar populations, which can then be systematically compared to the actual observed exoplanet data.

The heterogeneity of stellar populations, and their connection to galaxy evolution is an area of recent research that has has benefitted greatly from the advancements in simulation techniques (Crain & van de Voort 2023). Hydrodynamical simulations of galaxy evolution, in particular, have had great success in reproducing the properties of galaxies in the local universe (e.g. Navarro & White 1994; Vogelsberger et al. 2014; Schaye et al. 2015; Dubois et al. 2016), including MW-like galaxies. Amongst the MW-reproducing simulations, one successful approach has been to model the local (galaxy) environment through the use of constrained simulations. These simulations are able to generate MW analogues within environments that closely resemble the Local Group (Hoffman 2009; Gottloeber et al. 2010; Libeskind et al. 2020).

In this work, we aim to combine the recent advances in both galactic modelling and planetary population synthesis. We employ the HESTIA suite of constrained hydrodynamical galaxy formation simulations to study planet populations throughout the MW (Libeskind et al. 2020). HESTIA reproduces the Local Group, including the MW, Andromeda (M31), and surrounding satellite galaxies within the correct cosmographic context. As a zoom-in simulation, HESTIA creates a high-resolution MW-like galaxy with a smoothing length of approximately 220 pc and a star particle mass resolution of 10⁴M_⊙. The simulated galaxy closely matches the MW in terms of total mass, stellar mass, metallicity distribution, and size. We combine one of these simulated galaxies with the state-of-the-art New Generation Planet Population Synthesis (NGPPS) data, obtained from the Bern model for planet formation and evolution, which generates realistic planet populations for a wide range of host star masses and metallicities (Emsenhuber et al. 2021a). Its goal is to study the demographics of exoplanet throughout the MW.

This paper is organized as follows. Section 2 details the Bern model and NGPPS planet population, including the construction of a planet assignment model to associate planet occurrence rates with the mass and metallicity of the host stars. Section 3 describes the simulated MW analogue in HESTIA, comparing the key properties with those of the actual MW and presenting a decomposition of the simulated galaxy into bulge, thin disc, thick disc, and halo. This section also details the integration of the planet assignment model from Sec. 2 into the MW analogue. In Sec. 4, we present the resulting planet populations in different parts of the galaxy, focussing on Sun-like and low-mass stars. These results are discussed in Sec. 5, where we also highlight the limitations of our study. Finally, we summarise our findings in Sec. 6.

2 Bern planet model and planet populations

The third generation of the Bern global model for planet formation and evolution, NGPPS, is based on a semi-numerical approach to study the formation of planets and evolution of planetary systems (Emsenhuber et al. 2021a). The Bern model is based on the core-accretion paradigm (Pollack et al. 1996), which is thought to be the primary formation pathway for low-mass planets, and a majority of massive planets. Specifically, it does not include giant planet formation through gravitational instability (Boss 1997), which is thought to be a potential origin for very wide orbits (>10 AU, Rafikov 2005; Schib et al. 2021) and very massive planets (>4 M_Jupiter, Schlaufman 2018). The third generation of the model encompasses the growth of planets from (1–100) planet embryos with a mass of about 0.01 M_⊕ into planets with masses ranging from Earth-like planets to super-Jupiter giants. The planet embryos are immersed into fluid-like protoplanetary discs, which are described by gas and planetesimal fields from which the planets grow. These fields are described by continuous density fields, coupled to the planet embryos. In the current generation of the model, direct N-body interaction between the growing planets is also included, and Emsenhuber et al. (2021a) have found that embryo-embryo interactions increase the number of massive planets compared to isolated evolution.

The model is initialised at the stage when planet embryos have formed; the growth stages from dust to pebbles and planetesimals are not included. The initial number and distribution of the embryos are therefore free parameters. The simulation starts with a formation phase in which the planet embryos interact with the gas and planetesimal fields. This stage lasts for a fixed period of 20 Myr, after which the disc is considered dispersed. From this point on, the planets evolve isolated, affected by tidal migration and atmospheric escape. The assumed formation phase of 20 Myr is longer than typical protoplanetary disc lifetimes and therefore encompasses most of the growth stage. However, the late-stage giant impact phase that can affect low-mass planets is not fully simulated, which can lead to an underestimation of the masses for rocky planets. Emsenhuber et al. (2021a) found that for a minimum-mass solar nebula-like surface density of solids, planet formation is mostly complete within ~1 AU after 20 Myr, while more distant planets have not yet reached their final state. The model provides a wealth of data about the planet populations in the systems, including their mass, distance to the host star, radii, luminosities, and evolution tracks.

A wide range of physical processes are included in the model and, correspondingly, a large number of parameters need to be set. Some of these parameters are fixed to specific values (Emsenhuber et al. 2021b, see Table 1); however, some parameters are considered to be variable and a suite of simulations has been run varying these values. Four quantities are varied randomly and dubbed Monte Carlo variables.

Initial gas disc mass. The initial mass of the gas disc around the star, denoted as M_g, provides the building material for planet formation, with a more massive disc potentially leading to more numerous and larger planets. Emsenhuber et al. (2021b) opted for log-normal distributions (with median μ = −1.49, and standard deviation σ = 0.35) based on Very Large Array (VLA) observations reported by Tychoniec et al. (2018), which focussed on Class I discs in the Perseus star-forming region. They restricted the range of masses to 0.004–0.16 M_⊙ to guarantee that the discs have non-negligible mass and are self-gravitationally stable.

External photo-evaporation rate. The external photo-evaporation ${\dot{M}}_{\text{wind }}$ measurement is the main parameter determining the lifetime of the protoplanetary disc. While the inner part of the disc is evaporated by the internal radiation of the host star, in the outer regions (where most of the mass resides) evaporation is initiated primarily by the far-ultraviolet (FUV) radiation of nearby massive stars (Matsuyama et al. 2003). Emsenhuber et al. (2021b) assumed a fixed disc viscosity and a normal distribution (μ = −6, σ = 0.5) for log ${\dot{M}}_{\text{wind }}$ to reproduce observed protoplanetary lifetimes. With their parametrisation 50% of discs will have a lifetime of <4 Myr and virtually all discs are dispersed within 10 Myr.

Dust-to-gas ratio. The dust-to gas ratio, f_D/G = M_s/M_g, determines the amount of solid material in the disc, which is crucial for the growth of solid planets and the cores of gas giants. In the Emsenhuber et al. (2021a) model, the dust-to-gas ratio is assumed to be a function of stellar metallicity given by $\frac{f_{\mathrm{D}/\mathrm{G}}}{f_{\mathrm{D}/\mathrm{G},\odot }}={10}^{[\mathrm{F}\mathrm{e}/\mathrm{H}]}$, where [Fe/H] is the iron abundance and f_D/G⊙ = 0.0149 is the dust-to-gas ratio of the sun. Emsenhuber et al. (2021b) chose a normal distribution (μ = −0.02, σ = 0.22) for the [Fe/H] distribution truncated to the range [−0.6, 0.5]. This encompasses the vast majority of stars in the solar neighbourhood. However, in the galactic centre and stellar halo metal-poor stars exhibit significantly lower [Fe/H] abundances. This needs to be considered for the HESTIA planet population model.

Inner edge of the gas disc. The position the inner edge of the protoplanetary (gas) disc, r_in, is important for the final locations of close-in planets, since gas drag is a major driver of migration. Emsenhuber et al. (2021b) chose a log-normal distribution (μ = −1.26, σ = 0.2) for this quantity.

The mass of the planetesimal (i.e. solid material) disc, M_s, is not a Monte Carlo variable itself, but of central importance to planet formation since it provides the building blocks for rocky planets and giant cores. Within the Bern model, this is fully characterised by the mass of the gas disc and the dust-to-gas ratio, M_s = f_D/G · M_g. Furthermore, the embryos are initialised with a a mass of 10⁻² M_⊕, and distributed uniformly in the logarithm of distance between r_in and 40 AU. In addition to these quantities, two more parameters are varied systematically:

Initial number of embryos. The initial number of lunar-mass planetary embryos N_Embryo is crucial for the final number and properties of the planetary systems. Since the Bern model does not include the initial growth stage of dust and pebbles, this quantity cannot be predicted and needs to be set manually. Emsenhuber et al. (2021b) study the evolution for 1, 10, 20, 50 and 100 embryos. We will however only consider the range of 10–100 embryos, since they were created with a consistent sample of Monte Carlo variable and are thus more easily comparable.

Host star mass. The host star mass (M_*) is fixed to 1 M_⊙ in the work by Emsenhuber et al. (2021b). However Burn et al. (2021) extend the analysis to low mass stars with 10%, 30%, 50% and 70% of the sun’s mass. The host star’s mass affects the planet evolution in direct ways, for instance, through radiation, but is also correlated with other relevant quantities, such as the mass of the protoplanetary disc. In the work by Burn et al. (2021), only a N_Embryo = 50 run is considered, the Monte Carlo variables follow the same distribution as the for the Sun-like star but re-scaled through some fixed prescriptions, with M_g ∝ M_* and r_in ∝ $M_{\star }^{1/3}$. The metallicity and external photo-evaporation-rate distributions are assumed to be identical to the solar-mass case. Stars more massive than the sun were not studied, due to their relative rarity and short lifespan.

2.1 Original NGPPS sample and planet types

In the original set of publications, Emsenhuber et al. (2021a) have performed and analysed a set of population runs based on the parameter distributions described before. In their work, they fix the host star mass to M_* = 1 (Emsenhuber et al. 2021b), and run simulations for the different assumed initial embryo numbers. For each fixed value of N_Embryo, they draw 1000 random samples from the Monte Carlo variable distributions and run their simulations with these parameters. In a subsequent paper, Burn et al. (2021) perform a similar analysis for a fixed number of embryos N_Embryo = 50 and varying the host star mass between 0.1 and 1 M_⊙, using the same sample of Monte Carlo variables but re-scaled based on the host star mass.

Both Emsenhuber et al. (2021a) and Burn et al. (2021) classify the planets into categories based on their mass at the end of the formation phase. These encompass Earth-like, Super-Earth, Neptunian, Sub-giant, and giant, as detailed in Table 1. Figure 1 shows the distribution of the Monte Carlo variables and number of giant planets found in this system for the M_* = 1 and N_Embryo = 50 run.

Fig. 1 Distribution of Monte Carlo variables for original sample from Emsenhuber et al. (2021b), and effect on number of giant planets for the M* = 1 and NEmbryos = 50 model run. The colour marks the number of planets in the system after 20 Myr of evolution. Shown are the correlation between the gas mass, Mg, metallicity, [Fe/H], inner protoplanetary disc radius, rin, and external photo-evaporation rate, ${\dot{M}}_{\text{wind }}$.

2.2 Correlations between Monte Carlo variables and planet type abundances

Of primary concern for the study of planet populations throughout the galaxy is the relation between the number of planets in a given category and the chosen parameter, since those expected to vary systematically as a function of galactic environment. For example, Fig. 1 indicates that higher values of [Fe/H] and M_g lead to a larger number of giant planets. To get a feeling for the general relations, the correlation coefficients between the planet types and Monte Carlo variables are shown in Fig. 2 for the M_* = 1 and N_Embryos = 50 population. It is easy to see that the more massive planet types are positively correlated with a larger initial disc mass, while low-mass planets have a negative correlation. The same holds true for the metallicity. In contrast, correlations with the external photo-evaporation rate and inner disc edge are much weaker for the considered range of values. These results persist for initial values for N_Embryos between 10 and 100, and similarly for lower host star masses, although the correlations become weaker with decreasing host mass.

2.3 Planet assignment model

To populate the star particles in the HESTIA galaxy simulation with planets, we need a model that connects the properties of the host star (or stellar population) with the properties of the planets forming around the star. We constructed a simple assignment model that connects the occurrence rate of planets in each category with the Monte Carlo variable, where the occurrence rate is defined as: $$\text{Occurrence\hspace{0.17em}rate}=\frac{\text{Total\hspace{0.17em}number\hspace{0.17em}of\hspace{0.17em}planets}}{\text{Total\hspace{0.17em}number\hspace{0.17em}of\hspace{0.17em}stars}}.$$(1)

We note that this quantity differs from the fraction of planethosting stars, which compares the total number with systems with a planet of a given type compared to the total number of stars.

We determine the number of planets per category (as defined in Table 1) across the 1000 Monte Carlo runs (for a fixed N_Embryo and M_*). The results are smoothed and interpolated using a k-nearest neighbours algorithm with k = 30 (chosen via a ten-fold cross-validation on the original sample).

Considering that the photo-evaporation rate and inner disc radius exert minimal influence on planet classifications, we chose to further reduce the model by incorporating only the initial gas disc mass and [Fe/H] as occurrence rate predictors. This increases interpretability of the model at the cost of slightly reduced predictive power. The predicted occurrence rate of planets for the full model (incorporating all four Monte Carlo variables) and the reduced model (using only M_g and [Fe/H]) are shown in Table 2 for the M_* = 1 and N_Embryos = 50 population. The reduced model yields slightly fewer low-mass planets and more giant ones, as elevated photo-evaporation rates may shorten disc lifetimes and thus stop gas accretion early. These effects are however well within the variability of the model.

In Fig. A.1, the occurrence rates of planets are shown as a function of M_g and [Fe/H], for the M_* = 1 and N_Embryos = 50 population. Within the Bern model, occurrence rates of the different planet types are strongly dependent on the initial planetesimal mass in the disc (diagonal lines in Fig. A.1, since M_s = f_D/G · M_g), tend to peak around some value of M_s (dependent on the planet mass), and taper off to either side. The peak value increases as the characteristic mass of the planet type increases, suggesting an optimal range of available solid material to form planets within specific mass ranges. For Earth-like planets, the model indicates a value ~0.3 M_Jup, while Super-Earths are preferably formed in discs with ~0.7 M_Jup of available solid material. Sub-giant and giant planets in particular almost exclusively occur if the solid disc mass >1 M_Jup. Moreover, Fig. A.1 also shows the planet occurrence rates in relation to metallicity [Fe/H] alone. Notably, giant planets exhibit a much more pronounced metallicity-dependence in their occurrence rates compared to low-mass planets, in line with observations.

As part of the original publication series on NGPPS, Schlecker et al. (2021) performed a similar analysis on the predetermination of planet types from disc properties based on the Bern Model. Their approach differs from the one taken here in that they do not pre-define planet categories but use an unsupervised learning algorithm to categorise planets based on their mass, radius, and semi-major axis. Using a Gaussian Mixture Model, they identified four planet categories, which they call Neptunes, icy cores, giant planets, and (Super-)Earths. They found that the strongest predictor for a planet’s type is the location of its initial embryo, and emphasised that consistent treatment of planetary embryo formation is a priority in the future development of the Bern model. Given a fixed starting location, they found that giant planet formation is strongly dependent on the mass of the available solid material, and that the photo-evaporation has no strong correlation with planet type for the considered range of values. They interpret the latter finding as an indication that planet formation, on average, concludes on shorter time scales than the lifetime of the disc within the NGPPS population. Our results are therefore consistent with their findings despite the differences in the approach.

Fig. 2 Kendall correlation matrix between Monte Carlo variables and planet types for the M* = 1 and NEmbryos = 50 run (Emsenhuber et al. 2021b).

2.4 Influence of the number of embryos and stellar mass

One aspect that merits further interest is how the relation of planet populations to the systematically varied quantities. In Fig. 3, we show the total number of planet per category in all the simulation runs with a varying number of initial embryos (and fixed host star mass, M_* = 1M_⊙). The number of massive planets (Neptunians, Sub-giants, giants) are nearly constant for all runs with N_Embryo ≥ 10. The number of Earth-like planets and Super-Earths, on the other hand, increases significantly with the number of embryos. Figure 4 shows a similar figure for a fixed number of initial embryos (N_Embryo = 50) and varying host star mass. Here, the number of Earth-like planets does not depend strongly on the host star mass for masses between 0.1 and 1 M_⊙. The more massive planets however become more frequent for massive stars. Giant planets, in particular, only occur for host star masses >0.3 M_⊙. This agrees with our previous findings that occurrence rates depend primarily on the available solid mass. Burn et al. (2021) assume that the disc mass scales linearly with the mass of the host star. Low-mass stars therefore are less likely to possess the ~1M_Jup solid material needed to form giant planets within the Bern model.

Fig. 3 Variations of planet occurrence rates per category given different numbers of initial embryos, M* = 1 (Emsenhuber et al. 2021b).

Fig. 4 Variations of planet occurrence rates per category given different host star masses, NEmbryos = 50 (Burn et al. 2021).

2.5 Comparison with observations

The relationship between stellar metallicity [Fe/H] and occurrence rates of giant planets has been well-established over the years (Gonzalez 1997; Santos et al. 2001; Fischer & Valenti 2005; Johnson et al. 2010). Similarly, close-in (P < 10 days) small (<2 M_⊕) planets are also more likely to be found around metal-rich stars (Mulders et al. 2016; Lu et al. 2020), hinting at a metallicity-dependence for rocky planets as well, albeit much weaker than for giant planets. Petigura et al. (2018) find that warm Super-Earths with a period between 10 and 100 days have a nearly constant occurrence rate over metallicities in the range −0.4 < [Fe/H] < 0.4, while warm sub-Neptune occurrence rates double over the same range. The radii of planets have also been found to increase with host star metallicity (Narang et al. 2018; Swastik et al. 2022). This agrees well with our assignment model, where the metallicity-dependence is similarly weakened for low-mass planets (Fig. A.1). Giant planets are generally found to be rare around metal-poor stars (Buchhave et al. 2012; Thorngren et al. 2016), which is also reproduced in the assignment model.

Howard et al. (2012) and Mulders et al. (2015) find that the occurrence rate of Earth- to Neptune-sized close-in planets decreases for increasing stellar mass and temperature. He et al. (2021) and Yang et al. (2020) suggest that the fraction of stars with planets changes from ~0.3 for F-type stars to ~0.8 for K-types. Our assignment model, on the other hand, predicts an increase in both Super-Earth and Neptunian occurrence rates with increasing stellar mass (although the occurrence rate for Earth-like planets slightly decreases for host star masses between 0.7 and 1 M_⊙), therefore contradicting the observations. Van der Marel & Mulders (2021) suggest that gap formation in the protoplanetary disc caused by giant planet formation could impede inward migration of lower mass planets, leading to the observed anti-correlation. In the subsequent analysis, we therefore only compare planet population at fixed host star mass M_*, rather than across masses in order to avoid biases based on this spurious behaviour of the model.

3 Milky Way analogue in HESTIA

High-Resolution Environmental Simulations of The Immediate Area (HESTIA) is a suite of hydrodynamical cosmological simulations designed to replicate the Local Group in its correct cosmographic context (Libeskind et al. 2020). It is based on the Arepo code for gravitational N-body systems and magnetohydrodynamics, while the baryonic physics is handled by the Auriga (Grand et al. 2017) model of galaxy formation.

The initial conditions are derived from cosmographic observations to faithfully reconstruct the Local Group, both in terms of its internal structure and nearby large-scale structures. It accurately recreates prominent features of the local Universe such as the Virgo Cluster and Local Void, including their appropriate masses and relative locations to the Local Group (Tully et al. 2009; Karachentsev & Nasonova 2010). Centered on the Local Group midpoint, the simulation contains a high-resolution region with a radius of 3–5 Mpc. An ensemble of dark matter-only and low-resolution hydrodynamical simulations were run in order to construct a Local Group-like pair of MW and Andromeda analogues with accurate mass ranges and separations. The three most promising simulations have been re-run at a higher resolution. These analogues closely resemble their observational counterparts in terms of total mass, mass ratio, stellar disc mass, morphology separation, relative velocity, rotation curves, bulge-disc morphology, satellite galaxy stellar mass function, and satellite radial distribution. In terms of their formation history, the MW and Andromeda analogues grow slower than than comparable halos in unconstrained (McBride et al. 2009; Libeskind et al. 2020). In the highest-resolution runs, the attained mass resolutions of the individual particles are approximately m_DM ~ 10⁵, m_gas ~ 10⁴, and m_stars ~ 10⁴, while the spatial resolution is around 220 pc. It ought to be highlighted that the resolution is not sufficient to resolve individual stars, rather the star particles are treated as single stellar populations with uniform age and metallicity.

In this study, we focus on the MW analogue in the highest-resolution 37_11 run, which will be described in more detail in the following sections.

3.1 General comparison

As described in the previous section, the MW analogues in HESTIA have been purposefully crafted to mirror a galaxy that is MW-like in many properties relevant for galaxy formation and evolution. They were, however, not constructed to be a perfect copy of the MW, and some of the properties that are more relevant for planet formation were understandably not considered in the cosmological context they were built for. We therefore compare some of the more local and structural properties of our chosen MW analogue in more detail.

The virial mass and stellar masses of the analogue closely match those of the MW within observational constraints (see Table 3). This is unsurprising, given that the simulations were specifically designed to achieve this result. Assuming a Chabrier (2003) IMF, the stellar mass corresponds to roughly 135 billion stars in the analogue, of which 2.3 billion stars are main-sequence stars within 5% of the mass of the sun. Libeskind et al. (2020) analyzed the structural parameters of the analogue by fitting a surface brightness profile, which they modelled as the sum of a Sérsic profile describing the bulge and an exponential profile for the disc. The obtained effective bulge and disc radii are are in good agreements with those of the MW (Table 3). While the Sérsic index of the MW bulge lacks precise constraints, the simulation value of 1.80 is comparable to the observational value of 2.2 ± 0.3 for the Andromeda galaxy (Courteau et al. 2011), the closest comparable spiral galaxy.

Furthermore, as discussed in Sec. 2.3, a major factor in the planet assignment model is the metallicity of the host star. Within HESTIA, 9 elemental species are tracked throughout the simulation: H, He, C, N, O, Ne, Mg, Si, and Fe. From these, we can calculate the iron abundance [Fe/H] for any individual star particle. We assume all stars within a star particle (a single stellar population) share the same metallicity. The metallicity gradient of the MW analogue between 2 and 10 kpc from the galactic centre is −0.027 dex kpc⁻¹. This is a shallower slope compared to the values determined for the MW in observational studies (Lemasle et al. 2007; Luck & Lambert 2011; Genovali et al. 2014). Although, a more recent study by Lemasle et al. (2018) found MW values of −0.045 ± 0.007 dex kpc⁻¹ and −0.040 ± 0.002 dex kpc⁻¹ using F/10 double-mode Cepheids, closer to that of the analogue. Nevertheless, within the HESTIA analogue, an average solar neighbourhood metallicity value of [Fe/H] ~ −0.1 to 0 (Haywood 2001) is reached at a distance of ~7–8 kpc from the galactic centre, in good agreement with the distance of the solar neighbourhood from the MW centre (~8.2 kpc, McMillan 2017).

3.2 Galaxy components

Within the solar neighbourhood, the stars in the MW are generally categorised into three groups: the thin disc, the thick disc, and the stellar halo (Gilmore & Reid 1983). These populations are distinguished by their kinematics (Piffl et al. 2014; Sanders & Binney 2015), stellar age (Schuster et al. 2006), and chemical composition (Bensby et al. 2014; Masseron & Gilmore 2015; Hawkins et al. 2015). Thin disc stars have lower orbital velocity and higher velocity dispersion compared to stars in the thick disc. Chemically, thick disc stars are metal-poor and α-enhanced relative to thin-disc stars (Prochaska et al. 2000; Reddy et al. 2003, 2006). Halo stars, in contrast, have even larger velocity dispersion (Bond et al. 2010) and are metal-poor. Near the galactic centre, there is also a population of bulge stars characterised by high velocity dispersion (Portail et al. 2015) and a wide range of metallicities, spanning from [Fe/H] ~−1.5 to [Fe/H] ~+0.5 (Zoccali et al. 2008). These variations in formation history and chemistry are likely to influence the planet populations within these galactic components (Bashi & Zucker 2022; Nielsen et al. 2023). Therefore, we decompose the HESTIA analogue into four categories, broadly consistent with the four MW components.

While observational surveys benefit strongly from the separation of components based on combined kinematics and chemistry criteria (Bashi & Zucker 2019), we adopt a simpler classification approach in this work. We categorised objects based on the binding energy and circularity of the star particles, aligning with previous studies on galaxy formation simulations (Abadi et al. 2003; Yu et al. 2023) by employing the morphological decomposition code MORDOR by Zana et al. (2022). Face- and side-on projection maps of the components are shown in Fig. 5.

The determined galaxy components broadly align with anticipated properties and observations of the MW. The bulge, thin disc, and thick disc closely match those of the MW in terms of stellar mass, mean stellar age, and mean [Fe/H], as well as the minor-to-major axis ratio. Discrepancies are found for the stellar halo, which is more massive and metal-rich in the HESTIA analogue compared to the MW. This can likely be attributed to a misclassification of some bulge and disc stars based on the purely kinematic decomposition. The key attributes are summarised in Table 4, and some distributions are shown in Fig. A.2.

3.3 Integrating the planet model

To integrate the planet assignment model (Sec. 2.3) into the HESTIA analogue, we associate each star particle with a certain number of planets for each planet type, based on the properties of the star particle. With a mass of approximately 10⁴ M_⊙, we associate a star particle with a single stellar populations of uniform age and metallicity (Chevance et al. 2023).

The hydrogen and iron masses of the star particles are inherited from the gas cells from which the particles originate, meaning each particle is equipped with an [Fe/H] value at birth. This value tracks the chemical evolution of the galaxy and can be fed directly into the planet assignment model.

Since the resolution of the simulation is on the order of approximately 200 pc, processes that potentially influence the formation and evolution of the protoplanetary discs cannot be tracked. Instead, we associate every star particle with a random initial disc mass M_g sampled from the distribution provided by Emsenhuber et al. (2021b), and scaled to the host star mass according to the prescription given by Burn et al. (2021). This approach includes the assumption that the mass distribution of discs does not change significantly in different regions of the galaxy. The implications of this assumption are discussed in Sec. 5.3.

As discussion in Sec. 2.5, the assignment model predicts an increase of low-mass planet frequency with stellar mass, which is a behaviour not seen in observations. Since low-mass stars are common than massive ones, this could bias our conclusions if not controlled for. We therefore chose to analyse the planet populations for different host star masses, M_*, independently, and not across stellar masses. To calculate the number of eligible stars associated with the star particle, we integrate the IMF between M_* ± 5%, where M_* is the chosen (fixed) host star mass. Since planet formation is modeled only for main-sequence stars, we exclude stars with a lifetime shorter than the particle age.

In short, the planet prescription at fixed host star mass M_* is as follows. First, we choose an NGPPS population by fixing N_Embryos and M_*. Then, we obtain the age, [Fe/H], and mass of the star particle. Next, we sample the Monte Carlo variable M_g from the distribution given by Emsenhuber et al. (2021b). We integrate the (normalised) Chabrier (2003) IMF, ξ, in the interval M_* ± 5%, and multiply by the mass of the star particle to obtain the number of stars. Post-main sequence stars are removed based on the age of the star particle¹. We then assign the number of planets n per star using the planet assignment model for each category (Sec. 2.3). Finally, we multiply the number of stars by the number of planets per star².

Fig. 5 Side and face-on maps of the distribution of stars for the different galaxy components.

3.4 Metal-poor stars

The [Fe/H] distribution in the HESTIA analogue covers a wide range of values, with a majority of particles residing between [Fe/H] = −2.5 and [Fe/H] = 1. This interval is considerably broader than the range considered in the NGPPS database, which is limited to −0.6 < [Fe/H] < 0.5. This range was selected to cover the typical metallicities within the solar neighbourhood. In the HESTIA analogue, 80% of the star particles conform to this range. However, around 48% of the star particles in the thick disc and stellar halo fall beneath the NGPPS lower threshold.

Drawing from trends highlighted in the planet assignment model discussed in Sect. 2.3, considerably larger disc masses would be needed at such low metallicities to provide sufficient solid material for planet formation, especially for planets more massive than 2 M_⊕. In addition, observations suggest that protoplanetary discs disperse faster in low-metallicity environments (Yasui et al. 2009; Yasui 2021). From an observational standpoint, only 26 planets have been confirmed orbiting stars with [Fe/H] < −0.6 at the time of writing, with a minimum metallicity of [Fe/H] = −1.0, providing very limited information on the effects of low-metallicity environments on planet formation³.

Therefore, in our analysis, we consider two extreme scenarios. Firstly, we presume the planet frequencies to be unchanging at extremely low metallicities, adopting the same values as in the [Fe/H] = −0.6 case. Secondly, we assume that short disc lifetimes and a lack of solid material render planet formation highly unlikely to impossible, and assign an occurrence rate of zero to those particles.

Fig. 6 Occurrence rates of different planet types in the four galaxy components. The occurrence rates are calculated for the M* = 1M⊙ and NEmbryo = 50 run.

4 Results

4.1 Planets around Sun-like stars

Focussing initially on planets orbiting Sun-like stars, which (for our purposes) are stars within ±5% of the Sun’s mass, we observe distinctive patterns in occurrence rates across different planet categories. For the N_Embryo = 50 run, the occurrence rates in the different components are shown in Fig. 6. The planet populations in the bulge and thin disc are remarkably similar, the same is true for the thick disc and halo. This similarity is primarily attributed to these components covering similar metallicity regimes, which significantly influences planet occurrence rates in our model. A list of the occurrence rates can be found in Table A.1.

These occurrence rates of planets in the thin disc stellar population align closely with those produced by Emsenhuber et al. (2021b), which is consistent with their initial aim to model the nearby (primarily thin disc) planet population. The population in the bulge has similar occurrence rates, with the bulge favouring more massive planets slightly. The most common type of planets in these regions are Earth-like planets and Super-Earths with an occurrence rate of 3.5–4 planets per star. Sub-giants and giants have a frequency of 0.1–0.4, around one order of magnitude lower. However, the average age of a planetary system within these components differ. The average planet in the bulge is 7.6 Gyr old, while the average planet in the thin disc is younger at an age of 5.8 Gyr. Within a component, low-mass planets tend to be younger than more massive planets, due to the time scales of metal enrichment needed for giant planet formation (see Table 5).

The thick disc and halo bear a distinct planet population compared to the more metal-rich regions. The lower metallicity of the stars favour low-mass planet formation over more massive ones. The occurrence rate of Earth-like planets is a 1.5 times larger than for the bulge and thin disc component, at 5.5 planets per star, with a comparable drop in the occurrence rate of super-Earths to ~2.5 planets per star. More striking is the difference in occurrence rates for Neptunian and giant planets. Neptunian planets are five times more common in the thin disc and bulge compared to the thick disc and halo, with frequencies around 0.6. This ratio increases with increasing planet mass, the occurrence rates of giant planets with masses >300 M_⊕ are 20 times higher in the metal-rich components compared to the thick disc and halo. The planets in these regions are also older on average, at 7.7 Gyr for the thick disc and 8.2 Gyr for the halo (see also Table 5).

Taking a higher level view, we analysed the distribution of planet types across the whole galaxy. We show radial profiles in Fig. 7, while 2D distribution maps are shown in Fig. 8. Giant planets have a higher concentration in the centre of the galaxy, correlating with the region’s higher metallicity. Figure 7 shows the normalised planet density as a function of distance from the galactic centre. The spatial density of planets depends on the stellar density and planet occurrence rates. Thus, unsurprisingly, the highest density of planets is found towards the galactic centre due to the large stellar density in this region. Moving away from the galactic centre, the number of giant planets drops faster than for lower mass planets, due to the higher metallicity-dependence of the giant planet occurrence rates. In total, 90% of giant planets and 84% of Earth-like planets can be found within 8 kpc from the galactic centre, roughly the distance of the Sun from the MW centre.

Further, we find that the total number of planets is dependent on the assumed number of initial planet embryos. While the dependence is relatively weak for massive planets – echoing the conclusions of Emsenhuber et al. (2021a) – a more pronounced dependence is observed for lower mass planets. We calculated the total number of planets within R = 42 kpc of the HESTIA galaxy for N_Embryo = 10, 20, 50, and 100 and tabulate them in Table 6. Depending on the assumed number of embryos, this yields between 9 × 10⁹ and 25 × 10⁹ planets around Sun-like stars in the HESTIA analogue. It is probable that the number of embryos is not universal for protoplanetary systems and depends, amongst other things on the metallicity (see Sect. 5.3.2). The N_Embryo = 10 and N_Embryo = 100 cases can however be treated as a lower and upper estimate for the total number of planets.

Fig. 7 Spatial density of planets as a function of distance from the galactic centre (normalised to R = 0). The points correspond to the estimated values per radial bin, while lines are smoothing splines to show trends. The grey, dash-dotted line corresponds to the normalized density of stars.

Fig. 8 Maps of distribution of planets for each planet Type in the HESTIA MW analogue for stars with M* = 1 M⊙. The assignment is based on the NEmbryo = 50 NGPPS population.

4.2 Low-metallicity environments

Our analysis in the previous section assumed the occurrence rates are unchanging for metallicities below −0.6, since the original NGPPS sample is constrained to this region of parameter space. However, there are observational and theoretical considerations suggesting the occurrence rates differ in this regime (see Sec. 3.4), and the flat extrapolation is likely to overestimate the planet occurrence rates, especially in the thick disc and halo.

To quantify this effect, we considered an alternative metallicity model, where we assigned zero planets to stars below the metallicity threshold [Fe/H] = −0.6. We compared the occurrence rates to the previous prescriptions. Figure 9 shows the ratio between the occurrence rate for the fiducial and this truncated model for the different planet types and galaxy components. The variations are minimal for the bulge and thin disc, since most stars in these regions are above the metallicity threshold. For the low-metallicity environments, the occurrence rates of low-mass planets in the halo is up to 50% lower when assuming this more restrictive metallicity assumption.

Interestingly, the occurrence rate for giant planets does not significantly differ between the two models, not even in the low-metallicity galaxy components such as the halo, even though the occurrence rates of these planets have the highest metallicity dependence. This stems from the fact that, within our model, giant planets already show near-zero occurrence rates at a metallicity of [Fe/H] = −0.6. Consequently, the variation in the metallicity models does not significantly change their estimated frequency. In contrast, the occurrence rates for Earth-like planets show severe reduction. Given that the occurrence of lower mass planets is only weakly dependent on metallicity, and considering that approximately half of the star particles in the halo have metallicities below our cutoff, the occurrence rate for Earth-like planets is effectively halved in this scenario.

It is important to note that the actual MW halo is more metal-poor than that of the HESTIA analogue. The effect demonstrated here is therefore likely to be even more pronounced in reality.

Fig. 9 Ratio between the occurrence rates in the fiducial metallicity model (flat extrapolation) and the alternative case where we assign zero planets. The ratios are calculated for M* = 1 M⊙ and NEmbryo = 50.

4.3 Planets around intermediate- and low-mass dwarfs

Within the Bern model and NGPPS populations, the average number of planets per star decreases with host star mass. This is at odds with some observations which find higher occurrence rates around colder and less massive stars (Howard et al. 2012; Mulders et al. 2015; Yang et al. 2020; He et al. 2021). Further, while the NGPPS populations created by Burn et al. (2021) span the host mass range from 0.1 to 1 M_⊙, to reduce computing time only N_Embryo = 50 initial embryos have been considered. We therefore are unable to analyse the sensitivity of occurrence rates on N_Embryo for low-mass stars. Consequently, care should be taken when comparing the occurrence rates and absolute numbers for these stars directly to the ones discussed for Sun-like stars in previous section. Nonetheless, the relative frequencies still contain valuable information about the planet populations around cooler stars.

For lower mass stars, the correlation between planet occurrence rates with metallicity becomes weaker, although the average metallicity to form a given planet type is larger (Burn et al. 2021). Consequently, the planet populations in the different regions of the galaxy are more homogeneous compared to those around Sun-like stars. The occurrence rates are shown in Fig. 10 and Table A.1.

Around intermediate-mass dwarfs (M_* = 0.5 M_⊙), giant planets are extremely rare in the thick disc and halo with occurrence rates of 0.001–0.002. These values increase by a factor of 10 for the metal-rich regions. The occurrence rate of Earth-like planets for these stars is rather insensitive to the metallicity and stays between four and five planets per star. Super-Earths on the other hand are more metallicity-dependent with occurrence rates around one planet per star for the metal poor region, three times lower than for the metal rich regions.

For small M-type dwarfs with masses M_* = 0.3 M_⊙, no giant planets with masses >300 M_⊕ are found in any of the regions of the galaxy. The accretion discs around these stars simply do not provide enough solid material to form massive planetary cores within the lifetimes of the gas disc. Earth-like planets and Super-Earth occurrence rates are more uniform across the different regions, with occurrence rates of 4.5–5 and 1.7–2.1, respectively.

If the occurrence rates across planet masses are taken at face values, Earth-like planets are most common around intermediate mass stars M_* = 0.5 M_⊙ for metal-rich regions, while they are most common around Sun-like stars in the thick disc and halo. Neptunians and giants are most common around more massive stars, irrespective of galactic environment.

The weaker metallicity dependence also affects the trends in the occurrence rates with stellar mass (see Table A.1). For massive planets (Neptunians, sub-giants, giants), the occurrence rates consistently increase with stellar mass in all components. In contrast, the occurrence rate of Earth-like planets decreases with stellar mass in the metal-rich regions (bulge, thin disc) but increases in the metal-poor regions (thick disc, halo). For Super-Earths, the occurrence rates have a minimum at M_* = 0.5 M_⊙ in the metal-poor regions.

Fig. 10 Occurrence rates of planet around host star with masses of M* = 0.5 M⊙ and M* = 0.3 M⊙, for the NEmbryo = 50 run.

5 Discussion

5.1 Planet populations in the galactic context

Our results show that galactic chemical evolution and variations in stellar populations profoundly affect planet demographics throughout the galaxy. For Sun-like stars, we observe that metal-rich regions, specifically the bulge and thin disc of the MW, harbour the largest number of planets with masses >M_⊕, in both absolute numbers and relative frequencies. In contrast, the metal-poor regions, namely, the thick disc and halo, show a significant suppression in giant planet formation, which is attributed to the lack of solid material building blocks, which are crucial for forming the massive planet cores before the dispersal of the protoplanetary gas discs. Specifically, the occurrence rate of giant planets around Sun-like stars is an order of magnitude higher in the metal-rich regions. Earth-like planets and Super-Earths are more insensitive to the metallicity variations, but still affected. Earth-like planets around Sun-like stars are actually most frequent in the low metallicity environments of the thick disc and halo, being around twice as common as in the thin disc. Super-Earths show the opposite trend and are about half as frequent in the low metallicity regions.

For lower mass stars, with M_* = 0.5 M_⊙ and 0.3 M_⊙, we observe that planet populations are less dependent on metallic-ity, leading to a more uniform distribution across the galaxy. The lower amounts of solid material available in the protoplanetary discs around these stars generally results in host stars having a higher average metallicity compared to the Sun-like case for the same planet type. Thus, giant and sub-giant planets remain significantly rarer in low-metallicity regions. No giant planets with masses greater >300 M_⊕ are found around any stars of masses 0.3 M_⊙, regardless of the galactic environment.

Occasionally, the concept of a galactic habitable zone is discussed in the literature (Gonzalez et al. 2001; Lineweaver 2001; Prantzos 2008; Gowanlock et al. 2011), suggesting that life-supporting conditions are more likely in certain galactic regions, including some of our previous work involving hydrodynamical galaxy simulations (Forgan et al. 2017). The key idea is that the galactic centre, with its extreme environmental conditions, is less conducive to life due to higher radiation levels and more frequent stellar encounters. Conversely, the outer regions are thought to be too metal-depleted for habitable planet formation. While we have not explored the radiation-related arguments in the inner region, we find no evidence of an outer limit. Low-mass, rocky planets remain abundant throughout the galaxy far into the stellar halo. Moreover, some of these models suggest that giant planets might be detrimental to habitability or the formation of Earthlike planets (Lineweaver 2001; Forgan et al. 2017). Under this assumption, stars in the far-out metal-poor regions could be considered favourable, since giant planet formation is suppressed in these regions.

5.2 Empirical occurrence rates

The robustness of our results applied to the real MW depends on the accuracy of the NGPPS occurrence rates compared to true occurrence rates. Matuszewski et al. (2023) conducted a planet yield analysis for the PLATO mission (Rauer et al. 2014), and compared the occurrence rate estimates around FGK stars from various sources. This comparison includes the NGPPS N_Embryo = 100 run, and two empirical models based on Kepler data – one by Hsu et al. (2019) and the other by Kunimoto & Matthews (2020). Unlike the mass-based categorization in NGPPS, these models focus on planet radius, aligning with the observational constraints of transit missions.

The NGPPS model predicts an overall occurrence rate of 5.62 planets per star within an orbital period of 500 days, slightly higher than the 5.02 planets per star estimated by Hsu et al. (2019). In contrast, Kunimoto & Matthews (2020) report a significantly lower rate of approximately 1 planet per star, though this model is limited to orbital periods shorter than 400 days.

As shown in Fig. 11, the distribution of planet sizes predicted by these models show some notable differences. The Hsu et al. (2019) model predicts a substantially higher fraction of planets with R_p < 1R_⊕ compared to NGPPS and Kunimoto & Matthews (2020) models, especially for periods between 60 and 500 days. On the other hand, NGPPS has a larger fraction of planets above R_p > 1 R_⊕, with a marked peak around R_p ~ 5 R_⊕. This primarily stems from NGPPS predicting a large number of cold Neptune-sized planets with orbital periods >100 days, which remain largely unobserved to date.

The Kunimoto & Matthews (2020) model aligns with the Hsu et al. (2019) predictions for planets R_p > 3 R_⊕ but shows far fewer low-mass planets. The disagreement between the three models for smaller planets stems from the fact that observational constraints on Earth-sized planet occurrence rates are still limited. The upcoming PLATO mission aims to address this gap, and will improve these estimates.

Translating these findings to the galactic context, the true nature of low-mass planet occurrence rates and their dependence on metallicity is an important factor in understanding exoplanet demographics throughout the galaxy. Massive planets are more comprehensively understood and our findings more robust. However, other systematic effects have the potential to alter the giant planet population that we have not addressed in this study. We discuss some of these limitations in the next section.

Fig. 11 Occurrence rates of planets around FGK stars within an orbital period of 500 days, as estimated by Matuszewski et al. (2023), based on the works of Hsu et al. (2019), Kunimoto & Matthews (2020) and Emsenhuber et al. (2021b), as a function of planet radius.

5.3 Model assumptions and limitations

5.3.1 Disc lifetimes

The lifetime of protoplanetary discs, primarily governed by external photoevaporation, plays a crucial role in planet formation. We adopt the range of photoevaporation rates from the Emsenhuber et al. (2021b) NGPPS sample, which are chosen to reproduce typical disc lifetimes in the solar neighbourhood. Efficient massive planet formation (>10 M_⊕) in the Bern model require disc lifetimes around 3–6 Myr, indicating that early disc dispersal is not a dominant mechanism with respect to halting giant formation (Schlecker et al. 2021). This is reflected in the weak correlation coefficient between photoevaporation rate and planet type (Fig. 2).

However, it is important to note that most studies of disc lifetimes, including those by Emsenhuber et al. (2021a) and Ansdell et al. (2017), have focussed on regions with relatively low stellar densities. In contrast, the galactic centre presents a drastically different environment, where higher stellar densities lead to more intense photoevaporation and, consequently, shorter disc lifetimes (Winter & Haworth 2022). The Central Molecular Zone (CMZ) in particular, exhibits gas surface density of ~1000 M_⊙ pc⁻² (Henshaw et al. 2016), exceeding that of the solar neighbourhood by a factor of 100 (Kruijssen 2012), which is connected to larger and denser stellar clusters (Farias et al. 2024). This heightened density can significantly amplify photo-evaporation, and thus reduce disc lifetime. Winter et al. (2020) constructed an analytical model for the disc lifetime as a function of galactic environment, and find that virtually all discs in the CMZ have lifetimes <3 Myr, severely hindering giant planet formation.

While the CMZ is an exceptional region within the MW, with a diameter of approximately 0.6 kpc, it is not sufficiently resolved in our simulations. Nonetheless, it is plausible that the shortened disc lifetimes in high-density regions towards the galactic centre could lower the occurrence rates of giant planets.

5.3.2 Initial embryo distribution

The formation of planet embryos is not explicitly modeled in the Bern model. Instead, a fixed number of embryos are distributed logarithmically at the beginning of each run. However, an analysis of the connection between initial parameter and resulting planet population for the NGPPS population performed by Schlecker et al. (2021) found the location of a planet embryo is the single strongest predictor for the final planetary type and this distribution might vary significantly across different galactic environments.

Voelkel et al. (2021) have studied the connection between planetesimal formation and planet embryo formation using a one-dimensional evolution model coupled to the N-body code LIPAD. They find that the mean orbital separation of planet embryos tends towards a constant value when expressed in terms of the embryo Hill radii. This value (−10–20 R_Hill) is largely insensitive to planetesimal disc mass and density profile, and leads directly to a logarithmic distribution of embryos. This is consistent with previous studies (Kokubo & Ida 1998), and with assumptions made for the NGPPS runs.

In terms of the number of planet embryos, they find that the total number of embryos increases with increasing planetesimal disc mass and density slope. For low solid disc masses, specifically 6 M_⊕ between 0.5 and 5 AU, around 20–40 embryos are formed, depending on the assumed density profile. At a higher planetesimal disc mass of 27M_⊕, around 60 embryos are formed, irrespective of the assumed profile. Translated into the context of our planet population model, which assumes universal distribution of gas disc masses the galaxy, these findings imply that stars in metal-rich regions (boasting higher solid disc masses) would, on average, form a greater number of planet embryos. Since the occurrence rates of low-mass planets in the model increases with the initial number of embryos (Table 6), this would indicate that Earth-like planets and Super-Earths are more common in metal-rich regions. The frequencies of giant planets on the other hand are insensitive to NEmbryo and are not strongly affected.

5.4 Relation to previous works

Few studies have been dedicated to planet populations in the galactic context, primarily due to the difficulty in observing distant planets. Despite these challenges, some statistical studies have been performed to study the variations of these populations.

Bashi & Zucker (2022) analysed a sample of 506 planet candidates orbiting around 369 FGK host stars, categorising them into thin- and thick-disc components. They found that the occurrence rate of close-in super-Earths in the thin disc is approximately 1.3–3.6 larger than in the thick disc, albeit with low number statistics as only 8 of their planets are associated with thick disc stars. This is consistent with our findings, where the lower average [Fe/H] in the thick disc leads to an increased occurrence rate of Earth-like planets, with a corresponding drop in the frequency of super-Earths by a factor of 1.5. Bashi & Zucker (2022), however, attribute this variation primarily to the dynamical history of the stars, including potential disruptive encounters, and gravitational perturbations among planets in typically older systems, rather than differences in metallicity.

Koshimoto et al. (2021) estimated the planet-hosting probability of stars as a function of galactic radius by analysing gravitational microlensing events. Unlike transit and radial velocity detections, which are typically limited to within 1 kpc from the Sun, microlensing can probe much larger distances. Their analysis of 28 planetary microlensing events has shown no strong dependence of occurrence rates on galactic radius, consistent with our findings. However, larger microlensing surveys are needed to study statistical differences between the bulge and disc populations. The Nancy Grace Roman Space Telescope (Spergel et al. 2015) is expected to detect more than 50 000 microlensing events and is expected to find around 1400 planets in the galactic bulge (Penny et al. 2019), allowing for a robust estimation of planet occurrence rates around bulge stars.

Nielsen et al. (2023) constructed an analytical planet formation model to study the occurrence rates of planets around alpha-poor (broadly representing the thin disc), alpha-rich (broadly representing the thick disc), and halo stars. Their findings show minimal dependence of Earth-like planets on the galactic component, a slight increase in the occurrence of Super-Earths, and a strong dependence of giant planets on these galactic components. Similarly to our finding, this is primarily attributed to the metallicity-dependence of the occurrence rates. Earth-like planet frequency in their model is largely independent of metallicity, while giant planets strongly favour higher metallicities. Their overall findings are therefore similar to ours.

Furthermore, the works of Bashi & Zucker (2019), Bashi et al. (2020), Chen et al. (2022), and Nielsen et al. (2023) have expanded the scope of study to include the [α/Fe] fraction in addition to the [Fe/H] fraction of the host star and its implications for planet formation. This aspect is potentially important, as thick disc stars are α-enhanced compared to thin disc. We have not included this quantity, as the NGPPS populations have no direct α element dependence. The findings of Nielsen et al. (2023) however suggest that the α-enhancement of the host stars plays a secondary role compared to the overall metallicity, as traced by [Fe/H]. Yet, Swastik et al. (2022) find that planet masses increase marginally with decreasing [α/Fe].

6 Summary and conclusions

In this study, we combined a simulated MW analogue for HESTIA suite of high-resolution galaxy formation simulations (Libeskind et al. 2020) and the NGPPS dataset from the Bern model for planet formation and evolution (Emsenhuber et al. 2021a) to study the planet demographics throughout the galaxy. Using the stellar ages and metallicities provided by HESTIA, we were able to map the occurrence rate of various planet types as a function of galactic environment and quantify variations in the frequency of Earth-like planets, Super-Earths, Neptunians, and giant planets (defined through their planetary masses). Our key findings are as follows:

• 1.

  The occurrence rates of specific planet types within the framework of the NGPPS population model depend primarily on the available solid mass in the planetesimal disc, determined by the initial gas disc mass and host star metallicity. For a fixed initial gas mass, planet masses increase with increasing [Fe/H] (Sec. 2.3);

• 2.

  Based on this, the overall planet populations in the metal-rich galactic bulge and thin disc differ significantly from the populations found in the metal-poor thick disc and stellar halo. Planets with masses >300 M_⊕ are 10–20 times more common around stars in the thick disc than compared to the thin disc. Earth-mass planets around Sun-like stars are most abundant in the thick disc, with occurrence rates 1.4 times higher than in the thin disc (Sec. 4.1);

• 3.

  The planet populations around Sun-like stars are youngest in the thin disc, with Earth-like planets having an average age of 6.1 Gyr in the simulated MW analogue. The oldest Earth-like planets are found in the halo, at an average age of 8.3 Gyr. The average age decreases with increasing planet mass (Sec. 4.1);

• 4.

  The planet populations in different galactic environments vary more strongly around Sun-like stars than those around lower mass stars (M_* = 0.3–0.5 M_⊙). This is caused by stronger correlation between metallicity and occurrence rate for massive stars in the NGPPS sample. Earth-like planets are most common around intermediate mass dwarfs with masses 0.5–0.7 M_⊙ (Sects. 2.3 and 4.3).

The study of planet formation as a function of galactic environment is still in its infancy. Currently, transit and RV surveys are restricted mostly to stars within 1 kpc of the Sun. Future exoplanet surveys will however increase our understanding of planet formation and the effects of the environment. The upcoming PLATO mission (Rauer et al. 2014) will observe around 200 000 stars covering a large area of the sky and is likely to detect a sample of exoplanets around thick-disc stars. In addition, the Nancy Grace Roman Space Telescope (Spergel et al. 2015) is expected to find a statistical sample of around 1400 planets in the galactic bulge (Penny et al. 2019) via microlensing. Understanding the exoplanet demographics in different environments will yield new insights into processes involved in planet formation and evolution, while helping improve planet formation models.

Acknowledgements

C.B. thanks the Young Academy Groningen for their generous support through an interdisciplinary PhD fellowship. P.D. and M.T. acknowledge support from the NWO grant 016.VIDI.189.162 (“ODIN”). P.D. warmly thanks the European Commission’s and University of Groningen’s CO-FUND Rosalind Franklin program. K.R. is grateful for support from the UK STFC via grant ST/V000594/1. We want to thank Alexandre Emsenhuber for providing us with the Monte Carlo variable dataset used for the original NGPPS sample, and the Data & Analysis Center for Exoplanets (dace.unige.ch) for a seamless access to the NGPPS and observational exoplanet data. Access to the HESTIA data (Libeskind et al. 2020) was kindly provided by the CLUES collaboration (https://www.clues-project.org), and can be requested from them. The tables and further data presented in this work are available upon request to boettner@astro.rug.nl. We also would like to thank the referee for insightful feedback, which helped to improve this work. Software: The analysis has been performed using Python, and the scientific data analysis libraries numpy (Harris et al. 2020), pandas (The pandas development team 2023), statsmodels (Seabold & Perktold 2010) and scikit-learn (Pedregosa et al. 2011). The HESTIA simulation data was accessed and processed using yt (Smith 2015). Visualisation were made using matplotlib (Caswell et al. 2023), seaborn (Waskom 2021) and yt. The galaxy decomposition was performed using MORDOR (Zana et al. 2022).

Appendix A Additional material

Fig. A.1 Predicted number of planets by type from the planet model (M* = 1M⊙, NEmbryo = 50), based on initial gas disc mass and [Fe/H]. The white diagonal lines represent contours of equal solid disc mass, measured in Jupiter masses, showing a pronounced correlation between solid disc mass and planet count. The histograms on top illustrate the marginal distribution for [Fe/H], while dots mark the original Monte Carlo variable sample from Emsenhuber et al. (2021b).

Fig. A.2 Normalised distributions of age and metallicity distribution of the galaxy components.

References

All Tables

All Figures

Fig. 1 Distribution of Monte Carlo variables for original sample from Emsenhuber et al. (2021b), and effect on number of giant planets for the M* = 1 and NEmbryos = 50 model run. The colour marks the number of planets in the system after 20 Myr of evolution. Shown are the correlation between the gas mass, Mg, metallicity, [Fe/H], inner protoplanetary disc radius, rin, and external photo-evaporation rate, ${\dot{M}}_{\text{wind }}$.

In the text

	Fig. 2 Kendall correlation matrix between Monte Carlo variables and planet types for the M* = 1 and NEmbryos = 50 run (Emsenhuber et al. 2021b).
In the text

	Fig. 3 Variations of planet occurrence rates per category given different numbers of initial embryos, M* = 1 (Emsenhuber et al. 2021b).
In the text

	Fig. 4 Variations of planet occurrence rates per category given different host star masses, NEmbryos = 50 (Burn et al. 2021).
In the text

	Fig. 5 Side and face-on maps of the distribution of stars for the different galaxy components.
In the text

	Fig. 6 Occurrence rates of different planet types in the four galaxy components. The occurrence rates are calculated for the M* = 1M⊙ and NEmbryo = 50 run.
In the text

	Fig. 7 Spatial density of planets as a function of distance from the galactic centre (normalised to R = 0). The points correspond to the estimated values per radial bin, while lines are smoothing splines to show trends. The grey, dash-dotted line corresponds to the normalized density of stars.
In the text

	Fig. 8 Maps of distribution of planets for each planet Type in the HESTIA MW analogue for stars with M* = 1 M⊙. The assignment is based on the NEmbryo = 50 NGPPS population.
In the text

	Fig. 9 Ratio between the occurrence rates in the fiducial metallicity model (flat extrapolation) and the alternative case where we assign zero planets. The ratios are calculated for M* = 1 M⊙ and NEmbryo = 50.
In the text

	Fig. 10 Occurrence rates of planet around host star with masses of M* = 0.5 M⊙ and M* = 0.3 M⊙, for the NEmbryo = 50 run.
In the text

	Fig. 11 Occurrence rates of planets around FGK stars within an orbital period of 500 days, as estimated by Matuszewski et al. (2023), based on the works of Hsu et al. (2019), Kunimoto & Matthews (2020) and Emsenhuber et al. (2021b), as a function of planet radius.
In the text

Fig. A.1 Predicted number of planets by type from the planet model (M* = 1M⊙, NEmbryo = 50), based on initial gas disc mass and [Fe/H]. The white diagonal lines represent contours of equal solid disc mass, measured in Jupiter masses, showing a pronounced correlation between solid disc mass and planet count. The histograms on top illustrate the marginal distribution for [Fe/H], while dots mark the original Monte Carlo variable sample from Emsenhuber et al. (2021b).

In the text

	Fig. A.2 Normalised distributions of age and metallicity distribution of the galaxy components.
In the text