G Pascoli
email: pascoli@u-picardie.fr

A comparative study of MOND and MOG theories versus the κ-model: An application to galaxy clusters

Keywords: dark matter, MOND, modified gravity, κ-model, galaxies, galaxy clusters

. A new model, called κ-model, based on very elementary phenomenological considerations, has recently been proposed in the astrophysics field. This model shows that the presence of dark matter can be considerably minimized with regard to the dynamics of galaxies (Pascoli, 2022 a,b). The κ-model belongs to the general family of theories descended from MOND. Under this family of theories, there is no need to develop a highly uncertain dark matter sector of physics to explain the observations.

Introduction

The dark matter (DM) paradigm is considered to simply explain the dynamics in individual galaxies and galaxy clusters. Astrophysicists have been searching for DM evidence for years. However, the quantity of dark matter required is immense, and the ratio of the dark matter (DM) to the baryonic component (B), (DM/B), could largely exceed 10 in a few galaxy clusters, raising serious doubts about the validity of this hypothesis. A rapid explanation to remedy this problem is to say that a large quantity of invisible baryonic matter is not counted in the galaxies and galaxy clusters. Unfortunately, this immediate solution is may be acceptable up to a factor of 2, but fully irrealistic up to a factor of 10. What is then the nature of dark matter ? To date, no particle of dark matter has never been detected in lab (DAMA/LIBRA Collaboration, 2022), even though researchers have built larger and more sensitive detectors (Xenon Collaboration, 2023). This is rather an intriguing status and the 1 DM paradigm sounds suspiciously similar to the phlogistic and the aether in the nineteenth century. Additionally, some agnostic physicists believe that DM does not exist at all and have instead proposed alternative models of gravity. As a result, all sorts of theories have been built in order to remove dark matter in all astrophysical systems. There are two theories that stand out because they have been applied to numerous concrete situations. The first is modified Newtonian dynamics (MOND). MOND diverges from the standard Newton's laws at extremely low accelerations, which are a characteristic of the outer regions of galaxies (Milgrom, 1983(Milgrom, , 2015)). MOND postulates a modification to Newton's second law, such that the force applied on a particle is no longer proportional to the acceleration a, but rather to its square a 2 , when the accelerations are smaller than the critical limit a 0 = 10 -10 ms -2 (Milgrom, 1983). This model effectively explains the dynamics of individual galaxies without dark matter (Famey and Mc Gaugh, 2012). In modified gravity (MOG) (STVG or scalar-tensor-vector gravity), the approach is very different. The structure of space-time is described by the usual metric tensor g µν , complemented by a vector field, ϕ µ , and two scalar fields, G and µ, which represent a dynamical version of the Newtonian gravitational constant and the mass of the vector field, respectively (Moffat, 2006;2020). The vector field part produces a Yukawa-like modification of the gravitational force due to a point source. This model accurately explains not only the dynamics of galaxies but also the dynamics in galaxy clusters without dark matter (Moffat, 2020). The common and main objective of MOND and MOG theories is to eliminate in whole or in part the dark matter in the Universe. The agreement between these two models and the observational data are very remarkable for galactic dynamics; however, the situation is distinct for galaxy clusters where MOG appears to have an advantage over its competitor MOND (Brownstein and Moffat, 2006). MOG achieves to eliminate all dark matter content in these objects, but MOND still needs to consider some form of invisible matter in galaxy clusters (Mc Gaugh, 2015). However, this does not mean that MOND has been falsified, because the MONDian world is very rich and there are numerous extensions, such as the promising extended MOND (EMOND) (Hodson and Zhao, 2017a, b) 1 . Moreover, we can speculate whether MOG does not reintroduce DM in a disguised manner (the vector field ϕ µ is massive). Following this simple statement, only MOND is truly free of DM, and the path followed by EMOND is potentially preferable. The third model studied here is the κ-model; it is based on a phenomenological and MON-Dian procedure whose main aim is to simultaneously explain the dynamics of individual galaxies by minimizing the dark matter content and by maintaining the formal aspect of the Newtonian law of gravitation. It is based on a relational consideration; the mean density of matter is estimated at a very large scale, and the surroundings of a given observer influence the measurements made for the determination of the velocities and accelerations. In this regard, the observations depend not only on the reference frame of the observer but also on his environment. Thus, the κ-model uses a holistic or Machian approach. The velocities and accelerations that are environment-dependent are renormalized, and the observer measures apparent quantities depending on a coefficient denoted κ. An empirical (and universal) relationship is provided between the coefficient κ and the mean density. The coefficient κ intervenes in the acceleration term of the dynamics equation and imparts a MOND-type appearance to the κ-model. However, a physical support is now provided; this is the environment of the observer, which distorts the measurements. Thus, a naive analogy is that of an observer placed in a medium of given refraction index, who sees a magnification of both the size and the velocity of any object. However, this very simplistic comparison is not to be taken at face value because in a medium of given refraction index, the light is attenuated and the images can be extremely blurred. On the contrary, in the κ-model framework, no such medium is existing; the light is not attenuated and it still propagates in straight line with the speed c, constant and independent of frequency, as measured in vacuum by every observer (Pascoli, 2022a) 2 . The great interest of the κ-model is that no arbitrary parameter is introduced, contrary to MOG, which has two outer free parameters (these two models having a universal relationship for the fits). Another advantage of the κ-model is that it would permit to pass directly for any type of galaxies from the data on the spectroscopic velocity measurements to the mean densities (and inversely) without the ambiguous transition through the mass-to-light ratio knowledge. The relationship between the spectroscopic velocities and the mean densities is direct in the κ-model. Strongly contrasting with this view, DM can fit almost all rotational curves with any mean baryonic density profiles due to its very large flexibility, and for this reason the DM paradigm has unfortunately no predictive value. Eventually the κ-model can easily help to understand the weaker "DM effect" in the regions where the mean density is high (globular clusters, core of the galaxies) and, conversely, the stronger "DM effect" in the regions where the mean density is weaker (outer regions of the galaxies, low brightness galaxies and galaxy clusters). More specifically, the κ-model predicts that when the mean density ρ in a large-scale object (galaxy or galaxy cluster) is smaller than a critical value ∼ 4 10 -24 g cm -33 the measured velocities appear magnified compared to the estimated newtonian velocities. The κ-model has been already applied to various types of galaxies, such as large or small low surface brightness galaxies (LSBs) and high surface brightness galaxies (HSBs) (Pascoli, 2022 a,b). The κ-model fits and the observational curves has been shown in fairly good agreement. Even though this study needs to be extended to a larger database, such as SPARC (Lelli, McGaugh and Schombert;2017), these initial results and deductions provided by the κ-model are highly encouraging. While the κ-model appears to succeed at explaining the dynamics of numerous individual galaxies (Pascoli, 2022), to date, this model has not been tested on the galaxy cluster field, and this is the goal of our present study. We show that the κ-model can greatly lower the major part of the dark matter content in the galaxy clusters. The current ratio DM/B amounts around 10 in the outskirts of these objects, and the κ-model can reduce this ratio to approximately 0 -1. The situation is more difficult in the inner regions of the galaxy clusters, but by lowering the gas temperature in these regions, the problem can be easily solved. The κ-model provides a good predication for the ratio DM/B in the outer regions of galaxy clusters, without introduction of numerous free parameters other than those linked to the mean density of baryonic matter (Pascoli, 2022 a,b). Likely, this result is not a mere coincidence, and supports the consideration of our new proposal. Here, the κ-model is applied to a sample of galaxy clusters (paragraph 2) and then, eventually, to the Bullet Cluster (paragraph 3).

Galaxy clusters

The gas density profile for a galactic cluster can be approximately fitted by the following function (Cavaliere, and Fusco-Femiano, 1976):

ρ (r) = ρ M 1 + r r c 2 -3β 2 (1)
where ρ (r) is the intracluster medium (ICM) mass density profile and ρ M is the maximal value taken by ρ (r). r c and β are fit parameters for the distribution of the density. Due to this isotropic density distribution, the gas mass contained in a sphere of radius r is as follows:

M gas (r) = 4π r 0 drr ′ 2 ρ (r ′ ) = 2 F 1 (1.5, 1.5β, 2.5, -0.000017 r 2 ) (2)
where 2 F 1 is the hypergeometric function.

When r ≫ r c and β < 1, we can approximate this formula by the more manipulable relationship, as follows:

M gas (r) = 4πρ M r 3 c 3(1 -β) ( r r c ) 3(1-β) (3)
Unfortunately this relationship is clearly divergent when r -→ ∞. A cut-off for the distribution of gas, necessary for finite spatial extent, needs to be introduced. Following Brownstein and Moffat (2006), this cut-off, let rout, is chosen equal to the radius at which the density, drops to 10 -28 g cm -3 or 250 times the mean cosmological density of the baryons. On the other hand, assuming that the cluster is in hydrostatic (isothermal) equilibrium, the Newtonian dynamical mass is as follows (Brownstein and Moffat, 2006, eq. 19):

M N (r) = 3β k B T µm p G r 3 r 2 + r 2 c (4)
where T is the temperature, k B is Boltzmann's constant, µ ≈ 0.609 is the mean atomic weight and m p is the proton mass.

In eq.4 the quantities β, k B , T , µ, m p and G are κ-invariant, both r and r 2 c4 are transformed as r -→ κ κ E r. Then the κ-mass profile for a cluster is given by the following relationship:

M κ (r) = κ κ E M N (r) (5)
where M N (r) is the Newtonian mass evaluated at the radius r. Then the mean density ρ (r) is inserted in the relationship (κ, ρ) (see the appendix A eq. 18 placed at the end of the article) and this leads to the magnification ratio, as follows:

κ E κ = 1 + Ln ρ E ρ = 1 + Ln ρ E ρ M ρ M ρ = 1 + Ln ρ E ρ M + 3β 2 Ln 1 + r r c 2 (6)
with the mean mass density ρ E near the Sun estimated to 4 10 -24 g cm -3 . We used the sample of galaxy clusters from the paper of Brownstein and Moffat (2006). The relevant properties are listed in Table 1. By convention column (6) is the position, r out , at which the density drops to ρ out ≃ 10 -28 g cm -3 , or 250 times the mean cosmological density of the baryons.

Figure 1 shows the results of our analysis applied to the COMA cluster. The κ-curve (shown in amber) has approximately the same profile as MOND (shown in dashed-dotted cyan) when r > 100 kpc, but the κ-curve is much closer to the observational curve and even merges with the latter one in the outer regions of the cluster. Even though the κ-curve is not fully merged with the ICM gas curve inside the inner regions; however, the apparent gravitational mass has been largely lowered (by a factor in the range 7 -10 along the curve). Finally a residual gap has to be filled in between the κ-curve and the observational one in the inner regions of the COMA cluster. To do this, we propose to decrease the temperature in the inner regions. This proposal can be applied in the same way to the other galaxy clusters (Figure 2). Notably, in most of the cases presented on this array of figures, MOG theory leads exactly to the same prediction as the temperature profiles after comparing the curves (we can refer to A0085, A0133, NGC 507 and A062 as examples). Assuming a non-isothermal temperature profile, the Newtonian mass has to be re-calculated as follows (Brownstein and Moffat, 2006, eq. 18):

M N (r) = r k B T (r) µm p G 3βr 2 r 2 + r 2 c - dLnT (r) dLn(r) (7)
Then we used the following easy-to-manipulate temperature profile:

T (r) = T out exp[-α( r out -r r out )] (8)
where T out designates the temperature in the outer regions of the cluster. This parameter is provided in column (2) of Table 1. The coefficient α = 1 is used when M N ≲ 10 M gas . This is the case for most situations under study. For a few cases where this ratio is far beyond 10, then α = 2 is used. The temperature profile (eq.8) has been selected such that the outer temperature, T out (column (7) of Table 1) is that provided by Brownstein and Moffat (2006), as directly resulting from observations [START_REF] Reiprich | Cosmological Implications and Physical Properties of an X-Ray Flux-Limited Sample of Galaxy Clusters[END_REF].

M N (r) = r k B T (r) µm p G 3βr 2 r 2 + r 2 c - αr r out (9)
The green curve on the Figure 1 is the κ-profile with a non-isothermal temperature profile given by eq. 8. The κ-curve associated to a non-isothermal temperature profile becomes virtually superimposed to the ICM gas profile when r < 500 kpc (let us note that for r > 500 kpc, the κ-curve, associated to an isothermal temperature profile and that associated to a non-isothermal temperature profile flank the ICM gas profile). Thus, a prediction of the κ-model is that the temperature is lowered by a factor two (at r out /2) in the inner regions compared to that in the outskirts of the COMA cluster. More precisely, the most appropriate description is that of a non-isothermal core for r < 1000 kpc, which is itself surrounded by a quasi-isothermal shell when 1000 < r < 1954 kpc. A vey similar situation is encountered in the other clusters collected in Table 1. Thus, since the κ-model has no outer parameters, an optimal fit can be achieved by finely adjusting the temperature profile and, in return, can eventually predict this profile (whereas this refinement and prediction are not possible within the ad hoc, and definitely not falsifiable, DM paradigm, which is easily adapted to fit the ICM gas curve with any temperature profile indicating no predictive value). Figure 1 COMA cluster profile. The horizontal axis is the radius in kpc and the vertical axis is mass in units of the solar mass M ⊙ . The red long dashed curve is the ICM gas mass derived from X-ray observations (compilation of [START_REF] Reiprich | Cosmological Implications and Physical Properties of an X-Ray Flux-Limited Sample of Galaxy Clusters[END_REF][START_REF] Reiprich | Cosmological Implications and Physical Properties of an X-Ray Flux-Limited Sample of Galaxy Clusters[END_REF]; the short dashed blue curve is the Newtonian dynamic mass; the dashed-dotted cyan curve is the MOND dynamic mass; the solid black curve is the MSTG dynamic mass (Browstein and Moffat, 2006). Our contribution is displayed as the amber curve, showing the κ-model dynamic mass with the temperature T = 8.38 keV . The solid green curve displays the κ-model dynamic mass, assuming a non-isothermal temperature profile with α = 1 in eq. 8.

A perfect superposition with the ICM gas mass curve in COMA cluster can still be achieved by taking the temperature profile displayed in Figure 2. We note an increase of the mean temperature from the inner regions up to 1500 kpc and then a slow decrease toward the outskirts. However, a comparison of this profile with observational data is not very conclusive.

The main reason is that the COMA cluster modeled here, similar to other research (Browstein and Moffat, 2006), is in the form of a spherical and well-relaxed distribution of gas. The reality is much more complex. The density in the inner regions is relatively high, and the cooling through thermal bremsstrahlung emission must effectively be much more efficient than in the outskirts. However, the heating by active galactic nuclei in the inner regions can compensate for the cooling processes. Thus, the energetics are very complex and the temperature profile is not easily predicted (Bykov et al, 2015). In a general manner for the galaxy clusters, all temperature profiles found in the literature are model-dependent with a very large amount of dark matter. Moreover, even limiting to the outskirts, which are directly accessible, the physics is poorly understood (Walker et al, 2019). Thus, the temperature map of the outskirts of the COMA cluster simultaneously exhibits cool and hot regions with various substructures (Watanabe et al, 1999).

Figure 2 Predicted mean temperature profile in the COMA cluster in the framework of the κ-model We have suggested here to solve the problem of the very large excess of (newtonian) mass in the galaxy clusters by using a two-stage procedure. First, by significantly reducing the apparent attractive mass with the κ-model. With help of this operation the ratio M DM /M gas is reduced to rather more credible values, i.e. ∼ 0 in the outer regions and to ∼ 2 in the inner regions. Then the second step consists in adapting the temperature profile with cancellation of the residual excess of mass. It seems that MOG was facing the same situation, even though Brownstein and Moffat (2006) have not tried to perform this second step (for instance for A0085, ...). Another issue for this secong step would still be to follow the MOND galaxy cluster analysis of Banik and Zhao (2022), that is to say to add a dense core composed of sterile neutrinos (see Giunti and Lasserre (2019) for a review on these hypothetical particles). The observational data not being perfect it is very likely that a two-stage solution is also needed in the framework of other theories (for instance MOND) which have been proposed to eliminate the dark mater (taken into account of the measurement uncertainties of inclination, thickness, mass-to-light ratio for the individual spiral galaxies; density/temperature profiles and clumpiness in the case of galaxy clusters, and so on). The only model which runs directly in one step is dark matter given its undue flexibility (compared to the κ-model with no flexibility).

Figure 3 Plot of the radial mass profile for clusters of the sample in Table 1. For the caption see the Figure 1. The solid red curve displays the κ-model dynamic mass, with a non-isothermal temperature profile of α = 2.

Figure 3 Continued galaxy cluster mass profiles Table 1 Galaxy cluster properties

Note -This compilation is issued from Brownstein and Moffat (2006). We have added a column for the mass Mκ.

Column (1) Galaxy cluster name Column ( 6) radius where gas ≃ 10 -28 g cm -3

Column (2) X-ray temperature Column (7) ICM gas mass integrated to r out (3) 1, M κ is smaller than M M ST G when the temperatures are higher than 5 -6 keV , whereas the reverse is true when the temperatures are lower than 5 -6 keV . On the other hand, for T ∼ 5 -6 keV , then M κ ∼ M M ST G . Eventually, when the temperatures are smaller than 1 keV , then M κ are larger than M M ST G . By contrast, convergent MOND dynamic mass, M M ON D , are systematically located too high.

β (4) rc kpc (5) rout kpc (6) Mgas 10 14 M ⊙ (7) M N (8) M M ST G (9) M M ON D (
Even though MOND substantially decreases the dark matter content for convergent MOND dynamic mass, this theory is not able to completely remove dark matter in the galaxy clusters. By comparing the columns 10 (M M ON D ) and 7 (M gas ) of Table 1, M M ON D is sometimes superior by a mean factor of the order of 2 to M gas . Galaxy cluster stability cannot be fully explained by the current MOND formulation alone. MOND, at least its initial form (Milgrom, 1983), still needs to use a residual content of invisible matter for galaxy clusters. However, owing to the great success of MOND for the individual galaxies, a proper direction is to build a multiscale MOND, which is able to simultaneously explain both individual galaxies and galaxy clusters. Fortunately, this model has been built. This is the main motivation of EMOND; EMOND assumes that there is an increase in the fundamental parameter a 0 of the MOND paradigm in galaxy clusters compared to the value selected for this coefficient in the case of individual galaxies (Zhao and Famaey, 2012;Hodson and Zhao, 2017 a, b). This rescaling adapts the parameter a 0 to the size of the object under consideration and appears to be quite natural. In this case, MOND can eventually and adequately fit the ICM gas mass integrated to r out . Note, in the context of the κ-model a physical interpretation of this rescaling in MOND is supplied. The mean density in a galaxy cluster is much smaller than the corresponding one in an individual galaxy by a factor ∼ 40 (center of the cluster) and ∼ 200 -1000 (in the outskirts of the cluster). In our study, the rescaling is given by an estimate of κ, and κ G /κ cluster ∼ 5 -10.

Despite these promising considerations for the κ-model, MOG appears to be slightly better as for the fitting of the observational curves. We might argue that MOG has two free outer parameters, whereas the κ-model has none; therefore, it is easier to fit a curve with a sufficient number of outer free parameters than without any free parameter. However, on closer examination of MOG, we can see that a multifit procedure is used throughout the calculations. MOG relies on two parameters M 0 and r 0 , but these parameters are not constant across the series of galaxy clusters. Moreover, M 0 depends on the ICM gas mass for each individual clusters (Brownstein and Moffat, 2006, paragraph 4). Thus it seems that in MOG the results are already included at the start of the calculations, and the method appears somewhat post hoc. Thus, MOG is not a fully ab initio procedure. Due to these reasons, the κ-model and EMOND are more efficient than MOG as for their predictive power.

Bullet Cluster

The Bullet Cluster (1E 0657-56) is very often presented as a clear proof of the existence of dark matter (Figure 3). The Bullet Cluster is composed of two colliding clusters: a main cluster (Mc) and a small or sub-cluster (Sc) (the bullet per se) (Bradač el al, 2005(Bradač el al, , 2006)).

Figure 4 To determine how the κ-model reinterprets the observational data for the Bullet Cluster, we initially start with the known (but apparent) surface densities for both the hot gas and the visible galaxies (Figure 5). (2006). In addition, Table 2 provides the masses for the different components. The assumed DM content is adjusted such that the total mass ratio DM/B is ∼ 6.

Figure 5 Figure 5 Abscissa : distances, unit 500 kpc. Ordinate: surface densities, unit 3.1 10 3 M ⊙ pc -2 . The abscissa axis passes through the galaxy cluster centers. Figure 5a: scaled plot of hot gas density, Figure 5b: scaled plot of apparent galaxy density observed from Earth (Main cluster + Bullet) and Figure 5c: scaled plot of real galaxy density measured by a hypothetical observer located inside the Main or the Sub group of visible galaxies (represented by blue disks in Figure 4). The mass is invariable by the transform: apparent ←→ real surface densities, i.e. dxdy Σ s (x, y) = dxdy Σ ′ s (x, y). In the usual treatment of the galaxy clusters (paragraph 2), only the gaseous component are considered (the stellar fraction is negligible in galaxy clusters). In the specific case of the Bullet Cluster, the situation is very different. We have two clearly identified parts: a very massive gaseous component and a galactic component; however, these two components are separated. Based on the lensing diagram, only the gaseous component needs to be considered in the calculation of κ. Thus, if we assume that the Main and Sub groups of visible galaxies (the disks displayed in blue in Figure 6) have a low content of gas, the lensing due to the κ effect is much higher than that in the hot gas (the disk displayed in red in Figure 6). The usual relationship of the κ-model is as follows:

κ M κ = 1 + Ln Σ M Σ gas ( 10
)
where Σ M ∼ 0.075 × 3.1 10 3 M ⊙ pc -2 = 232.5 M ⊙ pc -2 is the maximal value taken by the surface density of hot gas (Figure 5a). Figure 5a can be used to provide the following relationship:

(Σ M ain gas ) outer , (Σ Sub gas ) outer ∼ 0.7 × 232.5

M ⊙ pc -2 = 162.7 M ⊙ pc -2 (11)
The latter quantity represents the surface density of the hot gas located in the immediate environment of the Main and Sub groups of visible galaxies; each of them being displayed by a blue disk in Figure 4. However, following the current interpretation, both the Main and Sub galaxy clusters were stripped from the hot gas during the collision. The amount of the gas remaining inside the groups of visible galaxies after this collision needs to be determined. Clearly the hot gas initially located inside the Main and Sub groups of visible galaxies is removed by the strong shock during the collision, and each of these groups is now located in a subdense bubble with a low content of hot gas. A pre-analysis of the α-diagram in the κ-model context, compared to that provided by observations (Bradač el al, 2005(Bradač el al, , 2006)), enables to predict that the amount of hot gas inside the subdense bubbles is ∼ 15% of the amount of hot gas surrounding them, i.e. (Σ M ain gas ) inner or (Σ Sub gas ) inner = 24.4 M ⊙ pc -2 . With equal temperatures, the gas pressure in the groups of visible galaxies (blue regions in the Figure 4) is thus predicted to be lower than that in the immediate region surrounding them (red region). Evidently, the system is not in hydrostatic equilibrium, and this situation cannot indefinitely persist. The required time to re-equalize the densities is approximately Rc 3cs ∼ 100 M years, where R c ∼ 0.5 M pc is the mean characteristic size of the Main or Sub groups of visible galaxies and c s is the speed of sound in the hot gas (T e ∼ 10 8 K)5 . With the relationship (10) and the aforementioned values of (Σ M ain gas ) inner or (Σ Sub gas ) inner , we can determine the amplification factor κ inside the two groups of visible galaxies as follows:

κ M κ M ain gas = κ M κ Sub gas ∼ κ M κ galactic gas = 3 (12)
Figure 6 Basic illustration of the combined action of both gravitational and kappa lensings.

The lensing by the hot gas is used as a reference for the two figures a and b.

In Figure 6a, the combined product of both the gravitational and kappa lensing effects is clearly centered on the Main and Sub groups of visible galaxies. In Figure 6b, as expected, the lensing is clearly centered on the hot gas (which contains 90% of the baryonic mass against 10% for the visible galaxies in a galaxy cluster). However, the reality shown in Figure 6b cannot be viewed by a unique observer (for instance, a terrestrial observer surrounded by its own environment). Specifically, two distinct observers are needed to determine the effect. The part located in the blue area (groups of visible galaxies) is perceived by a hypothetical observer situated inside any region where the mean gas density is the same as that in the Main or Sub group of visible galaxies; the part located in the red area (hot gas) is perceived by a hypothetical observer situated inside any region where the mean density is the same as that in the hot gas. If the κ-model paradigm is expected to be on the right track, then the Bullet Cluster is an illustrative (but rare) example that the perception of the objects in the Universe is observer-dependent. Beyond these qualitative considerations, calculations are evidently needed to ascertain our purpose. The procedure is well known. The field equations of general relativity can be linearized if the gravitational field is weak. Then, the deflection angle of a set of masses is simply the vectorial sum of the deflections due to individual lenses. The plane of the sky is (x, y). The deviation angle α can be written using the thin lens approximation as follows(Bartelmann and Schneider, 2001):

α (x, y) = S dx ′ dy ′ 4G c Σ(x ′ , y ′ ) r -r ′ |r -r ′ | 2 (13)
where G is the gravitational constant and c is the speed of light. The density distribution Σ is integrated over all the surface S of the cluster system. The total surface density Σ (x, y) is the sum of the gas and stellar surface densities.

Σ (x, y) = Σ g (x, y) + Σ ′ s (x, y) (14)
where Σ g (x, y) and Σ ′ s (x, y) are the fits of the distributions represented in Figures 4a and4c, respectively; these are assumed to be approximately circular and Gaussian. We know that the κ-model does not change the local physics, apart from a magnification factor locally applied. Thus, in the κ-model, the relationship (14) eventually becomes the following:

6 α (x, y) = S dx ′ dy ′ 4G c Σ g (x, y) + κ M κ galactic gas Σ ′ s (x, y) r -r ′ |r -r ′ | 2 (15)
Figure 7 shows the lensing diagrams |α| of both dark matter and κ-model and Figure 8 shows the components α x and α y . Even though the physical interpretation is very different, the κ-model and DM diagrams are very similar and both show that the lensing is centered on the visible galaxies and not on the hot gas, which is centered at (0, 0). The essential difference between the two models is that the κ-model does not need dark matter to explain the same facts. The method applied to the Bullet Cluster can be extended to other similar cases such as the Train Wreck Cluster (Abell 620). In the latter situation the lensing is centered on the hot gas. Even though the morphology of the Train Wreck Cluster is much more complex than that of the Bullet Cluster (Jee et al, 2014), a natural explanation in the framework of the κ-model is that the intergalactic gas filling rate of the bubbles containing the visible galaxies is different. With a filling rate of 60% instead of 15% (Bullet Cluster) we have been able to ascertain that the lensing is no longer centered on the visible galaxies.

Conclusion

Effectively tested MOND and MOG models were designed to understand the dynamics of the Universe without dark matter. By contrast, the dark matter paradigm is an ad hoc concept, where the dark matter content has to be adapted to each situation and consequently, has no predictive value. However, and paradoxically enough, the DM paradigm constitutes a near unanimous recognition among the astrophysicists. The main reason is that the DM paradigm is very easy to use and effectively works at all times due to its extreme flexibility.

On the other hand, in the framework of the κ-model, the single issue of the baryonic mass needs to be sufficient to understand the dynamics of galaxies and that of the galaxy clusters, without dark matter or artificial ingredients, such as the introduction of new parameters into the calculations. Specifically, the κ-model aims to determine a one-to-one relationship directly linking the sole observational data, i.e. the estimated mean density, to either the spectroscopic velocities in the galaxies or X-ray temperatures in galaxy clusters. The κmodel is a MOND-type model, and a behavior very similar to MOND needs to be found for the galaxy clusters. In reality, even though the shape of curves are effectively similar, the κ-model curves are strongly parallel-displaced and systematically move closer to the observational curves than MOND. The agreement between the κ-model prediction and the observational data for the total mass of the hot gas in a galaxy cluster is satisfactory when the mass ratio DM/B is less than 10, but for values exceeding this ratio the observational total mass of gas cannot be effectively predicted in the framework of an isothermal model. In this case, a lower temperature for the hot gas is predicted in the inner regions of the galaxy clusters.

Finally, the ordered series is as follows: the κ-model has no outer parameter (except the internal parameters relative to the system type (individual galaxies or galaxy clusters), i.e. mean density and temperature); MOND is similar with just a single outer parameter; MOG has two outer parameters (unfortunately, it is in fact a disguised gas-mass dependent multifit parameter for the galaxy clusters, i.e. the results are already included in the hypotheses), and DM is an artificial and ad hoc procedure that works at all times. Our proposed logical program was developed to fit the observational curves with only the observational parameters in a self-consistent manner (i.e. through the triplet mean densities, spectroscopic velocities, and temperatures), and this ultimate goal is now possible in the framework of the κ-model.

A Motion of an individual particle

This appendix is a summary of the preceding papers (Pascoli and Pernas, 2020;Pascoli, 2022 a, b). In a current manner for a terrestrial observer E the motion of any particle of mass m taken in an extended system of particles is simply determined by the following dynamic equation:

d(mκ E σ) dt = f E ( 16
)
where σ is the bare velocity of the particle (here, the derivative of the bare position given by the vector σ) and f E is the sum of forces applied to it 7 , as evaluated by the considered observer (the coefficient κ E is defined below). Additionally, the κ-model assumes that at a very large (galactic) scale, the immediate environment of a particle plays a role in the determination of its motion. For instance, the particle "feels" the presence of the other particles in such a way that its velocity is modified by a scale factor κ that is proportional to the mean density ρ; more exactly, this factor is proportional to the logarithm of the mean density. In a more concrete manner, all structures of the Universe are usually inserted in a homogeneous and isotropic base space Σ with coordinates σ. A phase space Π(σ, σ) is attached to Σ (the point on the letter designates the time derivative). However, every observer O has his own measuring gauge, which is dependent on his environment. Thus, this observer applies an apparent isotropic and homogeneous dilation to the phase space Π(σ, σ) 8

(σ, σ) -→ (κσ, κ σ) ( 17
)
where κ is a renormalization coefficient (for the terrestrial observer κ is denoted κ E ). For any observer O, the coefficient κ is considered to be independent of the point, and any other observer O ′ , placed in a distinct environment, follows the same reasoning, but with the essential difference that this one applies a coefficient κ ′ ̸ = κ. The basic idea is that any observer does not directly access to the phase space Π but rather visualizes an apparent homothetic replica κΠ. The κ-model is reduced to this sole operation. Consequently, this model is not a "theory" by itself. An underlying theory is needed, such as the Newtonian mechanics or the general relativity.

For any extended objects (a galaxy or a galaxy cluster) with a mean density profile with both definite upper ρ M and lower ρ m bounds, such as ρ M > ρ m , there exists a universal relationship for κ; this factor is linked to the local mean density ρ, as follows:

κ M κ = 1 + Ln( ρ M -ρ m ρ -ρ m ) ( 18
)
7 The force f is the resultant force of the gravitational, electrostatic and magnetic forces, and also of the centrifugal and Coriolis forces if the frame of reference is not inertial.

8 Once this operation is achieved, the quantity κσ forms an inseparable unit R, where κ and σ are no longer separately measurable for any observer (similarly for κ σ, which forms an inseparable unit Ṙ). Within the κ-model framework the Newtonian law is formally maintained, as follows: m dP dt = -∇ R [(Φ(R)] where P = m Ṙ is the impulsion of the particle and Φ is the gravitational potential measured in situ. For all these reasons the variation of κ is hidden to any observer, and the space Σ of vectors σ is not reachable. However, the ratios of the type κ ′ κ are always measurable for any observer, an operation which helps to exchange the information with other observers (Pascoli, 2022).

For a galaxy or a cluster of galaxies ρ M ≫ ρ m ≃ 0, this relationship simplifies in the form κ M κ = 1+Ln ( ρ M ρ ) (the function ρ has a definite upper bound ρ M and ρ m ≃ 0) 9 . The quantity ρ M is a reference value for the mean density. For any local observer the density measured, ρl , is equal to the following:

ρ 1 + Ln ρ M ρ 3 (19)
The density ρl , measured in situ, can be considered as real, whereas the density ρ measured by a terrestrial observer is only an "apparent" quantity. Eventually, the mean density around a point, labelled by σ, can be obtained by integrating over a suitably sized volume ω surrounding this point as follows:

ρ(σ) = ω d 3 σ ′ w(σ -σ ′ )ρ(σ ′ ) (20)
For simplicity, we can assume that the spread function w (σ) is isotropic and Gaussian, i.e.

w (σ) (πδ 2 ) (-3/2) e -σ 2
δ 2 . This spatial averaging is used to smooth the strongly varying density distribution of matter in a galaxy over distances of a few parsecs around each point. 9 However, in order to perform the analysis of the CMB in the framework of the κ-model, the complete form has to be used because ρ

M ≃ ρ m ( ρ M -ρ m ρ m ∼ 10 -5 ).
With these conditions in view, the dynamic equation ( 16) becomes 10

(mκ σ) dt = f (25)
or, the mass m is invariable, as follows:

m d(κ σ) dt = f (26)
As for the second member of the dynamic equation, the real force f is obtained from the apparent force f E (evaluated by a terrestrial observer) by applying the same transform (17), i.e. that which is measured by a local observer, where the force is really acting. In the following, and in order to simplify the writing, ρ will be replaced by the symbol ρ.

B Case of an extended set of points

As discussed in the preceding paragraph, in the κ-model, the dynamic equation, formally written for an infinitesimal element of matter of mass dm and subjected to a force df , is as follows:

10 Equation ( 16) is evidently not correct, and the well known consequence of its misuse constitutes the missing mass problem. The solution acknowledged today by a common consensus is to add dark matter. Thus, we have the following relationship:

m d dt (κ E σ) = f E + f DM = (1 + α)f E (21)
The observations provides α ∼ 1 -10, following the type of extended objects under consideration (the coefficient α, or more generally the profile of the function α when α varies in the objects, is predominantly adapted in an ad hoc manner). For instance, α is small in dense systems, such as globular clusters and, at the opposite side, high in diffuse low-density galaxies.

In the framework of the κ-model formalism, equation ( 16) can be rewritten as follows:

m d dt (κ σ) = κ E κ 2 f E = f (22)
where the second equality is consistent with the following relationship:

f C-→P = -GM m κCP κ 3 CP 3 = κ E κ 2 f E,C-→P (23)
This relationship applies for a test particle P of mass m under the action of a central attractor C of mass M (σ = CP ). The force is evaluated where the particle P resides (the real force is measured in situ). For a particle placed on a circular orbit and subjected to a central attractor of mass M , the coefficient κ is quasi-constant, and the following useful formula can be used:

v 2 spectro = (κ σ) 2 = κ E κ ( GM κ E σ ) ( 24
)
Two very distinct interpretations of the magnification factor κ E κ are then possible: an apparent magnification of the gravitational constant G (as in MOG) or an apparent magnification of the attractive mass M (as in DM with the identification 1 + α = κ E κ ). Note that v 2 spectro is observer-independent (we can replace κ E with any value). where κ is a local normalization coefficient applied to the spatial lengths and df is measured in situ, i.e. where the element of matter resides. However, this equation has no meaning without a reference frame predefined in advance. First, the internal motions of any particle (a star, a nebula) in a galaxy needs to be studied in the reference frame R A , in which the barycenter of the galaxy is at rest. Such a reference frame can be built by considering a collection S {A i } of open sets U A i which entirely covers the galaxy. A couple (A i , κ A i ), composed of an observer A i and a normalization coefficient applied to the spatial lengths, κ A i , is associated with each open set U A i . The observers A i are assumed to be at rest relatively to each other, as follows:

v(A i ) A j = 0 (28)
Then, the collection S {A i } composes a reference frame R A . The observers A i do not use the same unit of length along R A , but they are at rest relatively to each other as imposed by eq. ( 28). Next, for any point P in the set U A i , the following relationship is used:

v(P ) A i = κ A i σ(P ) A i (29)
Second, the motion of any galaxy as a whole has to be treated as if this one was a point confounded with its barycenter. Any galaxy evolves in a cluster of galaxies. A second collection S {B l } of open sets U B l is considered and once again covers the entire galaxy cluster.

The observers B l are assumed to be at rest relatively to each other, as follows:

  Column (3) ICM central mass density Column (8) Newtonian dynamic mass integrated to rout Column (4) model parameter Column (9) MSTG dynamic mass integrated to rout Column (5) model core radius parameter Column (10) convergent MOND dynamic mass Column (11) Mκ integrated to rout

  Figures 5a and 5b are reproduced from Brownstein and Moffat

Figure 7 Figure 8

 78 Figure 7 Comparative plots of the lensing diagram. The crosses indicate the positions of the centers of the different distributions (hot gas in red and groups of visible galaxies in blue)

  Figure 9

  By comparing columns 9 and 11 in Table1, MOG and κ-model provide reasonably close values to each other for the masses integrated to rout, respectively M M ST G and M

											Mκ
										10)	(11)
	A0085	6.90	0.34	0.532	58.5	2,241	1.48	9.02	1.15	1.83	0.77
	A0119	5.60	0.03	0.675	352.8	1728	0.73	6.88	0.73	1.76	0.60
	A0133	3.80	0.42	0.530	31.7	1417	0.37	3.13	0.28	0.55	0.27
	NGC507	1.26	0.23	0.444	13.4	783	0.05	0.48	0.02	0.04	0.04
	A0262	2.15	0.16	0.443	29.6	1334	0.26	1.39	0.11	0.13	0.12
	A0399	7.00	0.04	0.713	316.9	1791	0.90	9.51	1.07	3.07	0.82
	FORNAX	1.20	0.02	0.804	122.5	387	0.009	0.373	0.011	0.102	0.026
	NGC1550	1.43	0.15	0.554	31.7	632	0.034	0.548	0.024	0.086	0.047
	A1060	3.24	0.09	0.607	66.2	790	0.07	1.69	0.10	0.50	0.21
	A1367	3.55	0.03	0.695	269.7	1234	0.27	3.19	0.26	0.75	0.40
	MKW4	1.71	0.57	0.440	7.7	948	0.09	0.78	0.05	0.08	0.10
	ZwCl1215	5.68	0.05	0.819	303.5	1485	0.59	7.15	0.72	2.5	0.91
	NGC4636	0.76	0.33	0.491	4.2	216	0.001	0.088	0.001	0.019	0.011
	A3526	3.68	0.29	0.495	26.1	1175	0.20	2.35	0.17	0.45	0.24
	A3266	8.00	0.05	0.796	397.2	1,915	1.22	12.82	1.56	4.79	1.12
	A3395s	5.00	0.03	0.964	425.4	1223	0.32	5.77	0.49	2.34	0.50
	COMA	8.38	0.06	0.654	242.3	1954	1.13	11.57	1.38	3.81	0.99
	A2065	5.50	0.04	1.16	485.9	1302	0.49	8.01	0.76	3.83	0.69
	A2142	9.70	0.27	0.591	108.5	2537	2.39	15.93	2.32	4.36	1.37
	A2244	7.10	0.23	0.607	88.7	1773	0.84	8.36	0.92	2.45	0.71
	UGC03957	2.58	0.09	0.740	100.0	764	0.08	1.57	0.09	0.47	0.13
	S636	1.18	0.01	0.752	242.3	742	0.06	0.65	0.03	0.09	0.05
	M49	0.95	0.26	0.592	7.7	177	0.001	0.109	0.002	0.041	0.009
	A1689	9.23	0.33	0.690	114.8	1898	1.23	13.21	1.61	5.14	1.13
	A1800	4.02	0.04	0.766	276.1	1284	0.34	4.14	0.36	1.15	0.38
	A1914	10.5	0.22	0.751	162.7	1768	1.08	15.21	1.79	7.44	1.31
	NGC5813	0.52	0.18	0.766	17.6	166	0.001	0.072	0.001	0.021	0.006
	NGC5846	0.82	0.47	0.599	4.9	152	0.001	0.082	0.001	0.031	0.007
	A2151w	2.40	0.16	0.564	47.9	957	0.12	1.42	0.09	0.25	0.12
	TRIANGULUM	9.60	0.1	0.61	196.5	2385	1.98	15.22	2.11	4.48	1.32
	OPHIUCHUS	10.3	0.13	0.747	196.5	1701	0.91	14.11	1.59	6.91	1.21
	ZwC174	5.23	0.1	0.717	163.4	1354	0.43	5.49	0.50	1.79	0.47
	A3888	8.84	0.1	0.928	282.4	1455	0.71	12.61	1.33	7.14	1.08
	HGC94	3.45	0.11	0.514	60.6	1237	0.24	2.40	0.19	0.43	0.21
	RXJ2344	4.73	0.07	0.807	212.0	1222	0.34	4.97	0.43	1.78	0.43

κ . By comparing columns 7 and 11 in Table

1

, the agreement between M κ and M gas is fairly good, and it is clear that, in most cases, the κ-model does not necessitate dark matter in the outer regions of the galaxy clusters. In Table

The most complete theory regarding the MONDian world is that of Solznik and Skordis (2021), but the formalism EMOND is much easier to manipulate (the theory of Solznik and Skordis belongs to the large class of TeVeS theories, such as MOG with many free parameters and some disguised aspect of DM).

Further information is summarized in the two appendices A and B placed at the end of this paper. Initially, we start with an isotropic and homogeneous base space Σ. In this space the light propagates in straight line. However, any observer is located at the center of a homogeneous and isotropic universe κΣ, homothetic to the base space Σ; the coefficient κ is dependent on the environment of this observer (the mean density of matter surrounding the observer). Unfortunately, it is very difficult to have a global view of the situation with a unique R

-type space. Rather, the image needs to be that of a fiber bundle Σ × κ with Σ as the base and κ as the fiber (specifically, each observer projects all structures present in the Universe on his own stratum labelled by κ). Within the framework of this fiber bundle, DM effects can be re-interpreted as a κ-lensing, whose the Bullet Cluster is an illustrative example (paragraph 3).3 This value corresponds to the mean density of the baryonic matter detected in the solar neighborhood (or, correspondingly, ∼ 70 M ⊙ pc -2 , assuming a vertical thickness ∼ 1 kpc(Famaey, and McGaugh, Fig.19, 2012)).

The temperature T is measured in situ and is observer-independent, likewise for the dispersion velocities σ r , σ θ , σ ϕ (for instance in eq. 9 of Brownstein and Moffat, 2006 : σ 2 r = (k B T )/µm p is observer-independent, as is the gravitational potential Φ which is also measured in situ).

After this short period of time, the κ-model predicts that the lensing diagram will no longer be centered on the groups of visible galaxies, but on the hot gas and will likely resemble Figure7a. Thus, we see an instantaneous phase of a rapidly evolutionary process.

A multiplicative factor, κ E κ M , needs to be applied to α (x, y); however, the factor is a global and constant coefficient independent of x, y, which does not affect the relative position of the peaks on the lensing diagram.

Data availability statement: The author confirms that the data supporting the findings of this study are available within the article and the reference list.

Conflicts of Interest:

The author declares no conflict of interest.

Batik, I. & Zhao, H., 2022, Symmetry, 14, 1331Bartelmann, M., & Schneider, P., 2001, Physics Reports, 340, 291 Bradač, M., Schneider, P., Lombardi, M., & Erben, T., 2005, A & A, 436, 39 Bradač, M., Clowe, D., Gonzalez, A.H., Marshall, P., Forman, W., Jones, C., Markevitch, M., 2006, ApJ, 652, 937 The law of addition of velocities is applied, as follows:

with the following: v(P

Based on these relationships, the expression (25) can now be properly rewritten (df i represents the internal forces, df inertial , df e are the inertial forces and the external forces, respectively, applied on the mass dm), as follows:

This equation treats the case of the internal motions within the galaxy (for instance, the rotation). We can integrate eq. 34 over the volume Ω of the galaxy, as follows:

For an isolated (F e = 0, F inertial = 0) and symmetric (F i asym = 0) galaxy, the following relationship is obtained:

where the constant K can eventually considered nul.

The dynamic equation for a symmetric galaxy studied as a whole and covered by a set B l is as follows:

where M is the total mass of the galaxy (M = Ω dm)). The equation ( 38) can be used to treat the motion of a galaxy (viewed as a point here) in a galaxy cluster. The dynamic equation (34) relative to any point P of mass m in the galaxy becomes:

The force f e -m M F e is the tidal force produced by an outer mass. Equation ( 40) can be used to treat the internal motions in a galaxy.

11 For an isolated and asymmetric galaxy for which F i asym ̸ = 0, this equation becomes the following:

In this case, the galaxy can be auto-accelerated.

12 Let P contained in the set A i and P ′ contained in the set A i ′ The internal forces are calculated as follows:

The force needs to be measured in situ (i.e. respectively at P for f P ′ -→P and at P ′ for f P -→P ′ ) and, in this case, the measured force is a well-defined quantity and considered observer-independent.


Keywords: Storage of the ℓ-BFGS matrices, 8 5 Numerical results 9, 9 5, 2 Performance of the ℓ-BFGS preconditioners, 11

M1CG1 -A solver of symmetric linear systems by conjugate gradient iterations, using/building a BFGS/ℓ-BFGS preconditioner Version 1.2 (

Goal of the software

The module m1cg1 solves for x ∈ R n the linear system

Ax = b, (1.1)
by a possibly preconditioned Fletcher-Reeves conjugate gradient (CG) algorithm [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF][START_REF] Fletcher | Function minimization by conjugate gradients[END_REF][START_REF] Bonnans | Numerical Optimization -Theoretical and Practical Aspects[END_REF]5]. In this system, the matrix A ∈ R n×n and the vector b ∈ R n are supposed to be known. It is assumed that the matrix A is symmetric. Therefore, solving (1.1) is equivalent to finding a stationary point of the minimization problem inf

x∈R n 1 2 x T Ax -b T x . (1.2)
To be sure to be able to solve the linear system (1.1) [resp. the minimization problem (1.

2)] for an arbitrary b ∈ R n , the matrix A must be nonsingular [resp. positive definite]. When A is positive definite, (1.1) and (1.2) are equivalent. If the CG algorithm encounters a quasi-negative curvature direction, which is a conjugate direction v

1 satisfying v T Av ν v 2 , (1.3)
for some given threshold ν 0, m1cg1 just stops. In (1.3) and below, • denotes the ℓ 2 or Euclidean norm, i.e., x := ( i x 2 i ) 1/2 . The way A is actually stored is irrelevant for m1cg1, since the algorithm requires only products Av of A times various vectors v and since this operation is supposed to be done in a user-supplied subroutine. This gives the possibility to use the code even when A is not known explicitly.

M1cg1 has the nice feature of being able to build itself a preconditioning matrix while it solves the linear system (1.1). This matrix is an approximation of A -1 and can then be used to precondition a subsequent linear system with the same matrix A or a matrix close to A, and a different right-hand side. The preconditioning matrix is formed by BFGS updates [START_REF] Dennis | Numerical Methods for Unconstrained Optimization and Nonlinear Equations[END_REF][START_REF] Bonnans | Numerical Optimization -Theoretical and Practical Aspects[END_REF]5]. It is also possible to form an ℓ-BFGS preconditioning matrix [10,[START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF]8], so that the solver is adapted to problems of very large dimension. M1cg1 can also improve a previously built preconditioner; this feature is however presently limited by the fact that the method used to improve the preconditioner (i.e., full BFGS updates or ℓ-BFGS updates with a given selection procedure) must be the same as the used to build the previous preconditioner.

Another feature of m1cg1 is to be able to solve in parallel two linear systems: (1.1) and Ax = b.

(1.4)

The latter system uses the same matrix A as in (1.1), but another right-hand side b. The second linear system (1.4) is solved using the conjugate directions generated by the CG iterations for solving (1.1). This means that the two systems are solved without increasing the number of matrix-vector products. This feature is useful for the truncated SQP algorithm in constrained optimization [START_REF] Chauvier | A truncated SQP algorithm for solving nonconvex equality constrained optimization problems[END_REF].

The package

The module m1cg1 includes the following subroutines:

-m1cg1: the main subroutine, which has to be called and that calls m1cg1a after having structured the available memory; -m1cg1a: the actual solver; dbfgsi, dysave, dbfgsl, deuclid, dsyev, l2bfgs: subroutines used to build a BFGS preconditioning matrix.

On his part, the user has to provide:

-a subroutine, called mvprod by m1cg1, which computes matrix-vector products Av for arbitrary vectors v ∈ R n ; the method of storage for the matrix is left to the user convenience.

3 Usage

The optimization subroutine m1cg1

The subroutine definition statement is the following.

subroutine m1cg1 (mvprod, n, x, b, a, xx, bb, two, epsneg, restol, absrel, iter, imp, io, imode, mode, w, nw, bfgsp, pmat0, npmat0, m0, ilm0, nilm0, wlm0, nwlm0, bfgsb, pmat1, npmat1, m1, ilm1, nilm1, wlm1, nwlm1, select, izs, rzs, dzs)

In the description that follows, an argument flagged with (I) means that it is an input variable, which has to be initialized before calling m1cg1, an argument flagged with (O) means that it is an output variable, which has a meaning on return from m1cg1 only, and an argument flagged with (IO) is an input-output argument, which has to be initialized and has a meaning after the call to m1cg1. Arguments of the type (O) and (IO) are usually modified in m1cg1. Therefore, they should not be fortran constants! mvprod: Generic name of the user-supplied subroutine (see Section 3.2). This name must be declared external in the program calling m1cg1.

n (I): integer variable. Dimension n of the problem. It must be 1.

x (IO): double precision array of dimension n.

On entry: the initial guess for the solution x of (1.1).

On return: the approximate solution of ( • If absrel = 0, restol is assumed to give the required precision on the absolute residual. In other words, if during the CG iterations, the absolute residual satisfies Ax -b restol, the algorithm stops, since (1.1) is considered to be solved by the current x with enough precision.

• If absrel = 1, restol is assumed to give the required precision on the relative residual. In other words, if during the CG iterations, the relative residual satisfies

Ax -b Ax 1 -b restol,
the algorithm stops, since (1.1) is considered to be solved by the current x with enough precision.

When two = .true., there is no stopping test based on the behavior of the conjugate direction algorithm on (1.4). If (1.4) is solved before finding the approximate solution of (1.1), x is modified by adding terms that are zero in exact arithmetic and that should be close to zero on a computer. On return: when m1cg1 stops with mode = 0,

restol = Ax -b if absrel = 0 Ax-b Ax 1 -b if absrel = 1,
where x is the final iterate.

absrel (I): integer variable, which can take the value 0 or 1. It specifies how M1cg1 has to understand the value of the argument restol. See the description of that argument for the role of absrel.

iter (IO): integer variable. On entry: maximal number of iterations authorized. On return: actuel number of iterations performed During the run, iter is used as iteration counter, so that the user of m1cg1 can know the iteration number each time mvprod is called. It is important, however, not to change the value of iter while m1cg1 is running.

imp (I): Integer variable that controls the printings via the I/O channel io (it is another argument). = 0: nothing is printed. 1: printings on entry and return; m1cg1 prints also error messages. 2: one line is also printed at each CG iteration. It gives iter: iteration counter; cost: the cost 1 2 x T Ax -b T x; -|r| or |r|/|r1|: the absolute (if absrel = 0) or relative (if absrel = 1) norm of the residual (controlled by restol); does not decrease monotonically; -<(r,Pr): angle between the residual r = Ax -b and the preconditioned residual P r; should be 0 if bfgsp = 0; -(Av,v)/|v|^2: Rayleigh quotient of A along the current conjugate direction v; conj: cosine of the angle between 2 succesive residuals (for the scalar product defined by the preconditioning matrix P ); should be very close to 0; alpha: step-size; always positive; -|x|: norm of the current approxinate solution x; should increase if bfgsp = 0; -S: selection flag (only present if an ℓ-BFGS preconditioner is built, i.e., if bfgsb = 2); 'S' means that the (Av, v)-pair has been selected, '-' means that it has not been selected. If two = .true., cost and (|r| or |r|/|r1|) is also printed for the linear system (1.4). io (I): integer variable. Output channel for printings. imode (I): integer variable. Input mode of m1cg1.

imode(1):

• If imode(1) = 0, the value of x is set to 0 by m1cg1 and the CG algorithm starts with x 1 = 0; • If imode(1) = 1, the value in x is taken as starting point x 1 by the CG algorithm.

imode(2):

• If imode(2) = 0, the value of xx is set to 0 by m1cg1 and the conjugate direction algorithm starts with x1 = 0; • If imode(2) = 1, the value in xx is taken as starting point x1 by the conjugate direction algorithm.

imode(3): Specifies the meaning of the parameter epsneg.

• If imode(3) = 0, epsneg is considered to be the value of ν in (1.3).

• If imode(3) = 1, epsneg is considered to be the smallest authorized value for the ratio between the minimum and the maximum Rayleigh quotient encountered during the CG iterations. In other words, the inverse of epsneg is the maximum condition number that the CG iterations are supposed to support. In fact, ν in (1.3) is set to epsneg multiplied by the maximum Rayleigh quotient encountered during the CG iterations. In this case, ν depends on the iteration counter.

mode (0): integer variable. Output mode of m1cg1. = 0: normal terminaison (stop on restol); = 1: a dimension argument has a wrong value; = 2: (not used); = 3: A is probably indefinite (Rayleigh quotient less than ν); = 4: maximum number of iterations has been reached; = 5: (not used); = 6: the initial residual is zero; hence x 1 is solution of (1.1); = 7: error in the subroutine dysave, this may be due to a wrong value of select, a scalar product y T s < 0 (rounding error?), or a negative value of m1, = 8: significant increase of the cost, due to rounding errors. = 0: no preconditioning information; = 1: a full symmetric positive definite preconditioning n×n matrix is available in the structure (pmat0, npmat0); = 2: a limited memory BFGS preconditioning matrix is available in the ℓ-BFGS structure (m0, ilm0, nilm0, wlm0, nwlm0); this structure can be built up by a preceding call to m1cg1. pmat0, npmat0 (I): If bfgsp = 1, this structure must contain preconditioning information. It is not used otherwise.

pmat0: is a double precision array of dimension npmat0. It must contain an n × n symmetric positive definite matrix P approximating the inverse of A. The (i, j)th element P ij of P is supposed to be in pmat0(i,j), so that the leading dimension of pmat0 must be n. npmat0: integer variable. It is the dimension of pmat0 and must be n 2 . m0, ilm0, nilm0, wlm0, nwlm0 (I): If bfgsp = 2, this structure must contain ℓ-BFGS preconditioning information. It is not used otherwise.

m0: integer variable. It is the number of (Av, v)-pairs stored in the structure. ilm0: integer array of dimension nilm0. nilm0: integer variable, giving the dimension of ilm0. There must hold nilm0 m0 + 4.

wlm0: double precision array of dimension nwlm0 that contains the (Av, v)pairs. nwlm0: integer variable, giving the dimension of wlm0. There must hold nwlm0 m0(2n+1) + 1. bfgsb (I): integer variable. It specifies the type of BFGS preconditioning matrix that m1cg1 has to build while it solves (1.1). = 0: M1cg1 will not build any preconditioning matrix. = 1: M1cg1 has to build a full BFGS preconditioning matrix in the BFGS structure (pmat1, npmat1). If pmat1(1, 1) 0, the updates start from scratch, otherwise pmat1 is updated. = 2: M1cg1 has to build a limited memory BFGS preconditioning matrix in the ℓ-BFGS structure (m1, ilm1, nilm1, wlm1, nwlm1). If one of the following two conditions

ilm1(1) 0, ilm1 (1)
> 0 and select = 1, is satified, the updates start from scratch, otherwise the ℓ-BFGS structure (m1, ilm1, nilm1, wlm1, nwlm1) is updated, provided the selection procedure specified by the argument select is the same as the one that was used to built that previous preconditioner.

pmat1, npmat1 (O): This structure is filled in if bfgsb = 1 and not touched otherwise. Same meaning as for the structure (pmat0, npmat0). m1, ilm1, nilm1, wlm1, nwlm1 (O): This structure is filled in if bfgsb = 2 and not touched otherwise. Same meaning as for the structure (m0, ilm0, nilm0, wlm0, nwlm0). select (I): integer variable. In the case when bfgsb = 2, this argument monitors the selection of the (Av, v)-pairs used to build the ℓ-BFGS preconditioning matrix. The meaning is as follows. = 0: FIFO (First In First Out) selection: the oldest pair is discarded and the new one is saved. = 1: Mexican selection [9]: the pairs are distributed uniformly according to the iteration counter of the current CG run. = 2: Rayleigh quotient selection: it is tried to have pairs with a Rayleigh quotient distributed as uniformly as possible on the logarithmic scale. In this version of the software, if it is desired to improve a previously built preconditioner, stored in the structure (m1, ilm1, nilm1, wlm1, nwlm1), the selection procedure specified by select must be the same as the one that was used to build that previous preconditioner. izs, rzs, dzs: Variables or arrays, whose type is respectively integer, real, and double precision. These variables/arrays are not used by m1cg1. They are just passed to the user-supplied subroutine mvprod, so that the user can use them as wished.

The user-supplied subroutine mvprod

Mvprod is the user-supplied subroutine that computes the matrix-vector products Av for arbitrary vectors v, which are actually conjugate directions. M1cg1 assumes that this subroutine is defined by the following statement.

subroutine mvprod (n, a, v, av, izs, rzs, dzs)

Meaning of the arguments:

n, a, izs, rzs, dzs: same meaning as the arguments of m1cg1 with the same name.

v (I): double precision array of dimension n. A vector v ∈ R n given to multiply the matrix A.

av (O): double precision array of dimension n. Vecteur Av, product of the matrix A by the vector v, that mvprod has to compute.

Implementation details 4.1 Calling sequence

Necessarily, the instructions preceding the call to m1cg1 must include the following items:

1. the declaration in external of the name of the user-supplied subroutine doing matrix-vector products (named mvprod by m1cg1);

2. the calculation of the starting point x 1 (and possibly x1 );

3. the initialization of the input arguments of m1cg1;

4. the call to m1cg1.

Storage of the ℓ-BFGS matrices

The ℓ-BFGS matrices are stored in the structures formed of the variables mX, ilmX, nilmX, wlmX, and nwlmX, where the X stands for 0 (preconditioning matrix) or 1 (constructed preconditioner). A few has been said about these variables in the description of the subroutine m1cg1. Here are more details.

mX: integer variable. It is the maximal number of (Av, v)-pairs that can be stored in the structure.

ilmX(1): integer variable. It is the number of (Av, v)-pairs that have been selected for storage in the ℓ-BFGS structure, possibly replacing older pairs. ilmX(2): integer variable, specifying the selection procedure that has been used to select the (Av, v)-pairs stored in the ℓ-BFGS structure; see the variable select.

ilmX(3): integer variable. It is a cyclic pointer to the oldest (Av, v)-pair stored in the structure. Ones it has been initialized, its value is between 1 and mX.

The word cyclic means here that, when its value is mX and it is incremented, its new value is 5 Numerical results

Performance of the BFGS preconditioners

This section aims at experimenting the performance of the BFGS preconditioner built by M1cg1. The experiments are done for solving a randomly generated linear system In the first experiment, a BFGS preconditioner is formed during a single run made of n k = 50k (k = 0, 1, . . . , 10) unpreconditioned CG iterations. This preconditioner is then used to solve the same linear system with preconditioned CG iterations. The number of iterations to get an absolute precision of 10 -10 (restol=1.d-10 and absrel=0). BFGS preconditioner that was previously built using the number n k of (Av, v)-pairs given in abscissa. The solid blue curve refers to the first experiment: the preconditioner is generated by a single CG run. The dashed red curve refers to the second experiment: the preconditioner is generated/improved by preconditioned CG runs of 50 iterations. is shown in figure 5.1 as a function of n k (solid blue curve): there is a clear decrease of the number of preconditioned CG iterations when n k increases. Observe, however that for n k = 100 (the order of the linear system), the number of preconditioned CG iterations to solve the linear system is not 1, but still 180. This is not in concordance with the theory (see [5] for example) and, by lack of a deeper analysis, we attribute this to rounding errors (to back this opinion, note that a solution to the linear system by the unpreconditioned CG algorithm requires here 565 iterations to reach the required precision, and not 100 as claimed by the theory, which assumes that the computation is done in exact arithmetics).

In the second experiment, a BFGS preconditioner is also formed using n k = 50k (k = 0, 1, . . . , 4) CG iterations, but these iterations are not taken from a single unpreconditioned CG run, but out of k preconditioned CG runs, each of which being formed of 50 iterations. To be more precise, a first preconditioner is obtained using the first 50 (Av, v)-pairs generated by an unpreconditioned CG run. Next, the preconditioner is improved by the first 50 (Av, v)-pairs generated by a CG run preconditioned by the matrix obtained after the first run. And so on, k times. The final preconditioner is then used to solve the same linear system. The required number of iterations as a function of n k is shown by the dashed red curve in figure 5.1. We quote two facts: as for the solid blue curve, there is a clear decrease of the number of preconditioned CG iterations when n k increases; in addition the number of iterations is significantly less than the one obtained in the first experiment with the same n k .

Performance of the ℓ-BFGS preconditioners

This section aims at experimenting the performance of the ℓ-BFGS preconditioners built by M1cg1, those obtained when bfgsb is set to 2 and select is set to 0, 1, or 2. The experiments are done for solving the same randomly generated linear system (1.1) as in section 5.1, in particular, the order of the system is n = 100 In the first experiment, we compare the performance of the BFGS and ℓ-BFGS preconditioners, when these are formed using m = 10k updates (for k = 0, 1, . . . , 10). In this experiment, the ℓ-BFGS preconditioner is allowed to store 100 (Av, v)-pairs, so that the pair selection option, monitored by select, is ineffective. The comparison is illustrated in figure 5.2. A first observation is that the ℓ-BFGS preconditioner a BFGS (dashed blue curve) or ℓ-BFGS (solid blue curve) preconditioner that was previously built using the number of (Av, v)-pairs given in abscissa.

(solid blue curve) improves significantly when the number of updates increases; this observation was already observed in section 5.1 for the BFGS preconditioner, but the dashed blue curve shows here that this decrease is not monotone. A second observation is that the ℓ-BFGS preconditioner performs significantly better than the BFGS preconditioner. We attribute this to the diagonal scaling included in the ℓ-BFGS preconditioner, which can be rescaled at each update, but not in the BFGS preconditioner. Note, however, that when m is large, the use of the ℓ-BFGS preconditioner may require more computations than the BFGS preconditioner.

In the second experiment, we investigate the performance of the 3 possible selection procedures, corresponding to the values 0, 1, or 2 of the argument select. For this, we fix to m = 20 the number of (Av, v)-pairs that the ℓ-BFGS preconditioner can store and let M1cg1 select the pairs on an unpreconditioned CG run made of n k = mk iterations, for k = 1, . . . , 100/m, which is the abscissa of figure 5 figure gives the number of iterations required by M1cg1 using the built preconditioner, to solve the linear system to the absolute precision 10 -10 . The results are disappointing, showing that the selection procedures are inefficient in selecting the appropriate pairs. Among the selection procedures proposed by M1cg1, the one based on the Rayleigh quotient (select = 2), looks better than the others. Setting the ℓ-BFGS size to m = 5 or m = 20 provides similar results.

In the third and last experiment, we want to test the ability of M1cg1 to improve an ℓ-BFGS preconditioner by incorporating new (Av, v)-pairs during successive CG runs. For this, an ℓ-BFGS preconditioner is built during m = 20 unpreconditioned CG iterations. Then this preconditioner is updated k -1 times (k = 1, . . . , 5), each time by m preconditioned CG iterations (the preconditioner being the one so obtained before the current run). Finally, the linear system is solved with the resulting preconditioner. underlying the Mexican selection procedure is not adapted to this experiment). In ordinate, one finds the number of iterations required to solve the linear system with the preconditioner built with mk updates (k = 0, . . . , 5), a number given in abscissa.

These last two experiments clearly indicate that more research is needed to improve the selection procedures. In the present development of the software, if it is desired to improve the ℓ-BFGS preconditioner, it is much better to increase its number m of (Av, v)-pairs (provided the memory space allows such an increase) than to make it collecting information on a number of CG iterations larger than m.

  w (I): double precision array of dimension nw. Working zone. nw (I): integer variable. Give the dimension of w. One must have: nw 5n + m0 if two = .false. nw 6n + m0 if two = .true. bfgsp (I): integer variable. It specifies what type of preconditioning information is made available to m1cg1 to precondition the current CG run.

( 1 . 1 )

 11 of order n = 100 with a positive definite matrix A and a right hand side b. The experiments are implemented as test 10 in the program quad.f located in the directory examples of the standard distribution. It can be executed by typing make quad quad 10 which generates Matlab codes that can be run to obtain the curves in figure 5.1 below.

  Figure 5.1: Number of iterations needed to solve a randomly generated 100×100 linear system with a

  . The experiments are implemented as tests 11, 12, and 13 in the program quad.f located in the directory examples of the standard distribution. They can be run by typing make quad quad 11 quad 12 quad 13 which generates Matlab codes that can be launched to obtain the curves in the figures 5.2, 5.3, and 5.4 below.

Figure 5 . 2 :

 52 Figure 5.2: Number of iterations needed to solve a randomly generated 100 × 100 linear system with

Figure 5 . 3 :

 53 Figure 5.3: Number of iterations needed to solve a randomly generated 100 × 100 linear system with an ℓ-BFGS preconditioner, made of m = 20 (Av, v)-pairs, previously built by selecting these pairs on the n k iterations given in abscissa. The curves correspond to a different selection procedure: select = 0 for the solid blue curve, select = 1 for the dashed red curve, and select = 2 for the dotted-dashed green curve.

Figure 5 . 4 :

 54 Figure 5.4: Number of iterations needed to solve a randomly generated 100 × 100 linear system with an ℓ-BFGS preconditioner, made of m = 20 (Av, v)-pairs, previously built by selecting these pairs during 5 runs made of 20 preconditioned CG iterations. The curves correspond to a different selection procedure: select = 0 for the solid blue curve and select = 2 for the dotted-dashed green curve.

  1. ilmX(4): integer variable. It is a cyclic pointer to the newest (Av, v)-pair stored in the structure. Ones it has been initialized, its value is between 1 and mX. ilmX(4+i), for i = 1 to mX: integer variables; ilmX(4+i) is the index of the (Av, v) pair with the ith smallest Oren-Luenberger factor.

nilmX: integer variable; dimension of ilmX; = mX + 3. wlmX(1): double precision variable. It is the first element of an array of dimension (n,mX), containing the Av part of the mX (Av, v)-pairs. wlmX(n*mX+1): double precision variable. It is the first element of an array of dimension (n,mX), containing the v part of the mX (Av, v)-pairs. wlmX(2*n*mX+1): double precision variable. It is a scalar preconditioner. wlmX(2*n*mX+2): double precision variable. It is the first element of an array of dimension mX, containing the Oren-Luenberger factors of the (Av, v)-pairs. nwlmX: integer variable; dimension of wlmX; it must be mX*(2*n+1)+1.


Luciano Da
email: luciano@ifsc.usp.br

F Costa

A similarity approach to estimating the source of geometric and topologic random walks

This work addresses the issue of identifying the source of random walks in geometric and topologic spaces by using a methodology based on the Jaccard and coincidence similarity indices. Uniform random random walks are performed by a set of agents departing from the same source node and performing the same number of steps. The final positions of these agents are then taken into account while attempting to estimate the respective origins. In the case of geometrical spaces, the described methodology involves calculating the Jaccard similarity mean of the coordinates of the final points (also including respective outliers), which is found to have potential for providing accuracy higher than that allowed by the more traditional estimation of the center of mass (arithmetic mean). Topological random walks are addressed by identifying the maximum average coincidence similarity between topological the features (distances to every other node) of the termination nodes and corresponding features of every other node in the network. The respective results are shown to be more accurate than randomly choosing the source nodes. Though the reported developments focus on uniform random walks taking place in Watts-Strogatz networks, extensions to other types of random walks and complex networks provide interesting prospects for further investigations.

Introduction

The physical world can be understood as a complex system composed of several mutually interacting modules, which are often complex systems themselves. Diverse types of dynamics can take place, frequently in a simultaneous manner, within these modules, including linear dynamics such as diffusion and wave propagation, as well as non-linear dynamics, including integrate-and-fire and organized criticality.

A complex system can be decomposed into two parts: (a) the geometric/topologic structure of the interconnection between its parts; and (b) the type of dynamics taking place on (a). An issue of particular importance is then implied, which consists in trying to understand how these two aspects interrelate one another (e.g. [START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF]).

From the above perspective, one problem that has attracted substantial interest because of its theoretical and applied importance consists in, given a specific complex system composed by its structure and considered dynamics, to try to identify the initial state of an observed dynamics. For instance, in the case of a complex system defined by a discrete random walk (e.g. [START_REF] Zinn-Justing | From random walks to random matrices[END_REF][START_REF] Lovász | Random walks on graphs[END_REF][START_REF] Masuda | Random walks and diffusion on networks[END_REF][START_REF] Noh | Random walks on complex networks[END_REF]) taking place on a geometric space or on a complex network, it would be interesting to identify the initial point from which this diffusion dynamics emanates (e.g. [START_REF] Comin | Identifying the starting point of a spreading process in complex networks[END_REF]).

In the present work, a mean of similarities approach is described for estimating the original point from which a set of non-preferential random walks emanate. All random walks are supposed to start simultaneously and to last an identical number of steps N t . Two versions of the suggested similarity-based methodology are presented respectively to random walks taking place in: (i) geometrical random walks taking place in 2D geometrical spaces (extensions to higher dimensional spaces are immediate); and (ii) topological random walks unfolding on complex networks.

In the case of 2D geometrical spaces, the recently described concept of similarity between a data element and a reference set of data elements [START_REF] Da | Similarity means: Seeking for robustness to outliers[END_REF][START_REF] Da | Toward robust data classification[END_REF] is employed in order to determined the similarity mean of the positions of agents which performed random walks with the same length after having departed from a same source node. The basic idea underlying this methodology is that the similarity mean will tend to result near the starting point of the random walks. For comparison purposes, this methodology is compared with taking the center of mass (arithmetic mean) of the same final positions.

The suggested methodology to be applied to topological random walks involves comparing features (topological distances to every other node) of the final node reached by a set of same length random walks emanating from the same source and the corresponding features of every other node in the network. The node yielding the maximum average of maximum similarities observed between the available terminations and the network nodes is taken as the estimated respective source of the random walks. This rationale of this methodology assumes that neighboring nodes in a network tend to share most of their topological characteristics, here quantified in terms of distances to the other nodes in the network, so that the topological features of the termination nodes will tend to be maximally similar, in the average, to the features of the source node. Observe that not every type of network will present this property, which is more typically found in lattices and uniformly random complex networks.

The present work starts by describing some basic concepts and the adopted similarity-based methodology for estimation of the sources of random walks, and proceed by illustrating these methodologies respectively to some case examples. The obtained results support, at least for the considered data and parameter configurations, the potential of the similarity-based approaches for estimation of the sources of random walks.

Basic Concepts

All situations considered in the present work assume that the similarity comparisons take place between scalars and vectors composed exclusively by positive integer or real values. Extensions to generic real values (i.e. including negative values) can be found in [START_REF] Da | Coincidence complex networks[END_REF][START_REF] Da | On similarity[END_REF][START_REF] Da | Multiset neurons[END_REF] and references herein.

The Jaccard similarity between two positive-valued vectors x x x and y y y can be defined [START_REF] Da | Coincidence complex networks[END_REF][START_REF] Da | On similarity[END_REF][START_REF]Jaccard index[END_REF] in terms of the multiset (e.g. [START_REF] Da | Multisets[END_REF]) operations of union (max) and intersection (min) as follows:

J(x x x, y y y) = N i=1 min {x i , y i } N i=1 max {x i , y i } (1)
with 0 < J(x x x, y y y) ≤ 1. The interiority similarity index (also known as overlap, [START_REF] Vijaymeena | A survey on similarity measures in text mining[END_REF]) between the two above vectors can be expressed as:

I(x x x, y y y) = N i=1 min {x i , y i } min N i=1 x i , N i=1 y i (2)
with 0 < I(x x x, y y y) ≤ 1. The coincidence similarity index (e.g. [START_REF] Da | Further generalizations of the Jaccard index[END_REF][START_REF] Da | On similarity[END_REF][START_REF] Da | Multiset neurons[END_REF]) between the two vectors x x x and y y y with positive entries can now be calculated as:

C(x x x, y y y) = [J(x x x, y y y)]

D I(x x x, y y y).

where D > 0 corresponds to a parameter controlling how strict the implemented comparison is. The higher the value of D, the more selective and sensitive the similarity comparison becomes.

As with the two previous indices, we again have that 0 < C(x x x, y y y) ≤ 1.

Given that the Jaccard index cannot account for how much one of the multisets is contained (interior) into the other [START_REF] Da | Further generalizations of the Jaccard index[END_REF], the quantification of this property provided by the interiority index can be incorporated into the Jaccard index through respective multiplication as in the above equation. Therefore, the coincidence similarity index becomes capable of performing more selective and sensitive comparisons, incorporates intrinsic data values noramlization, and providies good resilience to data perturbations and outliers [START_REF] Da | Further generalizations of the Jaccard index[END_REF][START_REF] Da | Coincidence complex networks[END_REF][START_REF] Da | On similarity[END_REF][START_REF] Da | Multiset neurons[END_REF][START_REF] Da | Similarity means: Seeking for robustness to outliers[END_REF].

More recently [START_REF] Da | Similarity means: Seeking for robustness to outliers[END_REF], the concept of similarity between vectors has been extended to address the comparison of a scalar data element x and a reference set of N data elements y y y by using the following expression:

J(x, y y y) = N i=1 min {x, y i } N i=1 max {x, y i } (4)
Given a set S of samples of a random variable X, the respective similarity mean has been defined [START_REF] Da | Similarity means: Seeking for robustness to outliers[END_REF] as the value x that maximizes the respective similarity between the latter value and the values in the set S.

The above definition of similarity between a scalar data element and a reference set of scalar elements can be generalized to M × 1 vector data elements [START_REF] Da | Similarity means: Seeking for robustness to outliers[END_REF][START_REF] Da | Toward robust data classification[END_REF] by using the following expression:

J(x x x, y y y) = M j=1 N i=1 min {x j , y i,j } M j=1 N i=1 max {x j , y i,j } (5)
where y y y corresponds to a matrix containing each of the reference 1 × M vectors as subsequent lines.

In the case of the present work, this equation is employed to obtain the Jaccard similarity mean of the coordinates (x, y) involved in the geometric random walks.

Methodology

The methodology to perform estimation of the starting point of geometric random walks consists in identifying, among all possible points in the 2D continuous region R of interest (represented by a respective finite spatial quantization), the point corresponding to the Jaccard similarity mean of the coordinates of the observed termination points. As suggested in [START_REF] Da | Similarity means: Seeking for robustness to outliers[END_REF], this type of mean corresponds to taking the value in the domain of the data that maximizes its similarity with the points (with the same dimension) in a given reference set.

The Jaccard similarity between a vector and a reference set of vectors is here performed by using Equation 5with M = 2 (the x and y coordinates of the points in the 2D geometrical space).

In the present work, the 2D region R is spatially sampled, so that each of the N discrete points can be considered in this calculation. Each point is therefore represented by a respective 2 × 1 vector [x, y]

T . The reference set consists of a matrix y y y whose each line corresponds to one of the terminations of the random walks, while the columns corresponds to the respective x and y coordinates in the 2D domain of the random walks.

For the sake of comparison, the center of mass of the coordinates of the terminations are also calculated. Though it is known [START_REF] Révész | Random walk in random and non-random environments[END_REF][START_REF] Spitzer | Principles of random walk[END_REF] that the center of mass of 2D random walks tends to the origin of the respective walks, the present work considers practical situations in which additional points may also also present in the reference set which correspond to errors, terminations implied by random walks initiating at other points, noise, or outliers.

In the case of the estimation of the source of random walks in graphs or complex networks, we adopt a methodology which is also based on similarity, more specifically the coincidence similarity index. Because there is no coordinate system typically associated to a graph, the nodes of the given network of interest are henceforth represented in terms of its topological distances to the other N -1 network nodes. Each of the network nodes [ [ [x] [i] is thus represented in terms of a feature vector of the type:

x x x [i] = [d i,1 , d i,2 , . . . , d i,N ] T (6)
where d i,j is the shortest topological distance from node i to node j in the given network.

Once the network nodes have been represented as described above, N a uniformly distributed random walks (non-preferential) are performed in the network, starting at the same node s and successively choosing, at each time step, one of the possible outgoing links of the node where one of the agent is, which is done with identical probabilities. For instance, if the node where one of the agents is currently present contains n outgoing links, the next position of the agent will correspond to the destination node of one of these links, taken with probability 1/n. A total of N t steps is considered for all random walks. Therefore, a set of N a termination nodes are obtained, each represented by the respective vector of distances as described above.

The estimation of the possible source of the random walks is then performed by identifying which of the net-work nodes has the maximum average coincidence similarity of respective feature vectors considering all the available termination nodes. For each of the N candidate nodes, this is done by determining its coincidence similarity, by using Equation 3, with each of the available termination nodes, and then taking the respective sum (or average) to be maximized.

Geometric Random Walks

Henceforth, we consider uniform random walks taking place in a 2D geometrical space. After departing from the same source, at the same time, a set of N a moving agents perform N t steps and then stop. The termination nodes, described by respective coordinates, are then considered as subsidy for inferring the respective starting point.

Figure 1 illustrates one realization of the experiments that have been performed as part of the present work in order to test the above described methodology. We have a set of N a = 30 termination points resulting after 50 steps of uniform random walks initiation from the center (x c , y c ) of the depicted region. A set of outlier points is also incorporated in order to quantify the resilience of the proposed methodology respectively to the presence of outliers or noisy data. The error of the estimation of the source nodes can be quantified in terms of the Euclidean distance between that point and the point (x c , y c ). Figure 2 presents these errors obtained by employing the center of mass (arithmetic average of the termination points coordinates, including outliers) and the coincidence similarity mean. It can be readily verified that the coincidence similarity mean provided the most accurate identification of the sources of the random walks. This has been achieved thanks to the enhanced robustness of this similarity indes [START_REF] Da | On similarity[END_REF][START_REF] Da | Multiset neurons[END_REF] to data perturbations.

Topological Random Walks

In this work, the methodology for estimation of the sources of topological random walks is illustrated respectively lattice Watts-Strogatz networks (e.g. [START_REF] Watts | Collective dynamics of 'small-world'networks[END_REF][START_REF] Barabási | Network Sience[END_REF][START_REF] Newman | Networks: An introduction[END_REF][START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF]) with N = 144 nodes, average degree 12, and rewiring probability p r = 0.05. Figure 3 illustrates the network adopted in the henceforth described experiments.

The respectively obtained coincidence complex network obtained by considering the shortest distance feature vectors is illustrated in Figure 4. The coincidence similarity was calculated considering D = 31.

All the N a random walks initiate at the same node s, and unfold during N t = 5 steps. The thus obtained N a termination nodes are then considered as the reference set for the estimation.

In a similar manner with geometrical random walks, the error of the estimation of the source node can be quantified in terms of its topological distance to the known original starting point. Figure 5 presents these errors obtained for the above describe experiment.

It can be readily verified that the obtained errors are markedly smaller than the original distances between any pair of nodes in the given network, and also substantially smaller than the total number N t = 5 of steps in the

Concluding Remarks

The study of the interrelationship between the topological and dynamical properties of complex systems has motivated continuing research in respective scientific research areas. A particular related problem of special importance regards the identification of the source of random walks while considering a set of termination states reached after a given number of steps as performed by a number of random walk agents.

The present work has addressed the identification of the sources of geometrical random walks by considering Jaccard mean, related to recent developments [START_REF] Da | Similarity means: Seeking for robustness to outliers[END_REF] in which the concept of similarity between a given data element and a set of reference data elements has been proposed. In the case of topological random walks, averages (or sums) of coincidence similarities between each of the the network nodes and the set of reference nodes (including outliers) have been considered.

The obtained results indicate that, at least for the con-sidered data sets and parameter configurations, the similarity methods described in the present work present good potential for estimation of the starting points in geometrical and topological random walks. These results are mostly a consequence of the enhanced selectivity and sensitivity of the Jaccard and coincidence similarity index when compared to several other quantifiers of similarity (or distances). In addition, these two indices also present good potential for comparing data in presence of data perturbations, outliers, and noise. The reported concepts and developments pave the way to a number of related further research. In particular, it would be important to further evaluate their performance respectively other types of random walks, higher dimensional geometrical spaces, parameter configurations, as well as other types, sizes and parameters of considered complex networks. In addition, it would be interesting to investigate how the observed performance varies with the number of agents and number of time steps taken during the random walks.

Figure 1 :

 1 Figure 1: Estimation of the origin of random walks while taking into account the position of 30 agents after 50 time steps (shown in black) in the presence of a group of outlier points, shown in gray. All agents started at central position (50, 50), shown as a blue circle. The original positions estimated by the arithmetic mean and Jaccard mean methods are shown in green and red, respectively. The latter approach resulted more accurate.

Figure 2 :

 2 Figure 2: Histograms of the errors obtained from 100 estimations of the origin of the random walks as in Fig. 1 by using the arithmetic mean (a) and coincidence similarity mean (b) methodologies. The errors correspond to the Euclidean distance from the position of the origin of the random walks and the estimated origins. The latter approach resulted markedly more accurate.

Figure 3 :

 3 Figure 3: Watts-Strogatz complex network considered in the caseexample, containing N = 144 nodes, average degree 12 and rewiring probability of 0.05.

Figure 4 :

 4 Figure 4: Coincidence complex networks obtained from the WS network in Fig. 3 by adopting D = 31. The edge widths are proportional to the respective coincidence values between the topological distances (see text) from each compared pair of nodes.

Figure 5 :

 5 Figure 5: Results from the estimation of the source nodes after topological random walks with length 5 taking place on the WS network in Fig. 3. The histogram of the topological distances between the nodes in the original network is shown in (a), and the histogram of the errors (topological distance between the original and estimated sources) is presented in (b). The fact that the obtained distances are substantially smaller in (b) than in (a) suggests that the similaritybased estimation of the sources performed substantially better than selecting the sources randomly.

Acknowledgments Luciano da F. Costa thanks CNPq (grant no. 307085/2018-0) and FAPESP (grant 15/22308-2).

Observations

As all other preprints by the author, the present work contains preliminary work subject to further revision and validation. Respective modification, commercial use, or distribution of any of its parts are not possible, as this work has author copyright. Many of the preprints by the author are also available in Hal and/or arXiv. This work can also be cited by using the DOI number or article identification link. It is also interesting to check for new versions of this work, as they may become available. Thanks for reading.