text:
# Spectra of “Real–World” Graphs: Beyond the Semi-Circle Law

## I Introduction

Random graphs have long been used for modelling the evolution and topology of systems made up of large assemblies of similar units. The uncorrelated random graph model – which assumes each pair of the graph’s vertices to be connected with equal and independent probabilities – treats a network as an assembly of equivalent units. This model, introduced by the mathematicians Paul Erdős and Alfréd Rényi, has been much investigated in the mathematical literature. However, the increasing availability of large maps of real-life networks has indicated that real networks are fundamentally correlated systems, and in many respects their topology deviates from the uncorrelated random graph model. Consequently, the attention has shifted towards more advanced graph models which are designed to generate topologies in line with the existing empirical results . Examples of real networks, that serve as a benchmark for the current modelling efforts include the Internet , the World-Wide Web , networks of collaborating movie actors and those of collaborating scientists , the power grid  and the metabolic network of numerous living organisms .

These are the systems that we will call “real-world” networks or graphs. Several converging reasons explain the enhanced current interest in such real graphs. First, the amount of topological data available on such large structures has increased dramatically during the past few years thanks to the computerization of data collection in various fields, from sociology to biology. Second, the hitherto unseen speed of growth of some of these complex networks – e.g., the Internet – and their pervasiveness in affecting many aspects of our life has created the need to understand the topology, origin and evolution of such structures. Finally, the increased computational power available on almost every desktop has allowed, for the first time, to study such systems in unprecedented detail.

The proliferation of data has lead to a flurry of activity towards understanding the general properties of real networks. These efforts have resulted in the introduction of two classes of models, commonly called small-world graphs and the scale-free networks. The first aims to capture the clustering observed in real graphs, while the second reproduces the power law degree distribution present in many real networks. However, until now, most analyses of these models and data sets have been confined to real-space characteristics, that capture their static structural properties: e.g., degree sequences, shortest connecting paths and clustering coefficients. In contrast, there is extensive literature demonstrating that the properties of graphs and the associated adjacency matrices are well characterized by spectral methods, that provide global measures of the network properties. In this paper, we offer a detailed analysis of the most studied network models using algebraic tools intrinsic to large random graphs.

The paper is organized as follows. Section II introduces the main random graph models used for the topological description of large assemblies of connected units. Section III lists the – analytical and numerical – tools which we used and developed to convert the topological features of graphs into algebraic invariants. Section IV contains our results concerning the spectra and special eigenvalues of the three main types of random graph models: sparse uncorrelated random graphs in Section IV A; small-world graphs in Section IV B; and scale-free networks in Section IV C. Section IV D gives simple algorithms for testing the graph’s structure, and Section IV E investigates the variance of structure within single random graph models.

## II Models of random graphs

### A The uncorrelated random graph model and the semi-circle law

#### 1 Definitions

Throughout this paper, we will use the term “graph” for a set of points (vertices) connected by undirected lines (edges); no multiple edges and no loops connecting a vertex to itself are allowed. We will call two vertices of the graph “neighbors”, if they are connected by an edge. Based on Ref. , we shall use the term “uncorrelated random graph” for a graph, if (i) the probability for any pair of the graph’s vertices being connected is the same, $`p`$; (ii) these probabilities are independent variables.

Any graph $`G`$ can be represented by its adjacency matrix, $`A(G)`$, which is a real symmetric matrix: $`A_{ij}=A_{ji}=1`$, if vertices $`i`$ and $`j`$ are connected, or $`0`$, if these two vertices are not connected. The main algebraic tool which we will use for the analysis of graphs will be the spectrum – i.e., the set of eigenvalues – of the graph’s adjacency matrix. The spectrum of the graph’s adjacency matrix is also called the spectrum of the graph.

#### 2 Applying the semi-circle law for the spectrum of the uncorrelated random graph

A general form of the semi-circle law for real symmetric matrices is the following. If $`A`$ is a real symmetric $`N\times N`$ uncorrelated random matrix, $`A_{ij}0`$ and $`A_{ij}^2=\sigma ^2`$ for every $`ij`$, and with increasing $`N`$ each moment of each $`|A_{ij}|`$ remains finite, then in the $`N\mathrm{}`$ limit the spectral density – i.e., the density of eigenvalues – of $`A/\sqrt{N}`$ converges to the semi-circular distribution:

$`\rho (\lambda )=\{\begin{array}{cc}(2\pi \sigma ^2)^1\sqrt{4\sigma ^2\lambda ^2},\hfill & \text{if }|\lambda |<2\sigma \text{;}\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}`$ (1)
This theorem is also known as Wigner’s law , and its extensions to further matrix ensembles have long been used for the stochastic treatment of complex quantum mechanical systems lying far beyond the reach of exact methods . Later, the semi-circle law was found to have many applications in statistical physics and solid-state physics as well .

Note, that for the adjacency matrix of the uncorrelated random graph many of the semi-circle law’s conditions do not hold, e.g., the expectation value of the entries is a non-zero constant: $`p0`$. Nevertheless, in the $`N\mathrm{}`$ limit, the rescaled spectral density of the uncorrelated random graph converges to the semi-circle law of Eq. (1). An illustration of the convergence of the average spectral density to the semi-circular distribution can be seen on Fig. 1. It is necessary to make a comment concerning figures here. In order to keep figures simple, for the spectral density plots we have chosen to show the spectral density of the original matrix, $`A`$, and to rescale the horizontal ($`\lambda `$) and vertical ($`\rho `$) axes by $`\sigma ^1N^{1/2}=[Np(1p)]^{1/2}`$ and $`\sigma N^{1/2}=[Np(1p)]^{1/2}`$.

Some further results on the behavior of the uncorrelated random graph’s eigenvalues, relevant for the analysis of real-world graphs as well, include the following: The principal eigenvalue (the largest eigenvalue, $`\lambda _1`$) grows much faster than the second eigenvalue: $`lim_N\mathrm{}(\lambda _1/N)=p`$ with probability 1, whereas for every $`ϵ>1/2`$, $`lim_N\mathrm{}(\lambda _2/N^ϵ)=0`$. (see Refs.  and Fig. 1.). A similar relation holds for the smallest eigenvalue, $`\lambda _N`$: for every $`ϵ>1/2`$, $`lim_N\mathrm{}(\lambda _N/N^ϵ)=0`$. In other words: if $`k_i`$ denotes the average number of connections of a vertex in the graph, then $`\lambda _1`$ scales as $`pNk_i`$, and the width of the “bulk” part of the spectrum – the set of the eigenvalues $`\{\lambda _2,\mathrm{},\lambda _N\}`$ – scales as $`\sigma \sqrt{N}`$. Lastly, the semi-circular distribution’s edges are known to decay exponentially, and the number of eigenvalues in the $`\lambda >𝒪(\sqrt{N})`$ tail has been shown to be of the order of $`1`$.

### B Real-world graphs

The two main models proposed to describe real-world graphs are the small-world model and the scale-free model.

#### 1 Small-world graphs

The small-world graph is created by randomly rewiring some of the edges of a regular ring graph. The regular ring graph is created as follows. First draw the vertices $`1`$, $`2`$, …, $`N`$ on a circle in ascending order. Then, for every $`i`$, connect vertex $`i`$ to the vertices lying closest to it on the circle: vertices $`ik/2`$, …, $`i1`$, $`i+1`$, …, $`i+k/2`$, where every number should be understood modulo $`N`$ ($`k`$ is an even number). Fig. 9 shows, that this algorithm creates a regular graph indeed, because the degree of any vertex is the same number, $`k`$. Next, starting from vertex $`1`$ and proceeding towards $`N`$, perform the rewiring step. For vertex $`1`$, consider the first “forward connection”, i.e., the connection to vertex $`2`$. With probability $`p_\mathrm{r}`$, reconnect vertex $`1`$ to another vertex chosen uniformly at random and without allowing multiple edges. Proceed towards the remaining forward connections of vertex $`1`$, and then, perform this step for the remaining $`N1`$ vertices also. For the rewiring, use equal and independent probabilities. Note, that in the small-world model the density of edges is $`p=k_i/(N1)k/N`$. Throughout this paper, we will use only $`k>2`$.

If we use $`p_\mathrm{r}=0`$ in the small-world model, the original regular graph is preserved, and for $`p_\mathrm{r}=1`$, one obtains a random graph which differs from the uncorrelated random graph only slightly: every vertex has a minimum degree of $`k/2`$. Next, we will need two definitions. The separation between vertices $`i`$ and $`j`$ – denoted by $`L_{ij}`$ – is the number of edges in the shortest path connecting them. The clustering coefficient at vertex $`i`$ – denoted by $`C_i`$ – is the number of existing edges among the neighbors of vertex $`i`$ divided by the number of all possible connections between them. In the small-world model, both $`L_{ij}`$ and $`C_i`$ are functions of the rewiring probability, $`p_\mathrm{r}`$. Based on the above definitions of $`L_{ij}(p_\mathrm{r})`$ and $`C_i(p_\mathrm{r})`$, the characteristics of the small-world phenomenon – which occurs for intermediate values of $`p_\mathrm{r}`$ – can be given as follows: (i) the average separation between two vertices, $`L(p_\mathrm{r})`$, drops dramatically below $`L(p_\mathrm{r}=0)`$, whereas (ii) the average clustering coefficient, $`C(p_\mathrm{r})`$, remains high, close to $`C(p_\mathrm{r}=0)`$. Note, that the rewiring procedure is carried out independently for every edge, therefore the degree sequence and also other distributions in the system – e.g., path length and loop size – decay exponentially.

#### 2 The scale-free model

The scale-free model assumes a random graph to be a growing set of vertices and edges, where the location of new edges is determined by a preferential attachment rule. Starting from an initial set of $`m_0`$ isolated vertices, one adds $`1`$ new vertex and $`m`$ new edges at every time step $`t`$. (Throughout this paper, we will use $`m=m_0`$.) The $`m`$ new edges connect the new vertex and $`m`$ different vertices chosen from the $`N`$ old vertices. The $`i`$th old vertex is chosen with probability $`k_i/_{j=1,N}k_j`$, where $`k_i`$ is the degree of vertex $`i`$. (The density of edges in a scale-free graph is $`p=k_i/(N1)2m/N`$.) In contrast to the small-world model, the distribution of degrees in a scale-free graph converges to a power-law when $`N\mathrm{}`$, which has been shown to be a combined effect of growth and the preferential attachment. Thus, in the infinite time or size limit, the scale-free model has no characteristic scale in the degree size .

#### 3 Related models

Lately, numerous other models have been suggested for a unified description of real-world graphs . Models of growing networks with aging vertices were found to display both heavy tailed and exponentially decaying degree sequences  as a function of the speed of aging. Generalized preferential attachment rules have helped us better understand the origin of the exponents and correlations emerging in these systems. Also, investigations of more complex network models – using aging or an additional fixed cost of edges  or preferential growth and random rewiring  – have shown, that in the “frequent rewiring, fast aging, high cost” limiting case, one obtains a graph with an exponentially decaying degree sequence, whereas in the “no rewiring, no aging, zero cost” limiting case the degree sequence will decay as a power-law. According to studies of scientific collaboration networks  and further social and biological structures , a significant proportion of large networks lies between the two extremes. In such cases, the characterization of the system using a small number of algebraic constants could facilitate the classification of real-world networks.

## III Tools

### A Analytical

#### 1 The spectrum of the graph

The spectrum of a graph is the set of eigenvalues of the graph’s adjacency matrix. The physical meaning of a graph’s eigenpair (an eigenvector and its eigenvalue) can be illustrated by the following example. Write each component of a vector $`\stackrel{}{v}`$ on the corresponding vertex of the graph: $`v_i`$ on vertex $`i`$. Next, on every vertex write the sum of the numbers found on the neighbors of vertex $`i`$. If the resulting vector is a multiple of $`\stackrel{}{v}`$, then $`\stackrel{}{v}`$ is an eigenvector, and the multiplier is the corresponding eigenvalue of the graph.

The spectral density of a graph is the density of the eigenvalues of its adjacency matrix. For a finite system, this can be written as a sum of delta functions

$`\rho (\lambda ):={\displaystyle \frac{1}{N}}{\displaystyle \underset{j=1}{\overset{N}{}}}\delta (\lambda \lambda _j),`$ (2)
which converges to a continuous function with $`N\mathrm{}`$ ($`\lambda _j`$ is the $`j`$th largest eigenvalue of the graph’s adjacency matrix).

The spectral density of a graph can be directly related to the graph’s topological features: the $`k`$th moment, $`M_k`$, of $`\rho (\lambda )`$ can be written as

$`M_k={\displaystyle \frac{1}{N}}{\displaystyle \underset{j=1}{\overset{N}{}}}(\lambda _j)^k={\displaystyle \frac{1}{N}}\mathrm{Tr}\left(A^k\right)=`$ (3)
$`={\displaystyle \frac{1}{N}}{\displaystyle \underset{i_1,i_2,\mathrm{},i_k}{}}A_{i_1,i_2}A_{i_2,i_3}\mathrm{}A_{i_k,i_1}.`$ (4)
From the topological point of view, $`D_k=NM_k`$ is the number of directed paths (loops) of the underlying – undirected – graph, that return to their starting vertex after $`k`$ steps. On a tree, the length of any such path can be an even number only, because these paths contain any edge an even number of times: once such a path has left its starting point by chosing a starting edge, no alternative route for returning to the starting point is available. However, if the graph contains loops of odd length, the path length can be an odd number, as well.

#### 2 Extremal eigenvalues

In an uncorrelated random graph the principal eigenvalue, $`\lambda _1`$, shows the density of edges and $`\lambda _2`$ can be related to the conductance of the graph as a network of resistances. An important property of all graphs is the following: the principal eigenvector, $`\stackrel{}{e}_1`$, of the adjacency matrix is a non-negative vector (all components are non-negative), and if the graph has no isolated vertices, $`\stackrel{}{e}_1`$ is a positive vector. All other eigenvectors are orthogonal to $`\stackrel{}{e}_1`$, therefore they all have entries with mixed signs.

#### 3 The inverse participation ratios of eigenvectors

The inverse participation ratio of the normalized $`j`$th eigenvector, $`\stackrel{}{e}_j`$, is defined as

$`I_j={\displaystyle \underset{k=1}{\overset{N}{}}}\left[(e_j)_k\right]^4.`$ (5)
If the components of an eigenvector are identical – $`(e_j)_i=1/\sqrt{N}`$ for every $`i`$, – then $`I_j=1/N`$. For an eigenvector with one single non-zero component – $`(e_j)_i=\delta _{i,i^{}}`$ – the inverse participation ratio is $`1`$. The comparison of these two extremal cases illustrates that with the help of the inverse participation ratio, one can tell whether only $`𝒪(1)`$ or as many as $`𝒪(N)`$ components of an eigenvector differ significantly from $`0`$, i.e., whether an eigenvector is localized or non-localized.

### B Numerical

#### 1 General real symmetric eigenvalue solver

To compute the eigenpairs of graphs below the size $`N=5,000`$, we used the general real symmetric eigenvalue solver of Ref. . This algorithm requires the allocation of memory space to all entries of the matrix, thus to compute the spectrum of a graph of size $`N=20,000`$ ($`N=1,000,000`$) using this general method with double precision floating point arithmetic, one would need $`3.2`$ GB ($`8`$ TB) memory space and the execution of approximately $`30N^2=1.2\times 10^{10}`$ ($`3\times 10^{13}`$) floating point operations . Consequently, we need to develop more efficient algorithms to investigate the properties of graphs with sizes comparable to real-world networks.

#### 2 Iterative eigenvalue solver based on the thick-restart Lanczos algorithm

The spectrum of a real-world graph is the spectrum of a sparse real symmetric matrix, therefore the most efficient algorithms which can give a handful of the top $`n_\mathrm{d}`$ eigenvalues – and the corresponding eigenvectors – of a large graph are iterative methods. These methods allow the matrix to be stored in any compact format, as long as matrix-vector multiplication can be carried out at a high speed. Iterative methods use little memory: only the non-zero entries of the matrix and a few vectors of size $`N`$ need to be stored. The price for computational speed lies in the number of the obtained eigenvalues: iterative methods compute only a handful of the largest (or smallest) eigenvalues of a matrix. To compute the eigenvalues of graphs above the size $`N=5,000`$, we have developed algorithms using a specially modified version of the thick-restart Lanczos algorithm. The modifications and some of the main technical parameters of our software are explained in the following paragraphs.

Even though iterative eigenvalue methods are mostly used to obtain the top eigenvalues of a matrix, after minor modifications the internal eigenvalues in the vicinity of a fixed $`\lambda =\lambda _0`$ point can be computed as well. For this, extremely sparse matrices are usually “shift-inverted”, i.e., to find those eigenvalues of $`A`$ that are closest to $`\lambda _0`$, the highest and lowest eigenvalues of $`(A\lambda _0I)^1`$ are searched for. However, because of the extremely high cost of matrix inversion in our case, for the computation of internal eigenvalues we suggest to use the “shift-square” method with the matrix

$`B=\left[\lambda ^{}/2(A\lambda _0I)^2\right]^{2n+1}.`$ (6)
Here $`\lambda ^{}`$ is the largest eigenvalue of $`(A\lambda _0I)^2`$, $`I`$ is the identity matrix and $`n`$ is a positive integer. Transforming the matrix $`A`$ into $`B`$ transforms the spectrum of $`A`$ in the following manner. First, the spectrum is shifted to the left by $`\lambda _0`$. Then, the spectrum is “folded” (and squared) at the origin such that all eigenvalues will be negative. Next, the spectrum is linearly rescaled and shifted to the right, with the following effect: (i) the whole spectrum will lie in the symmetric interval $`[\lambda ^{}/2,\lambda ^{}/2]`$ and (ii) those eigenvalues that were closest to $`\lambda _0`$ in the spectrum of $`A`$, will be the largest now, i.e., they will be the eigenvalues closest to $`\lambda ^{}/2`$. Now, raising all eigenvalues to the $`(2n+1)`$st power increases the relative difference, $`1\lambda _i/\lambda _j`$, between the top eigenvalues $`\lambda _i`$ and $`\lambda _j`$ by a factor of $`2n+1`$. This allows the iterative method find the top eigenvalues of $`B`$ more quickly. One can compute the corresponding eigenvalues (those being closest to $`\lambda _0`$) of the original matrix, $`A`$: if $`\stackrel{}{b}_1`$, $`\stackrel{}{b}_2`$, …, $`\stackrel{}{b}_{n_\mathrm{d}}`$ are the normalized eigenvectors of the $`n_\mathrm{d}`$ largest eigenvalues of $`B`$, then for $`A`$ the $`n_\mathrm{d}`$ eigenvalues closest to $`\lambda _0`$ will be – not necessarily in ascending order – $`\stackrel{}{b}_1A\stackrel{}{b}_1`$, $`\stackrel{}{b}_2A\stackrel{}{b}_2`$, …, $`\stackrel{}{b}_{n_\mathrm{d}}A\stackrel{}{b}_{n_\mathrm{d}}`$.

The thick-restart Lanczos method uses memory space for the non-zero entries of the $`N\times N`$ large adjacency matrix, and $`n_\mathrm{g}+1`$ vectors of length $`N`$, where $`n_\mathrm{g}`$ ($`n_\mathrm{g}>n_\mathrm{d}`$) is usually between $`10`$ and $`100`$. Besides the relatively small size of required memory, we could also exploit the fact that the non-zero entries of a graph’s adjacency matrix are all 1’s: during matrix-vector multiplication – which is usually the most time-consuming step of an iterative method – only additions had to be carried out instead of multiplications.

The numerical spectral density functions of large graphs ($`N5,000`$) of this paper were obtained using the following steps. To compute the spectral density of the adjacency matrix, $`A`$, at an internal $`\lambda =\lambda _0`$ location, first, the $`n_\mathrm{d}`$ eigenvalues closest to $`\lambda _0`$ were searched for. Next, the distance between the smallest and the largest of the obtained eigenvalues was computed. Finally, $`\rho (\lambda _0)`$ this distance was multiplied by $`N/(n_\mathrm{d}1)`$, and was averaged using $`n_{\mathrm{av}}`$ different graphs. We used double precision floating point arithmetic, and the iterations were stopped, if (i) at least $`n_{\mathrm{it}}`$ iterations had been carried out and (ii) the lengths of the residual vectors belonging to the $`n_\mathrm{d}`$ selected eigenpairs were all below $`\epsilon =10^{12}`$.

## IV Results

### A Sparse uncorrelated random graphs: the semi-circle law is not universal

In the uncorrelated random graph model of Erdős and Rényi, the total number of edges grows quadratically with the number of vertices: $`N_{\mathrm{edge}}=Nk_i=Np(N1)pN^2`$. However, in many real-world graphs edges are “expensive”, and the growth rate of the number of connections remains well below this rate. For this reason, we also investigated the spectra of such uncorrelated networks, for which the probability of any two vertices being connected changes with the size of the system using $`pN^\alpha =c=\mathrm{const}`$. Two special cases are $`\alpha =0`$ (the Erdős-Rényi model) and $`\alpha =1`$. In the second case, $`pN\mathrm{const}`$ as $`N\mathrm{}`$, i.e., the average degree remains constant.

For $`\alpha <1`$ and $`N\mathrm{}`$, there exists an infinite cluster of connected vertices (in fact, it exists for every $`\alpha 1`$ ). Moreover, the expectation value of any $`k_i`$ converges to infinity, thus any vertex is almost surely connected to the infinite cluster. The spectral density function converges to the semicircular distribution of Eq. (1), because the total weight of isolated subgraphs decreases exponentially with growing system size. (A detailed analysis of this issue is available in Ref. .)

For $`\alpha =1`$ and $`N\mathrm{}`$ (see Fig. 2), the probability for a vertex to belong to a cluster of any finite size remains also finite . Therefore, the limiting spectral density contains the weighted sum of the spectral densities of all finite graphs. The most striking deviation from the semi-circle law in this case is the elevated central part of the spectral density. The probability for a vertex to belong to an isolated cluster of size $`s`$ decreases exponentially with $`s`$, therefore the number of large isolated clusters is low and the eigenvalues of a graph with $`s`$ vertices are bounded by $`\sqrt{s1}`$ and $`\sqrt{s1}`$. For these two reasons, the amplitudes of delta functions decay exponentially, as the absolute value of their locations, $`|\lambda |`$, increases.

The principal eigenvalue of this graph converges to a constant: $`lim_N\mathrm{}(\lambda _1)=pN=c`$, and $`\rho (\lambda )`$ will be symmetric in the $`N\mathrm{}`$ limit. Therefore, in the limit, all odd moments ($`M_{2k+1}`$), and thus, the number of all loops with odd length ($`D_{2k+1}`$) disappear. This is a salient feature of graphs with tree structure (because on a tree every edge must be used an even number of times in order to return to the initial vertex), indicating that the structure of a sparse uncorrelated random graph becomes more and more tree-like. This can also be understood by considering that the typical distance (length of the shortest path) between two vertices on both a sparse uncorrelated random graph and a regular tree with the same number of edges scales as $`\mathrm{ln}(N)`$. So except for a few short-cuts a sparse uncorrelated random graph looks like a tree.

### B The small-world graph

#### 1 Triangles are abundant in the graph

For $`p_\mathrm{r}=0`$, the small-world graph is regular and also periodical. Because of the highly ordered structure, $`\rho (\lambda )`$ contains numerous singularities, which are listed in Section VI A (see also Fig. 3). Note that $`\rho (\lambda )`$ has a high third moment. (Remember, that we use only $`k>2`$.)

If we increase $`p_\mathrm{r}`$ such that the small-world region is reached, i.e., the periodical structure of the graph is perturbed, then singularities become blurred and are transformed into high local maxima, but $`\rho (\lambda )`$ retains a strong skewness (see Fig. 3). This is in good agreement with results of Refs. , where it has been shown that the local structure of the small-world graph is ordered, however, already a very small number of shortcuts can drastically change the graph’s global structure.

In the $`p_\mathrm{r}=1`$ case the small-world model becomes very similar to the uncorrelated random graph: the only difference is that here, the minimum degree of any vertex is a positive constant, $`k/2`$, whereas in an uncorrelated random graph the degree of a vertex can be any non-negative number. Accordingly, $`\rho (\lambda )`$ becomes a semi-circle for $`p_\mathrm{r}=1`$ (Fig. 3). Nevertheless, it should be noted, that as $`p_\mathrm{r}`$ converges to $`1`$, a high value of $`M_3`$ is preserved even for $`p_\mathrm{r}`$ close to $`1`$, where all local maxima have already vanished. The third moment of $`\rho (\lambda )`$ gives the number of triangles in the graph (see Section III A 1); the lack of high local maxima – i.e., the remnants of singularities – shows the absence of an ordered structure.

From the above we conclude, that – from the spectrum’s point of view – the high number of triangles is one of the most basic properties of the small-world model, and it is preserved much longer, than regularity or periodicity, if the level of randomness, $`p_\mathrm{r}`$, is increased. This is in good agreement with the results of Ref. , where the high number of small cycles is found to be a fundamental property of small-world networks. As an application, the high number of small cycles results in special diffusion properties on small-world graphs .

### C The scale-free graph

For $`m=m_0=1`$, the scale-free graph is a tree by definition and its spectrum is symmetric. In the $`m>1`$ case $`\rho (\lambda )`$ consists of several well distinguishable parts (see Fig. 4). The “bulk” part of the spectral density – the set of the eigenvalues {$`\lambda _2,\mathrm{},\lambda _N`$} – converges to a symmetric continuous function which has a triangle-like shape for the normalized $`\lambda `$ values up to $`1.5`$ and has power-law tails.

The central part of the spectral density converges to a triangle-like shape with its top lying well above the semi-circle. Since the scale-free graph is fully connected by definition, the increased number of eigenvalues with small magnitudes cannot be accounted to isolated clusters, as before. As an explanation, we suggest, that the eigenvectors of these eigenvalues are localized on a small subset of the graph’s vertices. (This idea is supported by the high inverse participation ratios of these eigenvectors, see Fig. 7).

#### 1 The spectral density of the scale-free graph decays as a power-law

The inset of Fig. 4 shows the tail of the bulk part of the spectral density for a graph with $`N=20,000`$ vertices and $`100,000`$ edges (i.e., $`pN=10`$). Comparing this to the inset of Fig. 1 – where the number of vertices and edges is the same as here, – one can observe the power law decay at the edge of the bulk part of $`\rho (\lambda )`$. As shown later, in Section IV D, the power-law decay in this region is caused by eigenvectors localized on vertices with the highest degrees. The power law decay of the degree sequence – i.e., the existence of very high degrees – is, in turn, due to the preferential attachment rule of the scale-free model.

#### 2 The growth rate of the principal eigenvalue shows a crossover in the level of correlations

Since the adjacency matrix of a graph is a non-negative symmetric matrix, the graph’s largest eigenvalue, $`\lambda _1`$, is also the largest in magnitude (see, e.g., Theorem 0.2 of Ref. ). Considering the effect of the adjacency matrix on the base vectors $`(b_i)_j=\delta _{ij}`$ ($`i=1,2,\mathrm{},N`$), it can be shown that a lower bound for $`\lambda _1`$ is given by the length of the longest row vector of the adjacency matrix, which is the square root of the graph’s largest degree, $`k_1`$. Knowing that the largest degree of a scale-free graph grows as $`\sqrt{N}`$ , one expects $`\lambda _1`$ to grow as $`N^{1/4}`$ for large enough systems.

Fig. 5 shows a rescaled plot of the scale-free graph’s largest eigenvalue for different values of $`m`$. In this figure, $`\lambda _1`$ is compared to the length of the longest row vector, $`\sqrt{k_1}`$ on the “natural scale” of these values, which is $`\sqrt{m}N^{1/4}`$. It is clear that if $`m>1`$ and the system is small, then through several decades (a) $`\lambda _1`$ is larger than $`\sqrt{k_1}`$ and (b) the growth rate of $`\lambda _1`$ is well below the expected rate of $`N^{1/4}`$. In the $`m=1`$ case, and for large systems, (a) the difference between $`\lambda _1`$ and $`\sqrt{k_1}`$ vanishes and (b) the growth rate of the principal eigenvalue will be maximal, too. This crossover in the behavior of the scale-free graph’s principal eigenvalue is a specific property of sparse growing correlated graphs, and it is a result of the changing level of correlations between the longest row vectors (see Section VI B).

#### 3 Comparing the role of the principal eigenvalue in the scale-free graph and the $`\alpha =1`$ uncorrelated random graph: a comparison of structures

Now we will compare the role of the principal eigenvalue in the $`m>1`$ scale-free graph and the $`\alpha =1`$ uncorrelated random graph through its effect on the moments of the spectral density. On Figs. 4 and 5, one can observe that (i) the principal eigenvalue of the scale-free graph is detached from the rest of the spectrum, and (ii) as $`N\mathrm{}`$, it grows as $`N^{1/4}`$ (see also Sections IV C 2 and VI B). The inset of Fig. 4 indicates, that in the $`N\mathrm{}`$ limit, the bulk part will be symmetric, and its width will be constant (Fig. 4 rescales this constant width merely by another constant, namely $`Np[(1p)]^{1/2}`$). Because of the symmetry of the bulk part, in the $`N\mathrm{}`$ limit, the third moment of $`\rho (\lambda )`$ is determined exclusively by the contribution of the principal eigenvalue, which is $`N^1(\lambda _1)^3N^{1/4}`$. For each moment above the third (e.g., for the $`l`$th moment), with growing $`N`$, the contribution of the bulk part to this moment will scale as $`𝒪(1)`$, and the contribution of the principal eigenvalue will scale as $`N^{1+l/4}`$. In summary, in the $`N\mathrm{}`$ limit, the scale-free graph’s first eigenvalue has a significant contribution to the fourth moment; the 5th and all higher moments are determined exclusively by $`\lambda _1`$: the $`l`$th moment will scale as $`N^{1+l/4}`$.

In contrast to the above, the principal eigenvalue of the $`\alpha =1`$ uncorrelated random graph converges to the constant $`pN=c`$ in the $`N\mathrm{}`$ limit, and the width of the bulk part also remains constant (see Fig. 2). Given a fixed number, $`l`$, the contribution of the principal eigenvalue to the $`l`$th moment of the spectral density will change as $`N^1c^l`$ in the $`N\mathrm{}`$ limit. The contribution of the bulk part will scale as $`𝒪(1)`$, therefore all even moments of the spectral density will scale as $`𝒪(1)`$ in the $`N\mathrm{}`$ limit, and all odd moments will converge to $`0`$.

The difference between the growth rate of the moments of $`\rho (\lambda )`$ in the above two models (scale-free graph and $`\alpha =1`$ uncorrelated random graph model) can be interpreted as a sign of different structure (see Section III A 1). In the $`N\mathrm{}`$ limit, the average degree of a vertex converges to a constant in both models: $`lim_N\mathrm{}k_i=pN=c=2m`$. (Both graphs will have the same number of edges per vertex.) On the other hand, in the limit, all moments of the $`\alpha =1`$ uncorrelated random graph’s spectral density converge to a constant, whereas the moments $`M_l(l=5,6,\mathrm{})`$ of the scale-free graph’s $`\rho (\lambda )`$ will diverge as $`N1+l/4`$. In other words: the number of loops of length $`l`$ in the $`\alpha =1`$ uncorrelated random graph will grow as $`D_l=NM_l=𝒪(N)`$, whereas for the scale-free graph for every $`l3`$, the number of these loops will grow as $`D_l=NM_l=𝒪(N^{l/4})`$. From this we conclude that in the limit, the role of loops is negligible in the $`\alpha =1`$ uncorrelated random graph, whereas it is large in the scale-free graph. In fact, the growth rate of the number of loops in the scale free graph exceeds all polynomial growth rates: the longer the loop size ($`l`$) investigated, the higher the growth rate of the number of these loops ($`N^{l/4}`$) will be. Note that the relative number of triangles (i.e., the third moment of the spectral density, $`M_l/N`$) will disappear in the scale-free graph, if $`N\mathrm{}`$.

In summary, the spectrum of the scale-free model converges to a triangle-like shape in the center, and the edges of the bulk part decay slowly. The first eigenvalue is detached from the rest of the spectrum, and it shows an anomalous growth rate. Eigenvalues with large magnitudes belong to eigenvectors localized on vertices with many neighbors. In the present context, the low number of triangles, the high number of loops with length above $`l=4`$, and the buildup of correlations are the basic properties of the scale-free model.

### D Testing the structure of a “real-world” graph

To analyze the structure of a large sparse random graph (correlated or not), here we suggest several tests, that can be performed within $`𝒪(N)`$ CPU time, use $`𝒪(N)`$ floating point operations, and can clearly differentiate between the three ”pure” types of random graph models treated in Section IV. Furthermore, these tests allow one to quantify the relation between any real-world graph and the three basic types of random graphs.

#### 1 Extremal eigenvalues

In Section III A 2 we have already mentioned that the extremal eigenvalues contain useful information on the structure of the graph. As the spectra of uncorrelated random graphs (Fig. 1) and scale-free networks (Fig. 4) show, the principal eigenvalue of random graphs is often detached from the rest of the spectrum. For these two network types, the remaining bulk part of the spectrum – i.e., the set {$`\lambda _2`$, …, $`\lambda _N`$} – converges to a symmetric distribution, thus the quantity

$`R:={\displaystyle \frac{\lambda _1\lambda _2}{\lambda _2\lambda _N}}`$ (7)
measures the distance of the first eigenvalue from the main part of $`\rho (\lambda )`$ normalized by the extension of the main part. ($`R`$ can be connected to the chromatic number of the graph.)

Note, that in the $`N\mathrm{}`$ limit the $`\alpha =0`$ sparse uncorrelated random graph’s principal eigenvalue will scale as $`k_i`$, whereas both $`\lambda _2`$ and $`|\lambda _N|`$ will scale as $`2\sqrt{k_i}`$. Therefore, if $`k_i>4`$, the principal eigenvalue will be detached from the bulk part of the spectrum and $`R`$ will scale as $`(\sqrt{k_i}2)/4`$. If, however, $`k_i4`$, $`\lambda _1`$ will not be detached from the bulk part, and it will converge to $`0`$.

The above explanation and Fig. 6 show that in the $`k_i>4`$ sparse uncorrelated random graph model and the scale-free network, $`\lambda _1`$ and the rest of the spectrum are well separated, which gives similarly high values for $`R`$ in small systems. In large systems, $`R`$ of the sparse uncorrelated random graph converges to a constant, while $`R`$ in the scale-free model decays as a power law function of $`N`$. The reason for this drop is the increasing denominator on the rhs. of Eq. 7: $`\lambda _2`$ and $`\lambda _N`$ are the extremal eigenvalues in the lower and upper long tails of $`\rho (\lambda )`$, therefore, as $`N`$ increases, the expectation values of $`\lambda _2`$ and $`\lambda _N`$ grow as quickly as that of $`\lambda _1`$. On the other hand, the small-world network shows much lower values of $`R`$ already for small systems: here, $`\lambda _1`$ is not detached from the rest of the spectrum, which is a consequence of the almost periodical structure of the graph.

On Fig. 6, graphs with the same number of vertices and edges are compared. For large ($`N10,000`$) systems, for sparse uncorrelated random graphs $`R`$ converges to a constant, whereas for scale-free graphs and small-world networks it decays as a power-law. The latter two networks significantly differ in the magnitude of $`R`$. In summary, the suggested quantity $`R`$ has been shown to be appropriate for distinguishing between the following graph structures: (i) periodical or almost periodical (small-world), (ii) uncorrelated non-periodical and (iii) strongly correlated non-periodical (scale-free).

#### 2 Inverse participation ratios of extremal eigenpairs

Figure 7 shows the inverse participation ratios of the eigenvectors of an uncorrelated random graph, a small-world graph with $`p_\mathrm{r}=0.01`$ and a scale-free graph. Even though all three graphs have the same number of vertices ($`N=1,000`$) and edges ($`5,000`$), one can observe rather specific features (see also the inset of Fig. 7).

The uncorrelated random graph’s eigenvectors show very little difference in their level of localization, except for the principal eigenvector, which is much less localized than the other eigenvectors; $`I(\lambda _2)`$ and $`I(\lambda _N)`$ are almost equal. For the small-world graph’s eigenvectors, $`I(\lambda )`$ has many different plateaus and spikes; the principal eigenvector is not localized, and the $`2`$nd and $`N`$th eigenvectors have high, but different $`I(\lambda )`$ values. The eigenvectors belonging to the scale-free graph’s largest and smallest eigenvalues are localized on the “largest” vertices. The long tails of the bulk part of $`\rho (\lambda )`$ are due to these vertices. All three investigated eigenvectors ($`\stackrel{}{e}_1`$, $`\stackrel{}{e}_2`$ and $`\stackrel{}{e}_N`$) of the scale-free graph are highly localized. Consequently, the inverse participation ratios of the eigenvectors $`\stackrel{}{e}_1`$, $`\stackrel{}{e}_2`$ and $`\stackrel{}{e}_N`$ are handy for the identification of the $`3`$ basic types of random graph models used.

### E Structural variances

#### 1 Relative variance of the principal eigenvalue for different types of networks: the scale-free graph and self-similarity

Figure 8 shows the relative variance of the principal eigenvalue – i.e., $`\sigma (\lambda _1)/E(\lambda _1)`$ – for the three basic random graph types.

For non-sparse uncorrelated random graphs ($`N\mathrm{}`$ and $`p=\mathrm{const}`$) this quantity is known to decay at a rate which is faster than exponential. Comparing sparse graphs with the same number of vertices and edges, one can see, that in the sparse uncorrelated random graph and the small-world model the relative variance of the principal eigenvalue drops quickly with growing system size. In the scale-free model, however, the relative variance of the principal eigenvalue’s distribution remains constant with an increasing number of vertices.

In fractals, fluctuations do not disappear as the size of the system is increased, while in the scale-free graph, the relative variance of the principal eigenvalue is independent of system size. In this sense, the scale-free graph resembles self-similar systems.

## V Conclusions

We have performed a detailed analysis of the complete spectra, eigenvalues and the eigenvectors’ inverse participation ratios in three types of sparse random graphs: the sparse uncorrelated random graph, the small-world model and the scale-free network. Connecting the topological features of these graphs to algebraic quantities, we have demonstrated that (i) the semi-circle law is not universal, not even for the uncorrelated random graph model; (ii) the small-world graph is inherently non-correlated and contains a high number of triangles; (iii) the spectral density of the scale-free graph is made up of three, well distinguishable parts (center, tails of bulk, first eigenvalue), and as $`N\mathrm{}`$, triangles become negligible and the level of correlations changes.

We have presented practical tools for the identification of the above mentioned basic types of random graphs and further, for the classification of real-world graphs. The robust eigenvector techniques and observations outlined in this paper combined with previous studies are likely to improve our understanding of large sparse correlated random structures. Examples for algebraic techniques already in use for large sparse correlated random structures are analyses of the Internet and search engines and mappings of the World-Wide Web. Besides the improvement of these techniques, the present work may turn out to be useful for analyzing the correlation structure of the transactions between a very high number of economical and financial units, which has already been started in e.g., Refs. . Lastly, we hope to have provided quantitative tools for the classification of further “real-world” networks, e.g., social and biological networks.

ACKNOWLEDGEMENTS

We thank D. Petz, G. Stoyan, G. Tusnády, K. Wu and B. Kahng for helpful discussions and comments. This research was partially supported by a HNSF Grant No. OTKA T033104 and NSF PHY-9988674.

After submitting our manuscript we were made aware of a manuscript of K.-I. Goh, B. Kahng and D. Kim , which also investigates the eigenvalues of random graph models.

## VI Appendices

### A The spectrum of a small-world graph for $`p_\mathrm{r}=0`$ rewiring probability

#### 1 Derivation of the spectral density

If the rewiring probability of a small-world graph is $`p_\mathrm{r}=0`$, then the graph is regular, each vertex is connected to its $`k`$ nearest neighbors and the eigenvalues can be computed using the graph’s symmetry operations. Rotational symmetry operations can be easily recognized, if the vertices of the graph are drawn along the perimeter of a circle (see Fig. 9): let $`P^{(n)}`$ ($`n=0,1,\mathrm{},N1`$) denote the symmetry operation that rotates the graph by $`n`$ vertices in the anticlockwise direction. Being a symmetry operation, each $`P^{(n)}`$ commutes with the adjacency matrix, $`A`$, and they have a common full orthogonal system of eigenvectors.

Now, we will create a full orthogonal basis of $`A`$. (We will treat only the case when $`N`$ is an even number; odd $`N`$’s can be treated similarly.) It is known, that the eigenvalues of $`A`$ are real, however, to simplify calculations, we will use complex numbers first. The eigenvectors of every $`P^{(n)}`$ are $`\stackrel{}{e}_1,\stackrel{}{e}_2,\mathrm{},\stackrel{}{e}_N`$:

$`(e_l)_j=\mathrm{exp}\left(2\pi i{\displaystyle \frac{jl}{N}}\right),`$ (8)
where $`l=0,2,\mathrm{},N1`$ and $`i=\sqrt{1}`$. The eigenvalue of $`P^{(n)}`$ on $`\stackrel{}{e}_l`$ is

$`s_l^{(n)}=\mathrm{exp}\left(2\pi i{\displaystyle \frac{nl}{N}}\right).`$ (9)
By adding these values pairwise, one can obtain the $`N`$ eigenvalues of the graph:

$`\lambda _l=2{\displaystyle \underset{j=1}{\overset{k/2}{}}}\mathrm{cos}\left(2\pi {\displaystyle \frac{jl}{N}}\right).`$ (10)
In the previous exponential form the rhs. is a summation for a geometrical series, therefore

$`\lambda _l={\displaystyle \frac{\mathrm{sin}[(k+1)l\pi /N]}{\mathrm{sin}(l\pi /N)}}1.`$ (11)
In the $`N\mathrm{}`$ limit, this converges to

$`\lambda (x)={\displaystyle \frac{\mathrm{sin}[(k+1)x]}{\mathrm{sin}(x)}}1,`$ (12)
where $`x`$ is evenly distributed in the interval $`[0,\pi [`$.

#### 2 Singularities of the spectral density

The spectral density is singular in $`\lambda =\lambda (x)`$, if and only if $`\frac{d\lambda }{dx}(x)=0`$, which is equivalent to

$`(k+1)\mathrm{tan}(x)=\mathrm{tan}[(k+1)x].`$ (13)
Since $`k`$ is an even number, both this equation and Eq. (12) are invariant under the transformation $`x\pi x`$, therefore only the $`xϵ[0,\pi /2]`$ solutions will give different $`\lambda `$ values. If $`k=10`$ (see Fig. 3), Eq. (13) has $`k/2+1=6`$ solutions in $`[0,\pi /2]`$, which are $`x`$ $`=`$ $`0`$, $`0.410`$, $`0.704`$, $`0.994`$, $`1.28`$ and $`\pi /2`$. Therefore – according to Eq. (12) – in the $`N\mathrm{}`$ limit, the spectral density will be singular in the following points:

$`\lambda _i=3.46,2.19,2,\mathrm{\hspace{0.17em}0.043},\mathrm{\hspace{0.17em}0.536},\mathrm{and}k=10.`$ (14)
### B Crossover in the growth rate of the scale-free graph’s principal eigenvalue

The largest eigenvalue is influenced only by the longest row vector if and only if the two longest row vectors are almost orthogonal:

$`\stackrel{}{v}_1\stackrel{}{v}_2|\stackrel{}{v}_1||\stackrel{}{v}_2|.`$ (15)
For $`m>1`$, the lhs. of (15) is the number of simultaneous $`1`$’s in the two longest row vectors, and the rhs. can be approximated with $`|\stackrel{}{v}_1|^2=k_1`$, the largest degree of the graph. It is known, that for large $`j`$ ($`j>i`$), the $`j`$th vertex will be connected to vertex $`i`$ with probability $`P_{ij}=m/(2\sqrt{ij})`$. Thus, we can write (15) in the following forms:

$`{\displaystyle \underset{t=1}{\overset{N}{}}}P_{1t}^2{\displaystyle \underset{t=1}{\overset{N}{}}}P_{1t},`$ (16)
or

$`{\displaystyle \frac{\sqrt{N_\mathrm{c}}}{\mathrm{ln}N_\mathrm{c}}}{\displaystyle \frac{m}{4}}`$ (17)
where $`N_\mathrm{c}`$ is the critical system size.


format: text

text:
# Spin chains and combinatorics: twisted boundary conditions

## Abstract

The finite XXZ Heisenberg spin chain with twisted boundary conditions was considered. For the case of even number of sites $`N`$, anisotropy parameter $`1/2`$ and twisting angle $`2\pi /3`$ the Hamiltonian of the system possesses an eigenvalue $`3N/2`$. The explicit form of the corresponding eigenvector was found for $`N12`$. Conjecturing that this vector is the ground state of the system we made and verified several conjectures related to the norm of the ground state vector, its component with maximal absolute value and some correlation functions, which have combinatorial nature. In particular, the squared norm of the ground state vector is probably coincides with the number of half-turn symmetric alternating sign $`N\times N`$ matrices.

The Hamiltonian of the periodic XXZ Heisenberg spin chain with the number of sites $`N`$ and the anisotropy parameter $`\mathrm{\Delta }`$ has the form

$$H^{(N)}=\underset{i=1}{\overset{N}{}}\left[\sigma _i^x\sigma _{i+1}^x+\sigma _i^y\sigma _{i+1}^y+\mathrm{\Delta }\sigma _i^z\sigma _{i+1}^z\right].$$

where we assume that

$$\sigma _{N+1}^x=\sigma _1^x,\sigma _{N+1}^y=\sigma _1^y,\sigma _{N+1}^z=\sigma _1^z.$$

In terms of the operators $`\sigma _i^+=(\sigma _i^x+\mathrm{i}\sigma _i^y)/2`$, $`\sigma _i^{}=(\sigma _i^x\mathrm{i}\sigma _i^y)/2`$ and $`\sigma _i^z`$ one has the following expression for the Hamiltonian

$$H^{(N)}=2\underset{i=1}{\overset{N}{}}\left[\sigma _i^+\sigma _{i+1}^{}+\sigma _i^{}\sigma _{i+1}^+\right]\mathrm{\Delta }\underset{i=1}{\overset{N}{}}\sigma _i^z\sigma _{i+1}^z,$$
(1)
and the boundary conditions take the form

$$\sigma _{N+1}^+=\sigma _1^+,\sigma _{N+1}^{}=\sigma _1^{},\sigma _{N+1}^z=\sigma _1^z.$$
(2)
It was proved in paper  that for $`\mathrm{\Delta }=1/2`$ and for an odd $`N`$ the Hamiltonian $`H^{(N)}`$ has the eigenvalue $`3N/2`$, see also . In our paper  we analysed the explicit expressions for the corresponding eigenvectors and made the conjecture, which was verified for $`N17`$, that the eigenvalue $`3N/2`$ corresponds to the ground state of the system. This is in agreement with the earlier numerical results reported in paper . In the same paper  we made and verified for $`N17`$ a lot of conjectures concerning the norm of the ground state, the largest components of the ground state vector and some correlation functions. These conjectures have a combinatorial character and are related, in particular, to the number of alternating sign matrices (see, for example, ).

Several variations of the XXZ Heisenberg chain with a special ground state energy are also known . Very recently Batchelor, de Gier and Nienhuis considered these variations along with the corresponding O$`(n)`$ loop model at $`n=1`$ . They found the generalization for two of our conjectures for these systems and extended a list of related combinatorial objects.

First of all the authors of paper  reported about XXZ Heisenberg chain with twisted boundary conditions. The Hamiltonian in this case has again form (1), but instead of (2) one uses the boundary conditions

$$\sigma _{N+1}^+=\mathrm{e}^{\mathrm{i}\varphi }\sigma _1^+,\sigma _{N+1}^{}=\mathrm{e}^{\mathrm{i}\varphi }\sigma _1^{},\sigma _{N+1}^z=\sigma _1^z.$$

Explicitly one can write

$$H^{(N)}=2\underset{i=1}{\overset{N1}{}}[\sigma _i^+\sigma _{i+1}^{}+\sigma _i^{}\sigma _{i+1}^+]2\mathrm{e}^{\mathrm{i}\phi }\sigma _N^+\sigma _1^{}2\mathrm{e}^{\mathrm{i}\phi }\sigma _N^{}\sigma _1^+\mathrm{\Delta }\underset{i=1}{\overset{N}{}}\sigma _i^z\sigma _{i+1}^z.$$
(3)
In paper  on the base of numerical calculations it was conjectured that for the special value of the twisting angle $`\varphi =2\pi /3`$ and $`\mathrm{\Delta }=1/2`$ the ground state energy is equal to $`3N/2`$. The proof of the existence of such an eigenvalue of the Hamiltonian was given in paper . It was the case that was considered by the authors of paper . They formulated two conjectures concerning the sum of the components of the supposed ground state vector and the sum of its squared components enjoying remarkable combinatorial character.

In the present paper we also consider the case $`\varphi =2\pi /3`$, $`\mathrm{\Delta }=1/2`$ and study the eigenvectors of the Hamiltonian (3) corresponding to the eigenvalue $`3N/2`$ with the main purpose to find conjectures on correlation functions. We start with exposition of the explicit form of the eigenvectors for $`N=2,4,6`$. We hope that this information will help someone for further treating of the problem. As a matter of fact we found explicitly the eigenvectors for $`N=8,10`$ and 12 as well. But they have 70, 252 and 924 components, respectively, so we decided against their exposition. Let us fix the notations.

We use the basis in the state space formed by the normed common eigenvectors of the operators $`\sigma _i^z`$. Here we associate $``$ with the eigenvalue $`+1`$ and $``$ with the eigenvalue $`1`$. Thus, the basis vectors are marked by words consisting of $`N`$ letters $``$ and $``$. The set of such words is denoted $`W^{(N)}`$ and we write for a general state vector

$$|\mathrm{\Psi }=\underset{AW^{(N)}}{}\mathrm{\Psi }_A|A.$$

Furthermore, we use such representation for the operators $`\sigma _i^x`$ and $`\sigma _i^y`$ that their matrix elements coincide with usual Pauli matrices. Finally, the eigenvector of the Hamiltonian (3) with the eigenvalue $`3N/2`$ is denoted $`|\mathrm{\Psi }^{(N)}`$.

The eigenvector $`|\mathrm{\Psi }^{(2)}`$ has the form

$$|\mathrm{\Psi }^{(2)}=\frac{\mathrm{i}}{2}(1+\mathrm{i}\sqrt{3})+\frac{\mathrm{i}}{2}(1\mathrm{i}\sqrt{3}).$$

Let us discuss the used normalization condition. Consider the operator

$$R=\sigma _1^x\mathrm{}\sigma _N^x$$

which reverse the $`z`$-action projection of the spins and the antiunutary operator $`K`$ of complex conjugation acting on a general state vector as

$$K|\mathrm{\Psi }=\underset{A}{}\mathrm{\Psi }_A^{}|A.$$

It is clear that the Hamiltonian (3) commutes with the operator $`R^{}=RK`$111The Hamiltonian (3) also commutes with the product of the operator reversing the direction of the chain and the operator $`K`$ and such symmetry can be also used for the fixation of $`|\mathrm{\Psi }^{(N)}`$.. Moreover, we have $`R^2=1`$. We have found that for even $`N12`$ the eigenvalue $`3N/2`$ is nondegenerate. It follows from these facts that multiplying $`|\mathrm{\Psi }^{(N)}`$ by an appropriate phase factor one can satisfy the condition

$$R^{}|\mathrm{\Psi }^{(N)}=|\mathrm{\Psi }^{(N)}.$$
(4)
This condition is invariant under the multiplication of $`|\mathrm{\Psi }^{(N)}`$ by a real number. To fix this freedom we divide $`|\mathrm{\Psi }^{(N)}`$ by the least absolute value of its components. After this the vector $`|\mathrm{\Psi }^{(N)}`$ is defined up to a sign which is fixed by the requirement that the sum of the components of $`|\mathrm{\Psi }^{(N)}`$ be positive222Note that this sum is real due to the condition (4)..

Under the above normalization conditions the vector $`|\mathrm{\Psi }^{(4)}`$ looks as

$`|\mathrm{\Psi }^{(4)}={\displaystyle \frac{1}{2}}(1`$ $``$ $`\mathrm{i}\sqrt{3})+{\displaystyle \frac{1}{2}}(3\mathrm{i}\sqrt{3})+`$
$`+`$ $`+{\displaystyle \frac{1}{2}}(3+\mathrm{i}\sqrt{3})+{\displaystyle \frac{1}{2}}(1+\mathrm{i}\sqrt{3}).`$
For the case $`N=6`$ we obtain

$`|\mathrm{\Psi }^{(6)}`$ $`=`$ $`\mathrm{i}+{\displaystyle \frac{1}{2}}(\sqrt{3}5\mathrm{i})`$
$`+`$ $`(\sqrt{3}2\mathrm{i})+{\displaystyle \frac{1}{2}}(\sqrt{3}\mathrm{i})+(\sqrt{3}2\mathrm{i})`$
$`+`$ $`{\displaystyle \frac{5}{2}}(\sqrt{3}\mathrm{i})+{\displaystyle \frac{1}{2}}(3\sqrt{3}\mathrm{i})+{\displaystyle \frac{1}{2}}(3\sqrt{3}\mathrm{i})`$
$`+`$ $`{\displaystyle \frac{1}{2}}(3\sqrt{3}+\mathrm{i})+{\displaystyle \frac{1}{2}}(\sqrt{3}+\mathrm{i})+{\displaystyle \frac{1}{2}}(\sqrt{3}\mathrm{i})`$
$`+`$ $`{\displaystyle \frac{1}{2}}(3\sqrt{3}\mathrm{i})+{\displaystyle \frac{1}{2}}(3\sqrt{3}+\mathrm{i})+{\displaystyle \frac{1}{2}}(3\sqrt{3}+\mathrm{i})`$
$`+`$ $`{\displaystyle \frac{5}{2}}(\sqrt{3}+\mathrm{i})+(\sqrt{3}+2\mathrm{i})+{\displaystyle \frac{1}{2}}(\sqrt{3}+\mathrm{i})`$
$`+`$ $`(\sqrt{3}+2\mathrm{i})+{\displaystyle \frac{1}{2}}(\sqrt{3}+5\mathrm{i})+\mathrm{i}.`$
Due to the complex character of the components of the eigenvector $`|\mathrm{\Psi }^{(N)}`$, our data is less than for the periodic case . Nevertheless, we are able to formulate and verify analogues of conjectures given in paper . The first conjecture associates the eigenvalue $`3N/2`$ with the ground state.

###### Conjecture 1

The ground state of Hamiltonian (3) for an even $`N`$, $`\varphi =2\pi /3`$ and $`\mathrm{\Delta }=1/2`$ has the energy $`E=3N/2`$ and the total spin $`z`$-axes projection $`S_z=0`$.

Consider now the maximal and the minimal absolute value of the components of $`|\mathrm{\Psi }^{(N)}`$. Note that the minimal absolute value is taken, in particular, by the ferromagnetic component $`\mathrm{\Psi }_{\mathrm{}\mathrm{}}^{(N)}`$ and the maximal absolute value is taken by the antiferromagnetic component $`\mathrm{\Psi }_{\mathrm{}}^{(N)}`$. For the absolute value of the ratio $`\mathrm{\Psi }_{\mathrm{}}^{(2n)}/\mathrm{\Psi }_{\mathrm{}\mathrm{}}^{(2n)}`$ we have the sequence

$$1,\sqrt{3},5,14\sqrt{3},198,1573\sqrt{3},\mathrm{}$$
(5)
Remind that for the periodic boundary conditions and an odd number of sites we have obtained  that the similar ratio $`\mathrm{\Psi }_{\mathrm{}}^{(2n+1)}/\mathrm{\Psi }_{\mathrm{}\mathrm{}}^{(2n+1)}`$ coincides with the number $`A_n`$ of the alternating sign $`n\times n`$ matrices. The numbers $`A_n`$ are given by the formula

$$A_n=\underset{i=0}{\overset{n1}{}}\frac{(3i+1)!}{(n+i)!}.$$

and the sequence $`A_n`$ goes as

$$1,2,7,42,429,7436,\mathrm{}$$
(6)
Dividing sequence (5) by sequence (6), we come to

$$1,\frac{\sqrt{3}}{2},\frac{5}{7},\frac{\sqrt{3}}{3},\frac{6}{13},\frac{11\sqrt{3}}{413},\mathrm{}$$
(7)
Let us now we pay the special attention to the largest prime divisors in the denominators, 7 for $`n=3`$ and 13 for $`n=5`$. Their locations prompt to consider the auxiliary sequence

$$1,4,47,4710,471013,47101316,\mathrm{}$$
(8)
and to multiply sequence (7) by sequence (8). This leads to

$$1,2\sqrt{3},2^25,\frac{2^357\sqrt{3}}{3},2^4357,2^55711\sqrt{3},\mathrm{}$$
(9)
The largest prime divisors in the numerators hint at another auxiliary sequence

$$1,3,35,357,3579,357911,\mathrm{}$$
(10)
Dividing sequence (9) by sequence (10), we arrive at

$$1,\frac{2}{\sqrt{3}},\left(\frac{2}{\sqrt{3}}\right)^2,\left(\frac{2}{\sqrt{3}}\right)^3,\left(\frac{2}{\sqrt{3}}\right)^4,\left(\frac{2}{\sqrt{3}}\right)^5,\mathrm{}$$

As the result we obtain

###### Conjecture 2

The absolute value of the ratio of a component of $`|\mathrm{\Psi }^{(2n)}`$ with the largest absolute value to the component with the smallest absolute value is given by the formula

$$\left|\mathrm{\Psi }_{\mathrm{}}^{(2n)}/\mathrm{\Psi }_{\mathrm{}\mathrm{}}^{(2n)}\right|=\left(\frac{2}{\sqrt{3}}\right)^{n1}\frac{135\mathrm{}(2n3)(2n1)}{147\mathrm{}(3n5)(3n2)}A_n.$$

Consider now the squared norm of the vector $`|\mathrm{\Psi }^{(N)}`$:

$$𝒩_N^2=\mathrm{\Psi }^{(N)}|\mathrm{\Psi }^{(N)}=\underset{AW^{(N)}}{}|\mathrm{\Psi }_A^{(N)}|^2.$$

Fixing the normalization in the way discussed above we see that the sequence formed by the numbers $`𝒩_{2n}^2`$ goes as

$$2,10,140,5544,622908,198846076,\mathrm{}$$
(11)
For the case of periodic boundary conditions the corresponding sequence of squared norms $`𝒩_{2n+1}^2`$ is

$$1,3,25,588,39204,7422987,\mathrm{}$$
(12)
and we conjectured in paper  that

$$𝒩_{2n+1}^2=\left(\frac{3}{4}\right)^n\left[\frac{25\mathrm{}(3n1)}{13\mathrm{}(2n1)}\right]^2A_n^2.$$

Dividing sequence (11) by the sequence formed by the numbers $`A_n^2`$, we obtain

$$2,\frac{5}{2},\frac{2^25}{7},\frac{211}{7},\frac{2^211}{13},\frac{1117}{2^213},\mathrm{}$$

Using the same method as for conjecture 2, we come to the following analogue of conjecture 3 formulated in paper .

###### Conjecture 3

The squared norm of the vector $`|\mathrm{\Psi }^{(2n)}`$ is given by the formula

$$𝒩_{2n}^2=\frac{25\mathrm{}(3n1)}{14\mathrm{}(3n2)}A_n^2.$$

Let us remark that joining the sequences (11) and (12) we come to the sequence

$$1,2,3,10,25,140,588,5544,39204,622908,7422987,\mathrm{}$$

This is the sequence formed by the number of half-turn symmetric alternating sign $`N\times N`$ matrices $`A_N^{HT}`$ (see  and references therein). Therefore we can combine our conjecture 3 and conjecture 3 of paper  into the statement that $`𝒩_N^2=A_N^{HT}`$.

Batchelor, de Gier and Nienhuis discovered another analogue of the conjecture 3 of paper .

###### Conjecture 3 (Batchelor, de Gier and Nienhuis )

The sum of the squared components of the vector $`|\mathrm{\Psi }^{(2n)}`$ is given by

$$\underset{AW^{(2n)}}{}[\mathrm{\Psi }_A^{(2n)}]^2=A_n^2.$$

The same authors gave the analogue of conjecture 4 from paper .

###### Conjecture 4 (Batchelor, de Gier and Nienhuis )

The sum of the components of the vector $`|\mathrm{\Psi }^{(2n)}`$ is given by

$$\underset{AW^{(2n)}}{}\mathrm{\Psi }_A^{(2n)}=3^{n/2}A_n.$$

Note that Hamiltonian (3) has a modified shift invariance. Indeed, it commutes with the operator

$$T^{}=T\mathrm{exp}\left[\frac{\mathrm{i}\varphi }{2}\left(\sigma _N^z\frac{1}{N}\underset{i=1}{\overset{N}{}}\sigma _i^z\right)\right],$$

where $`T`$ is the operator of the right shift by a site. Our choice of the operator of modified shift invariance is dictated by the requirement $`T^N=1`$. Note that the operator $`T^{}`$ commutes with the operator $`R^{}`$. It follows from these facts and from requirement (4) that in the case under consideration either $`T^{}|\mathrm{\Psi }^{(N)}=|\mathrm{\Psi }^{(N)}`$ or $`T^{}|\mathrm{\Psi }^{(N)}=|\mathrm{\Psi }^{(N)}`$. It appears that for even $`N12`$ one has

$$T^{}|\mathrm{\Psi }^{(N)}=|\mathrm{\Psi }^{(N)}.$$
(13)
It is natural to suppose that it is so for any even $`N`$. It can be deduced from (13) that for any $`k=0,1,\mathrm{}`$

$$\underset{AW^{(N)}}{}[\mathrm{\Psi }_A^{(N)}]^{3(2k+1)}=0.$$

Let us proceed to the consideration of correlation functions. For any operator $`O`$ we denote

$$O_N=\frac{\mathrm{\Psi }^{(N)}|O|\mathrm{\Psi }^{(N)}}{\mathrm{\Psi }^{(N)}|\mathrm{\Psi }^{(N)}}.$$

It is evident that $`\sigma _{i_1}^z\mathrm{}\sigma _{i_k}^z_{2n}`$ does not depend on the order of the operators $`\sigma _i^z`$. From (4) we obtain

$$\sigma _{i_1}^z\mathrm{}\sigma _{i_{2k+1}}^z_{2n}=0$$

for any $`k=0,1,\mathrm{}`$, and the equality (13) implies that

$$\sigma _{i_1}^z\mathrm{}\sigma _{i_k}^z_{2n}=\sigma _{i_1+1}^z\mathrm{}\sigma _{i_k+1}^z_{2n}.$$

To formulate our conjectures on correlation functions introduce the operators

$$\alpha _i=(1+\sigma _i^z)/2.$$

The analogue of conjecture 5 of paper  is

###### Conjecture 5

$$\alpha _i\alpha _{i+2}_{2n}=\frac{71n^419n^2+20}{16(4n^21)^2}.$$

Unlike the periodic case, this conjecture is a simple consequence of

###### Conjecture 6

There is a simple formula for the correlation functions called the Probabilities of Formation of Ferromagnetic String :

$$\frac{\alpha _1\alpha _2\mathrm{}\alpha _{k1}_{2n}}{\alpha _1\alpha _2\mathrm{}\alpha _k_{2n}}=\frac{(2k2)!(2k1)!(2n+k1)!(nk)!}{(k1)!(3k2)!(2nk)!(n+k1)!}.$$

Note that the conjectured formulas for correlation functions have the same thermodynamic limit both for the twisted and periodic boundary conditions .

Acknowledgments The work was supported in part by the Russian Foundation for Basic Research.


text:
# Detecting barriers to transport: A review of different techniques

## 1 Introduction

Transport processes play a crucial role in many natural phenomena. Among the many examples, we just mention particle transport in geophysical flows which is of great interest for atmospheric and oceanic issues. The most natural framework for investigating such processes is the Lagrangian viewpoint, in which the tracer trajectory $`𝐱(t)`$ is advected by a given Eulerian velocity field $`𝐮(𝐱,t)`$ according to the differential equation

$$\frac{d𝐱}{dt}=𝐮(𝐱,t)$$
(1)
Despite of the simple formal relation (1), the problem of connecting the Eulerian properties of $`𝐮(𝐱,t)`$ to the Lagrangian properties of the trajectories $`𝐱(t)`$, and viceversa, is a very difficult task. In the last $`2030`$ years the scenario has become even more complex by the recognition of the ubiquity of Lagrangian chaos (chaotic advection). Even very simple Eulerian fields can generate unpredictable Lagrangian trajectories which are practically indistinguishable from those obtained in a complex, turbulent, flow .

In the following we will restrict our attention to the case of two-dimensional incompressible velocity field $`𝐮(𝐱,t)`$, with $`𝐱=(x,y)`$. Incompressibility is automatically satisfied by introducing the stream function $`\psi =\psi (𝐱,t)`$ and consequently defining the velocity field in terms of partial derivatives as $`𝐮=(\psi _y,\psi _x)`$. The evolution equation (1) then becomes

$$\frac{dx}{dt}=\psi _y,\frac{dy}{dt}=\psi _x.$$
(2)
Formally (2) is a Hamiltonian system with Hamiltonian $`\psi (𝐱,t)`$. Chaotic trajectories $`𝐱(t)`$ typically appear as a consequence of the time dependence of $`\psi `$ .

Many geophysical flows, when observed at sufficiently large scale, are within this class . Moreover time dependence can be often considered a perturbation over a given steady flow, i.e. (2) represents a quasi-integrable Hamiltonian system. It is well known that in quasi-integrable Hamiltonian systems chaos can be quite non uniform in the phase space, due to the presence of the invariant KAM tori with a chaotic layer around them . The presence of these regular islands (also called coherent structures in the context of geophysical flows) is of particular importance for the dispersion process because they act as barriers to transport. The sensitivity of different diagnostics of transport to the presence of barriers will be the main topic of our investigation. In particular, we will consider the Finite Time Lyapunov Exponent (FTLE), the Okubo-Weiss (OW) and Hua-Klein criteria (HK), and the local Finite Scale Lyapunov Exponent (FSLE).

We discuss these different methods by considering two examples. First, we study transport properties in a kinematic meandering jet model, formerly introduced for describing the Gulf stream . Second, we analyze a large number of stratospheric isoentropic trajectories, computed according to (1) from assimilated wind fields, describing Lagrangian motion around the polar vortex. In both situations our results show that the existence of barriers limiting the dispersion across the stream is well described by the FSLE but it is completely missed by the OK and HK criteria, which can only depict the landscape of alternating unstable hyperbolic and stable elliptic points of the flow. The FTLE will be discussed in relation to the OW and HK exponents and to the FSLE.

The remaining of this paper is organized as follows. In Section 2 we introduce and discuss the different diagnostics for characterizing dispersion. Section 3 is devoted to the analysis of the meandering jet and of the polar vortex. Finally, in Section 4 we present some conclusions.

## 2 Characterization of local transport properties

In presence of Lagrangian chaos, two close trajectories typically separate exponentially in time . Thus the natural statistics we adopt for describing chaotic particle spreading is the relative dispersion statistics. Relative separation between two particles $`𝐑(t)=𝐱^{^{}}(t)𝐱(t)`$ evolves according to the velocity difference

$$\frac{d𝐑}{dt}=𝐮(𝐱(t)+𝐑(t),t)𝐮(𝐱(t),t).$$
(3)
As far as particle separation remains much smaller than the typical length scale $`l_E`$ of the velocity field, we can linearize (3) around the trajectory $`𝐱`$ and, for a generic time dependent flow, we expect exponential growth of the separation, i.e.

$$R(t)R(0)e^{\lambda t}$$
(4)
where $`\lambda `$ is the Lagrangian Lyapunov Exponent (LE).

In the opposite limit, $`Rl_E`$, the two particles feel uncorrelated velocity fields and one recovers the standard diffusive regime, i.e.

$$R^2(t)2Dt$$
(5)
where the average is taken over many particle pairs and where $`D`$ is the diffusion coefficient.

It is important to remark that, in most realistic situations, both asymptotic regimes, i.e. very small $`R(t)`$ for (4) and very large $`R(t)`$ for (5) cannot be attained. From one side, particles separation can be not sufficiently small to justify the linearization leading to (4). In the opposite limit, large separations cannot be reached in presence of boundaries at scales comparable with $`l_E`$. As a consequence, the asymptotic quantities as $`\lambda `$ and $`D`$ cannot be computed and a non-asymptotic characterization of transport is needed .

Let us discuss, now, some techniques one can use to characterize local dynamical properties of a system, in particular the relative dispersion rate as a function of the initial position.

### 2.1 The Finite Time Lyapunov Exponent

Let us start from the definition of the Lagrangian Lyapunov Exponent:

$$\lambda =\underset{t\mathrm{}}{lim}\underset{R(0)0}{lim}\frac{1}{t}\mathrm{ln}\frac{R(t)}{R(0)}$$
(6)
In (6) one basically assumes that the linearization of the perturbation $`R(t)`$ on a generic reference trajectory holds for an infinite time. This is correct only if the perturbation can be considered infinitesimal at any time. The characteristic time naturally associated to the LE is known as the predictability time $`T_\lambda =\lambda ^1`$, which is the characteristic time at which one can predict the position of the tracer in the future. The FTLE is obtained by avoiding the limit $`t\mathrm{}`$ in (6). This gives the instantaneous growth rate over a finite interval $`\tau `$ as

$$\gamma _\tau (t)=\underset{R(t)0}{lim}\frac{1}{\tau }\mathrm{ln}\frac{R(t+\tau )}{R(t)}$$
(7)
which, at variance with $`\lambda `$, depends on the initial point $`𝐱(t)`$. The FTLE is distributed around a mean value which is nothing but the LE, $`<\gamma _\tau >=\lambda `$, where the average is computed over a virtually infinite number of $`\tau `$ intervals along the trajectory .

In principle, even at very small $`R(t)`$, one must wait a certain time interval, $`T_w`$, needed to drive the perturbation along the Lyapunov direction . This is necessary if we want to measure intrinsic properties of the system. In presence of many degrees of freedom, the possibility that the waiting time $`T_w`$ is of the same order of the predictability time $`T_\lambda `$ cannot be excluded (see  for a review on the predictability time in extended systems).

Let us now discuss some practical shortcomings arising when we want to analyze realistic situations described by experimental or model data. First, as we already rescaled, in quasi-integrable Hamiltonian systems, different regions in the phase space can display different behaviors . As a consequence one has a non trivial spatial distribution of Lyapunov exponents: zero if the trajectory lies in a regular island, positive if it diffuses across the stochastic layer. In special cases, when for instance structures of the velocity field are “localized” persistent features, at least within time intervals considerably longer than the characteristic Lagrangian time, a more refined description in terms of finite-time Lyapunov exponent can be more appropriate.

In order to measure $`\gamma _\tau (t)`$ at a point $`𝐱(t)`$ one can make use of the following procedure. Backward integration in time for an interval $`T_{}`$ bring the trajectory at the point $`𝐱(tT_{})`$. An infinitesimal perturbation $`\delta 𝐱(tT_{})`$ is switched on and it is integrated forward to $`\delta 𝐱(t)`$. It $`T_{}`$ is sufficiently long, i.e. $`T_{}T_w`$, the perturbation $`\delta 𝐱(t)`$ will be aligned along the Lyapunov eigenvector and a further integration to $`\delta 𝐱(t+\tau )`$ will give the FTLE according to (7) . In general, $`T_w`$ is not known a priori and can vary with the initial conditions. More serious problems arise from the limits of resolution in the knowledge of experimental or simulated Lagrangian data, which can easily disrupt the linear approximation scenario.

When $`R(t)`$ attains finite sizes, i.e. is of the order of the characteristic lengths of the system, the so-called Finite-Scale Lyapunov Exponent (FSLE) give an appropriate description of the intrinsic physical properties of dispersion at different scales of motion. We will discuss this point in Section 2.3.

### 2.2 The Okubo-Weiss and Hua-Klein criteria

In two-dimensional turbulent flows, the stirring properties of initially small material lines are related to the combined effect of eddy and jet features of the velocity field. In cases when continuous velocity fields are given, a popular way used to characterize the local rate of separation of initially close trajectories is the Okubo-Weiss (OW) criterion , based on the computation of the eigenvalues of the velocity gradient tensor. If explicit time dependence cannot be neglected, Hua and Klein (HK)  have proposed a generalization of the OW criterion, based on the computation of the eigenvalue of the acceleration gradient tensor, related to the distribution of the pressure field.

Let us recall the two criteria in the case of 2D incompressible velocity field with Lagrangian evolution given by (2). The evolution of an infinitesimal separation $`𝐑(t)`$ is given in the tangent space as

$$\frac{d𝐑}{dt}=\left(\begin{array}{cc}\psi _{xy}\hfill & \psi _{yy}\hfill \\ \psi _{xx}\hfill & \psi _{xy}\hfill \end{array}\right)𝐑\mathrm{𝐌𝐑},$$
(8)
where the Jacobian $`𝐌`$ has the property $`𝐌^2=\lambda _0\mathrm{𝟏}`$ with $`\lambda _0=det(𝐌)`$. At small times, the solution of (8) is

$$𝐑(t)=\left[\mathrm{𝟏}+𝐌t+\frac{1}{2}𝐍t^2\right]𝐑(0)+O(t^3)$$
(9)
where $`𝐍=\lambda _0\mathrm{𝟏}+d𝐌/dt`$. The Okubo-Weiss criterion consists in computing $`\lambda _0`$, i.e. the product of the eigenvaluesΦΦeigenvalues of $`𝐌`$. We recall that the quantity $`\lambda _0`$ can be written in terms of the square strain $`\sigma ^2`$ and the square vorticity $`\omega ^2`$ as

$$\lambda _0=\frac{1}{4}(\sigma ^2\omega ^2)$$
(10)
If $`\lambda _0`$ is positive, the two eigenvalues of $`𝐌`$ are real, the velocity field is locally hyperbolic and strain overcomes rotation. For negative $`\lambda _0`$, we have imaginary eigenvalues and a predominance of rotation over strain.

Of course, the Okubo-Weiss criterion may not be sufficient to determine the local strain-vorticity balance in a time dependent flow. In this respect, the Hua-Klein criterion, being based on the sign of the largest eigenvalue of the $`𝐍`$ matrix $`\lambda _+=\lambda _0+\lambda _1`$, with

$$\lambda _1=\sqrt{\frac{d\psi _{xy}}{dt}^2\frac{d\psi _{xx}}{dt}\frac{d\psi _{yy}}{dt}}$$
(11)
gives a “more Lagrangian” description. Both of them provide a picture of the distribution of stable elliptic points and unstable hyperbolic points in the flow. Let us observe that in the case of stationary velocity field, $`d𝐌/dt=0`$ and one has $`\lambda _+=\lambda _0`$.

Let us remark the relationship existing between the HK criterion and the FTLE. The “instantaneous” Lyapunov exponent $`\gamma _\tau `$ estimates the growth rate of a typical perturbation within a finite time interval $`\tau `$, after the perturbation has aligned along the most unstable direction, so that it is an intrinsic property shared by all the trajectories (except for a set of zero probability measure). The HK eigenvalue $`\lambda _+`$ estimates the local maximum strain rate, regardless any alignment time of the perturbation, so that, from the Lagrangian point of view, it may be biased by transient behaviors. In other words, measuring $`\lambda _+`$ corresponds to measuring $`\gamma _\tau `$ at very small $`\tau `$ (the integration time step) starting with a perturbation always aligned along the locally most unstable direction. In practical situations it can happen that the HK eigenvalue and the FTLE give similar local descriptions but, from a theoretical point of view, starting with a perturbation along the local most unstable direction is as arbitrary as choosing any other direction: after a transient, the time evolution will drive the (infinitesimal) perturbation definitely along the most unstable Lyapunov eigenvector.

### 2.3 The local Finite Scale Lyapunov Exponent

In most cases of interest, the linear regime during which the exponential growth of the inter-particle distance occurs can be not resolved. Particle spreading is generally observed on large spatial scales, of the order of the characteristic lengths of the system, and appropriate non linear techniques must be employed to quantify relative dispersion rates (see  for a review about non asymptotic properties of transport and mixing).

Let us consider a very small (infinitesimal) initial perturbation $`R(0)`$ on a trajectory $`𝐱`$. For a chaotic system, after the initial transient, $`R(t)`$ typically grows exponentially in time according to (4). We wait until $`R(t)`$ reaches a certain threshold $`\delta _i`$, at a given time $`t_i`$. Let $`𝐱_i`$ the position of the trajectory $`𝐱`$ at time $`t_i`$. At a later time $`t_f`$, $`R(t)`$ will reach a larger threshold $`\delta _f=r\delta _i`$ with an assigned $`r>1`$. We now define the $`r`$amplification time of $`R(t)`$, relatively to the initial position $`𝐱_i`$, as $`\tau (\delta _i,𝐱_i,r)=t_ft_i`$. From this quantity, the local Finite Scale Lyapunov Exponent is defined as:

$$\lambda _r(\delta _i,𝐱_i)=\frac{1}{\tau (\delta _i,𝐱_i,r)}\mathrm{ln}r$$
(12)
The exponent $`\lambda _r(\delta _i,𝐱_i)`$ is a measure of the local amplification rate of a perturbation of size $`\delta _i`$ on a trajectory passing by the point $`𝐱_i`$. In the case of periodic time dependence in the equation of motion, the local FSLE will depend also on the initial phase.

The above described prescription is necessary in order to characterize the perturbation growth as an intrinsic property of the system. Of course, in realistic situations, we have to consider the following problems: the flow is not simply periodic and thus $`\lambda _r`$ may depend explicitly on time. Moreover, we cannot observe distances below a certain threshold because of finite resolution, and thus the initial perturbation cannot be considered infinitesimal. Another important remark is that in practical applications it may not be possible to have information on a uniform distribution of initial conditions. In particular, let us consider a 2-D time dependent velocity field (e.g. the surface circulation of a sea or an isoentropic layer of the stratosphere) and the relative Lagrangian dispersion of trajectories. Usually one cannot set the initial distance between two particles to an arbitrarily small size, and wait until the relative position vector is aligned with the most unstable direction and starts expanding with a typical rate measured by the LE. Finite resolution imposes a lower limit to physically reasonable distance sizes, say $`\delta `$, and we can only hope to take into account all the possible realizations of the local dispersion rate by taking the average of $`\lambda _r(\delta ,𝐱)`$ ($`r>1`$) over a large number of directions around the initial point $`𝐱`$, i.e. on a sphere with radius $`\delta `$. Furthermore, if $`\delta `$ is not very small if compared to the characteristic lengths of the system, the linear regime of instability is already expired. The FSLE, by its nature, is a non linear indicator of trajectory instability so it can measure relative dispersion rates at finite scales of motion. How long can we follow two trajectories so that their FSLE is still meaningful as a local diagnostics ? The answer depends essentially on how much rapidly the velocity field varies in time relatively to the Lagrangian characteristic time. For instance, if $`T_E`$ is the time scale within which the Eulerian structures, e.g. current systems, change their geometrical and physical aspects, local FSLE’s are meaningful only if they are observed on times $`T_E`$, i.e. in an almost frozen field approximation.

Let us now briefly discuss the main objective of the present paper, i.e. the detection of barriers in particle transport. We seed the flow with an uniform distribution of Lagrangian tracers and compute the FSLE for any trajectory according to the prescription given above.

If $`\mathrm{}_E`$ represents some characteristic length in the flow, e.g. the typical size of the eddies around the Gulf Stream in the North Atlantic Ocean or around the polar jet current in the stratosphere, then particles trapped inside vortices or traveling down jet streams may never separates beyond the scale $`\mathrm{}_E`$, giving $`\lambda _r(\delta )=0`$ for $`\delta \mathrm{}_E`$. On the contrary, particles located in the chaotic layer will give a positive FSLE. Thus at an appropriate value of the threshold $`\delta `$, the map of FSLE (12) can be used as an efficient indicator of the presence of barriers in transport.

We will see how a simple periodically perturbed meandering jet can be a significant test to show that a more appropriate diagnostics is required when macroscopic barriers are under investigation. An interesting result related to geophysical data analysis concerns the barrier effects of the jet current of the stratospheric polar vortex (southern hemisphere). This technique is being used also to study local mixing properties in ocean systems.

## 3 Results

### 3.1 Numerical experiments

We first discuss a simple, but not trivial, kinematic model in which a barrier to motion is known to exist for certain values of the parameters. The transition to the barrier-breakdown occurs by variation of some parameters. We use this model to show what kind of information can be extracted from the different techniques previously discussed. The system, formerly introduced as a model of transport across the Gulf Stream , consists a time periodic streamline pattern forming an oscillating meandering (westerly) current with recirculations along its boundaries:

$$\mathrm{\Psi }=tanh\left[\frac{yBcoskx}{\sqrt{1+k^2B^2sin^2kx}}\right]+cy$$
(13)
where $`k`$ is the spatial wave number of the structure, $`c`$ is the retrograde velocity in the “far field”, $`B`$ is the amplitude of the meanders which varies periodically in time as

$$B=B_0+ϵcos(\omega t+\varphi )$$
(14)
The system can be fully mixing, i.e. any portion of the domain is definitely visited from any initial condition, in a certain portion of the parameter space $`(ϵ,\omega )`$ where in particular one has cross stream transport. This system has two separatrices, with a spatial periodic structure, see Fig. 1, coinciding with the borderlines of the current. At very small $`ϵ`$ the chaotic layer is restricted to a limited region around the separatrices, and no cross stream mixing occurs. In order to have large-scale chaotic mixing, i.e. particles jumping across the jet from a northern recirculation to a southern one, and vice-versa, one needs the overlap of the resonances , when $`ϵ`$ and $`\omega `$ are larger than certain critical thresholds.

In Figure 3 we report the OW indicator $`\lambda _0(𝐱)`$ as function of the initial position. Incidentally, the HK indicator $`\lambda _+(𝐱`$ shows no significant differences from $`\lambda _0`$ and it is not shown. It is important to observe that both the OW and HK criteria are not able to detect the existence or not of a dynamical barrier, i.e. the two figures are practically indistinguishable.

The finite time Lyapunov exponent $`\gamma _\tau (𝐱)`$ is shown in Figure 4. The computation is done according to (7) with initial separation $`R(t)=\delta `$ for all the particle pairs, without the “waiting procedure” (i.e. $`T_w=0`$). This is because we want to mimic a realistic situation of data analysis in which arbitrarily small separations cannot be attained.

We note that, although the indicator is able to detect the jet core (the low FTLE value filament inside the chaotic current), the transition between the confined chaos regime and the large scale mixing regime is not observed (compare Figure 4a and Figure 4b). The asymmetry of the FTLE map is due to the dependence of this indicator on the initial phase of the periodic flow.

Let us consider the local FSLE as $`\lambda _r(𝐱)`$, computed on the same trajectories and with the same initial separation of Figure 4. Figure 5 contains the results of the $`\lambda _r`$ maps, before and after the overlap of the resonances, in what we can call the Melnikov  and the Chirikov  regimes, respectively. The amplification factor $`r`$ and the lower threshold $`\delta `$ are chosen such that the upper threshold $`r\delta `$ is of the order of the jet width. In one case, black regions of zero FSLE values, i.e. particle pairs which never reach the upper limiting separation, are located in the jetcore and in the centers of the recirculation: when chaos is still confined in the vicinity of the separatrices, no cross-stream transport is allowed. In the other case, after the separatrix breaking has occurred, no zero FSLE values are present, indicating that particles can spread apart over any distance from any initial condition.

Let us conclude this section by remarking that a Lagrangian diagnostic based on the FSLE shows is major skill in detecting dynamical barriers in the flow.

### 3.2 Geophysical data

We now consider a geophysical example regarding Lagrangian motion on an isoentropic layer (i.e. at constant potential temperature) at a low stratospheric level, in presence of the winterly polar vortex , characterized by a quasi-zonal robust jet stream.

The typical flow pattern is usually represented by means of stereographic maps of isoentropic Potential Vorticity (PV), see Figure 6, and the modulus of its gradient shown in Figure 7. In winter, the stratospheric PV can be considered as a quasi-conserved quantity over a timescale of about 2-3 weeks. It is also widely accepted in literature that the outer border of the polar vortex, usually identified by the maximum horizontal gradient of isoentropic PV , can act as a strong barrier to meridional cross stream transport.

The kinematics is remarkably similar to that of the previously discussed meandering jet model, if we imagine the latter as closed on itself in a circular geometry. The Lagrangian data set we have used for the analysis account of about $`10^4`$ trajectory pairs, initially uniformly distributed over the southern hemisphere on the 475K isoentropic surface (lower stratosphere). Trajectories are calculated by means of the University of L’Aquila Trajectory Model  using analyzed wind, pressure, and temperature fields from the U.K. Meteorological Office (UKMO)  provided by the British Atmospheric Data Centre (BADC). Latitude coverage goes from poles to about tropics. The trajectories run from June 30th, 1997 up to a maximum observation time of 20 days. The initial distance between pair particles is $`R_i10Km`$.

Being the trajectory evolution simulated with wind fields relative to the southern hemisphere winterly season, we are observing a situation of stable polar vortex regime .

The local properties of the particle relative dispersion are obtained by computing the Okubo-Weiss eigenvalue $`\lambda _0`$ and the local FSLE $`\lambda (r)`$. Spatial derivatives of the velocity fields used for the calculations are estimated as finite differences over spatial grid steps of the order of $`100Km`$.

The results are presented in the two-dimensional maps in Figures 8 and 9. A further analysis (not shown) demonstrates that $`\lambda _0`$ does not change substantially in time and, as a consequence, it reproduces the essential features of the HK exponent $`\lambda _+`$. Neither of the exponents, as in the previous example, is able to detect the presence of a barrier to transport.

On the other hand, the FSLE map (see Figure 9) detects the dynamical barrier as the region of vanishing FSLE values. The location of the barrier is in good agreement with the definition of polar vortex border based on geophysical considerations, e.g. the potential vorticity gradient shown in Figure 7 .

## 4 Conclusions

We have discussed several techniques proposed for describing dispersion in two-dimensional flows. In particular our analysis has been focused on the capability of these techniques to detect the presence of barriers to transport. By means of two examples of geophysical relevance, we have shown that Eulerian-based techniques, such as the Okubo-Weiss criterium and its generalization proposed by Hua and Klein, are not sensible to the presence of barriers. The Lagrangian finite-time Lyapunov exponent is, in principle, useful for describing space variations of the chaotic properties, e.g. in a quasi-integrable Hamiltonian system, but it is limited to small-scale properties of dispersion. A recent non-linear Lagrangian diagnostics, based on the Finite Scale Lyapunov Exponent, is found to give the correct description of the presence of large-scale barriers. As final remark, we notice that the OW criterion has been recently shown to give poor information also in the case of fully developed turbulence : the probability distribution function of $`\lambda _0`$, $`P(\lambda _0)`$, for a typical 2D turbulent field is not sensitive to the presence of coherent structures, i.e. $`P(\lambda _0)`$ is the same as for a Gaussian field.

From a general point of view, it is not a surprise that purely Eulerian statistics, such as the OW quantity, are unable to predict the behavior of Lagrangian tracers. The presence of dynamical barriers is a fundamental information about the transport properties of the flow and thus can be considered as a good discriminatory for the diagnostics. It would be interesting to check the performance of the proposed methods on other geophysical flows.

We thank the British Atmospheric Data Center for the UKMO data. We are grateful to B. Joseph, D. Iudicone, B. Legras, R. Santoleri and G. Visconti for useful discussions. This work is partially supported by MURST (contract N. 9908264583).