Action calculation and analysis

Once the model of the galaxy is appropriately centred and rotated, we utilise AGAMA (Vasiliev, Reference Vasiliev2019), a publicly available software library for a broad range of applications in the field of stellar dynamics. A smooth approximation of the galactic potential is generated and used to calculate three standard action coordinates for star particles. These are the vertical action $J_z$ , the radial action $J_R$ , and the azimuthal action $J_\phi$ . Combined, they fully describe the orbit of the star particle: $J_R$ and $J_z$ describe oscillations in the radial and vertical directions, respectively, while the azimuthal action represents the z-component of the angular momentum $L_z$ in an axisymmetric potential. We will adopt the notation of $L_z$ for the azimuthal action for the remainder of this work.

Additionally, the angular momentum of the particle $L_\mathrm{c}(E)$ is calculated with total energy E, for a circular orbit (e.g. Abadi et al., Reference Abadi, Navarro, Steinmetz and Eke2003). Then, the circularity of the orbit can be calculated as $\xi = L_z/L_\mathrm{c}(E)$ , for each star particle. Circularity, defined this way, is just another way of assessing the eccentricity of the orbit. For example, particles with $\xi \geq 0.9$ are on nearly circular orbits, while the ones with lower circularity values have more eccentric orbits.

Figure 5 shows the distributions of these parameters – actions, circularity, and galactocentric distance, at two different times (initial and final). These distributions do not change significantly throughout the simulation, which is expected in a stable axisymmetric model of the Milky Way. Minor disk thickening is noticeable on distributions that include vertical action $J_z$ . The apparent immutability of distributions, a sign of the stability of the system as a whole, does not imply that individual stars do not migrate. It merely indicates that most migrations are balanced out on a larger scale. As the migrations of individual star particles are of the most interest for this work, we will focus on them once we state key observations from Figure 5. The majority of particles are in highly circular orbits, as is both intuitively expected and observationally grounded, as far as the Milky Way stellar population is concerned (e.g. Cubarsi, Stojanović, & Ninković, Reference Cubarsi, Stojanović and Ninković2021). The ones whose orbits are significantly eccentric (with lower circularity values) inhabit, almost exclusively, inner to central regions of the galaxy. This includes a broad Solar neighbourhood, in its entirety. Hence, even without looking into individual migrations, it is expected that the motions of stars affect the boundaries of the GHZ, at least to some extent. Stars that are not on nearly circular orbits (stars with circularities $\xi<0.9$ , e.g. Beraldo e Silva et al., Reference Beraldo e Silva, Debattista, Nidever, Amarante and Garver2021), account for $39.5\%$ of our stars subsystem, and the majority of them are confined within inner $9.65$ kpc initially (at $t=0.5$ Gyr), while the values change to $45.3\%$ and $10.17$ kpc by the end of the simulation (at $t=10$ Gyr). These stars, as expected, have lower values of angular momentum $L_z$ , while both their radial $J_R$ and vertical $J_z$ action span over a broad range. Interestingly, the $J_R-J_z$ distribution in Figure 5, shows explicit anti-correlation: stars with higher radial action have lower values of vertical counterpart and vice versa. This means that stars whose orbits deviate from nearly circular tend to oscillate either in a radial or in a vertical direction, but not simultaneously in both of these directions if the oscillations are drastic.

For individual stars, we calculate the absolute changes of relevant parameters (galactocentric distance R, circularity $\xi$ , and actions $J_R$ , $L_z$ , and $J_z$ ) as:

(1) \begin{equation} \Delta X = X(t=10\;\mathrm{Gyr})-X(t=0.5\;\mathrm{Gyr})\end{equation}

where X is any of the parameters. We show two-dimensional probability density distributions $R\;(t=0.5\;\mathrm{Gyr})-\Delta X$ , where $R\;(t=0.5\;\mathrm{Gyr})$ is an initial galactocentric distance, in Figure 6 to determine the magnitude of these changes and which regions of the disk are the most prone to them.

Figure 4. Evolution of the global properties of the stars subsystem, top to bottom: median galactocentric distance $\langle R \rangle$ , (arithmetic) mean height $z_\mathrm{rms}$ and velocity dispersion in cylindrical coordinates: radial $\sigma_{v_R}$ , circular $\sigma_{v_\phi}$ , and vertical $\sigma_{v_z}$ . Different simulations are represented with different line colours and styles, as indicated by the legend.

Angular momentum change, $\Delta L_z$ , is the only considered parameter whose distribution appears almost perfectly symmetrical around zero. In total, $51.2\%$ of stars satisfy $\Delta L_z<0$ , i.e., they experience angular momentum loss (it is implied that the rest have angular momentum gain). The most extreme angular momentum absolute changes (of any sign) correspond to stars initially located in the central region. Similarly, the most extreme absolute changes in the radial action $J_R$ correspond to stars initially located in inner to central regions, although the probability density distribution is not as symmetric around zero, and more stars experience positive change (with only $37.8\%$ of the sample satisfying $\Delta J_R<0$ ). The majority of stars are evolving towards orbits that oscillate more radially. On the contrary, both vertical action $J_z$ and circularity $\xi$ have the most extreme absolute changes in the innermost regions, and both probability density distributions are asymmetrical around zero. Across all galactocentric distances, vertical action predominantly increases (in total, $67.3\%$ of stars satisfy $\Delta J_z>0$ ), which is in line with previously discussed disk thickening. Possibly as a consequence, circularity predominantly decreases (in total, 68% of stars satisfy $\Delta \xi<0$ ). However, at higher distances, the absolute changes in circularity are minor, indicating that most stars that inhabit the outer regions do not have significant deviations from their initial circularity or orbital eccentricity. Interestingly, there are stars whose circularity increases. However small a fraction it is, this sub-sample of stars appears to be the most interesting from the habitability point of view: these stars evolve towards circular orbits, which means that their galactic environment conditions are likely to be more stable in the long term. We will consider this sub-sample in particular when analysing and discussing habitability.

Finally, the parameter of the most interest, the galactocentric distance, has the probability density distribution of its absolute change asymmetrical around zero. In total, about 46% of stars migrate inwards (i.e. satisfy $\Delta R<0$ ), while the rest migrate outwards (i.e. $\Delta R>0$ ), resulting in the net radial migration outwards of about 8%Footnote ^c. This result is in line with the report of, e.g., Roškar et al. (Reference Roškar, Debattista, Quinn, Stinson and Wadsley2008, Reference Roškar, Debattista, Quinn and Wadsley2012), who noticed that the outward radial migration is larger than its inward counterpart and our previous observation that the minor global radial migration outward is expected. The asymmetry of this probability density distribution is insightful. Interestingly, stars with the most extreme outward radial migrations are initially clustered around $R\simeq 4.44\;\mathrm{kpc}$ , while the ones with significant inward migrations originate from higher galactocentric distances, $R\simeq 6.84\;\mathrm{kpc}$ . This hints at a sort of radial mixing in the central disk parts and should, in principle, affect the traditionally constrained GHZ.

It is important to highlight that the radial mixing is typically associated with spiral patterns (e.g. Sellwood & Binney, Reference Sellwood and Binney2002) and that the previous studies, mentioned earlier, include transient spiral structure. The general agreement of our results implies that the spiral patterns are not a necessary condition for radial migrations and mixing and that similar outcomes can happen in axisymmetric models with sufficient, mild levels of dynamical heating. However, the inclusion of a transient spiral structure in our model would, most likely, result in slower radial migration outwards or a lower net outwards trend at the expense of even more efficient radial mixing in the central parts of the disk.

As mentioned in the introductory section, the aim is to separate different source causes of stellar radial migrations. For this purpose, it is investigated if there is a correlation between the relative changes in actions and the radial changes and of what significance and strength. More specifically, the Spearman’s correlation coefficient $\rho$ is calculated between $\Delta X/X_0$ and $\Delta R$ , where X corresponds to any of the three actions ( $J_R$ , $L_z$ or $J_z$ ), $\Delta X$ is defined with Equation 1 and $X_0$ represents appropriate initial value at $t=0.5\;\mathrm{Gyr}$ . The results of this test are listed in Table 2. Since the test is sensitive to outliers, a robust linear regression is performed (thus, de-weighting outliers) with 95% confidence interval. No difference is found, and both tests lead to essentially the same conclusions.

Figure 5. Probability density distributions, R is galactocentric distance, $J_R$ , $L_z$ , and $J_z$ are radial, azimuthal, and vertical actions, respectively, and $\xi$ is orbital circularity, for two different times represented with different line colours, as indicated by the legend.

Figure 6. Two-dimensional probability density distributions where x-axis, $R\;(t=0.5\;\mathrm{Gyr})$ , is the initial galactocentric distance and y-axis corresponds to the absolute change of parameters defined and presented in Figure 5. The absolute change is calculated as: $\Delta X = X(t=10\;\mathrm{Gyr})-X(t=0.5\;\mathrm{Gyr})$ , where X is any of the parameters.

Relative changes in angular momentum and the radial action are correlated with the radial change, but the relative change in the vertical action appears independent. In particular, the correlation is moderate for angular momentum and very weak for the radial action. This is expected and corroborates the conclusion of Frankel et al. (Reference Frankel, Sanders, Ting and Rix2020) that diffusion in angular momentum dominates in the secular orbit evolution and that, in comparison, radial heating is much weaker.

Vertical motions appear stochastic and show no explicit correlation with radial migrations. They might be irrelevant to the present study, as its focus is on the radial boundaries of the GHZ. This is manifestly not the case, however, when the habitability of disk galaxies is considered on a larger scale. Vertical motions directly affect the number of galactic plane crossings, which can be hazardous as the probability of close encounters peaks in mid-plane (e.g. Medvedev & Melott, Reference Medvedev and Melott2007; Bailer-Jones, Reference Bailer-Jones2009; Melott et al., Reference Melott, Bambach, Petersen and McArthur2012; Sloan, Alves Batista, & Loeb, Reference Sloan, Batista and Loeb2017), as do tidal stresses associated with galactic disk. Moreover, stellar orbits with higher vertical oscillations contribute to these stellar systems experiencing vastly different environments, which might not satisfy, in terms of stability, continuous habitability conditions. Any attempt to unify various habitability conditions and constraints and consider multiple factors towards a complex galactic habitability model should, thus, include vertical motions in some way. In such a case, it is essential to separate kinematically distinct (among other, non-dynamical differences) thin and thick disk components (see, e.g. Vieira et al., Reference Vieira, Carraro, Korchagin, Lutsenko, Girard and van Altena2022, and references therein), as their highest impact on the results should be related to the vertical motions. The single-disk model from this work, which should not be solely considered as representative of a thin disk, is suitable for this study, which aims to make a first step towards exploring and quantifying dynamical-related effects on the radial boundaries of the GHZ. However, when vertical motions are taken into consideration, the single-disk model becomes less applicable as multiple-component disk models should uncover insightful trends.

Comparison with LRMs

For the sake of completeness, as previously mentioned, radial migration trends are briefly analysed for all models and compared with the main results. As expected, in the LRM_UES simulation, radial migration trends are roughly the same, with a net outwards trend of roughly 8%, and the asymmetry of $R-\Delta R$ distribution around zero being the same. This could confirm the previous remark that LRMs could still give reliable results if numerical parameters, such as softening lengths of different particles, are carefully chosen.

On the contrary, simulation LRM considerably differs from the main results. The whole disk expands significantly in the vertical direction, and a net outwards trend of radial migrations is slightly higher (around 10%). The absolute magnitudes of radial migrations, i.e. $|\Delta R|$ values, are typically larger than in our main simulation, suggesting that the bar might be able to induce migrations of a longer range. Stars in the inner region of the disk, where the artificially formed bar is located, predominantly migrated inwards (around 62% of stars satisfy $\Delta R < 0$ ). Median magnitudes of radial migrations in this region are also angle-dependent: the largest migrations outwards are aligned with the bar’s major axis, whereas their inwards counterparts align with the minor axis. Unsurprisingly, as the stars captured by the bar evolve towards radial orbits, the negative changes of circularity are, on average, higher in magnitude than in our main simulation and aligned with the bar’s major axis. In the parts of the disk towards the central region, however, we did not notice signs of efficient radial mixing, as the $R-\Delta R$ distribution appears roughly symmetrical around zero (contrary to LRM_UES and HRM simulations).

While the results of the LRM simulation are not considered to be reliable, due to numerical effects, there is a striking agreement between the brief results presented here and the previous more rigorous works on the topic (see, e.g. Di Matteo et al., Reference Di Matteo, Haywood, Combes, Semelin and Snaith2013; Filion et al., Reference Filion, McClure, Weinberg, D’Onghia and Daniel2023, and references therein). The lack of efficient radial mixing in the presence of a bar should have been expected. As Di Matteo et al. (Reference Di Matteo, Haywood, Combes, Semelin and Snaith2013) point out, in non-axisymmetric galaxy models, migrations and mixing do not occur at the same time, and radial mixing can only be established when the phase of significant radial migrations, caused by the bar, is over. Angle-dependent trends of radial migrations and the inwards trend in the inner disk region are also in agreement with the results of Filion et al. (Reference Filion, McClure, Weinberg, D’Onghia and Daniel2023). However, despite these similarities, we cannot use this model in our habitability considerations due to the previously mentioned numerical effects that contaminate the results. For the results to be reliable and robust enough, our model has to have numerical effects as low as possible and the bar formation should arise from an actual (realistic) instability.

Table 2. Spearman’s correlation coefficient $\rho$ between relative change $\Delta X/X_0$ and absolute change $\Delta R$ in the galactocentric distance, where X corresponds to any of the three actions, $\Delta X$ is defined with Equation 1 and $X_0$ represents initial value at $t=0.5\;\mathrm{Gyr}$ .

The rate of close encounters

Aside from the metallicity condition, which defines the availability of material to form complex forms of life, another critical condition for habitability is the continuity of habitable-friendly conditions (Vukotić et al., Reference Vukotić, Steinhauser, Martinez-Aviles, Ćirković, Micic and Schindler2016, for example, call this notion “habitable time”). From a strictly dynamical point of view, this continuity will be disrupted by the close stellar encounters, which are catastrophic enough to change planetary orbits initially in circumstellar habitable zones. The rate of such encounters $\Gamma$ is given by:

(2) \begin{equation}\Gamma = n_\star \langle v \rangle \langle \sigma \rangle\end{equation}

where $n_\star$ is stellar number density, $\langle v \rangle$ relative stellar velocity, and $\langle \sigma \rangle$ cross-section of the encounter having adverse astrobiological consequences. Evidently, we can directly calculate $\langle v \rangle$ , and we will adopt $\langle \sigma \rangle \sim 100\;\mathrm{AU}^2$ (as suggested by, e.g. Laughlin & Adams, Reference Laughlin and Adams2000). Stellar number density estimate is particularly challenging since we can only calculate mass density $\rho_\star$ in our galaxy model and do not have information on the number of stars. For the sake of simplicity, we will assume that all stars have Solar masses, i.e., that $n_\star = \rho_\star/\mathrm{M}_\odot$ . Since our galaxy model is axisymmetric, it is safe to calculate the rate of close encounters $\Gamma$ as a function of galactocentric distance R, which we show in Figure 7.

Figure 7. The rate of close encounters $\Gamma$ , defined with Equation 2, as a function of galactocentric distance R.

Represented this way, the rate of close encounters is a declining function of galactocentric distance. This means that imposing a certain threshold to limit the boundaries of the GHZ is only possible for the inner one. The outer boundary is, rather inconveniently for this study, limited by the metallicity. Generally speaking, one can argue that the outer boundary of the GHZ is also limited by the environment where the galaxy resides (e.g. part of a cluster or group, having satellite galaxies, etc.), given the outside-in nature of tidal effects. High-resolution cosmological or zoom-in simulations should be the best currently available tool for studying these environmental effects and limiting the outer boundary of the GHZ based on those considerations, in addition to metallicity.

To constrain the inner boundary of the GHZ, a critical value for the rate of close encounters needs to be assumed, for example, $\Gamma_\mathrm{crit} = 1 \;(\mathrm{age\;of\;the\;Earth})^{-1} \simeq 0.22\;\mathrm{Gyr}^{-1}$ . The calculated rate of encounters is smaller than this critical value over the entire range, which implies that the inner GHZ boundary cannot be constrained in this manner. Apart from the relevance of the inner boundary determination, the precision of the estimate can also be scrutinised. We will address the latter part first. For the Solar system, $\Gamma (\mathrm{R_\odot}) \sim 10^{-6}\;\mathrm{Gyr}^{-1}$ , but Sloan, Alves Batista, & Loeb (Reference Sloan, Batista and Loeb2017) give an estimate of $3 \times 10^{-8}\;\mathrm{Gyr}^{-1}$ for this value (thus, even smaller). While rigorously studying the effects of stellar encounters on habitability, Bojnordi Arbab & Rahvar (Reference Bojnordi Arbab and Rahvar2021) found the expected number of threatening stellar encounters in the Solar neighbourhood is $\Gamma \sim 10^{-4}\;\mathrm{Gyr}^{-1}$ . Thus, the estimate from this work, falls exactly in the middle and is in line with previous findings. Despite this agreement, our estimate might be too crude for the inner parts of the disk. More specifically, our assumption for the stellar number density is most likely incorrect in these regions, as the stars may have, on average, lower mass than the assumed Solar mass. Accounting for this may increase the rate of close encounters by an order of magnitude in the central region, which is still not enough of an increase to constrain the GHZ outside the central galactic region. Another argument that can be made is that our adopted value for cross-section is too conservative and that higher values should be used (see, e.g. Li & Adams, Reference Li and Adams2015; Brown & Rein, Reference Brown and Rein2022, for cross-section discussion). The closest encounters of the Solar system within the present $10^7$ yr are found to occur at a distance of $10^4$ AU from Gaia Data Release 3 (Bailer-Jones, Reference Bailer-Jones2022). As adopting this value for the cross-section in our calculation might imply an overly restrictive view of habitability, we will consider its implications in relation to previously presented results. Since the cross-section is de facto a constant, adopting $\langle \sigma \rangle \sim 10^4\;\mathrm{AU}^2$ instead would not change the shape of the declining $\Gamma (R)$ function shown in Figure 7, but rather shift it upwards by two orders of magnitude. For the Solar system, we would then get $\Gamma (\mathrm{R_\odot}) \sim 10^{-4}\;\mathrm{Gyr}^{-1}$ , exactly the value given by Bojnordi Arbab & Rahvar (Reference Bojnordi Arbab and Rahvar2021), and the inner boundary of the GHZ could be placed at roughly $R=1\;\mathrm{kpc}$ .

Bojnordi Arbab & Rahvar (Reference Bojnordi Arbab and Rahvar2021) argue that the rate of disruptive stellar encounters is the most dependent on the stellar number density. This should be intuitively understandable, as discussed by Stojković, Vukotić, & Ćirković (Reference Stojković, Vukotić and Ćirković2019a), just by examining Equation 2: cross-section should be treated as a constant, leaving us with two variables, one of which, relative stellar velocity $\langle v \rangle$ , does not vary as much. Naturally, this leads to conclusion that the rate of close encounters very well may be irrelevant, or at least redundant, on larger scales (i.e., for constraining the GHZ) while still applicable on smaller scales (e.g., when considering substructures and local over-densities).

Despite not being of only dynamical significance, the Galactic gas reservoir with its redistribution and its state, should also be considered. In particular, close encounters with gaseous local over-densities such as giant molecular clouds (GMCs), regions of higher star-formation primarily located in spiral arms (Hou, Han, & Shi, Reference Hou, Han and Shi2009), can have severe effects on the orbits of stars (e.g. Fujimoto, Inutsuka, & Baba, Reference Fujimoto, Inutsuka and Baba2023), and moving through GMCs can even be hazardous for the planetary systems orbiting those stars (Kokaia & Davies, Reference Kokaia and Davies2019). Undoubtedly, as previous works suggest, gas can affect the motions of stars but only if local over-densities are present, or rather specific distribution of gas and its flows, typically related to non-axisymmetric features such as bars or spirals (e.g. Berentzen et al., Reference Berentzen, Shlosman, Martinez-Valpuesta and Heller2007; Yu et al., Reference Yu, Xu, Ho, Wang and Kao2022). In our axisymmetric galaxy model, we have little to no reason to assume that the inclusion of gas would affect radial migrations or mixing, as the distribution of gas would also be axisymmetric and uniform (unless we specifically model it differently, which is beyond the scope of this paper). However, in a more complex galaxy model, which includes one or more non-axisymmetric features, the possible effects of gas should be modelled or, at least, carefully discussed.

Stars evolving from non-circular to circular orbits

As previously mentioned, stars whose orbital circularity increases are the most interesting from an astrobiological point of view, as these stars evolve towards circular orbits making their environmental conditions stable. Once they evolve to a circular orbit, such stars are likely to permanently reside at a given galactocentric radii making a stable population that can be used to constrain the confinements of the GHZ. We will consider the sub-sample of stars whose circularity increases significantly enough (i.e. from non-circular to nearly circular orbits), regardless of galactocentric distance, since as earlier discussed the GHZ boundaries cannot be effectively constrained with the rate of close encounters. As the previous analysis has shown, these stars (initially in non-circular orbits) exclusively inhabit the inner to central regions of the disk, and this sub-sample represents merely a $3.8\%$ of the stars in the whole sample. The best candidates from this sub-sample should be the stars that stabilise on nearly circular orbits after they have migrated outwards. Thus, we show the initial versus final galactocentric distance of the stars from this sub-sample with a $y=x$ dashed black line to facilitate visual separation of inward from outward migrators in Figure 8. Additionally, we show a linear regression line with 95% confidence interval. It is clear that the regression line is not a perfect fit, nor is linear or any other particular functional relationship expected, as the most stars are scattered around. However, it is included for the purpose of an easier assessment of whether there is a global trend and how it deviates from the dashed black line separating migrators.

Figure 8. Initial versus final galactocentric distance of stars evolving from non-circular to nearly circular orbits, plotted with linear regression line with 95% confidence interval. The dashed black line corresponds to the $y=x$ line.

Figure 9. Probability density function of galactocentric distance R for stellar populations. The blue line represents the final distribution of the initial stellar population of the broad Solar neighbourhood, and the orange one represents the origin of its final stellar population (as indicated by the legend). The Solar neighbourhood, as well as percentages and directions of migrators, are clearly indicated.

Interestingly, there are both inward and outward migrators in this sub-sample, but outwards migrations are slightly more prevalent (with net outwards migrations of $11.5\%$ ) than in the general population (8%). However, at the smaller galactocentric distances, more stars tend to migrate outwards, while it is the opposite case for the larger ones (although there are fewer stars there falling into this sub-sample). This agrees with previously noticed asymmetry in $R-\Delta R$ distribution (Figure 6), indicative of efficient radial mixing in the region. Median galactocentric distances are $4.168\;\mathrm{kpc}$ and $4.358\;\mathrm{kpc}$ for the initial and the final ones, respectively. Mean values for both are slightly larger than the median, indicating that the distribution is skewed to the right (i.e., higher end). This implies (in a less obvious manner than by examining Figure 8 for example) that migrations outwards in this sub-sample are, on average, of a higher magnitude.

The most fascinating stars from this sub-sample are the ones scattered above the regression line. These stars migrate farther away before stabilising in nearly circular orbits in the Solar neighbourhood and beyond. Their region of origin can vary a lot, going as low as $R(t=0.5\;\mathrm{Gyr})<3\;\mathrm{kpc}$ . The result is particularly interesting given that the observations of the Galactic stellar populations support this, including the outward motion of the Solar System (Wielen, Fuchs, & Dettbarn, Reference Wielen, Fuchs and Dettbarn1996). Naturally, this highlights the need to explore the Solar neighbourhood, its stellar population, and the migrations of stars in detail.

Broad Solar neighbourhood

Given that the presented galaxy model mimics the Milky Way, the traditional, annular GHZ shape can be considered. The conservative estimate from Lineweaver, Fenner, & Gibson (Reference Lineweaver, Fenner and Gibson2004) place this region between $R=7\;\mathrm{kpc}$ and $R=9\;\mathrm{kpc}$ . This region is further referred to as a broad Solar neighbourhood and specifically, the orbits of its initial and final stellar population are explored. Figure 9 presents the probability density function of galactocentric distance R for stellar populations that, at some point, reside in the broad Solar neighbourhood. The blue line is the final distribution of the initial stellar population – the present-day distribution of stars born in the Solar neighbourhood, i.e., after radial migrations and mixing. The orange line represents the origin of the final stellar population – the distribution of birth radii of stars which represent the present-day stable population of the Solar neighbourhood, i.e., before migrations and mixing. Since the distributions are not perfectly symmetrical, we list appropriate basic properties (mainly median value, first $Q_1$ and third $Q_3$ quartile, and interquartile range $\mathrm{IQR} = Q_3 - Q_1$ ) in Table 3. In addition, we calculate the expected boundaries of these distributions: lower $\mathrm{MIN}= Q_1-1.5\;\mathrm{IQR}$ and upper $\mathrm{MAX}= Q_3+1.5\;\mathrm{IQR}$ , thus excluding outliers, i.e., extreme migrators.

Table 3. Properties of distributions represented in Figure 9: median (essentially second quartile $Q_2$ ), first $Q_1$ and third $Q_3$ quartile, and interquartile range $\mathrm{IQR} = Q_3 - Q_1$ . The expected boundaries calculated as $\mathrm{MIN}= Q_1-1.5\;\mathrm{IQR}$ and $\mathrm{MAX}= Q_3+1.5\;\mathrm{IQR}$ are also listed.

^a Migrations of initial Solar neighbourhood

^b Origin of final Solar neighbourhood

While we refer to this assumed GHZ as the broad Solar neighbourhood, it should be apparent that the zone is quite narrow when radial migrations are considered. It is also a crowded zone, with a majority of stars moving in and out (thus efficiently mixing in an even broader radial zone), and less than half of the total stellar population inhabits this region at all times. The final stellar population of this zone came primarily from the inner parts of the disk ( $38.5\%$ of the population, as opposed to $17.8\%$ that came from outer parts). This is in good agreement with previous studies. The simulations performed by Roškar et al. (Reference Roškar, Debattista, Quinn, Stinson and Wadsley2008) indicate that about half the stars in the Solar neighbourhood came primarily from the inner parts of the galaxy. Hayden et al. (Reference Hayden, Recio-Blanco, de Laverny, Mikolaitis, Guiglion, Hill and Gilmore2018), analysing the data from the Gaia-ESO Survey on metallicity and kinematics of stars from the Solar neighbourhood, also implied that it is likely that stars of higher metallicity have migrated to the Solar neighbourhood from the inner parts of the galactic disk. While studying the influence of satellites on radial mixing in the galactic disk, Bird, Kazantzidis, & Weinberg (Reference Bird, Kazantzidis and Weinberg2012) mention that 20% of stars that end up in the Solar annulus started at $R<6\;\mathrm{kpc}$ . Interestingly, this is similar to the results of this work: 25% of the present-day stable population of the Solar neighbourhood originate from smaller galactocentric radii, $R\leq6.33\;\mathrm{kpc}$ (parameter $Q_1$ , Table 3).

The initial stellar population of this zone also migrates out in both directions. It is important to highlight that inward migrations are more prevalent ( $30.8\%$ ) than outward counterparts ( $25.2\%$ ), opposite to the global trends. However, it is expected due to radial mixing, as shown in our previous analysis and observations, in particular because of the asymmetry in $R-\Delta R$ distribution (Figure 6). Not only do the stars born in this region have a higher possibility of inward migrations, but they are also more likely to migrate farther (compared to their outward counterparts). This brings us to the asymmetry of the two distributions around the initially assumed GHZ. Despite the global net positive radial migrations outwards, these distributions do not have a symmetric spread around the broad Solar neighbourhood – the spread inwards is of a higher magnitude. If we would assume some initial metallicity gradient, this would imply that the metal-poor stars, that are born in the Solar neighbourhood, predominantly migrate to the inner regions, and vice versa for the metal-rich stars, which effectively flattens the metallicity gradient (e.g. Roškar et al., Reference Roškar, Debattista, Quinn, Stinson and Wadsley2008; Minchev & Famaey, Reference Minchev and Famaey2010; Minchev et al., Reference Minchev, Famaey, Combes, Di Matteo, Mouhcine and Wozniak2011; Vincenzo & Kobayashi, Reference Vincenzo and Kobayashi2020).

Moreover, the zone of efficient radial mixing is much wider than the assumed GHZ, as should be evident from the mid-spread (i.e., IQR in Table 3) of two presented distributions, as well as from their expected boundaries. Hence, our results suggest that the region roughly between $R=3\;\mathrm{kpc}$ and $R=12\;\mathrm{kpc}$ represents the zone of radial mixing. Within such a region, it is possible to further constrain the limits of the GHZ based on other (e.g., metallicity) grounds. The stars, which at some point reside in the broad Solar neighbourhood, either originate from or migrate to this wider zone, making it challenging to distinguish habitable from non-habitable systems solely on this account if no other additional constraints are applied. Not all stars in stable orbits residing in the GHZ should be considered habitable systems. For example, some systems might come from the outer regions of the Disk and might not be sufficiently metal-rich to form Earth-like planets. Generally, it might be the case that the star was born in a too-metal-rich (or even too-metal-poor) environment, thus a non-friendly one, before migrating to a more stable orbit in the GHZ.

On the concept of a GHZ

While a typical earth-like planet in a local universe is more likely to be found in a spheroid-dominated galaxy (rather than a disk-dominated one), as Zackrisson et al. (Reference Zackrisson, Calissendorff, González, Benson, Johansen and Janson2016) reported, we will focus closely on previous studies that examined our own Galaxy, since our galaxy model represents the Milky Way. As already mentioned in the introductory section, Lineweaver, Fenner, & Gibson (Reference Lineweaver, Fenner and Gibson2004) constrained the GHZ to an annular zone between $7\;\mathrm{kpc}$ and $9\;\mathrm{kpc}$ that spreads outwards, while Prantzos (Reference Prantzos2008) argued that the GHZ annulus is more likely to expand, possibly making the whole disk suitable for hosting life. Interestingly, Gowanlock et al. (Reference Gowanlock, Patton and McConnell2011) reported that the most habitable regions are located in the inner disk ( $\sim2.5\;\mathrm{kpc}$ ) around the mid-plane (i.e., above or below). On the contrary, Vukotić et al. (Reference Vukotić, Steinhauser, Martinez-Aviles, Ćirković, Micic and Schindler2016) found the outskirts to be the most habitable, constraining the GHZ annulus between $10\;\mathrm{kpc}$ and $15\;\mathrm{kpc}$ . Spinelli et al. (Reference Spinelli, Ghirlanda, Haardt, Ghisellini and Scuderi2021) studied the history of the GHZ while considering the effects of gamma-ray bursts and supernovae and found that, indeed, the outskirts are more habitable during early times, but the GHZ shifts inwards, ending as an annulus between $2\;\mathrm{kpc}$ and $8\;\mathrm{kpc}$ at present. While studying the effects of interactions with GMCs, Kokaia & Davies (Reference Kokaia and Davies2019) gave two estimates of the GHZ annulus for the thin disk (between $5.8\;\mathrm{kpc}$ and $8.7\;\mathrm{kpc}$ ) and the thick disk (between $4.5\;\mathrm{kpc}$ and $7.7\;\mathrm{kpc}$ ). Taking only dynamical constraints into account, present results place the whole of the broad Solar neighbourhood into the GHZ annulus, extending it towards the inner galactocentric radii – a birthplace of our own planetary system.

All these different results suggest that habitability estimates are far too sensitive to initial assumptions, fine-tuning of the models, experiment design, and habitability determining factors explored (for a review of different factors influencing galactic habitability, see, e.g. Stojković, Vukotić, & Ćirković, Reference Stojković, Vukotić and Ćirković2019a). Considering the effect of stellar migrations on galactic habitability, this work has demonstrated that, in contrast to traditional GHZ studies, the boundaries of such an annulus are not likely to be firmly constrained. However, our results also did not imply that the whole disk is likely to be significantly populated with habitable systems, at least when considering close encounters and the secular evolution of stellar orbits.

Good agreement of our results, employing the axisymmetric galaxy model with mild levels of dynamical heating, with previous works that include models with transient spiral structure, should not come as a surprise. In particular, Brunetti, Chiappini, & Pfenniger (Reference Brunetti, Chiappini and Pfenniger2011) emphasise the importance of taking into account the bar while studying the impact of radial migration on the chemical evolution of the Milky Way. Similarly, Minchev et al. (Reference Minchev, Famaey, Quillen, Di Matteo, Combes, Vlajić, Erwin and Bland-Hawthorn2012) mention that the bars are the most effective drivers of radial migration. Filion et al. (Reference Filion, McClure, Weinberg, D’Onghia and Daniel2023) found that the stars in the outer zones likely originated at smaller radii and had their orbits evolve outwards, which has particularly significant implications for habitability. It suggests that, in the presence of bars, the enrichment of the galactic outskirts with potentially habitable stellar systems can happen at a faster pace. However, it should not be considered that the presence of spiral arms is negligible in a dynamical sense. Even from a strictly dynamical perspective, spiral arms are an important driver of radial mixing in the central region, which expands the GHZ (or blurs its boundaries). We also have to acknowledge that they represent local over-densities and influence the rate of close encounters, which would be higher than in the inter-arm regions. Thus, the complexity of robust habitability considerations is evident.

The next generation of numerical GHZ models should effectively consider the motion of the stellar component at the mass resolution of individual stars in order to constrain the parts of the Disk that are populated with stellar systems that are most likely to be habitable. However, running a full-scale version of such numerical models will require a substantial amount of computing resources. We have demonstrated that a test particle approach gives fiducial results and might be used as a path-finder for developing such large-scale projects. Thus, in the meantime, more reliable conclusions will be obtained with advanced models that include the effects of GMCs and quasi-periodical spiral-arm crossings on the galactic stellar orbits, as well as more nuanced habitability considerations. In addition, the effects of the galaxy-dwelling environment should also be accounted for. The presence of satellites and minor mergers predominantly affect the orbits of stars in the outskirts directly (Bird, Kazantzidis, & Weinberg, Reference Bird, Kazantzidis and Weinberg2012; Carr et al., Reference Carr, Johnston, Laporte and Ness2022). They also cause perturbations that lead to the formation of non-axisymmetric features and drive radial migrations and mixing indirectly. Long-lived non-axisymmetric structures can also form due to quick galaxy flybys (e.g. Pettitt & Wadsley, Reference Pettitt and Wadsley2018; Mitrašinović & Micic, Reference Mitrašinović and Micic2023), in which case the influence on the outskirts should be brief unless these interactions happen on inclined orbits and lead to the formation of warps (Kim et al., Reference Kim, Peirani, Kim, Ann, An and Yoon2014). Cosmological and zoom-in simulations are possibly the best tool for studying these environmental effects and their influence on habitability. They are already utilised to explore the habitability of the entire zoo of galaxies (Forgan et al., Reference Forgan, Dayal, Cockell and Libeskind2017; Stanway et al., Reference Stanway, Hoskin, Lane, Brown, Childs, Greis and Levan2018; Stojković et al., Reference Stojković, Vukotić, Martinović, Ćirković and Micic2019b), many of which display a diversity of shapes and substructures. However, we can also employ isolated N-body or hydrodynamical simulations and other similar methods, essentially utilising various tools at our disposal. Thus, we obtain richer, more detailed, and more provocative insight into the underlying physics of habitability.