\section{Introduction} \label{sec:introduction
Loop models are interesting examples of statistical models of
extended objects. They are related to
the O($n$) spin model \cite{Nienhuis:exact,DoMuNiSc:duality},
a surface growth model \cite{KoHe::four-coloring},
the self-avoiding walk \cite{DuSa:exact},
the protein folding problem \cite{ChDi:protein},
and so on. It includes the fully packed loop model \cite{BlNi:honeycomb}
and the Hamiltonian
cycle problem \cite{OrItDo:hamiltonian,Suzuki:regular,Higuchi:aspectratio}
as particular limits.

The partition function of an O($n$) loop model on a lattice
with $N$ sites at the inverse temperature $x$ is given by
\begin{equation}
  Z_{\text{loop}}(n,x^{-1}) = \sum_{c\in\mathcal{C}}
    x^{\mathcal{N}_\mathrm{S}(c)-N} n^{\mathcal{N}_\mathrm{L}(c)}
    \quad \quad(n,x\in\mathbb{R}).
  \label{loop-partition}
\end{equation}
The summation is taken over the set $\mathcal{C}$
of all the non-intersecting loop
configurations drawn along links of the lattice.
The number of loops and that of sites visited by them
are denoted by
$\mathcal{N}_\mathrm{L}(c)$ and $\mathcal{N}_\mathrm{S}(c)$, respectively.

One may hope to study the model \eqref{loop-partition} by the transfer
matrix approach. For $n\in\mathbb{Z}_+$,
this is done in a simple way; one introduces
link variables whose values are either occupied states with one of $n$
colors or an unoccupied state and lets them interact on sites.
A transfer matrix is written
as a product of vertex weights straightforwardly.

For  $n\not\in\mathbb{Z}_+$, however,
the partition sum \eqref{loop-partition} cannot be rewritten in terms
of local degrees of freedom such as link variables in a simple way.
It is not trivial to have a \emph{local} transfer matrix%
\footnote{
  The use of the connectivity basis is discussed in subsection
  \protect\ref{sec:connectivity}.
  }.
I say a transfer matrix is \emph{local} when its component is written
as a product of weights each of which
is determined by the local state configuration around a lattice site.

It is surprising that, in two dimensions,
$n\not\in\mathbb{Z}_+$ models
admit a mapping onto a state sum model with a local vertex weight and
thus have local transfer matrices
\cite{Baxter:coloring,Nienhuis:CoulombGas,BaBlNiYu:packedloop}.
In fact, by choosing $s\in\mathbb{C}$ satisfying $n=s+s^{-1}$,
$Z_{\text{loop}}(n,x^{-1})$ can be written as
\begin{equation}
\!\!\!\!\!\!
Z_{\text{loop}}(n,x^{-1})
=\sum_{c\in\mathcal{C}} x^{\mathcal{N}_\mathrm{S}(c)-N}
                              (s+s^{-1})^{\mathcal{N}_\mathrm{L}(c)}
= \sum_{c\in\overline{\mathcal{C}}}  x^{\mathcal{N}_\mathrm{S}(c)-N}
\prod_{L\in\overline{\mathcal{L}}(c)} s^{\pm1},
\label{directed}
\end{equation}
where $\overline{\mathcal{C}}$ is the set of loop configurations
with a direction associated with each loop.
The set $\overline{\mathcal{L}}(c)$ consists
of all the directed loops in a configuration $c$.
A loop with the (counter-)clockwise direction is given a weight
$s^{+1}$ $(s^{-1})$.
This weight can be realized by associating
$s^{+1/4}$($s^{-1/4}$) with each right(left)-turn site
and the model can be regarded as a state sum model with a local vertex
weight.
This trick has made the study of two-dimensional loop models very fruitful.

Physics of loops in three dimensions is very attractive.
It is  realistic in the context of condensed matter physics.
There has been a continuous suspicion that two-dimensional ones
have missed some important ingredient in real physics,
\textsl{e.g.} the protein folding problem.
Three-dimensional loops have also rich mathematical structures.
For instance, loops can be knotted or linked
in three dimensions\cite{Nechaev:entangled}.
It is noted that a number of attractive proposals
have been made to generalize the loop model to higher dimensions
\cite{WiKa:geometric,WiKa:coloredmanifolds}.

The analysis of loop models and their generalizations in higher
dimensions is, however,
extremely hard to perform.
Needless to say, the number of configurations increases considerably.
For fugacity $n\not\in\mathbb{Z}_+$, which includes the interesting
case of the self-avoiding walk ($n=0$),
no way of constructing local transfer matrices is known.
This is because specialties of two dimensions cannot be used
to simplify problems any more.
The mapping \eqref{directed} makes use of the fact that the
a directed loop in two dimensions turns around just once either clockwisely or
counter-clockwisely.
It appears that this kind of tricks  never works in
higher dimensions.

In this article, I propose a model
which generalizes \eqref{loop-partition} in a fashion
specific to three dimensions.
It is furnished with loop fugacity that depends on the global
three-dimensional shape of loops.
I show that, despite this generalization which makes
the model even more non-local,
a local transfer matrix for the system  can be constructed
for a number of choices of fugacity.
These choices include the ones that give zero or non-integer weight to
loops.

This article is organized as follows.
In section \ref{sec:fugacity}, I define a loop model in three
dimensions generalizing \eqref{loop-partition}.
Its local transfer matrix is constructed for a family of points in the
parameter space in section \ref{sec:transfer_matrix}.
In section \ref{sec:entropy},
this transfer matrix is numerically diagonalized to yield an estimate
of the site entropy in the fully packed limit $x^{-1}=0$.
In the last section \ref{sec:discussions}, I discuss my results and
their relation to combinatorial problems.
In an appendix, a technical issue on the block-diagonalization of the
transfer matrix is addressed.

\section{Generalized fugacity}\label{sec:fugacity
I define a statistical model of loops
on the three-dimensional simple cubic lattice
$\mathbb{Z}^3=
\{\sum_{i=1}^3 m_i \mathbf{e}_i\in\mathbb{R}^3|m_i\in\mathbb{Z}\}$ ,
$\mathbf{e}_i \cdot\mathbf{e}_j=\delta_{ij}$.
The partition function is given by
\begin{equation}
  Z_{\text{loop}}[n](x^{-1})
= \sum_{c\in\mathcal{C}}
x^{\mathcal{N}_\mathrm{S}(c)-N}
\prod_{L\in\mathcal{L}(c)} n(A(L)),
\label{rloop-partition}
\end{equation}
where $\mathcal{L}(c)$ is the set of loops in a configuration $c$.
The loop fugacity $n$ is now promoted to a function
which
depends on the shape of $L\in\mathcal{L}(c)$
through a quantity $A(L)\in\mathbb{R}$ defined below.

To define $A(L)$,
one begins with associating
a closed trajectory on the unit sphere with each loop $L$.
One picks a direction for $L$.
On every point $\mathbf{x}\in L \subset \mathbb{R}^3$ except
for sites where $L$ makes a turn,
there is a unit tangential vector $\mathbf{v}(\mathbf{x})$ to $L$;
it is either of $\pm \mathbf{e}_i, \;i=1,2,3$.
One may regard $\mathbf{v}(\mathbf{x})$
as a mapping from $L\setminus \text{(`turn-sites')}$
to the unit sphere $\mathrm{S}^2$.

As one walks along $L$,
$\mathbf{v}(\mathbf{x})$  jumps from a point to another on $\mathrm{S}^2$.
One can naturally interpolate these points to
define a continuous trajectory $\mathbf{v}: L\rightarrow \mathrm{S}^2$.
One has only to declare that $\mathbf{v}(\mathbf{x})$ moves
along the geodesic (of length $\tfrac12\pi$) on $\mathrm{S}^2$
at each turn-site.
This is equivalent with smoothing a loop in neighborhoods of turn-sites
keeping it within the plane (Fig.\ref{fig:smoothing}).
\begin{figure}[tb]
  \begin{center}
    \leavevmode
    \includegraphics[scale=0.5]{smoothing.eps}
    \caption{
      The definition of $A(L)$.
      (a) The original loop $L$ with a direction
      associated. (b) The loop smoothed at turn-sites.
      (c) Trajectory of the tangent vector
      on the unit sphere.
      An area of $A(L)=2\pi$ is enclosed.
      \label{fig:smoothing}
      }
  \end{center}
\end{figure}
Then one defines $A(L)$ to be
the oriented area encircled in the right of the trajectory
$\mathbf{v}(\mathbf{x})$.
On the lattice $\mathbb{Z}^3$, $A(L)$ takes values
\begin{equation}
  A(L)=\frac12m\pi \; (m\in\mathbb{Z}) \label{A-unit}.
\end{equation}
In two dimensions, the quantity $A(L)$ takes values $\pm2\pi$
and this signature corresponds to that of $s^{\pm1}$ in \eqref{directed}.
Therefore \eqref{rloop-partition} incorporates an essential
ingredient of three-dimensional loops and is regarded
as a natural generalization of \eqref{directed}.

As is evident from the above construction, there is certain ambiguity
for the value of $A(L)$.
First, because the trajectory is drawn on a closed surface of area $4\pi$,
$A(L)$ is well-defined up to $4\pi$.
Second, the signature of $A(L)$ is changed when the picked
direction of $L$ is reversed.
I require that  $n(\cdot)$ in \eqref{rloop-partition} absorbs this ambiguity.
Hence, it should satisfy
\begin{gather}
  n(A)=n(A+4\pi)     \label{A-modfourpi},\\
  n(-A)=n(A)         \label{A-signature}.
\end{gather}
Eqs. \eqref{A-unit}, \eqref{A-modfourpi},\eqref{A-signature}
imply that the fugacity function $n(\cdot)$
can be specified
by five parameters $n(A), A=0,\frac12\pi,\pi, \frac32\pi,2\pi$.

In spite of the above restriction on $n(A)$, the model \eqref{rloop-partition}
includes many interesting cases. Consider, for example, fugacity
\begin{equation}
  n(A)=n_0 \delta^{(4)}_0(A),   \label{delta-n}
\end{equation}
with $n_0\in\mathbb{R}$ and
\begin{equation}
  \delta^{(a)}_b(A)=
  \begin{cases}
    1 & \text{if\ }  \left(\frac{A}{\pi}-b\right)\equiv0 \bmod a,\\
    0 & \text{otherwise}.
  \end{cases}
  \label{periodic_delta}
\end{equation}
The sum in \eqref{rloop-partition} is then restricted
to configurations which consist only of loops
with the oriented area $A \equiv 0 \mod 4\pi$.
It should be interesting to compare the site entropy with that of the
model with  $n(A)=n_0$.
It is also tempting to ask whether such an additional constraint
changes the critical behavior or not.
The present case reminds one of
the fully packed loop model in two dimensions.
Its universality class differs from that of densely packed
loop phase when
the additional constraint that the loop length
must be even is imposed \cite{BaBlNiYu:packedloop,Jacobsen:universality}.

\section{Transfer matrices from local vertex weights}%
\label{sec:transfer_matrix}
In order to construct a local layer-to-layer transfer matrix for
the loop model \eqref{rloop-partition},
I define a vertex model and show that it is equivalent with
\eqref{rloop-partition}.

The local degree of freedom $z$ of the vertex model lives on each link
$\langle \mathbf{r}, \mathbf{r}\pm\mathbf{e}_i\rangle$,
$\mathbf{r}\in\mathbb{Z}^3$.
It takes one of three values $\leftarrow$, $\rightarrow$,  and  $-$(empty).
On each site,  six neighboring link variables interact
by the vertex weight $W$ defined in Fig. \ref{fig:vertex_weight},
where $s(\omega)$ is a function that satisfies
\begin{equation}
 s(A_1)\times s(A_2)=s(A_1+A_2)
 \label{s-representation}
\end{equation}
and is specified further below.
The partition function of the vertex model is
\begin{equation}
  Z_{\text{vertex}}[s](x^{-1})
=\sum_{z=\leftarrow, \rightarrow,-}
\;\;
\prod_{\mathbf{r}\in\mathbb{Z}^3}
W( \{z(\langle \mathbf{r}, \mathbf{r}\pm\mathbf{e}_i\rangle)\}).
\label{vertex-partition}
\end{equation}
Evidently, the partition function \eqref{vertex-partition} has a
local transfer matrix which is a product of $W$'s.
\begin{figure}[tb]
  \begin{center}
    \leavevmode
    \begin{center}
    \includegraphics[scale=0.5]{vertex_weight.eps}
    \end{center}
    \caption{Vertex weights $W$ for $Z_{\text{vertex}}[s](x^{-1})$.
      \label{fig:vertex_weight}
      }
  \end{center}
\end{figure}

Now I  show that \eqref{vertex-partition} for an appropriate $s$
is equivalent with \eqref{rloop-partition}.
Because the weight $W$ is nonzero only when there is one incoming and
one outgoing arrow, contribution to the partition sum
\eqref{vertex-partition} comes only from
the set $\overline{\mathcal{C}}$ of directed loop configurations
\begin{equation}
  Z_{\text{vertex}}[s](x^{-1})
  =\sum_{c\in\overline{\mathcal{C}}} x^{\mathcal{N}_\mathrm{S}(c)-N}
  \prod_{L\in\overline{\mathcal{L}}(c)}
\left[ \;\;
 \prod_{\mathbf{r}\in L \cap \mathbb{Z}^3} s(\omega(\mathbf{r}))
\right]
\label{vertex-dloop-partition}
\end{equation}
where $s(\omega(\mathbf{r}))$ is the weight for  the vertex
at $\mathbf{r}\in\mathbb{Z}^3$ in Fig.\ref{fig:vertex_weight}.

The factor in the square bracket in \eqref{vertex-dloop-partition}
is associated with a directed loop component $L$.
To evaluate this quantity, it is crucial to observe that
\begin{equation}
A(L)=\sum_{\mathbf{r}\in L\cap\mathbb{Z}^3} \omega(\mathbf{r}).
\label{area-local}
\end{equation}
Actually, the weight system in Fig.\ref{fig:vertex_weight} is designed
to have this property
in refs.\cite{AsItKo:randomwalk,AsItKo:diffusion,Itoi:smooth}
in the context of random walk with a spin factor %
\cite{Polyakov:transmutation,IsItMu:geometric}.
Combining the properties \eqref{area-local} and \eqref{s-representation},
one finds that the factor in the square bracket is simply $s(A(L))$.

Using the same trick as that used in \eqref{directed}, one can further
write the partition function \eqref{vertex-dloop-partition}
as a sum over undirected loop configurations
\begin{equation}
  Z_{\text{vertex}}[s](x^{-1})
=\sum_{c\in\mathcal{C}} x^{\mathcal{N}_\mathrm{S}(c)-N}
\prod_{L\in\mathcal{L}(c)} ( s(A(L))+s(-A(L))).
\end{equation}
Therefore,
if $n(A)$ in the loop model  \eqref{rloop-partition} can be written as
\begin{equation}
 n(A)=s(A)+s(-A)
 \label{fugacity-decomposition}
\end{equation}
with a function $s$ that satisfies \eqref{s-representation},
then
the vertex model partition function $Z_{\text{vertex}}[s](x^{-1})$
defined above is equal to \eqref{rloop-partition}.

The requirement \eqref{s-representation} together with
the restriction
\eqref{A-unit},
\eqref{A-signature} and  \eqref{A-modfourpi}
on $n(A)$ forces $s(A)$ to have a simple form
\begin{equation}
  s(A)= e^{iJ A}
\end{equation}
with $J\in\mathbb{Z}/2$.
It is enough \cite{AsItKo:randomwalk} to consider the cases $J=0,
\tfrac12, 1, \tfrac32$ and $2$ because of \eqref{A-unit}.
Hereafter, one introduces a shorthand notation
\begin{equation}
  Z_J(x^{-1}) =
Z_{\text{loop}}[ n(A)= e^{iJ A}+e^{-iJ A} ](x^{-1}).
\end{equation}

The vertex weights at $J=0$ and $2$ enjoy a special property
$s(\omega)=s(-\omega)$.
This enables one to define a model with
only two microscopic states $\leftrightarrow$ and $-$:
\begin{align}
  Z'_{\text{vertex}}[s](x^{-1})
=&\sum_{z=\leftrightarrow, -}
\prod_{\mathbf{r}\in\mathbb{Z}^3}
W( \{z(\langle \mathbf{r}, \mathbf{r}\pm\mathbf{e}_i\rangle)\})\nonumber\\
=&\sum_{c\in\mathcal{C}} x^{\mathcal{N}_\mathrm{S}(c)-N}
\prod_{L\in\mathcal{L}(c)} s(A(L)),
\end{align}
which I denote by $Z_{0'}$ and $Z_{2'}$.

The fugacity functions corresponding to
$J=0, \tfrac12, 1, \tfrac32,2,0'$ and $2'$ are listed in Table
\ref{tab:generators}.
They indeed give $n(A)\le0$ or $n(A)\not\in\mathbb{Q}$ for some loops.
\begin{table}[tb]
  \begin{center}
    \leavevmode
    \begin{tabular}{|l|l|}
      \hline
      $J$ &   $n(\cdot)$          \\
      \hline\hline
      ${0}$        &  $2\delta^{(\tfrac12)}_0(\equiv2)$   \\
      ${\tfrac12}$  &
$2(\delta^{(4)}_0-\delta^{(4)}_2)
+\sqrt{2}
(\delta^{(4)}_{\tfrac12}
+\delta^{(4)}_{\tfrac72}
-\delta^{(4)}_{\tfrac32}
-\delta^{(4)}_{\tfrac52})
 $\\
      $1$           & $2(\delta^{(2)}_0 - \delta^{(2)}_1)$ \\
      ${\tfrac32}$  &
$2(\delta^{(4)}_0-\delta^{(4)}_2)
-\sqrt{2}
(\delta^{(4)}_{\tfrac12}
+\delta^{(4)}_{\tfrac72}
-\delta^{(4)}_{\tfrac32}
-\delta^{(4)}_{\tfrac52})
 $ \\
      ${2}$        &  $2(\delta^{(1)}_0 - \delta^{(1)}_{\tfrac12})$  \\
\hline
      ${0'}$        &  $\delta^{(\tfrac12)}_0(\equiv1)$   \\
      ${2'}$        &  $\delta^{(1)}_0 - \delta^{(1)}_{\tfrac12}$ \\
\hline
    \end{tabular}
  \end{center}
  \caption{ The generalized fugacity $n$
    which generates the semigroup of allowed models.
    The function   $\delta^{(a)}_b(\cdot)$  is defined in
    eq.\protect\eqref{periodic_delta}.
    \label{tab:generators}
    }
\end{table}

Although only finite number of vertex models
$J=0,\tfrac12,1,\tfrac32,2,0'$ and $2'$ have been constructed  above,
it is possible  to construct an infinite number of ones by taking
the direct sum of the space of their microscopic states.
More precisely,
one generalizes the link variable to
take one of $2q+1$ $(q\in\mathbb{Z}_+)$ states:
$\leftarrow_k, \rightarrow_k$ with the $k$ labeling colors
$k=1,\ldots,q$ and an uncolored empty state $-$.
Introducing parameters
  $J_k\in\{0,\tfrac12,1,\tfrac32,2\}$,
the vertex weight assignments in Fig.\ref{fig:vertex_weight}
are supplemented by additional rules:
\begin{itemize}
\item If the both two arrows have the $k$-th color,
  then  $W=e^{i J_k \omega}$.
\item If the two colors do not agree, $W=0$.
\end{itemize}
The cases $J_k=0', 2'$ are handled in the obvious manner.

The fugacity of the
`direct sum' model $Z_{J_1\oplus J_2 \oplus \cdots \oplus J_q}$
is simply the sum:
\begin{gather}
n(A) = \sum_{k=1}^q \left[
 (e^{iJ_k A} + e^{-iJ_k A})\times
 B(J_k)\right],\\
B(J)=
\begin{cases}
  1       & ( J=0,\frac12,1,\frac32,2), \\
  \frac12 & ( J=0',2').
\end{cases}
\end{gather}
One immediately notices  that
\begin{equation}
  Z_{0'\oplus0'} = Z_{0}, \quad\quad  Z_{2'\oplus2'} = Z_{2}.
\end{equation}
Thus the fugacity functions expressible via  vertex models
form an infinite semigroup under addition
\footnote{
  This direct sum operation may be used for the lattice
  construction
  \protect
  \cite{AsItKo:randomwalk,AsItKo:diffusion,Itoi:smooth}
  of higher-spin three-dimensional field theories.
}
generated by
$J=0',\tfrac12,1,\tfrac32$ and $2'$.

One can also take the `tensor product' of the space of microscopic states of
$Z_{J_1}$  and $Z_{J_2}$ to define a model $Z_{J_1\otimes J_2}$.
Let the link variable take five values
$(z_1,z_2)=
  \uparrow\uparrow,   \uparrow\downarrow,
  \downarrow\uparrow, \downarrow\downarrow,$ and $||$.
The vertex weight $W$ is defined to be the product of $W$'s with
$J=J_1$ and $J_2$.
Then the loop fugacity becomes
\begin{equation}
n(A) =
\left[
 (e^{iJ_1 A} + e^{-iJ_1 A})\times
 B(J_1)\right] \times
\left[
 (e^{iJ_2 A} + e^{-iJ_2 A})\times
 B(J_2)\right].
\end{equation}
However, a new fugacity function cannot be realized because
the partition function $Z_{J_1\otimes J_2}$ is equivalent with
an appropriate direct sum
\begin{equation}
  Z_{J_1\otimes J_2} = Z_{\overline{J}_1\oplus\cdots\oplus\overline{J}_q}
\end{equation}
corresponding to the decomposition rule of the representation of SU(2).

\section{Entropy estimates}\label{sec:entropy
I numerically diagonalize the transfer matrices constructed in section
\ref{sec:transfer_matrix}.
Throughout this section, I concentrate on the fully packed limit $x^{-1}=0$
where all the sites are visited by a loop.
This simple case is in fact very interesting case;
in two dimensions, this limit yields a new universality class with a
shifted central charge
on several bipartite lattices and has been attracting much attention
\cite{Jacobsen:universality,%
  Higuchi:decorated,%
  DiGuKr:fully_packed_eulerian,%
  GuKrNi:hamiltonian_eulerian}.
It would be interesting to look at the limit
where the two strong constraints are combined:
the fully-packing constraint $x^{-1}=0$ and the constraint
\eqref{delta-n} on the shape of loops.

The site entropy in the thermodynamic limit is defined by
\begin{equation}
  f[n](\infty)
= \lim_{N\rightarrow\infty} \frac{1}{N} \log Z_{\text{loop}}[n](x^{-1}=0).
\label{thermodynamic-entropy}
\end{equation}
I evaluate this quantity on
quasi-one-dimensional geometry $L_1\times L_2 \times L_3,
L_3\rightarrow\infty $ while $L_1,L_2$ kept finite
by  calculating the largest eigenvalue of the transfer
matrix in an appropriate sector.

Let $T$ be the layer-to-layer transfer matrix in
$+\mathbf{e}_3$ (vertical)
direction for the vertex model defined in \eqref{vertex-partition}.
Then $T$ acts on linear combinations of
arrays of $L_1\times L_2$ vertical (colored) arrows.
One can take either hard-wall or periodic boundary condition
in the horizontal directions.

It is important to note that the transfer matrix $T$ commutes
with the operator giving the net flow of arrows of $k$-th color in $+\mathbf{e}_3$ direction
\begin{equation}
 d_k = (\#\uparrow_k) - (\#\downarrow_k), \label{baxter}
\end{equation}
which is understood as
\begin{equation}
   d_k = (\#\updownarrow_k) \bmod 2
\end{equation}
for $J_k=0'$ and $2'$.
Thus  $T$ is block-diagonalized as
\begin{equation}
  T=\bigoplus_{\mathbf{d}} T_{\mathbf{d}},
  \quad \mathbf{d}=(d_1,\ldots,d_q).
\label{block-decomposition}
\end{equation}

The quantity \eqref{thermodynamic-entropy} is obtained as
\begin{equation}
  f[n](\infty) =
\lim_{L_1, L_2\rightarrow\infty}
\frac{1}{L_1L_2} \log |\lambda_{\mathbf{0}}^{0}(L_1,L_2)|,
\end{equation}
where $\lambda^i_{\mathbf{d}}(L_1,L_2)$ is the $i$-th largest
eigenvalue of $T_{\mathbf{d}}(L_1,L_2)$.
The condition $\mathbf{d}=\mathbf{0}$ excludes unwanted configurations
that have unbalanced arrows traveling along the infinite direction.
Shown in Tables \ref{tab:irreducible} and \ref{tab:reducible}
are the finite-$L_1,L_2$ results%
\footnote{I have also measured several leading  eigenvalues of
  $T_\mathbf{d}$ with $\mathbf{d}=(d_1,\ldots,d_q)$, $d_k=0,\pm1$.
  These are related to correlation length of operators in the
  theory.
  These results will be reported elsewhere.
}.
The asymmetric Lanczos algorithm is utilized for the present sparse
eigenproblem.
\begin{table}[tb]
  \begin{center}
    \leavevmode
    \begin{tabular}{|l|l|l|l|l|l|l|}
        \hline
$J$
      & $2\times2$(h) & $3\times3$(h) & $3\times4$(h)
      & $2\times2$(p) & $3\times3$(p) & $3\times4$(p) \\
      \hline\hline
      ${0}$        & 0.54202495 & 0.59145447 & 0.63524092
                   & 1.0585126 & 0.83841678 & 0.83340128  \\
      ${\tfrac12}$ & 0.27123680 & 0.33576248 & 0.35050951
                   & 0.79451346 & 0.55900063 & 0.55895924  \\
      $1$           & 0.51585927 & 0.50234791 & 0.55646223
                    & 0.97170402 & 0.69812631 & 0.69832061 \\
      ${\tfrac32}$  & 0.35592318 & 0.35908194 & 0.37136801
                    & 0.79451346 & 0.57435935 & 0.49496387 \\
      ${2}$        & 0.49499647 & 0.45468972 & 0.51615498
                   & 1.0406166 & 0.66200716 &  0.67440981 \\
\hline
      ${0'}$        & 0.46298939 & 0.55650697 & 0.60072954
                    & 0.91847381 & 0.79631788 & 0.80135760 \\
      ${2'}$        & 0.38697370 & 0.37695844 & 0.42272584
                    & 0.87898824 & 0.55608904 & 0.58497412 \\
\hline
\hline
     & $4\times4$(p) &\multicolumn{5}{|l}{}\\
\cline{1-2}\cline{1-2}
$0'$ & 0.81947983 &\multicolumn{5}{|l}{}\\
$2'$ & 0.60931946 &\multicolumn{5}{|l}{}\\
\cline{1-2}
    \end{tabular}
  \end{center}
  \caption{
    The site entropy estimated numerically.
    $L_1\times L_2$ is the size of a layer while (p) and (h) mean
    periodic and hard-wall boundary conditions in a layer.
    \label{tab:irreducible}
    }
\end{table}

\begin{table}[tb]
  \begin{center}
    \begin{tabular}{|l|l|l|l|l|l|} \hline
    $\bigoplus_k J_k$  & $n$ & $2\times2$(h) & $3\times3$(h)
                             & $2\times2$(p) & $3\times3$(p) \\
\hline\hline
      ${0'\oplus2'}$ & $2\delta^{(1)}_0$ & 0.52330515 & 0.57390934
      &  1.0502400 & 0.81503626 \\
      ${0\oplus0}$ & 4 &  0.64498133 & 0.65020710
                       & 1.2354255 & 0.90747958 \\
      ${0\oplus2}$ & $4\delta^{(1)}_0$  & 0.63331428 & 0.62765833
                   & 1.232054 & 0.88008432 \\
      ${0'\oplus2'\oplus1}$ & $4 \delta^{(2)}_0$ & 0.63216104 & 0.61664018
                            & 1.2110124 & 0.86121744 \\
      ${0\oplus0\oplus0\oplus0}$ & $8$ & 0.76911076 &
                                 & 1.4501153 & 1.0108338   \\
      ${0\oplus0\oplus2\oplus2}$ & $8\delta^{(1)}_0$ & 1.3280134 &
                                 &  1.4490048 & 0.98982172  \\
      ${0\oplus1\oplus1\oplus2}$ & $8\delta^{(2)}_0$& 0.75985245 &
                                 & 1.4364447  & 0.97947213  \\
      ${0'\oplus2'\oplus\frac12\oplus1\oplus\frac32}$ & $8\delta^{(4)}_0$ & 0.36464934 &
            & 1.3863922  &  0.88604978 \\
      ${\frac12\oplus\frac32}$ &$4(\delta^{(4)}_0-\delta^{(4)}_2)$
& 0.47382047 & 0.41770663
& 1.0397208 &    0.63630800 \\
      ${0\oplus\frac12}$ & $\not\in\mathbb{Q}$  & 0.40013109 & 0.54670471
             & 1.1686609 & 0.84652598 \\
      ${0\oplus\frac32}$ & $\not\in\mathbb{Q}$ & 0.44457252 & 0.55636026
 & 1.1686609 & 0.84694306 \\
      ${0'\oplus0\oplus\frac12}$ & $\not\in\mathbb{Q}$
    & 0.48094569 & 0.57807813
    & 1.2520659 & 0.62795992  \\
     \hline
\hline
\cline{1-4}
        &       & $3\times4$(h) &$3\times4$(p) & \multicolumn{2}{|l}{}\\
\cline{1-4}
$0'\oplus2'$ & $2\delta^{(1)}_0$
& 0.61935581 & 0.81580656 &\multicolumn{2}{|l}{}\\
\cline{1-4}
    \end{tabular}
  \end{center}

  \caption{
    The site entropy estimated numerically.
    $L_1\times L_2$ is the size of a layer while (p) and (h) mean
    periodic and hard-wall boundary conditions in a layer.
    \label{tab:reducible}
    }
\end{table}
The obstacle in making $L_1$ and $L_2$ large in the actual
numerical work is of course the exponential growth of the
dimensionality of the transfer matrix.
The selection of $\mathbf{d}=\mathbf{0}$ sector helps
to reduce the dimensionality, though the improvement is polynomial.
For a direct sum $Z_{J_1 \oplus \cdots \oplus J_q}$,
the dimensionality is estimated to be
\begin{equation}
\dim (\mathbf{d}=\mathbf{0} \text{ sector})
\sim
(1 + 2u + p )^{L_1L_2}\times 2^{-p}
\left(\frac{3}{4\pi L_1L_2}\right)^{u/2},
\end{equation}
where $p=\#\{k| J_k=0' \text{\ or } J_k=2'\}$, $u=q-p$.
One notices that it is better to recast $Z_{2'\oplus 2'}$ as $Z_2$
if one is interested in its $\mathbf{d}=\mathbf{0}$ sector.

The exponential growth is severe even after
restricting to  the $\mathbf{d}=\mathbf{0}$ sector.
In order to increase $L_1L_2$
as much as possible within the available computer resources,
I have further decomposed $T_\mathbf{0}$ with respect to the
eigenvalue of shift (lattice momentum) operator
for the periodic boundary case where the translational symmetry is present.
I have looked at the zero-momentum sector ${T_\mathbf{0}}^{(0,0)}$
as described in the appendix \ref{sec:zero_momentum}.
It is quite natural to expect that the largest eigenvalue lies there.
By this decomposition, the dimensionality of the eigenproblem
is reduced, at most,  by $(L_1 L_2)^{-1}$.
As a drawback,
the matrix ${T_\mathbf{0}}^{(0,0)}$ becomes less sparse than
the original $T_\mathbf{0}$.
With both effects combined, some improvements in the memory
usage and the CPU time are observed.
Thus the analysis of larger systems becomes possible for the  periodic
boundary case, as seen in Tables \ref{tab:irreducible} and \ref{tab:reducible}.

\section{Discussions} \label{sec:discussions
In the comparison between the periodic and the hard-wall boundary
conditions,
one notices that the periodic case has always larger site entropy.
This is because many of the loops that wind  non-trivially in the
horizontal directions satisfy $A(L)=0$ and the loops with $A(L)=0$
contribute to every partition sum with a positive fugacity.

The numerical works in the present study have been carried out on
modest workstations.  Unfortunately, information in the thermodynamic
limit $L_1,L_2\rightarrow\infty$ is out of reach in the present
analysis.
For the study of criticality, ref.\cite{AsItKo:randomwalk},
where random walks with the weight in Fig.\ref{fig:vertex_weight} are studied,
is quite suggestive. It is reported that
Euclidean symmetry is not always recovered even in the continuum limit.

I discuss relations with combinatorial problems below.

\subsection{Even and odd number of loops of a specific type}
For most allowed values of $J$,
the loop fugacity takes both positive and negative values.
Some interesting combinatorial information is encoded in
these models.
For example, the linear combinations
$\tfrac12( Z_{0'} \pm Z_{2'} )$ counts the number
of loop configurations
such that there are
even (odd) number of loops for which  $\frac{2}{\pi}A(L)\equiv1\mod2$,
\textsl{e.g.},
\begin{equation}
  \frac12( Z_{0'} \pm Z_{2'} ) =
  \begin{cases}
   \sum_{c\in\mathcal{C}_{\text{even},\frac12}} 1,\\
   \sum_{c\in\mathcal{C}_{\text{odd},\frac12}} 1,
  \end{cases}
  \label{evenodd}
\end{equation}
where
$ \mathcal{C}_{\text{even(odd)},a}$ is the subset of $\mathcal{C}$
and is defined by the following properties.
\begin{itemize}
\item $c\in \mathcal{C}_{\text{even(odd)},a}$  contains even (odd)
  number of loops with $A(L)/a\pi \equiv 1 \mod 2$.
\item All other loops in $c\in \mathcal{C}_{\text{even(odd)},a}$
  satisfy  $A(L)/a\pi \equiv 0 \mod 2$.
\end{itemize}
Similarly,
the quantities
$  \frac12( Z_{0'\oplus2'} \pm Z_{1} )$ and
$  \frac12({Z_{0'\oplus2'\oplus1} \pm Z_{\frac12\oplus\frac32}})$
are interpreted as the sums over
$\mathcal{C}_{\text{even(odd)},1}$
and
$\mathcal{C}_{\text{even(odd)},2}$.

Of the two $Z$'s in \eqref{evenodd},
the one with the larger leading eigenvalue
dominates the sum in the limit $L_3\rightarrow\infty$
studied in section \ref{sec:transfer_matrix}.
In finite geometries, both terms contribute to yield an exact number.

\subsection{Self-avoiding walk}
The partition function \eqref{loop-partition} in
the limit $n\rightarrow0$
corresponds to the enumeration of self-avoiding walks.
Self-avoiding walks in three-dimensions have mainly been studied by
the exact enumeration method due to the lack of transfer matrix
formalism as pointed out in the introduction.

The model I propose in this paper can be regarded
as a step forward to overcome this difficulty;
in the models $Z_{0'\oplus2'}, Z_{0'\oplus2'\oplus1}$ and
$Z_{0'\oplus2'\oplus\frac12\oplus1\oplus\frac32}$,
the fugacity is set zero for families of loops.
This is, however, achieved by paying the cost of
having larger $n$ for another family of loops.
Within the present construction, loops with $A(L)\equiv0\mod 4\pi$
cannot have weights different from the number of possible link states.
Thus the partition function listed above serve only as a very loose
upper bound for the entropy of self-avoiding walks.

The problem of construction of a local transfer matrix to enumerate
self-avoiding walks on three-dimensional lattices
still remains open.

\subsection{Mapping to ribbon configurations}
The oriented area defined in \eqref{rloop-partition}
has a nice geometric interpretation as holonomy.
The tangent vector $\mathbf{v}(\mathbf{x})$
moves along a trajectory on $\mathrm{S}^2$.
Let the unit vector tangential to
this trajectory at $\mathbf{v}(\mathbf{x_0})$ be $\mathbf{u_0}$
(Fig.\ref{fig:holonomy}).
Consider the parallel transport (in the sense of riemannian
geometry) of $\mathbf{u}_0$
along the trajectory $\mathbf{v}(\mathbf{x})$ on $\mathrm{S}^2$.

When $\mathbf{u}_0$ is transported back to
$\mathbf{v}(\mathbf{x}_0)$, it gains some holonomy
( the angle $\theta$ in Fig.\ref{fig:holonomy}).
This holonomy angle is given by the integration of the scalar curvature of
$\mathrm{S}^2$ over the domain encircled by the trajectory and is
nothing but the oriented area $A(L)$.
\begin{figure}[tb]
  \begin{center}
    \leavevmode
    \begin{center}
    \includegraphics[scale=0.8]{holonomy.eps}
    \end{center}
    \caption{When $\mathbf{u}_0$, the tangent vector at
      $\mathbf{v}(\mathbf{x}_0)$, is parallel-transported
      along the trajectory on $\mathrm{S}^2$,
      it receives holonomy whose magnitude (angle $\theta$) is equal to the
      oriented area $A(L)$ the path encloses.
      \label{fig:holonomy}}
  \end{center}
\end{figure}
In the real space $\mathbb{Z}^3\subset\mathbb{R}^3$,
the holonomy described above is nicely kept track of
by broadening the loop segment to a `ribbon' with the
distinction of the right and the reversed sides.
The parallel transportation can be recasted
as a rule of bending ribbons on sites, which is
shown  in Fig.\ref{fig:ribbon}.
\begin{figure}[tb]
  \begin{center}
    \leavevmode
    \begin{center}
    \includegraphics[scale=0.8]{ribbon.eps}
    \end{center}
    \caption{
      Examples of allowed ways of bending ribbons at sites
      (except for the ones indicated by gray arrows)
      The holonomy is accumulated at the sites indicated  by gray arrows.
      The right and the reverse sides of ribbons are represented by
      white and black colors.
      \label{fig:ribbon}
      }
  \end{center}
\end{figure}
Among the loop configurations shown in Fig.\ref{fig:ribbon},
the partition sum $Z_{\text{loop}}$ with fugacity
$n=\delta^{(2)}_0$, $\delta^{(1)}_0$ and $\delta^{(\frac12)}_0$
receives contribution from
$\{\text{(a)}\}$,
$\{\text{(a),(b)}\}$,
and
$\{\text{(a),(b),(c)}\}$, respectively.

For $n=4\delta^{(2)}_0$, the sum in \eqref{rloop-partition} is
over ribbon loop configurations without mismatch.
In that interpretation,  the coefficient four is naturally regardes
as the number of directions the right side can face.
Therefore, $Z_{\text{loop}}[n=4\delta^{(2)}_0](x^{-1})$ is nothing but the
generating function of the number of allowed ribbon configurations%
\footnote{%
  This observation clearly indicate a simple way of
  constructing transfer matrices for this  fugacity function.
  The link variable represents the ribbon which can be placed in four ways.
  The rule of bending is implemented in the vertex weight.
  The construction in section \ref{sec:transfer_matrix}, however, has
  an advantage; the size of the matrix can be reduced more by the choice of the
  sector $\mathbf{d}=\mathbf{0}$.
}.

Similarly, the partition sum for $n=2\delta^{(1)}_0$ can be
interpreted as the sum over the configurations of ribbons without the
distinction between the right and the reverse sides,
while in the case $n=\delta^{(\tfrac12)}_0$, the loop segment is just a
chord.
This interpretation suggests that
\eqref{rloop-partition} may be regarded as a model of polymers with
various partially broken axial symmetry by, for instance, the presence of
side chains.

\subsection{Comparison with the connectivity basis}
\label{sec:connectivity}
The connectivity basis\cite{ScHiKl:compact,BlNi:square} is very powerful
in that one can always write a transfer matrix for a loop model with
respect to it.
It has been very useful for numerical calculation in two-dimensions.

Nevertheless, I have avoided the use of the connectivity basis
in the present work. The reason is the following.
First, its fundamental degrees of freedom are not the link variables
and the transfer matrix with respect to it is not local.
Local transfer matrices have  merits even in two dimensions.
Namely, it paved the way to the Bethe ansatz
solution\cite{Baxter:exactly,BaBl:honeycomb}
and the conformal field theoretic description
\cite{Cardy:geometrical,JaKo:fieldtheory}
via Coulomb gas representation.
Second, in three dimensions and higher, the size of connectivity
basis grows considerably because of the lack of the planarity constraint.
It is not clear if it is effective to perform numerical calculation in this
basis.
In two dimensions, the present basis is as good as the connectivity one
\cite{BaBlNiYu:packedloop,Higuchi:decorated}.

I suppose it is very important to see how useful the connectivity basis in
three dimensions is and to try to improve the efficiency of the
calculation in that basis.

\section*{Acknowledgments}
\pdfbookmark{Acknowledgements}
I thank S. Tanimura for useful discussions on geometric phases.
I gratefully acknowledge useful conversations with
M.~Asano, E.~Guitter, C.~Itoi,  S.~Hikami, K.~Minakuchi, and J.~Suzuki.
I thank T.~Iwai and Y.~Uwano for discussions and for hospitality at
Kyoto University, where a part of this work was done.

This work was supported by the Ministry of Education, Science
and Culture under Grant 08454106 and
\href{http://rice.c.u-tokyo.ac.jp/~hig/kaken.html}{10740108}
and by \href{http://www.jst.go.jp}{Japan Science and Technology Corporation}
under
\href{http://www2.jst.go.jp/crest/}{CREST}.  \medskip


\section{Introduction}

The Standard Model (SM) of strong and electroweak interactions possesses only
one source of CP violating phases: the phase in the CKM matrix. While the
physical origin of this phase remains unknown, it appears to satisfactorily
account for the observed CP violation in the K meson system, and future data
from B factories and high energy accelerators will be able to test its validity
for B meson phenomena. It was early on realized that supersymmetric (SUSY)
extensions of the Standard Model allowed for an array of new CP violating
phases, and that these phases could give large contributions to the electric
dipole moments [EDMs] of the electron and neutron \cite{redm1}, thus violating
the known experimental bounds
\cite{de,dn}. Several resolutions to this problem have been proposed,  e.g. one
might assume the phases are quite small i.e.
$\stackrel{<}{\sim}O(10^{-2})$\cite{redm1,redm3}, or suppress the diagrams by
assuming the relevant SUSY particles are very heavy \cite{redm4}. Both these
suggestions a priori appear unsatisfactory, i.e. the first possibility would
appear to require a significant amount of fine tuning, and the second would
imply a SUSY spectrum in the TeV domain, possibly beyond the reach of even the
LHC. More recently, it has been realized that some cancellations can occur in
different contributions to the EDMs in certain parts of the SUSY parameter
space \cite{nath1}
which has led to considerable examination of this possibility
\cite{nath1,nath2,nath3,falk1,bk,falk2,kane1,kane2,bartl,pokorski}. In the
minimal supersymmetric extension of the SM, the MSSM, (which we define as the
low energy SUSY extension of the SM obtained by supersymmetrizing the particle
spectrum with two Higgs doublets) this does appear to reduce the amount of fine
tuning needed, allowing for larger phases. The existence of such phases could
also have effects on other predictions of the MSSM
\cite{demir,carm,falk3,falk4,nath4,pila},
thus allowing future experimental checks of this idea.

In this paper we consider these questions from the viewpoint of supergravity
(SUGRA) grand unified models (GUTs). In particular we consider gravity mediated
models with R-parity invariance, where supersymmetry is broken in a hidden
sector at some scale above the GUT scale $M_G$ (presumably at the string or
Planck scale $M_{Pl}$), the breaking being transmitted to the physical sector by
gravity \cite{acn}. (For previous work from this viewpoint see
Refs.\cite{nath1,nath2,falk1,bk,kane2,bartl}.) Such models are much
more constrained than the MSSM. Thus low energy properties are now determined by
running the renormalization group equations (RGE) from $M_G$ to the electroweak
scale. CP violating phases that in the MSSM are arbitrary parameters now get
correlated by the RGE, and no longer can be assigned independently to satisfy
the EDM experimental bounds. Further, the radiative breaking of $SU(2)\times
U(1)$ puts additional constraints on the CP violating phases. We find that as a
consequence of these constraints, one must still fine tune some of the CP
violating phases at the electroweak scale to be quite small, except for small
regions of the parameter space with tan$\beta$ close to its minimum value, or
when SUSY masses are very  large. Further, except for the lowest tan$\beta$ a
serious fine tuning develops at the GUT scale. These phenomena appear to be true
both for the minimal mSUGRA model with universal soft breaking masses, and with
nonuniversal soft breaking extensions. Hence the required smallness of the
phases also tends to reduce the possibility of collateral evidence in other
phenomena for their existence in SUGRA models.

In Sec.2, we discuss the solutions of the RGE for mSUGRA models and explain
which phases that are large at the GUT scale get naturally suppressed at the
electroweak scale (and which do not), and what the constraints of electroweak
breaking imply both at the tree and 1-loop level. We then examine the amount of
fine tuning needed to satisfy the existing EDM data, and what would be needed if
those data were improved by a factor of 10 without finding any EDM. Sec.3
considers some effects of nonuniversal soft breaking on the EDM results.

In gravity mediated SUGRA models, the structure of the soft breaking parameters
is to be deduced from the nature of the Kahler potential and gauge kinetic
function at the GUT scale \cite{acn,cfgv,wb}. Thus what might be reasonable
sizes for CP violating phases are presumably Planck scale physics questions.
While at present there is no theory of such phenomena, one can still examine the
general framework to see under what circumstances small or large phases might
occur ``naturally'' at the GUT scale. This is analyzed in Sec.4, and a model
with ``naturally'' small CP violating phases is discussed.

Concluding remarks are given in Sec.5.

\section{EDM For mSUGRA Models} Supergravity GUT models with universal soft
breaking of supersymmetry, mSUGRA, depend upon five  parameters at the GUT
scale: $m_{1/2}$ (the universal gaugino mass), $A_0$ (the cubic soft breaking
mass, $B_0$ (the quadratic soft breaking mass), $\mu_0$ (the Higgs mixing
parameter) and $m_0$ (the universal squark and slepton mass). Of these, the
first four may be complex. However, it is always possible to make a phase
transformation on the gaugino fields to make $m_{1/2}$ real, and since the
reality is preserved by the RGE at 1-loop order, we will assume for now on that
$m_{1/2}$ is real. We parameterize the remaining phases at $M_G$ by :
\begin{eqnarray} A_0&=&|A_0|e^{i \alpha_{0A}};\, B_0=|B_0|e^{i
\theta_{0B}};\,\mu_0=|\mu_0|e^{i
\theta_{0\mu}}.
\label{amb}
\end{eqnarray} The RGE allow one to determine the low energy parameters in terms
of those of Eq.(\ref{amb}) and $m_{1/2}$ and $m_0$. Thus at the electroweak
scale, one obtains a different A parameter for each fermion ($A_u$, $A_d$,
$A_t$, $A_b$,
$A_e$ and
$A_\tau$), real gaugino masses $\tilde m_3$, $\tilde m_2$ and $\tilde m_1$ for
the $SU(3)_C$, $SU(2)_L$ and $U(1)_Y$ sectors, and complex $B$ and $\mu$
parameters.

The EDM, $d_f$ for fermion f, appears in the effective Lagrangian $L_f$ as
\begin{eqnarray}  L_f&=&-{i\over 2}d_f\bar f\sigma_{\mu\nu}\gamma^5 f F^{\mu\nu}
\label{lag}
\end{eqnarray} The basic diagrams giving rise to $d_f$ are shown in
Fig.\ref{fig1}, and involve neutralino ($\tilde\chi^0_i$, i=1,4), chargino
($\tilde\chi^{\pm}_i$, i=1,2) and gluino ($\tilde g$) loops with squarks ($\tilde
q$) and sleptons ($\tilde l$).

For the neutron dipole moment $d_n$ one must also take into account the gluonic
operators
\begin{eqnarray}  L^G&=&-{1\over
3}d^Gf_{abc}G_{a\mu\alpha}G_{b\nu}^{\alpha}\tilde{G}_c^{\mu\nu}
\label{lg}
\end{eqnarray} and
\begin{eqnarray}  L^C&=&-{i\over 2}d^C\bar
q\sigma_{\mu\nu}\gamma^5T^aqG_a^{\mu\nu}
\label{lc}
\end{eqnarray} where $\tilde{G}^{\mu\nu}_c= {1\over
2}\epsilon^{\mu\nu\alpha\beta}G_{c\alpha\beta}$, $f_{abc}$ are the SU(3)
structure constants and $T^a={1\over 2}\lambda^a$, where $\lambda^a$ are the
SU(3) Gell-Mann matrices. In addition to the diagrams of Fig.\ref{fig1}
with $\gamma$ replaced by $g$,
operators of Eqs.(\ref{lag}-\ref{lc}) receive contributions from the two loop
Barr-Zee diagrams of Fig.\ref{fig2}\cite{ckp}, and the two loop Weinberg type
diagram of Fig.\ref{fig3} \cite{dai}. For the neutron dipole moment one must use
the QCD factors
$\eta^{ED}$, $\eta^G$, $\eta^C$ to evolve the results at the electroweak scale
down to 1 GeV \cite{aln}.

The calculation of $d_n$ suffers from QCD uncertainties from several sources.
These include:
\par\noindent (1) How to relate the quark moments $d_u$, $d_d$ to the neutron
moment
$d_n$. We use here the non relativistic quark model relation. (For other
approaches see Ref.\cite{bartl}.)
\begin{eqnarray}  d_n&=&{1\over 3}(4d_d-d_u)
\label{qm}
\end{eqnarray} (2) How to relate color and gluonic contributions to the electric
dipole moment. We use here the ``naive dimensional analysis'' of Ref.\cite{mon}.
\par\noindent (3) Uncertainty in $m_s$, the strange quark mass, which affects the
determination of $m_u$ and $m_d$. While the quark mass ratios are well known
\cite{lr},
\begin{eqnarray}  {m_u\over m_d}&=&0.553\pm 0.043;\,\,\,\, {m_s\over
m_d}=18.9\pm 0.8
\label{qr}
\end{eqnarray} one has $m_s=(175\pm 25)$MeV from QCD sum rules and $m_s=(100\pm
20 \pm 10)$ MeV  (2 GeV) from the quenched lattice calculation\cite{ms1} (Using
unquenched  approximation one expects an even smaller value \cite{ms2}).

Thus the values of $d_n$ calculated below have a significant uncertainty. In
spite of this, it is useful to see the effect on the SUSY parameter space of
simultaneously imposing the experimental constraints on both
$d_n$ and $d_e$ and we will discuss below how the uncertainty in $m_s$ effects
these constraints.

To specify our phase conventions, we give the mass matrices for the particles
entering in Fig.\ref{fig1}. With a convenient choice of phases, the chargino and
the neutralino mass matrices are:
\begin{eqnarray} M_{\chi^{\pm}}&=&\left(\matrix{
 \tilde m_2                & \sqrt 2 M_W sin\beta  \cr
  \sqrt 2 M_W cos\beta            &-|\mu|e^{i\theta} }\right)\label{char}
\end{eqnarray}

\begin{eqnarray} M_{\chi^0}&=&\left(\matrix{
 \tilde m_1  &0            &a   &b\cr
  0          &\tilde m_2   &c   &d  \cr
  a          &c            &0   &|\mu|e^{i\theta}\cr
  b          &d            &|\mu|e^{i\theta}   &0  \cr}\right)\label{neut}
\end{eqnarray}
  where $a=-M_Z sin\theta_W cos\beta$, $b=M_Z sin\theta_W sin\beta$,
$c=-cot\theta_W a$, $d=-cot\theta_W b$,
$tan\beta=v_2/v_1$ ($v_{1,2}=\mid <H_{1,2}>\mid$) and $\theta_W$ is the weak
mixing angle. The phase $\theta$ is given by
\begin{eqnarray}
\theta=\epsilon_1+\epsilon_2+\theta_\mu
\end{eqnarray} where at the electroweak scale,
$<H_{1,2}>=v_{1,2}e^{i\epsilon_{1,2}}$, and
$\mu=|\mu|e^{i\theta_\mu}$. The squark mass matrices may be written as
\begin{eqnarray} M_{\tilde q}^2&=&\left(\matrix{
 m^2_{q_L}               & e^{-i\alpha_q}m_q(|A_q|+|\mu| R_q
e^{i(\theta+\alpha_q)})
\cr
  e^{i\alpha_q}m_q(|A_q|+|\mu| R_q e^{-i(\theta+\alpha_q)})
&m^2_{q_R} }\right)\label{sqrk}
\end{eqnarray} where $m_q$, $e_q$ are the quark mass and
charge,\begin{eqnarray} m_{q_L}^2&=&m^2_Q+m_q^2+(1/2-e_q sin^2\theta_W)M_Z^2
cos2\beta\\ m_{q_R}^2&=&m^2_U+m_q^2+e_q sin^2\theta_W M_Z^2 cos2\beta\label{qlqr}
\end{eqnarray}
$A_q=|A_q|e^{i\alpha_q}$, $m_Q^2$ and $m_U^2$ are given in Ibanez et al.
\cite{iba} and $R_q=cot\beta(tan\beta)$ for u(d) quarks. Similar expressions
hold for slepton masses, with phases $\alpha_l$. Our sign conventions for $A_q$
and $\mu$ are those of Ref.\cite{bop}.

The condition for electroweak symmetry breaking is obtained by minimizing the
effective potential $V_{eff}$ with respect to $v_1$, $\epsilon_1$, $v_2$ and
$\epsilon_2$. The Higgs sector of $V_{eff}$ is
\begin{eqnarray}
V_{eff}&=&m_{1}^2v_1^2+m_{2}^2v_2^2+2|B\mu|cos(\theta+\theta_B)v_1v_2+
{g^2_2\over 8}(v_1^2+v_2^2)^2+{{g^{\prime}}^2\over
8}(v_2^2-v_1^2)^2+V_1\label{veff}
\end{eqnarray} where $V_1$ is the one loop contribution, $m^2_i=\mu^2+m^2_{H_i}$
and
$m_{H_{1,2}}^2$ are the $H_{1,2}$ running masses. For particles of spin j one has
\begin{eqnarray}V_1={1\over {64 \pi^2}}\sum_a C_a(-1)^{2 j_a} (2 j_a+1) m_a^4
(ln{m_a^2\over Q^2}-{3\over 2})\label{v1}
\end{eqnarray}  where $C_a$ is the color degree of freedom of the $a^{th}$
particle. In the following we choose the low energy scale $Q$ to be $m_t$ (175
GeV), and
include the full third generation states, $t$, $b$ and $\tau$ in $V_1$. This
allows an examination of the large tan$\beta$ domain.

We view the minimization equations of $V_{eff}$ as equations to determine the
Higgs VEVs i.e. $v_1,\, v_2,\, \epsilon_1,\, \epsilon_2$. Thus in the tree
approximation, the extrema equations $\partial V_{eff}/\partial
\epsilon_i=0$
yield $2|B\mu| sin(\theta+\theta_B)=0$, and hence the minimum of $V_{eff}$
requires
\begin{eqnarray}\theta=\pi-\theta_B\label{theta}
\end{eqnarray} (the choice of $\theta=-\theta_B$ leads to a maximum). The one
loop corrections depend only on the mass eigenvalues, and hence from
Eq.(\ref{sqrk}) only on the phase $\theta+\alpha_q$ and $\theta+\alpha_l$. Thus,
at the one loop level, one gets a correction to Eq.(\ref{theta}) of the
form\cite{demir,carm}
\begin{eqnarray}\theta=\pi-\theta_B+ f_1(\pi-\theta_B+\alpha_q,\pi-\theta_B+\alpha_l
)\label{thetaf}
\end{eqnarray}  where $f_1$ is the one loop correction with $\theta$
approximated by its tree  value, Eq.(\ref{theta}). As we will see, this
correction can become significant for large tan$\beta$. The EDMs however depend
on both the mass eigenvalues and the rotation matrices. Hence they will depend
on $\theta+\alpha_q$ and $\alpha_q$ separately, with $\theta$ determined by the
electroweak symmetry breaking condition in Eq.(\ref{thetaf}).

The current experimental 90$\%$ C.L. upper bounds on $d_n$\cite{dn} and
$d_e$\cite{de} are
\begin{eqnarray}(d_n)_{exp}&<&6.3\times 10^{-26}ecm;\,\,\,\,(d_e)_{exp}<4.3\times
10^{-27}ecm\label{dnde}
\end{eqnarray} In discussing the data, it is convenient to define the quantity
\begin{eqnarray}  K=log_{10}\mid{d_f\over {(d_f)_{exp}}}\mid
\label{k}
\end{eqnarray} Thus $K=0$ corresponds to a theoretical value which saturates the
current experimental bound, while $K=-1$ would represent the case, should the
experimental bounds be reduced by a factor of 10. As pointed out in
Refs.\cite{nath1,nath2,nath3,falk1,bk,falk2,kane1,kane2,bartl,pokorski},
various
cancellations can occur among the different contributions to the EDM. This is
illustrated in Fig.\ref{fig4} where $K$ is plotted vs. $m_0$ for the electron
dipole moment (eEDM).
\noindent We note that eventually for very large $m_0$, the curves fall below the
$K=0$ bound (as expected). Further the allowed range of $m_0$ (so that
$K\leq 0$) decreases with increasing tan$\beta$ (and the allowed range would
become very small should $K=-1$, i.e. the experimental bounds be reduced by a
factor of 10). In addition, the position of the dips moves to lower $m_0$ with
increasing tan$\beta$. This happens due to the fact that the regions on the right of
the dips are dominated by the contribution from the  chargino diagram.
 As tan$\beta$ increases, the contribution from the  chargino diagram increases
much more than the contribution from the neutralino diagram. However, a decrease
in $m_0$ increases the neutralino  diagram much more than the chargino diagram.
As a result the dips shift towards the smaller $m_0$ value for larger
tan$\beta$. We use negative values (arising from the factor $\pi$ in the Eq.
\ref{theta}) of $\mu$ to satisfy the experimental constraints
on the BR of $b\rightarrow s\gamma$ \cite{bsg1,bsg2}.

We begin our analysis by examining the RGE to determine what GUT scale parameters
lead to acceptable CP violating phases at the electroweak scale. We have used
Ref.\cite{bop} for the RGE relating the GUT scale parameters to the
parameters
at the electroweak scale. In general, the RGE must be solved numerically
and all
results given below are obtained by accurate numerical integration. However,
approximate analytic solutions can be found for low and intermediate tan$\beta$
(neglecting b and
$\tau$ Yukawa couplings) or for the SO(10) limit of very large tan$\beta$
(neglecting the $\tau$ Yukawa coupling). These analytic solutions give some
insight into the nature of the more general numerical results.

The solution for the $A_t(t)$ parameter in the low and intermediate tan$\beta$
case  can be cast in the form \cite {landau}
\begin{eqnarray}A_0&=&{A_R(t)\over D_0}+{H_3(t)\over F(t)}m_{1/2}\label{a0}
\end{eqnarray} where $A_R$ is the residue at the t-quark Landau pole,

\begin{eqnarray}A_R&=&A_t+m_{1/2}(H_2(t)-{H_3(t)\over {F(t)}})\label{ar}
\end{eqnarray}  and
\begin{eqnarray}D_0&=&1-6{F(t)\over E(t)}Y(t)\label{d0}.
\end{eqnarray} Here $t=2ln(M_G/Q)$, the form factors $E$, $F$, $H_2$ and $H_3$
are real and are defined in
\cite{iba} and $Y(t)=h_t/{16 \pi^2}$, where $h_t$ is the t-quark Yukawa coupling
constant. $D_0$ vanishes at the t-quark Landau pole (for Q=$m_t$,
$D_0\cong 1-(m_t/{200 sin\beta})^2)$ and is generally small (i.e.
$D_0\stackrel{<}{\sim} 0.2$ for $m_t=175$ GeV). The imaginary part of
Eq.(\ref{a0}) gives:
\begin{eqnarray}D_0|A_0|sin\alpha_{0A}&=&|A_t|sin\alpha_t\label{d0a0}.
\end{eqnarray} Thus even if $\alpha_{0A}=\pi/2$, $\alpha_t$ at the electroweak
scale will be generally suppressed due to the smallness of $D_0$, i.e. the RGE
naturally make the phase $\alpha_t$ small due to the nearness of the Landau
pole.
(This result has been previously observed, for low tan$\beta$ in
\cite{falk1}.) The approximate SO(10) solution with large tan$\beta$
(where
$Y_t\cong Y_b$) gives a similar analytic form for $A_t(t)\cong A_b(t)$  with the
factor 6 replaced by 7 in Eq.(\ref{d0}). Thus the suppression effect on
$\alpha_t$ occurs over the entire tan$\beta$ domain. These effects thus allow
$\alpha_{0A}$ to be large without a priori violating the experimental EDM bounds.

In contrast there is no analogous RGE suppression effect in the first generation
$A$ parameters $A_u$, $A_d$ and $A_e$. From the RGE, one finds that
$A_u$ has the form of Eq.(\ref{a0}) for low and intermediate tan$\beta$ with the
factor 6 replaced by 1 in Eq.(\ref{d0}). Now $D_0\simeq 1$ and no suppression
effect occurs for $\alpha_u$ at the electroweak scale in the analogue of
Eq.(\ref{d0a0}). However, the EDMs are not very sensitive to the first
generation $\alpha_u$, $\alpha_d$, $\alpha_e$ and one can have large values of
these parameters without violating the experimental bounds in a reasonable
region of the SUSY parameter space.

From the RGE, one sees that the phase of the $\mu$ parameter does not run at the
1-loop level, i.e.
$\theta_{0\mu}=\theta_{\mu}$, where $\theta_{\mu}$ is the phase at the
electroweak scale. Thus a large GUT scale phase will lead to a large electroweak
scale phase. However, as seen in Eqs.(\ref{char}-\ref{sqrk}),
$\theta_\mu$ enters only in the combination $\theta$, and $\theta$ is determined
in terms of $\theta_B$, $\alpha_l$ and $\alpha_q$ by the minimization of the
effective potential, as given in Eq.(\ref{thetaf}). Thus (with our choice of the
phases in the mass matrices) $\theta_{\mu}$ does not enter in the EDMs.

The EDMs, however are highly sensitive to $\theta_B$, which by
Eq.(\ref{thetaf}), enters into all the diagrams, and over most of the parameter
space, $\theta_B$ must be small to satisfy the EDM experimental bounds. This is
illustrated in Fig.\ref{fig5} where we have plotted the eEDM $K$ vs.
$\theta_B$. We see that $\theta_B$ can be large i.e. $\theta_B\simeq 0.08$, for
low $\tan\beta=3$, though it becomes smaller for higher
$\tan\beta$. However, even for tan$\beta=3$, the range $\Delta\theta_B$ of
$\theta_B$ so that $K\leq 0$ is very small i.e.
$\Delta\theta_B\stackrel{<}{\sim} 0.02$ (and $\theta_B$ would become quite fine
tuned if the experimental bound were reduced by a factor of 3 i.e.
$K=-0.5$). Figs.6 show that $\Delta\theta_B$ remains small for
$m_{1/2}\stackrel{<}{\sim} 350$ GeV ($m_{\tilde g}\stackrel{<}{\sim} 1$ TeV) and
only relaxes somewhat when $m_{1/2}\simeq 700$ GeV ($m_{\tilde
g}\stackrel{<}{\sim}2$ TeV). In this figure we have plotted the $K\le 0$
regions
as a function of $\theta_B$ and
$m_{1/2}$ for three different values of tan$\beta$. Since the contribution from
the chargino diagram increases with tan$\beta$, the larger tan$\beta$ requires
a smaller phase to produce the necessary cancellation between the chargino
and the
neutralino diagrams.

The RGE allow us to examine the significance of the above results at the GUT
scale. The RGE for $B$ and $A_t$ can be solved in the low and intermediate
tan$\beta$ regime to obtain
\begin{eqnarray}B&=&B_0-{1\over 2}(1-D_0)A_0-\Phi m_{1/2}\label{b}
\end{eqnarray}  where
\begin{eqnarray}\Phi&=&-{1\over 2}(1-D_0){H_3\over F}+ [3h_2+{3\over
5}h_1]{\alpha_G\over{4\pi}}
\label{phi}
\end{eqnarray}
$h_i=t/(1+\beta_it)$, $\alpha_G\cong 1/24$ is the GUT gauge coupling  constant
and $\beta_i$ are the MSSM beta functions. The real and the imaginary parts of
Eq.(\ref{b}) give: \begin{eqnarray}|B|sin\theta_B&=&|B_0|sin\theta_{0B}-{1\over
2}(1-D_0)|A_0|sin\alpha_{0A}\label{bsin}
\end{eqnarray}
\begin{eqnarray}|B|cos\theta_B&=&|B_0|cos\theta_{0B}-{1\over
2}(1-D_0)|A_0|cos\alpha_{0A}-\Phi m_{1/2}\label{bcos}
\end{eqnarray} Eqs.(\ref{bsin},\ref{bcos}) can be viewed as determining the
electroweak scale values of $|B|$ and $\theta_B$ in terms of the GUT input
parameters. Thus
$\theta_B$ depends upon both the initial phase $\theta_{0B}$ and the $A_0$ phase
$\alpha_{0A}$. Alternatively, one may use Eqs.(\ref{bsin},\ref{bcos}) to impose
electroweak scale phenomenological constraints on the allowed GUT parameters.
The requirement that the GUT theory gives rise to electroweak symmetry breaking
gives
\begin{eqnarray}|B|&=&{1\over 2}sin{2\beta}{m_3^2\over|\mu|}\label{bmag}
\end{eqnarray} where $m_3^2=\mu^2_1+\mu^2_2$ and
$\mu_i^2=|\mu|^2+m_{H_i}^2+\Sigma_i$. (Here
$m_{H_i}^2$ are the running Higgs masses at the electroweak scale, $|\mu|$ is
determined by electroweak breaking and $\Sigma _i$ are the one loop corrections
\cite{ap}.) In addition, we have seen that the experimental  bounds on the EDMs
restrict $\theta_B$ to be small, i.e.
$|\theta_B|\stackrel{<}{\sim} 0.1$ (and usually much smaller) and the allowed
range,
$\Delta\theta_B$, that satisfies the EDMs is quite small. These conditions put
severe constraints on the GUT scale theory which we now discuss.

Consider first the case where tan$\beta$ is very close to its minimum value
e.g. tan$\beta$=2. Here $1/2 sin2\beta$ is not small and by   Eq.(\ref{bmag}),
$|B|$ is of normal size. Eq.(\ref{bsin}) then implies that $\theta_B$ and
$\theta_{0B}$ are of roughly the same size (even for $\alpha_{0A}=\pi/2$) as
previously noted in Ref.\cite{falk1}. Thus a reasonable GUT theory
requires one
only to justify the size of
$\theta_{0B}$. However, as tan$\beta$ increases, $1/2 sin2\beta$
 decreases, and unless $\mu$ becomes anomalously small or $m_3^2$ anomalously
large, Eq.(\ref{bmag}) implies that $|B|$ becomes small. Thus the lefthand side
of Eq.(\ref{bsin}) becomes small (and to first approximation can be neglected),
and Eq.(\ref{bsin}) implies that $\theta_{0B}$ is determined by $\alpha_{0A}$
mainly, i.e. $\theta_{0B}$ is large if
$\alpha_{0A}$ is large. This result is illustrated in Fig.\ref{fig7} where
$\theta_{0B}$ is plotted vs. $\theta_B$. We see that since $\theta_B$ is small
for all tan$\beta$, $\theta_{0B}$ is in general large, e.g. even for tan$\beta
=3$ one has $\theta_{0B}\simeq 0.8$.

Returning to Eq.(\ref{bsin}), we see, however, that $\Delta\theta_{0B}$, the
range of values of $\theta_{0B}$ that will satisfy the EDM constraints, is very
small. Thus for fixed $A_0$, Eq.(\ref{bsin}) implies for large tan$\beta$ that
\begin{equation}
\Delta\theta_{0B}\cong{|B|\over {|B_0|}}\Delta\theta_B\ll\Delta\theta_B
\label{DthetaB}
\end{equation}
and since $\Delta\theta_B$ is small, $\Delta\theta_{0B}$ will be
very small. This is illustrated in Fig.\ref{fig8} where the allowed values of
$\theta_{0B}$ for $K\leq 0$ are plotted vs. $m_{1/2}$ for $d_e$ and $d_n$ for
tan$\beta$=3, 10, 20. We see that the phenomenological constraints at the
electroweak scale imply that $\theta_{0B}$ is both large and its value is
sharply fine tuned (unless $\alpha_{0A}$ is small, SUSY masses are large or
tan$\beta$ is small). One may alternatively view this from the ``top down'': if
Planck physics determines $\alpha_{0A}$ and $\theta_{0B}$ to be large and fixed,
then if the model is to achieve electroweak symmetry breaking with EDMs below
existing bounds, there will be a fine tuning of other GUT scale parameters
(unless tan$\beta$ is small or the SUSY masses are large).

As discussed above, there are a number of uncertainties in the calculation of
the neutron EDM. Fig.\ref{fig9} illustrates the effects of varying the d-quark
mass in a plot of the region allowed by the experimental EDM bounds ($K\leq 0$)
in the $m_0-m_{1/2}$ plane for $\theta_B=0$, $|A_0|=300$ GeV,
$\alpha_{0A}=\pi/2$, tan$\beta$=3. The bound on $d_e$ already excludes all
regions below the lower solid curve, while the $d_n$ bound excludes regions to
the left of the upward running curves for $m_d=5$ MeV (dotted), $m_d=8$ MeV
(solid), $m_d=12$ MeV (dashed). Thus the combined exclusion region increases from
$m_{1/2}\stackrel{>}{\sim} 260$ GeV to $m_{1/2}\stackrel{>}{\sim} 440$ GeV as
$m_d$ increases, a significant change. In the other figures in this paper we
have used the middle value of $m_d=8$ MeV corresponding to $m_s\simeq 150$ MeV.

In Fig.\ref{fig9}, we have chosen $\theta_B=0$. As one increases
$\theta_B$, one will eventually arrive at a neutralino-chargino cancellation
region, and this will reduce the values of $m_0$ and $m_{1/2}$ which are
excluded. Further, for low tan$\beta$ this cancellation can occur at relatively
large $\theta_B$ (i.e. $\theta_B\simeq 0.2$) as shown in Ref.\cite{falk1}.
Alternatively, increasing tan$\beta$ can also cause this cancellation, as the
loop corrections, Eq.(\ref{thetaf}), become large with large tan$\beta$ and the
contribution from the  chargino diagram increases with tan$\beta$.  This is
illustrated for tan$\beta$=20 in Fig.\ref{fig10} where the $d_n$ curve for $K=0$
(solid) bends downward below $m_0$=500 GeV when 200 GeV$<m_{1/2}<400$ GeV,
showing in this exceptional type situation that one can satisfy the EDM constraints
with a light particle spectrum and large tan$\beta$. Note that the cancellation
can occur for a relatively wide band of $m_{1/2}$, and persists even for
$K=-0.5$.

The regions of the SUSY parameter space that get eliminated by a joint
imposition of the experimental bounds on the $d_e$ and $d_n$ are sensitive, of
course, to the choice of the parameters. Thus in Figs.\ref{fig6}  one sees that
for $m_0=100$ GeV, one requires $m_{1/2}\stackrel{>}{\sim} 350$ GeV ($m_{\tilde
g}\stackrel{>}{\sim} 1$TeV) to jointly satisfy the $d_e$ and $d_n$ bounds, and
this requirement is roughly independent on tan$\beta$. Fig.\ref{fig11} shows
that if $m_0$ is increased to 250 GeV one requires now
$m_{1/2}\stackrel{>}{\sim} 160$ GeV ($m_{\tilde g}\stackrel{>}{\sim} 450$ GeV)
to satisfy both bounds, i.e. if $m_0$ is increased, the lower bound on $m_{1/2}$
is decreased as expected. However, if $|A_0|$ is increased, it is not necessary
that the allowed domain of $m_{1/2}$ increases. This is illustrated in
Fig.\ref{fig12} for the choice $|A_0|$=800 GeV, $m_0$=100 GeV. The increase of
$|A_0|$ raises the allowed values of $\theta_B$ for both $d_e$ and $d_n$,
but the raising of $d_e$ is larger causing the joint
allowed domain to require $m_{1/2}\stackrel{>}{\sim}  550$ GeV
($m_{\tilde{g}}\stackrel{>}{\sim} 1.5$ TeV).

In order to exhibit the size of $\mid\mu\mid$, we have plotted
$\mid\mu\mid$ as a function of $m_{1/2}$ for $m_0$=100, 300 and 700 GeV in
Fig.\ref{fig13}. The EDM constraints have not been imposed in
Fig.\ref{fig13}, but from Fig.6a, one sees that they are
satisfied at least for $m_0$=100 GeV for $d_e$ for essentially the entire
$m_{1/2}$ range (which allows $\theta_B\simeq 0.08$). Note that
$\mid\mu\mid$ is large (i.e. $\mid\mu\mid^2/M_Z^2\gg 1$) over the full
$m_{1/2}$ domain so that one is generally in the gaugino scaling domain.

\section{Nonuniversal soft breaking}

There are three types of nonuniversalities that might be considered in SUGRA
models: nonuniversal gaugino masses at $M_G$, nonuniversal scalar masses at
$M_G$ and generation off diagonal cubic soft breaking matrices
$A_{ij}$ and scalar mass matrices $m^2_{ij}$. We will discuss here briefly the
first two possibilities, a more detailed analysis will be given in a later paper
\cite{aad}. (The third type of non universality has the possibility of
generating $\epsilon$ and $\epsilon^{\prime}/\epsilon$ CP violations independent
of the CKM phase \cite{epsilon}.)

In general, the U(1), SU(2) and SU(3) gaugino masses at $M_G$ can have the form
\begin{eqnarray} m_{1/2i}&=&|m_{1/2i}|e^{i\phi_i};\, i=1,2,3
\label{mg}
\end{eqnarray} A convenient phase choice is to set the SU(2) phase to zero:
$\phi_2=0$, and hence the chargino mass matrix, Eq.(\ref{char}), is
unchanged. The remaining phases produce effects in the RGE for $A_t(t)$:
\begin{eqnarray} A_t&=&D_0A_{0}+\sum_i\Phi_{ti}|m_{1/2i}|e^{i\phi_i}
\label{At}
\end{eqnarray} and hence the imaginary part gives
\begin{eqnarray}
|A_t|sin\alpha_t&=&|A_0|D_0sin\alpha_{0A}+\sum_i\Phi_{ti}|m_{1/2i}|sin\phi_i
\label{At2}
\end{eqnarray} Thus the Landau pole factor $D_0$ still suppresses any large phase
$\alpha_{0A}$ at the GUT scale, though $\alpha_t$ may become large due to the
$\phi_i$ phase which does not run at one loop order and hence is not suppressed.
This phenomena can affect $d_n$.

A second effect occurs in the neutralino mass matrix which now reads
\begin{eqnarray} M_{\chi^0}&=&\left(\matrix{
 |\tilde m_1|e^{i\phi_1}  &0            &a   &b\cr
  0          &\tilde m_2   &c   &d  \cr
  a          &c            &0   &|\mu|e^{i\theta}\cr
  b          &d            &|\mu|e^{i\theta}   &0  \cr}\right)\label{neut2}
\end{eqnarray} where $\theta =\epsilon_1+\epsilon_2+\theta_\mu$. Thus $\phi_1$
effects any neutralino-chargino cancellations in $d_e$ and $d_n$. An example is
illustrated in Fig.\ref{fig14}. One sees that $\theta_B$ can be large but is
still tightly constrained to satisfy the experimental EDM bound (i.e.
$\Delta\theta_B\simeq 0.02$). One is again led to a fine tuning at the GUT
scale.

As an example of nonuniversal sfermion masses, we consider an SU(5)-type model
where the GUT group possesses an SU(5) subgroup with matter embedded in the
usual way in $10+\bar 5$ representations. We also assume the first two
generations remain universal to suppress flavor changing neutral currents. The
Higgs masses have the form at $M_G$
\begin{eqnarray}
m_{H_1}^2&=&m_0^2(1+\delta_1);\,\,\,\,m_{H_2}^2=m_0^2(1+\delta_2)
\label{scanonuni}
\end{eqnarray} while the third generation $10$ representation (containing
$\tilde t_L$,
$\tilde t_R$,$\tilde b_L$, $\tilde\tau_R$) and the $\bar 5$ (containing
$\tilde b_R$,
$\tilde\nu_\tau$, $\tilde\tau_L$) masses are parametrized by
\begin{eqnarray} m_{10}^2&=&m_0^2(1+\delta_{10});\,\,\,\,m_{\bar
5}^2=m_0^2(1+\delta_5)
\label{scanonuni1}
\end{eqnarray} Here the $\delta_i\ge-1$ represent the deviations from the
universality. In addition the third generation cubic soft breaking parameters are
$A_{0t}$ and $A_{0b}=A_{0\tau}$.

One may use the non universal sfermion and Higgs masses to soften the effects of
the experimental EDMs. This is illustrated in Fig.\ref{fig15} where the K=0
constraint is imposed in the $m_0-m_{1/2}$ plane for a choice of
$\delta_i$ parameters. For $\delta_i=0$, the $d_n$ curve would continue to rise
for small $m_{1/2}$ as seen in Fig.\ref{fig9}, rather than turn over as in
Fig.\ref{fig15}. For example, for $m_{1/2}$=220 GeV and tan$\beta$=5,
$K(d_n)<0$ occurs for $m_0>580$ GeV
in the non universal case,
whereas $K(d_n)<0$ occurs for $m_0>750$ GeV in the universal case.
The particular set of
$\delta_i$s in the Fig.\ref{fig15} is chosen to reduce the stop and the sbottom
mass and to
satisfy the requirement of the radiative electroweak symmetry breaking.
 The lower stop and sbottom mass increases the Weinberg type diagram which has a subtractive
 effect in the net nEDM magnitude.

\section{SUGRA Model Of Small Phases}

While choosing the CP violating phases to be small, i.e. $\phi_i\simeq 10^{-2}$,
to satisfy the EDM constraints appears to be artificial in the low energy MSSM,
it is possible that a natural model of this type can arise in SUGRA GUT models.

SUGRA models depend on three functions of the chiral fields
${\Phi_{\alpha}}$: $f_{\alpha\beta}(\Phi_{\alpha} )$ the gauge kinetic energy,
$K(\Phi_{\alpha} ,\Phi_{\alpha}^{\dag} )$ the Kahler potential and
$W(\Phi_{\alpha} )$ the superpotential. One assumes a hidden sector exists where
some fields, e.g. moduli or dilaton, $\Phi_i$ grow VEVs of Planck size to break
supersymmetry i.e.
\begin{eqnarray} x_i&=&\kappa <\Phi_i>=O(1)
\label{xi}
\end{eqnarray} where $\kappa^{-1}=M_{Pl}=2.4\cdot 10^{18}$ GeV. Thus one can
write
$\{\Phi_{\alpha}\}=\{\Phi_i, \Phi_a\}$ where $\Phi_a$ are the physical
sector
fields. One might expand the Kahler potential in a power series of the physical
fields:
\begin{eqnarray} K&=&\kappa^{-2}c^{(0)}+(c_{ab}^{(2)}\Phi_a\Phi_b+{1\over M}
c_{abc}^{(3)}\Phi_a\Phi_b\Phi_c +...) +(\tilde c_{ab}^{(2)}\Phi_a\Phi_b^{\dag}
+{1\over M}
\tilde c_{abc}^{(3)}\Phi_a^{\dag}\Phi_b\Phi_c+...)
\label{kahler}
\end{eqnarray} where $M$ is a large mass. The $c^{(i)}, \tilde c^{(i)}$ are
dimensionless and we assume them to be $O(1)$. The first parenthesis is
holomorphic, and can be transferred to the superpotential by a Kahler
transformation giving rise to the $\mu$ term:
\begin{equation} W\rightarrow W +
\kappa^2W(c_{ab}^{(2)}\Phi_a\Phi_b+{1\over
M}c_{abc}^{(3)}\Phi_a\Phi_b\Phi_c+...)
\label{W}
\end{equation} The leading additional terms on the right occur when $W$ there is
replaced by its
VEV ($\kappa^2<W>\simeq m_{3/2}$) and one of the fields in the cubic term,
e.g.
$\Phi_c$, has a GUT scale VEV. Then $W\rightarrow W+\mu_{ab}\Phi_a\Phi_b$,
where
\begin{eqnarray}
\mu_{ab}&=&(c_{ab}^{(2)}+{M_G\over M}c_{abc}^{(3)})m_{3/2}
\label{muab}
\end{eqnarray} In the perturbative heterotic string we expect $M\simeq M_{Pl}$
and hence
$M_G/M\simeq O(10^{-2})$. If one assumes $c^{(2)}$ is real and $c^{(3)}$ is
complex with arbitrary size phase, then $\mu$ naturally has a phase of size
$\theta_{0\mu}\simeq O(10^{-2})$. A similar analysis can be done for the other
parameters (A, B, $\tilde m$) yielding automatically small phases. One assumes
only that the leading (renormalizable) term in the matter expansion is real,
while the higher terms (presumably arising from integrating out the tower of
massive states) can have arbitrary phases.

Other string scenarios can have $M\simeq O(M_G)$. In that case one expects the
nonrenormalizable terms to give rise to large phases, unless the low lying
members of the tower of heavy states do not couple to the fields of the physical
sector.

\section{Conclusions}

In SUGRA GUT models, the basic parameters are given at the GUT scale. The
renormalization group equations then carry this information to the electroweak
scale, allowing for experimental tests of the model. Thus what represents a
``natural'' choice of parameters is presumably a GUT scale question.

The experimental values of the electric dipole moments put strong constraints on
the low energy SUSY parameter space. However, for mSUGRA models at least, the
RGE suppress the value of the $A_t$ phase at the electroweak scale, $\alpha_t$,
so that any size of phase, $\alpha_{0A}$, even $\pi/2$ at
$M_G$ will generally lead to acceptable EDMs. This effect is due to the nearness
of the t-quark to its Landau pole.

The EDMs, however, are very sensitive to the B phase at the electroweak scale,
$\theta_B$. In the domain of light SUSY mass spectrum (i.e.
$\stackrel{<}{\sim}$1 TeV) and $\alpha_{0A}\simeq O(1)$, the EDMs allow a
large
$\theta_B$ (i.e.
$O(10^{-1}))$ only for low tan$\beta$, i.e. tan$\beta\stackrel{<}{\sim}3$.
Further, for fixed SUSY parameters in this domain, the allowed range
$\Delta\theta_B$ of
$\theta_B$ needed to satisfy the EDM constraints is very small for
tan$\beta\stackrel{>}{\sim}3$.

The above results strongly affect the GUT theory. Thus for $\alpha_{0A}$ large
and a light mass spectrum one finds at the GUT scale that $\theta_{0B}$ is
generally large, i.e. $O(1)$ even for low tan$\beta$, and the allowed range
$\Delta\theta_{0B}$ for satisfying the EDMs is exceedingly small for
tan$\beta\stackrel{>}{\sim} 3$. Thus for this situation, $\theta_{0B}$ must be
chosen large and be very precisely determined as tan$\beta$ gets large, leading
to a new fine tuning problem at the GUT scale. The origin of this difficulty
resides in the requirement that the GUT models give rise to parameters at the
electroweak scale that simultaneously satisfy the experimental EDM constraints
and give rise to radiative electroweak breaking.

The above discussion suggests that reasonable mSUGRA models with large phases
and mass spectrum below 1 TeV can only be constructed for small tan$\beta$ i.e.
tan$\beta\stackrel{<}{\sim}3$. The current LEP data combined with
cosmological constraints already require
tan$\beta\stackrel{>}{\sim} 2$ \cite{lep}, and LEP, the Tevatron and the LHC
 will be
able to explore the domain tan$\beta\stackrel{<}{\sim} 50$\cite{wag}. Thus the value of tan$\beta$ (or a
lower bound) is something that will be determined experimentally in the not too
distant future. Should tan$\beta$ turn out to be large, it is possible that only
models with small GUT phases $\phi_{0i}\simeq O(10^{-2})$ are reasonable, and it
was shown above that such models can arise naturally with $\phi_{0i}\simeq
O(M_G/M_{Pl})$.

Our analysis has included the loop corrections to the radiative breaking
condition. While these corrections are small, they grow with tan$\beta$
and are competitive to the $\theta_B$ contribution to the EDMs. Thus for
large tan$\beta$ (e.g. tan$\beta\simeq 20$) they can produce a new
cancellation phenomena (with $\theta_B$=0), allowing the EDM constraints to be satisfied for a
light mass spectra (e.g. $m_0\simeq$300 GeV, $m_{1/2}\simeq$ 200 GeV)
which for lower tan$\beta$ would be excluded by the EDM data.

The combined restrictions of the experimental EDM bounds for the neutron
and electron put strong constraints on the SUSY parameter space. However,
since the theoretical cancellations needed to satisfy experiment are
delicate, as one varies the parameters one can get rather different
excluded regions. Thus, while we have found that increasing $m_0$ gives an
allowed region of smaller $m_{1/2}$ (as one might have expected),
increasing $\mid A_0\mid$ can require $m_{1/2}$ to be increased to satisfy
the EDM bounds. However, it should be stressed that there are significant
theoretical uncertainties in the calculation of $d_n$, and we have seen,
for example, that the current uncertainty in the quark masses $m_u$ and
$m_d$ (arising from the uncertainty in $m_s$) can lead to a significant
variation in the domain the neutron EDM data can exclude.

Much of the above discussion holds also for non minimal models with non
universal gaugino and sfermion soft breaking masses. A more detailed discussion
of the effects of non universalities will be given elsewhere
\cite{aad}.

\section{Acknowledgement}

This work was supported in part by National Science Foundation Grant No.
PHY-9722090. One of us (R.A.) is pleased to acknowledge the hospitality of the
Korea Institute for Advanced Study where part of this paper was written.
 We would also like to thank Toby Falk and Keith Olive for valuable discussions.