Aurélien Beuf

Jean-Noël Gence

Philippe Carrière

Florence Raynal
email: florence.raynal@ec-lyon.fr

Chaotic mixing efficiency in different geometries of Hele-Shaw cells

Keywords: Chaotic advection, Mixing, Hele-Shaw flows

Flows in Hele-Shaw cells are generally laminar, and can be in a first approximation considered as quasi two-dimensional. Therefore, the only way to mix efficiently inside the cell is to generate a periodic flow and thus create chaos. Here we study the influence of the geometry of the cell on mixing efficiency: to this aim, we use three different geometries of the cell (circle, square and rectangle) and we numerically test them by using two protocols of mixing: the first protocol uses syringes, the second uses pumps. In a view of simplicity, this Hele-Shaw flow is modeled as two-dimensional, allowing us to analytically determine the velocity field corresponding to the different shapes. The results are then analyzed in terms of Poincaré sections (appearance, filling rate, homogeneity), and the Lyapunov exponents are calculated. We numerically show that the rectangular geometry leads to a better mixing, but also that the aspect ratio of the rectangle plays unexpectedly no important role on mixing.

Nomenclature c

distance from the center to a corner (square and rectangle) D diffusion coefficient h height of the Hele-Shaw cell h K Kolmogorov entropy i natural base of complex numbers K(k) Jacobi function L side of a square chamber P t=t 0 Poincaré section for initial time t 0 q flow-rate R radius of the circular cell S surface of the horizontal section of the cell t time T period of the flow-field T * dimensionless period of the flow-field w parameter measuring the distance between a point and the center w c (z) complex velocity field at complex point z x, y coordinates in a plane z complex coordinate, z = x + iy α dimensionless pulse volume β parameter for the position a source with respect to the center γ aspect ratio parameter for the rectangle λ Lyapunov exponent µ fraction of covered area τ (w) function for the time needed for a point to go from source to sink ξ intermediate parameter

Introduction

Microfluidic devices are widely used in nowadays' biotechnology, for applications ranging from transport of material in microchannels, sorting of cells, DNA analysis, etc. When chemical reactions come into play (e.g. in many such systems), efficient mixing is required in order to homogenize the solution, so as to favour the meeting of reactants. However, the Reynolds number of flows at the microscale are small, generally of the order of unity or smaller, so that Stokes flows should be commonly expected [START_REF] Stone | Engineering Flows in Small Devices: Microfluidics Toward Lab-on-a-Chip[END_REF].

It is now well-known that chaotic advection is the best way to mix efficiently in laminar flows [START_REF] Aref | Stirring by chaotic advection[END_REF][START_REF] Wiggins | Global Bifurcations and Chaos: analytical methods[END_REF][START_REF] Ottino | The Kinematics of Mixing: Stretching, Chaos and Transport[END_REF][START_REF] Toussaint | A numerical eulerian approach to mixing by chaotic advection[END_REF], all the more reason for microfluidic flows.

Stroock et al. proposed a famous chaotic mixer adapted to microdevices, the so-called staggered herringbone mixer (SHM) [START_REF] Stroock | Chaotic Mixer for Microchannels[END_REF], dedicated to mixing in microchannels. A typical microchannel has height and width of the order of 50 -100 µm, and length of a few centimeters, i.e. two small dimensions and the third one large compared to the others. In a smooth channel the mixing of material between streams in the flow is purely diffusive, so that the length of channel required for mixing can be prohibitely long. The SHM consists in adding specific bas-relief structures (ridges) on the floor of the channel, leading to chaotic trajectories; the mixing length is therefore greatly reduced. In this article we deal with mixing in a different microfluidic configuration, the Hele-Shaw flows, where the height is small (typically 25 -50 µm) and the length and width are "large", of the order of the centimeter.

Hele-Shaw flows can, in a first approximation, be considered as twodimensional. In those conditions, the only way to achieve chaotic advection is to make the flow non stationary; indeed, time-periodicity is enough.

We previously proposed two time-periodic protocols of mixing, one using two pairs of syringes, the other operating with pumps, and showed that they could exhibit chaotic trajectories in a square geometry [START_REF] Raynal | Study of a chaotic mixing system for DNA chip hybridization chambers[END_REF][START_REF] Raynal | Towards better DNA chip hybridization using chaotic advection[END_REF]. Other chaotic protocols in a Hele-Shaw cell were proposed by two other teams: MC Quain et al. [START_REF] Mcquain | Chaotic mixer improves microarray hybridization[END_REF] proposed a stirring protocol using pulsed pairs of syringe-driven source-sink systems in a rectangular chamber and tested it in a view of DNA chip hybridization: they showed that hybridization was significantly improved; then Stremler et al. [START_REF] Stremler | Designing for chaos: applications of chaotic advection at the microscale[END_REF][START_REF] Stremler | A maximum entropy approach to optimal mixing in a pulsed source-sink flow[END_REF] carried out the 2D mathematical analysis of mixing in their flow while the experimental analysis was done by Cola et al. [START_REF] Cola | A pulsed sourcesink fluid mixing device[END_REF]. Another work of importance was done by Hertzsch et al. [START_REF] Hertzsch | DNA microarrays : Design principles for maximizing chaotic, ergodic mixing[END_REF] who showed that flows generated by pulsed source-sink pairs could be studied as linked twist maps; therefore they could relate the flow to mathematically precise notions of chaotic mixing, and proposed a new design to generate a well-mixed flow; they compared our protocol A, together with the protocol by Stremler et al., with their new protocol and showed that their protocol, which has mathematical properties in favor of better mixing, was indeed best in a circular geometry (see also Sturman et al. [14] for more details on the mathematical background).

The idea of this article is not to compare all the protocols of stirring proposed in the literature for Hele-Shaw cells, but to test whether or not the shape of the chamber influences the quality of mixing. To this end, we use our two previous protocols [START_REF] Raynal | Study of a chaotic mixing system for DNA chip hybridization chambers[END_REF][START_REF] Raynal | Towards better DNA chip hybridization using chaotic advection[END_REF], which have been shown to exhibit very different mixing abilities. Fluid flow is created thanks to four vertical pipes located near the corners, which inject or extract fluid, and therefore act as sources or sinks. Indeed, those protocols can be considered as variants of the pulsed source-sink chaotic mixer first proposed by Aref and Jones [START_REF] Jones | Chaotic advection in pulsed source-sink systems[END_REF].

1. The first protocol (called protocol A in the following) relies on periodic injection of fluid thanks to four reversible syringes : during one period T , four steps of same duration occur (figure 1(a)).

2. The second protocol (protocol B) operates with two mono-directional pumps switched alternately. For this protocol, only two steps of same duration occur during T (figure 1(b)).

In all the following, like in our previous works [START_REF] Raynal | Study of a chaotic mixing system for DNA chip hybridization chambers[END_REF][START_REF] Raynal | Towards better DNA chip hybridization using chaotic advection[END_REF], the flow rate q, the height h, and the horizontal section S are chosen such that q/hS = 0.18 s -1 ; the shape of the cell is whether a circle, a square of side L (like in our previous works), or a rectangle of lateral sides γL and L/γ, where we took different values for γ > 1, so that the total surface is, in all cases, equal to S.

Numerical tools 2.1 Numerical model

The velocity field is modeled as purely two-dimensional in a first approximation, using the depth-averaged horizontal velocity field, which, in a Hele-Shaw flow, satisfies the Euler equation. Indeed, the aim of this work is to test the influence of the shape of the domain, and for this purpose such a 2D model is much more flexible than the calculation of the whole 3D-flow, for which a change of geometry requires a new computational grid, together with heavy calculations. In our case the flow between a source and a sink is calculated analytically in a circular geometry using the Milne-Thomson theorem [START_REF] Milne-Thomson | Theoretical Aerodynamics[END_REF]. It is then very easy to obtain the flow in the square or rectangular geometry using a Schwarz-Christoffel transformation [START_REF] Milne-Thomson | Theoretical Aerodynamics[END_REF]. All the analytical details of the 2D calculations are given in appendix A. The streamlines obtained with our model for a source-sink flow are shown in figure 2 for different geometries of the chamber. The trajectories are calculated thereafter using a standard fourth order Runge-Kutta integration. Like in the 3D-simulations, our calculations have been performed without taking into account molecular diffusion.

In our 2D-model, we impose that a particle entering in a sink reenters the flow domain on the same streamline. The underlying aim is to avoid adding artificial mixing other than that imposed by chaotic advection. With this hypothesis, a fluid particle at a given location at the end of an injection/extraction stage of protocol A, which is swallowed inside a pipe during the next stage comes back at exactly the same location when the sink is turned again into a source. In a laboratory experiment or in the 3D calculations however, the flow-field is not expected to be exactly zero in the non-active pipes: because of the flow inside the cell, friction causes in those pipes the formation of a large sequence of eddies, each exponentially weaker than its predecessor as height grows [START_REF] Van Dyke | An Album of Fluid Motion[END_REF][START_REF] Taneda | Visualization of separating stokes flows[END_REF]. Of course this effect is very weak; but since the streamlines are very tight inside an active pipe, a change in a fluid particle location during the "inactive" stage, leads the fluid particle to come back nearly, but not exactly at the same location. Even though this effect is small, it may affect regular regions visible on the 2D Poincaré sections after integration of fluid particle trajectories on very long times, as we will see later. Therefore, because our 2D-model neglects the weak mixing inside the non-active pipes, it is likely to mix less than in reality.

Note that long time integration of trajectories is also more easily available with the 2D-model than with a 3D flow on a computational grid, since the probalility that a fluid particle ends on a boundary is very low in the former case: first the flow-field is known analytically, and therefore satisfies exactly the incompressibility constraint, and second the velocity, which is a solution of the Euler equation, is not zero on the boundaries.

Correspondence between parameters from different authors

It is clear that the efficiency of mixing does not depend on the value of the Reynolds number inside the pipes and inside the chamber, as long as it remains small. Therefore, for a given geometry and a given protocol, the important non dimensional parameter is the dimensionless pulse volume α,

α = qT hS , (1)
which represents the volume of fluid displaced during one period compared to the volume of the chamber. Although α is defined for the 3D flow, one can easily switch from volumes to surfaces by dividing by the -constant -height h of the chamber. Therefore, α can be used identically as its 2D equivalent, the dimensionless pulse area, equal to the surface of fluid displaced during one period compared to the surface of the chamber.

In a view of comparison with our previous works where we used the dimensional period of the flow T instead of α, note that those two parameters are linked by the relation α = 0.18 T (T in seconds) with our set of parameters.

While Stremler and Cola also used α, Hertzsch et al. used a non dimensional period -denoted here T * -defined as T * = qT /2πa 2 , where the distance between the source and the sink is equal to 2a; in their simulations, the radius of the circular chamber is equal to R, chosen as R = 2a, which leads to the relation α = T * /2.

Tools of chaos

Mixing efficiency is then analyzed with the usual tools of chaos: Poincaré sections and Lyapunov exponents. In a time-periodic flow, the Poincaré section of a given trajectory is obtained by plotting on the same graph all the positions of the initial point taken every period. A good mixing is obtained when the Poincaré section of a single trajectory covers the whole surface, in the most homogeneous way and the smallest period possible.

Mixing can also be analyzed quantitatively using the Lyapunov exponent λ, expressed in seconds -1 , which measures the asymptotic exponential stretching between two nearby trajectories (rate of separation). Classically in a 2D-incompressible flow, the sum of the two Lyapunov exponents is zero, so that we only need to calculate the largest Lyapunov exponent (in opposition to our previous works in 3D-flows where we determined the three Lyapunov exponents [START_REF] Ph | Lyapunov spectrum determination from the FEM simulation of a chaotic advecting flow[END_REF]). We check that the calculation is sufficiently converged by taking different initial points inside the stochastic region; the simulation is stopped when all Lyapunov exponents have converged towards the same value.

In the square and rectangular geometry, it is also possible to calculate the "filling rate" of the Poincaré section, i.e. the fraction of surface covered, so as to obtain a quantitative rather than visual information (for a discussion on how to measure the filling rate of an ensemble of points, see [START_REF] Jay | Lagrangian particle calculations of distributive mixing: Limitations and applications[END_REF]). We can therefore investigate the use of the topological or Kolmogorov entropy h K through the same estimation as Stremler and Cola [START_REF] Stremler | A maximum entropy approach to optimal mixing in a pulsed source-sink flow[END_REF], as the product of the fraction of the phase space covered by a single chaotic trajectory multiplied by the corresponding Lyapunov exponent (see [START_REF] Lichtenberg | Regular and Chaotic Dynamics[END_REF], section Definitions and basic concepts page 305). Mixing is supposed to be all the more efficient as the topological entropy is high.

Validation of the 2D model

It is important first to validate the 2D model properly. Indeed, our 2D model uses the depth-averaged horizontal velocity field, whereas the horizontal velocity profile is almost everywhere parabolic. This implies that a fluid particle located at middle-height moves faster than what we compute, while another one located for example near the bottom can be much slower. Therefore, we could fear some important discrepancies between our 2D results and our 3D calculations.

Note also that, since a Poincaré section accumulates the positions of a fluid particle after each period, and since a fluid particle can change depth, the 3D Poincaré section is a 3D set of points, and a regular region is likely to be hidden by points located above and below this region (corresponding to points with different velocities), which is not possible in the 2D case.

On figure 3 are compared 3D (left) and 2D (right) Poincaré sections, together for protocol A for three values of α, and for protocol B for a very small α -the comparison for protocol B and large values of α is of little interest, since points are almost randomly spread in the whole domain.

Concerning protocol B (figures 3(a) and (b)), the matching is very good. For protocol A α = 0.36 (figures 3 (c) and (d)), the situation is also very similar: the peripheral region of the domain, more dense, corresponds to a chaotic zone, while the region inside is a nearly regular region For protocol A, α = 0.72 (figures 3(e) and (f )) , the chaotic area (area full of points) and large regular island (empty area or area having "elliptical trajectories") have comparable shape, size, and location. The regular islands located around the large one are very thin, and therefore can hardly be detected in the 3D calculations, as explained before.

For larger values of α (for instance α = 1.8 on figures 3(g) and (h)), some very large regular regions are present in the 2D case, that are not visible in the 3D case. At first sight, the correspondence may appear quite bad. However, a careful analysis shows that the 2D Poincaré section is indeed informative in the present case: let us consider the large island which repeats identically at three occurrences in the section (bottom, right and top, denoted by regions A, B and C on figure 12). A point in one of these islands is in fact a periodic point, and enters and leaves the different pipes eight times before it comes back to its initial position seven periods later, see appendix B (we remind that a complete period is made of four stages). As explained before, some mixing occurs inside the pipe, even when non active. This explains why these zones are not visible in the 3D section.

The trapezoidal-shaped island on the left is, on the other hand, never swallowed inside a pipe, but moves in the domain with time (this is indeed the same regular region as the large island seen on figure 3(f ) and at the middle of figure 3(d), which grows in size and is slightly displaced towards the left with increasing period). Therefore, we should expect a fluid particle trapped inside this region to remain there for many periods before leaving; this is indeed the case, as can be seen on figure 3(g), where a curve corresponding to the trapezoidal area is visible. It is also noticeable that both this region and the central region are less visited than the rest of the domain (less points in the section), in accordance with the corresponding regular regions in the 2D section (figure 3(h)). Therefore, even in this case, we recover special features of the 2D section in the 3D calculations. Note that in a real experiment, because of molecular diffusion, the mixing will be even more rapid than in the 3D calculations.

Since we present Lyapunov exponents as a quantitative tool, we could wonder whether the hypothesis of a quasi 2D flow does not alter quantitatively this parameter. On figure 4 are presented with symbols the Lyapunov exponents obtained with our model flow versus the dimensionless pulse area α. As a view of comparison are plotted with lines the results obtained with the full 3D calculations in a square chamber [START_REF] Raynal | Towards better DNA chip hybridization using chaotic advection[END_REF]; the dotted line corresponds to protocol A and is to be compared with the empty square symbols (2D flow), while the solid line (protocol B) shall be compared with the full square symbols: somewhat surprisingly, the comparison shows almost perfect matching, except for very small values of α, when the horizontal dispersion by the flow is weak, so that 3D effects cannot be neglected anymore.

We therefore conclude that this 2D-velocity field model is adapted for a complete study of the problem.

Results

Poincaré sections and Lyapunov exponents

Poincaré sections are usually the most straightforward tools for characterizing chaos, because of the visual information they provide. In figures 5 and 6 are shown Poincaré sections of protocols A and B respectively, for three values of α, namely α = 0.18, α = 0.72 and α = 1.8, (corresponding to the periods T = 1 s, T = 4 s and T = 10 s from our previous work).

The circular, square or rectangular domains have been plotted with the same area so as to be able to compare them visually -in the case of the circle, although this is not obvious, we indeed have S = πR 2 = L 2 , with R the radius of the circle, which implies that the diameter is 2R/L = √ 4π 1.13.

As expected [START_REF] Raynal | Towards better DNA chip hybridization using chaotic advection[END_REF], we recover with our 2D-model that mixing is much more efficient with protocol B compared to protocol A in the square geometry; this result can obviously be extended to different geometries of the chamber.

A more surprising result is the influence of the geometry of the chamber: the fact is not straightforward concerning protocol B, since mixing is quite good, even if not perfect in the square geometry (some small regular regions still persist, figures 6 (e) and (h), but which are much too small to be seen in a 3D computation or an experiment). However, mixing is as good in the rectangle as in the circle. The results are much more striking concerning protocol A, since mixing is not very good, neither in the circle nor in the square: the very robust regular region for protocol A, α = 0.72, seen in the 3D or 2D calculations (figures 3(e) and (f ))

of the square geometry, also present in the circular geometry (figure 5 (d)), has totally disappeared in the rectangular geometry. For the largest value of α studied, α = 1.8 (also protocol A), the very large islands seen in the circle and the square (figures 5(g) and (h)) have a much more limited extent in the rectangular geometry. In fact, we performed many 2D calculations, but we could hardly find cases with regular regions (except for protocol A and large values of α, where a fluid particle at the center is always swallowed into a sink after a source/sink phase, and always comes back to its initial position when the sink is turned into a source again)although we will show two particular cases (one for each protocol) where some small islands are indeed still present-.

We wanted to know whether the fact that mixing was best in the rectangular geometry was robust with respect to the aspect ratio of the rectangle, or else if the 2:1 aspect ratio was a particular case: on figure 7 are shown different Poincaré sections, for the same case protocol A, α = 0.72, from the square geometry to rectangles with aspect ratii ranging from 3:2 to 4:1. In all rectangular Poincaré sections, the regular island is not present (although all of them present the same star-shaped region where the repartition of points is not completely homogeneous), as if the topology of chaos did depend much on the geometry (circular, square or rectangular) but not on the aspect ratio of the rectangle either. A similar conclusion can be drawn when looking at figure 8, where the small regular islands seen both for protocol A with α = 1.44 or protocol B with α = 1.8 are topologically similar in the 2:1 and 3:1 rectangles.

We could calculate the Lyapunov exponent for both protocols and each geometry (circle, square, rectangle 2:1). The results are shown on figure 9, where the black symbols represent protocol B, while the empty symbols are for protocol A: surprisingly, although the Lyapunov exponent depends much on the protocol, it seems not to depend much on the geometry -it is slightly above in the rectangular case for α < 1.44, but is also slightly below for α > 1.44, although the Poincaré sections always show a much better chaos-. Moreover, the maximum of the Lyapunov exponent (the smallest value of α tested, α 0.045 for protocol B, and α 0.09 for protocol A) is not representative of a good mixing: the stirring is important, but in a restricted chaotic domain. We can conclude that the Lyapunov exponent is quite robust with respect to the protocol; however, as expected the Lyapunov exponent is not the relevant quantity to decide between different geometries.

Topological entropy

As mentioned at the end of the previous section, there is a need to find a different quantitative parameter which could help decide on whether the rectangular geometry indeed leads to more efficient mixing than square or circular shapes. We propose to try to use the topological entropy, calculated as the fraction of surface covered by a single trajectory, multiplied by the corresponding Lyapunov exponent [START_REF] Lichtenberg | Regular and Chaotic Dynamics[END_REF]. As seen before, the Lyapunov exponent depends poorly on the geometry of the chamber, so that the topological entropies for different geometries should follow more or less the behavior of the filling rate of the sections. However, for very small values of α the filling rate is low, even for protocol B, while for the same protocol the Lyapunov exponent is a decreasing function of α (figure 9); therefore we could wonder whether the product of the filling rate by the Lyapunov exponent would present a maximum that could be interpreted as an indicator of the optimum of α.

On figure 10(a) is shown the filling rate of the Poincaré section as a function of the dimensionless pulse area α, for both protocols in the square and rectangle 2:1, calculated using a single fluid particle followed for 10.000 periods. In accordance with the aspect of the Poincaré section, the filling rate is nearly always higher in the rectangular geometry. When multiplied by the Lyapunov exponent, the topological entropy is obtained, as shown in figure 10(b): concerning protocol B, no maximum was found, neither in the square nor in the rectangle; as for the Lyapunov exponent, the quantity is decreasing with increasing α, greater for the rectangle than for the square for α < 1.44, and lower thereafter. It is clear though that the smallest values of α are not those associated with the best chaos, as seen on figures 6(a), (b) and (c) for protocol B, α = 0.18. In the case of protocol A however, a maximum exists for the topological entropy, located at α 0.36 for the square (which corresponding Poincaré section can be seen on figure 3(d)), and α 0.72 for the rectangle (figure 5(f )). This may indeed correspond to an optimal period for mixing in the rectangle, but inspection of figure 3(d) reveals that this is not obvious in the case of the square. Moreover, the topological entropy does not differ much from square to rectangle, although the filling rate proves that chaos is much more extended in the rectangle for this protocol. Therefore, the use of the topological entropy, even if giving some interesting quantitative information, is not totally convincing. . . .

Discussion

The protocols have approximately the same symmetries for the different geometries: in the case of protocol A the equations of the trajectories in the 2D or 3D flow are invariant under the change of variables t -→ -t, x -→ x, y -→ -y, z -→ z. This implies that the Poincaré sections P t=0 , corresponding to the accumulated periodic positions of particles, starting at t = 0, are invariant under reflection in the plane y = 0. In the square, a translation in time t -→ t + T /4 leaves the flow unchanged after rotation of angle -π/2, and the Poincaré section P t=T /4 is obtained from the Poincaré section P t=0 by a rotation of angle -π/2; this is of course not true anymore in the rectangle, but if now we consider the translation in time t -→ t+T /2, two consecutive rotations of angle -π/2 are equivalent to a rotation of angle π, which is also a point symmetry. One can check that this last symmetry is also present in the rectangular geometry. In the case of protocol B, the invariance of particle trajectories under t -→ -t, x -→ -x, y -→ y, z -→ z implies that the Poincaré sections P t=0 are invariant with respect to reflection in the plane x = 0. The translation in time t -→ t + T /2 leaves the flow unchanged after reflexion in the plane y = 0, which implies that the Poincaré section P t=T /2 is, in all the geometries studied here, obtained from the Poincaré section P t=0 by reflection in y. Therefore, the improvement of results seen in the Poincaré sections in the rectangle does not result from a breaking of symmetry in the combined effect of protocol and new geometry.

Recently, Hertzsch et al. showed that flows generated by source-sink pairs could be studied as "Linked Twist Map" (LTM) [START_REF] Hertzsch | DNA microarrays : Design principles for maximizing chaotic, ergodic mixing[END_REF]. The central idea of their article was to construct crossed periodic flows having properties as close as possible to those of an ideal model, i.e., (I) crossing of streamlines as orthogonal as possible, and (II) velocity field as close to monotonous as possible. We propose to test properties (I) and (II) on our different geometries.

On figure 11 are shown superimposed streamlines for the two source-sink pairs for the different geometries. In all cases the streamlines are close to tangential near the boundaries. In the central region however, the greatest transversality is found in the circular and square geometry, where streamlines are close to orthogonality, whereas mixing is much more efficient in the rectangle. In the case of the rectangle 3:1 (figure 11(d)), the situation is even worse, since the streamlines in the central region are nearly tangential. Loss of transversality works here in favor of mixing, rather than against it.

We now try to find out whether the flow-field in the rectangular geometry is more monotonous than in the square or the circle. It is quite difficult to compare velocity profiles since the flow is not unidirectional, but we can compare the time needed to go from the source to the sink on each streamline instead. To this end, let us come back to figure 2 where streamlines of a source-sink pair are shown for the different geometries: for a given phase with a source-sink pair, let us consider the segment joining the two other corners, that we will chose as our axis of reference. We graduate it from -c to +c, where c is the distance from the center to one corner, and we denote by τ (w) the time needed to join the sink starting from the source (-c ≤ w ≤ c).

In the circle and the square, the streamlines are symmetric with respect to any line joining a source and a sink; therefore τ (w) is a symmetric function, and is absolutely not monotonic. In the off diameter source-sink model proposed by Hertzsch et al. in a circle [START_REF] Hertzsch | DNA microarrays : Design principles for maximizing chaotic, ergodic mixing[END_REF], τ (w) is dissymetric, and much closer to a monotonic function. In the rectangle however, although the two symmetries with respect to the lines joining a source-sink pair do not hold, the streamlines are clearly symmetric with respect to the center. This implies that τ (w) is, as for the circle or the square, a symmetric function, therefore not monotonic either... We conclude that the only reason for the improvement of mixing in the rectangular geometry is the loss of symmetry of the streamlines themselves. The line joining the source and the sink is no more a streamline of the corresponding flow. Instead, the streamline that goes through the center point wanders from one side of this line to the other, all the more as the rectangle has a large aspect ratio. This last point could explain also why, although superimposed streamlines are even more tangential when the aspect ratio of the rectangle is increased, mixing remains good, and the Poincaré section mostly unchanged.

Summary and conclusion

In this article we tested the influence of the geometry of Hele-Shaw cells on the quality of chaotic mixing, using two different mixing protocols that relied on chaotic advection. In a view of fast design, we used a 2D model consisting in the depth-average velocity field: it could be easily calculated analytically for a source-sink flow, first in a circle, and in a square or rectangle thereafter with a Schwarz-Christoffel transformation. The model was validated using various Poincaré sections in a square, that qualitatively matched our previous 3D calculations. As a more quantitative tool, 2D and 3D Lyapunov exponents were also compared, with good agreement. We thus could show that mixing in the rectangular geometry was nearly always best (better filling of the Poincaré sections, usually higher Lyapunov exponents, higher topological entropy), whatever the protocol used. The result did not depend much on the aspect ratio of the rectangle.

We tried to find out whether the improvement in mixing in the rectangle could be relied to loss of symmetry in the protocols compared to squared or circular geometry, with no success. Then we searched for properties proposed by Hertzsch et al. in order to approach more ergodic systems, namely increase in transversality of superimposed streamlines, and monotonicity of the velocity profile, but this study revealed that the profile was no more monotonic in the rectangle than in the square or the circle, and that, on the contrary, the streamlines were even less perpendicular. We then concluded that the only reason for this improvement was the loss of symmetry in the streamlines for a source-sink flow in the rectangular geometry: whereas in the square or the circle the streamlines have two axes of symmetry, in the rectangle only a point symmetry remains. More specifically, the straight line joining the source and the sink is no more a streamline in the rectangle.

This improvement in mixing in the rectangular geometry is all the more important in practice as a straightforward application of mixing in Hele-Shaw cells is the enhancement of reliability of DNA chips technology, and that most DNA chips are manufactured on microscope lids, which have a rectangular shape.

A Velocity field

In order to obtain the flow in the squared or rectangular domain, the velocity field is first easily calculated in a circular chamber using the Milne-Thomson theorem. Let v x (x, y) and v y (x, y) the components of the velocity-field at point (x, y). The complex velocity-field w c (z) is

w c (z) = v x (x, y) -i v y (x, y) (2)

Figure 1 :

 1 Figure 1: (a) : Alternating injection scheme for protocol A. Each quarterperiod step of the protocol involves two opposite syringes (the black ones in step 1 and 3, the grey ones in steps 2 and the two other ones being inactive. (b) : Same for protocol B. Each half-period step of the protocol involves one of the pumps, the other being inactive. One pump always pushes the fluid in the same direction, or else is inactive. Note : the DNA chip is symbolized by the central black square.

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 234567891011 Figure 2: Streamlines for a source-sink flow in the different geometries of the chamber.

Figure 12 :

 12 Figure 12: Sketch of the velocity field for the Poincaré section corresponding to protocol A, α = 1.8 in the square: during stage 1 (respectively 2, 3 or 4), the sink is the hole located at cross I (respectively II, III or IV ), the source is the hole located at the cross at the opposite corner.

where i is the natural base of complex numbers, defined as i 2 = -1, and z = x + i y. In a circular cavity of radius R, the flow of rate q created by one pair source/sink in opposition at length βR from the center of the cavity is:

Then a Schwarz-Christoffel function g, adapted to the transformation of a circle into a square or rectangle, is calculated [START_REF] Milne-Thomson | Theoretical Aerodynamics[END_REF][START_REF] Milne-Thomson | Theoretical Hydrodynamics[END_REF]: let k and ξ be solutions of equation:

where K is a Jacobi function [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] and a and b the lengths of sides of the squared (or rectangular) cavity. Then g is defined as:

where sn, cn and dn are elliptical functions [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF].

Finally, the velocity field w r (z) inside a square (or rectangular) shape is defined as:

B Study of the period-7 periodic point in figure 3(d)

We describe here briefly how the regular region that wanders in the whole squared domain in figure 3(d) returns to its initial position after 7 periods, while only three identical regions are seen on the Poincaré section. Consider figure 12: crosses I, II, III and IV represent the pipes. Protocol A consists in four stages per period; during the first one, denoted by stage 1 (resp. 2, 3 and 4), the active couple is I and its opposite hole, the sink being I (resp. II, III and IV ). For this set of parameter, a fluid particle located in regular region A (respectively D) is mapped into B (respectively C) during this stage, and a fluid particle located at B or C is swallowed into pipe I. We remind that a fluid particle swallowed into a pipe comes back at the same location when the sink is turned into a source again. We can therefore follow a fluid particle initially in regular region A, knowing that each period is composed of the four stages 1, 2, 3 and 4.


Luciano Da
email: luciano@ifsc.usp.br

Luciano Da Fontoura Costa

Autorrelation and Cross-Relation of Graphs and Networks

research documents, whether they are published or not. The documents may come

Autorrelation and Cross-Relation of Graphs and Networks

Introduction

Thanks to continuing advances from graph theory (e.g. [START_REF] Da Fontoura Costa | Introduction to graph theory[END_REF][START_REF] Bollobás | Modern graph theory[END_REF]) and network science (e.g. [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][START_REF] Da | Analyzing and modeling real-world phenomena with complex networks: a survey of applications[END_REF]), graphs and complex networks now underly a wide range of concepts and applications, extending from systems biology to transportation systems. In this work, we will understand complex networks, or networks for short, as corresponding to particularly intricate types of graphs (e.g. [START_REF] Da | What is a complex network?[END_REF][START_REF] Costa | Cost-based approach to complexity: A common denominator? Revista Brasileira de Ensino de Física[END_REF]).

One of the most interesting characteristics of graphs and complex networks is their ability to represent virtually every discrete system and phenomenon. Indeed, from the point of view of data structures (e.g. [START_REF] Knuth | The Art of Computing[END_REF][START_REF] Wirth | Algorithms & data structures[END_REF][START_REF] Cormen | Introduction to Algorithms[END_REF]), graphs can be understood as providing the most general possible manner to represent other structures, including lists, trees, lattices, etc.

The intrinsic representational generality of graphs and networks has motivated a directly related interesting issue, namely how to mathematically handle, combine and compare the represented structures. This bears a direct analogy with the manipulation of mathematical structures such as functions, fields, vectors and matrices by using arithmetic and algebraic concepts and methods including addition, product, correlation, etc.

Among the several applications operations between mathematical structures underlain by respective data structures (e.g. vectors and matrices), the operations of self-and cross-correlation, as well as the closely related self-and cross-convolution counterparts (e.g. [START_REF] Brigham | Fast Fourier Transform and its Applications[END_REF][START_REF] Da | Shape Classification and Analysis: Theory and Practice[END_REF][START_REF] Vetterli | Foundations of signal processing[END_REF][START_REF] Haykin | Neural Networks And Learning Machines[END_REF][START_REF] Gallager | Stochastic Processes: Theory for Applications[END_REF]), play a key role in several important areas, including signal and image processing and analysis, mathematical physics, computer vision, control theory, and electronic engineering, to name but a few possibilities. For instance, the blurring of an image is frequently obtained by convolving the image with a gaussian kernel (e.g. [START_REF] Da | Shape Classification and Analysis: Theory and Practice[END_REF][START_REF] Gonzalez | Digital Image Processing[END_REF]), and the recognition of an instance of a pattern within a signal can be approached in terms of the respective crosscorrelation (i.e. template matching, e.g. [START_REF] Da | Shape Classification and Analysis: Theory and Practice[END_REF][START_REF] Horn | Robot Vision[END_REF][START_REF] Davies | Machine Vision[END_REF]).

While self-and cross-correlations are well defined respectively to vectors and matrices, their extension to graphs and networks is not so straightforward. Related approaches include the concept of graph Fourier transform, which involves the spectral analysis (eigenvalues and eigenvectors) of the Laplacian matrix of a given graph or network (e.g. [START_REF] Shuman | A windowed graph fourier transform[END_REF][START_REF] Sardellitti | On the graph fourier transform for directed graphs[END_REF][START_REF] Ricaud | Fourier could be a data scientist: From graph fourier transform to signal processing on graphs[END_REF]). The term Fourier transform here is used as an analogy, not literally. Another approach related to the relationship between portions of a graph is correlation clustering, aimed at finding cuts that optimize the allocation of edges with large weights to result within respective clusters (e.g. [START_REF] Demaine | Correlation clustering in general weighted graphs[END_REF][START_REF] Hero | Hub discovery in partial correlation graphs[END_REF][START_REF] Bansal | Correlation clustering[END_REF]). Statistical-based approaches to correlation between graphs have also been suggested (e.g. [START_REF] Berg | Correlated random networks[END_REF][START_REF] Fujita | Correlation between graphs with an application to brain network analysis[END_REF]).

The present work sets out at investigating the possibility to use multiset similarities, in particular the coincidence similarity index [START_REF] Da | Further generalizations of the Jaccard index[END_REF][START_REF] Da | On similarity[END_REF][START_REF] Da | Coincidence complex networks[END_REF], as a means for obtaining analogues of the self-and cross-correlation operations respectively to graphs and networks. The main justification for this approach is the fact that the coincidence similarity provides a selective, sensitive, and stable/robust means for comparing any mathematical structures (e.g. [START_REF] Da | On similarity[END_REF][START_REF] Da | Comparing cross correlationbased similarities[END_REF][START_REF] Da | Multiset neurons[END_REF][START_REF] Da | Coincidence complex networks[END_REF]).

In particular, the coincidence operation is here employed in direct analogy to the inner product between two signals as the means for gauging the similarity between two graphs or networks. Then, in order to implement the operation analogous to the relative shifting of one signal respectively to the other, we consider the mean coincidence similarity between each node of the networks and the nodes belonging to successively respective neighborhoods, or hierarchies (e.g. [START_REF] Travençolo | Hierarchical spatial organization of geographical networks[END_REF][START_REF] Costa | Hierarchical characterization of complex networks[END_REF][START_REF] Costa | What are the best concentric descriptors for complex networks[END_REF]), at topological distance δ from the reference node. Thus, in this work the variable d plays an analogous role to the lag variable in traditional correlation analysis.

The resulting operations between two graphs are here called autorrelation and corss-relation so as to indicate their analogy with the classic autocorrelation and cross-correlation, but at the same time distinguishing the counterparts because they are based on different concepts: while the latter two are based on the inner product between functions or vectors, the autorrelation and corss-relation are founded on the coincidence similarity between features representing the topological (or other types) of characteristics of the considered network nodes.

In addition to presenting the thus motivated self-and cross-relation between two graphs or network (a graph with itself in the case of autorrelation), we illustrate the potential of the suggested concepts and methods respectively to the identification of lag-based relationships between model-theoretic and real-world networks. In addition, we also address the possibility to derive coincidence networks from the individual autorrelation signatures, with remarkable results that can reflect in a particularly effective manner the intricacy of the topology around the nodes of the original networks.

The obtained results indicate that the autorrelation and cross-relation between networks reflect several global and local properties of the analyzed networks, therefore constituting an interesting resource for network characterization.

We start this work by presenting the several involved concepts, presenting the suggested autorrelation and cross-relation operations, and illustrating their potential for the characterization of global and local properties of several model-theoretical and real-world networks.

Basic Concepts

Given two vectors x and y with N elements, their inner product can be expressed as:

x, y = N i=1 x[i][i] (1)
If the magnitudes | x| and | y| of the two vectors are kept fixed, the respective inner product will provide a measurement of the similarity between these two vectors.

The cross-correlation between two vectors (or signals) x and y with sizes N can be written as:

( x ⊗ y)[ ] = N -1 i=1 x[i] y[i -] (2)
where vector y is suitably padded with zeros and = 0, 1, . . . , N -1, so that the result has size N -1.

The convolution between two vectors (or signals) x and y with sizes N can be written as:

( x * y)[ ] = N -1 i=1 x[i] y[ -i] (3)
where vector y is suitably padded with zeros and = 0, 1, . . . , N -1, so that the result has size N -1.

Let g(t) be a function, with respective Fourier transform G(f ). Its autocorrelation (e.g. [START_REF] Brigham | Fast Fourier Transform and its Applications[END_REF][START_REF] Da | Shape Classification and Analysis: Theory and Practice[END_REF]) can be calculated as:

g(t) ⊗ g(t) ⇐⇒ G(f ) * G(f ) = |G(f )| 2 (4)
where |G(f )| 2 is often called the power spectrum of f (t). Given a undirected graph g, the degree of each of its nodes corresponds to the respective number of links. The degree of node, in general, possibly corresponds to the most important topological information about the local topology around that node. The clustering coefficient or transitivity of a node corresponds to the ratio between the total number of links existing between the neighbors of that node and the total maximum number of links that would be possible (e.g. [START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF][START_REF] Barabási | Network Sience[END_REF][START_REF] Newman | Networks: An introduction[END_REF]) Given a graph (or network) g, and one of its nodes i, the N δ (i) neighborhood of i corresponds to the set of all nodes that are at topological distance smaller or equal to δ. These successive neighborhoods can be understood as defining a respective hierarchy respectively to i.

Given two networks g and h with N nodes each, we say they are aligned whenever for each node i of g there corresponds an associated node i in h. Observe that the two networks can have distinct interconnections.

The multiset coincidence index has been suggested as a manner to generalize the Jaccard index (e.g. [START_REF] Jaccard | Distribution de la flore alpine dans le bassin des dranses et dans quelques régions voisines[END_REF][START_REF] Jaccard | Étude comparative de la distribution florale dans une portion des alpes et des jura[END_REF][START_REF] Samanthula | Secure multiset intersection cardinality and its application to jaccard coefficient[END_REF][START_REF]Jaccard index[END_REF][START_REF] Schubert | A note on the Jaccardized Czekanowski similarity index[END_REF]) in order to take into account the interiority between the two compared mathematical structures [START_REF] Da | Further generalizations of the Jaccard index[END_REF], being based on multiset theory (e.g. [START_REF] Hein | Discrete Mathematics[END_REF][START_REF] Knuth | The Art of Computing[END_REF][START_REF] Blizard | Multiset theory[END_REF][START_REF] Blizard | The development of multiset theory[END_REF][START_REF] Mahalakshmi | Properties of multisets[END_REF][START_REF] Singh | Complementation in multiset theory[END_REF]). The real-valued coincidence has been found to provide more strict comparison between two multisets than the Jaccard or overlap indices taken separately, which has paved the way to enhanced selectivity, sensitivity and robustness (e.g. [START_REF] Da | On similarity[END_REF][START_REF] Da | Comparing cross correlationbased similarities[END_REF][START_REF] Da | Multiset neurons[END_REF]). More recently, the average of the coincidence similarity established between each node of the network and all the others was suggested [START_REF] Da | On the self-coincidence structure of networks[END_REF] as an feature that can be used to effectively complement to the characterization of the topology around each node.

It is henceforth assumed that all network measurements are non-negative. Otherwise, the real-valued version of the coincidence index (e.g. [START_REF] Da | Further generalizations of the Jaccard index[END_REF][START_REF] Da | On similarity[END_REF]) should be considered.

If x and y are two vectors taking non-negative real values with size N , their multiset Jaccard similarity (e.g [START_REF] Da | Further generalizations of the Jaccard index[END_REF][START_REF] Da | On similarity[END_REF][START_REF] Da | Coincidence complex networks[END_REF]) can be defined as:

J ( x, y) = N i=1 min {x[i], y[i]} N i=1 max {x[i], y[i]} (5)
Observe that 0J ( x, y) ≤ 1. Their respective multiset interiority index (also overlap, e.g. [START_REF] Vijaymeena | A survey on similarity measures in text mining[END_REF]) can be expressed as;

I( x, y) = N i=1 min {x[i], y[i]} min {S x , S y } (6)
with 0 ≤ I( x, y) ≤ 1.

where:

S x = N i=1 x[i]; S y = N i=1 y[i] (7)
The multiset coincidence similarity between the vectors x and y corresponds to the product of the respective realvalued Jaccard and interiority indices, i.e.: Their respective real-valued interiority index (also overlap, e.g. [START_REF] Vijaymeena | A survey on similarity measures in text mining[END_REF]) can be expressed as;

C( x, y) = J ( x, y) I( x, y) (8)
with 0 ≤ C( x, y) ≤ 1.

The real-valued coincidence similarity between the two vectors x and y can be understood as being analogous the Pearson correlation coefficient between the two vectors.

Autorrelation and Cross-Relation between Two Graphs

Having revised the main basic concepts, we are now in position to mathematically formalize the concepts Autorrelation and Cross-Relation between two graphs or networks.

Given two aligned networks g and h with N nodes, as well as a set of M topological features x j , j = 1, 2, . . . , M , we can associate each node i = 1, 2, . . . , N of those networks to the respective features, therefore obtaining the following feature vectors:

X g,i = [x g,1 [i] x g,2 [i] . . . x g,M [i] ] T (9) X h,i = [x h,1 [i] x h,2 [i] . . . x h,M [i] ]
T [START_REF] Wirth | Algorithms & data structures[END_REF] where x g,j [i] means the value of the feature x j of node i in graph g.

The respective cross-relation of those two aligned vectors can now be defined as:

(g h)[i, δ] = 1 |N δ (i)| k∈N δ (i) C( X g,i , X h,k ) (11)
where δ = 1, 2, . . . , ∆ and |N δ (i)| is the number of neighbors of node i at topological distance δ.

In case g = h, the above expression becomes the autorrelation of a node i of the network g:

(g g)[i, δ] = 1 |N δ (i)| k∈N δ (i) C( X g,i , X g,k ) (12)
Observe that (g h)[i, δ] with δ = 1, 2, . . . , ∆ can be understood as a kind of signature of each of the network nodes i = 1, 2, . . . , N , corresponding to how much that node is similar in the average, as far as the coincidence index is concerned, to the nodes at each successive neighborhood N δ (i). The basic underlying assumption is that the system of progressive neighborhoods around each node i can establish a coordinate axis associated to that node, and with origin at that node when δ = 0. Therefore, N coordinate axes are established in this manner, allowing one of the networks to be shifted respectively to the other in this highly dimensional space.

As such, the autorrelation of a network g can be understood in analogy to the autocorrelation of a function or signal, while a similar relationship can be established between the cross-relation between two networks g and h and the cross-correlation between two signals or functions.

As defined above the autorrelation and cross-relation between two graphs will involve the coincidence comparison between the features of each node with those of successive neighborhoods, which tend to increase steadily as each node in the neighborhoods radiates more links, and then decrease as a consequence of reaching the border of the network. This implies that the average coincidence is taken respectively to varying numbers of nodes in each successive neighborhood. There are several ways in which the autorrelation and cross-relation between two graphs can be instantiated that can be used to keep this number more commensurate or even constant.

One such manner consists in comparing between each node in one network and a walk along the nodes in the other network. These walks can be specified in a vast number of ways, including the more traditional random walks and self-avoiding random walks (e.g. [START_REF] Lovász | Random walks on graphs[END_REF][START_REF] Noh | Random walks on complex networks[END_REF]).

The autorrelation (cross-relation) presents potential for characterizing how much one (two) graphs are similar respectively to a sequence of respective topological lags δ. In the remainder of this work, we illustrate the application of the autorrelation and cross-relation concepts respectively to several model-theoretical and real-world sections.

An important aspect of the autorrelation and crossrelation between graphs as suggested above concerns the fact that these operations are relative to the set of adopted node features, which can be of topological or complementary type. This flexibility to define specific set of features allow the lag-based comparison of networks to be investigated from the point of view of several choices of measurements, aimed at reflecting network properties of specific interest in each research program. For simplicity's sake, the measurements of node degree and clustering coefficient will be adopted henceforth.

Results and Discussion

In this section, the concepts of autorrelation and crossrelation between two graphs are illustrates respectively to model-theoretical and real-world examples. In all cases, networks that happen to be directed are first symmetrized, and only the largest connected component is taken into account. In addition, given that maximum autorrelation is always observed for δ = 0, these values are not shown for simplicity's sake. Observe that this property does not hold for the cross-relation.

We start the study of the application of the autorrelation and cross-relation concepts respective to three modeltheoretical complex networks (e.g. [START_REF] Da | Characterization of complex networks: A survey of measurements[END_REF][START_REF] Barabási | Network Sience[END_REF][START_REF] Newman | Networks: An introduction[END_REF]), namely: Erdős-Rényi -ER, Barabási-Albert -BA, and Watts-Strogatz -WS in its ring configuration. Consequently, we have the possibility to consider uniformly random interconnectivity (ER), scale free degree distribution (BA), as well as the model WS which varies, in terms of a respective rewiring probability, from a perfect lattice to an ER configuration. The ER, BA and WS with high reconnection probability are all characterized by having relatively small average shortest path distances.

Figure 1(a) presents the autorrelation obtained for an ER network with N = 300 nodes and average degree k = 4, considering topological lags δ = 1, 2, . . . , 10. Shown are the several autorrelation signatures obtained for each of the 300 nodes (several random colors), as well as the respective mean signature (shown in black). Inter-estingly, the autorrelation function of this network can be found to decrease monotonically with two negative slopes: one smaller taking place from δ equal to 1 to 4, and another steeper slope thereon. The autorrelation function for this network can then be subdivided into two portions, the latter starting when the neighborhoods reach the border of the network. Indeed, the lower curves in this figures correspond to some of the border nodes, implying individual autorrelation that is initially small (the border elements tend not to be similar to the inner nodes) and then increases as the neighborhoods reach the core of the network.

The autorrelation function of a BA network also with N = 300 nodes and k = 4 is presented in Figure 1(b). Compared to the ER case, a completely different result has been obtained. Now, the autorrelation increases progressively, and then undergoes a steep decrease. Similarity to the ER, though, the latter regime again coincides with the neighborhoods hitting the network border. The autorrelation increase, however, is a phenomenon specific to this type of network. It follows from the fact that, in a BA network, although only a few nodes are hubs, they are soon reached by the neighborhood progressions starting from the other nodes. Because the hubs are not similar to the other nodes, while the other nodes are mutually similar, the autorrelation tends to increase emphafter hitting the hubs, until the border is reached. The lower individual signatures are associated to the hubs of the network, which have markedly distinct topological features respectively to the other nodes.

Figure 1(c) depicts the autorrelation signatures obtained for a WS network (ring) with similar parameters as before, obtained by using a rewiring probability of p w = 0.05. A third type of autorrelation function can be observed for this type of network, characterized by an initial plateau, followed by the a steep decrease which again is a consequence of the finite size of the network. The obtained plateau derives from the intrinsic mutual similarity between most nodes of this network. Interestingly, the nodes corresponding to the shortcuts implemented by the rewirings can be readily observed as corresponding to the lower signatures in this figure .  In order to explore in more detail the autorrelation structure of the WS networks, Figure 2 presents the results respectively obtained for WS networks with N = 300 nodes and rewiring probabilities p w = 0.01, 0.02, and 0.03.

It can be readily observed that the increase of the rewiring probability not only reduces the range of lags characterized by non-zero autorrelations, but also influences the sharpness of the transition from core to border, as well as the height and extent of the plateaux. Also of interest are the lower signatures, which correspond to the nodes corresponding to the network rewirings.

A noticeably interesting possibility is to illustrate the autorrelations between the nodes in a given network in terms of a respective coincidence network [START_REF] Da | Coincidence complex networks[END_REF]. Examples of this type of coincidence network -obtained for the above ER, BA and WS networks, is shown in Figure 3 (a), (b), and (c), respectively. The size of the nodes reflect the degree of the nodes in the original network.

As could be expected, a relatively uniform autorrelation network has been obtained for the ER network, other than by statistical fluctuations. In addition, most of the nodes are shown in cyan, reflecting the small transitivity of the original respective reference nodes.

A more interesting structure can be observed in the case of the BA structure, with the hubs either resulting in one of the extremities of the larger connected component, which are followed by nodes of signatures obtained for medium and small degree original nodes. This type of structure reflects the fact that, in a BA network, the hubs are different from the other nodes regarding their topological properties along the successive neighborhoods, but are mutually similar. As could be expected, most original nodes have very small transitivity (most nodes are in cyan).

A strongly uniform autorrelation network has been obtained for the WS case, reflecting the predominant regularity of the original network. Interestingly, some nodes in cyan can be seen around the uniform mass of node, most of which with significant transitivity values (in magenta). The former nodes correspond to the nodes the nodes establishing shortcuts in the original network.

All in all, the autorrelation signatures obtained for the three model-theoretical networks confirm the ability of this function to reflect several of the critical properties of the respective networks, as well as their parameter configurations. The individual signatures also provide means for identifying special types of nodes.

We now proceed to verifying the types of autorrelation functions as obtained from some real-world networks.

Figure 4(a) presents the autorrelation signature obtained for the Zachary karate club network [START_REF] Zachary | An information flow model for conflict and fission in small groups[END_REF][START_REF] Csardi | The igraph software package for complex network research[END_REF]. Despite the relatively small size of this network, the respective autorrelation signatures suggest that it has hubs, possibly possessing a scale-free node distribution. The obtained autorrelation signatures indicate similarity with the ER model, being characterized by monotonically decreasing autorrelation values.

The autorrelation signatures obtained for the US airports network [START_REF] Csardi | The igraph software package for complex network research[END_REF] is presented in Figure 4(d). The average signature, which involves small values, indicates a relatively more complex autorrelation structure for this network, involving the three following regimes: (a) an initial plateaux that indicates mutual similarity between each node and those at the nearest respective neighborhoods; (b) a following peak of autorrelation values, possibly related to presence of hubs in this network; and (c) a moderate decrease implied by the neighborhood progressions hitting the border of this network. An impressive variety of individual autorrelation signatures can be identified, some characterized by particularly small autorrelation values, indicating the presence of several interesting outlier nodes.

Figure 5 depicts the coincidence network obtained by considering the individual autorrelation signatures obtained for the US airports network as input to the coincidence methodology for transforming datasets into respective networks [START_REF] Da | Coincidence complex networks[END_REF]. The resulting coincidence network has been thresholded at T = 0.87 for the sake of improved visualization.

Several interesting aspects can be observed from the obtained coincidence network. First, we have that the signatures have been organized along a sequence of interconnections related to the degree of the original nodes. The hubs, being most different from the other nodes, but tending to be similar each other, resulted at one of the extremities of the main obtained structure. As the degree of the nodes becomes smaller, the observed structure initially undergoes a kind of mostly sequential unfolding, up to a point where several branching and more intricate structures starts to emerge.

Of particular interest is the separation of this initially string-like structure into two major branches, each of which being presenting complex interconnections, including ring-like substructures. At the same time, several smaller disconnected components involving strongly interconnected nodes with smaller degree degree can be also observed. Also of particular interest are the small string component observed at the top left-hand side of the visualization in Figure 5, as well as the nodes that remained disconnected. These nodes, the vast majority of which being characterized by low transitivity, can therefore be understood as outliers in this network as far as their hierarchical similarity structure is concerned respectively to those of the other network nodes.

It is also interesting to compare the coincidence network obtained for the US airports individual autorrelation signatures with that of a BA network, as illustrated in Figure 3(b). Though these two networks share some important properties, such as the sequential connection in terms of the degree of the original nodes, they also present important differences, including the markedly more intricate structure of the US airports network for smaller values of original node degree, possibly as a consequence of the more local and geographical nature of the interconnections between smaller airports, thus suggesting that this network involves a mixture of structural organization principles.

In order to verify the effect of the extension of the considered range of lags taken into account while building the coincidence network from the individual autorrelation signatures, we shown in Figure 6 this type of network derived from the US airports individual signatures considering successive lags. The effect of taking into account longer lag extensions contributes to obtaining progressively more detailed coincidence networks. Interestingly, the complex topological structure of this network could not be discerned from the coincidence networks considering narrow lag extensions. At the same time, the effect of incorporating additional lags tends to be larger for smaller lag extensions, with progressively more similar coincidence networks being obtained as the maximum lag approaches the borders of the network. Interestingly, the enhanced selectivity and sensitivity of the coincidence similarity index provides allows the intricate level of details of the autorrelation networks to be effectively revealed.

The above example illustrates the potential of the autorrelation approach for revealing impressively detailed information about the network being analyzed from the perspective of its hierarchical similarity structure.

Examples of cross-relations between two graphs are de- picted in Figure 7, referring to the cross-relation between original and rewired versions of a same network. First, we have the ER network in Figure 1(a) before and after 500 rewirings. Interestingly, a noticeable effect on the average signature can be observed, which is a consequence of random fluctuations given the relatively small size of this network. The cross-relation between networks and their rewired versions is strongly affected by the necessarily implied reduction of the average shortest path length, which substantially narrows the extent of the signatures along the lag δ. Another interesting study of cross-relations between two graphs consists of considering the original graph and a respective version with a given percentage of hubs removed from the network. Figure 8 illustrates this type of analysis respectively to the removal of 5% of the nodes with the highest degrees respectively to the ER, BA and US airports networks.

As could be expected, this modification had almost no effect in the case of the ER network, which does not tend to have nodes with substantially higher degrees (hubs) except for statistical fluctuations. However, in the case of the BA network, several alterations can be verified as a consequence of the removal of only 5 hubs. We can observe a reduction of the extent of the signatures, as well as a substantial shape of the average signature, which nows corresponds mostly to a plateau instead of a peak. The effects on the US airports network are also strong, implying longer signatures as well as a reduction of the peak and change of the shape of most signatures. This result indicates that the hubs in this network have a strong role in defining its respective topological properties.

Figure 9 presents the cross-relation analysis of networks under another type of modification, now involving the complete permutation of the labels of all nodes respec-tively to the original network. The effect of this alteration was similar in the case of the ER, BA, and US airports networks, consisting in making the signatures more similar, therefore revealing the fact that the respective networks were perceived as being substantially more regular when compared under nodes permutation.

The above cross-relation examples above illustrate the potential of this approach for comparing between two networks while considering respective displacements along the distance lags. The obtained results indicate that the cross-relation average and individual signatures can provide interesting comprehensive information about the topological properties of the analyzed structures, especially regarding the presence of hubs and the uniformity of the topology around the nodes at successive topological scales.

To complement this section, we estimate an analogue of the "power spectrum" (see Equation 4) of the mean au- torrelation signature of the US airport network by using . The magnitude of the Fourier transform of the average autorrelation function obtained for the US airports network is shown in Figure 10.

Concluding Remarks

An important class of relationship between a pair of mathematical structures regards how much they are interrelated in terms of a relative displacement, controlled by respective lags. The important of this type of relationship is reflect in the key roles these operations play in areas including signal and image processing and analysis, control systems, computer vision, neuronal networks and AI at large.

In the present work, we aimed at developing operations that are analogue to the autocorrelation and crosscorrelation between functions, fields, vectors and matrices, but adapted to graphs and networks. In order to do so, we resourced to the multiset coincidence similarity, which implements a strict and stable quantification of the similarity between any two mathematical structures.

More specifically, we defined the order of neighborhoods around each network node as being analogous to the lag parameter in the classic correlation. Therefore, it became possible to compare the topological features each reference node with those of nodes in successive respective neighorhoods. This operation yields a function of the distance from the successive neighborhoods, which acts as the lag parameter. The possibility to chose among sev- eral possible topological, or other types of features accounts for additional flexibility in the implemented analysis, which can be customized to the aspects of specific interest in each research approach. Once these measurements have been obtained, the estimation of the autorrelation and cross-relation involves relatively little computational cost. In the case of the present work, which adopted the node degree and transitivity as features, the overall cost resulted particularly low given that these measurements are local and easy to be calculated.

As with the traditional autocorrelation and crosscorrelation between two functions, the autorrelation and cross-relation have their maximum value observed for lag zero. As it has been indicated by several examples involving several model-theoretical and real-world networks, the average and individual autorrelation and cross-relation signatures reflect several aspects of the analyzed networks, including the presence of hubs, communities, the overall regularity, scale-free behavior, presence of outliers, as well as the extent of topological scales along which the network is more or less strongly mutually similar. Of particular interest is the intricate coincidence similarity obtained from the individual autorrelation signatures which, in the case of the US airports network revealed highly intricate interrelationships between involved signatures, confirming the particularly complex organization of this network.

The reported concepts and methods paves the way to several interesting further developments, including the consideration of other measurements, other modeltheoretical networks (especially heterogeneous modular networks), the investigation of other parameter configurations (e.g. average degree and number of nodes). Other promising possibilities include the application to other real-world networks, as well as further relating the autorrelation and cross-relation signatures with other interesting properties of networks such as their spectra and Fourier transforms.

Acknowledgments.

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

Note:

As all other preprints by the author, this work is possibly being considered by a scientific journal. Respective modification, commercial use, or distribution of any of its parts are not possible. Many of the preprints by the author are also available in HAL and arXiv. This work can also be cited by using the DOI number or article identification link. Thanks for reading.

Figure 1 :

 1 Figure 1: The autorrelation functions in terms of δ = 1, 2, . . . , 30 for ER, BA and WS networks with N = 300 nodes and k = 4. The WS network was obtained by using rewiring probability pw = 0.05.

Figure 2 :

 2 Figure 2: The autorrelation functions of an WS networks with N = 300 nodes, p = 4 and pw = 0.01(a), 0.02(b), and 0.03 (c), respectively to topological lags δ = 1, 2, . . . , 30.

Figure 3 :Figure 4 :

 34 Figure 3: The coincidence networks obtained from the individual autorrelation signatures of the ER (a), BA (b), and WS (c) networks, shown after thresholding by T = 0.87 (T = 0.6 in the BA case) and α = 0.5. The size of the nodes corresponds to the degree of the original network nodes providing the reference for the autorrelation analysis. The width of the edges is proportional to the respective coincidence value between the respective individual autorrelation signatures. The colors (from cyan to magenta) indicate the transitivity of the original nodes.

  The individual and average autorrelation signatures obtained for the BA network are shown in Figure 7(b). Interestingly, the obtained average crossrelation signature is characterized by a less steep initial increase, which is a direct consequence of the partial loss of the scale free property, as a consequence of the hubs being rewired into more uniform configurations, which is corroborated by the associated reduction of the number of lower individual signatures. The rewiring effects on the US airports network, shown in Figure 7(c), included making the intermediate peak less salient, reducing the extent of the signatures at the same time as they became less dispersed, as well as making the final transition sharper.

Figure 5 :

 5 Figure 5: The coincidence network obtained from the individual autorrelation signatures of the US airports network, shown after thresholding by T = 0.87 and α = 0.5. The size of the nodes corresponds to the degree of the original network nodes providing the reference for the autorrelation analysis. The width of the edges is proportional to the respective coincidence value between the respective individual autorrelation signatures. The colors (from cyan to magenta) indicate the transitivity of the original nodes.

Figure 6 :

 6 Figure 6: The coincidence networks obtained from the individual autorrelation signatures of the US airports network considering successive lags varying fro 1 to 9. By supplying additional information about the topology of the original network, networks derived from individual signatures with larger lag values yield more detailed coincidence networks.

Figure 7 :

 7 Figure 7: The average (in black) and individual (random colors) cross-relation signatures obtained between the ER network in Fig. 1(b) and the respectively network obtained by 500 edge rewirings (a), the original BA network and respective version after 4000 rewirings (b), as well as the US airports network and respective network obtained for 3000 rewirings (c).

Figure 8 :

 8 Figure 8: The cross-relation functions between the original ER, BA, and US airports networks and their respective versions after removal of the 5% of nodes with the highest degree (hubs).

Figure 9 :

 9 Figure 9: The cross-relation functions between the original ER, BA, and US airports networks and their respective versions after permutation of nodes.

Figure 10 :

 10 Figure10: The Fourier transform of the average autorrelation function obtained for the US airports network. Though this should not be taken for the respective power spectrum, the obtained Fourier transform has special interest as it highlights eventual oscillations present in the autorrelation function, which would be related to periodical properties eventually found in the original network topology.