© ESO 2018

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0;), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

With the aim of using an improved method to study the dynamics of galaxy pairs, we constructed two new pair catalogs using well-defined criteria and improved observational data. The first is a catalog of ≈1000 pairs with high accuracy radial velocities (Chamaraux & Nottale 2016) extracted from Nilson’s Uppsala Galaxy Catalog (UGC; Nilson 1973), which has the advantage that it is complete in apparent diameter. The second is a catalog of ≈13 000 pairs (Nottale & Chamaraux 2018, the largest presently available to our knowledge), which makes use of the huge recent increase in astronomical data. We used the HyperLEDA database (2016, Makarov et al. 2014) to identify galaxy pairs and extract their parameters from the large surveys (e.g., SDSS; Alam et al. 2015 and 2MASS).

To improve our understanding of physics at extragalactic scales, the goal of this paper is to elaborate better methods of statistical analysis of galaxy pair data (which can also be applied to any type of pairs of astronomical objects). This is the goal of the present paper. We need to obtain statistical information on the physical characteristics of these pairs, particularly their masses through Kepler’s third law and the possible existence of anomalous dynamics in these pairs. Such a new dynamics is, according to various proposals, attributed to missing mass or dark matter (Bergstrom 2000), to modification of gravity (Sanders 2002), or to a new dark potential (Nottale 2011; Chavanis 2017a,b).

In order to compute those physical quantities, we have to know the 3D velocity difference between the pair members and their 3D interdistances. But we have at our disposal only one component of the velocity difference (along the line of sight) and two components of the interdistance (projection on the sky-plane). Therefore one has to find statistical methods to obtain those 3D quantities from the projected quantities. We call these methods statistical deprojection. Moreover, we note that for each galaxy pair we can determine the various 3D parameters for only one point of each galaxy orbit and only one instant.

Because of these limitations, various methods of analysis have been devised (Chengalur et al. 1996; Peterson 1979; Faber & Gallagher 1979). However, as shown by Faber & Gallagher (1979), these methods remain unsatisfactory from a mathematical viewpoint. In particular, recovering the mass remains very uncertain. For these reasons, we propose new methods of deprojection, which we show here to yield more precise results for recovering the 3D unprojected parameters.

The paper is organized as follows. In Sect. 2 we first analyze the statistical projection process and derive the probability distributions expected for the various relevant projected quantities, namely, intervelocity, v_z; interdistance, r_p; their ratio, r_p ∕v_z, which can be used as signature of circular orbits; and the product, , which is involved in mass determination.

In Sect. 3, we develop statistical methods of deprojection of these parameters, that is, intervelocity, interdistance, and mass, based on the previous analysis. We start from a histogram of the projected data organized in various ways, i.e., constant bins, moving bins, and variable moving bins. Then we apply inversion formulae, which amount in some cases to matrix inversion, to obtain the deprojected distributions. Finally we test these methods using numerical simulations, which also allow us to determine the uncertainties involved.

2 Statistical analysis

Let us use a cylindrical coordinate system whose axis z is oriented from the observer to the galaxy pair. In this case the (x, y) plane is the plane of the sky. The position vector r is defined as (1)

while the velocity vector v is defined as (2)

The measured quantities are the position vector projected on the plane of the sky, r_p = (x, y), and the radial velocity v_z. Then r_p, v_z, and φ are known while z, v_p, and φ_v are unknown variables.

2.1 Probability distribution of radial velocity

Consider a randomly oriented 3D vector (x, y, z) of fixed length r. Let us determine the probability distribution of any of its projections, say z. The part of the sphere of radius r which projects between z and z + dz has a surface (3)

Since z = rcosθ, then d z = −rsinθ dθ, and we finally find that dS∕dz = 2πr. The probability distribution of the projection on any axis of a randomly oriented 3D vector of length r can then be derived (4)

and it is therefore constant when r is constant. Since z can vary from − r to + r, we verify that this probability distribution is correctly normalized.

This can be easily applied to the velocity vector. Given a fixed pair configuration (r, v) with random orientation, the probability distribution of the values of the radial velocity differences between the pair members is (5)

We note that here no member of the pair is priviledged, such that we consider only positive differences, which then vary between 0 and v.

We performed a numerical simulation of this projection of a given pair with random orientation. The obtained distribution of v_z confirms this expectation (Fig. 1).

2.2 Probability distribution of projected distance

We consider a pair of objects, the 3D interdistance between which is r, fixed. Let be the distance between the pair members projected on the plane of the sky. The part of the sphere corresponding to projected distances lying in the interval [r_p, r_p + dr_p] is made of two rings of width rdθ and radius r_p. Therefore its surface isdS = 4πr_prdθ. Now, since r_p = rsinθ, we find (6)

The differential probability distribution of r_p is d P = dS∕(4πr²). Then the normalized probability density of r_p values p(r_p), projected from a given r value, is written as (7)

We performed a numerical simulation of such a projection. The obtained distribution is in fair agreement with this theoretical expectation (Fig. 2).

More generally, we consider now a set of pairs whose distances between their members are distributed with a probability P_r (r) when r lies in the interval [r₁, r₂]. When the pair orientation is random, the expected probability distribution of the projected distances on the plane of the sky is now given by (8)

2.3 Probability distribution of ratio ζ = r_p∕v_z

This ratio can be used as signature of circular orbits. We give in Appendix A the probability distribution of a product and of the ratio of two variables, of which the individual probability distributions are known, and also that of the inverse of a variable. From these basic formulae, one can infer the probability distribution functions (PDFs) of various combinations of the variables observed for pairs [v_z and (x, y), involving r_p].

Since with 0 < r_p < 1 (for r = 1) and p(v_z)= 1 with 0 < v_z < 1 (for v = 1), we obtain for the ratio ζ = r_p∕v_z [more generally, ζ = (r_p∕r)∕(v_z∕v)] (9)

which is written as (10)

Two cases must now be considered. When ζ > 1, the integration interval is reduced to [0, 1∕ζ]. One finds (11)

When ζ < 1, the integration interval is again [0, 1] and one finds (12)

We performed a numerical simulation of randomly oriented pairs with uncorrelated randomly oriented velocity differences. As can be seen in Fig. 3, the theoretically expected distribution agrees very well with the simulation.

In the case of circular orbits, once r_p is given, the possible values of v_z, instead of being uniformly distributed between 0 and v, are constrained to be smaller than v × (r_p∕r). As a consequence, the left part of the PDF of ζ = (r_p∕r)∕(v_z∕v) below ζ = 1 is expected to be empty, which achieves a possible statistical signature of circular orbits.

The probability distribution of the reverse function, χ = v_z∕r_p is easy to derive from these expressions and from Eq. (A.3). We find

It is therefore constant up to χ = 1, after which it decreases quickly toward 0. In the case of circular orbits, only the constant part of the PDF remains (0 < χ < 1), which means that the projection of v_z∕r_p becomes similar to that of v_z.

2.4 Probability distribution of product

This product is involved in mass determination. Indeed, the orbits of isolated galaxy pairs are subjected to Kepler’s third law, (15)

where M = M₁ + M₂, a is the semimajor axis of the orbit and T the period. Defining a characteristic velocity V = 2πa∕T, it yields the total mass of the system (16)

The perimeter of an ellipse varies between 4a (radial free fall) and 2πa (circular orbit). This perimeter is given by elliptic integrals, which can be approximated by a power series (17)

where e is the eccentricity of the ellipse, and therefore the average velocity is such that .

When the orbit is circular, r = a and v = V, such that in this case the total mass is given by (18)

Instead of r we have only access to r_p, and, instead of v, to v_z. This leads us to look for the probability distribution of .

Let us set , then , while we recall that when the radius is normalized to r = 1.

In the case of uncorrelated values of r and v, the general formula for the probability distribution of a product Eq. (A.1) then yields (19)

The integration limits are determined by the fact that 0 < v_z < 1 (for a normalized velocity v = 1) and hence U < 1; U² − η² must be > 0 due to the square root such that U > η. This integral can be integrated in terms of the hypergeometric function ₂ F₁ (a, b;c;z) as follows: (20)

This formula is in excellent agreement with the result of a numerical simulation of the projections of r to r_p and v to v_z (Fig. 4).

3 Methods of statistical deprojection

The knowledge of the expected statistical distribution of the various variables or of their combination allows one to construct methods of deprojection from the observed subset of variables.

3.1 Deprojection of x, y, and v_z

3.1.1 Theoretical deprojection of v_z

The simplest method of deprojection deals with the value of a vector projected on a single axis. This is the case in particular for v_z (radial velocity difference). The various methods described in this work are also valid for deprojection of the individual observed variables x and y (projections of r on the plane of the sky along right ascension and declination).

We have seen that, in this case, the expected probability distribution for a random orientation and a given value of the unprojected variable is constant. Therefore, if P_v (v) is the probability distribution of the 3D velocity v, the probability distribution of v_z is given by Nottale (2011) (21)

This distribution can be easily inverted. One obtains a deprojection formula, that is, (22)

and similar formulae for x and y separately; the deprojection of their combination in is considered in the following.

3.1.2 General algorithm of deprojection of v_z

The simplest way to achieve statistical deprojection of a single variable, such as the radial velocity v_z, consists of directly implementing Eq. (22). This formula means that if N_v pairs have a true velocity difference v, their radial velocities are uniformly distributed between 0 and v. The existence of these N_v objects at velocity v creates a jump δN = N_v in the observed distribution of v_z, N(v_z) = N_totP(v_z). Therefore, the contribution of a particular true intervelocity v to the observed PDF of v_z is given by the surface of the rectangle of sides v × δN (see Fig. 5).

3.1.3 Various methods of deprojection of v_z

Differences on adjacent constant bins. The simplest way to implement Eq. (22) consists of

• (1)

  constructing the histogram of radial (projected) velocities V _r in bins [V _i−1, V _i] of given width δV;

• (2)

  computing the differences between the numbers in successive bins;

• (3)

  multiplying by the rank i = V _i∕δV of the bin.

We obtainedin this way the deprojected numbers, N = V (δN^p∕δV), i.e., in terms of the elements of the histogram, . This numbershould be assigned to a bin of width δV and of mean deprojected velocity . Indeed, the projected velocity V _i = i δV is the limit between the bins of ranks i and i + 1 in the projected histogram, and the resulting N_i should be assigned to this velocity V _i in the mean, i.e., to a bin centered on this velocity and therefore shifted by δV∕2 with respect to the initial histogram.

We performed a numerical simulation in which we calculated the mean and standard deviations of 100 realizations of the deprojection (see Fig. 6). We recovered in these simulations the overall shape of the initial distribution (chosen to be a Gaussian of mean 190 km s⁻¹ and standard deviation 90 km s⁻¹). In particular the existence and position of the main peak is recovered in a satisfactory way. We can see in Fig. 6 the improvement obtained when going from a catalog of ≈ 1000 pairs to ≈ 10 000, regarding the decrease of both the optimal bin width and dispersion of the various PDFs. A given realization of a catalog deprojection is expected to be contained (with a high probability) between the two ± 1σ lines, as supported by an example given in Fig. 7 (left part).

However, Eq. (22) presupposes a strictly monotonic decreasing distribution. Otherwise it may lead to obtainingnegative numbers. The problem is that there are fluctuations that may locally break this expected monotony.

For large enough bin sizes, it is clear that the monotony is preserved. Therefore, a possible way to solve this difficulty without losing too much resolution in the deprojection is to choose the smallest of these large enough bins. This value depends on the total number of pairs. In our numerical simulations, the optimal bin width was found to be ≈ 30 km s⁻¹ for ≈1000 pairs and ≈ 20 km s⁻¹ for ≈10 000 pairs.

Another possible method consists of correcting the negative numbers by substracting these numbers to the adjacent bins. This method is possible only for small deviations and reveals to be equivalent to the previous method.

Differences between constant intervals separated by two bins. Actually, the previous method where the difference is taken between two adjacent bins is not optimized and it can therefore be improved. It is more efficient (as in finite difference methods) to take differences between two intervals separated by one bin, . This improvement is based on thefact that f(x + dx) − f(x) = f′(x)dx + O(dx²), while (f(x + dx) − f(x −dx))∕2 = f′(x)dx + O(dx³).

We checked the method by a numerical simulation (see right Fig. 7). We find that the ± 1σ standard deviation on the PDF is half that obtained with the previous method (differences ), thus achieving a significant improvement.

Differences on constant moving bins. The method that uses constant bins has drawbacks: the result isdigitalized too much, has low resolution, and is too dependent on limits between bins.

Then we devised another method that deals with more information about the original radial velocity PDF and is more continuous. It amounts to performing a histogram of the projected velocities on a large enough moving bin w_bin shifted by a low value, e.g., δv = 1 km s⁻¹. We obtain a function that is both monotonous and quasi-continuous. An example of such a function is given in Fig. 8: it is obtained from the projection of ≈ 12 000 velocities having a Gaussian distribution of mean 150 km s⁻¹ and standard deviation 50 km s⁻¹.

Then we applied on this function the deprojection formula (Eq. (22)), but with thesmall shift dv_r = δv. Then we performed a moving average of the resulting distribution using the original bin width w_bin.

Numerical simulations were achieved to validate the method. We first defined a given PDF of the 3D velocity v, performed a random realization of this PDF on N_v values v_i, randomly projected them to v_ri, and then we deprojected the distribution obtained using the moving bins method. The process was repeated N_s times to obtain the convenient statistics of the method.

We give in Fig. 9 the results obtained for an initial v distribution showing one unique peak described by a Gaussian for two different peak widths. The quality of the result depends on the width ofthe initial peak: for a large enough peak (left Fig. 9 and Fig. 10), the position, amplitude, and width of the peak are correctly recovered. When the peak is narrower (right Fig. 9), its position is correctly recovered by the deprojection, but its amplitude is slightly too low and correspondingly its width becomes too large (but only within one sigma).

Differences on varying moving bins. In order to account for this bias, we devised another more accurate method. The problem is that, when there is a peak in the PDF of the 3D true velocity v, it manifestsas a high value of the slope in the PDF of the radial (projected) velocity v_r. If the binwidth is too large, this slope is decreased by the smoothing out effect of the binning. As a consequence, the amplitude of the deprojected peak is too low. Thus we corrected this bias by using a moving bin of variable width, decreasing when the slope increases. A balance should be found for the variation of the binwidth, since a smaller bin also increases the fluctuations and therefore the final dispersion. This method gives very good results, since it allows us to recover multiple peaks present in the initial 3D distribution, as can be seen in Fig. 11.

3.2 Deprojection of r_p

3.2.1 Theoretical deprojection

Although one can obtain the probability distribution of r from x and y separately, more information is contained in their combination , from which a better deprojection is expected to be achieved. This expectation is confirmed by the distribution of r_p for a given value of r, which shows a strong peak at r_p = r (Fig. 2). Assuming a probability distribution P_r(r) of the unprojected variable r, one expects for r_p a probability distribution, i.e., (23)

For example, if the unprojected probability distribution is constant between 0 and r₀, the projected distribution is given by (see Fig. 12) (24)

It does not seem possible to invert this formula analytically. However, it is very possible to construct an algorithm to perform this inversion numerically.

3.2.2 Deprojection method for r_p

Let i_max be the total number of bins in the histogram of r_p values. We denote the observed number of projected values r_p in the bin of rank i. This number is the sum of contributions of 3D values r ≥ r_p from bins of rank j ≥ i – we assume a same number of bins for the r and r_p histograms (see Fig. 12).

Let us derive the probability law for the various contributions. The probability distribution of r_p for a given r value has been shown (Sect. 2.2) to be (25)

for x = r_p∕r. We assume, as an approximation valid for a small enough bin width (i.e., a large enough number of bins), that the probability distribution of r remains constant in any bin. Then the probability density of given r values projected in the bin of rank i and relative width b is written as (26)

Concerning the r values pertaining to the bin of rank j, the normalized x values (i.e., r_p ≤ r) belong to only a total number of j bins, such that the relative bin width (normalized to a total range x = 0 to 1) is b = 1∕j.

Let us now show that the probability law of the contributions of the various r values (given by index j) can be recovered from that of the projected values r_p (given by index i).

Indeed, let us introduce the matrix (27)

for i≤ j. The other elements of this matrix (i > j) are zero. Let P_ij be the transpose of this matrix, i.e., P = π^T.

If the column vector represents the various projected numbers in the bin of rank i, they are obtained from the initial 3D column vector by the matrix product (28)

Therefore the initial unprojected probability distribution is recovered from a mere matrix inversion, i.e., (29)

For example, for three bins the relation between the probability distributions is written as

The deprojected reverse relation is, in this case,

This method is illustrated in Fig. 13, where we deproject a projected distribution N_p = arccosr_p, which is the expected function for a constant initial distribution N_r = cst (in the null bin width limit). The deprojected distribution is constant as expected, except for a small bias mainly involving the first and last bins (which is therefore easy to correct and disappears when the number of bins is increased). Another example is given in the numerical simulation of Fig. 14, where the initial 3D distribution is randomly drawn from a Gaussian probability density (10 000 points). This distribution is nicely recovered from the matrix inversion method. In this case the bias involving the extreme points is unobservable, since these values are almost vanishing.

3.3 Deprojection of mass

3.3.1 Theoretical deprojection

The total mass M = M₁ + M₂ of a pair of objects in relative Keplerian motion is given by (32)

where G is Newton’s constant of gravitation, a is the semimajor axis of the orbit of one object around the other, and T is the periodand V = 2πa∕T.

However, there are two drawbacks that deteriorate the available information on mass, since we have no direct access to a and V: first, we deal only with instantaneous interdistance r and velocity v; and second these 3D values are themselves projected to r_p and v_z, i.e., . The first drawback vanishes for circular orbits, but it can lead to some additional uncertainty in the elliptical case, since the relation between M and rv² is written as (33)

where e is the eccentricity and ξ the parameter of the orbit, which is such that r = a(1 − ecosξ) and t = (T∕2π)(ξ − esinξ). The value of rv²∕GM fluctuates between 1 − e and 1 + e while its time average varies from 1 to 0.5 when e varies from 0 to 1, respectively; it can be approximated by the relation ⟨rv² ⟩∕GM = cos²(πe∕4) up to some percents.

The mass distribution should therefore be statistically deprojected from the observed products . We established in Sect. 2.4 the expected PDF of the projection of a given value of rv² (see Eq. (20) and Fig. 4).

From this formula, one can theoretically deproject any distribution of . Its analytical integration does not seem to be possible, but, as in the case of the deprojection of r_p, it is very possible to construct an algorithm to perform this inversion numerically. However, this deprojection is expected to be difficult, since the individual projection function vanishes for the original unprojected value rv² (see Fig. 4), while it is constant for v_z (easier deprojection) and shows a divergent peak for r_p (best deprojection).

3.3.2 Deprojection method for

As in the deprojection method for r_p, let us introduce a projection matrix constructed from the projection function p_η (η) given in Eq. (20) and Fig. 4. We divide the η range into N bins. The projection matrix is written as (for its non-null coefficients) (34)

where j = 1 to N and i = 1 to N + 1 − j (and A_ji = 0 for the remaining coefficients).

Then the deprojection matrix is obtained from this matrix by the following transformation: (35)

For example, for N = 3, the projection matrix is written as (36)

and the resulting deprojection matrix is written as (37)

We tested the method by numerical simulations in which we built a probability distribution for v and r (e.g., Gaussian distributions, see Fig. 15), then we projected these to v_z and r_p and computed the resulting distribution.Whatever the initial 3D distribution, the projected distribution is strongly decreasing, as expected from the shape of the projection function p_η(η) (see Fig. 4). This property is the reason why it is so difficult to estimate the pair mass (Faber & Gallagher 1979). However, despite this problem, our new method allows us to recover the 3D PDF of rv² with a fair accuracy (see Figs. 15 and 16).

This yields the total mass M of the system for circular orbits, but in the general case, where r and v are instantaneous values on elliptical orbits, there is, as previously specified, an additional uncertainty to go from the rv² distribution to the M distribution, which may be estimated from an evaluation of the eccentricity distribution. This will be studied in a forthcoming work.

4 Conclusions

In this paper, we developed statistical methods of deprojection of the main physical parameters of pairs of astronomical objects (in particular, galaxy pairs). These methods are needed since, in the extragalactic domain, only one component of the intervelocityand two components of the interdistance are available on only one point of the orbit instead of the six (x_k, v_k) coordinates on the full orbit.

We analytically determined the various PDFs of the projected variables: radial velocity v_z, projected interdistance on the sky-plane r_p, their ratio r_p ∕v_z, which can be used as signature of circular orbits, and , which intervenes in the calculation of the pair mass. These analytical solutions were validated by numerical simulations.

Then we described in detail the deprojection methods obtained by inversion of these projection functions. These methods were conceived in a digitalized way to deal with the real data that will be available in the form of histograms.

In this paper, we prepared the application to effective extragalactic galaxy pair data, and in particular to the two galaxy pair catalogs that we recently constructed, one of which contains 13 000 pairs (Nottale & Chamaraux 2018). These numerical simulations support the validity of our deprojection methods, and they also allow us to determine their error bars.

In a forthcoming work, we will apply these deprojection methods to the study of the physics of galaxy pairs, in particular to their dynamics. Such a study will then benefit from twin improvements to catalogs (membership criteria, data quality, and size) and to methods aiming at statistically recovering missing information.

Appendix A Various probability distributions

A.1 Probability distribution of a product of two variables

Consider two random variables x and y whose probability distribution p(x, y) is known. The probability distribution of the product η = xy is given by Rohatgi (1976) and Glen et al. (2004) (A.1)

It becomes in the uncorrelated case, where we have p_xy (x, y) = p_x(x) p_y(y), (A.2)

A.2 Probability distribution of the inverse of a variable

Let x be a random variable of probability distribution p_x (x). The probability distribution of its inverse X = 1∕x is written as (A.3)

A.3 Probability distribution of the ratio of two variables

From these two relations we easily derive the probability distribution of the ratio of two random variables, ζ = y∕x. We find (A.4)

Appendix B Circular orbits

We give here a criterion for identifying circular orbits (and more generally, all configurations when the vectors v and are perpendicular). We consider the scalar product of the position and velocity vectors, i.e., (B.1)

where ϕ is the angle between the two vectors. This relation is written in cylindrical coordinates as (B.2)

where ϕ_p is the angle between the projected vectors v_p and r_p on the plane of the sky. Since and , the ratio of the two observables v_z and r_p is written as (B.3)

When the velocity vector v is perpendicular to the position vector r, we have cosϕ = 0 and therefore this formula takes the simplified form (B.4)

The statistical distribution of the ratio v_z∕r_p is therefore expected to be very different between the orthogonal case and the general case, providing us with a statistical signature for circular orbits.

Let us establish the theoretical expectation of this distribution.

In the case when r and v are perpendicular, which corresponds to circular orbits and to special positions on elliptic orbits, the distributions of v_z and r_p become highly correlated. In particular, let us show that, in this case, |v_z |∕v ≤ r_p∕r.

Let us set α = (xv_y)∕(yv_x), we have , which may be written as α + 1∕α ≥ 2, i.e., (B.5)

The orthogonality of r and v writes xv_x + yv_y + zv_z = 0, such that (B.6)

Accounting for this relation, the inequality (B.5) becomes (B.7)

i.e., |v_z| r ≤ r_p v, QED.

This means that for given r, v and r_p, the possible values of v_z are no more uniform between 0 and v, but are limited to v_z ≤ v r_p∕r.

References

All Figures