\section{Introduction}

In recent years investigation of the
unusual behaviour of random heteropolymers with
a random sequence of different kinds of monomers
has attracted considerable attention (see \cite{3}-\cite{5} and
references therein).
With temperature decreasing such systems experience two
types of phase transitions
at which polymer either goes to freezing
state or has a microphase separation
that is inherent in  protein biomolecules type.

At folding process heteropolymer chooses out
of infinite set of possible states the  unique frozen
non-equilibrium but stationary state  with fixed
spatial conformation and monomer sequence (type of DNA).
The problem in considering of such transitions is to build up more
simple method of averaging thermodynamic values over
the quenched disorder.
Here, the methods of statistical physics first used in spin glass
investigation are taken as  basic ones \cite{6}.

Certain of polymers type of block copolymers, undergo the
so-called  microphase separation at which blocks of
different chemical composition segregate to microdomain structures to
gain a spatially inhomogeneous composition throughout the pattern.
The chemical  bonds between the blocks prevent the usual
macroscopic phase separation observed in homopolymer blends.
A variety of microphases,
such as one-dimensional lamellar structures,
hexagonal arrays of cylinders, body-centered-cubic
arrays of spheres and so on have been observed
for simple diblock copolymer melts.

At first, microphase separation had been studied
within mean-field approximation \cite{7} in the case of  A-B block
copolymers with arbitrary composition $f$.
It turns out that
 at $f\not= 0.5$ with temperature decrease the system
undergoes a sequence
of the microphase separations  of the first order,
where period
$2\pi/k_0$ of space structure is of order of block length
and does not depend on temperature.
In \cite{Braz} it was shown a
principal role of the fluctuations in phase transition.
By virtue of condition $k_0\not=
0$  the order parameter fluctuations give divergent
contribution to thermodynamic values, so that continuous
phase transition transforms to the first order weak transition
\cite{8}.

At passage to random heteropolymers both the microphase separation
and the freezing
remain, but here the space period acquires
 strong dependence on temperature \cite{3,9}.
According to field considerations \cite{10,11} the fluctuations
suppress both mentioned transitions in random copolymers,
but copolymer melts in disordered media can suffer these transformations.

Similar to spin glass \cite{6}, the main approach
of statistical mechanics of
random heteropolymers is the replica method
 \cite{3} that is not only pragmatic,
but very non-descriptive and artificial.
Apart from the replica  trick, the transfer matrix method
\cite{12}, kinetic  approach \cite{13} and others
(see \cite{14})
are used in random heteropolymers theory as well.
However, it is a well known
in theory of spin glass that within the framework of
Sherrington-Kirkpatrick model \cite{15} both replica method and
supersymmetry (SUSY) field approach  give identical results
\cite{16}.
On the other hand, the advantage of the SUSY field usage is
that  its components have
explicit physical meaning \cite{17}.

By recent moment, SUSY method in polymer theory had been proposed
by Vilgis \cite{18}, but did not obtain further development.
In this work we intend to complete this problem.
Our article is organized as follows.
In Section 2 we build up an effective motion equation
of disorder heteropolymer to describe the
different monomer alternation along the heteropolymer chain
\cite{20}.
The form of effective Hamiltonian of the problem under consideration
is obtained in Section 3 to avoid discrepancies of different methods.
Section 4 is devoted to consideration of a SUSY scheme \cite{21}
that is optimal for study of the disorder heteropolymer
in weak segregation limit \cite{22}.
In Section 5 a SUSY correlator is considered to examine
memory and non-ergodicity effects (Section 6).
Finally, Section 7 is devoted to discussion of obtained results.

\section{Effective motion equation}

Usually, the SUSY method is based on the dynamic
Langevin-type motion equation  \cite{21}.
But covalent bonds in the polymer chain
make such way inefficient because dynamic theory in polymers is
much more complicated, than statistical mechanics of usual many-body
systems \cite{23}.
Therefore, it is essential to
obtain the effective motion equation
instead of dynamic  one.

Let us initiate with directed polymer that is represented by
 the Gaussian
chain with probability $\Psi({\rm{\bf R}},N)$ to find the
end-point $N$ at coordinate ${\rm{\bf R}}$. As is known the
function $\Psi({\rm{\bf R}},N)$ obeys the Schrodinger-type
equation with imagine time $-{\rm i}N$ \cite{23}

\begin{equation}
\partial\Psi/\partial N=
\left (D\partial^2/\partial{\rm{\bf R}}^2-U({\rm{\bf R}},N)\right)\Psi,
\label{1}
\end{equation}
where number $N\gg 1$, $D\equiv b^2/6$ is effective diffusion coefficient
determined by the Kuhn segment length $b$, $U({\rm{\bf
R}},N)$ is external field.
In the limit $N\to\infty$ a solution of
Eq.(1) can be represented in the form of functional integral over
dependence ${\rm{\bf r}}(n)$ of chain coordinate on
number of internal monomers:

\begin{equation} \Psi({\rm{\bf
R}},N)=\int\exp(-S_{{\rm{\bf R}}N}\{{\rm{\bf
r}}(n)\}/2D)\delta{\rm{\bf r}}(n).
\label{2} \end{equation}
Here conventional action $S({\rm{\bf R}},N)\equiv S_{{\rm{\bf R}}
N}\{{\rm{\bf r}}(n)\}=\int_0^N L_0({\rm{\bf r}}(n)){\rm d}n$
corresponds to the fixed end-points ${\rm{\bf r}}(0)={\rm{\bf 0}}$,
and ${\rm{\bf r}}(N)={\rm{\bf R}}$ is  determined by Lagrangian
of the Euclidian field theory \cite{21,23}

\begin{equation} L_0={1\over
2}\left ({{\rm d}{\rm{\bf r}}(n)\over {\rm d}n}\right
)^2+2DU({\rm{\bf r}},n),
\label{3} \end{equation}
where the effective
kinetic energy, for which the continuum limit ${\rm{\bf
r}}(n+1)-{\rm{\bf r}}(n)\to{\rm d}{\rm{\bf r}}(n)/{\rm d}n$ is taken
into account, presents the covalent bonds between
monomers of polymer chain \cite{23}.
Inserting of ansatz (\ref{2}) into Eq.(\ref{1})
gives  the Jacobi-type equation

\begin{equation} {\partial
S\over\partial N}= D{\partial^2 S\over\partial{\rm{\bf R}}^2}-{1\over
2}\left ({\partial S\over\partial {\rm{\bf R}}}\right )^2+2DU.
\label{4} \end{equation}
After introducing generalized momentum
${\rm{\bf p}}\equiv\partial S/\partial{\rm{\bf R}}$ and total
derivative ${\rm d}{\rm{\bf p}}/{{\rm d}N}\equiv\partial{\rm{\bf
p}}/\partial N+({\rm{\bf p}}\partial/\partial{\rm{\bf R}}){\rm{\bf
p}}$, Eq.(4) takes the linear form of Burgers equation:

\begin{equation} {\rm d}{\rm{\bf
p}}/{\rm d}N=D\left (\partial^2{\rm{\bf p}}/\partial{\rm{\bf
R}}^2+2\partial U/\partial{\rm{\bf R}}\right ).  \label{5}
\end{equation}
The above pointed relation of
Eqs.(\ref{1})--(\ref{5}) is a well-known fact in
the theory of directed polymers,
stochastic growth and kinetic roughening phenomena (see \cite{24}).

The basic observation  for our purpose is that the
Schrodinger-type equation (\ref{1}) takes the form of the
Fokker-Planck equation \cite{25}

\begin{equation} {\partial P\over\partial
N}=\left (D{\partial^2\over\partial{\rm{\bf
R}}^2}-{\partial\over\partial{\rm{\bf R}}}{\rm{\bf F}}\right )P
\label{6} \end{equation}
for distribution probability

\begin{equation}
P({\rm{\bf R}},N)=\Psi({\rm{\bf R}},N)\exp\{-V({\rm{\bf R}})/2D\},
\label{7}
\end{equation}
which dependence on ${\rm{\bf R}}$ is determined by effective
potential

\begin{equation}
V\equiv -\int{\rm{\bf F}}{\rm d}{\rm{\bf R}}.
\label{8}
\end{equation}
Corresponding force ${\rm{\bf F}}$  relates with  the
initial potential $U$ in Eq.(\ref{1}) as follows:

\begin{equation}
U={1\over 4D}{\rm{\bf F}}^2+{1\over 2}{\partial{\rm{\bf
F}}\over\partial{\rm{\bf R}}}.  \label{9}
\end{equation}
As is well known \cite{26},
Eq.(\ref{6}) determines the probability
to realize solution  of corresponding Langevin-type equation

\begin{equation}
\partial{\rm{\bf
R}}/\partial N={\rm{\bf F}}({\rm{\bf R}},N)+{\rm{\bf \zeta}}(N)
\label{10}
\end{equation}
for stochastic variable ${\rm{\bf R}}={\rm{\bf R}}(N)$.
Here  the separated
stochastic force ${\rm{\bf\zeta}}$ is  fixed by the
white-noise conditions

\begin{equation}
\langle\zeta(N)\rangle=0,\qquad
\langle\zeta(N)\zeta(N')\rangle=2D\delta (N-N'), \label{11}
\end{equation}
where the angular brackets denote average with respect
to the distribution (\ref{7}).

In going from the above considered case of directed polymer to the
main object of our interest, the random heteropolymer AB, the
coordinate ${\rm{\bf R}}$ of end-point $N$ turns into a stochastic
Ising variable $\theta(n)$, where $\theta(n)=1$ if $n$-th segment
is of type A and $\theta(n)=-1$ otherwise.
Corresponding to the quenched disorder in fixed
sequence of different type segments, the law
$\theta(n)$ of monomer alternation along the chain
is described by master equation
type of that determines the Glauber dynamics \cite{27}.
Corresponding sequence correlator
$\overline{\sigma(n)\sigma(n')}$ for effective spin
$\sigma(n)\equiv\theta(n)-\overline{\theta(n)}$, being
deviation of microscopic value $\theta(n)$ from average
$\overline{\theta(n)}$, takes the form \cite{12}

\begin{eqnarray}
&&\overline{\sigma(n)\sigma(n')}=C_2\exp(-|n-n'|/l),\nonumber\\
&&C_2\equiv 4f(1-f),\qquad f\equiv(1/2)(1+\overline{\theta(n)}),
\label{12}
\end{eqnarray}
where overbar denotes the averaging over composition
(quenched disorder).
Here two sequence characteristics appear
being correlation length $l$ and fraction $f$ of type-A
monomers.

Stochastic variable $\sigma(n)$ that possesses of the
correlator (\ref{12}) with exponential form,
is governed by the effective motion equation

\begin{equation}
{\rm d}\sigma/{\rm d}n=-\sigma/l+s(n), \label{13} \end{equation}
where the stochastic source $s(n)$ is the white noise:

\begin{equation}
\overline{s(n)}=0,\qquad \overline{s(n)s(n')}=2C_2l^{-1}\delta(n-n').
\label{14}
\end{equation}
Relation between the microscopic value $\sigma(n)$ and stochastic
$\delta$-correlated variable $s(n)$ is given by equation

\begin{equation}
\sigma(n)=\int_0^n e^{-(n-m)/l}s(m){\rm d}m,
\label{15}
\end{equation}
that is the solution of Eq.(\ref{13}).
On the contrary to the colored noise
$\sigma(n)$, the white noise $s(n)$ possesses of the
Gauss distribution function

\begin{equation}
P\{s(n)\}=(4C_2\pi/l)^{-1/2}\exp\left\{-{l\over
4C_2}\int_0^N s^2(n){\rm d}n\right\},
\label{16} \end{equation}
that determines quenched disorder with intensity
$4C_2l^{-1}$.
Respectively, the local averaged  field

\begin{equation} \eta({\rm{\bf r}},n)\equiv (4C_2)^{-1/2}
\overline{\sigma(n)\delta({\rm{\bf r}}-{\rm{\bf r}}(n))}
\label{17} \end{equation}
represents the order parameter (here and below the
monomer volume takes to be equal unity).

Apart from usual terms type of that in Eq.(\ref{10}),
effective motion equation for the field (\ref{17}) must contain  the inhomogeneity
contribution $D\partial^2\eta/\partial{\rm{\bf
r}}^2$ (as the first term in r.h.s. of Eq.(\ref{4})) that takes
Fourier transform $-Dk^2\eta_{{\rm{\bf k}}}$.
As a result,
the motion equation for the Fourier transformation

\begin{equation}
\eta_{{\rm{\bf k}}}(n)=N^{-1/2}\int \eta({\rm{\bf r}},n)e^{-{\rm
i}{\rm{\bf k}}{\rm{\bf r}}}{\rm d}{\rm{\bf r}} \label{18}
\end{equation}
of the field (\ref{17}) takes the Langevin form:

\begin{equation}
\partial \eta_{{\rm{\bf k}}}/\partial
n=-(ak)^2 \eta_{{\rm{\bf k}}}- \partial
{\cal H}/\partial\eta^*_{{\rm{\bf k}}}+\zeta_{\rm{\bf k}}.
\label{19} \end{equation}
Here, as above, the continuum approximation
for effective time $n\gg 1$ is used,
the characteristic distance
$a\equiv D^{1/2}=6^{-1/2}b$ is
determined by the Kuhn segment length $b$, and
determination of effective force $f_{{\rm{\bf k}}}=-{\partial
{\cal H}/\partial\eta^*_{{\rm{\bf k}}}}$ (cf. Eq.(\ref{8})) are taken
into account.
The white noise $\zeta_{{\rm{\bf k}}}=\zeta_{{\rm{\bf
k}}}(n)$ is defined by conditions type of Eqs.(\ref{11}):

\begin{equation}
\langle\zeta_{{\rm{\bf k}}}\rangle=0,\qquad \langle\zeta_{{\rm{\bf
k}}}^*(n)\zeta_{{\rm{\bf k'}}}(n')\rangle=\delta_{{\rm{\bf
k}}{\rm{\bf k'}}}\delta(n-n'), \label{20} \end{equation}
where angular brackets denote average over thermal disorder.

\section{Effective Hamiltonian}

To obtain the effective Hamiltonian
${\cal H}\{\eta\}\equiv\Omega'\{m\}$ \cite{28}
in Eq.(\ref{19}) the thermodynamic
potential $\Omega'\{m\}$, being the average with
respect to both the conformation and the sequence scattering,
should be
determined as a function of the averaged order parameter

\begin{equation}
m({\rm{\bf r}})\equiv\sum_n\langle\eta({\rm{\bf
r}},n)\rangle=(4C_2)^{-1/2}\sum_n\overline{\langle\sigma(n)\delta
({\rm{\bf r}}-{\rm{\bf r}}(n))\rangle}.
\label{21} \end{equation}
With this aim, let us write the
partition function in the form of the functional integral \cite{21}

\begin{eqnarray}
&&\hspace{-2cm} Z=\int \delta m({\rm{\bf r}})\exp\left\{C_2\chi\int
m^2({\rm{\bf r}}){\rm d}{\rm{\bf r}}\right\}\times \nonumber\\
&&\times\overline{\left\langle\delta\left\{\sum_{n}\delta({\rm{\bf
r}}-{\rm{\bf
r}}(n))-1\right\}\delta\left\{\sum_{n}
(4C_2)^{-1/2}\sigma(n)\delta({\rm{\bf
r}}-{\rm{\bf r}}(n))-m({\rm{\bf r}})\right\}\right\rangle}.
\label{22}
\end{eqnarray}
Here $\chi >0$ is the composition Flory parameter, the first
$\delta$-function takes into account incompressibility condition, the
second one reduces to determination (\ref{21}) of order parameter
$m({\rm{\bf r}})$. Further, one follows to represent the
$\delta$-functions as functional Laplace expansions over auxiliary
fields $J_\rho$, $J_m$. Then the average expression in
Eq.(\ref{22}) takes the exponential form with the exponent
$\int(J_{\rho}+J_m m){\rm d}{\rm{\bf
r}}-\Omega\{J_{\rho},J_m\}$, where the last term, being
averaged over sequence and configuration sets, gets conventional free
energy at given fields $J_{\rho}$, $J_m$. The steady-state
magnitudes $\bar J_{\rho}$, $\bar J_m$ of these fields are
determined by conditions $\delta \Omega/\delta J_{\rho}=0$, $\delta
\Omega/ \delta J_m=-m$.  Inserting expression for $\bar
J_{\rho}$, $\bar J_m$ to functional $\Omega\{J_{\rho},J_m\}$
gets the thermodynamic potential defined by equation $Z=\int\delta
m_{{\rm{\bf k}}}\exp(-\Omega'\{m_{{\rm{\bf k}}}\})$, where (see more
detailed Refs.\cite{3,12,29})

\begin{eqnarray} &&\Omega'=\sum_{{\rm{\bf
k}}}\tau_{{\rm{\bf k}}}|m_{{\rm{\bf k}}}|^2+{1\over 2}\sum_{{\rm{\bf
k}}{\rm{\bf k}}'}w_{{\rm{\bf k}}{\rm{\bf k}}'}|m_{{\rm{\bf
k}}}|^2|m_{{\rm{\bf k}}'}|^2+\int v({\rm{\bf r}}){\rm d}{\rm{\bf r}},
\label{23}\\ &&\tau_{{\rm{\bf k}}}\equiv \tau+(ak)^2,\qquad
\tau\equiv l^{-1}-C_2\chi. \nonumber \end{eqnarray}
Here the
inhomogeneity contribution is taken into account in the quadratic term,
averaging procedure with respect to the quenched disorder (\ref{15})
results in appearance of $l^{-1}$-terms,
double ${\rm{\bf k}}$-dependence of the kernel
$w_{{\rm{\bf k}}{\rm{\bf k'}}}=4N^{-1}(la)^{-2}({\rm{\bf
k}}^2+{\rm{\bf k'}}^2)^{-1}$ is stipulated by averaging over the
distribution (\ref{2}).
It is worth mentioning that  expression (\ref{23}) agrees with
Refs.\cite{3,12}, whereas
in Refs.\cite{10,30}, where replica approach was used,
the second summand of Hamiltonian
(\ref{23})  with the opposite sign was obtained.
To correspond the
self-action effects the integrand in last term has the
usual expansion form \cite{3,12}

\begin{eqnarray}
&&v=-(\mu/3!)m^3+(\lambda/4!)m^4;\nonumber\\
&&\mu\equiv 12C_3C_2^{-1/2}l^{-1},\qquad
\lambda\equiv24(1+5C_3^2/C_2)l^{-1},
\label{24}\\
&&C_2\equiv 4f(1-f),\qquad
C_3\equiv|1-2f|.  \nonumber
\end{eqnarray} \\

As mentioned above, this method is based on the Eqs.(\ref{12}) that
allows to express quenched disorder correlator
within the transfer matrix approach.
If one uses the more popular
replica method, then the  field
$J_m$ and the order parameter $m$ must take a replica index $\alpha$
under which in the Hamiltonian (\ref{23}) the summation
from $1$ to $n\to 0$ should be made \cite{6}.
Then the quadratic contribution takes the form

\begin{equation}
{1\over 2}\sum_{{\bf
k}\alpha}A_{\alpha\alpha}({\bf k})
|m_\alpha ({\bf k})|^2+
{1\over 2}\sum_{{{\bf k}\atop\alpha\not=\beta}}
A_{\alpha\beta}({\bf k})
m_\alpha({\bf k})m_\beta(-{\bf k}),
\label{3.36}
\end{equation}\noindent\\
where in the limit $n\to 0$ for coincided replica indices is
$A_{\alpha\alpha}({\bf k})\to 2\tau_{\bf k}$.
As it was found out in spin glasses, the peculiarity of the systems
with quenched disorder is governed by the hierarchy of
the phase space  that is characterized by random overlapping
of the  different replicas \cite{6}.
Therefore overlapping parameter
$A_{\alpha\beta}({\bf k})$
in the second term  (\ref{3.36}) is the stochastic variable
under which the average should be made.
Let us take the corresponding distribution in the simplest
Gauss form

\begin{eqnarray}
{\cal P}\{A_{\alpha\beta} ({\bf k})\}\propto \exp\left\{ -{1\over 8}
\sum_{{\bf k}_1{\bf k}_2 \atop\alpha\ne\beta}  u_{{\bf k}_1{\bf
k}_2}^{-1} A_{\alpha\beta} ({\bf k}_1)A_{\alpha\beta} (-{\bf k}_2)
\right\} \label{3.37}
\end{eqnarray}\noindent\\
where dispersion
$ u_{{\bf k}_1{\bf k}_2}\equiv
\sigma^2(la)^{-2}N^{-1}({\bf k}_1^2+{\bf k}_2^2)^{-1}$
is given by parameter $\sigma$ (see \cite{10}).
Then after averaging of the partition function
$Z=\int Dm_{{\bf k}}\exp(-\Omega'\{m_{{\bf k}}\})$
the second term in Eq.(\ref{3.36}) takes the form

\begin{eqnarray}
-{1\over 2} \sum_{{\bf k}_1{\bf k}_2 \atop \alpha\ne\beta}
u_{{\bf
k}_1{\bf k}_2} m_{\alpha}({\bf k}_1) m_{\beta}(-{\bf
k}_1)m_{\alpha}(-{\bf k}_2) m_{\beta}({\bf k}_2). \label{3.38}
\end{eqnarray}\noindent\\
As a result the thermodynamic potential  (\ref{23})
in the replica form is

\begin{eqnarray}
&&\Omega'=\sum_{{\bf k}\alpha}\tau_{\bf k}|m_\alpha({\bf
k})|^2+{1\over 2}\sum_{{\bf k}_1{\bf k}_2\atop \alpha} w_{{\bf
k}_1{\bf k}_2}
|m_\alpha ({\bf k}_1)|^2|m_\alpha ({\bf
k}_2)|^2+ \nonumber\\
&&\label{3.39}\\
&&+\sum_\alpha\int v(m_\alpha){\rm d}{\bf r}-
{1\over 2} \sum_{{\bf k}_1{\bf k}_2\atop \alpha\ne\beta} u_{{\bf
k}_1{\bf k}_2} m_{\alpha}({\bf k}_1) m_{\beta}(-{\bf
k}_1)m_{\alpha}(-{\bf k}_2) m_{\beta}({\bf k}_2).
\nonumber \end{eqnarray}\\
The Eq.(\ref{3.39}) removes the above mentioned contradiction
in the choice of the effective Hamiltonian form:
the positive contribution of the second term
is governed by the  intra-replica interaction (see \cite{12,3}),
whereas negative contribution in \cite{30, 10} is determined
by replica overlapping.
The
distinction of mentioned terms makes itself evident in the fact that
the former leads to renormalization of the value  $\tau_{\bf k}$,
whereas the latter  corresponds to the memory and non-ergodic effects.

In order to carry out such renormalization (see \cite{28})
using of the mean-field approximation it needs to replace one of the
multiplier $|m_\alpha({\bf k})|^2$ in the second term of
Eq.(\ref{3.39}) by the bare Green function $G_{{\bf k}0}$
that corresponds to  $v=0$, $ u_{{\bf k}_1{\bf k}_2}=0$ and is
determined by equation

\begin{equation} G_{{\rm{\bf k}}0}^{-1}=r+2a^2(k-k_0)^2,
\label{25}
\end{equation}
independent on replica number $\alpha$.
To determine
the parameters $r$, $k_0$ one needs to substitute the
Eq.(\ref{25}) into the corresponding Dyson equation

\begin{equation} G_{{\rm{\bf k}}0}^{-1}=\tau_{{\rm{\bf
k}}}+\sum_{{\rm{\bf k'}}}w_{{\rm{\bf k}}{\rm{\bf k'}}} G_{{\rm{\bf
k'}}0}.  \label{26} \end{equation}
Then after integration over
wave vector ${\rm{\bf k'}}$ one obtains

\begin{equation}
r=\tau+(3/4\pi)l^{-2}(2r)^{-1/2},\qquad
k_0^{-1}=2\pi^{1/2}l(2r)^{1/4}a.
\label{27}
\end{equation}
According to the first equation the positive determined parameter $r$
increases smoothly with growth of bare parameter $\tau$
($r\sim\tau^{-2}$ at $\tau<0$, $|\tau|\gg 1$ and $r\sim\tau$ at
$\tau\gg 1$ (see Fig.1a)).
It means that without   self-action effect, random
heteropolymer is stable with respect to the microphase separation
\cite{8}.
The second equation (\ref{27}) means that space period
$\lambda\equiv 2\pi/k_0$ depends
on renormalized thermodynamic parameter $r$,
as it is inherent in random copolymers \cite{3}.
According to Fig.1b with increasing of
$\tau$  the value
$\lambda$ grows monotonously from $0$ to $\infty$.
In so doing the
greater the correlation length $l$, the faster change of
the period $\lambda$ near the point $\tau=0$.

The  final form of the effective Hamiltonian of
random heteropolymer   follows from the thermodynamic
potential (\ref{3.39}) renormalized under acting of
the fluctuations

\begin{eqnarray}
&&{\cal H}=\sum_{{\bf k}\alpha}\tau_{\bf k}|\eta_\alpha({\bf k})|^2
+ \sum_\alpha\int v(\eta_\alpha){\rm d}{\bf r}-
{1\over 2} \sum_{{\bf k}_1{\bf k}_2\atop \alpha\ne\beta}
u_{{\bf
k}_1{\bf k}_2} \eta_{\alpha}({\bf k}_1) \eta_{\beta}(-{\bf
k}_1)\eta_{\alpha}(-{\bf k}_2) \eta_{\beta}({\bf k}_2),
\label{28}
\\ &&r_{{\rm{\bf
k}}}\equiv r+2a^2(k-k_0)^2, \quad  u_{{\bf k}_1{\bf k}_2}\equiv
\sigma^2(la)^{-2}N^{-1}({\bf k}_1^2+{\bf k}_2^2)^{-1} .\nonumber
\end{eqnarray}

\noindent  Here the kernel
$v(\eta_\alpha)$ is determined by Eqs.(\ref{24}), where $m$ is
replaced by $\eta_\alpha$.

\section{Supersymmetric scheme}

To build a SUSY scheme on a basis of
effective motion equation (\ref{19})
let us introduce the generation functional \cite{21}

\begin{equation}
Z\{\eta_{\rm{\bf
k}}\}=\left\langle\delta\left ({\partial \eta_{\rm{\bf
k}}\over\partial n}+{\delta {\cal H}\over \delta \eta_{\rm{\bf
k}}^*}-\zeta_{\rm{\bf k}}\right )\det\left|{\delta\zeta_{\rm{\bf
k}}\over \delta \eta_{\rm{\bf k}}}\right|\right\rangle,\qquad
{\delta {\cal H}\over \delta \eta_{\rm{\bf k}}^*}\equiv{\partial
{\cal H}\over \partial \eta_{\rm{\bf
k}}^*}+2a^2(k-k_0)^2\eta_{\rm{\bf k}}, \label{33} \end{equation}
averaged over noise $\zeta_{{\rm{\bf k}}}(n)$, where
$\delta$-function considers the motion equation (19), the
determinant is Jacobian of transfer from $\zeta_{{\rm{\bf k}}}$ to
$\eta_{{\rm{\bf k}}}$. Then, the functional Laplace representation is
used for $\delta$-function that introduces a ghost field
$\varphi_{{\rm{\bf k}}}(n)$. To attach exponential form for the
determinant in Eq.(\ref{33}) Grassmann
conjugate fields $\psi_{{\rm{\bf k}}}(n)$, $\bar\psi_{{\rm{\bf
k}}}(n)$ should be used, that obey the conditions \cite{21}

\begin{eqnarray} &&
\{\psi,\psi\}=\{\psi,\bar\psi\}=\{\bar{\psi},\bar{\psi}\}=0,
\nonumber\\
&&\int\delta\psi=\int\delta\bar{\psi}=0,
\qquad\int\bar{\psi}\psi\delta^2\psi=1,\qquad
\delta^2\psi\equiv\delta\psi\delta\bar{\psi},
\label{34} \end{eqnarray}
where figure brackets denote anticommutator. Then, assuming that
the averaging over $\zeta_{\rm{\bf k}}(n)$
in Eq.(\ref{33}) is determined by Gauss distribution
with variance 1 (see Eqs.(\ref{20})), the standard form
of the partition function is obtained:

\begin{eqnarray}
&& Z\{\eta\}=\int
P\{\eta,\varphi;\psi,\bar{\psi}\}\delta\varphi\delta^2\psi,
\nonumber\\ &&
P\{\eta,\varphi;\psi,\bar{\psi}\}=\exp(-S\{\eta,\varphi;
\psi,\bar{\psi}\}),\qquad S=\int_0^N L{\rm d}n,
\label{35}\\
&& L=\int[(\varphi\dot
\eta-\bar{\psi}\dot\psi-\varphi^2/2)+({\cal
H}'\{\eta\}\varphi-\bar{\psi}{\cal H}''\{\eta\}\psi)]{\rm d}{\rm{\bf
r}}.  \nonumber \end{eqnarray}
Here point denotes derivative over
"time" $n$, the prime denotes the functional derivation with respect
to the field (\ref{17}), the last expression is taken in ${\rm{\bf
r}}$-representation.

The last expression in Eqs.(\ref{35}) takes the simplest form
\cite{22}

\begin{equation}
L=\int \Lambda (\Phi){\rm
d}^2\vartheta,\qquad \Lambda\equiv\sum_{\rm{\bf k}}
(\bar{{\cal D}}\Phi_{\rm{\bf
k}}^*)({\cal D}\Phi_{\rm{\bf k}})+{\cal H}\{\Phi_{\rm{\bf
k}}\}, \qquad {\rm d}^2\vartheta\equiv{\rm
d}\vartheta{\rm d}\bar{\vartheta}
\label{36}
\end{equation}
if one introduces SUSY generators

\begin{equation}
{\cal D}\equiv{\partial\over\partial\bar{\vartheta}}-
2\vartheta{\partial\over\partial n},
\qquad \bar{{\cal D}}\equiv{\partial\over\partial
\vartheta}
\label{37}
\end{equation}
and SUSY field

\begin{equation}
\Phi=\eta+\bar{\psi}\vartheta+\bar{\vartheta}\psi+
\bar{\vartheta}\vartheta\varphi,
\label{38}
\end{equation}
where Grassmann coordinates $\vartheta,\bar{\vartheta}$ obey the
same relations (\ref{34}) as for the fields $\psi,\bar{\psi}$.
Here functional ${\cal H}\{\Phi\}$ has the same form as the effective
Hamiltonian (\ref{28}), where order parameter $m_{{\bf k}\alpha}$ is
replaced by superfield $\Phi_{\bf k}(\vartheta)$, Eq.(\ref{38}).
In this case replica index $\alpha$ is removed by Grassmann variable
$\vartheta$ that is the formal reason of the replica and
supersymmetry methods identity.
The advantage of the latter method is in the Grassmann conditions
(\ref{34}) that correspond
to the limit $n\to 0$ in the replica
method.

According to consideration \cite{17,22},
the physical meaning of the components
of the SUSY field (\ref{38}) is as follow:
$\varphi$ is the most probable value of
fluctuations of the field conjugated to the order parameter $\eta$,
and combination $\bar{\psi}\psi$
determines density of sharp interphases.
So, using the
4-component SUSY field (\ref{38}) corresponds to the strong
segregation limit \cite{4,5}.
Further we shall consider more simple case of
the weak segregation limit where $\bar{\psi}\psi\equiv 0$.
Then, the SUSY field (\ref{38}) is reduced to the
2-component form

\begin{equation} \phi=\eta+{\theta}\varphi,
\label{39} \end{equation}
where self-conjugate nilpotent
variable ${\theta}\equiv\bar{\vartheta}\vartheta$ is introduced.
Respectively, Lagrangian (\ref{36}) takes the form

\begin{eqnarray}
&& L=\int\Lambda(\phi){\rm d}{\theta}, \qquad \Lambda\equiv
\sum_{\rm{\bf k}}\phi_{\rm{\bf k}}^*D\phi_{\rm{\bf k}}+
{\cal H}\{\phi_{\rm{\bf k}}\}; \nonumber\\
&& D=-{\partial\over\partial\theta}+\left
(1-2\theta{\partial\over\partial\theta}\right )
{\partial\over\partial n}.  \label{40} \end{eqnarray}
The motion equation
for the nilpotent field (\ref{39})
that corresponds to Lagrangian (\ref{40}) reads

\begin{equation} D\phi_{\rm{\bf k}}=-\delta
{\cal H}/\delta\phi_{\rm{\bf k}}^*. \label{41} \end{equation}
In the component form it leads to the equations for the order
parameter $\eta(n)$ and the amplitude of the more probable
fluctuation $\varphi(n)$ of the conjugate field  (see \cite{22}).

\section{Correlators}

Now let us consider SUSY correlator

\begin{equation}
C_{\rm{\bf k}}(n,{\theta};n',\theta ')\equiv
\left\langle\phi_{\rm{\bf k}}^*(n,\theta)\phi_{\rm{\bf
k}}(n',\theta ')\right\rangle. \label{42}
\end{equation}
Multiplying Eq.(\ref{41}) by value $\phi_{\rm{\bf
k}}^*$ and averaging  within zeroth approximation
($v=u=0$ in Eqs.(\ref{28})) one gets \cite{22}

\begin{equation}
C_{\nu{\rm{\bf k}}}^{(0)}(\theta,\theta ')={1+(r_{\rm{\bf
k}}-{\rm i}\nu)\theta+(r_{\rm{\bf k}}+{\rm i}\nu)\theta
'\over r_{\rm{\bf k}}^2+\nu^2}.  \label{43} \end{equation}
Here
conventional frequency $\nu$ denotes Fourier transformation over
"time" being the monomer number $n$, for example

\begin{equation}
C_\nu=\int_0^N C(n)e^{{\rm i}\nu n}{\rm d}n.\label{44}
\end{equation}
The expression (\ref{43}) has characteristic structure
with respect to combination
of the nilpotent variables $\theta$, $\theta '$, that is
inherent in, obviously, not only zeroth approximation but arbitrary
supercorrelator.
In this connection it is convenient to introduce
basis  supervectors

\begin{equation} A(\theta,\theta ')=\theta, \qquad
B(\theta,\theta ')=\theta ',\qquad T(\theta,\theta ')=1, \label{45}
\end{equation}
that have functional production

\begin{equation}
X(\theta,\theta ')=\int Y(\theta,\theta '')Z(\theta '',\theta){\rm
d}\theta '' \label{46} \end{equation}
for any vectors ${\bf X}$, ${\bf Y}$, ${\bf Z}$.
It is easy to see, that the basis  SUSY vectors (\ref{45})
obey the following multiplication rules: ${\bf A}^2={\bf A}$,
${\bf B}^2={\bf B}$, ${\bf B}{\bf T}={\bf T}$, ${\bf T}{\bf
A}={\bf T}$, other products are zero.
Because the set of
vectors ${\bf A}$, ${\bf B}$, ${\bf T}$ is complete, it
is convenient to expand any SUSY correlator over this basis:

\begin{equation}
{\bf C}=G_-{\bf A}+G_+{\bf B}+S{\bf T}. \label{47}
\end{equation}
Here and below subscripts ${\rm{\bf k}}$, $\nu$ are suppressed for
brevity.
Using Eqs.(\ref{39}), (\ref{42}), one obtains for coefficients
of the expansion (\ref{47})

\begin{equation}
G_-=\langle\eta\varphi^*\rangle,\qquad
G_+=\langle\eta^*\varphi\rangle,\qquad S=\langle|\eta|^2\rangle.
\label{48} \end{equation}
So, $G_{\pm}$ represent advanced and
retarded Green functions and $S$ is the structure factor.
In accordance with Eqs.(\ref{43}), (\ref{45}), (\ref{47}),
these functions
within zeroth approximation take the form

\begin{equation}
G_\pm^{(0)}=(r\pm{\rm i}\nu)^{-1}, \qquad
S^{(0)}=G_+^{(0)}G_-^{(0)}=(r^2+\nu^2)^{-1}. \label{49}
\end{equation}

The Dyson equation for SUSY correlator (\ref{42})
is as follows \cite{19,22}

\begin{equation}
{\bf C}^{-1}=({\bf C}^{(0)})^{-1}-u{\bf C}-{\bf \Sigma}.
\label{50}
\end{equation}
Here ${\bf\Sigma}$ is self-energy function,
$u=\sigma^2 (2N)^{-1}(lak_0)^{-2}$ is the typical value
of the  interreplica overlapping $u_{\bf kk'}$ at ${\bf k}=
{\bf k'}={\bf k}_0$.
By analogy with
the expansion (\ref{47})  the self-energy superfunction
${\bf \Sigma}$ that describes the self-action effects is

\begin{equation} {\bf \Sigma}=\Sigma_-{\bf A}+\Sigma_+{\bf
B}+\Sigma{\bf T}. \label{51} \end{equation}
Then using  Eqs.(\ref{49}),
the Dyson equations for component $G_\pm$, $S$ take
the form

\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}{a}}
\setcounter{equation}{51} \begin{equation}
G_\pm^{-1}+uG_\pm=(r\pm{\rm i}\nu)-\Sigma_\pm, \label{52a}
\end{equation} \def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}{b}}
\setcounter{equation}{51} \begin{equation} S=(1+2\pi
C_2l^{-1}\delta(\nu)+\Sigma)G_+G_-(1-uG_+G_-)^{-1}.  \label{52b}
\end{equation}
\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}}\setcounter{equation}{52}

\noindent Here $\delta$-term is  caused by quenched disorder.

To complete the system (52) it needs to express  the components
of the self-energy superfunction by supercorrelators.
Using
the SUSY perturbation theory with accounting of the cubic and
quartic anharmonicities (\ref{24}) for matrix
element of the self-energy function one  obtains  \cite{21}

\begin{equation}
\Sigma(z,z')=
{\mu^2\over 2!}\left (C(z,z')\right )^2
+{\lambda^2\over 3!}\left (C(z,z')\right )^3.
\label{53} \end{equation}
Here $z\equiv\{{\rm{\bf r}},n,\theta\}$ denotes the set of coordinate
${\rm{\bf r}}$, "time" $n$, and nilpotent coordinate $\theta$.
It is important that the multiplication rule in
Eq.(\ref{53}) is not the same as for the functional product
(\ref{46}).
Here the ordinary product should be used
according to the following multiplication rules \cite{16,22}:
$AT=TA=A$, $BT=TB=B$, $T^2=T$, other products are zero.
Then coefficients of the SUSY expansion (\ref{51})
take the form

\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}{a}} \setcounter{equation}{53}
\begin{equation}
\Sigma_\pm(n)=
\left(\mu^2+{\lambda^2\over 2}S(n)\right)S(n)G_\pm(n), \label{54a}
\end{equation}
\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}{b}} \setcounter{equation}{53}
\begin{equation}
\Sigma(n)={1\over 2}\left(\mu^2+{\lambda^2\over
3}S(n)\right)S^2(n),\label{54b} \end{equation}
\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}}\setcounter{equation}{54}

\noindent where ${\rm{\bf r}}$-representation  for
macroscopically homogeneous system is used.

\section{Memory and non-ergodicity effects}

Following to Edwards and Anderson \cite{31}
let us introduce the composition  memory parameter
$q\equiv\langle\eta(n=N)\eta(n=0)\rangle$,
which value determines the long-range correlation
in the different monomer alternation
along the whole heteropolymer chain.
Moreover, the non-ergodicity parameter
$\Delta\equiv g_0-g$ that is difference between the isothermal
susceptibility $g_0\equiv G_-(\nu=0)$ and thermodynamic value
$g\equiv G_-(\nu\to 0)$ is used.
Then, the main correlators acquire the elongated form:

\begin{equation}
G_\pm(\nu)=\Delta+G_{\pm 0}(\nu),\qquad S(n)=q+S_0(n), \label{55}
\end{equation}
where index $0$ denotes the components corresponding to ergodic
system without memory.
Substitution of Eqs.(\ref{55}) into Eqs.(54) gives

\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}{a}} \setcounter{equation}{55}
\begin{eqnarray}
&&\Sigma_\pm(n)=\left (\mu^2+{\lambda^2\over 2}q\right )q\left
(\Delta+G_{\pm 0}(n)\right )+\Sigma_{\pm 0}(n),
\nonumber\\
&&\Sigma_{\pm 0}(n)\equiv(\mu^2+\lambda^2q)S_0(n)G_{\pm 0}(n)
+{\lambda^2\over 2}S_0^2(n)G_{\pm 0}(n); \label{56a}
\end{eqnarray}
\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}{b}} \setcounter{equation}{55}
\begin{eqnarray}
&&\Sigma(n)={1\over 2}\left (\mu^2+{\lambda^2\over 3}q\right
)q^2+\left (\mu^2+{\lambda^2\over 2}q\right )qS_0(n)+\Sigma_0(n),
\nonumber\\
&&\Sigma_0(n)\equiv{1\over
2}(\mu^2+\lambda^2q)S_0^2(n)+{\lambda^2\over 6}S_0^3(n).
\label{56b}
\end{eqnarray}
\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}} \setcounter{equation}{56}

\noindent Here the terms, being non-linear with respect to
correlators $G_{\pm 0}$, $S_0$, are included
into summands $\Sigma_{\pm
0}$, $\Sigma_0$, the terms
that contain production $S_0\Delta\approx 0$ have been
dropped and, finally, one separates out the first summands that
disappear when parameter $q$ goes to zero.
It is
characteristic that to determine the self-energy functions (56)
the "time"-representation was used, whereas the Dyson equations
(52) require "frequency" Fourier transformations of non-linear
expressions (56). To avoid this difficulty
the fluctuation-dissipation theorem is used \cite{22}

\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}{a}} \setcounter{equation}{56}
\begin{equation}
S_0(n\to 0)=G_{\pm 0}(\nu\to 0)\equiv g,
\label{57a}
\end{equation}
\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}{b}} \setcounter{equation}{56}
\begin{equation}
\Sigma_{\pm 0}(\nu\to 0)=\Sigma_0(n\to 0)\equiv{1\over
2}(\mu^2+\lambda^2q)g^2+{\lambda^2\over 6}g^3,
\label{57b} \end{equation}
\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}}\setcounter{equation}{57}

\noindent
where  Eq.(\ref{56b}) is taken into account in the last equation.

Inserting Fourier transformations of Eqs.(\ref{55}), (56) to
Dyson equation (\ref{52b}), within the $\nu$-representation
one obtains

\begin{equation}
q_0\left[1-ug_0^2-{1\over 2}\left(\mu^2+{\lambda^2\over
3}q_0\right)q_0g_0^2\right]=C_2l^{-1}g_0^2,
\label{58} \end{equation}
\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}{a}} \setcounter{equation}{57}
\begin{equation}
S_0={{(1+\Sigma_0)G_+G_-}\over{1-[u+(\mu^2+\lambda^2q/2)q]G_+G_-}}.
\label{58a}
\end{equation}
\def\arabic{equation}}\setcounter{equation}{58{\arabic{equation}}\setcounter{equation}{58}

\noindent The first of these equations corresponds to $\delta$-term
being caused by memory effects, the second one - to "frequency"
$\nu\not=0$. At $\nu\to 0$ the characteristic production is
$G_+G_-\to g^2$, so that the pole of structure factor (\ref{58a})

\begin{equation}
u+\left(\mu^2+{\lambda^2\over 2}q\right)q=g^{-2}
\label{59} \end{equation}
determines the point of ergodicity breaking.
By analogy,
substituting of Eqs.(\ref{56a}), (\ref{57b}) into the Dyson equation
(\ref{52a})  and taking into account $g\equiv G_-(\nu\to 0)$
for microscopic susceptibility and memory parameter we have

\begin{equation}
1-rg+ug^2+{\mu^2\over 2}g\left[(g+q)^2-q^2\right]+{\lambda^2\over
6}g\left[(g+q)^3-q^3\right]=0.
\label{60} \end{equation}
The thermodynamic behaviour of random heteropolymer in the
vicinity of the ergodicity breaking point is described  by the system
of equations (\ref{58})-(\ref{60}).
By analogy with
spin glass \cite{6}, Eqs.(\ref{58}), (\ref{60}) play role of
Sherrington-Kirkpatrick equations, and Eq.(\ref{59}) determines the
point of de Almeida-Thouless instability.
On further  analyzing
of these equations one should distinguish  macro- and
microscopic values $q_0$, $q$ of the memory parameter and
corresponding susceptibilities $g_0$, $g$.
The peculiarity of such
a hierarchy is that microscopic values, which conform to the limit $\nu
\to 0$, are the usual thermodynamic parameters and depend on
temperature (Flory parameter $\chi$).
The macroscopic values $q_0$, $g_0$ conform to the point $\nu =0$ and
only depend on quenched disorder parameter $l$.
In the non-ergodic area
macroscopic quantities take the values at the point of ergodicity
breaking.

Because the system of three equations (\ref{58})--(\ref{60})
is insufficient for determination of four quantities
$q_0$, $q$, $g_0$, $g$,
it  must be completed by equation

\begin{equation}
q\left[1-ug^2-{1\over 2}\left(\mu^2+{\lambda^2\over
3}q\right)qg^2\right]=C_2l^{-1}g^2,
\label{61}
\end{equation}
that is microscopic analogue of the equation (\ref{58})
obtained as a result of separating out of the singular $\delta$-terms
 for the structure factor in the Dyson equation.
As is known from the spin glass theory \cite{6},
the hierarchy of such singularities, which correspond to
a set of infinitely decreasing "frequencies" $\nu \to 0$,
is a key point of the non-ergodic systems.
In our case, the equations (\ref{58}), (\ref{61})
correspond to the point $\nu=0$ and a minimal of these "frequencies"
respectively.
In turn, at $\nu=0$ the equation (\ref{59}) reads:

\begin{equation}
u+\left(\mu^2+{\lambda^2\over 2}q_0\right)q_0=g_0^{-2}.
\label{62}
\end{equation}
Eqs.(\ref{58}), (\ref{60})--(\ref{62}) are the complete system
for determination of the quantities $q_0$, $q$, $g_0$, $g$.

\section{Discussion}

According to Eqs.(\ref{58}), (\ref{62}) the macroscopic memory
parameter is given by cubic equation

\begin{equation}
(\mu^2/2+\lambda^2q_0/3)q_0^2=C_2l^{-1}.
\label{63} \end{equation}
The characteristic form of the dependence of the value $q_0$
on intensity  of quenched disorder $l$ are depicted in Fig.2.
For copolymers close to symmetric composition $f=0.5$ ($C_3\ll
C_2$) the first term in brackets of Eq.(\ref{63}) is negligible, and
the dependence $q_0\propto l^{1/3}$ takes place. In
the opposite case of
dilute copolymer, where $f\ll 1$ ($C_2\ll C_3$), one obtains the
dependence $q_0\propto f~l^{1/2}$.

Simultaneous solution of Eqs.(\ref{58}), (\ref{60}) and (\ref{62})
gives the point of ergodicity breaking $\chi_0$,
which dependence on correlation length $l$
is depicted in Fig.3 (solid curve).
A characteristic of this phenomenon is that non-zeroth value of
$\chi_0$ appears above a critical value of correlation length $l$,
and with further growth of $l$ dependence reaches its maximum and
then  monotonously falls down.
In so doing, ergodic region is located under curve $\chi_0(l)$
and contracts  with correlation length growth.
Condition ${\rm d} g/{\rm d}\chi =-\infty$ added by Eqs.(\ref{60}),
(\ref{61}) gives the equation

\begin{equation} u+\mu^2(g_f+q)+{\lambda^2\over 2}(g_f+q)^2=
g_f^{-2},
\label{64} \end{equation}
that defines $\chi_f$ value of Flory parameter in the
freezing point, under which the microscopic susceptibility
$g$ takes zeroth value (see Fig.4).
Corresponding dependence  $\chi_f(l)$ on the correlation length $l$
is depicted in Fig.3 (thin curves).
A characteristic is that the dependence $\chi_f(l)$ is
under ergodicity
breaking curve $\chi_0 (l)$ and has the same form.
The influence of the composition on the values of  $\chi_0$
and $\chi_f$ is shown in Fig.3a.
As we recede from the composition $f=0.5$,
the growth of the above mentioned parameters is observed.
The more complicated behaviour is realized at the growth of the
interreplica overlapping  $\sigma$ (see Fig.3b).
At small magnitudes of the correlation length $l$,
 the growth of both $\chi_0$ and $\chi_f$
values is observed with $\sigma$ increasing,
whereas at large magnitudes $l$,
parameters $\chi_0$ and $\chi_f$ decrease.

Dependencies of the macroscopic $g_0$ and microscopic $g$
susceptibilities on the parameter $\chi$ are depicted in Fig.4.
Under point of  ergodicity breaking $\chi_0$, these susceptibilities
and corresponding  memory parameters $q$, $q_0$ coincide.
The dependence $g(\chi)$ has
a cut-off at the point $\chi=\chi_f$
(below it, susceptibility $g$ takes zeroth value
that  corresponds to the freezing state).
Above the ergodicity breaking point $\chi_0$,
the macroscopic susceptibility $g_0$ is constant and
the microscopic one has smooth decrease
(latter can be obtained by simultaneous solution
of Eqs.(\ref{60}), (\ref{61})).
According to Fig.4a, with moving away from the composition $f=0.5$,
the values of the susceptibilities $g_0$, $g$
in the freezing point and in the point of ergodicity breaking
decrease but corresponding values of
parameters $\chi_0$ and $\chi_f$ increase.
The dependence on the correlation length $l$
is shown in Fig.4b: with increasing of $l$ ergodic area contracts
as it should be.
The influence of the interreplica  overlapping is shown in Fig.4c.
With increasing of the corresponding parameter  $\sigma$, the values
of the susceptibilities $g(\chi)$ and $g_0(\chi)$ decrease and hence
the interreplica  overlapping  prevents the
heteropolymer freezing.

The influence of the thermodynamic parameter $\chi$ on the
microscopic memory parameter $q$ is depicted in Fig.5.
The absence of memory below $\chi_f$
is the characteristic feature of the freezing region.
Non-zeroth value
of $q$ appears in the freezing point $\chi_f$ and with further
growth of the parameter $\chi$, the memory parameter monotonously increases.
The step-like behaviour of the parameter $q$ is
inherent in the first order transition.
Evidently, the
physical reason of the mentioned  behaviour is the fluctuation
contribution into the thermodynamic potential of the heteropolymer.
According to Fig.5a, with moving away from the composition $f=0.5$,
the dependence $q(\chi)$ becomes more slight.
On the
contrary, the growing of the  correlation length $l$ results in
more abrupt growth of the memory parameter (see Fig.5b).
At last, one can see from Fig.5c that influence of the
interreplica overlapping  $\sigma$ above and below of the ergodicity
breaking point occurs to be opposite.

According to the Fig.6, the non-ergodicity parameter
$\Delta (\chi)$  grows monotonously with  increasing $\chi$ from
the   ergodicity breaking point $\chi_0$.
The deflection from the composition $f=0.5$, the decrease
of the correlation length $l$ and the growth of the interreplica
overlapping $\sigma$ result in  weakening of the non-ergodicity effects
(see Fig.6).

In the search for a new polymers with predetermined features,
the  phase diagram  plays the basic
role  that defines possible thermodynamic states
at different values of the Flory parameter $\chi$ and
composition $f$.
According to Fig.7, such diagram has a concave form
for both the freezing point $\chi_f(f)$ and the point
of ergodicity breaking $\chi_0(f)$.
A region of the large $\chi$ adjacent
to the composition $f=0.5$ corresponds to the non-ergodic
unfreezing state. With decreasing of the values $\chi$ and $|f-1/2|$,
at first the system goes to the ergodic state  and  then it freezes.
The comparison of the Fig.7a and Fig.7b shows that the increasing of
the correlation length $l$ results in expansion of the unfreezing
and non-ergodic phases.
On the contrary, from Fig.7a and Fig.7c one can see that interreplica
overlapping leads to its contraction.

\newpage


\section{Introduction}
The detailed properties of multiplicity distributions of gluon and quark jets
are studied in experiment nowadays \cite{pep,del,opa}. Theoretical
description of the data in the framework of perturbative QCD is rather
successful qualitatively, often at the accuracy  of 10--15$\%$
or better. Such an accuracy, obtained sometimes only in low-order
approximations, is somewhat surprising because the QCD expansion parameter is
rather large at present energies (about 0.5), so higher order
contributions should be estimated. It happens that different physical
quantities are sensitive in a different way to higher order corrections. A
thorough discussion of average multiplicities of gluon and quark jets is
given in [4-8] with respect to their dependence on these corrections. Here we
extend these studies to the widths and higher moments of multiplicity
distributions of gluon and quark jets, calculating them in the
next-to-next-to-next-to leading order approximation. Such an accuracy becomes
possible because the equations for the generating functions admit a
perturbative expansion. Besides finding out the appropriate coefficients, we
also discuss some general problems, revealed by such an analysis of average
multiplicities, widths, and higher moments.

\section{Widths of multiplicity distributions in perturbative QCD}
Any moment of the parton multiplicity distribution of gluon and quark jets in QCD can be obtained from the equations for the generating functions
\begin{eqnarray}
&G_{G}^{\prime }&= \int_{0}^{1}dxK_{G}^{G}(x)\gamma _{0}^{2}[G_{G}(y+\ln x)G_{G}
(y+\ln (1-x)) - G_{G}(y)] \nonumber \\
&+&n_{f}\int _{0}^{1}dxK_{G}^{F}(x)\gamma _{0}^{2}
[G_{F}(y+\ln x)G_{F}(y+\ln (1-x)) - G_{G}(y)] ,   \label{50}
\end{eqnarray}
\begin{equation}
G_{F}^{\prime } = \int _{0}^{1}dxK_{F}^{G}(x)\gamma _{0}^{2}[G_{G}(y+\ln x)
G_{F}(y+\ln (1-x)) - G_{F}(y)] ,                                   \label{51}
\end{equation}
where  $G_i$ are the generating functions of the multiplicity distributions
$P_n^{(i)}$ of gluon ($i=G$) and quark ($i=F$) jets, defined by
\begin{equation}
G_i(y,z)=\sum _{n=0}^{\infty }(z+1)^nP_n^{(i)}=\sum _{q=0}^{\infty }
\frac {z^q}{q!}\langle n_i\rangle ^qF_q^{(i)}.   \label{gen}
\end{equation}
In these expressions,
$\langle n_i\rangle =\sum _{n=0}^{\infty }nP_n^{(i)}$ is the average
multiplicity, $z$ is an auxiliary variable, $y=\ln (p\Theta/Q_0)$ is the
evolution variable, $p, \Theta $ are the momentum and the opening angle of a
jet, $Q_0$=const ,
 $G^{\prime }(y)=dG/dy ,$ and $ n_f$ is the number of active flavors. Moreover,
\begin{equation}
\gamma _{0}^{2} (y)=\frac {2N_{c}\alpha _S(y)}{\pi } ,  \;\;
\alpha _S(y)=\frac {2\pi }{\beta_0y}\left ( 1-\frac {\beta _1\ln 2y}{\beta _0^2y}\right ),               \label{52}
\end{equation}
\begin{equation}
\beta _0=\frac {11N_c-2n_f}{3}, \;\;  \beta _1=\frac {17N_c^2-n_f(5N_c+3C_F)}{3},
\end{equation}
$\alpha _S$ is the running coupling strength, $N_c$ is the number of colours,
and $C_{F}=  (N_{c}^{2}-1)/2N_{c}=4/3$ in QCD. The argument of $\gamma _0^2$
in the integrals is chosen to be $y+\ln x(1-x)$, as determined
 by the transverse momentum of partons at the splitting vertex. The kernels of the equations are
\begin{equation}
K_{G}^{G}(x) = \frac {1}{x} - (1-x)[2-x(1-x)] ,    \label{53}
\end{equation}
\begin{equation}
K_{G}^{F}(x) = \frac {1}{4N_c}[x^{2}+(1-x)^{2}] ,  \label{54}
\end{equation}
\begin{equation}
K_{F}^{G}(x) = \frac {C_F}{N_c}\left[ \frac {1}{x}-1+\frac {x}{2}\right] .
\label{55}
\end{equation}
The normalized factorial moment of any rank $q$ can be obtained by
differentiation
\begin{equation}
F_q^{(i)}=\frac {1}{\langle n_i\rangle ^q}\frac {d^qG_i}{dz^q}\vert _{z=0}, \label{fq}
\end{equation}
or, equivalently, by using the series (\ref {gen}) and collecting the terms
with equal powers of $z$ on both sides of the equations (\ref{50}), (\ref{51}).

The normalized second factorial moment $F_2$ defines the width of the
multiplicity distribution, and is related to its dispersion
$D^{2}=\langle n^2\rangle -\langle n\rangle ^2$ by the formula
\begin{equation}
D^{2}=(F_{2}-1)\langle n\rangle ^{2}+\langle n\rangle =K_{2}\langle n\rangle ^{2}
+\langle n\rangle ,     \label{disp}
\end{equation}
where $K_2$ is the second cumulant.

The second factorial moments normalized to their own average
multiplicities squared are
\begin{equation}
F_2^G =\frac {\langle n_G(n_G-1)\rangle }{\langle n_G\rangle ^2}, \;\;\;\;
F_2^F =\frac {\langle n_F(n_F-1)\rangle }{\langle n_F\rangle ^2}. \label{fgff}
\end{equation}

Let us write down their perturbative expansions up to $\gamma _0^3$-terms as
\begin{equation}
F^G_2=\frac{4}{3}(1-f_1\gamma_0-f_2\gamma_0^2-f_3\gamma _0^3),  \label{f2g}
\end{equation}
\begin{equation}
F^F_2=(1+\frac{r_0}{3})(1-\phi_1\gamma_0-\phi_2\gamma_0^2-\phi _3\gamma _0^3), \label{f2f}
\end{equation}
where $r_0=N_c/C_F$   determines the asymptotical value of the ratio of
multiplicities in gluon and quark jets (see (\ref{ra}) below). It is equal
to 9/4 in QCD  and to 1 in SUSY QCD. Actually, the asymptotical values of
$F_2^G$ and $F_2^F$ in front of the brackets in (\ref{f2g}), (\ref{f2f}) are
found out from the equations below by equating the leading terms on both
sides. However, we have inserted  their explicit expressions directly here to
simplify further notations.

Using the Taylor series expansion of $G$'s as proposed in \cite{dr}, one can
rewrite the Eqs. (\ref{50}), (\ref{51}) as
\begin{eqnarray}
\frac {1}{\gamma _0^2}[\ln G_G]^{''}=G_G-1-2G_G^{'}{\rm v}_1+G_G^{''}{\rm v}_2+
0.5G_G^{'''}{\rm v}_3+
\nonumber
\\
\left (\frac {G_G^{'2}}{G_G}\right )^{'}{\rm v}_{12}+4B\gamma _0^2[(G_G-1){\rm v}_1-
G_G^{'}{\rm v}_2]+
\nonumber
\\
\frac {n_f}{4N_c} [ \left (\frac {G_F^2}{G_G}\right )^{'}{\rm v}_4+
2\left (\frac {G_FG_F^{'}}{G_G}\right )^{'}{\rm v}_5
+\left (\frac {G_FG_F^{''}}{G_G}\right )^{'}{\rm v}_6+\left (\frac
{G_F^{\prime 2}}{G_G}
\right  )^{'}{\rm v}_{13}-
\nonumber
\\
2B\gamma _0^2\left  ( \left (\frac {G_F^2}{G_G}-1\right ){\rm v}_4+2\left (
\frac {G_FG_F^{'}}{G_G}+
\left ( \frac {G_F^2}{G_G}\right )^{'}\right ){\rm v}_5\right ) ] ,    \label{gg''}
\\
\frac {r_0}{\gamma _0^2}[\ln G_F]^{''}=G_G-1-G_G^{'}{\rm v}_7-G_G^{''}{\rm v}_8-
0.5G_G^{'''}{\rm v}_9-
\nonumber
\\
\left (\frac {G_GG_F^{'}}{G_F}\right )^{'}{\rm v}_{10}-0.5\left
(\frac {G_GG_F^{''}}{G_F}\right )^{'}{\rm v}_{11}-
\left (\frac {G_G^{'}G_F^{'}}{G_F}\right )^{'}{\rm v}_{14}+
\nonumber
\\
2B\gamma _0^2[(G_G-1){\rm v}_7+2G_G^{'}{\rm v}_8+
(G_G^{'}+\frac {G_GG_F^{'}}{G_F}){\rm v}_{10}].      \label{gf''}
\end{eqnarray}
The terms up to the third derivative of $G$ are kept everywhere because
each derivative gives rise to the factors $\gamma $ or $\gamma _F$ as seen
from their definition below in
Eq. (\ref{an}), and we are interested in corrections up to $\gamma _0^3$.
We use $B=\beta _0/8N_c, \; B_1=\beta _1/4\beta _0N_c$.
The corresponding expansions for the anomalous QCD dimensions $\gamma $ and
$\gamma _F$ defined as
\begin{equation}
\langle n_G\rangle \propto \exp (\int ^y\gamma (y')dy'), \;\;
\langle n_F\rangle \propto \exp (\int ^y\gamma _F(y')dy' )     \label{an}
\end{equation}
are
\begin{equation}
\gamma=\gamma _0(1-a_1\gamma _0-a_2\gamma _0^2-a_3\gamma _0^3)   \label{gam}
\end{equation}
with $\gamma _F=\gamma -r'/r$, where
\begin{equation}
r=\langle n_G\rangle /\langle n_F\rangle =r_0(1-r_1\gamma _0-r_2
\gamma _0^2-r_3\gamma _0^3).    \label{ra}
\end{equation}
All the coefficients $a_i, r_i$ have been calculated and tabulated in
\cite{cd}. The integrals ${\rm v}_i$ and terms in the right-hand sides of
equations (\ref{gg''}), (\ref{gf''}) proportional to
$\langle n_G\rangle ^2z^2/2$ are given in the Appendix. The corresponding
terms in the left-hand sides can be written as
\begin{equation}
[\ln G_G]''=\frac {\gamma _0^2\langle n_G\rangle ^2z^2}{2}\frac {4}{3}[1+
\sum _{1}^{3}M_n\gamma _0^n];  \;\; r_0[\ln G_F]''=\frac {\gamma _0^2\langle
n_G\rangle ^2z^2}{2}
\frac{4}{3}[1+\sum _{1}^{3}N_n\gamma _0^n],    \label{mn}
\end{equation}
where
\begin{eqnarray}
M_1=-(2a_1+0.5B+4f_1),
\nonumber
\\
M_2=a_1^2-2a_2+Ba_1-4f_2+4f_1(2a_1+1.5B),
\nonumber
\\
M_3=2a_1a_2-2a_3+0.5B(3a_2-B_1)-4f_3+2f_2(4a_1+5B)
\nonumber
\\
-4f_1(a_1^2-2a_2+2Ba_1+0.75B^2),
\nonumber
\\
N_1=-2a_1-0.5B-(1+\frac{3}{r_0})\phi _1+2r_1,
\nonumber
\\
N_2=a_1^2-2a_2+Ba_1-(1+\frac {3}{r_0})\phi_2 +(1+\frac {3}{r_0})\phi _1
(2a_1+1.5B-2r_1)+
\nonumber
\\
2r_2+r_1(3r_1-4a_1-3B),
\nonumber
\\
N_3=2a_1a_2-2a_3+0.5B(3a_2-B_1)-(1+\frac {3}{r_0})\phi _3 +
(1+\frac {3}{r_0})\phi _2(2a_1+
\nonumber
\\
2.5B-2r_1)-(1+\frac {3}{r_0})\phi _1[a_1^2-2a_2+2Ba_1+0.75B^2+
\nonumber
\\
2r_2+r_1(3r_1-4a_1-5B)]
+2r_3+r_2(6r_1-4a_1-5B)+
\nonumber
\\
r_1(4r_1^2-7.5Br_1-6a_1r_1+4Ba_1+2a_1^2-4a_2+1.5B^2).  \label{mnc}
\end{eqnarray}
The terms with the same power of $\gamma _0$ in expressions on both sides
should be equal. Therefore, one gets
\begin{eqnarray}
f_1=\frac{1}{3}[4{\rm v}_1-2a_1-B/2
-\frac{n_f}{4N_c}{\rm v}_4\left(1-\frac{5}{r_0}+\frac{6}{r_0^2}\right) ],
\\
\phi_1=\frac{1}{1+3/r_0}(f_1+2r_1+2{\rm v}_7-2a_1-B/2),
\\
f_2 = \frac{2}{3}f_1(4a_1-2{\rm v}_1+3B) +
     \frac{1}{3}( a_1^2 - 2a_2 - 4{\rm v}_2 - 4a_1 {\rm v}_1 + B a_1 - 4 B {\rm v}_1 )
\nonumber
\\
 +\frac{n_f}{4N_c}\left\{
 {\rm v}_4 [\phi_1\frac{1}{r_0}(\frac{1}{3}+\frac{1}{r_0}) - \frac{2}{3} f_1
 + 4\frac{r_1}{r_0}(\frac{1}{3}-\frac{1}{r_0})
                + (a_1+B)(\frac{1}{3}-\frac{5}{3r_0} + \frac{2}{r_0^2}) ]
\right.
\nonumber
\\
\left.
- 4\frac{{\rm v}_5}{r_0^2}(1-\frac{r_0}{3})
\right\},
\\
\phi_2 = \phi_1( 2a_1 + \frac{3}{2}B-2r_1 ) +
     \frac{1}{1+3/r_0}(
     f_2 + 2 r_2 + 3 r_1^2 - 2 a_2 - 4 a_1 r_1
\nonumber
\\
+ a_1^2 - 3 B r_1 + B a_1 - 2{\rm v}_7 (f_1+a_1+B) + 4 {\rm v}_8 + \frac{4}{r_0}
{\rm v}_{10}) \ \end{eqnarray}

\begin{equation}
f_3=\frac {1}{3}[d_f+\sum _{i=1}^{6}d_i{\rm v}_i+d_{12}{\rm v}_{12}+d_{13}{\rm v}_{13}],
  \label{f3}
\end{equation}
where $d_f=M_3+4f_3$,
\begin{equation}
\phi _3=\left (1+\frac {3}{r_0}\right )^{-1}[f_3+d_{\phi }+\sum _{i=7}^{11}d_i
{\rm v}_i+d_{14}{\rm v}_{14}],  \label{p3}
\end{equation}
where $d_{\phi }=N_3+(1+\frac {3}{r_0})\phi _3$.  The coefficients $d_i$ are
given in the Appendix.
The numerical values  of $f_i, \, \phi_i$ for different number of active
flavors are shown in the Table 1.

\bigskip

Table 1

\bigskip

\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
$n_f$ & $f_1$ & $f_2$&$f_3$ & $\phi_1$ & $\phi_2$&$\phi _3$\\
\hline
3 & 0.364 & -0.0279 &  0.795 & 0.637 & -0.276 &  2.12      \\
\hline
4 & 0.358 & -0.0457 & 0.740 & 0.631 & -0.286 &  2.04  \\
\hline
5 & 0.352 & -0.0629 & 0.689 & 0.625 & -0.295 &  1.95        \\
\hline
S & 0.313 & 0.310 &  -0.120    & 0.313 & 0.310 & -0.120  \\
\hline
\end{tabular}

\bigskip

Herefrom it is easy to see that the asymptotical $(\gamma _{0}\rightarrow 0)$
values of $F_{2}^{G}$ and $F_{2}^{F}$ are different as has been known since
long ago \cite{and}:
\begin{equation}
F_{2, as}^{G}=\frac {4}{3}, \;\;\;\; F_{2, as}^{F}=1+\frac {r_0}{3}=
\frac {7}{4}.  \label{fas}
\end{equation}
The experimental values \cite{opa} for 41.8 GeV gluon jets, $F_2^G=1.023\pm
0.008\pm 0.011$,
and for 45.6 GeV uds quark jets, $F_2^F=1.0820\pm 0.0006\pm 0.0046$, are much
lower than the above asymptotical limits. If one accepts, however, the
effective value of $\alpha _S$, averaged over all the energies of partons
during the jet evolution, to be 0.2, then one gets
$F_2^G(NLO)\approx 1.039$ and $F_2^F(NLO)\approx 1.068$,
by taking into account only
the first correction proportional to $\gamma _0$. This is quite close to the
experimental results \cite{opa}. For $\alpha _S=0.12$, one gets about 10$\%$
higher values of $F_2$'s. Thus we conclude that the NLO-approximation
describes the widths at $Z^0$-energies within 10--15$\%$ accuracy.

At lower energies the widths should be slightly smaller due to the slow
increase of $\alpha _S$ and somewhat smaller effective values of $n_f$
leading to larger $f_1, \phi _1$. Using the values of $\alpha _S$ given
in the PDG-data \cite{pdg} and $n_f=4$, we predict the energy dependence of
NLO-values of second factorial moments shown by the solid
curves in Fig. 1.

\noindent
The additional dots at the ends of these curves demonstrate the effect due
to possible change of effective values of $n_f$ to 3 at $Q=10$ GeV and to
5 at $Q=90$ Gev, where $Q$ is the total energy as in \cite{pdg}. They
show how small is the indefiniteness imposed by effective number of
active flavours which is the only free parameter in such an approach .
The curves indicate that the widths are closer to
Poissonian ones at lower energies. Qualitatively, it agrees with experimental
trends observed by the DELPHI collaboration \cite{del}. The coupling strength
$\alpha _S$ changes within the interval of $Q$ shown here from 0.18 to 0.12.
The choice of the energy scale for jets is, however, highly nontrivial
(see, e.g., \cite{del,opa}). Therefore we do not plot experimental values
here, just claiming the qualitative agreement within 15$\%$ accuracy.

The cut-off of the integration region at $\varepsilon=e^{-y}\approx e^{-2\pi /
\beta _0\alpha _S}$ from below and at $1-\varepsilon $ from above is not very
crucial at present energies as seen from the dashed curves in Fig. 1. It
diminishes the correction terms and, therefore, slightly increases the
widths and flattens their energy dependence. Thus, the role of the power
corrections is not very important.

Unfortunately, the higher order corrections do not improve our estimates. On
the contrary, the 2NLO-term is positive and tends to violate slightly the
agreement with experiment while 3NLO corrections are negative and so large
that lead even to sub-Poissonian widths of distributions ($F_2<1$) for
$\alpha _S=0.2$. These terms are approximately equal to NLO corrections due
to large values of $f_3$ and $\phi _3$. The inception of such large values
can be traced to rather large contributions of integrals containing $\ln ^2x$,
i.e. to the region of very soft gluons. Thus the cut-off at $e^{-y}$ and
$1-e^{-y}$ becomes more important for these terms.

The 3NLO-corrections are overestimated due to the adopted Taylor series
expansion with the assumption $y\gg \vert \ln x\vert $ which is invalid for
soft gluons. For example, the $k_t$-dependence of the coupling strength is
transformed so that
\begin{equation}
\alpha _S\propto \frac {1}{y+\ln x(1-x)}\approx \frac {1}{y}\left (1-
\frac {\ln x(1-x)}{y}\right ),  \label{yx}
\end{equation}
and the second term becomes infinitely large at $x\rightarrow 0$. The cut-off
at $x=e^{-y}$ leads to a factor of 2 only. Thus the above expansion implies
some special presumption about the coupling strength behavior in the
nonperturbative region as well as its modification at the limits of the
perturbative one.
 The series (\ref{f2g}), (\ref{f2f}) are with the sign-changing and increasing
(in modulus) terms. From the values of $f_i, \phi _i$ in Table 1,
one concludes that higher order
terms are more important for the width of a quark jet compared with
a gluon jet.

The slopes of the widths are especially sensitive to these higher order terms
because each of them is enlarged by the factor  $n$ when differentiating
$\gamma _0^n$. Thus 3NLO contribution is about 3 times larger than the NLO
term in the slopes of widths. It demonstrates that any precise quantitative
estimates of slopes become impossible. In particular, at present energies one
cannot trust NLO estimates of these slopes as being small: $F_2^{G'}(NLO)
\approx 0.04; \; F_2^{F'}(NLO)\approx 0.092$ at $\alpha _S=0.2$. However, one
can predict the asymptotical value of the ratio of slopes as
\begin{equation}
\frac {(F_{2}^{G})^{'}_{as}}{(F_{2}^{F})^{'}_{as}}=\frac {16f_1}{21\phi _1}
\approx 0.43,   \label{rfpr}
\end{equation}
which surely coincides with their NLO ratio. It demonstrates that the second
factorial moment of quark jets approaches its asymptotical value faster than
for gluon jets.

Let us stress that all slopes and curvatures in pQCD are related to the running
property  of the QCD coupling constant since they are proportional to its
derivatives which are equal to zero for a fixed coupling constant. It is
interesting to note that the two-loop term in $\alpha _S$ proportional to
$\beta _1$ contributes only to the left-hand side of Eqs (\ref{mn}), namely,
to the coefficients $M_3$ and $N_3$ in (\ref{mnc}), and its role is very mild
there (about 1--2 $\%$). Thus it can be accounted with high accuracy
considering it only in the expressions for the expansion parameter $\gamma _0$.

One can also easily check that all the relations of SUSY QCD (where $n_f =N_c =
C_F$) are valid for all the coefficients shown above (e.g.,
$F_{2}^{G}=F_{2}^{F}$ etc.).
The SUSY values of $f_i=\phi _i$ are  shown in the lower line of the Table 1
marked by S. The asymptotical SUSY values of $F_2$ are equal to 4/3.

\section{Third moments of the multiplicity distributions}

The system of equations (\ref{gg''}), (\ref{gf''}) valid up to 3NLO-terms
of $\gamma _0^3$ was also applied by us to calculation of third moments of
the multiplicity distributions. It was done by equating the terms proportional
to $z^3$ on both sides of the equations and writing down the third factorial
moments (defined by Eq. (\ref{fq}) at $q=3$) as
\begin{equation}
F_3^G=h_0(1-\sum _{i=1}^{3}h_i\gamma _0^i); \;\;
F_3^F=g_0(1-\sum _{i=1}^{3}g_i\gamma _0^i).  \label{f32}
\end{equation}
Proceeding in the same way as done above, one gets $h_0=9/4, \; g_0=1+r_0+
r_0^2/4$ and the values of the coefficients $h_i, g_i$ in (\ref{f32}) shown
in the Table 2. The asymptotic limit of the third moment of quark jets is
about twice larger than that of gluon jets.
The analytic expressions are too lengthy to be presented here.

\bigskip

Table 2

\bigskip

\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
$n_f$ & $h_1$ & $h_2$ & $h_3$ & $g_1$ & $g_2$ & $g_3$ \\
\hline
3     & 0.986  & -0.342 & 2.49 & 1.61 & -1.58 & 7.74   \\
\hline
4     & 0.972 & -0.380 & 2.36 & 1.60 & -1.59 & 7.54 \\
\hline
5     & 0.957 & -0.417 & 2.25 & 1.59 & -1.60 & 7.34 \\
\hline
S     & 0.844 & 0.722 & -1.09 & 0.844 & 0.722 & -1.09 \\
\hline
\end{tabular}

\bigskip

Comparing $f_i, \phi _i$ with $h_i, g_i$ one concludes that the corrections
increase for higher moments even in the NLO-approximation. Moreover, at present
values of $\gamma _0\approx 0.5$ they are rather large.
Let us note the similarity in the structure of corrections for widths and third
moments. They alternate in sign, and third coefficients are larger than the
first ones. It is an indication on the sign-alternating asymptotic series,
and Borel summation can be effective here. The increase of the coefficients
originates from the terms containing the integrals of the type $\int _0^1
\ln ^nxdx \propto n!$. The termination of the cascade at $\varepsilon =e^{-y}$
leading to power corrections becomes more important below $Z^0$. This reminds of
the situation with renormalons (see, e.g., \cite{zak}).

In SUSY QCD the asymptotical values of third moments are equal to 9/4. The
first correction given by $h_1$(SUSY)=$g_1$(SUSY)=0.844 is almost as large as
for ordinary gluon jets. It is similar to the correction for second moments.
However, NNLO and 3NLO-terms for moments in SUSY QCD differ drastically from
those for ordinary jets both in absolute values and signs as seen from the
Tables 1 and 2. It demonstrates their sensitivity to the value of $r_0$
which is drastically different in the two cases.

The similar procedure can be used for higher rank moments as well.

\section{General discussion and conclusions}
The equations (1), (2) are dealing with probabilities and, therefore, are of
classical nature. However, the quantum-mechanical interference has been
accounted in the angular ordering effect. Nevertheless, there is no proof of
their validity at all orders. The approach advocated above treats
these equations as the kinetic equations in QCD for partonic
cascade processes. Implicitly we have assumed that these equations describe
the cascade down to extremely low energies of partons by imposing the limits
0 and 1 of integration in the shares of energy. Thus the non-perturbative
region of soft partons was assumed to be described in the same manner as the
perturbative cascade. Probably, one should include into the consideration
only the perturbative region by the requirement that the evolution parameters
under the integral in (1), (2) stay always positive i.e. by introducing the
cut-off at $\varepsilon $
and $1-\varepsilon $ where $\varepsilon =e^{-y} $. Such a modification would
lead to the power-like corrections which can be neglected asymptotically but
contribute at present energies. Their role was considered for average
multiplicities in \cite{cd}. They are not very important for the moments in
the NLO-approximation as shown above.
However, more thorough treatment is needed,
especially, in view of rather large contribution of soft gluons to $f_3$ and
$\phi _3$.

Other  analytic solutions of these equations  different from the perturbative
one as well as nonperturbative modifications of the equations
\cite{imd,los,ekg} can be looked for.
Especially interesting would be to learn more about
the singularities of the generating functions in the $z$-plane which are
known up to now for the leading order solution only.

Leaving this program for future studies, we can now compare the results
obtained from perturbative solutions of the equations for average
multiplicities  [3-7] and for widths of the multiplicity distributions of
gluon and quark jets up to 3NLO approximation of pQCD. Leading order
predictions for any quantity are quite far from present experimental data.
NLO corrections are always pointing in the right direction of closer
agreement with experiment.  In particular, the energy dependence of the gluon
jet average multiplicity, of the ratio of its slopes for gluon and quark jets
and the values of widths at $Z^0 $-energies can be fitted with rather high
accuracy. However, it still fails in the ratio of average multiplicities of
gluon and quark jets differing from experiment by about 30$\%$ in absolute
values and even more in its energy dependence slope \cite{cd}. Here,  2NLO
and 3NLO terms improve the situation so that the ratio $r$ differs by about
10--15$\%$ only. They do not spoil good qualitative features of NLO in the
energy dependence \cite{dg}.
However, these corrections become very large for the widths and for the third
moments as shown above, as well as for the slope and curvature in the energy
dependence of average multiplicities \cite{cd}.

Thus we conclude that present perturbative QCD results can describe the
experimental data within 10--15$\%$ accuracy. The perturbative series breaks
down, however, at different orders for different quantities.
There seems to be no standard way to improve the analytic results and
unanimously predict where one should truncate the expansion.
Nevertheless, the general trends obtained from the
perturbative approach are steadily indicating qualitative convergence of
theory and experiment. Moreover, the computer results \cite{lo} show that
the exact solutions of QCD equations can be even closer to experiment.

\acknowledgements
This research is supported in part by the Natural Sciences and
Engineering Research Council of Canada, and the Fonds pour la Formation de
Chercheurs et l'Aide \`a la Recherche of Qu\'ebec, and also by INTAS
and Russian Foundation for Basic Research.


\section{\@startsection {section}{1}{\z@}%
                                   {-4.3ex \@plus -.2ex \@minus -.2ex}%
                                   {2.1ex \@plus.2ex}%
                                   {\normalfont\underline}}
\renewcommand\subsection{\@startsection{subsection}{2}{\z@}%
                                     {-2.1ex\@plus .1ex \@minus -.1ex}%
                                     {2.1ex \@plus .2ex \@minus -.1ex}%
                                     {\normalfont\underline}}
 \makeatother
 \renewcommand\figurename{FIGURE}
 \renewcommand\tablename{TABLE}
 \setlength{\parskip}{0pt}
 \setlength{\belowcaptionskip}{-5mm}
 \setlength\captionindent{\parindent}%

\renewcommand{\thefootnote}{\fnsymbol{footnote}}

\begin{document}
 \pagestyle{fancy}
 \setlength{\headrulewidth}{0pt}
 \lhead{}
 \rhead{\begin{picture}(0,0)\put(120,70){\it\thepage}\end{picture}}
 \chead{}
 \cfoot{}
 \hphantom{.}\vskip10mm
\begin{center}
{Linear and non-linear dielectric spectroscopy\\ on ammonium doped Rochelle salt }
\end{center}
\vskip10.5mm
\begin{center}
 U. Schneider, P. Lunkenheimer, J. Hemberger and A. Loidl\\
 Experimentalphysik V, Universit\"{a}t Augsburg, D-86135 Augsburg, Germany
\end{center}
\vskip10.5mm

 \noindent
\underline{Abstract} In this work we present a thorough investigation of the mixed system
of Rochelle salt doped with ammonium, with $0\leq x\leq 1$ for
[Na$($K$_{1-x}($NH$_4)_x)$C$_4$H$_4$O$_6\cdot 4$H$_2$O]. This interesting system is known
to exhibit a rich phase diagram with a variety of polar and non-polar phases. However,
most of the earlier investigations date back many decades and the experimental advances
made since then should enable new insight into the physical nature of its phase diagram.
We studied single crystals of 18 doping concentrations using linear and {\it
non}\,-linear frequency-dependent dielectric response spectroscopy. In the low field
measurements a variety of relaxation features are detected and analyzed in detail. The
results for two concentrations of $x=0$ and $x=0.15$ are described qualitatively with a
simple model. Possibly the most important outcome of the present work is the detection of
a new electrically ordered phase for concentrations $0.18<x<0.89$.
 \vskip 10mm

\noindent \underline{Keywords}\quad ferroelectric material; Rochelle salt; dielectric
spectroscopy; hysteresis
 \vskip 14mm

\noindent \underline{PACS Numbers}\quad 77.22.Ej, 77.22.Gm, 77.80.Bh
 \vskip 14mm

\section{INTRODUCTION}
Rochelle salt (NaK\,C$_4$H$_4$O$_6\cdot{}4$H$_2$O) was first synthesized in 1655 by the
pharmacist {\it Pierre Seignette} in La Rochelle (France). Therefore the name {\it
Seignette salt} is also in use in the literature. Although it enjoyed a rapid
distribution for medical and chemical purposes \cite{Mac77} from the very beginning, its
unusual physical properties, especially the large piezoelectric coefficient, were
discovered in the late 18$^{\rm th}$-century only \cite{Poc06}. Its extraordinary
ferroelectric properties along the crystallographic {\it a}\,-axis were revealed even
later by Valasek \cite{Val21} who transferred the terminology of magnetic phase
transitions ({\it Curie}\,-temperature) into the physics of dielectric materials. In his
work he found two second-order phase transition temperatures at 297\,K and 255\,K which
indicate a narrow temperature range where ferroelectric (FE) order exists.

As reported in literature \cite{kourtschatov36,takagi58,makita58,jona62,sandy67}, on
doping the pure compound, NaK-tartrate\footnote{The name {\it tartrate} is derived from
french {\it cr\`eme de tartre} for the crystalline furring deposit in containers of
wine.}, with ammonium, the temperature range of the FE phase is reduced till the order is
fully suppressed at approximately $x=0.025$. Mixed crystals up to $x \approx 0.18$ remain
paraelectric (PE) down to lowest temperatures. For ammonium concentrations $0.18<x<0.89$
there is a FE phase below a second order transition which was reported to prevail down to
lowest temperatures \cite{jona62}. Between $x=0.89$ and the pure NaNH$_4$-tartrate
compound this FE phase is no longer present and a new phase appears with a first order
phase transition. This phase shows a polar order that cannot be reversed by electrical
fields $E \leq 20$\,kV/cm but by applying additional mechanical stress.

The crystal structure of the non-polar phases of pure NaK-tartrate was determined to be
orthorhombic ($P2_12_12$) \cite{beevers41}. The structure of the polar phase was found to
be monoclinic ($P2_111$) \cite{frazer54}. Despite numerous theoretical approaches there
is no overall consensus concerning the explanation of the strange re-entrant behaviour.
To enhance our understanding of this puzzling feature we were especially interested in
investigating the dipolar reorientation in the low temperature PE phase below the second
phase transition of the pure NaK-tartrate and the doping range of ammonium with
$0.025<x<0.18$ between the two FE phases. We applied both ordinary broadband dielectric
spectroscopy to examine the frequency-dependent linear response and high-field
spectroscopy to detect the non-linear hysteretic response. Various samples of 18
different doping concentrations, were investigated covering all characteristic regions of
the phase diagram.

Finally we will show that the simple model of Blinc {\it et al.} \cite{blinc72} can be
used to describe the doping range of $0<x<0.18$ qualitatively.

 \chead[U. Schneider, P. Lunkenheimer, J. Hemberger and A. Loidl]
 {Linear and non-linear dielectric spectroscopy...}

\section{EXPERIMENTAL DETAILS}

The crystals were grown from an aqueous solution. The crystallographic axes were
identified according to the geometrical properties reported in literature ({\it e.g.}
\cite{jona62,mason50}). The average size of the raw crystals was about 1\,cm$^3$. They
were polished down to a thickness along the {\it a}\,-axis $\approx{}1$\,mm with an area
in the {\it b-c}\,-plane of $\approx 100-150$\,mm$^2$. For applying dielectric
spectroscopy to these samples, both faces were covered with silver paint.

We used commercially available closed-cycle $^4$He refrigerators for cooling. The
frequency-dependent linear dielectric measurements were performed with an autobalance
bridge (HP4284) using pseudo 4-point technique. The frequency range covered is
$20$\,Hz\,$\leq\nu\leq$\,$1$\,MHz. The high-field non-linear measurements were made with
a home-made setup \cite{hem94} designed according to Sawyer and Tower \cite{sawyer30}.

\section{RESULTS AND DISCUSSION}

Frequency-dependent linear measurements were performed in the whole ammonium doping range
of $0\leq x\leq 1$ for the filed along the crystalline $a$ axis. The various PE and FE
phases showing up in the phase diagram are discussed in detail below. In the following,
the presentation of the spectral results will be restricted to certain concentrations,
which exhibit dielectric properties, typical for the respective phases. In
Figure~\ref{fig:lin000}
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig1.eps}
\caption{Temperature dependence of real (a) and imaginary part (b) of the dielectric
permittivity of NaK-tartrate for three frequencies. The inset shows a {\it Curie}\,-plot
of $\varepsilon$'$(T)$ for $\nu = 20$\,Hz. The lines are drawn to guide the eye.}
\label{fig:lin000}
\end{center}
\end{figure}
we present the temperature dependent data of both, the real and the imaginary part of the
dielectric function, $\varepsilon^*=\varepsilon$'$-i\varepsilon$'', of pure NaK-tartrate
for three frequencie s. Clearly the two second-order phase transitions at
$T_{c1}=296.6$\,K and $T_{c2}=254.6$\,K appear as sharp anomalies in $\varepsilon$'$(T)$
(Figure~\ref{fig:lin000}a). In the temperature range of the ferroelectric phase for both,
$\varepsilon$' and $\varepsilon$'' (Figure~\ref{fig:lin000}b) there are indications of
relaxational behaviour which, however, cannot be analyzed further due to the underlying
strong temperature dependence in the vicinity of the phase transitions. Moreover, for
$\varepsilon$'' one can observe unusually high losses due to domain wall movement. The
inset in Figure~\ref{fig:lin000}a shows the inverse of $\varepsilon$' for $\nu = 20$\,Hz.
This plot reveals Curie-Wei{\ss} behaviour at $T<T_{c1}$ and $T>T_{c2}$ while at
$T_{c1}<T<T_{c2}$ its detection is hampered by the relatively close vicinity of both
phase-transition temperatures and the relaxation behaviour mentioned above. The
re-entrant phase transitions, PE$\rightarrow$FE$\rightarrow$PE, prevail up to an ammonium
concentration of $0.025$.

As an example for the doping range of $0.025<x<0.18$ we present the data for $x=0.15$ in
Figure~\ref{fig:lin015}.
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig2.eps}
\caption{Temperature dependence of the real part of the dielectric function of
Na$($K$_{0.85}($NH$_4)_{0.15})$-tartrate for three frequencies. The inset shows the
imaginary part below 100\,K. The lines are drawn to guide the eye.} \label{fig:lin015}
\end{center}
\end{figure}
In this range there is no phase transition. The FE correlations, which produce the
ordered phase for doping levels below $2.5$\,\%, are reduced below a critical limit (see
section~\ref{sec:theory}). In Figure~\ref{fig:lin015} the real part $\varepsilon$'$(T)$
shows a relatively flat, broad peak at about 217\,K, but no FE ordering is detected in
the nonlinear results (not shown), which reveal simple PE loss-loops. At low
temperatures, between 10 and 90\,K there is a pronounced relaxation feature. It shows up
as sigmoidally shaped curves in $\varepsilon$'$(T)$ and as peaks in $\varepsilon$''$(T)$
(inset of Figure~\ref{fig:lin015}), both shifting to higher temperatures with increasing
frequency. The temperature evolution of the corresponding relaxation times will be
discussed below.

The doping range $0.18<x<0.89$ was reported to exhibit a ferroelectrically ordered phase
\cite{kourtschatov36,takagi58,makita58,jona62,sandy67}. All samples investigated in this
range show qualitatively the same characteristics.
Figure~\ref{fig:lin030}
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig3.eps}
\caption{Temperature dependence of real (a) and imaginary part (b) of the dielectric
permittivity of Na$($K$_{0.70}($NH$_4)_{0.30})$-tartrate for three frequencies. The inset
shows a {\it Curie}\,-plot of $\varepsilon$'$(T)$ for $\nu = 1.16$\,kHz. The numbers (I),
(II), and (III) denote the three observed dominant relaxation r\'egimes. The lines are
drawn to guide the eye.} \label{fig:lin030}
\end{center}
\end{figure}
shows the data of a typical crystal (30\,\% doping) within this range. Like the pure
NaK-tartrate this sample shows a sharp peak in $\varepsilon$'$(T)$
(Figure~\ref{fig:lin030}a) due to a second-order transition at 205\,K. On closer
inspection of the data a shoulder on the low temperature side of this peak can be
identified at 187\,K, which indicates a second phase transition. This feature is even
more pronounced in the inverse representation shown in the inset. This additional
transition is detected in all data on crystals in this doping range. The imaginary part
in Figure~\ref{fig:lin030}b shows at least three relaxation processes, namely in the
temperature ranges of $10-90$\,K (I), $90-180$\,K (II) and $180-200$\,K (III)
superimposed to the transition anomalies. The three relaxations can be clearly marked by
their dispersive behaviour which will be discussed in detail below.

For $x>0.89$ a new phase transition appears. The data on pure NaNH$_4$-tartrate are
presented in Figure~\ref{fig:lin100}.
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig4.eps} \caption{Temperature dependence of
the real part of the dielectric function of NaNH$_4$-tartrate for three frequencies. The
lines are drawn to guide the eye. The inset shows a magnified view of the first order
phase-transition region.} \label{fig:lin100}
\end{center}
\end{figure}
The transition at 112\,K from the PE to an ordered phase is of first order as can be
deduced from the hysteretic step in $\varepsilon$'$(T)$ (inset of
Figure~\ref{fig:lin100}).

We will now focus on the variety of relaxation processes mentioned above. At the
temperature of the peak maximum in $\varepsilon$''$(T)$ the relaxation rate $\nu_p$ can
be estimated to be equal to the measuring frequency. The temperature dependent $\nu_p$,
determined in this way for the collection of evaluated processes is shown in
Figure~\ref{fig:relax1}
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig5.eps}
\caption{Collection of relaxation rates in Arrhenius representation of the processes,
which are observed in the whole doping range; part (b) shows an enlarged view of the high
temperature data of (a).} \label{fig:relax1}
\end{center}
\end{figure}
in an {\it Arrhenius} representation. One may classify three groups of separate processes
in the inverse temperature ranges of $17-30$ (I), $6-12$ (II) and $4-6$ (III) in units of
1000/K (with the exception of one relaxation in the $x=0.2$ doped samples). Within the
error bars the various data sets can be fitted\footnote{The fits are not shown to improve
readability of the plots.} by straight lines (nearly parallel within one group),
indicating thermally activated behaviour according to
 \begin{equation}
\nu_{p}(T)=\nu_0 \exp\left(-\frac{E_b}{k_BT}\right)\;  .\label{eq:Arrhenius}
\end{equation}
Here $\nu_0$ is the attempt frequency and $E_b$ denotes the energy barrier. The fits
result in similar parameters for each group. In
Figure~\ref{fig:relax2}
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig6.eps}
\caption{Parameters $E_b$ (a) and $\nu_0$ (b) of the {\it Arrhenius}\,-fits of $\nu_p(T)$
(cf. Equation~\protect\ref{eq:Arrhenius}) versus doping concentration $x$. The lines
suggest interpolated doping dependences of the plotted parameters for the three ranges
(I), (II) and (III) classified according to the corresponding temperature ranges.}
\label{fig:relax2}
\end{center}
\end{figure}
we plot the values of these parameters versus the doping concentration $x$. The lines are
drawn to demonstrate the continuous development of the parameters within the three groups
(I), (II) and (III) up to $x\approx 0.6$. The lines, connecting the Arrhenius parameters
of the three concentrations above $x=0.8$, also belong to the group (III) in terms of the
absolute values of the relaxation times, while their Arrhenius parameters deviate
considerably from those of the other members of this group. The involved errors for the
values shown in Figure~\ref{fig:relax2} are about 25\,\% resulting from the error bars
given in Figure~\ref{fig:relax1}. Within these limits all parameters follow the behaviour
sketched by the lines in Figure~\ref{fig:relax2} reasonably well. Presumably the observed
relaxations result from intra-molecular motions, e.g. of the carboxyle groups or the
hydrogen bonds. This finding is (for group I) also corroborated by the relatively small
doping dependence of the plotted parameters. From the extremely high values of the
attempt frequencies of group (III), deviations from Arrhenius behaviour can be expected
at high temperatures, indicating a collective nature of the probed motion. Assuming an
increase of disorder upon doping the width of the associated random bonds distribution is
proportional to a characteristic "freezing" temperature $T_f$ which thus shifts to higher
temperatures for higher doping concentrations. The Arrhenius fit therefore yields a
higher value for $E_b$ in the same temperature interval and an increased overestimation
of $\nu_0$. For all relaxation processes observed, the inspection of
$\varepsilon$''$(\nu)$ (not shown) reveals a half width of the loss peaks much broader
than expected for Debye relaxation. Such behaviour is usually ascribed to a distribution
of relaxation times.

In the following we will discuss the high field results of two samples with the doping
concentrations $x=0.2$ and $x=1$. In
Figure~\ref{fig:nonlinear020}
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig7.eps}
\caption{Nonlinear response of Na$($K$_{0.80}($NH$_4)_{0.20})$-tartrate for various
temperatures.} \label{fig:nonlinear020}
\end{center}
\end{figure}
the field dependent polarization for the 20\,\% doped sample is presented for various
temperatures. At temperatures $T>200$\,K linear behaviour is observed. At 200\,K a
hysteresis loop starts to open up, typical for a transition to FE order. Below 160\,K a
frequency independent wasp waisted contraction of the curves near the origin is observed.
This effect can be attributed to an overlay of the present FE order with {\it
anti}\,ferroelectric (AFE) correlations \cite{grindlay70,Haertling87,Xu91}. Another
possible reason for this phenomenon found also in magnetically ordering materials might
be a reorientation of domains about 90$^\circ$ as proposed by Weber {\it et al.} for
cobalt films on copper \cite{Weber96}. However, the shoulder, observed in
$\varepsilon$'$(T)$ as discussed above (Figure~\ref{fig:lin030}a) speaks in favour of the
first interpretation. It clearly indicates a second-order transition. Below about 140\,K
the polarization becomes significantly reduced, which can be ascribed to growing
coercitive fields with decreasing temperature leading to an incomplete reorientation of
the polarization.

Figure~\ref{fig:nonlinear100}
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig8.eps}
\caption{Nonlinear response of NaNH$_4$-tartrate for various temperatures.}
\label{fig:nonlinear100}
\end{center}
\end{figure}
shows the high field response of pure NaNH$_4$-tartrate: even with the highest applied
fields (just below the threshold limit) no electrical ordering can be introduced. An
abrupt change to ellipsoidally shaped $P(E)$ loops can be observed around the temperature
of the first-order phase transition at 112\,K indicating an enhanced dielectric loss in
this region.

Finally we will describe the modified phase diagram for the full ammonium concentration
range of $0<x<1$ in Figure~\ref{fig:phases}.
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig9.eps}
\caption{Phase diagram of Na$($K$_{1-x}($NH$_4)_x)$C$_4$H$_4$O$_6\cdot 4$H$_2$O. The
filled symbols indicate transition temperatures determined in the present work except for
the range of $0.025<x<0.18$ where the dotted line follows the temperature of the flat,
broad maximum observed for these samples. The open symbols are taken from the
literature.} \label{fig:phases}
\end{center}
\end{figure}
All previously published phase transitions are reproduced and refined on the basis of our
results. With decreasing temperature we observe a succession of
PE$\rightarrow$FE$\rightarrow$PE phases for $0<x<0.025$, denoted as PE I, FE I and PE II.
Above $x=0.025$ there is no indication for any ordering: the FE I phase vanishes
successively on increasing doping concentration; the dotted line in the phase diagram
indicates the position of a flat, broad maximum in $\varepsilon$'$(T)$. A possible
microscopic explanation for this feature is discussed in the following section. In the
range of $0.18<x<0.89$ however, our data justify a modification of the phase diagram:
while there is indeed a ferroelectric phase transition as reported before, this phase (FE
II) is followed by another r{\'e}gime with {\it anti}\,-ferroelectric correlations at
lower temperatures (AFE). We therefore introduce a new phase line in
Figure~\ref{fig:phases}. This result is supported by the fact that this behaviour was
reproduced in different samples and in measurements using two different experimental
setups. Above $x=0.89$ we observe a "polar" phase (PE III) where no order can be induced
electrically by fields below the threshold limit. The detailed physical nature of this
phase is still unresolved.

\section{MODEL CALCULATIONS}\label{sec:theory}

There has been a tremendous effort
\cite{blinc72,mitsui58,takagi72,blinc74,takagi75,sannikov79,kozlov88,okada89a,okada89b}
to explain the reason of the puzzling behaviour of the pure NaK-tartrate compound, namely
showing two second order transitions and a ferroelectric phase between them, taking into
account dipole-dipole and dipole-charge interaction between sublattices. From vibrational
spectroscopy of soft-mode behaviour \cite{kamba95,volkov99} there is strong evidence in
favour of the model of coupled oscillators of Reference~\cite{kozlov88,volkov99}. In the
present paper we apply the simple model of Blinc and \v{Z}ek\v{s} \cite{blinc72,blinc74}
which considers the interaction of the hydroxyl groups. The model, which is sometimes
called {\it pseudo-spin model}, is based on an asymmetric double-well potential with two
intertwined sublattices being mirror images of each other. Omitting a tunneling parameter
for the present case, because of the high temperatures at which the transition takes
place, the Hamiltonian of this arrangement may be expressed as follows:
\begin{align*}
 {\cal H}= &\frac{1}{2} \sum_{ij}\left[K_{ij}(S_{i1}^zS_{j1}^z+S_{i2}^zS_{j2}^z) + L_{ij}S_{i1}^zS_{j2}^z\right] \nonumber \\
           & - \Delta\sum_j(S_{j1}^z-S_{j2}^z)-2\mu\sum_j\left[E_{j1}(t)S_{j1}^z + E_{j2}(t)S_{j2}^z \right]
\end{align*}
where the $S_{ij}$ are {\it pseudo}\,spin-$\frac{1}{2}$ operators with the indices 1 and
2 denoting the two sublattices; the tensors $K_{ij}$ and $L_{ij}$ stand for the effective
interaction constants of the dipoles in the same or in the other sublattice. $\Delta$ is
a measure for the symmetry of the double wells and $\mu$ is the dipole moment interacting
with an external field. These assumptions and a mean-field approximation lead to a
self-consistent system of two equations. For a detailed derivation the reader may refer
to \cite{blinc74}. From this model we get the scalar parameters $K,L$ and $\Delta$ which
determine that the system is stable for a certain set of parameters
 \begin{center}
\begin{tabular}{rcl}
 for     &$ T<T_{c,1}$ or $T>T_{c,2}$ & with $\langle S_1^z\rangle = -\langle S_2^z\rangle$ \\[1ex]
 and for &$T_{c,1}\leq T\leq T_{c,2}$ & with $\langle S_1^z\rangle \neq -\langle S_2^z\rangle$
\end{tabular}
\end{center}
Here $L$ is the dominant term which leads to an antiparallel alignment of the sublattices
for all temperatures and all applied external fields. The asymmetry term $\Delta$ is also
responsible for an antiparallel tendency of the sublattices. The parameter $K$ describes
an antiparallel influence of overnext neighbours resulting in a ferroelectric tendency
for the whole system in the temperature range around the value of $K$.

Describing the physics of the model qualitatively one may explain the mechanism as
follows. The parts of the molecule now called "particles" which are responsible for the
polar order are located in double well potentials. There is no net polarization in the
ground state ({\it i.e.} the dipoles of both sublattices are opposite and equal); and
there are excited states with a finite dipole moment. Furthermore there may be an
interaction between the dipoles of the excited states which favours a parallel alignment.

At $T=0$ all molecules are in the ground state. Therefore there is no total polarization.
If the dipole-dipole interaction is large enough, the excitation energy of the particles
can be compensated by the gain in total energy due to the alignment of the dipoles. In
this scenario the state of ferroelectric order is favoured with respect to the total
energy of the system.

Now, the dipole-dipole interaction may be just below the critical value for this
scenario. For $T>0$ we have to take both the energy and the entropy into account: at a
finite temperature there is a probability for a certain number of particles to be in the
excited state ({\it i.e.} resulting in a net polarization). From this follows, that the
effective magnitude of the polarization is equal to the original magnitude multiplied by
this probability. On the other hand, the influence of the dipole-dipole interaction is
decreasing with rising temperature. Therefore, at $T>0$, a ferroelectric phase is only
possible if the probability for the excited states increases more rapidly than the
influence of the order decreases. At high temperatures the order will vanish in any case.

We will now apply the model to data of two ammonium concentrations: the sample of pure
NaK-tartrate ($x=0$) and the sample with $x=0.15$, the latter being in the range where no
FE ordering occurs any more. We take the data for the lowest frequency measured, namely
20\,Hz, to be close to the static limit of the theoretical model.
Figure~\ref{fig:theo}
\begin{figure}[htbp]
\begin{center}
\includegraphics[clip,width=9.5cm]{fig10.eps}
\caption{Comparison of the measured data for the ammonium concentrations of $x=0$ (a) and
$x=0.15$ (b) with calculated curves according to the model of Blinc and \v{Z}ek\v{s}
\protect\cite{blinc72} in (c) and (d). The lines in (a) are drawn to guide the eye.}
\label{fig:theo}
\end{center}
\end{figure}
shows that the theoretical result is in good qualitative agreement with our measurements.
We yield the parameters shown in Table~\ref{tab:1}.
\begin{table}[h]
\begin{center}
\begin{tabular}{|r|r|r|} \hline
 \,\,parameter\,\, & \,\,$x=0$\,\,        & \,\,$x=0.15$\,\,     \\ \hline
    \,\,K\,\,      & \,\,$271.275$\,K\,\, & \,\,$213.195$\,K\,\, \\ \hline
    \,\,L\,\,      & \,\,$5877.18$\,K\,\, & \,\,$4618.87$\,K\,\, \\ \hline
 \,\,$\Delta$\,\,  & \,\,$1701$\,K\,\,    & \,\,$1342$\,K\,\,    \\ \hline
\end{tabular}
\caption{Parameters for modeling the dielectric spectra for the two doping concentrations
$x=0$ and $x=0.15$ with the {\it pseudo-spin model} described in the text.} \label{tab:1}
\end{center}
\end{table}
From this table it is obvious that the set of parameters does not have to be changed very
much to describe the r{\'e}gime of suppressed FE order. For the description of the doped
sample the effective parameters $K$ and $L$, as well as the effective asymmetry parameter
$\Delta$ are lowered by about $27$\,\%.

\section{CONCLUSIONS}

We have presented results from low and high field frequency-dependent dielectric
measurements on 18 different ammonium doping concentrations of the compound
Na$($K$_{1-x}($NH$_4)_x)$C$_4$H$_4$O$_6\cdot 4$H$_2$O. Characteristic spectra of
concentrations representative for the various phases in the complex phase diagram were
shown. In the low field dielectric measurements a richness of relaxational features was
found over the whole doping range and systematically analyzed. The detected relaxation
processes were classified roughly into three groups. Presumably two of these groups
represent single particle intramolecular motions, whereas the third one is of collective
nature. From the hysteresis loops of the high field measurements we have provided
striking evidence for a new phase  showing anti-ferroelectric correlations in the doping
range of $0.18<x<0.89$. This finding is corroborated by the detection of an anomaly in
the temperature dependent dielectric permittivity obtained from the low field
measurements. A new phase boundary is introduced into the established phase diagram
showing the systematic doping dependence of the new phase. Finally we demonstrated for
the first time that the simple model by Blinc and \v{Z}ek\v{s} can be used to describe
not only the data on pure Rochelle salt, but the whole range up to an ammonium content of
$x=0.18$ at least qualitatively. This is a remarkable fact and may improve the
understanding of the re-entrant behaviour of the pure Rochelle salt vanishing upon
doping.

\section*{ACKNOWLEDGEMENTS}

We want to express our gratitude towards A. Maiazza at the Technische Universit\"{a}t
Darmstadt for growing the high-quality Rochelle salt crystals this work is based upon. We
are indebted to H. Ries for help with the high field non-linear measurements. One of the
authors (U.S.) also wants to acknowledge fruitful discussions with J. Petzelt and S.
Kamba at the Institute of Physics in Prague and A. Volkov at the General Physics
Institute in Moskow.