Gravitational wave cosmology with extreme mass-ratio inspirals
Danny Laghi, Nicola Tamanini, Walter Del Pozzo, Alberto Sesana, Jonathan Gair, Stanislav Babak
MMNRAS , 1–20 (2021) Preprint 4 February 2021 Compiled using MNRAS L A TEX style file v3.0
Gravitational wave cosmology with extreme mass-ratioinspirals
Danny Laghi (cid:63) , Nicola Tamanini , , Walter Del Pozzo ,Alberto Sesana , Jonathan Gair , Stanislav Babak Dipartimento di Fisica “Enrico Fermi”, Università di Pisa, and INFN Sezione di Pisa, Pisa I-56127, Italy Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, Am Mühlenberg 1, 14476 Potsdam-Golm, Germany Laboratoire des 2 Infinis - Toulouse (L2IT-IN2P3), Université de Toulouse, CNRS, UPS, F-31062 Toulouse Cedex 9, France Dipartimento di Fisica “G. Occhialini”, Università di Milano - Bicocca, Piazza della Scienza 3, 20126 Milano, Italy Laboratoire Astroparticule et Cosmologie, CNRS, Université de Paris, 10 rue Alice Domon et Léonie Duquet, 75013 Paris, France
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The Laser Interferometer Space Antenna (LISA) will open the mHz frequency window of the gravitational wave (GW)landscape. Among all the new GW sources expected to emit in this frequency band, extreme mass-ratio inspirals(EMRIs) constitute a unique laboratory for astrophysics and fundamental physics. Here we show that EMRIs canalso be used to extract relevant cosmological information, complementary to both electromagnetic (EM) and otherGW observations. By using the loudest EMRIs (SNR > H can reach ∼ ∼ Λ CDM parameters fixed by other observations, we further show that in ourbest (worst) case scenario ∼ ∼ w , theDE equation of state parameter. Besides being relevant in their own right, EMRI measurements will be affected bydifferent systematics compared to both EM and ground-based GW observations. Cross validation with complementarycosmological measurements will therefore be of paramount importance, especially if convincing evidence of physicsbeyond Λ CDM emerges from future observations.
Key words:
Gravitational waves astrophysics – Cosmology – Extreme mass ratio inspirals
The first direct detection of gravitational waves (GWs) in2015 (Abbott et al. 2016a) ended a long experimental questand opened a new observational window onto the Universe.Since then, the three ground-based interferometers currentlyoperated by the LIGO-Virgo Collaboration, have announceda total of 50 candidate GW events emitted from binary blackhole (BBH) mergers, binary neutron star (BNS) mergers, andpossible binary neutron star-black hole mergers (Abbott et al.2019a, 2020a,b). These observations have shed new light onthe physics of compact objects (Abbott et al. 2019a,c) andhave allowed tests of General Relativity in the dynamicalstrong field regime for the first time (Abbott et al. 2016b,2019b,d, 2020c). With observational runs currently ongoingand improvements in the sensitivity of the interferometersplanned for the forthcoming years, many more detections areexpected with a progressive increase in the parameter esti-mation accuracy. These new observations are quickly consol-idating GWs into a new observational science, namely GW (cid:63) e-mail:[email protected] astronomy (Barack et al. 2019). Although the direct scopeof GW observations consists in recovering the intrinsic andastrophysical properties of GW sources, interesting cosmolog-ical information can be extracted from the data as well, es-pecially if associated electromagnetic (EM) counterparts canbe identified. Cosmological inference can thus be considered asubpart of GW astronomy, which we can call GW cosmology . The poster example of how GW observations can be usedto infer cosmological constraints is given by GW170817, thefirst BNS merger ever detected (Abbott et al. 2017b). Theobservation of an EM counterpart to the GW signal (Abbottet al. 2017d) allowed the identification of the host galaxy ofthe source, and thus the use of the event as a cosmic dis-tance indicator. In this context such GW sources are usuallyreferred to as standard sirens (Schutz 1986; Holz & Hughes2005; Nissanke et al. 2010). The simultaneous measurement ofboth the luminosity distance (from the GW waveform) andthe redshift (from EM observations) of a GW source pro-vides data points to fit the so-called distance-redshift relation © 2021 The Authors a r X i v : . [ a s t r o - ph . C O ] F e b D. Laghi et al. (Peebles 1993; Weinberg 2008), which links the luminositydistance to the redshift of each point in the Universe and isa function of the cosmological parameters characterizing thecosmic background expansion. At low redshift this relationbecomes the
Hubble law , that only depends on the Hubbleconstant H . Given the low redshift ( z (cid:39) . ) of GW170817,this event provided constraints on H only. The results arein general agreement, though not competitive, with previousmeasurements (Abbott et al. 2017c). Future observations ofsimilar events will reduce the uncertainty on H (Dalal et al.2006; Nissanke et al. 2013; Chen et al. 2018), and possiblywill help solving the tension on its measured value betweenlocal and CMB observations (e.g. Ade et al. (2016); Aghanimet al. (2018); Riess et al. (2016, 2019); Mörtsell & Dhawan(2018); Feeney et al. (2019)).An EM counterpart is not strictly necessary to use compactbinary mergers such as BBHs and BNSs as standard sirens.By matching the sky localisation region of GW sources –which can be inferred from the GW measurements – withgalaxy catalogs, one might in fact be able to extract com-plementary information on the redshift of the sources, with-out the need of spotting an EM counterpart. The idea wasoriginally proposed by Schutz (1986) and it has subsequentlybeen used and developed in different analyses (Holz & Hughes2005; MacLeod & Hogan 2008; Del Pozzo 2012; Petiteau et al.2011; Chen et al. 2018; Gray et al. 2020). It has already beentested with real data collected by the LIGO and Virgo detec-tors (Fishbach et al. 2019; Soares-Santos et al. 2019; Abbottet al. 2020d), though the constraints obtained so far with this“statistical” method are not competitive with the ones derivedfrom GW170817 and its EM counterpart, mainly because ofthe poor spatial resolution of the current network of ground-based interferometers. Future observations, taken with an en-larged network of ground-based GW detectors, will allow forbetter cosmological measurements (Chen et al. 2018), mainlythanks to the improved sky localisation accuracy. Other com-plementary methods, which analogously do not require theidentification of an EM counterpart, might yield interestingresults as well (Taylor & Gair 2012; Del Pozzo et al. 2017;Oguri 2016; Mukherjee & Wandelt 2018; Mukherjee et al.2020c,d,a,b; Farr et al. 2019; Ezquiaga & Holz 2020). Theera of precise cosmological measurements with GWs will how-ever start only with next generation interferometers, suchas the Einstein Telescope (ET) (Punturo et al. 2010; Mag-giore et al. 2020; Sathyaprakash et al. 2010; Belgacem et al.2019a) and the
Cosmic Explorer (Abbott et al. 2017a; Reitzeet al. 2019a,b) on the Earth, or
TianQin (Mei et al. 2020),
Taiji (Luo et al. 2020), and the
Laser Interferometer SpaceAntenna (LISA) (Amaro-Seoane et al. 2017) in space. Thelatter instrument is the focus of the present investigation. Inwhat follows, we will briefly introduce LISA and review pre-vious studies of LISA’s capability to do cosmological analysesusing standard sirens. More details on how to extract cosmol-ogy from GWs by statistically matching with galaxy catalogswill be given in Sec. 2.
LISA is a space mission designed to detect GWs. It has beenselected by the European Space Agency (ESA) for the L3slot of its
Cosmic Vision program, with launch expected inthe early 2030s (Amaro-Seoane et al. 2017). By using inter- ferometric technology already tested in space (Armano et al.2016, 2018), with a much larger arm-length baseline ( . × km) compared to present ground-based detectors ( ∼ few km),LISA aims at measuring GWs in the mHz frequency band,which is expected to be populated by GWs emitted by manydifferent sources.Expected GW sources include: massive black hole binary(MBHB) mergers (Klein et al. 2016), with masses rangingfrom M (cid:12) to M (cid:12) and detectable up to redshift ∼ ∼
100 M (cid:12) , the merger of which will be detectable by ground-based interferometers; extreme mass-ratio inspirals (EMRIs)Babak et al. (2017), which are BBH systems formed by amassive black hole (MBH) and a stellar-origin black hole(SOBH); Galactic and Local-Group binaries (Breivik et al.2018; Korol et al. 2017, 2018; Lau et al. 2019), i.e., compactstellar binaries in the Milky Way and nearby galaxies; andstochastic backgrounds of GWs (Caprini et al. 2016; Bartoloet al. 2016; Caprini & Figueroa 2018), of both cosmologicaland astrophysical origin. The measurement and analysis ofall these sources will yield unprecedented astrophysical andfundamental physics information, and will allow us to testGeneral Relativity in as yet unprobed regimes.The GW sources that LISA will observe at cosmologicaldistances can be used as standard sirens. These include MB-HBs, EMRIs, and SOBHBs. Unfortunately, only for the firstof these types of sources are EM counterparts plausibly ex-pected to be produced and observed by future EM facilities(Tamanini et al. 2016). MBHBs are in fact expected to emita large amount of EM radiation in different bands at mergeror during long-lasting ( ∼ weeks/months) afterglows (see, e.g.,Palenzuela et al. 2010; Dotti et al. 2012; Giacomazzo et al.2012; Moesta et al. 2012), and possibly even through pre-merger signals (Kocsis et al. 2008; O’Shaughnessy et al. 2011;Kaplan et al. 2011; Haiman 2017; Dal Canton et al. 2019). Ifsufficiently accurate sky localisation can be attained from theGW parameter estimation analysis and if the EM counterpartis sufficiently powerful to be spotted by EM telescopes, thenwe expect to identify the host galaxy of up to a few LISAMBHB mergers per year (Tamanini et al. 2016; Tamanini2017). These golden sources can then be used as high-redshiftstandard sirens to map the expansion of the Universe up to z ∼ . Although the low number of expected EM counter-parts and the high redshift of MBHB mergers are not idealto test standard cosmological models such as Λ CDM or toplace constraints on late-time dark energy (DE) (Tamaniniet al. 2016; Tamanini 2017; Belgacem et al. 2019b), they canefficiently be used to probe deviations from Λ CDM at earliercosmological epochs, specifically in the interval (cid:46) z (cid:46) (Caprini & Tamanini 2016; Cai et al. 2017; Belgacem et al.2019b; Speri et al. 2020).At lower redshift, LISA will provide other GW sources thatcan be used as standard sirens. SOBHBs will be mainly de-tected at redshifts z (cid:46) . , while EMRIs might be observedup to z ∼ , with a broad peak around z ≈ . Unfortunately,the most widely accepted formation channels for these typesof sources do not predict associated EM counterparts (see,e.g., Belczynski et al. 2002; Amaro-Seoane et al. 2007; Ro-driguez et al. 2016), implying that statistical matching withgalaxy catalogs will be necessary in order to extract cosmo-logical information from them. A few investigations have al- MNRAS000
100 M (cid:12) , the merger of which will be detectable by ground-based interferometers; extreme mass-ratio inspirals (EMRIs)Babak et al. (2017), which are BBH systems formed by amassive black hole (MBH) and a stellar-origin black hole(SOBH); Galactic and Local-Group binaries (Breivik et al.2018; Korol et al. 2017, 2018; Lau et al. 2019), i.e., compactstellar binaries in the Milky Way and nearby galaxies; andstochastic backgrounds of GWs (Caprini et al. 2016; Bartoloet al. 2016; Caprini & Figueroa 2018), of both cosmologicaland astrophysical origin. The measurement and analysis ofall these sources will yield unprecedented astrophysical andfundamental physics information, and will allow us to testGeneral Relativity in as yet unprobed regimes.The GW sources that LISA will observe at cosmologicaldistances can be used as standard sirens. These include MB-HBs, EMRIs, and SOBHBs. Unfortunately, only for the firstof these types of sources are EM counterparts plausibly ex-pected to be produced and observed by future EM facilities(Tamanini et al. 2016). MBHBs are in fact expected to emita large amount of EM radiation in different bands at mergeror during long-lasting ( ∼ weeks/months) afterglows (see, e.g.,Palenzuela et al. 2010; Dotti et al. 2012; Giacomazzo et al.2012; Moesta et al. 2012), and possibly even through pre-merger signals (Kocsis et al. 2008; O’Shaughnessy et al. 2011;Kaplan et al. 2011; Haiman 2017; Dal Canton et al. 2019). Ifsufficiently accurate sky localisation can be attained from theGW parameter estimation analysis and if the EM counterpartis sufficiently powerful to be spotted by EM telescopes, thenwe expect to identify the host galaxy of up to a few LISAMBHB mergers per year (Tamanini et al. 2016; Tamanini2017). These golden sources can then be used as high-redshiftstandard sirens to map the expansion of the Universe up to z ∼ . Although the low number of expected EM counter-parts and the high redshift of MBHB mergers are not idealto test standard cosmological models such as Λ CDM or toplace constraints on late-time dark energy (DE) (Tamaniniet al. 2016; Tamanini 2017; Belgacem et al. 2019b), they canefficiently be used to probe deviations from Λ CDM at earliercosmological epochs, specifically in the interval (cid:46) z (cid:46) (Caprini & Tamanini 2016; Cai et al. 2017; Belgacem et al.2019b; Speri et al. 2020).At lower redshift, LISA will provide other GW sources thatcan be used as standard sirens. SOBHBs will be mainly de-tected at redshifts z (cid:46) . , while EMRIs might be observedup to z ∼ , with a broad peak around z ≈ . Unfortunately,the most widely accepted formation channels for these typesof sources do not predict associated EM counterparts (see,e.g., Belczynski et al. 2002; Amaro-Seoane et al. 2007; Ro-driguez et al. 2016), implying that statistical matching withgalaxy catalogs will be necessary in order to extract cosmo-logical information from them. A few investigations have al- MNRAS000 , 1–20 (2021) osmology with EMRIs ready assessed the potential of LISA to test the Hubble law bystatistically matching the sky localisation region of SOBHBswith galaxy catalogs (Kyutoku & Seto 2017; Del Pozzo et al.2018). It was found that constraints on the Hubble constantcan reach at best a few %, though uncertainties on the ex-pected sensitivity of LISA at high-frequency could well under-mine this result (Moore et al. 2019). No thorough investiga-tion has so far been performed considering EMRIs as possiblestandard siren sources. The only analysis that can be foundin the literature (MacLeod & Hogan 2008) suggests that ∼ z ∼ . could be used to constrain theHubble constant at the 1% level. However, beside consideringan old configuration of LISA, this study was highly idealized:the authors assumed only linear cosmic expansion neglectingthe acceleration of the Universe, they employed a simplifiedstatistical framework and did not perform parameter estima-tion over the GW signals, using approximate relations only.A complete cosmological investigation with LISA EMRIs iscurrently missing in the literature. The scope of the present investigation is to provide an in-depth and up-to-date analysis of the prospects of using EM-RIs detected by LISA as GW standard sirens. As alreadystressed and cited (above), the standard capture scenario togenerate EMRIs does not predict observable EM counter-parts. We thus employ a statistical method, assigning to eachgalaxy within the LISA 3D error volume a probability of be-ing the host of the EMRI within a Bayesian framework. Wewill start by presenting the statistical methodology that weuse to infer cosmological parameters by cross-matching GWsky localisation regions with galaxy catalogs (Sec. 2). We willthen describe the catalogs of EMRIs observed by LISA thatwe will use in our analysis (Sec. 3), and outline the procedurefor matching sky localisations with galaxy catalogs (Sec. 4).We will subsequently present the results of our investigation(Sec. 5), discuss them (Sec. 6), and finally we will draw ourconclusions (Sec. 7).
To infer the value of the cosmological parameters, we operatewithin the framework of Bayesian inference (Jaynes 2003).The starting point of our analysis is: i) a list of GW EMRIobservations D , combined with ii) a catalog of galaxies thatare potential hosts of each individual EMRI. We will describethe EMRI and galaxy catalogs in Sec. 3 and 4, respectively.Here we provide the mathematical details of the analysis,given those populations.Our inference model is constructed starting from Del Pozzo(2012), with updates from Del Pozzo et al. (2017); Del Pozzoet al. (2018) to account for the uncertainty on the redshift ofpotential counterparts to GW events. We will, however, usea different notation from both references. In the Bayesianframework, all knowledge about the cosmological parameterswe are interested in, Ω ≡ { H , Ω m , w , w a } , is summarisedby the posterior probability distribution, p (Ω | D H I ) = p (Ω | H I ) p ( D | Ω H I ) p ( D | H I ) , (2.1) where H is the cosmological model, that defines the relationbetween redshift, distance, and cosmological parameters, I represents all background information relevant for the infer-ence of Ω , and D ≡ { D , . . . , D N } is the set of GW ob-servations, with D i the data from the i ’th EMRI event. Ina realistic setting, this data would be the strain time seriescorresponding to the event observed by the LISA detector, sodifferent EMRI events can overlap in time and therefore sharethis time series. For the current analysis, we will howeverrepresent the observed data instead by D i = { ˆ d L , ˆ θ gw , ˆ φ gw } i ,that is, point estimates of the event luminosity distance andsky position in ecliptic coordinates, with associated uncer-tainties, and we assume these estimates for each event areindependent.The terms on the right hand side of Eq. (2.1) are the prior probability distribution p (Ω | H I ) , the likelihood func-tion p ( D | Ω H I ) , and the evidence p ( D | H I ) . Since we arenot interested in doing model selection, at this stage the ev-idence is considered as a normalization constant for the pos-terior, so we will only define the prior and the likelihood,eventually renormalizing Eq. (2.1) at the end of the compu-tation.To start with, we formalize our prior information and cos-mological model in the following propositions: I : “A GW event can be hosted by only one galaxy. Notall the galaxies within the comoving volume are visible. Ob-served and non-observed galaxies obey the same cosmol-ogy. Galaxies are highly clustered with each other. The hostgalaxy catalog is considered to be complete”. H : “The cosmological model obeys a flat ( k = 0 ) Friedmann-Lemaítre-Robertson-Walker (FLRW) metric (Weinberg1972) with negligible energy contribution from radiation( Ω r = 0 ), predicting that the luminosity distance d L isa function of the redshift z gw and of the cosmologicalparameters Ω : d L = d (Ω , z gw ) (Peebles 1993; Hogg 1999;Weinberg 2008)”.With assumption I , we are stating that the catalogue ofgalaxies is representative of the whole galaxy distribution,and it is equivalent to assuming that the EMRI lies in oneof the galaxies in the catalogue. EMRIs hosted in galaxiesnot in the catalogue are assumed to be sufficiently close togalaxies in the catalogue that this assumption holds. The like-lihood for every event is determined by the same catalogueof galaxies, but for ease of computation we only use a sub-set of possible galaxies to analyse each event. Possible hostsare any galaxies that lie within the 2 σ credible region for theGW direction and distance, for at least one value of the cos-mological parameters within the prior range. This selects a3D co-moving volume in direction and redshift (referred toas an “error box”) containing N g,i possible hosts, each withsky position ( θ j , φ j ) and redshift z j .We will further specify our model H investigating two dif-ferent cosmological scenarios , that will define the subset of Ω which we will be interested in. In this work we consider twomodels:(i) A Λ CDM scenario characterized by a parameter space ( h, Ω m ) to be explored ( h ≡ H / km − s Mpc, while Ω Λ is determined through the boundary condition Ω m + Ω Λ = 1 ,which holds assuming spatial flatness). We consider a priorrange h ∈ [0 . , and Ω m ∈ [0 . , . , with fiducial values MNRAS , 1–20 (2021)
D. Laghi et al. dictated by the Millennium run (Springel et al. 2005) (i.e. h =0 . , Ω m = 0 . , Ω Λ = 0 . ).(ii) A dark energy (DE) scenario in which we assume ( h, Ω m , Ω Λ ) to be pre-determined by other probes at thevalues of the Millennium run (Springel et al. 2005) and wesearch on the parameters ( w , w a ) defining the DE equationof state (Linder 2003) w ( z ) = w + z/ (1 + z ) w a . We drawfrom the prior range w ∈ [ − , − . , w a ∈ [ − , , withfiducial values w = − and w a = 0 , corresponding to thecosmological constant Λ .We assume each GW event is statistically independent fromany other event. Astrophysically, every EMRI event occursindependently of the others since we will most likely neverobserve multiple EMRIs occurring in the same galaxy. Thereis, however, coupling between parameter estimation for theevents because they will be overlapping in the LISA dataset and hence subject to the same noise fluctuations. Weexpect the observed EMRIs to be essentially orthogonal toeach other and therefore they will be susceptible to differentcomponents of the instrumental noise. So, mutual indepen-dence should be a good approximation. Thus, the likelihoodfunction simplifies to the product of the likelihoods for eachindividual observation, p ( D | Ω H I ) = N (cid:89) i =1 p ( D i | Ω H I ) . (2.2)Therefore, we only need to determine how to construct thelikelihood for a single GW event, whose main ingredients weintroduce in what follows. The relevant quantities for the in-ference of the cosmological parameters are the GW luminos-ity distance d L (which is directly measured) and redshift z gw (which is inferred through the cosmology priors) . Assumingthat the correlation between sky localisation and other pa-rameters is negligible, a multiple application of marginalisa-tion and product rule lead to a single-event likelihood writtenas: p ( D i | Ω H I ) = (cid:90) d d L d z gw p ( d L | z gw Ω H I ) · p ( z gw | Ω H I ) p ( D i | d L z gw Ω H I ) , (2.3)where the integrals on d L and z gw go, in principle, from 0 to ∞ .Once we know the redshift z gw along with the values ofthe cosmological parameters (assuming their priors) and thecosmological model H , we have (Del Pozzo 2012) p ( d L | z gw Ω H I ) = δ ( d L − d (Ω , z gw )) , (2.4)where under our background information I the function d (Ω , z gw ) is given by (Hogg 1999): d (Ω , z gw ) = c (1 + z gw ) H (cid:90) z gw d z (cid:48) E ( z (cid:48) ) , (2.5)with the function E ( z (cid:48) ) given by: E ( z (cid:48) ) = (cid:112) Ω m (1 + z (cid:48) ) + Ω Λ g ( z (cid:48) , w , w a ) , (2.6) Although these values of cosmological parameters are outdated,this is irrelevant for the purpose of testing our statistical method-ology for constraining cosmology with standard sirens. To ease the notation, when there is no possibility of confusionwe will be dropping the EMRI index i , as in the quantities d L and z gw . where g ( z (cid:48) , w , w a ) = (1 + z (cid:48) ) w + w a ) e − w z (cid:48) z (cid:48) , (2.7)and as already noted in the definition of H we neglect theradiation density Ω r , since Ω r (cid:28) Ω m , Ω Λ .Our prior information I prescribes that each GW is hostedby a galaxy; we are assuming that the distribution of z gw isthe same as that of the host galaxies. Thus we have p ( z gw | Ω H I ) ∝ N g,i (cid:88) j =1 w j exp (cid:40) − (cid:18) z j − z gw σ z j (cid:19) (cid:41) , (2.8)where we have introduced relative weights w j for each galaxyhost (this is different from MacLeod & Hogan (2008), whereall galaxies within the error-box were assumed to be equallylikely hosts). The uncertainty σ z j includes the contributiondue to the peculiar velocity of the galaxy j , and the sum goesover all the possible galaxy hosts N g,i associated to the i -thevent. As in Del Pozzo et al. (2018), we assume a typical red-shift error of ∆ z = 0 . corresponding to a rms peculiarvelocity of 500 km/s. We assign the relative weights w j toeach galaxy from the marginal distribution over the sky posi-tion angles computed in each galaxy, as defined by Eq. (4.4)and described in detail in Sec. 4. We note that the product w j p ( D i | d L z gw Ω H I ) should equal the full likelihood evalu-ated at the sky location of the galaxy, which is assumed tobe known perfectly. Therefore w j depends on both d L andon the observed data, which is not explicit in the preced-ing equation. Here we approximate w j as a constant, givenby marginalising the posterior over luminosity distance atthe sky location of the galaxy, and similarly assume that theprojected distance posterior does not depend on sky loca-tion. This is equivalent to assuming that the joint likelihoodon sky location and distance factorises. Assigning weights inthis way we are treating all galaxies as a priori equally likelyhosts of observed EMRIs, but weightings can also be used toreflect the relative probability that a galaxy is the GW host,for example based on the galaxy luminosity. The expressionslook the same, but the weights are then proportional to theproduct of the sky-location contribution and the pre-assignedweight.The remaining term in Eq. (2.3) is: p ( D i | d L z gw Ω H I ) = p ( D i | d L z gw H I ) , (2.9)which is the likelihood for the GW data. As described above,we approximate the GW data as point estimates of the rele-vant parameters, with Gaussian uncertainties, and so replacethis term by a quasi-likelihood specified as a multivariateGaussian, with covariances estimated from the Fisher Matrix(see Sec. 3), and also accounting for the weak lensing uncer-tainty σ WL,i as modelled in Tamanini et al. (2016), see Eq.(7.3) therein (see also Hirata et al. (2010); Cusin & Tamanini(2020)).The final expression for the single-event likelihood is: p ( D i | Ω H I ) = 12 π (cid:90) d z gw,i N g,i (cid:88) j =1 w j σ z j (cid:113) σ d L,i + σ WL,i · exp (cid:40) − (cid:34) (cid:0) z j − z gw,i (cid:1) σ z j + (cid:0) ˆ d L,i − d (Ω , z gw,i ) (cid:1) σ d L,i + σ WL,i (cid:35)(cid:41) . (2.10)As mentioned above, ( ˆ d L , cos ˆ θ gw , ˆ φ gw ) i are the point esti-mates of the parameters, and the uncertainties, σ d L,i etc.,
MNRAS , 1–20 (2021) osmology with EMRIs are determined from the Fisher matrix evaluated at the trueparameter values. We have also assumed that the joint like-lihood on sky location and distance factorises as mentionedabove.Until now, we have assumed that all EMRIs that occur dur-ing the observation period are included in the analysis. How-ever, in practice our detectors have limited sensitivity and wewill only include in the analysis the events that are success-fully “detected” by our analysis pipelines, that is, those forwhich some ranking statistic computed from the data is abovesome predetermined threshold. Whether or not an event isfound is a property of the data only, and assuming that thetotal number of events does not convey any information theselection effects can be accounted for (Mandel et al. 2019) byreplacing p ( D i | Ω H I ) with ˆ p ( D i | Ω H I ) = p ( D i | Ω H I ) (cid:82) D> threshold p ( D i | Ω H I ) d D , (2.11)where the integral is over all data sets that would give de-tection statistics above the threshold and hence be includedin a cosmological analysis. In practice, the detectability ofa GW event is primarily determined by its signal-to-noiseratio (SNR), which in turn is primarily determined by lumi-nosity distance. In the quasi-likelihood model used here wecould thus use a cut on the observed luminosity distance, ˆ d L > d crit , as a proxy to represent selection effects. As thechoice of cosmological parameters is varied, the centres ofthe Gaussian distributions in the likelihood will move insideand outside the selection cut. Although a small number of theGaussian distributions span the boundary, the selection func-tion normalisation can be well approximated by the numberof host galaxies that remain consistent with being inside theluminosity distance horizon. In our analysis we found thatthis number changed by very little over the range of cos-mological parameters consistent with our priors, and so weapproximated this normalisation as a constant. However, wealso verified that our results were insensitive to that approx-imation.This completes the definition of our inference problem. Wewill further discuss the limitations and approximations of ourframework in Sec. 3, 7, and Appendix A.We explore the posterior distribution in Eq. (2.1) us-ing cpnest (Veitch et al. 2020), a parallel nested sam-pling algorithm implemented in Python . We make use ofthe
LALCosmologyCalculator library in
LAL (LIGO ScientificCollaboration 2018) and of numpy (van der Walt, Stéfan et al.2011), scipy (Virtanen et al. 2020), matplotlib (Hunter2007), cython (Behnel et al. 2011), and the plotting utilitiesin Foreman-Mackey (2016). The inference code utilised in thispaper is publicly available in Del Pozzo & Laghi (2020).
The catalogs of EMRIs detected by LISA which we use inour work are based on the analysis of Babak et al. (2017).Here we review how these catalogs have been constructedand outline their main properties. For more information thereader is referred to Babak et al. (2017) (see also Babak et al.2015; Gair et al. 2017).Because of their small mass-ratio, EMRIs present a very slow inspiral, producing many cycles ( - ) within LISA’ssensitivity band. For the same reason, the detailed dynam-ics of these systems strongly differ from equal mass compactbinaries. The motion of the stellar BH can in fact be approx-imated by geodesics of the massive BH, with small but rele-vant corrections due to its own self-force (for a review, see,e.g., Poisson et al. 2011). Unfortunately ongoing perturbativecalculations, which exploit the extreme difference in massesof these systems, are not yet at the level needed to calculatethe emitted GW waveform with the accuracy required forLISA observations (i.e., by keeping track of the phase overthe entire inspiral). The analysis of Babak et al. (2017) con-sidered thus an approximate analytical model to estimate thewaveform generated by EMRIs, the so-called analytical kludge model (Barack & Cutler 2004). This approximation was usedto produce waveforms under two different endpoints for thedynamical evolution of the system: until the Schwarzschildlast stable orbit (more pessimistic) and until the Kerr laststable orbit (more optimistic). In our study we will only con-sider results obtained with the latter of these assumptions(denoted “AKK” in Babak et al. (2017)).In order to simulate the response of LISA and estimate theuncertainty on the parameters of the GW waveform gener-ated by an EMRI, the investigation made by Babak et al.(2017) employed a Fisher matrix approach, useful to quicklyanalyse many different events. Their results were obtainedsetting a threshold of SNR=20 for LISA detection and con-sidering a LISA sensitivity curve as specified by the LISAWhite Paper written in response to ESA’s call for the L3mission slot (Amaro-Seoane et al. 2017). The results of thispaper will thus be based on that sensitivity curve as well. Wenote here that although the laser stability requirement hasbeen relaxed since, this has an impact to the detector per-formance at f > . Hz, which is unlikely to significantlyaffect EMRIs.As pointed out by Babak et al. (2017), the main uncer-tainties in forecasting how many EMRIs will be detected byLISA are of astrophysical origin. The expected rate of EMRIsdepends in fact on several astrophysical assumptions, includ-ing: • the MBH population in the accessible LISA mass range(from M (cid:12) to M (cid:12) ), the redshift evolution of their massfunction, and their spin distribution; • the fraction of MBHs hosted in dense stellar cusps, whichconstitute the nurseries for the formation of EMRIs; • the EMRI rate per individual MBH, and the mass andeccentricity distribution of the inspiralling compact object.In the analysis of Babak et al. (2017), by considering 12 differ-ent combinations of prescriptions for the assumptions listedabove, the authors produced forecasts for 12 different scenar-ios. Among all 12 scenarios, the rates for LISA detections us-ing the AKK model always fall in between 10% to 20% of thetotal EMRI population, and in absolute numbers they span arange from one to a few thousands of detections per year, withsimilar properties between the 12 EMRI populations corre-sponding to each scenario. In particular, irrespective of theastrophysical scenario, EMRIs detected by LISA will comefrom MBHs with masses from × M (cid:12) to × M (cid:12) , overa redshift range that is broadly peaked between . < z < ,with tails usually reaching z ∼ .In this paper we concentrate on three EMRI models, span- MNRAS , 1–20 (2021)
D. Laghi et al. ning the bulk of the uncertainty range in the number ofEMRI detections, corresponding to models M1, M5, and M6of Babak et al. (2017) (cf. their Table I): • our fiducial model is based on M1, which depends onrather standard assumptions for the MBH mass function andEMRI properties; • our pessimistic model is based on M5, featuring a down-turn of the MBH mass function at M < M (cid:12) , whichstrongly suppresses the occurrence of EMRIs; • our optimistic model is based on M6. This is similar toM1, but ignores the effect of cusp erosion following galaxymergers, which boosts the EMRI rate by a factor of ∼ .Although there are more optimistic and pessimistic EMRImodels in Babak et al. (2017), those are based on rather ad-hoc assumptions, in particular about the relative occurrenceof plunges and EMRIs. In fact, when a BH is scattered on analmost radial orbit, chances are that it plunges directly ontothe MBH, rather than forming an EMRI. N-body simula-tions of realistic stellar cusps suggest that there are few suchdirect plunges for each EMRI (Merritt 2015). The modelslisted above assume a rather conservative ratio of 10 plungesper EMRI. We note that, as the number of useful events N increases further, one can reasonably expect that our final re-sults, i.e., the constraints that we obtain on the cosmologicalparameters, scale approximately with √ N (cf. Sec. 6).Each EMRI model predicts different events, thus it can bethought of as a different data set D . For each data set, wewill adopt the cosmological inference models H detailed inSec. 2. Because of the complexity of their waveform, a relatively largeSNR is required for EMRI detection. As already noted, thethreshold is customarily set to SNR = 20 , which results in sev-eral hundreds-to-thousands of EMRIs detected over 10 yearsfor the three models considered here, as reported in the sec-ond column of Table 1.Although such an abundance of sources is in principle anasset for our analysis, a number of considerations have to bemade when designing a practical implementation of the algo-rithm. In fact, the speed of computation of Eq. (2.10) stronglydepends on the dimension of the EMRI catalog. This imposesa limitation on the number of events that we can analyse ina reasonable time. At the current state, catalogs of hundredsof EMRIs are prohibitive, as they require O ( months ) to beanalysed (especially with the limited computational resourcesavailable to this project). This is due to two main reasons:i) the increase in the number of single-event likelihoods thathave to be computed and ii) the large number of hosts N g,i that some events will have, resulting in a significant slowdownof each likelihood evaluation.Point ii) suggests some useful guidelines for devising anEMRI selection strategy. As shown for our fiducial model M1by the green histograms in Fig. 1, the majority of EMRIswill be observed at z > . As z increases, so do the errorsin the estimate of their parameters and the associated LISA3D localisation volume ∆ V . When ∆ V gets too large, thenumber of candidate galaxy hosts within it, N g,i , increasesdramatically, washing out the information enclosed in their clustering properties, which is of paramount importance forthe success of our technique.This suggests two approaches to the selection of the events.One possible strategy, that we indicate as our first selectionprocedure, is to pre-select EMRIs with a good LISA 3D skylocation estimate up to some maximum redshift. Followingthis route, we select only EMRIs with sky localisation error ∆Ω = 2 π sin θ gw (cid:112) ∆ θ gw ∆ φ gw − (Σ θ gw φ gw ) < deg , being θ gw and φ gw the latitude and longitude of the source in eclip-tic coordinates and Σ θ gw φ gw their correlation extracted bythe 3D marginalised correlation matrix of the Fisher analysis.We further require a relative error in the luminosity distancedetermination: ∆ d L /d L < . . Moreover, an additional con-straint stems from the fact that, in order to keep the numberof candidate hosts within a maximum of a few thousands,we use the mock Millennium sky to select galaxy hosts up to z = 1 . This implies that we keep only those events for whichthe furthest candidate host compatible with the adopted pri-ors in the cosmological parameters is at z < . In fact, thecut in our galaxy catalogue at z = 1 implies that we discardsome of the events at z (cid:38) . for which the furthest candidatecompatible with the range of cosmological priors lies at z > .We note that this cut depends on the adopted cosmologicalmodel, and yields a slightly different selection for the Λ CDMand the DE cases. The number of events selected accordingto this procedure for the Λ CDM and DE models are given inthe third and fourth columns of Tab. 1.Although legitimate, this selection procedure is rather con-trived, as it combines a number of ad hoc choices in the selec-tion process. A more straightforward alternative is to observa-tionally select events by simply imposing an SNR threshold.This is illustrated for our fiducial model M1 by the orangehistograms in Fig. 1. In fact, high-SNR events are expected tobe those yielding a better determination of the EMRI param-eters (distance and sky location) and are also preferentially atrelatively low redshifts: this is evident looking at panels a), b),and c). This combination of properties automatically limitsthe number of candidate hosts – see panel d) – making themthe most informative sources for our analysis. We found em-pirically that imposing a LISA SNR threshold SNR >
100 re-turns only events with N g (cid:46) out to z ≈ . , as shown inpanel a). The imposition of an SNR threshold, together withthe cut at z = 1 dictated by the galaxy catalog that is subse-quently used (which is sensitive to the assumed cosmologicalmodel), yields the number of events reported in the fifth andsixth columns of Table 1. The first selection procedure pro-duces a larger number of EMRIs (see the third and fourthcolumns of Tab. 1), among which we found several eventswith relatively low SNR and N g > . Tests performed onthese samples produced reasonable estimates of the cosmolog-ical parameters, but also displayed some undesired behaviourrelated to the limitations of our likelihood implementation,which became particularly evident when the majority of theselected events had low SNR. These aspects will be furtherdiscussed later on and in particular in Appendix A2. There-fore, unless otherwise stated, results reported in this workfollow the SNR >
100 selection criterion.
MNRAS , 1–20 (2021) osmology with EMRIs detected(SNR > ∆ d L /d L < . SNR > & ∆Ω < deg Λ CDM DE Λ CDM DE
M1 ( fid ) 2941 180 202 32 33M5 ( pess ) 472 20 21 6 6M6 ( opt ) 4788 260 281 72 73 Table 1.
Number N of EMRIs observed by LISA in 10 yearsof operation. For each model (first column) we report the totalnumber of detections (second column), those additionally local-ized within ∆ d L /d L < . < deg , referred in the textas first selection procedure, condition that depends on the chosencosmological model (third and fourth columns), and those withSNR > (fifth and sixth columns). In both selection procedureswe require no potential host galaxy above z = 1 . Figure 1.
Relevant properties of the events used in the analysis,taken from model M1. Green histograms in panels a), b), and c)show the redshift, distance error and sky location error distribu-tions for the entire population of EMRIs detected by LISA in 10years of operations. For the purpose of our analysis we select onlysystems with SNR > σ er-ror volume, averaged over the three realisations of M1 (see section4 for details). Selecting only high-SNR events allows us to avoid some lim-itations and complications related to a number of featuresthat are not (as of yet) modeled in our pipeline.First, we do not fold into our analysis the possible incom-pleteness of the galaxy catalogs. As we will see in the nextsection, galaxy catalogs are constructed starting from theflux-limited mock sky realizations of Henriques et al. (2012)based on the Millennium simulation (Springel et al. 2005). Asreported in our prior information I , we are going to assumethat each GW event is hosted by one of the galaxies withinthe reconstructed co-moving volume, regardless of whetherthe galaxy is luminous enough to be listed in the mock cat-alog or not. We assume that EMRIs are hosted by galaxies with stellar mass M ∗ > M (cid:12) , for which we found thatthe Millennium sky maps are complete out to z ≈ . . Thelarge majority of events we select are at z < . (see Fig. 1).Thus, the mass threshold imposed in constructing our cata-logs should not yield a relevant Malmquist bias.Another possible source of bias may be given by evolutioneffects, i.e, the fact that the host galaxy population evolveswithin the allowed redshift range for the host. Also in thiscase, the high-SNR cut preferentially selects low- z events withprecise distance measurement, thus limiting the redshift ofcandidate hosts within a window that is narrow enough thatevolution effects can be neglected.Last, this selection leaves out events with a large numberof candidate hosts, that are not expected to make significantcontributions to the inference. When testing our pipeline, wefound that less informative, low SNR events not only do notadd much information to the inference, but tend to producea bias in the estimate of cosmological parameters, especiallythose that are weakly constrained. More details about thisbias are given in App. A2. Its origin is still under investigationand will be a subject of future work. Although some level ofbias can also be seen in our main results (see the posterior for w a in Fig. 5), it does not affect the estimate of the constrainedparameters. Having selected a sample of “useful” EMRIs, the next step isto simulate realistic distributions for candidate host galaxies.To this end, we closely follow the procedure outlined in DelPozzo et al. (2018), with some improvements, as we now sum-marize.We use the simulated observed sky of Henriques et al.(2012) constructed from the Millennium run (Springel et al.2005) which, critically for this work, reproduces the ob-served clustering properties of galaxies. The catalog, thatwe use up to z = 1 , lists a number of galaxy properties,including mass, observed redshift z obs and cosmological red-shift z cos , which differ due to the extra redshift contribu-tion imprinted on the former by galaxy peculiar velocities, ∆ z v p ≡ z obs − z cos = (1 + z cos ) v p /c ( v p (cid:28) c ) , v p being thepeculiar velocity.The simulation assumes a Λ CDM cosmology with h =0 . , Ω m = 0 . , Ω Λ = 0 . , which will thus correspond tothe fiducial values for our analysis as well. Although thoseparameters are outdated (Ade et al. 2016), this is irrelevantfor the illustrative purpose of our analysis.The LISA 3D error volume is approximated by a multivari-ate Gaussian distribution: p ( D i | Θ ) = 1 (cid:112) (2 π ) | Σ | exp (cid:26) − ˜Θ Ti Σ − ˜Θ i (cid:27) , (4.1)where ˜Θ i ≡ Θ − ˆ Θ i ( D i ) , and the 3D parameter vector in-cludes the source luminosity distance, the cosine of its decli-nation and its right ascension, i.e. Θ = ( d L , cos θ gw , φ gw ) .The vector ˆΘ i ( D i ) = ( ˆ d L , cos ˆ θ gw , ˆ φ gw ) i defines the bestmeasured values of the parameters (i.e. the center of theGaussian), which depend on the observed data. We assumethat these are the only quantities measured from the data MNRAS , 1–20 (2021)
D. Laghi et al.
Figure 2.
Each row of three panels represents the properties of the galaxy host candidates for a selected event extracted from modelM1_2. The left panels show a scatter plot of candidate hosts in the θ − φ plane; the shades of blue, from lighter to darker, are in orderof increasing probability of hosting the EMRI, according to the LISA sky localisation in the celestial sphere. The red dot is the truehost, which is generally not in the center of the sky localisation area. The middle panels show the redshift probability distributions of thecandidates hosts (black histograms), i.e. the distribution of galaxies, each weighted by its own probability of being the host according toits sky location. The dark-green vertical line is the fiducial redshift measured by LISA and the two light-green lines bracket the 2 σ redshiftinterval implied by the LISA measurement error in the luminosity distance alone, assuming the true cosmology. The vertical yellow lines atthe edges bracket the redshift interval uncertainty when the cosmological parameters are allowed to vary within the assumed priors (here Λ CDM is assumed). The redshift of the true host is marked by the vertical red line. Finally, the right panels represent the deviation ofthe x quantities, θ (magenta), φ (blue), and d L (black), from the best measured value x t normalized to the respective LISA measurementerrors. The vertical dotted lines are the values of the true hosts. Note that the d L of candidate galaxies can deviate from the best measuredvalue by many σ , as measured by LISA, due to the prior uncertainties on the cosmological parameters. and drop the explicit dependence on D i in subsequent equa-tions. The measurement uncertainties and correlations aredescribed by the 3D correlation matrix Σ , which we ex-tract from the full fifteen-dimensional correlation matrix ofthe EMRI parameter estimation carried out in Babak et al.(2017), by marginalising over all other parameters. In prac-tice, these posterior widths could depend on both the trueparameters of the EMRI and on the particular realisationof the noise in the LISA data. For the well localised EM-RIs we select in our analysis, these dependencies are likely to be weak. Therefore, we assume here that the measure-ment uncertainties are fixed and known for each event. Thefirst task is to associate a true host to each EMRI, consis-tent with these errors. To pick the true host, we first leave (cos θ gw , φ gw ) free, we compute d L for all galaxies in the Mil-lennium sky and then evaluate the marginalised luminosity MNRAS000
Each row of three panels represents the properties of the galaxy host candidates for a selected event extracted from modelM1_2. The left panels show a scatter plot of candidate hosts in the θ − φ plane; the shades of blue, from lighter to darker, are in orderof increasing probability of hosting the EMRI, according to the LISA sky localisation in the celestial sphere. The red dot is the truehost, which is generally not in the center of the sky localisation area. The middle panels show the redshift probability distributions of thecandidates hosts (black histograms), i.e. the distribution of galaxies, each weighted by its own probability of being the host according toits sky location. The dark-green vertical line is the fiducial redshift measured by LISA and the two light-green lines bracket the 2 σ redshiftinterval implied by the LISA measurement error in the luminosity distance alone, assuming the true cosmology. The vertical yellow lines atthe edges bracket the redshift interval uncertainty when the cosmological parameters are allowed to vary within the assumed priors (here Λ CDM is assumed). The redshift of the true host is marked by the vertical red line. Finally, the right panels represent the deviation ofthe x quantities, θ (magenta), φ (blue), and d L (black), from the best measured value x t normalized to the respective LISA measurementerrors. The vertical dotted lines are the values of the true hosts. Note that the d L of candidate galaxies can deviate from the best measuredvalue by many σ , as measured by LISA, due to the prior uncertainties on the cosmological parameters. and drop the explicit dependence on D i in subsequent equa-tions. The measurement uncertainties and correlations aredescribed by the 3D correlation matrix Σ , which we ex-tract from the full fifteen-dimensional correlation matrix ofthe EMRI parameter estimation carried out in Babak et al.(2017), by marginalising over all other parameters. In prac-tice, these posterior widths could depend on both the trueparameters of the EMRI and on the particular realisationof the noise in the LISA data. For the well localised EM-RIs we select in our analysis, these dependencies are likely to be weak. Therefore, we assume here that the measure-ment uncertainties are fixed and known for each event. Thefirst task is to associate a true host to each EMRI, consis-tent with these errors. To pick the true host, we first leave (cos θ gw , φ gw ) free, we compute d L for all galaxies in the Mil-lennium sky and then evaluate the marginalised luminosity MNRAS000 , 1–20 (2021) osmology with EMRIs distance likelihood given by: p ( d L | D i ) = exp (cid:40) − (cid:0) d L − ˆ d L ( D i ) (cid:1) σ d L + σ WL (cid:41)(cid:113) π (cid:0) σ d L + σ WL (cid:1) , (4.2)at the luminosity distance of each galaxy. A host galaxyis then randomly selected from the catalogue with prob-ability proportional to this marginalised likelihood. Leav-ing cos θ gw and φ gw free ensures that the pool of galax-ies from which we select the host is large enough that theselection process reflects the clustering properties of galax-ies at the measured distance of the event, while the lumi-nosity distance of the selected host, d L,H , is fully consis-tent with the LISA measurement error distribution by con-struction . The host is described by the parameter vec-tor Θ H = ( d L,H , cos θ gw,H , φ gw,H ) . We then fix d L = d L,H in Eq. (4.1) and draw a pair (cos θ (cid:48) gw , φ (cid:48) gw ) from the 2Dslice of the p ( Θ ) distribution. The LISA error volume isthen re-positioned so that cos ˆ θ gw = cos θ gw,H − cos θ (cid:48) gw and ˆ φ gw = φ gw,H − φ (cid:48) gw , which, together with ˆ d L , redefine thevector ˆΘ . This procedure ensures that the true host is drawnconsistently with the clustering properties of galaxies, andthat its position relative to the best LISA measured values ˆΘ is drawn according to the 3D volume error described byEq. (4.1). Now that we have identified a true host and haveplaced the LISA error volume appropriately with respect toit, we can truncate the galaxy catalogue by sub-selecting allgalaxies g j so that z obs , j is consistent with the LISA measure-ment, including the uncertainties in the cosmological parame-ters. In principle, the sum over galaxies in Eq. (2.10) includesall galaxies in the catalogue up to z = 1 . However, galaxieslying in the tail of the Gaussian distribution for all choices ofcosmological parameters in the prior are exponentially sup-pressed in the likelihood. Therefore, it is more efficient toremove these from the catalogue. We retain those galaxiesfor which z − < z obs , j < z + , where z − and z + are implicitlygiven by d L ± σ d L = (1 + z ± ) d (cid:90) z ± d z (cid:48) F ± ( z (cid:48) ) . (4.3)Here d = c/H is the Hubble distance, and F ( z ) = hE ( z ) ,with E ( z ) given in Eq. (2.6). F ± ( z ) are the realizations of F ( z ) that minimize and maximize the d L − z conversionwithin the assumed prior on the cosmological parameters.In practice Eq. (4.3) extends the ∆ z due to the σ errorin the LISA measurement of d L as much as allowed by theprior range on the cosmology. Finally, we need to keep inmind that due to peculiar velocities, the true cosmologi-cal redshift of each galaxy z cos (cid:54) = z obs , the difference be-tween the two being ∆ z v p . We thus consider all galaxies with z obs ∈ [ z − − ∆ z − v p , z + + ∆ z + v p ] ≡ ∆ z tot , where ∆ z − v p , ∆ z + v p are simply ∆ z v p calculated at z + , z − , taking a characteristic This is necessary because the EMRI catalogues used in this workwere generated independently from the Millennium sky and itsclustering properties. By allowing hosts to be selected within thewhole / π [( d L + 2 σ d L ) − ( d L − σ d L ) ] , we ensure that at leastthe hosts probability reflects the clustering of the Millennium skywithin that luminosity distance range. value of v p = 500 km s − , consistent with the standard de-viation of the galaxy radial peculiar velocity distribution inthe Millennium run (Springel et al. 2005).We compute the weights w j , appearing in Eq. (2.10), foreach galaxy g j within ∆ z tot by marginalization of Eq. (4.1)over d L (assuming uniform priors): p ( D i | cos θ j φ j ) = (cid:90) d d L p ( D i | Θ ) ≡ w j , (4.4)evaluated at (cos θ j , φ j ) (where w j must not be confused withthe DE equation of state parameters w and w a ). We doso for all galaxies falling into a sky region ∆Ω tot covering σ in the determination of cos θ gw and φ gw . In practice weassign probabilities only based on the marginalized 2D skylocation error, discarding further information included in thecorrelation of cos θ j and φ j with d L , which is a conservativeassumption.In summary, to assess the power of EMRI cosmology withLISA, we take three EMRI population models (M1, M5, andM6) from Babak et al. (2017), and we consider the two cos-mology scenarios presented in Sec. 2: Λ CDM and DE. We se-lect only EMRIs with SNR >
100 and associate to each eventcandidate hosts according to the procedure outlined above.For each model, the procedure of drawing the true host isrepeated three times for each EMRI, so that we have threeindependent realizations of the galaxy host candidates, mak-ing our analysis robust with respect to statistical fluctuationsin the galaxy distribution. We identify the independent real-izations with a number following the name of the model, sothat for model M1 we refer to M1_1, M1_2, and M1_3, thesame being for M2 and M3. Finally, for each model, we in-vestigate the impact of the mission lifetime by consideringobservation of EMRIs carried out in 4 years and 10 years ofLISA operations.An example of the outcome of the procedure outlined aboveis shown in Fig. 2, where we depict the properties of the se-lected galaxies within the error box ∆ z tot × ∆Ω tot for threeEMRI events selected from model M1_2 under the assump-tion of a Λ CDM cosmology. The cos θ − φ correlation is vis-ible in the left column panels, together with the clusteringpattern of the galaxies in the sky. Note that the true host,which differently from MacLeod & Hogan (2008) we do notremove from the group of possible hosts and is marked bythe red dot in the left panels and the red vertical line in themiddle panels, is generally offset from the center of the errorbox by construction, which mimics a realistic situation. Thevarious elements of the ∆ z tot construction procedure can beappreciated in the central panels. The dark-green line marksthe redshift corresponding to the best LISA distance mea-surement ˆ d L . The light-green lines bracket the ∆ z intervalassociated to the σ d L LISA error assuming the true (i.e.,theMillennium run) cosmology. Finally, the yellow vertical linesbracket the interval ∆ z tot , allowed by the prior range on thecosmological parameters and accounting for the errors dueto peculiar velocities (which are generally subdominant, un-less z (cid:28) . ). The black histograms show the probabilitydistributions of observed redshifts z obs of the hosts in theerror volume, where the probability of each individual host, w j , is given by Eq. (4.4). As it should be, z obs of the truehost generally falls within the light-green lines, although thisis not always strictly the case due to peculiar velocities. Fi-nally, the right panels show sanity checks of the offset distri- MNRAS , 1–20 (2021) D. Laghi et al. butions of d L , cos θ , and φ of the candidate hosts from thebest-measured LISA values normalized to the respective mea-surement errors. The distributions for d L (black histograms)extend to several σ due to the extra contribution allowed byvarying the cosmological parameters within the prior range.The top panels in Fig. 2 show a “very good EMRI”, forwhich the putative hosts cluster in the ∆ z range allowed bythe true cosmology. The middle panels show a “non infor-mative EMRI”, displaying an essentially flat (although verynoisy) probability distribution of the hosts across the wholeredshift range allowed by the cosmological prior. Finally, thebottom panels show a “bad EMRI”. In this latter case, thetrue host was picked in a small group of galaxies at z ≈ . ,and by chance, in the same sky region, there is a much largergroup of galaxies at z ≈ . , which skews the probabilitydistribution of the hosts towards a cosmology different fromthe true, although still within the allowed prior. Events ofthis type are expected to be subdominant, and by stackingthe posteriors of several events, the true cosmological param-eters naturally emerge from the analysis, as we now show. In this section we report the results of our investigation fol-lowing the procedure described in the previous sections. Wewill first show the results for the Λ CDM scenario (Sec. 5.1)and then for the DE scenario (Sec. 5.2). For both scenarioswe will consider 4 and 10 years of LISA observational time.Now we describe how we obtain the posteriors reported inthis section. Concerning the 4-year results, in order to accu-mulate enough statistical evidence to produce reliable resultsfor a 4-year LISA mission, for each EMRI model we gener-ated three 4-year catalogs from each of the three realisationsof the 10-year catalog described in Sec. 3.2 and 4. The num-ber of events in each 4-year catalog was obtained as follows:i) we randomly picked an integer n by drawing from aPoisson distribution with mean equal to the 4/10 of the totalnumber of events (see Table 1); ii) we then selected randomly n events from the original catalog; iii) finally we appliedthe SNR >
100 cut to obtain the final set of EMRIs for theanalysis. In this way we have a total of nine 4-year reali-sations per EMRI scenario (for example, in case of M1 weproduced three realisations for M1_1, three for M1_2, andthree for M1_3), each of which has been analysed separately.The posteriors for the cosmological parameters are then av-eraged (Del Pozzo et al. 2018) over the different realisations.This procedure is applied for both the Λ CDM and the DEcosmological inference models. Our results with 10 years ofLISA mission lifetime for each EMRI scenario are simply ob-tained averaging the analysis of the three different catalogrealisations (again, for M1 it will be the average of M1_1,M1_2, and M1_3). Such catalogs contain the exact numberof EMRI events as given in the fifth and sixth columns ofTable 1, depending on the assumed cosmological model.For illustrative purposes in Fig. 3 we show a distance-redshift diagram as obtained from our analysis in the Λ CDMscenario for the M1_2 realisation only of our fiducial EMRImodel, assuming 10 years of LISA mission lifetime. The cred-ible regions pictured in Fig. 3 for each EMRI event are theposterior distributions resulting from combining the LISAdistance measurement with the combined information col- lected from all galaxies within the sky localisation region as-sociated to that particular event. For each EMRI, to comparethe prior range for z gw with the corresponding range of galaxyhost redshifts, in Fig. 3 we show the redshift z j of each galaxy j within the sky localisation region of each EMRI event, dis-played as w j -weighted dots at a fixed luminosity distance set,for convenience, equal to ˆ d L . Note also how measurementsat low redshift are more accurate due to the lower numberof galaxies contained within the reconstructed volume, whileevents at higher redshift possess on average less constrainingpower due to the large number of galaxies within their skyerror boxes. Λ CDM
Results for the Λ CDM scenario are reported in Fig. 4 for ourthree selected EMRI models: M5 (pessimistic), M1 (fiducial),and M6 (optimistic). We show corner-plots with 2D posteri-ors in the h - Ω m parameter space, together with marginal-ized 1D posteriors on both h and Ω m . Results with both 4years (upper-row) and 10 years (lower-row) of LISA missionoperation are shown. The constraints reported in the figureindicate 90% confidence levels around the measured medianvalue. In what follows we will comment in detail on these twoscenarios. As described above, the joint posteriors reported in Fig. 4are averaged over the nine independent realizations of eachEMRI scenario.As reported in Fig. 4, the constraints we obtain on h with 4 years of LISA observations range from a maximumof 0.026 (3.6%), within the M5 scenario, to a minimum of0.012 (1.6%), within the M6 scenario. The constraints on Ω m go instead from 0.150 (60.0%) in the M5 scenario, to 0.081(32.0%) in the M6 scenario. In general, the predominance oflow-redshift events in the distance-redshift diagram is suchthat they are more efficient at constraining H rather than Ω m . Note that the posterior on Ω m in the M5 (pessimistic)case is actually dominated by our prior, implying that in thisscenario the LISA data do not yield statistically relevant in-formation on Ω m . Also note that, as expected, our resultsobtained with M1 (our fiducial scenario) are superseded bythose obtained with M6 (optimistic), as in this case the totalnumber of useful EMRIs is more than doubled (cf. Table 1). As shown by the lower-row plots in Fig. 4, the constraintson h now range from a worst-case result of 0.019 (2.6%) toa best-case result of 0.008 (1.1%), while the constraints on Ω m range from 0.107 (42.8%) to 0.058 (23.2%). Differentlyfrom the 4-year case, the posterior on Ω m in the M5 scenario(pessimistic) is not dominated by the prior, meaning that inthis scenario 10 years of mission duration could start yieldingmeaningful information on Ω m .The results obtained with M5, our pessimistic scenario,might seem too optimistic if one considers the fact that only6 EMRI events are counted in the 10-year catalog that weused (see Table 1). However, we have selected the loudest MNRAS , 1–20 (2021) osmology with EMRIs z d L / G p c c r e d i b l e l e v e l Figure 3. d L − z regression line for M1_2, one of the three realisations of our fiducial scenario M1 in the 10-year mission case: median(solid black) and 68% and 90% credible regions in light seagreen and light gray, respectively. The dashed line corresponds to the MillenniumSimulation fiducial values of h and Ω m . The coloured regions show the posterior distribution for z gw and d L for each EMRI event. Thehorizontal dots show the redshift of each candidate galaxy host for that particular EMRI. For illustrative purpose, we assigned to eachgalaxy a luminosity distance equal to ˆ d L . The dots are also colour-coded from light to dark for increasing values of the weights w j . events, which should be the most informative from a param-eter estimation point of view. An analogous considerationcould be made for the 4-year results, where the number ofEMRI events in the employed M5 catalogs are on averageonly n = 2 . .We have moreover checked that the 10-year M5 catalogcontains a low-redshift event which in the three versions ofthe 10-year catalog happens to have only 1, 3, and 8 galax-ies within its 3D localisation error box. In the catalog wherethis EMRI has a single host, the inference of its redshift isevidently optimal. In the catalog where it has three galax-ies, two of these galaxies are spot-on the true redshift value,while the third one is less relevant, because it lies at the mar-gins of the sky localisation region and thus, according to howour analysis is implemented (cf. Sec. 2), it is weighted witha lower probability with respect to the others. This eventalone therefore provides a highly accurate measurement ofthe Hubble constant at low redshift, which dominates theresults reported in the first column of Fig. 4. We note thatalthough the detection of such an advantageous event is aquestion of luck, the actual probability of finding a similarevent in the few EMRI events included in our catalogs is ac-tually non-negligible, in that our analysis automatically se-lects for well-localized events at high SNR. In general, anyEMRI event of this kind will substantially drive our results,while less well-localized, low-SNR events are not expected tocontribute (but see App. A2 for a test-case in which we only use the faintest events of our fiducial model). Results for the DE cosmological scenario are reported inFig. 5, again for all three EMRI models. We show corner-plotswith 2D posteriors in the w - w a parameter space, togetherwith marginalized 1D posteriors on both w and w a , for both4 years (upper-row) and 10 years (lower-row) of LISA mis-sion operations. Again constraints shown in Fig. 5 define 90%confidence levels around median values. In what follows wewill address in more details the two observational scenarios.As already noted in Sec. 4, the total number of useful eventsfor each EMRI model slightly varies according to the cosmo-logical model under consideration (cf. Table 1); thus the DEcatalogs have a slightly different total number of useful EM-RIs with respect to the ones used for Λ CDM.
From the top-row panels of Fig. 5 we can notice that con-straints on w range from a maximum of 0.123 (12.3%) toa minimum of 0.068 (6.8%) for the M5 (pessimistic) andM6 (optimistic) models, respectively. On the other hand, themeasurements on w a are largely uninformative as they yield90% confidence intervals of ∼ . (corresponding to ∼ of the prior) irrespective of the EMRI cosmological scenario.This implies that although LISA EMRIs might be able totell us something interesting on the current value of the DEequation of state ( w ), they are likely to not have the statisti-cal power to tell us anything about its current time evolution( w a ). MNRAS , 1–20 (2021) D. Laghi et al. Λ CDM, 4yr
Model M5 (pessimistic)
Model M1 (fiducial)
Model M6 (optimistic) h = 0.720 +0.0260.023 .
66 0 .
72 0 .
78 0 . h . . . . m . . . . mm = 0.279 +0.1440.150 h = 0.728 +0.0180.020 .
65 0 .
70 0 .
75 0 .
80 0 . h . . . . m . . . . mm = 0.254 +0.1120.095 h = 0.725 +0.0110.012 .
70 0 .
72 0 .
74 0 . h . . . . m . . . . mm = 0.291 +0.0810.070 Λ CDM, 10yr
Model M5 (pessimistic)
Model M1 (fiducial)
Model M6 (optimistic) h = 0.721 +0.0190.017 .
66 0 .
69 0 .
72 0 .
75 0 . h . . . . m . . . . mm = 0.292 +0.1070.103 h = 0.729 +0.0110.013 .
70 0 .
72 0 .
74 0 . h . . . . m . . . . mm = 0.251 +0.0710.058 h = 0.726 +0.0080.008 .
705 0 .
720 0 .
735 0 . h . . . . . m .
16 0 .
24 0 .
32 0 .
40 0 . mm = 0.290 +0.0580.052 Figure 4.
Corner plots of the posteriors for the parameters h and Ω m in the Λ CDM scenario. In each plot the lower-left panel reportsthe contours of the joint posterior, while the upper and right panels show the same posterior after marginalization over one parameter.In each panel, the cyan lines mark the fiducial values and the dashed lines indicate the median and 90% credible interval extracted fromthe marginalized posterior.
Results are reported in the bottom row of Fig. 5. We cansee that with 10 years of observations LISA EMRIs will yieldconstraints on w ranging from 0.083 (8.3%) to 0.059 (5.9%),slightly improving over the 4-year results. Results for w a areagain largely uninformative, as they still produce 90% con-fidence constraints around 0.6-0.7, which still constitute anon-negligible fraction of the prior. However, we notice theapparent railing of w a in case of M6, which becomes promi-nent in the 10-year case. We will comment on this in App. A1. Overall the results we obtained from our investigation showthat EMRIs detected by LISA will have an interesting cos-mological potential. Assuming the Λ CDM model, the Hub-ble constant will be probed at the percent level, while inthe evolving DE scenario the equation of state of DE will be tested with an accuracy better than 10%. Some constraint,even though not particularly strong ( ∼ ), can also beput on Ω m . The constraints on these three parameters, forall three EMRI models and both 4 and 10 years of LISAobservations, have been summarised in Fig. 6. On the otherhand, only a marginal gain will be achieved on w a (for DE),which usually is better constrained by measurements at highredshift, which are not so numerous in our catalogs. How do our results compare with current and future EMobservations? Current EM measurements of the Hubble con-stant have reached percent levels of accuracy with observa-tions collected from the CMB (Ade et al. 2016; Aghanim et al.2018) and from local-Universe distance indicators (Riess et al.2016, 2019), in particular type-Ia supernovae (SNIa). How-ever, as mentioned in the Introduction, these two techniquesyield different values of the Hubble constant which are now
MNRAS000
MNRAS000 , 1–20 (2021) osmology with EMRIs DE, 4yr
Model M5 (pessimistic)
Model M1 (fiducial)
Model M6 (optimistic) w = 1.058 +0.1100.123 . . . . w . . . . . w a . . . . . w a w a = 0.051 +0.6480.703 w = 0.987 +0.0750.074 .
25 1 .
00 0 .
75 0 . w . . . . . w a . . . . . w a w a = 0.029 +0.6550.676 w = 1.015 +0.0680.064 . . . . . w . . . . . w a . . . . . w a w a = 0.145 +0.5920.697 DE, 10yr
Model M5 (pessimistic)
Model M1 (fiducial)
Model M6 (optimistic) w = 1.053 +0.0830.079 .
20 1 .
05 0 . w . . . . . w a . . . . . w a w a = 0.057 +0.6520.721 w = 1.000 +0.0620.063 .
12 1 .
04 0 .
96 0 .
88 0 . w . . . . . w a . . . . . w a w a = 0.085 +0.6120.676 w = 1.023 +0.0590.043 .
12 1 .
04 0 .
96 0 . w . . . . . w a . . . . . w a w a = 0.396 +0.4320.640 Figure 5.
Corner plots of the posteriors for the parameters w and w a in the DE scenario. Same as Fig. 4. See the Appendix for adiscussion on the M6 model. in tension with each other at more than the 4 σ statisticallevel (Riess et al. 2019). Whether this tension is due to sys-tematics in one of the two measurements or to new physicsbeyond Λ CDM, it will be only decided by future observations,including the ones performed with GW standard sirens (Chenet al. 2018). The advantage of standard sirens in this respectis that their measurements present completely different sys-tematics with respect to both CMB and local EM measure-ments, and will thus deliver completely independent infor-mation on the Hubble constant, which will hopefully pointtowards its correct value. Similarly, the interesting aspect ofthe estimates that we obtained on H from our analysis isnot that we will reach a level of accuracy comparable withcurrent and future EM probes, but that EMRIs detected withLISA will provide yet another complementary measurementof H , useful to corroborate results from both EM observa-tions and ground-based GW cosmological measurements. Inthe likely scenario in which by the time LISA flies we willalready have solved the Hubble tension, the measurementsof H with LISA EMRIs will help consolidate our controlover all systematics, including the ones affecting GWs mea- surements with ground-based interferometers, and our con-fidence over the results obtained, especially in the case inwhich physics beyond Λ CDM is discovered.We can also compare our results for the dynamical DE cos-mological scenario with current measurements and forecastsfor the future. At present, distance measurements taken withSN Ia alone, i.e., not combined with CMB observations, con-strain w at the ∼ . ( ∼ Ω m ,but not for w a ). Future EM observations collected with theEuclid mission are expected to improve the constraints up to ∆ w (cid:39) . and ∆ w a (cid:39) . (Amendola et al. 2018). Theresults we obtained for w in Sec. 5 reach the same level of ac-curacy of present EM observations in the worst case scenario(M5), while they match the expected results from the Euclidmission in the best case scenario (M6). This implies that cos-mological measurements with LISA EMRIs will be of greatinterest to probe the DE equation of state, not only becausethey will reach constraints comparable with EM probes, butmost importantly because, as stressed above, they will pro-vide this level of accuracy with completely different system- MNRAS , 1–20 (2021) D. Laghi et al.
Summary plots h M M M M M M h N = N = N = N = N = N = Ω m M M M M M M m N = N = N = N = N = N = w M M M M M M w N = N = N = N = N = N = Figure 6.
Summary of results. We show 90% (black) and 68%(red) percentiles, together with the median (red dot) for both h , Ω m , and w for each EMRI model (M5, M1, M6) and LISA ob-servational scenario (4 and 10 years) considered. The blue dashedhorizontal line denotes the true cosmology. For each data point,we also report the average number N of EMRIs considered in theanalysis, see Table. 1. atics, increasing our confidence on such measurements (espe-cially if the cosmological constant value w = − appears tobe ruled out).The constraints obtained on the other two parameters con-sidered in our analysis, namely Ω m for Λ CDM and w a forDE, seem not to be comparable with the level of precision of current and future EM observations (Scolnic et al. 2018;Amendola et al. 2018). This means that although they willbe affected by completely different systematics, most likelythey will not be able to provide useful additional informationwith respect to EM observations. In other words, the situa-tion will be similar to the H constraints that we have todaywith GWs (Abbott et al. 2017c): even if they constitute acomplementary measurement of H , affected by completelydifferent systematics, the level of precision of current GW ob-servations does not offer a result comparable to EM measure-ments, implying that interesting cosmological insight beyondbroad consistency between different measurements cannot beobtained.Better constraints on cosmological parameters beyond H and w might be obtained by extending our analysis to EM-RIs with lower SNR, since ideally we would have additionalpoints to place on our d L − z plot. However, as noticed inSec. 3 and as discussed in App. A2, the statistical identi-fication method presents some difficulties when applied tohigh-redshift EMRI events having a large number of poten-tial hosts. In our present framework the modelling and infer-ence of galaxy evolution features as well as the inclusion oflow-SNR EMRIs, which could allow to better probe furthercosmological parameters, are not taken into account and willbe investigated in future studies. A more promising strategynevertheless may be to combine measurements from differentGW sources detected by LISA, as discussed hereafter. As we mentioned in Sec. 1, EMRIs are not the only LISAGW sources that can be used as standard sirens. Other in-vestigations considering SOBHBs and MBHBs detected byLISA have already produced cosmological estimates. How doour results compare with those?LISA SOBHBs will be mainly detected at low redshift( z (cid:39) . ) and thus will be only useful to constrain H . Recentanalyses (Kyutoku & Seto 2017; Del Pozzo et al. 2018) sug-gest that SOBHBs could be used to probe H at the few %level, similarly to the results we obtained with EMRIs. SOB-HBs can however reach this level of accuracy only if the detec-tion rate will fall in the optimistic range, which will cruciallydepend on the yet uncertain sensitivity that LISA will beable to achieve at high frequencies (Moore et al. 2019). Notealso that an “archival” search for SOBHBs will be possible:detecting the black hole merger with the third generation ofground-based GW detectors and then search those sources inthe (archival) LISA data with the reduced prior (Ewing et al.2021). On the other hand, EMRIs rely on the LISA sensitiv-ity around the mHz (middle-band) which is guaranteed to beat the level considered in our analysis (Amaro-Seoane et al.2017), if not better. Nevertheless as we have seen, EMRI ex-pected rates are affected by other uncertainties of astrophys-ical and theoretical nature. It is important to notice thatLISA will fail in deliver compelling low-redshift cosmologicalresults on H only if both the LISA high-frequency sensitiv-ity does not perform as expected and the true astrophysicalpopulation of EMRIs falls on the very pessimistic side, thusnot providing enough detections. In all other scenarios wemight expect LISA to deliver useful low-redshift constraintson H , from either SOBHBs, EMRIs, or both. MNRAS , 1–20 (2021) osmology with EMRIs LISA MBHB mergers, which are expected to produce ob-servable EM counterparts, will instead probe the expansionof the Universe at much higher redshift, providing in this waycomplementary constraints to the ones obtained at low red-shift with both SOBHBs and EMRIs. LISA MBHBs mightyield constraints on H at the few % level (Tamanini et al.2016; Tamanini 2017; Belgacem et al. 2019b), comparable toSOBHB and EMRI results. This implies that LISA will de-liver accurate measurements of H from data sets at differentredshift ranges, possibly providing further insights into thecurrent Hubble tension between local measurements (low- z )and CMB observations (high- z ). LISA MBHBs will more-over be useful to probe other cosmological parameters whichcannot be constrained at low redshift, for example Ω m for Λ CDM. The integration of low- and high-redshift standardsirens will thus allow LISA to test the expansion of the Uni-verse with GWs from early-to-late cosmological times, with-out the need to combine it with other probes. This makesLISA a unique cosmological observatory (Tamanini 2017),whose thorough implications are currently being investigatedand will be presented in a future study (Tamanini et al. 2021).Finally, we can compare our cosmological results with theones forecast for future Earth-based GW interferometers. Thecurrent LIGO-Virgo network of detectors, with further im-provements and the addition of more interferometers, is ex-pected to produce constraints on H at the few % level onlyif the rate of multi-messenger detections will lie on the opti-mistic side Chen et al. (2018). Future third-generation detec-tors, such as ET (Punturo et al. 2010; Maggiore et al. 2020),are instead expected to precisely probe the expansion of theUniverse up to z ∼ (Sathyaprakash et al. 2010; Taylor &Gair 2012; Del Pozzo et al. 2017; Belgacem et al. 2019a). Theywill provide complementary results with respect to LISA inapproximately the same time period (2030s), implying thatcross-checking GW cosmological measurements from spaceand from the ground will increase our confidence in the ob-tained results and will help us to get a handle on the expectedsystematics. Because of the reasons mentioned in Sec. 3, in our cosmolog-ical investigation we selected and analysed only three out ofthe twelve EMRI models considered in Babak et al. (2017).In first approximation, we can expect that our results canbe extrapolated to the other EMRI models according to thesquare-root of the number of detected events. In other words,the constraints we found on h and w should improve by afactor inversely proportional to √ N (Schutz 1986), where N is the total number of EMRIs useful for cosmology, reportedin Sec. 3 for each model.Unfortunately, it is difficult, if not impossible, to estimatethis number for the other EMRI models of Babak et al. (2017)without performing the detailed analysis we carried out forthe models M1, M5, and M6. For this reason, in what followswe use the total number of detected EMRI events, as givenby Table I in Babak et al. (2017) (AKK column), as a roughproxy for this number. We must keep in mind however thatdifferent EMRI scenarios will of course yield a population ofEMRIs with different properties, which will certainly trans-late in a different number of cosmologically useful events afterone applies the threshold we defined in Sec. 3. The estimates that follow should hence only be taken as simple indicativenumbers useful to quickly explore other possible EMRI sce-narios, rather than statistically robust results.The majority of the EMRI scenarios considered in Babaket al. (2017) present a number of detections comparable toM1, our fiducial model. We can thus expect that for all thesescenarios we would achieve cosmological constraints compara-ble to the ones we obtained with M1. Model M8 yields a num-ber of detections similar to M5, and thus for M8 we expectroughly the same cosmological results as the more pessimisticones we found in our analysis. The most pessimistic modelof Babak et al. (2017), namely model M11, provides only onedetected event per year, way below the detections achievedwith M5, our pessimistic scenario. Needless to say, we wouldnot expect any useful cosmological results with M11.On the other hand, there are three models which predictmore detections than those obtained in our optimistic sce-nario M6. We can thus extrapolate results starting from thedetections of M6 by using our simple √ N estimate. ModelM3 presents roughly twice the detections obtained with M6,which translates into improved 4-year (10-year) constraints(90% C.I.) as ∆ h (cid:39) . (0.0062) and ∆ w (cid:39) . (0.045),corresponding to relative errors of 1.3% (0.84%) and 5.2%(4.5%), respectively. The two most optimistic scenarios arehowever models M7 and M12, which deliver an order of mag-nitude more EMRI detections with respect to M6. WithinM7 we estimate improved 4-year (10-year) constraints as ∆ h (cid:39) . (0.034) and ∆ w (cid:39) . (0.025), correspond-ing to 0.69% (0.46%) and 2.9% (2.5%) relative uncertainties,respectively. For M12 instead we find ∆ h (cid:39) . (0.0027)and ∆ w (cid:39) . (0.020), corresponding to 0.55% (0.37%)and 2.3% (2.0%).According to these rough estimates, we could reach sub-percent constraints on H and two-percent constraints on DE( w ) if the rate of EMRIs detected by LISA is at the mostoptimistic end of expectations. From a cosmological perspec-tive these would constitute outstanding results, which evenEM probes might not be able to achieve in the near future.Having said that, we must stress again that these are onlysimple estimates which rely on extremely optimistic EMRIscenarios and do not take into account all the complexities ofthe true detected population of EMRIs. The analysis performed in the present study can be im-proved in several ways. First of all, the numbers and prop-erties of the EMRI population detected by LISA could bebetter characterized by considering enhanced astrophysicalmodelling and by refining the LISA response. In the errorestimation we have used the long-wavelength approximationto the response which takes into account the amplitude andthe phase (Doppler) modulation of the signal due to LISA’smotion around the Sun. This modulation is what gives usthe sky position used in this paper. The full response also in-cludes the sky-dependent time delays due to GW propagationacross satellites in the constellation. Taking the full response,therefore, could improve the sky localisation for EMRIs with
M < few × M (cid:12) . In addition, one could decisively improvethe data analysis treatment by starting from better EMRIwaveforms, possibly including inputs from self-force calcula-tions, and by performing a full Bayesian parameter estima- MNRAS , 1–20 (2021) D. Laghi et al. tion with a better model for the instrumental noise. Althoughsuch improvements would require further theoretical develop-ments and costly numerical computations, they would permitto consider full non-Gaussian posteriors in both distance andsky localisation, which will be important in making our anal-ysis more realistic, especially for low-SNR events.Another important addition that might be taken into con-sideration is modelling galaxy time evolution in order to char-acterize departures from homogeneity and uniformity in theemployed galaxy catalog (see Sec. 2). This would be a morerealistic representation of the Universe, but the price to payis the need to perform simultaneous inference of both the cos-mological parameters and the properties of the galaxy evo-lution, which on the one hand is computationally costly andmight degrade cosmological results, but on the other handwill provide complementary astrophysical insight.Another improvement would be to account for incomplete-ness of the galaxy catalogue. In the current analysis we haveignored this and assumed that out-of-catalogue hosts willbe close to hosts that are included in the catalogue. Thisbecomes increasingly problematic for more distant events,where a greater fraction of potential hosts will be missing.Incompleteness can be accounted for in the analysis, for ex-ample by weighting galaxies in the catalogue by the numberof nearby galaxies that are missing, or by adding an appro-priate number of missing galaxies into the assumed redshiftdistribution (Abbott et al. 2020d). The impact of such cor-rections should be explored in the future.One can also think of using galaxy observational features,such as luminosity, mass, metallicity, and others, to weightgalaxies differently within each GW sky localisation region.Moreover, one could use empirical relations, for example be-tween the mass of the MBH located at the center of thegalaxy, which is inferred from EMRI parameter estimation,and the mass or luminosity of galaxies in order to identifymore likely host galaxy candidates. Although these meth-ods would reduce the number of possible host galaxies foreach EMRI event, they would also introduce a dependencyupon astrophysical modelling into the analysis, possibly in-troducing new systematics, or the need to marginalise overadditional astrophysical parameters.We conclude this section by noting a few aspects of ourinference model that could be improved. First of all, we de-cided to neglect the intrinsic EMRI population evolution inthe prior assignment over z gw . The possibility of some biasin the estimation of the cosmological parameters due to arapid evolution of the rate of EMRIs with redshift, whichmeans that the weighting within each box should not be uni-form in z , has been already pointed out in MacLeod & Hogan(2008), although they did not account for this effect in theiranalysis. Just like for the cosmological evolution of the galaxypopulation, inclusion of this additional feature might degradethe overall inference over Ω . However, in scenarios where thenumber of detected EMRIs is large, this term would im-pose additional constraints on the whole EMRI populationthrough the time dependence of the merger rate, thus po-tentially increasing the amount of information contained ineach posterior. Moreover, we would be able to relax the ar-bitrary SNR cutoff of , using more faintest and possiblyfarther events, thus exploring a larger co-moving volume ofthe Universe.We stress that the present investigation constitutes a first simple attempt at assessing the cosmological potential ofLISA EMRIs. Further studies, improving the analysis alongthe lines outlined above, will be needed in order to providemore reliable forecasts and to prepare for the developing ofthe pipelines needed to analyse the expected data from LISA. In this article we investigated the cosmological potential ofEMRIs detected by LISA. By statistically matching the skylocalisation region of the loudest EMRI events (as given bythe analysis of Babak et al. (2017)) with the position and red-shift of galaxies within a given catalog (in our case based onthe Millennium run), we extracted constraints on the param-eters characterizing the background evolution of two cosmo-logical models: Λ CDM and a dynamical DE scenario. Our re-sults show that interesting cosmological insight can be gainedfrom EMRIs as standard sirens. Over three different EMRImodels and two LISA mission durations, constraints on H and w (the equation of state of DE) can respectively reachthe ∼ ∼ ∼ ∼ H and w to ∼ ∼ H can be reached with asfew as 20 EMRIs at redshift z < . . However, that analy-sis assumed parameter estimation accuracies appropriate fora more optimistic LISA design. It was shown in Gair et al.(2017) that 7/1/8 EMRI events satisfying these optimisticparameter constraints were found in two years of observationfor models M1/M5/M6, respectively. Therefore, we would ex-pect to achieve comparable precision with ten years of obser-vation of models M1 and M6, as we find here. The fact thatwe do not do better, even with the slightly larger numberof well-localized events and the inclusion of additional lesswell-localized events in the analysis, is most likely due to thesimplifications employed in the analysis in MacLeod & Hogan(2008), such as the use of a linear cosmic expansion modelwith H as the parameter.Finally, as in MacLeod & Hogan (2008), we let each EMRIsource contribute on an equal footing to the likelihood, while,differently from MacLeod & Hogan (2008), we do not assumeequal probability over the whole sky localisation region, butweight galaxies by the marginalised likelihood, see Eqs. (2.8)and (4.4).As a general final remark, we stress that LISA, a spacemission dedicated to GW science, will reveal itself as a uniquecosmological probe, through which we will be able to mapthe expansion of the Universe using different GW sources asstandard sirens at different redshifts, including, as thoroughlydemonstrated by our study, EMRIs. Future more-in-depthinvestigations may deliver cosmological forecast analyses with MNRAS , 1–20 (2021) osmology with EMRIs LISA EMRIs extended at higher redshift, thus allowing us toexplore the high-redshift Universe with dark sirens.
ACKNOWLEDGEMENTS
D.L. and W.D.P. thank Stefano Rinaldi for discussions.N.T. acknowledges partial support from the COST ActionCA16104 “Gravitational waves, black holes and fundamen-tal physics” (GWverse), supported by COST (European Co-operation in Science and Technology). A.S. is supported bythe European Research Council (ERC) under the EuropeanUnion’s Horizon 2020 research and innovation program ERC-2018-COG under grant agreement No 818691 (B Massive).
APPENDIX A: DISCUSSION OF CAVEATS ANDSYSTEMATICS
Here we discuss two aspects of the results presented in thepaper: (i) the railing of the parameter w a observed in the 10-year M scenario; (ii) the consequences of imposing an cutoffin SNR on the Λ CDM and DE analyses.
A1 Railing w a The right panels of Fig. 5 show the results for w and w a inthe 4 and 10-year optimistic scenarios (M6). It can be notedthat while w is correctly measured (being the true valuealways well within the 90% credible intervals), w a , althoughuninformative, shows a railing against the upper prior bound.This is particularly evident in the 10-year analysis (bottomright panel of Fig. 5, where the median value is 40% off thetrue value and the correlation between w and w a tends topush for w < − . It is worth to focus on this specific scenarioand see what is the possible cause of this behaviour.We use the same set of 73 loudest EMRI events used inthat scenario, this time associating to each event only onegalaxy host, chosen to be the nearest in redshift to the EMRI.This analysis may be viewed as representative of the physi-cal scenario in which all our events have an EM counterpart,so it can be used also to investigate what the DE paradigmwould predict in that case. In our likelihood, Eq. (2.10), thischoice corresponds to assigning w j = 1 to the nearest hostand w j = 0 to the others. The scope of this test is to seewhether the railing seen for w a is due to the particular na-ture of the EMRI catalog or to our formulation of the prob-lem. The results, shown in Fig. A1, seem to indicate thatindeed this is imputable to the latter. The railing is now ab-sent, showing posteriors for w and w a that are fully consis-tent with the fiducial values. Thus, cross-matching the EM-RIs with their nearest-in- z hosts gives substantially unbiasedposteriors. This seems to point to some limitations on ourformulation concerning the way in which the galaxy hostsare treated.Looking at Fig. A1, it is also interesting to note that theCLs on the measure of w are comparable to those of thegeneral case – where we include all galaxy hosts – reportedin Fig. 5: this shows that even if one were to observe EMcounterparts to these EMRIs, the inference would not sub-stantially improve. w = 1.003 +0.0570.056 . . . . . w . . . . . w a . . . . . w a w a = 0.021 +0.5090.518 Figure A1.
M6 DE 10 years, obtained analysing the same set ofevents used to produce the analogous corner plot in Fig. 5, butnow keeping only the nearest-in- z host for each EMRI. As can beseen, the railing behaviour of w a has now disappeared. A2 Low-SNR events
Even though we limited our analysis to high-SNR events, itis interesting to estimate the contribution of low-SNR eventsto the inference of the cosmological parameters under theassumption of likelihood (2.10). As an illustrative example,let us consider the EMRI catalogs for our fiducial model M1 inthe 10-year mission case. We will first adopt our first selectionprocedure detailed in Sec. 3, which impose constraints on themeasurement precision of the sky position and the luminositydistance (still associating to each EMRI hosts drawn from ourgalaxy catalog up to z = 1 ). Then, to limit ourselves to low-SNR events, we impose an upper cutoff SNR <
40 and analysethe events assuming our two cosmological models. In case of Λ CDM, the number of events is 30, while for DE there are39 events, see Table A1.The posteriors for the cosmological parameters are shownin the upper-row panels of Fig. A2: using only the quietestEMRIs does not yield uninformative posteriors. In particular,the correlation between the two parameters is mainly lost, asexpected for low-SNR events, while we observe a preferencefor high values of the parameters.To exclude the possibility of sampling issues or that thelow-SNR events show some unaccounted-for systematic dif-ference from the rest of the EMRI population, we analysedthe same set of low-SNR events, but keeping only the nearesthost in redshift to the EMRI, in a manner analogous to whatwe did in Sec. A1 for the loudest events of M6 in the DEscenario. As already noted in that section, such an analysismay be interpreted as the scenario in which we observe anEM counterpart to all our low-SNR EMRIs. The resultingposteriors are shown in the bottom-row panels of Fig. A2.The correlations between the parameters is restored and therailing has now disappeared, resulting in fairly informative
MNRAS , 1–20 (2021) D. Laghi et al. detected(SNR > ∆ d L /d L < . SNR < & ∆Ω < deg Λ CDM DE Λ CDM DE
M1 ( fid ) 2941 180 202 30 39 Table A1.
Number N of EMRIs observed by LISA in 10 years ofoperation imposing an upper cut in SNR, as done for the resultspresented in Fig. A2. The upper cutoff SNR < is imposed inorder to study the effect of the faintest events on the inferenceproblem. posteriors for h and w , a mildly informative posterior for Ω m , and a substantially uniform posterior for w a . The poste-riors indicate that the low-SNR population of EMRIs is con-sistent with the general one and that the inclusion of multiplepotential hosts in the analysis is likely the culprit. We havethus another indication that our treatment of multiple galaxyhosts is conditioned on our main assumptions, which thusdeserve further investigation. Hence, one should be aware ofthe caveats to this likelihood in a general analysis where alsolow-SNR events are included. It is also true that the faintestevents are expected to be less relevant to the inference prob-lem, so they should not contribute in a relevant way. This isdue to the large number of hosts that are typically associ-ated to each EMRI. The results presented in this work arebased on a different selection criterion (see Sec. 5) which donot suffer from any relevant bias in the estimate of the mainparameters we are interested in, that is, h , Ω m , and w . Weplan to further investigate the inclusion of low-SNR events inthe cosmological inference elsewhere. REFERENCES
Abbott B. P., et al., 2016a, Phys. Rev. Lett., 116, 061102Abbott B. P., et al., 2016b, Phys. Rev. Lett., 116, 221101Abbott B. P., et al., 2017a, Classical and Quantum Gravity, 34,044001Abbott B. P., et al., 2017b, Phys. Rev. Lett., 119, 161101Abbott B. P., et al., 2017c, Nature, 551, 85Abbott B. P., et al., 2017d, Astrophys. J., 848, L13Abbott B. P., et al., 2019a, Phys. Rev., X9, 031040Abbott B. P., et al., 2019b, Phys. Rev. Lett., 123, 011102Abbott B. P., et al., 2019c, Astrophys. J., 882, L24Abbott B. P., et al., 2019d, Phys. Rev., D100, 104036Abbott B. P., et al., 2020a, preprint ( arXiv:2001.01761 )Abbott R., et al., 2020b, preprint ( arXiv:2010.14527 )Abbott R., et al., 2020c, preprint ( arXiv:2010.14529 )Abbott B. P., et al., 2020d, preprint ( arXiv:1908.06060 )Ade P. A. R., et al., 2016, Astron. Astrophys., 594, A13Aghanim N., et al., 2018, preprint ( arXiv:1807.06209 )Amaro-Seoane P., Gair J. R., Freitag M., Miller M. C., Mandel I.,Cutler C. J., Babak S., 2007, Classical and Quantum Gravity,24, R113Amaro-Seoane P., et al., 2017, preprint ( arXiv:1702.00786 )Amendola L., et al., 2018, Living Rev. Rel., 21, 2Armano M., et al., 2016, Phys. Rev. Lett., 116, 231101Armano M., et al., 2018, Phys. Rev. Lett., 120, 061101Babak S., Gair J. R., Cole R. H., 2015, Fund. Theor. Phys., 179,783Babak S., et al., 2017, Phys. Rev., D95, 103012Barack L., Cutler C., 2004, Phys. Rev., D69, 082005Barack L., et al., 2019, Class. Quant. Grav., 36, 143001Bartolo N., et al., 2016, JCAP, 1612, 026 Behnel S., Bradshaw R., Citro C., Dalcin L., Seljebotn D. S., SmithK., 2011, Computing in Science & Engineering, 13, 31Belczynski K., Kalogera V., Bulik T., 2002, ApJ, 572, 407Belgacem E., Dirian Y., Foffa S., Howell E. J., Maggiore M.,Regimbau T., 2019a, JCAP, 08, 015Belgacem E., et al., 2019b, JCAP, 1907, 024Breivik K., Kremer K., Bueno M., Larson S. L., Coughlin S.,Kalogera V., 2018, Astrophys. J., 854, L1Cai R.-G., Tamanini N., Yang T., 2017, JCAP, 1705, 031Caprini C., Figueroa D. G., 2018, Class. Quant. Grav., 35, 163001Caprini C., Tamanini N., 2016, JCAP, 1610, 006Caprini C., et al., 2016, JCAP, 1604, 001Chen H.-Y., Fishbach M., Holz D. E., 2018, Nature, 562, 545Cusin G., Tamanini N., 2020, preprint ( arXiv:2011.15109 )Dal Canton T., Mangiagli A., Noble S. C., Schnittman J., Ptak A.,Klein A., Sesana A., Camp J., 2019, Astrophys. J., 886, 146Dalal N., Holz D. E., Hughes S. A., Jain B., 2006, Phys. Rev. D,74, 063006Del Pozzo W., 2012, Phys. Rev., D86, 043011Del Pozzo W., Laghi D., 2020, cosmolisa, https://github.com/wdpozzo/cosmolisa
Del Pozzo W., Li T. G. F., Messenger C., 2017, Phys. Rev. D, 95,043502Del Pozzo W., Sesana A., Klein A., 2018, Mon. Not. Roy. Astron.Soc., 475, 3485Dotti M., Sesana A., Decarli R., 2012, Adv. Astron., 2012, 940568Ewing B., Sachdev S., Borhanian S., Sathyaprakash B. S., 2021,Phys. Rev. D, 103, 023025Ezquiaga J. M., Holz D. E., 2020, preprint ( arXiv:2006.02211 )Farr W. M., Fishbach M., Ye J., Holz D., 2019, Astrophys. J. Lett.,883, L42Feeney S. M., Peiris H. V., Williamson A. R., Nissanke S. M.,Mortlock D. J., Alsing J., Scolnic D., 2019, Phys. Rev. Lett.,122, 061105Fishbach M., et al., 2019, Astrophys. J., 871, L13Foreman-Mackey D., 2016, The Journal of Open Source Software,24Gair J. R., Babak S., Sesana A., Amaro-Seoane P., Barausse E.,Berry C. P. L., Berti E., Sopuerta C., 2017, J. Phys. Conf. Ser.,840, 012021Giacomazzo B., Baker J. G., Miller M. C., Reynolds C. S., vanMeter J. R., 2012, Astrophys. J., 752, L15Gray R., et al., 2020, Phys. Rev. D, 101, 122001Haiman Z., 2017, Phys. Rev., D96, 023004Henriques B. M. B., White S. D. M., Lemson G., Thomas P. A.,Guo Q., Marleau G.-D., Overzier R. A., 2012, MNRAS, 421,2904Hirata C. M., Holz D. E., Cutler C., 2010, Phys. Rev. D, 81, 124046Hogg D. W., 1999, arXiv e-prints, pp astro–ph/9905116Holz D. E., Hughes S. A., 2005, Astrophys. J., 629, 15Hunter J. D., 2007, Computing in Science & Engineering, 9, 90Jaynes E. T., 2003, Probability Theory: The Logic of Science. Cam-bridge University Press, doi:10.1017/CBO9780511790423Kaplan D. L., O’Shaughnessy R., Sesana A., Volonteri M., 2011,Astrophys. J., 734, L37Klein A., et al., 2016, Phys. Rev., D93, 024003Kocsis B., Haiman Z., Menou K., 2008, Astrophys. J., 684, 870Korol V., Rossi E. M., Groot P. J., Nelemans G., Toonen S., BrownA. G. A., 2017, Mon. Not. Roy. Astron. Soc., 470, 1894Korol V., Koop O., Rossi E. M., 2018, Astrophys. J., 866, L20Kyutoku K., Seto N., 2017, Phys. Rev., D95, 083525LIGO Scientific Collaboration 2018, LIGO Algorithm Library -LALSuite, free software (GPL), doi:10.7935/GT1W-FZ16Lau M. Y. M., Mandel I., Vigna-Gómez A., Neijssel C. J., Steven-son S., Sesana A., 2019, preprint ( arXiv:1910.12422 )Linder E. V., 2003, Phys. Rev. Lett., 90, 091301Luo Z., Guo Z., Jin G., Wu Y., Hu W., 2020, Results in Physics,16, 102918MNRAS , 1–20 (2021) osmology with EMRIs Λ CDM, 10yr (SNR <
40) DE, 10yr (SNR < Model M1 Model M1 h = 0.790 +0.0440.053 . . . . . h . . . . m . . . . mm = 0.413 +0.0630.109 w = 0.568 +0.1720.267 . . . . . w . . . . . w a . . . . . w a w a = 0.361 +0.4410.688 Λ CDM, 10yr (SNR <
40) DE, 10yr (SNR < Model M1 (nearest host) Model M1 (nearest host) h = 0.727 +0.0370.033 . . . . h . . . . m . . . . mm = 0.264 +0.1350.126 w = 1.002 +0.0930.089 . . . . w . . . . . w a . . . . . w a w a = 0.002 +0.6560.660 Figure A2.
Corner plots of the posteriors on the parameters h , Ω m , w , and w a in the Λ CDM and DE scenario, respectively, for ourfiducial model M1 using the faintest events, selected imposing an upper threshold SNR < .MacLeod C. L., Hogan C. J., 2008, Phys. Rev., D77, 043512Maggiore M., et al., 2020, JCAP, 03, 050Mandel I., Farr W. M., Gair J. R., 2019, Monthly Notices of theRoyal Astronomical Society, 486, 1086–1093Mei J., et al., 2020, Progress of Theoretical and ExperimentalPhysicsMerritt D., 2015, ApJ, 814, 57Moesta P., Alic D., Rezzolla L., Zanotti O., Palenzuela C., 2012,Astrophys. J., 749, L32Moore C. J., Gerosa D., Klein A., 2019, Mon. Not. Roy. Astron.Soc., 488, L94Mukherjee S., Wandelt B. D., 2018, preprint ( arXiv:1808.06615 )Mukherjee S., Wandelt B. D., Nissanke S. M., Silvestri A., 2020a, preprint ( arXiv:2007.02943 )Mukherjee S., Wandelt B. D., Silk J., 2020b, preprint( arXiv:2012.15316 ), doi:10.1093/mnras/stab001Mukherjee S., Wandelt B. D., Silk J., 2020c, Phys. Rev. D, 101,103509Mukherjee S., Wandelt B. D., Silk J., 2020d, Mon. Not. Roy. As-tron. Soc., 494, 1956Mörtsell E., Dhawan S., 2018, JCAP, 1809, 025Nissanke S., Holz D. E., Hughes S. A., Dalal N., Sievers J. L., 2010,Astrophys. J., 725, 496Nissanke S., Holz D. E., Dalal N., Hughes S. A., Sievers J. L.,Hirata C. M., 2013, preprint ( arXiv:1307.2638 )O’Shaughnessy R., Kaplan D. L., Sesana A., Kamble A., 2011,MNRAS , 1–20 (2021) D. Laghi et al.
Astrophys. J., 743, 136Oguri M., 2016, Phys. Rev. D, 93, 083511Palenzuela C., Lehner L., Liebling S. L., 2010, Science, 329, 927Peebles P. J. E., 1993, Principles of physical cosmology. PrincetonUniversity PressPetiteau A., Babak S., Sesana A., 2011, Astrophys. J., 732, 82Poisson E., Pound A., Vega I., 2011, Living Rev. Rel., 14, 7Punturo M., et al., 2010, Class. Quant. Grav., 27, 194002Reitze D., et al., 2019a, preprint ( arXiv:1903.04615 )Reitze D., et al., 2019b, preprint ( arXiv:1907.04833 )Riess A. G., et al., 2016, Astrophys. J., 826, 56Riess A. G., Casertano S., Yuan W., Macri L. M., Scolnic D., 2019,Astrophys. J., 876, 85Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Phys. Rev. D,93, 084029Sathyaprakash B. S., Schutz B. F., Van Den Broeck C., 2010, Class.Quant. Grav., 27, 215006Schutz B. F., 1986, Nature, 323, 310Scolnic D. M., et al., 2018, Astrophys. J., 859, 101Sesana A., 2016, Phys. Rev. Lett., 116, 231102Soares-Santos M., et al., 2019, Astrophys. J., 876, L7Speri L., Tamanini N., Caldwell R. R., Gair J. R., Wang B., 2020,preprint ( arXiv:2010.09049 )Springel V., et al., 2005, Nature, 435, 629Tamanini N., 2017, J. Phys. Conf. Ser., 840, 012029Tamanini N., Caprini C., Barausse E., Sesana A., Klein A., Pe-titeau A., 2016, JCAP, 1604, 002Tamanini N., et al., 2021, in preparation.
Taylor S. R., Gair J. R., 2012, Phys. Rev. D, 86, 023502Veitch J., Del Pozzo W., Cody Pitkin M., 2020, cp-nest, doi:10.5281/zenodo.835874, https://doi.org/10.5281/zenodo.835874
Virtanen P., et al., 2020, Nature Methods, 17, 261Weinberg S., 1972, Gravitation and Cosmology: Principles and Ap-plications of the General Theory of Relativity. John Wiley andSons, New YorkWeinberg S., 2008, Cosmology. Cosmology, OUP Oxfordvan der Walt, Stéfan Colbert S. C., Varoquaux G., 2011, Comput-ing in Science & Engineering, 13, 22This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000