A new parameter space study of cosmological microlensing
Abstract
Cosmological gravitational microlensing is a useful technique for understanding the structure of the inner parts of a quasar, especially the accretion disk and the central supermassive black hole. So far, most of the cosmological microlensing studies have focused on single objects from 90 currently known lensed quasars. However, present and planned allsky surveys are expected to discover thousands of new lensed systems. Using a graphics processing unit (GPU) accelerated rayshooting code, we have generated 2550 magnification maps uniformly across the convergence () and shear () parameter space of interest to microlensing. We examine the effect of random realizations of the microlens positions on map properties such as the magnification probability distribution (MPD). It is shown that for most of the parameter space a single map is representative of an average behaviour. All of the simulations have been carried out on the GPUSupercomputer for Theoretical Astrophysics Research (gSTAR).
keywords:
gravitational lensing: micro – accretion, accretion discs – quasars: general1 Introduction
Understanding the properties of the supermassive black holes located at the centres of most galaxies (Ferrarese & Ford, 2005), is essential to our understanding of the evolution of the Universe. The energy output observed in quasars, L (Peterson, 1997), is believed to stem from the interaction of a central supermassive black hole with a surrounding accretion disc. By studying accretion discs it is then possible to infer properties of the supermassive black holes, which in turn could constrain theories on their origin and evolution.
While plausible accretion disc models exist (see Abramowicz & Fragile, 2013, for a review), current observations are not sufficient to fully test them. This is mainly because the angular diameter of the accretion disc at typical quasar distances, arcsec at redshift z2, is orders of magnitude smaller than the resolution of our current best telescopes. However, by detecting the imprints of gravitational microlensing over such cosmological scales we can extract useful information from the observations that can constrain accretion disc and, ultimately, black hole properties (Rauch & Blandford, 1991; Morgan et al., 2010).
Light propagating from a distant quasar can be gravitationally lensed by a foreground galaxy resulting in the creation of luminous arcs, Einstein rings, multiple images, magnifications, etc. In the case of an image of a background source being seen through the bulge of a foreground galaxy, individual stars that lie close to the line of sight can further deflect the light rays, acting as gravitational microlenses. Their effect can be observed as uncorrelated variability in the light curve of a pair of images (Paczynski, 1986), or as anomalous flux ratios between the multiple images (Witt et al., 1995; Schechter & Wambsganss, 2002). Quasar microlensing was originally proposed by Chang & Refsdal (1979), and later detected by Irwin et al. (1989). For a recent review of the progress in this field see Schmidt & Wambsganss (2010).
The collective effect of microlensing by individual stars in a foreground galaxy is most often investigated using a magnification, or caustic, map; a 2dimensional pixellated map showing the magnification on a finite region of the source plane (Fig. 1). The relative velocity of the observer, the source and the lens galaxy can take a background source from regions of low to high magnification on this map, inducing variability in the light curve. The rate and amplitude of this variability is related to the size of the source with respect to the Einstein radius, , of the microlenses in the lensing galaxy; small sources will display larger and more rapid variations while extended ones will vary more smoothly (Wambsganss et al., 1990; Schmidt & Wambsganss, 2010). Since quasar accretion discs are much smaller than the corresponding Einstein radii in all known systems ( cm compared to cm, see Mosquera & Kochanek, 2011), they will be subject to microlensing essentially all the time.
Lightcurveanalysis methods require continuous monitoring of a system for periods from a few months to several years (Mosquera & Kochanek, 2011). Such monitoring has been practical at some level for 20 to 30 systems (for some examples see Eigenbrod et al., 2007; Morgan et al., 2010, and references therein). Alternatively, because gravitational microlensing is achromatic, single epoch multiwavelegth observations (snapshots) can estimate the accretion disc size in various wavelengths i.e. the temperature profile of the disc: (Blaes, 2004; Abramowicz & Fragile, 2013).
The difficulties involved in both these methods of observing microlensed systems e.g. monitoring over large periods of time, coordinating observations at various wavelengths, measuring the lens and source redshifts, measuring timedelays etc., have kept the number of microlensed systems low; only 23 systems have been studied in detail to date (Bate & Fluke, 2012, hereafter BF12). Most studies have focused on single objects in the effort to constrain accretion disc properties, and only recently have collections of objects been considered (Morgan et al., 2010; Blackburne et al., 2011; Mosquera et al., 2011; Sluse et al., 2012). However, future planned allsky surveys, like the Large Synoptic Survey Telescope (LSST; LSST Science Collaborations et al., 2009), are expected to discover hundreds, or even thousands, of multiply imaged systems (Oguri & Marshall, 2010).
Studying the parameter space of cosmological microlensing will be useful in understanding how the microlensing characteristics of a system can help constrain underlying accretion disc models. Of particular interest is understanding how the properties of magnification maps, such as the magnification probability distribution (MPD), depend on the physical parameters convergence () and shear (). However, there are a number of additional parameters related to properties of the microlenses (mass, positions) and map characteristics (resolution, size, accuracy). The imminent transition from single objects to statistically meaningful samples of microlensed quasars makes a parameter space approach to quasar microlensing timely (BF12).
1.1 Explorations of parameter space
A first exploratory study of microlensing parameter space was conducted by Wambsganss (1992, hereafter W92). Although the resolution of the maps was low [e.g. to allow for an investigation of the power law behaviour of the MPD, , for high magnifications (Schneider, 1987; Witt, 1990)], W92 provided some early evidence that the MPD is independent of the mass distribution of the microlenses. Lewis & Irwin (1995, hereafter LI95) performed a similar parameter space exploration confirming this result. A summary of the parameter values used in W92, LI95 and the present parameter space study (hereafter VF13) is shown in Table 1.
Since W92 and LI95, no attempt to systematically explore parameter space has been made, mainly due to the lack of enough computational power to do so (W92). Moreover, taking into account the low number of systems that exhibit microlensing, there has been little motivation in carrying out parameter space studies. As a result, systems have been modelled individually. However, the advent of supercomputers based on graphics processing units (GPUs) has enabled the systematic study of magnification maps across the entire microlensing parameter space (BF12).
The assumption of using a single map per combination of parameters, the usual approach in most microlensing studies, has not been examined in detail with respect to the random microlens positions. Varying the microlens positions is an important test for the consistency of the existing microlensing techniques and for the robustness of the resulting accretion disc constraints. Future highresolution parameter space explorations, like the one proposed in BF12, should take into account any dependence of map properties on the microlens positions.
W92  LI95  VF13  

Method  tree code  contours  direct 
27  29  170  
16  1  1  
4  1  1  
5  5  1  
1  1  15  
Total  64  61  2550 
In this work, we determine the number of magnification maps needed to fairly sample statistical properties, such as the MPD, for a given set of microlensing parameters. We use the GPUD (Thompson et al., 2010) brute force rayshooting code to uniformly cover the largest part of the parameter space. We provide 170 combinations, producing 15 maps with different random microlens positions for each combination. In Section 2 we introduce the direct inverse ray shooting technique and the microlensing parameter space. In Section 3 we present our results, statistical properties of the MPD are extracted and an average behaviour per parameter space () pair is presented. Discussion and conclusions follow in Sections 4 and 5.
2 Method
At the heart of each microlensing technique lies the caustic magnification pattern created by compact objects and smooth matter in the lensing galaxy (Kayser et al., 1986; Schneider, 1987). Once provided with a suitable map, one can study the background source (with a given size and/or profile) by either generating model lightcurves or by calculating the MPD. Key parameters for generating a magnification map are the convergence, , and shear, . Convergence describes the focusing power of the microlenses and has two components, one from matter in compact objects (microlenses) , and another from smoothly distributed matter, . Shear describes the external distortion applied to the images, which is dominated by the mass distribution of the lens galaxy. The values of and arise from the macromodelling of the observed multiple images and the lensing galaxy (Keeton, 2001). One can then use the lens equation to obtain the map:
(1) 
where is in the image plane, is in the source plane, and are the mass and position of each microlens. The number of microlenses, , is related to through:
(2) 
for microlenses with mean mass , uniformly distributed over a surface S (see Bate et al., 2010, for an explanation of the choice of S). The magnification is , with A being the Jacobian of the transformation from to .
It was realized early on that analytic solutions of equation (1) would exist only for a limited number of lens models (Paczynski, 1986). In the case of quasar microlensing, can vary from tens to hundreds of thousands of microlenses, therefore, numerical methods had to be developed. Most of these methods are based on the inverse rayshooting technique introduced by Kayser et al. (1986) and Schneider (1987). Light rays are propagated backwards from the observer, through the lens plane where they are deflected using equation (1) and then mapped on to the source plane. The source plane is divided into a grid of pixels. By counting the number of rays reaching each pixel, , and comparing it to the number of rays that would have reached each pixel if there was no lensing, , one obtains an estimate of the magnification:
(3) 
i.e. a pixelated magnification map.
A number of numerical methods to generate magnification maps exist in the literature. Wambsganss (1990) used a hierarchical tree method to replace microlenses at a given distance from a light ray with a higher mass pseudolens, effectively reducing the final number of microlenses in the calculation. Kochanek (2004) allowed for the lens and source planes to be periodic and based his method on Fourier techniques (similar to algorithms used in Nbody problems), which reduced the number of microlenses as well. Mediavilla et al. (2006, 2011) present an efficient method of inverse polygon mapping from the lens to the source plane which leads to the same accuracy as the previous approaches using far less rays (same accuracy achieved in less than 1 per cent of the computational time, Mediavilla et al., 2006). Finally, Thompson et al. (2010) made use of the high degree of parallelization inherent in the lens equation with a bruteforce solution on a modern GPU. Bate et al. (2010) compared compute times for this direct GPU approach and the more sophisticated treecode approach and found it to be highly competitive (for a single CPU core), especially in regions of the parameter space where the number of microlenses is lower than .
2.1 Microlensing map parameters
The computation of the caustic structure in a microlensing magnification map depends on eight main parameters, which can be categorized into three groups:

Macromodel (external) parameters: and with physical interpretation as explained in Section 2, and the smooth matter fraction, , which describes the smooth matter component.

Map characteristics and statistics: the map size (i.e. side length) should be large enough to enable simultaneous accretion disc and broad emission line region modelling. The map resolution^{1}^{1}1we are following the notation in our previous papers and use the imaging term of resolution: the total number of pixels in an image dimension. should be high enough to resolve the inner parts of an accretion disc (see BF12 for a thorough discussion). The average number of rays per map pixel, , is a measure of the accuracy of a map and in general it is method dependent. For example, the polygon mapping technique (Mediavilla et al., 2011) uses more rays near the caustics and fewer away from them, where extra accuracy is not needed.

Parameters of the microlenses: the microlens mass, , appears explicitly in the lens equation (equation 1), while the mean microlens mass, , is involved in the calculation of the total number of microlenses (equation 2). Microlensing studies can be performed once a mass function is chosen for stars in the lens galaxy. The microlens positions, in equation (1), cannot be measured by observations. Therefore, random realizations of the microlens positions have been used in generating magnification maps.
In most existing studies, only a single set of microlens positions is used under the assumption that different microlens configurations produce statistically equivalent maps. In the present study we investigate the validity of this assumption by generating 15 maps with different sets of random microlens positions on the parameter space. We now describe our choice of the remaining parameter values.
A complete treatment of the stellar mass function introduces several additional parameters (initial mass function, minimum mass cutoff, maximum mass etc). Such a treatment is more practical in the case of single object studies; including it in a parameter space approach would result in additional computations. Nevertheless, it is expected that the MPD will be relatively insensitive to the choice of the stellar mass function over most of parameter space (W92; LI95; Wyithe & Turner 2001), except in some specific cases e.g. a highly magnified negative parity macroimage, as the one explored in Schechter et al. (2004) for . Therefore, we adopt the simplest treatment of the stellar mass function viz. a constant mass = 1 M for all microlenses.
The total number of rays, for which we are calculating the deflections in equation (1), is . This results in an average number of rays per pixel, , among the 2550 maps, equal to . Bate et al. (2010) and BF12 consider this number of rays to be sufficient for the statistical accuracy of the maps. We performed a test of varying the map accuracy in selected regions of the parameter space, presented in Appendix A. The results indicate that the suggestion of Bate et al. (2010) and BF12 is a reasonable choice.
We choose a size of 24 for our maps and a resolution of 4096 pixels, following the strategy proposed in BF12. This means that we can resolve sources down to 0.006 in size. If the size of a map is small compared to the typical size of a caustic (), for fixed resolution, then a particular caustic structure may dominate the map and the corresponding MPD (see fig. 1 in Wambsganss et al., 1990, and panels 6164 of fig. 1 in W92). Therefore, varying the microlens positions on smaller maps is expected to have a larger effect.
In our present results we have assumed all matter to be in the form of compact objects i.e. for all maps. For fixed , is maximized for ; i.e. for any given , adding a contribution from smooth matter will reduce which is proportional to the number of microlenses (equation 2). In this case, the statistical equivalence of maps with different microlens positions can be matched to the cases without smooth matter and a lower .
2.2 The parameter space
The location in space of existing macromodels of multiply lensed systems and microlensing parameter space studies is presented in Fig. 2. The grey circles show the results from the analysis of 85 macromodels of 23 known multiply imaged systems, as compiled and presented in BF12. Each of these 85 models corresponds to two or four pairs, depending on whether the modelled system has a doubly or quadruply imaged background quasar. Parameter space values investigated by W92 and LI95 are shown as red triangles and yellow squares respectively.
The macromagnification, , for a given pair is:
(4) 
with when . This is the critical line (black line in Fig. 2), which divides the parameter space into three distinct regions where macroimages can form: at the minima (); saddlepoints (); and maxima () of the Fermat traveltime surface (Blandford & Narayan, 1986). Images found in these regions of parameter space exhibit different characteristics (Schechter & Wambsganss, 2002), e.g. minima are always magnified (), while saddlepoints and maxima can be demagnified (). None of the 23 systems in the BF12 compilation has an observed maximum image, also called a central or odd image, as can be seen in Fig. 2. Central images occur near the centre of the lensing galaxies, where the optical depth is high () and are thus subject to high demagnifications. Although this makes their detection harder, central images may prove important for quasar microlensing (Dobler et al., 2007).
Fig. 2 shows the present status of macromodelling of the 23 multiply imaged systems that are suitable for microlensing studies. Once future synoptic allsky surveys start to discover and monitor new multiply imaged quasars, new points will be added to this figure. Moreover, alternative macromodels of existing or new systems can also increase the coverage of parameter space e.g. SDSS J0924+0219 is represented with 13 models (Keeton et al., 2006; Morgan et al., 2006; Mediavilla et al., 2009).
3 Results
We have used the GPUD code of Thompson et al. (2010) to generate 15 maps for each of 170 , pairs: , , with . Our choice of parameter space fully includes the critical line where the computations are more demanding. As there have been no previous attempts at wholeparameterspace studies, we focused our effort on the more complex and challenging part of the parameter space along the critical line. From the compilation of 85 lens macromodels presented in BF12, 95 per cent are found within our selected ranges of .
All our simulations were carried out at gSTAR, the GPU Supercomputer for Theoretical Astrophysics Research, based at Swinburne University of Technology. gSTAR made available to us a total of 50 compute nodes, each containing 2 NVIDIA Tesla C2070 GPUs which perform at greater than 1 Tflop/s (single precision). Our dataset consists of a total of 2550 magnification maps at 4096 resolution, corresponding to 0.16 TB of data and 5100 GPU hours on gSTAR’s C2070 cards.
3.1 Probability distributions
One way of comparing different realizations of a magnification map is through the MPD. By eye, the caustic structure in the three magnification maps shown in Fig. 1 appears subtly different even though they all share the same macromodels. E.g. the map on the left is almost completely covered by caustics, while the caustics on the map on the right appear clustered in horizontal bands. The argument that different microlens positions lead to statistically equivalent maps means that although the maps in Fig. 1 appear different, their MPDs should be similar. In Fig. 3 we show the MPDs for 15 individual maps generated for , where each map has different random microlens positions. As we are covering only a limited range of (due to the finite number of pixels in a map), our MPDs are normalized as:
(5) 
where is the probability to find a pixel of magnification in the map, and is the total number of discrete values. Each of the 15 MPDs is scaled by and then binned in 100 logarithmic bins (grey points). All 15 of the MPDs are similar for high magnifications, while one of them, shown in magenta, is different for low magnifications (log). This means that at least one of the 15 maps for the same may not be statistically equivalent to the rest. Because the MPD is used in deriving accretion disc constraints (e.g. anomalous flux ratios), using a single map may lead to biased results. In the following sections, we examine the statistical equivalence of the maps across the parameter space in more detail.
The thick black line in Fig. 3 is the average over the 15 individual distributions, and we refer to it in the following as the mean MPD or . The orange shaded area shows the and the light blue shaded area the limits, where is the standard deviation of the in each bin. The vertical dashed line indicates the limit for no magnification, i.e. .
3.2 Statistical properties
In order to tackle the problem of statistical equivalence between the MPDs at given values for different microlens positions, we now examine their statistical properties. We calculate the mean, the mode, the standard deviation, the median and the skewness for each of the 15 MPDs, for each pair.
The panels on the lefthand column of Fig. 8 show the mean (noted as angle brackets enclosing the property under consideration) of each statistical property, while the righthand column shows the corresponding standard deviation, s, among the 15 MPDs at each position. In Figs 8(a), (b), (c) and (d) we show the mean of the mean, , the mode, Mo, the standard deviation, , and the median, Md, of the MPDs, measured in units of the theoretical magnification (equation 4). The corresponding standard deviation of these properties are shown in Figs 8(a’), (b’), (c’) and (d’) measured in units of as well. Figs 8(e) and (e’) show the logarithm of the mean and the standard deviation of the skewness, , which is a dimensionless quantity.
In Fig. 8(a) we use a reference value of (white on the color scale). Greater values () are shown in red and lower () in blue. We do the same for the mode, Fig. 8(b), but only the blue part of the colorbar is used because the mode is always found to be lower than . In Figs 8(c), (d) and (e) low values are shown in white and high values in red, but the range differs in each case. Skewness was found to be always positive, with per cent of the cases having , quite higher than the rest, therefore, we are showing the logarithm of the skewness and the corresponding standard deviation in Figs 8(e) and (e’). A consistent color scale (white to red) is used for all the standard deviation plots (panels on the right column of Fig. 8), starting always from zero and having a varying maximum value.
3.3 Kolmogorov–Smirnov tests
Along with calculating standard statistical properties of the MPDs, we also perform Kolmogorov–Smirnov (KS) tests between pairs of MPDs and between MPDs and the corresponding .
In the first case, the null hypothesis is that the MPDs are sampled from the same underlying distribution. For each value we have 15 MPDs, which combine into 105 unique pairs. We perform the KS test for each pair; if the resulting pvalue is less than 0.05, then that particular pair of MPDs failed the test i.e. the MPDs are from different underlying distributions. In Figs 10(a), (b) and (c) we show the percentage of pairs that failed the KS test, , across the parameter space, for MPDs binned in 100, 150 and 200 bins respectively. The color scale ranges from 0 (white) to 100 (black) per cent.
In the second case, the null hypothesis is that each MPD and the corresponding are originally from the same underlying distribution. If the pvalue of the test is less than 0.05, then the individual MPD failed the test with the . The results for MPDs and binned in 100, 150 and 200 bins are shown in Figs 10(a’), (b’) and (c’). The color code indicates the actual number of MPDs that failed the KS test with the mean, ranging from 0 (white) to 15 (black).
We perform the KS test between pairs of MPDs while also varying the side length of the map area from which the MPDs are calculated (with the number of bins fixed to 100). From our 24 maps, we are selecting central areas of widths equal to 20 , 16 and 12 . The number of pixels per stays the same (170) in each case as in our original maps. In Figs 11(a), (b) and (c) we perform the KS test between pairs of MPDs for areas of our maps with widths equal to 20 , 16 and 12 respectively. The color scale ranges from 0 (white) to 100 (black) per cent. In Figs 11(a’), (b’) and (c’) we perform the KS test between the MPDs of the smallersized maps and the of the 24 maps, which is our best approximation of the MPD for each value. The color code indicates the actual number of the smallersize map MPDs that failed the KS test with the of the 24 maps, and it ranges from 0 (white) to 15 (black).
Finally, we have chosen 5 pairs of values across the parameter space (the same which appear in Fig. 12) and generated maps with reduced width and resolution. We have set the width and resolution equal to 20 and 3413 pixels, 16 and 2730 pixels and 12 and 2048 pixels (the number of pixels per has been kept the same as our original 24 and 4096pixel maps), and generated 15 maps with different microlens positions in each case (a total of 225 maps). We find that the MPDs have the same form as the central regions of reduced width of our original 24 and 4096pixel maps, and display similar variation with respect to the different microlens positions. Therefore, we can use our original maps across the parameter space to examine the effect of map width and microlens positions on the MPDs.
system  image(s)  model(s)  reference 

QJ01584325  A  f=0.1,0.2,0.3  Morgan et al. (2008) 
MG0414+0534  A1,A2,B  b=1.5  Witt et al. (1995) 
SDSSJ0924+0219  all  f=0.0,0.1,0.2  Morgan et al. (2006) 
HE11041805  A  f=0.1  Morgan et al. (2008) 
PG1115+080  A1,A2,C (B)  b=1.5  Witt et al. (1995) 
RXJ11311231  B,C (A)  f=0.1  Dai et al. (2010) 
HE1413+117  B,C (A)  b=1.5  Witt et al. (1995) 
Q22370305  all  “best fit”  Wambsganss & Paczynski (1994) 
4 Discussion
We produced microlensing magnification maps uniformly covering most of the parameter space of interest, using 15 different sets of random microlens positions at each value. The effect of changing the random microlens positions on the MPD varies across the parameter space and depends on the map width. This was identified in W92, who gave an example of four different map widths for the same and microlens positions (W92, cases 6164 of table 1). However, using our systematic parameter space approach we can estimate the number/width required for magnification maps to show an average behavior.
4.1 Statistical properties of the maps and the mean distribution
Having 15 maps per point makes it possible to examine average statistical properties of the MPDs across the parameter space. We calculate the mean, the mode, the standard deviation, the median and the skewness of the MPDs per point, as shown in Fig. 8. We also calculate the among the 15 cases per value and use it as an average MPD at each point of parameter space.
The mean magnification of the MPDs, in Figs 8(a) and (a’), stays almost always within 5 per cent of the theoretical value (equation 4), in agreement with what was found in W92, which is explained by the finite size of the maps (24 ). However, has a systematic behaviour when compared to in the parameter space: in the minima, in the maxima, while in the saddlepoint region the situation is unclear. This behaviour agrees with the expectation of minima being always magnified (), while saddlepoints and maxima can be demagnified () as well (Schechter & Wambsganss, 2002).
Throughout the parameter space, the appears to be dominated by two peaks and matches those produced by W92 and LI95 for the cases they studied. The highest peak corresponds to the mode of the MPDs [Figs 8(b) and (b’)], the maximum probability, and it is located at as expected. The secondary peak has been identified previously in the works of Rauch et al. (1992), W92 and LI95, and its origin has been attributed to the creation of very bright pairs of microimages at caustics. It is predicted that higher order pairs would appear at even higher magnification, effectively creating additional weaker peaks, which, however, require higher resolution to be seen. Our resolution is sufficient to confirm this by observing a third peak in at least six cases: with and with (Fig. 4).
There seems to be an interplay between the two peaks of the across the parameter space. In the minima region we have a relatively narrow peak located just below , together with a second wider peak that appears at higher magnification. This explains the location of the mode of the MPDs below , and the high values of the standard deviation, the median and the skewness (Fig. 8); positive skewness means assymetric MPDs with tails expanding into high magnifications. As we approach the critical line, the secondary peak approaches and gradually dominates as the primary peak vanishes. Close to the critical line, only one peak is present and located at . The mode is closest to , the standard deviation, the median and the skewness get their minimum values; the MPDs become narrow and symmetric. Crossing into the saddlepoint region, for only one peak remains while for we have two prominent peaks, one in low and one in high magnification, which cause the values of the standard deviation and the median to increase. As we approach again the critical line at higher optical depth, the low magnification peak approaches and gradually dominates, while the high magnification peak vanishes. Reaching the critical line at from the saddlepoint region, there is again only one peak located at and the MPDs are again narrow and symmetric. This peak persists as we go from near the critical line into the maxima region, but becomes wider and reaches into demagnifications the further we go. Using the results described above we can conclude that:

Along the critical line, the MPDs are centered around , extend in a narrow area around it, , and appear symmetric (in log space).

Away from the critical line, the MPDs generally become wider and assymetric, mostly due to the appearance of a high magnification tail and/or a secondary peak.
There seems to be a smooth transition in the shape of the across the parameter space. This could mean that a fit of an analytical formula to the mean MPD could be possible. In other words, an empirical MPD at any given point in parameter space could be constructed. Moreover, posessing a makes it possible to calculate the slope at high , similar to W92, but in all points in parameter space and compare it to theoretical predictions (Schneider, 1987; Witt, 1990). Performing these tasks is beyond the scope of this study and we will deal with them in the future.
A narrow strip of parameter space for is worthy of a special mention. The standard deviation appears significantly lower in this strip than its adjacent neighbours at . This does not happen for the median and the skewness, which get their maximum values in this strip. The explanation is that the MPDs located at are very narrow and have an extended high magnification tail. This tail is dominated by the appearance of a secondary peak in the MPD.
4.2 Interpreting the KS tests
A known systematic of the KS test is related to the choice of the nubmer of bins for the distributions being compared. A large number of bins picks up the fine features of each distribution, leading frequently to a pvalue lower than 5 per cent and the test failing. On the other hand, a small number of bins smooths out the differences among the distributions. Because the binning of the MPDs is arbitrary, we performed the KS test for 3 different numbers of bins, namely 100, 150 and 200. The anticipated behaviour of more MPDs failing the test as the number of bins grows is seen in Figs 10(a), (b) and (c). A number of 100 bins seems to describe the MPD in enough detail for the purposes of accretion disc modelling.
In Fig. 10(a), we show the percentage of all possible MPD pairs at each value that failed the KS test, . is zero or lower than 510 per cent for most of the minima and the saddle point region and increases to around 2030 per cent in the maxima region. In an extreme case, where particular sets of random microlens positions lead to significantly different MPDs, the presence of one such map in our sample of 15 would result in 14 failed KS tests ( per cent) at the given value. In practice, it is possible to have a single MPD for a given value that behaves different to the rest (e.g. magenta line in Fig. 3), but it will fail the test with some and not all of the remaining 14 MPDs. This is what we generally observe in our results, with around 1020 per cent being mostly due to 2 or 3 particular MPDs. However, all the 15 maps, and the corresponding MPDs, represent equally valid choices of random microlens positions, and the fact that some of them appear different to the rest is most likely due to random clustering of microlenses on the image plane. Therefore, in the case of a pair of MPDs failing the KS test, we do not have any criterion yet to distinguish which of the two MPDs should be considered different; they are both equally “correct”. A way out is to perform the KS test between each of the 15 MPDs and the computed for each value. In this case, the is more “correct” than a simple MPD; it represents an average behaviour of 15 MPDs. If the KS test fails in this case, it is because that particular MPD is different to all the remaining 14 that are represented by their . The results of such a test are shown in Figs 10(a’), (b’) and (c’). In the maxima region we have at least 1 out of 15 maps per value, whose MPD is different to the rest. This also occurs on the critical line and for . Increasing the number of bins leads to more MPDs failing the test with the in larger areas of parameter space, similarly to the results of the KS test among pairs of distributions.
The combined effect of random microlens positions and map width can be seen in Fig. 11 throughout the parameter space. For small map widths, 12 , the MPDs significantly vary for different random microlens positions. As the map width is increased, the MPDs are converging. Our maps have a 24 side length, with atypical maps appearing mainly in the maxima area, at a rate of 1 or 2 out of 15. Therefore, to achieve an average behaviour in the maxima area more than one map is required, or alternatively, maps wider than 24 should be generated.
It should be noted here, that the regions of parameter space where there is a higher probability to have an atypical map correlate with regions of the highest number of microlenses, (equation 2); in general, increases with and gets its highest values for combinations, where in equation 4 approaches infinity i.e. along the critical line. The contrary behaviour was anticipated: maps with low numbers of microlenses should be more prone to lens position systematics (e.g. random clustering in the lens plane) while for higher numbers () the purely smooth matter limit should be reached (Garsden et al., 2012). The highest number of microlenses in our maps does not exceed (for and ), therefore, further investigation is warranted to determine the cause for atypical maps appearing in this intermediate region. Unfortunately, this is exactly where the computations are the most demanding.
5 Conclusions
We have performed a uniform systematic study of microlensing parameter space in preparation for future upcoming allsky surveys. We have shown that at 24 map width, a single map should display an average behaviour, except in the maxima region where 2 or 3 maps should be considered. For maps smaller than 24 , the effect of the microlens positions is stronger, leading to a more frequent occurence of atypical maps, even in the minima and saddlepoint regions. For maps wider than 24 , the corresponding MPDs are expected to display an average behaviour independently of the microlens positions.
We have introduced the in Fig. 4 as a representative MPD of each point in the parameter space. The MPD of any magnification map, independently of the method used to generate it, can be compared to the shown in Fig. 4.
It is of great interest to investigate how atypical maps effect derived quasar accretion disc constraints. In general, observations are effected mainly by high magnification events. For example, when convolving a magnification map with a given source profile, the result will depend more on the highly magnified pixels i.e. the high magnification tail of the MPD. Upon constraining accretion disc models, part of the uncertainty of the results could be attributed to the uncertainty of the MPD. As demonstrated, in certain parts of the parameter space the maps are almost statistically equivalent. Therefore, for such values the MPDs are welldefined and an improvement on accretion disc constraints by using more than one map is not anticipated. It would be of interest to examine to what extent an accretion disc model will differ when using MPDs that tend to fail the KS test with the corresponding , however, this is out of the scope of this paper. In Table 2 we identify which of the macromodels of the known lensed systems shown in Fig. 2 fall into regions of space where we have more than 10 per cent of failed pairs, or at least one MPD that failed the test with the (for 24 maps). These are also the systems that will get more effected by microlens positions systematics if a map smaller than 24 is used. In conclusion, in some parts of the parameter space, using several maps with different microlens positions is expected to lead to more robust constraints in accretion disc models. However, the magnitude of this improvement, and consequently the gain of following a manymap approach, is something to be determined in practice. In any case, a map larger than 24 most likely shows an average behaviour.
The results presented above constitute the most complete covering of the convergence and shear parameter space accomplished to date. The different properties of the MPDs across the , parameter space should be considered in designing a microlensing parameter survey, like the one proposed in BF12.
Acknowledgements
The authors would like to thank Nick Bate and Darren Croton. This work was performed on the gSTAR national facility at Swinburne University of Technology. gSTAR is funded by Swinburne and the Australian Government’s Education Investment Fund. Our referee provided very helpful suggestions, which improved the quality of the wrok.
References
 Abramowicz & Fragile (2013) Abramowicz M. A., Fragile C. P., 2013, Living Rev. Relativity, 16, 1
 Bate & Fluke (2012) Bate N. F., Fluke C. J., 2012, ApJ, 744, 90
 Bate et al. (2010) Bate N. F., Fluke C. J., Barsdell B. R., Garsden H., Lewis G. F., 2010, New Astronomy, 15, 726
 Blackburne et al. (2011) Blackburne J. a., Pooley D., Rappaport S., Schechter P. L., 2011, ApJ, 729, 34
 Blaes (2004) Blaes O. M., 2004, preprint (astroph/0211368)
 Blandford & Narayan (1986) Blandford R., Narayan R., 1986, ApJ, 310, 568
 Chang & Refsdal (1979) Chang K., Refsdal S., 1979, Nature, 282, 561
 Dai et al. (2010) Dai X., Kochanek C. S., Chartas G., KozÅowski S., Morgan C. W., Garmire G., Agol E., 2010, ApJ, 709, 278
 Dobler et al. (2007) Dobler G., Keeton C. R., Wambsganss J., 2007, MNRAS, 377, 977
 Eigenbrod et al. (2007) Eigenbrod A., Courbin F., Meylan G., 2007, A&A, 465, 51
 Ferrarese & Ford (2005) Ferrarese L., Ford H., 2005, Space Sci. Rev., 116, 523
 Garsden et al. (2012) Garsden H., Bate N. F., Lewis G. F., 2012, MNRAS, 420, 3574
 Irwin et al. (1989) Irwin M. J., Webster R. L., Hewett P. C., Corrigan R. T., Jedrzejewski R. I., 1989, AJ, 98, 1989
 Kayser et al. (1986) Kayser R., Refsdal S., Stabell R., 1986, A&A, 166, 36
 Keeton (2001) Keeton C. R., 2001, preprint (astroph/0102340)
 Keeton et al. (2006) Keeton C. R., Burles S., Schechter P. L., Wambsganss J., 2006, ApJ, 639, 1
 Kochanek (2004) Kochanek C. S., 2004, ApJ, 605, 58
 Lewis & Irwin (1995) Lewis G., Irwin M., 1995, MNRAS, 276, 103
 LSST Science Collaborations et al. (2009) LSST Science Collaborations et al. 2009
 Mediavilla et al. (2011) Mediavilla E., Mediavilla T., Muñoz J. a., Ariza O., Lopez P., GonzalezMorcillo C., JimenezVicente J., 2011, ApJ, 741, 42
 Mediavilla et al. (2009) Mediavilla E., Muñoz J. a., Falco E., Motta V., Guerras E., Canovas H., Jean C., Oscoz A., Mosquera a. M., 2009, ApJ, 706, 1451
 Mediavilla et al. (2006) Mediavilla E., Munoz J. a., Lopez P., Mediavilla T., Abajas C., GonzalezâMorcillo C., GilâMerino R., 2006, ApJ, 653, 942
 Morgan et al. (2008) Morgan C. W., Eyler M. E., Kochanek C. S., Morgan N. D., Falco E. E., Vuissoz C., Courbin F., Meylan G., 2008, APJ, 676, 80
 Morgan et al. (2006) Morgan C. W., Kochanek C. S., Morgan N. D., Falco E. E., 2006, ApJ, 647, 874
 Morgan et al. (2010) Morgan C. W., Kochanek C. S., Morgan N. D., Falco E. E., 2010, ApJ, 712, 1129
 Mosquera & Kochanek (2011) Mosquera A. M., Kochanek C. S., 2011, ApJ, 738, 96
 Mosquera et al. (2011) Mosquera A. M., Muñoz J. A., Mediavilla E., Kochanek C. S., 2011, ApJ, 728, 145
 Oguri & Marshall (2010) Oguri M., Marshall P. J., 2010, MNRAS, 405, 2579
 Paczynski (1986) Paczynski B., 1986, ApJ, 301, 503
 Peterson (1997) Peterson B., 1997, An introduction to Active Galactic Nuclei. Cambridge Univ. Press, Cambridge
 Rauch & Blandford (1991) Rauch K., Blandford R., 1991, ApJ, 381, L39
 Rauch et al. (1992) Rauch K. P., Mao S., Wambsganss J., Paczynski B., 1992, ApJ, 386, 30
 Schechter & Wambsganss (2002) Schechter P. L., Wambsganss J., 2002, ApJ, 580, 685
 Schechter et al. (2004) Schechter P. L., Wambsganss J., Lewis G. F., 2004, ApJ, 613, 77
 Schmidt & Wambsganss (2010) Schmidt R. W., Wambsganss J., 2010, Gen. Relativ. and Gravit., 42, 2127
 Schneider (1987) Schneider P., 1987, ApJ, 319, 9
 Sluse et al. (2012) Sluse D., Hutsemékers D., Courbin F., Meylan G., Wambsganss J., 2012, A&A, 544, A62
 Thompson et al. (2010) Thompson A. C., Fluke C. J., Barnes D. G., Barsdell B. R., 2010, New Astronomy, 15, 16
 Wambsganss (1990) Wambsganss J., 1990, Phd thesis, Munich Univ.
 Wambsganss (1992) Wambsganss J., 1992, ApJ, 386, 19
 Wambsganss & Paczynski (1994) Wambsganss J., Paczynski B., 1994, AJ, 108, 1156
 Wambsganss et al. (1990) Wambsganss J., Paczynski B., Schneider P., 1990, ApJ, 358, L33
 Witt (1990) Witt H., 1990, A&A, 236, 311
 Witt et al. (1995) Witt H., Mao S., Schechter P., 1995, ApJ, 443, 18
 Wyithe & Turner (2001) Wyithe J. S. B., Turner E. L., 2001, MNRAS, 320, 21
Appendix A Map accuracy and the direct inverse rayshooting technique
Each pixel of a magnification map must have a number of rays large enough to accurately estimate the magnification of the enclosed area of the source plane. In order to achieve this even for the low magnification pixels, we typically have to shoot billions of rays. The total number of rays, , which we are using for each map is determined through:
(6) 
where is an input parameter, same as per equation 6 of Bate et al. (2010). For all our maps we have set , leading to . In this appendix we demonstrate that this is a reasonable trade off between accuracy and computational time.
We have chosen 5 pairs of values across the parameter space, for which we calculate a magnification map with the same microlens positions and set to 500, 1000, 1500 and 2000. The choice of values aims at representing distinct parameter space regions i.e. the minima, saddlepoint, maxima and critical line regions. The corresponding MPDs in each case are shown in Fig. 12 (the axis ranges and units are the same as in Figs 3 and 4). We notice that changing the value of has an effect only towards low magnifications (). For (dotted line) the appearance of dips on the MPD for all values indicates that some magnification values are missing; map pixels corresponding to lowmagnifications, which would have been reached by a small number of rays, were not reached by any because we did not use enough. For (solid line), (dashed line) and (dotteddashed line) the MPDs are essentially indistinguishable above .
Next, we examine the effect of changing the random microlens positions at increased map accuracy. In Fig. 13 we show the MPDs of 15 maps for and , calculated with exactly the same random microlens positions as the 15 maps in Fig. 3, but with twice as many rays per map (and twice the computational time). We see that the microlens position effect persists and only a minor smoothing of the MPDs at low magnification occurs. Because the computational time scales linearly with (equations 7 and 8 in Thompson et al., 2010), we find to be the best choice for our calculations.
Finally, we examine parameter space behaviour of (see Section 2.1). The calculated among all 2550 maps is , however, in Fig. 14 we show its behaviour across the parameter space. We calculate the mean and the standard deviation of the average number of rays per pixel of each map, , among the 15 maps for each value across the parameter space. In Fig. 14(a) we compare the results with a value of 300 (white). We see that Figs 14(a) and (a’) are identical to Figs 8(a) and (a’) where we show the mean magnification, , of the MPDs. This should be anticipated due to a simple relation between and the mean magnification of a map:
(7) 
where is the total number of pixels and we have used equation 3 in the second step. The systematic behaviour seen in distinct parts of the parameter space i.e. overestimating in the minima region and underestimating it in the maxima with respect to , is because the rayshooting area used in the GPUD approach depends on the size and shape of the minimum ellipse (Bate et al., 2010), which in turn depends on the ratio .
We give an estimate of the efficiency of the GPUD approach, with respect to the total number of rays used, , across the parameter space. After we perform the rayshooting and count the rays in each pixel, we are left with a number of rays for which the whole series of calculations was performed, but they were excluded from the final map, , or as a fraction of the total number of rays, . We can write the average number of rays per pixel of the map as:
(8) 
where we used the fact that the number of rays in a map is equal to Replacing from equation 6 and using the calculated , we get or an efficiency of per cent, averaged over the parameter space. If instead of the calculated global we use the individual values across the parameter space, we essentially end up with the same parameter space behaviour of as the one shown in Fig. 14 for ; the GPUD code is more efficient in the minima than in the maxima region.
The low efficiency of the GPUD code directly reflects the brute force approach followed. While it is possible to make use of more sophisticated algorithms (tree code; Wambsganss 1990, inverse polygon mapping; Mediavilla et al. 2006, 2011) to increase the GPUD efficiency, it may prove to be a quite complex task. Instead, we argue that a brute force GPU approach in combination with a GPU cluster like gSTAR can lead more quickly to practical results.