\section{\label{sec:intro}Introduction}

The effect of inhomogeneity on the critical behaviour of magnetic systems has been considered in various contexts (e.g. disorder, coupling randomness, quasiperiodic structures); in particular, discrete-spin models defined on fractal topologies possess critical properties significantly different and richer than those found for translationally invariant systems \cite{gefen1,gefen2,gefen3,milan,cassi,done}.

The interest in fractal structures is not purely theoretical: many condensed-matter systems display strong nonuniformity on all length scales and can therefore be characterized as fractal objects; examples include the backbone of percolation clusters, aggregates obtained from diffusion-limited growth processes, and absorbent surfaces.

One of the most known fractals is the Sierpinski gasket (SG), which, due to its exact decimability, allows analytical approaches; in particular, by means of renormalization group techniques, it was proved that the Ising model on the SG exhibits phase transition only at zero temperature, while at any finite temperature the system breaks into domains and loses long-range order \cite{gefen1,gefen2,gefen3,donetti}. While this result was found in the thermodynamic limit, at the mesoscopic sizes peculiar
and interesting thermodynamic properties arise \cite{liu}. More precisely, the Ising model defined on a finite SG exhibits critical-like
 features at nonzero, low temperatures and the solution found in the thermodynamic limit turns out to be a poor approximation. This anomalous behavior has been investigated from a thermodynamic point of view and derives from long-range, slowly decaying correlations at low temperatures \cite{luscombe}.

The way such a behavior is reflected by the evolution of spin configurations is an item so far overlooked (even if, for rectangular lattices, the idea of studying cluster dynamics may be traced back to Peierls and Griffith \cite{peierls, griffith}). Given the importance of the dynamics of Ising-like clusters in many research areas, from condensed matter to biological systems \cite{appl1,appl2,appl3}, the definition and the development of proper tools for this kind of analysis would be very useful.
Moreover, it would be particularly intriguing for inhomogeneous substrates, due to the emerging non-trivial thermodynamic behavior; on the other hand, it is just on such structures that the definition of a proper metrization or evolutionary dynamics can be more awkward.

The interest in cluster mobility actually extends to an extremely wide class of models: indeed the Ising model on the SG may be seen as a particular realization of cellular automata on graphs, i.e. discrete time dynamical systems assigning to each node of a graph $\mathbf G$ a value chosen in a alphabet $\mathbb K $, along a rule depending only on a finite neighborhood at previous time
\cite{wolfram}.

In this work we aim to introduce and develop proper algorithms for the quantitative characterization of cluster dynamics on graphs, and we use this approach for the Ising model on a finite SG, meant as a prototype of cellular automata on graphs. The procedure requires a projection of evolving configurations into an appropriate {\sl partition space}, where an information-based metrics (Rohlin distance) and  a method measuring the effective emergence of configurational novelty (reduction process and amplification parameter) may be naturally defined and worked out in order to focus the changing and the stable components of configurations. The  algorithmic implementations of Rohlin distance and related quantities are deeply affected by the topological features of the substrate \cite{nostro}. For instance, in previous implementations designed specifically for automata on regular lattices, the very passage from one to two dimensions yields a much higher order of complexity, and could not be exported on different substrates \cite{parti1,parti2,soc,entro}.
On the contrary, the algorithm developed here can be directly applied to generic cellular automata for which the metric characterization of cluster dynamics gets feasible; a brief description is given in the appendix $B$.

Our investigations on the Ising SG highlight the existence of two ``critical'' temperatures, $T_I$ and $T_{II},$ demarcating three main regimes which recover, both qualitatively and quantitatively, the results of the previous thermodynamic analysis; in addition, within the above mentioned regimes, we obtain a more detailed dynamical and geometrical characterization. In particular, at very low temperatures
($T < T_I$), a long-range order is established, and the few small-sized clusters display poor overlap from one time step to the next one: this makes the distance close to zero, and the reduction ineffective. As the spot sizes start to increase, the slowness of the evolution is still such that both non-similarity and overlap between successive configurations rapidly grow.
At greater temperatures ($T_I < T < T_{II}$), large-scale correlations start to decay, clusters of all sizes appear, and overlap gets easier; then, for $T> T_{II}$, any trace of order has vanished and, even if overlaps are very important, the complexity of magnetic pattern is irreducible. The progression of such different kinds of non-similarity and dynamical overlap is well described by our parameters (distance, amplification and intersection, see below).
We will also evidence the existence of different scales controlling the disappearance of local and correlated order.
Thus, the phenomenology provided by our method, confirms with a deeper geometric insight the peculiar critical-like behavior exhibited by the system

We finally notice that, in view of future extensions to non-equilibrium situations, we will adopt a microcanonical dynamics \cite{ACV2009}, which allows implementations even in the presence of temperature gradients.

In the next two sections, after recalling basic notations on graphs, we review some facts concerning the thermodynamic properties of the Ising model on the Sierpinski gasket (Sec.~\ref{sec:SG}), and the microcanonical dynamics working as evolutionary dynamics (Sec.~\ref{sec:Micro}). Then, we introduce general procedures for cluster identification and reduction (Sec.~\ref{sec:metrica}) and we show our results on the SG (Sec.~\ref{sec:numerics}). Finally, we present our conclusions and perspectives (Sec.~\ref{sec:con}). Technical remarks about partitions on graphs can be found in the appendices.

\section{\label{sec:SG}Ising model on the Sierpinski gasket}

A generic graph $\mathbf{G}$ is mathematically specified by the
pair $\{\Lambda, \Gamma \}$ consisting of a non-empty, countable
set of points $\Lambda$ joined pairwise by a set of links
$\Gamma$.
The cardinality  $|\Lambda| = N$ of $\Lambda$ represents the number of sites making up the graph, i.e. its
volume. From an algebraic point of view, a graph
is completely described by its
adjacency matrix $\mathbf{A}$. Every entry of this off-diagonal, symmetric
matrix corresponds to a pair of sites, and it equals $1$ if and
only if this couple is joined by a link, otherwise it is $0$.

Here we consider the Sierpinski gasket which can be built recursively with the following procedure: the initial state ($\mathbf{G}_0$) is a triangle and the $g$th stage $\mathbf{G}_g$ is obtained joining two of the three external corners of three $\mathbf{G}_{g-1}$ to form a bigger triangle (see Fig.~\ref{fig:sierp}). In this way, the volume of $\mathbf{G}_{g}$ is $3(3^g-1)/2$. The Gasket is obtained as the limit for $g \rightarrow \infty$ of this procedure.

\begin{figure}[tb] \begin{center}
\includegraphics[width=.30\textwidth]{sierp.eps}
\caption{\label{fig:sierp} (Color on line) Sierpinski gasket of generation
$2$ (a), $3$ (b) and $4$ (c), with volume $N=6$, $N=15$ and $N=42$,
respectively.}
\end{center}
\end{figure}

\subsection{\label{ssec:termo_SG}Thermodynamic properties of the Ising model on the SG}
The Ising model on a generic graph $\mathbf{G}=\{\Lambda, \Gamma \}$ is defined associating the spin variable $s_i = \pm 1$ to every site $i$ of the graph, and considering a nearest-neighbours interaction between points $i$ and $j$, such that $A_{ij}=1$. The Hamiltonian is therefore
\begin{equation}
H_{\mathbf{G}}(\mathbf{s},J) = - J \sum_{i=1}^N \sum_{j>i}^N A_{ij} s_i s_j,
\end{equation}
where $\mathbf{s}=\{s_i\}_{i \in \Lambda}$ denotes the magnetic configuration of the system and the coupling $J$ is assumed to be the same for any couple; in the following we will set $J \equiv 1$.

As it is well known, a magnetic model defined on a finite lattice cannot exhibit critical behaviour at nonzero temperatures; critical features can only emerge when the underlying lattice becomes infinitely large, i.e. in the thermodynamic limit. Another necessary condition in order to have a nonzero critical temperature concerns the topology of the (infinite) substrate: for Euclidean structures it was rigorously shown that the dimension must be larger than 1. Analogously, it has been shown that the discrete symmetry Ising model on finitely ramified fractals cannot have a nonzero critical temperature \cite{gefen3,luscombe,bhattacharya}.
However, the critical behaviour of one-dimensional systems and finitely ramified fractals can be markedly different, since for the latter it is further governed by additional geometric aspects such as ramification, lacunarity, and connectivity \cite{gefen1}.

The thermodynamic properties of the Ising model on the SG were studied in details in \cite{gefen1,gefen2,gefen3,milan,liu,luscombe,stosic}, where it was shown that its scale-invariant, fractal structure leads to highly cooperative correlations and, at sufficiently low temperatures, the correlation length $\xi$ becomes extremely large and slowly decaying (compared with the one-dimensional case), so that any system with size smaller than $\xi$ displays long-range order.
Indeed, one can define an apparent magnetic transition temperature $T_1$ as the point where $\langle m^2 \rangle =0.5$; for a system of generation $g$ ($g \gg 1$),
$T_1 \approx 4J / \ln(4g)$, for example, for generation $g=5$ and $g=6$ one finds $T_1 \approx 1.335$ and $T_1 \approx 1.259$, respectively.
Hence, as the system size is enlarged, $T_1$ diminishes slowly; more precisely, being $N=3(3^g+1)/2$ the number of nodes for a gasket of generation $g$, one has $T_1 \sim 1/ \ln[\ln(N)]$.

On the other hand, the specific heat $c$ does not display any anomaly associated with long-range order, yet it exhibits a peak at a temperature $T_2 \approx 2J$  \cite{liu}. Conversely, the ``reduced'' specific heat $\sim ~ c T^2$, (basically the derivative of the magnetic energy per link with respect to the inverse temperature \cite{luscombe}), evidences a qualitative difference between the SG and a one-dimensional system, since for the former it exhibits a peak, while for latter it grows monotonically as the temperature is increased. Indeed, for the SG, as $T$ is increased from small values, the energy decreases rather slowly as a result of a relatively large cost due to the large (bulk) coordination number; at larger temperatures it becomes progressively easier to reduce the energy and this reflects the fact that the large fluctuations begin to develop; finally, since the energy must ultimately vanish, the rate of change gets smaller and smaller \cite{luscombe}.

In conclusion, the SG displays a non-trivial thermodynamic behavior which can be summarized as follows.
At $T< T_1$ a long-range order is established; for $T_1 < T < T_2$ a short-range order is still present with large  fluctuations on all length scales. At larger temperatures $T>T_2$ any trace of order has disappeared and a paramagnetic state is approached.

\section{\label{sec:Micro}Microcanonical Dynamics}
When studying transport properties, microcanonical dynamics are usually chosen as they allow to describe an isolated system, or its isolated bulk, without any assumption on
the equilibrium state between the system and the surrounding.

Here we adopt a recently introduced microcanonical dynamics \cite{ACV2009}, which features a high degree of flexibility, being ergodic in any temperature range and implementable on a generic structure, even in the presence of disorder. On regular lattices (e.g. cylinder, torus) such a dynamics has already been shown to be able to lead the system to thermalized states compatible with those expected from a canonical dynamics and to allow the study of out-of-equilibrium properties \cite{ACV2009,future}.

Although we will focus only on equilibrium regimes, the reason for choosing this dynamics is twofold: First we test its reliability on an inhomogeneous structure; second, we pave the ground for the study of transport properties on such a substrate.

In the following, we briefly resume how it works having in mind as substrate a generic graph $\mathbf{G}$ described by the adjacency matrix $\mathbf{A}$.

For each pair of connected sites, namely each link $i, j$, such that $A_{ij}=1$, besides the magnetic energy $E^m_{ij} = J s_i s_j$, we introduce a local {\it kinetic energy} $E_{ij}>0$ which is, in principle, unbounded.
Now, the dynamical rule proceeds as follows:
\begin{enumerate}
\item Start from a (discrete) distribution of energies $\{E_{ij}\}_{i,j=1,...,N}$;
\item choose randomly a link  $i \sim j$;
\item extract one over the possible four spin-configurations for the couple of sites $i,j$, and evaluate the magnetic energy variation $\Delta E^m_{ij}$ induced by the move;
\item if $\Delta E_{ij}^m \leq 0$, accept the move and increase the link energy $E_{ij}$ of $\Delta E_{ij}^m$. When $\Delta E^m_{ij} >0$, the move is accepted only if $E_{ij} \geq \Delta E_{ij}^m$ and the link energy is consequently decreased of $\Delta E_{ij}^m$; otherwise  the move is not accepted and the link energy is not updated.
\end{enumerate}

It is worth remarking that $\Delta E_{ij}^m$ allows for energy variations occurred on the link  $i \sim j$ as well as on those pertaining to links adjacent to sites $i$ or $j$:
$$
\Delta E_{ij}^m = 2 s_i \delta_i  \sum_{k \in \Lambda} A_{ik} s_k + 2 s_j \delta_j  \sum_{k \in \Lambda} A_{jk} s_k - 4 s_i s_j \delta_i \delta_j,
$$
where $\delta_i = 1$ if the $i$th spin has undergone a spin flip, otherwise it is zero.
It is therefore clear that, due to the discreteness of the system and to the fact that $J$ is constant over all links, both $E^m_{ij}$ and $E_{ij}$ are discrete variables: the former can only assume two different values corresponding to the aligned and non aligned configurations of the adjacent spins $i$ and $j$; the latter can only assume integer values deriving by proper combinations of the pertaining $\Delta E_{ij}^m$.

We also notice that $E_{ij}$ works as an additional degree of freedom and the above dynamics conserves the total energy given by the following Hamiltonian function
\begin{equation}\label{ham}
H_{\mathbf{G}}(\mathbf{s},\{ E_{ij} \})=\sum_{i,j \in \Lambda}A_{ij}  \left( s_i s_j + E_{ij}\right).
\end{equation}

As shown in \cite{ACV2009}, the magnetic and kinetic energies result to be non-correlated: this allows a natural definition of temperature at equilibrium, which, in a very natural way, depends only on the average kinetic energy. In fact, the link energy satisfies the Boltzmann distribution $\exp(-\beta E_{ij})$, and the fitted constant $\beta=1/T$ just corresponds to the expected inverse temperature of the system.

Finally, we stress that the possible coupling with thermostats set at a temperature $\bar{T}$ can be realized straightforwardly by selecting a subset of links $\Gamma' \subset \Gamma$ (or, analogously a subset of nodes) and by extracting the pertaining kinetic energies according to the Boltzmann distribution, being $\beta = 1/ \bar{T}$  \cite{ACV2009}.

\section{Metrization}\label{sec:metrica}
As foresaid, an interesting characterization of the Ising model on the SG may be accomplished by a configurational analysis defined in the wider context of Cellular Automata on graphs.
Such an analysis can be realized by a particular metrization referring non directly to the configuration space, whose points are the states $\mathbf a, \mathbf b,\mathbf c,.. $ of the system, but to a peculiar {\sl partition space} containing - among other elements - the cluster distributions of the system. The mathematical framework, which  is summarized in the Appendix $A.1$, requires that the graph is endowed with the structure of a probability space. Precisely, we consider the triple $(\mathbf G,\mathcal G, \mu ) $, where the measure $\mu$ on the subset algebra $\mathcal G $ of the gasket $\mathbf G$ is simply given by the normalized number of nodes in every subset. A {\sl configuration} (or {\sl state}) on $\mathbf G$ is a function assigning to each node a value in an alphabet
$\mathbb K $. The set of all possible configurations will be denoted as $\mathcal C (\mathbf G)$, the configuration space.
Since we consider the Ising model on $\mathbf G $, the alphabet is binary, but all we are going to say is independent of the number  $ \mid \mathbb K \mid $. The adjacency matrix, combined with the list of values on nodes, allows an easy definition
of clusters on an arbitrary structure (the procedure is identical to the recognization of connected subset in graph colouring): two homogeneous nodes belong to the same cluster if they are connected through a path of nodes sharing the same value. Thus, every state $\mathbf a$ determines in a natural way  a {\sl partition}  of $\mathbf G $, i.e. an exhaustive collection $\alpha = \Phi(\mathbf a )  \equiv \{ A_1, A_2,...,A_n \} $ of disjoint subsets $ A_i$, each connected and homogeneous, commonly called {\sl atoms} of the partition. The set of all partitions of $\mathbf G$ constitutes the {\sl partition space} $\mathcal Z \equiv  \mathcal Z (\mathbf G )$. The application $ \Phi : \mathcal C (\mathbf G) \to \mathcal Z (\mathbf G )$ is many-to-one, since, for instance, permutations in $\mathbb K $ produce different states but the same partition. Obviously, $\mathcal Z (\mathbf G ) $  contains much more partitions than those derived from clusterization, e.g. since, in general, atoms do not require to be connected sets, as clusters are. In the present case,  $\mathcal Z (\mathbf G)$ is discrete and finite because $\mathbf G$ is such, but the formalism applies in abstract probability spaces (see \cite{bill,AA,rohlin,sinai,rohlin2}).

Clearly, when a dynamics is defined on the graph, this determines configuration orbits  $\{ \mathbf a(t_n) \} $ starting from any initial state $ \mathbf a(t_0) $, and the corresponding partition orbits $\{ \alpha(t_n) \}$, where $ \alpha(t_n)  \equiv  \Phi (\mathbf a(t_n) )$

Basic operations between two arbitrary partitions $\alpha$ and $\beta$  are the minimal common multiple
$ \gamma = \alpha \vee  \beta$, and the maximal common factor (m.c.f.) or {\it intersection} $\sigma = \alpha \wedge \beta$  (see Fig.~\ref{fig:example} and Appendix $A.1$ for details).
The entropy $H({\sigma})$ of the intersection $\sigma$, when calculated between partitions at next steps along an orbit, is an index of the relevance of the (instantaneous) non evolving part.

The metrization of the partition spaces is  based on the Rohlin distance, which describes the non-similarity of two arbitrary partitions $\alpha$ and $\beta$. This distance, requiring the Shannon conditional entropy $H(\alpha | \beta) $ of measurable partitions (see  Eq.~\ref{condiz}), is given by the functional of Eq.~\ref{rohlin}, we anticipate here:
\begin{equation}
d_R (\alpha,\beta)
=H(\alpha|\beta)+H(\beta|\alpha).
\end{equation}
In order to amplify non similarity, in the Appendix $A.2$ we present also a method, referred to as ``reduction process'' and denoted $\pi$, which acts on couples of partitions and uses both operations $\vee$ and $\wedge$. More precisely, given two partitions, say $\alpha$ and $\beta$, their reduction is obtained by first defining their intersection $\sigma = \alpha \wedge \beta$ and by keeping from both partitions only those subfactors $ \widehat{\alpha}_k$ and $ \widehat{\beta}_k$, prime with $\sigma$, namely such that $\widehat{\alpha}_k \wedge \sigma = \widehat{\beta}_k \wedge \sigma = \nu$. Then, the reduced partitions are given by, $\widehat{\alpha} = \vee_j \widehat{\alpha}_j$ and $\widehat{\beta} = \vee_j \widehat{\beta}_j$ respectively (see Fig.~\ref{fig:All}, lower panel). The process $\pi$ gives evidence of the essentially different sub-partitions of any couple in $ \mathcal Z (\mathbf G )$, and therefore amplifies their distance: $d_R(\alpha ,\beta) \leq d_R(\pi(\alpha ,\beta))$. Hence, by comparing the distance between reduced and non-reduced couples, it is possible to introduce an {\sl amplification ratio}, that is $d_R( \pi[\alpha, \beta])/d_R( \alpha, \beta)$, which provides further information about the cluster distribution and mobility. However, all this analysis has to be performed in correlation with other observables. As explained in Appendix $A.3$, the reduction process is effective, namely gives rise to a large amplification ratio, whenever one of the two partitions, say $\alpha$, displays at least one cluster which in $\beta$ is exactly decomposed into smaller ones. Now, in the case under study $\alpha$ and $\beta$ are partitions defined by cluster configurations at two successive steps, and the existence of a common cluster is furthered by the special topology of the SG: clusters corresponding to subgraphs which are (combinations of) gaskets of generation $<g$ are rather stable (the border is made by two vertices only), nonetheless internal fluctuations may occur and hence decompose the cluster itself.

It is worth underlining that while the    $ \gamma = \alpha \vee  \beta$ operation is rather trivial, the $\sigma = \alpha \wedge \beta$ operation and the reduction process $\pi$ are quite tricky. Reduction in particular is the main algorithmical obstacle in handling large graphs ($g > 6$).

Of course, the Rohlin distance is deeply different from the well known Hamming distance $d_H$ in the configuration space, i.e.
\begin{equation} \label{eq:dH}
d_H (\mathbf a,\mathbf b) = {1\over N}\sum_1^{N} \rho (a_i,b_i),
\end{equation}
where $a_i$ and $b_i$ are the values of the $i-$th node, $\rho$ is a distance functional in the alphabet. The simplest distance in $ \mathbb K $ is $\rho(x,y)= 1 - \delta_{xy}$, leading to Hamming distance 0.5 for purely random configurations. Looking for instance to Fig.~\ref{fig:All} (upper panel), we would get the maximal Hamming distance for the configurations $\mathbf{a}$ and $\mathbf{b}$,  and a minor distance for  $\mathbf{a}$ and $\mathbf{c}$,  while the Rohlin metrics on the corresponding partitions gives null distance in the former case and a high distance in the latter.
In other terms, $d_R$ and $d_H$ are  deeply different in principle
not only because they refer to different objects (defined in  $\mathcal Z(\mathbf G) $ and $\mathcal C (\mathbf G)$, respectively), but also because the former has do to with mutual distribution of clusters, which could involve geometrical features and long range correlations, while the latter is the sum of strictly local differences. Therefore, even if the correspondence $\Phi : \mathcal C \to \mathcal Z $ is many-to-one, the loss of information should be compensated by the fact that by $d_R$ we get a global estimate on cluster distributions, instead of a bare sum of uncorrelated differences.
Finally, more details concerning the implementation of the algorithm for measure of the above mentioned metric observables can be found in Appendix $A.4$.

\section{Numerical Experiments and Results}\label{sec:numerics}

We now focus on the Ising model on the SG, with the microcanonical dynamics described in Sec.~\ref{sec:Micro}.
The equilibrium regime is ensured  by coupling the system to thermostats set at the same  temperature; in this particular case the simplest way is to consider, for a given generation $g$, the external triangle, i.e the $3 \times 2^{g-1}$ links defining the perimeter of the gasket; if then we exclude the six angular links, we get $3$ separated thermostats.
This could be useful in the future, as it consistutes a quite simple way of imposing temperature gradients to the system. However, it should be noticed that, due to the fact that the contact between the thermostats and the system gets vanishingly small with respect to the bulk as $g$ is increased, the time for thermalization is expected to grow with the size of the system. Anyway, as mentioned, we are now interested only in the equilibrium behavior.

Before proceeding, we underline that our checks strongly confirm that the thermalized states reached by the microcanonical dynamics are consistent with those expected from a canonical dynamics, e.g. based on the Metropolis algorithm. In particular, we verified that, for a given temperature, macroscopic observables like the magnetization and the energy measured with the two kinds of dynamics are indistinguishable.

For any given temperature $T$, the geometric observables have been calculated as time series starting after a thermalization time, and lasting an observation time $t_{\mathrm{max}}$, where time is measured in units of Monte Carlo (MC) steps; i.e. $N$ elementary moves.
In particular, we consider finite segments of trajectories $\{ \mathbf a(t_k) \} $ in $\mathcal C$ or $\{ \alpha(t_k) \equiv \Phi (\mathbf a(t_k)  )\} $ in $\mathcal Z$, for $t_k = 0,1,2,...,t_{\mathrm{max}}$, as well as the related
intersections  $\sigma(t_k) = \alpha(t_k) \wedge \alpha(t_{k+1})$ and the couples of reduced partitions at successive time steps
$(\widehat \alpha ,\widehat \beta ) = \pi(\alpha ,\beta ) $.
From such equilibrium trajectories we obtain segments of time series for the following  quantities:

\begin{enumerate}
\item the entropy $H(\cdot)$, applied either to the orbits $\{ \alpha(t_k)\}$ or $\{ \sigma(t_k)\}$;
\item the Rohlin distance $d_R(t_k) \equiv  d_R( \alpha(t_{k-1}), \alpha(t_k))$, (see Eq.~\ref{rohlin});
\item the amplified Rohlin distance $ \widehat{d_R}(t_k) \equiv d_R( \pi[\alpha(t_{k-1}), \alpha(t_k)])$;
\item the amplification ratio $\widehat{d_R}(t_k)  / d_R(t_k)$;
\item the Hamming distance $d_H(t_k) = d_H(\mathbf{a}(t_{k-1}),\mathbf{a}(t_k))$, (see Eqs.~\ref{eq:dH} and \ref{hamming}).
\end{enumerate}

Here we focus to the case of unitary time steps ($t_{k+1}- t_k$ equals one MC step), leaving the study of the role of the time gap to future investigations.

After setting the experimental parameters (size, thermalization time, confidence length of trajectories $t_{\mathrm{max}}$), we calculate time averages, variances etc. for each of the series above. Results found for different choices of size (we especially focused on gaskets of generation $g=4$, $g=5$ and $g=6$, corresponding to $N=120 $, $N=366$ and $N=1095$ sites respectively) and of $t_{\mathrm{max}}$ are qualitatively in very good agreement. Moreover, to approach the thermalized state, the initial configuration is taken ferromagnetic; this minimizes the likelihood of pinning effects during the evolution \cite{ale}. In the following we report only the essential information, dropping redundant numerical outputs.

In general, the observables analyzed highlight the existence of different regimes, demarcated by remarkable temperatures $T_I\approx 1.3$ and $T_{II} \approx 2.0$. More precisely, $T_I$ corresponds to
the flex in the Rohlin distance and to the peak in the amplification ratio, while $T_{II}$ corresponds to the peak in the variance of Rohlin distance and to the crossover in the intersection entropy (see Figs.~$2-4$). Interestingly, $T_I$ and  $T_{II}$ recover the ``critical'' temperatures $T_1$ and $T_2$
evidenced by thermodynamic analysis (see Sec.~$2$ and \cite{gefen1,gefen2,liu}). Indeed, consistently with thermodynamic results, the highlighted regimes correspond to a long-range order region and to a disordered  region with a critical-like transition region in between. More precisely:

\begin{figure}\label{fig:d_r}
\includegraphics[width=1.0\textwidth]{Rohlin.eps}
\caption{(Color on line) Mean (left panel)  and variance (right panel) for Rohlin's distance $d_R$ as a function of the temperature; three different sizes are compared.}
\end{figure}

\begin{figure}\label{fig:d_h}
\includegraphics[width=1.0\textwidth]{Hamming.eps}
\caption{(Color on line) Mean (left panel)  and variance (right panel) for Hamming's distance $d_H$ as a function of the temperature; three different sizes are compared.}
\end{figure}

\begin{itemize}
\item for $ T < T_{I} $, $\langle m^2 \rangle$ is close to $1$ (see \cite{liu}), while the distances $d_H$ and $d_R$ are approximately $0$ (see Figs.~$2$ and $3$); in fact, at such small temperatures a ferromagnetic order is established over large length-scales: clusters are constituted by few single spots in the large ``sea'' of equally oriented spins and  spin-flips are rather unlikely to happen. So, configurations - and partitions - at consecutive steps
differ for such small spots that distances are  extremely small. Moreover, during a MC step, the rare spin flips occurring are yet able to change the atoms in such a way that overlaps between consecutive partitions are quite improbable. This inhibits the reduction, and the amplification is close to 1.
For $T$ approaching $T_I$,  such spots get larger, but, due to the slowness of the evolution at low temperature, their borders can remain sufficiently unchanged for several steps. Hence, the emergence of spots within such clusters allows the reduction, processing couples of consecutive partitions, to get more effective, and the time average of the amplification is manifestly increased. The growth of the spots can be retrieved by the intersection entropy increase, as shown in the inset of Fig.~$4$ (left panel).

\item for $T \approx T_I$, $\langle m^2 \rangle$ is close to $0.5$ and both $d_H$ and $d_R$ start to be significantly larger than $0$. This  is a consequence of the fact that spin flips are getting more frequent. More interestingly, the amplification ratio reaches a maximum.
We have seen that, at the middle of the previous regime, two conditions cooperated to start the growth of the amplification ratio: non empty intersection $\sigma$, and the fact that at consecutive steps there exist uniform large clusters which are decomposed internally, yielding an effective reduction (as explained in Appendix $A.3$). Here, the peak in the amplification ratio means that cluster sizes and the speed of the dynamics optimally fit such conditions; moreover, fixed clusters are more likely. The exact determination of the peak temperature is tricky, due to the growing complexity of the configuration dynamics, but we argue that $T_I$
(or $T_1$) is a good approximation. This behavior is consistent with the apparent transition occurring at $T_1$ as a long-range order breakdown.

\item for $T_I < T < T_{II}$,  $\langle m^2 \rangle$ is approaching zero. Clusters get more intricate due to the growing temperature, and  $H({\sigma})$, which measures the relevance of overlapping between successive configurations, exhibits a rapid growth (see Fig.~$4$, left panel). The coexistence of fragmentation  of large clusters, which is a signature of decay for long range correlations, and overlapping may be related to ``critical slowing down'' effects \cite{CSD}.  Clearly, fragmented clusters have a higher probability to overlap, but are unlikely to include small spots; As a consequence the amplification ratio is still larger than $1$, though rapidly decreasing.
Note that $H({\sigma})$ is rather  far from saturation (namely, uncorrelated chaos) and this suggests the persistence of
a local short-length order.
\item for $ T \sim T_{II}$, the amplification ratio is practically 1: such a breakdown of the reduction, as forementioned,  has a completely different meaning than in the case of very low temperatures, where overlapping was insignificant for the smallness of
sparse clusters in the ``big sea'' of dominant magnetic orientation. Here, due to the smallness of the $\sigma$ atoms, it is unlikely that a sufficiently large cluster is internally decomposed at the next step (see also the Appendix $A.3$).

At this temperature the variance of $d_R$ has a maximum (whose value scales like $1/\sqrt{N}$), and this constitutes a signature that, along the trajectory, the fluctuation  of the distance is particularly important. Indeed, the complexity of clusters shape and their mobility can give rise to wide fluctuations of distances in time.
Also, $H({\sigma})$ exhibits a crossover: from $T_{II}$, the growth is due only to fragmentation, while the ``critical
slowing down'' is over. Both the persistent growth of $H({\sigma})$ and the value of $d_H $ (neatly below $0.5$) indicate that the complete chaos is far from being established at this
temperature and for these sizes. Moreover, the large fluctuations make the system more susceptible to spin-flips, i.e. energy changes, and this is consistent with the peak in the specific heat.

\item for $T > T_{II}$ , the Rohlin variance decreases roughly as $\sim 1/N$ toward an asymptotic value corresponding to the uncorrelated chaos, whose onset may be recognized by $d_H \simeq 0.5$. We argue that the fragmentation and the mobility of clusters stabilize the behavior of the time series for global quantities like distances.

\begin{figure}\label{fig:ampl}
\includegraphics[width=1.0\textwidth]{Hsigma2.eps}
\caption{(Color on line) Intersection entropy $H(\sigma)$ (left panel) and amplification ratio $\widetilde{d_R}/d_R$ (right panel) as a function of the temperature. For the latter, accurate numerical data are presented only for generation $g=5$, due to time consuming calculations; nonetheless we have checked that for other values of $g$ the qualitative behavior is robust. The inset of the left panel shows a magnification of the departure from zero of the entropy; notice the semilogarithmic scale.}
\end{figure}

We note however that the maximum for the $d_H$-variance occurs at $T\approx 3$, well above the corresponding maximum of the
$d_R$-variance at $T\approx 2$. We deduce that in this temperature interval big clusters of all sizes, typical of the previous fractal configurations, have started their fragmentation, so that the global shape diversity between subsequent configurations is
 stabilizing and this explains why the $d_R$-variance decreases. However, the fragmentation process is still slow enough that the local matching between different or equal spins at next steps is highly unstable in time, allowing for the growth of the $d_H$-variance. Indeed, spins belonging to the inner part of a cluster (even if of small size)  result from the dynamics to be more stable than those belonging to its periphery. For $T> 3$ the fragmentation is such that both local overlaps and medium length correlations tend to stabilize the time behavior. There is therefore a temperature interval exhibiting a subtle interplay between correlation length and time stability. In other terms, in the way to the chaos we recognize two time scales: one in terms of global similitude (shape) of clusters along the orbits, which starts to stabilize just at the end of the ``critical slowing down''
($T\approx 2$); the other in terms of local overlaps, whose stabilization requires a higher temperature  ($T\approx 3$). The direct visual inspection could be awkward for such properties, undetectable also from the mean magnetization.

\end{itemize}

\section{Conclusions and perspectives}\label{sec:con}

In this work we have performed a quantitative characterization of cluster dynamics for the Ising model defined on a finite Sierspinski Gasket, whose thermodynamic behavior is known to be non-trivial.
The analysis is based on a set of geometric observables, such as Hamming and Rohlin distances, and on reduction operations among partitions which allow to detect the evolving and the stable components of clusters.

The phenomenology evidenced by previous thermodynamic analysis is qualitatively and quantitatively confirmed by the present metric approach, which provides in addition a geometric characterization of the anomalous behavior of the system. Indeed, we highlight first the existence of two ``critical'' temperatures: $T_I$ corresponding to a peak in the amplification ratio, meaning that
there an abrupt change in the clusters behavior occurs, and $T_{II}$, corresponding to a peak in
the Rohlin variance and other crossovers, evidencing the loss of short range order.
More precisely, the amplification ratio gives detailed information about the development of clusters, starting
from little independent spots (amplification ratio equal to 1) to first overlaps (amplification ratio larger than 1 and growing to
the maximum) while the next phase of decreasing ratio indicates the complex effect of decreasing efficiency of the reduction
due to the fragmentation, up to $T_{II}$ when the fragmentation is such that,
notewithstanding the large intersection, the reduction process is inhibited by the extreme improbability of clusters with a decomposable inner part (see  \ref{ampli}).

Another new information we obtain in the high temperature regime is the existence of  different scales in the destruction of local and correlated order, as evidenced by the Hamming and Rohlin variances, with an anticipated peak for the latter quantity.

Thus, the analysis above enabled us to distinguish different behaviors and phases demarcated by $T_I$ and $T_{II}$ as main milestones; we cannot speak of course of ``phase transitions'', since the onset of these distinct behaviors seems to be a rather smooth process. However, from such a detailed description of the collective motion we get an insight on the interplay between dynamic and statistical features, especially the decay of spatial and temporal correlations.

This work opens a lot of possible extensions and further insights.
The next step is the metric characterization of clusters dynamics in condition of non-equilibrium;
a kind of analysis allowed by the chosen dynamics even in the presence of a disordered coupling pattern \cite{ACV2009}.

In particular, it could be interesting to deepen the effect of the double scale in the approach to uncorrelated chaos on the conductivity. In the same context, another interesting item to explore is the relevance of topology for the cluster diffusion. Indeed, the whole set of operations performed on the SG, from cluster identification to  reduction, immediately applies to any automaton on arbitrary, connected graphs: the process depends only on the adjacency matrix, on a generic alphabet and on a dynamics, working as a proper external engine generating a succession of configurations.
Therefore, by the bare substitution of the adjacency matrix, the algorithm is ready to fit a great variety of statistical models (e.g. Pott's model on lattices of arbitrary dimension or graphs), or other network problems where a dynamics is defined. If necessary, nodes could also be weighted, defining alternative probability measures.
True algorithmic problems and non trivial extensions could only arise from alternative definitions of partitions (assuming e.g. atoms of a different kind with respect to the clusters) or from a different factorization, modifying the reduction process.


\section{Introduction and summary}

The first inflationary model widely discussed in literature, dubbed
$R^2$-inflation \cite{starobinsky}, works in a very economic way,
exploiting one and the same interaction---gravity---to accomplish both
inflation and subsequent reheating. It is tempting to exploit gravity
somewhat further, addressing two other important issues: dark matter
production and generation of the baryon asymmetry of the Universe.

To this end we consider in this work the universal mechanism of particle
production operating in $R^2$-inflationary model: scalaron decay. We
have found that free scalars heavier than 10~keV and fermions heavier
than $10^7$~GeV are forbidden in this scenario, as they would
overclose the Universe. The fermion of $10^7$~GeV is a viable
dark matter candidate in this model, while light scalars could
contribute to the {\it hot} dark matter component at best.

We have also found that with two additional right-handed sterile
neutrinos one can explain the neutrino oscillations (via standard
seesaw mechanism) and baryon asymmetry of the Universe (via standard
non-thermal leptogenesis). Curiously, the amount of baryon
  asymmetry available in the model is strongly constrained from above,
  so that the observed amount is only one order of magnitude below the
  model upper limit.

Both the dark matter fermions and sterile neutrinos are
produced in post-inflationary Universe in scalaron decays. This
suggests that SM, supplemented with $R^2$-term and three right-handed
sterile neutrinos, forms a kind of {\em naturally complete }
theory. This theory explains different phenomena beyond
the SM---inflation, reheating, dark matter, baryon asymmetry of the
Universe---involving {\em one and the same mechanism }
based on the peculiarities of
scalaron interactions with itself and other fields. These phenomena
together with neutrino oscillations are the main observational facts
pointing at incompleteness of the SM. All of them can be explained
within the proposed model.

\section{Gravitational production of Dark Matter in $R^2$-inflation}

We start with the following Lagrangian in the Jordan frame:
\begin{equation}
\label{T1*}
S^{JF}=-\frac{M_P^2}{2}\int \!\! \sqrt{-g}\,d^4x\, \left( R-\frac{R^2}{6\,\mu^2} \right) +
S^{JF}_{matter}\;,
\end{equation}
where we use the reduced Planck mass $M_P$ related to the Planck mass $M_{Pl}$
as $M_P=M_{Pl}/\sqrt{8\pi}=2.4\times 10^{18}$~GeV;
$S^{JF}_{matter}$ includes the action of the Standard Model
and other new fields.
In particular, for free
scalar $\varphi$ and Dirac fermion $\psi$
one has\footnote{These fields are free in the sense of particle
  physics: in both the Jordan (Eqs.~\eqref{T1**} and \eqref{T1+})
and Einstein frames (Eqs.~\eqref{T2*} and \eqref{T2+}),
they are free   at $M_P\to\infty$.}
\begin{align}
\label{T1**}
S^{JF}_\varphi&=\int \!\! \sqrt{-g}\,d^4x\,
\left(
\frac{1}{2} \, g^{\mu\nu}\partial_\mu\varphi\partial_\nu\varphi
-\frac{1}{2} \, m_\varphi^2\varphi^2
\right)\;,\\
\label{T1+}
S^{JF}_\psi&=\int \!\! \sqrt{-g}\,d^4x\,
\left(
i\bar\psi
\hat{\cal D}
\psi
-m_\psi\bar\psi\psi
\right)\;,
\end{align}
where $\hat{\cal D}$ is the covariantly generalized Dirac operator,
see e.g.~\cite{Weinberg3}.
It is convenient to go to the Einstein frame by the conformal
transformation
\[
g_{\mu\nu}\to \tilde{g}_{\mu\nu}=\chi\, g_{\mu\nu}\;, ~~~~~~
\chi= {\rm exp}\left( \sqrt{2/3}\,\phi/M_P\right)\;.
\]
Scalar and fermion fields are rescaled then as
\[
\varphi\to\tilde \varphi = \chi^{-1/2}\,\varphi\;,~~~~~~
\psi\to \tilde\psi = \chi^{-3/4}\,\psi\;,~~~~~~  \hat{\cal D} \to
\hat{\tilde{{\cal D}}} =
\chi^{-1/2}
\, \hat{\cal D}\;.
\]
In the Einstein frame the gravity action takes the
Einstein--Hilbert form, but additional scalar degree of freedom $\phi$
emerges and couples to all matter fields. Thus the original action
\eqref{T1*} transforms into
\begin{equation}
\label{T1++}
S^{EF}=\int \!\!\sqrt{-\tilde g}\,d^4x\,
\left[ -\frac{M_P^2}{2}\,\tilde R + \frac{1}{2}\,
\tilde{g}^{\mu\nu} \partial_\mu \phi \partial_\nu\phi
-\frac{3\,\mu^2 M_P^2}{4}\, \left( 1-\frac{1}{\chi\left( \phi\right)} \right)^2
\right] + S^{EF}_{matter}\;,
\end{equation}
where $S^{EF}_{matter}$ includes interactions with the field
$\phi$. In particular, now the actions \eqref{T1**} and \eqref{T1+} read
\begin{align}
\label{T2*}
S^{EF}_\varphi&=\int \!\! \sqrt{-\tilde g}\,d^4x
\left(
\frac{1}{2} \, \tilde g^{\mu\nu}\partial_\mu\tilde\varphi\partial_\nu\tilde \varphi
-\frac{1}{2\,\chi}\, m_\varphi^2\tilde \varphi^2
+\frac{\tilde\varphi^2}{12\,M_P^2}\, \tilde
g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi
+\frac{\tilde\varphi}{\sqrt{6}\, M_P}\,
\tilde g_{\mu\nu} \partial_\mu \tilde\varphi \partial_\nu\phi
\right)\;,\\
\label{T2+}
S^{EF}_\psi&=\int \!\! \sqrt{-\tilde g}\,d^4x
\left(
i\bar{\tilde \psi}
\hat{\tilde{ {\cal D}}}
\tilde\psi
-\frac{m_\psi}{\sqrt{\chi}} \, \bar{\tilde\psi}\tilde\psi
\right)\;.
\end{align}
At small values of $\phi$ this new degree of freedom decouples from
all other fields and both frames become identical, as in this limit
$R$ is also small.

Cosmology of the homogeneous and isotropic Universe described by the
action \eqref{T1*} has inflationary
stage \cite{starobinsky} at large values of $R$. In the Einstein frame
this stage is realized as large-field inflation in the slow-roll
regime taking place at
super-Planckian values of the field $\phi$ which serves as  the inflaton. The
equivalent scalar mode in the action \eqref{T1*}
was named scalaron \cite{starobinsky} and we use both names in what
follows.

Inflation gives rise to primordial scalar perturbations,
whose amplitude normalization to the observed CMB anisotropy and
large-scale structure yields the estimate~\cite{Faulkner:2006ub}
\[
\mu= 1.3\times 10^{-5}\; M_P\;.
\]
The spectral index of the scalar perturbations and
parameters of the generated tensor perturbations are
 consistent with
observational constraints~\cite{Komatsu:2008hk}.

When the slow-roll conditions get violated, inflation terminates
and the inflaton $\phi$ starts to oscillate rapidly, with frequency
equal to scalaron mass $\mu$. This drives the
Universe expansion like at matter-dominated stage.
The
intermediate stage naturally ends up with inflaton decays into
ordinary particles due to {\em universal } interactions as in
\eqref{T2*}, \eqref{T2+}. In particular, scalaron decay rates into a
pair of sufficiently light scalars and into a pair of sufficiently
light fermions are\footnote{These estimates are in agreement
with similar ones in the Jordan frame originally obtained in
\cite{starobinsky,Vilenkin:1985md}, while \eqref{T4+} and \eqref{T4++}
are two times larger and four times
smaller, respectively, than the estimates in Ref.~\cite{Faulkner:2006ub}.}
\begin{align}
\label{T4+}
\Gamma_{\phi\to \varphi\varphi}&=\frac{\mu^3}{192\pi\,M_P^2}
\;,\\
\label{T4++}
\Gamma_{\phi\to\bar\psi\psi}&=\frac{\mu\,m_\psi^2}{48\pi\,M_P^2}\;,
\end{align}
respectively.
In this {\em universal } way the energy flows from the inflaton to
ordinary particles. There is no any amplification of the energy drain to
the SM particles due to coherent effects (like ones in
Ref.~\cite{Kofman:1994rk}), since the produced particles
interact strongly enough.

Formulas \eqref{T4+}, \eqref{T4++} generally mean that scalar
particles
play the major role in the scalaron decay process.  The dominant contribution
to the decay rate comes from coupling to kinetic term in
Eq.~\eqref{T2*}. There is no similar term in the fermion case, see
Eq.~\eqref{T2+}, because of the conformal invariance\footnote{This
  implies quite a similar situation for a scalar {\it conformally }
  coupled to gravity in the Jordan frame. We do not consider this
  case, but note that the coherent processes at preheating stage we
  alluded to above may be relevant then. }  of the fermion kinetic term.  The
energy transfers to the relativistic scalars by the time
$t_{reh}\gtrsim 1/\sum_s\Gamma_{\phi\to \varphi\varphi}$. We
define the reheating temperature $T_{reh}$ as an effective temperature
of produced relativistic matter at the time of equality between
inflaton and matter energy densities. Then in the
absence of coherent effects mentioned above,
for $N_s$ scalar components
contributing to the inflaton decay the effective
temperature of produced relativistic matter is
\begin{equation}
\label{T5+}
T_{reh}\approx 4.5\times 10^{-2}\times g_*^{-1/4}\cdot
\left( \frac{N_s\,\mu^3}{M_P}\right)^{1/2}\;,
\end{equation}
where $g_*$ denotes the number of relativistic species.
In case of the Standard Model with $g_*=106.75$ and $N_s=4$,
Eq.~\eqref{T5+} gives numerically
\begin{equation}
\label{T5*}
T_{reh}\approx 3.1\times 10^9~{\rm GeV}\;,
\end{equation}
which is well in the region where the gauge interactions of the
Standard Model are in equilibrium. So, the value
\eqref{T5*} is the {\em maximum } temperature of the primordial plasma
or the {\em reheating } temperature, indeed.

With low reheating temperature \eqref{T5*} the post-inflationary
matter-dominated stage lasts long enough, so that the subhorizon scalar
perturbations growing at this stage become nonlinear (for a recent
study see e.g. \cite{Jedamzik:2010dq}). One expects that scalarons
then start to form self-gravitating clumps of linear sizes much smaller than
horizon, so the Universe still expands as at matter domination. The
scalaron overdensity in clumps can be estimated similar to analysis of
dark matter halos and clumps in the
late-time Universe, see e.g. \cite{LB}. The scalaron density in clumps
is not high enough to initiate scalaron scatterings due to
self-interaction. Therefore, relativistic (SM) particles production can
be described as decays of (non-relativistic) scalarons and the estimate of
reheating temperature \eqref{T5*} remains intact.

The decay rate formulas \eqref{T4+}, \eqref{T4++} show that the
scalaron produces ordinary particles in a universal way, so that
all scalars and fermions of the theory will eventually populate the
expanding Universe. With all non-gravitational
interactions between ordinary particles switched off, the resulting
abundances would be determined mostly by spin and, for fermions,
by mass of the particles.

Now let us consider a new field which is either free in the Jordan
frame or couples to other fields very weakly, so that these new
particles never equilibrate in the primordial plasma\footnote{ In
  particular, the Planck-scale suppressed nonrenormalizable couplings
  are allowed.}. If stable at cosmological time scale, this particle
is a good candidate to be dark matter.

Given the discussion above, the particle mass is the only free
parameter and hence its value is fixed by the requirement of
comprising all dark matter whose relative contribution to the present
energy density $\rho_c$ is~\cite{Komatsu:2008hk} $\Omega_{DM}\approx 0.223$.
The mass of dark matter particles $m_{DM}$ and their number density at
present $n_{DM,0}$ are related to  $\rho_c$ and $\Omega_{DM}$ as follows
\begin{equation}
\label{T6+}
m_{DM}=\frac{\Omega_{DM}\,\rho_c}{s_0}\,\frac{s_0}{n_{DM,0}}\;,
\end{equation}
where we introduced the present entropy density $s_0$.

Let us estimate the entropy-to-dark-matter ratio $s/n_{DM}$
by the time of reheating. The dark matter production after inflation
can be described as the decay of non-relativistic scalarons of mass
$m_\phi$ whose number density $n_\phi$ evolves with scale factor
$a=a\left( t\right)$ at $t\gtrsim t_{reh}$ as
\[
n_\phi\left( a\right)= \frac{\mu}{2}\,\phi_{reh}^2\,\left( \frac{a_{reh}}{a}\right)^3\;.
\]
Here $\phi_{reh}$, $a_{reh}$ refer to the values of inflaton field and
scale factor at reheating. We again neglect any coherent effects
related to the new fields and also treat this new decay channel of
scalaron as subdominant one, so the reheating temperature is still given by
Eq.~\eqref{T5*}. Both points are discussed below in due course: we
will see that the
obtained results justify our choice.
Assuming two-body decays of the inflaton to
dark matter particles with decay rate $\Gamma_{\phi\to DM}$ one
writes down the Boltzmann equation for the dark matter density,
\[
\frac{d}{dt}\,\left( n_{DM}\,a^3\right) = 2\, n_\phi\!\left( a\right) \, \Gamma_{\phi\to
  DM} \, a^3\;.
\]
This equation has a solution
\[
n_{DM}\left( t_{reh}\right) =\frac{\rho_{reh}}{\mu}\, \Gamma_{\phi\to DM} \, t_{reh}\;,
\]
where $\rho_{reh}$ is the total energy density at the time of reheating,
  when the corresponding contributions of inflaton and relativistic
species coincide. The energy density is related to the Hubble parameter at
reheating $H_{reh}$ by
\[
\rho_{reh}=2\, g_*\,\frac{\pi^2}{30}\, T^4_{reh}=3\, M_P^2 \,
H_{reh}^2 \;.
\]
Taking for the numerical estimate of the reheating
time
\[
t_{reh}\simeq \frac{1}{\sqrt{3}\,H_{reh}}\;,
\]
we finally obtain for the entropy-to-dark-matter ratio at reheating
\begin{equation}
\label{T8+}
\frac{s}{n_{DM}}\left( T_{reh} \right) =
\frac{2\pi\,\sqrt{g_*} }{3\,\sqrt{15}} \,
\frac{\mu}{\Gamma_{\phi\to DM}}\, \frac{T_{reh}}{M_P}\;.
\end{equation}
Since this ratio remains intact at the further hot stages, we have the
estimate
\[
\frac{s}{n_{DM}}\simeq\frac{s_0}{n_{DM,0}}\;.
\]
Eqs.~\eqref{T8+}, \eqref{T6+} imply the following equation for the mass
of dark matter particles universally produced by the scalaron decay,
\begin{equation}
\label{T8*}
m_{DM}=\frac{\Omega_{DM}\,\rho_c}{s_0}\,
\frac{2\pi\,\sqrt{g_*} }{3\,\sqrt{15}} \,
\frac{\mu}{\Gamma_{\phi\to DM}}\, \frac{T_{reh}}{M_P}\;.
\end{equation}
Then in the cases of scalar \eqref{T4+} and fermion \eqref{T4++} one
obtains the numerical estimates
 \begin{align}
\label{T9+}
m_\varphi&\approx 6.9\;{\rm keV}\!\times\!\left(\! \frac{1.3\!\times\!
  10^{-5}\,M_P}{\mu}\!\right)^{\!\!1/2} \left(\! \frac{N_s}{4}\!\right)^{\!\!1/2}
\left( \frac{g_*}{106.75} \right)^{\!\!1/4} \left(\! \frac{\Omega_{DM}}{0.223}\!\right)
\left(\! \frac{\rho_c}{0.52\!\times\! 10^{-5}\,{\rm GeV}/{\rm
    cm}^3}\!\right), \\
\label{T9*}
m_\psi&\approx1.2\!\times\!10^7\;{\rm GeV}\!\times\! \left( \!\frac{\mu}{1.3\!\times\!
  10^{-5}\,M_P}\!\right)^{\!\!1/2} \left( \!\frac{N_s}{4}\!\right)^{\!\!1/6}
\left( \frac{g_*}{106.75} \right)^{\!\!1/12} \left( \!\frac{\Omega_{DM}}{0.223}\!\right)^{\!\!1/3}
\left( \!\frac{\rho_c}{0.52\!\times\! 10^{-5}\,{\rm GeV}/{\rm
    cm}^3}\!\right)^{\!\!1/3}.
\end{align}
Here we have corrected numerical coefficients (by tens of percent)
to match the results obtained
with numerical integration of the Boltzmann equations which consistently
describe both the energy transfer from scalaron to all relativistic
scalars and the dark matter production.

The dark matter particles are produced highly relativistic with 3-momenta
$\sim \mu/2$, which exceed $T_{reh}$ and hence
momenta of particles in the plasma by about four orders of magnitude,
\[
\mu \sim 10^4\times T_{reh}\;.
\]
So, the free scalars of 10\,keV mass would contribute to {\em hot dark
  matter } component only. This component is certainly subdominant,
hence free scalars cannot solve the dark matter problem in the case of
$R^2$-inflation. The viable choice of dark matter is free fermions of
mass $m_{DM}\approx 10^7$\,GeV. These particles would form {\em cold
  dark matter. }

We end up the discussion of universal dark matter in
$R^2$-inflationary model with several remarks. First,
for the considered case of free boson particles coherent
effects \cite{Kofman:1994rk} during scalaron oscillations could change
the obtained results. However, it only makes worse the situation with
scalar dark matter, and we do not consider these processes.

Second, if at a later stage of the Universe expansion
the entropy gets produced (say, due to decays of nonrelativistic
particles out of equilibrium, etc.), the estimates \eqref{T9+},
\eqref{T9*} have to be corrected to account for the corresponding dilution
factor $r_s=s_{new}/s_{old}$. Thus, r.h.s. of \eqref{T9+} has to
be divided by $r_s$, while r.h.s. of \eqref{T9*} has to be divided
by $r_s^{1/3}$. In realistic models with moderate (if any) entropy
production this does not save scalar dark matter, and only mildly
changes the prediction for the fermion dark matter mass \eqref{T9*}.

Third, in the extensions of particle physics with additional scalars
coupled to the SM fields,
both reheating temperature \eqref{T5+} and would-be scalar dark
matter mass increase. With very large number of new scalars
this makes the scalar dark matter viable.

Fourth, if there are other sources of dark matter in the
model---other stable particles or other mechanisms of out-of-equilibrium
production of heavy fermions---the estimate \eqref{T9*} implies the
{\em upper } limit on the mass of universal dark matter in
$R^2$-inflationary model.

Fifth, in case of the Majorana fermions the mass is larger by a factor of
$2^{1/3}$ as compared to the Dirac case \eqref{T9*}, since the decay
rate of scalaron into the Majorana particles is two times smaller.

Sixth, the estimates \eqref{T9+}, \eqref{T9*} mean that {\em heavier free
particles are forbidden } in models with $R^2$-inflation, otherwise the
Universe would be overclosed. This conclusion is true for the
particles lighter than the scalaron, when these estimates of particle
production rates are applicable. Clearly, production rate of heavier
particles by means of scalaron oscillations is suppressed as compared
to \eqref{T4+}, \eqref{T4++}. The heavier the particles, the stronger
the suppression. Hence, free particles of mass $m\gg \mu$ are
not forbidden in $R^2$-inflation. Moreover, with strongly (most
probably exponentially) tuned value of mass $m\sim\mu$ they can form
viable dark matter. We do not study this situation here.

\section{Leptogenesis in $R^2$-inflation}

The study presented in the previous Section revealed that heavy
free (in the Jordan frame) fermion is a viable dark matter candidate
in $R^2$-inflationary model. These fermions are naturally produced in
the post-inflationary Universe by scalaron decays, so the right amount
of dark matter is achieved with the mass of about $10^7$\,GeV
\eqref{T9*}. One can further
ask whether it is possible to make use of this universal production
mechanism to unravel another cosmological problem which SM fails to
solve: the baryon asymmetry of the Universe. The answer is positive,
and we illustrate it with example of nonthermal leptogenesis via
decays of heavy sterile
neutrinos~\cite{Fukugita:1986hr} universally produced
by scalaron decays. This particular example is strongly motivated by
observed oscillations of active neutrinos, phenomena which also
lack an explanation within SM. Heavy sterile neutrinos give masses to
active neutrinos via seesaw mechanism. This mechanism is responsible for the
hierarchy between the neutrino masses and electroweak scale. In this
Section we show that both the correct neutrino masses and successful
leptogenesis can be realized
in the $R^2$-inflationary model.

The modification of SM we consider consists of two new Majorana
fermions $N_I$, $I=1,2$, which are right singlets with respect
to the SM gauge group. The most general renormalizable lagrangian for
these fermions is
\begin{equation}
\label{c2+}
{\cal L}= i \bar N_I \gamma^\mu\partial_\mu N_I - y_{\alpha I} \bar L_\alpha
N_I \tilde \Phi - \frac{M_I}{2}\bar N_I^cN_I + h.c.\;,
\end{equation}
where $y_{\alpha I}$ are new Yukawa couplings, $\Phi$ is the SM Higgs
doublet and $\tilde \Phi=\epsilon\,\Phi^*$ with $\epsilon$ being
antisymmetric $2\times2$ matrix.

When electroweak symmetry breaks down, $\Phi$ acquires vacuum
expectation value, $\Phi^T=\left( 0, v/\sqrt{2}\right)$, $v=246$~GeV. Then
\eqref{c2+} yields active-sterile mixing in the neutrino mass
matrix. Assuming $yv\ll M_I$, this mixing gives for the active
neutrino masses:
\begin{equation}
\label{c2*}
{m_{\nu}}_{\,\alpha\beta}= - \sum_Iy_{\alpha I}\frac{v^2}{2\,M_I}y_{\beta I}\;.
\end{equation}
Hence, $m_\nu\ll v$, that is the seesaw
mechanism\footnote{Note that with only two seesaw sterile neutrinos we
end up with one out of three active neutrino being massless.}.
The formula
\eqref{c2*} together with atmospheric neutrino mass
splitting, $\Delta m_{atm}^2\simeq 3\times
10^{-3}$~eV$^2$~\cite{Amsler:2008zzb}, implies an order-of-magnitude
upper limit on the {\em lightest sterile neutrino } mass:
\begin{equation}
\label{c3*}
M_{lightest}\lesssim \frac{v^2}{2\sqrt{\Delta m_{atm}^2}}\simeq 5\times
  10^{14}~{\rm GeV}\;.
\end{equation}
The smaller the Yukawa couplings $y_{\alpha I}$ and the heavier the
active neutrinos, the smaller $M_{lightest}$, see Eq.~\eqref{c2*}.

Lagrangian \eqref{c2+} added to SM explicitly violates the lepton
number. The source of CP-violation, which is one of
the Sakharov's conditions~\cite{Sakharov:1967dj}
for successful baryo- and leptogenesis, is the complex Yukawa
couplings $y_{\alpha I}$. The third Sakharov's condition, departure
from thermal equilibrium for lepton number violating processes, is
achieved due to decay of non-relativistic sterile
neutrinos, which is evidently an out-of-equilibrium process.

Let us discuss the production of these heavy sterile neutrinos in the
post-inflationary Universe with oscillating scalaron. Neutrinos are
produced via Yukawa-type coupling to scalaron, see
Eq.~\eqref{T2+}. The production is effective only for neutrinos
lighter than scalaron. Indeed, the amplitude of scalaron oscillation
never exceeds $M_P$. At the same time the scalaron coupling to a
fermion is proportional to the fermion mass, so the latter remains
almost intact. Therefore, an enhancement of heavy fermion production
by Yukawa source due to periodic decrease of fermion effective
mass~\cite{Giudice:1999fb} does not work here.

Thus, we consider the models with sterile neutrinos lighter than
scalaron, $M_I<m_\phi$. By the time of reheating, $t\simeq t_{reh}$,
the ratio of number density of produced in scalaron decays heavy
neutrinos to entropy is given by Eq.~\eqref{T8+}. Substituting
relevant numbers in \eqref{T8+} and \eqref{T4++} (recall that for
Majorana fermions the decay rate is two times lower than for the Dirac
fermions) one obtains
\begin{equation}
\label{c3+}
\frac{n_{N_I}}{s}\left( T_{reh}\right)=2.9\times 10^{-6}\times \left(
\frac{M_I}{5\times 10^{12}~{\rm GeV}}\right)^2\;,
\end{equation}
where we have included a correction of about tens of percent in order to match
the solution of the relevant Boltzmann equations (see comment below
Eqs.~\eqref{T9+}, \eqref{T9*}).

These sterile neutrinos, however, decay quite rapidly due to the same
Yukawa couplings in \eqref{c2+} responsible for the seesaw
mechanism. Their total decay rates are
\[
\Gamma_{N_I}=\frac{M_I}{8\pi}\sum_\alpha \left| y_{\alpha I} \right|^2 \;.
\]
Disregarding any hierarchy in the neutrino Yukawa sector one takes
for the order-of-magnitude numerical estimate of Yukawa
couplings relation \eqref{c2*} with $m^2_\nu\sim \Delta m^2_{atm}$
and finds for the lightest sterile
seesaw neutrino
\[
\Gamma_{N_I}\sim \frac{\sqrt{\Delta m^2_{atm}}}{4\pi}
\frac{M_I^2}{v^2}\;.
\]
Thus, sterile neutrino of mass above $10^{10}$~GeV decays before
reheating\footnote{Eqs.~\eqref{c3+} and \eqref{T5*} mean that sterile
  neutrino decays contribute not much to the Universe reheating.}.
This is a very process which generates lepton asymmetry in
the early Universe. To simplify formulas we assume further, that only
the lightest neutrino, $N_1$, contributes and  $M_1\ll M_2$. Then
the value of lepton asymmetry $\Delta_L=n_L/s$ is
\begin{equation}
\label{c4+}
\Delta_L= \delta_L\cdot \frac{n_{N_1}}{s}\;,
\end{equation}
with~\cite{Fukugita:1986hr} (for reviews see
e.g. \cite{Buchmuller:2005eh})
\begin{equation}
\label{c4*}
\delta_L=-\frac{3\,M_1}{8\,\pi v^2 } \cdot
\frac{1}{\sum_{\alpha,\, I} \left| y_{\alpha I}\right|^2}
\sum_{\alpha\beta} \Im \left( y_{\alpha 1} y_{\beta 1} y_{\alpha 2}^*
\, \frac{v^2}{2M_2}\,y_{\beta 2}^*  \right)
\sim \frac{3M_1\,\sqrt{\Delta m^2_{atm}}}{8\pi v^2}\;,
\end{equation}
where we have used the order-of-magnitude estimate \eqref{c2*}
with $m^2_\nu\sim \Delta m^2_{atm}$. Finally,
from \eqref{c4+}, \eqref{c4*} and \eqref{c3+} we obtain
\begin{equation}
\label{c4++}
\Delta_L\sim 1.5\times 10^{-9}\times
\left( \frac{M_1}{5\times 10^{12}~{\rm GeV}}\right)^3\;.
\end{equation}

The asymmetry is generated mostly at late times just before reheating,
when the ratio of sterile neutrino (if stable) density to entropy
becomes maximum, see Eq.~\eqref{T8+}.
The sterile neutrinos are quite heavy so the lepton
asymmetry is not washed out by inverse decay processes, which are out
of equilibrium. Likewise, another
important lepton violating process, scatterings $\Phi L
\leftrightarrow \Phi \bar L$ through virtual sterile neutrinos (with
cross section $\sigma \sim \frac{m_\nu^2}{8\pi\, v^4}$)
are out-of-equilibrium at the epoch of interest. Thus, lepton
asymmetry \eqref{c4++} transfers later on into baryon asymmetry by sphaleron
processes \cite{Khlebnikov:1996vj},
\[
\Delta_B\simeq\Delta_L/3 \sim 0.5\times 10^{-9}\;.
\]
Given the present baryon asymmetry of the Universe
$\Delta_{B,0}=0.88\times 10^{-10}$~\cite{Komatsu:2008hk} and
the estimates above, we conclude that heavy seesaw
sterile neutrinos of mass $10^{12}$-$10^{13}$\,GeV, produced in scalaron
decays, in $R^2$-inflationary model can be responsible for the baryon
asymmetry of the Universe. Since the sterile neutrino are lighter than
scalaron, the available room is not large, that one
could interpret as {\it an order-of-magnitude prediction of the baryon
  asymmetry, } which is otherwise (and usually) appears as an
accidental number.

The estimates above are an order-of-magnitude only. In the
particular model under discussion with only two seesaw sterile
neutrinos one can relate all the Yukawa couplings to the measured active
neutrino mass differences and mixing angles, so that only few model
parameters remain free. In this way one can refine our estimates.
Likewise the numerical estimates may be changed in models with a particular
hierarchy of Yukawa couplings in the neutrino sector.

Note in passing that second, heavier sterile neutrino $N_2$ if not too
heavy can also
contribute to lepton asymmetry generation, its decay products do not
equilibrate.
There is also some contribution from sterile neutrino decays at
earlier times, when they are relativistic. The inverse processes are
rare because the particle density is small.

\section{Discussion}

To conclude, we have shown that in $R^2$-inflation the scalaron decay
responsible for reheating can also take care of dark matter production
and generation of the baryon asymmetry of the Universe. The latter
occurs through nonthermal leptogenesis by introducing seesaw
neutrinos. Thus, neutrino oscillations are also explained in this
model, which pretends to be {\it phenomenologically self-contained. }

We end up with several comments. First, at relatively low energy
scale,
the model we discussed can be further modified in gravity sector
(attempting to explain the late-time accelerated expansion of the
Universe, as e.g.
in Refs.~\cite{Nojiri:2003ft,Starobinsky:2007hu}) or in particle
physics sector (say, by adding $PQ$-axion \cite{Peccei:1977hh}
to solve the strong CP-problem). All these and similar modifications can be
safely adopted provided the early time cosmology remains intact.

Second, we do not discuss here the gauge hierarchy problem in particle
physics. Certainly, $R^2$-term does not contribute to this problem:
all couplings of scalaron to other particles are suppressed by the
Planck scale. With the absence of new scales in {\em particle physics,}  one
cannot discard that still unknown complete quantum theory at gravity scale
(quantum gravity)
is responsible for the
cancellation of dangerous quantum corrections (if any) to the Higgs boson
mass. In fact, a nontrivial space structure that shows up at Planck scales
could also invalidate the strong CP-problem (see
e.g. \cite{Khlebnikov:1987zg}). Thus, quantum gravity could
solve all problems, indeed. Otherwise, one can think about other
solutions usually considered. In particular, our logic can be adopted
for supersymmetric extensions of the SM and models with axion. In all
these cases free fermion again serves as dark matter candidate, though
it is somewhat heavier than in the case of SM because of larger
number $N_s$ of scalars responsible for reheating and larger
number of degrees of
freedom $g_*$, see Eq.~\eqref{T9*}.

Third comment, which is related to the second one.
The heavy sterile neutrinos (of
mass $M_N$) we used contribute to gauge hierarchy problem making it
worse as compared to the SM. Indeed, their coupling to the SM Higgs
boson gives rise to both divergent and finite, of order $M_N$,
contributions to the Higgs boson mass. This worsening can be avoided if sterile
neutrinos coupled to the Higgs boson have masses of order electroweak
scale or smaller, as e.g. in $\nu$MSM model~\cite{Asaka:2005an}.
Indeed, heavy sterile
neutrinos can be replaced by the light ones in $R^2$-model we
discuss. In that case, however, sterile neutrinos, required for
leptogenesis, should be produced not by the scalaron decay (which, as we saw, requires
large mass to be efficient enough), but via neutrino oscillations in
the early Universe~\cite{Akhmedov:1998qx}.
The $R^2$-model with one heavy fermion (free in
the Jordan frame) and two light seesaw sterile neutrinos can be {\it both
phenomenologically and theoretically self-contained } upto the quantum
gravity effects. At the same time, this modification has also one
more advantage as compared to what we discuss in the main text:
sterile neutrino sector here can be directly
tested~\cite{Gorbunov:2007ak}  in
particle physics experiments.

We thank S. Odintsov, V. Rubakov, M. Shaposhnikov, A. Starobinsky,
I. Tkachev and M. Vasiliev for
valuable discussions.  The work is supported in part by the grant of
the President of the Russian Federation NS-5525.2010.2 (government
contract 02.740.11.0244).  D.G. thanks ITPP EPFL and the organizers of
the long-term workshop in Yukawa Institute YITP-T-10-01 for
hospitality.  The work of D.G.  is supported in part by the Russian
Foundation for Basic Research (grants 08-02-00473a), by FAE program
(government contract $\Pi$520) and by SCOPES program. The work of
A.P. is supported in part by the "Dynasty" Foundation (awarded by the
Scientific board of ICPFM), by FAE program (government contract
$\Pi$2598) and by the grant of the President of the Russian Federation
MK-4317.2009.2.