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

The limitation of the length of a genome by the replication accuracy
per nucleotide ($q$) has led to a deadlock in the theories of the
origin of life based on the evolution of competing self-replicating
polynucleotides.
According to Eigen's quasispecies model \cite{eigen,review}, which may
serve as a paradigm here, polynucleotides have to replicate with high
accuracy in order to reach a certain length, a requirement that is
impossible to fulfill without the aid of specialized catalysts. However, to
build those catalysts a blueprint is necessary that amounts to a large
genome (the nucleotide sequence), which itself cannot be maintained without
the catalysts. In particular, for polynucleotides of fixed length $L$, the
quasispecies model predicts the existence of a minimal replication accuracy
per genome $Q_c =q_c^L$, below which the genetic information is
irreversibly lost. This information crisis has been termed error threshold
transition. Above $Q_c$ the population is composed of a master copy
together with a cloud of structurally similar mutants (quasispecies)
\cite{eigen,review}. Equally important is the finding that, except in a
trivially degenerate case, two or more quasispecies cannot coexist
\cite{swetina}, thus precluding the coexistence of templates
(i.e. polynucleotides) sufficiently different from each other to code for
any useful set of catalysts. Although it has been claimed that the
information crisis is not really a fundamental issue, since the error
threshold transition appears only in some pathological, discontinuous
replication landscapes \cite{wiehe}, the coexistence problem seems to be
more pervasive, as it is associated to the form of the growth functions in
the chemical kinetics equations \cite{tree,wills}.

In order to circumvent the aforementioned limitations of the quasispecies
model, Eigen and Schuster proposed the hypercycle \cite{hyper}, that is, a
catalytic feedback network whereby each template helps in the replication
of the next one, in a regulatory cycle closing on itself. This model has
gain plausibility when the ability of polynucleotides to help propagate
each other was established experimentally through the study of the
catalytic activity of the RNA (ribozymes) \cite{doudna,bio}. Interestingly,
though the error threshold phenomenon has traditionally been considered the
main motivation for the proposal of the hypercycle (see \cite{maynard}, for
instance), most of the seminal works in this field have dealt with the
coexistence issue only, as they assume perfect replication accuracy for the
hypercycle elements \cite{hyper,hofbauer}. In this case an arbitrary number
of templates permanently coexist in a dynamical equilibrium state; if $n>4$,
however, the template concentrations vary with time \cite{hyper},
periodically  decreasing to very small values. In practice,
large hypercycles are therefore susceptible to extinction via fluctuations,
see e.g.\ \cite{nuno0}, hence the information gain due to the coexistence
of different templates in the hypercycle may not be very impressive after
all. Furthermore, we will argue in this paper that coexistence in the
absence of a stable equilibrium can also be achieved by a simpler
arrangement, namely, the free chains, in which the cyclic order of the
catalysts is interrupted.

The effect of error-prone replication (mutation) in the hypercyclic
organization was investigated by introducing a mutation field as a
perturbation of the error-free kinetic equations \cite{stadler}. This
approach, however, is not very appropriate to study the error threshold
phenomenon, since the results obtained cannot be easily compared with those
of the quasispecies model.  In this sense, a better approach is to assume
the existence of a special class of templates with no catalytic activity,
so-called error-tail, that appear as a consequence of the replication
errors of the hypercycle elements \cite{garcia,nuno1,nuno2}. However, the
particular catalytic network investigated extensively within that framework
was not the hypercycle, except for a short discussion in \cite{nuno3}, but
the fully connected network in which each element helps the replication of
all other elements of the network \cite{nuno1,nuno2}.  (Clearly, in the
case of $n=2$ elements, these two networks become identical \cite{garcia}.)
Such a network is more robust than the hypercycle since the mal-functioning
or extinction of one of its elements does not compromise the whole network.
Nevertheless, besides its aesthetic appeal, the cyclic coupling of the
hypercycle seems to be more realistic \cite{bio}.

The goal of this paper is to investigate analytically the
steady-states of a deterministic system comprised of two parts,
namely, a hypercycle made up of $n$ self-replicating templates $I_1,
I_2, \ldots, I_n$ and its error tail $I_e$. These parts are coupled
such that any erroneous copy of the hypercycle elements will belong to
the error tail.  The focus of the present analysis is on the location
in the parameters space of the model (i.e. replication accuracy per
template, non-catalyzed and catalyzed productivity values, and
hypercycle size) of the regions of stability of the diverse
possibilities of coexistence between the templates composing the
hypercycle. In particular, we give emphasis to the characterization of
the critical parameters at which the hypercycle becomes unstable
against the error tail.

The remainder of the paper is organized as follows. In Sec.\
\ref{sec:level2} we present the chemical kinetics equations that
govern the time evolution of the system and motivate the specific
choice of the parameters used throughout the paper. The fixed-points
of the kinetic equations are obtained analytically in Sec.\
\ref{sec:level3} and their stability discussed in Sec.\
\ref{sec:level4}. The phase-diagrams showing the regions of stability
of the diverse coexistence states are presented and analyzed in Sec.\
\ref{sec:level5}.  Finally, some concluding remarks are presented in
Sec.\ \ref{sec:level6}.

\section{The model}\label{sec:level2}

\begin{figure}
\psfig{file=fig1.ps,width=0.47\textwidth}
\caption{The system composed of a hypercycle of size $n=3$
and its error tail $I_e$. The thin arrows represent the
non-catalyzed self-replication reactions and the thick arrows
represent the self-replication catalytically assisted by
the neighbor template}
\label{Fig.1}
\end{figure}

We consider a system composed of a hypercycle made up of $n$ elements $I_1,
\ldots, I_n$ and its error tail $I_e$, as illustrated in
Fig.~\ref{Fig.1}. In contrast to the so-called elementary hypercycle
\cite{hyper}, we assume that the templates are capable of self-replication
with productivity values $A_i ~\left( i=1, \ldots, n \right )$ and
$A_e$. Moreover, as usual, the growth promotion of template $I_i$ as a
result of the catalysis from template $I_{i -1}$ is measured by the kinetic
constants $K_i$. The key ingredient in the modeling is that in both
processes of growth of template $I_i$ the probability of success is given
by the parameter $Q \in [0,1]$, so that an erroneous copy, which will then
belong to the error tail, is produced with probability $1-Q$.  Hence the
concentrations $x_i~\left ( i=1, \ldots,n \right )$ of the hypercycle
elements and the concentration $x_e$ of the error-tail evolve in time
according to the kinetic equations
\begin{equation}\label{ode_h}
\dot x_{i} = x_{i} \left ( A_{i}Q + K_i x_{i-1}Q - \Phi \right )
~~~~~~~ i=1,...,n
\end{equation}
and
\begin{equation}\label{ode_e}
\dot x_{e} = x_{e} \left ( A_{e} - \Phi \right ) + \left ( 1-Q \right )
\sum_{i=1}^{n} x_{i} \left ( A_{i}  + K_i x_{i-1} \right )
\end{equation}
where $x_0 \equiv x_n$ and
\begin{equation}\label{flux}
\Phi= \sum_{i=1}^n x_i \left ( A_i + K_i x_{i-1} \right ) + A_e x_e
\end{equation}
is a dilution flux that keeps the total concentration constant, i.e.,
$ \sum_{i=1}^{n} \dot x_{i} + \dot x_{e}  = 0$. As usual,
the dot denotes a time derivative. Henceforth we will assume
that \cite{nuno1,nuno2}
\begin{equation}\label{CC}
 \sum_{i=1}^{n} x_{i} +x_{e} = 1 .
\end{equation}
Clearly, this formulation is equivalent to considering polynucleotides
of length $L \rightarrow \infty$ whose replication accuracy per
nucleotide $q$ goes to 1 such that the replication accuracy per genome
is finite, i.e.  $q^L \rightarrow Q$. In this limit, the
back-mutations from the error-tail elements to the templates that
compose the hypercycle, as well as the mutations between those
templates, can be safely neglected.  Hence, mutations can only
increase the concentration of the elements in the error-tail.  The
advantage of working in this limit is that the error threshold
transition can be precisely located by determining the value of $Q$ at
which the concentration of a relevant template vanishes. For finite
$L$, as well as for finite population sizes, the characterization of
this transition is more involved, being achieved through the use of
finite-size scaling techniques \cite{prac}.

In this work we consider the single-sharp-peak replication landscape
\cite{eigen,review}, in which we ascribe the productivity value
$A_{1}=a>1$ to the so-called master template $I_1$ and $A_{i} = A_e=1$
to the $n-1$ other elements of the hypercycle as well as to the
error-tail. Also, for the sake of simplicity we set $K_i = K$ for all
$i$. The motivation for this particular choice of parameters is the
observation that the emergence of the hypercycle requires both spatial
and temporal coexistence of the templates forming the network, and
this can be achieved by a quasispecies distribution, which guarantees
the coexistence of the master template and its close mutants, despite
the purely competitive character of the quasispecies model
\cite{emergence}.  Once the coexistence is established, the appearance
of catalytic couplings between the templates is not a very unlike
event. Of course, as soon as those cooperative couplings become
sufficiently strong to balance the competition imposed by the constant
concentration constraint, the mutants will certainly depart from the
master template due to the relentless pressure of the mutations, so
that no trace will remain of the original quasispecies distribution.

\section{Fixed points}\label{sec:level3}

Let us distinguish between surviving templates $x_i>0$ and extinct
templates $x_i=0$.  A survivor $I_j$ is said isolated if
$x_{j-1}=x_{j+1}=0$. Hence,
\begin{equation}\label{j}
\dot x_j = x_j \left ( Q - \Phi \right )  ~~~~ j > 1
\end{equation}
and
\begin{equation}\label{e}
\dot x_e = x_e(1\!-\!\Phi)\!+\!(1\!-\!Q)\!
\left[x_1\left(a\!-\!1\right)\!+\!1\!+\!
K\!\!\sum_{i\neq j,j+1}\!\!\! x_{i}x_{i-1}\right].
\end{equation}
In the steady-state regime, Eq.\ (\ref{j}) yields $\Phi = Q$ which,
for $Q < 1$, is incompatible with Eq.\ (\ref{e}), since the term
within brackets in this equation is positive. Therefore, all isolated
survivors with the exception of the master template are unstable
against the error tail.  Next consider the following chain of
surviving templates:
\begin{equation}
\dot x_i     = x_i( Q - \Phi )
\end{equation}
\begin{equation}
\dot x_{i+1} = x_{i+1}( Q +KQx_{i} - \Phi )
\end{equation}
\begin{equation}
\dot x_{i+2} = x_{i+2}( Q +KQx_{i+1} - \Phi )
\end{equation}
\begin{equation}
\vdots
\end{equation}
\begin{equation}
\dot x_{k}   = x_{k}( Q +KQx_{k-1} - \Phi )
\end{equation}
which does not contain $x_1$. Again, in the steady-state regime the
first equation yields $\Phi = Q$, implying $KQx_{i}=0$, i.e.,
$x_{i}=0$. So, there is no fixed point corresponding to such a
chain. Any chain of survivors therefore must start with template $n$
or $1$. In the first case we get $ \Phi = Q$ from $\dot x_n=0$ and
$aQ+KQx_n=Q$ from $\dot x_1=0$, yielding $x_n=(1-a)/K<0$, which rules
out this possibility.  The equilibria of interest for our study thus
are either the interior equilibrium in which all templates survive, or
a fixed point that corresponds to a chain of survivors beginning with
$I_1$.

Accordingly, we define a $m$-coexistence state as the $n$-component
template vector ${\bf x} = \left( x_1, x_2, \ldots, x_n \right )$ in
which the first $m$ components are strictly positive and the rest are
equal to zero. Clearly, given the template vector ${\bf x}$ the
concentration of the error tail $x_e$ is determined by the constraint
(\ref{CC}).  In the following we solve analytically the kinetic
equations in the steady-state regime $\dot x_i = 0 ~\forall i$.

The simplest fixed point is the zero-coexistence state ($m=0$) which
corresponds to the solution $x_1=\ldots=x_n = 0$ and $x_e= 1$,
existing for the complete range of parameter values.

In the case of chains, i.e. $0 < m < n$, the steady-state solutions of
Eqns. (\ref{ode_h}) and (\ref{ode_e}) are straightforward. In fact,
since $x_n = 0$ by definition, we get $\Phi = aQ$ from $\dot x_1 = 0$
which then yields
\begin{equation}
x_1 = x_2 = \ldots = x_{m-1} = \frac{a -1}{K}.
\end{equation}
Next we insert this result in Eq.\ (\ref{flux}) to obtain
\begin{equation}\label{xm}
x_m = \frac{Qa-1}{a-1} - \left( m-1 \right) \frac{a-1}{K}.
\end{equation}
However, since $x_i \in \left ( 0,1 \right )~\forall i$, this solution
is physically meaningful in the region $K > a-1$ and $Q > Q_m$ where
\begin{equation}\label{Qm}
Q_m = \frac{1}{a} + \frac{1}{Ka}~ \left ( m-1 \right ) \left ( a-1 \right )^2 .
\end{equation}
We note that the $1$-coexistence state (quasispecies) is obtained by
setting $m=1$ in Eq.\ (\ref{xm}) and its region of existence is simply
$Q > 1/a$, since the other condition, namely $K > a-1$, is derived by
considering the other templates in the chain. In fact, this very
simple result quantifies nicely the notion that the cooperative
couplings must reach a certain minimum strength so as to balance the
competition between templates.

The analysis of the hypercycle, i.e. $m=n$, is a little more involved.
>From $\dot x_2 = 0$ we get $\Phi - Q = K Q x_1$ which, inserted in the
equations $\dot x_3 = \ldots = \dot x_n=0$, yields $x_1 = x_2 = \ldots
= x_{n-1}$ and
\begin{equation}\label{co_n}
x_n = x_1 - \frac{a-1}{K} .
\end{equation}
Finally, using these results in Eq. (\ref{flux}) we find that $x_1$
is given by the roots of the quadratic equation
\begin{equation}\label{quad}
n K x_1^2 - \left ( K Q + a -1 \right ) x_1 + 1 - Q = 0 .
\end{equation}
For $ K < \left ( a -1 \right )^2/4n$, this equation has real roots
for all $Q \geq 0$, otherwise it has real roots for $Q \geq Q_h$ where
$Q_h$ is the unique positive root of the equation
\begin{equation}\label{Qh}
K^2 Q_h^2  + 2 K  \left(  a -1 + 2 n \right ) Q_h +
\left ( a -1 \right )^2 - 4 n K = 0 .
\end{equation}
In particular, for large $K$ we find $Q_h \approx 2 \sqrt{n/K}$.
Furthermore, it can be easily seen from Eqns.\ (\ref{co_n}) and
(\ref{quad}) that $x_n$ vanishes at $Q = Q_n$ with $Q_n$ given in Eq.\
(\ref{Qm}). To understand the role of $Q_h$ and $Q_n$ (we note that
$Q_n \geq Q_h$) in delimiting the region of existence of the
$n$-coexistence state we must look at the behavior of the two real
roots of Eq.\ (\ref{quad}).  Let us denote them by $x_1^+$ and $x_1^-$
with $x_1^+ \geq x_1^-$, which, according to Eq.\ (\ref{co_n}),
correspond to $x_n^+$ and $x_n^-$, respectively.  Of course, these
roots become identical at $Q=Q_h$ and so the two solutions for $x_n$
will vanish simultaneously only at the value of $K=K_h$ at which $Q_h$
equals $Q_n$. Explicitly, we obtain
\begin{equation}\label{Kh}
K_h = \left ( a-1 \right ) \left [ n \left ( a+1 \right ) - 1 \right ] .
\end{equation}
by inserting Eq.\ (\ref{Qm}) into Eq.\ (\ref{Qh}).  Although both roots
$x_1^+$ and $x_1^-$ are in the simplex $\left ( 0,1 \right )$, this is not
so for $x_n^+$ and $x_n^-$. In particular, for $K < K_h$ both
concentrations are negative within the range $Q_h \leq Q < Q_n$.  However,
while $x_n^+$ becomes positive for $Q > Q_n$ (it vanishes at $Q_n$),
$x_n^-$ remains always negative.  Since $K_h > \left ( a -1 \right )^2/4n$
the same conclusion holds in the range $K < \left ( a -1 \right )^2/4n$ as
well, provided we define $Q_h = 0$ in this region.  The situation is
reversed for $K > K_h$: both concentrations are positive within the range
$Q_h \leq Q < Q_n$, but now it is $x_n^-$ that vanishes at $Q_n$ and
becomes negative for $Q > Q_n$ while $x_n^+$ remains always positive.
Despite the small region in the parameters space where the root $x_1^-$
yields concentrations inside the simplex, the linear stability analysis
discussed in the sequel indicates that this solution is always unstable, so
we only need to consider the root $x_1^+$.  Thus the range of existence of
the hypercycle fixed point $m=n$ is $Q \geq Q_n$ if $K \leq K_h$ and $Q
\geq Q_h$ if $K > K_h$.

In models without error-tail, i.e., pure replicator equations a much
stronger statement on coexistence is possible. The ``time average theorem''
\cite{hs-book} states that if there is a trajectory along which a certain
subset of templates $J$ survives, then there is a fixed point with exactly
the $J$-coordinates non-zero. While we have not been able to prove the
``time average theorem'' in full generality for Eqns.~(\ref{ode_h}) and
(\ref{ode_e}), it is easily verified for free chains. Hence, if there is
no $m$-coexistence equilibrium, then there is no trajectory at all along
which the templates $I_1$ through $I_m$ survive.

Hitherto we have determined the ranges of the parameter $Q$ where the
$m$-coexistence states are physically meaningful, in the sense that
$x_i \in \left ( 0,1 \right )~\forall i$. The next step is to find the
regions where these states are locally stable.

\section{Stability analysis}\label{sec:level4}

In order to perform a standard linear stability analysis of the fixed
points obtained in the previous section, it is convenient to rewrite
the kinetic equations (\ref{ode_h}) and (\ref{ode_e}) as follows
\begin{equation}
\dot x_i = x_i F_i \left ( {\bf x} \right ) ~~~~i=1, \ldots,n
\end{equation}
where
\begin{equation}\label{F}
F_i({\bf x}) =  A_i Q + K Q x_{i-1} - A_e -\!\sum_j x_j
\left(A_j\!-\!A_e\!+\!K x_{j-1}\right)
\end{equation}
and we have used the constraint (\ref{CC}) to eliminate $x_e$.  The
stability of a fixed point is ensured provided that the real parts of
all eigenvalues of the $n \times n$ Jacobian ${\cal J}$ are negative.
In our case the elements of the Jacobian are given by
\begin{equation}\label{M}
J_{ij} = \delta_{ij} F_i + x_i \frac{\partial{F_i}}{\partial{x_j}}
~~~~i,j=1,\ldots,n .
\end{equation}
The evaluation of the eigenvalues is simple only for the
zero-coexistence state, since in this case the Jacobian is diagonal
with elements $J_{11} = a Q - 1$ and $J_{ii} = Q-1,~ i > 1$. Therefore
this steady state becomes unstable for $Q > 1/a$, which coincides with
the lowest replication accuracy required for the existence of the
$1$-coexistence state. However, for a general $m$-coexistence state we
have to resort to a numerical evaluation of the Jacobian eigenvalues.

Fortunately, in the case of chains $0 < m < n$ there is an alternative
way to look at the stability of the fixed points, as hinted by the
stability analysis of the zero-coexistence state, which becomes
unstable due to the emergence of the $1$-coexistence state. In fact,
it can be easily seen that any perturbation of the $m$-coexistence
fixed point which makes the concentration $x_{m+1}$ non-zero will be
amplified if $A_{m+1}Q + K Q x_m - \Phi$ is positive. For $m>0$ we use
$\Phi = aQ$ and $A_{m+1} = 1$ together with the value of $x_m$ given
in Eq.\ (\ref{xm}) to obtain the following (necessary) condition for
the stability of the $m$-coexistence state,
\begin{equation}\label{stab_c}
Q < Q_{m+1}  ~~~~ m > 0,
\end{equation}
with $Q_m$ given in Eq.\ (\ref{Qm}).  Hence the maximum value of $Q$
allowed for the stability of the $m$-coexistence state coincides with
the minimum $Q$ required for the existence of the $(m+1)$-coexistence
state.  Interestingly, though for $m=0$ we have $\Phi = 1$, $A_1 = a$
and $x_0 = 0$, condition (\ref{stab_c}) holds true in this case too.

At this point two caveats are in order. First, the entire argument
leading to the stability condition (\ref{stab_c}) is flawed if the
$(m+1)$-coexistence state happens to be unstable. Therefore, we must
guarantee via the numerical evaluation of the Jacobian eigenvalues
that the $l$-coexistence state is stable before using that condition
to study the stability of chains with $m < l$.  In particular, we have
carried out the numerical analysis for the hypercycle solution $l=n$
and found the following general results:
\begin{itemize}
\item[(i)] For $ n \leq 4$, it  is always stable;
\item[(ii)] for $ n =5$, it is stable in a very small range of $Q$ above
$Q_5$; and
\item[(iii)] for $n \geq 6$, it is always unstable.
\end{itemize}
Second, the derivation of the stability condition (\ref{stab_c}) is
based on the analysis of a single eigenvalue of the Jacobian and so it
does not yield a sufficient condition for the stability of the fixed
points. Nevertheless, we have verified through the numerical
evaluation of all $n$ eigenvalues that, provided the
$(m+1)$-coexistence state is stable, the eigenvalue associated to
fluctuations leading to an increase of the chain length is the first
one to become positive.

\section{Discussion}\label{sec:level5}

Combining the existence and the stability results derived in the
previous sections we can draw the phase-diagrams in the plane $(K,Q)$
for fixed $a$ and $n$. In particular, for $n \leq 4$ the
$m$-coexistence state is stable within the interval
\begin{equation}\label{chain}
Q_m <Q<Q_{m+1} ~~~~~ m < n
\end{equation}
with $Q_m$ given in Eq.\ (\ref{Qm}), provided that $K > a-1$.
Interestingly, for fixed $Q$, Eq.\ (\ref{chain}) shows that the
increment $\delta K$ in the catalytic coupling needed to incorporate a
new template into the chain is
\begin{equation}\label{DK}
\delta K = \frac{ \left ( a - 1 \right)^2 }{aQ - 1} ,
\end{equation}
regardless the number of elements in the chain.  The case $ K < a -
1$, for which no chains with $m>1$ are allowed, does not require any
special consideration. In fact, we find that the only stable states
are the zero-coexistence ($Q < 1/a$) and the $1$-coexistence ($1/a
\leq Q \leq 1$) states. However, since $Q_2 \leq 1$ only for $K \geq a
-1$, this result is consistent with Eq.\ (\ref{chain}).

\begin{figure}
\psfig{file=fig2.ps,width=0.47\textwidth}
\caption{Phase-diagram in the space $\left ( K,Q \right )$ for $n=4$
and $a=10$ showing the regions of stability of the diverse coexistence
states.  The numbers between parentheses indicate the number of
coexisting templates.  Regions of bistability appear for $K >
K_h$. The thin lines are (from bottom to top) $Q_1, Q_2, Q_3$ and
$Q_4$. The thick line is $Q_h$.}
\label{Fig.2}
\end{figure}

The $n$-coexistence state (i.e., the hypercycle solution) is stable
for $Q > Q_n$ if $K \leq K_h$ and for $ Q > Q_h$, otherwise, where
$K_h$ and $Q_h$ are given by Eqns.\ (\ref{Kh}) and (\ref{Qh}),
respectively.  We define the error threshold of the hypercycle as the
value of the replication accuracy $Q$ that delimits the region of
stability of the $n$-coexistence state. The phase diagram for $a=10$
and $n=4$ shown in Fig.~\ref{Fig.2} illustrates the major role played by $K_h$
in the hypercyclic organization: only for $ K > K_h$ the hypercycle
becomes more stable than a chain of the same length.  Another
important quantity is the value of $K$, denoted by $K_c$, at which
$Q_h$ equals $1/a$, the minimal replication accuracy of the
quasispecies. It is given by
\begin{equation}\label{Kc}
K_c = a \left ( a - 1 \right ) \left [ 2n - 1 + \sqrt{n \left ( n - 1
\right)} \right ] .
\end{equation}
Beyond this value the error threshold of the hypercycle $Q_h$ is
smaller than that of the quasispecies.  Moreover, as mentioned before,
for large $K$, it vanishes like $1/\sqrt{K}$.  A rather frustrating
aspect of $K_h$ and $K_c$ is that both are of order $a^2$, indicating
then that the productivity of catalytically assisted self-replication
is much larger than that of non-catalyzed self-replication. While this
is obviously true for biochemical catalysis, it is difficult to argue
for the existence of such efficient catalysts in prebiotic conditions.
On the other hand, we can take a different, more optimistic viewpoint
and argue that modern biochemical catalysts (enzymes) are so efficient
because their precursors had to satisfy the stringent conditions
imposed by surpassing $K_h$.

\begin{figure}
\psfig{file=fig3.ps,width=0.47\textwidth}
\caption{Same as Fig.~\ref{Fig.2} but for $n=5$. There are no stable fixed
points above the dashed curve. The solid curves are (from bottom to
top) $Q_1, Q_2, Q_3,Q_4,$ and $Q_5$.}
\label{Fig.3}
\end{figure}

In Fig.~\ref{Fig.3} we present the phase-diagram for $n=5$. The main
difference from the previous figure is that the $5$-coexistence state is
stable only within the thin region between $Q_5$ and the dashed curve,
obtained through the numerical evaluation of the Jacobian eigenvalues. As
these curves intersect at some $K \leq K_h$, the $5$-membered hypercycle is
not very interesting, since it has the same characteristics of a chain of
length $m=5$.  To confirm this result we have carried out the numerical
integration of the kinetic equations using the ninetieth-order Runge-Kutta
method. The results are shown in Fig.~\ref{Fig.4}, which illustrates the
time evolution of the concentrations $x_i~ \left ( i=1, \ldots, 5 \right )$
inside and outside the region of stability. Although the behavior pattern
in the region of instability seems periodic, we have not explored
completely the space of parameters to discard the existence of chaotic
behavior. Hence we use the term complex dynamics to label this region in
Fig.~\ref{Fig.3}. We note that the phase-diagram shown in this figure
describes also the regions of stability of hypercycles and chains of size
$n \geq 5$, since $m$-coexistence states with $m >5$ are always unstable.

\begin{figure}
\psfig{file=fig4.ps,width=0.47\textwidth}
\caption{Time evolution of the five concentrations of the templates
composing a hypercycle of size $n=5$ for $a=10$, $Q= 1$, and (a) $K=
37$ (inside the region of stability) and (b) $K= 40$ (outside the
region of stability).  The initial state is $x_i (0) = 0.2 ~\forall
i$.}
\label{Fig.4}
\end{figure}

An interesting limiting case which deserves special attention is the
symmetric hypercycle ($a=1$).  According to the argument put forward
in the beginning of Sec.\ \ref{sec:level3}, the only fixed points in
this case are the zero-coexistence state and the hypercycle, i.e.,
chains are not allowed. Moreover, Eq.\ (\ref{co_n}) yields $x_1 = x_2
= \ldots = x_n$ where $x_1$ is given by Eq.\ (\ref{quad}) with $a$
replaced by $1$.  The analysis of the roots of that quadratic equation
and the numerical evaluation of the Jacobian eigenvalues yield that
the symmetric hypercycle is stable for
\begin{equation}
K > \frac{4 n}{Q^2}~\left ( 1 - Q \right ) ,
\end{equation}
provided that $n \leq 4$. The region of stability observed in Fig.~\ref{Fig.3}
for the $5$-coexistence state does not appear in the symmetric case
$a=1$, so it must be a consequence of the asymmetry in the
productivity values of the non-catalyzed self-replication reaction. We
note that, differently from the asymmetric case ($a > 1$), the
zero-coexistence state is always stable.

For the sake of completeness, we present some results on the
elementary hypercycle ($A_i = 0, ~i=1, \ldots,n$) coupled to an error
tail ($A_e = 1$) via the imperfect catalytically assisted
self-replication. Inserting these parameters into Eq.\ (\ref{F}) and
setting $F_i = 0~\forall i$ yields $x_1 = \ldots = x_n$ with $x_1$
given by the larger root of the quadratic equation
\begin{equation}
K n x_1^2 - \left ( n + K Q \right ) x_1 + 1 = 0 ,
\end{equation}
since we have verified that the smaller root is always unstable.  As
in the symmetric case discussed above, for $ n \leq 4$ the stability
condition coincides with the condition for real $x_1$, namely,
\begin{equation}\label{ele}
Q \geq 2 \sqrt{\frac{n}{K}} - \frac{n}{K}.
\end{equation}
Thus the term in the right hand side of this inequality yields the
error threshold of the elementary hypercycle.

\begin{figure}
\psfig{file=fig5.ps,width=0.47\textwidth}
\caption{Time evolution of the five  concentrations of the templates
composing a free chain of size $n=5$. The parameters and initial state
are the same as for Fig.~\ref{Fig.4}(b).}
\label{Fig.5}
\end{figure}

Before concluding this section, we must note that the chains of size $m$
considered hitherto are bonded in a hypercycle of size $n > m$.  We could
study {\it free} chains of size $n$ as well by simply setting $x_0 = 0$ in
the kinetic equations (\ref{ode_h}) and (\ref{ode_e}).  Not surprisingly,
the results are essentially the same as for bonded chains, with $Q_m$
playing a similar fundamental role in delimiting the regions of stability
of the shorter chains ($m < n$). Although a full discussion of the
stability of the complete chain ($m=n$) is beyond the scope of this paper,
we point out that free chains with $n>4$ are always unstable. Moreover, as
illustrated in Fig.~\ref{Fig.5}, the oscillatory behavior pattern of the
template concentrations, which ensues a dynamic coexistence among all
templates in the chain, is similar to that observed in the hypercycle
(compare with Fig.~\ref{Fig.4}).  In this sense, the free chains seem as
good as the hypercycle to attain that kind of coexistence. However, we must
emphasize that, for sufficiently large $K$ (i.e. $K > K_h$), the fixed
points describing the coexistence of $n \leq 4$ templates in a hypercycle
are much more robust against replication errors than their counterparts in
the free chains.

\section{Conclusion}\label{sec:level6}

Our study of the steady states of an asymmetric hypercycle composed of
$n$ error-prone self-replicating templates indicates that, for $n \leq
4$, the error threshold of the hypercycle ($Q_h$) becomes smaller than
that of the quasispecies ($Q_1 = 1/a$) for catalytic couplings ($K$)
of order of $a^2$, where $a$ is the productivity value of the master
template. In particular, $Q_h$ vanishes like $\sqrt{n/K}$ for large
$K$. Perhaps, even more important is our finding that the asymmetry in
the non-catalyzed self-replication reaction ($a > 1$) entails the
existence of chains of size $n \leq 5$. We note that these chains are
unstable in the symmetric hypercycle as well as in the fully connected
network \cite{nuno1,nuno2}.

Adding to the scenario for the emergence of the hypercycle described
in Sec.\ \ref{sec:level2}, which starts with an isolated quasispecies,
our results indicate that the chains may well be the next step in this
complex evolutionary process.  In fact, according to Eq.\ (\ref{DK})
the strengths of the catalytic couplings needed to form a chain are of
order $a$, while the hypercycle only acquires its desirable stability
characteristics for strengths of order $a^2$ (see Eq.\ (\ref{Kh})).
Although we realize that an evolutionary step leading from chains to
hypercycles is still a major one, it is certainly much more plausible
than a direct transition from quasispecies to hypercycle. In any
event, we think that the emergence of the hypercycle can be explained
as a series of plausible smooth transitions, without need to
postulating the hypercycle as an unique event in prebiotic
evolution. In this vein, this work represents a modest first step to
tackle this fundamental problem within a firm basis.

To conclude we must mention an alternative resolution for the
information crisis in prebiotic evolution which has received some
attention recently \cite{book}, namely, the stochastic corrector model
\cite{Szat}. This model builds on ideas of the classical group
selection theory for the evolution of altruism \cite{altru}, since it
considers replicative templates competing inside replicative
compartments, whose selective values depend on their template
composition.  However, the chemical kinetics equations governing the
dynamics of the templates inside the compartments display a
non-physical (non-integer exponents) dependence of growth on template
concentrations. It seems to us that this basic assumption of the
stochastic corrector model must be relaxed, or at least justified,
before it can be considered an important alternative to the more
traditional approach based on the hypercycle and its variants.

\bigskip

P.R.A.C.\ thanks Prof.\ P.\ Schuster and Prof.\ P.\ F.\ Stadler for
their kind hospitality at the Institut f\"ur Theorestische Chemie, where
part of his work was done, and Prof.\ P.\ E.\ Phillipson for illuminating
discussions.  The work of J.\ F.\ F.\ was supported in part by Conselho
Nacional de Desenvolvimento Cient{\'\i}fico e Tecnol{\'o}gico (CNPq).
The work of P.F.S.\ was supported in part by the Austrian
Fonds zur F{\"o}rderung der Wissenschaftlichen Forschung, Proj.\ No.\
13093-GEN. P.R.A.C.\ is supported by FAPESP.


\section{Introduction}

Since the discovery of high-$T_{\rm c}$ superconductors,
there have been revealed many novel superconducting properties
which can not be expected for conventional superconductors.
One of these is the $d$-wave symmetry of the pair potential
giving rise to
novel properties of phase coherence which has been confirmed
by both
experimental and theoretical studies
in tunneling and Josephson effects
\cite{Sigrist1,Wollman,Hu,Kashiwaya2,Tanaka1,
Iguchi1,Tanaka94,Sigrist2,Tanaka3,Barash2}.
The strong two dimensionality and the anomalous transport
along $c$-axis direction have brought about another
important feature of the existence of the intrinsic Josephson
effects in Bi$_{2}$Sr$_{2}$CaCu$_{2}$O$_{8}$(BSCCO)
\cite{Kleiner1,Kleiner2}.
It was reported that a small BSCCO single crystal behaved like a
series of Josephson junctions.
It is natural to consider a BSCCO superconductor as a stack
of superconducting sheets consisting of CuO bilayers separated by
BiO and SrO layers acting as insulators.
After the discovery of intrinsic Josephson effect by Kleiner
\cite{Kleiner1,Kleiner2},
there are several evidences which support the existence of
Josephson coupling between adjacent CuO bilayers
\cite{Itoh,Matsuda1,Matsuda2,Suzuki,Latyshev}.
Up to now there are several theoretical works about this problem
from the view point of the layered structure \cite{Tachiki}.
However, existing theories \cite{Klemm,Radtke,Chen}
do not treat the two important facts, i.e.,
short coherence length and the strong correlation
on the equal footing.
Taking account of the short coherence length effect,
it is suitable to express BSCCO by the lattice model.
One of the authors  calculated the intrinsic Josephson current
in layered attractive Hubbard model within mean field approximations
where the $s$-wave conventional paring is considered\cite{Tanaka}.
Later, Schmitt $et$ $al.$ \cite{Hanke}
calculated the Josephson current
in the same model
based on the Quantum Monte Carlo
Calculation and established the existence of
intrinsic Josephson current from the numerical calculations
with much more accuracy.
However, in these two papers, strong correlation effect
and the resulting $d$-wave symmetry were not considered.
To overcome this problem,
in this paper,
we study the intrinsic Josephson current using
the layered two-dimensional (2D) $t$-$J$ model,
which includes the important features above mentioned.
\par
The $t$-$J$ model \cite{Zhang1}
is one of the promising models which explain the
low-energy excitations in high-$T_{\rm c}$ superconductors.
Although the
analytic solution of this model have not yet been obtained,
phase diagrams as a function of doping rate $\delta$ and a
superexchange interaction $J$ are numerically studied
at $T=0$
for one-dimension \cite{Ogata1} and two-dimension
\cite{Putikka,Dagotto1,Dagotto2,Yokoyama}.
Especially, in the 2D $t$-$J$ model,
the obtained phase diagram
as a function of doping is consistent with actual
high-$T_{\rm c}$ superconductors \cite{H.F.}.
\par
In this letter, we study the intrinsic Josephson effect in the layered
two-dimensional $t$-$J$ model.
The intrinsic Josephson current which flows perpendicular to the $t$-$J$
planes is calculated and the current-phase relations are obtained for
several doping rates.
The Josephson current
obtained
as a function of the doping concentration $\delta$
has a  maximum near the optimum doping region.
\par

The Hamiltonian of the $m$-th $t$-$J$ plane is written as
\begin{eqnarray}
{\cal H}_{t-J}^{m} &=& -t\sum_{\langle i,j \rangle \sigma}
( c^{m\dag}_{i\sigma} c^{m}_{j\sigma} + {\rm h.c.} )
\nonumber \\
& & + \sum_{\langle i,j\rangle }
(J\mib{S^{m}_{i}}\cdot \mib{S^{m}_{j}}
 - \frac{J_{N}}{4}n^{m}_{i}n^{m}_{j} )
\label{tj}
\end{eqnarray}
where the Hilbert space is defined on the subspace without double
occupancy, $i, j$ refer to planar sites on a square lattice,
and $c_{i\sigma}^{m\dag}, c_{i\sigma}^{m}$ represent
fermion operators within the $m$-th layer.
The spin and the number operators in this layer are defined as
\[
\mib{S^{m}_{i}} = \sum_{\alpha,\beta}c_{i\alpha}^{m\dagger}
(\frac{1}{2}
\mib{\sigma} )_{\alpha\beta}c^{m}_{i\beta}  ~,
~~n^{m}_{i} = \sum_{\sigma}c^{m\dag}_{i\sigma} c^{m}_{i\sigma},
\]
respectively.
There seems to be an agreement in the conclusion that the $c$-axis conductivity
in the normal state generally cannot be explained as a coherent interlayer
transport.
We assume, however, that the coherency among
the layers is restored in the superconducting state.
Thus we allow for the electron hopping between the $m$-th and the $m+1$-th
layers,
\begin{equation}
{\cal H}_{z}^{m} =
-t_{z}( c^{m+1\dag}_{i\sigma} c^{m}_{i\sigma} + {\rm h.c.} ).
\end{equation}
Now the total Hamiltonian of the layered $t$-$J$ model is written as
\begin{equation}
{\cal H} = \sum_{m}({\cal H}_{t-J}^{m} + {\cal H}_{z}^{m})
-\mu\sum_{m,i,\sigma} c_{i\sigma}^{m\dag}c_{i\sigma}^{m}
\end{equation}
where $\mu$ is the chemical potential of the whole system.

For this Hamiltonian, we consider Gutzwiller-type variational wave functions
$P_{G}|\Phi\rangle$ with
$P_{G} = \Pi_{m,i}(1 - n_{i\uparrow}^{m}n_{i\downarrow}^{m})$
being the Gutzwiller projection operator which excludes the double occupancy,
and $|\Phi\rangle$ being a one-body mean-field wave function.
Due to the Gutzwiller projection,
it is usually difficult to calculate the variational energy.
Thus, we use a Gutzwiller approximation in which the effect of the projection
is represented by statistical weight factors.
For the two-dimensional $t$-$J$ model, it is known that
this approximation and the variational Monte Carlo simulation give very
similar variational energies \cite{Yokoyama}.

In the Gutzwiller approximation,
the expectation values of the terms in ${\cal H}_{t-J}^{m}$ are estimated as
\begin{eqnarray}
\langle c^{m\dag}_{i\sigma}c_{j\sigma}^{m} \rangle
&=& g_{t}\langle c^{m\dag}_{i\sigma}c_{j\sigma}^{m}\rangle_{0},
\nonumber \\
\langle  \mbox{\boldmath $S_{i}^{m}\cdot S_{j}^{m}$} \rangle
&=& g_{s}\langle \mbox{\boldmath $S_{i}^{m}\cdot S_{j}^{m}$}\rangle_{0},
\nonumber \\
\langle n_{i}^{m}n_{j}^{m}\rangle
 &=& g_{n}\langle n_{i}^{m}n_{j}^{m}\rangle_{0},
\end{eqnarray}
where $\langle\cdots\rangle$ and $\langle\cdots\rangle_{0}$
represent the expectation values in terms of
$P_{G}|\Phi\rangle$ and $|\Phi\rangle$, respectively.
The renormalization factors $g_{t}, ~g_{s}$ \cite{Zhang2} and $g_{n}$
\cite{Yokoyama} are
determined by the ratios of the probabilities of the corresponding
physical processes in the states $P_{G}|\Phi \rangle$ and $|\Phi \rangle$:
\begin{equation}
g_{t} = \frac{2\delta}{1+\delta},
~g_{s} = \frac{4}{(1+\delta )^{2}},
~g_{n} = 1,
\end{equation}
where $\delta = 1 - n$ is the doping rate.
Since the double occupancy is excluded in each layer, we can naturally
assume that the interlayer hopping term is also estimated as
\begin{equation}
\langle c^{m\pm 1\dag}_{i\sigma}c_{j\sigma}^{m} \rangle
= g_{t}\langle c^{m\pm 1\dag}_{i\sigma}c_{j\sigma}^{m}\rangle_{0}.
\end{equation}

Next, we consider the relation among the superconducting phases of the layers.
We introduce order parameters in the $m$-th layer as
\begin{equation}
 \Delta_{\tau}^{m} = \langle c_{i+\tau,\uparrow}^{m\dag}
 c_{i,\downarrow}^{m\dag}\rangle_{0},
\end{equation}
with $\tau = x$ and $y$, $i+\tau$ denotes the nearest neighbor of $i$ in
the $\tau$ direction.
Since we consider the $d_{x^{2}-y^{2}}$-wave symmetry,
$\Delta_{x}^{m} = -\Delta_{y}^{m}$.
When the Josephson current flows along the $z$-direction, the order parameters
satisfy the following phase relation:
\begin{equation}
\Delta_{\tau}^{m} = \Delta_{\tau}^{0}\exp (im\phi).
\end{equation}
Thus, the fermion operator is subject to the following condition:
\begin{equation}
 c_{i\sigma}^{m\pm 1\dag} = \exp(\pm\frac{i}{2}\phi)c_{i\sigma}^{m\dag}.
\end{equation}

Using the Gutzwiller approximation and eq.(9), we can obtain the effective
mean-field Hamiltonian of (3) in momentum space as
\begin{eqnarray}
{\cal H}_{\rm eff} &=& -\sum_{k,\sigma}
T_{k}c_{k\sigma}^{\dag}c_{k\sigma} + N_{s}\epsilon_{0} \nonumber \\
& & - \sum_{k}\left[ \Delta_{k_{\rm 2D}}c_{-k\downarrow}c_{k\uparrow}
 + \Delta_{k_{\rm 2D}}^{\ast}
c_{k\uparrow}^{\dag}c_{-k\downarrow}^{\dag} \right],
\end{eqnarray}
where $k = (k_{x}, k_{y}, k_{z})$ and
$k_{\rm 2D} = (k_{x}, k_{y})$
are the three- and
two-dimensional wavevectors, respectively,
$k_{\rm 2D} = (k_{x}, k_{y})$ is the two-dimensional wavevectors,
$N_{s}$ is the total number of sites,
\begin{eqnarray*}
T_{k} &=& (g_{t}t+J_{1}\xi)\gamma_{k_{\rm 2D}}
 - 2g_{t}t_{z}\cos(k_{z}a_{z}-\frac{\phi}{2})
\nonumber \\
& & + J_{N}n + \mu, \nonumber \\
\Delta_{k_{\rm 2D}} &=& J_{2}|\Delta_{\tau}^{0}|\eta_{k_{\rm 2D}}, \nonumber \\
\epsilon_{0} &=& 4J_{2}|\Delta_{\tau}^{0}|^{2} + 4J_{1}\xi^{2}
+ \frac{J_{N}}{2}n^{2}, \nonumber \\
\gamma_{k_{\rm 2D}} &=& 2(\cos(k_{x}a_{x}) + \cos(k_{y}a_{y}) ), \nonumber \\
\eta_{k_{\rm 2D}} &=& 2(\cos(k_{x}a_{x}) - \cos(k_{y}a_{y}) ), \nonumber \\
J_{1} &=& \frac{3}{4}g_{s}J - \frac{1}{4}J_{N},
~~~J_{2} = \frac{3}{4}g_{s}J + \frac{1}{4}J_{N}, \nonumber
\end{eqnarray*}
and $a_{x}, a_{y}, a_{z}$ are the lattice constants.

This Hamiltonian can be diagonalized using a standard Bogoliubov transformation:\begin{eqnarray}
{\cal H}_{\rm eff} &=& \sum_{k}\left( E_{k} - 2g_{t}t_{z}\sin k_{z}a_{z}\sin
\frac{\phi}{2}\right)\left( A_{k}^{\dag}A_{k} + B_{-k}^{\dag}B_{-k} \right)
\nonumber \\
& & + \sum_{k}( \xi_{k} - E_{k} ) + N_{s}\epsilon_{0},
\end{eqnarray}
where
\begin{equation}
E_{k} = \sqrt{\xi_{k}^{2} + |\Delta_{k_{\rm 2D}}|^{2} },
~~~\xi_{k} = -\frac{1}{2}( T_{k} + T_{-k} )
\end{equation}

The parameters $\Delta_{\tau}^{0}, ~\xi, ~\mu$ are determined by
the following self-consistent equations:
\begin{eqnarray}
1 &=& \frac{J_{2}}{4N_{s}}\sum_{k}\frac{\eta_{k_{\rm 2D}}^{2}}
{2E_{k}},
\\
\xi &=& -\frac{1}{4N_{s}}\sum_{k}\frac{\xi_{k}\gamma_{k_{\rm 2D}}}
{2E_{k}},
\\
\delta &=&  \frac{1}{N_{s}}\sum_{k}\frac{\xi_{k}}{E_{k}}.
\end{eqnarray}
Then we obtain the ground state energy per site as the function of $\phi$;
\begin{eqnarray}
E_{0}(\phi)/N_{s} &=& \frac{1}{N_{s}}\sum_{k}(\xi_{k} - E_{k}) \nonumber \\
& & + 4J_{2}|\Delta_{\tau}^{0}|^{2} + 4J_{1}\xi^{2} + \frac{J_{N}}{2}n^{2}.
\end{eqnarray}
In the following calculation, we take $J/t = J_{N}/t = 0.2$  and
~$t_{z}/t = 0.1$.

First, we look at the $\phi$-dependence of $E_{0}(\phi)/N_{s}$.
In order to see the variation in the energy clearly, we define
${\it{\Delta}} E(\phi) = ( E_{0}(\phi) - E_{0}(0) )/N_{s}$.
In Fig.2, $\it{\Delta} E(\phi)$ is plotted against the phase $\phi$ for
various values of the doping rate $\delta$.
We can see that $\it{\Delta} E(\phi)$ is the monotonic increasing
 function of $\phi$,
and that $\it{\Delta} E(\phi)$ increases with $\delta$ for arbitrary
 values of $\phi$.
When $\delta\geq 0.25$ ( the line (e) ), however, $E(\phi)$ exceeds
the energy of the normal state, i.e., the superconducting state becomes unstable
in the region $\phi\geq 0.3\pi$.
Let us discuss this situation a little more.
From eq.(6), we can interpret $g_{t}t_{z}$ as the effective amplitude
of the interlayer hopping.
Thus, as $\delta$ increases, the coupling strength between layers increases.
On the other hand, the gap amplitude $\Delta_{\tau}^{0}$
linearly decreases with
increasing of $\delta$ \cite{Zhang2}.
If the gap amplitude $\Delta_{\tau}^{0}$
is smaller than the effective interlayer
hopping amplitude $g_{t}t_{z}$, the quasiparticle energy
\begin{equation}
 \sqrt{\xi_{k}^{2}+|\Delta_{k_{\rm 2D}}|^{2}}
 - 2g_{t}t_{z}\sin k_{z}a_{z}\sin\frac{\phi}{2}
\end{equation}
becomes
negative for $\phi\geq\phi_{c}$, where $\phi_{c}$ is a critical phase.
Consequently, the superconducting state is no more stable
as compared to the normal state.
\begin{figure}
\vspace{20pt}
\begin{center}
\leavevmode
 \epsfxsize=80mm
 \epsfbox{fig1.eps}
\end{center}
\caption{
The $\phi$ dependence of ${\it{\Delta}}(\phi)$
for (1) $t_{z}/t=0.1$.
We take (a) $\delta=0.05$, (b) $\delta=0.1$, (c) $\delta=0.15$,
(d) $\delta=0.2$ and (e) $\delta=0.25$.
}
\end{figure}

Next, we calculate the Josephson current defined as
\begin{equation}
j(\phi) = \frac{2e}{\hbar}\frac{\partial}{\partial\phi}E_{0}(\phi).
\end{equation}
Figure 2 shows the Josephson current-phase relations
obtained for
various values of $\delta$.
The currents show the sinusoidal $\phi$-dependence
as expected
and have their maximum values $j_{\rm max}(\delta)$ at
$\phi = \frac{\pi}{2}$.
When $\delta\leq0.2$, $j_{\rm max}$ increases with $\delta$.
This is because the effective mass of a Cooper pair, that is,
its localization tendency in a layer decreases in rough proportion to
$1/|\Delta_{\tau}^{0}|$ \cite{Tanaka},
with $|\Delta_{\tau}^{0}|$ being the decreasing function of $\delta$.
On the contrary, when $\delta\geq0.25$, $j_{\rm max}$
decreases with increasing of $\delta$ because the Josephson current
vanishes for $\phi\geq\phi_{c}$.
Thus $j_{\rm max}$ reaches its maximum in the region
$0.2\leq\delta\leq0.25$.
In order to see this situation in more detail,
the doping dependence of $j_{\rm max}$ is plotted explicitly in Fig.3.
\begin{figure}
\vspace{20pt}
\begin{center}
\leavevmode
 \epsfxsize=80mm
 \epsfbox{fig2.eps}
\end{center}
\caption{
The $\phi$ dependence of the Josephson current $j(\phi)/t$
for $t_{z}/t=0.1$.
We take (a) $\delta=0.05$, (b) $\delta=0.1$, (c) $\delta=0.15$,
(d) $\delta=0.2$ and (e) $\delta=0.25$.
}
\end{figure}
\begin{figure}
\vspace{20pt}
\begin{center}
\leavevmode
 \epsfxsize=80mm
 \epsfbox{fig3.eps}
\end{center}
\caption{
The maximum Josephson current $j_{\rm max}(\delta)/t$ plotted as
a function of $\delta$ for $t_{z}/t=0.1$.
}
\end{figure}
We can see that $j_{\rm max}$ has its maximum value
at $\delta\simeq 0.23$ in the present case. \par
In the recent experiment\cite{Matsuda2},
it has been reported that the intrinsic
Josephson current increases with hole-doping (in our notation, $\delta$)
when the system is in the under-doped region.
We note here that, if we do not assume
the doping dependence of interlayer hopping
as in Eq. (10), $i.e.$,
the interlayer hopping amplitude has no doping dependence,
$j(\phi)$ becomes a decreasing function of $\delta$.
Since the renormalization of interlayer hopping  in
eq.(10) is a natural assumption to treat the layered $t$-$J$ model
based on the Gutzwiller approximation,
the doping dependence of the intrinsic Josephson current observed
in actual experiments can be used
to check the validity of the Gutzwiller approximations
of the layered $t$-$J$ model.
Of course there is only a few experiments concerning a doping dependence of the
Josephson current so far, so that more extensive
comparison between the theory and experiments is necessary.

In this letter, we have calculated the d.c. Josephson
 current which flows perpendicular to the
CuO$_{2}$ planes in high $T_{C}$ superconductors using the
layered $t$-$J$ model within the mean field theory and Gutzwiller
approximation. The maximum Josephson current has its maximum
near the optimum doping rate.
This is due to the competition between the enhancement of the
effective transfer energy along the $c$-axis and the
decrease of the magnitude of the order parameter
with the increase of the doping rate.
We hope such a behavior will be observed in experiments near future.

\section*{Acknowledgment}
The authors wish to thank Y. Matsuda and M. Ogata for useful discussions.
This work is supported by
the Core Research for
Evolutional Science and technology (CREST) of the Japan Science
and Technology Corporation (JST).
One of the author Y.T. is supported by
Grant-in Aid for Scientific Research from the
Ministry of Education, Science, Sports and Culture.
The computational aspect of this work has been performed at the
facilities of the Yukawa Institute.


\section{Introduction}

One of the major unresolved puzzles in the Standard Model (SM) is the
observed hierarchy in the masses and mixings of quarks and leptons.  The SM readily
accommodates the measured values, but at the
expense of introducing Yukawa couplings which span in strength from $10^{-6}$ (for the
electron) to $1$ (for the top quark).  An understanding of the small couplings
involved is highly desirable and has been much sought after.

In this Letter we propose a new approach towards resolving the fermion
mass hierarchy puzzle.  We show that an extremely good fit to all the masses and
mixings can be obtained by using higher dimensional operators
involving the relevant fermion fields and
successive powers of the SM Higgs doublet field.
We do not employ any singlet scalar field.
These non--renormalizable
operators will have inverse mass dimensions, but the dimensionless couplings
which multiply them can all be of order one.  Thus the small Yukawa couplings
of the SM will emerge naturally in our approach.

Our proposal differs from most other approaches \cite{other} which address the fermion mass
hierarchy puzzle in one crucial way:  Since we do not utilize vacuum expectation
values (VEVs) of scalar fields that are singlets of SM, the scale of flavor
physics cannot be much above the weak scale.  In fact, from our fits to the fermion
masses and mixings, we find that the scale $M$ where flavor symmetries becomes
manifest, is in the range $(1-2)$ TeV.  While such a scale is high enough to be
consistent with all the experimental constraints, it is within reach
of planned accelerators in the
near future for direct detection.  Had we used VEVs of singlet scalar
fields to explain the mass hierarchy as in most other approaches, the scale of
flavor physics could be (and usually is) near the Planck mass.

In our approach, although the spectrum of the theory below a TeV is that of
the minimal Standard Model, there are a variety of new signals associated with
the Higgs boson that are distinct from the SM.  As noted, the effective Lagrangian involves
higher powers of the SM Higgs doublet field.  As a consequence, the couplings of the
Higgs boson to light fermions are enhanced by a flavor--dependent
numerical factor relative to that of the SM.  These enhancement factors are:
$(3,5,7)$ for ($b,s,d)$ and $(\tau, \mu,e)$; and $(1,3,7)$ for $(t,c,u)$.
These numbers are related to the dimension of the operators from which a particular
fermion receives its mass.  The enhanced $\bar{b}bh^0$ coupling (by a factor 3)
would have a significant impact on Higgs detection at hadron colliders; the
enhancement in the muonic coupling implies that the Higgs production cross section
at a muon collider \cite{muon} will be increased by a factor of 25.

A further consequence of using successive higher dimensional operators
is that there are flavor changing neutral current processes
mediated by the SM Higgs boson.  The mass matrix and the Yukawa coupling
matrix do not diagonalize simultaneously, although there is only a single
Higgs doublet field.  We show that all the low energy constraints, such as from
$K^0-\bar{K^0}$ mass difference, are well satisfied.  The model predicts
$D^0-\bar{D^0}$ mixing to be close to the present experimental limit.  An
improvement by a factor of $(10-20)$ could test the model, if the Higgs
boson mass is not much above $300$ GeV.

The proposed scenario predicts an interesting new discovery channel for the Higgs boson
at the upgraded Tevatron and the LHC.  It uses the flavor changing Higgs vertex $\bar{t}ch^0$
which has a strength  comparable to that of the
flavor conserving $\bar{b}bh^0$ vertex.
If the mass of $h^0$ is below 150 GeV, then in the $\bar{t}t$ production at colliders,
the decay of the top quark
can lead to Higgs boson signals with $h^0 \rightarrow \bar{b}b$ providing the
invariant mass of the Higgs while the $W+\bar{b}$ decay of the $\bar{t}$ tagging the event.
If the Higgs mass is above $m_t$ ($m_{h^0} \approx
(200-350)$ GeV), the production
process $gg \rightarrow h^0$ followed by $h^0 \rightarrow \bar{t}c+\bar{c}t$
can lead to a new signal at LHC.

In the following we will describe in more detail the scheme and its predictions
just outlined.  We will also comment on some possible ways of realizing the effective
Yukawa Lagrangian from renormalizable theories with flavor symmetries broken at the TeV scale.

\section{Fermion mass hierarchy from higher dimensional operators with SM Higgs}

We assume that some flavor symmetries prevent the direct Yukawa coupling of the SM Higgs doublet
($H$) to the light fermions.  These flavor symmetries are spontaneously broken at a scale
$M$, which will be seen to be in the $(1-2)$ TeV range.  The effective theory below
$M$ is the SM with one Higgs doublet but with non--renormalizable terms in the
Higgs Yukawa couplings.  For the quark sector the effective Yukawa Lagrangin
is taken to be (in an obvious notation for the quark fields and with
$\tilde{H} \equiv -i \tau_2 H^*$):
\begin{eqnarray}
{\cal L}^{\rm Yuk} &=& h_{33}^u Q_3 u_3^c \tilde{H} + \left({H^\dagger H \over M^2}\right)
\left(h_{33}^d Q_3 d_3^c H + h_{22}^u Q_2 u^c_2 \tilde{H}+h_{23}^u Q_2 u_3^c \tilde{H}
+ h_{32}^u Q_3 u_2^c \tilde{H}\right) \nonumber \\
&+& \left({H^\dagger H \over M^2}\right)^2
\left(h_{22}^dQ_2 d_2^c H + h_{23}^d Q_2 d_3^c H + h_{32}^d Q_3 d_2^c H+ h_{12}^uQ_1u_2^c
\tilde{H} + h_{21}^uQ_2u_1^c \tilde{H} \right. ) \nonumber \\
&+& \left.  h_{13}^uQ_1 u_3^c \tilde{H} + h_{31}^u Q_3 u_1^c
\tilde{H} \right)
+ \left({H^\dagger H \over M^2}\right)^3 \left(h_{11}^u Q_1 u_1^c \tilde{H}
+ h_{11}^dQ_1 d_1^c H  \right.  \nonumber \\
&+& \left. h_{12}^dQ_1 d_2^c H + h_{21}^d Q_2 d_1^c H + h_{13}^d Q_1 d_3^c H
+ h_{31}^d Q_3 d_1^c H \right ) + H.C.
\end{eqnarray}
Only the top quark has a renormalizable Yukawa coupling, other couplings are
suppressed by successive powers of $(H^\dagger H/M^2)$.
This provides the
small expansion parameter needed to explain the light fermion masses.  We will show that
an good fit to all the fermion masses and mixing angles can be obtained
with the dimensionless couplings ($h_{ij}^{u,d}$) in Eq. (1) taking values
of order one.  Note that $SU(2)_L$ invariance
prevents even number of Higgs doublet fields in the expansion.  In each term
of Eq. (1) there is a unique $SU(2)$ group contraction.

In most other attempts to explain the mass hierarchy, singlet scalar
fields ($S_i$) are employed and the expansion parameter is
$\left\langle S_i \right \rangle/M$.  In such cases, both $\left\langle S_i
\right \rangle$ and $M$ both can be large, for eg., near the Planck scale.  As long as
$\left\langle S_i\right\rangle/M$ is finite (i.e., not extremely small), a fit
to the fermion masses can be obtained.
There is no a priori
reason in such cases
for the flavor symmetry to be broken near the weak scale.  The low energy
theory will be identical to the SM with no modification in the Higgs boson
interactions.  In contrast, in our
approach, since no singlet fields are employed,
the expansion parameter is $\left\langle H^\dagger H \right\rangle/M^2$,
which implies that $M$ cannot be much above the weak scale.  It also results
in new interactions of the SM Higgs boson, which can be directly tested.

We will comment on possible ways of deriving the effective Lagrangian of Eq. (1)
from renormalizable theories towards the end of the paper.  We remark here
that the coefficients $h_{ij}^{u,d}/M^2$ could be thought of as background fields
which carry flavor quantum numbers.  In writing Eq. (1), we have allowed for
the possibility that the underlying flavor structure might be parity invariant.
Identify $h_{ij} = h_{ji}^*$ in that case.  The  dimension of the
various operators in Eq. (1) are
determined by examining the magnitudes of the
quark masses and mixings.  We have allowed for this to be as low as
it can be, consistent with our demand that the dimensionless coefficients be
order one.  This is the most general case.
Interesting special cases which might arise due to the constraints of
flavor symmetry such
as $h_{13}^{u,d}=h_{31}^{u,d} = 0$ or $h_{ij}^d = 0 (i \neq j),$ can be readily
incorporated into our analysis.  The effective Lagrangian inducing the charged lepton
masses is taken to have a form identical to that for the down--type quarks (replace
the couplings $h_{ij}^d$ by $h_{ij}^\ell$ in Eq. (1) for the leptons).

Writing $H = (h^0/\sqrt{2}+v,~~0)^T$ in unitary gauge with $v \simeq 174$ GeV,
and defining a small parameter
\begin{equation}
\epsilon \equiv { v \over M},
\end{equation}
Eq. (1) leads to the following mass matrices for the up--type and the down--type
quarks:
\begin{eqnarray}
M_u = \left(\matrix{h_{11}^u \epsilon^6 & h_{12}^u \epsilon^4 & h_{13}^u \epsilon^4 \cr
h_{21}^u\epsilon^4 & h_{22}^u \epsilon^2 & h_{23}^u \epsilon^2 \cr
h_{31}^u \epsilon^4 & h_{32}^u \epsilon^2 & h_{33}^u}\right)v, ~~~~~
M_d = \left(\matrix{h_{11}^d \epsilon^6 & h_{12}^d \epsilon^6 & h_{13}^d \epsilon^6 \cr
h_{21}^d\epsilon^6 & h_{22}^d \epsilon^4 & h_{23}^d \epsilon^4 \cr
h_{31}^d \epsilon^6 & h_{32}^d \epsilon^4 & h_{33}^d \epsilon^2}\right)v~.
\end{eqnarray}
The charged lepton mass matrix is obtained from $M_d$ by replacing
$h_{ij}^d \rightarrow h_{ij}^\ell$.  The couplings $h_{ij}$ are complex
in general, but $\epsilon$ can be taken to be real without loss of generality.

The masses of the quarks and leptons can be read off from Eq. (3) in the
approximation $\epsilon \ll 1$.  They are:
\begin{eqnarray}
\{m_t, m_c, m_u\}  &\simeq& \{|h^u_{33}|,~ |h^u_{22}|\epsilon^2,~ |h^u_{11} -
h^u_{12}h_{21}^u/h_{22}^u|\epsilon^6\}v ,\nonumber \\
\{m_b, m_s, m_d\} &\simeq& \{|h_{33}^d| \epsilon^2, ~|h_{22}^d| \epsilon^4,~
|h_{11}^d|\epsilon^6\}v, \nonumber \\
\{m_\tau, m_\mu, m_e\} &\simeq& \{|h_{33}^\ell| \epsilon^2,~ |h_{22}^\ell| \epsilon^4,~
|h_{11}^\ell|\epsilon^6\}v.
\end{eqnarray}
The quark mixing angles are found to be:
\begin{eqnarray}
|V_{us}| &\simeq& \left|{h_{12}^d \over h_{22}^d}-{h_{12}^u \over h_{22}^u}\right| \epsilon^2,
\nonumber \\
|V_{cb}| &\simeq& \left|{h_{23}^d \over h_{33}^d}-{h_{23}^u \over h_{33}^u}\right|
\epsilon^2 ,\nonumber \\
|V_{ub}| &\simeq& \left|{h_{13}^d \over h_{33}^d}-{h_{12}^u h_{23}^d \over h_{22}^uh_{33}^d}
- {h_{13}^u \over h_{33}^u}\right|\epsilon^4.
\end{eqnarray}

In order to see how well our scheme explains the observed
mass and mixing hierarchy, we choose an illustrative
set of input values for the quark masses \cite{gasser}: $m_u(1~{\rm GeV}) = 5.1~
{\rm MeV},~m_c(m_c) = 1.27~{\rm GeV},~ m_t^{\rm phys} = 175~{\rm GeV},~m_d(1~{\rm GeV})
= 8.9~{\rm MeV}, ~m_s(1~{\rm GeV}) = 175~{\rm MeV},~m_b(m_b) = 4.25~{\rm GeV}$, along
with the well known charged lepton masses.  We extrapolate all masses to a common
scale, conveniently chosen as $m_t(m_t) \simeq 166$ GeV using three loop QCD and
one loop QED beta functions.  With $\alpha_s(M_Z) = 0.118$, the values of the running
masses at $\mu = m_t(m_t)$ are found to be:
\begin{eqnarray}
\{m_t, m_c, m_u\} &\simeq& \{166,~ 0.60,~ 2.2 \times 10^{-3}\}~{\rm GeV}, \nonumber \\
\{m_b, m_s, m_d\} &\simeq& \{2.78,~ 7.5 \times 10^{-2},~ 3.8 \times 10^{-3}\}~{\rm GeV},
\nonumber \\
\{m_\tau, m_\mu, m_e\} &\simeq& \{1.75,~ 0.104,~ 5.01 \times 10^{-4}\}~{\rm GeV}.
\end{eqnarray}
An optimal choice of the expansion parameter is $\epsilon = 1/6.5$.  ($\epsilon \simeq
1/7-1/6$ gives reasonable fits.)
This corresponds to a relatively low scale of flavor symmetry breaking: $M \simeq 1.1$
TeV.  With $\epsilon = 1/6.5$, we determine the dimensionless coefficients $h_{ij}^{u,d,\ell}$
from the fermion masses quoted in Eq. (6).  They are:
\begin{eqnarray}
\{|h_{33}^u|, |h_{22}^u|, |h_{11}^u-h_{12}^u h_{21}^u/h_{22}^u|\} &\simeq& \{0.96, 0.14,
0.95\}, \nonumber \\
\{|h_{33}^d|, |h_{22}^d|, |h_{11}^d|\} &\simeq& \{0.68, 0.77, 1.65\}, \nonumber \\
\{|h_{33}^\ell|, |h_{22}^\ell|, |h_{11}^\ell|\} &\simeq& \{0.42, 1.06, 0.21\}.
\end{eqnarray}
We see that, remarkably, all couplings are of order unity.  The largest deviation
from one is the charm Yukawa coupling $|h_{22}^u| \simeq 0.14$.  This small
fluctuation actually goes in the right direction to explain the magnitude of the
Cabibbo angle.  From the expression for $|V_{us}|$ (Eq. (5)), with $|h_{22}^u| \simeq
0.14$, and $h_{12}^{u,d}$ of order one, the correct value of
$|V_{us}|\simeq 0.2$ follows naturally.  (In our scheme, $|V_{us}|$ is a linear
combination of
${\cal O}(\sqrt{m_u/m_c})$ and  ${\cal O} (m_d/m_s)$ terms.)
Similarly, $|V_{ub}| \sim 7 \epsilon^4 \simeq 0.004$, where
the factor $7$ is due to $1/h_{22}^u$ enhancement in
the second term of $|V_{ub}|$ (see Eq. (5)).
$|V_{cb}| \simeq 0.04$ can be fit for example, by choosing
$|h_{23}^u| \simeq 1.4$ (if that term dominates) or with $|h_{23}^d| \simeq 0.84$
(if it dominates).  We thus see that the fit to all the quark and lepton masses
as well as the mixing angles is
extremely good with all couplings taking order one values.
We regard the success of the
fit to be a strong enough motivation to take
the proposed scheme  seriously and study its
experimental consequences.  This is what we do in the discussions that follow.

\section{Enhanced Yukawa couplings and Higgs mediated FCNC}

The Yukawa coupling matrices of the SM Higgs boson to the up--type and
the down--type quarks that follow from Eq. (1) are:
\begin{eqnarray}
Y_u = \left(\matrix{7h_{11}^u \epsilon^6 & 5h_{12}^u \epsilon^4 & 5h_{13}^u \epsilon^4 \cr
5h_{21}^u\epsilon^4 & 3h_{22}^u \epsilon^2 & 3h_{23}^u \epsilon^2 \cr
5h_{31}^u \epsilon^4 & 3h_{32}^u \epsilon^2 & h_{33}^u}\right), ~~~~~
Y_d = \left(\matrix{7h_{11}^d \epsilon^6 & 7h_{12}^d \epsilon^6 & 7h_{13}^d \epsilon^6 \cr
7h_{21}^d\epsilon^6 & 5h_{22}^d \epsilon^4 & 5h_{23}^d \epsilon^4 \cr
7h_{31}^d \epsilon^6 & 5h_{32}^d \epsilon^4 & 3h_{33}^d \epsilon^2}\right),
\end{eqnarray}
with the charged lepton Yukawa coupling matrix $Y_\ell$ obtained from $Y_d$
by replaing $h_{ij}^d \rightarrow h_{ij}^\ell$.

There are two striking features in Eq. (8).  One is that the diagonal couplings
are enhanced relative to the respective SM Higgs boson Yukawa couplings by a numerical factor.
These enhancement factors are $(1,3,5)$ for $(t,c,u)$ and $(3,5,7)$ for
$(b,s,d)$ as well as $(\tau, \mu, e)$.  In the usual SM, these factors are all
equal to one.  In two Higgs doublet models or in supersymmetric models, while
the couplings of $h^0$ might be enhanced, they are not flavor dependent, nor do they
take these specific values.
The use of higher dimensional operators in our approach
make them flavor--dependent and different from one.  The second feature is related
to the flavor dependence which has
an important consequence: The Yukawa coupling matrix and the corresponding
mass matrix do not diagonalize simultaneously.  This will lead to
Higgs mediated flavor changing neutral current processes, even though there
is only a single electroweak Higgs boson in the theory \cite{2higgs}.

The enhanced fermionic Yukawa couplings has a strong influence on the strategy for discovering
the Higgs boson.  We will discuss this in the next Section.  First we will address
the FCNC processes and establish that they are consistent with low energy
data such as $K^0-\bar{K^0}$ mass difference.

The FCNC couplings of the SM Higgs boson $h^0$ to fermions can be obtained from Eqs. (8) and
(3).  In the quark sector in terms of the mass eigenstates they are given by:
\begin{eqnarray}
{\cal L}^{\rm FCNC} &\simeq& {h^0 \over \sqrt{2}}(2h_{12}^d \epsilon^6 ds^c+
2h_{21}^d \epsilon^6 sd^c+4h_{13}^d \epsilon^6
db^c + 4h_{31}^d \epsilon^6 bd^c + 2h_{23}^s \epsilon^4 sb^c+2h_{32}^d \epsilon^4 bs^c)
\nonumber \\
&+& {h^0 \over \sqrt{2}}(2h_{12}^u \epsilon^4 uc^c+
2h_{21}^u \epsilon^4 cu^c+4h_{13}^d \epsilon^4
ut^c + 4h_{31}^u \epsilon^4 tu^c + 2h_{23}^u \epsilon^2 ct^c+2h_{32}^u \epsilon^2 tc^c) + H.C.
\end{eqnarray}
For FCNC Higgs couplings in
the charged lepton sectors, replace $h_{ij}^d \rightarrow h_{ij}^{\ell}$
in Eq. (9).
The factor in the $ds^c$ vertex is obtained as $7-5=2$, where $5$ is the
coefficient of the (2,2) entry in $Y_d$ and $7$ that of the (1,2) entry.

There is a tree--level contribution mediated by the Higgs boson for $K^0-\bar{K^0}$
mass difference.  We estimate it in the vacuum saturation approximation for the
hadronic matrix element \cite{soni}.
The new contribution, $\Delta m_K^{\rm Higgs}$, is given by
\begin{eqnarray}
\Delta m_K^{\rm Higgs} &\simeq& {4 \over 3} {f_K^2 m_K B_K \over m_{h^0}^2}\epsilon^{12}[
\{ {1 \over 6} {m_K^2 \over (m_d+m_s)^2} + {1 \over 6}\} {\rm Re}\left[\left({h_{12}^d+
h_{21}^{d*} \over \sqrt{2}}\right)^2\right] \nonumber \\
&-& \{ {11 \over 6} {m_K^2 \over (m_d+m_s)^2} + {1 \over 6}
\} {\rm Re}\left[\left({h_{21}^d-h_{12}^{d*} \over \sqrt{2}}\right)^2\right]]~.
\end{eqnarray}
Here $B_K$ is the bag parameter.
The first term arises from the scalar piece of the matrix element, while
the second term is from the pseudoscalar piece.  Note that the pseudoscalar
contribution is larger than the scalar one by about an order of magnitude.  On the
other hand, if the underlying theory is parity invariant, $h_{12}^d = h_{21}^{d*}$,
so only the scalar piece will be significant.  We shall also allow for this possibility.
Using $B_K = 0.75, f_K \simeq 160$ MeV, $\epsilon \simeq 1/6.5$ $m_s = 175$ MeV,
$m_d = 8.9$ MeV,
and with a parity non--invariant choice $h_{12}^d = 1,  h_{21}^d = 0.5$,
we obtain $\Delta m_K^{\rm Higgs} \simeq 3.1 \times 10^{-17}$ GeV, for $m_{h^0} = 100$ GeV.
This is two orders of magnitude below the experimental value.  If parity invariance is
assumed, with $h_{12}^d \simeq h_{21}^{d*} = 1$, the same choice of parameters
gives $\Delta m_K^{\rm Higgs} \simeq 6 \times 10^{-16}$ GeV, which is a factor of 6 below the
experimental value.  We see broad consistency with
data.

As for the CP violating parameter $\epsilon_K$, it receives a
new contribution from the Higgs exchange which can be significant and can even dominate over the
the SM CKM contribution.  For example, in the choice of parameters with
$|h_{12}^d| = 1, |h_{21}^d| = 0.5$, but with their phases being of order
one, $\epsilon_K$ arising from the Higgs exchange can explain the
observed value entirely.  This possibility will be tested at the $B$ factory.
New contributions to $\epsilon'$ are negligible.
However, since the prediction for $\epsilon_K$ is modified, the
standard model fit to the CKM parameter
Im($V_{td}^*V_{ts})$ will be modified by a factor of order one (depending on
the relative strength of the CKM and the Higgs exchange contribution to $\epsilon_K$).
Since $\epsilon'$ is directly proportional to Im($V_{td}^* V_{ts})$, its prediction
will be altered by a factor of order one \cite{buras}.

Electric dipole moments (EDM) of the neutron and the electron in our scheme is much
larger than the SM prediction.  The dominant source of the neutron EDM ($d_n$)
is from the two--loop
Barr-Zee \cite{barr} diagram involving the SM Higgs boson and the $Z$ boson.
In our scheme, $h^0$ has a scalar as well as a pseudoscalar coupling to the
$u$--quark.  In the physical basis of the $u$--quark the pseudoscalar coupling
has a strength of order Im($h_{11}^u) (9 \epsilon^6)$.  From this, we estimate the neutron
EDM to be in the range $(10^{-26}-10^{-27})$ ecm for the phase of order unity
and Higgs mass of order 100 GeV.  The EDM of the electron ($d_e$) will arise from an
analogous diagram, but its magnitude is suppressed by an additional power of
$\epsilon^2$.  We estimate $d_e \approx 10^{-27}-10^{-28})$ ecm.  Both $d_e$
and $d_n$ are within reach of future experiments.

All other constraints from low energy flavor changing processes are easily
satisfied in the model.  For example, the process $K_L \rightarrow \mu^+ \mu^-$
will proceed via tree--level Higgs exchange.  The four--Fermion interaction
for the decay is given by:
\begin{equation}
{\cal L}^{\rm eff} \simeq {5 \over 2} h_{22}^\ell (h_{21}^d-h_{12}^{d*})
{\epsilon^{10} \over m_{h^0}^2} (\bar{d}\gamma_5 s)(\bar{\mu}\mu)+ H.C.
\end{equation}
This amounts to an amplitude of strength $\sim 10^{-7}G_F$,
which implies that the branching ratio for
$K_L \rightarrow \mu^+\mu^-$ will be of order $10^{-14}$ from Higgs exchange.
This is well below the experimental limit.  Other processes such as
$K_L \rightarrow \mu e$, $K \rightarrow \pi \bar{\nu} \nu$, $\mu \rightarrow
e \gamma$, $\mu \rightarrow 3e$, $B_d-\bar{B_d}$ mixing, etc are all orders of magnitude
below the corresponding experimental limits.

$D^0-\bar{D^0}$ mixing, on the other hand, is predicted to be near the present
experimental limit in our scenario.  Note that the FCNC $uc^ch^0$ vertex is
enhanced by a factor of $\epsilon^2$ compared to $ds^ch^0$ vertex (see Eq. (9)).
The new contribution  $\Delta m_D^{\rm Higgs}$ is given by an expression
analogous to Eq. (10).  Using $f_D \simeq 200$ MeV, $B_D \simeq 0.75$,
$m_D^2/(m_c+m_u)^2 \simeq 2$, and $h_{12}^u =1, h_{12}^u = 0.5$, we estimate
$\Delta m_D^{\rm Higgs} \simeq 4.7 \times 10^{-14}$ GeV for a Higgs mass of 200
GeV, to be compared with the present experimental limit of $ \Delta m_D \le
1.6 \times 10^{-13}$
GeV \cite{pdg}.  Even if parity invariance is assumed, the new contribution does not
diminish for the $D^0$ system.  Allowing for reasonable order one uncertainties in the hadronic
matrix element and the FCNC couplings, and varying $\epsilon$ in the
range $1/6-1/7)$, we conclude that $D^0-\bar{D^0}$ mixing
should be not more than a factor of $(10-20)$ below the present limit, provided that
the Higgs boson mass is below about 300 GeV.  This prediction should be
testable in the near future.  We should remark that the SM long distance contribution
to $D^0-\bar{D^0}$ mixing \cite{pakvasa} is expected to be about three orders of magnitude below
the current limit, so its discovery should be clear--cut signal for new
short--distance physics, such as the one we propose.

\section{New signals for the Higgs boson}

There are a variety of new signals associated with the production and decay
of the Higgs boson in our scenario.  While the tree-level Higgs couplings
to the gauge bosons are identical to that of the SM, the Yukawa interactions
are modified from that of the SM.  The consequences are significant for the
strategy to discover the Higgs boson at colliders \cite{higgs}.  We list
below  the processes
where the differences from the SM are most striking.

1.  The Higgs boson couplings to light fermions are enhanced by
a flavor--dependent numerical factor.  At the NLC and perhaps the LHC, the Yukawa couplings
to $(b,\tau,t)$ will be measured \cite{dawson}.  The scenario can thereby be directly
tested.  The enhanced coupling to $\mu$ implies that
the Higgs production cross section will increase by a factor of 25 at
a muon collider \cite{muon}, relative to that of SM.

2.  The branching fraction
for Higgs decay into $\bar{b}b$ will increase by a factor of 9 relative to SM.
In the SM, the $\bar{b}b$ and the $W^*W^*$ branching ratios become
comparable for a Higgs mass of about 135 GeV \cite{spira}.  Because of the difference
in its $b$--quark coupling, this cross--over occurs in our scheme
for $m_{h^0} \simeq 155$ GeV.  This difference should be incorporated into the strategy
for discovering the Higgs boson at the upgraded Tevatron \cite{han}.

3.  The $h^0\gamma\gamma$ coupling, which arises primarily from the
$W^\pm$ loop, is unaffected in our scheme.  Since the $\bar{b}b h^0$
vertex is enhanced by a factor of $9$, the branching ratio for $h^0
\rightarrow \gamma \gamma$ will decrease by a factor of approximately 9, if
the Higgs boson mass is below about 155 GeV.

3.  At hadron colliders such as the upgraded Tevatron and the LHC,
gluon fusion is a dominant source for Higgs production.  Since the
top--quark coupling $\bar{t}th^0$ in our scheme is identical to
that of the SM, the production cross section via gluon fusion is
little affected.  However, there is some difference since the
$b$--quark contribution is not negligible. In the SM, for Higgs
mass of $(100,150,200)$ GeV, the $b$--quark loop contribution to
$gg \rightarrow h^0$ decreases the production cross section by
about $(9,5,3)\%$ (for $m_b(m_b) = 4.25$ GeV). In our scenario,
Higgs production production is decreased from the SM value by
$(17,10,6)\%$.  This should be measurable at
the LHC.

4.  There is an exciting new discovery channel for the Higgs boson at the
upgraded Tevatron.  This utilizes the flavor--changing $\bar{t}ch^0$ vertex.
Note that in our scenario this vertex has a strength of order $2\epsilon^2$,
which is similar in magnitude to that of the flavor conserving $\bar{b}bh^0$
coupling.

If $m_{h^0}$ is less than $m_t$, the decays $t \rightarrow ch^0$
and $\bar{t} \rightarrow \bar{c}h^0$ can provide a new channel for
Higgs discovery.  Once produced, $h^0$ will decay into $\bar{b}b$ with
a significant branching ratio.  The width for $t \rightarrow c h^0$ in our scheme is
given by:
\begin{equation}
\Gamma(t c h^0) \simeq {G_F m_t^3 \over 8\pi \sqrt{2}}\left(1-{m_{h^0}^2 \over
m_t^2}\right)^2\left[2 \epsilon^4\{\left|{h_{23}^u \over h_{33}^u}\right|^2+
\left|{h_{32}^u \over
h_{33}^u}\right|^2\}\right]~.
\end{equation}
For $m_{h^0}= 100$ GeV, and with $|h_{23}^u/h_{33}^u| \simeq |h_{32}^u/h_{33}^u| \simeq
1$, the branching ratio is $Br(t\rightarrow ch^0) \simeq 1.1 \times 10^{-3}$.
With an integrated luminosity of $20 {\rm fb}^{-1}$, about $1.4 \times 10^{5}$~
$t\bar{t}$
pairs are expected at Tevatron running at $\sqrt{s} = 2$ TeV.  This would lead
to about 300 Higgs events from $t$ and $\bar{t}$ decays.  The invariant mass of
the $b\bar{b}$ jets will provide the Higgs signal.  The QCD background can be
brought under control by tagging the $t$ (or $\bar{t})$ by its decay into $W +b$.
This process could be useful to discover a Higgs boson of mass as large as about
150 GeV at the Tevatron (at which point the kinematic suppression becomes significant).
This reach can be as large as 170 GeV at LHC.

For $m_{h^0}$ between $m_t$ and $2m_t$, there is another way to look for the Higgs boson
at hadron colliders.  The cross section for $h^0$ production via $gg$ fusion is
about 15 pb for $m_{h^0} = 200$ GeV at the LHC running at $\sqrt{s}=14$ TeV.
The Higgs will decay dominantly into
$W$ pairs, but the branching fraction into $\bar{t}c+\bar{c}t$ is not negligible.
For $h_{23}^u \simeq h_{32}^u \simeq 1$, the branching ratio, allowing for the
kinematic suppression, is found to be $Br(h^0\rightarrow
\bar{t}c+\bar{c}t) \simeq 1\times 10^{-3}$.  With $100 {\rm fb}^{-1}$
of data, this will result in 1500 $\bar{t}c + \bar{c} t$ events.
The signal will be thus a $b$--jet, a charm and a $W$.
It might be possible
to reconstruct the invariant mass of the Higgs boson by studying the leptonic
decay of the $W$ from the top quark.  (Although there is a neutrino involved, its
four--momentum can be reconstructed by measuring the charged lepton momentum, up to
a two--fold ambiguity.)  The background from SM single top production can be
substantially reduced by the invariant mass requirement.  Two jet $+W$ production
(where a jet is misidentified as a $b$) and $\bar{b}bW$ where one $b$ is missed
are other dominant backgrounds.  The presence  of a top quark in the signal but
not in these background events can be utilized to provide further cuts\footnote{We thank
John Conway for very helpful discussions on this signal.}.

5.  Since the scale of flavor physics is identified to be $M \approx
(1-2)$ TeV from fermion masses and mixings, new phenomena associated with
flavor physics will show up at experiments performed with energies
greater than $(1-2)$ TeV.  This will happen at the LHC.  One concrete example
is the unraveling of the effective vertices in Eq. (1).  For example, at
$s \sim M^2$, the process $p p \rightarrow \bar{b}b (3h^0)$ will proceed without
much suppression in the coupling.
The relevant dimensionful coupling
is $h_{33}^d/M^2 \simeq m_b/v^3$.  For $m_{h^0} \sim 100$ GeV, we estimate
the cross section for this process at LHC to be in the fb range.
The signature will be quiet dramatic: 8 $b$--jets with
3 pairs adding to to the same invariant Higgs mass.

Before concluding we wish to remark briefly on possible ways of inducing the
effective Yukawa Lagrangian on Eq. (1).  It is quite conceivable that it has
a dynamical origin \cite{nambu}.  In order to show that Eq. (1) can be derived from
perturbative physics at the TeV scale we have constructed models which
should serve as proof of existence.  Suppose there are vector--like
fermions at the TeV scale in doublet and singlet representations.  They could
have flavor--symmetric Yukawa couplings to the SM fermions and the SM Higgs
doublet.  The masses of these vector fermions break the flavor symmetry
softly (or spontaneously).  Once these vector fermions are integrated out,
the effective Lagrangian of Eq. (1) will result.  Since the extra matter is
vector--like under the SM, no significant departures from the SM predictions
are expected in the precision
electroweak observables.

In summary, we have suggested a novel scenario to explain the observed
hierarchy of the fermion masses and mixings.  By using successive higher dimensional
operators involving the SM Higgs doublet field we obtained a
good fit to the fermion masses and mixings without putting small numbers
into the theory.  We pointed out several new features associated with the
Higgs boson that are predicted by this scheme.  While low energy constraints
are well satisfied, there is one prediction in $D^0-\bar{D^0}$ mixing that
will test the model.  The flavor changing $\bar{t}ch^0$ coupling, as large in
strength as the flavor conserving $\bar{b}hh^0$ coupling, leads to new possibilities
to discover the Higgs boson at the upgraded Tevatron and the LHC.

We wish to thank Mike Berger, John Conway and Sally Dawson for discussions.
This work is supported in part by  the Department of Energy Grant No. FG03-98ER41076
and by funds provided by the Oklahoma State University.