Steve: A hierarchical Bayesian model for Supernova Cosmology

S. R. Hinton11affiliationmark: , T. M. Davis11affiliationmark: , A. G. Kim22affiliationmark: , D. Brout33affiliationmark: , C. B. D’Andrea33affiliationmark: , R. Kessler44affiliationmark: 55affiliationmark: , J. Lasker44affiliationmark: 55affiliationmark: , C. Lidman66affiliationmark: , E. Macaulay77affiliationmark: , A. Möller88affiliationmark: 66affiliationmark: , M. Sako33affiliationmark: , D. Scolnic55affiliationmark: , M. Smith99affiliationmark: , R. C. Wolf1010affiliationmark: , S. Allam1111affiliationmark: , J. Annis1111affiliationmark: , S. Avila77affiliationmark: , E. Bertin1212affiliationmark: 1313affiliationmark: , D. Brooks1414affiliationmark: , D. L. Burke1515affiliationmark: 1616affiliationmark: , A. Carnero Rosell1717affiliationmark: 1818affiliationmark: , M. Carrasco Kind1919affiliationmark: 2020affiliationmark: , J. Carretero2121affiliationmark: , M. Childress99affiliationmark: , C. E. Cunha1515affiliationmark: , L. N. da Costa1818affiliationmark: 2222affiliationmark: , C. Davis1515affiliationmark: , J. De Vicente1717affiliationmark: , D. L. DePoy2323affiliationmark: , P. Doel1414affiliationmark: , T. F. Eifler2424affiliationmark: 2525affiliationmark: , B. Flaugher1111affiliationmark: , P. Fosalba2626affiliationmark: 2727affiliationmark: , J. Frieman1111affiliationmark: 55affiliationmark: , J. García-Bellido2828affiliationmark: , E. Gaztanaga2626affiliationmark: 2727affiliationmark: , D. W. Gerdes2929affiliationmark: 3030affiliationmark: , R. A. Gruendl1919affiliationmark: 2020affiliationmark: , J. Gschwend1818affiliationmark: 2222affiliationmark: , G. Gutierrez1111affiliationmark: , W. G. Hartley1414affiliationmark: 3131affiliationmark: , D. L. Hollowood3232affiliationmark: , K. Honscheid3333affiliationmark: 3434affiliationmark: , E. Krause2424affiliationmark: , K. Kuehn3535affiliationmark: , N. Kuropatkin1111affiliationmark: , O. Lahav1414affiliationmark: , M. Lima3636affiliationmark: 1818affiliationmark: , M. A. G. Maia1818affiliationmark: 2222affiliationmark: , M. March33affiliationmark: , J. L. Marshall2323affiliationmark: , F. Menanteau1919affiliationmark: 2020affiliationmark: , R. Miquel3737affiliationmark: 2121affiliationmark: , E. Morganson2020affiliationmark: , R. L. C. Ogando1818affiliationmark: 2222affiliationmark: , A. A. Plazas2525affiliationmark: , E. Sanchez1717affiliationmark: , V. Scarpine1111affiliationmark: , R. Schindler1616affiliationmark: , M. Schubnell3030affiliationmark: , S. Serrano2626affiliationmark: 2727affiliationmark: , I. Sevilla-Noarbe1717affiliationmark: , M. Soares-Santos3838affiliationmark: , F. Sobreira3939affiliationmark: 1818affiliationmark: , E. Suchyta4040affiliationmark: , G. Tarle3030affiliationmark: , D. Thomas77affiliationmark: , V. Vikram4141affiliationmark: , and Y. Zhang1111affiliationmark:
(DES Collaboration)
School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA The Research School of Astronomy and Astrophysics, Australian National University, ACT 2601, Australia Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK ARC Centre of Excellence for All-sky Astrophysics (CAASTRO) School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK Graduate School of Education, Stanford University, 160, 450 Serra Mall, Stanford, CA 94305, USA Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA CNRS, UMR 7095, Institut d’Astrophysique de Paris, F-75014, Paris, France Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d’Astrophysique de Paris, F-75014, Paris, France Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450, Stanford University, Stanford, CA 94305, USA SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA Department of Astronomy/Steward Observatory, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain Instituto de Fisica Teorica UAM/CSIC, Universidad Autonoma de Madrid, 28049 Madrid, Spain Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA Department of Physics, The Ohio State University, Columbus, OH 43210, USA Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain Brandeis University, Physics Department, 415 South Street, Waltham MA 02453 Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 13083-859, Campinas, SP, Brazil Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831 Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA

We present a new Bayesian hierarchical model (BHM) named Steve for performing type Ia supernova (SN Ia) cosmology fits. This advances previous works by including an improved treatment of Malmquist bias, accounting for additional sources of systematic uncertainty, and increasing numerical efficiency. Given light curve fit parameters, redshifts, and host-galaxy masses, we fit Steve simultaneously for parameters describing cosmology, SN Ia populations, and systematic uncertainties. Selection effects are characterised using Monte-Carlo simulations. We demonstrate its implementation by fitting realisations of SN Ia datasets where the SN Ia model closely follows that used in Steve. Next, we validate on more realistic SNANA simulations of SN Ia samples from the Dark Energy Survey and low-redshift surveys (DES Collaboration et al., 2018). These simulated datasets contain more than SNe Ia, which we use to evaluate biases in the recovery of cosmological parameters, specifically the equation-of-state of dark energy, . This is the most rigorous test of a BHM method applied to SN Ia cosmology fitting, and reveals small -biases that depend on the simulated SN Ia properties, in particular the intrinsic SN Ia scatter model. This -bias is less than on average, less than half the statistical uncertainty on . These simulation test results are a concern for BHM cosmology fitting applications on large upcoming surveys, and therefore future development will focus on minimising the sensitivity of Steve to the SN Ia intrinsic scatter model.

Subject headings:
cosmology: supernovae

1. Introduction

Two decades have passed since the discovery of the accelerating universe (Riess et al., 1998; Perlmutter et al., 1999). Since that time, the number of observed type Ia supernovae (SN Ia) has increased by more than an order of magnitude, with contributions from modern surveys at both low redshift (Bailey et al., 2008; Freedman et al., 2009; Hicken et al., 2009a; Contreras et al., 2010; Conley et al., 2011), and higher redshift (Astier et al., 2006; Wood-Vasey et al., 2007; Frieman et al., 2008; Balland et al., 2009; Amanullah et al., 2010; Chambers et al., 2016; Sako et al., 2018). Cosmological analyses of these supernova samples (Kowalski et al., 2008; Kessler et al., 2009b; Conley et al., 2011; Suzuki et al., 2012; Betoule et al., 2014; Rest et al., 2014; Scolnic et al., 2017) have been combined with complementary probes of large scale structure and the CMB. For a recent review, see Huterer & Shafer (2018). While these efforts have reduced the uncertainty on the equation-of-state of dark energy () by more than a factor of two, it is still consistent with a cosmological constant and the nature of dark energy remains an unsolved mystery.

In attempts to tease out the nature of dark energy, active and planned surveys are continually growing in size and scale. The Dark Energy Survey (DES, Bernstein et al., 2012; Abbott et al., 2016) has discovered thousands of type Ia supernovae, attaining both spectroscopically and photometrically identified samples. The Large Synoptic Survey Telescope (LSST, Ivezic et al., 2008; LSST Science Collaboration et al., 2009) will discover tens of thousands of photometrically classified supernovae. Such increased statistical power demands greater fidelity and flexibility in modelling supernovae for cosmological purposes, as we will require reduced systematic uncertainties to fully utilise these increased statistics (Betoule et al., 2014; Scolnic et al., 2017).

As such, considerable resources are aimed at developing more sophisticated supernova cosmology analyses. The role of simulations mimicking survey observations has become increasingly important in determining biases in cosmological constraints and validating specific supernova models. First used in SNLS (Astier et al., 2006) and ESSENCE analyses (Wood-Vasey et al., 2007), and then refined and improved for SDSS (Kessler et al., 2009b), simulations are a fundamental component of modern supernova cosmology. Betoule et al. (2014) quantise and correct observational bias using simulations, and more recently Scolnic & Kessler (2016) and Kessler & Scolnic (2017) explore simulations to quantify observational bias in SN Ia distances as a function of multiple factors to improve bias correction. Approximate Bayesian computation methods also make use of simulations, trading traditional likelihoods and analytic approximations for more robust models with the cost of increased computational time (Weyant et al., 2013; Jennings et al., 2016). Bayesian Hierarchical models abound (Mandel et al., 2009; March et al., 2011; March et al., 2014; Rubin et al., 2015; Shariff et al., 2016; Roberts et al., 2017), and either use simulation-determined distance-corrections to correct for biases, or attempt to find analytic approximations for effects such as Malmquist bias to model the biases inside the BHM itself.

In this paper, we lay out a new hierarchical model that builds off the past work of Rubin et al. (2015). We include additional sources of systematic uncertainty, including an analytic formulation of selection efficiency which incorporates parametric uncertainty. We also implement a different model of intrinsic dispersion to both incorporate redshift-dependent scatter and to increase numerical efficiency, allowing our model to perform rapid fits to supernovae datasets.

Section 2 is dedicated to a quick review of supernova cosmology analysis methods, and Section 3 outlines some of the common challenges faced by analysis methods. In Section 4 we outline our methodology. Model verification on simulated datasets is given in Section 5, along with details on potential areas of improvement. We summarise our methodology in Section 6.

2. Review

Whilst supernova observations take the form of photometric time-series brightness measurements in many bands and a redshift measurement of the supernova (or its assumed host), most analyses do not work from these measurements directly. Instead, most techniques fit an observed redshift and these photometric observations to a supernova model, with the most widely used being that of the empirical SALT2 model (Guy et al., 2007, 2010). This model is trained separately before fitting the supernova light curves for cosmology (Guy et al., 2010; Betoule et al., 2014). The resulting output from the model is, for each supernova, an amplitude (which can be converted into apparent magnitude, ), a stretch term , and color term , along with a covariance matrix describing the uncertainty on these summary statistics (). As all supernovae are not identical, an ensemble of supernovae form a redshift-dependent, observed population of , and , where the hat denotes an observed variable.

This ensemble of , and represents an observed population, which – due to the presence of various selection effects – may not represent the true, underlying supernova population. Accurately characterising this underlying population, its evolution over redshift, and effects from host-galaxy environment, is one of the challenges of supernova cosmology. Given some modelled underlying supernova population that lives in the redshift-dependent space (absolute magnitude of the supernova, traditionally in the Bessell band), , and , the introduction of cosmology into the model is simple – it translates the underlying population from absolute magnitude space into the observed population in apparent magnitude space. Specifically, for any given supernova our map between absolute magnitude and apparent magnitude is given by the distance modulus:


where is the mean absolute magnitude for all SN Ia given , is the stretch correction (Phillips, 1993; Phillips et al., 1999), and is the color correction (Tripp, 1998) that respectively encapsulate the empirical relation that broader (longer-lasting) and bluer supernovae are brighter. refers to the host-galaxy mass correlation discussed in Section 4.4.3. The distance modulus is a product of our observations, however a distance modulus can be precisely calculated given cosmological parameters and a redshift. The ‘other corrections’ term often includes bias corrections for traditional analyses. Bias corrections can take multiple forms, such as a redshift-dependent function (Betoule et al., 2014) or a 5D function of , , , and (Kessler & Scolnic, 2017; Scolnic et al., 2017).

2.1. Traditional Cosmology Analyses

Traditional analyses such as that found in Riess et al. (1998); Perlmutter et al. (1999); Wood-Vasey et al. (2007); Kowalski et al. (2008); Kessler et al. (2009b); Conley et al. (2011); Betoule et al. (2014), minimise the difference in distance modulus between the observed distance modulus attained from equation 1, and the cosmologically predicted distance modulus. The function minimised is


where is an uncertainty matrix that combines statistical and systematic uncertainties (see Brout et al. (2018) for a review of these uncertainties for the DES supernova analysis). The predicted is defined as


where is the luminosity distance for redshift given a specific cosmology, is the current value of Hubble’s constant in and , , and represent the energy density terms for matter, curvature and dark energy respectively.

The benefit this analysis methodology provides is speed – for samples of hundreds of supernovae or less, efficient matrix inversion algorithms allow the likelihood to be evaluated quickly. The speed comes with several limitations. Firstly, formulating a likelihood results in a loss of model flexibility by assuming Gaussian uncertainty. Secondly, the method of creating a covariance matrix relies on computing partial derivatives and thus any uncertainty estimated from this method loses information about correlation between sources of uncertainty. For example, the underlying supernova color population’s mean and skewness are highly correlated, however this correlation is ignored when determining population uncertainty using numerical derivatives of population permutations. Whilst correlations can be incorporated into a covariance matrix, it requires human awareness of the correlations and thus methods that can automatically capture correlated uncertainties are preferable. Thirdly, the computational efficiency is dependent on both creating the global covariance matrix, and then inverting a covariance matrix with dimensionality linearly proportional to the number of supernovae. As this number increases, the cost of inversion rises quickly, and is not viable for samples with thousands of supernovae. A recent solution to this computational cost problem is to bin the supernovae in redshift bins, which produces a matrix of manageable size and can allow fitting without matrix inversion at every step. Whilst binning data results in some loss of information, recent works tested against simulations show that this loss does not result in significant cosmological biases (Scolnic & Kessler, 2016; Kessler & Scolnic, 2017).

Selection efficiency, such as the well known Malmquist bias (Malmquist K. G., 1922) is accounted for by correcting the determined from the data, or equivalently, adding a distance bias to the prediction. Specifically, Malmquist bias is the result of losing the fainter tail of the supernova population at high redshift. An example of Malmquist bias is illustrated in Figure 1, which simulates supernovae according to equation (1). Simulations following survey observational strategies and geometry are used to calculate the expected bias in distance modulus, which is then subtracted from the observational data. When using traditional fitting methods such as that found in Betoule et al. (2014), these effects are not built into the likelihood and instead are formed by correcting data. This means that the bias uncertainty is not captured fully in the distribution, and subtle correlations between cosmological or population parameters and the bias correction is lost. Recent developments such as the BBC method (Kessler & Scolnic, 2017) incorporate corrections dependent on and , improving their capture of uncertainty on bias corrections in the likelihood.

Figure 1.— An example of the effects of Malmquist bias. Here are shown 1000 simulated supernovae redshifts and distance modulus given fiducial cosmology. The simulated survey is magnitude limited, and all supernovae brighter than magnitude 24 are successfully observed (shown as blue dots), and all dimmer than 24th magnitude are not successfully observed (shown as red dots). By binning the supernovae along redshift, and taking the mean distance modulus of the supernovae in each bin, we can see that at higher redshift where Malmquist bias kicks in, the population mean drops and becomes biased. This source of bias must either be corrected by adjusting the data (such as subtracting the found bias) for by incorporating Malmquist bias explicitly in the cosmological model.

2.2. Approximate Bayesian Computation

To avoid the limitations of the traditional approaches, several recent methods have adopted Approximate Bayesian Computation, where supernova samples are forward modelled in parameter space and compared to observed distributions. Weyant et al. (2013) provides an introduction into ABC methods for supernova cosmology in the context of the SDSS-II results (Sako et al., 2014) and flat CDM cosmology, whilst Jennings et al. (2016) demonstrates their superABC method on simulated first season Dark Energy Survey samples. In both examples, the supernova simulation package SNANA (Kessler et al., 2009a) is used to forward model the data at each point in parameter space.

Simulations provide great flexibility and freedom in how to treat the systematic uncertainties and selection effects associated with supernova surveys. By using forward modelling directly from these simulations, data does not need to be corrected, analytic approximations do not need to be applied, and we are free to incorporate algorithms that simply cannot be expressed in a tractable likelihood such as those found in traditional analyses from Section 2.1. This freedom comes with a cost – computation. The classical method’s most computationally expensive step in a fit is matrix inversion. For ABC methods, we must instead simulate an entire supernova population – drawing from underlying supernova populations, modelling light curves, applying selection effects, fitting light curves and applying data cuts. This is an intensive process.

One final benefit of ABC methods is that they can move past the traditional treatment of supernovae with summary statistics (, , and ). Jennings et al. (2016) presents two metrics, which are used to measure the distance between the forward modelled population and observed population, and are minimised in fitting. The first metric compares forward modelled summary statistic populations (denoted the ‘Tripp’ metric) and the second utilises the observed supernova light curves themselves, moving past summary statistics. However, we note that evaluation of systematic uncertainty was only performed using the Tripp metric.

2.3. Bayesian Hierarchical Models

Sitting between the traditional models simplicity and the complexity of forward modelling lies Bayesian hierarchical models (BHM). Hierarchical models utilise multiple layers of connected parameters, with the layers linked via well defined and physically motivated conditional probabilities. For example, an observation of a parameter from a population will be conditioned on the true value of the parameter, which itself will be conditioned on the population distribution of that parameter. We can thus incorporate different population distributions, and parameter inter-dependence which cannot be found in traditional analyses where uncertainty must be encapsulated in a covariance matrix. Unlike ABC methods, which can model arbitrary probability distributions, BHM methods are generally constrained to representing probabilities using analytic forms.

With the introduction of multiple layers in our model, we can add more flexibility than a traditional analysis whilst still maintaining most of the computational benefits that come from having a tractable likelihood. Mandel et al. (2009); Mandel et al. (2011); Mandel et al. (2017) construct a hierarchical model that they apply to supernova light-curve fitting. March et al. (2011) derive a hierarchical model and simplify it by analytically marginalising over nuisance parameters to provide increased flexibility with reduced uncertainty over the traditional method, but do not incorporate bias corrections. March et al. (2014) and Karpenka (2015) improve upon this by incorporating redshift-dependent magnitude corrections from Perrett et al. (2010) to remove bias, and validate on 100 realisations of SNLS-like simulations. The recent BAHAMAS model (Shariff et al., 2016) builds on this and reanalyses the JLA dataset (using redshift dependent bias corrections from Betoule et al., 2014), whilst including extra freedom in the correction factors and , finding evidence for redshift dependence on . Ma et al. (2016) performed a reanalysis of the JLA dataset within a Bayesian formulation, finding significant differences in and values from the original analysis from Betoule et al. (2014). Notably, these methods rely on data that is bias corrected or the methods ignore biases, however the UNITY framework given by Rubin et al. (2015) incorporates selection efficiency analytically in the model, and is applied to the Union 2.1 dataset (Suzuki et al., 2012). The assumption made by the UNITY analysis is that the bias in real data is perfectly described by an analytic function. They validate their model to be free of significant biases using fits to thirty realisations of supernova datasets that are constructed to mimic the UNITY framework. The well known BEAMS (Bayesian estimation applied to multiple species) methodology from Kunz et al. (2007) has been extended and applied in several works (Hlozek et al., 2012), most lately to include redshift uncertainty for photometric redshift application as zBEAMS (Roberts et al., 2017) and to include simulated bias corrections in Kessler & Scolnic (2017). For the latter case, by inferring biases using Bayesian models, sophisticated corrections can be calculated and then applied to more traditional models.

Whilst there are a large number of hierarchical models available, none of them have undergone comprehensive tests using realistic simulations to verify each models’ respective bias. Additionally, testing has generally been performed on supernovae simulations with either CDM cosmology or Flat CDM cosmology. However, quantifying the biases on CDM cosmology simulations with realistic simulations is becoming critically important as precision supernovae cosmology comes into its own, and focus shifts from determination of to both and .

The flexibility afforded by hierarchical models allows for investigations into different treatments of underlying supernova magnitude, color and stretch populations, host-galaxy corrections, and redshift evolution, each of which will be discussed further in the outline of our model below. Our model is designed to increase the numerical efficiency of prior works whilst incorporating the flexibility of hierarchical models. We reduce our dependence on an assumed scatter model in simulations by not utilising bias-corrections in an effort to provide a valuable cross-check on analysis methodologies which utilise scatter-model-dependant bias corrections.

3. Challenges in Supernova Cosmology

The diverse approaches and implementations applied to supernova cosmology are a response to the significant challenges and complications faced by analyses. In this Section, we outline several of the most prominent challenges.

Forefront among these challenges is our ignorance of the underlying type Ia intrinsic dispersion. Ideally, analysis of the underlying dispersion would make use of an ensemble of time-series spectroscopy to characterise the diversity of type Ia supernovae. However this data is difficult to obtain, and recent efforts to quantify the dispersion draw inference from photometric measurements. The underlying dispersion model is not a solved problem, and we therefore test against two dispersion models in this work. The first is based on the Guy et al. (2010, hereafter denoted G10) scatter model, the second is sourced from Chotard et al. (2011, hereafter denoted C11). As the SALT2 model does not include full treatment of intrinsic dispersion, each scatter model results in different biases in , , and when fitting the SALT2 model to light curve observations, and results in increased uncertainty on the summary statistics that is not encapsulated in the reported covariance . These two scatter models are currently assumed to span the possible range of scatter in the underlying supernova population. We have insufficient information to prefer one model over the other, and thus we have to account for both possible scatter models.

The underlying supernova population is further complicated by the presence of outliers. Non-type Ia supernovae often trigger transient follow-up in surveys and can easily be mistaken for type Ia supernovae and represent outliers from the standardised SN Ia sample. This contamination is not just a result of non-SN Ia being observed, but can also arise from host galaxy misidentification causing incorrect redshifts being assigned to supernovae. Different optimizations to the host-galaxy algorithm can result in misidentification of the host at the 3% to 9% level (Gupta et al., 2016), resulting in a broad population of outliers. For spectroscopic surveys, where both supernova type and redshift can be confirmed through the supernova spectra, this outlier population is negligible. However, for photometric surveys, which do not have the spectroscopic confirmation, it is one of the largest challenges; how to model, fit, and correct for contaminants.

Finally, one of the other persistent challenges facing supernova cosmology analyses are the high number of systematics. Because of the rarity of SN Ia events, datasets are commonly formed from the SN Ia discoveries of multiple surveys in order to increase the number of supernovae in a dataset. However, each survey introduces additional sources of systematic error, from sources within each survey such as band calibration, to systematics introduced by calibration across surveys. Peculiar velocities, different host environments, and dust extinction represent additional sources of systematic uncertainty which must all be modelled and accounted for.

4. Our Method

We construct our hierarchical Bayesian model Steve with several goals in mind: creation of a redshift-dependent underlying supernova population, treatment of an increased number of systematics, and analytic correction of selection effects, including systematic uncertainty on those corrections. We also desire Steve more computationally efficient than prior works, such that cosmological results from thousands of supernovae are obtainable in the order of hours, rather than days. As this is closest to the UNITY method from Rubin et al. (2015, hereafter denoted R15), we follow a similar model scaffold, and construct the model in the programming language Stan (Carpenter et al., 2017; Stan Development Team, 2017). The primary challenge of fitting hierarchical models is their large number of fit parameters, and Stan, which uses automatic differentiation and the no-U-turn Sampler (NUTS, a variant of Hamiltonian Monte Carlo), allows us to efficiently sample high dimensional parameter space.

At the most fundamental level, a supernova cosmology analysis is a mapping from an underlying population onto an observed population, where cosmological parameters are encoded directly in the mapping function. The difficulty arises both in adequately describing the biases in the mapping function, and in adding sufficient, physically motivated flexibility in both the observed and underlying populations whilst not adding too much flexibility, such that model fitting becomes pathological due to increasing parameter degeneracies within the model. In the following sections, we will describe these layers, mapping functions, and occurrences of these fatal pathologies. Summaries of observables and model parameters are shown in Table 1 for easy reference.

4.1. Observed Populations

4.1.1 Observables

Like most of the BHM methods introduced previously, we work from the summary statistics, where each observed supernova has a brightness measurement (which is analogous to apparent magnitude), stretch and color , with uncertainty on those values encoded in the covariance matrix . Additionally, each supernova has an observed redshift and a host-galaxy stellar mass associated with it, , where the mass measurement is converted into a probability of being above solar masses. We also have a probability of each supernova being a type Ia, . Our set of observables input into the Steve is therefore given as , as shown in the probabilistic graphical model (PGM) in Figure 2.

As we are focused on the spectroscopically confirmed supernovae for this iteration of the method, we assume the observed redshift is the true redshift such that . Potential sources of redshift error (such as peculiar velocities) are taken into account not via uncertainty on redshift (which is technically challenging to implement) but instead uncertainty on distance modulus. Similarly, we take the mass probability estimate as correct, and do not model a latent variable to represent uncertainty on the probability estimate. One of the strengths of Steve (and the R15 analysis) is that for future data sets where supernovae have been classified photometrically, and we expect some misclassification and misidentification of the host galaxies, these misclassifications can naturally be modelled and taken into account by introducing additional populations that supernovae have a non-zero probability of belonging to.

Parameter Description
Global Parameters
Matter density
Dark energy equation of state
Stretch standardisation
color standardisation
Scale of the mass-magnitude correction
Redshift-dependence of mass-magnitude correction
Systematics scale
Mean absolute magnitude
Survey Parameters
Selection effect deviation
Mean stretch nodes
Mean color nodes
Skewness of color population
Population magnitude scatter
Population stretch scatter
Population color scatter
Extra color dispersion
Redshift-dependence of extra color dispersion
Supernova Parameters
True flux
True stretch
True color
True redshift
Derived absolute magnitude
Derived distance modulus
Input DataaaNot model parameters but shown for completeness.
Measured flux
Measured stretch
Measured color
Covariance on flux, stretch and color
Observed redshift
Determined mass probability
Determined Ia probability
Table 1Model Parameters

4.1.2 Latent Variables for Observables

The first layer of hierarchy is the set of true (latent) parameters that describe each supernova. In contrast to the observed parameters, the latent parameters are denoted without a hat. For example, is the true color of the supernova, whilst is the observed color, which, as it has measurement error, is different from .

For the moment, let us consider a single supernova and its classic summary statistics , , . For convenience, let us define . A full treatment of the summary statistics would involve determining , where represents the observed light curves fluxes and uncertainties. However, this is computationally prohibitive as it would require incorporating SALT2 light curve fitting inside our model fitting. Due to this computational expense, we rely on initially fitting the light curve observations to produce a best fit along with a covariance matrix describing the uncertainty on . Using this simplification, our latent variables are given by


As discussed in Section 3, the SALT2 model does not include full treatment of intrinsic dispersion, and thus this approximation does not encapsulate the full uncertainty introduced from this dispersion.

Figure 2.— Probabilistic graphical model for our likelihood. Double-lined nodes represent observed variables, diamond nodes represent deterministic variables, single-lined ellipse nodes represent fit variables. The SN box represents observed variables and latent variables for each individual supernova, whilst the survey box represents survey-specific variables, which in general describe the supernova population for the survey and the systematics associated with it. Variables that appear outside both boxes represent top level model parameters.

4.2. Underlying Population

4.2.1 Type Ia population

Unlike many previous formalisms which utilise as a singular number and model magnitude scatter on the apparent magnitude , we incorporate this scatter into the underlying rest-frame population by having a population in absolute magnitude space. To denote this difference, we refer to the mean of our absolute magnitude population with .

In addition to absolute magnitude, the underlying supernova population is also characterised by distributions in color and stretch. We follow the prior work of R15 and model the color population as an independent redshift-dependent skew normal distribution for each survey. For the stretch population, we adopt a redshift-dependent normal distribution, and magnitude dispersion is modelled as a normal distribution. We also tested a skew-normal approach for these parameters, reverting to the normal distributions as they are computationally easier to evaluate and we found no reduction in cosmological bias with the skew-normal distributions for stretch and magnitude. Following R15 we allow the mean color and stretch to vary over redshift, anchoring four equally spaced redshift nodes spanning the redshift range of each survey, linearly interpolating between the nodes to determine the mean stretch and color for a given redshift. These nodes are represented as and . Both the color and stretch means are modelled with normal priors. Initial versions of our model adopted a fully covariant multivariate skew normal (with skewness set to zero only for the magnitude component) to capture correlations between and , however pathological fitting complications required us to simplify our treatment. We parametrise the color skewness by sampling which itself is given a uniform prior that allows to span positive values in the realm of physical plausibility as determined from constraints in Scolnic & Kessler (2016). We sample in log space for efficiency in sampling close to zero. The width of the population, represented by the vector is subject to Cauchy priors.

The only constant between survey populations is the absolute magnitude . We model the colour skewness and the redshift-dependent means and width of the colour and skew distributions individually for each survey. The probability for a supernova to have true values , , given the underlying population is thus given as


where for legibility.

4.2.2 Outlier populations

For future use with photometrically classified surveys, we include a simplistic outlier population model. We follow R15 (and therefore Kunz et al., 2007) by implementing a Gaussian mixture model. For surveys with low rates of contamination, it is not possible to fit a contamination population, and the mean of the outlier population has been fixed to the SN Ia population in prior works. However, with the increased number of contaminants expected in the DES photometric sample, we seek a more physically motivated outlier population. We find that an acceptable parametrisation is to model the outlier population with a mean absolute magnitude of , where is constrained to be positive, or even to be greater than a small positive number to reduce degeneracy between the two populations. We note that this represents the mean brightness of outliers, and so outliers could both be brighter and dimmer than the mean SN Ia absolute magnitude. We set the population width to in , and . The probability of each supernova falling into either population is determined by the observed type Ia probability discussed in the previous section. For the spectroscopic survey, we set this to unity. For the photometric proof-of-concept we provide an accurate probability estimate. Further investigation on the effect of inaccurate estimates will be left for future improvements during the analysis of the DES photometric sample.

4.3. Correcting biased summary statistics

With the fitted summary statistics being biased and their uncertainty under-reported, we face a significant challenge utilising these statistics naively in supernova cosmology. We must either correct the observables to remove the biases introduced by the intrinsic dispersion of the underlying population, or incorporate this dispersion into our model. We should also avoid assuming a specific dispersion model – either the G10 or C11 model, or utilise the results of computing the bias from both models.

We model the extra dispersion only in color, and do so by adding independent uncertainty on the color observation . We note that extra dispersion in magnitude (from coherent scatter) is absorbed completely by the width of the underlying magnitude population (discussed in Section 4.2.1) without introducing cosmological bias, which is not true of the color term, hence the requirement for modelling additional color dispersion. Tests on incorporating extra dispersion on stretch as well show that stretch is less biased than color, and causes negligible bias in cosmology.

As shown in (Kessler et al., 2013b), the extra color dispersion shows heavy redshift dependence, increasing with redshift. This is an artifact of different filters, however as we may be subject to similar effects in our observational data, we decide to incorporate redshift dependence in our extra uncertainty. We thus add to our observed color uncertainty (in quadrature). The parameters are highly degenerate with the width of the intrinsic color population . We subject them to Cauchy priors centered on zero and with width , where is bounded between and . We pick this maximum value to allow extra dispersion without completely subsuming the intrinsic population widths due to the severe degeneracy, where this maximum value easily encapsulates the determined dispersion according to the results of Kessler et al. (2013b). As such, our combined covariance on the observation is given by .

Fully covariant extra dispersion on (rather than just dispersion on ) was also tested, by modelling the dispersion as a multivariate Gaussian, but it showed negligible improvement in recovering unbiased cosmology over just color dispersion, and was far more computationally inefficient. We note here that we model dispersion in magnitude, but this is done at the level of underlying populations, not observed populations. This magnitude dispersion is modelled with redshift independence.

4.4. Mapping function

4.4.1 Cosmology

We formulate our model with three different cosmological parameterisations; Flat CDM, Flat CDM, and standard CDM. is given the prior , was treated with and the equation of state was similarly set to a flat prior . For calculating the distance modulus, we fix . If the Hubble constant has a different value, the absolute magnitude is with the other cosmological parameters unaffected.

4.4.2 Supernova Standardisation Parameters

With increasingly large datasets and more nuanced analyses, the choice of how to handle and becomes an important consideration when constructing a model. R15 employs a broken linear relationship for both color and stretch, where different values of and are adopted depending on whether and are positive or negative (although the cut could be placed at a location other than 0). Shariff et al. (2016) instead model as redshift-dependent, testing two phenomenological models; and a second model which effects a rapid but smooth change in at a turnover redshift .

We tested two models with varying against simulated supernova sets; and . See Section 5.2 for details on simulation generation. We found for both models that non-zero values for are preferred even with constant used in simulation, due to severe degeneracy with selection effects. This degeneracy resulted in a significant bias in recovered cosmology. Due to the recovery of non-zero , we continue to adopt the constant and found in traditional analyses. As such, our calculation of distance modulus mirrors that found in Equation (3).

4.4.3 Host Galaxy Environment

There are numerous results showing statistically significant correlations between host-galaxy environment and supernova properties (Kelly et al., 2010; Lampeitl et al., 2010; Sullivan et al., 2010; D’Andrea et al., 2011; Gupta et al., 2011; Johansson et al., 2013; Rigault et al., 2013). The latest sample of over 1300 spectroscopically confirmed type Ia supernovae show evidence for correlation between host mass and luminosity (Uddin et al., 2017). The traditional correction, as employed in analyses such as Suzuki et al. (2012) and Betoule et al. (2014), invokes a step function such that , where is the Heaviside step function, is the galaxy mass in solar masses and represents the size of the magnitude step. The scale of this step function varies from analysis to analysis, and is treated as a fit parameter. In this work we adopt the model used in R15, which follows the work from Rigault et al. (2013), such that we introduce two parameters to incorporate a redshift-dependent host galaxy mass correction:


where represents the correction at redshift zero, and a parameter allowing the behaviour to change with increasing redshift. We take flat priors on and .

4.4.4 Uncertainty Propagation

The chief difficulty with including systematic uncertainties in supernova analyses is that they have difficult-to-model effects on the output observations. As such, the traditional treatment for systematics is to compute their effect on the supernova summary statistics – computing the numerical derivatives , , , where represents the th systematic.

Assuming that the gradients can be linearly extrapolated – which is a reasonable approximation for modern surveys with high quality control of systematics – we can incorporate into our model a deviation from the observed values by constructing a matrix containing the numerical derivatives for the systematics and multiplying it with the row vector containing the offset for each systematic. By scaling the gradient matrix to represent the shift over of systematic uncertainty, we can simply enforce a unit normal prior on the systematic row vector to increase computational efficiency.

This method of creating a secondary covariance matrix using partial derivatives is used throughout the traditional and BHM analyses. For each survey and band, we have two systematics — the calibration uncertainty and the filter wavelength uncertainty. We include these in our approach, in addition to including HST Calspec calibration uncertainty, ten SALT2 model systematic uncertainties, a dust systematic, a global redshift bias systematic, and also the systematic peculiar velocity uncertainty. A comprehensive explanation of all systematics is given in Brout et al. (2018); see Table 4 for details. This gives thirteen global systematics shared by all surveys, plus two systematics per band in each survey. With , our initial conditional likelihood for our observed summary statistics shown in Equation (6) becomes


4.4.5 Selection Effects

One large difference between traditional methods and BHM methods is that we treat selection effects by incorporating selection efficiency into our model, rather than relying on simulation-driven data corrections. We describe the probability that the events we observe are drawn from the distribution predicted by the underlying theoretical model and that those events, given they happened, are observed and pass cuts. To make this extra conditional explicit, we can write the likelihood of the data given an underlying model, , and that the data are included in our sample, denoted by , as


As our model so far describes components of a basic likelihood , and we wish to formulate a function that describes the chance of an event being successfully observed, we rearrange the likelihood in terms of those functions and find


where the denominator represents an integral over all potential data . This is derived in Appendix A.1 As represents the vector of all model parameters, and represents a vector of all observed variables, this is not a trivial integral. Techniques to approximate this integral, such as Monte-Carlo integration or high-dimensional Gaussian processes failed to give tractable posterior surfaces that could be sampled efficiently by Hamiltonian Monte-Carlo, and post-fitting importance sampling failed due to high-dimensionality (a brief dismissal of many months of struggle). We therefore simplify the integral and approximate the selection effects from their full expression in all of -space, to apparent magnitude and redshift space independently (not dependent on or ), such that the denominator of equation (11), denoted now for simplicity, is given as


where can be expressed by translating the underlying , , and population to given cosmological parameters. A full derivation of this can be found in Appendix A.2.

We now apply two further approximations similar to those made in R15 – that the redshift distribution of the observed supernovae reasonably well samples the distribution, and that the survey color and stretch populations can be treated as Gaussian for the purposes of evaluating . We found that discarding the color population skewness entirely resulted in highly biased population recovery (see Figure 12 to see the populations), and so we instead characterise the skew normal color distribution with a Gaussian that follows the mean and variance of a skew normal; with mean given by and variance . This shifted Gaussian approximation for color completely removes the unintended bias when simply discarding skewness. This shift was not required for the stretch population, and so was left out for the stretch population for numerical reasons. The impact of this approximation on the calculated efficiency is shown in Figure 3, and more detail on this shift and resulting population recovery can be found in Appendix A.3.

Figure 3.— Testing the correctness of our normal approximation to the skewed color distribution. The ‘correct’ line (shown in black) represents the exact integral where is an error function (following our high-redshift surveys) and , calculated numerically. The -axis is analogous to is cosmological context. As expected, all efficiencies drop towards zero as shift increases (as objects get fainter). The unshifted normal approximation shows significant discrepancy in the calculated efficiency as it transitions from 1 to 0, whilst the shifted normal approximation shows negligible error to the correct solution. From these plots, further refinement of the normal approximation (such as including kurtosis or higher powers) are unnecessary.

The population becomes , where


What then remains is determining the functional form of . For the treatment of most surveys, we find that the error function, which smoothly transitions from some constant efficiency down to zero, is sufficient. Formally, this gives


where is the complementary cumulative distribution function and and specify the selection function. The appropriateness of an error function has been found by many past surveys (Dilday et al., 2008; Barbary et al., 2010; Perrett et al., 2012; Graur et al., 2013; Rodney et al., 2014). However, for surveys which suffer from saturation and thus rejection of low-redshift supernovae, or for groups of surveys treated together (as is common to do with low-redshift surveys), we find that a skew normal is a better analytic form, taking the form


The selection functions are fit to apparent magnitude efficiency ratios calculated from SNANA simulations, by taking the ratio of supernovae that are observed and passed cuts over the total number of supernovae generated in that apparent magnitude bin. That is, we calculate the probability we would include a particular supernova in our sample, divided by the number of such supernovae in our simulated fields. To take into account the uncertainty introduced by the imperfection of our analytic fit to the efficiency ratio, uncertainty was uniformly added in quadrature to the efficiency ratio data from our simulations until the reduced of the analytic fit reached one, allowing us to extract an uncertainty covariance matrix for our analytic fits to either the error function or the skew normal. This is mathematically identical to fitting the efficiency ratio with a second ‘intrinsic dispersion’ parameter which adds uncertainty to the efficiency ratio data points.

We thus include into our model parametrised selection effects by including the covariance matrix of selection effect uncertainty. Formally, we include deviations from the determined mean selection function parameters with parameter vector , and apply a normal prior on this parameter as per the determined uncertainty covariance matrix. Whilst this uncertainty encapsulates the potential error from the simulations not matching the analytic approximations, it does not cover potential variations of the selection function at the top level — varying cosmology or spectroscopic efficiency. Tests with changing the intrinsic scatter model used in the selection efficiency simulations show that the uncertainty introduced is negligible.

With the well-sampled redshift approximation we can remove the redshift integral in Eq (12) and replace it with a correction for each observed supernova. For the error function (denoted with the subscript ‘CDF’) and skew normal selection functions respectively (denoted with a subscript ‘Skew’), the correction per SN Ia becomes


and is incorporated into our likelihood. This is illustrated in Figure 4. Our corrections for the DES spectroscopic data utilise the CDF functional form, with the combined low redshift surveys being modelled with the skew normal efficiency. Further details on this choice are given in Section 5.2.

Figure 4.— The efficiency for supernova discovery at an arbitrary redshift. Shown in both panels in dashed blue is the SN Ia population distribution, which takes the form of a normal distribution. The top panel shows a CDF based survey efficiency (green dotted line), whilst the bottom panel shows a skew normal based survey efficiency (red dotted line), as functions of apparent magnitude. The survey efficiency, given the SN Ia population, is shown as a solid line in both panels, and the probability of observing a SN Ia is found by integrating over the population detection efficiency as described in equation (12), and has been shown by shading the area integrated. This area is what is analytically given by equations (17) and (18).

5. Model Verification

In order to verify our model we run it through stringent tests. First, we validate on toy models, verifying that we recover accurate cosmology when generating toy supernovae data constructed to satisfy the assumptions of the BHM construction. We then validate our model on SNANA simulations based on a collection of low redshift surveys and the DES three-year spectroscopic sample, termed the DES-SN3YR sample

5.1. Applied to Toy Spectroscopic Data

Figure 5.— Population distributions shown in redshift and uncorrected absolute magnitude for 1000 supernovae in both high-redshift and low-redshift surveys. Selection effects are visible in both samples, where red supernovae are often cut as redshift increases, creating a skewed color population. The color of the data points is representative of the supernova color itself, a negative color value showing bluer supernovae, with positive color values representing redder supernovae.
Figure 6.— Maximal posterior points for 100 realisations of supernova data with the Flat CDM model, with a representative contour from a single data realisation shown for context. Even a large supernova sample when treated robustly is insufficient to provide tight constraints on either or separately due to the severe degeneracy between the parameters.
Figure 7.— Maximal posterior points for 100 realisations of supernova data with the Flat CDM model, with a representative contour from a single data realisation shown for context. The well known banana shaped contour is recovered, with the marginalised distributions in and providing misleading statistics due to the non-Gaussian nature of the posterior surface. The recovered posteior maximums show the same degeneracy direction as the representative surface, and scatter around the truth values input into the simulation, which are shown in dashed lines.
Model (scatter) (scatter)
Flat CDM
Flat CDM

Note. – Cosmological parameters determined from the surfaces of 100 fits to independent realisations of toy supernova data. As described in the main text, each dataset comprised 1000 low-redshift supernovae and 1000 high-redshift supernovae. For each chain, we record the mean and standard deviation, and then show the average mean and average standard deviation in the table. The scatter introduced by simulation variance (the standard deviation of the 100 mean parameter values) is shown in brackets. Model bias would appear as shifts away from the simulation values of , .

Table 2Cosmological Parameters for toy supernova data

We generate simple toy data to validate the basic premise of the model. The data generation algorithm is described below:

  1. Draw a redshift from a power law distribution. For the low redshift survey this is , and for the DES-like survey this is .

  2. Draw a random mass probability from and calculate the mass-brightness correction using , , and equation (8).

  3. Draw an absolute magnitude, stretch and color from the respective distributions , , . Calculate the true absolute magnitude using Equation (1).

  4. Calculate given the drawn redshift and cosmological parameters , under Flat CDM cosmology. Use this to determine the true apparent magnitude of the object .

  5. Determine if the SN Ia is detected using detection probability for the low redshift survey (numeric values obtained by fitting to existing low redshift data). For the DES-like survey, accept with probability . Repeat from step one until we have a supernova that passes. We use realistic values for the selection probability to ensure our model is numerically stable with highly skewed selection functions.

  6. Add independent, Gaussian observational error onto the true , , using Gaussian widths of , , respectively (following the mean uncertainty for DES-like SNANA simulations). Add extra color uncertainty in quadrature of , where .

The selection functions parameters (a skew normal for low-redshift and a complimentary error function for high-redshift) are all given independent uncertainty of (mean and width for the CDF selection function, and mean, width and skewness for the skew normal selection function). Draw from each survey simulation until we have 1000 low- supernovae and 1000 DES-like supernovae, representing a statistical sample of greater power than the estimated 350 supernovae for the DES-SN3YR sample. Sample data for 1000 high and low redshift supernovae are shown in Figure 5, confirming the presence of strong selection effects in both toy surveys, as designed.

We test four models: Flat CDM, Flat CDM, CDM, and Flat CDM with a prior , with the latter included to allow sensitive tests on bias for . To achieve statistical precision, we fit 100 realisations of supernovae datasets. Cosmological parameters are recovered without significant bias. Combined posterior surfaces of all 100 realisations fits for CDM are shown in Figure 6 and fits for Flat CDM are shown in Figure 7. By utilising the Stan framework and several efficient parametrisations (discussed further in Appendix B), fits to these simulations of 2000 supernovae take only on order of a single CPU-hour to run.

To investigate biases in the model in fine detail, we look for systematic bias in in the Flat CDM cosmology test, and bias in for the Flat CDM test with strong prior . This allows us to investigate biases without the investigative hindrances of non-Gaussian or truncated posterior surfaces. The strong prior on cuts a slice through the traditional ‘banana’ posterior surface in the - plane of Figure 7. Without making such a slice, the variation in is larger due to a shift along the degeneracy direction of the ‘banana’. By focusing the slice at an almost fixed , we can see the variation in the mean value of approximately perpendicular to the lines of degeneracy, instead of along them. The results of the analysis are detailed in Table 2, and demonstrate the performance of our model in recovering the true cosmological parameters. As we are using 100 independent realisations, the precision of our determination of the mean simulation result is a tenth of the quoted scatter. The deviation from truth values is below this threshold, no significant bias is detected in either the Flat CDM model or the Flat CDM model. With this simple data, we also correctly recover underlying supernova populations, which can be seen in Figure 12.

5.2. DES SN data validation

Many BHM methods have previously been validated on data constructed explicitly to validate the assumptions of the model. This is a useful consistency check that the model implementation is correct, efficient and free of obvious pathologies. However, the real test of a model is its application to realistic datasets that mimic expected observational data in as many possible ways. To this end, we test using simulations (using the SNANA package) that follow the observational schedule and observing conditions for the DES and low- surveys, where the low- sample is based on observations from CfA3 (Hicken et al., 2009a; Hicken et al., 2009b), CfA4 (Hicken et al., 2012) and CSP (Contreras et al., 2010; Folatelli et al., 2010; Stritzinger et al., 2011). Simulation specifics can be found in Kessler et al. (2018).

Additionally, prior analyses often treated intrinsic dispersion simply as scatter in the underlying absolute magnitude of the underlying population (Conley et al., 2011; Betoule et al., 2014), but recent analyses require a more sophisticated approach. In our development of this model and tests of intrinsic dispersion, we analyse the effects of two different scatter models via simulations, the G10 and C11 models described in Section 3. The G10 models dispersion with 70% contribution from coherent variation and 30% from chromatic variation whilst the C11 model has 25% coherent scatter and 75% from chromatic variation. These two broadband scatter models are converted to spectral energy distribution models for use in simulations in Kessler et al. (2013a).

In addition to the improvements in testing multiple scatter models, we also include peculiar velocities for the low- sample, and our full treatment of systematics as detailed in Brout et al. (2018). Our simulated populations are sourced from Scolnic & Kessler (2016, hereafter SK16) and shown in Table 3. Initial tests were also done with a second, Gaussian population with color and stretch populations centered on zero and with respective width and , however cosmological parameters were not impacted by choice of the underlying population and we continue using only the SK16 population for computational efficiency. The selection effects were quantified by comparing all the generated supernovae to those that pass our cuts, as shown in Figure 8. It is from this simulation that our analytic determination of the selection functions for the low- and DES survey are based. We run two simulations to determine the efficiency using the G10 and C11 scatter models and find no difference in the functional form of the Malmquist bias between the two models.

SK16 low- & &

Note. – The SK16 low- stretch distribution is formed as sum of two bifurcated Gaussians, with the mean and spread of each component given respectively.

Table 3Supernova population distributions
Figure 8.— Fitting the selection function for both the DES 3YR spectroscopically confirmed supernova sample and the low- sample. Blue errorbars represent the efficiency calculated by determining the ratio of discovered to generated supernovae in apparent magnitude bins for SNANA simulations. The black line represents the best fit analytic function for each sample, and the light grey lines surrounding the best fit value represent random realisations of analytic function taking into account uncertainty on the best fit value.

Note. – Standardisation parameters and base intrinsic scatter parameter results for the 100 fits to G10 and C11 simulations. We show the average parameter mean and average standard deviation respectively, with the simulation scatter shown in brackets, such that each cell shows . The width of the intrinsic scatter () does not have an input truth value as it is determined from the scatter model. Model G10 Stat + Syst G10 Stat C11 Stat + Syst C11 Stat

Table 4Realistic simulation standardisation parameters
Figure 9.— Maximum posterior points for 100 realisations of supernova data for two intrinsic dispersion models - the G10 model for the top panel and the C11 model for the bottom panel. Points are shown for parameters , , and , with the other fit parameters being marginalised over. As we are unable to fully correct observed summary statistics, a step required by the lack of intrinsic scatter in the SALT2 model, we expect to see an offset in and . This in turn effects cosmology, resulting in small biases in .

Note. – Investigating the combined 100 fits to G10 and C11 simulations, fitting with both statistics only and also when including systematics. The quoted value for represents the average mean of the fits, with the average uncertainty being shown in square brackets and the simulation scatter (the standard deviation of the mean of 100 fits) shown in standard brackets. The bias significance represents our confidence that the deviation in the mean away from is not due to statistical fluctuation. Model (scatter) -Bias G10 Stat + Syst G10 Stat C11 Stat + Syst C11 Stat

Table 5Realistic simulation determination of

Note. – Correlations determined from the combined 100 simulation fits. Correlations for the low- band systematics and the latent parameters representing selection function uncertainty are not shown but have negligible correlation. Zero superscripts indicate the Dark Energy Survey, and a superscript one represents the low- survey. Parameter G10 Stat+Syst C11 Stat+Syst 0.19 0.21 0.17 0.20 0.29 0.23 0.68 0.66 0.00 0.00 0.00 0.00 0.04 0.07 0.23 0.18 0.04 0.03 0.05 0.01 0.01 0.11 0.08 0.04 0.04 0.04 0.03 0.01 0.10 0.05 0.20 0.17 0.05 0.01 0.01 0.01 0.01 0.05 0.02 0.02 0.04 0.04 0.03 0.06 0.06 0.06 0.04 0.02 0.04 0.04 0.08 0.03 0.05 0.12 0.11 0.03 0.11 0.06 0.14 0.04 0.11 0.11 0.15 0.08 0.12 0.13 0.12 0.06 0.05 0.05 0.01 0.02 0.10 0.09 0.03 0.03 0.08 0.09 0.01 0.02 0.05 0.07 0.11 0.10 0.01 0.02 0.02 0.02 0.03 0.02 0.07 0.07 0.00 0.01 0.01 0.00 0.05 0.11 0.16 0.10 0.16 0.18 0.26 0.26 0.16 0.20 0.05 0.06 0.00 0.01 0.09 0.07

Table 6Reduced parameter correlations with .
Figure 10.— Parameter correlations for the combined fits to the 100 G10 scatter model simulations. We see that the primary correlations with enter through , and , as shown in Table 6. Whilst is generally thought to be a nuisance parameter, we find cosmological correlation. We note that, by fixing in our distance modulus calculation, absorbs any cosmological uncertainty on this term. Additionally also effects the selection efficiency, which was computed from simulations with a fixed value, introducing a second plausible source of correlation. Also visible in this figure are several other interesting relationships. is strongly anti-correlated with intrinsic dispersion for both surveys (DES-like and low-), with showing strong anti-correlation with . This relationship is indeed expected — as grows larger (more unexplained dispersion on the color observation), the width of the supernova population in apparent magnitude space increases. As the fit prefers it to conform to the observed width of the distribution, the extra width in color causes the inherent magnitude smearing amount to decrease. And with extra freedom on the observed color from , shifts in response. The other striking feature in the plot is the strong correlation blocks in the bottom right and the anti-correlation stripes on the edges. These too are expected, for they show the relationship between the color distribution’s mean value, its width and its skewness. As skewness or population width increases, the effective mean of the population shifts (see Appendix A.3 for details), creating anti-correlation between skewness and the (Gaussian) mean color population. Strong anti-correlation between and with reveals the strong population degeneracy, and – for the C11 simulation results – a constrained positive value shows that a finite non-zero extra color dispersion is indeed preferred by our model.

Each realisation of simulated SN Ia light curves contains the SALT2 light-curve fits and redshifts to 128 low- supernovae, and 204 DES-like supernovae, such that the uncertainties found when combining chains is representative of the uncertainty in the DES-SN3YR sample. As our primary focus is Dark Energy, we now focus specifically on the Flat CDM model with matter prior.

Points of maximum posterior for 100 data realisations are shown in Figure 9. The parameter bounds and biases for are listed in Table 5, and secondary parameters are shown in Table 4.

Table 5 shows that the G10 model is consistent with , whilst the C11 model show evidence of bias on , scattering high. However, their deviation from the truth value represents a shift of approximately when taking into account the uncertainty on fits to . The bias is sub-dominant to both the size of the uncertainty for each fit, and the scatter induced by statistical variance in the simulations. We also note that the simulations do not vary cosmological parameters nor population. As our model does include uncertainty on those values, the simulation scatter is expected to be less than the model uncertainty, and represents a minimum bound on permissible uncertainty values.

Table 4 shows a clear difference in both and across the G10 and C11 simulations. As expected, the C11 simulations recover a far smaller intrinsic magnitude scatter, giving results of approximately when compared to the result of for the G10 simulations. The extra smearing of the C11 model does not result in a significantly biased value compared to the average uncertainty on , with recovery of close to the input truth value of , however the recovery for the G10 simulations is biased high, finding with an input of . Interestingly, -bias is only found for the C11 simulations. A measure of the significance of the parameter bias can be calculated by comparing the bias to a tenth of the scatter (as our Monte-Carlo estimate uncertainty is of the scatter). From this, we can see that most biases are detected with high statistical significance due to the large number of simulations tested against.

We investigate the cosmological bias and find its source to be a bias in the observed summary statistics (i.e. the , , output from SALT2 light curve fitting), in addition to incorrect reported uncertainty on the summary statistics. To confirm this, we run two tests. The first of which, we replace the SALT2-fitted , and with random numbers drawn from a Gaussian centered on the true SALT2 , and values with covariance as reported by initial light curve fits. With this test, both the G10 and C11 fits recover exactly. Our second test aims to test our model whilst allowing biases in the summary statistics not caused from intrinsic scatter through. To this end, we test a set of 100 simulations generated using an intrinsic dispersion model of only coherent magnitude scatter, and also find , showing that the source of the biases in summary statistics is the underlying intrinsic scatter model. From this, the main challenge of improving our methodology is to handle the fact that observational uncertainty reported from fitting the SALT2 model to light curves is incorrect, non-Gaussian and biased. Our current model and techniques can quantify the effect of different scatter models on biasing the observed summary statistics, but being unable to constrain the ‘correct’ (simulated) scatter model in our model fit means we cannot fully correct for the bias introduced by an unknown scatter model.

Unfortunately, adding extra fit parameters to allow for shifting observables washes out our ability to constrain cosmology, and applying a specific bias correction requires running a fiducial simulation (assuming cosmology, population and scatter model), which presents difficulties when trying to account for correlations with population and scatter model. This is compounded by the fact that bias corrections do not in general improve fits (increase the log posterior), and so are difficult to fit parametrically. Works such as Kessler & Scolnic (2017) show that bias corrections can be applied to supernovae datasets that can robustly handle multiple intrinsic scatter models, and future work will center on uniting these methodologies — incorporating better bias corrections that separate intrinsic scatter bias and non-Gaussian summary statistic bias from Malmquist bias, without having to precompute standardisation parameters and populations.

Difficulty in providing an adequate parametrisation for realistic intrinsic dispersion, and the simplification of Malmquist bias to only apparent magnitude also leads to biases in the population parameters. Tests when fixing the population parameters to true does not resolve the cosmological bias observed in the C11 simulations. As the population parameters recovered using the simplistic toy supernova data in the previous section do not exhibit significant bias, future work will focus on intrinsic dispersion and Malmquist bias rather than alternate parametrisations of the underlying supernova population.

Table 6 lists the fit correlations between our model fit parameters (excluding the low- band systematics, and Malmquist bias uncertainty parameters which had negligible correlation), showing (in order) cosmological parameters, standardisation parameters, population width and skewness parameters, intrinsic dispersion parameters, mass-step parameters, population mean parameters, SALT2 model systematics, dust systematic, global HST calibration systematic, peculiar velocity systematic, global redshift systematic and DES band magnitude and wavelength systematics. Figure 10 show the full correlations between all non-systematic model parameters. Other interesting correlations are shown and discussed in Figure 10. The band systematics for DES filters , and also show significant correlation with , highlighting the importance of minimising instrumental uncertainty.

For the sample size of the DES low- supernova samples (332 supernova), the bias from intrinsic scatter models is sub-dominant to the statistical uncertainty, as shown in Figure 9. For our full systematics model, the bias represents a deviation between to depending on scatter model, and given that they remain sub-dominant, we will leave more complicated treatment of them for future work.

5.3. Uncertainty Analysis

Note. – Error budget determined from analysing uncertainty on simulation data whilst progressively enabling model features. We start from the top of the table, only varying cosmological parameters and , and then progressively unlock parameters and let them fit as we progress down the table. The cumulative uncertainty shows the total uncertainty on on the fit for all, where the term is derived by taking the quadrature difference in cumulative uncertainty as we progress. Feature Parameters Cumulative Cosmology , 0.051 0.051 Standardisation , , , , 0.046 0.068 Intrinsic scatter , 0.020 0.071 Redshift-independent populations , , , 0.022 0.074 Redshift-dependent populations , 0.030 0.080 Systematics , 0.054 0.096

Table 7 error budget

With the increased flexibility of Bayesian hierarchical models over traditional models, we expect to find an increased uncertainty on parameter inference. To characterise the influence of the extra degrees of freedom in our model, we analyse the uncertainty on averaged across 10 nominal simulations of the DES-SN3YR sample with various model parameters allowed to either vary or stay locked to a fixed value. By taking the difference in uncertainty in quadrature, we can infer the relative contribution for each model feature to the uncertainty error budget.

The error budget detailed in Table 7 shows that our uncertainty is dominated by statistical error, as the total statistical uncertainty is on is . With the low number of supernovae in the DES-SN3YR sample, this is expected. We note that the label ‘Systematics’ in Table 7 represents all numerically computed systematics (as discussed in Section 4.4.4) and systematic uncertainty on the selection function.

5.4. Methodology Comparison

We compare the results of our model against those of the BBC+CosmoMC method (Kessler & Scolnic, 2017). BBC+CosmoMC has been used in prior analyses, such as the Pantheon sample analysis of Scolnic et al. (2017) and is being used in the primary analysis of the DES-SN3YR sample (DES Collaboration et al., 2018). The BBC method is a two-part process, BBC computes bias corrections for observables, and then the corrected distances are fit using CosmoMC (Lewis & Bridle, 2002). For shorthand, we refer to this combined process as the BBC method hereafter in this paper, as we are concerned with the results of cosmological parameter inference. As a leading supernova cosmology method, it provides a good consistency check as to the current levels of accuracy in recovering cosmological parameters.

To this end, we take the results of the BBC method which were also run on the same set of 200 validation simulations and compare the recovered values to those of our method. The results are detailed in Table 8, and a scatter plot of the simulation results is presented in Figure 11.

As shown in Brout et al. (2018), the BBC method recovers cosmological parameters without bias so long as the intrinsic scatter model is known. As we do not know the correct intrinsic scatter model, the BBC method averages the results when using bias corrections from G10 and from C11. As such, we expect the BBC method to have a -bias in one direction for G10 simulations and the other direction for C11 simulations. These results are consistent with those displayed in Table 8. Both BBC method and Steve are sensitive to the intrinsic scatter model, finding differences of and respectively in when varying the scatter model. The BBC method finds biased low for G10 and biased high for C11 (by about ), so taking the average result only results in a small bias of in . Our method shows a small improvement in the insensitivity to the intrinsic scatter model (having a decrease in difference in between the G10 and C11 models), finding no bias for G10 but a biased high for C11. This decrease in error is not statistically significant as we have statistical uncertainty of for simulation realisations. The average bias over the two scatter models is , representing a larger bias than the BBC method.

When comparing both the G10 and C11 set of simulations independently, our model differs from BBC in its average prediction of by and respectively. For the G10 model this difference is a result of bias in the BBC results, however for the C11 simulations this is a result of both bias from BBC, and a larger bias from our method. These results also allow us to state the expected values for when run on the DES-SN3YR sample. When using Planck priors our uncertainty on is reduced compared to using our simulation Gaussian prior on , shrinking the average -difference from to . After factoring this into our uncertainty, we expect our BHM method to, on average, recover .

Having established that our method exhibits similar shifts in the recovery of compared to BBC, future work will focus on improving the parametrisation of intrinsic scatter model into our framework, with the goal of minimising the effect of the underlying scatter model on the recovery of cosmology.

Note. – We characterise the bias on using the 100 simulations for the G10 scatter model and 100 simulations for C11 scatter model. We also show the results when combining the G10 and C11 models into a combined set of 200 simulations. The mean value for our method and BBC are presented, along with the mean when averaging the difference between our method and BBC for each individual simulation. Averages are computed giving each simulation sample the same weight. In the model, represents Steve - BBC. The final row shows the scatter between Steve and BBC for the different simulations. Model G10 C11 (G10 + C11) Steve BBC aaThis value is computed with each simulation having the same weight. It disagrees with the value provided in Brout et al. (2018, Table 10, row 3) which uses inverse variance weighted averages. We do not utilise this weight because the variance is correlated with the value of due to the prior applied in the fitting process. We note that if inverse variance weighting is applied to both datasets, they both shift by , and thus the predicted difference between the BBC method and Steve remains the same.

Table 8 bias comparison
Figure 11.— Recovered for the 200 validation simulations with full treatment of statistical and systematic errors. Uncertainty on the recovered value is shown for every second data point for visual clarity.

6. Conclusions

In this paper we have outlined the creation of a hierarchical Bayesian model for supernova cosmology. The model takes into account selection effects and their uncertainty, fits underlying populations and standardisation parameters, incorporates unexplained dispersion from intrinsic scatter color smearing and incorporates uncertainty from peculiar velocities, survey calibration, HST calibration, dust, a potential global redshift offset, and SALT2 model uncertainty. Furthermore, our uncertainties in standardisation, population, mass-step and more, being explicitly parametrised in our model, are captured with covariance intact, an improvement on many previous methods. The model has been optimised to allow for hundreds of supernovae to be modelled fully with latent parameters. It runs in under an hour of CPU time and scales linearly with the number of supernovae, as opposed to polynomial complexity of matrix inversion of other methods.

The importance of validating models using high-precision statistics gained by performing fits to hundreds of data realisations cannot be overstated, however this validation is lacking in many earlier BHM models for supernova cosmology. We have validated this model against many realisations of simplistic simulations with well-known and well-defined statistics, and found no cosmological bias. When validating using SNANA simulations, we find evidence of cosmological bias which is traced back to light curve fits reporting biased observables and incorrect covariance. Allowing fully parametrised corrections on observed supernovae summary statistics introduces too many degrees of freedom and is found to make cosmology fits too weak. Allowing simulation based corrections to vary in strength is found to give minor reductions in bias, however the uncertainty on the intrinsic scatter model itself limits the efficacy of the bias corrections. For the data size represented in the DES three-year spectroscopic survey, the determined biases should be sub-dominant to other sources of uncertainty, however this cannot be expected for future analyses with larger datasets. Stricter bias corrections calculated from simulations are required to reduce bias. Ideally, this would include further work on the calculation of intrinsic dispersion of the type Ia supernova population such that we can better characterise this bias.

With our model being validated against hundreds of simulation realisations, representing a combined dataset over more than simulated supernovae, we have been able to accurately determine biases in our model and trace their origin. With the current biases being sub-dominant to the total uncertainty, we now prepare to analyse the DES three-year dataset.


Plots of posterior surfaces and parameter summaries were created with ChainConsumer (Hinton, 2016).

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under Grant Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.


Appendix A Selection Effect Derivation

a.1. General Selection Effects

When formulating and fitting a model using a constraining dataset, we wish to resolve the posterior surface defined by


which gives the probability of the model parameter values () given the data. Prior knowledge of the allowed values of the model parameters is encapsulated in the prior probability . Of primary interest to us is the likelihood of observing the data given our parametrised model, . When dealing with experiments that have imperfect selection efficiency, our likelihood must take that efficiency into account. We need to describe the probability that the events we observe are both drawn from the distribution predicted by the underlying theoretical model and that those events, given they happened, are subsequently successfully observed. To make this extra conditional explicit, we write the likelihood of the data given an underlying model, , and that the data are included in our sample, denoted by , as:


A variety of selection criteria are possible, and in our method we use our data in combination with the proposed model to determine the probability of particular selection criteria. That is, we characterise a function , which colloquially can be stated as the probability of a potential observation passing selection cuts, given our observations and the underlying model. We can introduce this expression in a few lines due to symmetry of joint probabilities and utilising that :