Trans-Dimensional Bayesian Inference for Gravitational Lens Substructures

Mon. Not. R. Astron. Soc. 000 , 000–000 (0000) Printed 9 Octob er 2018 (MN L A T E X style file v2.2) T rans-Dimensional Ba y esian Inference for Gra vitational Lens Substructures Brendon J. Brew er 1 ? , Da vid Huijser 1 , Gerain t F. Lewis 2 1 Dep artment of Statistics, The University of Auckland, Private Bag 92019, A uckland 1142, New Ze aland 2 Sydney Institute for Astr onomy, Scho ol of Physics, A28, The University of Sydney, NSW 2006, Austr alia ABSTRA CT W e in tro duce a Ba yesian solution to the problem of inferring the densit y profile of strong gra vitational lenses when the lens galaxy ma y contain m ultiple dark or faint substructures. The source and lens mo dels are based on a sup erp osition of an un- kno wn num b er of non-negativ e basis functions (or “blobs”) whose form was c hosen with speed as a primary criterion. The prior distribution for the blobs’ properties is sp ecified hierarchically , so the mass function of substructures is a natural output of the metho d. W e use reversible jump Marko v Chain Mon te Carlo (MCMC) within Dif- fusiv e Nested Sampling (DNS) to sample the posterior distribution and ev aluate the marginal likelihoo d of the mo del, including the summation ov er the unknown num b er of blobs in the source and the lens. W e demonstrate the method on t wo sim ulated data sets: one with a single substructure, and one with ten. W e also apply the metho d to the g-band image of the “Cosmic Horsesho e” system, and find evidence for more than zero substructures. How ever, these ha ve large spatial extent and probably only p oin t to misspecifications in the mo del (such as the shap e of the smo oth lens comp onen t or the p oin t spread function), whic h are difficult to guard against in full generality . Key w ords: gra vitational lensing: strong — methods: data analysis — metho ds: statistical 1 INTR ODUCTION Galaxy-galaxy gravitational lensing is a pow erful astrophys- ical to ol for studying the distribution of matter, including dark matter, in galaxies (T reu 2010). One promising appli- cation of lensing is to study and measure the prop erties of dark matter substructures in the lens galaxy (Ko opmans 2005). In recent y ears, this promise has been realised in sev eral lens systems whic h are thought to con tain at least one dark substructure (V egetti et al. 2012, 2010; V egetti, Czosk e, & Ko opmans 2010). In these studies, the lens w as mo delled as a sup erposition of a smo oth comp onen t (such as an elliptical pow er-law matter distribution) plus a pixel- lized “non-parametric” correction term. The estimated spa- tial structure of the correction term pro vides clues ab out the lo cations of p ossible substructures. In a second mo delling step, the lens is assumed to b e a smo oth ov erall component plus a compact “blob” near any lo cations suggested by the correction term. This kind of inference has also been done with point-lik e multiple images of quasars (F adely & Keeton 2012). ? T o whom correspondence should be addressed. Email: bj.brewer@auckland.ac.nz In this pap er, we in tro duce an approach for solving this same problem (inferring the n umber and prop erties of dark substructures giv en an image) in a more direct fashion, by fitting a model with an unknown num b er of substructures directly to the image data. How ever, with appropriate mo d- ifications, it may also prov e useful in other galaxy-galaxy lens modelling areas. One example is time-delay cosmog- raph y , whic h requires realistically flexible lens models, and reliable quantification of the uncertain ties (Suyu et al. 2013, 2014; Grillo et al. 2015). Lens modelling has been a topic of muc h interest ov er the last tw o decades. Most approaches are based on Ba yesian inference (Sivia & Skilling 2006; O’Hagan and F orster 2004), maxim um lik eliho od (Millar 2011) or v ariations thereof. These approaches may b e categorized according to the fol- lo wing criteria: (i) Whether to use a simply parameterized (e.g. a S´ ersic profile), mo derately flexible (e.g. a mixture model, Brew er et al. 2011) or free-form (e.g. pixellated) mo del for the surface brigh tness profile of the source; (ii) Whether to use a simply parameterized (e.g. Singular Isothermal Ellipsoid [SIE]) or flexible model (e.g. pixellated) for the mass profile of the lens; and (iii) Ho w to compute the results (e.g. optimization meth- c  0000 RAS 2 Br ewer, Huijser and L ewis o ds, Mark ov Chain Monte Carlo [MCMC] with source pa- rameters marginalized out analytically , or MCMC with the source parameters included). Recen t sophisicated approaches inv estigating different prior assumptions and computational approaches include Coles, Read, & Saha (2014), T agore & Jackson (2015), and Birrer, Amara, & Refregier (2015). Eac h approach inv olves tradeoffs b et ween conv enience and realism. Simply parameterized models, such as S´ ersic surface brigh tness profiles for the source, and singular isothermal ellipsoid plus external shear (SIE+ γ ) mo dels for the lens (Kormann, Sc hneider, & Bartelmann 1994), are v ery con venien t. They only ha ve a few adjustable parameters, and they capture (to “first order”) relev ant prior information that we hav e ab out the source and lens profiles. Of course, they are clearly simplifications, and can pro duce misleading results if the actual profile is v ery different from an y mem b er of the assumed family . On the other hand, pixellated mo dels for the source (e.g. Suyu et al. 2006) or the lens (e.g. Coles, Read, & Saha 2014) can in principle represent “an y” source surface brightness profile or lens pro jected densit y profile. Ho wev er, the prior distribution ov er pixel v alues is often chosen to b e a multi- v ariate gaussian for mathematical reasons, so that the source can b e analytically marginalized out (W arren & Dy e 2003). Unfortunately , a m ultiv ariate gaussian prior o ver pixel inten- sities usually corresp onds to a p oor mo del of our prior b eliefs ab out the source. It assigns virtually zero prior probability to the hypothesis that the source actually lo oks lik e a galaxy , and very high prior probability to the h yp othesis that the source lo oks like noise (or blurred noise). These priors also assign nonzero probability to negative surface brightness or densit y v alues; in fact, the marginal prior probability that an y pixel is negative is typically 0.5. Brew er et al. (2011) argued that an ideal mo delling ap- proac h lies somewhere betw een simply-parameterised and pixellated mo dels, which is also a motiv ation behind the in- v estigation of shapelets by T agore & Jackson (2015). One w ay of achieving this is with mixture mo dels. The source and the lens can be built up from a mixture of a mo derate n umber (from a few to a few h undred) of simply param- eterised comp onen ts. F or the source, this allows us to in- corp orate prior knowledge ab out the local correlations (the surface brigh tness at any particular point is lik ely to be sim- ilar to that at a nearby p oin t) and the fact that most of the sky is dark (Brew er & Lewis 2006). The Brew er et al. (2011) mo del also allow ed for multi-band data, and allow ed for our exp ectation that the image may or ma y not b e similar in differen t bands. The present pap er is similar in approach to Brew er et al. (2011). The main differences are: i) w e use a sim- pler and faster set of basis functions; ii) we apply it to the lens as well as the source, allowing for the p ossibil- it y of substructure; and iii) the implementation is based on a C++ template library developed by Brewer (2014) whic h allows for hierarchical priors and trans-dimensional MCMC to be implemented in a relatively straightforw ard manner. How ever, in the present paper w e do not account for m ulti-band data; this will be reserved for a future con- tribution. The source co de for this pro ject is a v ailable at http://www.github.com/eggplantbren/Lensing2 . If w e are in terested in detecting and measuring the prop erties of p ossible dark substructures in a lens galaxy , a sup erp osition of an unkno wn num b er of “blobs” is the most natural mo del. In addition, if we w ant to constrain the mass function of these blobs (e.g. V egetti & Koopmans 2009; V egetti et al. 2014), we will need a hierarchical model which sp ecifies the prior probabilit y distribution for the masses conditional on some hyperparameters. Inferring the mass function of the blobs then reduces to calculating the p oste- rior distribution for the hyperparameters. This is related to the general principle that we should (ideally) construct our inference metho ds so that the quantities we infer directly an- sw er our scientific questions (Alsing et al. 2015; Schneider et al. 2015; Pancoast, Brewer, & T reu 2011). Here, we are referring to the mass function of the substructures within a single lens galaxy . An additional level of hierarch y would b e needed to answer questions about a sample of lens galaxies. 2 THE MODEL W e now describe the details of our mo del and its parame- terization. The motiv ation for most of our mo delling c hoices is a compromise betw een computational efficiency , realism, and ease of implementation. None of the choices we hav e made are final in any sense; rather, this mo del should b e considered a proof of concept. W e encourage exploration of other c hoices. Since w e can compute marginal likelihoo ds, Ba yesian mo del comparison b et ween our choices here and an y proposed alternatives should b e straightforw ard, if a union of our hypothesis space and another can b e consid- ered a reasonable mo del of prior uncertain ty . 2.1 The Source The surface brightness profile of the source is assumed to b e comp osed of a sum of a finite num b er of “blobs”, or basis functions, in order to allo w some flexibilit y while incorpo- rating prior information ab out the non-negativity of surface brigh tness, and the spatial correlation exp ected in real sur- face brightness profiles. F or computational sp eed, we c ho ose the follo wing functional form for the surface brightness pro- file of a single blob centered at the origin: f ( x, y ) = ( 2 A πw 2  1 − r 2 w 2  , r 6 w 0 , otherwise (1) where r = p x 2 + y 2 , w is the width of the blob, and A is the total flux of the blob (i.e. the in tegral of the surface brigh tness ov er the en tire domain). These basis functions are in verted paraboloids, which are faster to ev aluate than gaus- sians, since they do not contain an exp onen tial function. In addition, the finite supp ort means that each blob will ev alu- ate to zero ov er a large fraction of the domain whic h confers an additional sp eed adv antage. More conv entional choices suc h as gaussians are p ossible, as are alternativ e compact basis functions, p erhaps inspired by the smoothed particle h ydro dynamics literature (Dehnen & Aly 2012). Alterna- tiv e choices, such as S ´ ersic profiles, may b e more realistic despite their computational cost. This is esp ecially likely in the case of an early-t yp e source galaxy (e.g. Auger et c  0000 RAS, MNRAS 000 , 000–000 3 al. 2011). Suc h modifications are p ossible without great ef- fort using the RJObject library , and once implemented, the marginal lik eliho od can be compared across differen t choices for the source mo del. If our mo del con tains N src suc h blobs, positioned at { ( x src i , y src i ) } with widths { w i } and total fluxes { A i } , the o verall surface brigh tness profile is: f ( x, y ) = N src X i =1 ( 2 A i πw 2 i  1 − r 2 i w 2 i  , r i 6 w i 0 , otherwise (2) where r i = p ( x − x src i ) 2 + ( y − y src i ) 2 . Under these assump- tions, the source can b e describ ed in its entiret y by the fol- lo wing parameters: n N src , { θ src i } N i =1 o (3) where θ src i denotes the parameters for source blob i : θ src i = { x i , y i , A i , w i } . (4) The dimensionality of the parameter space for the source de- p ends on the minimum and maximum v alues of N src , which w e set to 0 and 100 respectively (more detail ab out priors is given in Section 3). Therefore the source is described b y 1 – 401 parameters. 2.2 The Lens The surface mass densit y profile of the lens is mo delled as a sup erposition of a singular isothermal ellipsoid plus external shear (SIE+ γ ) and N circular lens “blobs”. The SIE+ γ is a simple and widely used lens model with analytically a v ail- able deflection angles, and is in tended to account for the bulk of the lensing effect due to the smooth spatial distribution of (visible and dark) matter in the lens galaxy . Unfortunately , along with simplicity comes a lack of realism, and generalis- ing the smo oth lens mo del is an imp ortant future step. The softened pow er la w elliptical mass densit y or SPEMD mo del is a popular generalization of the SIE where the deflection angles can still b e computed quic kly using the appro xima- tions giv en by Bark ana (1998). Ultimately , concentric mixtures of smo oth comp onen ts ma y b e the most useful for realistic inference of the den- sit y profile of the halo. A common approac h in the lensing comm unity , that is similar in spirit, is a sup erposition of an elliptical p o wer law, external shear, and a pixellated p oten- tial correction. How ever, the direct use of blobs is closer to the scientific question at hand when inv estigating p osten tial dark substructures. The SIE+ γ has nine free parameters: the (circularized) Einstein radius b , axis ratio q , central p osition ( x c , y c ), ori- en tation angle θ , the external shear γ , and the orientation angle of the external shear, θ γ . The N lens blobs are in tended to model p ossible dark or fain t substructures in the pro jected mass profile of the lens. A lensing blob with mass M and width v centered at the origin has the follo wing surface mass density profile: ρ ( x, y ) = ( 2 M πv 2  1 − r 2 v 2  , r 6 v 0 , otherwise (5) whic h is the same as the surface brightness profile of a source blob. The deflection angles for a single blob are: α x ( x, y ) =   2 − r 2 /v 2  M x πv 2 , r 6 v M x πr 2 , otherwise (6) α y ( x, y ) =   2 − r 2 /v 2  M y πv 2 , r 6 v M y πr 2 , otherwise (7) where r = p x 2 + y 2 . These do not dep end on an y compu- tationally exp ensiv e functions such as square ro ots or exp o- nen tials. F or N such blobs the deflection angles are summed o ver all blobs. Therefore the lens can b e entirely describ ed b y the following parameters:  θ SIE , N lens , n θ lens i o N i =1  (8) where θ lens i denotes the parameters for lens blob i : θ lens i = n x lens i , y lens i , M i , v i o , (9) and θ SIE describ es the parameters for the SIE+ γ comp o- nen t: θ SIE = n b, q , x SIE c , y SIE c , θ , γ , θ γ o . (10) Since the model is not of fixed dimension (the num- b er of source and lens “blobs” is unkno wn), we used a trans-dimensional MCMC sampler based on rev ersible jump MCMC (Green 1995). This framework has been used successfully in many astronomical inference problems (e.g. Jones, Kash yap, & v an Dyk 2014; Umst¨ atter et al. 2005). W e implemented the MCMC using the RJObject soft ware (Brewer 2014), a C++ library for implementing trans-dimensional MCMC with hierarchically-specified pri- ors when the mo del is a “mixture model” or similar, as is the case here. The RJObject library uses the Diffusiv e Nested Sampling algorithm (DNS; Brew er, P´ arta y , & Cs´ an yi 2011) for its sampling, but the trans-dimensionalit y is handled by the MCMC mov es. Therefore, the marginal li kelihoo d we ob- tain is one that inv olves a sum ov er the hypothesis space for N src and N lens . In other w ords, we do not need separate runs with different trial v alues of N src and N lens . Previous astro- nomical applications of RJObject include Huppenkothen et al. (2015) and Brew er & Dono v an (2015). Proposal mov es for the parameters (p ositions and masses of the source and lens blobs), and birth/death mov es to add or remov e blobs to the source or lens, are handled internally b y RJObject and fit in to the general framew ork describ ed in Brewer (2014). 3 PRIOR PR OBABILITY DISTRIBUTIONS The prior probability distributions for all h yp erparameters, parameters, and the data, are given in T able 1, and a fac- torization of the joint distribution is display ed in Figure 1. With these priors, we aim to express a large degree of prior uncertain ty (hence the lib eral use of uniform, log-uniform, and Cauc hy distributions). Ho wev er, we also sp ecify some priors hierarchically to allow (for example) blobs to clus- ter around a typical cen tral p osition, rather than implying a high probabilit y for substructure p ositions b eing spread uni- formly (in a frequency sense) o ver the sky . W e also applied hierarc hical priors to the substructure masses (and source fluxes) so that, while the typical order of magnitude of the c  0000 RAS, MNRAS 000 , 000–000 4 Br ewer, Huijser and L ewis masses is unknown, the masses themselv es are likely to be roughly the same order of magnitude. The h yp erparameters of these distributions can also be interpreted as straight- forw ard answers to questions about the substructure mass function. W e assign somewhat informative (but heavy-tailed) Cauc hy priors to the cen tral positions of the source and the lens. One migh t ob ject to the Cauch y priors for the central p ositions (and argue instead for gaussians) on the basis that they do not hav e rotational symmetry . How ever, the priors w e are assigning are prior only to the v alues of the data pixel in tensities only , and not other facts ab out the data, suc h as its dimensions in arc seconds, or its rectangular shap e. Given these, there is no reason to insist on rotational symmetry . The Cauch y priors are intended to enhance the plausibility (relativ e to what a uniform prior would imply) that the lens and source are somewhere near the centre of the image, but in a cautious wa y . The “circular” conditional prior for the blob positions has the following density: ρ ( x, y ) dx dy = 1 2 π W 1 r exp  − r W  dx dy (11) where r = p ( x − x c ) 2 + ( y − y c ) 2 . This results from an ex- p onen tial distribution (with exp ectation W ) for the radial co ordinate and a uniform distribution for the angular coor- dinate in plane polar co ordinates. This w as chosen (as op- p osed to a more “obvious” choice suc h as a normal distribu- tion) for computational reasons: the RJObject co de requires a function to transform from Uniform(0, 1) distributions to this distribution and bac k (making use of cumulativ e distri- bution functions and their inv erses), and the normal w ould ha ve required the use of sp ecial functions for this. 3.1 Conditional prior for the data The conditional prior for the data giv en the parameters is a pro duct of indep enden t gaussian distributions, one for each pixel. The mean of the gaussian is given by the “mo c k” noise-free image w e would expect based on the parameters. The standard deviation of the gaussian is a combination of three terms; the first (denoted s ij ) is a “noise map” which is loaded from a file, the second is an unkno wn constan t σ 0 whic h applies to the whole image, and the third is prop or- tional to the square ro ot of the mock image, with propor- tionalit y constan t σ 1 . This is usually called the “sampling distribution”, or sometimes just the “likelihoo d”. Ho wev er, sampling distri- bution is misleading since no physical frequency distribution exists whic h is b eing sampled from. The term lik eliho od is also usually used to refer to the scalar function of the param- eters obtained when the data are known. F or a discussion of the view that this ob ject is really a prior distribution, see Catic ha (2008, pp. 33–35). It is common to provide a “v ariance image” to lens mo d- elling softw are, which sp ecifies the standard deviations used in the lik eliho od function. How ever, telescopes do not pro- vide v ariance images. It also do es not make sense to sp ecify v ague priors prior to the image pixel v alues but p osterior to the v ariance image (which usually resembles the image itself, and therefore would b e highly informativ e). Allowing the standard deviation of the gaussian in the lik eliho od to dep end on the mock image brightness is a more principled w ay of mo delling our actual inferential situation. Nev erthe- less, we allo w for an input v ariance image as well, which can b e used for ad-ho c masking of troublesome regions (such as unmo delled flux from other sources). 4 COMPUT A TION Computing the p osterior distribution o ver the parameters of suc h a model requires that we can implemen t Marko v Chain Mon te Carlo (MCMC) ov er the space of possible sources and lenses. T o compute the p osterior distribution for the parameters (of which there are 24–464, depending on the v alues of N src and N lens ), we use Diffusive Nested Sampling (Brew er, P´ arta y , & Cs´ anyi 2011), a form of Nested Sampling (Skilling 2006) that uses the Metropolis algorithm to mov e around the parameter space. The proposal distributions for the blob parameters (b oth source and lens) are handled by the RJObject library (Brew er 2014). This includes birth and death proposals that increase or decrease either N src , or N lens , as well as prop os- als that mo ve the blobs (in their parameter space) while k eeping the n umber of blobs fixed. RJObject also facilitates prop osals that change the h yp erparameters (either for the lens or the source) while keeping the actual blobs in place, as well as prop osals that c hange the h yp erparameters and shift all of the relev ant blobs in so they are appropriate for the new v alues of the h yp erparameters. Ev aluating the likelihoo d function requires that w e com- pute a “mo c k” image from the current setting of the param- eters. This mo c k image is calculated using standard ray- tracing metho ds with a uniform grid of n × n ra ys fired p er image pixel. How ever, certain kinds of prop osals do not af- fect the image in any wa y , such as those whic h c hange the noise parameters. F or efficiency we do not recompute the mo c k image in these cases. 5 “EASY” SIMULA TED D A T A T o demonstrate the metho d, we generated a simulated dataset where the lens w as an SIE+ γ profile with a sin- gle additional substructure. The image is shown in Figure 2, and consists of 100 × 100 pixels. The p oin t spread func- tion (PSF) was compact and was defined on a 5 × 5 grid. The data was created by firing only one ray p er pixel, and the inference was also carried out under the same lev el of appro ximation. Since the mo del assumptions are all correct for this image, the demonstration here is purely to illustrate the computational tractability of the problem, and the kind of outputs the metho d can pro duce. The SIE+ γ parameters w ere, to three significan t figures, { b, q, x SIE c , y SIE c , θ , γ , θ γ } = { 5 . 25 , 0 . 323 , 0 . 191 , − 0 . 341 , 0 . 833 , 0 . 0322 , 0 . 497 } , where the length units in the image plane are arbitrary , and rotation angles are measured in radians. The single substructure is lo cated at ( − 2 . 57 , 2 . 20), ab out halfwa y b et ween the center and the upp er left image. Its mass is 2.50 units, using the con ven tion where the critical densit y is 1. F or comparison, the SIE mass within its critical ellipse is π b 2 = 86 . 6 units. W e executed the code to generate 5,000 samples from the “mixture of constrained priors” distribution of DNS. Our c  0000 RAS, MNRAS 000 , 000–000 5 Quantit y Meaning Prior Num b ers of Blobs N src Number of source blobs Uniform(0 , 1 , ..., 100) N lens Number of lens blobs Uniform(0 , 1 , ..., 10) Source hyperparameters ( α src ) ( x src c , y src c ) Typical p osition of source blobs iid Cauch y(lo cation=(0 , 0), scale=0.1 × imageWidth ) R src Typical distance of blobs from ( x src c , y src c ) LogUniform(0.01 × imageWidth , 10 × imageWidth ) µ src Typical flux of source blobs ln( µ src ) ∼ Cauch y(0, 1) T ( − 21 . 205 , 21 . 205) W src max Maximum width of source blobs LogUniform(0.001 × imageWidth , imageWidth ) W src min Minimum width of source blobs Uniform(0, W src max ) Lens hyperparameters ( α lens ) ( x lens c , y lens c ) Typical p osition of lens blobs iid Cauch y(lo cation=(0 , 0), scale=0.1 × imageWidth ) R lens Typical distance of blobs from ( x lens c , y lens c ) LogUniform(0.01 × imageWidth , 10 × imageWidth ) µ lens Typical mass of lens blobs ln( µ lens ) ∼ Cauch y(0, 1) T ( − 21 . 205 , 21 . 205) W lens max Maximum width of lens blobs LogUniform(0.001 × imageWidth , imageWidth ) W lens min Minimum width of lens blobs Uniform(0, W lens max ) Source Blob Parameters ( θ src i ) ( x src i , y src i ) Blob position Circular(location=( x src c , y src c ), scale= R src ) A i Blob flux Exponential(mean= µ src ) w i Blob width Uniform( W src min , W src max ) Lens Blob Parameters ( θ lens i ) ( x lens i , y lens i ) Blob position Circular(location=( x lens c , y lens c ), scale= R lens ) M i Blob mass Exponential(mean= µ lens ) v i Blob width Uniform( W lens min , W lens max ) Smo oth Lens P arameters ( θ SIE ) b SIE Einstein Radius LogUniform(0.001 × imageWidth , imageWidth ) q Axis ratio Uniform(0, 0.95) ( x SIE c , y SIE c ) Central p osition iid Cauc hy(location=(0 , 0), scale=0.1 × imageWidth ) θ Orientation angle Uniform(0 , π ) γ External shear Cauch y(0 , 0 . 05) T (0 , ∞ ) θ γ External shear angle Uniform(0 , π ) Noise Parameters ( σ ) σ 0 Constant comp onent of noise variance ln( σ 0 ) ∼ Cauch y(0, 1) T ( − 21 . 205 , 21 . 205) σ 1 Coefficient for v ariance increasing with flux ln( σ 1 ) ∼ Cauch y(0, 1) T ( − 21 . 205 , 21 . 205) Data ( D ) D ij Pixel intensities Normal( m ij , s 2 ij + σ 2 0 + σ 1 m ij ) T able 1. The prior distribution for all hyperparameters, parameters, and the data, in our model. Uniform( a, b ) is a uniform distribution between a and b . LogUniform( a, b ) is a log-uniform distribution (with density f ( x ) ∝ 1 /x , sometimes erroneously called a Jeffreys prior) between a and b . The notation T ( α, β ) after a distribution denotes truncation to the in terv al [ α, β ]. The constan t imageWidth is the geometric mean of the image dimensions in the x and y directions. thinning factor w as 10 5 , so 5 × 10 8 MCMC iterations w ere ac- tually p erformed, taking approximately 48 hours on a mo d- est desktop PC 1 . After resampling these samples to reflect the p osterior distribution, w e w ere left with 500 (equally w eighted) p osterior samples. The p osterior distributions for complex mo dels, suc h as the mixture mo dels used here, are often challenging and un- in tuitive to summarise. One effective wa y to visualise the 1 The computer was purchased in 2012 and has a second- generation in tel i7 processor, and the process was run on 8 threads. uncertain ty in the inferences is to play a movie where each frame is a sample from the p osterior distribution. The de- gree to which the frames differ from each other con veys the uncertain ty remaining after taking the data in to account. F or the purp oses of a pap er, static summaries are more con venien t than movies. One useful summary is based on the concept of empirical measure. The empiric al me asur e of the substructure p ositions is a function that tak es the actual substructure p ositions ( x lens i , y lens i ) and pro duces a “densit y function” ov er tw o dimensions, comp osed of delta functions c  0000 RAS, MNRAS 000 , 000–000 6 Br ewer, Huijser and L ewis Source Blobs i = 1 , ..., N src Lens Blobs i = 1 , ..., N lens α lens σ N src θ SIE θ src i θ lens i α src D N lens Figure 1. A probabilistic graphical model (PGM) of the dep endence structure of the prior information, produced using D AFT ( www.daft-pgm.org ). The prior for the source and lens blob parameters are sp ecified conditional on hyperparameters. The purp ose of this is to induce dep endence in the prior distribution for the blob parameters, implying (for example) that the mass of one blob is slightly informative about the mass of another, and that the lo cations might b e clustered around a certain typical lo cation. Sim ulated Dataset Figure 2. A simulated image of a simple “galaxy” lensed by an SIE+ γ lens close to the center of the image, indicated by the white circle, plus a single substructure close to the top left image (indicated by a star symbol). The image was blurred b y a PSF and had some noise added. In an arbitrary system of units, this image extends from -8 to 8 in the x and y directions. at the p ositions themselves: n ( x lens i , y lens i ) o N lens i =1 = ⇒ N lens X i =1 δ 2  x − x lens i , y − y lens i  . (12) In tuitively , the empirical measure is a mathematical ob- ject that is like an “infinite resolution histogram”, in this case a t wo dimensional histogram, of the substructure p osi- tions. Being a function of the actual substructure p ositions, the empirical measure is not av ailable to us since we do not kno w those p ositions with certaint y . How ever, we ha ve sam- ples from the posterior distribution for those positions, and can use these (trivially) to create samples from the posterior distribution for the empirical measure. W e can also sum- marise this p osterior, for example, by taking its exp ected v alue. The p osterior exp ected v alue of the empirical measure is: Z p ( θ | D ) N lens X i =1 δ 2  x − x lens i , y − y lens i  d θ (13) where θ denotes all parameters and hyperparameters (in- cluding the N s) and “ R d θ ” is an in tegral and summation o ver the en tire parameter space. Since w e can approximate p osterior expectations using Mon te Carlo, w e can obtain the expected v alue of the em- pirical measure using: 1 n n X k =1 N k lens X i =1 δ 2  x − x lens ik , y − y lens ik  (14) where n is the num b er of p osterior samples. The resulting function is an “image” with a point mass wherever a sub- structure occurred. F or visualisation purposes the image can b e blurred, or calculated at a lo wer resolution by replacing c  0000 RAS, MNRAS 000 , 000–000 7 80 81 82 83 84 85 86 87 88 89 SIE mass (within critical ellipse) 0 5 10 15 20 T otal substructure mass Figure 3. The joint posterior distribution for the SIE mass (in- tegrated within its critical ellipse, which is not the critical curve of the lens ov erall), and the total mass in substructures. the Dirac-delta function with a discrete version whic h re- turns a nonzero constan t if a substructure app ears in a pixel or zero otherwise. The masses of the smo oth and substructure comp onen ts of the lens are usually of interest. Since the total mass of an SIE lens is infinite, the question needs to be redefined, so we ask about the mass within some ap erture of finite area. F or an SIE, the mass (in dimensionless units based on the critical densit y) within the critical ellipse is simply π b 2 . Ho wev er, the blobs hav e finite total masses { M i } . T o obtain the posterior distribution for the lens masses, one m ust b e clear ab out exactly which mass they are talking ab out, and man y definitions are p ossible, although some might b e more meaningful or well constrained than others. In the present pap er we do not address the question of exactly which quan tities related to the densit y profile of the lens are most scientifically in teresting. Nevertheless, we can verify that the results of the inference do b eha ve in un- derstandable wa ys. F or example, in Figure 3, we plot the join t p osterior distribution for the SIE mass within its crit- ical ellipse and the total substructure mass o ver the entire domain. As one w ould exp ect, there is a strong dep endence b et ween these t wo quantities in the p osterior distribution, as mass in the SIE can b e traded off with mass in substructures to some extent, while the mo del remains (loosely sp eaking) “consisten t with the data”. The p osterior distribution for N lens is also clearly of in- terest, and is display ed in Figure 4. The prior for this param- eter w as uniform from 0 to 10 (inclusiv e), and the p ossibilit y N = 0 has b een (loosely sp eaking) “ruled out” by the image data. The true solution ( N = 1) has the highest probability . Ho wev er, the possibilities with N > 1 are still fairly plausi- ble, since it’s possible that a very small substructure exists, and it’s also p ossible that tw o substructures near each other could mimic the effect of one. The degree to which these p ossibilities are plausible is related to the choice of prior for the blob amplitudes and p ositions. The estimated marginal likelihoo d of our mo del (av- eraged o ver all parameters including N lens and N src ) is ln [ p ( D )] ≈ − 14141 . 4, and the information (Kullback- Leibler divergence from the prior to the p osterior) is H ≈ 99 . 4 nats. The information represen ts the degree of com- 0 1 2 3 4 5 6 7 8 9 10 N lens 0 10 20 30 40 50 60 70 80 90 Num b er of P osterior Samples Figure 4. The marginal p osterior distribution for N lens , the num- ber of substructures in the lens. The prior was uniform and the true v alue used to generate the data was 1. Ho wev er, the data is only informative enough to suppress the probability of v alues above 1 slightly , since it is p ossible (given this data) that lo w mass substructures might exist somewhere, or that what we think is one substructure migh t actually b e t wo close together, and other suc h possibilities. pression of the p osterior distribution with respect to the prior, and can b e interpreted quite literally as how muc h w as learned ab out the parameters from the data. It is also straigh tforward to estimate from Nested Sampling (Skilling 2006). Its definition is: H = Z p ( θ | D ) ln  p ( θ | D ) p ( θ )  d θ (15) and we can build in tuition ab out its meaning based on some simple examples. One such example is a uniform prior ov er a volume V 0 and a p osterior which is uniform ov er a smaller v olume V 1 con tained within V 0 . In this case H = ln( V 0 /V 1 ) nats. Therefore, a v alue of H = 100 nats implies the p oste- rior distribution o ccupies roughly e − 100 of the prior volume. F or the sampled parameter v alues representing the pos- terior distribution, there was no structure in the residuals (not sho wn), when normalised by the noise standard de- viation in eac h pixel. This is unsurprising since the mo del assumptions are entirely correct for this dataset. 6 “HARDER” SIMULA TED D A T A T o further test the algorithm, we created a “harder” sim- ulated image, where the lens consisted of an SIE+ γ profile plus 10 additional substructures. The image is shown in Fig- ure 6. With the easy dataset, there was little hope of mea- suring the substructure mass function parameter µ lens since there was only a single substructure which w ould provide little information about µ lens (only constraining its general order of magnitude), even if its mass were measured p er- fectly . As with the ‘easy’ dataset, we generated 5,000 sam- ples from the DNS target distribution and thinned b y a factor of 10 5 , so 5 × 10 8 MCMC iterations w ere actually p erformed, taking approximately 60 hours. The slow er run- time in this case was because N lens = 10, more CPU time w as spent exploring parts of the hypothesis space c  0000 RAS, MNRAS 000 , 000–000 8 Br ewer, Huijser and L ewis Substructure P ositions Figure 5. The ‘easy’ sim ulated data, with substructure p ositions ov erlaid (a Mon te Carlo approximation to the p osterior exp ecta- tion of the empirical measure). The density of black points in an y region is prop ortional to the expected num b er of substructures whose centers lie within that region. In this case, there is strong evidence for a substructure close to the top-left image (where one was actually placed). Sim ulated Dataset Figure 6. The ‘harder’ sim ulated data, consisting of a more com- plex source profile, and a lens with an SIE+ γ (whose cen tre is indicated b y a white circle) and ten substructures. Only fiv e of the substructures (star sym b ols) occured within the square image boundary , but the other five still hav e an effect on the image. where N lens is high. The run resulted in 557 equally- w eighted posterior samples. The posterior distribution for the SIE mass and the total substructure mass is giv en in Figure 7. The samples can b e view ed as a mo vie at www.youtube.com/watch?v=o3ppfKSk248 . T o demonstrate the claim that we can infer something ab out the mass function of substructures directly from the image data, we ha ve plotted the p osterior distribution for µ lens , the hyperparameter for the substructure masses, in Figure 8. The true v alue w as 1, and the posterior distribu- tion indicates a fairly wide uncertaint y range that is at least 60 65 70 75 80 85 SIE mass (within critical ellipse) 5 10 15 20 25 30 T otal substructure mass Figure 7. Same as Figure 3, but for the ‘harder’ simulated dataset. 0 . 0 0 . 5 1 . 0 1 . 5 2 . 0 2 . 5 3 . 0 3 . 5 4 . 0 4 . 5 µ lens 0 10 20 30 40 50 Num b er of p osterior samples Figure 8. The p osterior distribution for µ lens , the hyperparam- eter which determines the prior expected v alue of substructure masses, given the ‘harder’ simulated data. The true value (in ar- bitrary units) was 1, which is indicated by the vertical line. consisten t with the true v alue in some sense. Given that the image only con tained a few substructures whose masses could b e measured with any accuracy , the large uncertaint y is not surprising. The estimated marginal likelihoo d of our model, for the ‘harder’ sim ulated data, is ln [ p ( D )] ≈ − 5671 . 9, and the information (Kullbac k-Leibler divergence from the prior to the p osterior) is H ≈ 155 . 7 nats. As with the easy dataset, the standardised residuals of the sampled mo dels resembled noise. 7 THE COSMIC HORSESHOE As a further demonstration the method, w e apply it to the g- band data of the Cosmic Horsesho e J1004+4112 (Belokurov et al. 2007; Dye et al. 2008) taken with the Isaac Newton T elescop e (INT). The image is sho wn in Figure 11. F or this system, more data is a v ailable (i and U band data, as well as more recent Hubble Space T elescop e (HST) imaging), al- though the g-band image is the highest signal-to-noise of c  0000 RAS, MNRAS 000 , 000–000 9 0 1 2 3 4 5 6 7 8 9 10 N lens 0 50 100 150 200 250 Num b er of P osterior Samples Figure 9. Same as Figure 4, but for the ‘harder’ simulated dataset. The true num b er of substructures was 10, which also happens to be the p osterior mo de. Only five substructures were located within the domain of the image, but N lens = 10 is prob- able b ecause the spread-out p ositions of the substructures imply that R lens is probably large. Substructure P ositions Figure 10. Same as Figure 5, but for the ‘harder’ simulated dataset. F or all of the substructures located close to the ring, the p ositions were well constrained. There was weak evidence for some substructures that did not in fact exist. the INT images. The source redshift is z s = 2 . 379, and the lens is a massiv e luminous red galaxy at z l = 0 . 4457. Unfortunately our current implementation do esn’t allow for m ulti-band data (unlike Brewer et al. (2011)), and is quite slo w when running on the larger HST image. Therefore this section should b e considered a further demonstration of the tec hnique, and not a thorough study of this system. W e generated 5,000 samples from the DNS target dis- tribution and thinned by a factor of 10 5 , so 5 × 10 8 MCMC iterations were actually performed, taking approximately 40 hours. After resampling, this resulted in 718 equally w eighted posterior samples. Although this image is smaller than the simulated data (46 × 46 pixels), we fired 2 × 2 rays p er pixel for greater accuracy . The joint p osterior distribution for the SIE mass total substructure mass for the Cosmic Horsesho e is given in Fig- -7.66 0 7.66 x (arcseconds) -7.66 0 7.66 y (arcseconds) Cosmic Horsesho e Figure 11. The g-band INT image of the Cosmic Horseshoe, with the lens galaxy subtracted. ure 12. As with the simulated data (Figure 3), we see an exp ected negative correlation betw een these tw o quantities. Ho wev er, the p oin ts are not as smoothly distributed in this case (we discuss this issue further in Section 7.1). An upp er “limit” on the SIE mass, with 95% p osterior probabilit y , is 5 . 20 × 10 12 solar masses. The algorithm has also found some p ossibilities where the substructure mass is m uch greater than this. In these cases, the substructures are far from the image — it is the environmen t that is b eing mo delled. The p osterior distribution for N lens (Figure 13) also sho ws some evidence for more than zero substructures. In particular, the prior probabilit y for N lens > 0 was 10/11 ≈ 0.91, whereas the posterior probability is appro ximately 0.997. This corresp onds to a Ba yes F actor of ab out 300 in fa vour of N lens > 0 v ersus the alternative N lens = 0. Despite w eak evidence for the existence of substruc- ture, when w e examine the exp ected v alue of the empirical measure of substructure p ositions (Figure 14) we find no consistency in their positions, unlike for the simulated data (Figure 5). This situation is not uncommon. F or example, Brew er & Dono v an (2015) found evidence for a large n umber of Keplerian signals in a time series, but the num b er of such signals with w ell constrained properties w as m uch low er. Re- lated to this, we can compute p osterior probabilities for any h yp othesis ab out the substructure masses. With 95% prob- abilit y the substructure mass is less than 1.53 × 10 13 solar masses and with 75% probabilit y it is less than 5.10 × 10 12 solar masses. As the Einstein ring itself implies a mass of ab out 5 × 10 12 solar masses, these summaries are affected b y the p ossibilit y of blobs outside the main Einstein ring. As with any inference, the results presented here may b e sensitiv e to many of the input mo delling assumptions, and a slightly differen t (yet still reasonable) set of c hoices migh t yield differen t results. F or example, if the densit y pro- file of the lens was smo oth but not of the SIE form, and the data were informative enough to show this, the curren t mo del would only b e able to explain the data b y adding substructures. Hence, an alternative explanation for these c  0000 RAS, MNRAS 000 , 000–000 10 Br ewer, Huijser and L ewis 0 1 2 3 4 5 6 7 SIE mass (within critical ellipse) ( M  ) × 10 12 0 1 2 3 4 5 6 7 T otal substructure mass ( M  ) × 10 12 Figure 12. Same as Figure 3, but for the Cosmic Horseshoe. The joint p osterior distribution for the SIE mass (in tegrated within its critical ellipse, which is not the critical curve of the lens ov erall), and the total mass in substructures. The units are defined by the critical density , so that a mass of π units would have an Einstein radius of one arcsecond. T o calculate the masses in solar masses, we assumed a flat cosmology with Ω m = 0 . 3, Ω Λ = 0 . 7, and H 0 = 70 km/s/Mpc. 0 1 2 3 4 5 6 7 8 9 10 N lens 0 20 40 60 80 100 120 140 160 180 Num b er of P osterior Samples Figure 13. Same as Figure 4, but for the Cosmic Horseshoe. The posterior distribution for the num b er of lens substructures in the Cosmic Horsesho e system. There is moderate evidence in fav our of the hypothesis that N lens 6 = 0. results is that, rather than containing a substructure, the densit y profile of the lens is simply not within the SIE+ γ family . In fact, some of the p osterior samples obtained con- tain massiv e substructures outside of the image, whic h could b e mo delling higher-order effects of the en vironment b ey ond what is captured by the simple external shear model. This is one wa y of violating the SIE+ γ assumption, and another is simply to ha ve a different pro jected densit y profile. One w ay of further inv estigating this is to use a different family of smooth lens models and doing mo del selection based on the marginal likelihoo ds. Another option whic h we defer to future work is to implement a more flexible model (suc h as Substructure P ositions Figure 14. Same as Figure 5, but for the Cosmic Horseshoe. The positions of substructures encountered in the MCMC sampling for the Cosmic Horseshoe; technically a Mon te Carlo approximation to the p osterior exp ected v alue of the empirical measure of sub- structure p ositions. Unlike Figure 5, for the Cosmic Horsesho e there is little consistency in the positions of the substructures. a mixture of concentric elliptical lenses). This result is con- sisten t with that of Dye et al. (2008), who found a slight preference for an elliptical p o wer la w profile ov er the SIE (whic h is a sp ecific case of an elliptical p o wer law). Another p oten tial inadequacy of our mo del is the as- sumption that the PSF is kno wn. In practice, the PSF is often estimated from the image of nearby star(s), or from theoretical kno wledge of the telescope (esp ecially in the case of space based imaging). Ho wev er, if the PSF is misspecified, or in fact v aries across the image, this could induce subtle effects in the imaging which our current mo del could only explain using substructure. The estimated marginal likelihoo d of our mo del is ln [ p ( D )] ≈ − 6948 . 5, and the information (Kullback-Leibler div ergence from the prior to the posterior) is H ≈ 123 . 1 nats. Three example mo dels (lens and source) sampled from the p osterior are shown in Figure 15. The mo del w as able to fit the data do wn to the noise level. 7.1 Repro ducibilit y of the results The “effectiv e sample size” returned b y DNS, whic h we hav e describ ed as the num b er of p osterior samples, takes into ac- coun t the fact DNS’s target distribution is not the p oste- rior. How ever, it do es not tak e auto correlation in to account, and th us can presen t an optimistic picture of the accuracy of an y Monte Carlo approximations to p osterior quantities. Most standard diagnostic techniques used in MCMC can b e applied here. The simplest of these is a c hec k of repro ducibil- it y . If different runs yield substantially different results, the MCMC output should be treated with caution. F or exam- ple, in Figure 12, there is a correlation b et ween the mass attributed to the SIE comp onen t and that attributed to substructures. How ever, the distribution also app ears some- what “lump y” or multimodal. This may not b e a feature of the actual p osterior distribution, but could arise due to c  0000 RAS, MNRAS 000 , 000–000 11 Source Densit y in substructures Image Figure 15. Three example sources, the density profile in substructures, and the corresp onding lensed, blurred images, representativ e of the posterior distribution. The angular scales are 7.5 arcseconds for the sources and 15 arcseconds for the images. imp erfect sampling. Whereas a standard Metropolis sampler w ould mov e around very slowly in the hypothesis space, DNS naturally sp ends a non-negligible fraction of the time sam- pling the prior. Therefore, a particle exploring the parameter space can “forget” its go od-fitting position, mov e somewhere completely different, and find another go od-fitting mo del in a different lo cation. This is a natural feature of the algo- rithm (and is also present in related algorithms suc h as par- allel temp ering). Despite this, the patc hy nature of Figure 12 suggests that posterior exploration is still challenging in this problem. The results for the simulated data (Figure 3) w ere less complex because of the definite existence of one substruc- ture. Its mass and p osition provide some information ab out the hyperparameters α lens , restricting the probability that high mass substructures exist far from the image. 8 CONCLUSIONS W e ha ve developed a trans-dimensional Bay esian approac h motiv ated directly by the question of whether dark substruc- tures exist in a lens galaxy , given image data. By making use of the Diffusive Nested Sampling algorithm (Brew er, P´ artay , & Cs´ an yi 2011) and the RJObject library (Brewer 2014), we outsource the difficulties asso ciated with choosing Metrop o- lis proposals for a hierarchical mo del of non-fixed dimension. The mo del allo ws for source and lens “blobs” to app ear as needed to explain the data. The prior for the blobs’ prop- erties is specified hierarc hically , to more realistically mo del sensible prior beliefs, and to tie the mo del parameters di- rectly to questions of scientific interest such as the mass function of substructures. As a pro of of concept, w e demonstrated the successful reco very of a single substructure from simulated data for whic h all of the mo del assumptions were true. W e then ap- plied the method to an image of the Cosmic Horsesho e sys- tem and found mo derate evidence in fav our of the existence of a substructure, or at the v ery least, a departure from a simple SIE+ γ lens profile. The main output of our method is a set of p osterior samples, each representing a plausible scenario for the lens and the source giv en the data and the assumed prior information. These samples can be display ed as a movie, whic h is a very useful metho d for intuitiv ely un- c  0000 RAS, MNRAS 000 , 000–000 12 Br ewer, Huijser and L ewis derstanding the remaining uncertain ty . Any p osterior sum- maries of in terest can b e produced from the samples, but the in terpretation of many of these summaries is not necessarily straigh tforward. W e ha ve suggested that the posterior ex- p ectation of the empirical measure of substructure positions is one of the most helpful summaries. The mo del assumptions used in this pap er are fairly simple. In future w ork we intend to generalize the mo del to multi-band data, and extend the smooth lens mo del be- y ond the ov erly simplistic SIE+ γ assumption. Applications to other systems are also forthcoming, as well as v ariants on the mo del using different conditional priors for the sub- structure prop erties giv en the hyperparameters. A CKNOWLEDGEMENTS It is a pleasure to thank Phil Marshall (Stanford), T om- maso T reu (UCLA), Ross F adely (NYU), Alan Heav ens (Ed- in burgh), and the referee for v aluable discussion. Thomas Lumley (Auc kland) also deserves credit for ending BJB’s us- age of the jet colormap. This w ork was funded by a Marsden F ast Start grant from the Roy al Society of New Zealand. REFERENCES Alsing J., Heav ens A., Jaffe A. H., Kiessling A., W andelt B., Hoffmann T., 2015, arXiv, Auger M. W., T reu T., Brewer B. J., Marshall P . J., 2011, MNRAS, 411, L6 Bark ana R., 1998, ApJ, 502, 531 Belokuro v V., et al., 2007, ApJ, 671, L9 Birrer S., Amara A., Refregier A., 2015, arXiv, Brew er B. J., Lewis G. F., 2006, ApJ, 637, 608 Brew er B. J., P´ arta y L. B., Cs´ anyi G., 2011, Statistics and Computing, 21, 4, 649-656. Brew er B. J., Lewis G. F., Belokurov V., Irwin M. J., Bridges T. J., Ev ans N. W., 2011, MNRAS, 412, 2521 Brew er, B. J., 2014, preprint. ArXiv: 1411.3921 Brew er B. J., Donov an C. P ., 2015, MNRAS, 448, 3206 Catic ha, A. 2008. Lectures on Probability , En trop y , and Statistical Ph ysics. ArXiv e-prints Coles J. P ., Read J. I., Saha P ., 2014, MNRAS, 445, 2181 Dehnen W., Aly H., 2012, MNRAS, 425, 1068 Dy e S., Ev ans N. W., Belokurov V., W arren S. J., Hewett P ., 2008, MNRAS, 388, 384 F adely R., Keeton C. R., 2012, MNRAS, 419, 936 Green, P . J., 1995, Reversible Jump Marko v Chain Monte Carlo Computation and Ba yesian Mo del Determination, Biometrik a 82 (4): 711–732. Grillo C., et al., 2015, ApJ, 800, 38 Hupp enk othen D., et al., 2015, ApJ, 810, 66 Jones D. E., Kashy ap V. L., v an Dyk D. A., 2014, arXiv, Ko opmans L. V. E., 2005, MNRAS, 363, 1136 Kormann R., Schneider P ., Bartelmann M., 1994, A&A, 284, 285 Millar, R. B., Maximum likelihoo d estimation and infer- ence: with examples in R, SAS and ADMB. V ol. 111. John Wiley & Sons, 2011. O’Hagan, A., F orster, J., 2004, Bay esian inference. London: Arnold. P ancoast A., Brewer B. J., T reu T., 2011, ApJ, 730, 139 Sc hneider M. D., Hogg D. W., Marshall P . J., Dawson W. A., Mey ers J., Bard D. J., Lang D., 2015, ApJ, 807, 87 Sivia, D. S., Skilling, J., 2006, Data Analysis: A Bay esian T utorial, 2nd Edition, Oxford Universit y Press Skilling, J., 2006, “Nested Sampling for General Bay esian Computation”, Bay esian Analysis 4, pp. 833-860 Suyu S. H., Marshall P . J., Hobson M. P ., Blandford R. D., 2006, MNRAS, 371, 983 Suyu S. H., et al., 2014, ApJ, 788, L35 Suyu S. H., et al., 2013, ApJ, 766, 70 T agore A. S., Jackson N., 2015, arXiv, T reu T., 2010, ARA&A, 48, 87 Umst¨ atter R., Christensen N., Hendry M., Mey er R., Simha V., V eitch J., Vigeland S., W oan G., 2005, PhRvD, 72, 022001 V egetti S., Ko opmans L. V. E., 2009, MNRAS, 400, 1583 V egetti S., Lagattuta D. J., McKean J. P ., Auger M. W., F assnach t C. D., Ko opmans L. V. E., 2012, Natur, 481, 341 V egetti S., Czoske O., Koopmans L. V. E., 2010, MNRAS, 407, 225 V egetti S., Ko opmans L. V. E., Bolton A., T reu T., Ga v azzi R., 2010, MNRAS, 408, 1969 V egetti S., Ko opmans L. V. E., Auger M. W., T reu T., Bolton A. S., 2014, MNRAS, 442, 2017 W arren S. J., Dye S., 2003, ApJ, 590, 673 c  0000 RAS, MNRAS 000 , 000–000