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

Recent advances in submicrometer technology have made it possible to
fabricate relatively large arrays of ultrasmall SIS
(superconductor to insulator to superconductor) Josephson
junctions\cite{moij,gee,herre,herre2,tighe,chalmers}. These arrays
consist of Josephson junctions made of superconducting islands separated by
an insulating barrier\cite{moij}. The areas of these junctions can vary
from a few microns to submicron sizes. In the latter case the effective
capacitance of the junctions can be smaller than  femtoFarads (fF=$10\times
10^{-15}$ Farads).
Under these circumstances the long range phase coherent properties of the
Josephson junction arrays (JJA)
depend crucially on the interplay between the Josephson ($E_J$) and
charging energies ($E_C$). In this paper we will be concerned with
studying this interplay for the specific parameter ranges of the experiments
carried out in Delft\cite{herre2}. The charging energy associated with adding
a single charge to an island is $E_{C_{\rm s}}=\frac{e^2}{2C_{\rm s}}$,
whereas the corresponding energy necessary to transfer a charge from an
island to a nearby one is $E_{C_{\rm m}}=\frac{e^2}{2C_{\rm m}}$. Here $e$
is the electronic charge. For the Delft samples, the
self capacitances are  typically $C_{s}\sim 3\times 10^{-18}$F,
while the mutual capacitances are on the order
of $C_{m}\sim 1\times 10^{-15}$F. This means that $C_m$ is
three orders of magnitude larger than $C_s$. Prior to the advent of these
types of JJA, most theoretical studies assumed
that the dominant term was the self-capacitive contribution, since the
calculations are easier and also no significant changes in the
results were expected from having the extra $C_m$ contribution. In the
Delft experiments
a phase diagram was obtained as a function of temperature vs the quantum
parameter
\begin{equation}
\label{alfa}{\alpha\equiv \frac{E_C}{E_J}}.
\end{equation}
The phase diagram that includes the experimental results plus our theoretical
results  is shown in Fig. 1. When $\alpha=0$, the array is modeled by the
two-dimensional XY model and it has been found to be a faithful
representation \cite{bkt2} of the
Berezinskii-Kosterlitz-Thouless (BKT) scenario\cite{bkt,jkkn}. As $\alpha$
increases (in the Delft experiments the values of $\alpha$ cover the
range $0.13$ to $4.55$, while in the Harvard samples the $\alpha$ can be as
large as 33),
the critical temperature decreases monotonically. As $\alpha$ increases
further, at fixed temperature, one is expected to change from a
superconducting phase characterized by quasi-long range phase coherence (SC)
to a charge dominated insulating phase (I). For small and fixed $\alpha$,
as temperature increases, we move from a SC phase to a resistive or
conducting one. For larger values of $\alpha$ the transition is expected to
be between an insulating and a conducting phase.
In the insulator normal phase boundary,
the idea is that in the case when the $C_s=0$, a mutual capacitance model
would map to a two-dimensional Coulomb gas model and thus a charge unbinding
BKT transition would ensue\cite{yoshi,moij2,fazio}. This scenario has not
been found experimentally since in practice, although $C_s$ can be
three orders of magnitude smaller than $C_m$, the electrostatic screening
length is much
smaller than the lattice size and the possible BKT transition is masked.
This means that experimentally the boundary I to N is strictly speaking not
a thermodynamic  phase boundary but a cross over region\cite{tighe,herre2}.
The situation is different in the SC to N boundary
since that is a true thermodynamic phase boundary. An important question
is, however, if the models that have been proposed to describe the JJA
are faithful representations of the experimental systems  studied.
This question is non
trivial since there is evidence for having extraneous elements in the
arrays, e.g. stray charges, that have not been accounted for properly in the
models so far considered. It is for this reason
that we have thought it important to carry out detailed {\bf quantitative}
calculations of the SC to N phase boundary  and to make detailed
comparisons with experiments.
We can already note that our results shown in Fig.1, coming from a
WKB renormalization group calculation,
fit the experimental results quite well in the small $\alpha$ range. On the
other hand, a rather  good fit to the experimental results
comes from carrying out a nonperturbative quantum Monte Carlo  calculation.
In this paper we describe mainly the WKB-RG calculation, and we will briefly
mention the QMC results, which will be discussed in more detail
elsewhere\cite{jens,cristian}.

An important result emerging from the WKB-RG calculation is the possible
existence of a low temperature instability. This $QUIT$ (QUantum fluctuation
Induced Transition) was originally found in a self-capacitive
model\cite{jos84,jacobs}. Here we show that the results are also true in the
case when $C_m$ dominates. This extension is non trivial since the
self-capacitive model {\bf has no} I to N phase boundary with a BKT type charge
unbinding transition, even theoretically. One can then wonder if the
competition between the insulating and
superconducting phases at sufficiently low temperatures may quench the QUIT.
Within the context of the  WKB-RG calculation this is, in fact,
not the case as we discuss in this paper.
{}From the WKB-RG analysis the QUIT instability may be conjectured to be
a phase boundary from a SC to a normal phase\cite{jos84}.

There was early support for the existence of the QUIT in the $C_m=0$,
$C_s\neq 0$ case from nonperturbative quantum Monte Carlo
simulations\cite{jacobs}.
In the $f=0$ case it was found that the helicity modulus had a discontinuity
between two states with {\bf finite} superfluid densities each. This lead to
the surmise that the QUIT is in fact a transition between two superfluid
states, one dominated by thermal and the other by quantum fluctuations.
Later studies considered the case of $f=1/2$, the fully
frustrated limit, which gave a larger discontinuity in the superfluid density
and at higher $T_{QUIT}(f=1/2)$ temperature\cite{jacobs}. As we shall see
from the  WKB-RG calculation described
below, the $T_{QUIT}(f=0)\sim T_{QUIT}(f=1/2)$. However, what was found
in the QMC calculations was that $T_{QUIT}(f=1/2)\sim 10\times T_{QUIT}(f=0)$.
In a recent calculation which uses an improved QMC algorithm, that
reduces the slow approach to equilibrium\cite{novo}, the large
hysteresis loop seen before in the $f=1/2$ case
was suggested to be directly related
to the long lived metastable states seen in the
classical fully frustrated XY model. The authors, however, did not carry out
the calculations to sufficiently low temperature to test the $f=0$, or
$f=1/2$ QUIT. The question then remains as to the physical nature of the QUIT,
if it does exist. From the WKB-RG analysis there is clearly a low
temperature QUIT instability. The known QMC results for the $f=0$,
 $C_m=0$ case suggest a QUIT between two different superconducting states.
Further studies are needed to ascertain the nature of the putative QUIT.
As we also mention later in this paper there is an experimental low
temperature instability that could very well be a manifestation of the QUIT.

Here we shall report on our main WKB-RG and QMC results and sketch how they
were obtained while a fuller more detailed description will appear
elsewhere\cite{jens,cristian}.

The outline of the rest of this paper is the following:
In Sec. II we define  the model studied here and discuss
the path integral representation of the partition function used in our
analysis. In Section III we outline our WKB-RG analysis in the case where the
Josephson energy dominates. We obtain generalized RG equations from which
we obtain the phase diagram, which is then compared with the experimental
results\cite{herre2}, as well as with the Monte Carlo simulation
results\cite{jens}. As we mentioned above, if
$\xi _E$ is large in the insulating phase the charges interact with a
logarithmic potential and a BKT type transition must
follow\cite{yoshi,moij2,fazio}. In section IV we address the
question as to what happens when we perturb the Coulomb gas model ($C_s=0$)
with the Josephson term in the hamiltonian; Is there a QUIT in the charge
dominated phase? We have carried out a WKB-RG perturbative analysis
in the insulating region and found that {\bf there is no analog}
to the QUIT transition found in the SC phase,
at least to the lowest order in the perturbative analysis.

\section{\bf {The model}}

A Josephson junction array can be modelled by a periodic lattice of
superconducting islands separated by insulating barriers. Each island is
characterized by a Ginzburg-Landau order parameter $\Psi (\vec r_i)=
|\Psi _0(\vec r_i)|e^{i \phi (\vec r_i)}$, where $\vec r_i$ is a
two-dimensional vector denoting the position of each island.
Each one of the islands becomes independently superconducting about the
bulk transition temperature $T_{c0}$. When the temperature is lowered further,
the magnitude of the order parameter, $|\Psi _0(\vec r_i)|$, is nonfluctuating
and the onset of long range phase coherence is responsible for the
zero resistance in the arrays. The onset of phase coherence is due to the
tunneling of Josephson currents between the islands and it is characterized
by a sharp drop to the zero resistance state. The onset temperature
in the arrays can be significantly  modified by making the junctions
ultrasmall. When the junction's capacitance is small the charging energy,
i.e. the energy necessary to transfer a Cooper pair between the islands,
can be large to the point where no Cooper pairs can tunnel any more,
and thus the Josephson current can be reduced to the point
where it is completely quenched. The competition between the Josephson
tunneling and the charging energy, in zero field, can be modeled by
the hamiltonian,
\begin{eqnarray}
\label{ha}{\hat {\cal H}}=H_C+H_J
&=&
{{q^2}\over{2}}\sum_{\vec r_1,\vec r_2}\hat n(\vec r_1){\bf {\bf C}^{-1}
(\vec r_1,\vec r_2)}\hat n(\vec r_2)+
E_J \sum_{<\vec r_1,\vec r_2>}[1-
\cos(\hat \phi(\vec r_1)-\hat \phi(\vec r_2))]
\end{eqnarray}
where $q=2e$.  Here $\hat \phi (\vec r_i)$ is the quantum phase operator
while $\hat n (\vec r_i)$
is its canonically conjugate number operator, which measures the excess number
of Cooper pairs in the island placed at $\vec r_i$. These canonically
conjugate operators satisfy the commutation relations
\begin{equation}
[\hat n(\vec r_1),\hat\phi(\vec r_2)]=-i\delta_{\vec r_1,\vec r_2},
\end{equation}
which imply that in the phase representation we can write $\hat n(\vec r_i)=
-i\frac{\partial}{\partial {\phi(\vec r_i)}}$. A Josephson junction is
essentially a capacitor with an insulating barrier. The geometric capacitance
matrix ${\bf C (\vec r_1,\vec r_2)}$ is assumed to include only the island
capacitance with respect to the ground plane, $C_{\rm s}$, and the
mutual capacitance between nearest neighbor islands, $C_{\rm m}$.
In the square array case, of interest here,
we can write the capacitance matrix as
\begin{equation}
\label{cap}{\bf
C(\vec r_1,\vec r_2)} = (C_{\rm s}+4C_{\rm m})\delta_{\vec
r_1,\vec r_2} + C_{\rm m} \sum_{\vec d} \delta_{\vec r_1,\vec r_2+\vec
d}\,\, ,
\end{equation}
where the sum in the second term is over the nearest neighbor
vectors $\vec d$. The electrostatic potential associated with this capacitance
matrix model leads to the screening length
$\xi_E=\sqrt{\frac{C_{\rm m}}{C_{\rm s}}}$.
This means that when $C_{\rm m}=0$ the charges are maximally screened.
In the arrays fabricated so far the screening length, as calculated from
Eq.~(\ref{cap}) is between $15$ and $20$ plaquettes\cite{herre2,tighe}.
This means that screening
effects in these arrays are important and that including $C_{\rm m}$
in the analysis is essential. As mentioned in the introduction, however,
prior to the fabrication of the arrays with
ultrasmall capacitances, many theoretical studies concentrated either on
mean field approaches\cite{simanek,fishman}, which neglect the special
nature of the BKT ordering in two-dimensions, or in 2-D studies with
$C_{\rm m}=0$ \cite{jos84,jacobs}. Initial MC studies of
the off diagonal problem have also appeared\cite{jens2}. Since the
samples fabricated have a typical ratio $C_{\rm m}\sim
 {C_{\rm s}}10^{3}$, it is clear that if one wants to make a direct
 quantitative comparison to
experimental results we need to extend the analysis to include
$C_{\rm m}$ explicitly. This is the purpose of this paper.

We should stress that the model defined by Eqs.(2-4) is not the most general
model one could consider, for it assumes that all the electric fields are
confined to the two-dimensional plane of the array. This may not be quite
the case for real samples but one would expect that in the regime where
$\xi _E$ is sufficiently large this model should be appropriate.
As we shall discuss below the fact that one is able to make successful
comparisons between experiment and theory in the {\it superconducting
to normal} regime leads us to believe that in that regime at least
the model gives a correct representation of the experimental system.

There are other possibly important elements missing in the model studied
in this paper. To wit: (i) We do not consider here the influence of a constant
external magnetic field. The physics in that case, even in the classical
limit, is highly nontrivial and requires a special treatment\cite{ramirez}.
(ii) We do not include the self-induced magnetic fields, that can also lead to
interesting new physics\cite{daniel}. The reason is that  the
critical currents in the
arrays are found to be too low to lead to any significant magnetic field
screening effects. (iii) We do not include quasi-particle
dissipation. Here again it has been
found experimentally that in the SC to N phase these effects are not important
\cite{moij2}. A low temperature study that includes
quasi-particle dissipation for the diagonal capacitance matrix model
has shown that dissipation counteracts the charging effects so as to
strengthen the Josephson tunneling\cite{choi}. (iv) We could also
consider the effect of disorder, in particular the one due to random stray
charges in the array. These random frustration charges may also be more
important in the insulating phase than in the superconducting one, which is
the one we concentrate on in this paper.

Here we are interested in calculating the thermodynamic properties of the
model defined by $\hat {\cal {H}}$. The corresponding partition function is
defined by
\begin{equation}
Z = {\rm Tr}\left\{ e^{-\beta\hat {\cal H}}\right\}\label{zeta}.
\end{equation}
The trace is taken over  $\hat \phi$ or $\hat n$. It is convenient, for
calculational purposes, to change from the operator expression to the
imaginary-time Feynmann path integral representation of $Z$. We accomplish
this by using
\begin{equation}
           <n(\vec r_1)|\phi(\vec r_2)> = \delta_{\vec r_1,\vec r_2}
           \frac{\exp\{ i n(\vec r_1)\phi(\vec r_1)\}}{\sqrt{2\pi}}.
\end{equation}
For this problem we follow steps parallel as those described in Ref.
{}~\cite{jacobs}.
This means that we discretize the imaginary-time axis into
$L_\tau$ time slices separated by a distance $\epsilon =
\frac{\beta \hbar}{L_\tau}$. We  recover, in principle, the fully
quantum results in the limits $L_\tau\to\infty$, $\epsilon\to 0$ with
$\beta \hbar$ kept fixed. Following Ref.~\cite{jacobs} we get  the
partition  function expression (up to order $O(1/L_{\tau})$)
\begin{equation}
Z = \prod_{\tau=0}^{L_{\tau}-1} \sqrt{{\rm det}[ C ] }
    \prod_{\vec r} \int_{0}^{2\pi}\!\!\sqrt{\frac{L_{\tau}}
    {2\pi\beta q^2} } d\phi(\vec r,\tau)\!\!\!\!\sum_{m(\vec r,\tau)
     = -\infty}^{\infty}\!\!\!\!\exp\bigg[-\frac{1}{\hbar}S[\{\phi\},\{m\}]
     \bigg]
     \label{eq(2.3)}
\end{equation}
where we defined the action
\begin{eqnarray}
\frac{1}{\hbar}S[\{\phi\},\{m\}] = \frac{ \beta}{L_{\tau}}
     \sum_{\tau = 0 }^{L_{\tau}-1} \Bigg[& &H_J(\{\phi(\vec r,\tau)\})
     -\frac{L_{\tau}}{2\beta q^2}\!\!\!\sum_{<\vec r_1,\vec r_2>}
       [\phi(\vec r_1,\tau\! +\! 1) - \phi(\vec r_1,\tau)
       + 2\pi m(\vec r_1,\tau) ] \times \nonumber \\
       & &\times{\bf C}(\vec r_1,\vec r_2)[\phi(\vec r_2, \tau\! +\! 1)
       -\phi(\vec r_2,\tau) + 2\pi m(\vec r_2,\tau)]\Bigg] +
               \nonumber \\
               & &+ O(1/L_{\tau}).
               \label{eq(2.4)}
\end{eqnarray}
Here $H_J(\{\phi(\vec r,\tau)\})$ is the Josephson hamiltonian in terms
of the phase variable $\phi(\vec r,\tau)$.
The quantum nature of this c-number functional integral comes from imposing
the periodic boundary condition
\begin{equation}
\phi(\vec r,L_\tau) = \phi(\vec r,0).       \label{eq(2.5)}
\end{equation}
These are the basic equations used in the  WKB-RG and Monte Carlo
calculations described in the next section.

\section{\bf {WKB-RG study about the SC to N phase boundary}}

Our approach here is to perturb the physics described by
$H_J(\{\phi(\vec r,\tau)\})$, in the limit when the charging energy is
small. When $\alpha=0$, $H_J$ describes the physics of the classical
2-D XY model\cite{bkt,jkkn}. In this case we have the BKT scenario that
depends on the
thermal nucleation of vortex-antivortex pairs (VAP). The density of
VAP increases exponentially as the temperature rises until they unbind
at the critical temperature
$T_{\rm BKT}(\alpha=0)=T_{\rm BKT}^{(0)}= (\pi E_J/2k_{\rm B})$.
The BKT scenario is best understood in terms of a renormalization
group (RG) analysis\cite{bkt,jkkn}. The RG
flow diagram is obtained from a perturbation expansion in powers of the
vortex pair density $y=y_0e^{-\frac{\pi ^2}{2}K}$. Here $K=\beta E_J$ and
$y_0$ is the initial condition for the bare vortex pair density.
In the standard BKT picture there is a line of fixed points for $0\leq T\leq
T_{\rm BKT}^{(0)}$, with algebraically decaying correlation functions.
In the self-capacitive model, at $T=0$, one can map the problem to an
anisotropic three-dimensional XY model, which must have a standard phase
transition at a critical value of $\alpha ^c_s$. Around the $T=0$
critical point $\alpha_s^c$, one expects to have exponentially decaying
correlation functions while at $\alpha_s^c$ the correlations must decay
algebraically. The important question
is then: How do we go  from the BKT regime, with algebraic decaying correlation
functions to the very low temperature exponentially decaying one? There must
be a discontinuity in going from one limit ($\alpha _s=0$, $T\neq 0$) to the
other one ($\alpha _s\neq 0$, $T=0$). Below we discuss
the evidence we have found, including a possible experimental candidate,
that there may indeed be a QUIT at low temperatures.

The situation  when $C_s=0$ and $C_m\neq 0$  is actually
quite different. In that case the $T=0$ limit can be approximately
represented by two
coupled three-dimensional XY models, one describing the phase degrees of
freedom and the other the charges. As a function of $\alpha_m$ we can go
from a phase dominated region, with a 3-D type XY model critical properties
to one dominated by a 3-D Coulomb gas.
There is not much known about the critical properties,
and in particular the correlation function behavior, of the two coupled
3D XY models. So, strictly speaking, we can not state what kind of crossover
we should expect when going from $T=T_{BKT}(\alpha _m,\alpha_s)$ to $T=0$.
Some understanding of the physics in this limit can be obtained by using the
Villain transformation, both for the charging energy term and for the phase
contribution\cite{fazio}. The Villain approximation is, however, valid only in
a restricted range of $\alpha$ values which do not cover the full experimental
range. One could conjecture, however, that the properties of the
Villain approximated models is in the same universality class as the full
coupled XY models, which is in fact the case when $\alpha=0$ \cite{jkkn},
but this needs to be explicitly shown.
Furthermore, the general case treated here where both $C_s\neq 0$ and
$C_m\neq 0$ is more complicated since the effective Coulomb gas in the
insulating phase has a finite screening length. All these issues need to be
studied further.

To find  the corrections to the BKT scenario due to the charging effects
we carry out a semiclassical or WKB analysis of the model. This was
originally done for the self-capacitive model in Ref.~\cite{jos84}.
Here we follow a similar approach, except that technically the problem is
more demanding.

To evaluate the partition function in the SC to N, or small $\alpha$ regime, we
notice that we can extend the range of integration of the phases in
Eq.~(\ref{eq(2.3)}) from
[0,$2\pi$] to [$-\infty,\infty$], while at the same time all but one of
the summations over the set of $\{m\}$'s can be eliminated.
The resulting expression for the partition function, in the
$L_\tau\rightarrow\infty$ limit, then reads
\begin{eqnarray}
           Z = \sqrt{{\rm det}[{\bf C}]} \int_{0}^{2\pi}
               \prod_{\vec r}
               \sqrt{\frac{L_{\tau}}{2\pi \beta q^2}} d\phi(\vec r,0)
               \!\!\!\!\sum_{m(\vec r )=-\infty }^{\infty} \int_{-\infty}^
               {\infty} & &\prod_{\tau = 0}^{L_{\tau}-1}\!\!\sqrt{{\rm det}
               {\bf [C]}}
               \prod_{\vec r} \sqrt{\frac{L_{\tau}}{2\pi \beta q^2}}
               d\phi(\vec r,\tau) \times \nonumber \\
               & &\times\exp\bigg[-\frac{1}{\hbar}\int_{0}^{\beta\hbar}
               d\tau L_E \bigg],
           \label{eq(3.1)}
\end{eqnarray}
where the euclidean lagrangian is
\begin{equation}
L_E = \frac{1}{2}\left(\frac{\hbar}{q}\right)^2\sum_{<\vec r_1,\vec r_2>}
\frac{d\phi}{d\tau}(\vec r_1,\tau){\bf C}(\vec r_1,\vec r_2)\frac{d\phi}
{d\tau}(\vec r_2,\tau) + H_J(\{\phi\}).
\label{eq(3.2)}
\end{equation}
The boundary condition, Eq.~(\ref{eq(2.5)}), now reads
\begin{equation}
\phi(\vec r,\beta\hbar) = \phi(\vec r,0) + 2\pi m(\vec r),
\label{eq(3.3)}
\end{equation}
where the $\{m(\vec r)\}$'s are the winding numbers. We note that since
the lagrangian is invariant under the transformation
$\phi(\vec r,0)\rightarrow\phi(\vec r,0) + 2\pi l(\vec r)$
for all integers $\{l(\vec r)\}$, we can extend the limits of
integration over $\phi(\vec r,0)$ to $[-\infty,\infty]$,
the difference coming  only  from an overall multiplicative constant.
Now that the limits of Eq.~(\ref{eq(3.1)}) are all from $[-\infty,\infty]$ we
can make the following change of variable\cite{schulman}
\begin{equation}
\phi(\tau,\vec r) = \frac{2\pi}{\beta\hbar}m(\vec r)\ \tau\ +\
                    \overline{\phi}(\vec r)+\phi_f(\vec r,\tau).
\label{eq(3.4)}
\end{equation}
Here $\phi_f(\vec r,\tau)$ represents the quantum fluctuations of
the path about its mean value $\overline{\phi}(\vec r)$. These
quantum fluctuations become larger than the thermal ones as $\alpha$
increases or as the temperature decreases. This
means that we would need to take higher order harmonics in the Fourier
series into account when the quantum effects are not relatively small.
Because of  the periodicity in Eq~(\ref{eq(3.3)}),
 $\phi_f(\vec r,\tau)$  can be expanded in the Fourier series
\begin{equation}
\phi_f(\tau,\vec r) = (\beta\hbar)^{-1/2} \sum_{k=1}^{\infty}[\phi_k(\vec r)
e^{i\omega_k\tau}+C.C.],
\label{eq(3.5)}
\end{equation}
where the $\omega_k=2\pi k/\beta\hbar$ are the Bose Matsubara
frequencies. Substituting Eqs.~(\ref{eq(3.4)}) and ~(\ref{eq(3.5)}) in
Eq.~(\ref{eq(3.1)}),  expanding the Josephson term up to second order in
$\phi_k(\vec r)$, i.e. up to  order $O(q^2)$ or equivalently to $O(\alpha)$,
we obtain an effective  action for the classical
variables $\overline\phi(\vec r)$ after evaluating the integrals. In obtaining
the effective action we note that once the integrations over the
$\phi_k(\vec r)'s$ are carried out, the partition function still includes
a summation over the $m(\vec r)$'s.
In the semiclassical limit the contributions to the partition function
from configurations with $m(\vec r)$ different from zero are exponentially
small, so that we can safely take $m(\vec r)=0$ for all $\vec r$.
A most important property of the Josephson hamiltonian
$H_J$ is that it is a periodic function of its argument. This implies that
in the expansion the second order derivative with respect to the argument
in $H_J$ is proportional to $H_J$ itself. Specifically, for the cosinusoidal
form of $H_J$ we have $H_J'' = -H_J$ + constant. This important property of
$H_J$ allows us to write the effective partition function as a 2-D classical
XY model with an effective  coupling constant. The effective partition
function to this order of approximation is\footnote{ We must note that in
obtaining the effective partition function we assumed that the phases
took values between $[-\infty,\infty]$, whereas in the classical XY model
the phases are constrained to lie in the $[0,2\pi]$ range. Following this
route makes the derivation of the effective action  more direct. However,
in the small $\alpha$ regime of interest here the differences between
the two ranges for the phases can be shown to be exponentially small.}
\begin{equation}
Z_{eff} = \int \prod_{\vec r} \frac{d\overline\phi(\vec r)}
{2\pi} \exp[-\beta_{\rm eff}H_J(\overline\phi)],
\label{eq(3.6)}
\end{equation}
where the effective inverse temperature is given by
\begin{equation}
\beta_{\rm eff} = \beta -\frac{(q\beta)^2)}{12}\left[{\bf C}^{-1}(0)-
{\bf C}^{-1}(\vec d)\right].
\label{eq(3.7)}
\end{equation}
Note that we have explicitly used the fact that  ${\bf C}(\vec r_1,\vec r_2)
={\bf C}(|\vec r_1 -\vec r_2|)$, valid for a periodic lattice,
so that we can  Fourier transform the capacitance matrix.

Once we have a hamiltonian which is just like the 2-D classical XY model,
we can write down the corresponding effective RG recursion relations
to lowest order in $x$ as
\begin{eqnarray}
\frac{dK}{dl} &=& 4\pi^3 K^2 {\tilde y}^2 \frac{(1-xK)^2}{(2Kx-1)},
\label{eq(3.9)}\\
\frac{d\tilde y}{dl} &=& [2-\pi K(1-xK)]\tilde y.
\label{eq(3.10)}
\end{eqnarray}
In writing these equations we defined the variables
\begin{eqnarray}
x &=& \frac{q^2}{12E_J}\left[{\bf C}^{-1}(0)-{\bf C}^{-1}(\vec d)\right],\\
\tilde y &=& \exp[-{\frac{\pi^2}{2}K(1-xK)}]\equiv
\exp[-{\frac{\pi^2}{2}K_{eff}}].
\label{eq(3.8)}
\end{eqnarray}
The variable $x$ is the $\alpha$ parameter when the capacitance matrix is
not just the self or the mutual capacitance.
The RG equations are solved using as initial conditions
$K_{eff}(l=0)\equiv K_{eff}^0$ and $\tilde y (l=0)\equiv \tilde y_0$.
As written, the RG equations are valid for an arbitrary ratio between the
self and the mutual capacitances.
We first notice that in the $x=0$ limit the RG equations reduce to
the standard Kosterlitz RG equations\cite{bkt,jkkn}, as they should.
The form of the vortex density $\tilde y$ is
most important. As
mentioned above, $\tilde y(x=0)$ grows exponentially with temperature.
When $x\neq 0$ and as a function of temperature, $\tilde y$
exhibits a low temperature minimum. This is shown in the  discontinuous line
in Fig. 2. This behavior for $\tilde y$ is easy to understand physically.
At high temperatures the difference between $\tilde y$ and $y$ is very small.
However, below the minimum the increase in the vortex pair density is
due to nucleation of VAP via quantum phase slips.
Of course, we need to remember that we have done a perturbative
calculation in $x$ and therefore
we may not be on safe ground when $\tilde y$ starts increasing again
at low temperatures. Nevertheless, as often happens with WKB derived results,
the fact that the perturbative analysis shows a low temperature instability
is likely to be true. In fact
we also have found numerical evidence for the low temperature instability in
our QMC calculations\cite{jacobs}.

The RG equations have two nontrivial fixed points, one that
corresponds to the effective BKT thermal fluctuations driven transition,
and the other that corresponds to the
QUantum fluctuations Induced transition (QUIT)\cite{jos84,jacobs}.

The RG level curves in the ($\tilde y, K$) phase space, to lowest
order in $x$, result from solving
Eq.~(\ref{eq(3.8)}) and
\begin{equation}
\pi x K -\pi {\rm ln}K -\frac{2}{K} + 2\pi^3 {\tilde y}^2 = A,
\end{equation}
for different initial conditions. Figure 2 shows the RG flows for
different initial conditions starting with different
values for ($\tilde y_0, K^0$)
along the discontinuous line in the figure. Each RG flow line corresponds
to a different temperature with the arrows indicating the direction of
increasing $l$. We clearly see from the figure that we can divide the
temperature axis into three different regions. In the region between
$[K^{-1}_{QUIT},K^{-1}_{BKT}]$, as the value of $l$ increases we eliminate
VAP, with the unusual property that the vortex density can initially
grow for a while
before tumbling to the critical line $\tilde y=0$. This means that in this
region at $l=\infty$ there are no VAP with infinite separation, i.e.
unbounded. Below $K^{-1}_{QUIT}$, as $l$ increases, the RG trajectories grow
away from the $\tilde y=0$ line, nonmonotonically, indicating that the
perturbation expansion in $\tilde y$ is no longer valid.
If we associate the instability in the perturbation theory
with the normal state behavior, one could say then that
this behavior is characteristic of a reentrant phase transition, i.e.
going from N to SC to N. This is certainly  the case in the
high temperature regime but not necessarily so at low temperatures. The
single line that divides the two types of behaviors mentioned above is the
separatrix that determines the critical temperature. This is the line with
 the highest temperature for which we can touch the $\tilde y=0$
line. The corresponding separatrix value of the constant $A=A_c$ is determined
from the condition that it passes through
the point $(K^{-1}_{BKT},\tilde y=0$). This leads to the result that the
critical point is obtained from solving the equation,
\begin{equation}
 \pi xK_c-\pi{\rm ln}K_c-\frac{2}{K_c}
 2\pi^3 \exp\{-\pi^2 K_c(1-xK_c) \} =A_c(x),
\end{equation}
where $A_c(x)$ is obtained from  $A_c(x)=2\pi x K_c-\pi -\pi
\rm {ln} K_c$, and $K_c$ is the solution to the equation
$2=\pi x K_c(1-xK_c)$. We'll come back to the problem of determining
$K^{-1}_{BKT}$ in the next section, where we make a comparison with the
experimental results. It is known, however,  that the determination of
$T_{BKT}$ using the RG equations is not quantitatively
exact, since the RG analysis is explicitly derived  for the Villain
action. In comparing with experiment in the next section we
will take this into account. Here we  present the corrections to
the classical results to the leading order in $x$, which give the correct
qualitative trends. Specifically, expanding in powers of $x$ we find that
$T_{BKT}$ and $T_{QUIT}$ are given by
\begin{eqnarray}
T_{\rm BKT}&\approx& T_{\rm BKT}^{(0)} - \frac{E_J}{k_{\rm B}}x + O(x^2),
\label{eq(3.11)}\\
T_{\rm QUIT}&\approx& \frac{E_J}{k_{\rm B}}x+O(x^2).
\label{eq(3.12)}
\end{eqnarray}

Notice that these equations are applicable not only in 2-D,
for if the system described by Eq.~(\ref{eq(3.6)}) has a
transition point at some $K_{\rm eff}^{\rm c}$ then the equation
$K_{\rm eff}^{\rm c}~=~K-xK^2$ has two solutions for $K$, which are the ones
implied in Eq.~(\ref{eq(3.10)}). Moreover, notice that the results to the
first order in $x$ are independent of the specific value
of $T_{\rm BKT}^{(0)}$.
This means that if we consider the finite magnetic field case,
the corresponding critical temperature will be $T_c(B)\approx
T_c^{(0)}(B)-(E_J/k_{\rm B})x+O(x^2)$. Furthermore, we notice that,
to the lowest order in $x$,
the $T_{QUIT}$ must be the {\bf same} with and without a field.
This fact will be compared with the  experiments in the next section.

The explicit leading order calculation of the correction to the BKT critical
temperature in the asymptotic limits in which either the self or the mutual
capacitance dominates results in
\begin{equation}
\frac{T_{\rm BKT}}{T_{\rm BKT}^{(0)}} = \left\{
\begin{array}{ll}
    1-\frac{4}{3\pi}\frac{E_{C_{\rm s}}}
    {E_J}+O\left[\left(\frac{E_{C_{\rm s}}}
    {E_J}\right)^2\right],
    & \mbox{if $C_{\rm s}\gg
     C_{\rm m}$}\\
    1-\frac{4}{3z\pi}\frac{E_{C_{\rm m}}}
    {E_J}+O\left[\left(\frac{E_{C_{\rm m}
    }}{E_J}\right)^2\right],
    & \mbox{if $C_{\rm s}\ll
    C_{\rm m}$}
     \end{array}
     \right.
     \label{eq(3.14)}
\end{equation}
here $z$ is the coordination number of the array, and for a square array in
two dimensions, $z=4$.

\section{\bf {Comparison to experiment}}

We now move on to a brief discussion of how  we obtained the results presented
in Fig. 1. As mentioned before, in trying to find quantitative correspondence
between experiment and theory it is important to ascertain the validity
of the theoretical models employed to study the arrays. We have carried
out two different types of checks. One based on the RG analysis described
in the previous sections and the other from a nonperturbative quantum
Monte Carlo calculation\cite{jens,cristian}.
We discuss  the RG analysis here and only
briefly mention the QMC results, with more details left for
forthcoming publications \cite{jens,cristian}.

As mentioned  in section III , it is known that the RG equations
do not lead to quantitatively exact results for the critical temperatures.
They do, of course, lead to the correct universal
critical  exponents. However, what has so far been measured experimentally
is the phase boundary between the SC and the N phases. We then need  a
consistent way to compare the RG results with the experimental results.

We first note that the phase diagram of Fig. 1 is plotted as a function of
$\alpha_m=\frac{E_{C_m}}{E_J}$, since experimentally $\alpha_s$ is
three orders of magnitude smaller. We can then write $K_{eff}=K(1-xK)$
with $C_S=0$ as
\begin{equation}
K_{eff} = K - \frac{\alpha_m}{6}K^2.
\end{equation}
Next we  set the critical temperature for the classical model to be the one
obtained in classical MC simulations (e.g. \cite{ramirez})
$1/K_c^{(0)}\approx 0.93$, so that the critical temperature $T_c(\alpha_m$)
of the actual model is given by the equation
\begin{equation}
k_BT_c(\alpha_m) = \frac{1}{2K_c^{(0)}} \left[1+\sqrt{1
-\frac{2\alpha_m}{3}K_c^{(0)}}\right].
\end{equation}
To lowest order  in $\alpha_m$ this equation gives
\begin{equation}
T_c(\alpha) = T_c(0) - \frac{\alpha_m}{6},
\end{equation}
whereas the maximum value of $\alpha_m$ for which there is a
physical solution is $\alpha_m = \frac{3}{2} k_BT_c(0)
\approx 1.4$. The results
obtained from this analysis are shown as a discontinuous line in Fig. 1. By
following this approach we see that for $\alpha_m\leq 1$ the RG result
is actually quite good when compared to the experimental and the
MC results.

We also have  extended our previous QMC calculations to the case when
the off-diagonal capacitance is dominant. The results are shown in Fig. 1
by the crosses, including their error bars. It is clear from these results
that the correspondence between experiment and QMC results is excellent,
up to nonperturbative values of $\alpha_m$. This leads us to the conclusion
that the model studied here does provide a good representation of the
experimental system, at
least in the SC to N regime. In the next section we will briefly discuss
what happens in the insulating region to normal region.

The discussion presented above dealt with the SC to N phase boundary.
What about evidence for a QUIT? In this regard we note that in the
experimental results of Ref.~\cite{herre2} there are results
for a sample with a nominal
$\alpha_m=1.67$ for which there is a double type of reentrant behavior.
The low temperature glitch seen in the resistance versus temperature diagram
occurs at $T'=40mK$. Moreover the latter instability is also characterized by
an increase in noise fluctuations in the IV characteristics measured
in this sample. If we assume that what is seen at $T'$ is related to the QUIT,
using Eq.(24) and the parameter values of the experiment we get a
$T_{QUIT}=33mK$, rather close to the experimental value. Furthermore
when the same experiment is repeated in a small magnetic field the low
temperature instability is found at he same temperature, i.e.;
$T'(f=0)=T'(f=0.08)$. This result is also consistent with
Eq. (24) which also leads
to a $T_{QUIT}$ independent of $f$ at leading order in $\alpha_m$.
These results may just be coincidental and more work needs to be done
to conclusively connect the $T_{QUIT}$ with the low temperature
instability already seen in the Delft experiments.

\section{\bf {Insulating to Normal cross over}}

As mentioned in the introduction, as $\alpha_m$ increases there is a SC to I
transition at finite temperature. In the $\alpha_s=0$ case,
the insulating phase has been modelled as a two-dimensional Coulomb
gas of charges with a possible BKT charge unbinding
transition\cite{yoshi,moij2,fazio}. This
situation has been studied extensively, in particular in Ref.~\cite{fazio}.
Here we ask if there is an equivalent QUIT in the
insulating phase at low temperatures. If we draw an analogy to
the quantum induced vortices in the SC phase one could
also imagine that the the number of free dipoles in the arrays could increase
due to quantum fluctuations. However, as we show below the vortices and
charges are not dual to each other in that sense.

The calculation described here aims at finding the leading correction
to the charging hamiltonian due to Josephson junction fluctuations.
 We then expand the Josephson contribution to $Z$ as
 \begin{equation}
\exp\left[-\frac{1}{\hbar}\int_{0}^{\beta\hbar}d\tau H_J(\tau)\right]\approx
1-\frac{1}{\hbar}\int_{0}^{\beta\hbar}d\tau H_J(\tau)+\frac{1}{2\hbar^2}
\int_{0}^{\beta\hbar}\int_{0}^{\beta\hbar}d\tau d\tau' H_J(\tau)H_J(\tau')
+\dots
\label{eq(4.1)}
\end{equation}
As before, we use  Eqs.~ (12) and (13), but
this time we integrate out both $\phi_f(\tau,\vec r)$
and $\overline{\phi(\vec r)}$, which leaves us with an
effective action for the $\{m\}$'s,
\begin{equation}
Z\approx Z_\phi \prod_{\vec r} \sum_{m(\vec r)=-\infty}^{\infty}
\exp\left[-\frac{1}{4\tilde K_{\rm eff}}\sum_{\vec r,\vec d}\left(m
(\vec r+\vec d) - m(\vec r)\right)^2\right],
\label{eq(4.2)}
\end{equation}
where we have assumed that $C_{\rm s}\ll C_{\rm m}$. The function
$Z_\phi$ does not contain the $\{m\}$ variables, and the effective
coupling constant $\tilde K_{\rm eff}$ is
\begin{equation}
\tilde K^{-1}_{\rm eff} = \tilde K^{-1}\left[1+\left(\frac{2\pi^2}{\alpha_m}
\right)^2 g(\tilde K)\right],
\label{eq(4.3)}
\end{equation}
given as a function of $\tilde K=(\beta E_{C_m})/2\pi$.
As in the vortex dominated case, we have ended up with a Coulomb gas problem
but with a renormalized coupling constant. The
function $g(\tilde K)$ determines the importance of the zero point
fluctuations of the phases on the charge dominated phase.
The function $g(\tilde K)$  is defined as
\begin{equation}
g(\tilde K)=\int_{0}^{1/2} dt[1-\cos(2\pi t)]\exp\left[-\frac{(2\pi)^2}{z}
\tilde Kt(1-t)\right].
\label{eq(4.4)}
\end{equation}
For a BKT type phase transition we have
 $\tilde  {K}_{eff}=K_{\rm BKT}^{c} = 2/\pi$ \cite{fazio}. It is important
 here  to see if, within this approximation, the system shows a
QUIT or a reentrant transition. We then need  to study the number of
solutions to this equation. From the fact that $K_{\rm BKT}^{c}\sim O(1)$
we can see that for any value of $1/{\alpha_m}$ the function $g(\tilde K)$
has  only one solution for $\tilde K_c$. This means that, to this order
of approximation, {\bf there is no}
QUIT in the charge dominated phase in this model.

A question that immediately arises is: Why is there a difference between the
vortex and  charge dominated phases, in particular in view of the duality
between the two phases extensively studied in \cite{fazio}? The reason is that
the duality is not exact since there is a term in the action, obtained using
the Villain approximation, that breaks this symmetry. If one includes this
term we then see that the cost of producing quantum fluctuations in the
vortices is bounded from above whereas the corresponding
cost in the charge dominated phase is unbounded. In our calculation we
have kept these contributions  intact.

\acknowledgments
We thank  H. van der Zant for many useful discussions and for providing us
with copies of his unpublished work prior to publication.
We also thank J. Houlrik and M. Novotny
for useful conversations and suggestions in the early stages of this work.
This work has been partially supported  by $NSF$ grant DMR-9211339.

redpajama_set_name: RedPajamaArXiv


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

Photoproduction of pion pairs on the proton has been studied extensively both in theory and experiment since the 1960s. In recent years the interest in the process has been revived by improved experimental data which were obtained by using polarized tagged photon beams at the MAMI synchrotron (Mainz, Germany) \cite{Daphne}, the GRAAL facility (Grenoble, France) \cite{GRAAL}, and the Jefferson laboratory (Newport News, USA) \cite{Jefferson}. At GRAAL, the high-energy photon beam is produced by Compton backscattering of laser light on a relativistic electron beam. These studies allow insights into the internal structure and excitation spectrum of the proton. A particular focus lies on polarization asymmetry measurements which are sensitive to interference cross terms. Various theoretical models based on effective Lagrangians for the pion-nucleon interaction have been developed (see \cite{Oset} and references therein).

When protons or other particles interact directly with intense laser beams (rather than with single photons from synchrotron or Compton backscattering sources), multiphoton processes may arise involving the simultaneous absorption of more than one photon. Multiphoton $e^+e^-$ pair creation was observed in ultrarelativistic electron-laser collisions at SLAC \cite{SLAC} and theoreticians have studied related $e^+e^-$ pair creation processes by relativistic proton impact on intense laser beams \cite{Report,MVG, Kaminski,  Milstein, Kuchiev, Baur}. In this setup, the laser frequency and field strength are largely Doppler-enhanced in the projectile rest frame. The process of $e^+e^-$ pair creation induced by high-energy neutrino impact on an intense laser beam has also been examined theoretically \cite{Tinsley}.

Inspired by the sustained progress in laser technology, very recently theoreticians started to study laser-induced $\mu^+\mu^-$ pair creation in proton-laser collisions \cite{Kuchiev, Deneke}. The apparent gap between the laser photon energy and the $\sim 100$\,MeV energy scale of the process can in principle be bridged by combining upcoming X-ray laser sources ($\hbar\omega_{\rm lab}\sim 10$\,keV) with the ultrarelativistic proton beam at CERN-LHC ($\gamma\approx 7000$) \cite{LHC}. The Doppler-upshifted photon energies $\omega\approx 2\gamma\omega_{\rm lab}\sim 100$\,MeV lie in the desired range. Large-scale \cite{XFEL} as well as table-top \cite{tabletopXFEL} free-electron lasers (FELs) are currently being developed aiming at the generation of intense coherent X-ray pulses. Coherent X-rays are also envisaged via high-order harmonics from oscillating plasma surfaces \cite{surface} or atomic gas jets \cite{Klaiber}. Such compact and portable X-ray sources hold the potential to be operated in conjunction with the LHC proton beamline.

In this paper, we consider $\pi^+\pi^-$ pair creation by multiphoton absorption in ultrarelativistic proton-laser collisions, i.e. the reaction $p+n\omega\longrightarrow p+\pi^+ +\pi^-$ with the photon number $n>1$. To this end, we combine the well-established approach to multiphoton processes in quantum electrodynamics (QED) \cite{strongfieldQED} with a simple phenomenological model to describe the pion-nucleon interaction. In general terms, the present study may be considered a first step towards an extension of the theory of laser-dressed QED into the realm of hadronic physics. More specifically, our main goals are (i) to provide order-of-magnitude estimates for $\pi^+\pi^-$ multiphoto-production rates, demonstrating the observability of the process; (ii) to compare with the corresponding rates for $\mu^+\mu^-$ production and show that a range of laser frequencies exists where $\pi^+\pi^-$ production dominates over the (direct) production of muons; and (iii) to discuss prospects why detailed investigations of multiphoto-production of pions might be useful.

We note that highly energetic reactions can also be induced when high-power laser pulses interact with solid targets. In the resulting plasma wakefields, electrons are accelerated to relativistic energies and emit secondary bremstrahlung $\gamma$-rays. These have led to the observation of photonuclear reactions \cite{LaserNuclear} and efficient $e^+e^-$ pair creation through the Bethe-Heitler effect \cite{Chen}. The setup also offers prospects for Bethe-Heitler creation of muon pairs \cite{Kaempfer} and single pion photoproduction through the reaction $\gamma p\to \pi^+n$ via the $\Delta(1232)$ resonance \cite{Karsch}. We point out that all these processes rely on single photo-absorption and do not exhibit multiphoton character. Particle reactions such as $\mu^+\mu^-$ and $\pi^+\pi^-$ production were also considered in an $e^+e^-$ plasma coupled to a photon field \cite{muonPLB}. Single pion production in strong magnetic fields \cite{ADP} and in collisions of laser-accelerated protons with nuclei have been studied theoretically as well \cite{pion}.

In the following, natural units with $\hbar = c = \epsilon_0 = 1$ are used.
\section{QED description of multiphoton muon pair creation}
\label{muon_creation}

We briefly review the laser-dressed QED approach to muon pair creation in the combined fields of an electromagnetic wave and an atomic nucleus \cite{Report,MVG,Kaminski, Kuchiev, Deneke}. Within the external-field approximation of QED, the total Lagrangian of the problem reads
\begin{eqnarray}
\label{L_QED}
{\mathcal L}={\overline \Psi}[i\gamma_\mu (\partial^\mu + ieA^\mu_L + ieA^\mu_C)-m]\Psi\,,
\end{eqnarray}
where $A_L^\mu$ and $A_C^\mu$ are the four-potentials of the laser wave and the nucleus, respectively, and $m$ is the muon mass. In the spirit of the Furry picture, one may split the total Lagrangian ${\mathcal L}={\mathcal L}_V+{\mathcal L}_C$ into an unperturbed part ${\mathcal L}_V={\overline \Psi}(i\gamma_\mu D^\mu - m)\Psi$ and the remaining interaction ${\mathcal L}_C=e{\overline \Psi}\gamma_\mu A^\mu_C\Psi$. Here, $D^{\mu}=\partial^{\mu}+ieA_L^{\mu}$ denotes the covariant derivative with respect to the laser field. The field theory of ${\mathcal L}_V$ is solved exactly by the so-called Dirac-Volkov states $\Psi_{p_\pm, s_\pm}$ which include the interaction of the leptons with the plane-wave laser field up to all orders~\cite{Itzykson}. The leptons are characterized by their four-momenta $p_\pm^\mu$ and spin projections $s_\pm$ outside the laser field.

The Dirac-Volkov solutions may be used as basis states in perturbative calculations with respect to the remaining interaction with the nuclear field. As a result, the leading-order $S$-matrix element for multiphoton muon pair creation on a proton at rest reads
\begin{eqnarray}
\label{Smuon}
S_{\mu^+\mu^-} = -i\int d^4x\, {\overline \Psi}_{p_-,s_-}(x) \gamma^0 \Psi_{p_+,s_+}(x) V_C(x)\,,
\end{eqnarray}
where
\begin{eqnarray}
\label{Coulomb}
V_C(x) = eA^0_C(x) = \dfrac{e^2}{4\pi |{\boldsymbol x}|}
\end{eqnarray}
denotes the Coulomb potential energy in this frame and $\gamma^0$ is a Dirac matrix.

To be specific, we assume a monochromatic laser wave of circular polarization. Choosing the propagation direction along the $z$ axis, the corresponding four-potential in the radiation gauge reads
\begin{equation}
\label{laser_pot}
A_L^{\mu}(x)= (0,{\boldsymbol A}_L)\,,
{\boldsymbol A}_L(x)=A_0[{\bf e}_x \cos(kx)+{\bf e}_y \sin(kx)]
\end{equation}
where $k^\mu=\omega(1,0,0,1)$ is the wave four-vector and $(kx)=\omega t-{\boldsymbol k}\cdot{\boldsymbol x}$ is the laser phase. Although the laser field is treated as a classical electromagnetic wave, photons arise from a mode expansion of the oscillatory parts in Eq.\,\eqref{Smuon} which contains multiphoton processes of arbitrary order \cite{Report,MVG, Kaminski, Kuchiev, Deneke}.

From the $S$-matrix, the $\mu^+\mu^-$ production rate is obtained as
\begin{eqnarray}
\label{R}
R_{\mu^+\mu^-} = \int \frac{d^3{\boldsymbol p}_+}{(2\pi)^3} \int \frac{d^3{\boldsymbol p}_-}{(2\pi)^3}\sum_{s_+,s_-} |S_{\mu^+\mu^-}|^2.
\end{eqnarray}
The cross section of the process may be obtained by dividing out the photon flux $j=\omega A_0^2/4\pi$. As an alternative to the $S$-matrix approach, the laser-dressed polarization operator can be employed to calculate total pair creation rates \cite{Milstein}.

A few additional remarks are in order here. (1) For a proper description of muon pair creation in the field of a heavy nucleus, the finite nuclear extension needs to be taken into account in general \cite{Deneke}. Here we restrict our consideration to protons as projectiles, which may be treated as pointlike [see Eq.\,\eqref{Coulomb}] to a good approximation. (2) Being interested in order-of-magnitude estimates for the production rates, the recoil suffered by the proton during the process is ignored for simplicity. (3) We will consider X-ray laser fields where the value of the Lorentz-invariant parameter $\xi\equiv eA_0/m$ is much smaller than unity. In this regime of laser-matter coupling, a process involving $n$ photons could, in principle, be calculated within $n$-th order of perturbation theory in the photon field. We find it more convenient, however, to work within the framework of laser-dressed QED.
\section{Effective hadronic model of multiphoton pion pair creation}
\label{pion_creation}
Various models of $\pi^+\pi^-$ photoproduction on the nucleon have been developed which are based on effective Lagrangians for the pion-nucleon interaction \cite{Oset}. It was found that a proper description of the process requires the inclusion of more than 40 Feynman graphs. They consist of one-point, two-point and three-point diagrams with respect to the number of hadronic vertices. In the present case of interest, where the pions are created by the absorption of more than one photon, the number of relevant Feynman diagrams increases tremendously, which renders a consideration on this level almost prohibitorily involved. However, order-of-magnitude estimates for $\pi^+\pi^-$ multiphoto-production rates may be obtained by applying the following simplified model.

We treat the pion as a scalar particle whose dynamics is governed by the covariant Lagrangian
\begin{eqnarray}
{\mathcal L}= \left(D^{\mu}\phi\right)^*\left(D^{\nu}\phi\right)\left[\delta_{\mu\nu}-\mathcal{K}_{\mu\nu}\right]-m_\pi^2\phi\phi^*+{\mathcal L}_{\rm had},
\label{ld}
\end{eqnarray}
where $\phi$ denotes the scalar pion field, $m_\pi$ the pion mass, and $\delta_{\mu\nu}$ the Kronecker symbol. The symmetric contribution
\begin{eqnarray}
\mathcal{K}_{\mu\nu}=-\dfrac{\lambda_e+\lambda_m}{m_\pi} A_0^2 k_\mu k_\nu
\end{eqnarray}
represents a measure for the strength of the induced electric and magnetic dipole moments of the pion by the external photon field \cite{Kruglov}. It characterizes, to some degree, their underlying structure from quarks. The experimental value $\lambda_e+\lambda_m\approx0.16\times 10^{-4}~\mathrm{fm^3}$ of the sum of the electric and magnetic pion polarizabilities  was recently measured in Primakoff scattering of high-energy pions at CERN-SPS~\cite{ref:Abbon}.

The electromagnetic coupling of the pion to the laser field is contained in the unperturbed part, ${\mathcal L}_0\equiv{\mathcal L}-{\mathcal L}_{\rm had}$, of the Lagrangian in Eq.~\eqref{ld}. Similarly as in Sec.II, the field theory of ${\mathcal L}_0$ can be solved exactly. The corresponding scalar states $\phi_p$ were derived in \cite{Kruglov}. They represent generalizations  to composite particles of the usual scalar Volkov solutions of the Klein-Gordon equation in the presence of an electromagnetic wave. In fact, all terms in Eq.~\eqref{ld} multiplied by the Kronecker symbol contribute to the usual Lagrangian of a pointlike spinless boson in an external laser field, which is solved by these Gordon-Volkov states \cite{Report}.

The term ${\mathcal L}_{\rm had}$ in Eq.~\eqref{ld} describes the hadronic interaction between the pion and the proton. For the sake of simplicity of our model, we assume a scalar Yukawa coupling of the form
\begin{eqnarray}
\label{g}
{\mathcal L}_{\rm had} = g{\overline\psi}\psi\phi\,,
\end{eqnarray}
where $\psi$ denotes the nucleon field. We consider the coupling constant $g$ as a free parameter of our model which will be adjusted to experimental data. Note that in reality the pion-nucleon coupling is more involved and consists of pseudoscalar and pseudovector interactions of Yukawa type \cite{Walecka}. Nevertheless, the simple form in Eq.~\eqref{g} already incorporates main features of the hadronic interaction.

The leading-order $S$-matrix element for $\pi^+\pi^-$ multiphoto-production on the proton in the presence of an intense laser field may be written as
\begin{eqnarray}
\label{Spion}
S_{\pi^+\pi^-} = -i\int d^4x\, \left(\phi_{p_-}^*\partial_t\phi_{p_+}-\phi_{p_+}\partial_t\phi_{p_-}^*\right) V_{\rm had}(\vec{x})
\end{eqnarray}
with the Yukawa potential
\begin{eqnarray}
\label{Yukawa}
V_{\rm had}(x)=\frac{g^2}{4\pi}
\frac{{\rm e}^{-m_\pi |{\boldsymbol x}|}}{|{\boldsymbol x}|},
\end{eqnarray}
arising from the interaction \eqref{g}. Formally, this potential is associated with a scattering-type Feynman diagram of second order in the hadronic coupling constant. In the spirit of a meson-exchange model, it reflects the short-range nature of the hadronic interaction which represents an essential difference with the long-range Coulomb potential \eqref{Coulomb}. The total $\pi^+\pi^-$ production rate, $R$, is obtained from the $S$-matrix as in Eq.\,\eqref{R}, with the spin sum omitted.

We determine the numerical value of the coupling constant $g$ appearing in Eq.~\eqref{Yukawa} by taking reference to the available experimental data on $\pi^+\pi^-$ production by single-photon absorption on the proton. By setting $g=4.3$, we reproduce the $\pi^+\pi^-$ photoproduction cross section $\sigma = R/j$ at 450 MeV as measured at MAMI \cite{Daphne}. With this choice, we find reasonable agreement with the experimental data in the energy range from threshold up to $\approx 480$\,MeV (see Fig.~\ref{TC}).

\begin{figure}[h]
\centering
\includegraphics*[width=20pc]{Dadi_Fig1.ps}
\caption{(Color online) Total cross sections for $\pi^+\pi^-$ production by single-photon absorption on the proton. Shown are the experimental data recorded with the DAPHNE detector at MAMI \cite{Daphne} and the results from the present phenomenological model. The coupling constant for the pion-nucleon interaction was set to $g=4.3$ in order to reproduce the measured cross section at 450 MeV.}
\label{TC}
\end{figure}

Due to the adjustment of $g$ to a measured cross section, the Yukawa potential in Eq.~\eqref{Yukawa} becomes an effective potential which mimics the contributions from all relevant Feynman graphs in an approximative way. Our phenomenolgical model thus effectively accounts for various effects which are not included explicitely. These comprise, in particular, contributions from nucleonic resonances such as the $\Delta(1232)$.
Note that the value of the coupling constant extracted from measurements of pion-nucleon scattering is significantly larger, $g_{\pi N}\approx 13$ \cite{Arndt}, which clearly demonstrates the phenomenological character of our model.

At higher photon energies, our predictions nevertheless would deviate significantly from the experimental data in \cite{Daphne}. In fact, within our model the cross section raises approximately like a power law, $\sigma\sim (\omega-2m_\pi)^2$, and thus cannot reproduce the saturation occuring above 700 MeV observed in experiment. Such high energies, however, are not crucial in view of $\pi^+\pi^-$ production by multiphoton absorption. Namely, when a pion pair is created by the absorption of two photons with a total energy of 700 MeV, the energy of each single photon already lies above the $\pi^+\pi^-$ threshold. Since the probability for an $n$-photon process generally scales like $\xi^{2n}$ in the perturbative coupling regime ($\xi\ll 1$), the two-photon process is strongly suppressed as compared with the dominant single-photon process. Hence, for our purposes the good agreement in the energy range up to 480 MeV may be considered sufficient to obtain order-of-magnitude estimates for $\pi^+\pi^-$ multiphoto-production.

\section{Results on multiphoton $\pi^+\pi^-$ production}
Based on Eq.\,\eqref{Spion}, we have calculated total rates and angular spectra for multiphoto-production of $\pi^+\pi^-$ pairs in ultrarelativistic proton-laser collisions. The value $g=4.3$ was used for the pion-nucleon coupling constant in the effective Yukawa potential~\eqref{Yukawa}. The threshold for $\pi^+\pi^-$ production by two-photon absorption lies at an X-ray laser frequency of $\omega_{\rm lab}\approx m_\pi/2\gamma\approx 10$\,keV in the laboratory frame, when the proton beam is counterpropagating the laser beam with a Lorentz factor of $\gamma=7000$. This frequency domain is aspired by the currently emerging X-ray laser facilities \cite{XFEL,tabletopXFEL,surface,Klaiber}. Below, we mainly concentrate on processes involving two photons.

Figure~\ref{TR} shows our results on the total production rates as a function of the photon energy, refering to the rest frame of the proton. In the laboratory frame, the rates are reduced to $R_{\rm lab}=R/\gamma$ due to time dilation.

Let us consider an example. For a photon energy of 200\,MeV in the proton frame, the $\pi^+\pi^-$ production rate by two-photon absorption amounts to $R\sim 10^3$\,s$^{-1}$ (see Fig.~\ref{TR}). Hence, in the collision of a proton beam containing $10^{11}$ particles \cite{LHC} with an X-ray pulse of 100\,fs duration \cite{XFEL}, the probability for production of one pion pair is $10^{-3}$, assuming perfect beam overlap.
An average $\pi^+\pi^-$ production rate of 10 events per second is obtained by taking the
envisaged X-ray pulse repetition rate of 10\,kHz into account. Note that the latter can be synchronized with the revolution frequency of the LHC proton beam.

We note that the effect of the electric and magnetic pion polarizabilities turns out to be immaterial for the present parameters. The reason is that the parameter $(\lambda_e+\lambda_m)(p_\pm k)^2/m_\pi\sim (\lambda_e+\lambda_m)m_\pi^3\sim 10^{-5}$, which effectively quantifies the coupling of the external laser field to the pion polarizabilities \cite{Kruglov}, is much smaller than its coupling $e^2\approx 0.1$ to the pion charge.

In order to detect the process of $\pi^+\pi^-$ production by two-photon absorption unambigously in experiment, it is necessary to discriminate it from the sequential process where first a single $\pi^+$ and a neutron are created, $\gamma p \to \pi^+ n$, followed by the reaction $\gamma n \to \pi^- p$. For the purpose of discrimination, the angular distribution of the pions might serve as a fingerprint (see below).

It is interesting to compare the $\pi^+\pi^-$ production rate by two-photon absorption with the corresponding rate for $\mu^+\mu^-$ production \cite{Deneke}. Although the muon mass is substantially smaller than the pion mass, $\pi^+\pi^-$ production is dominant in the frequency range between $\approx 150$--210\,MeV, as Figure~\ref{TR} shows. For example, at a photon energy of 200\,MeV, the production of pions exceeds the production of muons by three orders of magnitude. The dominance of $\pi^+\pi^-$ production is due to the much larger strength of the hadronic interaction. The purely electromagnetic contribution to $\pi^+\pi^-$ production [resulting from the Coulomb potential of the proton in Eq.\,\eqref{Coulomb}] is suppressed by four orders of magnitude and, thus, smaller than the $\mu^+\mu^-$ rate. For higher photon energies ($\omega>2m$), the channel of $\mu^+\mu^-$ production by a single photon opens and represents the dominant process.

Our results imply that in the frequency interval $150\,{\rm MeV}\lesssim\omega\lesssim 210$\,MeV muons are predominantly produced indirectly via two-photon $\pi^+\pi^-$ production and subsequent pion decay, $\pi^+\to\mu^+ + \nu_\mu$ and $\pi^-\to\mu^- + {\bar \nu}_\mu$.

\begin{figure}[h]
\centering
\includegraphics*[width=20pc]{Dadi_Fig2.ps}
\caption{(Color online) Total rates for $\pi^+\pi^-$ and $\mu^+\mu^-$ pair production on the proton by one- and two-photon absorption, as functions of the photon energy. The amplitude of the laser vector potential is $eA_0=5.1\,$keV, corresponding to lab-frame photon intensities $I_{\rm lab}=\omega_{\rm lab}^2 A_0^2/4\pi$ of the order of 10$^{22}$\,W/cm$^2$.}
\label{TR}
\end{figure}

Figure~\ref{DS} displays our results on the angular distribution $dR/d\cos\theta$ of one of the created pions by one-, two- and three-photon absorption. The respective photon energies were chosen such that the energy absorbed in total is always 360\,MeV. The polar emission angle $\theta$ is measured with respect to the incident laser wave vector $\boldsymbol{k}$.

The calculated angular spectrum for $\pi^+\pi^-$ production by one photon attains a maximum at $\theta\approx 29^{\circ}$. It is in reasonable agreement with experimental observation where a maximum around $\approx 30^{\circ}$ was found in the integrated energy interval from 350 to 500\,MeV \cite{spectrum}. For two- and three-photon absorption the position of the maximum is shifted towards larger angles ($\theta\approx 35^{\circ}$ and $\theta\approx 38^{\circ}$, respectively) and the distributions become slightly more narrow. These are characteristic multiphoton effects which have also been found in $e^+e^-$ pair creation by few-photon absorption \cite{Sarah}.

We draw again a comparison with two-photon production of $\mu^+\mu^-$ pairs. As Figure~\ref{DS} shows, here the angular distribution is peaked at smaller angles around $\theta_+\approx25^{\circ}$ and the contribution from large angles is suppressed. The differences are caused by the long-range nature of the Coulomb potential \eqref{Coulomb} as compared to the short-ranged Yukawa potential \eqref{Yukawa}. In the latter case, the particles are mainly produced at short distances which gives rise to large momentum transfers ${\boldsymbol q}$ corresponding to large emission angles. For the muons, small momentum transfers are relatively more important since the Fourier transform of the Coulomb potential is proportional to $1/{\boldsymbol q}^2$, as compared to $1/({\boldsymbol q}^2+m_\pi^2)$ for the Yukawa potential.

\begin{figure}[h]
\centering
\includegraphics*[width=20pc]{Dadi_Fig3.ps}
\caption{(Color online) Angular distributions of one of the particles in $\pi^+\pi^-$ and $\mu^+\mu^-$ photoproduction on the proton. The spectra refer to the proton frame and the angle is measured with respect to the laser beam direction. In all cases, an identical amount of energy is absorbed from the laser field ($n\omega=360$\,MeV). The spectra have been normalized to the same height in order to facilitate their comparison.}
\label{DS}
\end{figure}

\section{Conclusion and Outlook}
A simple phenomenological model of $\pi^+\pi^-$ production by multiphoton absorption on the proton has been presented. The electromagnetic interaction of the pions with the photon field was incorporated into scalar Volkov states, while an effective Yukawa potential was used to describe the hadronic pion-nucleon coupling. The model contains one free parameter which has been adjusted to a measured cross section of $\pi^+\pi^-$ single-photon production at 450 MeV. Multiphoto-production of pion pairs could, in principle, be realized by utilizing the 7 TeV proton beam at CERN-LHC in conjunction with a compact 10--15\,keV X-ray laser source.

Predictions for total rates of two-photon $\pi^+\pi^-$ production were provided. In particular, it was shown that the two-photon production of $\pi^+\pi^-$ pairs substantially dominates over the (direct) production of muon pairs for photon energies above 150 MeV up to the threshold for $\mu^+\mu^-$ single-photoproduction. In this range of photon energies, muons are thus predominantly created indirectly via pion decays. We note that a similar conclusion was drawn in \cite{Kaempfer} with respect to pion and muon production by single bremsstrahlung photons in a laser-generated plasma.

Angular pion spectra have also been calculated which reflect the short-range nature of the hadronic interaction. When the $\pi^+\pi^-$ pair is created by more than one photon, larger emission angles arise which correspond to larger momentum transfers in the process.

As an outlook, we briefly discuss some general issues for which the process of multiphoto-production of pions might prove to be especially useful.

(1) Due to their long wavelength, photons interact uniformly over the nuclear volume and can therefore provide quantitative information on the entire structure of a given nucleus. Multiphoton processes such as $\pi^+\pi^-$ production might be particularly sensitive to the inner nuclear structure since they involve larger momentum transfers. The enhanced contribution from large momentum transfers moreover implies that multiphoto-production of pions probes a different spatial region inside the nucleus. This might help to gain an improved understanding of the short-range repulsion and intermediate-range attraction in the nucleon-nucleon potential.

(2) A main focus of pion photo-production is presently lying on polarization observables \cite{GRAAL,Jefferson}. Pion pair production via multiphoton absorption from a circularly polarized laser beam could be of particular interest in this regard. Since all photons carry the same helicity, the simultaneous absorption of several of them might lead to characteristic signatures on the polarization properties of the produced pions.

(3) Finally, photoproduction of pions exhibits an interesting analogy in atomic physics. A formal relation between the strong-field processes of $e^+e^-$ pair production and atomic photoionization in intense laser fields has already been revealed (see, e.g., \cite{MVG}). In the same spirit, a connection between photoionization of atoms and photoproduction of pions may be established, taking into account that pions can be produced both as single particles and in pairs. An analogy with single and double ionization of helium by one- and two-photon absorption is therefore suggestive. These four different ionization processes have found a sustained interest in atomic physics because of their sensitivity to electron-electron correlations \cite{helium}. While the production of single pions and pion pairs by one photon has been examined for a long time, we presented here a first study of pion pair production by two-photon absorption. The picture consisting of all four channels may be completed by a consideration of single pion production by two-photon absorption.

We hope that our study may stimulate further research on multiphoto-production of pions on the nucleon.

\section*{Acknowledgements}
We thank B. Kopeliovich, H. J. Pirner, and B. Povh for helpful discussions on various aspects of hadronic interactions.


text:
\section{Non-Fermi liquid behavior}
\label{sec_nfl}

\begin{figure}[t]
\begin{minipage}{6cm}
\vspace*{1cm}
\includegraphics[width=5cm]{sdl1.eps}
\end{minipage}\begin{minipage}{9cm}
\includegraphics[width=9cm]{sdl2p.eps}
\end{minipage}
\caption{One loop corrections to the fermion propagator and the
quark gluon vertex in dense QCD. The solid squares denote hard
dense loop insertions.}
\label{fig_sdl}
\end{figure}

 At high baryon density the relevant degrees of freedom are
particle and hole excitations which move with the Fermi
velocity $v$. Since the momentum $p\sim v\mu$ is large,
typical soft scatterings cannot change the momentum by very
much and the velocity is approximately conserved. An effective
field theory of particles and holes in QCD is given by
\cite{Hong:2000tn}
\be
\label{l_hdet}
{\cal L} =\psi_{v}^\dagger \left(iv\cdot D
   - \frac{1}{2p_F}D_\perp^2 \right) \psi_{v}
   + {\cal L}_{HDL}
   -\frac{1}{4}G^a_{\mu\nu} G^a_{\mu\nu}+ \ldots ,
\ee
where $v_\mu=(1,\vec{v})$. The field $\psi_v$ describes particles
and holes with momenta $p=\mu(0,\vec{v})+l$, where $l\ll\mu$. We will
write $l=l_0+l_{\|}+l_\perp$ with $\vec{l}_{\|}=\vec{v}(\vec{l}
\cdot \vec{v})$ and $\vec{l}_\perp = \vec{l}-\vec{l}_{\|}$.
At energies below the screening scale $g\mu$ hard dense loops
have to be resummed. The generating functional for hard dense loops
in gluon $n$-point functions is given by \cite{Braaten:1991gm}
\be
\label{S_hdl}
{\cal L}_{HDL} = -\frac{m^2}{2}\sum_v \,G^a_{\mu \alpha}
  \frac{v^\alpha v^\beta}{(v\cdot D)^2}
G^b_{\mu\beta},
\ee
where $m^2=N_f g^2\mu^2/(4\pi^2)$ is the dynamical gluon mass and
the sum over patches corresponds to an average over the direction
of $\vec{v}$.

 The hard dense loop action describes static screening of electric
fields and dynamic screening of magnetic modes. Since there is no
screening of static magnetic fields low energy gluon exchanges are
dominated by magnetic modes. The resummed transverse gauge boson
propagator is given by
\be
\label{d_trans}
D_{ij}(k) = \frac{\delta_{ij}-\hat{k}_i\hat{k}_j}{k_0^2-\vec{k}^2+
i\eta |k_0|/|\vec{k}|} ,
\ee
where $\eta=\frac{\pi}{2}m^2$ and we have assumed that $|k_0|<|\vec{k}|$.
We observe that the gluon propagator becomes large in the regime
$|\vec{k}|\sim (\eta k_0)^{1/3}\gg k_0$. This leads to an unusual scaling
behavior of Green functions in the low energy limit. Consider a
generic Feynman diagram and scale all energies by a factor $s$.
Because of the functional form of the gluon propagator in the
Landau damped regime gluon momenta scale as $|\vec{k}|\sim s^{1/3}$.
This implies that the gluon momenta are much larger than the
gluon energies. The quark dispersion relation is $k_0\simeq
k_{||}+k_\perp^2/(2p_F)$. The only way a quark can emit a very
spacelike gluon and remain close to the Fermi surface is if the
gluon momentum is transverse to the Fermi velocity. We find
\be
k_0 \sim s, \hspace{0.5cm}
k_{||}\sim s^{2/3},\hspace{0.5cm}
k_\perp \sim s^{1/3},
\ee
and $k_0\ll k_{||}\ll k_\perp$. These scaling relations have many
interesting consequences. As an example we consider the fermion self
energy, see Fig.~\ref{fig_sdl}. The one-loop self energy is
\be
\Sigma(p) \sim \int dk_0
\int dk_\perp^2 \ \frac{k_\perp}{k_\perp^2+i\eta k_0}
 \int dk_{||} \  \frac{\Theta(p_0+k_0)}{k_{||}+p_{||}
          -(k_\perp^2+p_\perp^2)/(2p_F)+i\epsilon}
 \sim p_0\log(p_0).
\ee
A more careful calculation gives \cite{Ipp:2003cj}
\be
\label{sig}
\Sigma(\omega) = \frac{g^2}{9\pi^2} \left[
  \omega \log\left( \frac{4\sqrt{2}m}{\pi|\omega|} \right)
  +\omega + i\frac{\pi}{2}|\omega| \right].
\ee
There are no corrections of the form $g^{2n}\log^n(\omega)$
\cite{Schafer:2004zf}. Higher order corrections involve fractional
powers $(\omega/m)^{1/3}$. Eq.~(\ref{sig}) shows that cold quark matter
is not a Fermi liquid. The Fermi velocity vanishes on the Fermi surface
and the specific heat scales as $T\log(T)$.

\section{CFL Phase}
\label{sec_cfl}

 If the temperature is sufficiently low quark matter is
expected to become a color superconductor and quarks acquire
a gap due to pairing near the Fermi surface. In the following
we will concentrate on the color flavor locked (CFL) phase
which is the ground state of three flavor quark matter at very
high baryon density. Our starting point is the effective theory
of the CFL phase derived in \cite{Casalbuoni:1999wu,Kryjevski:2004jw}.
The effective lagrangian contains Goldstone boson fields $\Sigma$
and baryon fields $N$. The meson fields arise from chiral symmetry
breaking in the CFL phase. The leading terms in the effective theory
is
\be
\label{l_mes}
{\cal L} =  \frac{f_\pi^2}{4} {\rm Tr}\left(
  \nabla_0 \Sigma \nabla_0 \Sigma^{\dagger} - v_\pi^2
  \vec{\nabla} \Sigma \vec{\nabla}\Sigma^\dagger \right).
\ee
Baryon fields originate from quark-hadron complementarity. The
effective lagrangian is
\bea
\label{l_bar}
 {\cal L} &=&
 {\rm Tr}\left(N^\dagger iv^\mu D_\mu N\right)
 - D{\rm Tr} \left(N^\dagger v^\mu\gamma_5
               \left\{ {\cal A}_\mu,N\right\}\right)
 - F{\rm Tr} \left(N^\dagger v^\mu\gamma_5
               \left[ {\cal A}_\mu,N\right]\right)
  \nonumber \\
 & &  \mbox{} + \frac{\Delta}{2} \left\{
     \left( {\rm Tr}\left(N_LN_L \right)
   - \left[ {\rm Tr}\left(N_L\right)\right]^2 \right)
   - \left( {\rm Tr} \left(N_RN_R \right)
   - \left[ {\rm Tr}\left(N_R\right)\right]^2 \right)
     + h. c.  \right\},
\eea
where $N_{L,R}$ are left and right handed baryon fields in the adjoint
representation of flavor $SU(3)$, $v^\mu=(1,\vec{v})$ is the Fermi
velocity, and $\Delta$ is the superfluid gap. We can think of $N$
as being composed of a quark and a diquark field, $N_L\sim q_L
\langle q_Lq_L\rangle$. The interaction of the baryon field with
the Goldstone bosons is dictated by chiral symmetry. The covariant
derivative is given by $D_\mu N=\partial_\mu N +i[{\cal V}_\mu,N]$.
The vector and axial-vector currents are
\be
 {\cal V}_\mu = -\frac{i}{2}\left\{
  \xi \partial_\mu\xi^\dagger +  \xi^\dagger \partial_\mu \xi
  \right\}, \hspace{1cm}
{\cal A}_\mu = -\frac{i}{2} \xi\left(\partial_\mu
    \Sigma^\dagger\right) \xi ,
\ee
where $\xi$ is defined by $\xi^2=\Sigma$. The low energy constants
$f_\pi,v_\pi,D,F$ can be calculated in perturbative QCD. Symmetry
arguments can be used to determine the leading mass terms in the
effective lagrangian. Bedaque and Sch\"afer observed that $X_L=
MM^\dagger/(2p_F)$ and $X_R=M^\dagger M/(2p_F)$ act as effective
chemical potentials and enter the theory like the temporal components
of left and right handed flavor gauge fields \cite{Bedaque:2001je}.
We can make the effective lagrangian invariant under this symmetry
by introducing the covariant derivatives
\bea
\label{V_X}
 D_0N &=& \partial_0 N+i[\Gamma_0,N], \hspace{0.5cm}
 \Gamma_0 = -\frac{i}{2}\left\{
  \xi \left(\partial_0+ iX_R\right)\xi^\dagger +
  \xi^\dagger \left(\partial_0+iX_L\right) \xi
  \right\}, \\
\nabla_0\Sigma &=& \partial_0\Sigma+iX_L\Sigma-i\Sigma X_R.
\eea
Using eqs.~(\ref{l_bar}-\ref{V_X}) we can calculate the dependence
of the gap in the fermion spectrum on the strange quark mass. For
$m_s=0$ there are 8 quasi-particles with gap $\Delta$ and one
quasi-particle with gap $2\Delta$. As $m_s$ increases some of
the gaps decrease. The gap of the lowest mode is approximately
given by $\Delta=\Delta_0-3\mu_s/4$ where $\mu_s=m_s^2/(2p_F)$
and $\Delta_0$ is the gap in the chiral limit.

\begin{figure}
\begin{center}
\includegraphics[width=10cm]{kstar_mass_2_txt.eps}
\vspace*{-1cm}
\end{center}
\caption{Screening mass in the CFL phase as a function of the
effective chemical potential $\mu_s=m_s^2/(2p_F)$. The screening
mass is defined as the second derivative of the energy density
with respect to an isospin or hypercharge current. The blue
and black curves show the result in the gapless CFL phase with
(pKCFL) and without (gCFL) a supercurrent.  }
\label{fig_kst}
\end{figure}

 For $\mu_s>4\Delta_0/3$ the system contains gapless fermions
interacting with light or even massless Goldstone bosons. This
situation is superficially similar to the normal phase discussed
in Sect.~\ref{sec_nfl}, but this is not the case. The gapless
CFL phase is unstable with respect to the formation of non-zero
currents \cite{Huang:2004bg}. This can be seen by considering
the dispersion relation in the presence of a hypercharge or
isospin current. The dispersion relation of the lowest mode
is \cite{Kryjevski:2005qq,Schafer:2005ym}
\be
\label{disp_ax}
\omega_l = \Delta_0 +\frac{l^2}{2\Delta_0}-\frac{3}{4}
  \mu_s -\frac{1}{4}\vec{v}\cdot\vec{\jmath}_K,
\ee
where $l$ is the momentum relative to the Fermi surface and
$\jmath_K$ is the current. The energy relative to the CFL phase
is the kinetic energy of the current plus the energy of occupied
gapless modes
\be
\label{efct}
{\cal E} = \frac{1}{2}v_\pi^2f_\pi^2\vec{\jmath}_K^2
 + \frac{\mu^2}{\pi^2}\int dl \int
 \frac{d\Omega}{4\pi} \;\omega_l \theta(-\omega_l) .
\ee
The energy functional can develop a minimum at non-zero $\jmath_K$
because the current lowers the energy of the fermions near one of
the pole caps on the Fermi surface. Introducing the dimensionless
variables $x=\jmath_K/(a\Delta)$ and $h=(3\mu_s-4\Delta)/(a\Delta)$
we can write
\be
\label{efct_x}
{\cal E} = Cf_h(x), \hspace{0.5cm}
 f_h(x) = x^2-\frac{1}{x}\left[
   (h+x)^{5/2}\Theta(h+x) - (h-x)^{5/2}\Theta(h-x) \right] ,
\ee
where $C$ and $a$ are numerical constants. The functional eq.~(\ref{efct_x})
was analyzed in \cite{Son:2005qx,Kryjevski:2005qq,Schafer:2005ym}.
There is a critical chemical potential $\mu_s= (4/3+ah_{crit}/3)\Delta$
above which the groundstate contains a non-zero supercurrent $\jmath_K$.
This current is canceled by a backflow of gapless fermions. The
screening mass $m_V^2=(\partial^2{\mathcal E})/(\partial\jmath_K^2)$
is shown in Fig.~\ref{fig_kst}. Without the supercurent an instability
occurs for $\mu_s= 4\Delta/3$, but the instability is resolved by
a non-zero current. The new phase is analogous to $p$-wave pion
condensates at lower densities because the current is carried by
Goldstone kaons in the CFL phase and the instability is caused
by the $p$-wave interaction between kaons and fermions.