1. INTRODUCTION

As the gravitational fields are highly nonuniform during the formation of galaxies and their clusters, the distribution of the former plays an important role. Moreover, the character of the gravitational field distribution permits to discern the specific scenario of galaxy formation. The classical scenarios of such formation were proposed a long time ago (^{Peebles 1969}; ^{Zeldovich 1970}; ^{Sunyaew & Zeldovich 1972}; ^{Doroshkevich 1973}; ^{Shandarin 1974}; ^{Efstathiou & Silk 1983}; ^{Dekel 1985}) and dealt primarily with the mechanical instability of the so-called Zeldovich pancake model (^{Zeldovich 1970}), based on a large gravitating body (of pancake shape, hence the name), without any randomness in the gravitational fields of the constituents (^{Zeldovich 1970}; ^{Shandarin & Zeldovich 1989}; ^{Longair 2008)}. These classical scenarios have not lost significance since new scenarios are essentially modifications of the classical ones and can be classified according to them. The improved form of the above classical scenarios has been put forward more recently (see ^{Shandarin et al. 2012}; ^{Giahi-Saravani & Schäfer 2014} for relevant references). The presence or absence of the random gravitational field fluctuations contributes to the studies of the angular momentum acquisition of galaxies during their formation stages. This may in principle permit to specify the most probable (among many other) scenario of the emergence of large stellar objects. This is because the final test of the correctness of a a given scenario is the comparison of its predictions with observations. The investigations of the variations of the angular momenta of structures gives us the opportunity to do so. Note that different scenarios make different predictions about galaxy orientations, i.e. the angular momenta alignments in structures (^{Peebles 1969}; ^{Doroshkevich 1973}; ^{Shandarin 1974}; ^{Efstathiou & Silk 1983}; ^{Catelan & Theuns 1996}; ^{Li 1998}; ^{Lee & Pen 2000}, ²⁰⁰¹, ²⁰⁰²; ^{Navarro at al. 2004}; ^{Trujillio et al. 2006}; ^{Zhang et al. 2013}). For this reason the analysis of the orientation of the galaxy planes is regarded as a standard test of different scenarios of cosmic structure formation (^{Romanowsky & Fall 2012}; ^{Joachimi et al. 2015}; ^{Kiessling et al. 2015)}.

One of the natural sources of randomness of gravitational fields is the mass distribution of stellar objects. The simplest possible model of the mass distribution dates back to ^{Chandrasekhar (1943)}, where the a priori given distribution function of masses τ(M) has been considered. In this paper, all characteristics of the stellar ensemble are expressed through the different average powers of mass. The mass averaging was performed implicitly with the above function τ(M). In principle, such an approach can be generalized to averaging over quadrupole and higher multipole moments of galaxies in the spirit of ^{Stephanovich & Godłowski (2015)}. Although this effect may change some of the results quantitatively, we speculate that its overall influence will be rather weak. The next step was done by ^{Press & Shechter (1974)}, who considered the distribution function of stellar object masses within the model of a self-similar gravitational condensation. In this work, within the model of the expanding Universe in a Friedmann cosmology, the stellar ensembles were represented as a “gas” of self-gravitating masses, which could condense into aggregates with larger mass, forming finally very large clumpy objects. This model permits to derive the distribution of masses in the form

where A is a normalization constant, (see below). Note that the explicit expression of the Shechter function (1) has been listed in ^{Shechter (1976)}, a work devoted to the luminosity distribution of galaxies. Below we will use function (1) for the calculation of the amended (on the mass distribution) distribution function of the gravitational fields and angular momenta. Using this function, we assume that the mass is proportional to the first power of the luminosity: m ∼ L. Below we give arguments why a possible nonlinearity m ∼ L ^γ (γ 6= 1) will not change our results qualitatively.

In the present paper, we consider tidal interactions in the ensemble of galaxies and their clusters in a Friedmann-Lemaître-Robertson-Walker Universe with a Newtonian self-gravitating dust fluid (p = 0) containing both luminous and dark matter. The commonly accepted model of such a Universe is the spatially flat homogeneous and isotropic ΛCDM model. The clumpy objects like galaxies and their clusters are formed as a result of almost scale-invariant Gaussian fluctuations (^{Silk 1968}; ^{Peebles & Yu 1970}; ^{Sunyaew & Zeldovich 1970}). This assumption is the base of the so-called hierarchical clustering model (^{Doroshkevich 1970}; ^{Dekel 1985}; ^{Peebles 1969}). The models with non-Gaussian initial fluctuations have also been considered in ^{Bartolo et al. (2004)}. The non-Gaussian character of the distribution function was postulated there, rather than calculated. A calculation was presented in ^{Stephanovich & Godłowski (2015}, ²⁰¹⁷⁾, where the non-Gaussian distributions of gravitational fields and momenta were calculated using the method of ^{Chandrasekhar (1943)}. Here we generalize this calculation considering the mass distribution (1). Note that the calculations made in ^{Stephanovich & Godłowski (2015}, ²⁰¹⁷⁾ dealt with the equilibrium situation only. To consider a non-equilibrium situation, it is necessary to use differential equations of Fokker-Planck type with so-called fractional derivatives (^{Garbaczewski & Stephanovich 2009}, ²⁰¹¹⁾. In this case we can begin with an ubiquitous Gaussian distribution and arrive at a non-Gaussian one as a result of a primordial, fast time evolution. After it, the slower evolution, dictated by the ΛCDM scenario, takes place. Note that recently time evolution of intrinsic galaxy alignments has been found by ^{Schmitz et al. (2018)}.

In the hierarchical clustering approach, the large clumpy structures form as a result of gravitational interactions between smaller objects. In other words, the spin angular momenta of galaxies arise as a result of tidal interactions with their neighbors (^{Schäfer 2009)}. Note that, in the present paper, the angular momentum is the result of tidal interactions¨ with the entire environment, which occur via interaction transfer from close to distant galaxies (see below). That is to say, our approach is the generalization of ^{Schäfer & Merkel (2012)}; ^{Catelan & Theuns (1996}, ^1996a); ^{Lee & Pen (2002)},¨ where the average tidal interaction with the entire environment has been considered. In the present work we perform the theoretical and statistical analysis of the influence of tidal interactions between astronomical objects on the formation of larger (then initial) structures. We also perform the comparison of our model predictions with vast data arrays, derived from observational and numerically simulated data. It turns out that our theoretical results are in good coincidence with the above observational and numerical data. Our theoretical model includes an additional mass distribution, obeying the Shechter function (1). It turns out that the mass distribution (1) does not change our main result (^{Stephanovich & Godłowski 2015}, ²⁰¹⁷⁾ that in stellar systems with multipole (tidal) gravitational interactions, the distribution function of gravitational fields cannot be Gaussian. The crux of the matter here is the long-range character of the Newtonian (and derived multipole) interaction between stellar objects. Such a character implies that distant objects (like galaxies, their clusters and even dark matter haloes) still “feel each other”, which is not the case for a Gaussian distribution. The derived non-Gaussian distribution function allows us to calculate the distribution of virtually any observable (like angular momentum) of the astronomical structures (not only galaxy clusters, but also smooth components like halos, which massdominate the total mass of a cluster (see ^{Kravtsov & Borgani 2012)} in any Eulerian approach (linear or non-linear).

The relation between the angular momentum of galaxy clusters and their mass has also been investigated observationally. It is not difficult to analyze the distribution of angular momenta for the luminous matter. In the real universe, the luminous galaxies and their structures are surrounded by dark matter halos. These halos are often much more extended and massive than the luminous component of the structures. Unfortunately, direct observation of dark matter halos and their angular momenta is complicated. One should not forget, however, that there are observed correlations between luminous and dark matter (sub)structures. This implies the certain dependence between dark matter halos and luminous matter (real galaxies) orientations (^{Trujillio et al. 2006}; ^{Paz at al. 2008}; ^{Pereira et al. 2008}; ^{Bett et al. 2010}; ^{Paz et al. 2011}; ^{Kimm et al. 2011}; ^{Varela et al. 2012)}. Recently, the results of ^{Okabe et al. (2018)}; ^{Codis et al. (2018)} based on Horizon-AGN simulation show a similar dependence. This allows to conclude that the analysis of angular momentum of luminous matter gives also information about the angular momentum of the total structure (i.e. that with dark matter halos). As a result, the analysis of the angular momentum of “real” (luminous) galaxies and their structures, is still important as a test for a possible structure formation scenario. Note that investigations of galaxy orientation in clusters are also very important for the analysis of weak gravitational lensing (see ^{Heavens et al. 2000}; ^{Heymans et al. 2004}; ^{Kiessling et al. 2015}; ^{Stephanovich & Godłowski 2015}; ^{Codis et al. 2016} for more details).

As generally galaxy clusters do not rotate (^{Hwang & Lee 2007}; ^{Tovmassian 2015)}, the angular momentum of a cluster is primarily due to the spins of member galaxies. Unfortunately, usually we do not know the angular momenta of galaxies. So the orientations of galaxies are investigated instead (^{Oepik 1970}; ^{Hawley & Peebles 1975}; see ^{Romanowsky & Fall 2012}; ^{Pajowska et al. 2019} for recent reviews), since it is assumed that the rotational axes of galaxies are normal to their planes. Such an assumption seems to be quite reasonable at least for spiral galaxies. As a result, a stronger alignment of galaxies in a structure means a larger angular momentum of the latter.

The question is whether there is any relation between the alignment and mass of the structure. The general result of the previous papers is that there is insufficient evidence for galaxy alignment in less massive structures like groups and poor clusters. However, we do observe the alignment of galaxies in rich clusters; see ^{Godłowski (2011)}; ^{Pajowska et al. (2019)} for a review. The first results (^{Godłowski, Szydłowski, & Flin 2005}; ^{Aryal et al. 2007)} were qualitative only. Because of that, ^{Godłowski et al. (2010)} and ^{Godłowski (2012)} investigated quantitatively the orientation of galaxies in a sample of 247 rich Abell clusters, using the improved ^{Hawley & Peebles (1975)} method (see ^{Pajowska et al. 2019} for the latest review). In these papers, it was found that the alignment is present in the above sample. Moreover, galaxy orientation increased with richness of the cluster. However, the data were insufficient both to resolve the question about the exact form of this relationship and to confirm the hypothesis that the angular momentum of the structure increases with time. This is the reason why we decided to extend our sample and to compare observational results with those obtained from the simulations.

2. THE FORMALISM

Similar to the papers by ^{Stephanovich & Godłowski (2015}, ²⁰¹⁷⁾ we consider here the quadrupolar (tidal) interaction of the stellar objects

where G is the gravitational constant, Q_i and m_i are, respectively, the quadrupole moment and mass of i-th object, r_ij ≡ |r _ij |, r _ij = r _j − r _i is a relative distance between objects, and θ is the apex angle. The Hamiltonian function (2) describes the interaction of quadrupoles formed both from luminous and dark matter; see ^{Stephanovich & Godłowski (2017)} for details.

To account for the mass distribution (1), we begin with the expression for the characteristic function F(ρ) of the random gravitational fields distribution (^{Stephanovich & Godłowski 2015}, ²⁰¹⁷⁾.

In the spirit of the article of ^{Chandrasekhar (1943)}, we rewrite the expression (3) in the form

where n(r ,m) is the number density (concentration, proportional to probability, see below) of stellar objects (galaxies, their clusters, and also dark matter haloes) at the position r with a mass m. As the average density at large scales can be regarded as constant (slowly spatially fluctuating to be specific, see ^{Chandrasekhar 1943}; ^{Press & Shechter 1974)}, the number density n in equation (4) can be safely considered to be spatially uniform, i.e. n = n(m). In this case the expression (4) reads

where

is the quadrupolar field (^{Stephanovich & Godłowski 2015)}, m is the mass of a stellar object (like a galaxy or cluster) and R is its mean radius. We take the function n(m) ≡ f(m) in the form of the Shechter function (1), where m _∗ and α are adjustable parameters. We obtain the normalization constant A from the condition (see equation 4 of ^{Press & Shechter 1974)}

where n is our previous constant concentration (^{Stephanovich & Godłowski 2015}, ²⁰¹⁷⁾. Note that there is no problem to take any other dependence n(m) which will not complicate our consideration. As we mentioned above, here following ^{Shechter (1976)}, we assume that the luminosity is directly proportional to the first power of the mass ^L˜ ∼ m. But there is no problem to consider the higher powers in this relation, like L~ ∼ m ^k , k = 4. In this case, the argument of the function (1) will be m ^k instead of m. The explicit calculation gives

Here Γ(z) is Euler Γ - function (^{Abramowitz & Stegun 1972)}. Finally we have from (5)

where E(r ,m) is given by equation (6). It turns out that equation (8) can be reduced to equations (17) and (18) from ^{Stephanovich & Godłowski (2015)} but with a slightly renormalized coefficient before ρ ^3/4 . This is because, under the assumption that n does not depend on the coordinates (it depends only on the mass, see equation 1), the coordinates and the mass turn out to be effectively decoupled. To proceed, we perform first the integration over d ³ r in (8). This integration is exactly the same as that in ^{Stephanovich & Godłowski (2015)} (since the mass enters equation (6) through parameter E ₀ which is unimportant for coordinate integration) so that we have from (9)

The integral in (10) can be performed as follows

Finally

where

The expressions (12), (13) give the answer for the case when we have a Shechter distribution for the galaxy masses. The difference between the previous case (^{Stephanovich & Godłowski 2015}, ²⁰¹⁷⁾ of a single mass is that now the width of the distribution function of random gravitational fields depends on the fitting parameters m _∗ and α. Note that for any function n(m) the result for the characteristic function F(ρ) will be equation (12) but with a different coefficient κ.

3. CALCULATION OF THE MASS DEPENDENCE OF THE MEAN ANGULAR MOMENTUM

To derive the mass dependence of the mean angular momentum, we should first calculate the distribution function of gravitational fields f(E), then using a linear relation between angular momentum L and field E (here, without loss of generality, we consider the moduli of corresponding vectors, see ^{Stephanovich & Godłowski 2015}, ²⁰¹⁷ for details), we derive the distribution function f(L), from which we obtain the desired dependence.

The expression for the field distribution reads (^{Stephanovich & Godłowski 2015)}

where F(ρ) is the characteristic function (12). Function f(E) is normalized as follows

and coefficient κ is given by the expression (13). The distribution function of angular momenta can be expressed in the usual way from the f(E):

This gives explicitly (see ^{Stephanovich & Godłowski 2015}, ²⁰¹⁷⁾

where

Function L ₀(t) defines the model (CDM or ΛCDM) used by ^{Stephanovich & Godłowski (2017)}.

3.2. ΛCDM Model

Here, similar to ^{Stephanovich & Godłowski (2015)}, we should isolate the contribution from the time dependent functions f1,2(τ)(τ=t/t0, t0=2/(3H0ΩΛ), see equation (46) of ^{Stephanovich & Godłowski (2015)}; i = 1,2 numbers the orders (first and second respectively) of perturbation theory. Following ^{Stephanovich & Godłowski (2015)}, we have for the argument of the distribution function H(λ,t) (equation 44).

λ(τ)=LIκ4/31fi(τ). (27)

Similar to the above CDM model, the maximum of the distribution function λ_max = 0.602730263 generates the following relation

Lmax=Iκ4/3⋅0.6027fi(τ)≡ηΛCDMi(τ)m2R4n1/3,ηΛCDMi(τ)=fi(τ)0.6027Ψ4/3G2(α+1),  Ψ=2π⋅0.418⋅Γα+74Γα+1. (28)

The relation (28) is almost similar both in CDM and ΛCDM models. The only difference is in the functions f_i (τ), where in the CDM model

f1(τ)=(2/3)τ,  f2(τ)=(-4/3)τ1/3. (29)

It is seen that the substitution of functions f_1,2 (29) yields immediately the expressions (25) and (26) for the CDM model. In the ΛCDM model the functions f_1,2 (τ) should be taken from the solution of the differential equations (30) and (31) of ^{Stephanovich & Godłowski (2015)}.

As the expression (28) for L_max is formally equivalent to equations (25) and (26), the dependencies L_max (M) are the same as those defined by equations (25) and (26) except that we should use now η_ΛCDMi (τ).

The representative plots of the function H(λ,τ) = 2I(λ,τ)/(πλ) (^{Stephanovich & Godłowski 2015)} and the maximal value λ _max(τ) for the mass-dependent density (1) are reported in Figure 1. The left and middle panels show the results for the first (equation 30 of ^{Stephanovich & Godłowski 2015}) and second orders (equation 31) of the perturbation theory. First, it is seen that the distribution of masses does not make a qualitative difference in the shape of the distribution function. That is, the shape of the dashed black curves (those without mass distribution) and the red and blue ones is the same. At the same time, the distribution functions H(λ,τ) have substantially smaller amplitudes in the case of a mass distribution. This means that the presence of a mass distribution changes the functions H(λ,τ) quantitatively. The influence of the Shechter exponent α (^{Press & Shechter 1974}; Shechter 1976) is minimal - the curves H(λ,τ) as well as λ_max (τ) (right panel) are almost the same for very different values of α = 1.25 and 0.5. This means that while the distribution of masses by itself changes the distribution function quantitatively, the value of the constant α in that distribution is of minute influence. Our analysis shows that the above tendency persists for any time instant (in Figure 1 we have only two time instants τ = 0.5 and 1) and any reasonable α > 0. The maximal values λ_max are almost independent of the presence of a mass distribution. Really, it is seen from the right panel of Figure 1 that the curves for no mass distribution (black dashed lines) and those for it lie very close to each other. This means simply that the physics of the system under consideration is determined by the distribution of random gravitational fields, rather than that of the masses of stellar objects. Maybe the “complete” dependence n(r ,m) (4) (rather than the present simplified situation n(m) in each spatial point r (5)) will improve the situation. On the other hand, it is acceptable that the robust distribution of gravitational fields is simply not susceptible to small corrections like a mass distribution. The latter, in turn, may mean that the next important step in the physics of galaxies formation is to consider the short-range interaction between galaxies (due to the presence of dark matter, for instance) so that the real average angular momentum (and not the distribution function maximum, considered so far) will appear; see equation (48) of ^{Stephanovich & Godłowski (2015)} and equation (28) of ^{Stephanovich & Godłowski (2017)}.

Fig. 1 Time evolution of the distribution function H(λ,τ) in the ΛCDM model for the case of mass-dependent density (1). Left panel: first order perturbation theory. Middle panel: second order. Right panel: the dependence λ_max (τ) on the first and second orders of perturbation theory. In all panels, dashed lines (marked as “No distr.” in the right panel) correspond to the previous case of mass-independent density n. We consider two Shechter exponents . (red curves) and 1.25 (blue curves), coded by colors and explained in the legend in the right panel. Parameter of ΛCDM model α1=(1-ΩΛΩΛ)1/3=1. The color figure can be viewed online. 

4. OBSERVATIONAL DATA

The first part of our data is the sample of rich Abell clusters containing at least 100 member galaxies each (^{Pajowska et al. 2019}). The sample contains 247 clusters and was selected on the basis of the PF catalogue (^{Panko & Flin 2006}, see ^{Pajowska et al. 2019} for details). However, in the present paper we decided to restrict ourselves to 187 clusters which have explicit redshifts. As our PF cluster sample was not numerous enough to confirm the hypothesis that galaxy alignment decreases with redshift, we decided to enlarge our sample with the DSS base.

From the ACO Catalogue (^{Abell, Corwin & Olowin 1989)} we selected all Abell clusters with galactic latitude b > 40^◦ and richness class ≥ 1. This gave us 1238 structures of galaxies from which we selected only those with redshifts z < 0.2 (^{Struble & Rood 1999)}. Therefore, 377 clusters were left for analysis. From DSS we extracted the area covering 2Mpc × 2Mpc (h = 0.75, q ₀ = 0.5) around each cluster. We applied the FOCAS package (^{Jarvis & Tyson 1981)} to the extracted regions and obtained catalogues of galaxies, considering objects within the magnitude range (m ₃, m ₃ + 3), where m ₃ is the magnitude of the third brightest galaxy. The catalogues obtained automatically were visually corrected in order to reduce the possible incorrect star/galaxy classification. FOCAS calculates the catalogue parameters using the moments of pixel distribution in an object. There are three steps from the basic image to the object list in FOCAS: segmentation, area assembly and object evaluation. The various parameters characterizing the individual images in the segmented areas were calculated. In FOCAS, the object location is defined by centroids:

where M ₀₀ is the zero moment, which is equal to:

The summation over A means that the sum includes all pixels in the object-defining area A. I(x,y) is the intensity corresponding to the density at the location (x,y) in the digital plate image. I_s is the intensity corresponding to the average plate density at the object location. Shape information about the object is obtained from the higher central moments:

The object position angle is calculated using the above central moments:

The galaxy ellipticity is

where

Each catalogue contains information about the right ascension and declination of each galaxy, its coordinates x and y on the photographic plate, instrumental magnitude, object area, galaxy ellipticity and the position angle of the major axis of the galaxy image. The equatorial galaxy coordinates for the epoch 2000 were computed according to the rectangular coordinates of the DSS scans. We calculate the position angle and ellipticity of each galaxy cluster using the method described by ^{Carter & Metcalfe (1980)} which is also based on the first five moments of the observed distribution of galaxy coordinates x_i , y_i .

5. STATISTICAL STUDIES

^{Hawley & Peebles (1975)} proposed to analyze the distribution of galaxies angular momenta by that of the observed position angles of the galactic image major axes. The direction of the angular momentum is then believed to be perpendicular to that of the major galaxy axis. This means that in the original version of the method the face-on and nearly face-on galaxies must be excluded from the analysis. This method can also be extended to studies of the spatial orientation of galaxy planes (^{Flin & Godłowski 1986)}.

The idea of ^{Hawley & Peebles (1975)} is to use the statistical tests to investigate the position angle distribution. A high value of statistics means a greater deviation from an isotropic distribution i.e. a stronger alignment of galaxy angular momenta in the analyzed structures. Since the ^{Hawley & Peebles (1975)} paper, this method has become the standard tool for searching galactic alignments. A recent improvement and revision of this method was presented by ^{Pajowska et al. (2019)}.

In the present paper we follow the analysis from ^{Stephanovich & Godłowski (2017)}. The entire range of investigated angles is divided into n bins. As the aim of the method is to detect non-random effects in the galaxy orientation, we first check if the considered distribution deviates from an isotropic one. Following ^{Stephanovich & Godłowski (2017)}, in the present paper we use χ ² and Fourier tests. We also extend our analysis with first auto-correlation and KolmogorovSmirnov (K-S) tests (^{Hawley & Peebles 1975}; ^{Flin & Godłowski 1986}; ^{Godłowski et al. 2010}; ^{Godłowski 2012}; see ^{Pajowska et al. 2019} for last review). The statistic χ ² is:

where p_k are probabilities that the chosen galaxy falls into the k-th bin, N is the total number of galaxies in a sample (in a cluster in our case), N_k is the number of galaxies within the k-th angular bin and N_0,k = Np_k is the expected number of galaxies in the k-th bin. Note that the number of degrees of freedom of the χ ² test is n−1, mean value E(χ ²) = n−1 while the variance σ ²(χ ²) = 2(n−1). As in our analysis n = 36, we obtain the values E(χ ²) = 35 while σ ²(χ ²) = 70, i.e. σ(χ ²) = 8.367.

The first auto-correlation test quantifies the correlations between galaxy numbers in neighboring angle bins. The statistic C is given by

where N_n+1 = N ₁. When, as in the present paper, we analyze the distribution of the position angles, than all_√ N_k,0 =

N p_k are equal to each other and E(C) = −1 while D(C) ≈ n (i.e. σ(C) ≈n = 6), see ^{Godłowski (2012)}; ^{Pajowska et al. (2019)} for details.

If we assume that the deviation from isotropy is a slowly varying function, we can use the Fourier test:

In this test, statistically important are the amplitudes

(only the first Fourier mode is taken into account) or

where the first and second Fourier modes are analyzed together. We investigate the statistics Δ1/σ(Δ1)=(Δ112/σ2(Δ11)+Δ212/σ2(Δ21))1/2 and Δ/σ(Δ)=(Δ112/σ2(Δ11)+Δ212/σ2(Δ21)+Δ122/σ2(Δ12)+Δ222/σ2(Δ22))1/2 (see ^{Godłowski et al. 2010}; ^{Godłowski 2012}; ^{Pajowska et al. 2019} for details). Note that the expression means that the elements of ∆² should be divided by their errors rather than the total factor ∆² _j divided by its error. The expectation values of the total factors are EΔ1σ(Δ1)=1.2247 and EΔσ(Δ)=1.8708 while σ²(∆₁ /σ(∆₁)) = 1/2, and σ²(∆/σ(∆)) = 1/2 (i.e. errors of the total factors equal to √2 - see ^{Godłowski 2012}; ^{Pajowska et al. 2019} for details). In the case of the K-S test, the statistic under study is λ:

which is given by limiting Kolmogorov distribution, where

and F(x) and S(x) are the theoretical and observational distributions of the investigated angle, respectively. ^{Wang et al. (2003)} analyzing the limiting form of the D_n function, found that µ(λ) = 0.868731 while σ ²(λ) = 0.067773 i.e σ(λ) = 0.260333 (see also ^{Pajowska et al. 2019} for discussion).

Using the extended ^{Hawley & Peebles (1975)} method it is possible to analyze both the alignment dependence on a particular parameter, like richness of a galaxy cluster ^{Godłowski et al. (2010)}, and to quantitatively answer the question if an alignment is present in a sample (see ^{Pajowska et al. 2019} for last revision). In our previous papers (^{Godłowski et al. 2010}; ^{Stephanovich & Godłowski 2017)}, using a sample of 247 PF rich Abell clusters, it was shown that the alignment of galaxies in a cluster increases significantly with its richness. Unfortunately, the available data were insufficient for a persuasive conclusion about the correctness of the theoretically predicted dependence of the analyzed statistics on redshift z (^{Stephanovich & Godłowski 2015}, ²⁰¹⁷).

For this reason we performed the investigations of our samples of galaxy clusters checking if there is a significant dependence between the analyzed statistics and both richness and redshift of the clusters. As a first step we analyzed the linear model Y = aX + b. The Y are the values of the statistics i.e. χ ², ∆₁ /σ(∆₁), ∆/σ(∆), C and λ (see ^{Stephanovich & Godłowski 2017} for details) while X is the number of analyzed galaxies in each particular cluster or its redshift z respectively. Our null hypothesis H ₀ is that the statistic Y does not depend on X . This means that we should analyze the statistic t = a/σ(a), which is Student’s distribution with u−2 degrees of freedom, where u is the number of analyzed clusters. In other words, we test the H ₀ hypothesis that t < 0 against the H ₁ hypothesis that t > 0, where t > 0 corresponds to the case of dependence on the number of member galaxies in clusters and t < 0 to the case of dependence on redshift z. In order to reject the H ₀ hypothesis, the value of observed statistic t should be greater than t_cr which could be obtained from the tables. For example, for sample of 247 clusters analyzed in ^{Stephanovich & Godłowski (2017)} (our sample A) at the significance level α = 0.05, the value t_cr = 1.651. For our sample B (564 clusters) at the significance level α = 0.05, the value t_cr = 1.648.

The results of our analysis are presented in Table 1 and Figures 2, 3. One could discern from the Table 1 that if the analysis of sample A confirms that the alignment increases with cluster richness, then any test confirms the negative deviation of the linear regression parameter from zero in the case of the alignment depending on redshift. However, both the above dependencies can be confirmed from the analysis of the full sample (564 clusters). At first we analyzed the statistic i.e. alignment, which increases significantly with the richness of a cluster, confirming the result obtained in ^{Godłowski et al. (2010)}; ^{Stephanovich & Godłowski (2017)} as well as the theoretical prediction (^{Stephanovich & Godłowski 2015}). The details are also presented in Figures 2(a) - 2(e). Moreover, we could conclude that the alignment decreases with z, which means that it increases with time as predicted by ^{Stephanovich & Godłowski (2015}, ²⁰¹⁷⁾; ^{Schmitz et al. (2018)} (see Figures 3a - 3e).

Fig. 2 The dependence of the number N of galaxies in a cluster on different statistical parameters of the Sample B. The error boundaries, at a confidence level 95%, are presented as well. 

Fig. 3 The dependence of cluster redshift z on different statistical parameters of the Sample B. The error boundaries at a confidence level 95%, are presented as well. 

However, a closer look at the results shows that situation is not yet clear. This is because in the real data the cluster richness usually decreases also with redshift z. Quantitatively, in a linear model, the dependence between richness and redshift of a cluster is N(z) = az +b. In this model, we obtain the value of the t statistic t = −7.066. This is the reason why we repeated our analysis as a 3D model Y = a ₁ N + a ₂ z + b. Note that until now, such a 3D analysis has not been performed in galaxy alignment studies, but due to the above reason, we consider it to be necessary here. In this extended analysis, the test statistic t is given by formulae t ₁ = a ₁ /σ(a ₁) and t ₂ = a ₂ /σ(a ₂). From Table 2, we could confirm that the alignment increases significantly with cluster richness. One should note that all our tests show that the statistic is decreasing with z. However, the values of t are too small to be statistically significant (significance level α = 0.05).

6. SIMULATIONS

The ^{Illustris Project (2018)} was the simulation base for our present study. The Project uses the AREPO code for the hydrodynamic realizations of a (106.5 Mpc)³ cosmological volume (^{Springel 2010}). The simulation assumes a ΛCDM cosmology with Ω _m = 0.2726, Ω_Λ = 0.7274, Ω _b = 0.0456, σ ₈ = 0.809, n_s = 0.963, and H ₀ = 100·h·km·s⁻¹Mpc⁻¹ with h = 0.704. It contains multiple resolution runs with the highest resolution performed for Illustris - 1. Three different physical configurations have been applied: dark matter only, as well as non-radiative and full galaxy formation. In the first case (dark matter only), the mass was treated as collisionless in the simulations. The non-radiation configuration also adds gas hydrodynamics, but ignores radiative cooling and star formation processes. The full galaxy formation physics contains (in addition to the previously mentioned ones) also processes related to galaxy emergence, through a model described in ^{Vogelsberger et al. (2013)}. Illustris-1 consists of 136 runs for different redshifts z, where the initial conditions were generated at z = 127 for snapshot 0 and evolved to z = 0 for snapshot 135.

Illustris successfully follows the coevolution of dark and visible matter. Haloes, subhaloes, and their basic properties have been identified with the FOF and SUBFIND algorithms (^{Davis et al. 1985}; ^{Springel et al. 2001}; ^{Dolag et al. 2009)}, at every one of the 136 stored snapshots. We have added information from the supplementary catalog to the resulting directory of Haloes from ^{Zjupa & Springel (2017)}. The code was written in such a manner that it can run both as a postprocessing option to increase existing catalogues, or as part of the regular group finding.

From Illustris-1 we selected haloes at z = 0. We obtained 119 haloes with total mass 1013M⊙ and 1435 with total mass larger than 1012M⊙. The angular momentum parameter for the extracted haloes was taken from ^{Zjupa & Springel (2017)}.

Ilustris simulations give direct values of both the mass of the structures and their angular momentum. Present available data from Ilustris are evolved to z = 0, so it is possible to study the dependence of angular momentum as a function of cluster mass, but unfortunately not of redshift. However, as we know directly the cluster angular momentum, it is not necessary to assume a linear relation between angular momentum and mass. Since theoretical modeling predicts usually power law relations (see also ^{Stephanovich & Godłowski 2015} for review) we could study the model J = b · M ^a . The latter relation could easily be rendered as a linear model: ln J = ln b + a ln M.

The results of the analysis are presented in Table 3 and Figure 4. The analysis of Ilustris simulations confirms that the angular momentum of a cluster increases with its mass. In this case, the coefficient a = 1.807 ± 0.028 (Figure 4a). The analysis of a sample with only the most massive clusters (mass M > 10¹³ solar mass) left, gives a = 1.708 ± 0.114 (Figure 4b) which is close to the most popular theoretical prediction a = 5/3 ≈ 1.667 (see ^{Godłowski et al. 2010}; ^{Stephanovich & Godłowski 2015}, ²⁰¹⁷ for details).

Fig. 4 Relation between angular momentum and mass of a cluster derived from the Ilustris simulation for different total masses M. The error boundaries at a confidence level 95%, are presented as well. 

7. OUTLOOK

In the present paper we have shown that the distribution of mass of the stellar objects does not alter substantially the distribution function of their gravitational fields. This shows that the main contribution to the latter distribution function comes from the long-range Newtonian interaction between astronomical objects, rather than from the distribution of their masses. At the same time, the mass distribution alters the dependence of L_max on the concentration of astronomical objects n (see equation 25) and on the total cluster mass M (see equation 26), which is observationally important. To discern which of the dependencies (25) or (26) is realized in practice, additional observational data are needed. As the total interaction potential contains both luminous and dark matter components, one can ask a question about the alignment of sub-dominant galaxies, even though the majority of the angular momenta of galaxy clusters is related to the smooth dark matter halo component. This question becomes important in view of the fact that the mass distribution (1) alters the dependence on total cluster mass M (26). In the halo model (^{Schneider & Bridle 2010)} where the galaxies are embedded in a dark matter halo, the latter may mediate the intergalactic interaction, adding possible short-range terms to it. That is to say, to “see each other” in a dark matter halo, the galaxies should be closer than in an empty space.

Note that the observational results about a lack of alignment of galaxies for less clumpy (so called poor) clusters, as well as evidence for such an alignment in the clumpy (rich) ones (^{Godłowski, Szydłowski, & Flin 2005}; ^{Aryal et al. 2007}; see also ^{Godłowski 2011} for incremental study and relevant references) clearly show that the angular momentum of galaxy groups and clusters increases with their mass (richness). The generalized analysis, based on equation (4), where n = n(r,m) (i.e. the mass becomes spatially distributed), will improve the overall understanding, which can additionally be tested against observed galaxy shape distributions and alignments. The problem of angular momenta alignment due to their interactions, as well as those with dark matter haloes, has been simulated by ^{Hahn et al. (2007)}. The main effect there is the presence of a threshold cluster mass (richness) value. The latter is related to the mutual alignment of clusters and dark matter haloes axes. This fact can be analyzed on the base of a more general model (4), which accounts for a spatially inhomogeneous distribution of the number density of stellar objects, as well as for its mass dependence. We postpone the consideration of this interesting question for future publications.

Our formalism permits studying this effect (see ^{Stephanovich & Godłowski 2017)} as well as the nonequilibrium time evolution of luminous astronomical objects (with respect to dark matter haloes) within the ΛCDM model. The combination of stochastic dynamical approaches (^{Garbaczewski & Stephanovich 2009}, ²⁰¹¹⁾ along with a deterministic one defined by the ΛCDM model, may permit to answer (at least qualitatively) the question about the initial alignment of galaxies and their clusters) at the time of their formation. The questions about how dark matter haloes influence (mediate) the latter alignment can also be answered within the above dynamical approach.

Our statistical analysis of the Abell cluster sample shows that the alignment of the angular momenta of galaxies and their clusters increases substantially with cluster mass. This result is confirmed also by 3D analysis, consisting in the study of the dependence of galaxy alignment in a cluster both on its richness and redshift. We have also found that the alignment decreases with redshift i.e. increases with time, but 3D studies show that this effect is too faint to be confirmed statistically at the significance level α = 0.05. A probable reason is that the corresponding relaxation time is too long. So, for future investigations more extensive data containing larger numbers of galaxy clusters (and with higher redshifts) are required. The comparison of our theoretical results with those of the Ilustris simulation also confirms the increase of the angular momenta of galaxies with cluster mass. Moreover, the latter comparison confirms the power law relation with a coefficient very close to 5/3, which is the value favored by most popular theoretical predictions; see, for instance, ^{Stephanovich & Godłowski (2015}, ²⁰¹⁷⁾. This suggests that the approach of ^{Stephanovich & Godłowski (2015}, ²⁰¹⁷⁾ correctly reflects the main features of the formation of galaxies and their clusters.

Let us finally note that, as we have shown above, our results are in close conformity with the commonly preferred model of galaxy formation, i.e. the so-called hierarchic clustering model (^{Peebles 1969}), improved recently by taking into account a tidal torque scenario.