text:
\section{Introduction}
Born as a theoretical background for thermodynamics, statistical
mechanics provides nowadays a flexible approach to several
scientific problems whose depth and wideness increases
continuously. In fact, in the last decades statistical mechanics
has invaded fields as diverse as spin glasses \cite{MPV}, neural
networks \cite{amit}, protein folding \cite{huang}, immunological
memory \cite{immune}, social networks \cite{PL1}, theoretical
economy \cite{ton} and urban planning \cite{intro2}. As a
consequence, an always increasing need for models and proper
techniques must be fulfilled. Coherently, recently, several models
have been systematically tackled via the smooth cavity field by
the authors, namely the Curie-Weiss model \cite{abarra}, the fully
connected $p$-spin model \cite{abarra3}, the
Sherrington-Kirkpatrick model \cite{barra1}, its diluted counter
part Viana-Bray model \cite{bdsf} and the diluted ferromagnetic
model \cite{ABC}. All these models can just be seen as different
components of a more general class including models based on
binary agents with mean field interactions (Fig. \ref{fig:cubo}).
Now, in order to complete the analysis of the free energies for
the whole class, the X-OR-SAT (of the Random Optimization Theory
\cite{xorsat}) and the diluted ferromagnetic $p$-spin model, are
still missing; this paper is devoted to the study of the latter.

\begin{figure}[tb]\label{fig:cubo}
\includegraphics[height=100mm]{cubo3.pdf}
\caption{\label{fig:cubo} Schematic representation of the connections among different models based on mean
field interactions between variables endowed with discrete symmetry.}
\end{figure}

In a nutshell, the system is a ferromagnet in which the
interactions happen in $p$-plets, instead of more classical couples,
and the interacting agents live on a diluted random network, i.e.
the Erd\"{o}s-Renyi graph. In general, the graph can be specified by fixing the
 number of nodes $N$ and its ``connectivity'' $\alpha$, which represents the number of nearest neighbors per site.

As standard ferromagnets, the model is shown to exhibit two
phases, a paramagnetic one and a (replica symmetric) ferromagnetic
one, on the the other hand, as a difference with respect to the
standard ferromagnet, the phase transition does not display
criticality for $p>2$. The model is investigated by means of
cavity field technique and extensive numerical simulations.

We find an expression for the free energy as a function of $p$, of
the network connectivity  $\alpha$ and of the (inverse)
temperature $\beta$, showing that it is consistent with known
results. In particular, by properly tuning $p$ and $\alpha$ we
recover the Curie-Weiss model \cite{abarra}, the diluted Ising
model \cite{ABC} and the fully-connected $p$-spin model
\cite{abarra3}; moreover, regardless of the (finite) dilution, for
$p=2$ criticality is restored. Full agreement with Monte Carlo
simulations is obtained both on the absence of the critical
behavior and on the free energy structure.

The paper is organized as follows: In Section $2$ the model is
introduced and some of its properties worked out together with the
introduction of a proper statistical mechanics machinery, while in
section $3$ its equilibrium is solved via the smooth cavity field
technique. Section $4$ deals with the properties of the free
energy and its consistency with well known models, while in
section $5$ our numerical analysis is presented. Section $6$ is
left for a summary and outlook. Finally, Section  $7$, as an
Appendix, contains the detailed proofs of the theorems introduced.

\section{The diluted even-$p$-spin ferromagnet}
In this section we explore the properties of a diluted even-$p$-spin ferromagnet: we restrict ourselves only to even values of
$p$ for mathematical convenience as the investigation with the
cavities is much simpler. However, due to monotonicity of all the
observables in $p$, such restriction does not imply any loss of
generality, as confirmed also by numerical simulations performed
on both even and odd values of $p$.

Before proceeding, it is worth recalling some concepts concerning
the diluted random network where the magnetic system is set. Such
a network is an Erd\"{o}s-Renyi graph \cite{ER} defined as
follows: given a number $N$ of nodes, we introduce connections
between them in such a way that each pair of vertices $i,j$ has a
connecting link with independent probability equal to $\alpha/N$,
with $0 \leq \alpha \leq N$. As a result, the probability
distribution for the number of links per node (or coordination
number) is binomial with average $\alpha$. Hence, the parameter
$\alpha$ provides a measure of the ``degree of connectivity'' of
the graph itself: the smaller $\alpha$ the more diluted the
system; for $\alpha = 0$ and $\alpha =N$ the extreme cases of
fully disconnected and fully connected graphs, respectively, are
recovered. Notice that in the thermodynamic limit $N \rightarrow
\infty$ the binomial distribution converges towards the Poisson
distribution \cite{BCC}.

The ER graph can be algebraically described by the so-called adjacency
matrix $\mathbf{A}$ which is an $N \times N$ symmetric matrix
whose entry $A_{ij}$ is $1$ if $i \neq j$ and the two nodes are
connected together, otherwise it is zero.

We now associate to each node $i$ a binary variable $\sigma_{i}=\pm1, \ i \in
[1,N]$ and we introduce $p$ families
$\{i_{\nu}^1\},\{i_{\nu}^2\},...,\{i_{\nu}^p\}$ of i.i.d. random
variables uniformly distributed on the previous interval. Then, the
Hamiltonian is given by the following expression
\begin{equation}\label{ham}
H_{N}(\sigma,\gamma)=-\sum_{\nu=1}^{k_{\gamma N}}
\sigma_{i_{\nu}^1}\sigma_{i_{\nu}^2}...\sigma_{i_{\nu}^p},
\end{equation}
where $k$ represents the number of connected $p$-plets present in
the graph. Reflecting the underlying network, $k$ is a Poisson
distributed random variable with mean value $\gamma N$. The
relation among the coordination number $\alpha$ and $\gamma$ is
$\gamma \propto \alpha^{p-1}$: this will be easily understood a
few lines later by a normalization argument coupled with the high
connectivity limit of this mean field model.

The quenched expectation of the model is given by the composition
of the Poissonian average with the uniform one performed over the
families $\{i_{\nu}\}$
\begin{equation}
\textbf{E}[\cdot] = E_PE_i[\cdot] = \sum_{k=0}^{\infty}
\frac{e^{-\gamma N}(\gamma
N)^k}{k!N^p}\sum_{i_{\nu}^1....i_{\nu}^p}^{1,N}[\cdot],
\end{equation}
where the term $N^p \approx N!/(N-p)!$ accounts for the number of
possible ordered $p$-plets.

As they will be useful in our derivation, it is worth stressing
the following properties of the Poisson distribution: Let us
consider a function $g:\mathbb{N}\to\mathbb{R}$, and a Poisson
variable $k$ with mean $\gamma N$, whose expectation is denoted by
$\mathbb{E}$.

It is easy to verify that
\begin{eqnarray}\label{Pp1}
\mathbb{E}[k g(k)] &=& \gamma N \mathbb{E} [g(k-1)] \\ \label{Pp2}
\partial_{\gamma N}\mathbb{E}[g(k)] &=&
\mathbb{E}[g(k+1)-g(k)]\\ \label{Pp3}
\partial^2_{(\gamma N)^2}\mathbb{E}[g(k)] &=&
\mathbb{E}[g(k+2)-2g(k+1)+g(k)].
\end{eqnarray}

The Hamiltonian written as in eq.~(\ref{ham}), has the advantage that it
is the sum of (a random number of) i.i.d.\ terms. To see the
connection to a more familiar Hamiltonian written in terms of
adjacency matrix elements, we first notice that being $\alpha /N$ the probability that two nodes are connected, among the $N^p$ possible
$p$-plets, the number of connected $p$-plets is
Poisson-distributed with average $\alpha^{p-1} N + O(\sqrt{N})$ for
large $N$. We now define the adjacency tensor $A_{i_1,...,i_p} \equiv A_{i_1,i_2}A_{i_1,i_3}...A_{i_1,i_p}$ which equals $1$ whenever the $p$-plet  $i_1,...,i_p$ occurs to be connected;  $A_{i_1,...,i_p}$ is Poisson distributed and has mean $\gamma N/N^{p} \sim (\alpha / N)^{p-1}$.
Hence, we can write the following Hamiltonian which is thermodynamically equivalent
to $H_N(\sigma,\gamma)$ appearing in eq.~(\ref{ham}):
\begin{equation}\label{parag}
-H_N(\sigma;  \gamma) \sim - \hat{H}_N(\sigma;\textbf{A}) = \sum_{i_1,...,i_p}^N
A_{i_1,...,i_p}\sigma_{i_1}...\sigma_{i_p}.
\end{equation}

Then, it is enough to consider the streaming of  the following
interpolating free energy (whose structure proves the statement a
priori by its thermodynamic meaning), depending on the real
parameter $t\in[0,1]$
$$
\phi(t) = \frac{\mathbb{E}}{N}\ln\sum_{\sigma}e^{\beta
\large(\sum_{\nu=1}^k \sigma_{i^1_{\nu}}...\sigma_{i^p_{\nu}} +
\sum_{i_1,...,i_p}^N A_{i_1,...,i_p}
\sigma_{i_1}...\sigma_{i_p}\large)},
$$
where $k$ is a Poisson random variable with mean $\gamma N t$ and
$A_{i_1,...,i_p}$ are random Poisson  variables with mean
$(1-t)\gamma/N^{p-1}$. In this way the two separated models are
recovered in the two extremals of the interpolation (for $t=0,1$).
By computing the $t$-derivative, we get
\begin{eqnarray}
\frac{1}{\gamma}\frac{d \phi(t)}{dt} &=& \mathbb{E}\ln(1+
\Omega(\sigma_{i_0^1}...\sigma_{i_0^p})\tanh(\beta))
\\ \nonumber &-& \frac{1}{N^{p}}\sum_{i_1,...,i_p}^N\ln(1+\Omega(\sigma_{i_1}...\sigma_{i_p})\tanh(\beta))=0,
\end{eqnarray}
where the label $0$ in $i_0^k$ stands for a new spin, born in the
derivative, according to the Poisson property (\ref{Pp2}); as
the $i_0$'s are independent of the random site indices in the
$t$-dependent $\Omega$ measure, the equivalence is proved.

\bigskip

Following a statistical mechanics approach, we know that the
macroscopic behavior, versus the connectivity $\alpha$ and the
inverse temperature  $\beta=1/T$, is described by the following
free energy density (often called {\itshape quenched pressure})
\begin{eqnarray}
A(\alpha,\beta) &=& \lim_{N \to \infty} A_N(\alpha,\beta) \\
\nonumber &=& \lim_{N \to \infty}\frac1N\textbf{E}\ln
Z_N(\gamma,\beta),\
\end{eqnarray}
where
\begin{equation}
Z_N(\gamma,\beta) = \sum_{\{\sigma\}}e^{-\beta H_N(\sigma,\gamma)}
\end{equation}
is the partition function. Taken $g(\sigma)$ as a generic function, the Boltzmann state is therefore given
by
\begin{equation}
\omega(g(\sigma)) = \frac{1}{Z_{N}(\gamma,\beta)}
\sum_{\{\sigma_N\}} g(\sigma) e^{-\beta H_{N}(\sigma,\gamma)},
\end{equation}
with its replicated form
\begin{equation}
\Omega(g(\sigma)) = \prod_s \omega^{(s)}(g(\sigma^{(s)}))
\end{equation}
and the total average $\langle g(\sigma)\rangle$ is defined as
\begin{equation}
\langle g(\sigma)\rangle = \textbf{E}[\Omega(g(\sigma))].
\end{equation}
\medskip
\newline
Let us introduce further, as order parameters of the theory, the
multi-overlaps \begin{equation}
q_{1...n}=\frac1N\sum_{i=1}^{N}\sigma^{(1)}_{i}...\sigma^{(n)}_{i},
\end{equation} with a particular attention at the magnetization $m = q_1
=(1/N)\sum_{i=1}^{N}\sigma_{i}$ and to the two replica overlap
$q_{12} = (1/N)\sum_{i=1}^{N}\sigma_{i}^1\sigma_i^2$.

The normalization constant of the quenched pressure can be checked
by performing the expectation value of the cost function:
\begin{eqnarray}\nonumber
\textbf{E}[H] &=& -\gamma N m^p \\  \textbf{E}[H^2] -
\textbf{E}^2[H] &=& \gamma^2 N^2\Big[(q_{12}^p - m^p ) + O \left(
\frac{1}{N} \right) \Big],
\end{eqnarray}
by which it is easy to see that the model is well defined, in
particular it is linearly extensive in the volume $N$. Then, in
the high connectivity limit each agent interacts with all the
others ($\alpha \sim N$) and, in the thermodynamic limit, $\alpha
\to \infty$. Now, such a high-connectivity limit, i.e. a linear
divergence of $\alpha$, is properly recovered for any finite $p$,
$p<N$. In particular, if $p=2$ the amount of couples in the
summation scales as $N(N-1)/2$ and $\gamma = 2\alpha$; if $p=3$
the amount of triples scales as $N(N-1)(N-2)/3!$ and, with $\gamma
= 3!\alpha^2$.

Before starting our free energy analysis, we want to point out
also  the connection between this diluted version and the fully
connected counterpart.
Let us remember that the Hamiltonian of the fully connected  $p$-spin
model (FC) can be written as \cite{abarra3} \begin{equation} H^{FC}_{N}(\sigma)
= \frac{p!}{2 N^{p-1}}\sum_{1 \leq i_1 < ... < i_p \leq N}
\sigma_{i_1}\sigma_{i_2}...\sigma_{i_p}, \end{equation} and let us consider
the  trial function $\hat{A}(t)$ defined as follows \begin{equation} \hat{A}(t)
= \frac{1}{N}\mathbb{E}\ln \sum_{\sigma} \exp\Big[ \beta
\sum_{\nu}^{P_{\gamma N
t}}\sigma_{i_{\nu}^1}\sigma_{i_{\nu}^2}...\sigma_{i_{\nu}^p} +
(1-t)\frac{\beta' N}{2}m^p \Big], \end{equation} which interpolates between
the fully connected  $p$-spin model and the diluted one, such that
for $t=0$ only the fully connected survives, while the opposite
happens for $t=1$. Let us work out the derivative with respect to
$t$ to obtain \begin{eqnarray} \partial_t \hat{A}(t) &=&
(p-1)\alpha^{p-1}\ln\cosh(\beta) \\ \nonumber &-& (p-1)
\alpha^{p-1}\sum_{n}\frac{-1^n}{n}\theta^n \langle q_n^p\rangle -
\frac{\beta'}{2}\langle m^p \rangle, \end{eqnarray} by which we
see that the correct scaling, in order to recover the proper
infinite connectivity model, is obtained when $\alpha \to \infty$,
$\beta \to 0$ and $\beta' = 2(p-1) \alpha^{p-1}\tanh(\beta)$ is
held constant.

\begin{remark}
It is worth noting that for $p=2$ we recover the correct scaling
of the diluted Curie-Weiss model \cite{ABC}, furthermore the
dilute $p$-spin model reduces to the fully connected one, in the
infinite connectivity limit, uniformly in the size of the system.
\end{remark}

\section{The smooth cavity approach}

In this section we want to look for an iterative expression of the
free energy density by using a version of the cavity strategy
\cite{barra1,abarra} that we briefly recall: the idea behind the
cavity techniques \cite{guerra2,MPV}, which, for our purposes,
resembles the stochastic stability approach \cite{cg2,parisiSS},
is that information concerning the free energy density can be
extrapolated when looking at the incremental extensive free energy
given by the addition of a spin.
\newline
In diluted models, this additional spin changes also (infinitesimally
in the high $N$ limit) the connectivity and, in evaluating how the
free energy density varies conformingly with this, we are going to
prove that it can be written by a cavity function and such a
connectivity shift.
\newline
So the behavior of the system is encoded into these two parts. The
latter is simpler as it is made up only by stochastically stable
terms (a proper definition of these terms will follow in the current
section). Conversely, the former term needs to be expressed via these
terms and this must be achieved by iterative expansions.

At first we show how the free energy density can be decomposed via
these two parts (the cavity function and the connectivity shift).
Then, we analyze each term separately. We will see that they can
be expressed by the momenta of the magnetization and of the
multi-overlaps, weighted in a perturbed Boltzmann state, which
recovers the standard one in the thermodynamic limit.
\begin{theorem}\label{primolevi}
\textit{In the thermodynamic limit, the quenched pressure of the
even  $p$-spin diluted ferromagnetic model is given by the
following expression} \begin{equation} A(\alpha,\beta) = \ln2
-\frac{\alpha}{p-1}\frac{d}{d\alpha}A(\alpha,\beta) +
\Psi(\alpha,\beta,t=1), \end{equation}
\end{theorem}
where the cavity function $\Psi(t,\alpha,\beta)$ is introduced as
\begin{eqnarray}\label{cavity1}
&& \textbf{E}\Big[\ln\frac{\sum_{\{\sigma\}}
e^{\beta\sum_{\nu=1}^{k_{\tilde{\gamma}N}}
\sigma_{i_{\nu}^1}\sigma_{i_{\nu}^2}...\sigma_{i_{\nu}^p}}\;
e^{\beta\sum_{\nu=1}^{k_{2\tilde{\gamma}t}}
\sigma_{i_{\nu}^1}\sigma_{i_{\nu}^2}...\sigma_{i_{\nu}^{p-1}}}}
{\sum_{\{\sigma\}} e^{\beta\sum_{\nu=1}^{k_{\tilde{\gamma}N}}
\sigma_{i_{\nu}^1}\sigma_{i_{\nu}^2}...\sigma_{i_{\nu}^p}}}\Big]= \nonumber \\
&&  \textbf{E}\Big[\ln \frac{Z_{N,t}(\tilde{\gamma},\beta)}
{Z_{N}(\tilde{\gamma},\beta)}\Big] =
\Psi_N(\tilde{\gamma},\beta,t),
\end{eqnarray}
with \begin{equation}\label{cavity2} \Psi(\gamma,\beta,t) =
\lim_{N\rightarrow\infty}\Psi_N(\tilde{\gamma},\beta,t). \end{equation}
\medskip
\newline
For the sake of clearness and to avoid interrupting the paper with
long technical calculations, the proof of the theorem is reported
in the Appendix.
\newline
Thanks to the previous theorem, it is possible to figure out an
expression for the pressure by studying the properties of the
cavity function $\Psi(\alpha,\beta)$ and the connectivity shift
$\partial_{\alpha}A(\alpha,\beta)$.
\newline
Using the properties of the Poisson distribution (\ref{Pp1},
\ref{Pp2}), we can write
\begin{eqnarray}
\frac{d}{d\alpha}A(\alpha,\beta)
&=& \frac{(p-1)}{N}\alpha^{p-2}\frac{d}{d\gamma}\textbf{E}\Big[\ln Z_N(\gamma,\beta)\Big] = \nonumber \\
&=& (p-1)\alpha^{p-2}\textbf{E}\Big[\ln \sum_{\{\sigma\}}
e^{\beta\sum_{\nu=1}^{k+1}
\sigma_{i_{\nu}^1}...\sigma_{i_{\nu}^p}} -  \nonumber \\ \nonumber
&-& \ln \sum_{\{\sigma\}} e^{\beta\sum_{\nu=1}^{k}
\sigma_{i_{\nu}^1}...\sigma_{i_{\nu}^p}}\Big]. \nonumber
\end{eqnarray}
Now considering the relation (and definition)
\begin{eqnarray}
e^{\beta\sigma_{i_0^1}...\sigma_{i_0^p}} &=&
\cosh\beta + \sigma_{i_0^1}...\sigma_{i_0^p}\sinh\beta, \\
\theta &=& \tanh\beta,
\end{eqnarray}
we can write
\begin{eqnarray}
&& \frac{d}{d\alpha}A(\alpha,\beta) = \\ \nonumber &&
(p-1)\alpha^{p-2}\Big[\ln\cosh\beta + \textbf{E}[\ln(1 +
\omega(\sigma_{i_{\nu}^1}...\sigma_{i_{\nu}^p})\theta)]\Big].
\end{eqnarray}
At the end, expanding the logarithm, we obtain
\begin{eqnarray}
\frac{d}{d\alpha}A(\alpha,\beta) &=&
(p-1)\alpha^{p-2}\ln\cosh\beta - \\ \nonumber
 &-& (p-1)\alpha^{p-2}\sum_{n=1}^{\infty}\frac{(-1)^n}{n}\theta^n
\langle q_{1,...,n}^p \rangle.
\end{eqnarray}
\medskip
\newline
With the same procedure it is possible to show that
\begin{eqnarray}\label{Psit}
\frac{d}{dt}\Psi(\tilde{\alpha},\beta,t) &=&
2\tilde{\alpha}^{p-1}\ln\cosh\beta - \\ \nonumber &-&
2\tilde{\alpha}^{p-1}\sum_{n=1}^{\infty}\frac{(-1)^n}{n}\theta^n
\langle q_{1,...,n}^{p-1} \rangle_{\tilde{\alpha},t},
\end{eqnarray}
where
$$
\tilde{\alpha} = \alpha \Big[\frac{N}{N+1}\Big]^{\frac{1}{p-1}}
\qquad.
$$
Now, by eq.~(\ref{Psit}), we see that
even the cavity function, once the r.h.s. of eq.(\ref{Psit}) is
integrated back against $t$, can be expressed via all the order
parameters of the model:
\begin{equation}\nonumber
\Psi(\tilde{\alpha},\beta,t)=2\tilde{\alpha}^{p-1}\Big(\ln\cosh(\beta)-
\sum_{n=1}^{\infty}\frac{(-\theta)^n}{n}\int_0^t \langle
q_{1,...,n}^{p-1}\rangle_{\tilde{\alpha},t} \Big).
\end{equation}

So, as expected, we can understand the properties of the free
energy by analyzing the properties of the order parameters:
magnetization and overlaps, weighted in their extended Boltzmann
state $\tilde{\omega}_t$.
\newline
Further, as we expect that the order parameters is able to describe
thermodynamics even in the true Boltzmann states $\omega, \Omega$
\cite{landau}, accordingly to the following definitions, we are
going to show that {\em filled} order parameters (the ones
involving even numbers of replicas) are stochastically stable or,
in other words, are independent of the $t$-perturbation in the
thermodynamic limit, while the others, not filled, become filled,
again in this limit (such that even for them $\omega_t\to \omega$
in the high $N$ limit and thermodynamics is recovered). The whole
is explained in the following definitions and theorems of this
section.

\begin{definition}
We define the t-dependent Boltzmann state $\tilde{\omega}_t$ as
\begin{eqnarray}\label{dente}
&& \tilde{\omega}_t(g(\sigma)) = \\ \nonumber &&
\frac{1}{Z_{N,t}(\gamma,\beta)} \sum_{\{\sigma\}}g(\sigma)
e^{\beta\sum_{\nu=1}^{k_{\tilde{\gamma}N}}
\sigma_{i_{\nu}^1}...\sigma_{i_{\nu}^p} +
\beta\sum_{\nu=1}^{k_{2\tilde{\gamma}t}}
\sigma_{i_{\nu}^1}...\sigma_{i_{\nu}^{p-1}}}\label{dante},
\end{eqnarray}
where $Z_{N,t}(\gamma,\beta)$ extends the classical partition
function in the same spirit of the numerator of eq.(\ref{dente})
itself, and $\tilde{\gamma} = \gamma(1+N^{-1})$.
\end{definition}
We see that the original Boltzmann state of a $N$-spin system is
recovered as $t$ approached $0$, while, in the limit $t \to 1$ and
gauging the spins, it is possible to build a Boltzmann state of a
$N+1$ spins, with a little shift both in $\alpha,\beta$, which
vanishes in the $N \to \infty$ limit.

Now, coherently with the implication of thermodynamic limit (by which
$A_{N+1}(\alpha,\beta)-A_N(\alpha,\beta) =0$ for $N \to \infty$),
we are going to define the {\em filled} overlap monomials and show
their independence (stochastic stability) with respect to the
perturbation encoded by the interpolating parameter $t$.
These parameters are already `good" order parameters describing
the theory, while the others (the not-filled ones) must be
expressed via the formers, and this will be achieved by expanding
them with a suitably introduced streaming equation.

\begin{definition} We can split the class of monomials of the
order parameters in two families:
\begin{itemize}

\item We define {\itshape filled} or equivalently {\itshape stochastically
stable} those overlap monomials with all the replicas appearing an
even number of times (i.e. $q_{12}^2$,\ $m^2$,\
$q_{12}q_{34}q_{1234}$).

\item We define {\itshape non-filled}
those overlap monomials with at least one replica appearing an odd
number of times (i.e. $q_{12}$,\ $m$,\ $q_{12}q_{34}$).
\end{itemize}
\end{definition}
We are going to show three theorems that will play a guiding role
for our iteration: as this approach has been deeply developed in
similar contexts (as fully connected Ising and  $p$-spin models
\cite{abarra,abarra3}, fully connected spin glasses \cite{barra1}
or diluted ferromagnetic models \cite{ABC,BCC}, which are the
``boundaries'' of the model of this paper) we will not show all
the details of the proof, but we sketch them in the appendix as
they are really intuitive. The interested reader will found a
clear derivation in the appendix and can deepen this point by
looking at the original works.
\begin{theorem}\label{ciccia}
In the thermodynamic limit and setting $t=1$ we have \begin{equation}
\tilde{\omega}_{N,t}(\sigma_{i_1}\sigma_{i_2}...\sigma_{i_n}) =
\tilde{\omega}_{N+1}(\sigma_{i_1}\sigma_{i_2}...\sigma_{i_n}\sigma_{N+1}^n).
\end{equation}
\end{theorem}
\begin{theorem}\label{saturabili}
Let $Q_{ab}$ be a not-filled monomial of the overlaps (this means
that $q_{ab}Q_{ab}$ is filled). We have \begin{equation}
\lim_{N\rightarrow\infty}\lim_{t\rightarrow1} \langle Q_{ab}
\rangle_t = \langle q_{ab}Q_{ab} \rangle, \end{equation} (examples:
\newline for $N \rightarrow \infty$ we get $\langle m_1 \rangle_t
\rightarrow \langle m_1^2 \rangle,\quad \langle q_{12} \rangle_t
\rightarrow \langle q_{12}^2 \rangle$).
\end{theorem}
\begin{theorem}\label{saturi}
In the $N\rightarrow\infty$ limit, the averages
$\langle\cdot\rangle$ of the filled polynomials are t-independent
in $\beta$ average.
\end{theorem}
\medskip

\section{Properties of the free energy}
In this section we are going to address various points: at first
we work out the constraints that the model must fulfil, which are
in agreement both with a self-averaging behavior of the
magnetization and with the replica-symmetric behavior of the
multi-overlaps \cite{gulielmo}; then we write an iterative
expression for the free energy density and its links with known
models as diluted ferromagnets ($p\to2$ limit) and fully connected
$p$-spin models ($\alpha \to \infty$ limit).

With the following definition \begin{eqnarray} \tilde{\beta} &=&
2(p-1)\tilde{\alpha}^{p-1}\theta \\ \nonumber &=&
2(p-1)\alpha^{p-1}\frac{N}{N+1}\theta \quad
\stackrel{N\rightarrow\infty}{\longrightarrow}
2(p-1)\alpha^{p-1}\theta = \beta', \end{eqnarray} we show (and
prove in the Appendix) the streaming of replica functions, by
which not filled multi-overlaps can be expressed via filled ones.
\begin{proposition}\label{stream}
Let $F_s$ be a function of s replicas. Then the following
streaming equation holds
\begin{eqnarray}
\frac{\partial\langle F_s \rangle_{t,\tilde{\alpha}}}{\partial t}
&=& \tilde{\beta} \Big[\sum_{a=1}^s\langle F_s
m_a^{p-1}\rangle_{t,\tilde{\alpha}} - s \langle F_s
m_{s+1}^{p-1}\rangle_{t,\tilde{\alpha}}\Big] \quad
\\ \nonumber
&+& \tilde{\beta}\theta \Big[ \sum_{a<b}^{1,s}\langle F_s
q_{a,b}^{p-1} \rangle_{t,\tilde{\alpha}} - s\sum_{a=1}^s\langle
F_s q_{a,s+1}^{p-1}\rangle_{t,\tilde{\alpha}} \\ \nonumber &+&
\frac{s(s+1)}{2!}\langle F_s
q_{s+1,s+2}^{p-1}\rangle_{t,\tilde{\alpha}}\Big] + O(\theta^2).
\end{eqnarray}
\end{proposition}
\begin{remark}
We stress that, at the first two level of approximation presented
here, the streaming has the structure of a $\theta$-weighted
linear sum of the Curie-Weiss streaming ($\theta^0$ term)
\cite{abarra} and the Sherrington-Kirkpatrick streaming
($\theta^1$ term) \cite{barra1}, providing mathematical structures of disordered systems with a certain degree of
independence with respect to the kind of quenched noise (frustration or
dilution).
\end{remark}

\bigskip

It is now immediate to obtain the linear order parameter
constraints (often known as Aizenman-Contucci polynomials
\cite{ac,abarra,BCC}) of the theory: in fact, the generator of
such a constraint is the streaming equation when applied on each
filled overlap monomial (or equivalently it is possible to apply
the streaming on a not-filled one and then gauge the obtained
expression; for the sake of clearness both the methods will be
exploited, the former for $q_2$ and the latter for $m$).

As examples, dealing with the terms $m^{p-1}$ and $q_{2}^{p-1}$,
it is straightforward to check that \begin{eqnarray}\nonumber 0
&=& \lim_{N \to \infty} \frac{\partial\langle m_N^{p-1}
\rangle_{t,\tilde{\alpha}}}{\partial t} = \tilde{\beta} \Big(
\langle m_1^{2(p-1)} \rangle - \langle m_1^{p-1} \rangle^2 \Big)\\
\nonumber &+& \tilde{\beta}\theta \Big( \langle m_1^{p-1}q_2^{p-1}
\rangle - \langle m_1^{p-1}\rangle \langle  q_2^{p-1} \rangle
\Big) + O(\theta^3), \nonumber \end{eqnarray} then, by gauging the
above expression, in the thermodynamic limit, (as $\lim_{N \to
\infty}\langle m_N^{p-1}\rangle_t \to \langle m^p \rangle$), we
get
$$
\Big( (\langle m_1^{2p} \rangle - \langle m_1^p \rangle^2) +
\theta ( \langle q_{2}^{2p} \rangle - \langle q_{2}^p \rangle^2
)\Big)=0, \ \ \forall \theta \in \mathbb{\mathcal{R}}^+.
$$
The fact that the previous expression holds for every $\theta$
suggests self-averaging for the energy (by which all the linear
constraints can be derived \cite{BCC}) due to the first term, as
well as replica symmetric behavior of the two replica overlap due
to the last one. Analogously, the contribution of the $\langle q_{2}^2\rangle$
generator is  \begin{eqnarray}  \nonumber 0 &=& \Big( (\langle
q_{12}^{p-1} m_1^{p-1} \rangle + \langle q_{12}^{p-1} m_2^{p-1}
\rangle - 2 \langle q_{12}^{p-1} m_3^{p-1} \rangle) +  \\
\nonumber &+& \theta ( \langle q_{12}^{p-1} q_{12}^{p-1} \rangle -
4 \langle q_{12}^{p-1} q_{23}^{p-1} \rangle + 3 \langle
q_{12}^{p-1} q_{34}^{p-1} \rangle) \Big), \end{eqnarray} which
shows replica symmetric behavior of the magnetization by the first
term and the classical Aizenman-Contucci relations \cite{ac,BCC}
by the latter.

Furthermore, turning now our attention to the free energy, it is
easy to see that the streaming equation  allows to generate all
the desired overlap functions coupled to every well behaved $F_s$.
In this way, if $F_s$ is a not filled overlap, we can always
expand recursively it into a filled one, the only price to pay
given by the $\theta$ order that has to be reached or, which is
equivalent, the number of derivatives that have to be performed.

Let us now remember the t-derivative of the cavity function
(\ref{Psit}), showing explicitly the first two terms of its
expansion
\begin{eqnarray}\label{Psit2}
\frac{d}{dt}\Psi(\tilde{\alpha},\beta,t) &=&
2\tilde{\alpha}^{p-1}\ln\cosh\beta + \tilde{\beta}\langle
m_1^{p-1} \rangle_{\tilde{\alpha},t} - \\ \nonumber  &-&
\frac{\tilde{\beta}}{2}\theta\langle q_{12}^{p-1}
\rangle_{\tilde{\alpha},t} -
2\tilde{\beta}^{p-1}\sum_{n=3}^{\infty}\frac{-1^n\theta^n}{n}
\langle q_{1,...,n}^{p-1} \rangle_{\tilde{\alpha},t}.
\end{eqnarray}

As derivative of fillable terms involve filled ones, we can arrive
to an analytical form of $\Psi(\alpha,\beta)$ if we calculate it
as the $t$-integral of its $t$-derivative, together with the
obvious relation $\Psi(t=0) = 0$. Hence, if we apply the streaming
equation machinery to the overlaps constituting equation
(\ref{Psit2}), we are able to fill them and to remove their
$t$-dependence in the thermodynamic limit. In this way we are
allowed to bring them out from the final $t$-integral.
\newline
In fact, without gauging (so, not only in the ergodic regime,
where symmetries are preserved), we can expand the streaming of
$\langle m^{p-1}\rangle_t$:
\newline
\begin{eqnarray}
\frac{d\langle m_1^{p-1} \rangle_t}{dt} &=&
\tilde{\beta}\Big[\langle m_1^{2(p-1)} \rangle - \langle m_1^{p-1}m_2^{p-1} \rangle_t\Big] + \nonumber \\
&-&  \tilde{\beta}\theta\Big[\langle m_1^{p-1}q_{12}^{p-1}
\rangle_t - \langle m_1^{p-1}q_{23}^{p-1} \rangle_t\Big] +
O(\theta^2). \nonumber
\end{eqnarray}
We can note the presence of the filled monomial $\langle
m_1^{2(p-1)} \rangle$, whose $t$-dependence has been omitted
explicitly to underline its stochastic stability, while the overlaps
$\langle m_1^{p-1}m_2^{p-1} \rangle_t$ and $\langle
m_1^{p-1}q_{12}^{p-1} \rangle_t$ can be saturated in two steps of
streaming. This will be sufficient, wishing to have a fourth order
expansion for the cavity function.
\newline
We now derive these two functions and apply the same scheme to all
the overlaps that appear and that have to be necessary filled in
order to obtain the desired result.
\begin{eqnarray}
&& \frac{d\langle m_1^{p-1}m_2^{p-1} \rangle_t}{dt} = \nonumber \\
&& 2\tilde{\beta}\Big[
\langle m_1^{2(p-1)}m_2^{p-1} \rangle_t - \langle m_1^{p-1}m_2^{p-1}m_3^{p-1} \rangle_t\Big] + \nonumber \\
&& \theta \tilde{\beta} \Big[\langle
m_1^{p-1}m_2^{p-1}q_{12}^{p-1} \rangle -
4\langle m_1^{p-1}m_2^{p-1}q_{13}^{p-1} \rangle_t +  \nonumber \\
&&  3\langle m_1^{p-1}m_2^{p-1}q_{34}^{p-1} \rangle_t\Big],
\label{m1m2}
\end{eqnarray}
\begin{eqnarray}
&& \frac{d\langle m_1^{2(p-1)}m_2^{p-1} \rangle_t}{dt} = \\
\nonumber && 2\tilde{\beta}\Big[ \langle m_1^{2(p-1)}m_2^{2(p-1)}
\rangle_t\Big] +\tilde{\beta}\Big[\mbox{unfilled terms}\Big] +
O(\theta^2).
\end{eqnarray}
Integrating back in $t$ and neglecting higher order terms we have
\begin{equation}
\langle m_1^{2(p-1)}m_2^{p-1} \rangle_t =
\tilde{\beta}\Big[\langle m_1^{2(p-1)}m_2^{2(p-1)} \rangle\Big]t,
\end{equation}
and we can write
\begin{eqnarray}
&& \langle m_1^{p-1}m_2^{p-1} \rangle_t = \\ \nonumber &&
\tilde{\beta}\theta\langle m_1^{p-1}m_2^{p-1}q_{12}^{p-1} \rangle
t + \tilde{\beta}^2\langle m_1^{2(p-1)}m_2^{2(p-1)} \rangle t^2.
\end{eqnarray}
Let us take a look now at the other overlap $\langle
m_1^{p-1}q_{12}^{p-1} \rangle_t$:
\begin{eqnarray}\label{m1m2}
\frac{d\langle m_1^{p-1}q_{12}^{p-1} \rangle_t}{dt} &=&
\tilde{\beta}\Big[
\langle m_1^{2(p-1)}q_{12}^{p-1} \rangle_t - \langle m_1^{p-1}m_2^{p-1}q_{12}^{p-1} \rangle_t  \nonumber \\
&-& 2\langle m_1^{p-1}m_2^{p-1}m_3^{p-1}q_{12}^{p-1}
\rangle_t\Big] + O(\theta^2),
\end{eqnarray}
that gives
\begin{equation}
\langle m_1^{p-1}q_{12}^{p-1} \rangle_t = \tilde{\beta}\langle
m_1^{p-1}m_2^{p-1}q_{12}^{p-1} \rangle t + O(\theta^2).
\end{equation}

At this point we can write for $\langle
m_1^{p-1}\rangle_{t,\tilde{\alpha}}$ (and consequently for
$\langle q_{12}^{p-1}\rangle_{t,\tilde{\alpha}}$)
\begin{eqnarray}\nonumber
\langle m_1^{p-1}\rangle_{t,\tilde{\alpha}} &=&
\tilde{\beta}\langle m_1^{2(p-1)}\rangle t -
\frac{\tilde{\beta}^3}{3}\langle m_1^{2(p-1)}m_2^{2(p-1)}\rangle
t^3 \\ \nonumber &-&
\tilde{\beta}^2\theta\langle m_1^{p-1}m_2^{p-1}q_{12}^{p-1}\rangle t^2 + O(\theta^3), \nonumber \\
\langle q_{12}^{p-1}\rangle_{t,\tilde{\alpha}} &=&
\tilde{\beta}\theta\langle q_{12}^{2(p-1)}\rangle t +
\tilde{\beta}^2\langle m_1^{p-1}m_2^{p-1}q_{12}^{p-1}\rangle t^2 +
O(\theta^3).\nonumber
\end{eqnarray}
With these relations, eq. (\ref{Psit2}) becomes
\begin{eqnarray}\nonumber
&& \frac{d}{dt}\Psi_N(\alpha,\beta,t) = 2\alpha^{p-1}\ln\cosh\beta
+ \tilde{\beta}^2\langle m_1^{2(p-1)}\rangle t  \\ \nonumber && \
\ \ \  -\frac{\tilde{\beta}^2\theta^2}{2}\langle
q_{12}^{2(p-1)}\rangle t - \frac{3\tilde{\beta}^3\theta}{2}
\langle m_1^{p-1}m_2^{p-1}q_{12}^{p-1}\rangle t^2 \\ \nonumber &&
\ \ \ \  - \frac{\tilde{\beta}^4}{3}\langle
m_1^{2(p-1)}m_2^{2(p-1)}\rangle t^3 + O(\theta^5),\nonumber
\end{eqnarray}
which  ultimately allows us to write an iterated expressions for
$\Psi$ evaluated at $t=1$
\begin{eqnarray}
&& \Psi_N(\alpha,\beta,1) = \\ \nonumber &&
2\alpha^{p-1}\ln\cosh\beta + \frac{\tilde{\beta}^2}{2}\langle
m_1^{2(p-1)}\rangle  -
\frac{\tilde{\beta}^2\theta^2}{4}\langle q_{12}^{2(p-1)}\rangle - \\
&& \frac{\tilde{\beta}^3\theta}{2} \langle
m_1^{p-1}m_2^{p-1}q_{12}^{p-1}\rangle -
\frac{\tilde{\beta}^4}{12}\langle m_1^{2(p-1)}m_2^{2(p-1)}\rangle
t^3 + O(\theta^5). \nonumber
\end{eqnarray}
Overall the result we were looking for, namely a Landau-like polynomial
form for the free energy, reads off as
\begin{eqnarray}\label{maino}
A(\alpha,\beta) &=& \ln2 \:+\: \alpha^{p-1}\ln\cosh\beta +\\
&+& \frac{\beta'}{2}\Big(\beta'\langle m^{2(p-1)}\rangle - \langle
m^{p}\rangle\Big) + \nonumber \\ \nonumber &+&
\frac{\beta'\theta}{4}\Big(\beta'\theta\langle
q_{12}^{2(p-1)}\rangle - \langle q_{12}^{p}\rangle\Big) +
O(\theta^5). \nonumber
\end{eqnarray}

Now, several conclusions can be addressed from the expression
(\ref{maino}):
\newline
In fact, as we are going to see immediately through remarks, this
formula can bridge free-energies of quite different models
(diluted versus non-diluted, critical versus uncritical) and acts
as a general free energy expression close to the phase transition.
\begin{remark}
At first let us note that, by constraining the interaction to be
pairwise, critical behavior should arise \cite{landau}.
Coherently, we see that for $p=2$ we can write the free energy
expansion as
$$A(\alpha,\beta)_{p=2} = \ln 2 + \alpha
\ln\cosh(\beta) - \frac{\beta'}{2}(1-\beta')\langle m^2 \rangle -
\frac{\beta' \theta}{4} \langle q_2^2 \rangle,$$ which coincides
with the one of the diluted ferromagnet \cite{ABC} and displays
criticality at $2\alpha\theta=1$, where the coefficient of the
second order term vanishes, in agreement with previous results
\cite{ABC} and Landau theory \cite{landau}.
\end{remark}
\begin{remark}
The free energy density of the fully connected $p$-spin model
is \cite{abarra3} $A(\beta')= \ln 2 + \ln\cosh(\beta m^{p-1}) -
(\beta/2)m^p$, which coincides with the expansion (\ref{maino}) in
the limit of $\alpha \to \infty$ and $\beta \to 0$ with $\beta' =
2(p-1) \alpha^{p-1}\theta$.
\end{remark}
\begin{remark}
It is worth noting that the connectivity no longer plays a linear
role in contributing to the free energy density, as it does happen
for the diluted two body models \cite{ABC,gt2}. This is
interesting in applications to economic networks, where, for high
values of coordination number it may be interesting to develop
strategies with more than one coupling \cite{cicciona}.
\end{remark}

\section{Numerics}\label{sec:numerics}
We now analyze the system described in the previous section, from
the numerical point of view by performing extensive Monte Carlo
simulations. Within this approach it is more convenient to use the
second Hamiltonian introduced (see eq.(\ref{parag})):
\begin{equation}\label{my_hamiltonian}
\hat{H}_N(\sigma, \mathbf{A})=-\sum_{i_i}^{N}
\sigma_{i_1}\sum_{i_2<i_3<...<i_p=1}^{N} A_{i_1,...,i_p}
\sigma_{i_{2}}\sigma_{i_{3}}...\sigma_{i_{p}}.
\end{equation}
The product between the elements of the adjacency tensor ensures
that the $p-1$ spins considered in the second sum are joined by a
link with $i_1$.
\newline
The evolution of the magnetic system is realized by means of a
single spin-flip dynamics based on the Metropolis algorithm
\cite{barkema}. At each time step a spin is randomly extracted and
updated whenever its coordination number is larger than $p-1$. For
$\alpha$ large enough (at least above the percolation threshold,
as  obviously holds for the results found previously) and $p=3,4$
this condition is generally verified. The updating procedure for a
spin $\sigma_i$ works as follows: Firstly we calculate the energy
variation $\Delta E_i$ due to a possible spin flip, which for
$p=3$ and $p=4$ reads respectively
\begin{eqnarray}
\Delta E_i &=& 2 \sigma_i \sum_{j<k=1}^{N} A_{i,j}A_{i,k}
\sigma_{j}\sigma_{k},
\\
\Delta E_i &=& 2 \sigma_i \sum_{j<k<w=1}^{N} A_{i,j}A_{i,k}A_{i,w}
\sigma_{j}\sigma_{k}\sigma_{w}.
\end{eqnarray}
Now, if $\Delta E_i <0$, the spin-flip $\sigma_i \rightarrow -
\sigma_i$ is realized with probability $1$, otherwise it is
realized with probability $e^{-\beta \Delta E}$.

The cases $p=3, 4$ were studied in details, while  for $p=2$ we refer to
\cite{ABC}. Our investigations are aimed to evidence the existence of a phase transition and its nature and also to highlight a proper scaling for the temperature as the
parameter $\alpha$ is tuned.

As for the first point, we measured the so-called Binder cumulants
defined as follows:
\begin{equation}
G_N(T(\alpha)) \equiv 1 - \frac{ \langle m^4 \rangle_N }{3  \langle m^2 \rangle_N ^2},
\end{equation}
where $\langle \cdot \rangle_N$ indicates the statistical average obtained for a system of size $N$ and $T=\beta^{-1}$ \cite{binder}. The study of Binder
cumulants is particularly useful to locate and catalogue the
phase transition. In fact, in the case of continuous phase
transitions, $G_N(T)$ takes a universal positive value at the
critical point $T_c$, namely all the curves obtained for different
system sizes $N$ cross each other. On the other hand, for a
first-order transition $G_N(T)$ exhibits a minimum at $T_{min}$,
whose magnitude diverges as $N$. Moreover, a crossing point at
$T_{cross}$ can be as well detected when curves pertaining to
different sizes $N$ are considered \cite{vollmayr}. Now, $T_{min}$
and $T_{cross}$ scale as $T_{min}-T_c \propto N^{-1}$ and
$T_{cross}-T_c \propto N^{-2}$, respectively.

\begin{figure}[tb]
\begin{center}
\includegraphics[height=60mm]{Binder_doppio.pdf}
\caption{\label{fig:Binder} Binder cumulants $G_L(T)$ for systems with $p=3$ and different size
$N$, as shown in the legend, and connectivity
$\alpha=50$ (left panel) and $alpha=80$ (right panel).}
\end{center}
\end{figure}

In Fig.~\ref{fig:Binder} we show data for $G_N(T)$ obtained for
systems of different sizes ($N=400$, $N=500$, and $N=800$) but
equal connectivity ($\alpha=50$ and $\alpha=80$, respectively) as
a function of the temperature $T$. The existence of a minimum is
clear and it occurs at $T \approx 625$ and $T \approx 1600$.
Similar results are found also for $p=4$ and they all highlight
the existence of a first-order phase transition (hence lack of
criticality) at a temperature which depends on the connectivity
$\alpha$.

In order to deepen the role of connectivity in the evolution of
the system we measure the macroscopic observable $\langle m \rangle$ and its (normalized) fluctuations $\langle m^2 \rangle -
\langle m \rangle^2$, studying their dependence on $T$ and on $\alpha$. Data for different choices of
size and dilution are shown in Figure \ref{fig:p3} for $p=3$ and in Figure \ref{fig:p4} for $p=4$.

The profile of the magnetization, with an abrupt jump, and the
correspondent peak found for its fluctuations confirm the
existence of a first order phase transition at a well defined
temperature $T_c$ whose value depends on the dilution $\alpha$.
More precisely, by properly normalizing the temperature in
agreement with analytical results, namely $\tilde\beta \equiv
\beta \; \alpha^{p-1}$ we found a very good collapse of all
the curves considered. Hence, we have agreement among analytic and
numerics concerning the scaling of the temperature as
$\alpha^{p-1}$. Moreover our data provide a very clear hint
suggesting that the critical temperature can be written as $T_c =
f(p) \alpha^{p-1}$, where $f(p)$ is a monotonic decreasing function of $p$.

\begin{figure}[tb]\begin{center}
\includegraphics[height=60mm]{mag_Susc.pdf}
\caption{\label{fig:p3}Magnetization (main figure) and its
normalized fluctuations (inset) for $3$-spin systems of different sizes and
different dilution as a function of $\beta \; \alpha^{p-1}$. The
collapse of all the curves provides a strong evidence for the
scaling of the temperature.} \end{center}
\end{figure}

\begin{figure}[tb]\begin{center}
\includegraphics[height=60mm]{p4.pdf}
\caption{\label{fig:p4}Magnetization for $4$-spin systems of different sizes and
different dilution as a function of $\beta \; \alpha^{p-1}$. The
collapse of all the curves provides a strong evidence for the
scaling of the temperature.} \end{center}
\end{figure}

\section{Conclusions}
In this paper we performed an analysis of the ferromagnetic
diluted $p$-spin model via cavity field technique and numerical
simulations. Several questions have been addressed, including an expression
for the free energy, the self-averaging families for the order
parameters and a study of the phase transition among a
paramagnetic and a ferromagnetic regime. Despite a
rigorous picture for the lacking of replica symmetry breaking in
diluted ferromagnet is still unavailable, we supported strong
evidence toward a full replica symmetric behavior in the whole
phase diagram. In particular, we showed the vanishing of criticality for $p>2$ and we found a proper scaling for the transition temperature as a function of the system dilution, namely $T_c \sim  \alpha^{p-1}$.

Further development should be two-fold: from one side the same
analysis is still to be performed on the X-OR-SAT model which
constitutes another element making up the class of models based on
binary agents with mean filed interaction. On the other side, the
whole mathematical architecture still suffers a not exhaustive
development; in fact the difference among even and odd $p$ model,
at least for large $p$, is thermodynamically almost irrelevant,
while the lacking of the gauge symmetry in the latter rules out
the method at this stage. Moreover, it is highlighted the need to
develop a Hamilton-Jacobi technique \cite{BG} in order to handle
this kind of problem to avoid the iteration procedure implied by
the cavity method.

\section{Appendix: Analytical proofs}

In this section the proofs of al the Theorems and the Proposition
$1$ are reported.

\bigskip

\textbf{Proof of Theorem \ref{primolevi}}
\newline
Bridging a system made of by $N+1$ spins with one made of by $N$
spins implies the definition of rescaled $\gamma, \alpha$
parameters, accordingly to \cite{ABC}\cite{BCC}
\begin{eqnarray}
\tilde{\gamma} &=& \gamma \frac{N}{N+1} \qquad
\stackrel{N\rightarrow\infty}{\longrightarrow} \quad \gamma \\
\tilde{\alpha} &=& \alpha \Big[\frac{N}{N+1}\Big]^{\frac{1}{p-1}}
\qquad \stackrel{N\rightarrow\infty}{\longrightarrow} \quad
\alpha.
\end{eqnarray}
We have, in distribution, the Hamiltonian of a system made of
$N+1$ particles writable as
\begin{eqnarray}\label{Hsplit}
H_{N+1}(\sigma,\gamma) &=& -\sum_{\nu=1}^{k_{\gamma (N+1)}}
\sigma_{i_{\nu}^1}\sigma_{i_{\nu}^2}...\sigma_{i_{\nu}^p}\,\sim\,
\\ \nonumber &-& \sum_{\nu=1}^{k_{\tilde{\gamma} N}}
\sigma_{i_{\nu}^1}\sigma_{i_{\nu}^2}...\sigma_{i_{\nu}^p} -
\sum_{\nu=1}^{k_{2\tilde{\gamma}}}
\sigma_{i_{\nu}^1}\sigma_{i_{\nu}^2}...\sigma_{i_{\nu}^{p-1}}\sigma_{N+1},
\end{eqnarray}
that we may rewrite as \begin{equation} H_{N+1}(\sigma,\gamma) =
H_N(\sigma,\tilde{\gamma}) + \hat{H}_N(\sigma,2\tilde{\gamma}).
\end{equation}
\medskip
\newline
Following the above decomposition, let us consider the partition
function of the same $N+1$ spin model and let us introduce the
gauge transformation $\sigma_{i} \rightarrow \sigma_i\sigma_{N+1}$
which is a symmetry of the Hamiltonian known as
\textit{spin-flip}.
\begin{eqnarray}
Z_{N+1}(\gamma,\beta) &\sim& \sum_{\{\sigma_{N+1}\}} e^{-\beta
H_{N}(\sigma,\tilde{\gamma}) -
\beta\hat{H}_{N}(\sigma,\tilde{\gamma})\sigma_{N+1}} = \label{T11}
\\ \nonumber &=& \sum_{\{\sigma_{N+1}\}}
e^{\beta H_{N}(\sigma,\tilde{\gamma}) +
\beta\sum_{\nu=1}^{k_{2\tilde{\gamma }}}
\sigma_{i_{\nu}^1}...\sigma_{i_{\nu}^{p-1}}\sigma_{N+1}} =
\label{T12} \\ \nonumber &=& 2 \sum_{\{\sigma_{N}\}}
e^{\beta\sum_{\nu=1}^{k_{\tilde{\gamma}N}}
\sigma_{i_{\nu}^1}...\sigma_{i_{\nu}^p} +
\beta\sum_{\nu=1}^{k_{2\tilde{\gamma }}} \sigma_{i_{\nu}^1}...\sigma_{i_{\nu}^{p-1}}} = \label{T13} \\
&=& 2 Z_N(\tilde{\gamma},\beta)\tilde{\omega}(e^{-\beta
\hat{H}_{N}}),\nonumber
\end{eqnarray}
where  the new Boltzmann state $\tilde{\omega}$, and its
replicated $\tilde{\Omega}$, are introduced as \begin{eqnarray}
\tilde{\omega}(g(\sigma)) &=&
\frac{\sum_{\{\sigma_{N}\}}g(\sigma)e^{-\beta
H_N(\tilde{\gamma},\sigma)}} {\sum_{\{\sigma_{N}\}}e^{-\beta
H_N(\tilde{\gamma},\sigma)}}, \\  \qquad \tilde{\Omega}(g(\sigma))
&=& \prod_i\tilde{\omega}^{(i)}(g(\sigma^{(i)})). \end{eqnarray}
\newline
To continue the proof we now take the logarithm of both sides of
the last expression in eq. (\ref{T11}), apply the expectation
$\textbf{E}$ and subtract the quantity $\textbf{E}[\ln
Z_{N+1}(\tilde{\gamma},\beta)]$. We obtain
\begin{eqnarray}\nonumber
&& \textbf{E}[\ln Z_{N+1}(\gamma,\beta)] - \textbf{E}[\ln
Z_{N+1}(\tilde{\gamma},\beta)]=  \\  && \ln2 - \textbf{E}[\ln
\frac{Z_{N+1}(\gamma,\beta)}{Z_N(\tilde{\gamma},\beta)}] +
\Psi_N(\tilde{\gamma},\beta,1),
\end{eqnarray}
The left hand side gives
\begin{eqnarray}
\textbf{E}[\ln Z_{N+1}(\gamma,\beta)] &-& \textbf{E}[\ln
Z_{N+1}(\tilde{\gamma},\beta)] = \\ \nonumber &=&
(\gamma - \tilde{\gamma})\frac{d}{d\gamma}\textbf{E}[\ln Z_{N+1}(\gamma,\beta)]|_{\gamma=\tilde{\gamma}} = \nonumber \\
&=& \gamma\frac{1}{N+1}\frac{d}{d\gamma}\textbf{E}[\ln Z_{N+1}(\gamma,\beta)]|_{\gamma=\tilde{\gamma}} = \nonumber \\
&=& \gamma \frac{d}{d\gamma}A_{N+1}(\gamma,\beta).
\end{eqnarray}

Considering the $\alpha$ dependence of $\gamma$, we have
$$
\partial_{\gamma} \propto \frac{1}{(p-1)\alpha^{p-2}}\partial_{\alpha} \quad
\Rightarrow \quad \gamma\frac{d}{d\gamma}A \propto
\frac{\alpha}{p-1}\frac{d}{d\alpha}A,
$$
where the symbol $\propto$ instead of $=$ reflects the
arbitrariness by which we include  the $p!$ term, multiplying
$\alpha$, inside the definition of $\gamma$, or directly in
$\alpha$.
\newline
Performing now the thermodynamic limit, we see that at the right
hand side we have \begin{equation} \lim_{N \to \infty}\textbf{E}[\ln
\frac{Z_{N+1}(\alpha,\beta)}{Z_N(\tilde{\alpha},\beta)}]
\longrightarrow A(\alpha,\beta) \end{equation} and the theorem is proved
$\Box$.

\bigskip

\textbf{Proofs of Theorems
\ref{ciccia},\ref{saturabili},\ref{saturi}}
\newline
In this sketch we are going to show how to get Theorem
(\ref{ciccia}) in some details; It automatically has as a
corollary Theorem (\ref{saturabili}) which ultimately gives, as a
simple consequence when applied on filled monomials,
Theorem(\ref{saturi}).
\newline
Let us assume for a generic overlap correlation function $Q$, of
$s$ replicas, the following representation
$$
Q = \prod_{a=1}^s\sum_{i_l^a}\prod_{l=1}^{n^a}\sigma_{i_l^a}^a
I(\{i_l^a \})
$$
where $a$ labels the replicas, the internal product takes into
account the spins (labeled by $l$) which contribute to the a-part
of the overlap $q_{a,a'}$ and runs to the number of time that the
replica $a$ appears in $Q$. The external product takes into
account all the contributions of the internal one and the $I$
factor fixes the constraints among different replicas in $Q$; so,
for example, $Q=q_{12}q_{23}$ can be decomposed in this form
noting that $s=3$, $n^1=1,n^2=2$,
$I=N^{-2}\delta_{i_1^1,i_1^3}\delta_{i_1^2,i_2^3}$, where the
$\delta$ functions fixes the links between replicas $1,3
\rightarrow q_{1,3}$ and $2,3 \rightarrow q_{2,3}$. The averaged
overlap correlation function is
$$
\langle Q \rangle_t = \mathbf{E}\sum_{i_l^a}I(\{i_l^a
\})\prod_{a=1}^s \omega_{t}(\prod_{l=1}^{n^a}\sigma_{i_l^a}^a).
$$
Now if $Q$ is a fillable polynomial, and we evaluate it at $t=1$,
let us decompose it, using the factorization of the $\omega$ state
on different replica, as
$$ \langle Q \rangle_t = \mathbf{E}\sum_{i_l^a,i_l^b}I(\{i_l^a \}, \{i_l^b \})\prod_{a=1}^u
\omega_a ( \prod_{l=1}^{n^a}\sigma_{i_l^a}^a) \prod_{b=u}^s
\omega_b ( \prod_{l=1}^{n^b}\sigma_{i_l^b}^b),
$$
where $u$ stands for the number of the unfilled replicas inside
the expression of $Q$. So we split the measure $\Omega$ into two
different subset $\omega_{a}$ and $\omega_{b}$: in this way the
replica belonging to the $b$ subset are always in even number,
while the ones in the $a$ subset are always odds. Applying the
gauge $\sigma_i^a \rightarrow \sigma_i^a\sigma_{N+1}^a, \forall i
\in (1,N)$ the even measure is unaffected by this transformation
$(\sigma_{N+1}^{2n} \equiv 1)$ while the odd measure takes a
$\sigma_{N+1}$ inside the Boltzmann measure. \begin{eqnarray} &&
\langle Q \rangle = \\ \nonumber && \sum_{i_l^a,i_l^b}I(\{i_l^a
\}, \{i_l^b \}) \prod_{a=1}^u \omega ( \sigma_{N+1}^a
\prod_{l=1}^{n^a}\sigma_{i_l^a}^a) \prod_{b=u}^s \omega (
\sigma_{N+1}^b\prod_{l=1}^{n^b}\sigma_{i_l^b}^b). \end{eqnarray}
At the end we can replace in the last expression the index $N+1$
of $\sigma_{N+1}$ by $k$ for any $k \neq \{ i_l^a \}$ and multiply
by one as $1=N^{-1}\sum_{k=0}^N$. Up to orders $O(1/N)$, which go
to zero in the thermodynamic limit, we have the proof.
\medskip
\newline
It is now immediate to understand that Theorem (\ref{ciccia}) on a
fillable overlap monomial has the effect of multiplying it by its
missing part to be filled (Theorem \ref{saturabili}), while it has
no effect if the overlap monomial is already filled (Theorem
\ref{saturi}). $\Box$

\bigskip

\textbf{Proof of Proposition \ref{stream}}
\newline
The proof works by direct calculation:
\begin{eqnarray}
&& \frac{\partial\langle F_s \rangle_{t,\tilde{\alpha}}}{\partial
t} = \\ \nonumber && \frac{\partial \textbf{E}}{\partial t}  \Big[
\frac{\sum_{\{\sigma\}}F_s
e^{\sum_{a=1}^s(\beta\sum_{\nu=1}^{k_{\tilde{\gamma} N}}
\sigma_{i_{\nu}^1}^a...\sigma_{i_{\nu}^p}^a +
\beta\sum_{\nu=1}^{k_{2\tilde{\gamma}t}}
\sigma_{i_{\nu}^1}^a...\sigma_{i_{\nu}^{p-1}}^a)}}
{\sum_{\{\sigma\}}
e^{\sum_{a=1}^s(\beta\sum_{\nu=1}^{k_{\tilde{\gamma} N}}
\sigma_{i_{\nu}^1}^a...\sigma_{i_{\nu}^p}^a + \beta
\sum_{\nu=1}^{k_{2\tilde{\gamma}t}}
\sigma_{i_{\nu}^1}^a...\sigma_{i_{\nu}^{p-1}}^a)}}\Big] =
\\ \nonumber && 2\tilde{\alpha}^{p-1} \textbf{E}\Big[
\frac{\tilde{\Omega}_t (F_s
e^{\sum_{a=1}^s(\beta\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a)})}
{\tilde{\Omega}_t(e^{\sum_{a=1}^s(\beta\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a)})}\Big]
- 2\tilde{\alpha}^{p-1}\langle F_s \rangle_{t,\tilde{\alpha}} = \\
\nonumber &&  2\tilde{\alpha} \textbf{E}\Big[
\frac{\tilde{\Omega}_t (F_s \Pi_{a=1}^{s}(\cosh\beta +
\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a\sinh\beta))}
{\tilde{\Omega}_t (\Pi_{a=1}^{s}(\cosh\beta +
\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a\sinh\beta))}\Big]- \\
\nonumber &&
2\tilde{\alpha}^{p-1}\langle F_s \rangle_{t,\tilde{\alpha}} = \\
\nonumber &&   2\tilde{\alpha}^{p-1} (\textbf{E}\Big[
\frac{\tilde{\Omega}_t (F_s \Pi_{a=1}^{s}(1 +
\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a\theta))} {(1 +
\tilde{\omega}_t(\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a)\theta)^s}\Big]
- \langle F_s \rangle_{t,\tilde{\alpha}}),
\end{eqnarray}
Now noting that
\begin{eqnarray} \nonumber
\Pi_{a=1}^{s}(1 &+&
\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a\theta) = 1 +
\sum_{a=1}^{s}\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a\theta
\\ \nonumber &+& \sum_{a<b}^{1,s}\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a
\sigma_{i_{0}^1}^b...\sigma_{i_{0}^{p-1}}^b\theta^2
+ ...\nonumber \\
\frac{1}{(1 + \tilde{\omega}_t \theta)^s} &=& 1 -
s\tilde{\omega}_t \theta + \frac{s(s+1)}{2!}\tilde{\omega}_t^2
\theta^2 + ... \nonumber
\end{eqnarray}
\medskip
we obtain
\begin{eqnarray}
\frac{\partial\langle F_s \rangle_{t,\tilde{\alpha}}}{\partial t}
&=& 2\tilde{\alpha}^{p-1} \Big(\textbf{E}\Big[ \tilde{\Omega}_t
\Big(F_s(1 +
\sum_{a=1}^{s}\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a\theta +
\\ \nonumber &+& \sum_{a<b}^{1,s}\sigma_{i_{0}^1}^a...\sigma_{i_{0}^{p-1}}^a
\sigma_{i_{0}^1}^b...\sigma_{i_{0}^{p-1}}^b\theta^2
+ ...)\Big) \times \nonumber \\
&\times& \Big(1 - s\tilde{\omega}_t \theta +
\frac{s(s+1)}{2!}\tilde{\omega}_t^2 \theta^2 + ...\Big)\Big] -
\langle F_s \rangle_{t,\tilde{\alpha}}\Big), \nonumber
\end{eqnarray}
from which our thesis follows. $\Box$

\section*{Acknowledgment}

The authors are grateful to Francesco Guerra, Pierluigi Contucci
and Raffaella Burioni for interesting discussions.

\addcontentsline{toc}{chapter}{Bibliografia}


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

Many solids,
like alloys or doped semiconductors,
form crystals consisting of two or more chemical species.
In contrast to amorphous solids they possess, at least
approximately, a regular lattice whose sites are randomly occupied by
the different components.
Understanding their electronic structure is
an important task which we will address here through an
approach that allows for a comprehensive description of such
substitutionally disordered systems.
This approach, which we call the local distribution (LD) approach,
considers the local density of states (LDOS) $\rho_i(\omega)$,
which is a quantity of primary importance in systems
with prominent local interactions or scattering,
e.g. in the formation of local magnetic moments.
What makes the LD approach `non-standard' is that it directly deals with the
\emph{distribution} of the LDOS in the spirit that Anderson
introduced in his pioneering work.~\cite{anderson58}
While the LDOS is directly related to the amplitude of the electron's
wave-functions on a certain lattice site, its distribution captures the
fluctuations of the wave-functions through the system.

The LD approach has been originally developed
for a description of Anderson localization~\cite{abou73},
and furthermore been applied to systems with topological
disorder.~\cite{logan84,logan87}
In this article we like to demonstrate that it can also describe
systems with a bimodal disorder distribution, as the binary alloy model.
It has been noted several times
(see e.g. Refs.~\onlinecite{kirk72,taka74,souk92}), that
the physics of this model is only partially covered in a mean field
description as provided by the coherent potential approximation
(CPA).
Instead one has to account for the specific configuration of atoms in the
vicinity of a lattice site.~\cite{dean58,tsukada69}
To give a better description, a variety of
(cluster) extensions to CPA has been devised,
which explicitly treat correlations on finite clusters (see
e.g. Ref.~\onlinecite{taka74} or Ref.~\onlinecite{elliott74}).
We will show that the LD approach can serve as a conclusive extension of
CPA---indeed it contains
the CPA as a limit, see App.~\ref{appendix:LD}---which implicitly
contains these correlations.

Besides giving very precise results for non-interacting systems,
an appealing feature of the LD approach
is its possible application to interacting disordered systems.
Recently some progress in this direction has been made using the LD
approach in combination with dynamical mean field theory.~\cite{dobro98}
We could e.g. show that polaron formation in an
electron-phonon-coupled system is enhanced in the presence
of impurities, leading to polaron like defect states.~\cite{baf04}
The Mott transition in a binary alloy is another example,
where the LD approach might help to substantiate present
results.~\cite{byczuk04}

The outline of this article is as follows:
We will shortly introduce the binary alloy model and the associated
distributions, and then apply the LD approach.
As a limiting case we will consider the (quantum) percolation model,
and finally conclude. The two appendices contain the derivation of the
LD approach and its application to the Anderson localization problem.

\section{Model and distributions}\label{sec:modeldis}

The simplest model for an electron moving in a crystal with
substitutional disorder is given by the Hamiltonian
\begin{equation}\label{eq:hamilton}
H=\sum\limits_i \epsilon_i c^\dagger_i c^{\phantom{\dagger}}_i
- t
\sum\limits_{\langle i,j \rangle }
c^\dagger_i c^{\phantom{\dagger}}_j \; ,
\end{equation}
where $t$ denotes the tight-binding hopping integral between nearest
neighbour sites on a given lattice,
and the $\epsilon_i$'s are random on-site potentials.

The model is further specified through the probability
distribution of the random variables $\epsilon_i$.
We only consider models where the
$\epsilon_i$ are identically
independently distributed random variables
with a fixed distribution $p(\epsilon_i)$.
The Anderson model of localization, the binary alloy
model and the quantum percolation model are examples.

For given values of the $\epsilon_i$,
i.e. for a specific disorder configuration (`one specimen'),
the local density of states (LDOS)
\begin{equation}\label{def:LDOS}
\rho_i(\omega)=- \mathrm{Im}\; G_{ii}(\omega)/\pi
\end{equation}
 (expressed through the retarded
Green function $G_{ii}(\omega)$)
has a definite value on each lattice site $i$.
For the ordered system, translational symmetry implies that
$\rho_i(\omega)$ has the same value for every $i$.
In a disordered system, $\rho_i(\omega)$ varies with $i$,
so it takes on different values with a certain probability.
The corresponding probability distribution $p(\rho_i,\omega)$
captures the fluctuation of the LDOS through the system.
It is reasonable to assume that for an infinite system $p(\rho_i,\omega)$
does not depend on the explicit values of the $\epsilon_i$, but only
on the probability distribution $p(\epsilon_i)$.
We thus assume that the distribution
$p(\rho_i,\omega)$ is, in a slight `abuse of language',
self-averaging: It does not depend on the specific disorder configuration
looked at, but takes on a definite value in the thermodynamic limit.
Moreover  $p(\rho_i,\omega)$ has a second meaning.
If we look at a fixed lattice site $i$
but consider all possible values of the $\epsilon_i$,
the LDOS is a random variable in its own right, whose probability
distribution is equal to $p(\rho_i,\omega)$ as defined above.
So $p(\rho_i,\omega)$ gives
(i) the probability that the LDOS has a certain value on some lattice
site, when all $\epsilon_i$ are fixed, and (ii) the probability
that the LDOS has a certain value on a fixed lattice site $i$, when
the $\epsilon_i$ vary.
Note that $p(\rho_i,\omega)$ does not depend on $i$, in contrast to $\rho_i(\omega)$.

One basic physical quantity that is calculated from
$p(\rho_i,\omega)$ is the
\emph{arithmetically averaged density of states} (DOS)
\begin{equation}
 \rho(\omega) =  \int\limits_0^\infty
 \rho_i \; p(\rho_i,\omega) \; d\rho_i
\end{equation}
which counts the number of states at energy~$\omega$.
For the ordered system, when $\rho_i(\omega)$ does not depend on $i$,
$p(\rho_i,\omega)=\delta(\rho_i-\rho(\omega))$
is a $\delta$-peak at the DOS.
In the presence of disorder
$p(\rho_i,\omega)$ attains a certain width
and provides additional information on the
character of the electronic states in the system.
A narrow distribution corresponds to more or less homogeneous states,
when electron scattering is weak,
while a broad distribution reflects strong scattering leading to
very inhomogeneous states.
The essential information on the character of states at energy $\omega$ is thus
provided by the distribution $p(\rho_i,\omega)$, which gives
a more detailed description of the disordered system
than the DOS $\rho(\omega)$ alone.

In particular we can decide whether states at $\omega$ are localized
if we employ the precise definition of $p(\rho_i,\omega)$.
The retarded Green function $G_{ii}(\omega)$, hence the LDOS and its
distribution, is usually calculated for a complex energy
$\omega+\mathrm{i}\eta$ with a small positive imaginary part $\eta$,
followed by analytical continuation to the real axis, i.e. $\eta\to0$.
A finite $\eta$ gives a Lorentzian broadening with respect to
$\omega$, i.e. a finite energy resolution.
For $\eta\to0$, the resolution increases until peaks in the LDOS
(corresponding to poles of $G_{ii}(\omega)$) and bands
(corresponding to branch cuts of $G_{ii}(\omega)$) can be separated.
The behaviour of the distribution in the limit $\eta\to0$ is thus different depending
on the spectral properties of the system.
For a spectrum consisting only of discrete peaks---this corresponds
to localized states---the distribution becomes singular for $\eta\to0$.
Contrary to this we get a regular distribution if the spectrum is
continuous as for extended band (`Bloch') states.
We use this property in the study of Anderson localization (see
App.~\ref{appendix:anderson}), and will use it for the fragmented
spectra of the quantum percolation model.

In this article we obtain $p(\rho_i,\omega)$ directly through the
LD approach (for details see App. \ref{appendix:LD}).
It is natural to compare the results for $\rho(\omega)$
of the LD approach to the corresponding CPA results.
While the lattice enters the CPA calculations only through the lattice
DOS for the ordered system the LD approach is explicitly constructed for a
Bethe lattice~\footnote{
The (half-infinite) Bethe lattice is a loop-free tree
with a semicircular DOS
$\rho(\omega)=(4/\pi W^2)\sqrt{W^2-4\omega^2}$.~\cite{economou83}
If not stated otherwise, we use a Bethe lattice with
$K=2$ nearest neighbours to any lattice site.
The hopping matrix element $t$ is chosen to give a full
bandwidth $W=4 t \sqrt{K}=1$, so energies will be measured in units of $W$.
}.
Owing to the absence of loops the Bethe lattice is a kind of
mean-field approximation to (hyper-) cubic lattices.
Its particular geometry does not affect the
qualitative behaviour away from the localization transition, which
indeed is similar to a cubic lattice in three dimensions, see e.g. our
discussion of the quantum percolation model (Sec.~\ref{sec:perc}) or
the phase diagram obtained for the Anderson model
(App.~\ref{appendix:anderson}). Conversely the existence of different
phases and transitions in between is correctly described by the LD
approach for dimensions~$\ge 3$, similar to the success of dynamical
mean field theory for interacting systems in high dimensions.
Of course at the phase transition, when the correlation length
diverges, results for Bethe and cubic lattices differ
(see Ref.~\onlinecite{mirlin94} and App.~\ref{appendix:anderson}).
Weak localization in two dimensions, on the other hand,
is related to interference on closed loops
and thus not seen on a Bethe lattice.

\section{Binary Alloy Model}\label{sec:bam}

The binary alloy model has a bimodal disorder distribution
\begin{equation}
p(\epsilon_i)=c_A \delta(\epsilon_i-E_A)+(1-c_A)
\delta(\epsilon_i-E_B) \quad,
\end{equation}
corresponding to a crystal randomly composed of A-atoms (B-atoms) at energy
$E_A$ ($E_B$).
We set $E_A=-\Delta/2$, $E_B = \Delta/2$, i.e. $\Delta=E_B-E_A$.
Different aspects of the involved physics of the binary alloy model
have been discussed previously (see
e.g. Refs.~\onlinecite{kirk72,taka74,souk92,dean58,tsukada69}).
The local electron motion
strongly depends on the configuration of atoms at (and in the vicinity of) a
certain lattice site.
This `cluster effect', showing up as peaks in the DOS, is not obtained
by the CPA.
We will show that the LD approach can give a (more) complete picture with
fair ease. Especially the DOS becomes accessible and will
turn out to be important for a thorough understanding.

\subsection{Low A-concentration -- small separation energy}\label{subsec:bam1}
We first study the situation of low concentration
$c_A=0.1$ of A-atoms (which is below the classical percolation
threshold, cf. Sec. \ref{sec:perc}).
The A-atoms can be considered as the minority species in a bulk system
of B-atoms (doping a semiconductor is an example).

For small separation energy $\Delta=0.3$
the energy levels of the minority A-atoms
lie inside the B-band.
The DOS mainly consists of the B-band
centered at $E_B=+0.15$, but shows some additional structures at the lower
band edge,
which are absent in the CPA DOS
(cf. Fig.~\ref{fig:badosprob}).
\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.85\linewidth]{Fig1.eps}
\includegraphics[width=0.85\linewidth]{Fig2.eps}
\end{center}
\caption{
Binary alloy with $\Delta=0.3$, $c_A=0.1$.
The upper panel shows the  DOS, obtained within the LD approach and CPA.
The inset displays the conditional  DOS for the cluster configuration
$A-BB$, the dashed curve has been calculated using Eq.~(\ref{singleA}) with
broadening $\eta=10^{-2}$ (see text).
The lower panel shows the distribution of the LDOS at $\omega=0.6$
(curve on top).
The six other curves show the distribution of the conditional LDOS for
a specified cluster of three sites.}
\label{fig:badosprob}
\end{figure}
They originate from the strong fluctuations in the local
environment of an atom,
which can be clearly seen in the distribution of the LDOS (cf. Fig.~\ref{fig:badosprob}).
Note the various peaks of this distribution (in contrast to the
Anderson model, see App.~\ref{appendix:anderson}), and its large width
displaying the pronounced fluctuations of the electron's wave-functions.

To understand the consequences of different atoms
situated at a lattice site
we can look at the `conditional' LDOS
$\rho^{A/B}_i(\omega)$,
subject to the constraint that an $A$-
respectively $B$-atom is located at site $i$.
More generally we can specify a certain configuration of atoms on
a cluster of sites centered at $i$.
By considering larger and larger clusters
every peak in the distribution of the conditional LDOS (and in the DOS) can be attributed to a specific configuration
(cf. Fig.~\ref{fig:badosprob}).
As the simplest example the pronounced peak at
$\omega \approx -0.35$ results from a single A-atom surrounded by B-atoms
(cf. inset Fig.~\ref{fig:badosprob}).
Its approximate position and form
follow from the simple formula
\begin{equation}\label{singleA}
  G_A(\omega)=\left[\frac{\omega+E_B}{2}-E_A+
\sqrt{\left(\frac{\omega-E_B}{2}\right)^2-\frac{W^2}{16}\;}
\; \right]^{-1}
\end{equation}
for the local Green function of one A-impurity embedded in a
B-lattice.
The corresponding DOS, using a Lorentzian broadening
$\eta=0.02$ in the energy argument $\omega+\mathrm{i}\eta$ of the
Green function to mimic the effect of tunnelling between different
A-atoms, fits the conditional DOS very well (dashed curve in
inset of Fig.~\ref{fig:badosprob}).
The complementary situation of a B-atom surrounded by A-atoms contributes to
the `hump' in the B-band (cf. Fig.~\ref{fig:badosprob}),
which arises from B-atoms neighbouring to A-atoms.

It should be noted that Eq.~(\ref{singleA})
gives a $\delta$-peak outside the B-band,
corresponding to the impurity state at the A-atom,
only above a critical value of $\Delta$.
This is another evidence that the Bethe lattice should be understood as
an approximation to lattices in dimensions $d \ge 3$.

\subsection{Low A-concentration -- split band case}\label{subsec:bam2}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.85\linewidth]{Fig3.eps}
\end{center}
\caption{
  DOS for the binary alloy with $\Delta=2.0$, $c_A=0.1$.
  To broaden the peaks of the `A-states'
  we add an artificial imaginary part $\eta=10^{-3}$ to the energy
  argument $\omega$ of the Green function $G_{ii}(\omega)$ in
  Eq.~(\ref{def:LDOS}).
  Arrows mark the position of three peaks each in the A- and B-band.}
\label{fig:basplit}
\end{figure}

Increasing the separation energy to $\Delta=2.0$
(still with $c_A=0.1$)
leads to the split band case (see Fig.~\ref{fig:basplit}).
The majority B-atoms still form a band whose DOS shows some additional spikes.

The concentration of the minority A-atoms is below the classical percolation
threshold hence only finite A-clusters exist (cf. Sec.~\ref{sec:perc}).
Due to the large energy separation the `A-states' on these clusters
are strongly damped through scattering on B-atoms.
Accordingly the `A-states' do not form a band of (extended)
states but a fragmented set of peaks with varying
height and width, reminiscent of the percolation model (cf. Sec. \ref{sec:perc}).
Again these peaks emerge from A-clusters embedded in the B-lattice.
The central peak at $\omega=-1.031$
whose position can be again calculated with Eq.~(\ref{singleA})
corresponds to a single A-atom surrounded by B-atoms.
Similarly the two side peaks result from two adjacent A-atoms, and so forth.
The weight of these peaks decreases exponentially as $c_A^N$
for a $N$-atomic cluster.
The complementary configurations, reversing the role of A- and B-atoms,
yield spikes in the B-band.
Since the concentration of B-atoms is large these spikes do not form
peaks but merge with the B-band.

Owing to the different strength of scattering we expect that
  all `A-states' are localized, while the `B-band' remains extended.
In accordance with localization in the Anderson model
(App.~\ref{appendix:anderson}) we can extract the nature of states from
the behaviour of the distribution $p(\rho_i,\omega)$
for $\omega$ shifted by a small imaginary part $\eta$.
Indeed we find that, for $\eta\to0$,
the distribution is singular in the energy range of the
`A-states' but regular for the B-band,
which proves that `A-states' are localized and the `B-band' remains
extended.

\subsection{Phase diagram}\label{subsec:bam3}

In Fig.~\ref{fig:baphase} we show the phase diagram of the binary
alloy for low concentration $c_A=0.1$.
For $\Delta$ below a critical value $\Delta_c \approx 0.5 $ the spectrum
consists of a single band, which splits into two bands above $\Delta_c$.
The value of $\Delta_c$ is a little bit larger than obtained within CPA
($\Delta_c^\mathrm{CPA} \approx 0.45$), but significantly smaller than
$\Delta =1$, when two bands with bandwidth $W=1$ centered around
$\pm \Delta/2$ would no longer overlap.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.85\linewidth]{Fig4.eps}
\includegraphics[width=0.85\linewidth]{Fig5.eps}
\end{center}
\caption{\label{fig:baphase}
(Colour online)
Phase diagram of the binary alloy,
for low concentration $c_A=0.1$ (upper panel),
and equal concentration $c_A=c_B=0.5$ (lower panel).
  The solid filled curves show the DOS for various $\Delta$.
  The dashed curves show the CPA band edges,
  and the dotted lines mark $\omega=\pm \Delta/2 \pm W/2$.}
\end{figure}

The phase diagram shows the rich physics of the binary alloy As
discussed in the previous two sections:
(i) formation of many gaps, leading to a strongly fragmented minority `band',
(ii) the bands between the gaps show additional spikes,
(iii) for split bands and low concentrations,
minority states are localized impurity states,
while majority states are extended band states.
These three important effects go beyond the CPA,
which gives only the band splitting.
The CPA works reasonably well for the majority band,
where the local environment of a lattice site is less important,
but it cannot describe the fine structure of the minority band in even
a crude way.
These findings are in full agreement with Ref.~\onlinecite{souk92},
and we can additionally show the exact form of the DOS which
is one outcome of the LD approach.
So far we have dealt with a low concentration $c_A$ of the
A-species below the classical percolation threshold $p_c = 1/K$ for
the Bethe lattice.
Then, when only finite A-clusters exist, strong signatures and
fragmentation of the DOS could be observed.
Increasing $K$ with $c_A$ fixed,
these signatures are partially weakened as long as $c_A < p_c$
(see Fig.~\ref{fig:IncK}).
Above the percolation threshold ($c_A > p_c$), infinite A-clusters
exist, leading to a band which is no longer fragmented.
Nevertheless the DOS still has peaks.
For large $K$ (equivalent to $p_c \ll c_A$) the DOS approaches a
semicircular form, anticipating the CPA result in the limit $K=\infty$
(cf. the discussion in the appendices).

For equal concentration $c_A=c_B=0.5$ and $K=3$, both atom species are
in the percolating regime, which corresponds not to a
doped material but to a stoichiometric compound.
The phase diagram (Fig.~\ref{fig:baphase}) shows that
the signatures in the DOS become more pronounced with increasing
$\Delta$, when the binary alloy model approaches
the percolation model (see Sec.~\ref{sec:perc}).

\begin{figure}
\begin{center}
\includegraphics[width=0.85\linewidth]{Fig6.eps}
\end{center}
\caption{\label{fig:IncK}
  DOS of the impurity band for $\Delta=2.0$ and $c_A=0.2$,
  with increasing number of nearest neighbours $K$ of the Bethe
  lattice.
  For $K=2,4$, the concentration $c_A$ is below the classical
  percolation threshold $p_c=1/K$, for $K=6,12$ above.
  In contrast to Fig.~\ref{fig:basplit} we did not broaden the peaks
  with an additional $\eta$. }
\end{figure}

\section{Combined Anderson and binary alloy model}\label{sec:cab}
The binary alloy model shows a distinct tendency towards
peak and gap formation.
This behaviour is a generic feature of systems with a bimodal disorder
distribution.
For example consider the binary alloy
model with additional on-site box disorder $\gamma$, i.e.
\begin{equation}\label{eq:bap}
p(\epsilon_i) = \frac{c_A}{\gamma} \Theta(\gamma/2-|\epsilon_i-E_A|)
+ \frac{1-c_A}{\gamma} \Theta(\gamma/2-|\epsilon_i-E_B|) \; .
\end{equation}
The CPA suggests that the system is the combination of two rescaled
Anderson models.
This is consistent with its phase diagram concerning the
position of mobility edges.~\cite{plyu03}
The LD approach additionally shows that
structures similar to the `pure' binary alloy appear
(cf. Fig.~\ref{fig:bap}).
Evidently the system is not adequately described in terms
of rescaled Anderson models.
For the parameters considered $\gamma$ is on the same energy
scale as the hopping matrix element $t$
(here $t=W/\sqrt{32}$, i.e. $\gamma \approx 0.56 t$).
It is thus reasonable to assume that in a real alloy,
e.g. a doped semiconductor,
the peaked structure of the binary alloy DOS can be found.
Furthermore, if we extract localization properties of the system by the
$\eta\to0$-limit (cf. App.~\ref{appendix:anderson}), we find the A-band
to be entirely localized (see Fig.~\ref{fig:bap}).
Nevertheless $\gamma$ is much smaller than the critical disorder
($\gamma_c \simeq 3 \times \mathrm{bandwidth}$) for Anderson
localization of the A-subband (in agreement with Ref.~\onlinecite{plyu03}).
Therefore localization of the A-band is increased due to the strong
scattering occurring for the minority band in the binary alloy.
This suggests that impurity band states in a doped semiconductor will almost
always be localized.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.85\linewidth]{Fig7.eps}
\end{center}
\caption{\label{fig:bap}
The upper panel shows the DOS for the
binary alloy ($\Delta=2.0$, $c_A=0.1$) with additional on-site disorder
$\gamma=0.1$. (see Eq.~\ref{eq:bap}).
The bottom row displays the dependence on $\eta$ of the LDOS
distribution.
Lower left panel: In the center of the A-band at $\omega=-1$
the distribution becomes singular for $\eta\to0$,
i.e. states are localized.
Lower right panel: For the B-band at $\omega=+1$
the distribution is independent of $\eta$,
and has a peaked structure known from the binary alloy model
(cf. Fig.~\ref{fig:badosprob}). }

\end{figure}

\section{Quantum Percolation Model}\label{sec:perc}

In the limit $E_B \to \infty$, now keeping $E_A=0$ fixed, the
binary alloy model reduces to the percolation model
where A-sites are embedded in an impenetrable host medium~\footnote{On
  the Bethe lattice site and bond percolation are equivalent.}.
Certainly, in a real alloy with $\Delta < \infty$,
the electron has the chance to tunnel through the host barrier.
Nevertheless, the percolation model
shows very general features
which have already shown up in the split band case of the binary alloy
model (see above).

Let us first study the DOS.
If the concentration $c_A$ of A-atoms is below the classical
percolation threshold $p_c$ only finite clusters exist.
The corresponding spectrum is a pure point spectrum which densely fills
the energy interval $[-W/2,W/2]$ of the tight-binding band
(`Dirac comb').
An $N$-site cluster with occurrence probability
$c_A^N(1-c_A)^{N+1}$ contributes peaks at energies spread
over the full possible range.
Consequently the weight of a peak varies
in contrast to the Anderson model non-monotonically with $\omega$.
Different statistical properties of the finite clusters can be
straightforwardly calculated on the Bethe lattice.
For instance the weight
\begin{eqnarray}
  w_\mathrm{fin} &=& \sum\limits_{N=1}^\infty c_A^{N-1}(1-c_A)^{N+1}
  C_N   \\
  &=&
\begin{cases}
 1, \quad &    c_A \le 0.5 \\
 (1-c_A)^2/c_A^2, &   c_A > 0.5
\end{cases}
\end{eqnarray}
of all finite clusters~\footnote{Again on a Bethe lattice with
  coordination number $K=2$, cf. first footnote.}
directly follows
by help of the generating function for the Catalan numbers $C_N$, which give
the number of binary trees with $N$ sites.
For concentrations above $p_c \; (=0.5)$, when $w_\mathrm{fin}<1$,
an infinite percolating cluster exists.
This cluster can support extended states
which contribute to an continuous spectrum.
Since the `effective' dimension of the percolating cluster
is smaller than $K$ (especially close to $p_c$),
the bandwidth of the resulting band is smaller than $W$.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.85\linewidth]{Fig8.eps}
\end{center}
\caption{\label{fig:percdos}
(Colour online)
  DOS $\rho(\omega)$ for the percolation
  model,
  for various concentrations $c_A$.
  Upper panel: Full spectrum with Lorentzian broadening $\eta=10^{-3}$.
Lower left panel: Central region around $\omega=0$ without Lorentzian
broadening.
Lower right panel: Region around the second peak at $\omega = t =
1/\sqrt{32} \approx 0.177$  without Lorentzian
broadening. }
\end{figure}

The DOS shows strong signatures (peaks and dips),
which have already occurred in the binary alloy model (see top picture in
Fig.~\ref{fig:percdos}).
These signatures arise from both isolated finite clusters and
finite clusters attached to the backbone of the infinite
percolating cluster.
With increasing concentration their weight reduces,
and the signatures are washed out.
At $\omega=0$ a pronounced $\delta$-peak surrounded by a
dip in the DOS exists.
With increasing concentration the peak reduces in weight and the dip
narrows, eventually both merge to a spike in the band.

The origin of this central peak and the dip
can be understood on the same level of reasoning as for
Eq.~(\ref{singleA}).
If a single atom is attached to the backbone of the percolating
cluster the Green function is modified by the additional hopping to this atom.
The DOS then shows a $\delta$-peak at
 $\omega=0$, which is the energy of the state located at the atom,
and a dip around this peak which arises from damping of states on
the percolating cluster.
The same argumentation holds for any finite cluster instead of a
single atom.
Since larger clusters have lower probability of occurrence the central
peak is the most pronounced, and exists up to the largest
concentrations.

Nevertheless the percolating cluster is not made up of its backbone
plus one additional finite cluster, but many finite clusters attached
to it.
The Green function $G^c(\omega)$ of a (half-)infinite chain with one
additional atom attached to each site obeys the recursion relation
\begin{equation}
G^c(\omega)=(\omega - t^2 /\omega - t^2 G^c(\omega)  )^{-1}
\end{equation}
which is solved by $G^c(\omega)=G^0(\omega-t^2/\omega)$,
where $G^0(\omega)$ is the Green function of the Bethe lattice.
Of course $G^c(\omega=0)=0$ as before,
but the diverging real part of $1/\omega$ produces not a dip but a \emph{gap}
in the DOS around $\omega=0$.
The validity of this simple argumentation
can be tested within the LD approach since its energy resolution
is not limited by finite size effects.
Without Lorentzian broadening, i.e. $\eta\to0$,
the $\delta$-peaks arising from finite cluster states
do not contribute to the DOS
(cf. Sec.~\ref{sec:modeldis}),
and only the continuous spectrum from extended band states survives.
This spectrum indeed shows a gap around $\omega=0$ for sufficiently
small concentrations (cf. Fig.~\ref{fig:percdos}, bottom left).
Lowering the concentration, gaps open at the energies corresponding
to any finite cluster eigenstate.
The first additional gaps open at $\omega=\pm t$
(Fig.~\ref{fig:percdos}, bottom right), corresponding to
two site clusters, which have the second largest weight among all finite
clusters.
These gaps are filled by peaks from the dense spectrum of finite
cluster states, thus are absent in the DOS for any
finite energy resolution (cf. Ref.~\onlinecite{schubert05}).
The formation of gaps---which is a kind of `level repulsion'---is a
significant quantum feature with no counterpart in the classical
model.

Besides the DOS the nature of states is
important.
For concentrations $c_A < p_c$, when only finite clusters exist,
all states are localized and no electron transport is possible.
For concentrations above $p_c$ a classical electron can traverse
the system along the percolating cluster. But a quantum mechanical
electron scatters off all irregularities and is possibly localized.
This raises the question of a quantum percolation threshold $p_q$
\emph{above}~$p_c$.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.85\linewidth]{Fig9.eps}
\end{center}
\caption{\label{fig:perchisto}
LDOS distribution for the percolation model at $\omega=0.1$.
  Left panel: For $c_A=0.7$ and two values of regularization
  $\eta=10^{-8}$ and $\eta=10^{-20}$.
  For $\eta\to0$
 only extended states contribute to $p(\rho_i,\omega)$ at finite $\rho_i$.
  Right panel:
  Integrated probability
  \mbox{$p(\rho_i > \xi,\omega)$}
 in dependence on $\eta$,
  for various $c_A$ and one particular value $\xi=10^{-10}$.
  For $c_A \lesssim 0.54$ no states exist at $\omega=0.1$,
  and $p(\rho_i,\omega)=0$. }
\end{figure}

Close to $p_c$ states on the percolating cluster have only
exponentially small weight. It is thus difficult to extract the
information about those states from the background of finite cluster
contributions.
However, we can gain some information from the distribution of the
LDOS, if $c_A$ and $\eta$ are varied.
Again, if we let $\eta\to0$, the energy resolution increases until
contributions from extended states on
the percolating cluster are separated from the discrete peaks of
localized states (see left panel of Fig.~\ref{fig:perchisto}).
For fixed energy $\omega$ and suitably large concentration $c_A$
the distribution $p(\rho_i,\omega)$ has finite weight at
finite $\rho_i$ in the limit $\eta\to0$ (see right
panel of Fig.~\ref{fig:perchisto}), which
indicates that extended states exist.
Their weight can be estimated through the integrated probability
$p(\rho_i > \xi,\omega)=\int_\xi^\infty p(\rho_i,\omega) \;d\rho_i$,
with $\eta < \xi$.
The weight decreases with $c_A$,
and abruptly drops to zero at a certain concentration (for $c_A
\approx 0.54$ at $\omega=0.1$).
Then no states exist at $\omega$, and $\rho_\mathrm{ave}(\omega)=0$.
This behaviour has to be attributed to the fragmentation of the spectrum
discussed before.
At a given energy $\omega$ a gap will open for concentrations $c_A$
sufficiently close to $p_c$.
Before a gap opens at $\omega$,
states are extended.
After the gap has opened, no states at $\omega$ exist.
Between the gaps, which open at different
concentrations for different energies,
 extended states can survive even for very small
$c_A$, although all states are extremely damped and practically localized.
A definite localization transition from extended to localized states
does however not take place, and the only transition occurs when the
spectrum is fully fragmented, that is at the classical transition $c_A=p_c$.
We conclude, on the basis of these arguments and our numerical
results, that, for the percolation problem on the Bethe lattice, a
quantum percolation threshold above the classical one does not exist.
Note that scattering on the percolating cluster is of
a different type than for the Anderson model.
The finite clusters attached to the percolating backbone do not act
as coherent but incoherent scatterers.
So states on the backbone will not be localized for small $c_A$
(i.e. strong scattering), or even immediately localized, as in one
dimension.

\section{Conclusions}\label{sec:concl}

In this article we demonstrated how the
LD approach can be used to study
localization and percolative effects in alloys.
With very moderate computational demands
this scheme suffices to resolve the rich structures
in the DOS originating from comparably strong disorder
fluctuations.
Even for the extreme limit of the binary alloy model,
the percolation model, convincing results are easily obtained.
For instance the question whether gaps form in the percolation model
could be definitely answered,
and the possibility of a quantum percolation threshold above the
classical one could be almost definitely ruled out for the Bethe lattice.

We conclude that the LD approach is a convenient framework
for investigations of disorder and localization, and suggest
its application to interacting disordered systems.

redpajama_set_name: RedPajamaArXiv