Maximizing Survey Volume for Large-Area Multi-Epoch Surveys with Voronoi Tessellation

TL;DR

This paper introduces a Voronoi tessellation method to dissect large-area surveys into mini-surveys, enabling more accurate volume estimation and recovering about 20% more objects in mock data, especially useful for rare object analysis.

Contribution

The work presents a novel Voronoi tessellation approach to improve survey volume estimation and object recovery in proper motion-limited samples, addressing discrete survey property changes.

Findings

Recovered ~20% more objects than existing methods.

Demonstrated unbiased volume estimation with the new method.

Showed how to increase tessellation resolution with artificial points.

Abstract

The survey volume of a proper motion-limited sample is typically much smaller than a magnitude-limited sample. This is because of the noisy astrometric measurements from detectors that are not dedicated for astrometric missions. In order to apply an empirical completeness correction, existing works limit the survey depth to the shallower parts of the sky that hamper the maximum potential of a survey. The number of epoch of measurement is a discrete quantity that cannot be interpolated across the projected plane of observation, so that the survey properties change in discrete steps across the sky. This work proposes a method to dissect the survey into small parts with Voronoi tessellation using candidate objects as generating points, such that each part defines a `mini-survey' that has its own properties. Coupling with a maximum volume density estimator, the new method is demonstrated to be unbiased and recovered {\sim}20% more objects than the existing method in a mock catalogue of a white dwarf-only solar neighbourhood with Pan--STARRS 1-like characteristics. Towards the end of this work, we demonstrate one way to increase the tessellation resolution with artificial generating points, which would be useful for analysis of rare objects with small number counts.

Figures

Tables

Table 1. Table 1: Physical properties of the Galaxy used in the Monte Carlo simulation.

Parameter	Thin disc	Thick disc	Stellar halo
$⟨U⟩$/km s-1	$-{8.62}^a$	$-{11.0}^d$	$-{26.0}^d$
$⟨V⟩$/km s-1	$-{20.04}^a$	$-{42.0}^d$	$-{199.0}^d$
$⟨W⟩$/km s-1	$-{7.10}^a$	$-{12.0}^d$	$-{12.0}^d$
${\sigma }_U$/km s-1	${32.4}^a$	${50.0}^d$	${141.0}^d$
${\sigma }_V$/km s-1	${23.0}^a$	${56.0}^d$	${106.0}^d$
${\sigma }_W$/km s-1	${18.1}^a$	${34.0}^d$	${94.0}^d$
H/pc	250b	1000e	$\infty$
n/pc-3	0.00310c	0.00064c	0.00019c

Table 2. Table 2: Parameters for the sky background count.

Filter	Mean	Standard Deviation	Lower limit
	(photon s-1)	(photon s-1)	(photon s-1)
g	$39.972$	$12.355$	$16.761$
r	$135.728$	$47.981$	$50.712$
i	$262.725$	$72.499$	$108.612$
z	$257.930$	$90.175$	$92.511$
y	$272.195$	$72.493$	$88.358$

Equations

$$N=\sqrt{(F+d+s)\times t+r^2}$$

$$N_{i,j,k}(D)=\sqrt{\left({\mathcal{F}}_i(D)+d_{j,k}+s_{j,k}\right)\times t_{j,k}+r^2}$$

$${\mathcal{F}}_i(D)=F_i\times {\left(\frac{D_i}{D}\right)}^2.$$

$$\delta {\mathcal{F}}_{i,j,k}(D)=\frac{N_{i,j,k}(D)}{{\mathcal{F}}_i(D)}.$$

$${\sigma }_{i,j,k}(D)=\sqrt{\delta {\mathcal{F}}_{i,j,k}^2(D)+{\sigma }_s^2},$$

$$\underbrace{\left(\begin{array}[]{cc}\dfrac{1}{\sigma_{1}}&\dfrac{\Delta t_{1}}{\sigma_{1}}\\ \cdot&\cdot\\ \cdot&\cdot\\ \dfrac{1}{\sigma_{M(j)}}&\dfrac{\Delta t_{M(j)}}{\sigma_{M(j)}}\end{array}\right)}_{\mathbf{A}}\times\left(\begin{array}[]{c}\overline{\alpha}\\ \mu_{\alpha}\end{array}\right)=\left(\begin{array}[]{c}\dfrac{\Delta\alpha_{1}}{\sigma_{1}}\\ \cdot\\ \cdot\\ \dfrac{\Delta\alpha_{M(j)}}{\sigma_{M(j)}}\end{array}\right)$$

$$\mathbf{\mathsf{A}}^{\mathrm{T}}\mathbf{\mathsf{A}}=\left[\begin{array}[]{cc}\sum\limits_{k}\left(\dfrac{1}{\sigma_{k}}\right)^{2}&\sum\limits_{k}\left(\dfrac{1}{\sigma_{k}}\dfrac{\Delta t_{k}}{\sigma_{k}}\right)\vspace*{1mm}\\ \sum\limits_{k}\left(\dfrac{1}{\sigma_{k}}\dfrac{\Delta t_{k}}{\sigma_{k}}\right)&\sum\limits_{k}\left(\dfrac{\Delta t_{k}}{\sigma_{k}}\right)^{2}\end{array}\right]$$

$$\frac{1}{{\sigma }_{{\mu }_{\alpha }\cos \delta }^2}=\sum _k{\left(\frac{\Delta t_k}{{\sigma }_k}\right)}^2$$

$${\sigma }_{\mu }=\sqrt{{\sigma }_{{\mu }_{\alpha \cos \delta }}^2+{\sigma }_{{\mu }_{\delta }}^2}=\sqrt{2}{\sigma }_{{\mu }_{\alpha \cos \delta }}.$$

$$\begin{gathered} {\mathrm{V}}_{\max }=\sum _j{\Omega }_j{\int }_{D_{\min ,\mathrm{j}}}^{D_{\max ,\mathrm{j}}}\frac{\rho (D)}{{\rho }_{\odot }}\times D^2 \\ \text{hfill}\times \left[{\int }_{a(D)}^{b(D)}{P}_{l(j)}(v_{\mathrm{T}})\mathrm{d}{v}_{\mathrm{T}}\right]\mathrm{d}D \end{gathered}$$

$$a(D)=\max \left[v_{\min },4.74\times s\times {\sigma }_{\mu }(D)\times D\right]$$

$$\frac{\rho (D)}{{\rho }_{\odot }}=\exp \left(-\frac{\left|z\right|}{H}\right)=\exp \left(-\frac{\left|D\sin b\right|}{H}\right),$$

$$\Phi (M_{\mathrm{bol}})={\int }_{{\mathcal{M}}_l}^{{\mathcal{M}}_u}\frac{\mathrm{d}{t}_{\mathrm{cool}}}{\mathrm{d}{M}_{\mathrm{bol}}}\psi \left(t_0-t_{\mathrm{cool}}-t_{\mathrm{MS}}\right)\varphi \left(\mathcal{M}\right)\mathrm{d}\mathcal{M},$$

$$t_0=t_{\mathrm{cool}}\left(M_{\mathrm{bol}},\xi ({\mathcal{M}}_l)\right)+t_{\mathrm{MS}}\left({\mathcal{M}}_l,Z\right)$$

$$\xi \left({\mathcal{M}}_i\right)={\mathcal{M}}_f\left({\mathcal{M}}_i\right)=0.101{\mathcal{M}}_i+0.463{\mathcal{M}}_{\odot }.$$

Full text

††pagerange: Maximising Survey Volume for Large Area Multi-Epoch Surveys with Voronoi Tessellation–Maximising Survey Volume for Large Area Multi-Epoch Surveys with Voronoi Tessellation††pubyear: 2017

1 Introduction

A number of types of transient, variable and moving sources are not rare, but their detection requires repeated observations of the same part of sky. This was not possible to perform over a large sky area until the era of digital astronomy. The highly automated observing runs and efficient digital detectors allow efficient data collection, while faster processors and the automated data reduction pipelines allow the production of high volume of output. Earliest attempts for such automations digitise the photographic plates from large sky area surveys, where measurements are made objectively with computer, as opposed to measuring manually. Large scale projects of this kind include the PPM catalogue (Röser & Bastian, 1991; Bastian & Röser, 1993), Automated Plate Machine Project (Evans, 1992; Evans & Irwin, 1995), USNO A 1.0, A 2.0 and B 1.0 (Monet, 1996, 1998; Monet et al., 2003), SuperCOSMOS (Hambly et al., 2001a, b; Hambly et al., 2001c), UCAC 1, 2, 3 and 4 (Zacharias et al., 2000, 2004, 2010, 2013) and SUPERBLINK (Lépine et al., 2002, 2003; Lépine, 2005, 2008), all of which had several epochs and simple tiling strategies.

In the current era of digital astronomy, some surveys continue to use simple tiling patterns where multiple pawprints are combined immediately to produce full coverage over a sky cell, for example in the UKIDSS (Lawrence et al., 2007), four pawprints can cover a cell while VISTA employs six (Sutherland et al., 2015). In other cases, SDSS Stripe 82 had nine epochs on average (Bramich et al., 2008), ALLWISE has a coverage from 12 to over 200 frames (Kirkpatrick et al., 2014, 2016); and the Pan-STARRS1 (PS1) 3$\pi$ Steradian Survey (3SS) typically has 60 epochs (Magnier et al., 2013), The Dark Energy Survey will scan $\sim$5,000 deg2 10 times (The Dark Energy Survey Collaboration, 2005), Gaia will have on average 81 transits, with over 140 in the most-visited parts of the sky at the end of the 5 yr nominal mission (Gaia Collaboration et al., 2016), and LSST will provide close to 1,000 epochs for half of the sky towards the end of the 10 yr survey mission (LSST Science Collaboration et al., 2009)111https://github.com/LSSTScienceCollaborations/. The key survey characteristics (depth and epoch coverage) vary on small scales and in complex ways because tiling strategies/overlapping patterns for these surveys are extremely complicated to maximize coverage due to losses from, for example, chip gaps between CCDs and unfavourable observing conditions, unlike the situation with large-format photographic plates. This in turn complicates the analysis of any survey sample culled from them because optimal techniques like $V_{\mathrm{max}}$ require the precise survey characteristics. This problem becomes even more complex when different surveys are combined to expand the wavelength coverage and/or maximum epoch difference, for example when SDSS is combined with USNO-B 1.0 to derive the proper motions (Gould & Kollmeier, 2004; Munn et al., 2004), the survey typically has five epochs (four from USNO-B and one from SDSS).

Existing methods tackle inhomogeneity by limiting the analysis to the shallower part of the survey and by applying a global correction in order not to run into unaccountable incompleteness (e.g. Harris et al., 2006; Munn et al., 2017, hereafter M17). In view of this problem, when deriving the white dwarf luminosity function (WDLF) with a maximum volume density estimator, Rowell & Hambly (2011, hereafter RH11) measured the empirical photometric and astrometric uncertainty for different skycells as defined by the tiling of the photographic plates. Lam et al. (2015, hereafter LRH15) improved the completeness correction due to kinematic selections. The methods described in these works allow an analysis to probe deeper, but in order to maximize the use of the data, we propose a new method based on Voronoi tessellation that can further maximize the analytical survey volume within which completeness and other biases can be corrected. Voronoi tessellation has been employed in defining simulation grids, clustering analysis, visualization etc. However, its property that partitions sources into well defined grids has not been used in all domains of astrophysics. By dividing the sky through Voronoi tessellation into a number of cells that is equal to the number of candidate sources, we can treat each cell as a mini-survey that has well-defined local properties. Analysis can also be performed at lower or higher resolution if needed.

In Section 2, we will describe the mathematical framework and the construction of the simulated solar neighbourhood in Section 3. The new method is applied to the simulated data in Section 4 under different selection criteria and we describe one procedure through which the resolution can be increased. The bias due to the choice of model is briefly discussed. In the final section, we discuss one possible way to increase the resolution for the analysis and conclude this work.

2 Mathematical Framework of the Voronoi Method

The maximum volume density estimator (Schmidt, 1968) tests the observability of a source by finding the maximum volume in which it can be observed by a survey (e.g. at a different part of the sky at a different distance). It is proven to be unbiased (Felten, 1976) and easily can combine multiple surveys (Avni & Bahcall, 1980). In a sample of proper motion sources, we need to consider both the photometric and astrometric properties (see LHR15 for details). The number density is found by summing the number of sources weighted by the inverse of the maximum volumes. For surveys with small variations in quality from field to field and from epoch to epoch, or with small survey footprint areas, the survey limits can be defined easily. However, in modern surveys, the variations are not small; this is especially true for ground-based observations. Therefore, properties have to be found locally to analyse the data most accurately. Through the use of Voronoi tessellation, sources can be partitioned into individual 2D cells within which we assume the sky properties are defined by the governing source. Each of these cells has a different area depending on the projected density of the population.

An important assumption for using the Voronoi method is that the distributions of the observing parameters of the cells at different resolutions are very similar to each other, hence the integrated maximum volume is approximately equal to the exact solution. In the rest of the article, cell will be used to denote Voronoi cell and h-pixel for HEALPix pixel (Górski et al., 2005, see Section 2.3) and pixels for the ones on a detector.

2.1 Voronoi Tessellation

A Voronoi tessellation is made by partitioning a plane with $n$ points into $n$ convex polygons such that each polygon contains one point. Any position in a given polygon (cell) is closer to its generating point than to any other for the case of Voronoi tessellation using Euclidean distance $\left(D_{E}=\sqrt{(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}}\right)$. For use in astronomy, such a tessellation has to be done on a spherical surface (two-sphere).

In the following work, the tessellation is constructed with the SciPy package spatial.SphericalVoronoi, where each polygon is given a unique ID that is combined with the vertices to form a dictionary. The areas are calculated by first decomposing the polygons into spherical triangles with the generating points and their vertices222https://github.com/tylerjereddy/py_sphere_Voronoi and then by using L’Huilier’s Theorem to find the spherical excess. For a unit-sphere, the spherical excess is equal to the solid angle of the triangle. The sum of the constituent spherical triangles provides the solid angle of each cell.

2.2 Cell Properties

For a Voronoi cell $j$, the properties of the cell are assumed to be represented by generating source $i$. Both $i$ and $j$ are indexed from $1$ to $\mathcal{N}$, but since each source has to be tested for observability in each cell to calculate the maximum volume, $i$ and $j$ cannot be contracted to a single index. Furthermore, the cells do not need to be defined by only the sources. Arbitrary points can be used for tessellation such that $i$ and $j$ will not have a one-to-one mapping. The epoch of the measurement is labelled by $k$. When tested for observability, epoch-wise information is essential in calculating the photometric and proper motion uncertainties as functions of distance. The major difference in the following approach is that the proper motion uncertainty is found from the formal propagation of errors instead of measuring the empirical form as a function of magnitude, $\sigma_{\mu}($mag$)$, which limits the survey to the worst part of a tile (RH11). This new approach does not need to take into account the scatter in $\sigma_{\mu}($mag$)$, due to different local sky properties and different colours of the sources and their neighbours. Different types of source can differ by up to a few magnitudes in the optical/infrared colours, so two sources with similar magnitudes in one filter can have very different proper motion uncertainties if one is close to the detection limit in another filter. The modelling of the photometric uncertainties from CCD detectors is much simpler than that for photographic plates, because the photometric response of the modern detectors is much more linear at both the faint and bright ends. Thus, the uncertainties can be estimated with relatively simple equations.

2.2.1 Photometric Uncertainty

When a source is being tested for the observability, it is ‘placed’ at a different distance so the apparent brightness changes as a consequence. The background and other instrumental noises are constant, but the Poisson noise from the source changes with the measured flux, hence the photometric uncertainties are functions of distance. The total noise, $N$, of a photometric measurement can be estimated by

[TABLE]

where $F$ is the instrumental flux per unit time, $d$ is the dark current per unit time, $s$ is the sky background flux per unit time, $t$ is the exposure time and $r$ is the read noise. Among these quantities, $d$, $s$, $t$ and $r$ are fixed quantities in a given epoch, only the flux varies as a function of distance. We use $F$ as the measured flux and $\mathcal{F}(D)$ as the flux at an arbitrary distance $D$. Therefore, in a Voronoi cell $j$ at epoch $k$, the photometric noise of source $i$ is

[TABLE]

where the flux at $D$ is calculated from applying the inverse square law on the observed flux $F_{i}$ and observed distance $D_{i}$,

[TABLE]

The random photometric uncertainty of a source at an arbitrary distance in a given epoch is the inverse signal-to-noise ratio,

[TABLE]

The total photometric uncertainty of the source as a function of distance, combining with the systematic uncertainty, $\sigma_{s}$, coming from the absolute calibration of the detector, is therefore

[TABLE]

which represents the photometric uncertainty as a function of the distance to the source.

2.2.2 Astrometric Uncertainty

The least square solution of proper motion in one direction for source $i$ can be expressed in the following matrix form, the epoch is labelled by the subscript from $1$ to $M(j)$ where $M$ is the number of epochs in cell $j$,

[TABLE]

where $\Delta t_{k}$ is the time difference between the mean epoch and epoch $k$, $\Delta\alpha_{k}$ is the positional offset from the mean position, $\overline{\alpha}$, and proper motion, $\mu_{\alpha}$, in the direction of the right ascension. The associated uncertainties can be found from the diagonal terms of the normal matrix,

[TABLE]

so for each cell,

[TABLE]

and the total proper motion uncertainty is

[TABLE]

The uncertainties in the $\alpha$ and $\delta$ directions are symmetrical in four-parameter astrometric solution (two positions and two proper motions). In the case of five-parameter solution where parallax is solved and for the seven-parameter solution where, in addition, the acceleration terms in both directions are solved for, the uncertainties will not be symmetrical due to the parallactic term, so the off-diagonal terms have to be taken into account which would otherwise be negligible compared to the diagonal terms. However, for a variance-weighted mean epoch, the off-diagonal terms are exactly zero by definition.

2.3 Consequence to $V_{\mathrm{max}}$ Calculation

There is only one minor adjustment to the volume integral – the lower proper motion limit. Instead of finding the limit by measuring from a number of nearby sources that include mostly sources with different colours, the limit is defined by the properties of the Voronoi cell that comes only from the generating source of the cell. A different set of pixelization by HEALPix, denoted by $l$, is for the line-of-sight tangential velocity completeness correction (see RH11 for detailed description). For source $i$, the maximum volume has to be tested in each cell $j$, the expression is almost identical to that in LRH15, except for the $j$ and $l$ indexes

[TABLE]

where $\dfrac{\rho(D)}{\rho_{\odot}}$ is the density normalized by that at the solar neighbourhood, $P_{l(j)}$ is the tangential velocity distribution, $l(j)$ denotes the h-pixel mapped from cell $j$ with area $\Omega_{j}$, $v_{\mathrm{T}}$ is the tangential velocity, $D_{\mathrm{min}}$ and $D_{\mathrm{max}}$ are the minimum and maximum photometric distances, and $\sigma_{\mu}(D)$ is the proper motion uncertainty as a function of the distance to the source. This is to model the change in the proper motion uncertainties due to varying apparent magnitude with distance (i.e. at greater distance the proper motion uncertainty will be larger because the source becomes fainter which increases the single-epoch positional uncertainty). The lower tangential velocity limit in the inner integral, $a(D)$, is

[TABLE]

where the factor of 4.74 comes from the unit conversion from arcsec yr*-1* to km s*-1* at distance $D$, $v_{\mathrm{min}}$ is the global lower tangential velocity limit and $s$ is the significance of the proper motions. The expression is identical to that in LRH15, but $\sigma_{\mu}(D)$ is calculated in a completely different way.

The inner integral can vanish before reaching the distance limits so the integrator must use a small step size or the distances at which the inner integral vanish have to be calculated explicitly. The cell ID $j$ and h-pixel ID $l$ can be set as a one-to-one mapping by calculating the tangential velocity distribution for each of the Voronoi cell. However, the Voronoi tessellation is dependent on the sample, while the tangential velocity correction is fixed on the sky, using a precomputed look-up table because the tangential velocity correction can significantly reduce computation time.

3 Simulated Data Set

To demonstrate the power of the Voronoi method described in Section 2, we apply it to catalogues of simulations of the solar neighbourhood. This section details the construction of the Monte Carlo simulation.

We generated snapshots of white dwarf (WD)-only populations in the solar neighbourhood containing six dimensional phase space information. The procedure is very similar to that described in LRH15, however, we introduce changes to the noise model of the system and include epoch-wise information. The volume probed is assumed to be small such that the simulation is performed in a Cartesian space, instead of a plane polar system centred at the Galactic Centre. The Galaxy has three distinct kinematic components: a thin disc, a thick disc and a stellar halo, all of which we model with no density variations along the co-planar directions of the Galactic plane. The vertical structures of the discs follow exponential profiles, with scale height $H_{\mathrm{thin}}$ and $H_{\mathrm{thick}}$ such that

[TABLE]

where $z$ is the vertical distance from the Galactic plane and $b$ the Galactic latitude. None of the three components are tilted relative to each other. The velocity components, $U$, $V$ and $W$, of each WD are drawn from the Gaussian distributions described by the measured means and standard deviations of the three sets of kinematics that describe the three populations in the solar neighbourhood (Table 1).

Theoretical WDLFs are used as the probability distribution functions (pdfs) in the simulations. The normalizations of the pdfs are taken from the WD densities found in RH11. The input parameters for a WDLF are the star formation rate (SFR), initial mass function (IMF), MS evolution model and WD cooling model. The standard equation for modelling the WDLF with those four given inputs is

[TABLE]

where $\Phi(M_{\mathrm{bol}})$ is the number density of WDs at magnitude $M_{\mathrm{bol}}$. The derivative inside the integral is the characteristic cooling time of WDs, $\psi(t)$ is the SFR at time $t$ and $\phi$ is the IMF. The input parameters are assumed to be invariant with time and are summarized in Table 1. The integral also depends on the lifetimes of MS progenitors, $t_{\mathrm{MS}}$, as a function of mass and metallicity. We have adopted the stellar evolution tracks from the Padova group (PARSEC; Bressan et al. 2012) with a metallicity of $Z=0.019$ and $Y=0.30$ (Girardi et al., 2000). By assuming a fixed surface gravity $\log{g}=8.0$ and pure hydrogen atmosphere (DA)333http://www.astro.umontreal.ca/$\sim$bergeron/CoolingModels/, the WD cooling time, $t_{\mathrm{cool}}$, is a function of mass and luminosity (Holberg & Bergeron, 2006; Kowalski & Saumon, 2006; Bergeron et al., 2011; Tremblay et al., 2011), and $t_{0}$ is the total time since the onset of star formation. The integral is over all MS masses that have had time to produce WDs at the present day, with the magnitude-dependent lower limit, $\mathcal{M}_{l}$, corresponding to the solution of

[TABLE]

and the upper limit for WD production $\mathcal{M}_{u}\approx 8\mathcal{M}_{\odot}$. The initial–final mass relation, $\xi$, relates the MS progenitor mass to the mass of the WD, is adopted from Kalirai et al. (2009) without including globular clusters in the analysis where the final WD mass can be expressed as

[TABLE]

The thin and thick disc populations are assigned with constant SFR since look back time, $\tau=8$ Gyr and $\tau=10$ Gyr respectively, while the halo has a starburst of $1$ Gyr at $\tau=13$ Gyr. The IMF has an exponent of $-2.3$ in the mass range of interest (Kroupa, 2001).

From the pdfs of the kinematics, distance and bolometric magnitude, we calculate the apparent magnitudes in the PS1 gP1, rP1, iP1, zP1 and yP1 filters444http://panstarrs.stsci.edu/ (Tonry et al., 2012; Schlafly et al., 2012; Magnier et al., 2013; Chambers et al., 2016; Waters et al., 2016; Magnier et al., 2016a; Flewelling et al., 2016; Magnier et al., 2016b, c). The line-of-sight and the projected velocity can be derived from the given 3D kinematics and 3D position. The uncertainties in the five filters are calculated from the sky background flux, exposure time, dark current and read noise that are representative of the PS1 3SS at Processing Version 2 (PV2). The sky background flux is drawn from a Gaussian distribution measured from the 3SS in each of the filters (Table 2). The means and standard deviations555$1.4826$ times the median absolute deviation is used for robust estimation of the standard deviation. were measured from $100$ fields drawn randomly across the survey footprint area.

To simulate the variations in the observing properties of the sky, HEALPix is used to pixelate the sky using a resolution of $n_{\mathrm{side}}=256$, i.e. each has a size of $\sim 0.0534$ deg2, which is sufficiently small compared to the projected density of white dwarfs that is less than $1$ deg*-2* (e.g. RH11). The HEALPix resolution is on a much finer scale than the Voronoi tessellation used for determining the maximum volume in order to test the accuracy later (Section 4.4). One feature of this approach is that when the analysis is done at a higher resolution, the new set of cells will be guaranteed to land on a different set of h-pixels and as such will provide a self-consistency check. There is not a switching from Voronoi cells to HEALPix pixels in the analysis: there is simply a look-up table for matching a given cell with the h-pixel that the cell-generating source lands on. Each h-pixel is given a sky background noise for each epoch of the measurement. When the sky brightness is below the lower limit, a sky brightness that is fainter than $99.73\%$ of the measured values, the background noise is resampled until it is above the limit. An ADU of $1$ e*-* per photon is assumed. The 3SS has $12$ epochs on average in each of the filters, so the number of epochs for each source in the simulation is drawn from a distribution666This study only focuses on tessellation; the effect of non-detection is another huge step in the optimization of analysis. that follows $1+P(11)$, where $P(11)$ is Poisson distribution with a mean of $11$ and the epochs are drawn from a random distribution over a period of $3$ yr with $6$ h on either side of the source masked out to simulate seasonal observing. When sources are distributed over the sky, they will take the set of values defined by the nearest pixel. Using the treatments from Section 2.2.1, with a dark current of $0.2\,$e*-s-1*, exposure times of $43,40,35,30$ and $30$ s, zero-point magnitudes at $24.563$, $24.750$, $24.611$, $24.250$ and $23.320$ mag in the five filters, and a constant read noise of $5.5\,$e*-* (Metcalfe et al., 2013), each source is assigned with proper motion uncertainty using equation (9). These inputs produce an all-sky survey that has $10\sigma$ detections (and their standard deviations) in gP1, rP1, iP1, zP1 and yP1 at $21.98\pm 0.04$, $21.53\pm 0.05$, $21.12\pm 0.04$, $20.54\pm 0.05$ and $19.59\pm 0.04\mathop{}\!\,\mathrm{mag}$. They are similar to the PV2 values, but the distribution is much narrower because the noise model is itself noiseless (e.g. the sources are not affected by diffraction spikes, optical ghosts, cosmic rays or other effects that lead to larger scatter).

4 Application to WDLFs

This section describes the application of our survey dissection method (presented in Section 2) to simulated PS1-like WD catalogues generated using the recipe described in Section 3. The bright limits at all filters are set at $15\mathop{}\!\,\mathrm{mag}$. The faint limits are at $21.5$, $21.0$, $20.5$, $20.0$ and $19.5\mathop{}\!\,\mathrm{mag}$ in gP1, rP1, iP1, zP1 and yP1 filters respectively, which are the typical magnitudes at which the 3SS is complete. The lower proper motion limit is set to five times the proper motion uncertainties, $\sigma_{\mu}$, unless specified otherwise; and the upper proper motion limit is set at $0.08438$ and $0.4219$ arcsec yr*-1* for the cases using lower tangential velocity limits at $40$ and $200\mathop{}\!\,\mathrm{km\,s}^{-1}$, respectively. The two limits correspond to a minimum distance of $100$ pc; this is to avoid any bias coming from the unaccounted parallax signature of very nearby sources. The upper tangential velocity limits are different in each analysis. The photometric parallaxes were not derived, and the real distances and bolometric magnitudes are used. The volume and the maximum volume are found by integrating equation (10) from $D_{\mathrm{min}}$ to $D$, and from $D_{\mathrm{min}}$ to $D_{\mathrm{max}}$ respectively.

4.1 Comparison with the RH11 selection

The RH11 method increases the survey volume by restricting shallow survey depths only over areas that are severely limited by a small number of poor observations. In this section we illustrate how the Voronoi method can further increase the number of sources that can be recovered while rigorous completeness correction can be performed. In Fig. 1, under the selection criteria: $40\mathop{}\!\,\mathrm{km\,s}^{-1}<v_{\mathrm{tan}}<60\mathop{}\!\,\mathrm{km\,s}^{-1}$, proper motion less than half an arcsec year*-1* and a minimum distance of $100$ pc, the number of sources recovered by the Voronoi method is plotted as a solid line, RH11 method as a dashed line and the global $95$th percentile as a dotted line. The ratio between the RH11 and Voronoi methods is plotted as a thick black solid line, while that between a global lower proper motion limit and the Voronoi method is plotted as a thick grey solid line. With the Voronoi method, more sources can be recovered. The ratio between the RH11 and Voronoi methods slowly decreases as the absolute bolometric magnitudes increase; the ratio with the global limit simply plummets given that the upper proper motion limit is only $0.5$ arcsec yr*-1*. Fig. 2 shows that the ratios stay fairly constant until some faint limits. The sources lost with the RH11 treatment are due to a combination of (1) the reassignment of proper motion uncertainties based on empirical observations where $95$-per-cent of all sources are given larger uncertainty values than their measured ones and (2) the loss of the deepest areas of the survey.

Due to the simplistic noise model of the simulation (i.e. no bad pixels, saturation, diffraction, optical ghosts and other effects that affect the photometric and astrometric precision significantly), the distribution of the proper motion uncertainties in the simulation is typically narrower than real measurements. Nevertheless, at $5\sigma$ level the Voronoi method can recover $\sim 15\%$ more sources than the RH11 method and much more sources than applying a global lower proper motion limit (Fig. 2).

4.2 Thin Disc and Combined Discs

Study of the thin disc WDLF requires a selection of the low-velocity population in order to minimize the contaminations from older populations, which typically possess higher velocities. In this section, we show the observed WDLFs from a thin disc-only simulation and from a mixed thin disc, thick disc and halo simulation. The WDLF comparison plots are displayed with the WDLF in the top panel, differences between the input and calculated WDLFs and the $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle$ as a function of bolometric magnitude are in the middle panel and at the bottom panel respectively.

4.2.1 Thin disc-only sample

In an analysis selecting only thin-disc WDs, the observed WDLF agrees very closely with the input function down to $M_{\mathrm{bol}}\approx 15\mathop{}\!\,\mathrm{mag}$ when the number of sources drops significantly (Fig. 3). From the $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle(\mathrm{mag})$ distribution, the derived solution is very stable throughout, except at the brightest and the faintest ends where the $\mathop{}\!\,\frac{1}{\mathrm{V}_{\mathrm{max}}}$ method is known to become less reliable as the number of sources decreases. In theory, $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle=0.5$, because it is the expected value of a uniform distribution between [math] and $1$. Statistically, it is expected that only $\sim 60\%$ of the time the $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle$ lies within the error bar. The uncertainty in $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle$ is $\frac{1}{\sqrt{12N}}$. The small oscillation about the line at $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle(\mathrm{mag})=0.5$ is a good indication that the sample is unbiased over a large dynamic range of magnitudes. The outliers at the extreme ends result from the application of the density estimator to a small number of sources, and so likely do not represent the true values. Taking $40$ and $60\mathop{}\!\,\mathrm{km\,s}^{-1}$ as the lower and upper tangential velocity limits of the inner integral (equation 10 & 11 and the equivalent set of the upper limit), the total integrated number density of the work is $3.65\pm 1.03\times 10^{-3}$ pc*-3*, compared to the input $3.10\times 10^{-3}$ pc*-3*. The overall $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle=0.4944\pm 0.0031$, which is very close to $0.5$, indicating an unbiased sample.

4.2.2 Mixed Population ($40-60\mathop{}\!\,\mathrm{km\,s}^{-1}$)

The modification to the density estimator itself is small, only the lower limit of the inner integral is changed (equation 11); and it is not expected that the effect due to contamination should differ from the previous analysis in LRH15. The extra depth enabled by the new method could have led to a significant increase in the measured density due to a combination of two effects: (1) an increase in contamination fraction as the thin disc contribution drops rapidly with distance: at the Galactic poles, the thin disc and thick disc densities equate at $525$ pc; and (2) the kinematic completeness correction applied on contaminants, which are more common at fainter magnitudes.

The kinematics of the two discs are well measured; however, the relative density of WDs in them is much less studied – there is only one measurement on record (RH11). To understand the effects of contamination, a better understanding of the two populations is needed. Nevertheless, we can compare the WDLFs from the last section to a mixed population with the same set of upper and lower tangential velocity limits ($40$ and $60\mathop{}\!\,\mathrm{km\,s}^{-1}$; Fig. 4). The total integrated number density is $4.00\pm 1.03\times 10^{-3}$ pc*-3* as compared to $3.10\times 10^{-3}$ pc*-3* for the thin disc and $0.64\times 10^{-3}$ pc*-3* for the thick disc, which sum to $3.74\times 10^{-3}$ pc*-3*. If it is treated as a pure thin disc WDLF, there is a roughly constant overestimation of $0.1$ dex at all magnitudes. When both discs are considered, the small over-density ($0.26\times 10^{-3}$ pc*-3*) is due to contamination from the thick disc where the using of a thin disc scale-height on these sources will lead to an overestimation of the maximum volume. In this simulation, $16.0\%$ of the data are from the thick disc.

The $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle$ distribution is very similar to the clean sample, and for the entire sample $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle=0.4984\pm 0.0028$ which is within 1 standard deviation of the ideal value $0.5$. We believe this velocity range is a good choice for driving an upper limit of the thin disc white dwarf density in the solar neighbourhood.

4.2.3 Mixed Population ($40-120\mathop{}\!\,\mathrm{km\,s}^{-1}$)

In M17, $40$ and $120\mathop{}\!\,\mathrm{km\,s}^{-1}$ are used as the tangential velocity limits to study both discs together, with a scale-height of $300$ pc instead of $250$ pc. From H06, it is known that the effect of scale-height is larger at the bright end because sources can be seen at a larger distance hence the density-correction is larger in equation (10). The effect on the total normalization is small because faint sources dominate after density and completeness corrections. However, studies in, for example, star formation history (Rowell, 2013) or high energy exotic particles (Isern et al., 2008) are sensitive to the whole range of magnitudes. This cannot simply be assumed to be negligible. Fig. 5 investigates this effect by comparing the cases of $250$ and $300$ pc scaleheights. It shows that the choice of scale-height has almost no effect to the WDLFs except at the brightest magnitudes, and for the distribution of $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle$, the differences are negligible. However, the absolute normalizations from using the two scale-heights are consistently overestimated by $\sim 0.1$ dex. In this simulation, $23.1\%$ of the sources in the range of $40-120\mathop{}\!\,\mathrm{km\,s}^{-1}$ are from the thick disc and the halo; in comparison, only $15.5\%$ of sources are not from the thin disc in the $40-60\mathop{}\!\,\mathrm{km\,s}^{-1}$ selection.

4.3 Halo

The study of the halo WDLF requires a selection of high velocity population in order to minimise the contaminations from the thick disc. In this section, we will show the observed WDLFs from a halo-only simulation and from a mixed thin disc, thick disc and halo simulation.

4.3.1 Halo-only sample

In the halo-only simulation, the observed WDLF agrees very well with the input function (Fig. 6). From the $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle(\mathrm{mag})$ distribution, the derived solution is stable throughout, except at the faintest bin where there are only two sources. The small oscillation about the line at $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle=0.5$ is a good indication that the sample is unbiased over the entire range of magnitudes. The upper and lower tangential velocity limits are set at $200$ and $500\mathop{}\!\,\mathrm{km\,s}^{-1}$ which define the survey limits of the inner integral (equations 10 & 11 and the equivalent set of the upper limit), the total integrated number density of the work is $1.77\pm 0.10\times 10^{-4}$ pc*-3*, compared to the input $1.90\times 10^{-4}$ pc*-3* (the integrated density up to $15\mathop{}\!\,\mathrm{mag}$ is $1.68\times 10^{-4}$ pc*-3*) and $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle=0.5116\pm 0.0200$.

4.3.2 Mixed Population ($200-500\mathop{}\!\,\mathrm{km\,s}^{-1}$)

In 10 realizations of the halo-only simulation, under the selection of $200-500\mathop{}\!\,\mathrm{km\,s}^{-1}$, there is a mean contamination rate of $7.3\%$ (minimum at $3.6\%$ and maximum at $11.9\%$). In the thin disc analysis, where the contamination is over $15\%$, we do not observe any significant bias in the WDLF or $\mathop{}\!\,\left\langle\mathop{}\!\,\mathrm{V}/\mathrm{V}_{\mathrm{max}}\right\rangle$(mag) distribution. We believe these fractions of contaminations have little effect on the analysis so the samples should be representative of the halo, we therefore do not consider them further.

4.4 Sensitivity to resolution

The purpose of this new method is to tackle a complex survey that has small-scale variations that reduce the maximum survey volume that is available to a study. The way this method divides the sky has avoided the detailed treatment of the survey and approximates it with the source properties that go into the analysis, which raises a concern whether the resolution for a small sample would be sufficient, and not causing significant systematic bias. We suggest a method to increase the resolution by a factor of $\sim 3$, which can be repeated to further increase the resolution if needed, as it would be useful if there are only a few sources over the sky. A higher resolution can be achieved by using the vertices as new generating points (Fig. 7). However, an increased resolution requires much more computation time, because the time complexity of the Voronoi method is $\mathrm{O}\left(\mathcal{N}^{2}\right)$. There is a trade-off between the accuracy and the computing time. The properties at the new cells can be approximated by those carried by the nearest sources (or they can be extracted directly from the field from the raw data). Fig. 8 demonstrates with $10$ halo analyses that for a very well-behaved survey (small differences in survey depths), a resolution in the order of $\sim 30$ deg2 per cell compared to $\sim 100$ deg2 per cell will only lead to an increase of $<1\%$ in number density (top panel). The increase arises from the deeper parts of the survey that the standard resolution always underestimates (as it is statistically less likely to land on the deeper parts). To understand the effect at lower resolutions, we simulate this by using one in three sources to generate the Voronoi tessellation. The change in number density is only in the range of a few percentage points (bottom panel). Over a large number of simulations, the ratios should average to $1$. The asymmetric distribution comes from the inverse proportionality between the maximum volume and the number density.

5 Conclusion

In this work, we have demonstrated that the use of Voronoi tessellation can increase the survey volume to more optimally retrieve sources from a large sky area multi-epoch survey. The assumption behind this is not ideal, but it is not possible to take an average value of the number of epochs and their properties for each cell. Further subdivisions will take much more time to compute the volumes as this algorithm scales as $x\times\mathcal{N}^{2}$, where $x$ is the number of subdivision of each cell and $\mathcal{N}$ is the number of sources. Nevertheless, this method is one big step in approaching the optimal sampling of the survey footprint with limited computing power.

From a mixed population simulation, we find that under the framework of our galactic models, the new method recovers $\sim 10-15\%$ more sources than the RH11 method under a typical lower proper motion selection. When considering a restricted tangential velocity selection ($40-60\mathop{}\!\,\mathrm{km\,s}^{-1}$), we do not observe any bias in the WDLFs brighter than $15\mathop{}\!\,\mathrm{mag}$. Similar result is observed from the $40-120\mathop{}\!\,\mathrm{km\,s}^{-1}$ sample. However, this conclusion is valid only up to $\mathrm{M}_{\mathrm{bol}}=15\mathop{}\!\,\mathrm{mag}$, the result should not be extrapolated to fainter magnitudes where thin disc contribution to the number density drops significantly compared to the thick disc and the halo. This work has deliberately removed sources closer than $100$ pc to avoid bias caused by parallactic displacements, which as a consequence removed all the faint sources that can be seen only from a small distance. In the high velocity regime, at the given thick disc-to-halo density ratio, a $200-500\mathop{}\!\,\mathrm{km\,s}^{-1}$ selection will only contain a small fraction of contaminations so it is a good sample for studying the halo.

We have demonstrated one way to increase the resolution of the tessellation and it shows that for a well behaved survey, a low resolution only limits the volume by $\sim 1\%$. An adaptive way that only subdivides cells larger than a certain solid angle can provide a grid of cells that have similar areas should it be more useful in some certain scenarios. When applying this method to real surveys, careful treatment at the boundaries is needed because the area of the survey is important. Leaving the boundaries untreated one will always end up with $4\pi$ steradian of sky area. In order to have a correct boundary that defines the survey, artificial points have to be added to create a layer of bounding cells surrounding the survey area such that the boundaries of the second last layer of cells overlap the survey footprint. One can identify the artificial points by using the survey boundary as a cell boundary, and then locate a generating point that can produce the correct boundary geometry. However, one should note that a typical survey boundary is given by a small circle on the celestial sphere, while the Voronoi tessellation constructs cell boundaries along the great circles. The total area described by the Voronoi cells is always going to be slightly different from the true survey area (unless it is a full sky survey). However, the difference in area is very small, comparable to the unaccounted non-perfect fill-factor or dead pixels at the detector.

In future surveys, the complexity of the tiling pattern, scanning strategy as well as the detector arrays at the focal plane will only increase. There is an increasing need for a more optimal analytic tool to maximally use the available data. The Voronoi method can include the faintest objects that would otherwise be neglected because of unaccountable incompleteness.

Acknowledgments

The Pan–STARRS 1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan–STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE).

We thank the PS1 Builders and PS1 operations staff for construction and operation of the PS1 system and access to the data products provided. ML acknowledges financial support from the STFC Consolidated Grant of the Institute for Astronomy, University of Edinburgh. ML also wishes to thank Dr. Nigel Hambly and the referee Dr. Floor van Leeuwen for helpful comments and suggestions that have led to major improvements in the clarity and presentation of the article.