© ESO, 2017

1. Introduction

The development of new detectors with high throughput over large areas has revolutionized observational astronomy during recent decades. These technological advances, aided by a considerable increase of computing power, have opened the way to outstanding ground-based and space-borne all-sky or very large area imaging projects (e.g. the 2MASS Skrutskie et al. 2006; Cutri et al. 2003; SDSS Ahn et al. 2012, 2013; and WISE Wright et al. 2010; Cutri et al. 2014, surveys). These surveys have provided an essential astrometric and photometric reference frame and the first true digital maps of the entire sky.

As an illustration of this flood of data, the number of catalogue entries in the VizieR service at the Centre de Données astronomiques de Strasbourg (CDS) which was about 500 million in 1999 has reached almost 18 billion as on February 2016. At the 2020 horizon, European space missions such as Gaia and Euclid together with the Large Synoptic Survey Telescope (LSST) will provide a several-fold increase in the number of catalogued optical objects while providing measurements of exquisite astrometric and photometric quality.

This exponentially increasing flow of high quality multi-wavelength data has radically altered the way astronomers now design observing strategies and tackle scientific issues. The former paradigm, mostly focusing on a single wavelength range, has in many cases evolved towards a systematic fully multi-wavelength study. In fact, modelling the spectral energy distributions over the widest range of frequencies, spanning from radio to the highest energy gamma-rays has been instrumental in understanding the physics of stars and galaxies.

Many well designed and useful tools have been developed worldwide concurrently with the emergence of the virtual observatory. Most if not all of these tools can handle and process multi-band images and catalogues. When assembling spectral energy distributions using surveys obtained at very different wavelengths and with discrepant spatial resolution, one of the most acute problems is to find the correct counterpart across the various bands. Several tools such as TOPCAT (Taylor 2005) or the CDS XMatch Service (Pineau et al. 2011a; Boch et al. 2012) offer basic cross-matching facilities. However, none of the publicly available tools handles the statistics inherent to the cross-matching process in a fully coherent manner. A standard method for a dependable and robust association of a physical source to instances of it in different catalogues (cross-identification) and in diverse spectral ranges is still absent.

The pressing need for a multi-catalogue probabilistic cross-matching tool was one of the strong motivations of the FP7-Space European program ARCHES (Motch et al. 2016)¹. Designing a cross-matching tool able to process, in a single pass, a theoretically unlimited number of catalogues, while computing probabilities of associations for all catalogue configurations, using the background of sources, positional errors and eventually introducing priors on the expected shape of the spectral energy distribution is one of the most important outcomes of the project. A preliminary description of this algorithm was presented in Pineau et al. (2015). Although ARCHES was originally focusing on the cross-matching of XMM-Newton sources, the algorithms developed in this context are clearly applicable to any combination of catalogues and energy bands (see for example Mingo et al. 2016).

2. Going beyond the two-catalogue case

Computing probabilities of identifications when cross-correlating two catalogues in a given area can be quite straightforward (provided the area is small enough so that the density of sources can be considered more or less constant, but large enough to provide sufficient statistics). For each possible pair of sources (one from each catalogue), we compute the distance normalized by positional errors (called normalized distance, σ-distance, χ-distance or more generally in this paper Mahalanobis distance D_M). Then we build the histogram of the number of associations per bin of D_M. This histogram is the sum of two components (see Fig. 2): the “real” or “true” associations (T) for which the distribution p(D_M | T) follows a Rayleigh distribution; the spurious or “false” associations, for which the distribution p(D_M | F) follows a linear (Poisson) distribution. Knowing these two distributions and the total number of associations (n_{T + F}), we may fit the histogram with the function $\mathit{f}\mathrm{\left(}{\mathit{D}}_{\mathrm{M}}\mathrm{\right)}\mathrm{=}\mathrm{(}{\mathit{n}}_{\mathit{T}\mathrm{+}\mathit{F}}\mathrm{-}{\mathit{n}}_{\mathit{F}}\mathrm{)}\mathit{p}\mathrm{\left(}{\mathit{D}}_{\mathrm{M}}\mathrm{\right|}\mathit{T}\mathrm{)}\mathrm{+}{\mathit{n}}_{\mathit{F}}\mathit{p}\mathrm{(}{\mathit{D}}_{\mathrm{M}}\mathrm{\left|}\mathit{F}\mathrm{\right)}$(1)to estimate the number of spurious associations (n_F) and thus the number of good matches (n_T = n_{T + F}−n_F). Hence, we are able to attribute to an association with a given normalized distance the probability of being a good match: $\mathit{p}\mathrm{\left(}\mathit{T}\mathrm{\right|}{\mathit{D}}_{\mathrm{M}}\mathrm{)}\mathrm{=}\frac{\mathrm{(}{\mathit{n}}_{\mathit{T}\mathrm{+}\mathit{F}}\mathrm{-}{\mathit{n}}_{\mathit{F}}\mathrm{)}\mathit{p}\mathrm{\left(}{\mathit{D}}_{\mathrm{M}}\mathrm{\right|}\mathit{T}\mathrm{)}}{\mathit{f}\mathrm{\left(}{\mathit{D}}_{\mathrm{M}}\mathrm{\right)}}\mathrm{·}$(2)Dividing both the numerator and the denominator by n_{T + F} we recognize Bayes’ formula considering (n_{T + F}−n_F) /n_{T + F} as the prior p(T) and considering either n_F/n_{T + F} as the prior p(F) or f(D_M) = n_{T + F}p(D_M).

The present paper basically extends this simple approach to more than two catalogues. Instead of fitting histograms to find the number of spurious associations, we directly compute them from the input catalogues data and from geometrical considerations.

Previously, Budavári & Szalay (2008) developed a multi-catalogue cross-match. For a given set of n sources from n distinct catalogues, they compute a “Bayes’ factor” based on both astrometric and photometric data. The “Bayes’ factor” is then used as a score: a pre-defined threshold on its value is applied to select or reject the given set of n sources. We discuss the astrometric part of Budavári & Szalay (2008) “Bayes’ factor” and compare it to our selection criterion in Sects. 5.6 and 6.1.

Throughout the present paper we consider a set of n catalogues. We use a Chi-square criterion based on individual elliptical positional errors to select, in these catalogues, sets of associations containing at most one source per catalogue. We call this selection a χ-match. We then compute probabilities for each set of associations. To compute probabilities, we consider only the result set in which each set of associations contains exactly n sources (one per catalogue, see below for partial matches). For people familiar with databases, it can be seen as the result of inner joins, joining successively each catalogue using a Chi-square criteria. The probabilities we then compute are only based on positional coincidences. Although we show how it is possible to add likelihoods based on photometric considerations, the computation of such photometric likelihoods is beyond the scope of this paper.

As the result of a χ-match, two distinct sets of associations may have sources in common: a source having a large positional error in one catalogue may for example be associated to several sources with smaller errors in another catalogue. We do not take into account in our probabilities the “one-to-several” and the “one-to-one” associations paradigms defined in Fioc (2014): it becomes far too complex when dealing with a generic number of catalogues and it is not that simple when a source may be blended, etc. We use a several-to-several-(to-several-...) paradigm. In other words, we compute probabilities for a set of associations regardless of the fact that a source in the set can be in other sets of associations. So a same detection in one catalogue may have very high probabilities of associations with several (sets of) candidates in the other catalogues. We think it is the responsibility of the photometric part to disentangle such cases.

Requiring one candidate per catalogue for each set of associations (i.e. each tuple) is somewhat restrictive. But, if one or several catalogues do not contain any candidates for a tuple, then we compute the probabilities from the cross-match of the subset of catalogues providing one candidate to that tuple. Those probabilities are computed independently of the “full” n catalogues probabilities. For example, if we cross-match three catalogues and if a set of associations (a tuple) contains one source per catalogues (A, B and C), then we will compute five probabilities: one for each possible configuration (ABC, AB_C, A_BC, AC_B and A_B_C in which the underscore “_” separates the catalogue entries associated to different actual sources, see Sect. 6.2.2). Now, if one source from A has a candidate in B and no candidate in C, we will compute only two probabilities (AB and A_B, see Sect. 6.2.1) considering only the result of the cross-match of A with B. Likewise for A and C only and for B and C only. These four cross-matches will yield eleven distinct probabilities. It is possible to deal with “missing” detections when computing photometrically based likelihoods (taking into account limit fluxes, ...) but it is not the case in the astrometric part of this work.

When χ-matching n catalogues, the number of hypotheses to be tested, and thus the number of probabilities to be computed for a given set of associations, increases dramatically with n. This number is 203 for 6 catalogues and reaches 877 for seven catalogues (see Table 2 in Sect. 6.2.4). To be able to compute probabilities when χ-matching more than seven catalogues we may start by merging catalogues for which the probability of making spurious associations is very low (e.g. catalogues of similar wavelength and similar astrometric accuracy), and handle the merged catalogue as a single input catalogue.

In Sect. 3 we lay down the assumptions we use to work on a simplified problem. We then (Sect. 4) define the notations and the standards used throughout the paper and link them to the standards adopted in a few catalogues. We then describe in detail the candidate selection criterion (Sect. 5) before providing (Sect. 6) an exhaustive list of all hypotheses we have to account for to apply Bayes’ formula. In Sect. 5 we also show how the “Bayesian cross-match” of Budavári & Szalay (2008) may be interpreted as an inhomogeneous χ-match. Then (Sect. 7) we show how it is possible to estimate the rates of spurious associations and hence “priors”. In Sect. 8 we compute an integral which is related to the probability the selection criterion has to select a set of n sources for a given hypothesis. This integral is crucial to compute likelihoods defined in Sect. 9 and to normalize likelihoods in Sect. 10. Finally, after showing how to introduce the photometric data into the probabilities (Sect. 11), and before concluding (Sect. 14), we explain the tests we carried out on synthetic catalogues in Sect. 12. Since this paper is long and technical, we put a summary of the steps to follow to perform a probabilistic χ-match in Sect. 13.

3. Simplifying assumptions

Cross-correlating catalogues taking into account an accurate model of the sky on one hand, and the effects and biases due to the catalogue building process on the other hand is a daunting task. To make progress towards this objective, we have to start by making simplifying assumptions.

First of all, we assume that there are no systematic offsets between the positions of each possible pair of catalogues. It means that the positions are accurate (no bias). We also assume that positional errors provided in catalogues are trustworthy. It means that they are neither overestimated nor underestimated: for instance, no systematic have to be quadratically added or removed. The first point supposes an accurate astrometric calibration of all catalogues. This is somewhat the “dog chasing its tail” problem since a proper astrometric calibration should be based on secure identifications, themselves based on... cross-identification! Ideally the astrometric calibration and the cross-identification should be performed simultaneously in an iterative process. It will not be developed here but we point out that the present work can be used to calibrate astrometrically n catalogues at the same time from one reference catalogue, taking into account all possible associations in all possible catalogue sub-sets. However, carrying out careful identification of primary or secondary astrometric standards is only important when the density of bright astrometric references is very low, typically in deep small field exposures. Reliable cross-identification is also crucial when the wavelength band of the image to calibrate differs widely from that of the astrometric reference image. In most large scale surveys such as 2MASS (Skrutskie et al. 2006) or SDSS (Pier et al. 2003) the density of bright Tycho-2 (Høg et al. 2000) or UCAC (Zacharias et al. 2004) astrometric reference stars is high enough to ensure an excellent overall calibration without any ambiguity in the associations.

Although the idealized vision of an immutable and static sky is long gone, we ignore proper motions in this analysis. There are at least two ways of taking them into account: either we may force associations to include at least one source from a catalogue containing measured proper motions; or we may try to fit proper motions during the cross-match process. In this last case, if a set of n sources detected at different epochs in n distinct catalogues does not satisfy the candidate selection defined in Sect. 5.2, we may make the hypothesis that they nonetheless are from a same source but having a proper motion. We can then estimate the proper motion and the associated error based on positions, (Gaussian) positional errors and epochs (see Appendix B). From the n observed positions and associated errors and from the n theoretical estimated positions and associated errors we can compute a Mahalanobis distance which follows a χ distribution with 2(n−2) degrees of freedom. Similarly to the candidate selection criterion in Sect. 5.2 we can then reject the hypothesis “same source with proper motion” if the Mahalanobis distance is larger than a given threshold.

We neglect clustering effects. We suppose that in a given area Ω, source properties are homogeneous. This implies that the local density of sources, the positional error distributions and the associations priors (probabilities of true associations that in principle depend on the astrophysical nature of the sources and on the limiting flux) are uniform over the sky area considered. As usual we have to face the following dilemma: on the one hand, the larger the area Ω, the better the statistic; on the other hand, the larger the area Ω, the less probable the uniform density, errors distributions and priors hypothesis. In the ARCHES project, for instance, we grouped the individual XMM-Newton EPIC fields of view of ≈0.126 deg² each into installments of homogeneous exposure times and galactic latitude so as to ensure as much uniformity as possible. Each installment contained on the order of several hundred sources.

Finally, we neglect blending. If two sources are separated in one catalogue and blended in the other one, the position of the blended source will be something like the photocentre of the two sources. Either the blended source will not match any of the two distinct sources, or only one of the two distinct sources will match, the match likely being in the tail of the Rayleigh distribution, possibly leading to a low probability of identification. It will then not be problematic to consider the match as spurious since the observed flux is contaminated by the flux of the nearby source. Finally, if the positional accuracy of the blended source is well below that of the distinct sources, both distinct sources will match the blended source, leading to a non-unique association requiring further investigations to be disentangled.

4. Notations and links with catalogues

4.1. Notations

This article uses almost exclusively the notations defined in the ISO 80000-2:2009(E) international standard. Exceptionally we waive the notation detA for determinant and replace it by the equivalent but more compact notation | A |.

We consider n catalogues defined on a common surface of area Ω. We assume that each catalogue source has individual elliptical positional errors defined by a bivariate normal (or binormal) distribution. For this, we assimilate locally the surface of the sphere to its zenithal (or azimuthal) equidistant projection (see ARC projection in Calabretta & Greisen 2002), that is to its local Euclidean tangent plane. In this frame, the position of a point at distance d arcsec from the origin O (the tangent point) and having a position angle ϕ (east of north) is simply $\begin{array}{ccc}\mathit{x}& \mathrm{=}& \mathit{d}\sin \mathit{\varphi ,}\\ \mathit{y}& \mathrm{=}& \mathit{d}\cos \mathit{\varphi }\mathit{.}\end{array}$This approximation is acceptable since typical positional errors, distances and surfaces locally considered are small.

We note the binormal probability density function (p.d.f.) representing the position of a source S and its associated uncertainty: ${\mathrm{𝒩}}_{{\mu }\mathit{,}{V}}\mathrm{\left(}{p}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\sqrt{\det {V}}}\exp \left\{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{Q}\mathrm{\left(}{p}\mathrm{\right)}\right\}\mathrm{d}{p}\mathit{,}$(5)with

• μ = (μ_x,μ_y)^⊺ the position of the source S provided in a catalogue, that is the mean of the binormal distribution;

• V the provided variance-covariance – also simply called covariance – matrix which defines the error on the source position;

• p = (x,y)^⊺ any given two-dimensional position;

• Q(p) the quadratic form Q(p) = (p−μ)^⊺V^-1(p−μ), that is the square of the weighted distance between a given position p and the position of the source S$\begin{array}{ccc}{V}& \mathrm{=}& \left(\begin{array}{c}\\ & \\ & \end{array}\right)\mathit{,}\\ {{V}}^{-1}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{x}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{y}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}^{\mathrm{2}}\mathrm{)}}\left(\begin{array}{c}\\ & \\ & \end{array}\right)\mathit{,}\end{array}$where

  • σ_x is the standard deviation along the x-axis (i.e. the east axis);

  • σ_y is the standard deviation along the y-axis (i.e. the north axis);

  • ρ the correlation factor between σ_x and σ_y;

• $\det {V}\mathrm{=}{\mathit{\sigma }}_{\mathit{x}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{y}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}^{\mathrm{2}}\mathrm{)}$ the determinant of V;

• dp = dxdy.

A covariance matrix V represents a 1σ ellipse. The “real” position of the source S has ≈39% chances to be located inside this 1σ-ellipse. It must not be confused with the 1-dimensional 1σ-segment which contains a real “value” with a probability of ≈68%.

4.2. Classical positional errors in catalogues

In astronomical catalogues like the 2MASS All-Sky Catalog of Point Sources (2MASS-PSC, Cutri et al. 2003) positional errors are described by three parameters defining the 1σ positional uncertainty ellipse²: err_maj or a the semi-major axis, err_min or b the semi-minor axis, and err_ang or ψ the positional angle (east of north) of the semi-major axis. We give the formula to transform the ellipse into a covariance matrix (see Appendix A.2 of Pineau et al. 2011b and footnote 11 in Fioc 2014): $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{x}}& \mathrm{=}& \sqrt{{\mathit{a}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\psi }\mathrm{+}{\mathit{b}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\psi }}\mathit{,}\\ {\mathit{\sigma }}_{\mathit{y}}& \mathrm{=}& \sqrt{{\mathit{a}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\psi }\mathrm{+}{\mathit{b}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\psi }}\mathit{,}\\ \mathit{\rho }{\mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\mathit{y}}& \mathrm{=}& \cos \mathit{\psi }\sin \mathit{\psi }\mathrm{(}{\mathit{a}}^{\mathrm{2}}\mathrm{-}{\mathit{b}}^{\mathrm{2}}\mathrm{)}\mathit{.}\end{array}$In the AllWISE catalogue (Cutri et al. 2014), the coefficients of the covariance matrix are (almost) directly available. Instead of providing the unitless correlation factor ρ or the covariance ρσ_xσ_y (in arcsec²), the authors chose to provide the co-sigma (σ_αδ) because, as they state³, the latter is in the same units as the other uncertainties. We thus have $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{x}}& \mathrm{=}& {\mathit{\sigma }}_{\mathit{\alpha }}\mathrm{=}\mathit{sigra,}\\ {\mathit{\sigma }}_{\mathit{y}}& \mathrm{=}& {\mathit{\sigma }}_{\mathit{\delta }}\mathrm{=}\mathit{sigdec,}\\ \mathit{\rho }{\mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\mathit{y}}& \mathrm{=}& {\mathit{\sigma }}_{\mathit{\alpha \delta }}\mathrm{\times }\mathrm{\left|}{\mathit{\sigma }}_{\mathit{\alpha \delta }}\mathrm{\right|}\mathrm{=}\mathit{sigradec}\mathrm{\times }\mathrm{\left|}\mathit{sigradec}\mathrm{\right|}\mathit{.}\end{array}$In catalogues like the Sloan Digital Sky Survey since its eighth data release (SDSS-DR8, Aihara et al. 2011) positional errors contain two terms: the error on RA (raErr) and the error on Dec (decErr). In this case the parameters of the covariance matrix are simply $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{x}}& \mathrm{=}& \mathit{raErr,}\\ {\mathit{\sigma }}_{\mathit{y}}& \mathrm{=}& \mathit{decErr,}\\ \mathit{\rho }{\mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\mathit{y}}& \mathrm{=}& \mathrm{0.}\end{array}$In catalogues like the XMM catalogues (e.g. the 3XMM-DR5, Rosen et al. 2016) a single error is provided. Ideally, one would like to have access to the two one-dimensional errors, even if their respective values are often very close. The column named radecErr is the total error, so the quadratic sum of the two computed (but not provided) 1-dimensional errors, one computed on RA and one computed on Dec. If one uses σ_x = radecErr and σ_y = radecErr, the total error will be $\mathit{\sigma }\mathrm{=}\sqrt{{\mathit{\sigma }}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathit{y}}^{\mathrm{2}}}\mathrm{=}\sqrt{\mathrm{2}}\mathit{radecErr}$ instead of radecErr. In output of the astrometric calibration process, the XMM pipeline provides a systematic error sysErrCC which is quadratically added to radecErr to compute the “total radial position uncertainty”⁴posErr. As for radecErr, we must divide posErr by $\sqrt{\mathrm{2}}$ to obtain the 1-dimensional error. The appropriate errors to be used (including a systematic) are then $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{x}}& \mathrm{=}& {\mathit{\sigma }}_{\mathit{y}}\mathrm{=}\mathit{posErr}\mathit{/}\sqrt{\mathrm{2}}\mathit{,}\\ \mathit{\rho }{\mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\mathit{y}}& \mathrm{=}& \mathrm{0.}\end{array}$The factor $\sqrt{\mathrm{2}}$ has not been taken into account in Pineau et al. (2011b). It partly explains why the fit of the curve in the right panel of Fig. 3 mentioned in Sect. 5 of this paper does not lead to a Rayleigh scale parameter equal to 1.

Similarly to the XMM case, the error posErr provided in the GALEX All-Sky Survey Source Catalog (GASC) catalogue⁵ (which also includes the systematic) is a “total radial error”. It is thus the Rayleigh parameter σ which is the quadratic sum of two one-dimensional errors. As for XMM, the appropriate errors to be used are $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{x}}& \mathrm{=}& {\mathit{\sigma }}_{\mathit{y}}\mathrm{=}\mathit{posErr}\mathit{/}\sqrt{\mathrm{2}}\mathit{,}\\ \mathit{\rho }{\mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\mathit{y}}& \mathrm{=}& \mathrm{0.}\end{array}$In catalogues like the ROSAT All-Sky Bright Source Catalogue (1RXS, Voges et al. 1999) the error provided is the radius of the cone containing the real position of a source with a probability of ≈68.269% (the 1 dimensional 1σ). Authors like Rutledge et al. (2003; given the details provided in Voges et al. 1999, Sect. 3.3.3) call this radius the 1σ-radius. We note it r_68%. But, in the Rayleigh distribution, the scale parameter σ is defined such that the cone of radius r = σ contains the real position with a probability 100 × (1−exp(−1/2)) ≈ 39.347%. Adjusting such that $\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{\times }{\mathit{r}}_{\mathrm{68}\mathrm{\%}}^{\mathrm{2}}\mathit{/}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{)}\mathrm{=}\mathrm{0.6827}$ leads to $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{x}}& \mathrm{=}& {\mathit{\sigma }}_{\mathit{y}}\mathrm{=}\mathit{\sigma }\mathrm{=}\frac{{\mathit{r}}_{\mathrm{68}\mathrm{\%}}}{\sqrt{\mathrm{-}\mathrm{2}\ln \mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0.6827}\mathrm{)}}}\\ & \mathrm{=}& \frac{{\mathit{r}}_{\mathrm{68}\mathrm{\%}}}{\sqrt{\mathrm{2}\ln \mathrm{\left(}\mathrm{3.1515}\mathrm{\right)}}}\mathrm{\approx }\frac{{\mathit{r}}_{\mathrm{68}\mathrm{\%}}}{\mathrm{1.51517}}\mathrm{·}\end{array}$(21)Similarly if the provided error is the radius of the cone containing the real position with a probability of 90% (e.g. in the WGACAT, White et al. 2000) ${\mathit{\sigma }}_{\mathit{x}}\mathrm{=}{\mathit{\sigma }}_{\mathit{y}}\mathrm{=}\mathit{\sigma }\mathrm{=}\frac{{\mathit{r}}_{\mathrm{90}\mathrm{\%}}}{\sqrt{\mathrm{-}\mathrm{2}\ln \mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0.90}\mathrm{)}}}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{90}\mathrm{\%}}}{\sqrt{\mathrm{2}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}}}\mathrm{\approx }\frac{{\mathit{r}}_{\mathrm{90}\mathrm{\%}}}{\mathrm{2.14597}}\mathrm{·}$(22)The description (White et al. 1997) and the on-line documentation⁶ of the FIRST catalogue (Helfand et al. 2015a,b) provide an “empirical expression” to compute the semi-major and semi-minor axis of the 90% positional accuracy associated to each source: $\begin{array}{ccc}{\mathit{a}}_{\mathrm{90}\mathrm{\%}}& \mathrm{=}& \mathrm{fMaj}\left(\frac{\mathrm{RMS}}{\left(\mathrm{Fpeak}\mathrm{-}\mathrm{0.25}\right)}\mathrm{+}\frac{\mathrm{1}}{\mathrm{20}}\right)\mathit{,}\\ {\mathit{b}}_{\mathrm{90}\mathrm{\%}}& \mathrm{=}& \mathrm{fMin}\left(\frac{\mathrm{RMS}}{\mathrm{(}\mathrm{Fpeak}\mathrm{-}\mathrm{0.25}\mathrm{)}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{20}}\right)\mathit{,}\end{array}$in which fMaj (fMin) is the major (minor) axis of the fitted FWHM, RMS “is a local noise estimate at the source position” and Fpeak is the peak flux density. The position angle ψ of the accuracy equals the fitted FWHM angle fPA. We first obtain the 1σ accuracy ellipse by resizing the 90% ellipse axes divinding them by the same factor as for the WGACAT (i.e. $\sqrt{\mathrm{2}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}}$) $\begin{array}{ccc}\mathit{a}& \mathrm{=}& \frac{{\mathit{a}}_{\mathrm{90}\mathrm{\%}}}{\sqrt{\mathrm{2}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}}}\mathit{,}\\ \mathit{b}& \mathrm{=}& \frac{{\mathit{b}}_{\mathrm{90}\mathrm{\%}}}{\sqrt{\mathrm{2}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}}}\mathrm{·}\end{array}$After possibly adding systematics, the variance-covariance matrix is obtained applying the equations used for the 2MASS catalogue.

Errors in catalogues like the Guide Star Catalog Version 2.3.2⁷ (Lasker et al. 2007, GSC2.3) should not be used in the framework of this paper. As stated in Table 3 of Lasker et al. (2008): these astrometric and photometric errors are not formal statistical uncertainties but a raw and conservative estimate to be used for telescope operations.

Table 1 summarizes the transformation of catalogues positional errors into the coefficients of covariance matrices V.

5. Candidates selection: the χ-match

We make the hypothesis that n sources from n distinct catalogues are n independent detections of a same real source. With p the unknown position of the real source and μ_i the observed position of detection i, the probability for the n detections to be located at the observed positions is expressed by the joint density function: $\begin{array}{ccc}{\mathit{f}}_{\mathit{p}}\mathrm{\left(}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathit{,}\mathit{...}\mathit{,}{{\mu }}_{{n}}\mathrm{\right|}{p}\mathrm{)}& \mathrm{=}& \underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽙}}}{\mathrm{𝒩}}_{{{\mu }}_{{i}}\mathit{,}{{V}}_{{i}}}\mathrm{\left(}{p}\mathrm{\right)}\mathit{,}\\ & \mathrm{=}& \end{array}$(27)

5.1. Estimation of the real position given n observations

We introduce the notations μ_Σ and V_Σ for the weighted mean position of the n sources and its associated error respectively. The inverse of the covariance matrix V_Σ is ${{V}}_{{\Sigma }}^{-1}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{V}}_{{i}}^{-1}\mathit{,}$(28)leading to (see demonstration in Sect. A.1) ${{V}}_{{\Sigma }}\mathrm{=}\frac{\mathrm{1}}{\det {{V}}_{{\Sigma }}^{-1}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{{V}}_{{i}}}{\det {{V}}_{{i}}}$(29)which is used in the weighted mean position expression ${{\mu }}_{{\Sigma }}\mathrm{=}{{V}}_{{\Sigma }}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{V}}_{{i}}^{-1}{{\mu }}_{{i}}\mathit{.}$(30)Using both the weighted mean position and its error, the sum of quadratics in Eq. (27) can be divided into two parts and written as (see demonstration Sect. A.2)$\begin{array}{ccc}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{Q}}_{\mathit{i}}\mathrm{\left(}{p}\mathrm{\right)}& \mathrm{=}& \\ & \mathrm{=}& {\mathit{Q}}_{\mathit{p}}\mathrm{\left(}{p}\mathrm{;}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathit{,}\mathit{...}\mathit{,}{{\mu }}_{{n}}\mathrm{\right)}\mathrm{+}{\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}\mathrm{\left(}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathit{,}\mathit{...}\mathit{,}{{\mu }}_{{n}}\mathrm{\right)}\mathit{,}\end{array}$with $\begin{array}{ccc}{\mathit{Q}}_{\mathit{p}}\mathrm{\left(}{p}\mathrm{;}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathit{,}\mathit{...}\mathit{,}{{\mu }}_{{n}}\mathrm{\right)}& \mathrm{=}& {\mathrm{(}{p}\mathrm{-}{{\mu }}_{{\Sigma }}\mathrm{)}}^{⊺}{{V}}_{{\Sigma }}^{-1}\mathrm{(}{p}\mathrm{-}{{\mu }}_{{\Sigma }}\mathrm{)}\mathit{,}\\ {\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}\mathrm{\left(}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathit{,}\mathit{...}\mathit{,}{{\mu }}_{{n}}\mathrm{\right)}& \mathrm{=}& \end{array}$In the case of two catalogues the latter term can be written as in Eq. (51). Moreover, if both covariances are null, it takes the simple and common form ${\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}\mathrm{=}\frac{\mathrm{\Delta }{\mathit{\alpha }}^{\mathrm{2}}}{{\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathrm{1}}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathrm{2}}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{2}}}{{\mathit{\sigma }}_{{\mathit{\delta }}_{\mathrm{1}}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{{\mathit{\delta }}_{\mathrm{2}}}^{\mathrm{2}}}\mathrm{·}$(35)Back to the general case, the term Q_χ² can also be put in the more computationally efficient form (only one loop over i) $\begin{array}{ccc}{\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}\mathrm{\left(}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathit{,}\mathit{...}\mathit{,}{{\mu }}_{{n}}\mathrm{\right)}& \mathrm{=}& \end{array}$(36)From those formulae, it appears that the weighted mean position (Eq. (30)) is the maximum likelihood estimator of the “true” position of the source: the second term (Q_χ²) is constant with respect to p so the maximum of the likelihood function ℒ(p;μ₁,μ₂,...,μ_n) = f_p(μ₁,μ₂,...,μ_n | p) is obtained when the first term (Q_p) is null, so when p = μ_Σ. The error on this estimate is simply V_Σ, the inverse of the Hessian of the likelihood function.

5.2. Candidates selection criterion

For the candidate selection, we are interested in the probability the n sources have to be located at the same position. Let’s first rewrite Eq. (27) to exhibit a product of a binormal distribution by another multi-dimensional normal law: $\begin{array}{ccc}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽙}}}{\mathrm{𝒩}}_{{{\mu }}_{{i}}\mathit{,}{{V}}_{{i}}}\mathrm{\left(}{p}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\sqrt{\det {{V}}_{{\Sigma }}}}\exp \left\{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{Q}}_{\mathit{p}}\right\}\\ & & \end{array}$(37)When integrating Eq. (37) over all possible positions (i.e. over p) the first term integrates to 1, since it is the p.d.f. of a normal law in p, so we obtain $\mathrm{\int }\mathrm{\int }\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽙}}}{\mathrm{𝒩}}_{{{\mu }}_{{i}}\mathit{,}{{V}}_{{i}}}\mathrm{\left(}{p}\mathrm{\right)}\mathrm{d}{p}\mathrm{=}\sqrt{\frac{\det {{V}}_{{\Sigma }}}{\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽑}}}\det {{V}}_{{i}}}}\frac{\exp \mathrm{\{}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}\mathrm{\}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{1}}}\mathrm{·}$(38)We are supposed to integrate on the surface of the unit sphere. But the errors being small, we consider the infinity being at a relatively close distance, before effects of the sphere curvature become non-negligible.

In the previous equation, only the Q_χ² term remains. It can also be written (see demonstration Sect. A.3) ${\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}\mathrm{\left(}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathit{,}\mathit{...}\mathit{,}{{\mu }}_{{n}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\sum _{\mathit{j}\mathrm{=}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{n}}{\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{{\mu }}_{{j}}\mathrm{)}}^{⊺}{{V}}_{{i}}^{-1}{{V}}_{{\Sigma }}{{V}}_{{j}}^{-1}\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{{\mu }}_{{j}}\mathrm{)}\mathit{.}$(39)Equation (38) is equivalent to P(D | H) in Budavári & Szalay (2008) and Eq. (39) – multiplied by the $\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}$ factor in the exponential (Eq. (38)) – is the generalization for elliptical errors of Eq. (B12) in Budavári & Szalay (2008). In practice, we never use Eq. (38) since the number of terms to be computed increases with O(n(n−1)/2) while it increases with O(n) in Eq. (36) or in its iterative form (see Eq. (48) in Sect. 5.3). We use here the big O notation, to be read as “the order of”.

We can see Eq. (34) as the result of a 2n-dimensional weighted least squares in which the model is the “real” position of the source and the solution is μ_Σ (by similarity with Eq. (31)). Putting all positional errors matrices in a 2n × 2n block diagonal matrix M, Q_χ² is the square of the Mahalanobis distance ${\mathit{D}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{\left(}{\mu }\mathrm{\right)}$ defined by ${\mathit{D}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{\left(}{\mu }\mathrm{\right)}\mathrm{=}{\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}\mathrm{\left(}{\mu }\mathrm{\right)}\mathrm{=}{{v}}^{⊺}{{M}}^{-1}{v}\mathit{,}$(40)${v}\mathrm{=}\left(\begin{array}{c}\\ {{\mu }}_{{1}}\mathrm{-}{{\mu }}_{{\Sigma }}\\ {{\mu }}_{{2}}\mathrm{-}{{\mu }}_{{\Sigma }}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {{\mu }}_{{n}}\mathrm{-}{{\mu }}_{{\Sigma }}\end{array}\right)\mathit{,}{{M}}^{-1}\mathrm{=}\left(\begin{array}{c}\\ {{V}}_{{1}}^{-1}& \mathrm{0}& \mathit{...}& \mathrm{0}\\ \mathrm{0}& {{V}}_{{2}}^{-1}& {\begin{array}{c}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& \mathrm{0}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& {\begin{array}{c}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& {\begin{array}{c}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& \mathrm{0}\\ \mathrm{0}& \mathit{...}& \mathrm{0}& {{V}}_{{n}}^{-1}\end{array}\right)\mathit{,}$(41)which follows in our particular case a χ² distribution with 2(n−1) degrees of freedom, or equivalently, (n−1)χ² distributions with two degrees of freedom. Equation (40) is probably the Mahalanobis distance mentioned without giving its expression in Adorf et al. (2006).

If ${\mathit{D}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{\left(}{\mu }\mathrm{\right)}$ follows a ${\mathit{\chi }}_{\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}$ distribution, then its square root, the distance D_M(μ), follows a χ_{d.o.f. = 2(n−1)} distribution.

We perform a statistical hypothesis test on a set of n sources, defining the null hyothesis H₀ as follows: all sources in the set are detections of the same “real” source. The alternative hypothesis H₁ would thus be: not all sources in the set are detections of the same “real” source; in other words the set of n sources contains at least one spurious source; or, expressed differently, the n sources are n observations of at least two distinct real sources. We adopt Fisher’s approach, that is we will reject the null hypothesis if, the null hypothesis being true, the observed data is significantly unlikely.

From now on, we indifferently write x or D_M the Mahalanobis distance. Assuming the null hypothesis is true, the “theoretical” probability we had to get the actual computed (square of) Mahalanobis distance is given by a Chi(-square) distribution with 2(n−1) degrees of freedom: $\begin{array}{ccc}\mathit{p}\mathrm{(}\mathit{X}\mathrm{=}\mathit{x}\mathrm{)}& \mathrm{=}& {\mathit{\chi }}_{\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{(}\mathit{X}\mathrm{=}\mathit{x}\mathrm{)}\mathrm{d}\mathit{X,}\\ \mathit{p}\mathrm{(}\mathit{X}\mathrm{=}{\mathit{x}}^{\mathrm{2}}\mathrm{)}& \mathrm{=}& {\mathit{\chi }}_{\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}\mathrm{(}\mathit{X}\mathrm{=}{\mathit{x}}^{\mathrm{2}}\mathrm{)}\mathrm{d}\mathit{X}\mathit{.}\end{array}$The probability we had to get an actual computed (square of) Mahalanobis distance less than or equal to a given threshold (or critical value) ${\mathit{k}}_{\mathit{\gamma }}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ is given by the value of the cumulative distribution function of a the Chi(-square) at the given threshold $\mathit{\gamma }\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}}\mathit{p}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}}{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{d}\mathit{X}\mathrm{=}{\mathit{F}}_{{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\mathrm{\right)}\mathit{.}$(44)We can indifferently work on x with the χ distribution or on x² with the χ² distribution. The threshold k_γ we obtain on x is simply the square root of the threshold ${\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}$ we obtain on x². Although we find the Chi test more natural in the present case, most astronomers are familiar with the Chi-square test.

In the framework of statistical hypothesis tests, it is the complementary cumulative distribution (or tail distribution) function which is usually used by defining the p-value $\mathit{p}\mathrm{-}\mathrm{value}\mathrm{=}{\mathrm{\int }}_{{\mathit{x}}^{\mathrm{2}}}^{\mathrm{+}\mathrm{\infty }}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{d}\mathit{X}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{F}}_{{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}}\mathrm{\left(}{\mathit{x}}^{\mathrm{2}}\mathrm{\right)}\mathit{,}$(45)and a significance level α defined by $\mathit{\alpha }\mathrm{=}{\mathrm{\int }}_{{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}}^{\mathrm{+}\mathrm{\infty }}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{d}\mathit{X}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{F}}_{{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\mathit{\gamma }$(46)is fixed. The null hypothesis is then rejected if p - value <α. In the Neyman-Pearson framework α is the type I error, or the false positive rate, that is the probability the null hypothesis has to be rejected (positive rejection test) while it is true (wrong/false decision). In our case we fix γ (hereafter called completeness), the fraction of real associations we “theoretically” select over all real associations. The candidates selection criterion, or fail of rejection criterion, we use is then ${\mathit{D}}_{\mathrm{M}}\mathrm{\left(}{\mu }\mathrm{\right)}\mathrm{\le }{\mathit{k}}_{\mathit{\gamma }}$(47)in which ${\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\mathrm{=}{\mathit{F}}_{{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}}^{-1}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}$ or, equivalently, ${\mathit{k}}_{\mathit{\gamma }}\mathrm{=}{\mathit{F}}_{{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}}^{-1}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}$. This inequality is equivalent to p - value <α. It is important to write “fail of rejection” since nothing proves that if Eq. (47) is satisfied the null hypothesis is true: at this point the selected set of sources is nothing else than a set of candidates. Nevertheless we do call region of acceptance the set of D_M(μ) values satisfying Eq. (47). This region of acceptance will be useful to define the domain of integration used to normalize likelihoods when computing probabilities for each hypothesis from Sect. 7. Its volume (see e.g. Eq. (64)) is the volume of the 2n-ellipsoid defined by M (see Eq. (40)) divided by the error ellipse associated to the weighted mean position μ_Σ and defined by V_Σ (it thus is a volume in a 2(n−1) space).

In practice, the value ${\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}$ is computed numerically using Newton’s method to solve ${\mathit{F}}_{{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{-}\mathit{\gamma }\mathrm{=}\mathrm{0}$. The initial guess we use is the approximate value returned by Eq. (A.3) of Inglot (2010).

The value of γ we fix is independent of n, the number of candidates. In practice we often set this input parameter to γ = 0.9973. In one dimension this value leads to k_γ = 3, that is the famous 3σ rule. It means that for 10 000 real associations in a dataset, we theoretically miss 27 of them by applying the candidate selection criterion. From now on we call this cross-correlation a χ_γ-match, or simply a χ-match.

In the particular case of two catalogues D_M(μ) follows a χ distribution with 2 degrees of freedom – that is a Rayleigh distribution – and k_{γ = 0.9973} = 3.443935. This latter value is used in the two-catalogues χ-match of Pineau et al. (2011b).

5.3. Iterative form: catalogue by catalogue

Somewhat similarly to the Bayes factor in Budavári & Szalay (2008, Sect. 6) it is noteworthy that Q_χ² can be computed iteratively, summing (n−1) successive χ² with two degrees of freedom computed from (n−1) successive two-catalogues cross-matches.

After each iteration, the new position to be used for the next cross-match is the weighted mean of all already matched positions and the new associated error is the error on this weighted mean. The strict equality between Eq. (48) and the non iterative form, for example Eq. (34), proves that the result is independent of the successive cross-matches order.

The maximum number of cross-matches to be performed must be known in advance in order to put an upper limit on k_γ since it depends on the degree of freedom of the total χ². The iteration formula is simply ${\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{2}}^{\mathit{n}}{\mathrm{(}{{\mu }}_{{{\Sigma }}_{{i}{-}{1}}}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}}^{⊺}\mathrm{(}{{V}}_{{{\Sigma }}_{{i}{-}{1}}}\mathrm{+}{{V}}_{\mathit{i}}{\mathrm{)}}^{-1}\mathrm{(}{{\mu }}_{{{\Sigma }}_{{i}{-}{1}}}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}$(48)in which $\begin{array}{ccc}{{V}}_{{{\Sigma }}_{{i}{-}{1}}}^{-1}& \mathrm{=}& \\ {{\mu }}_{{{\Sigma }}_{{i}{-}{1}}}& \mathrm{=}& {{V}}_{{\mathrm{\Sigma }}_{\mathit{i}\mathrm{-}\mathrm{1}}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{i}\mathrm{-}\mathrm{1}}{{V}}_{{k}}^{-1}{{\mu }}_{{k}}\mathit{.}\end{array}$We find it from the 2-catalogues case, for which (see Sect. A.4) ${\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}\mathrm{=}{\mathrm{(}{{\mu }}_{{1}}\mathrm{-}{{\mu }}_{{2}}\mathrm{)}}^{⊺}\mathrm{(}{{V}}_{\mathrm{1}}\mathrm{+}{{V}}_{\mathrm{2}}{\mathrm{)}}^{-1}\mathrm{(}{{\mu }}_{{1}}\mathrm{-}{{\mu }}_{{2}}\mathrm{)}\mathit{.}$(51)We can demonstrate by direct calculation that $\begin{array}{c}\det \end{array}\mathrm{(}{{V}}_{{1}}\mathrm{+}{{V}}_{{2}}\mathrm{)}\begin{array}{c}\det \end{array}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{=}\begin{array}{c}\det \end{array}{{V}}_{{1}}\begin{array}{c}\det \end{array}{{V}}_{{2}}$(52)and so, iteratively, we find the general expression $\underset{\mathit{i}\mathrm{=}\mathrm{2}}{\overset{\mathit{n}}{\mathrm{􏽙}}}\begin{array}{c}\det \end{array}\mathrm{(}{{V}}_{{{\Sigma }}_{{i}{-}{1}}}\mathrm{+}{{V}}_{{i}}\mathrm{)}\mathrm{=}\underset{\mathit{i}\mathrm{=}\mathrm{2}}{\overset{\mathit{n}}{\mathrm{􏽙}}}\frac{\det {{V}}_{{{\Sigma }}_{{i}{-}{1}}}\det {{V}}_{{i}}}{\det {{V}}_{{{\Sigma }}_{{i}}}}\mathrm{=}\frac{\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽑}}}\det {{V}}_{{i}}}{\det {{V}}_{{\Sigma }}}$(53)which is consistent with Eq. (38). The volume of the acceptance region of the statistical hypothesis test is the volume of a 2(n−1) dimensional ellipsoid. More precisely, it is the product of the previous equation Eq. (53) by the volume of a 2(n−1)-sphere of radius k_γ. This will be crucial when computing the rate of spurious associations.

5.4. Iterative form: by groups of catalogues

Instead of iterating over catalogues one by one, we can also perform G sub-cross-matches, each associating n_g distinct sources such that $\sum _{\mathit{g}\mathrm{=}\mathrm{1}}^{\mathit{G}}{\mathit{n}}_{\mathit{g}}\mathrm{=}\mathit{n}$. We note Q_{χ², { g }} the square of the Mahalanobis distance associated with the group g: ${\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}\mathit{,}\mathrm{\left\{}\mathit{g}\mathrm{\right\}}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{2}}^{{\mathit{n}}_{\mathit{g}}}{\mathrm{(}{{\mu }}_{{{\Sigma }}_{{i}{-}{1}}}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}}^{⊺}\mathrm{(}{{V}}_{{\mathrm{\Sigma }}_{\mathit{i}\mathrm{-}\mathrm{1}}}\mathrm{+}{{V}}_{\mathit{i}}{\mathrm{)}}^{-1}\mathrm{(}{{\mu }}_{{{\Sigma }}_{{i}{-}{1}}}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}\mathit{.}$(54)We show that we can compute Q_χ² iteratively from the G weighted mean positions μ_Σ{g} and their associated errors ${{V}}_{{\mathrm{\Sigma }}_{\mathrm{\left\{}\mathit{g}\mathrm{\right\}}}}^{-1}$. The square of the Mahalanobis distance can be written $\begin{array}{ccc}{\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}& \mathrm{=}& \sum _{\mathit{g}\mathrm{=}\mathrm{2}}^{\mathit{G}}{\mathrm{(}{{\mu }}_{{{\Sigma }}_{{g}{-}{1}}}\mathrm{-}{{\mu }}_{{{\Sigma }}_{{\left\{}{g}{\right\}}}}\mathrm{)}}^{⊺}\mathrm{(}{{V}}_{{\mathrm{\Sigma }}_{\mathit{g}\mathrm{-}\mathrm{1}}}\mathrm{+}{{V}}_{{\mathrm{\Sigma }}_{\mathrm{\left\{}\mathit{g}\mathrm{\right\}}}}{\mathrm{)}}^{-1}\mathrm{(}{{\mu }}_{{{\Sigma }}_{{g}{-}{1}}}\mathrm{-}{{\mu }}_{{{\Sigma }}_{{\left\{}{g}{\right\}}}}\mathrm{)}\\ & & \end{array}$(55)In other words, the square of the Mahalanobis distance is the sum of the square of the intra-group Mahalanobis distances plus the inter-group iterative one. With k being an index defined inside each of the G groups { g }$\begin{array}{ccc}{{V}}_{{{\Sigma }}_{{\left\{}{g}{\right\}}}}^{-1}& \mathrm{=}& \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{g}}}{{V}}_{{k}}^{-1}\mathit{,}\\ {{\mu }}_{{{\Sigma }}_{{\left\{}{g}{\right\}}}}& \mathrm{=}& {{V}}_{{{\Sigma }}_{{\left\{}{g}{\right\}}}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{g}}}{{V}}_{{k}}^{-1}{{\mu }}_{{k}}\mathit{,}\\ {{V}}_{{{\Sigma }}_{{g}{-}{1}}}^{-1}& \mathrm{=}& \sum _{{\mathit{g}}^{\mathrm{\prime }}\mathrm{=}\mathrm{1}}^{\mathit{g}\mathrm{-}\mathrm{1}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{{\mathit{g}}^{\mathrm{\prime }}}}{{V}}_{{k}}^{-1}\mathit{,}\\ {{\mu }}_{{{\Sigma }}_{{g}{-}{1}}}& \mathrm{=}& {{V}}_{{{\Sigma }}_{{g}{-}{1}}}\sum _{{\mathit{g}}^{\mathrm{\prime }}\mathrm{=}\mathrm{1}}^{\mathit{g}\mathrm{-}\mathrm{1}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{{\mathit{g}}^{\mathrm{\prime }}}}{{V}}_{{k}}^{-1}{{\mu }}_{{k}}\mathit{.}\end{array}$In fact, it is a straightforward generalization of the G = 2 groups case for which $\begin{array}{ccc}{{V}}_{{{\Sigma }}_{{g}{=}{1}}}^{-1}& \mathrm{=}& \mathrm{(}{{V}}_{{{\Sigma }}_{{\left\{}{1}{\right\}}}}\mathrm{+}{{V}}_{{{\Sigma }}_{{\left\{}{2}{\right\}}}}{\mathrm{)}}^{-1}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{V}}_{{{\Sigma }}_{{i}}}\mathrm{=}{{V}}_{{\Sigma }}^{-1}\mathit{,}\\ {{\mu }}_{{{\Sigma }}_{{g}{=}{1}}}& \mathrm{=}& {{V}}_{{{\Sigma }}_{{g}{=}{1}}}\mathrm{(}{{V}}_{{{\Sigma }}_{{\left\{}{1}{\right\}}}}^{-1}{{\mu }}_{{{\Sigma }}_{{\left\{}{1}{\right\}}}}\mathrm{+}{{V}}_{{{\Sigma }}_{{\left\{}{2}{\right\}}}}^{-1}{{\mu }}_{{{\Sigma }}_{{\left\{}{2}{\right\}}}}\mathrm{)}\mathit{,}\\ & \mathrm{=}& {{V}}_{{\Sigma }}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{V}}_{{i}}^{-1}{{\mu }}_{{i}}\mathrm{=}{{\mu }}_{{\Sigma }}\mathit{.}\end{array}$Here again, $\underset{\mathit{i}\mathrm{=}\mathrm{2}}{\overset{\mathit{G}}{\mathrm{􏽙}}}\begin{array}{c}\det \end{array}\mathrm{(}{\mathit{V}}_{{\mathrm{\Sigma }}_{\mathit{g}\mathrm{-}\mathrm{1}}}\mathrm{+}{\mathit{V}}_{{\mathrm{\Sigma }}_{\mathit{g}}}\mathrm{)}\mathrm{=}\underset{\mathit{g}\mathrm{=}\mathrm{2}}{\overset{\mathit{G}}{\mathrm{􏽙}}}\frac{\det {\mathit{V}}_{{\mathrm{\Sigma }}_{\mathit{g}\mathrm{-}\mathrm{1}}}\underset{\mathit{k}\mathrm{=}\mathrm{1}}{\overset{{\mathit{n}}_{\mathit{g}}}{\mathrm{􏽑}}}\det {\mathit{V}}_{\mathit{k}}}{\det {\mathit{V}}_{{\mathrm{\Sigma }}_{\mathit{g}}}}\mathrm{=}\frac{\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽑}}}\det {\mathit{V}}_{\mathit{i}}}{\det {\mathit{V}}_{{\mathrm{\Sigma }}_{\mathit{n}}}}\mathrm{·}$(63)Again, k_γ depends on the number of degrees of freedom of the total χ², thus on the total number of cross-correlated tables. It means that to be complete, all sub-cross-correlations must use the candidate selection threshold k_γ(2(n−1)) computed from the total number of tables instead of k_γ(2(n_g−1)) computed from the number of tables in a group.

5.5. Summary and Interpretation

Equations (34), (36), (39), (40), (48) and (55) are all equivalent and they lead to the same value, that is to the same squared Mahalanobis distance. All sources are retained as possible candidates if Eq. (47) is verified, so if the Mahalanobis distance is smaller or equal to k_γ. This threshold is the inverse of the cumulative χ distribution function at the chosen completeness γ, for 2(n−1) degrees of freedom.

As this criterion is no other than a χ-test criterion (or χ²-test criterion if we work on squared Mahalanobis distances) we call the result of such a criterion a χ-match.

The χ-match criterion defines a region of acceptance which is a 2(n−1)-ellipsoid of radius k_γ. Its volume is computed from Eq. (53): ${\mathrm{𝒱}}_{\mathit{n}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}\mathrm{=}{\left[\frac{\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\mathrm{􏽑}}\det {{V}}_{{i}}}{\det {{V}}_{{\Sigma }}}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{{\mathit{\pi }}^{\mathit{n}\mathrm{-}\mathrm{1}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}}{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}}\mathit{,}$(64)with ${\mathit{\pi }}^{\mathit{n}\mathrm{-}\mathrm{1}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{/}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}$ the volume of a 2(n−1)-sphere of radius k_γ. It will be later used to compute the expected number of spurious associations.

5.6. Comment on the “Bayesian cross-match” of Budavári & Szalay (2008)

We mention in Sect. 6.1 what appears to be a conceptual problem in calling B (Eq. (65)) a Bayes factor for more than two catalogues in the astrometrical part of Budavári & Szalay (2008).

Performing a cross-match by fixing a lower limit L on the “Bayes factor” B defined in Eq. (18) of Budavári & Szalay (2008) is no other than performing a χ-match with a significance level which depends both on the number of sources n and on the volume of the 2(n−1)-ellipsoid of radius 1. In fact, using the factor B of Budavári & Szalay (2008) in which w_i is the inverse of the cirular error on the position of the source i and φ_ij is the angular distance between sources i and j, we have the equivalence $\begin{array}{ccc}\mathit{B}& \mathrm{=}& \\ & \mathrm{\iff }& \frac{{\sum }_{\mathit{i}\mathit{<}\mathit{j}}{\mathit{w}}_{\mathit{i}}{\mathit{w}}_{\mathit{j}}{\mathit{\phi }}_{\mathit{ij}}^{\mathrm{2}}}{\sum {\mathit{w}}_{\mathit{i}}}\mathrm{\le }\mathrm{2}\ln \left(\frac{{\mathrm{2}}^{\mathit{n}\mathrm{-}\mathrm{1}}}{\mathit{L}}\frac{􏽑\mathit{w}i}{\sum {\mathit{w}}_{\mathit{i}}}\right)\mathrm{·}\end{array}$(65)We showed that the quantity on the left side of the inequality is equal to Eq. (39) in the present paper and thus follows a χ² distribution for “real” associations. It means that the “Bayesian” candidate selection criterion B ≥ L is equivalent to a χ² test having a significance level equal to $\mathit{\alpha }\mathrm{=}{\mathrm{\int }}_{\mathrm{2}\ln \left(\frac{{\mathrm{2}}^{\mathit{n}\mathrm{-}\mathrm{1}}}{\mathit{L}}\frac{􏽑\mathit{w}i}{\sum {\mathit{w}}_{\mathit{i}}}\right)}^{\mathrm{+}\mathrm{\infty }}{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}\mathit{.}$(66)The larger the volume of the 2(n−1)-ellipsoid of radius 1 (∝$\sum {\mathit{w}}_{\mathit{i}}\mathit{/}{}^{\mathrm{􏽑}}\mathit{w}_{\mathit{i}}$), the more “real” associations are missed and the less spurious associations are retrieved. We could replace the criterion B ≥ L by $\mathit{x}\mathrm{\le }\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{(}\mathit{n,}{}^{\mathrm{􏽑}}\mathit{V}_{\mathit{i}}\mathit{/}{\mathit{V}}_{\mathrm{\Sigma }}\mathrm{)}$. This is somewhere between the fixed radius cone search and the fixed significance level χ-match. The rate of missed “real” associations is not homogeneous but depends on the positional errors. Only if positional errors are constant in all catalogues, then the B ≥ L constraint becomes equal to the χ-match which is equal to a fixed radius cross-match.

6. Hypotheses from combinatorial considerations

A χ-match output is made of sets of associations, each set of associations containing one source per catalogue. For each set of associations we want to compute the probability all sources of the set have to come from a same actual source. In this section, especially in Sect. 6.2 we make explicit the sets { h_i } of hypotheses we have to formulate to compute probabilities of identification when cross-correlating n catalogues.

6.1. Generalities

Given a set { h_k } of pairwise disjoint hypotheses whose union is the entire set of possibilities, the law of total probabilities for an observable x is $\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{k}}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathrm{\left)}\mathit{p}\mathrm{\right(}{\mathit{h}}_{\mathit{i}}\mathrm{)}\mathit{.}$(67)Leading to Bayes’ theorem $\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{j}}\mathrm{\right|}\mathit{x}\mathrm{)}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{j}}\mathrm{\left)}\mathit{p}\mathrm{\right(}{\mathit{h}}_{\mathit{j}}\mathrm{)}}{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{k}}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathrm{\left)}\mathit{p}\mathrm{\right(}{\mathit{h}}_{\mathit{i}}\mathrm{)}}\mathrm{·}$(68)We stress that Bayes’ factor (also called likelihood ratio) is defined only in cases involving two and only two hypotheses $\mathit{LR}\mathrm{=}\mathit{K}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathrm{1}}\mathrm{)}}{\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathrm{2}}\mathrm{)}}\mathit{,}$(69)and is used when no trustworthy priors p(h₁) and p(h₂) are available. We can transform any set of pairwise disjoint hypotheses into two disjoint hypotheses. In this case, using the negation notation ¬$\mathit{LR}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{j}}\mathrm{)}}{\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathrm{\neg }{\mathit{h}}_{\mathit{j}}\mathrm{)}}\mathit{,}$(70)with $\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathrm{\neg }{\mathit{h}}_{\mathit{j}}\mathrm{)}\mathrm{=}\frac{\sum _{\mathit{i}\ne \mathit{j}}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathrm{\left)}\mathit{p}\mathrm{\right(}{\mathit{h}}_{\mathit{i}}\mathrm{)}}{\mathit{p}\mathrm{\left(}\mathrm{\neg }{\mathit{h}}_{\mathit{j}}\mathrm{\right)}}\mathit{,}$(71)and $\mathit{p}\mathrm{\left(}\mathrm{\neg }{\mathit{h}}_{\mathit{j}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\ne \mathit{j}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\mathit{.}$(72)Such a likelihood ratio (Eq. (70)) is not interesting since it is not only computed from likelihoods, but also from priors.

The term B(H | K) in Budavári & Szalay (2008, Eq. (8)) is improperly called Bayes factor when dealing with more than two catalogues. As a matter of fact, the union of the two hypotheses – all sources are from the same real source and each sources is from a distinct real source – is only a subset of all possibilities so the law of total probabilities and hence Bayes’ formula are not valid.

In Pineau et al. (2011b) the term LR(r) in Eq. (9) is also improperly called likelihood ratio since a likelihood is a probability density function and so integrates to 1 over its domain of definition. It is obviously not the case of dp(r | spur) in Eq. (8). The built quantity is related to the ratio between the probability the association has to be “real” over the probability it has to be spurious, but formally it is not a likelihood ratio. The very same “abuse of term” is made in Wolstencroft et al. (1986; who, moreover, adds a prior in the likelihood ratio), in Rutledge et al. (2000), Brusa et al. (2007) and probably other publications.

6.2. Possible combinations and the Bell number

Let’s suppose we have selected one set of n distinct sources from n different catalogues, one source per catalogue. Those n sources possibly are n detections of k distinct real sources, with k ∈ [1,n]. The case k = 1 corresponds to the situation where all sources are n observations of the same real source and the case k = n corresponds to the situation where there are n distinct real sources detected independently, one in each catalogue.

We call A the source from catalogue number one, B the source from catalogue number two and so on.

6.2.1. Two-catalogues case: two hypotheses

The classical two-catalogues case is trivial. We formulate only two hypotheses:

• AB, the match is a real match, the two sources are two observations of a same real source, that is k = 1;

• A_B, the match is spurious, the two sources are two observations of two different real sources, that is k = 2.

6.2.2. Three-catalogues case: five hypotheses

For three sources A, B and C from three different catalogues, we formulate five hypotheses:

• ABC, all three sources come from a same real source, that is k = 1;

• AB_C, A and B are from a same real source and C is from a different real source, that is k = 2;

• AC_B, A and C are from a same real source and B is from a different real source, that is k = 2;

• A_BC, B and C are from a same real source and A is from a different real source, that is k = 2;

• A_B_C, all three sources are from three different real sources, that is k = 3.

6.2.3. Four-catalogues case: 15 hypotheses

For four sources A, B, C and D we have to formulate 15 hypotheses:

• ABCD, when k = 1;

• ABC_D, ABD_C, ACD_B and BCD_A, but also

• AB_CD, AC_BD and AD_BC for k = 2;

• AB_C_D, AC_B_D, AD_B_C, BC_A_D, BD_A_C and DC_A_B when k = 3;

• A_B_C_D when k = 4.

6.2.4. n-catalogues case: Bell number of hypotheses

We now generalize to n catalogues. For each possible value of k, the number of ways the set of n sources can be partitioned into k non-empty subsets – each subset correspond to a real source – is given by the Stirling number of the second kind denoted ${{{\mathrm{\{}}^{}}_{\mathit{k}}^{\mathit{n}}}^{\mathrm{\}}}$. The total number of hypotheses to be formulated is equal to the Bell number. The Bell number counts the number of partitions of a set and is given by ${\mathit{B}}_{\mathit{n}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{n}}\left\{\begin{array}{c}\mathit{n}\\ \mathit{k}\end{array}\right\}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{n}\mathrm{-}\mathrm{1}}{\mathit{C}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathit{k}}{\mathit{B}}_{\mathit{k}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{n}\mathrm{-}\mathrm{1}}\frac{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}}{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}\mathrm{)}\mathrm{!}\mathit{k}\mathrm{!}}{\mathit{B}}_{\mathit{k}}\mathit{.}$(73)Its seven first values are provided in Table 2 and a graphic illustration representing all possible partitions for five catalogues is provided in Fig. 1⁸.

We face a combinatorial explosion of the number of hypotheses to be tested when increasing the number of catalogues. Although the theoretical developments presented here deal with any number of catalogues, the exhaustive analysis may be in practice limited to a few catalogues (n< 10).

Hereafter we note h_i the hypothesis number i, we explicit it with letters for example h_AB, and we note h_{k = i} an hypothesis in which n observed sources are associated to i real sources.

7. Frequentist estimation of spurious associations rates and priors

We have defined a candidate selection criterion to perform χ-matches. We recall that we note x the Mahalanobis distance, and we note s the “event” x ≤ k_γ, that is a given set of sources satisfies the selection criterion.

In a first step we want to estimate the number of “fully spurious” associations we would expect to find in a χ-match output and derive the prior p(h_{k = n} | s) from this estimate. By “fully spurious” we mean that each candidate from each catalogue is actually associated with a different “real” source. A good such estimate is simply the mean sky area of the test acceptance region (see Eq. (64)) over all possible sources of all catalogues, multiplied by the number of sources in one of the catalogues and by the density of sources in the other ones. Written differently for n catalogues of n_i sources each, on a common surface area Ω, it leads to an estimated number of spurious associations equals to: $\mathit{n̂}{\mathrm{\Omega }}_{\mathrm{spur}}\mathrm{=}\frac{{\mathit{\pi }}^{\mathit{n}\mathrm{-}\mathrm{1}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}}{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}{\mathrm{\Omega }}^{\mathit{n}\mathrm{-}\mathrm{1}}}\sum _{{\mathit{i}}_{\mathrm{1}}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{1}}}\sum _{{\mathit{i}}_{\mathrm{2}}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{2}}}\mathit{...}\sum _{{\mathit{i}}_{\mathit{n}}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{n}}}{\left[\frac{\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽑}}}\det {{V}}_{{{i}}_{{j}}}}{\det {{V}}_{{\Sigma }}}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}$(74)Or, having histograms or more generally discretized positional error distributions: $\mathit{n̂}{\mathrm{\Omega }}_{\mathrm{spur}}\mathrm{=}\frac{{\mathit{\pi }}^{\mathit{n}\mathrm{-}\mathrm{1}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}}{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}{\mathrm{\Omega }}^{\mathit{n}\mathrm{-}\mathrm{1}}}\sum _{{\mathit{b}}_{\mathrm{1}}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{1}}}\sum _{{\mathit{b}}_{\mathrm{2}}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{2}}}\mathit{...}\sum _{{\mathit{b}}_{\mathit{n}}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{n}}}\underset{\mathit{k}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽙}}}{\mathit{c}}_{{\mathit{b}}_{\mathit{k}}}{\left[\frac{\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽑}}}\det {{V}}_{{{b}}_{{j}}}}{\det {{V}}_{{\Sigma }}}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(75)in which N_k are the numbers of bins in histograms – or number of points in a discrete distribution– and c_bk are number of counts in given bins of a histogram. The number of counts may be replaced by the value of the discrete distribution (or weight w_bk) times the number of elements: c_bk = n_kw_bk.

To perform quick estimations using only a one dimensional error histogram per catalogue, we approximate elliptical errors by circular errors of same surface area.

The remainder of this section explains how we can compute priors from the rate of “fully spurious” associations and the number of associations found in all possible sub-cross-matches.

7.1. Case of two catalogues

Let’s suppose that we have two catalogues A and B and each catalogue contains only one source in the common surface area Ω. We note μ_a1, V_a1 and μ_b1, V_b1 the position of the source and associated covariance matrix in A and B respectively. If we fix the position μ_a1 of the first source, the second source will be associated with the first one by a χ_γ-match if Eq. (47) is satisfied. So if the second source is located in an ellipse of surface area $\mathit{\pi }\sqrt{\det \mathrm{(}{{V}}_{{a}{1}}\mathrm{+}{{V}}_{{b}{1}}\mathrm{)}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}$ centred around the position of the first source. We temporarily waive the ISO 80000-2 notation detM and replace it by the equivalent and more compact notation | M |. We also replace V_a1 + V_b1 by V_1,1 to rewrite the last term in the pithier form $\mathit{\pi }\mathrm{|}{{V}}_{{1}{,}{1}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}$. We now suppose that both sources are unrelated and that μ_a1 and μ_b1 are uniformly distributed in Ω. Then, neglecting border effects, the probability that the two sources are associated by chance when performing a χ_γ-match is given by the ratio of the acceptance ellipse to the total surface area Ω: $\mathit{p}\mathrm{=}\frac{{\mathrm{\int }}_{\mathit{x}\mathrm{=}\mathrm{0}}^{\mathit{x}\mathrm{=}{\mathit{k}}_{\mathit{\gamma }}}\mathrm{d}\mathrm{(}{{\mu }}_{{a}{1}}\mathrm{-}{{\mu }}_{{b}{1}}\mathrm{)}}{\mathrm{\Omega }}\mathrm{=}\frac{\mathrm{|}{\mathit{V}}_{\mathrm{1}\mathit{,}\mathrm{1}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{k}}_{\mathit{\gamma }}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathit{x}\mathrm{d}\mathit{x}\mathrm{d}\mathit{\theta }}{\mathrm{\Omega }}\mathrm{=}\frac{\mathit{\pi }\mathrm{|}{\mathit{V}}_{\mathrm{1}\mathit{,}\mathrm{1}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}}{\mathrm{\Omega }}\mathrm{·}$(76)We now suppose that the second catalogue contains n_B sources uniformly distributed in Ω. And if all of them are unrelated to the source of the first catalogue, then the estimated number of spurious associations is simply the sum of the previous probability over the n_B sources of the second catalogue $\mathit{n̂}\mathit{A}\mathrm{\_}\mathit{B}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{B}}}{\mathit{p}}_{\mathrm{1}\mathit{,j}}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{B}}}\frac{\mathit{\pi }\mathrm{|}{\mathit{V}}_{\mathrm{1}\mathit{,j}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}}{\mathrm{\Omega }}\mathrm{·}$(77)We now suppose that the first catalogue contains n_A sources also uniformly distributed in Ω, all unrelated to catalogue B sources. Still neglecting border effects, the estimated number of spurious associations is simply the sum of the previous estimation over all catalogue A sources $\mathit{n̂}\mathit{A}\mathrm{\_}\mathit{B}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{A}}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{2}}}{\mathit{p}}_{\mathit{i,j}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{1}}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{B}}}\frac{\mathit{\pi }\mathrm{|}{\mathit{V}}_{\mathit{i,j}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}}{\mathrm{\Omega }}\mathrm{·}$(78)In practice, evaluating this quantity can be time-consuming since we have to compute and sum n_A × n_B terms. Fortunately, we can evaluate it exactly for circular errors and approximately for elliptical errors computing only n_A + n_B terms. In fact $\begin{array}{ccc}\mathrm{|}{\mathit{V}}_{\mathit{i,j}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}& \mathrm{=}& \mathrm{|}{\mathit{V}}_{\mathit{ai}}\mathrm{+}{\mathit{V}}_{\mathit{bj}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\\ & \mathrm{=}& \left(\mathrm{|}{\mathit{V}}_{\mathit{ai}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{|}{\mathit{V}}_{\mathit{bj}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)\sqrt{\mathrm{1}\mathrm{+}\frac{\mathit{C}}{\mathrm{\left(}\mathrm{\right|}{\mathit{V}}_{\mathit{ai}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{|}{\mathit{V}}_{\mathit{bj}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{)}}^{\mathrm{2}}}}\\ & \mathrm{\approx }& \left(\mathrm{|}{\mathit{V}}_{\mathit{ai}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{|}{\mathit{V}}_{\mathit{bj}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)\mathrm{\times }\\ & & \mathrm{\times }\left(\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{C}}{\mathrm{\left(}\mathrm{\right|}{\mathit{V}}_{\mathit{ai}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{|}{\mathit{V}}_{\mathit{bj}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{-}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{+}\mathit{...}\hspace{0.17em}\right)\end{array}$in which $\begin{array}{ccc}\mathit{C}& \mathrm{=}& \mathrm{(}{\mathit{\sigma }}_{{\mathit{x}}_{\mathit{i}}}{\mathit{\sigma }}_{{\mathit{y}}_{\mathit{j}}}\mathrm{-}{\mathit{\sigma }}_{{\mathit{y}}_{\mathit{i}}}{\mathit{\sigma }}_{{\mathit{x}}_{\mathit{j}}}{\mathrm{)}}^{\mathrm{2}}\\ & & \mathrm{+}\mathrm{2}{\mathit{\sigma }}_{{\mathit{x}}_{\mathit{i}}}{\mathit{\sigma }}_{{\mathit{y}}_{\mathit{i}}}{\mathit{\sigma }}_{{\mathit{x}}_{\mathit{j}}}{\mathit{\sigma }}_{{\mathit{y}}_{\mathit{j}}}\left(\mathrm{1}\mathrm{+}{\mathit{\rho }}_{\mathit{i}}{\mathit{\rho }}_{\mathit{j}}\mathrm{-}\sqrt{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathit{j}}^{\mathrm{2}}\mathrm{)}}\right)\mathit{,}\end{array}$(83)and thus $\mathit{C}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}{\mathit{\rho }}_{\mathit{i}}\mathrm{=}{\mathit{\rho }}_{\mathit{j}}\mathrm{=}\mathrm{0}\mathrm{;}\\ & \mathrm{iferrorsarecircular}\mathit{.}\end{array}$(84)For ordinary ellipses, that is ellipses having a position angle different from 0 and π/ 2, the approximation is valid if $\mathit{C}\mathrm{\ll }\mathrm{\left(}\mathrm{\right|}{\mathit{V}}_{\mathit{ai}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{|}{\mathit{V}}_{\mathit{bj}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{)}}^{\mathrm{2}}$. In the particular case of circular errors, Eq. (78) becomes $\mathit{n̂}\mathit{A}\mathrm{\_}\mathit{B}\mathrm{=}{\mathit{n}}_{\mathit{A}}{\mathit{n}}_{\mathit{B}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\frac{{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{+}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}}{\mathrm{\Omega }}\mathit{,}$(85)in which ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}$ and ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}$ are the mean surface area of all positional error ellipses in catalogues A and B respectively: ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{n}}_{\mathit{A}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{A}}}\mathit{\pi }\mathrm{|}{\mathit{V}}_{\mathit{ai}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{and}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{n}}_{\mathit{B}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{B}}}\mathit{\pi }\mathrm{|}{\mathit{V}}_{\mathit{bi}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathit{.}$(86)For simple circular errors σ_ai and σ_bi, this reduces to ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{n}}_{\mathit{A}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{A}}}\mathit{\pi }{\mathit{\sigma }}_{\mathit{ai}}^{\mathrm{2}}\mathrm{and}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{n}}_{\mathit{B}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{B}}}\mathit{\pi }{\mathit{\sigma }}_{\mathit{bi}}^{\mathrm{2}}\mathit{.}$(87)If errors are constant for all sources in each catalogue, this reduces to ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{=}\mathit{\pi }{\mathit{\sigma }}_{\mathit{a}}^{\mathrm{2}}\mathrm{and}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}\mathrm{=}\mathit{\pi }{\mathit{\sigma }}_{\mathit{b}}^{\mathrm{2}}\mathit{.}$(88)These estimates based on geometrical considerations have the advantage of being very fast to compute.

Theoretically, we should remove from the double summation in Eq. (78) the pairs (i,j) which are real associations. We have no mean to do this since we do not known in advance the result of the cross-identification. Fortunately this effect is negligible in common cases. Indeed, if the result of the cross-match of the two catalogues contains n_AB real associations – that is sources of both catalogues from a same real source – and supposing that the positional error distribution of sources having a counterpart is similar to the global error distribution we should modify Eq. (85) by$\mathit{n̂}\mathit{A}\mathrm{\_}\mathit{B}\mathrm{=}\mathrm{(}{\mathit{n}}_{\mathit{A}}{\mathit{n}}_{\mathit{B}}\mathrm{-}{\mathit{n}}_{\mathit{AB}}\mathrm{)}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\frac{{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{+}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}}{\mathrm{\Omega }}\mathrm{·}$(89)In practice this estimate will tend to be overestimated since the distribution of sources in a catalogue cannot be uniform because of the limited angular resolution preventing the detection of very close sources in a same image. This effect is usually deemed to be of negligible importance. However one can detect its presence in particular circumstances. For instance, if the actual counterpart is located in the wings of a much brighter nearby source it may not be detected. This effect probably accounts for the presence of a fraction of the stellar identifications in high Galactic latitude X-ray surveys, in particular those with a much higher F_x/F_opt flux ratios and harder X-ray spectra than normal for active coronae in which cases a faint AGN may be the correct identification (Watson 2012; Menzel et al. 2016). One way to account for this effect and to limit the overestimation is to remove from the surface area Ω small areas around each source. The value of those areas depends for example on the source brightness. In addition, again because of the angular resolution: for real associations in catalogues having similar positional errors, the chance a source has to be also associated with a spurious source is low. More precisely, the start of the Poisson distribution will be truncated. In extreme cases in which the Poisson distribution is truncated for x<k_γ, meaning that sources in a real association cannot be part of a spurious association, we should remove those sources from the estimate . We thus have to rewrite the previous equation Eq. (89) as $\mathit{n̂}\mathit{A}\mathrm{\_}\mathit{B}\mathrm{=}\mathrm{(}{\mathit{n}}_{\mathit{A}}\mathrm{-}{\mathit{n}}_{\mathit{AB}}\mathrm{)}\mathrm{(}{\mathit{n}}_{\mathit{B}}\mathrm{-}{\mathit{n}}_{\mathit{AB}}\mathrm{)}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\frac{{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{+}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}}{\mathrm{\Omega }}\mathrm{·}$(90)Knowing the total number of associations, n_T, resulting from the χ-match, we can estimate from Eq. (89) the number of spurious associations, and thus the number of real associations is estimated by $\mathit{n̂}\mathit{AB}\mathrm{=}\frac{{\mathit{n}}_{\mathit{T}}\mathrm{-}{\mathit{n}}_{\mathit{A}}{\mathit{n}}_{\mathit{B}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\frac{{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{+}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}}{\mathrm{\Omega }}}{\mathrm{1}\mathrm{-}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\frac{{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{+}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}}{\mathrm{\Omega }}}\mathrm{·}$(91)If mean error ellipses in both catalogues are very small compared to the total surface area – that is ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{+}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}\mathrm{\ll }\mathrm{\Omega }$ – we can use the approximation $\mathit{n̂}\mathit{AB}\mathrm{\approx }{\mathit{n}}_{\mathit{T}}\mathrm{-}{\mathit{n}}_{\mathit{A}}{\mathit{n}}_{\mathit{B}}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\frac{{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}}}\mathrm{+}{\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}}}}{\mathrm{\Omega }}\mathit{,}$(92)which is equivalent to using directly equation Eq. (85), that is without taking care of removing real associations. is but an estimate and nothing prevents it from being negative due to count statistics in cross-matches with very few real associations and a lot of spurious associations. In practice, we have to define a lower limit such as .

Hence we can estimate the priors in the sample of associations satisfying the selection criterion (s) $\begin{array}{ccc}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{AB}}\mathrm{\right|}\mathit{s}\mathrm{)}& \mathrm{=}& \\ \mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{A}\mathrm{\_}\mathit{B}}\mathrm{\right|}\mathit{s}\mathrm{)}& \mathrm{=}& \end{array}$After a first two-catalogues cross-match, we may compare the expected histogram of detV_i,j for spurious associations with the same histogram obtained from all associations. We may then derive the estimated distribution of this quantity (detV_i,j) for “real” associations and compute the two likelihoods p(detV_i,j | h_AB,s) and p(detV_i,j | h_{A_B},s).

Similarly we may build the histograms of the quantity detV_Σ for both the spurious and the “real” associations. This quantity is the determinant of the covariance matrix – that is the positional error – associated with the weigthed mean positions. We proceeded likewise in Pineau et al. (2011b) using the “likelihood ratio” (see our comment on the abuse of term likelihood ratio) quantity instead of positional uncertainties.

7.2. Case of three catalogues

We recall that for 3 catalogues, the output contains five components (see Sect. 6.2.2): ABC, AB_C, A_BC, AC_B, A_B_C. We would like to estimate the number of spurious associations, that is the number of associations in the four components other than ABC. To do so, we need to perform the three two-catalogue cross-matches A with B, A with C and B with C. We are thus able to estimate n_AB and n_{A_B}, n_AC and n_{A_C} and finally n_BC and n_{B_C} respectively. To compute n_{AB_C}, we proceed like in the previous section considering the two catalogues AB and C. AB is the result of the χ-match of A with B: the positions in catalogue AB are the weighted mean positions (μ_Σ, Eq. (30)) of associated A and B sources and the associated errors (or covariance matrix) are given by V_Σ (Eq. (29)). The only difference with the two-catalogues case is that for the first catalogue (AB) we replace the simple mean elliptical error surface ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{AB}}}$ over the n_TAB entries by the weighted mean accounting for the probabilities the AB associations have to be “real” (i.e. not spurious) ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{AB}}}\mathrm{=}\frac{\mathrm{1}}{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{{\mathit{T}}_{\mathit{AB}}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{AB}}\mathrm{\right|}\mathit{...}\mathrm{)}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{{\mathit{T}}_{\mathit{AB}}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{AB}}\mathrm{\right|}\mathit{...}\mathrm{\left)}\mathit{\pi }\mathrm{\right|}{\mathit{V}}_{{\mathrm{\Sigma }}_{\mathit{AB}}\mathit{i}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathit{,}$(95)in which p(h_AB | ...) is the probability the association has to be a real association knowing some parameters (“...”), and V_ΣABi is the covariance matrix of the error on the weighted mean position i (see Sect. 5.3, particulary Eq. (49)). Such a probability will be computed in the next sections. We then compute ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{C}}}$ and derive like in the two catalogues case replacing A by AB and B by C. Similarly to , we can estimate and .

We now want to estimate n_{A_B_C}, with a result which is independent from the cross-correlation order. Although we may use Eq. (74), it is possibly time consuming. Another solution is to use its discretized form Eq. (75). To do so quickly at the cost of an approximation we may circularize the errors by replacing the coefficients of the covariance matrix V by a single error equal to $\sqrt{\det {V}}$ and setting the correlation (or covariance) parameter equal to zero. It means that the new covariance matrix is diagonal and both diagonal elements are equal to $\sqrt{\det {V}}$. We choose this value to preserve the surface area of the two-dimensional error since the determinant (∝area) of the circular error equals the determinant (∝area) of the ellipse. This approximation is the same as the one made in the previous section. For each catalogue we then make the histogram of $\sqrt{\det {V}}$ values using steps of for example one mas and we apply Eq. (75). In this case – circular errors – we simplify the equation using $\begin{array}{c}\det \end{array}{{V}}_{{\Sigma }}\mathrm{=}\frac{\mathrm{1}}{{\left(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{1}}{\sqrt{\det {{V}}_{{i}}}}\right)}^{\mathrm{2}}}$(96)and thus $\frac{\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽑}}}\det {{V}}_{{i}}}{\det {{V}}_{{\Sigma }}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\underset{\mathit{j}\mathrm{=}\mathrm{1}\mathit{,j}\ne \mathit{i}}{\overset{\mathit{n}}{\mathrm{􏽙}}}\sqrt{\begin{array}{c}\det \end{array}{{V}}_{{j}}}\mathit{.}$(97)We do not use this last form but give it for comparison with the denominator of Eq. (17) in Budavári & Szalay (2008).

Another option is to compute the number of “fully” spurious associations three times by following what was done in the previous section (and in the beginning of this section), but computing ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}\mathrm{\_}\mathit{B}}}$ instead of ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{AB}}}$. Similarly to Eq. (95): ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}\mathrm{\_}\mathit{B}}}\mathrm{=}\frac{\mathrm{1}}{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{{\mathit{T}}_{\mathit{AB}}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{A}\mathrm{\_}\mathit{B}}\mathrm{\right|}\mathit{...}\mathrm{)}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{{\mathit{T}}_{\mathit{AB}}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{A}\mathrm{\_}\mathit{B}}\mathrm{\right|}\mathit{...}\mathrm{\left)}\mathit{\pi }\mathrm{\right|}{\mathit{V}}_{{\mathrm{\Sigma }}_{\mathit{AB}}\mathit{i}}{\mathrm{|}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathit{.}$(98)Computing ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{C}}}$ we derive . Similarly we can compute ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{A}\mathrm{\_}\mathit{C}}}$ and ${\overline{\mathrm{\Omega }}}_{{\mathit{e}}_{\mathit{B}\mathrm{\_}\mathit{C}}}$ and estimate the number of fully spurious associations taking the mean of , and .

Having the estimated number of associations being part of the components AB_C, AC_B, A_BC and A_B_C plus knowing the total number of associations n_T, we are able to estimate and to compute the priors, for example $\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{ABC}}\mathrm{\right|}\mathit{s}\mathrm{)}\mathrm{=}\frac{\mathit{n̂}\mathit{ABC}}{{\mathit{n}}_{\mathit{T}}}\mathrm{·}$(99)

7.3. Case of n catalogues

We can easily generalise the previous section using recursion. For n = 4 catalogues, we estimate the number of associations in component A_B_C_D knowing the number of associations in the result of the four-catalogue χ-match and estimating recursively (from the three-catalogue χ-matches) the number of associations in the 14 other components (AB_C_D, ...). So for n catalogues, the total number of distinct (sub)-cross-matches to be performed to compute all priors recursively is ${\mathit{N}}_{\mathit{\chi }\mathrm{-}\mathrm{match}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{2}}^{\mathit{n}\mathrm{-}\mathrm{1}}{\mathit{C}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathit{k}}$(100)in which terms ${\mathit{C}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathit{k}}$ are the binomial coefficients (n−1) ! /(k ! (n−1−k) !). For five catalogues, N_χ−match = 26 and for six catalogues, N_χ−match = 57.

8. Probability of being χ-matched under hypothesis h_i

In this section we compute p(s | h_i), the probability that n sources from n distinct catalogues have to satisfy the candidate selection criteria under hypothesis h_i. We will show in section Sect. 9 that the p.d.f. of the Mahalanobis distance for χ-match associations under hypothesis h_i is the p.d.f. of the Mahalanobis distance without applying the candidate selection criteria, normalized by the probability p(s | h_i) we compute in this section: $\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathrm{)}}{\mathit{p}\mathrm{\left(}\mathit{s}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathrm{)}}\mathrm{·}$(101)We show here that p(s | h_i) is proportional to the integral we note I_hi,n(k_γ) (see Eq. (107)) which is independent of positional uncertainties and which also plays a role in Sect. 10. We will see also that p(s | h_i) and I_hi,n(k_γ) can be simplified to p(s | h_k) and I_k,n(k_γ) respectively, that is the probability n sources from n distinct catalogues have to satisfy the candidate selection criteria knowing they are actually associated with k distinct real sources.

If k = 1, that is all sources are from a same real source, we – logically – find p(s | h_{k = 1}) = γ, the cumulative χ distribution function evaluated at the threshold k_γ.

If k = n, all sources are spurious, we – also logically (see Eq. (74)) – find for I_{k = n,n}(k_γ) the volume of a 2(n−1)-dimensional sphere of radius k_γ, and p(s | h_{k = n}) equals the volume of the 2(n−1)-dimensional ellipsoid defined by the test acceptance region divided by the common χ-match surface area raised to the power of the number of χ-matches (i.e. n−1).

We note x the total Mahalanobis distance, that is the square root of Eq. (36). The vectorial form x = (x₁,x₂,...x_n−1) denotes the n−1 terms, also Mahalanobis distances, which are summed in the catalogue by catalogue iterative form Eq. (48). We rewrite this equation with the new notations ${\mathit{x}}^{\mathrm{2}}\mathrm{=}{\mathit{x}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{x}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{+}{\mathit{x}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{2}}\mathit{.}$(102)So x is the radius of an hypersphere in the n−1 successive Mahalanobis distances space. The relation between x and x of dimension n−1 is the polar transformation F:Rⁿ⁻¹ → Rⁿ⁻¹, (x₁,x₂,...,x_n−1) = F(x,θ₁,...,θ_n−2), $\mathit{F}\mathrm{:}\left(\begin{array}{c}\\ \mathit{x}\\ {\mathit{\theta }}_{\mathrm{1}}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\end{array}\right)\mathrm{\to }\left(\begin{array}{c}\\ {\mathit{x}}_{\mathrm{1}}\mathrm{=}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{x,}{\mathit{\theta }}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{\right)}\\ {\mathit{x}}_{\mathrm{2}}\mathrm{=}{\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{x,}{\mathit{\theta }}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{\right)}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{x}}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{=}{\mathit{f}}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{\left(}\mathit{x,}{\mathit{\theta }}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{\right)}\end{array}\right)\mathit{,}$(103)with ${\mathit{f}}_{\mathit{j}}\mathrm{\left(}\mathit{x,}{\mathit{\theta }}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathit{j}\mathrm{=}\mathrm{1}\mathrm{;}\\ & \mathrm{\forall }\mathit{j}\mathit{>}\mathrm{1.}\end{array}$(104)The associated differential transform is $\mathrm{d}{\mathit{x}}_{\mathrm{1}}\mathrm{d}{\mathit{x}}_{\mathrm{2}}\mathit{...}\mathrm{d}{\mathit{x}}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{=}\mathrm{\left|}\begin{array}{c}\det \end{array}{\mathit{J}}_{\mathit{F}}\mathrm{\right(}\mathit{x,}{\mathit{\theta }}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{\left)}\mathrm{\right|}\mathrm{d}\mathit{x}\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathit{,}$(105)with J_F the determinant of the Jacobian of F which is for example computed in Stuart & Ord (1994, Chap. II “Exact sampling distributions”, p. 375): $\begin{array}{c}\det \end{array}{\mathit{J}}_{\mathit{F}}\mathrm{=}{\mathit{x}}^{\mathit{n}\mathrm{-}\mathrm{2}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{3}}{\mathrm{􏽙}}}{\cos }^{\mathit{n}\mathrm{-}\mathit{i}\mathrm{-}\mathrm{2}}{\mathit{\theta }}_{\mathit{i}}\mathit{.}$(106)We now define the I_k,n(k_γ) integral which will be crucial in the next sections ${\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{x}\mathrm{=}\mathrm{0}}^{\mathit{x}\mathrm{\le }{\mathit{k}}_{\mathit{\gamma }}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathit{k}}{\mathrm{􏽙}}}{\mathit{\chi }}_{\mathrm{2}}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\underset{\mathit{i}\mathrm{=}\mathit{n}\mathrm{-}\mathit{k}\mathrm{+}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{i}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathrm{d}{\mathit{x}}_{\mathit{i}}$(107)in which k denotes the hypothetical number of real sources and so ranges from 1 to n. Written this way, the integral is simpler than the equivalent form: ${\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{x}\mathrm{=}\mathrm{0}}^{\mathit{x}\mathrm{\le }{\mathit{k}}_{\mathit{\gamma }}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{{\mathit{n}}_{\mathrm{1}}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{\chi }}_{\mathrm{2}}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\mathrm{d}{\mathit{x}}_{\mathit{i}}\left[\underset{\mathit{g}\mathrm{=}\mathrm{2}}{\overset{\mathit{k}}{\mathrm{􏽙}}}\left(\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{g}\mathrm{-}\mathrm{1}}\mathrm{d}{\mathit{x}}_{\mathit{g}\mathrm{-}\mathrm{1}}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{{\mathit{n}}_{\mathit{g}}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{\chi }}_{\mathrm{2}}\mathrm{\left(}{\mathit{x}}_{\mathit{j}}\mathrm{\right)}\mathrm{d}{\mathit{x}}_{\mathit{j}}\right)\right]$(108)in which we have k groups containing each n_g sources associated to a same real source so that $\sum _{\mathit{g}\mathrm{=}\mathrm{1}}^{\mathit{k}}{\mathit{n}}_{\mathit{g}}\mathrm{=}\mathit{n}$ (see the iterative candidate selection by groups of catalogues in Sect. 5.4). Here Eq. (55) takes the form ${\mathit{x}}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{g}\mathrm{=}\mathrm{2}}^{\mathit{k}}{\mathit{x}}_{\mathit{g}\mathrm{-}\mathrm{1}}^{\mathrm{2}}\mathrm{+}\sum _{\mathit{g}\mathrm{=}\mathrm{1}}^{\mathit{k}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{g}}}{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}$(109)where x_g−1 are inter-group Mahalanobis distances and x_i are intra-group Mahalanobis distances. In this version of the formula, we suppose that we iteratively cross-correlate the catalogues by groups. We suppose that each group corresponds to one real source. So inside each group, we multiply Rayleigh distributions and when associating each group, we multiply by two-dimensional Poisson distributions.

We compute I_k,n using recursive integration by parts, leading to (see Sects. A.5.3 and A.6) ${\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathit{k}\mathrm{=}\mathrm{1}\mathrm{;}\\ & \mathrm{if}\mathrm{1}\mathit{<}\mathit{k}\mathit{<}\mathit{n}\mathrm{;}\\ & \mathrm{if}\mathit{k}\mathrm{=}\mathit{n}\mathit{.}\end{array}$(110)We provide the exhaustive list of values of I_k,n(k_γ) for n = 2,3,4 and 5 in Table 3. Remark: k_γ depends on the selected completeness γ and on the number of catalogues n. When we call I_k,n−k(k_γ), the k_γ to be used in the integral is always the k_γ computed for n catalogues. So I_1,n−k will no longer be equal to γ but to F_{χ2(n−k−1)}(k_γ). For example, we fix γ to 0.9973. Then for a n = 2 catalogues χ-match, k_γ ≈ 3.4 and I_1,2 = γ. But for a n = 3 catalogues χ-match, k_γ ≈ 4.0, I_1,3 = γ and I_{1,3−1 = 2} ≠ γ.

We call p(s | h_k) the marginalized probability of observing a Mahalanobis distance less than or equal to k_γ in a group of n sources knowing they are actually associated to k real sources. $\mathit{p}\mathrm{\left(}\mathit{s}\mathrm{\right|}{\mathit{h}}_{\mathit{k}}\mathrm{)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathit{k}\mathrm{=}\mathrm{1}\mathrm{;}\\ & \mathrm{if}\mathit{k}\mathit{>}\mathrm{1.}\end{array}$(111)We obtain this equality by replacing 2πx_g−1 by 2πx_g−1det(V_Σg−1 + V_Σg)/Ω in Eq. (108) and then applying Eq. (63). The factor 1/Ω comes from the normalisation of the Poisson distribution so it integrates to one over the common surface area of the cross-matched catalogues. For the particular case in which all sources are spurious, we logically find the summed terms in Eq. (74).

And the distribution (p.d.f.) associated to the probability p(x | h_k) of observing a given Mahalanobis distance x knowing h_k is simply given by the derivative of p(s | h_k)dx, so is proportional to dI_k,n(x)dx.

9. Simple Bayesian probabilities

In this section, we compute Bayesian probabilities which depend on the Mahalanobis distance only.

9.1. General formula

Given a set of n candidates from n distinct catalogues satisfying the candidate selection criterion, we know

• x, the Mahalanobis distance (Eq. (40)), or χ value, which is a real value;

• s, the result of the selection criterion x ≤ k_γ, that is a boolean always equals to true for the sets of associations we keep, so p(s) = 1;

• { h_i }, i ∈ [1,B_n], the set of B_n (Eq. (73)) hypotheses to be formulated for each set of association.

We then note h_k hypotheses in which the n sources are associated with k “real” sources.

For a given set of n candidates from n distinct catalogues, the probabilities associated with the various hypotheses are given by Bayes’ formula $\begin{array}{ccc}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}\mathit{x,s}\mathrm{)}& \mathrm{=}& \frac{\mathit{p}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}\mathit{x}\mathrm{\left|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{\right)}}{\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{B}}_{\mathit{n}}}\mathit{p}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}\mathit{x}\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{\right)}}\mathit{,}\\ & \mathrm{=}& \end{array}$In this formula, p(h_i | s) are priors (considering only χ-matches, hence only s = true) and correspond to the number of associations satisfying the candidate selection (χ-matches) and hypotheses h_i over the total number of associations satisfying the candidate selection. We can transform the likelihood p(x | h_i,s) in $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathit{p}\mathrm{\left(}\mathit{x,}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{\right)}}{\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{\right)}}\mathit{,}\\ & \mathrm{=}& \frac{\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathrm{\left)}\mathit{p}\mathrm{\right(}\mathit{s}\mathrm{\left|}\mathit{x,}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}}{\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\mathit{p}\mathrm{\left(}\mathit{s}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathrm{)}}\mathit{,}\\ & \mathrm{=}& \frac{\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathrm{)}}{\mathit{p}\mathrm{\left(}\mathit{s}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathrm{)}}\mathit{,}\end{array}$because we keep in our sample only associations satisfying the candidate selection criteria we have p(s | x,h_i) = 1. In other words, the likelihood we use is a classical likelihood normalized so it integrates to one over the χ test acceptance region (defined by x ≤ k_γ).

We easily compute priors from the numbers estimated in Sect. 7. And likelihoods are simply computed from Sect. 8$\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\mathrm{d}\mathit{p}\mathrm{(}{\mathit{x}}^{\mathrm{\prime }}\mathit{<}\mathit{x}\mathrm{|}{\mathit{h}}_{\mathit{k}}\mathrm{)}\mathrm{d}\mathit{x}}{\mathit{p}\mathrm{\left(}\mathit{s}\mathrm{\right|}{\mathit{h}}_{\mathit{k}}\mathrm{)}}\mathrm{=}\frac{\mathrm{d}{\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathrm{·}$(117)We make explicit this result in the next section for the case of two, three and four catalogues.

9.2. Likelihoods p(x|h_i,s)

In this section, we compute the likelihoods p(x | h_i,s), that is the p.d.f. of the Mahalanobis distance of χ-matches under hypothesis h_i.

9.2.1. Case of two catalogues

For a set of two sources from two distinct catalogues, we have only two hypotheses, hence the two likelihoods: $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{1}}\mathrm{)}& \mathrm{=}& \mathit{\chi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}\mathrm{;}\\ \mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{2}}\mathrm{)}& \mathrm{=}& \frac{\mathrm{2}\mathit{\pi }\sqrt{\det \mathrm{(}{\mathit{V}}_{\mathrm{1}}\mathrm{+}{\mathit{V}}_{\mathrm{2}}\mathrm{)}}\mathit{x}\mathrm{d}\mathit{x}}{\mathrm{4}\mathit{\pi }}\mathrm{·}\end{array}$Knowing that the selection criterion is satisfied, we have to normalize so the integral of each likelihood over the domain defined by the selection criteria equals one (likelihoods are p.d.f.): $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{1}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathit{\chi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{k}}_{\mathit{\gamma }}}\mathit{\chi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}\mathrm{=}\frac{\mathit{\chi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{\mathit{\gamma }}\mathrm{=}\frac{\mathrm{d}{\mathit{I}}_{\mathrm{1}\mathit{,}\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathrm{1}\mathit{,}\mathrm{2}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathrm{;}\\ \mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{2}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathit{x}\mathrm{d}\mathit{x}}{{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{k}}_{\mathit{\gamma }}}\mathit{x}\mathrm{d}\mathit{x}}\mathrm{=}\frac{\mathrm{2}}{{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}}\mathit{x}\mathrm{d}\mathit{x}\mathrm{=}\frac{\mathrm{d}{\mathit{I}}_{\mathrm{2}\mathit{,}\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathrm{2}\mathit{,}\mathrm{2}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathrm{·}\end{array}$All constant terms, that is terms independent of x, vanish with the normalisation. The likelihoods are plotted in Fig. 2.

Remark: p(x | h_{k = 2}) mixes the derivative of the surface area of an ellipse in the Euclidean plane and the surface area of the sphere. It is an approximation valid as long as $\sqrt{\det \mathrm{(}{{V}}_{\mathrm{1}}\mathrm{+}{{V}}_{\mathrm{2}}\mathrm{)}}\mathit{x}$ is small enough so effects of the curvature of the sphere are negligible.

9.2.2. Case of three catalogues

For a set of three sources A, B, and C from three distinct catalogues, we have five hypotheses (see Sect. 6.2.2). In the five hypotheses, the number of “real” sources can be either one, two or three. We have as many likelihoods as possible distinct number of “real” sources. There are three ways of performing an iterative cross-match, leading to the same Mahalanobis distance x:

• cross-match A with B and then with C,${\mathit{x}}^{\mathrm{2}}\mathrm{=}{\mathit{x}}_{\mathit{AB}}^{\mathrm{2}}\mathrm{+}{\mathit{x}}_{\mathit{ABC}}^{\mathrm{2}}$;

• cross-match B with C and then with A, ${\mathit{x}}^{\mathrm{2}}\mathrm{=}{\mathit{x}}_{\mathit{BC}}^{\mathrm{2}}\mathrm{+}{\mathit{x}}_{\mathit{BCA}}^{\mathrm{2}}$;

• cross-match A with C and then with B, ${\mathit{x}}^{\mathrm{2}}\mathrm{=}{\mathit{x}}_{\mathit{AC}}^{\mathrm{2}}\mathrm{+}{\mathit{x}}_{\mathit{ACB}}^{\mathrm{2}}$.

We denote x_AB the Mahalanobis distance between A and B and we denote x_ABC the Mahalanobis distance between C and the weighted mean position of A and B.

x₁ and x₂ are used to designate without distinction x_AB, x_BC or x_AC and x_ABC, x_BCA or x_ACB respectively.

Although it may be tempting to write $\begin{array}{ccc}\mathit{p}\mathrm{\left(}{\mathit{x}}_{\mathit{AB}}\mathit{,}{\mathit{x}}_{\mathit{ABC}}\mathrm{\right|}{\mathit{h}}_{\mathit{ABC}}\mathrm{)}& \mathrm{=}& \frac{\mathit{\chi }\mathrm{\left(}{\mathit{x}}_{\mathit{AB}}\mathrm{\right)}\mathit{\chi }\mathrm{\left(}{\mathit{x}}_{\mathit{ABC}}\mathrm{\right)}\mathrm{d}{\mathit{x}}_{\mathit{AB}}\mathrm{d}{\mathit{x}}_{\mathit{ABC}}}{{\mathit{I}}_{\mathrm{1}\mathit{,}\mathrm{3}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathit{,}\\ \mathit{p}\mathrm{\left(}{\mathit{x}}_{\mathit{AB}}\mathit{,}{\mathit{x}}_{\mathit{ABC}}\mathrm{\right|}{\mathit{h}}_{\mathit{AB}\mathrm{\_}\mathit{C}}\mathrm{)}& \mathrm{=}& \frac{\mathit{\chi }\mathrm{\left(}{\mathit{x}}_{\mathit{AB}}\mathrm{\right)}\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{ABC}}\mathrm{d}{\mathit{x}}_{\mathit{AB}}\mathrm{d}{\mathit{x}}_{\mathit{ABC}}}{{\mathit{I}}_{\mathrm{2}\mathit{,}\mathrm{3}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathit{,}\\ \mathit{p}\mathrm{\left(}{\mathit{x}}_{\mathit{AC}}\mathit{,}{\mathit{x}}_{\mathit{ACB}}\mathrm{\right|}{\mathit{h}}_{\mathit{AC}\mathrm{\_}\mathit{B}}\mathrm{)}& \mathrm{=}& \frac{\mathit{\chi }\mathrm{\left(}{\mathit{x}}_{\mathit{AC}}\mathrm{\right)}\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{ACB}}\mathrm{d}{\mathit{x}}_{\mathit{AC}}\mathrm{d}{\mathit{x}}_{\mathit{ACB}}}{{\mathit{I}}_{\mathrm{2}\mathit{,}\mathrm{3}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathit{,}\\ \mathit{p}\mathrm{\left(}{\mathit{x}}_{\mathit{BC}}\mathit{,}{\mathit{x}}_{\mathit{BC}\mathrm{\_}\mathit{A}}\mathrm{\right|}{\mathit{h}}_{\mathit{A}\mathrm{\_}\mathit{BC}}\mathrm{)}& \mathrm{=}& \frac{\mathit{\chi }\mathrm{\left(}{\mathit{x}}_{\mathit{BC}}\mathrm{\right)}\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{BC}\mathrm{\_}\mathit{A}}\mathrm{d}{\mathit{x}}_{\mathit{BC}}\mathrm{d}{\mathit{x}}_{\mathit{BCA}}}{{\mathit{I}}_{\mathrm{2}\mathit{,}\mathrm{3}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathit{,}\\ \mathit{p}\mathrm{\left(}{\mathit{x}}_{\mathit{AB}}\mathit{,}{\mathit{x}}_{\mathit{ABC}}\mathrm{\right|}{\mathit{h}}_{\mathit{A}\mathrm{\_}\mathit{B}\mathrm{\_}\mathit{C}}\mathrm{)}& \mathrm{=}& \frac{\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{AB}}\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{ABC}}\mathrm{d}{\mathit{x}}_{\mathit{AB}}\mathrm{d}{\mathit{x}}_{\mathit{ABC}}}{{\mathit{I}}_{\mathrm{3}\mathit{,}\mathrm{3}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathit{,}\end{array}$we cannot directly compute probabilities p(h_i | x) from those likelihoods since infinitesimals (dx_AB, dx_AC, ...) are not the same and so do not vanish when applying Bayes’ formula.

It seems that the only measurement one can use to obtain coherent (and symmetrical) probabilities is the total Mahalanobis distance x. So we have to integrate the above probabilities over the domain defined by ${\mathit{x}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{x}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\le }{\mathit{x}}^{\mathrm{2}}$ and then evaluate their derivatives for x. We obtain the following likelihoods represented in Fig. 3. $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{1}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathrm{d}{\mathit{I}}_{\mathrm{1}\mathit{,}\mathrm{3}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathrm{1}\mathit{,}\mathrm{3}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{\chi }}_{\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{=}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{\mathit{\gamma }}\mathit{,}\\ \mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{2}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathrm{d}{\mathit{I}}_{\mathrm{2}\mathit{,}\mathrm{3}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathrm{2}\mathit{,}\mathrm{3}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{2}\mathit{x}\mathrm{(}\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}{\mathit{x}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{\left)}\mathrm{\right)}\mathrm{x}}{{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\mathrm{-}\mathrm{2}\mathrm{(}\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{\left)}\mathrm{\right)}}\mathit{,}\\ \mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{3}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathrm{d}{\mathit{I}}_{\mathrm{3}\mathit{,}\mathrm{3}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathrm{3}\mathit{,}\mathrm{3}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{4}{\mathit{x}}^{\mathrm{3}}\mathrm{x}}{{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{4}}}\mathit{,}\end{array}$in which h_{k = 1} is the hypothesis h_ABC, h_{k = 3} is the hypothesis h_{A_B_C} and h_{k = 2} is either the hypothesis h_{AB_C} or h_{AC_B} or h_{A_BC}.

9.2.3. Case of four catalogues

For a set of four sources A, B, C and D from four distinct catalogues, we have fifteen hypotheses (see Sect. 6.2.3). In the fifteen hypotheses, the number of “real” sources can be either one, two, three or four. We have as many likelihoods as possible distinct numbers of “real” sources. They are represented in Fig. 4: $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{1}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathrm{d}{\mathit{I}}_{\mathrm{1}\mathit{,}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathrm{1}\mathit{,}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{\chi }}_{\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{=}\mathrm{6}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{\mathit{\gamma }}\mathrm{;}\\ \mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{2}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathrm{d}{\mathit{I}}_{\mathrm{2}\mathit{,}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathrm{2}\mathit{,}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{;}\\ \mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{3}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathrm{d}{\mathit{I}}_{\mathrm{3}\mathit{,}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathrm{3}\mathit{,}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{;}\\ \mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{k}\mathrm{=}\mathrm{4}}\mathit{,s}\mathrm{)}& \mathrm{=}& \frac{\mathrm{d}{\mathit{I}}_{\mathrm{4}\mathit{,}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathrm{4}\mathit{,}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{6}{\mathit{x}}^{\mathrm{5}}\mathrm{d}\mathit{x}}{{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{6}}}\mathrm{·}\end{array}$In which h_{k = 1} is the hypothesis h_ABCD; h_{k = 4} is the hypothesis h_{A_B_C_D}; h_{k = 2} is either the hypothesis h_{A_BCD} or h_{B_ACD} or h_{C_ABD} or h_{D_ABC} or h_{AB_CD} or h_{AC_BD} or h_{AD_BC}; and h_{k = 3} is either the hypothesis h_{AB_C_D} or h_{AC_B_D} or h_{AD_B_C} or h_{BC_A_D} or h_{BD_A_C} or h_{CD_A_B}.

9.3. Advantage and limits

The main advantage of using p(h_i | x,s) is that the likelihoods it is based on do not depend on the positional errors: the only input parameter is the Mahalanobis distance x. Although it is true that x is computed from positional errors, once the χ-match has been performed we do not need the errors anymore: the distributions we use relies only on x. Changing positional errors modifies the priors, not the likelihoods. So we can easily add independent likelihoods based on magnitudes or other parameters.

There are two main problems. The first problem is precisely that the likelihoods depend only on x. It means that a set of very close sources with very accurate positions may have the same probability than a set of distant sources with large positional errors, even if intuitively the risk the first set have to contain spurious association should be far lower than in the second case. The second limitation is due to the fact that likelihoods are the same for hypotheses considering the same number of “real” sources. In the three catalogues case, p(x | h_{AB_C},s) = p(x | h_{A_BC},s) = p(x | h_{AC_B},s). The priors being constants, if p(h_{AB_C} | s) >p(h_{A_BC} | s) >p(h_{AC_B} | s), we always obtain posterior probability p(x | h_{AB_C},s) >p(h_{A_BC} | s) >p(h_{AC_B} | s).

10. Bayesian probabilities with positional errors

In this section, we compute Bayesian probabilities which include explicitly positional errors.

10.1. Warning about the non independence of positional uncertainties

In surveys providing individual uncertainties, positions of unsaturated bright sources are often more precise than positions of faint sources. The reason has to do with the higher photometric signal-to-noise ratio of bright sources compared to faint sources, while the FWHM is similar. An example of computation of positional uncertainties based on photon statistics can be found for example in the documentation of the SExtractor software (Bertin & Arnouts 1996). As mentioned in the documentation, the photon statistics based error is a lower value estimate.

It means that we cannot blindly assume that the positional uncertainties and photometric quantities like apparent magnitudes are independent. Moreover, if the positional errors of sources are related to their magnitudes and if the magnitudes of the sources in different catalogues are also related, then positional uncertainties in the different catalogues are related too. It means that we cannot blindly assume that the positional uncertainties of matching objects in different catalogues are independent from each other, at least not for h_k<n, that is the hypothesis in which at least two sources are from a same actual source.

One has to keep this in mind when using the naive independent hypothesis to simplify Bayes probabilities.

10.2. Probability using the Mahalanobis distance

To (at least partly) solve the first issue mentioned in Sect. 9.3 one possibility is to introduce likelihoods based for example on the volume V of the Chi test acceptance region writing: $\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}\mathit{x,V,s}\mathrm{)}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}\mathit{x}\mathrm{\left|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{\left(}\mathit{V}\mathrm{\right|}\mathit{x,}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{)}}{\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{B}}_{\mathit{n}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}\mathit{x}\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{\left(}\mathit{V}\mathrm{\right|}\mathit{x,}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{)}}\mathrm{·}$(134)From Sect. 7, it is easy – even though it may be time consuming – to build the estimated histogram n_{A_B}p(V + ΔV | h_{A_B}) (in which ΔV is the width of the histogram’s bars). Given this histogram and the result of a two-catalogues cross-match, we can also build an estimated histogram n_ABp(V + ΔV | h_AB). And so on for multiple catalogues, performing all possible sub-cross-matches.

If V and x are independent for all hypotheses, and knowing (having estimates of) n_AB and n_{A_B}, we have all the ingredients to compute $\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}\mathit{x,V}\mathrm{+}\mathrm{\Delta }\mathit{V,s}\mathrm{)}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}\mathit{x}\mathrm{\left|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{(}\mathit{V}\mathrm{+}\mathrm{\Delta }\mathit{V}\mathrm{|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{)}}{\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{B}}_{\mathit{n}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}\mathit{x}\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{(}\mathit{V}\mathrm{+}\mathrm{\Delta }\mathit{V}\mathrm{|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{)}}\mathit{,}$(135)even if it is not elegant to introduce a somewhat arbitrary slicing in V histograms.

10.3. Putting aside the Mahalanobis distance

We also consider the alternative form which puts aside the Mahalanobis distance and relies on the full sets of positions μ and associated errors V$\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}{\mu }\mathit{,}{V}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}{V}\mathrm{\left|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{\left(}{\mu }\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathit{,}{V}\mathit{,s}\mathrm{)}}{\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{B}}_{\mathit{n}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}{V}\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{\left(}{\mu }\mathrm{\right|}{\mathit{h}}_{\mathit{k}}\mathit{,}{V}\mathit{,s}\mathrm{)}}\mathit{,}$(136)in which the probabilities explicitly depend on the “configuration” of each position and on the associated errors. It also depends on the distribution of positional errors for a given hypothesis. Although p(V | h_i,s) can be estimated performing all possible sub-cross-matches, it is not trivial since it is a joint distribution in a space of dimension equal to the number of “actual” sources considered in h_i (using the circular error approximation).

10.3.1. Likelihoods p(μ |h_i, V, s)

We make the hypothesis that n_g sources are n_g detections of a same true source having a given position p. The probability to observe the set of positions μ_{g} = { μ₁,μ₂,...,μ_ng }, knowing p and the set of errors V_{g} = { V₁,V₂,...,V_ng } is $\mathit{p}\mathrm{\left(}{{\mu }}_{{\left\{}{g}{\right\}}}\mathrm{\right|}{p}\mathit{,}{{V}}_{{\left\{}{g}{\right\}}}\mathrm{)}\mathrm{=}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{{\mathit{n}}_{\mathit{g}}}{\mathrm{􏽙}}}{\mathrm{𝒩}}_{{{\mu }}_{{i}}\mathit{,}{\mathit{V}}_{\mathit{i}}}\mathrm{(}{p}\mathrm{)}\mathrm{d}{p}\mathrm{d}{{\mu }}_{{1}}\mathit{...}\mathrm{d}{{\mu }}_{{{n}}_{{g}}}\mathit{.}$(137)In practice we do not know the position of the real source p. So the probability to observe the set of positions μ_{g} knowing the set of errors V_{g} is obtained by integrating over all possible positions $\begin{array}{ccc}\mathit{p}\mathrm{\left(}{{\mu }}_{{\left\{}{g}{\right\}}}\mathrm{\right|}{{V}}_{{\left\{}{g}{\right\}}}\mathrm{)}& \mathrm{=}& \left(\mathrm{\int }\mathrm{\int }\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{{\mathit{n}}_{\mathit{g}}}{\mathrm{􏽙}}}{\mathrm{𝒩}}_{{{\mu }}_{{i}}\mathit{,}{\mathit{V}}_{\mathit{i}}}\mathrm{\left(}{p}\mathrm{\right)}\mathrm{d}{p}\right)\mathrm{d}{{\mu }}_{{1}}\mathit{...}\mathrm{d}{{\mu }}_{{{n}}_{{g}}}\mathit{,}\\ & \mathrm{=}& \sqrt{\frac{\det {{V}}_{{{\Sigma }}_{{\left\{}{g}{\right\}}}}}{\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{{\mathit{n}}_{\mathit{g}}}{\mathrm{􏽑}}}\det {{V}}_{{i}}}}\frac{\exp \mathrm{\{}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathit{,}\mathit{...}\mathit{,}{{\mu }}_{{{n}}_{{g}}}\mathrm{\right)}}\mathrm{\}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{{\mathit{n}}_{\mathit{g}}\mathrm{-}\mathrm{1}}}\\ & & \mathrm{\times }\mathrm{d}{{\mu }}_{{1}}\mathit{...}\mathrm{d}{{\mu }}_{{{n}}_{{g}}}\mathit{.}\end{array}$This result is the same as Eq. (38) in Sect. 5.2 but applied here to the sub-set of positions { μ₁,μ₂,...,μ_ng }. The difference is that in Sect. 5.2 we wanted to estimate the probability the sources had to be at the same location whereas here, knowing (making the hypothesis) they are at the same location, we compute the probability we had to observe this particular outcome. The particular case n_g = 1 leads to p(μ_{g} | V_{g}) = dμ₁.

We now consider G groups and the selection criteria s (x ≤ h_γ). Each of the n input sources is part of one, and only one group. Given the G groups, the errors on the positions and the candidate selection criteria, the probability to observe positions μ is $\begin{array}{ccc}\mathit{p}\mathrm{\left(}{\mu }\mathrm{\right|}{\mathit{h}}_{\mathit{G}}\mathit{,}{V}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}{\mu }\mathrm{\right|}{\mathit{h}}_{\mathit{G}}\mathit{,}{V}\mathrm{)}}{\mathit{p}\mathrm{\left(}\mathit{s}\mathrm{\right|}{\mathit{h}}_{\mathit{G}}\mathit{,}{V}\mathrm{)}}\mathrm{=}\frac{\underset{\mathit{g}\mathrm{=}\mathrm{1}}{\overset{\mathit{G}}{\mathrm{􏽑}}}\mathit{p}\mathrm{\left(}{{\mu }}_{{\left\{}{g}{\right\}}}\mathrm{\right|}{{V}}_{{\left\{}{g}{\right\}}}\mathrm{)}}{{\mathrm{\int }}_{\mathit{x}\mathrm{\le }{\mathit{k}}_{\mathit{\gamma }}}\underset{\mathit{g}\mathrm{=}\mathrm{1}}{\overset{\mathit{G}}{\mathrm{􏽑}}}\mathit{p}\mathrm{\left(}{{\mu }}_{{\left\{}{g}{\right\}}}\mathrm{\right|}{{V}}_{{\left\{}{g}{\right\}}}\mathrm{)}}\mathrm{·}& & \end{array}$(140)The denominator ensures that the likelihood integrates to one over its domain of definition, domain delimited by the candidate selection criteria, that is the region of acceptance of the χ² test.

Let us compute the integral in the denominator. The differential of the substitution x = yV^-1y transforms as $\mathrm{d}{y}\mathrm{=}\sqrt{\det {V}}\mathit{x}\mathrm{d}\mathit{x}\mathrm{d}\mathit{\theta }$. y can be the difference between two positions (e.g. μ_i−μ_Σi−1) and (x,θ) the polar coordinates of y in the basis defined by the eigenvectors of V and reduced by its eigenvalues. V can for example be V_Σi−1 + V_i. Using the iterative form of Sect. 5.3 and Eq. (53), we can rewrite p(μ_{g} | V_{g})$\mathit{p}\mathrm{\left(}{{\mu }}_{{\left\{}{g}{\right\}}}\mathrm{\right|}{{V}}_{{\left\{}{g}{\right\}}}\mathrm{)}\mathrm{=}\frac{\exp \left\{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{g}}\mathrm{-}\mathrm{1}}{\mathit{x}}_{\mathit{i}}\right\}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{{\mathit{n}}_{\mathit{g}}\mathrm{-}\mathrm{1}}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{{\mathit{n}}_{\mathit{g}}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{x}}_{\mathit{i}}\mathrm{d}{\mathit{x}}_{\mathit{i}}\mathrm{d}{\mathit{\theta }}_{\mathit{i}}\mathit{.}$(141)Integrating over all θ_i we obtain $\begin{array}{ccc}\mathit{p}\mathrm{\left(}{{x}}_{{\left\{}{g}{\right\}}}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathrm{2}\mathit{\pi }\\ \mathrm{0}\end{array}\mathit{...}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathit{p}\mathrm{\left(}{{\mu }}_{{\left\{}{g}{\right\}}}\mathrm{\right|}{{V}}_{{\left\{}{g}{\right\}}}\mathrm{)}\mathit{,}\\ & \mathrm{=}& \exp \left\{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{g}}\mathrm{-}\mathrm{1}}{\mathit{x}}_{\mathit{i}}\right\}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{{\mathit{n}}_{\mathit{g}}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{x}}_{\mathit{i}}\mathrm{d}{\mathit{x}}_{\mathit{i}}\mathit{,}\\ & \mathrm{=}& \underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{{\mathit{n}}_{\mathit{g}}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\mathrm{d}{\mathit{x}}_{\mathit{i}}\mathit{.}\end{array}$This joint p.d.f. of the successive Mahalanobis distances is different from the p.d.f. of their quadratic sum which gives the total Mahalanobis distance $\mathit{p}\mathrm{\left(}{\mathit{x}}_{\mathrm{\left\{}\mathit{g}\mathrm{\right\}}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{\mathrm{0}}^{\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{p}\mathrm{\left(}{{x}}_{{\left\{}{g}{\right\}}}\mathrm{\right)}\mathrm{=}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathrm{(}{\mathit{n}}_{\mathit{g}}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{\left(}{\mathit{x}}_{\mathrm{\left\{}\mathit{g}\mathrm{\right\}}}\mathrm{\right)}\mathrm{d}{\mathit{x}}_{\mathrm{\left\{}\mathit{g}\mathrm{\right\}}}\mathit{,}$(145)in which ${\mathit{x}}_{\mathrm{\left\{}\mathit{g}\mathrm{\right\}}}^{\mathrm{2}}\mathrm{=}\mathrm{\left|}\mathrm{\right|}{{x}}_{{\left\{}{g}{\right\}}}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{g}}}{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}$.

Putting all together, the integral in the denominator is no other than the integral I_{k = G,n}(k_γ) defined in Sect. 8. Written explicitly: $\mathit{p}\mathrm{\left(}{\mu }\mathrm{\right|}{\mathit{h}}_{\mathit{G}}\mathit{,}{V}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\underset{\mathit{g}\mathrm{=}\mathrm{1}}{\overset{\mathit{G}}{\mathrm{􏽑}}}\mathit{p}\mathrm{\left(}{{\mu }}_{{\left\{}{g}{\right\}}}\mathrm{\right|}{{V}}_{{\left\{}{g}{\right\}}}\mathrm{)}}{{\mathit{I}}_{\mathit{k}\mathrm{=}\mathit{G,n}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathrm{·}$(146)

10.3.2. Classical two catalogues case

In the case of two catalogues, the probabilities are simply:

$\begin{array}{ccc}& & \mathit{p}\mathrm{\left(}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathrm{\right|}{\mathit{h}}_{\mathrm{1}}\mathit{,}{{V}}_{{1}}\mathit{,}{{V}}_{{2}}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\exp \mathrm{\{}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}\mathrm{\}}\mathrm{d}{{\mu }}_{{1}}\mathrm{d}{{\mu }}_{{2}}}{\mathrm{2}\mathit{\pi }\sqrt{\det \mathrm{(}{{V}}_{{1}}\mathrm{+}{{V}}_{{2}}\mathrm{)}}\mathit{\gamma }}\mathrm{;}\\ & & \mathit{p}\mathrm{\left(}{{\mu }}_{{1}}\mathit{,}{{\mu }}_{{2}}\mathrm{\right|}{\mathit{h}}_{\mathrm{2}}\mathit{,}{{V}}_{{1}}\mathit{,}{{V}}_{{2}}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\mathrm{d}{{\mu }}_{{1}}\mathrm{d}{{\mu }}_{{2}}}{\mathit{\pi }{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}}\mathrm{·}\end{array}$

If we compare the likelihood ratio LR_μ,V computed from those formulae with the likelihood ratio LR_x computed from the previous result (Sect. 9.2.1) we obtain: $\begin{array}{ccc}\mathit{L}{\mathit{R}}_{\mathit{x}}& \mathrm{=}& \frac{{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}{\mathit{x}}^{\mathrm{2}}\mathit{/}\mathrm{2}}}{\mathrm{2}\mathit{\gamma }}\mathrm{;}\\ \mathit{L}{\mathit{R}}_{{\mu }\mathit{,}{V}}& \mathrm{=}& \frac{{\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}{\mathit{x}}^{\mathrm{2}}\mathit{/}\mathrm{2}}}{\mathrm{2}\mathit{\gamma }\sqrt{\det \mathrm{(}{{V}}_{\mathrm{1}}\mathrm{+}{{V}}_{\mathrm{2}}\mathrm{)}}}\mathrm{·}\end{array}$Contrary to LR_x, LR_μ,V accounts for the size of positional errors. As we will see in the next section, the drawback is that we can hardly combine the likelihoods p(μ₁,μ₂ | h_k,V₁,V₂,s) with photometry based likelihoods.

11. Bayesian probabilities with photometric data

All probabilities of association discussed so far are based on the likelihood that the positions recorded in various catalogues are consistent with that of a unique astrophysical object. However, one may wish to make additional assumptions on the nature of the source (e.g. star, active galactic nucleus, etc.) that could help decrease or, conversely, increase the plausibility of a given association of catalogue entries. This is particularly important when one seeks to gather homogeneous samples of objects. Spectral energy distributions assembled from photometric catalogues can be usefully compared with templates and assigned a probability of being representative of the targeted class of objects. This procedure has been presented in Budavári & Szalay (2008) and recently used in Naylor et al. (2013) and Hsu et al. (2014). However, following Budavári & Szalay (2008) we underline that entering criteria of resemblance to a given class of objects in the computation of association probabilities is done at the expense of the capability to find scientifically interesting outliers.

Another possibility may consist in building colour-colour diagrams for random and χ-matched associations to derive colour-colour diagrams for real associations (more precisely, we have to derive each diagram for each possible hypothesis h_i). Those normalized diagrams are p.d.f. that can be interpreted as likelihoods (p(m | h_i,s) in the following equations). Smoothing those diagrams, one can see them as the likelihoods of the kernel density classification (Richards et al. 2004), replacing the object types by the hypothesis h_i, and using magnitudes from different catalogues instead of just one.

We detail below how photometric data can be folded into the output of the purely astrometric method discussed in this paper, without making additional assumptions on the nature of the source.

Suppose we note a, a vector containing all the astrometric information we have about a set of n candidates (a may contain the positions, the associated covariance matrices, ...). We note m the set of photometric informations we have about the same set of n candidates (m may contain magnitudes and/or colours, associated errors, ...). Then we can write the Bayes formula: $\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}{a}\mathit{,}{m}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}{a}\mathrm{\left|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{\left(}{m}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathit{,s,}{a}\mathrm{)}}{\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{B}}_{\mathit{n}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}{a}\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{\left(}{m}\mathrm{\right|}{\mathit{h}}_{\mathit{k}}\mathit{,s,}{a}\mathrm{)}}\mathrm{·}$(151)If a and m are independent (naive hypothesis which is not granted, see Sect. 10.1) $\mathit{p}\mathrm{\left(}{m}\mathrm{\right|}{\mathit{h}}_{\mathit{k}}\mathit{,}{a}\mathit{,s}\mathrm{)}\mathrm{=}\mathit{p}\mathrm{(}{m}\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{\right)}\mathit{,}$(152)and Eq. (151) becomes $\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}{a}\mathit{,}{m}\mathrm{)}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}{a}\mathrm{\left|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{\left(}{m}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{)}}{\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{B}}_{\mathit{n}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}\mathit{s}\mathrm{\left)}\mathit{p}\mathrm{\right(}{a}\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{\right)}\mathit{p}\mathrm{\left(}{m}\mathrm{\right|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{)}}\mathrm{·}$(153)Let’s imagine we perform a cross-match taking into account astrometric data only. We compute probabilities p(h_i | a,s) for all possible hypotheses. Then if a and m are independent, and if we are able to compute likelihoods based on photometric data only p(m | h_i,s), then we can compute the probabilities p(h_i | a,m,s) in a second step from the probabilities computed in the astrometric part: $\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}{a}\mathit{,}{m}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right|}{a}\mathit{,s}\mathrm{\left)}\mathit{p}\mathrm{\right(}{m}\mathrm{\left|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{\right)}}{\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{B}}_{\mathit{n}}}\mathit{p}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}{a}\mathit{,s}\mathrm{\left)}\mathit{p}\mathrm{\right(}{m}\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathit{,s}\mathrm{\right)}}\mathit{,}$(154)which is equivalent to Eq. (153).

Unfortunately, positional errors and magnitudes are not necessarily independent. So one should not use p(h_i | μ,V,s) without any due caution in Eq. (154). However, one can use the Mahalanobis distance x which is independent of the photometry, that is probabilities p(h_i | x,s) (Eq. (113)).

12. Tests on synthetic catalogues

In the context of the ARCHES project, we developed a tool implementing the statistical multi-catalogue cross-match described in this paper. We added to the tool the possibility to generate synthetic catalogues that can be cross-matched like real tables. It has been allowing us to perform tests and to check both the software and the theory.

We present here such a test and provide the associated script (see Appendix C) so anybody can try it independently, possibly changing the input values. Currently the tool is accessible both via a web interface and an HTTP API⁹. Future plans are discussed in the conclusion Sect. 14.

We generate three synthetic catalogues, setting the numbers of sources they contain and have in common: we call n_ABC the number of common sources in the three catalogues A, B and C; n_AB the number of common sources in A and B only; n_A the number of sources in catalogue A only; and so on. Knowing a priori common and distinct sources in the catalogues, we can track the associations which are real and the spurious ones in the cross-match output. We can also check for missing associations.

The error associated to each individual position is a random value which follows a user define distribution. A different distribution is used for each catalogue. For catalogue A, we choose a constant value equal to 0.4′′; for catalogue B, the positional error distribution follows a linear function between 0.8′′ and 1.2′′; for catalogue C, the positional errors follow a Gaussian distribution of mean 0.75′′ and standard deviation 0.1′′ truncated to the 0.5–1′′ range.

We set the input sky area to be a cone of radius 0.42 degrees. Each position is randomly (uniform distribution) placed in this cone. For each catalogue in which the source is included, we randomly pick a positional error that we associate to the source and we blur the position using its error.

We first compute $\overline{{{V}}_{{\Sigma }}}$ for pairs AB, AC and BC. Given the chosen error distributions, the mean errors are equal to the median errors and the mean of the inverse of the errors is quite close to the inverse of the mean errors. So, for this particular case, we use the inverse of the mean errors 0.4, 1 and 0.75 instead of the means of the inverse. Given this approximation and using Eq. (96), we obtain $\begin{array}{ccc}\overline{\sqrt{\begin{array}{c}\det \end{array}{{V}}_{{{\Sigma }}_{{AB}}}}}& \mathrm{=}& \frac{\mathrm{1}}{\frac{\mathrm{1}}{\overline{\sqrt{\det {{V}}_{{A}}}}}\mathrm{+}\frac{\mathrm{1}}{\overline{\sqrt{\det {{V}}_{{B}}}}}}\mathit{,}\\ & \mathrm{=}& \frac{{\mathrm{0.4}}^{\mathrm{2}}\mathrm{\times }{\mathrm{1.0}}^{\mathrm{2}}}{{\mathrm{0.4}}^{\mathrm{2}}\mathrm{+}{\mathrm{1.0}}^{\mathrm{2}}}\mathrm{=}\mathrm{0.138}\mathrm{;}\\ \overline{\sqrt{\begin{array}{c}\det \end{array}{{V}}_{{{\Sigma }}_{{AC}}}}}& \mathrm{=}& \frac{\mathrm{1}}{\frac{\mathrm{1}}{\overline{\sqrt{\det {{V}}_{{A}}}}}\mathrm{+}\frac{\mathrm{1}}{\overline{\sqrt{\det {{V}}_{{C}}}}}}\mathit{,}\\ & \mathrm{=}& \frac{{\mathrm{0.4}}^{\mathrm{2}}\mathrm{\times }{\mathrm{0.75}}^{\mathrm{2}}}{{\mathrm{0.4}}^{\mathrm{2}}\mathrm{+}{\mathrm{0.75}}^{\mathrm{2}}}\mathrm{=}\mathrm{0.125}\mathrm{;}\\ \overline{\sqrt{\begin{array}{c}\det \end{array}{{V}}_{{{\Sigma }}_{\mathit{BC}}}}}& \mathrm{=}& \frac{\mathrm{1}}{\frac{\mathrm{1}}{\overline{\sqrt{\det {{V}}_{{B}}}}}\mathrm{+}\frac{\mathrm{1}}{\overline{\sqrt{\det {{V}}_{{C}}}}}}\mathit{,}\\ & \mathrm{=}& \frac{{\mathrm{1.0}}^{\mathrm{2}}\mathrm{\times }{\mathrm{0.75}}^{\mathrm{2}}}{{\mathrm{1.0}}^{\mathrm{2}}\mathrm{+}{\mathrm{0.75}}^{\mathrm{2}}}\mathrm{=}\mathrm{0.36.}\end{array}$And, to estimate the number of “fully” spurious associations, we have to compute the mean of the square root of Eq. (97) over all possible source trios, which can be approximated in this specific case by $\begin{array}{ccc}\overline{{\left(\frac{\det {{V}}_{\mathit{A}}\det {{V}}_{\mathit{B}}\det {{V}}_{{C}}}{\det {{V}}_{{{\Sigma }}_{\mathit{ABC}}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}}& \mathrm{\approx }& {\mathrm{1.0}}^{\mathrm{2}}\mathrm{\times }{\mathrm{0.75}}^{\mathrm{2}}\mathrm{+}{\mathrm{0.4}}^{\mathrm{2}}\mathrm{\times }{\mathrm{0.75}}^{\mathrm{2}}\\ & & \mathrm{+}{\mathrm{0.4}}^{\mathrm{2}}\mathrm{\times }{\mathrm{1.0}}^{\mathrm{2}}\mathit{,}\\ & \mathrm{\approx }& \mathrm{0.8125.}\end{array}$We note that in this particular case, the error distribution of sources A involved in AB, AC and ABC associations is the same. Idem for the error distribution of B and C sources.

We are now able to compute all components, depending on the size of histograms bins (step). To do this we note $\begin{array}{ccc}{\mathit{n}}_{\mathit{AB}\mathrm{\ast }}& \mathrm{=}& {\mathit{n}}_{\mathit{ABC}}\mathrm{+}{\mathit{n}}_{\mathit{AB}}\mathit{,}\\ {\mathit{n}}_{\mathit{AC}\mathrm{\ast }}& \mathrm{=}& {\mathit{n}}_{\mathit{ABC}}\mathrm{+}{\mathit{n}}_{\mathit{AC}}\mathit{,}\\ {\mathit{n}}_{\mathit{BC}\mathrm{\ast }}& \mathrm{=}& {\mathit{n}}_{\mathit{ABC}}\mathrm{+}{\mathit{n}}_{\mathit{BC}}\mathit{,}\\ {\mathit{n}}_{\mathit{A}\mathrm{\ast }}& \mathrm{=}& {\mathit{n}}_{\mathit{ABC}}\mathrm{+}{\mathit{n}}_{\mathit{AB}}\mathrm{+}{\mathit{n}}_{\mathit{AC}}\mathrm{+}{\mathit{n}}_{\mathit{A}}\mathit{,}\\ {\mathit{n}}_{\mathit{B}\mathrm{\ast }}& \mathrm{=}& {\mathit{n}}_{\mathit{ABC}}\mathrm{+}{\mathit{n}}_{\mathit{AB}}\mathrm{+}{\mathit{n}}_{\mathit{BC}}\mathrm{+}{\mathit{n}}_{\mathit{B}}\mathit{,}\\ {\mathit{n}}_{\mathit{C}\mathrm{\ast }}& \mathrm{=}& {\mathit{n}}_{\mathit{ABC}}\mathrm{+}{\mathit{n}}_{\mathit{AC}}\mathrm{+}{\mathit{n}}_{\mathit{BC}}\mathrm{+}{\mathit{n}}_{\mathit{C}}\mathit{,}\end{array}$to finally obtain $\begin{array}{ccc}\mathit{n̂}\mathit{ABC}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \mathit{step}\mathrm{\times }{\mathit{n}}_{\mathit{ABC}}\mathrm{\times }{\mathit{\chi }}_{\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{=}\mathrm{4}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}\\ \mathit{n̂}\mathit{AB}\mathrm{\_}\mathit{C}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{\approx }& \mathit{step}\mathrm{\times }{\mathit{n}}_{\mathit{AB}\mathrm{\ast }}\mathrm{\times }{\mathit{n}}_{\mathit{C}\mathrm{\ast }}\frac{\mathrm{0.138}\mathrm{+}{\mathrm{0.75}}^{\mathrm{2}}}{\mathit{\pi }\mathrm{(}\mathrm{0.42}\mathrm{\times }\mathrm{3600}{\mathrm{)}}^{\mathrm{2}}}\\ & & \mathrm{\times }\mathrm{2}\mathit{\pi x}\left(\mathrm{1}\mathrm{-}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}\right)\right)\mathit{,}\\ \mathit{n̂}\mathit{AC}\mathrm{\_}\mathit{B}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{\approx }& \mathit{step}\mathrm{\times }{\mathit{n}}_{\mathit{AC}\mathrm{\ast }}\mathrm{\times }{\mathit{n}}_{\mathit{B}\mathrm{\ast }}\frac{\mathrm{0.125}\mathrm{+}{\mathrm{1.0}}^{\mathrm{2}}}{\mathit{\pi }\mathrm{(}\mathrm{0.42}\mathrm{\times }\mathrm{3600}{\mathrm{)}}^{\mathrm{2}}}\\ & & \mathrm{\times }\mathrm{2}\mathit{\pi x}\left(\mathrm{1}\mathrm{-}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}\right)\right)\mathit{,}\\ \mathit{n̂}\mathit{BC}\mathrm{\_}\mathit{A}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{\approx }& \mathit{step}\mathrm{\times }{\mathit{n}}_{\mathit{BC}\mathrm{\ast }}\mathrm{\times }{\mathit{n}}_{\mathit{A}\mathrm{\ast }}\frac{\mathrm{0.36}\mathrm{+}{\mathrm{0.4}}^{\mathrm{2}}}{\mathit{\pi }\mathrm{(}\mathrm{0.42}\mathrm{\times }\mathrm{3600}{\mathrm{)}}^{\mathrm{2}}}\\ & & \mathrm{\times }\mathrm{2}\mathit{\pi x}\left(\mathrm{1}\mathrm{-}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}\right)\right)\mathit{,}\\ \mathit{n̂}\mathit{A}\mathrm{\_}\mathit{B}\mathrm{\_}\mathit{C}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{\approx }& \mathit{step}\mathrm{\times }{\mathit{n}}_{\mathit{A}\mathrm{\ast }}{\mathit{n}}_{\mathit{B}\mathrm{\ast }}\mathrm{\times }{\mathit{n}}_{\mathit{C}\mathrm{\ast }}\frac{\mathrm{0.8125}}{{\mathit{\pi }}^{\mathrm{2}}\mathrm{(}\mathrm{0.42}\mathrm{\times }\mathrm{3600}{\mathrm{)}}^{\mathrm{4}}}\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}{\mathit{x}}^{\mathrm{3}}\mathit{,}\\ \mathit{n̂}\mathit{Tot}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{\approx }& \mathit{n̂}\mathit{ABC}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{n̂}\mathit{A}\mathrm{\_}\mathit{B}\mathrm{\_}\mathit{C}\mathrm{\left(}\mathit{x}\mathrm{\right)}\\ & & \mathrm{+}\mathit{n̂}\mathit{AB}\mathrm{\_}\mathit{C}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{n̂}\mathit{AC}\mathrm{\_}\mathit{B}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{n̂}\mathit{BC}\mathrm{\_}\mathit{A}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{.}\end{array}$Each component equals step × n × p(s | h_i)p(x | h_i,s), see Eqs. (111) and (117). We simplify the expression since p(s | h_i) ∝ I_k,n(k_γ) and p(x | h_i,s) ∝ 1 /I_k,n(k_γ). The normalized histograms associated to each component are distributed according to each likelihood p(x | h_i,s). Both histograms made from the data and theoretical curves are plotted in Fig. 5. The theoretical results fit very well the result of the Chi-square cross-match based on simulated data.

We also verify that the number of “good” ABC matches we obtain as output of the cross-match is coherent with n_ABC times the input completeness γ.

When cross-matching real catalogues, the number of sources n_ABC, etc. are not known. But the previous “theoretical” curves can be built after a χ-match from the number of sources estimated to compute priors in Sect. 7.

13. Summarized recipe

In this section, we give the main steps and equations to perform a χ-match and to compute for each association the probability it has to be a good match (or any other possible hypothesis).

For a small and compact sky area, project all the sources of all catalogues on an Euclidian plane using for example the ARC projection (Calabretta & Greisen 2002).

To select matching candidates, for each possible set of n sources from n distinct catalogues:

• compute their weighted mean position (Eq. (30)) andthe associated error (Eq. (29));

• derive x, the Mahalanobis distance defined by the square root of Eq. (36);

• fix a constant threshold on all Mahalanobis distances, that is

  • set the fraction α of real associations it is acceptable to miss – the type I error – and

  • derive numerically the threshold k_γ inverting Eq. (46) based on the Chi-square distribution with 2(n−1) degrees of freedom (the result is the same computing the threshold from the Chi distribution with 2(n−1) degrees of freedom);

keep all sets of n sources having a Mahalanobis distance less than the threshold k_γ (Eq. (47)) as possibly being n observations of a same real source.

To compute Bayes’ probabilities, as many hypotheses as the number of possible partitions of the set of n sources (see Eq. (73), Table 2 and Fig. 1) have to be formulated. Depending on whether one wants to be able to account for photometry in a second step or not, a set of likelihoods may be chosen among several such sets.

In the first case, the likelihood associated to each hypothesis (knowing the selection criteria is fulfilled) depends only on the Mahalanobis distance and on the number of real sources k in the hypothesis h_i (see e.g. Fig. 4). The likelihoods are (Eq. (117)) $\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right|}{\mathit{h}}_{\mathit{i}}\mathit{,s}\mathrm{)}\mathrm{=}\frac{\mathrm{d}{\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}{{\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}{\mathit{k}}_{\mathit{\gamma }}\mathrm{\right)}}\mathit{,}$(175)with I_k,n(k_γ) given in Eq. (110). The formulae of I_k,n(k_γ) and dI_k,n(x) are provided for n ≤ 5 in Tables 3 and 4 respectively.

In the second case (no photometry to be taken into account), one can use likelihoods defined by Eq. (146).

Finally, to apply Bayes’ formula, priors p(h_i | s) are needed. This is more tricky and the steps detailed in Sect. 7 have to be considered. A pre-requisite is to work on an area (Ω) uniformly covered by all catalogues. For two catalogues, the number of spurious associations can be estimated by computing for each catalogue the mean area covered by the error ellipses for a radius equal to the threshold k_γ (so the mean area of the 1σ error ellipses times ${\mathit{k}}_{\mathit{\gamma }}^{\mathrm{2}}$). The two means are summed and the result is divided by Ω to obtain the mean probability to spuriously associate two unrelated sources. It is then multiplied by the product of the number of sources in both catalogues to finally obtain the mean expected number of spurious associations (Eq. (85)). Knowing the number of associations in the cross-correlation output, both the probability that one such association is spurious and the complementary probability of having a real association (the two priors of the two-catalogues cross-match) can be estimated (see Eqs. (93) and (94)). Similarly, performing all possible sub-χ-matches, all priors needed for a n-catalogues cross-match can be derived.

All needed ingredients to compute the probabilities associated with each hypothesis are thus available. Those probabilities can be computed applying Eq. (113).

14. Conclusions

In this paper we developed a comprehensive framework for performing the cross-correlation of multiple astronomical catalogues, in one pass. The approach employs a classical χ²-test to select candidates. We computed two sets of likelihoods based on positions, individual elliptical positional errors and the χ²-test region of acceptance: one that can be mixed without any caution with other parameters such as photometric values; and one for which the naive hypothesis of independence between positional uncertainties and magnitudes has to be tested. We also presented a way to estimate “priors” from the region of acceptance of the χ²-test. Probabilities for each possible hypothesis can thus be computed from those likelihoods and “priors”.

In practice the number of hypotheses, and thus the number of “priors”, increases dramatically with the number of catalogues. To be able to cross-match more than six or seven catalogues, it is necessary to simplify the problem. One possibility consists of merging two catalogues of similar astrometric accuracy and similar wavelength range, considering all matches as non-spurious matches. Doing so we would effectively reduce the number of input catalogues by one.

A large part of the statistical work carried out here depends on the simplifying assumptions made in Sect. 3: perfect astrometrical calibration (no systematic offsets), no proper motions, no clustering and no blending. In real life, the “normalized” distance between two detections of a same source present in two distinct catalogues hardly follows a Rayleigh distribution. The “actual distribution” (in practice it is not easy to build such a distribution since it requires secure identifications) often has a broader tail (see for example Rosen et al. 2016, Fig. 5) and a log-normal distribution may better fit it than the Rayleigh distribution. This is probably due to a combination of causes like small proper motions, imperfect reduction, systematics or bias from the calibration process, under or overestimated errors, etc.

In practice this means that the number of associations missed by the candidate selection criteria (based on Rayleigh) is larger that the chosen theoretical value (γ). We could for example add larger systematics to positional errors. The risk is then to distort (even more) the Rayleigh distribution. We could also try to re-calibrate locally the set of catalogues we want to cross-match, but we need secure identifications to do it properly; for each catalogue, all sources in the local area must have been calibrated at once (to possibly correct for a locally uniform systematic using four simple parameters Δα, Δδ, scale, θ). Those two constraints (having secure identifications and at once calibration) are in practice quite hard to satisfy.

If we re-calibrate using a “secure” population (i.e. a population of objects having no proper motions like QSOs) we introduce a bias since QSOs are fainter than most stars in the optical and thus have errors larger than the global population of objects. And adding stars, we introduce noise due to proper motions.

For these reasons, we believe that in case of “old” optical surveys based on photographic plates, a classical fixed radius cross-match may be more efficient that the χ-match to select candidates. We are nonetheless conviced that the equations we derived in this paper can help in building new catalogues, based for example on both multi-band and multi-epoch observations, and can be used to assess and improve the quality of coming surveys.

We generated and processed synthetical catalogues, which meet the simplifying assumptions, in the tool we developed for the ARCHES project. The consistency between the theoretical results derived in this paper – completeness of the candidate selection criterion, likelihoods and priors – and the outputs of the tool has allowed us to cross-validate both the method and its implementation. The tool has also been used to generate ARCHES products which were used in the scientific work packages of the project. Currently the CDS XMatch Service (Pineau et al. 2011a; Boch et al. 2012; Pineau et al. 2015) provides a basic but very efficient facility to cross-correlate two possibly large (>1 billion sources) catalogues. It is planned to include the ARCHES tool into the CDS XMatch. This paper will be the basic reference for the extension of the latter to multi-catalogue statistical χ-match.

---

Acknowledgments

A large part of this work was supported by the ARCHES project. ARCHES (No. 313146) was funded by the 7th Framework of the European Union and coordinated by the University of Strasbourg. All figures (except Fig. 1) were made using the ctioga2 plotting program developped by Vincent Fourmond. F. J. Carrera also acknowledges financial support through grant AYA2015-64346-C2-1-P (MINECO/FEDER).

References

Appendix A: Demonstrations

Appendix A.1: From V$\begin{array}{c}-1\\ {\Sigma }\end{array}$ to V_Σ

From Eq. (28), we compute for 2 × 2 symmetric square matrices: $\begin{array}{ccc}{\mathit{V}}_{\mathrm{\Sigma }}& \mathrm{=}& \mathrm{(}{{V}}_{{\Sigma }}^{-1}{\mathrm{)}}^{-1}\mathit{,}\\ & \mathrm{=}& \frac{\mathrm{1}}{\det {{V}}_{{\Sigma }}^{-1}}\mathrm{adj}\mathrm{\left(}{{V}}_{{\Sigma }}^{-1}\mathrm{\right)}\mathit{,}\\ & \mathrm{=}& \frac{\mathrm{1}}{\det {{V}}_{{\Sigma }}^{-1}}\mathrm{adj}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{V}}_{{i}}^{-1}\mathit{,}\\ & \mathrm{=}& \frac{\mathrm{1}}{\det {{V}}_{{\Sigma }}^{-1}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{adj}{{V}}_{{i}}^{-1}\mathit{,}\\ & \mathrm{=}& \frac{\mathrm{1}}{\det {{V}}_{{\Sigma }}^{-1}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{adj}\frac{\mathrm{adj}{{V}}_{{i}}}{\det {{V}}_{{i}}}\mathit{,}\\ & \mathrm{=}& \frac{\mathrm{1}}{\det {{V}}_{{\Sigma }}^{-1}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{{V}}_{{i}}}{\det {{V}}_{{i}}}\mathit{,}\end{array}$in which adjA is the adjugate matrix of A, that is the transpose of the cofactor matrix of A.

Appendix A.2: Sum of quadratics canonical form

First expanding and then factoring: $\begin{array}{ccc}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{Q}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{(}{x}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}}^{⊺}{{V}}_{\mathit{i}}^{-1}\mathrm{(}{x}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}\mathit{,}\\ & \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{x}}^{⊺}{{V}}_{{i}}^{-1}{x}\mathrm{-}\mathrm{2}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{x}\\ & & \mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{\mu }}_{{i}}\mathit{.}\end{array}$We use $\begin{array}{ccc}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{\Sigma }}}^{⊺}{{V}}_{{i}}^{-1}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}\mathit{,}\\ {{{\mu }}_{{\Sigma }}}^{⊺}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{V}}_{{i}}^{-1}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}\mathit{,}\\ {{{\mu }}_{{\Sigma }}}^{⊺}{{V}}_{{\Sigma }}^{-1}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}\mathit{,}\\ {{{\mu }}_{{\Sigma }}}^{⊺}& \mathrm{=}& \left(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}\right){{V}}_{{\Sigma }}\mathit{,}\end{array}$that we introduce in the previous equations to finally find $\begin{array}{ccc}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{Q}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{x}}^{⊺}{{V}}_{{i}}^{-1}{x}\mathrm{-}\mathrm{2}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{\Sigma }}}^{⊺}{{V}}_{{i}}^{-1}{x}\mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{\Sigma }}}^{⊺}{{V}}_{{i}}^{-1}{{\mu }}_{{\Sigma }}\\ & & \mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{\Sigma }}}^{⊺}{{V}}_{{i}}^{-1}{{\mu }}_{{\Sigma }}\mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{\mu }}_{{i}}\mathit{,}\\ & \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{(}{x}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}}^{⊺}{{V}}_{{i}}^{-1}\mathrm{(}{x}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}\\ & & \mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{{{\mu }}_{{i}}}^{⊺}{\mathit{V}}_{\mathit{i}}^{-1}{{\mu }}_{{i}}\mathrm{-}{{{\mu }}_{{\Sigma }}}^{⊺}{\mathit{V}}_{\mathrm{\Sigma }}^{-1}{{\mu }}_{{\Sigma }}\mathit{,}\end{array}$or, $\begin{array}{ccc}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{Q}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{(}{x}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}}^{⊺}{{V}}_{{i}}^{-1}\mathrm{(}{x}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}\\ & & \mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{\mu }}_{{i}}\mathrm{-}\mathrm{2}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{\mu }}_{{\Sigma }}\\ & & \mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{\Sigma }}}^{⊺}{{V}}_{{i}}^{-1}{{\mu }}_{{\Sigma }}\mathit{,}\\ & \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{(}{x}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}}^{⊺}{{V}}_{{i}}^{-1}\mathrm{(}{x}\mathrm{-}{{\mu }}_{{i}}\mathrm{)}\\ & & \mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{{\mu }}_{{\Sigma }}\mathrm{)}}^{⊺}{\mathit{V}}_{\mathit{i}}^{-1}\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{{\mu }}_{{\Sigma }}\mathrm{)}\mathit{.}\end{array}$

Appendix A.3: Expanding the Q_χ² term

Noting that $\mathit{I}\mathrm{=}{{V}}_{{\Sigma }}{{V}}_{{\Sigma }}^{-1}\mathrm{=}{{V}}_{{\Sigma }}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{V}}_{{j}}^{-1}\mathit{,}$(A.17)we can write $\begin{array}{ccc}{{V}}_{{i}}^{-1}\mathrm{-}{{V}}_{{i}}^{-1}{{V}}_{{\Sigma }}{{V}}_{{i}}^{-1}& \mathrm{=}& {{V}}_{{i}}^{-1}\mathrm{(}{I}\mathrm{-}{{V}}_{{\Sigma }}{{V}}_{{i}}^{-1}\mathrm{)}\mathit{,}\\ & \mathrm{=}& {{V}}_{{i}}^{-1}{{V}}_{{\Sigma }}\sum _{\mathit{j}\mathrm{=}\mathrm{1}\mathit{,j}\ne \mathit{i}}^{\mathit{n}}{{V}}_{{j}}^{-1}\mathit{,}\end{array}$and thus, with square symmetric matrices: $\begin{array}{ccc}{\mathit{Q}}_{{\mathit{\chi }}^{\mathrm{2}}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{\mu }}_{{i}}\mathrm{-}{{{\mu }}_{{\Sigma }}}^{⊺}{{V}}_{{\Sigma }}^{-1}{{\mu }}_{{\Sigma }}\mathit{,}\\ & \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{\mu }}_{{i}}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{V}}_{{\Sigma }}{{V}}_{{i}}^{-1}{{\mu }}_{{i}}\\ & & \mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\sum _{\mathit{j}\mathrm{=}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{n}}\mathrm{2}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{V}}_{{\Sigma }}{{V}}_{{j}}^{-1}{{\mu }}_{{j}}\mathit{,}\\ & \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{V}}_{{\Sigma }}\sum _{\mathit{j}\mathrm{=}\mathrm{1}\mathit{,j}\ne \mathit{i}}^{\mathit{n}}{{V}}_{{j}}^{-1}{{\mu }}_{{i}}\\ & & \mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\sum _{\mathit{j}\mathrm{=}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{n}}\mathrm{2}{{{\mu }}_{{i}}}^{⊺}{{V}}_{{i}}^{-1}{{V}}_{{\Sigma }}{{V}}_{{j}}^{-1}{{\mu }}_{{j}}\mathit{,}\\ & \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\sum _{\mathit{j}\mathrm{=}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{n}}{\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{{\mu }}_{{j}}\mathrm{)}}^{⊺}{{V}}_{{i}}^{-1}{{V}}_{{\Sigma }}{{V}}_{{j}}^{-1}\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{{\mu }}_{{j}}\mathrm{)}\mathit{.}\end{array}$

Appendix A.4: Sum of 2 quadratics: Q_χ term

We first develop: $\begin{array}{ccc}{{{\mu }}_{{{\Sigma }}_{{2}}}}^{⊺}{{V}}_{{{\Sigma }}_{{2}}}^{-1}{{\mu }}_{{{\Sigma }}_{{2}}}& \mathrm{=}& \mathrm{(}{{{\mu }}_{{1}}}^{⊺}{{V}}_{{1}}^{-1}\mathrm{+}{{{\mu }}_{{2}}}^{⊺}{{V}}_{{2}}^{-1}\mathrm{)}{\mathit{V}}_{{\mathrm{\Sigma }}_{\mathrm{2}}}\\ & & \mathrm{\times }\mathrm{(}{{V}}_{{1}}^{-1}{{\mu }}_{{1}}\mathrm{+}{{V}}_{{2}}^{-1}{{\mu }}_{{2}}\mathrm{)}\mathit{,}\\ & \mathrm{=}& \mathrm{(}{{{\mu }}_{{1}}}^{⊺}{{V}}_{{1}}^{-1}\mathrm{+}{{{\mu }}_{{2}}}^{⊺}{{{V}}_{{2}}^{-1}}^{\mathrm{)}}\frac{\mathrm{1}}{\mathrm{|}\frac{{{V}}_{{1}}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}\mathrm{+}\frac{{{V}}_{{2}}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}\mathrm{|}}\\ & & \mathrm{\times }\left(\frac{{{V}}_{{1}}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}\mathrm{+}\frac{{{V}}_{{2}}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}\right)\mathrm{(}{{V}}_{{1}}^{-1}{{\mu }}_{{1}}\mathrm{+}{{V}}_{{2}}^{-1}{{\mu }}_{{2}}\mathrm{)}\mathit{,}\\ & \mathrm{=}& \mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}[{{{\mu }}_{{1}}}^{⊺}\frac{{{V}}_{{1}}^{-1}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{\mu }}_{{1}}\mathrm{+}{{{\mu }}_{{2}}}^{⊺}\frac{{{V}}_{{2}}^{-1}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{\mu }}_{{2}}\\ & & \mathrm{+}{{{\mu }}_{{1}}}^{⊺}\frac{{{V}}_{{2}}^{-1}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{\mu }}_{{2}}\mathrm{+}{{{\mu }}_{{2}}}^{⊺}\frac{{{V}}_{{1}}^{-1}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{\mu }}_{{1}}]\\ & & \mathrm{+}\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}[{{{\mu }}_{{1}}}^{⊺}\frac{{{V}}_{{1}}^{-1}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{\mu }}_{{2}}\mathrm{+}{{{\mu }}_{{2}}}^{⊺}\frac{{{V}}_{{2}}^{-1}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{\mu }}_{{1}}\\ & & \mathrm{+}{{{\mu }}_{{1}}}^{⊺}{{V}}_{{1}}^{-1}\frac{{{V}}_{{2}}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{V}}_{{1}}^{-1}{{\mu }}_{{1}}\\ & & \mathrm{+}{{{\mu }}_{{2}}}^{⊺}{{V}}_{{2}}^{-1}\frac{{{V}}_{{1}}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{V}}_{{2}}^{-1}{{\mu }}_{{2}}]\mathit{.}\end{array}$Computing ${{V}}_{{{\Sigma }}_{{2}}}{{V}}_{{{\Sigma }}_{{2}}}^{-1}$: $\begin{array}{ccc}{{V}}_{{{\Sigma }}_{{2}}}{{V}}_{{{\Sigma }}_{{2}}}^{-1}& \mathrm{=}& \mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}\left(\frac{{{V}}_{{1}}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}\mathrm{+}\frac{{{V}}_{{2}}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}\right)\mathrm{(}{{V}}_{{1}}^{-1}\mathrm{+}{{V}}_{{2}}^{-1}\mathrm{)}\mathit{,}\\ {I}& \mathrm{=}& \frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{I}\mathrm{+}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{V}}_{{1}}{{V}}_{{2}}^{-1}\\ & & \mathrm{+}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{V}}_{{2}}{{V}}_{{1}}^{-1}\mathrm{+}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{I}\mathit{,}\\ \left(\mathrm{1}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}\right){I}& \mathrm{=}& \frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{V}}_{{1}}{{V}}_{{2}}^{-1}\mathrm{+}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{V}}_{{2}}{{V}}_{{1}}^{-1}\end{array}$We can write: $\begin{array}{ccc}{{V}}_{{2}}^{-1}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}\right)\mathit{I}& \mathrm{=}& \frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{V}}_{{2}}^{-1}{{V}}_{{1}}{{V}}_{{2}}^{-1}\mathrm{+}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{V}}_{{1}}^{-1}\mathit{,}\\ \frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{V}}_{{2}}^{-1}{{V}}_{{1}}{{V}}_{{2}}^{-1}& \mathrm{=}& \left(\mathrm{1}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}\right){{V}}_{{2}}^{-1}\\ & & \mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{V}}_{{1}}^{-1}\mathit{,}\end{array}$and, similarly: $\begin{array}{ccc}{{V}}_{{1}}^{-1}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}\right){I}& \mathrm{=}& \frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{V}}_{{2}}^{-1}\mathrm{+}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{V}}_{{1}}^{-1}{{V}}_{{2}}{{V}}_{{1}}^{-1}\mathit{,}\\ \frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{V}}_{{1}}^{-1}{{V}}_{{2}}{{V}}_{{1}}^{-1}& \mathrm{=}& \left(\mathrm{1}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}\mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}\right){{V}}_{{1}}^{-1}\\ & & \mathrm{-}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{V}}_{{2}}^{-1}\mathit{.}\end{array}$We use the above 3 relations to develop ${{{\mu }}_{{1}}}^{⊺}{{V}}_{{1}}^{-1}{{\mu }}_{{1}}\mathrm{+}{{{\mu }}_{{2}}}^{⊺}{{V}}_{{2}}^{-1}{{\mu }}_{{2}}\mathrm{-}{{{\mu }}_{{{\Sigma }}_{{2}}}}^{⊺}{{V}}_{{{\Sigma }}_{{2}}}^{-1}{{\mu }}_{{{\Sigma }}_{{2}}}\mathit{,}$(A.34)which leads to ${\mathrm{(}{{\mu }}_{{1}}\mathrm{-}{{\mu }}_{{2}}\mathrm{)}}^{⊺}\left(\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{2}}\mathrm{\right|}}{{V}}_{{1}}^{-1}\mathrm{+}\frac{\mathrm{\left|}{{V}}_{{{\Sigma }}_{{2}}}\mathrm{\right|}}{\mathrm{\left|}{{V}}_{{1}}\mathrm{\right|}}{{V}}_{{2}}^{-1}\right)\mathrm{(}{{\mu }}_{{1}}\mathrm{-}{{\mu }}_{{2}}\mathrm{)}\mathit{.}$(A.35)Using Eq. (52) we found that ${{{\mu }}_{{1}}}^{⊺}{{V}}_{{1}}^{-1}{{\mu }}_{{1}}\mathrm{+}{{{\mu }}_{{2}}}^{⊺}{{V}}_{{2}}^{-1}{{\mu }}_{{2}}\mathrm{-}{{{\mu }}_{{{\Sigma }}_{{2}}}}^{⊺}{{V}}_{{{\Sigma }}_{{2}}}^{-1}{{\mu }}_{{{\Sigma }}_{{2}}}$(A.36)is equal to ${\mathrm{(}{{\mu }}_{{1}}\mathrm{-}{{\mu }}_{{2}}\mathrm{)}}^{⊺}\mathrm{(}{{V}}_{{1}}\mathrm{+}{{V}}_{{2}}{\mathrm{)}}^{-1}\mathrm{(}{{\mu }}_{{1}}\mathrm{-}{{\mu }}_{{2}}\mathrm{)}\mathit{.}$(A.37)

Appendix A.5: χ and χ² distributions

Appendix A.5.1: Definition

The χ and χ² distributions with k = 2(n−1) degrees of freedom are defined as $\begin{array}{ccc}{\mathit{\chi }}_{\mathit{k}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \frac{{\mathrm{2}}^{\mathrm{1}\mathrm{-}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}}\mathit{,}\\ {\mathit{\chi }}_{\mathit{k}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \frac{{\mathrm{2}}^{\mathrm{-}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{x}}^{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{x}}{\mathrm{2}}}\mathit{,}\end{array}$with the gamma function ∀l ∈N, Γ(l) = (l−1) ! it leads to $\begin{array}{ccc}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \frac{{\mathrm{2}}^{\mathrm{2}\mathrm{-}\mathit{n}}}{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{!}}{\mathit{x}}^{\mathrm{2}\mathit{n}\mathrm{-}\mathrm{3}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}}\mathit{,}\\ {\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \frac{{\mathrm{2}}^{\mathrm{1}\mathrm{-}\mathit{n}}}{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{!}}{\mathit{x}}^{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{x}}{\mathrm{2}}}\mathit{.}\end{array}$So for $\begin{array}{cc}\mathit{n}& \mathrm{=}\mathrm{2}\\ \mathrm{=}\mathit{x}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}}& {\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\\ \mathit{n}& \mathrm{=}\mathrm{3}\\ \mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{3}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}}& {\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{4}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\\ \mathit{n}& \mathrm{=}\mathrm{4}\\ \mathrm{=}\frac{\mathrm{1}}{\mathrm{8}}{\mathit{x}}^{\mathrm{5}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}}& {\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{6}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\\ \mathit{n}& \mathrm{=}\mathrm{5}\\ \mathrm{=}\frac{\mathrm{1}}{\mathrm{48}}{\mathit{x}}^{\mathrm{7}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}}& {\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{8}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\end{array}$and so on.

Appendix A.5.2: Sum of two χ functions

We show here that χ_{k = 2p}(x₁) + χ_{k = 2q}(x₂) = χ_{k = 2(p + q)}(x). The distribution we are looking for is the density function of χ_{k = 2p}(x₁)χ_{k = q}(x₂)dx₁dx₂ given $\mathit{x}\mathrm{=}\sqrt{{\mathit{x}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{x}}_{\mathrm{2}}^{\mathrm{2}}}$. We use polar coordinates so, dx₁dx₂ = xdxdθ, x₁ = xcosθ and x₂ = xsinθ: $\begin{array}{ccc}\mathit{\chi }& \mathrm{=}& {\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{)}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x,}\\ & \mathrm{=}& \int \begin{array}{c}\frac{\mathit{\pi }}{\mathrm{2}}\\ \mathrm{0}\end{array}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathit{p}}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathrm{\right)}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathit{q}}\mathrm{\left(}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}\mathit{x}\mathrm{d}\mathit{x}\mathrm{d}\mathit{\theta ,}\\ & \mathrm{=}& \int \begin{array}{c}\frac{\mathit{\pi }}{\mathrm{2}}\\ \mathrm{0}\end{array}\frac{{\mathrm{2}}^{\mathrm{1}\mathrm{-}\mathit{p}}}{\mathrm{\Gamma }\mathrm{\left(}\mathit{p}\mathrm{\right)}}{\mathit{x}}^{\mathrm{2}\mathit{p}\mathrm{-}\mathrm{1}}{\cos }^{\mathrm{2}\mathit{p}\mathrm{-}\mathrm{1}}\mathit{\theta }\\ & & \mathrm{\times }\frac{{\mathrm{2}}^{\mathrm{1}\mathrm{-}\mathit{q}}}{\mathrm{\Gamma }\mathrm{\left(}\mathit{q}\mathrm{\right)}}{\mathit{x}}^{\mathrm{2}\mathit{q}\mathrm{-}\mathrm{1}}{\sin }^{\mathrm{2}\mathit{q}\mathrm{-}\mathrm{1}}\mathit{\theta }{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}}\mathit{x}\mathrm{d}\mathit{x}\mathrm{d}\mathit{\theta ,}\\ & \mathrm{=}& \frac{{\mathrm{2}}^{\mathrm{2}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{)}}}{\mathrm{\Gamma }\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{\Gamma }\mathrm{\left(}\mathit{q}\mathrm{\right)}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{)}\mathrm{-}\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}}\mathit{B}\mathrm{\left(}\mathit{p,q}\mathrm{\right)}\mathit{x}\mathrm{d}\mathit{x,}\\ & \mathrm{=}& \frac{{\mathrm{2}}^{\mathrm{1}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{)}}}{\mathrm{\Gamma }\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{)}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}}{\mathrm{2}}}\mathrm{d}\mathit{x,}\end{array}$in which B(p,q) is the beta function $\begin{array}{ccc}\mathit{B}\mathrm{\left(}\mathit{p,q}\mathrm{\right)}& \mathrm{=}& \mathrm{2}{\mathrm{\int }}_{\mathrm{0}}^{\frac{\mathit{\pi }}{\mathrm{2}}}{\cos }^{\mathrm{2}\mathit{p}\mathrm{-}\mathrm{1}}\mathit{\theta }{\sin }^{\mathrm{2}\mathit{q}\mathrm{-}\mathrm{1}}\mathit{\theta }\mathrm{d}\mathit{\theta ,}\\ & \mathrm{=}& \frac{\mathrm{\Gamma }\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{\Gamma }\mathrm{\left(}\mathit{q}\mathrm{\right)}}{\mathrm{\Gamma }\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{)}}\mathrm{·}\end{array}$

Appendix A.5.3: χ_{k = 2(n−1)} cumulative distribution function

Directly integrating for k = 2${\mathit{F}}_{{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{x}}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\left[\mathrm{-}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}\right]\begin{array}{c}\mathit{x}\\ \mathrm{0}\end{array}\mathrm{=}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}}$(A.53)For k = 4, we integrate by parts with $\begin{array}{c}\\ \mathit{u}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}& & {\mathit{u}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}{\mathit{x}}^{\mathrm{\prime }}\\ \mathit{v}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& {\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}& & {\mathit{v}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}{\mathit{x}}^{\mathrm{\prime }}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}\end{array}$(A.54)so $\begin{array}{ccc}{\mathit{F}}_{{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{4}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{x}\\ \mathrm{0}\end{array}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{4}}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}\\ & \mathrm{=}& {\mathit{F}}_{{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\left[\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}\right]\begin{array}{c}\mathit{x}\\ \mathrm{0}\end{array}\mathit{,}\\ & \mathrm{=}& {\mathit{F}}_{{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}}\mathit{.}\end{array}$For k = 6, also integrating by parts, we note $\begin{array}{c}\\ \mathit{u}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{8}}{\mathit{x}}^{\mathrm{\prime }\mathrm{4}}& & {\mathit{u}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{3}}\\ \mathit{v}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& {\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}& & {\mathit{v}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}{\mathit{x}}^{\mathrm{\prime }}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}\end{array}$(A.58)and so $\begin{array}{ccc}{\mathit{F}}_{{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{6}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{x}\\ \mathrm{0}\end{array}{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{6}}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}\\ & \mathrm{=}& {\mathit{F}}_{{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{4}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\left[\mathrm{-}\frac{\mathrm{1}}{\mathrm{8}}{\mathit{x}}^{\mathrm{\prime }\mathrm{4}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}\right]\begin{array}{c}\mathit{x}\\ \mathrm{0}\end{array}\mathit{,}\\ & \mathrm{=}& {\mathit{F}}_{{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{4}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{8}}{\mathit{x}}^{\mathrm{4}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}}\mathit{.}\end{array}$We deduce the general form for k = 2(n−1): ${\mathit{F}}_{{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}}\sum _{\mathit{i}\mathrm{=}\mathrm{2}}^{\mathit{n}}\frac{{\mathrm{2}}^{\mathrm{2}\mathrm{-}\mathit{i}}}{\mathrm{(}\mathit{i}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{!}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{i}\mathrm{-}\mathrm{2}\mathrm{)}}\mathit{.}$(A.62)

Appendix A.6: Computing the I_k,n(x) integral

Let us first expand a few notations for a better readability $\begin{array}{ccc}\begin{array}{c}\det \end{array}{\mathit{J}}_{\mathit{F}}& \mathrm{=}& {\mathit{x}}^{\mathit{n}\mathrm{-}\mathrm{2}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{3}}{\mathrm{􏽙}}}{\cos }^{\mathit{n}\mathrm{-}\mathit{i}\mathrm{-}\mathrm{2}}{\mathit{\theta }}_{\mathit{i}}\mathit{,}\\ & \mathrm{=}& {\mathit{x}}^{\mathit{n}\mathrm{-}\mathrm{2}}{\cos }^{\mathit{n}\mathrm{-}\mathrm{3}}{\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathit{n}\mathrm{-}\mathrm{4}}{\mathit{\theta }}_{\mathrm{2}}\mathit{...}\cos {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\mathit{.}\end{array}$We are supposed to use | detJ_F | but x is positive and all angles θ are ∈ [0,π/ 2], so detJ_F is always positive. Let us also expand $\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{x}}_{\mathit{i}}\mathrm{=}{\mathit{x}}^{\mathit{n}\mathrm{-}\mathrm{1}}{\cos }^{\mathit{n}\mathrm{-}\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathit{n}\mathrm{-}\mathrm{3}}{\mathit{\theta }}_{\mathrm{2}}\mathit{...}\cos {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\sin {\mathit{\theta }}_{\mathit{i}}\mathit{.}$(A.65)Multiplying both expressions leads to $\begin{array}{ccc}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{x}}_{\mathit{i}}\begin{array}{c}\det \end{array}{\mathit{J}}_{\mathit{F}}& \mathrm{=}& {\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\mathit{...}\\ & & \mathrm{\times }\cos {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\sin {\mathit{\theta }}_{\mathit{i}}\end{array}$(A.66)

Appendix A.6.1: Case I_k=1,n(x)

$\begin{array}{ccc}{\mathit{I}}_{\mathit{k}\mathrm{=}\mathrm{1}\mathit{,n}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\le }\mathit{x}\\ {\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}\end{array}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{\chi }}_{\mathrm{2}}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathrm{d}{\mathit{x}}_{\mathit{i}}\mathit{,}\\ & \mathrm{=}& \int \begin{array}{c}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\le }\mathit{x}\\ {\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}\end{array}\left(\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{x}}_{\mathit{i}}\right){\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}\\ & & \mathrm{\times }\left({\mathit{x}}^{\mathrm{\prime }\mathit{n}\mathrm{-}\mathrm{2}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{3}}{\mathrm{􏽙}}}{\cos }^{\mathit{n}\mathrm{-}\mathit{i}\mathrm{-}\mathrm{2}}{\mathit{\theta }}_{\mathit{i}}\right)\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\mathrm{d}{\mathit{\theta }}_{\mathit{i}}\mathit{,}\\ & \mathrm{=}& \int \begin{array}{c}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\le }\mathit{x}\\ {\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}\end{array}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{\mathrm{0}}^{\frac{\mathit{\pi }}{\mathrm{2}}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathit{i}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathit{i}}\sin {\mathit{\theta }}_{\mathit{i}}\mathrm{d}{\mathit{\theta }}_{\mathit{i}}\mathit{,}\\ & \mathrm{=}& \int \begin{array}{c}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\le }\mathit{x}\\ {\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}\end{array}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}}}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\frac{\mathit{B}\mathrm{\left(}\mathit{i,}\mathrm{1}\mathrm{\right)}}{\mathrm{2}}\mathit{,}\\ & \mathrm{=}& \int \begin{array}{c}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\le }\mathit{x}\\ {\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}\end{array}{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}\\ & \mathrm{=}& {\mathit{F}}_{{\mathit{\chi }}_{\mathit{k}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{.}\end{array}$The exact solution is given by Eq. (A.62).

We note the particular case I_k,n(k_γ) = γ if k_γ computed for this particular value of n.

Appendix A.6.2: Case I_k=n,n(x)

$\begin{array}{ccc}{\mathit{I}}_{\mathit{k}\mathrm{=}\mathit{n,n}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\le }\mathit{x}\\ {\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}\end{array}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{i}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathrm{d}{\mathit{x}}_{\mathit{i}}\mathit{,}\\ & \mathrm{=}& \mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{\int }}_{{\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{{\mathit{x}}^{\mathrm{\prime }}\mathrm{\le }\mathit{x}}\left(\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{x}}_{\mathit{i}}\right)\\ & & \mathrm{\times }\left({\mathit{x}}^{\mathrm{\prime }\mathit{n}\mathrm{-}\mathrm{2}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{3}}{\mathrm{􏽙}}}{\cos }^{\mathit{n}\mathrm{-}\mathit{i}\mathrm{-}\mathrm{2}}{\mathit{\theta }}_{\mathit{i}}\right)\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\mathrm{d}{\mathit{\theta }}_{\mathit{i}}\mathit{,}\\ & \mathrm{=}& \mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{\int }}_{{\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{{\mathit{x}}^{\mathrm{\prime }}\mathrm{\le }\mathit{x}}{\mathit{x}}^{\mathrm{\prime }\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\frac{\mathit{B}\mathrm{\left(}\mathit{i,}\mathrm{1}\mathrm{\right)}}{\mathrm{2}}\mathit{,}\\ & \mathrm{=}& \mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{1}}\frac{\mathrm{1}}{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}\frac{{\mathrm{2}}^{\mathrm{2}\mathrm{-}\mathit{n}}}{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{!}}\mathit{,}\\ & \mathrm{=}& \frac{{\mathit{\pi }}^{\mathit{n}\mathrm{-}\mathrm{1}}}{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{,}\end{array}$which is the volume of an hypersphere of dimension 2(n−1), also called 2(n−1)-sphere.

Appendix A.6.3: Intermediate case I_k,n(x)

For k> 1 and k<n: $\begin{array}{ccc}{\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}\mathit{X}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathit{k}}{\mathrm{􏽙}}}{\mathit{\chi }}_{\mathrm{2}}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\underset{\mathit{i}\mathrm{=}\mathit{n}\mathrm{-}\mathit{k}\mathrm{+}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{i}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathrm{d}{\mathit{x}}_{\mathit{i}}\mathit{,}\\ & \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{k}\mathrm{-}\mathrm{1}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}\\ & & \mathrm{\times }\sin {\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\sin {\mathit{\theta }}_{\mathrm{2}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\mathit{...}\\ & & \mathrm{\times }\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\cos {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}}\\ & & \mathrm{\times }\mathrm{d}\mathit{x}\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathit{,}\\ & \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\left(\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{k}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathit{B}\mathrm{\left(}\mathit{i,}\mathrm{1}\mathrm{\right)}\right)\\ & & \mathrm{\times }{\mathit{\pi }}^{\mathit{k}\mathrm{-}\mathrm{1}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}\sin {\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\\ & & \mathrm{\times }\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}{\cos }^{\mathrm{2}\mathit{k}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}}\\ & & \mathrm{\times }\mathrm{d}\mathit{x}\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathit{.}\end{array}$In a first step, we integrated by parts using $\begin{array}{ccc}\mathit{u}\mathrm{\left(}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathrm{\right)}& \mathrm{=}& {\cos }^{\mathrm{2}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathit{,}\\ {\mathit{u}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\mathrm{2}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathit{,}\\ {\mathit{c}}^{\mathrm{2}}& \mathrm{=}& {\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}\mathit{,}\\ \mathit{v}\mathrm{\left(}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathrm{\right)}& \mathrm{=}& {\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}}\mathit{,}\\ {\mathit{v}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathrm{\right)}& \mathrm{=}& {\mathit{x}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\cos {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}}\mathit{,}\\ \left[\mathit{uv}\right]\begin{array}{c}\mathit{\pi }\mathit{/}\mathrm{2}\\ \mathrm{0}\end{array}& \mathrm{=}& \mathrm{-}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}}\mathit{,}\\ \mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\frac{\mathit{\pi }}{\mathrm{2}}}{\mathit{u}}^{\mathrm{\prime }}\mathit{v}& \mathrm{=}& \int \mathrm{2}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\\ & & \mathrm{\times }{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathit{.}\end{array}$We thus find that $\begin{array}{ccc}{\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}\mathit{X}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{(}\mathit{k}\mathrm{+}\mathrm{2}\mathrm{)}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{(}\mathit{k}\mathrm{+}\mathrm{2}\mathrm{)}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}}{\mathit{\pi }}^{\mathit{k}\mathrm{-}\mathrm{1}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}}\\ & & \mathrm{\times }\sin {\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}{\cos }^{\mathrm{2}\mathit{k}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}\\ & & \mathrm{\times }\left[\mathrm{\int }{\mathit{u}}^{\mathrm{\prime }}\mathit{v}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}}\right]\mathrm{d}\mathit{x}\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}\mathit{.}\end{array}$(A.88)We finally find the recurrence formula ${\mathit{I}}_{\mathit{k,n}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathit{k,n}\mathrm{-}\mathrm{1}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{-}\mathrm{2}\mathit{\pi }{\mathit{I}}_{\mathit{k}\mathrm{-}\mathrm{1}\mathit{,n}\mathrm{-}\mathrm{1}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathit{,}$(A.89)since $\begin{array}{ccc}{\mathit{I}}_{\mathit{k,n}\mathrm{-}\mathrm{1}}\mathrm{\left(}\mathit{X}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}{\mathrm{􏽙}}}{\mathit{\chi }}_{\mathrm{2}}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\underset{\mathit{i}\mathrm{=}\mathit{n}\mathrm{-}\mathit{k}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{i}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\mathrm{d}{\mathit{x}}_{\mathit{i}}\\ & \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{k}\mathrm{-}\mathrm{1}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}}\\ & & \mathrm{\times }\sin {\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\sin {\mathit{\theta }}_{\mathrm{2}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{4}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\mathit{...}\\ & & \mathrm{\times }\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\cos {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}}\\ & & \mathrm{\times }\mathrm{d}\mathit{x}\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\\ & \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\left(\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{k}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\mathit{B}\mathrm{\left(}\mathit{i,}\mathrm{1}\mathrm{\right)}\right)\\ & & \mathrm{\times }{\mathit{\pi }}^{\mathit{k}\mathrm{-}\mathrm{1}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}}\\ & & \mathrm{\times }\sin {\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}\\ & & \mathrm{\times }{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}}\mathrm{d}\mathit{x}\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{k}}\end{array}$and $\begin{array}{ccc}{\mathit{I}}_{\mathit{k}\mathrm{-}\mathrm{1}\mathit{,n}\mathrm{-}\mathrm{1}}\mathrm{\left(}\mathit{X}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathit{k}}{\mathrm{􏽙}}}{\mathit{\chi }}_{\mathrm{2}}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\underset{\mathit{i}\mathrm{=}\mathit{n}\mathrm{-}\mathit{k}\mathrm{+}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\mathrm{2}\mathit{\pi }{\mathit{x}}_{\mathit{i}}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\mathrm{d}{\mathit{x}}_{\mathit{i}}\\ & \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{k}\mathrm{-}\mathrm{2}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}}\\ & & \mathrm{\times }\sin {\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\sin {\mathit{\theta }}_{\mathrm{2}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{4}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\\ & & \mathit{...}\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\cos {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\\ & & \mathrm{\times }{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}}\mathrm{d}\mathit{x}\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{3}}\\ & \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{\le }\mathit{X}\\ \mathit{x}\mathrm{=}\mathrm{0}\end{array}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\left(\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{k}\mathrm{-}\mathrm{2}}{\mathrm{􏽙}}}\mathit{B}\mathrm{\left(}\mathit{i,}\mathrm{1}\mathrm{\right)}\right)\\ & & \mathrm{\times }{\mathit{\pi }}^{\mathit{k}\mathrm{-}\mathrm{2}}{\mathit{x}}^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}}\\ & & \mathrm{\times }\sin {\mathit{\theta }}_{\mathrm{1}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{(}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\\ & & \mathrm{\times }{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}}\mathrm{d}\mathit{x}\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{1}\mathrm{-}\mathit{k}}\mathit{.}\end{array}$Knowing that ${\mathrm{\int }}_{\mathit{\theta }\mathrm{=}\mathrm{0}}^{\mathit{\theta }\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\sin \mathit{\theta }{\cos }^{\mathrm{2}\mathit{m}\mathrm{-}\mathrm{1}}\mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{B}\mathrm{\left(}\mathit{m,}\mathrm{1}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{1}}{\mathit{m}}\mathit{,}$(A.96)we integrate the differents parts: $\begin{array}{ccc}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathit{k}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathit{k}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}\mathit{...}{\mathrm{\int }}_{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{=}\mathrm{0}}^{{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}}& & \sin {\mathit{\theta }}_{\mathrm{2}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\sin {\mathit{\theta }}_{\mathrm{2}}{\cos }^{\mathrm{2}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\\ & & \mathit{...}\sin {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\cos {\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{\theta }}_{\mathit{n}\mathrm{-}\mathrm{2}}\\ & \mathrm{=}& \underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{k}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\frac{\mathit{B}\mathrm{\left(}\mathit{i,}\mathrm{1}\mathrm{\right)}}{\mathrm{2}}\mathit{,}\\ & \mathrm{=}& \frac{{\mathrm{2}}^{\mathrm{1}\mathrm{-}\mathit{k}}}{\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}}\mathrm{·}\end{array}$

Appendix B: Proper motion estimation and testing

Appendix B.1: Estimating proper motions

In this section, we show how it is possible to estimate the proper motion v of a source if the simplifying assumption of null proper motion made in Sect. 3 is not met. We neglect the parallax and the long term effect of the radial motion of the source, but those extra parameters could also be fitted provided we have enough catalogue measurements. In our simple case, the position p of a source at any time t can be computed from its position p₀ at a reference epoch (e.g. 2000): ${p}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{p}}_{\mathrm{0}}\mathrm{+}{v}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\left(\begin{array}{c}\\ \mathit{x}\\ \mathit{y}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ {\mathit{x}}_{\mathrm{0}}\mathrm{+}{\mathit{v}}_{\mathit{x}}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\\ {\mathit{y}}_{\mathrm{0}}\mathrm{+}{\mathit{v}}_{\mathit{y}}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\end{array}\right)\mathit{.}$(B.1)We assume we have n observations of the source at various epochs t_i. We want to estimate the proper motion, so to estimate the 4 unknowns (v_x,v_y,x₀,y₀). To do so we use the maximum-likelihood estimate which consists in maximizing the likelihood$\mathit{L}\mathrm{=}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{\mathrm{􏽙}}}{\mathrm{𝒩}}_{{{\mu }}_{{i}}\mathit{,}{{V}}_{{i}}}\mathrm{\left(}{p}\mathrm{\right(}{\mathit{t}}_{\mathit{i}}\mathrm{\left)}\mathrm{\right)}\mathit{,}$(B.2)therefore minimizing ${\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{p}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{p}\mathrm{(}{\mathit{t}}_{\mathit{i}}\mathrm{\left)}\mathrm{\right)}}^{⊺}{\mathit{V}}_{\mathit{i}}^{-1}\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{p}\mathrm{(}{\mathit{t}}_{\mathit{i}}\mathrm{\left)}\mathrm{\right)}\mathit{,}$(B.3)by solving the system of equations $\{\begin{array}{c}\\ \frac{\mathit{\partial }{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{p}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}}{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}& \mathrm{=}& \mathrm{0}& \mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{,}\\ \frac{\mathit{\partial }{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{p}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}}{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}& \mathrm{=}& \mathrm{0}& \mathrm{\left(}\mathit{b}\mathrm{\right)}\mathit{,}\\ \frac{\mathit{\partial }{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{p}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}}{\mathit{\partial }{\mathit{x}}_{\mathrm{0}}}& \mathrm{=}& \mathrm{0}& \mathrm{\left(}\mathit{c}\mathrm{\right)}\mathit{,}\\ \frac{\mathit{\partial }{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{p}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}}{\mathit{\partial }{\mathit{y}}_{\mathrm{0}}}& \mathrm{=}& \mathrm{0}& \mathrm{\left(}\mathit{d}\mathrm{\right)}\mathit{.}\end{array}$(B.4)To do so, we compute the derivative of χ²(p(t)) according to a parameter a_k$\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{p}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}}{\mathit{\partial }{\mathit{a}}_{\mathit{k}}}& \mathrm{=}& \mathrm{-}\mathrm{2}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}}[\frac{\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}\mathrm{-}{\mathit{p}}_{\mathit{x}}\mathrm{)}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathit{x}}}{\mathit{\partial }{\mathit{a}}_{\mathit{k}}}\mathrm{+}\frac{\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}\mathrm{-}{\mathit{p}}_{\mathit{y}}\mathrm{)}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathit{y}}}{\mathit{\partial }{\mathit{a}}_{\mathit{k}}}\\ & & \end{array}$(B.5)and the Jacobian matrix of p(t | v,p₀)$\begin{array}{ccc}{{J}}_{\mathrm{p}}\mathrm{\left(}{v}\mathit{,}{{p}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\left(\begin{array}{c}\\ \\ \end{array}\right)& \mathrm{=}& \left(\begin{array}{c}\\ & & & \\ & & & \end{array}\right)\mathit{,}\\ & \mathrm{=}& \left(\begin{array}{c}\\ & & & \\ & & & \end{array}\right)\mathit{.}\end{array}$So we have to solve the following system of equations, noting Δt_i = t_i−t₀$\{\begin{array}{c}\\ \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}}\left[\frac{\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}\mathrm{-}\mathit{x}\mathrm{)}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}\mathrm{-}\mathit{y}\mathrm{)}\right]\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}& \mathrm{=}& \mathrm{0}& \mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{,}\\ \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}}\left[\frac{\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}\mathrm{-}\mathit{y}\mathrm{)}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}\mathrm{-}\mathit{x}\mathrm{)}\right]\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}& \mathrm{=}& \mathrm{0}& \mathrm{\left(}\mathit{b}\mathrm{\right)}\mathit{,}\\ \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}}\left[\frac{\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}\mathrm{-}\mathit{x}\mathrm{)}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}\mathrm{-}\mathit{y}\mathrm{)}\right]& \mathrm{=}& \mathrm{0}& \mathrm{\left(}\mathit{c}\mathrm{\right)}\mathit{,}\\ \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}}\left[\frac{\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}\mathrm{-}\mathit{y}\mathrm{)}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\mathrm{(}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}\mathrm{-}\mathit{x}\mathrm{)}\right]& \mathrm{=}& \mathrm{0}& \mathrm{\left(}\mathit{d}\mathrm{\right)}\mathit{.}\end{array}$(B.8)We can for example use Cramer’s rule to solve the general problem ${AX}\mathrm{=}{\Lambda }$(B.9)with, in our case, and using the notation ${\mathit{\xi }}_{\mathit{i}}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}$$\begin{array}{ccc}{A}& \mathrm{=}& \left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)\mathit{,}\\ & \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{\xi }}_{\mathit{i}}\left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)\mathit{,}\end{array}$${X}\mathrm{=}\left(\begin{array}{c}\\ {\mathit{x}}_{\mathrm{1}}\\ {\mathit{x}}_{\mathrm{2}}\\ {\mathit{x}}_{\mathrm{3}}\\ {\mathit{x}}_{\mathrm{4}}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ {\mathit{v}}_{\mathit{x}}\\ {\mathit{v}}_{\mathit{y}}\\ {\mathit{x}}_{\mathrm{0}}\\ {\mathit{y}}_{\mathrm{0}}\end{array}\right)\mathit{,}$(B.12)and ${\Lambda }\mathrm{=}\left(\begin{array}{c}\\ {\mathit{e}}_{\mathrm{1}}\\ {\mathit{e}}_{\mathrm{2}}\\ {\mathit{e}}_{\mathrm{3}}\\ {\mathit{e}}_{\mathrm{4}}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\xi }}_{\mathit{i}}\left[\frac{{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\right]\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}\\ \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\xi }}_{\mathit{i}}\left[\frac{{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\right]\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}\\ \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\xi }}_{\mathit{i}}\left[\frac{{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\right]\\ \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\xi }}_{\mathit{i}}\left[\frac{{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\right]\end{array}\right)\mathit{,}$(B.13)leading to the solution ${\mathit{x}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{\left|}{{A}}_{{i}}\mathrm{\right|}}{\mathrm{\left|}{A}\mathrm{\right|}}\mathit{,}$(B.14)where $\begin{array}{ccc}{{A}}_{{1}}\mathrm{=}\left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)\mathit{,}{{A}}_{{2}}\mathrm{=}\left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)\mathit{,}& & \\ {{A}}_{{3}}\mathrm{=}\left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)\mathrm{and}{{A}}_{{4}}\mathrm{=}\left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)\mathit{.}& & \end{array}$

Appendix B.2: Estimating the error on the proper motion estimate

The covariance matthat isV_f on the estimated proper motion parameters is provided by the inverse of the Hessian matrix Hf of $\ln \left(\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{n}}{{}^{\mathrm{􏽑}}}}{\mathrm{𝒩}}_{\mathit{b}}\mathrm{\left(}{p}\mathrm{\right(}{\mathit{t}}_{\mathit{i}}\mathrm{\left)}\mathrm{\right)}\right)$, evaluated with the estimated parameters, that is by the matrix $\mathrm{-}{{H}}_{{f}}^{-1}\mathrm{=}{{\left(\begin{array}{c}\\ \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{v}}_{\mathit{x}}\mathit{\partial }{\mathit{v}}_{\mathit{y}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{v}}_{\mathit{x}}\mathit{\partial }{\mathit{x}}_{\mathrm{0}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{v}}_{\mathit{x}}\mathit{\partial }{\mathit{y}}_{\mathrm{0}}}\\ \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{v}}_{\mathit{y}}\mathit{\partial }{\mathit{v}}_{\mathit{x}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{v}}_{\mathit{y}}\mathit{\partial }{\mathit{x}}_{\mathrm{0}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{v}}_{\mathit{y}}\mathit{\partial }{\mathit{y}}_{\mathrm{0}}}\\ \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{x}}_{\mathrm{0}}\mathit{\partial }{\mathit{v}}_{\mathit{x}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{x}}_{\mathrm{0}}\mathit{\partial }{\mathit{v}}_{\mathit{y}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{x}}_{\mathrm{0}}^{\mathrm{2}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{x}}_{\mathrm{0}}\mathit{\partial }{\mathit{y}}_{\mathrm{0}}}\\ \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{y}}_{\mathrm{0}}\mathit{\partial }{\mathit{v}}_{\mathit{x}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{y}}_{\mathrm{0}}\mathit{\partial }{\mathit{v}}_{\mathit{y}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{y}}_{\mathrm{0}}\mathit{\partial }{\mathit{x}}_{\mathrm{0}}}& \frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{y}}_{\mathrm{0}}^{\mathrm{2}}}\end{array}\right)}^{-1}|}_{\mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{x}}\end{array}\mathit{,}\begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{y}}\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{x}}_{\mathrm{0}}\end{array}\mathit{,}\begin{array}{c}\mathrm{˜}\\ {\mathit{y}}_{\mathrm{0}}\end{array}\mathrm{\right)}}$(B.15)${-}{{H}}_{{f}}^{-1}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{\xi }}_{\mathit{i}}\left(\begin{array}{c}\\ \frac{{\mathrm{\Delta }}^{\mathrm{2}}{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}& \mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}{\mathrm{\Delta }}^{\mathrm{2}}{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}& \frac{\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}& \mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\\ \mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}{\mathrm{\Delta }}^{\mathrm{2}}{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}& \frac{{\mathrm{\Delta }}^{\mathrm{2}}{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}& \mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}& \frac{\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\\ \frac{\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}& \mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}& \frac{\mathrm{1}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}& \mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\\ \mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}& \frac{\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}& \mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}& \frac{\mathrm{1}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\end{array}\right)$(B.16)${{V}}_{{f}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{\left|}{{H}}_{{f}}\mathrm{\right|}}{\mathrm{com}\mathrm{\left(}{{H}}_{{f}}\mathrm{\right)}}^{⊺}\mathit{.}$(B.17)

Appendix B.3: Simple case: no covariance

If all positional errors are circles (i.e. ∀i ∈ [1,n] ,ρ_i = 0), the simplifications leads to the classical formulae in which x and y are computed independently (see Press et al. 2007, p. 781, Sect. 15.2 “Fitting Data to a Straight Line”) $\begin{array}{c}\\ {\mathrm{\Delta }}_{\mathit{x}}& \mathrm{=}& {\mathit{S}}_{\mathit{x}}{\mathit{S}}_{{\mathit{t}}_{\mathit{x}}{\mathit{t}}_{\mathit{x}}}\mathrm{-}\mathrm{(}{\mathit{S}}_{{\mathit{t}}_{\mathit{x}}}{\mathrm{)}}^{\mathrm{2}}\mathit{,}& & {\mathrm{\Delta }}_{\mathit{y}}& \mathrm{=}& {\mathit{S}}_{\mathit{y}}{\mathit{S}}_{{\mathit{t}}_{\mathit{y}}{\mathit{t}}_{\mathit{y}}}\mathrm{-}\mathrm{(}{\mathit{S}}_{{\mathit{t}}_{\mathit{y}}}{\mathrm{)}}^{\mathrm{2}}\mathit{,}\\ {\mathit{v}}_{\mathit{x}}& \mathrm{=}& \frac{{\mathit{S}}_{{\mathit{t}}_{\mathit{x}}{\mathit{t}}_{\mathit{x}}}{\mathit{S}}_{{\mathit{\mu }}_{\mathit{x}}}}{{\mathrm{\Delta }}_{\mathit{x}}}\mathit{,}& & {\mathit{v}}_{\mathit{y}}& \mathrm{=}& \frac{{\mathit{S}}_{{\mathit{t}}_{\mathit{y}}{\mathit{t}}_{\mathit{y}}}{\mathit{S}}_{{\mathit{\mu }}_{\mathit{y}}}}{{\mathrm{\Delta }}_{\mathit{y}}}\mathit{,}\\ {\mathit{x}}_{\mathrm{0}}& \mathrm{=}& \frac{{\mathit{S}}_{\mathit{x}}{\mathit{S}}_{{\mathit{t}}_{\mathit{x}}{\mathit{\mu }}_{\mathit{x}}}\mathrm{-}{\mathit{S}}_{{\mathit{t}}_{\mathit{x}}}{\mathit{S}}_{{\mathit{\mu }}_{\mathit{x}}}}{{\mathrm{\Delta }}_{\mathit{x}}}\mathit{,}& & {\mathit{y}}_{\mathrm{0}}& \mathrm{=}& \frac{{\mathit{S}}_{\mathit{y}}{\mathit{S}}_{{\mathit{t}}_{\mathit{y}}{\mathit{\mu }}_{\mathit{y}}}\mathrm{-}{\mathit{S}}_{{\mathit{t}}_{\mathit{y}}}{\mathit{S}}_{{\mathit{\mu }}_{\mathit{y}}}}{{\mathrm{\Delta }}_{\mathit{y}}}\mathit{,}\end{array}$(B.18)where $\begin{array}{c}\\ {\mathit{S}}_{\mathit{x}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{1}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\mathit{,}& & {\mathit{S}}_{{\mathit{t}}_{\mathit{x}}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\mathit{,}& & {\mathit{S}}_{{\mathit{\mu }}_{\mathit{x}}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\mathit{,}\\ {\mathit{S}}_{{\mathit{t}}_{\mathit{x}}{\mathit{t}}_{\mathit{x}}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{\mathrm{\Delta }}^{\mathrm{2}}{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\mathit{,}& & {\mathit{S}}_{{\mathit{t}}_{\mathit{x}}{\mathit{\mu }}_{\mathit{x}}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{\mathit{\mu }}_{{\mathit{i}}_{\mathit{x}}}\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\mathit{,}& & & & \\ {\mathit{S}}_{\mathit{y}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{1}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\mathit{,}& & {\mathit{S}}_{{\mathit{t}}_{\mathit{y}}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\mathit{,}& & {\mathit{S}}_{{\mathit{\mu }}_{\mathit{y}}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\mathit{,}\\ {\mathit{S}}_{{\mathit{t}}_{\mathit{y}}{\mathit{t}}_{\mathit{y}}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{\mathrm{\Delta }}^{\mathrm{2}}{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\mathit{,}& & {\mathit{S}}_{{\mathit{t}}_{\mathit{y}}{\mathit{\mu }}_{\mathit{y}}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{\mathit{\mu }}_{{\mathit{i}}_{\mathit{y}}}\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\mathit{,}& & & & \end{array}$(B.19)and associated errors are $\begin{array}{c}\\ {\mathit{\sigma }}_{{\mathit{v}}_{\mathit{x}}}^{\mathrm{2}}& \mathrm{=}& \frac{{\mathit{S}}_{{\mathit{t}}_{\mathit{x}}{\mathit{t}}_{\mathit{x}}}}{{\mathrm{\Delta }}_{\mathit{x}}}\mathit{,}& & {\mathit{\sigma }}_{{\mathit{v}}_{\mathit{y}}}^{\mathrm{2}}& \mathrm{=}& \frac{{\mathit{S}}_{{\mathit{t}}_{\mathit{y}}{\mathit{t}}_{\mathit{y}}}}{{\mathrm{\Delta }}_{\mathit{y}}}\mathit{,}\\ {\mathit{\sigma }}_{{\mathit{x}}_{\mathrm{0}}}& \mathrm{=}& \frac{{\mathit{S}}_{\mathit{x}}}{{\mathrm{\Delta }}_{\mathit{x}}}\mathit{,}& & {\mathit{\sigma }}_{{\mathit{y}}_{\mathrm{0}}}& \mathrm{=}& \frac{{\mathit{S}}_{\mathit{y}}}{{\mathrm{\Delta }}_{\mathit{y}}}\mathit{,}\\ \mathit{\rho }{\mathit{\sigma }}_{{\mathit{v}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{x}}_{\mathrm{0}}}& \mathrm{=}& \frac{\mathrm{-}{\mathit{S}}_{{\mathit{t}}_{\mathit{x}}}}{{\mathrm{\Delta }}_{\mathit{x}}}\mathit{,}& & \mathit{\rho }{\mathit{\sigma }}_{{\mathit{v}}_{\mathit{y}}}{\mathit{\sigma }}_{{\mathit{y}}_{\mathrm{0}}}& \mathrm{=}& \frac{\mathrm{-}{\mathit{S}}_{{\mathit{t}}_{\mathit{y}}}}{{\mathrm{\Delta }}_{\mathit{y}}}\mathrm{·}\end{array}$(B.20)

Appendix B.4: Verifying the results

We have implemented and tested the result given by equation Eq. (B.14). We compare the results with a modified version of the Levenberg-Marquardt (LM) method (see Press et al. 2007, p. 801, Sect. 15.5.2 “Levenberg-Marquardt method”) we designed tohandle binormal distributions. The algorithm is the same except that we replace the term $\frac{\mathit{\partial }{\mathit{\chi }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{a}}_{\mathit{k}}}$ in β_k by Eq. (B.5) and α_kl by $\begin{array}{ccc}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}}& & [\frac{\mathrm{1}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}^{\mathrm{2}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathit{x}}}{\mathit{\partial }{\mathit{a}}_{\mathit{k}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathit{x}}}{\mathit{\partial }{\mathit{a}}_{\mathit{l}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}^{\mathrm{2}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathit{y}}}{\mathit{\partial }{\mathit{a}}_{\mathit{k}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathit{y}}}{\mathit{\partial }{\mathit{a}}_{\mathit{l}}}\\ & & \mathrm{-}\frac{{\mathit{\rho }}_{\mathit{i}}}{{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{x}}}{\mathit{\sigma }}_{{\mathit{i}}_{\mathit{y}}}}\left(\frac{\mathit{\partial }{\mathit{p}}_{\mathit{x}}}{\mathit{\partial }{\mathit{a}}_{\mathit{k}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathit{y}}}{\mathit{\partial }{\mathit{a}}_{\mathit{l}}}\mathrm{-}\frac{\mathit{\partial }{\mathit{p}}_{\mathit{y}}}{\mathit{\partial }{\mathit{a}}_{\mathit{k}}}\frac{\mathit{\partial }{\mathit{p}}_{\mathit{x}}}{\mathit{\partial }{\mathit{a}}_{\mathit{l}}}\right)]\mathrm{·}\end{array}$(B.21)We initialize the LM parameters with the approximate solutions given in Eq. (B.18). The results obtained using both methods (LM and Eq. (B.14)) are identical.

Appendix B.5: Testing the unique source hypothesis

When estimating the proper motion, we formulated the hypothesis H than our n observations come from a single underlying source. The Chi-square of equation Eq. (B.3) follows a Chi-square distribution with 2n−4 = 2(n−2) degrees of freedom. Therefore the criteria not to reject H is ${\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{p}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{p}\mathrm{(}{\mathit{t}}_{\mathit{i}}\mathrm{\left)}\mathrm{\right)}}^{⊺}{{V}}_{{i}}^{-1}\mathrm{(}{{\mu }}_{{i}}\mathrm{-}{p}\mathrm{(}{\mathit{t}}_{\mathit{i}}\mathrm{\left)}\mathrm{\right)}\mathrm{\le }{\mathit{k}}_{\mathit{\gamma }\mathrm{\approx }\mathrm{0.9973}}^{\mathrm{2}}\mathit{,}$(B.22)in which ${\mathit{k}}_{\mathit{\gamma }\mathrm{\approx }\mathrm{0.9973}}^{\mathrm{2}}\mathrm{=}{\mathit{F}}_{{\mathit{\chi }}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}}^{\mathrm{2}}}^{-1}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}$.

Appendix C: Synthetic catalogues generation script

Here is the script used to generate three synthetical tables and cross-match them with the online ARCHES XMatch Tool. The language of the script is specific to the tool. Both the tool and its documentation are available online¹⁰.

synthetic seed=1 nTab=3 prefix=true \
  geometry=cone ra=22.5 dec=33.5 r=0.42 \
  nA=40000 nB=20000 nC=35000 \
  nAB=6000 nAC=12000 nBC=18000 \
  nABC=10000 \
  poserrAtype=CIRCLE poserrAmode=formula paramA1=0.4 \
  poserrBtype=CIRCLE poserrBmode=function paramB1func=x \
                     paramB1xmin=0.8 paramB1xmax=1.2 \
                     paramB1nstep=100 \
  poserrCtype=CIRCLE poserrCmode=function \
                     paramC1func=exp(-0.5*(x-0.75)*(x-0.75)/0.01)/(0.1*sqrt(2*PI)) \
                     paramC1xmin=0.5 paramC1xmax=1 \
                     paramC1nstep=100
save prefix=simu3 suffix=.fits common=simu3.fits format=fits

cleartables

get FileLoader file=simu3A.fits
set pos ra=posRA dec=posDec
set poserr type=CIRCLE param1=ePosA param2=ePosB param3=ePosPA
set cols *

get FileLoader file=simu3B.fits
set pos ra=posRA dec=posDec
set poserr type=CIRCLE param1=ePosA param2=ePosB param3=ePosPA
set cols *

xmatch chi2 completeness=0.9973 nStep=1 nMax=2 join=inner
merge pos chi2
merge dist mec

get FileLoader file=simu3C.fits
set pos ra=posRA dec=posDec
set poserr type=CIRCLE param1=ePosA param2=ePosB param3=ePosPA
set cols *

xmatch chi2 completeness=0.9973 nStep=2 nMax=2 join=inner
merge pos chi2
merge dist mec

save simu3.ABC.fits fits   

All Tables

All Figures