Catalogue of wide binary, trinary and quaternary candidates from the Gaia data release 2 (region $\left\vert b\right\vert>25\,\deg$)

The occurrence of multiple stars, dominantly binaries, is studied using the \textit{Gaia}-ESA DR2 catalogue. We apply the optimized statistical method that we previously developed for the analysis of 2D patterns. The field of stars is divided into a mosaic of small pieces, which represent a statistical set for analysis. Specifically, data input is represented by a grid of circles (events) with radius $0.02\,\deg$ covering the sky in the field of galactic latitude $\left\vert b\right\vert>25\,\deg$. The criteria for selecting candidates for multiple stars are based on two parameters: angular separation and collinearity of proper motion. Radial separation, due to limited accuracy, is used only as a weaker supplementary constraint. Due attention is paid to the accurate calculation of the background, which is a necessary input for evaluating the quality of the candidates. Our selection algorithm generates the catalogue of candidates: $900,842$ binaries, $5,282$ trinaries and $30$ quaternaries.


INTRODUCTION
It is believed that statistics of binaries and multiple stars can provide deeper insights into the formation and evolution of galaxies. Wide binaries may serve as a sensitive probe of the Galactic gravitational potential. Some recent studies suggest that very wide binaries may provide data on the presence of dark matter in the Galaxy (Peñarrubia et al. (2016); Hernandez et al. (2019); Pittordis&Sutherland (2019)). At present, an extremely rich source of data on stars of our Galaxy is provided by the Gaia-ESA mission and collected in the catalogue DR2 (Gaia Collaboration (2018)). The release of the next and more expanded full version (DR3) is expected in the first half of 2022. Several authors have recently studied various aspects of binaries using the Gaia DR1 (Oelkers et al. (2017); Oh et al. (2017)), DR2 (Ziegler et al. (2018); El-Badry&Rix (2018); Godoy-Rivera (2018); Jiménez-Esteban et al. (2019); González-Santamaría et al. (2020); Sapozhnikov et al. (2020); Hartman&Lépine (2020); Tian et al. (2020)) and EDR3 (Early Data Release 3) (El-Badry et al. (2021); Gaia Collaboration (2021)) catalogues. Other important articles examining wide binaries before Gaia include, for example, Cabalero (2009) and Close et al. (1990).
In our previous study (Zavada&Píška (2018) we have developed and applied statistical methods for the identification of binaries and multiple stars with an accurate estimate of the background. The resulting catalogue is based on the 3D event analysis and contains about 8 × 10 4 binary candidates at a distance of up to 340 pc. This catalogue involves the binary stars whose total separation does not exceed the event sphere diameter that was 4 pc. Because determining the radial separation Z ij can have a significant error, we have been selecting the binary stars in the events based on the distribution of the projected distances ∆ ij , the error of which is only slightly affected by the parallax error. At the same time, the error of parallax and radial distance causes that the constraint on radial separation Z ij ≤ 4 pc can exclude a large number of true binary stars. In the present study, we solve this drawback of 3D events by a more general procedure. The procedure combines the analysis of angular separations d ij of sources inside 2D events with the information on the radial separation Z ij . Simultaneously, we apply the condition of approximate collinearity of proper motion here as before. The result of the procedure is the selection of candidates for the binary with the defined probability that it is the true binary and not a random background.
Our 2D analysis is described in Sec.2 and consists of several parts. The basic notions, which our method deals with are defined in section 2.1. Input data from the sector of the Gaia catalogue are defined in section 2.2. In the next section 2.3 we analyze the probability of the binary depending on the angular separation and collinearity of proper motion of its components. This analysis sets the rules for the selection of binary candidates for our new catalogue. Using parallax, we can define the distance of binary and the projection of the absolute separation, which is more important for physics than just the angular separation. In section 2.4 we show that distributions of the projection and correlations with proper motion allow us to obtain information about binary orbit. In section 2.5 we show the results of our analysis on the occurrence of trinaries and quaternaries. Candidates with a high degree of reliability are listed in the electronic catalogue. Its structure and content are described in Sec.3. The comparison with the other catalogues (El-Badry et al. (2021); Hartman&Lépine (2020); Jiménez-Esteban et al. (2019); Zavada&Píška (2020)) is done in Sec.4. The last Sec.5 is devoted to the summary and concluding remarks.

Methodology
The method of 2D analysis is described in detail in our previous paper Zavada&Píška (2020), so here we repeat only basic notions.
The data for analysis are represented by the grid of circles with patterns of stars covering a defined region of sky ( Fig.1). The input data for generating the grid are given in the galactic reference frame. So, the position L of a source is defined by spherical coordinates L, l and b (distance from the Sun, galactic longitude and latitude): In the centre of each circle we define local orthonormal frame defined by the basis: k l = (− sin l 0 , cos l 0 , 0), where k r = n 0 (b 0 , l 0 ) defines angular position of the circle centre. Unit vector k l is perpendicular to k r and has direction of increasing l. Unit vector k b is defined as k b = k r × k l and has direction of increasing b. Vector k r has radial direction, perpendicular vectors k b and k l lies in the transverse plain. The stars inside the circle of small radius ρ 2 satisfy |n i − n 0 | ≤ ρ 2 , i = 1, ...M. (3) where {x i , y i } are local rectangular coordinates defined by the basis (2): The circles in the grid are, in fact, spherical caps. However, their radius ρ 2 will be so small ( 3.5 × 10 −4 rad), that the caps can be reliably considered as flat circles. A set of stars defined by (4) is called the event. We define the pair separations It is useful to define the scaled separationξ where ρ 2 is the angular radius of the events. Distribution of separations of randomly distributed sources inside the circle does not depend on M and is given by relation which was proved in (Zavada&Píška (2020)). This curve is important for subtraction of random background. Before practical use, we add a few general remarks: i) Of course, the grid in Fig.1 can be rectangular only locally. In fact, we work with circles with the same radius aligned only in rows with a constant galactic latitude.
ii) The radius ρ 2 must be set to be significantly larger than a typical separation of true binary. At the same time, it must be so small that the distribution of stars within the event can be considered random and uniform. In Sec. 2.3 these conditions will be explained in more detail. We will also explain that the obtained results are not sensitive to the exact setting of ρ 2 .
iii) Events with too high multiplicity M , in which various higher dense structures may dominate, are excluded from processing.
iv) The circular shape of the events is chosen due to a simple but accurate formula (8) for calculating a random background. Another shape would lead to a more complex formula depending on other shape parameters (triangle, square, orientation ...). On the other hand, the grid of circle events only partially covers the sky. The complete set of binaries is obtained in Sec.3 from several shifted grids.

Input data
The parameters of the events we work with are listed in Tab.1. By symbol R we denote the union of both regions in table: In Fig.2 we have shown distribution of the distances of the stars in the table from the Sun (origin of the galactic reference frame). If the parallax p is given in angular units as (= 1 ), then the corresponding distance is .
2500 5000 7500 10 000 L[pc] 0 400 000 800 000 1.2 × 10 6 P(L) The figure demonstrates the scale of distances measured by Gaia. This histogram also suggests the convention valid for all histograms P(x) in this paper: ordinate represents the number of entries of variable x of given value into the bins of the width specified in the figure caption. Stars in the sky are divided into the circle regions (events) of angular radius ρ 2 and in Fig.3 we show distribution of the star multiplicities in these events. The distribution of stars in region R 1 is roughly homogenous, so the multiplicity distribution has a nearly Poisson shape Since region R 2 shows greater density fluctuations resulting in a broader distribution then Poisson, we fail to make its fit. This region is dominated by an inhomogeneous stellar density due to patchy inter-stellar extinction, the presence of open clusters, and the imprint of the Galactic spiral pattern. The proper motion of the stars in DR2 is defined by two angular velocities in directions of the right ascension and declination in the ICRS. So the corresponding transverse 2D velocity U is given as For the pair of stars we can define the angle between both transverse velocities: Our selection of binaries is based on The second condition was applied also in our former 3D analysis. Panel b shows distribution P (α ij ) in the band d ij ≤ 15 as, similarly in panel c we have distribution P (d ij ) for α ij ≤ 15 deg. Recall that the stars separated by d ij 0.5 as are missing due to current resolution limit in the DR2 data set as stated in (Arenou et al. (2018)). In Fig.5 we have shown distribution of radial separations Z ij defined as and obtained with the use of parallax (10). Panel a(b) shows distribution of pairs outside (inside) the window (15). The sharp peak at small radial separations in the second panel is expected for binaries, the tail with greater separations corresponds partly to the background pairs and partly to the true binaries having large errors of parallaxes. Alternative distributions P (∆L/L), where are shown in Fig.6a,b. Relation (10) implies In Fig.6c we show distribution P (∆p/p) obtained directly from the Gaia data, where the parallax of the source is accompanied by its estimated error. Our ∆L, L related to binaries are also obtained from parallaxes, however, our estimate of relative error of the parallax is independent and can serve as a cross-check. We observe a noticeable similarity between the distributions b,c. A smaller difference can arise because the window of binaries also contains a background with random radial separations.

Selection of binaries and calculation of background
In the next we will work with the numbers: n 2 − number of (true) binaries b 2 − number of background pairs N 2 − total number of pairs: These numbers are related to given sample of pairs defined by the corresponding cuts. One cannot decide if the pair is a binary or background, but we can find out the probability that the pair is binary (n 2 /N 2 ) or background (b 2 /N 2 ). Obviously, in the domain of the peak in Fig.4a the probability of binaries n 2 /N 2 is high. Radial separation Z ij can be used to further increase the probability of binaries in the peak, but due to low accuracy, the cut must be set judiciously. Too strict cut generates a cleaner sample of binaries (higher ratio n 2 /N 2 ) but more binaries are excluded. And conversely, too soft cut preserves more binaries, but at the price of the higher background and lower ratio n 2 /N 2 .
The area d × α of Fig.4a can be divided into four the windows A,B,C and D (like in Fig.8): which is the domain of the peak (15) with the high population of binaries. The remaining windows are Distribution P (d ij ) can be represented by its normalized form The random background is described by the normalized function q(ξ) defined by Eq.(8). Both distributions are shown in Fig.7a for region R 1 . Due to equal normalization, an excess in the peak of binaries in distribution P is compensated by the lack of pairs in the region of background. Assuming that for d ij ≥ 2/3d max orξ ≥ 2/3 distribution P involves only background, we renormalize q correspondingly: In panel b we have shown distribution Q(ξ) together with the renormalized background q r (ξ). The same distributions for region R 2 is shown in panel c. Obviously, the distribution of binaries is given by their difference It is now clear from these figures why the diameter of events 2ρ 2 is chosen much greater than typical separation of the binaries d bin (the parameterξ bin = d bin /2ρ 2 0.1 in the figures). The larger diameter allows us to more accurately determine the weight of the background curve to be subtracted. At the same time we observe the distribution of separations is perfectly random outside the region of binary peak. The probability that the pair in the event is the true binary and not a random background is given by the ratio where one can replaceξ → d ij = 2ρ 2ξ . Examples of the function β will be given below. The shape of the background distribution q depends onξ only. In the background, we can see an increase or peak at lowξ, which is a signal of the presence of binaries. An additional selection of pairs with small α ij or a small radial separation Z ij increases the dimension of the peak on the background, but the shape of the background curve does not change. For instance, Fig.7b,c corresponds to α ij ≤ 15 deg . The reason is simple: in the sky without binaries (only background can be observed) there is no correlation between angular separation d ij and the parameters α ij or Z ij . From the data we have known the numbers of pairs in four windows defined above: N A 2 , N B 2 , N C 2 and N D 2 . The number N A 2 is called the number of binary candidates. According to an assumption above, we have also numbers of background pairs obtained forξ ≥ 2/3 in domains If we denote then we can calculate the numbers of background pairs in the windows A-D as Therefore, the numbers of binaries in the respective windows read These numbers, obtained in both regions R 1 and R 2 under different conditions are given in Fig.8. The numbers in upper panels A, B (the set All events) follow from panels b,c in Fig.7. The peaks of binaries are evident, however the area under the background curve is also considerable. Correspondingly, the quality ratios β A = n A 2 /N A 2 (in the second row of the panels A), representing the probability of the true binary in window A are not satisfactory, particularly for the region R 2 . Parameter β A can be increased by applying the additional cut on radial separation ∆L. The condition ∆L ≤ ∆L max = 500 pc (30) applied in the region R 1 gives "enriched" sample with the quality ratio β A = 0.77 (left middle panel). In the region R 2 the situation is more complicated. The density of stars is higher, so the background grows. Moreover, the average density varies with longitude. So, we have divided the region into three subregions with different ∆L max cuts, which give a similar β A in these subregions. Their definition is shown in Tab.2. These cuts give resulting β A and other parameters related to R 2 listed in right middle panel. One could further squeeze the cuts to obtain a cleaner sample of binaries (higher ratio n A 2 /N A 2 ), however at the price that more true binaries are excluded (n A 2 is smaller). This is illustrated by the numbers in upper and middle yellow panels in the figure. The cuts exclude not only background pairs but also true binaries in window A. It is due to a rather poor accuracy of radial separation that we have shown in Fig.5b. If we had an accurate radial separation, then a suitable cut would suppress only background and preserve true binaries.
The situation is more favorable with brighter stars. The lower panels in Fig.8 show that for binaries of magnitude G ≤ 15 the quality ratio β A is very good even without any cut on radial separation. The quality of this sample is illustrated also by Fig.9, where we observe the high binary peaks with low background. In any case, our criterion is based on the resulting ratio β A , so we can accept candidates with greater parallax uncertainty if β A is greater than requested. In other words, if the signature from d ij and α ij is sufficient, then the parallax is not important.
Equivalently, abundance of binaries depending on their separation is defined by β function (25). In Fig.10a,b,c we have shown this function for some subsets of regions R 1 and R 2 defined by Fig.8. Figure Fig.10 demonstrates that higher β A could be also reached by squeezing the parameter d c in window A (20). Note, the effectively wider peak for magnitude G ≤ 15 is due to a lower level of background pairs, which can be seen from the comparison of Fig.9 with Fig.7b  We have shown that the function β is a key to determining the quality and quantity of binaries in the selected subset. It is, therefore, necessary to verify that this function does not change with the choice of the purely technical parameter ρ 2 . Fig.10d shows the function β calculated for the sub-sample of enriched exposition with event radius ρ 2 = 36as (half the radius of standard events). The functions in panels b,d agree perfectly. A similar agreement was also verified for Figs.12 -14 below.
What is the effect of input data uncertainties on the selection quality? The initial distribution P (d ij , α ij ) in Fig.4 consists of two parts corresponding to the true binaries and the background. The shape of the distribution P bin (d ij , α ij ) is affected by the measurement errors. In general, these errors cause a widening of this distribution, but its integral n 2 does not change, as suggested in Fig.11. However, the errors will affect the value of the quality ratio in the domain of the peak. As suggested in the figure, wider P bin means larger b 2 , which implies a smaller ratio β P . In this way, the measurement errors of the selection parameters directly affect the β function, larger errors mean worse β. At the same time, larger errors expand the space for selection. A similar result follows from the algorithm in (Halbwachs (1986)) for the selection of common proper motion pairs. At the end of this section, we summarize the selection algorithm of binary candidates: 1. In the distribution P (d ij , α ij ) in Fig.4 we observe the peak on a smooth background. The peak is a sign of the presence of binaries. We take the edges (15) of this peak as the limits of our selection of binary candidates.
2. The number of true binaries can be accurately determined by subtracting the random background as shown in Fig.7b,c. The background curve is given by the function (8). This curve corresponds to the random distribution of sources within a circle. That is why we work with circular events and with the distribution of angular separations of the pairs.
3. The background level in the peak domain can be further reduced by the cuts on radial separation ∆L max . η P bin,bg (η) Figure 11. Distributions of true binaries P bin (red and blue, which is wider, due to measurement errors) and background P bg (green). Parameter η represent some selection parameter (dij, αij, ...). The numbers of the true binaries n2 are given by the area under red and blue curves. Corresponding numbers of background pairs b2 are given by the yellow area and the whole area under the green curve. 4. The cuts α c and ∆L max reduce the background level, but does not change the shape of the corresponding background curve q(ξ). This distribution, after the renormalization (10), is used for calculation of the probabilistic function β, which defines the selection quality.

Projected absolute separations, periods and masses of binaries
In Fig.12 we have shown an important correlation in the peak of binaries, where angular separation d ij , as one would expect, strongly correlate with the projected distance ∆ ij that is calculated as Distribution of ∆ ij in the region of peak (windows A in Fig.8) is shown in Fig.13. In panels b,c we have shown distributions from subsets with higher rate of binaries. We observe that which confirms our former result from 3D analysis. In panel a we observe presence of pairs, which contradicts this constraint. It is due to high rate of background pairs in this set. Note that nearby binaries with wider separation can be outside of our acceptance window A. In fact only candidates with separation where L is average distance of the pair, are accepted. Similarly as before (Zavada&Píška (2020)), we can estimate the projection of the orbital velocity of the binary as From the plot in Fig.14a, we can roughly estimate averages of the orbital periods and total masses of the binary systems. The binaries are accumulated in a peak very similar to that obtained earlier in 3D analysis. Very similar plot (with much lower statistics) can be obtained also for the subset of binaries with magnitude G ≤ 15. If we take the sources roughly in the domain of half-width of the maximum (panels b,c) ∆ ij ≤ 0.01 pc, ν ij ≤ 2.5km/s (37) then in the approach described in (Zavada&Píška (2020)) we now obtain which is comparable with the previous approximate estimate (8 × 10 4 y, 0.8M ) obtained in the cited paper. These numbers can also be compared with the results of the sample of thirty CPM (Common Proper Motion) pairs, which were discussed in (Duquennoy&Mayor (1991)), Table 4. These pairs represent wide binaries in the solar neighborhood, for which the table lists the parameters log T i and masses M 1i and M 2i . We have calculated the mean values T = 10 log Ti ≈ 2.6 × 10 4 y, which are also well comparable with our estimates above.

Wide trinaries and quaternaries
So far we have assumed the excess of close pairs (window A) is due to only binaries, which is correct only in a first approximation. A more detailed analysis indicates also a limited presence of multiple systems, trinaries and quaternaries. To identify their candidates, we will generalize our criterion of binary candidates. For each pair in the candidate multiple system we require the same condition as for separation of binaries (20 which gives ≈ 1.6% in the set of enriched exposition. For each triplet of stars inside the event we define the triangle separation as If αβ are subscripts of the most separated pair, then the components of triangle separation are defined as Distributions of these parameters are shown in Fig.15 together with the distributions corresponding to the background, which is generated by randomly distributed sources in the event. The high level of the background under the slight trinary peak at d ijk 20 as (upper panel left) is due to primarily by binaries. Each binary generates also excess of triangle separations. It is illustrated in Figs Fig.8 gives the total number of binaries roughly 1.5 × 10 6 . However, as we discussed in (Zavada&Píška (2020)), part of them (mainly in windows B,C) may be an image of widening pairs that were less separated but weakly bound in the past.

CATALOGUE
In this section, we describe the catalogue of the multiple star candidates, which are selected with the use of the events defined by Tab.1. For the present version of the catalogue, we accept only the candidates from window A(R) of the enriched exposition defined in the middle part of Fig.8. So, the selected pairs satisfy the conditions of angular separation, collinearity and radial separation: where ∆L max is defined in relation (30)    part of the sky (corresponding to the fraction π/4). We also lose candidates between neighboring events when the pairs are split between the neighbors. To recover these losses, we work with modified coverage: i) The event circles of radius 72 as are replaced by squares of edge 144as with no gaps between them. In each square, we search for multiple star candidates.
ii) The procedure is repeated with the same squares centred in the corners and the edge centres of the former squares (we have 4 grids in total). Then the search results are merged, the summary data are given in Tab.4. The quality ratios β = n A m /N A m are taken from the analysis of circle events with well defined background ( Fig.8 and Tab.3). Now we can use them to estimate n A m and corresponding statistical errors It is seen the number of multiple systems decreases rapidly with their multiplicity: In Fig.19 we have shown distribution of distances of all stars accepted to the catalog. Note that selection criteria exclude preferably the most distant sources from distribution in Fig.2. Obviously, the candidates of higher quality (but lower quantity) can be obtained from the catalog by reselection with more strict cuts (43). The catalogue of selected candidates is represented by a matrix, which is defined as follows. Each row represents one star and there are the following data in the columns: 1-2: Group ID and Group size (n = 2, 3, 4) to match stars with the group they belong to. The quality of our catalogue can be easily verified for brighter binaries, where the DR2 data includes radial velocities. The catalogue contains 6469 such pairs, only 16 of them involves a star with magnitude G > 15. In Fig.20 we show the correlation of both velocities, which confirms the common radial motion of the pairs. This result agrees very well with the same correlation obtained from the DR2 catalogue in (Andrews et al. (2018)) and (Sapozhnikov et al. (2020)).

COMPARISON WITH OTHER CATALOGUES
The catalogues of wide binaries (El-Badry et al. (2021); Hartman&Lépine (2020); Jiménez-Esteban et al. (2019); Zavada&Píška (2020)) and our present one are based on the selection defined by three parameters measured by the Gaia: angular position (separation), proper motion and parallax. However, catalogues differ in the choice of algorithm, which processes these variables. To illustrate, we will compare the algorithm used in this study (A1) with the algorithm used for the recent catalogue (El-Badry et al. (2021)) containing a comparable amount of binaries (A2).
i) angular separation  This parameter is measured with a high accuracy. The presence of binaries generates a clear peak in its distribution at small separations. The random background in the events is exactly described by the function (8). The peak in angular separations is the basic signature of binaries in the algorithm A1.
From angular separation one can, with the use of parallax, calculate projected separation. Both parameters are strongly correlated. The cut on projected separations ∆ ≤ 1pc is applied in the algorithm A2. However, the question of the background is here more complex, see discussion below.
ii) proper motion In algorithm A2, the difference of velocities v ij (36) is limited by the cut The algorithm A1 is based on conditions (15), which define the domain of the peak with the high population of binaries (Fig.4). The algorithm includes an accurate calculation of the background under the peak. What are the velocities v ij inside and outside this peak? The answer is given in Fig.21. The left panel shows the distribution of velocities for all pairs within the event, the peak at small v ij corresponds to the binaries. The right panel shows the same distribution for the pairs in the peak domain. Corresponding pairs are mostly binaries, so the large velocities are strongly suppressed. The panel represents the distribution of the orbital motion projection. One can admit, the tail of distribution may still be contaminated with some background pairs. Alternatively, instead of our cut it would be possible to set the cut (46). In the present study, we preferred (47), because the angle α ij does not depend on the parallax, which has a large measurement error. For future analysis of the Gaia DR3 data, in which a more accurate parallax measurement is expected, we plan to try replacing (47) with (46), which can more effectively suppress the background. iii) parallax The parallax p or the distance L = 1/p are determined with a large error, ∆L is usually larger than the expected binary dimension. Large radial separation does not necessarily mean that the pair is not true binary. However, a smaller radial separation (or difference of parallaxes) increases the probability that the pair is true binary. Therefore the algorithm A1 requires limited radial separation defined by (8) and Tab.2. Similarly, the A2 requires limited difference of parallaxes.
Finally, both algorithms similarly reject denser clusters of sources: event multiplicity M > 25 (A1) and number of neighbors > 30 (A2). In general, algorithms have a similar philosophy but differ in technical details.
A very important technical step is to define random background (A1) or equivalently the chance alignments (A2). In this respect, the two algorithms differ significantly. Let's make a comparison.
A1: The distribution of separations of random sources inside a circle is given exactly by formula (8). The shape of this curve does not depend on other parameters such as magnitude, the direction of proper motion or the parallax. The binary peak can be separated from the background described by this curve, see Sec.2.3 and examples in Figs.10   and 10. The normalized background curve is q ξ ≈ 8ξ + ...
for smallξ (peak region). The probability of finding binary in a selected subset is given by function β defined by ratio (25), examples are shown in Fig.10. A2: From the Gaia EDR3 data input, there are produced two files involving pairs with seven parameters x = {angular separation, distance, parallax difference uncertainty,...}: a) Catalogue of candidates b) Catalogue of chance alignments (shifted catalogue) The densities of pairs N candidates (x) and N chance align (x) in the seven-dimensional parameter space are approximated by the Gaussian kernel density estimates. The ratio of these approximations is the parameter, which provides classification of the quality of candidates: the low R means a high probability that the candidate is true binary, as illustrated in Fig. 5 in El-Badry et al. (2021). In the right panel of this figure, we observe: the lower R (and higher quality of the candidate) means stronger suppression of more separated pairs. It is the result, which correlates with the shape of probabilistic function β in A1. However, the ratio R of both approximations does not strictly mean probability, which the authors admit. We will compare the contents of the catalogues listed in Tab.5. We define three types of candidates for any two compared catalogues Ai, Aj: 1) unique -the candidate appears in only one catalogue N uniq i 2) identical -the candidate appears in both catalogues N iden ij 3) indefinite -cannot be decided, for example, two candidates from two catalogues have only one common star, or the pair in one catalogue is part of a greater system in another one.
The results of the comparison are shown in Tab.6 In the section above the diagonal, there are the numbers N iden ij , below the diagonal we have the numbers N uniq i \N uniq j . In the catalogue A2, only such pairs can be used for comparison, where the IDs of both stars appear also in the Gaia DR2 data. We perform the comparison of A2 only with our catalogue (A1, |b| > 25 deg), so the other corresponding places in the table are empty (x). Unique candidates are based on different selection conditions in A1 and A2. For example, in A2, the unique part is generated mostly by candidates with greater angular separation (d > 15 as) than is accepted in A1. On the other hand, A1 imposes much weaker constraints on parallaxes, which generate its unique part. That is why in Fig.19 we observe many candidates more distant than 1000pc (upper limit in A2). The occurrence of trinaries and quaternaries is analyzed only in A1. So, the content of both independent catalogues is partly complementary and partly identical.
Our current analysis is limited to |b| > 25 deg because we have verified that in the dense region |b| ≤ 25 deg, the efficiency of selection (ratio β A ) based on conditions (15) from the DR2 catalogue is even lower than in the region R 2 . And additional cuts on radial separation (similar to Tab.2 ) are not sufficiently effective. With the expected DR3 data release, where higher accuracy of astrometric data (mainly of the parallax) is assumed, we plan to recalculate the selection of multiple stars in the full angular range (4π).
The content of our previous catalogue A4 has already been compared with the A5, which contains 3, 055 binaries with magnitude G ≤ 13 (Zavada&Píška (2020)). Comparison A1 with A5 implies that only 381 of the A5-candidates belong to A1-region |b| > 25 deg and only about 60 of which meet selection criterion d ≤ 15 as. As shown in Tab.6, similar comparisons were made also for other catalogues. In general, the only partial overlap of different catalogues is due to mainly different cuts in selections algorithms.
For example, the content of present catalogue A1 is compared with our previous version A4, where the binary candidates inside the surrounding cubic region (400 pc Fig.18e that is excluded for the new catalog, but still one pair can meet the criteria of the previous one. At the same time, within the part of the cubic region that is common to both catalogues, the number of binary candidates of the new catalogue is 79, 771. It is 79, 771/25, 604 ≈ 3 times more, than the candidates in the previous one. This ratio proves high efficiency of the optimized method applied in the present analysis.
Our new catalogue A1 covers the region |b| > 25 deg up to the distance 3, 000 pc and contains about 10 times more candidates than the previous A4, which covered full 4π, but only up to 340 pc. A total of almost 10 6 candidates are recorded in both merged catalogues. For more detailed comparison of the catalogues A1, A3, A4, A5 and A6, which are based on the Gaia DR2 data, we have created the merged catalogue. This catalogue is the list of pairs consisting of four items: order number of binary, mask 103456 denoting origin from A1,A3,A4,A5,A6 and two DR2 sources ID. The catalogues A1, A4 and the merged one are available in the csv form on the website https://www.fzu.cz/ ∼ piska/Catalogue/.

SUMMARY AND CONCLUSION
With the use of the optimized statistical method for analysis of 2D patterns we studied occurrence of wide multiple stars: binaries, trinaries and quaternaries. The candidates are selected using the astrometric data collected in the Gaia-DR2 catalogue. So, we have studied the pairs with angular separation wider than ≈ 0.5 as, which is present Gaia lower limit for resolution of two sources. Our present analysis covers the region of galactic latitude |b| > 25 deg and radial distance L 15, 000 pc. In this space we have analyzed about 1.3 × 10 7 circle events of diameter 144 as involving 7.5 × 10 7 sources. The circle shape is advantageous for calculation of the random background. The total number of processed sources with positive parallax is about 1.2 × 10 8 .
The analysis is focused on two basic parameters related to any pair of sources in the multiple systems: angular separation d and collinearity α of their proper motion. Distribution of these parameters is compared with distributions generated by the random background. The domain of the clear peak of binaries is limited by conditions (15) that serve as the cut for the selection of candidates. The exact knowledge of the background allows us to define probabilistic parameter β representing the quality of candidates. Additional condition, which is required for the radial separation ∆L, improves quality of candidates. After this selection (enriched exposition) the candidates are recorded in the attached catalogue. Total numbers of candidates of wide multiple stars in the region |b| > 25 deg are shown in Tab.4. The catalog is compared with some other catalogues of wide binaries selected from the Gaia data.
We have also shown that the results of the present 2D analysis are fully consistent with our previous 3D analysis of Gaia DR2 data. We confirm that the projection of wide binary orbit satisfies approximate relation (34). The average period and mass of the wide binary systems (38) are very similar to our previous estimates based on 3D analysis. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa. int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/ web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The work was supported by the project LTT17018 of the MEYS (Czech Republic). We are grateful to A.Kupčo for the critical reading of the manuscript and valuable comments. Further, we are grateful to J.Grygar for his deep interest and qualified comments and O.Teryaev for useful discussions with interesting ideas. We thank S.A. Sapozhnikov for providing the link to their catalogue that we used for comparison.