Evaluating galactic habitability using high-resolution cosmological simulations of galaxy formation

Abstract

We present the first model that couples high-resolution simulations of the formation of local group galaxies with calculations of the galactic habitable zone (GHZ), a region of space which has sufficient metallicity to form terrestrial planets without being subject to hazardous radiation. These simulations allow us to make substantial progress in mapping out the asymmetric three-dimensional GHZ and its time evolution for the Milky Way (MW) and Triangulum (M33) galaxies, as opposed to works that generally assume an azimuthally symmetric GHZ. Applying typical habitability metrics to MW and M33, we find that while a large number of habitable planets exist as close as a few kiloparsecs from the galactic centre, the probability of individual planetary systems being habitable rises as one approaches the edge of the stellar disc. Tidal streams and satellite galaxies also appear to be fertile grounds for habitable planet formation. In short, we find that both galaxies arrive at similar GHZs by different evolutionary paths, as measured by the first and third quartiles of surviving biospheres. For the MW, this interquartile range begins as a narrow band at large radii, expanding to encompass much of the Galaxy at intermediate times before settling at a range of 2–13 kpc. In the case of M33, the opposite behaviour occurs – the initial and final interquartile ranges are quite similar, showing gradual evolution. This suggests that Galaxy assembly history strongly influences the time evolution of the GHZ, which will affect the relative time lag between biospheres in different galactic locations. We end by noting the caveats involved in such studies and demonstrate that high-resolution cosmological simulations will play a vital role in understanding habitability on galactic scales, provided that these simulations accurately resolve chemical evolution.

Introduction

A useful paradigm in astrobiology is ‘the habitable zone’, a preferred location or region in space for planets to produce a biosphere (Huang 1959; Hart 1979). Most commonly, this term is applied to individual star systems. Radiative transfer calculations of the atmospheric response of an Earthlike planet to insolation (Kasting et al. 1993; Kopparapu et al. 2013, 2014) allow the identification of planetary orbits that permit surface liquid water. This is the principal definition of ‘the habitable zone’, and we henceforth refer to it as the stellar habitable zone (SHZ).

As a consequence of stellar evolution, the SHZ is time-dependent. The increasing luminosity of ageing stars tends to push the habitable zone further outwards (Rushby et al. 2013) – we can therefore define a continuously habitable zone as the region which remains inside the SHZ for the entirety of the star's hydrogen fusion phase, usually referred to as the Main Sequence.

Binary systems have especially non-trivial SHZs, as the contribution to insolation from both sources, as well as perturbations of the planet's orbit, must be accounted for (Eggl et al. 2012; Forgan 2012, 2014; Haghighipour & Kaltenegger 2013; Kaltenegger & Haghighipour 2013; Mason et al. 2013; Cuntz 2014). The interaction of the two stars also influences the planet's high energy radiation environment (Johnstone et al. 2015; Zuluaga et al. 2015).

Exomoons may also provide sites for biospheres, and their SHZs are complicated even further by the effects of tidal heating (Forgan & Kipping 2013; Heller & Barnes 2015), frequent eclipses (Heller 2012) and by infrared radiation from the host planet (Heller & Barnes 2013; Forgan & Yotov 2014). These SHZs implicitly assume that magnetic fields can shield the moon from high energy radiation, which bears its own caveats (Heller & Zuluaga 2013).

Despite these complications and uncertainties, the SHZ remains one of the principal tools used by astrobiologists to assess the potential of exoplanet systems for producing habitable worlds. To date, 29 planets with physical radii less than 2.5 times the Earth's radius, with orbits inside their local SHZ ¹ . This is of course only the first step towards determining if these worlds are indeed habitable, as their atmospheric chemistry and tectonic activity will play a crucial role (see Guedel et al. 2014 for a review).

Somewhat less discussed is the habitable zone at galactic scales, as was first proposed by Gonzalez, Brownlee & Ward (2001). In principle, we should be able to identify galactic habitable zones (GHZs), regions of galaxies that are more amenable to planet formation and the growth of a biosphere, while being relatively unmolested by hazardous events such as nearby supernovae (SNe), gamma ray bursts, close stellar encounters, etc.

This implies a balance. Terrestrial planet formation requires the formation of chemical elements up to iron and beyond. Astronomers use the concept of metallicity to describe the fraction of elements with higher atomic number than helium, so in these terms a minimum metallicity must be generated locally for terrestrial planet formation to proceed. Metallicity generation occurs in the cores of stars, so metallicity generation requires star formation and stellar death to release the requisite elements into the interstellar medium, so that they might be incorporated into planets. Equally, stellar death often results in violent mechanical and radiative events such asSNe. The local supernova rate must be reasonably low to prevent large doses of hazardous radiation sterilizing terrestrial planets. In analogy with the SHZ, the GHZ does not guarantee that planets inside it are habitable, but merely indicates more favourable conditions for habitability.

The above is a quite broad explanation of the primary factors that define the GHZ – in practice, the GHZ is troublesome to elucidate, as it depends on a host of secondary factors that are difficult to model, and as such there is a great deal of disagreement over the shape and evolution of the Milky Way (MW)'s GHZ.

Gonzalez et al. (2001)'s initial study outlines the key constraints on the MW's GHZ – the metallicity gradient, and the supernova rate distribution. This analysis results in the ‘classic’ picture of the GHZ as an annulus in the thin disc, which moves outwards in time as metallicity is generated in the Galaxy. They assume a minimum metallicity, of approximately half that of Solar, is required for terrestrial planet formation. In their picture, increasing metallicity increases the chance for planet formation (but their calculations of radiogenic heating show that elemental composition is as important as simply finding elements more massive than hydrogen and helium).

Lineweaver et al. (2004) would improve on this analysis with a more sophisticated model of the distribution of terrestrial planets in space and time (Lineweaver 2001), in particular employing the latest statistical data on the probability of planets given stellar metallicity, rather than assuming a simple metallicity threshold for planet formation. This work indicated that the MW's GHZ is an annulus between 7 and 9 kiloparsecs from the Galactic Centre, with the Sun close to the centre of this zone at approximately 8 kiloparsecs.

Prantzos (2007) was highly critical of these approaches, and noted the uncertainties that then surrounded the planet-metallicity relation. In particular, the planet-metallicity relation applied more strictly to giant planets, especially given the observational biases placed upon the data at that time. Indeed, it seemed that high metallicity could damage terrestrial planet formation through the migration of Hot Jupiters from the outer planetary system to very close orbits. Given this, and the uncertainties in the metallicity gradient, Prantzos (2007) proposed that the GHZ could encompass the entire MW disc!

In light of this, later works have been more critical of the GHZ limits they produce. One of the most sophisticated models of the MW's GHZ is that of Gowanlock et al. (2011), which attempts to map the GHZ by modelling individual star systems. As a result, they are better placed to describe habitability above and below the Galactic midplane, as well as the properties of individual planetary systems. They give a much more detailed account of both Type I and Type II SNe, as well as tidal locking of planets to their parent star.

In recent years, other galaxies have had their GHZs mapped, in particular elliptical galaxies. Suthar & McKay (2012), working principally from metallicity constraints, suggested that the elliptical M87, and the dwarf Galaxy M32 could possess relatively broad GHZs. Later works (Carigi et al. 2013) used one-dimensional (1D) Galaxy chemical evolution models to map the GHZ of Andromeda (M31), which incorporated more realistic physics, such as radial mixing (Spitoni et al. 2014) to improve the metallicity gradient.

It is therefore clear that each Galaxy's GHZ depends critically on its evolutionary history, in particular on its accretion history. Understanding this history is a crucial goal of Galaxy formation theory, which aims to explain the observed evolution of Galaxy properties as a function of time.

According to the standard paradigm, Dark Matter (hereafter DM) collapses to form halos with complex baryonic physics governing Galaxy formation within these potential wells (White & Rees 1978). However, the process of Galaxy formation is much more complicated, with DM halos growing both through anisotropic accretion and mergers and Galaxy formation proceeding by the combined action of clumpy and anisotropic gas accretion and mergers of dwarf galaxies. These processes result in clumpy and disky galaxies whose anisotropic shapes are only now beginning to be captured by the latest hydrodynamic cosmological simulations (e.g. Stinson et al. 2006; Libeskind et al. 2010; Kobayashi & Nakasato 2011; Knebe et al. 2011; Pilkington et al. 2012; Scannapieco et al. 2012; Vogelsberger et al. 2014; Graziani et al. 2015).

Most studies of the GHZ to date assume a relatively steady, monotonic growth of the MW from an infall model, or in the case of Gowanlock et al. (2011), the MW's stellar distribution is initialized in its present day configuration, and azimuthal symmetry is still assumed (their study is restricted to a 1° sector of the Galaxy to boost spatial resolution).

In this work, we take a different approach to GHZ modelling. We analyse cosmological simulations of Galaxy formation, which were designed to explore the formation history of the Local Group, whose principal components are the MW, Andromeda (M31) and Triangulum (M33) galaxies. By analysing the habitability of star particles in the simulation, we are able to investigate non-axisymmetric habitable zones, and their non-trivial evolution with redshift ² .

Method

Cosmological simulations

We use the CLUES (http://www.clues-project.org) simulations for the calculations presented in this work. We briefly summarize the simulations here and the interested reader is referred to Libeskind et al. (2010), Knebe et al. (2010), Libeskind et al. (2011) and Knebe et al. (2011) for complete details on the central galaxies, their discs and the properties of their z = 0 substructure population, as well as details regarding the simulation initial conditions, gas-dynamics and star formation.

The simulations used in this study were run with the PMTree-Smoothed Particle Hydrodynamics (SPH) code GADGET2 (Springel 2005) in a cosmological box of size 64 h ⁻¹ comoving Mpc (cMpc) ³ . The runs employed the standard cosmological model of the Universe, referred to as Lambda cold dark matter (ΛCDM) due to the presence of both cold dark matter and dark energy as the principal contributions to the Universe's energy density. This defines standard ΛCDM initial conditions for the particles and parameters for the relative contributions of each component of the Universe's energy density, which are taken from the Wilkinson Microwave Anisotropy Probe Release 3 (WMAP3) data release (Spergel et al. 2007) such that Ω_m = 0.24, Ω_b = 0.042, Ω_Λ = 0.76, h = 0.73, σ₈ = 0.73 and n = 0.95 (the reader is referred to Peacock (1999) for definitions of these terms).

The z = 0 density field is described by a reconstruction that uses observations of objects in the local Universe as constraints. Initial conditions are then selected using the Hoffman-Ribak method (Hoffman & Ribak, 1991). These constrained initial conditions force the z = 0 linear scales to resemble the local Universe; non-linear scales are unconstrained. To obtain a reliable LG comprising the MW) and Andromeda (M31), around 200 low-resolution (256³ particles) constrained runs were performed until a suitable candidate (defined in terms of the mass, relative distance and velocity of the MW and M31 haloes) was found. Once selected, particles in a sphere of radius 2 h ⁻¹ cMpc at the z = 0 location of this region were then traced back to the initial condition (at z = 100) and replaced with higher resolution particles (equivalent to having 4096³ particles spanning the full cosmological box) using the prescriptions given by Klypin et al. (2001). This refining results in producing two objects resembling the MW and M31 in the high-resolution region, each of which contains about 10⁶ particles within their virial radii at z = 0. Outside of this region, the simulation box was populated with lower-resolution (i.e. higher mass) particles to mimic the correct large scale environment of the LG. The DM and baryons (i.e. gas and star particles) in the refined region have masses of $2.1 \times 10^5 h^{ - 1} M_{\odot} $ and $4.4 \times 10^4 h^{ - 1} M_{\odot} $ , respectively.

Around 500 processors and 1.5 terabytes of memory were used to generate the initial conditions for these simulations. The high-resolution N-Body plus gas (SPH) simulations used 500 MPI processors with a typical computing time of the order of 1 million central processing unit (CPU) hours, yielding in excess of 6.1 terabytes of data (Gottloeber et al. 2010). The resulting gas distribution of the simulated LG can be seen in Fig. 1.

Fig. 1. The gas distribution of the Local Group in the CLUES simulations on large scales (left picture, about 2 Mpc h⁻¹ across, viewed from a distance of 3.3 Mpc h⁻¹) and the gas disks of the three main galaxies (right panels, about 50 kpc h⁻¹ across, from a distance of 250 kpc h⁻¹). For the zoomed pictures the colour mapping is shifted to higher densities (factor 10°.5) in order to enhance the spiral arm features of the gas disks. From the CLUES project, image courtesy of K. Riebe.

As for the baryonic modelling, the initial mass function (IMF) used is Salpeter between $0.1$ and $100M_{\odot} $ . This is less than ideal for our purposes, as a Salpeter IMF over-estimates the proportion of low mass stars in a given population, and as a result can distort the calculated supernova rates (cf Gibson 1997). However, our calculations must remain consistent with the simulation's own properties. The simulation employs the feedback rules of Springel & Hernquist (2003): hot ambient gas and cold gas clouds in pressure equilibrium form the two components of the interstellar medium. Gas properties are calculated assuming a uniform but evolving ultraviolet background generated by quasi stellar objects (QSOs) and Active Galactic Nuclei (AGNs) (Haardt & Madau 1996). Metal line (or molecular) dependent cooling below 10⁴ K is ignored. Star formation is treated stochastically, choosing model parameters that reproduce the Kennicutt law for spiral galaxies (Kennicutt 1983; Kennicutt 1998). This empirical relation gives the star formation rate surface density as a powerlaw of the surface density of molecular gas (or in some cases the surface density of cold molecular gas, which produces a slightly different powerlaw exponent).

Each gas particle was allowed to undergo two episodes of star formation, each time spawning a star particle of half the original mass (i.e. $2.2 \times 10^4 h^{ - 1} M_{\odot} $ ). The star particles have an effective radius of 150 pc, defining a resolution limit for our habitability calculations. The instantaneous recycling approximation is assumed: cold gas cloud formation (by thermal instability), star formation, the evaporation of gas clouds and the heating of ambient gas by supernova driven winds all occur at the same instant. We assume kinetic feedback in the form of winds driven by stellar explosions. Assuming instantaneous recycling limits our ability to make habitability arguments based on the nitrogen, carbon, hydrogen, oxygen, phosphorus and sulphur (NCHOPS) elements, as we cannot model the production of various chemical species over stellar lifetimes. Instead, we must make more general arguments as to the availability of elements above hydrogen and helium in the Periodic Table.

Haloes and subhaloes are identified using the publicly available Amiga Halo Finder (AHF; Knollmann & Knebe 2009). AHF uses an adaptive grid to find local over-densities; after calculating the gravitational potential, particles that are bound to the potential are assigned to the halo.

We identify all the particles within the virial radius of MW and M31 at z = 0 and follow them back in time in order to locate the LG progenitors at high-z. If bound to a structure, we identify the corresponding halo: each of these halos was then considered a progenitor of the LG at high-z and used in our calculations, as explained in what follows (see Dayal & Libeskind 2012; Dayal et al. 2013 for complete details of the halo identification procedure).

Calculating habitability

The cosmological simulations deliver a list of ‘star’ particles {i} for each snapshot. Each particle has a total mass $M_i = 3.1612 \times 10^4 M_{\odot} $ , equivalent to a moderately sized star cluster, so we should consider each particle as consisting of an ensemble of stars.

Each particle has an effective radius of around 500 parsecs. SNe at distances of around 8 parsecs or less could produce enough ionizing flux at the surface of a terrestrial planet to affect its ozone layer and disrupt the biosphere (Gehrels et al. 2003). As such, we do not expect neighbouring star particles to produce sufficiently strong SNe or ionizing flux to influence each other's habitability, and as such we can consider the habitability of each particle independently. We note that in reality, star clusters will dissolve and disperse on timescales of a few Myr, and diffuse azimuthally throughout the Galaxy on timescales of an orbital period, i.e. a few 100 Myr (see e.g. Moyano Loyola et al. 2015). The CLUES simulations cannot model this diffusion, so our calculations do not capture the consequences, e.g. the diffusion of supernova progenitors into the vicinity of other star particles.

Each particle has an associated star formation rate S _i , and all share the same IMF, i.e. the same distribution of stellar masses. The cosmological simulation assumes a Salpeter IMF, and therefore for consistency we also assume it here:

(1) $$P(M_{\ast} )dM_{\ast} \propto M_{\ast}^{ - \alpha},$$

where M _* is star mass, and the constant α = 2.35. By specifying a minimum and maximum mass (0.1 and 100 $M_{\odot} $ , respectively), integrating a normalized version of this expression fixes the number of stars in each star particle, as well as the proportions of stars of a given mass, which we use to calculate Type Ia and Type II supernova (SNe) rates. We are required to assume a Salpeter IMF, but this function overestimates the true proportion of low mass stars which can affect SNe rates by a factor of a few for a given stellar population.

Type II SNe progenitors are stars with masses above 8 $M_{\odot} $ , therefore we calculate the fraction of stars that are formed above this critical mass, and multiply this by the star formation rate to obtain a Type II SNe rate. This approach does not take into account the time delay between star formation and supernova. The longest lived stars that produce Type II SNe are the lowest mass, and they remain on the Main Sequence for approximately 40 million years.

Type Ia SNe progenitors are white dwarfs, with masses below the Chandrasekhar limit of 1.4 $M_{\odot} $ . However, not all stars below this mass will become Type Ia SNe – the process is believed to require the presence of extra material for the white dwarf to accrete, such as a companion star. As white dwarf progenitors span the 0.08–8 $M_{\odot} $ mass range, and it is estimated that 1% of the white dwarf population ignites as a Type 1a SNe (Pritchet et al. 2008), we therefore assume that the probability of any individual star below 8 $M_{\odot} $ exploding as a Type 1a SNe is 0.01.

Again, we assume that the Type Ia SNe rate is directly related to the production rate of progenitors, and does not incorporate a time delay due to stellar evolution, which will be many billions of years for the lowest mass objects. We return to this no-delay issue in the section Discussion. The above prescription fixes the ratio of Type II to Type 1a SNe to ~2, which is slightly lower than the measured ratios of 3–5 for disk galaxies (cf Mannucci et al. 2007).

Star particles that experience combined SNe rates that are too high will experience a reduced probability that the habitable planets they host will possess surviving biospheres (see later).

The simulation also delivers metallicities for each star particle. These metallicities allow us to calculate two probabilities related to planets:

1. 1. The probability that terrestrial planets may form in the SHZ, P _{plannet, t}

2. 2. The probability that close-in giant planets may form, P _{plannet, g}

We wish to compare with azimuthally symmetric studies, so we therefore follow Gowanlock et al. (2011) by suggesting that for a system to produce habitable worlds, it should produce terrestrial planets in the SHZ without producing close-in giant planets. Giant planets cannot be formed in situ close to the star; instead, they are formed at larger distances and migrate inwards. We assume that such events are destructive to habitable terrestrial planets ‘in the way’ of giant planet migration, but we note that simulations indicate that giant planet migration may improve the odds of terrestrial planet formation in the SHZ in some cases (Raymond et al. 2005; Fogg & Nelson 2007).

We assume terrestrial planet formation occurs when the metallicity Z rises above a threshold, following Prantzos (2007):

(2) $$P_{{\rm planet},\,t} = \left\{ {\matrix{ {0.4} & {\displaystyle{Z \over {Z_{\rm \odot}}} \gt 0.1} \cr {0.0} & {\displaystyle{Z \over {Z_{\odot}}} \le 0.1} \cr}} \right.$$

This is in line with observations that show the correlation between host star metallicity and the presence of a planet is much weaker below Neptune masses (Fischer & Valenti 2005; Buchhave et al. 2012; Wang & Fischer 2014; Schlaufman 2015). The exact value of the boundary we use here is not tailored to the latest Kepler results, but it does produce terrestrial planets for a large fraction of star systems, which generally speaking is thought to be the case given the Kepler dataset (Petigura et al. 2013). Above Neptune masses, the metallicity correlation is stronger, so we assume that the probability of giant planet formation is governed by:

(3) $$P_{{\rm planet},\,g} = \left\{ {\matrix{ {0.03 \times 10^{log\displaystyle{Z \over {Z_{\rm \odot}}}}} & {\displaystyle{Z \over {Z_{\odot}}} \gt 1} \cr {0.03} & {\displaystyle{Z \over {Z_{\odot}}} \le 1} \cr}} \right.$$

Therefore, the probability of a terrestrial planet being habitable is (1 − P _{plannet, t} ), and the number of habitable planets hosted by a star particle is

(4) $$N_{{\rm hab}} = N_{{\rm planet}} (1 - P_{{\rm planet},g} ) = N_{\ast} P_{{\rm planet},t} (1 - P_{{\rm planet},g} )$$

Habitable planets that exist in hazardous environments will be sterilized. We can therefore propose a probability of survival, related to the local SNe rate SNR (supernova rate):

(5) $$P_{{\rm survive}} = 1 - \displaystyle{{{\rm SNR}} \over {2{\rm SNR}_{\odot}}} $$

We propose this ansatz to allow the probability of survival to be at a maximum when the local SNe rate is zero, and to be zero if the SNR is twice the Solar rate. This is similar to many previous studies (Lineweaver et al. 2004; Carigi et al. 2013), although in some cases a step function is used, as opposed to a linear approximation. As we cannot track individual stellar positions, we cannot use Gowanlock et al. (2011)'s more detailed approach of tracking individual stellar distances to nearby SNe, and hence we are forced to use this ansatz, which is clearly oversimplified. This is a weak point in our analysis. While our results are not sensitive qualitatively on the above criterion, we should expect some quantitative effects on our calculations, but even a more appropriate, detailed calculation is still subject to general uncertainty in how SNe can sterilize habitable planets (see the section Discussion).

Once calculated, this survival probability can be multiplied by the number of habitable planets to give the number of planets that remain in clement conditions, N _survive:

(6) $$N_{{\rm survive}} = N_{{\rm hab}} P_{{\rm survive}} $$

We do not place age constraints on habitability, unlike other studies, which demand a minimum of 4 Gyr of clement conditions to produce complex life forms. The pros and cons of such constraints are addressed in the section Discussion.

Results

For brevity, we analyse only two of the LG objects – the MW, and M33. The object representing Andromeda (M31) was omitted as it did not form a well-defined disc – this can be seen in Fig. 3 of Libeskind et al. (2013) ⁴

M33

Figure 2 shows the evolution of metallicity, star formation rate and number density of stars in the M33 simulation. The main stellar disc truncates at around 10 kpc, with a broadly exponential surface density profile. The disc is sufficiently metal rich at early times to surpass the threshold for terrestrial planet formation in the inner regions. However, significant star formation at early times suppresses habitability, with the disc more or less fully formed by t = 9.5 Gyr. The age-metallicity relation has a large amount of scatter, which is to be expected as the metals are fully entrained in the gas, and not allowed to diffuse (see e.g. Pilkington et al. 2012).

Fig. 2. The physical properties of the M33 galaxy. We plot the local average stellar metallicity (top left), the local average star formation rate (top right) and the number density of stars (bottom left) at three instances in the simulation, and the age-metallicity relation for all stars in the simulation (bottom right) at z = 0.

Figure 3 shows how N _survive evolves over the course of the simulation. Each panel represents a 2D binning to give the average value of N _survive (bins with N _survive = 0 are not plotted). Note that the simulations have been rotated so that the z axis is aligned with the system's global angular momentum vector.

Fig. 3. The evolution of N _survive in the M33 galaxy. The 2D binned value of N _survive is shown for the x–y plane (left) and x–z plane (right).

The early snapshots (top panels of Fig. 3) show indications of two-armed spiral structure, but remain poorly resolved. A well-defined disc is formed around 4 Gyr from the present day (middle panels of Fig. 3), but it is clear that the GHZ is very non-axisymmetric before disc formation. Initially, the most habitable regions are to be found in clumps orbiting greater than 10 kpc from the centre of the Galaxy, and this persists to the present day (bottom panels). In particular, satellite galaxies in particular appear to have broad GHZs in this analysis.

As the disc forms, the habitability of the inner regions increases, but a halo of good habitability remains at the disc's edges. The radial profiles of M33 (right panels of Fig. 4) show that large numbers of planets remain unsterilized down to distances of 1 kpc from the centre. Presumably at distances within 1 kpc, the ionizing radiation produced by the Galaxy's supermassive black hole (not modelled in this analysis) is likely to render this region inhospitable.

Fig. 4. Left: The axial distributions of surviving planets in the disc of the M33 galaxy. In each snapshot, the number of surviving planets is binned along the x, y and z axes. The strong peak in the z-curve between –10 and –5 kpc is due to the presence of a satellite orbiting above the galactic plane. Right: the radial distribution of surviving planets.

However, we do see that at regions between around 4 and 10 kpc, the number of surviving planets increases slightly over the linear trend (in log-log space). Beyond the disc's edge, spikes in survivability belie the presence of satellites at these galactocentric radii.

As we are analysing 3D data, we have the advantage of comparing the radial habitability of M33 with its habitability along the x, y and z axes (left panel of Fig. 4). Initially, the simulation shows habitable planets forming in an azimuthally symmetric configuration. Over time, the formation of the disc narrows the range of habitability in z, and the x and y profiles are quite similar (but not identical). As satellites and streams merge into the disc, we can see asymmetries grow and disappear in the profiles. As we reach the present day, the GHZ becomes more uniform, but the simulations indicate that such uniformity is a recent phenomenon. Note that our vertical density distributions are quite similar to those shown by Jurić et al. (2008), suggesting that these galaxies follow exponential profiles both radially and vertically.

MW

We now consider the other Galaxy in our simulation data, which is selected as it partially resembles the MW. Figure 5 again shows the Galaxy's physical properties, which suggests a more punctuated, chaotic formation history compared with M33. As can be seen also in Fig. 6, this Galaxy does not form a well-defined disc, but exhibits a more triaxial structure. It has accreted more material ‘off-axis’ from its angular momentum vector, as the flow of material from streams and satellites has been more isotropic. This accretion history is hinted at by the Galaxy's large number of satellites, even at the present day, and the satellites in general contain more star particles. The total number of star particles at the present day is 545 349 which is nearly double that of M33. Star formation rates remain relatively high even at late times, reducing the number of habitable planets in the outer disc, and as a result favouring locations beyond the disc, especially its larger satellites, for surviving biospheres until much later times.

Fig. 5. As Fig. 2, but for the Milky Way.

Fig. 6. As Fig. 3, but for the Milky Way.

Figure 7 shows the radial and axial profiles of this Galaxy. The larger number of satellites for this galaxy produce larger asymmetries, particularly in the axial profiles. There is much less flattening of habitability in the z-direction than in the case of M33, and the radial profiles follow a more uniform trend – there is no ‘bump’ in habitability in the 4–10 kpc range as seen in M33. Once again, we assume that the inner 1 kpc of this system, while recording a large number of surviving planets by our analysis, should be discounted due to the effect of the Galaxy's central supermassive black hole.

Fig. 7. As Fig. 4, but for the Milky Way.

Discussion

General trends

Several features are present in both galaxies:

1. 1. The main Galaxy produces relatively broad GHZs,

2. 2. The satellites also produce broad GHZs,

3. 3. Streams of stars between satellites and galaxies produce relatively large numbers of habitable planets. The satellites benefit from low SNe rates, old stellar populations, less major mergers and the weak metallicity dependence of terrestrial planet formation.

4. 4. The inner regions of the galactic disc contain a large number of surviving worlds by virtue of high stellar density, but the outer regions produce higher survival probabilities for individual star particles.

These trends are in line with recent GHZ calculations, such as Gowanlock et al. (2011), who also show improved habitability above the midplane and at galactocentric radii of order a few kpc, much lower than the canonical limits given by Lineweaver et al. (2004). Given that we use some of Gowanlock et al.'s prescriptions for planet formation, etc., this is perhaps not surprising. However, their models do not account for tidal streams and satellites. Given our fairly loose constraints on terrestrial planet formation, and reduced metallicity suppressing giant planet formation, these structures on the outer limits of galaxies appear to be good sites for habitable planets. However, we caution that our model assumptions remain simplistic, and there may be other hazards in these regions that we do not simulate.

Table 1 gives the first and third quartiles of N _survive as a function of radius (equivalently, the 25th and 75th percentiles), and this shows some important differences between the M33 and MW simulations, where we have again ignored the inner 1 kpc. The M33 quartiles are initially well separated, and draw closer together with time to give a range between 3 and 10 kpc. The MW shows a less uniform evolution, where the quartiles are initially very close together, on the disc's outskirts (due to vigorous star formation in the interior), and then draw very far apart at intermediate times before settling to values of approximately 2 and 13 kpc.

Table 1. The first and third quartiles (Q1 and Q3) of N_survive as a function of radius for both galaxies

As with all GHZ studies, we must consider these results in the light of the assumptions made, and the physics still missing from the analysis, as we do in the following section.

Limitations of our results

This work is proof in principle that habitability can be determined from cosmological simulations – however, the simplicity of the analysis (like many GHZ analyses) leaves the conclusions we draw open to uncertainty.

Much of our issues stem from the fact that these simulations were not designed with GHZ calculation in mind. Perhaps most pressingly, the resolution of the simulations limit us from exploring habitability at scales below a few hundred parsecs. To date, no cosmological simulations exist that can resolve both the formation of individual stars and the global evolution of Galaxy groups. Our work pushes at the limits of what is currently possible with state-of-the art algorithms and computational power, but we are optimistic that future simulations will be able to address this.

Given that most sterilizing events take place within a few tens of parsecs from a planet, we cannot investigate how neighbouring star particles influence each other. While this has allowed a certain expediency in our calculations, it prevents us from calculating sterilization rates convincingly in dense stellar environments, especially once natal star clusters have dissolved (Moyano Loyola et al. 2015). It also prevents us from measuring the effect of stellar encounters on local habitability (Jiménez-Torres et al. 2013), which is impossible without better estimates of the kinematics of individual stars. Semi-analytic prescriptions for such properties could be employed, but we leave this for future work.

Our resolution limits also limit our ability to model individual sterilization events. Each star particle's star formation rate is calculated from semi-analytic approximations to reproduce empirical relations for star formation, in particular the Kennicutt-Schmidt relation (Kennicutt 1983, 1998). We estimate supernova rates based on these star formation rates, assuming no time lag between star birth and star death. However, Type I SNe progenitors are low mass stars with a relatively long lifetime, and as such we should expect supernova rates that are less intense at early times than we calculate, with more supernova activity at later times as lower mass stars begin to ignite. We should therefore note that our simulations underestimate habitability for young star particle ages, and overestimate it at older ages. On the other hand, our prescriptions for Type I and Type II SNe tend to underproduce Type I relative to Type II, which goes some way to correcting this issue. That being the case, this time lag issue is not the principal uncertainty in our habitability calculation, it being superseded by our uncertainty about the required radiation dose for sterilization, as we discuss later.

The cosmological simulations do track the metallicity of the gas assuming this supernova activity is occurring, but the individual abundances of chemical elements are not traced. This simplification is advantageous when attempting to determine broad trends in galaxy formation, but it leaves us unable to make more nuanced predictions about regions of favourable chemical composition. For example, the abundance of iron and oxygen (and to some extent silicon) may give strong indications where there are large reservoirs of material to build planets with similar bulk compositions to the Earth, but a biota (at least based on what we know of terrestrial life) requires the presence of NCHOPS as a mimimal set of elements colocalized at the scale of the organism and co-located with liquid water and an energy supply. These requirements cannot be easily predicted for any given planet.

In a similar vein, we do not consider the ‘metals’ as a separate entity to the gas in the simulations. While there are many density and pressure regimes where we expect volatile elements and compounds (‘gas’) to entrain the refractory materials (‘dust’) on similar dynamical trajectories, there are equally regimes where we expect these two phases to separate. This separation is clearly visible in optical observations of the ‘dust lanes’ of galaxies at various wavelengths, where dust obscures our view of the interior gas. This evolution will affect the local metallicity, but it also affects the optical depth of the medium to biohazardous radiation. One could imagine circumstances where planetary systems in regions of enhanced dust density are afforded a shield against ionizing radiation and cosmic rays. Equally, this biohazardous radiation may become biogenic radiation if it encounters dust/gas mixtures amenable to grain surface chemistry, enhancing the production of organic molecules and potentially even prebiotics (Hill & Nuth 2003; Stark et al. 2014). This will have significant effects on the calculated habitability of low density regions such as tidal streams and the outer edges of satellite galaxies, and shows that a fully matured GHZ calculation must account for the anisotropy of radiative feedback as it encounters asymmetric distributions of intervening material.

On the subject of strong radiation sources, while we have considered the effect of AGN implicitly by excluding the inner kiloparsec from our considerations of habitability, we do not explicitly calculate the interaction between the AGN and its environs. AGN go through quiet and active phases depending on the black hole's accretion rate and the galaxy's own accretion history, and we expect their effects to extend well beyond the 500 pc radius of our star particles. Indeed, it is thought that ‘AGN feedback’ may be an important mechanism for quenching star formation (Springel et al. 2005), although it might also have the opposite effect (Ishibashi & Fabian 2013). If the accretion rate of the central supermassive black hole – and consequently the fuelling of the central engine of the AGN – is included in the simulation, then the models could link mergers and accretion to AGN activity and subsequent star formation, which has consequences both for Galaxy evolution and galactic habitability.

Another powerful source of sterilizing radiation comes from gamma ray bursts (GRBs), powerful explosions generated by the mergers of compact objects, or the collapse of extremely massive stars (see Fishman & Meegan 1995; Berger 2014 for reviews). As GRBs tend to be more frequent in the past than in the present, and their sterilization radius is large compared with the Galactic disc, the GHZ at early times could be even more asymmetric and patchy, if it exists at all. Recent studies have characterized the GRB threat to Earth by studying the GRB luminosity function (Piran & Jimenez 2014), which confirms the trend that GRBs are more frequent in low metallicity environments. As a result, the GRB rate is higher at earlier times, which suggests that the MW and most other galaxies might have very little GHZ to speak of before z = 0.5. This sort of ‘global regulation mechanism’ has been proposed to answer Fermi's Paradox, or the observed ‘Great Silence’ of transmissions from, or other evidence of, intelligent technological civilisations (Annis 1999; Vukotic & Cirkovic 2007). More modelling in this area is essential.

Perhaps most importantly, this work has demonstrated the current limitations of the GHZ concept. Modelling biosphere formation using a single parameter, metallicity, is unlikely to yield the ‘ground truth’ of how habitable worlds are formed. Our understanding of the timescales required for metazoan lifeforms to evolve remains informed by Earth's geological history, but essentially incomplete due to gaps in our knowledge of comparative exoplanetology, for example, our understanding of plate tectonics on super Earths, which remains a source of debate (Stamenkovic et al. 2010; Stamenković & Breuer 2014). We also remain quite ignorant of the levels of hazardous radiation required to definitively exterminate a biosphere, evidenced by the fact that the only known biosphere – our own – continues to persist. A significant fraction of damaging cosmic rays are absorbed or scattered away from the Earth by the Sun's heliosheath, and the Earth's own protective magnetic field. Simply stating that a SN is ‘deadly within 8pc’ is an extremely glib statement, and belies the subtleties involved in determining the effect of hazardous interstellar radiation on planets within stellar magnetospheres (Gehrels et al. 2003; Martin et al. 2009; Thomas et al. 2015). Future GHZ modelling must address these concerns, first put forward by Prantzos (2007), and which remain unresolved.

Our principal contribution to this debate has been to demonstrate that the GHZ will be morphologically complex. We have used state-of-the-art simulations to reproduce galaxy assembly histories, but have been forced to use relatively simple chemical evolution modelling and planet formation criteria to leverage the simulation data in ways that it was not designed. Future work in this area must use better chemodynamical models to better trace the distribution of NCHOPS elements, as well as the stellar population (see e.g. Gibson et al. submitted), while still resolving the assembly history of the LG accurately in cosmological context.

Conclusions

For the first time, we have estimated galactic habitability in three dimensions by analysing high-resolution numerical simulations of the formation of the LG. We apply habitability metrics commonly found in other studies of 1D azimuthally symmetric GHZs.

Thanks to mergers and the accretion of satellite galaxies, we find that if the GHZ exists, it must be fundamentally asymmetric. These dynamical events produce tidal streams and spiral arms which produce elliptical and even triaxial regions of improved habitability. While the simulations indicate that large numbers of habitable, non-irradiated planets can be found in the inner regions of galaxies and their satellites, the probability that an individual planetary system contains habitable worlds increases with distance from the galactic centre, with its peak near the edges of galactic discs. Regions above the midplane are also favoured locations, but the density of stars in these areas reduces the total number of habitable planets.

However, the idea that the outer edges of galactic discs are ‘more habitable’ is not wholly confirmed by our analysis (see e.g. Table 1). For the MW, this appears to be true at early times, where the first and third quartiles of all habitable planets are initially at large distances, but over time this interquartile range increases greatly before settling to values of approximately 2 and 13 kpc. For M33, we see quite opposite evolution, but a similar final interquartile range between 3 and 10 kpc. These large changes suggests that Prantzos (2007)'s criticisms are vindicated – we do not see a GHZ that is well constrained in space/time, and habitable planet formation appears to be possible at many locations in a Galaxy. Furthermore, the divergence in evolution between both galaxies suggests that assembly history plays a crucial role in subsequent habitability.

However, we acknowledge that galactic habitability is still maturing as a conceptual tool. This paper has demonstrated that cosmological simulations can be extremely informative regarding both galactic habitability and the weaknesses of our assumptions about planet formation and sterilization. As these simulations continue to improve in spatial resolution and input physics, we believe that the GHZ will become as useful as the circumstellar habitable zone in assessing the astrobiological potential of our local cosmic volume.