Meissner response of anisotropic superconductors
V. G. Kogan
The response field of a half-space anisotropic superconductor is evaluated for an arbitrary weak external field source. Example sources of a point magnetic moment and a circular current are considered in detail. For the penetration depth $\lambda \ll L$ with $L$ being any other relevant distance (the source size, or the distance between the source and the superconductor), the major contribution to the response is the $\lambda$ independent field of the source image. It is shown that the absolute value of $\lambda$ cannot be extracted from the response field with a better accuracy than that for the source position. Similar problems are considered for thin films. %Low frequency polarized electromagnetic waves are %reflected by anisotropic superconductors being elliptically polarized.
cond-mat0303381.tex
\documentclass[twocolumn,showpacs,amsmath,amssymb,prb,graphics,graphicx
]{revtex4}
\usepackage{bm}
\usepackage{graphicx}
\usepackage{dcolumn}
\begin{document}
%\date{Dec.1999 - Aug.2000 - \today}
\title
{Meissner response of anisotropic superconductors}
\author{V. G. Kogan}
\affiliation{Ames Laboratory DOE and Physics Department ISU, Ames, IA 50011}

\begin{abstract}
 The response field of a half-space anisotropic superconductor is
evaluated for an arbitrary weak external field source. Example sources  of a
 point magnetic moment and a circular current are considered in
detail. For the penetration depth $\lambda \ll L$ with $L$ being any other
relevant distance (the source size, or the distance between the source and the
superconductor), the major contribution to the response  is 
the $\lambda$ independent field of the source image. It is shown that the
absolute value of $\lambda$ cannot be extracted from the response field with a
better accuracy than that for the source position. Similar problems are
considered for thin films. 
%Low frequency polarized electromagnetic waves are
%reflected by anisotropic superconductors being elliptically polarized.  
  \end{abstract}
\pacs{74.20.-z,74.78.-w,74.78.Fk,74.81.-g}
\maketitle

\section{Introduction}

 A new experimental technique, the scanning SQUID microscopy (SSM), has
recently been developed for measuring magnetic field distributions due to 
 vortices exiting superconducting samples.\cite{Kirtley} Knowing the
distributions  one can, in principle, extract the London penetration depth
$\lambda$ (either isotropic or anisotropic) and its temperature
dependence.\cite{KirtleyT} In other implementation of this method, one can 
  measure the Meissner response of a superconductor to a
weak external source of static magnetic field. Again, the response  may provide
information about the $\lambda$ temperature dependence and its anisotropy. Also,
one can use tunable external field sources to study elementary forces acting upon
vortices exiting the surface.\cite{Moler-Gardner} Similar problems are
encountered in the magnetic force microscopy applied to surfaces of anisotropic
superconductors. 

There are a few publications addressing these problems for 
  isotropic superconductors.\cite{Clem,Milosevic,Obukhov} For
anisotropic materials, however, the response field is asymmetric even when the 
source has certain symmetries, and one cannot use methods developed for  
isotropic materials. To deal with the problem, one can utilize the
two-dimensional (2D) Fourier transform with respect to coordinates $x,y$ of the
interface, provided the  equations for the field distributions inside and outside
the superconducting half-space are {\it linear}. This is the case if one adopts
the London description for the field inside the material. Then, one solves the
remaining system of ordinary differential equations in the variable 
$z$, normal to the interface. This approach has been developed in Ref.
\onlinecite{KSL} for the problem of vortices crossing the superconductor surface. 

Let us consider a source with known field distribution   ${\bm
h}^s $ in the absence of superconductor. 
In the presence of the superconductor occupying the  half-space
$z<0$, the total field in vacuum $z>0$ can be written as 
\begin{equation}
{\bm h} = {\bm h}^s + {\bm h}^r\,,
\end{equation}
 where ${\bm h}^r$ is the response field which satisfies div$\,{\bm
h}^r=$ curl$\,{\bm h}^r=0$ in vacuum. One can look for this field as
$\nabla\varphi^r$ with the potential $\varphi^r$ obeying the Laplace equation
and the zero boundary condition far from the surface. The general form of such
a potential is
\begin{equation}
\varphi^r ({\bm r},z)=\int\frac{d^2{\bm k}}{(2\pi)^2}\,\varphi^r ({\bm k})\,
e^{i{\bm k}\cdot{\bm r}-kz}\,.
\label{Laplace}
\end{equation}
  Here, ${\bm r}=(x,y)$ and $z$ is directed normal to the superconducting flat
surface at $z=0$; $\varphi^r ({\bm k})e^{-kz}$ is the  2D Fourier
transform with respect to variables $x,y$ at any fixed  $z>0$. The potential
(\ref{Laplace}) is defined only in the upper half-space; hence, there is no
problem of uniqueness which is in general associated with the description of the
magnetic field by a potential. 

The field inside the superconductor satisfies London equations which read
in general anisotropic case as\cite{K81} 
\begin{equation}
h_i-\lambda^2m_{lk}\,e_{lst}\,e_{kni}\,  h_{t,ns}= 0 \,,\quad i=x,y,z\,.
\label{general}
\end{equation}
Here $e_{ikl}$ is the unit antisymmetric tensor, $m_{ij}$  is the mass
tensor, $h_{t,ns}$ abbreviates $\partial^2h_t/\partial x_n\partial x_s$. 
The average penetration depth $\lambda =(\lambda_a\lambda_b\lambda_c)^{1/3}$ is
related to the actual penetration depth for the currents, e.g., along the
crystal direction $a$: $\lambda_a=\lambda\sqrt{m_a}$. The masses are normalized
so that $m_am_bm_c=1$. 

The problem, therefore, is to match the solutions for the field inside and outside
the superconductor with boundary conditions of the field continuity at the
interface. 

\section{Response field}

Usually in situations of interest, the sample surface is normal to one of the
principal crystal directions. We call this direction $c$ and choose  the frame   
$x,y,z$ as coinciding with $a,b,c$.   In this situation, the mass
tensor is diagonal ($m_{xx}= m_a,\quad m_{yy}=m_b,\quad  m_{zz}= m_c$),
and Eqs. (\ref{general}) reduce to
\begin{eqnarray}
h_x
&+&\lambda_b^2(h_{z,x}-h_{x,z})_{,z}+\lambda_c^2(h_{y,x}-h_{x,y})_{,y}=0\,
,\nonumber\\
h_y &+&\lambda_c^2(h_{x,y}-h_{y,x})_{,x}+\lambda_a^2(h_{z,y}-h_{y,z})_{,z}=0\,,
\label{London_system}\\
h_z
&+&\lambda_a^2(h_{y,z}-h_{z,y})_{,y}+\lambda_b^2\,(h_{x,z}-h_{z,x})_{,x}=0\,.
\nonumber
\end{eqnarray}
One can replace any of one of these equations by div${\bm
h}=h_{x,x}+h_{y,y}+h_{z,z}=0$. It is convenient to replace the third one, and
exclude $h_{z,z}$ from the first two:
\begin{eqnarray}
\noindent  h_x-\lambda_b^2h_{x,xx}&-&\lambda_c^2h_{x,yy} -
\lambda_b^2h_{x,zz}+(\lambda_c^2-\lambda_b^2)h_{y,xy} =0,\nonumber\\ 
\noindent  h_y-\lambda_c^2h_{y,xx}&-&\lambda_a^2h_{y,yy} -
\lambda_a^2h_{y,zz}+(\lambda_c^2-\lambda_a^2)h_{x,xy} =0,\nonumber\\
   h_{z,z}= &-&\,(h_{x,x}+h_{y,y}) \,.
\label{London_1} 
\end{eqnarray}
Note that the first two equations are decoupled from the third. 
 
We now apply  the $x,y$ Fourier transform to Eqs.
(\ref{London_1}):
\begin{eqnarray}
&&(1+\lambda_b^2k_x^2 +\lambda_c^2k_y^2)h_x -\lambda_b^2 h''_x
+(\lambda_b^2-\lambda_c^2)k_xk_yh_y  =0\,,\nonumber\\
&&(\lambda_a^2-\lambda_c^2)k_xk_yh_x +(1+\lambda_c^2k_x^2 +\lambda_a^2k_y^2)h_y
-\lambda_a^2 h''_y   =0\,,\nonumber\\
&&i(k_xh_x+k_yh_y)+h_z'=0\,. \label{FT}
 \end{eqnarray}
Here,  $ h_i $ are functions of $ k_x,k_y$ and $z$,
and the prime denotes the derivative with respect to $z$. Hence, we are left with
the linear system of ordinary 2nd order differential equations with respect to
the variable $z$ for $h_i({\bm k},z)$.   

The solutions are linear combinations of simple exponentials:
\begin{equation}
h_i({\bm k},z) = \sum_n H_i^{(n)}({\bm k})e^{q_nz} \,.
\label{exp}
\end{equation}
The parameters $q_n$ and their number are to be determined. Substituting each
term of Eq. (\ref{exp}) in the system (\ref{FT}) we obtain a linear
homogeneous system for $H_i^{(n)}$: 
\begin{eqnarray}
\noindent
&&H_x\,(1+\lambda_b^2k_x^2+\lambda_c^2k_y^2-\lambda_b^2q^2)+H_y\,(\lambda_b^2-\lambda_c^2)k_xk_y
 =0,\nonumber\\
\noindent &&H_x\,(\lambda_a^2-\lambda_c^2)k_xk_y
+H_y\,(1+\lambda_c^2k_x^2+\lambda_a^2k_y^2-\lambda_a^2q^2)=0, \nonumber\\
\noindent &&i\,(k_xH_x +k_yH_y)+qH_z=0\,. \label{FT1} 
\end{eqnarray}
for each $n$; the superscript $n$ is omitted for brevity.
As has been mentioned, the first two equations here are decoupled from the third;
they have a nonzero solution provided  their determinant is zero. This gives a
quadratic equation for $q^2$ which is readily solved: 
\begin{eqnarray}
q^2_{1,2}&=&\frac{ P\pm\sqrt{Q}}{2\lambda^2_a\lambda^2_b}
\,, \label{q's}\\ 
P&=&\lambda^2_a+\lambda^2_b+\lambda^2_c(\lambda^2_bk_x^2+\lambda^2_ak_y^2)+
\lambda^2_a\lambda^2_bk^2\,,\nonumber \\
Q&=& P^2
-4\lambda^2_a\lambda^2_b(1+\lambda^2_bk_x^2+\lambda^2_ak_y^2)(1+\lambda^2_ck^2)\,.
\nonumber
\end{eqnarray}
  Note that $Q<P^2$ and both $q^2_1$ and $q^2_2$ are positive;
therefore,  there are only two  positive $q$'s (i.e., $n=1,2$) which satisfy the 
requirement of vanishing fields at $z=-\infty$. 

The quantities $q$ determine how the field attenuates in the superconductor. We
now have to find the ``amplitudes" $H_i$ from Eqs. (\ref{FT1})
for each $q$. It is worth noting that solving the homogeneous system of linear
equations implies, in fact,  expressing some unknowns in terms of others. To
this end, one has to determine the rank of the matrix of coefficients for the
system (\ref{FT1}), choose a proper subsystem to solve, etc.  The actual
procedure might differ depending on the situation in question. 

Solving the system (\ref{FT1}) we obtain:
 \begin{eqnarray}
H_x&=&-i\frac{d-\lambda_a^2q^2}{k_xq\,(\lambda^2_a-\lambda^2_b)}\, H_z
,\quad H_y=i \,\frac{ d-\lambda^2_bq^2}{k_yq\,(\lambda^2_a-\lambda^2_b)}\,H_z
\,,\nonumber\\
d&=&1+\lambda_b^2k_x^2+\lambda_a^2k_y^2\,.\label{Hxy}
\end{eqnarray}
 for each $q$ of Eq. (\ref{q's}). 

Let us turn now to the field ${\bm h}^s$ of the source. As a consequence of
div$\,{\bm h}^s=$ curl$\,{\bm h}^s=0$ out of the source, the 2D Fourier
components of ${\bm h}^s$ are not independent.  As with the response field, we
can look for this field in the from ${\bm h}^s=\nabla\,\varphi^s$ such that
$\nabla^2\varphi^s=0$. In our situation,  the source is situated in the upper
half-space $z>0$, and we are interested in the field ${\bm h}^s$ ``under" the
source. The general solution of the Laplace equation which vanishes as
$z\rightarrow -\infty$ is:
\begin{equation}
\varphi^s ({\bm r},z)=\int\frac{d^2{\bm k}}{(2\pi)^2}\,\varphi^s ({\bm k})\,
e^{i{\bm k}\cdot{\bm r}+kz}\,.
\label{Laplace_s}
\end{equation}
The 2D Fourier components of the source field at the interface then read:
\begin{equation}
 h_{\alpha}^s = i\,k_{\alpha} \varphi^s ({\bm k})\quad (\alpha=x,y)\,,\quad   
h_z^s =k\varphi^s ({\bm k})\,.\label{2dh}
\end{equation}
Hence, the boundary conditions take the form:
\begin{eqnarray}
  i\,k_{\alpha}(\varphi^s+\varphi^r) &=& H_{\alpha}^{(1)} +
H_{\alpha}^{(2)} 
\quad (\alpha=x,y),\nonumber\\   
 k(\varphi^s - \varphi^r) &=& H_z^{(1)} +H_z^{(2)}\,.\label{bc1} 
\end{eqnarray} 
Since the components $H_{\alpha}$ are expressed in terms of $H_z$'s, we can solve 
the system  (\ref{bc1}) to obtain:
  \begin{eqnarray}
\varphi^r &=& \varphi^s \,\frac{ k(q_1+q_2-k)- 
q_1q_2(1-1/d)}{ k(q_1+q_2+k)+ q_1q_2(1-1/d)}\,,\nonumber \\  
H_z^{(1)}&=&\frac{k[\varphi^s(q_2-k)-\varphi^r(q_2+k)]}{q_2-q_1}\,,\label{phi_an}
\\
H_z^{(2)}&=&\frac{k[\varphi^s(k-q_1)+\varphi^r(q_1+k)]}{q_2-q_1}\,.\nonumber
\end{eqnarray}
Thus, the response outside and inside the superconductor is
expressed in terms of the source field $\varphi^s$. It is worth noting that since $\varphi^s({\bm k})$ can be replaced with
$ h_z^s({\bm k})/k$, the response field can be expressed in terms of the
$z$ component of the source field at the interface.  

 
 One can see, in particular, that the total flux in
$z$ direction "reflected" by the superconductor is equal and opposite in
sign to the incident flux of the source crossing the interface. Indeed,
this flux is 
\begin{equation}
h_z^r|_{k=0}=-k\varphi^r|_{k=0} = -k \varphi^s|_{k=0}=-h_z^s|_{k=0} \,.
\end{equation}


It is seen from the general formulas (\ref{Hxy}) that the case
$\lambda_a=\lambda_b$ is singular. The formal reason for this is that the 
first two equations of the system (\ref{FT1}) (which we have used to express
$H_x$ and $H_y$ in terms of $H_z$) are no longer  independent.  This
situation should be treated separately. 


\subsection{Small $\bm\lambda$ }

In many situations, the penetration depths are small relative to other relevant
lengths in the problem such as the distance $z_0$ between the source and the
interface. The characteristic $k$'s then satisfy $k\lambda \ll 1$. The relations
between the response and the source fields then simplify. In this
approximation, Eqs. (\ref{q's}) give $q_1=1/\lambda_a$ and $q_2=1/\lambda_b$ for
$\lambda_a<\lambda_b$. Then, Eqs. (\ref{Hxy}) yield $d=1$ and
 \begin{eqnarray}
H_x^{(1)}=0\,,\quad H_z^{(1)}=-ik_y\lambda_aH_y^{(1)}\,,\label{Hy1_sm}\\
H_y^{(2)}=0\,,\quad H_z^{(2)}=-ik_x\lambda_bH_x^{(2)}\,.\label{Hx2_sm}
\end{eqnarray}
Finally, the boundary conditions give:
  \begin{eqnarray}
\varphi^r &=& \varphi^s \,\Big(1-2\,\frac{\lambda_bk_x^2 +\lambda_ak_y^2}{k
}\Big),\label{f_sm}\\  
H_y^{(1)}&=&{k_y\over k_x}H_x^{(2)}=2ik_y\varphi^s \,\Big(1-\frac{\lambda_bk_x^2
+\lambda_ak_y^2}{k }\Big).\label{hx1_sm}    
\label{small_l}
\end{eqnarray}
Note that in this approximation $\lambda_c$ does not enter the result; in other
words, the currents along $z$ are small and can be disregarded. The results 
(\ref{Hy1_sm})-(\ref{hx1_sm}) can also be obtained   directly 
  starting with the London Eqs. (\ref{London_system}) and taking advantage of 
  $\partial/\partial z\gg \partial/\partial x_{\alpha}$. 

It is also worth observing that in zero order in   $k\lambda\ll 1$,
$\varphi^r=\varphi^s$. This means that the $x$ and $y$ components of the response
field are the same as those for the source, whereas $h_z^r=-h_z^s$, see Eqs.
(\ref{Laplace}) and (\ref{Laplace_s}). In other words, in this approximation the
response field is the mirror image of the source field. 

\subsection{Isotropic materials}

It is readily seen that in this case 
\begin{equation}
q_1=q_2=\sqrt{\lambda^{-2}+k^2}\,.\label{q_is}
\end{equation}
The first two of Eqs. (\ref{FT1}) turn identities, whereas the third gives one of
$H_i$'s in terms of two others, e.g., $H_z=(k_xH_x +k_yH_y)/iq$. The boundary
conditions then give:
\begin{equation}
\varphi^r =\frac{q-k }{q+k}\,\varphi^s \,,\quad H_{\alpha}=\frac{2iq
}{q+k}\,k_{\alpha}\,\varphi^s\,. 
\label{isotropic}
\end{equation}
%These formulas, in principle, give both the response field outside and the field
%inside the superconductor for any given source.

\subsection{${\bm \lambda}_a={\bm\lambda}_b< {\bm\lambda}_c$}

This is the case, e.g., of BSCCO crystal with the $ab$ plane being the surface. 
Equations (\ref{q's}) yield:
\begin{equation}
q_1=\sqrt{\lambda_{ab}^{-2}+k^2}\,,\qquad
q_2=\sqrt{\lambda_{ab}^{-2}+\gamma^2k^2}\,,\label{BSCCO1}
\end{equation}
where $\gamma=\lambda_c/\lambda_{ab}$ is the anisotropy parameter. Substituting
$q_1$ in the first two equations of the system (\ref{FT1}), we obtain two
identical results $H_x^{(1)}k_y-H_y^{(1)}k_x=0$ which together with the third
equation yield:
\begin{equation}
H_x^{(1)}=\frac{iq_1k_x}{k^2}\,H_z^{(1)}\,,\quad
H_y^{(1)}=\frac{iq_1k_y}{k^2}\,H_z^{(1)}\,.\label{BSCCO2}
\end{equation}
Doing the same for $q_2$, we obtain from the first two equations
$H_x^{(2)}k_x+H_y^{(2)}k_y=0$ which is compatible with the third only if
$H_z^{(2)}=0$. Hence, we have:
\begin{equation}
 H_z^{(2)}=0\,,\qquad
H_y^{(2)}=-\,\frac{ k_x}{k_y}\,H_x^{(2)}\,.\label{BSCCO3}
\end{equation}
The boundary conditions (\ref{bc1}) 
yield:
\begin{equation}
\varphi^r =\frac{q_1-k}{q_1+k}\,\varphi^s \,,\quad H_z^{(1)}=\frac{2k^2
}{q_1+k}\, \varphi^s\,,\quad {\bm H}^{(2)}=0\,.\label{BSCCO5}
\end{equation}
Note that for this  case, $q_2$ along with $\lambda_c$
drop off the result for {\it any} source. In other words, the response is as if
the superconductor were isotropic with the penetration depth $\lambda_{ab}$. 

\subsection{${\bm \lambda}_c={\bm\lambda}_a< {\bm\lambda}_b$}

This is the case of  the screening by the ``side surface" of a
uniaxial crystal like BSCCO. Our notation is, however, differs from the commonly
used (for standard notation, we should have replaced in our formulas   
$\lambda_c$ and $\lambda_a$ with 
$\lambda_{ab}$ and $\lambda_b\rightarrow\lambda_c$). We then obtain using Eq.
(\ref{q's}):
\begin{equation}
q_1=\sqrt{\lambda_{a}^{-2}+ k^2}\,,\quad 
q_2=\sqrt{1+\lambda_{b}^2k_x^2+\lambda_a^2k_y^2}\Big/\lambda_b\,.
\label{q_1_2}
\end{equation}
With $q=q_1$, Eq. (\ref{Hxy}) gives
\begin{equation}
H_x^{(1)}=i\,\frac {k_x}{q_1}\,H_z^{(1)}\,,\quad
H_y^{(1)}=i\,\frac{1+\lambda_{a}^2k^2_y}{\lambda_a^2k_yq_1}\,H_z^{(1)}\,,\label{H1}
\end{equation}
and for $q=q_2$
\begin{equation}
 H_x^{(2)}=i\,\frac{q_2}{k_x}\,H_z^{(2)}\,,\qquad
H_y^{(2)}=0\,.\label{H1}
\end{equation}
Finally, Eqs. (\ref{phi_an}) give $\varphi^r$,  $H_z^{(1)}$, and $H_z^{(2)}$ in
terms of the source field $\varphi^s$.  

\section{Examples of sources}

To apply formulas derived above for the response fields outside and inside
superconductor, one needs the Fourier transform $\varphi^s({\bm k})$ of the source
field at the interface. Below we provide  examples for which $\varphi^s({\bm
k})$ can be calculated analytically.

\subsection{Point magnetic moment}

Consider a magnetic moment ${\bm\mu}$ situated at the hight $z_0$ above the
superconductor. The corresponding potential (in the absence of the
superconductor) at $z=0$ is
\begin{equation}
\varphi^s=-\,\frac{{\bm\mu}\cdot   {\bm R}}{R^3} =  \frac{\mu_z\,z_0 -
{{\bm\mu}}\cdot{\bm r}}
 {(r^2+ z_0^2)^{3/2}}\,.
\label{phi_moment}
\end{equation}
Here, ${\bm R}=(x,y,z-z_0)$ is the radius-vector
originating at the source; the first minus sign  is due to the definition
${\bm h}^s=\nabla\,\varphi^s$.  The 2D Fourier transform is:\cite{Grad}
\begin{equation}
\varphi^s({\bm k})  = 2\pi  \,e^{-k z_0}\Big(\mu_z+i\,\frac{
{{\bm\mu}}\cdot{\bm k}}{k}\Big)\,.
\label{FT-mu}
\end{equation}
The response field can now be calculated with the help of Eqs. (\ref{phi_an}).
In general, this can be done numerically; for $\lambda \ll z_0$,
the analytic evaluation is possible.  

 As an example take the moment  directed along $z$ above the flat isotropic
superconducting surface; for the isotropic thin film this problem has been
considered in Ref. \onlinecite{Milosevic}. According to Eq. (\ref{f_sm}) for the
isotropic case
$\varphi^r =
\varphi^s  (1-2 k\lambda)$. Transforming back to real space we obtain:
 \begin{eqnarray}
\varphi^r (r,z) &=& \mu\Big( \frac{Z}{R^3} -2\lambda\,
\frac{2Z^2-r^2}{ R^5}\Big)\,,\nonumber\\
 R &=& \sqrt{r^2+Z^2}\,,\qquad Z=z+z_0\,.
\label{phi(r,z)is2}
\end{eqnarray}
Here, the first term is the field of a moment $- {\bm \mu}$ at
$z=-z_0$, i.e., of the image source. The second term is  the field of the 
proportional to $\lambda$ magnetic quadrupole at the same point. 

 For the magnetic force microscopy, the quantity of interest is the interaction
energy which is given by 
\begin{eqnarray}
{\cal E} &=& -\,{\bm\mu}\cdot {\bm h}^r(0,0,z_0)\label{E}\\ 
&=& -\int \frac{d^2{\bm k}}{(2\pi)^2}\,\varphi^r({\bm
k})e^{-kz_0}(-\mu_zk+i \mu_{\alpha}  k_{\alpha})\nonumber
\end{eqnarray}
Substituting here Eqs. (\ref{f_sm}) and (\ref{FT-mu}) and 
integrating we obtain:
\begin{eqnarray}
{\cal E} &=&\frac{\mu_z^2}{4z_0^3}\Big[1+\frac{3(\lambda_a+\lambda_b)}{4z_0}\Big]
+\frac{\mu_x^2+\mu_y^2}{8z_0^3}\nonumber\\
&+&\frac{3}{32z_0^4}\big[\mu_x^2(3\lambda_b+\lambda_a)+\mu_y^2(3\lambda_a+
\lambda_b)\big]\,.\label{E_tot}
\end{eqnarray} 
In addition to a repulsive force $-\partial{\cal E}/\partial z_0$, the
magnetic moment ${\bm \mu}$ experiences a torque, because the energy depends on
the moment orientation. It is readily seen that if the position and the value of
the magnetic moment are fixed, the minimum of ${\cal E}$ corresponds to  
${\bm \mu}$ situated in the plane $xy$  and parallel to the direction of largest
$\lambda$. The $z$ component of the torque is easily evaluated:
\begin{equation}
\tau_z=\frac{3\mu_{\perp}^2}{16 z_0^4}\,(\lambda_a-\lambda_b)\,\sin 2\beta\,,
\label{torque}
\end{equation}
 where $\mu_{\perp}$ is the in-plane part of the magnetic moment, and $\beta$ is
the angle between ${\bm \mu}_{\perp}$ and $\hat x$. For $\lambda_a=\lambda_b$,
the position of ${\bm \mu}_{\perp}$ in the $xy$ plane is arbitrary, still there
is a torque which tends to rotate ${\bm \mu}$ out of the $z$ direction and to
place it in the $xy$ plane.
 
\subsection{The source as a current loop}

Let the source be a circular current of a radius  $a$ situated in the plane
$z=z_0$. If $a$ and $z_0$ are of the same order of magnitude
(practically, they are both of a few-micron size),  modelling of the
source by a point-size magnetic moment does not suffice. The scalar
potential of the field created by a loop in the plane $z=0$ reads:\cite{LL}   
\begin{equation} 
\varphi^s ({\bm r}) = -\,\frac{I}{c}\int\frac{d{\bm S}\cdot {\bm R}}{R^3}\,.
\label{LL}
\end{equation}
The source current $I$ flows in the counterclockwise direction relative to the
$z$ axis so that an area element $d{\bm S} = dS\,{\hat {\bm z}} $. 
${\bm R}$ is the radius-vector from this element to the interface point $({\bm
r},0)$, and the integral is over the area of the current contour. The position of
the element
$d{\bm S}=d^2{\bm r}^{\prime}{\hat {\bm z}} $ in our situation is $({\bm
r}^{\prime},z_0)$ so that
$ {\bm R} = {\bm r} - {\bm r}^{\prime} - z_0{\hat z}$, and 
\begin{equation} 
\varphi^s ({\bm r},0) = \frac{Iz_0}{c}\int\frac{d^2{\bm r}^{\prime}  
 }{[({\bm r} - {\bm r}^{\prime})^2+z_0^2]^{3/2}}\,.
\label{LL1}
\end{equation}
For the circular loop of the radius $a$, the integral here is over the
circle area. Comparing this with Eq. (\ref{phi_moment}) we see that the
source field can be considered as created by magnetic moments
distributed uniformly over the loop area with the density $I{\hat {\bm z}}/c$ so
that the total moment of the loop is $\pi a^2I{\hat {\bm z}}/c$. 

The 2D Fourier transform of $z_0[({\bm r} - {\bm
r}^{\prime})^2+z_0^2]^{-3/2}$ with respect to ${\bm r}$ is $2\pi
e^{-kz_0}e^{-i{\bm k}\cdot{\bm r}^{\prime}}$, see Eq. (\ref{FT-mu}); therefore,
\begin{eqnarray} 
\varphi^s ({\bm k}) &=& \frac{I}{c}\,2\pi e^{-kz_0}\int d^2{\bm
r}^{\prime} e^{-i{\bm k}\cdot{\bm r}^{\prime}} \nonumber\\ 
&=&  \frac{4\pi^2Ia}{ck} e^{-kz_0}J_1(ka) \,.
\label{ft_loop}
\end{eqnarray}
According to (\ref{BSCCO5}) the 2D Fourier transform of the response potential
for isotropic superconductor is:
\begin{equation} 
\varphi^r ({\bm k}) =   
 \frac{4\pi^2Ia}{ck}\,\frac{q-k}{q+k}\, e^{-kz_0}J_1(ka) \,.
\label{ft_loop_r}
\end{equation}
The $z$-component of the response field follows:
\begin{eqnarray} 
h_z^r ({\bm r},z) =-\, 
 \frac{Ia}{c}\,\int  d^2{\bm k}\frac{q-k}{q+k}\, e^{-k(z+z_0)}J_1(ka)
e^{i{\bm k}\cdot{\bm r}}\nonumber\\
=-\, \frac{2\pi Ia}{c}\,\int_0^{\infty}  dk\,k\,\frac{q-k}{q+k}\,
e^{-k(z+z_0)}J_1(ka)
 J_0(kr)\,.
\label{ft_loop_hz}
\end{eqnarray}

 Further, one can evaluate the response flux through a flat probe  placed above
the superconductor. If the probe is a circular loop of a radius
$a_p$ with the center at ${\bm r}=0$ at the height $z_p\le z_0$, the flux is
given by 
\begin{equation}
\Phi_z^r = -\frac{4\pi^2 Iaa_p}{c}\,\int_0^{\infty}  dk\,
\frac{q-k}{q+k}\, e^{-k(z_0+z_p)}J_1(ka)J_1(ka_p)\,.\label{40}
\end{equation}
All formulas for the isotropic case have been worked out earlier by Clem and
Coffey making use of the cylindrical symmetry of the problem.\cite{Clem}

 
\subsection{Interaction with vortices}

It is of interest to evaluate the force acting on the vortex tip by the
screening currents in the sample created by a circular current as the field
source; experiments for which this  is relevant are described in Ref.
\onlinecite{Moler-Gardner}. A similar problem for isotropic films and the
magnetic moment as a source has been considered in Ref.
\onlinecite{Milosevic}. We do this for materials isotropic in the $xy$ plane, for
which the force depends only on the distance $r$ from the loop center along with
the loop hight $z_0$. One can calculate
$F_x(x,0;z_0)$ and in the result replace $x$ with $r$:
\begin{eqnarray} 
F_x(x,0) &=&  
 \frac{\phi_0}{c}\,\int_{-\infty}^0  dz j_y = 
\frac{\phi_0}{4\pi}\,\int_{-\infty}^0  dz \,(h_{x,z}-h_{z,x})\nonumber\\
&=&\frac{\phi_0}{4\pi}\,\int \frac{d^2{\bm
k}}{(2\pi)^2}\,e^{ik_xx}\Big(H_x^{(1)}-\frac{ik_x}{q_1}H_z^{(1)}\Big)\,.
\label{Fx1}
\end{eqnarray}
%This function is plotted in Ref. \onlinecite{Moler_Gardner}
With the help of Eqs. (\ref{BSCCO1}), (\ref{BSCCO2}), and (\ref{ft_loop}) we
  obtain
\begin{equation}
F_r(r) =-\frac{\phi_0 Ia}{c \lambda_{ab}^2}\,\int_0^{\infty} 
\frac{ dk\,k\,e^{-k z_0}}{  q_1(q_1+k) } \,J_1(ka) J_1(kr)\,.
\end{equation}
 
Using $1/\lambda^2=q_1^2-k^2$ and the substitution $t=ka$, we write:
\begin{eqnarray}
F_r  &=&-\frac{\phi_0 I}{2ca}\,\int_0^{\infty} dt\,G\Big({\lambda\over
a}\,t\Big) t\,e^{-z't} J_1(t) J_1( r't) \,,\nonumber\\
G&=& 1-{\lambda\over a}\,t\Big(1+{\lambda^2\over a^2}\,t^2\Big)^{-1/2} ,
\end{eqnarray}
where $z'=z_0/a$ and $r'=r/a$. Usually, the parameter $\lambda/a\ll 1$.
Besides, due to the factor $e^{-z't} J_1(t)$,  the region contributing to the
integral is $0<t<\min (1/z',1)$ because of the oscillating $J_1(t)$ at
large
$t$ (unless $r'\approx 1$). Then, one can expand $G$ in powers of $\lambda t/a$
and keep only the linear term:
\begin{eqnarray}
F_r  &=&-\frac{\phi_0 I}{2ca}\,\Big(f_0-{\lambda\over
a}\,f_1\Big)\,,\nonumber\\ 
f_0&=&\int_0^{\infty} dt\,t\,e^{-z't}J_1(t) J_1( r't) \,,\label{f12}\\ 
f_1&=&\int_0^{\infty} dt\,t^2\,e^{-z't} J_1(t)J_1( r't)\,,\nonumber  
\end{eqnarray}
The integrals here can be expressed in terms of the hypergeometric functions
convenient for the numerical evaluation, see Appendix A.

One can also define a potential energy $U(r)$ such that $F_r=-dU/dr$:
\begin{eqnarray}
U  &=&-\frac{\phi_0 I}{2c}\,\Big(u_0-{\lambda\over
a}\,u_1\Big)\,,\nonumber\\ 
u_0&=&\int_0^{\infty} dt\, e^{-z't}J_1(t) J_0( r't) \,,\label{u12}\\ 
u_1&=&\int_0^{\infty} dt\,t \,e^{-z't} J_1(t)J_0( r't)\,,\nonumber  
\end{eqnarray}


This energy is defined so that $U(\infty)=0$; its value at the origin is
\begin{equation}
U(0) =-\frac{\phi_0
I}{2c}\,\Big[1-\frac{z_0}{(z_0^2+a^2)^{1/2}} -\frac{\lambda
a^2}{(z_0^2+a^2)^{3/2}}\,\Big]\,.\label{U}
\end{equation}
 Thus, the source loop creates a potential well of the depth $|U(0)|$ or a barrier
of the  hight $|U(0)|$ for vortices underneath depending on the current and
vortex directions. This opens an interesting possibility for studying the
behavior of vortices (or antivortices) in tunable potentials.\cite{Moler-Gardner}
In a similar manner, one can evaluate interactions of vortices with other types
of sources.
 
\subsection{Accuracy of the SSM determination of $\lambda$}

Magnetic fluxes of the response field, in particular, $\Phi_z=\int d^2{\bm
r}\,h_z^r({\bm r},z_p)$   (the integral is over the area of a pick-up coil
placed at the hight $z_p$ above the superconducting surface) can be measured with
high accuracy for a given source, given geometry of the coil, and a known height
$z_p$. In principle, this leads to a possibility to measure  the
penetration depth. However, the accuracy of this determination in bound by the
accuracy with which the heights $z_0$ and $z_p$ are known. To demonstrate this
consider the response field of an isotropic material for small $\lambda$'s:
 \begin{eqnarray}
&&h_z^r({\bm r},z_p)  = -\int \frac{d^2{\bm k}}{(2\pi)^2}\,k\varphi^r({\bm
k})e^{i{\bm k}\cdot{\bm r}-kz_p} \,\nonumber\\
&&=-\int \frac{d^2{\bm k}}{(2\pi)^2}\,k\varphi^s({\bm
k} )(1+2\lambda k)e^{i{\bm k}\cdot{\bm r}-kz_p}.    
\end{eqnarray}
If $z_p$ varies by $\delta z_p$, the response field variation is
\begin{equation}
\delta h_z^r =\delta z_p \int \frac{d^2{\bm k}}{(2\pi)^2}\,k^2\varphi^s({\bm
k})(1+2\lambda k)e^{i{\bm k}\cdot{\bm r}-kz_p}\,.
\end{equation}
If only $\lambda$ varies, we have:
\begin{equation}
\delta h_z^r =-2\,\delta\lambda \int \frac{d^2{\bm
k}}{(2\pi)^2}\,k^2\varphi^s({\bm k})(1+2\lambda k)e^{i{\bm k}\cdot{\bm r}-kz_p}\,,
\end{equation}
in other words,
\begin{equation}
\frac{\delta h_z^r}{\delta\lambda} =-2\,\frac{\delta
h_z^r}{\delta z_p} \,.
\end{equation}

Therefore, the accuracy of extracting $\lambda$ from the data on $h_z^r$ cannot
be much better than the knowledge of $z_p$ (the same is true about the source
position $z_0$). The latter is usually known within a fraction of a micron. This
is a severe restriction upon the accuracy of the absolute determination of
$\lambda$.  Still, in principle, SSM allows one to determine accurately the
temperature dependence of $\lambda$ (for  fixed $z_p$ and $z_0$).  

\section{Thin films}

The Meissner response of superconducting thin films can be probed by the SSM
method in yet greater detail than that of bulk samples. The films often
have a large Pearl length $\Lambda=2\lambda^2/d$  ($d$ is the film
thickness) exceeding substantially the size of the SSM sensing loop, making the
SSM measurement to a local probe. Formally, the problem of a thin film in a
field of an external source is simpler than that of the superconducting
half-space, because in the film case there is no ``internal problem" to solve;
instead, the film provides a boundary condition for the outside field
distribution. 

Consider a film in the $x,y$ plane made of a
uniaxial material with the $c$ axis at an angle $\theta$ to the film normal $z$. 
We write the London equation (\ref{general}) for
$i=z$, $h_z-4\pi\lambda^2(m_{xl}j_{l,y}-m_{yl}j_{l,x})/c
=0$, and integrate it over the film thickness:
\begin{equation}
h_z-\frac{2\pi}{c}\Lambda (m_{xx}g_{x,y}-m_{a}g_{y,x}) =0\,,\quad
\Lambda={2\lambda^2\over d},
\label{London_film}
\end{equation}
where $\bm g$ is the sheet current. Note that only two components of the mass
tensor, $m_{xx}=m_a\cos^2\theta+m_c\sin^2\theta$ and
$m_{yy}=m_a$, determine the film anisotropy. 

Using the relation of sheet currents ${\bm g}$ to the tangential components
of the response field,  
\begin{equation}
{2\pi\over c} g_x =  -h_y^r(+0)\,,\,\,\,{2\pi\over c}g_y =h_x^r(+0)  \,, 
\label{e4}
\end{equation}
($+0$ denotes  the upper  film face; the
tangential components satisfy  $h_t^r(+0) =-h_t^r(-0) $), we obtain for $z=+0$: 
\begin{equation} 
  h_z^s+h_z^r +\Lambda (m_{xx}\, h_{y,y}^r  +m_a\, h_{x,x}^r 
 )  =0\,,
\label{L-film}
\end{equation}
   for more detail, see Ref.\onlinecite{KSL}. 
% Further, since $h_z=(\varphi^s+\varphi^r),z$ whereas the tangential components
% of the source field are continuous, e.g., $h_x^+-h_x^- = 2\varphi^r_{,x}$, the
% 2D Fourier transform of the last equation reads:
Further, since $h_z^s+h_z^r
=k(\varphi^s-\varphi^r)$ 
%$k(\varphi^s-\varphi^r)-2\Lambda (m_{xx}\,k_y^2+m_ak_x^2)=0$. This yields
we obtain:
\begin{equation}
\varphi^r =\frac{k\varphi^s}{k+ \Lambda
(m_{xx}\,k_y^2 +m_ak_x^2)}  \,,\quad z=+0.
\label{potential-film}
\end{equation}
In particular, for $\theta=0$  we have  
\begin{equation}
\varphi^r
=\frac{\varphi^s}{1+ \Lambda_a
 k}  \,,\qquad \Lambda_a=m_a\Lambda ={2\lambda_{ab}^2\over d}\,.
\label{c||z}
\end{equation}

This result holds also for isotropic materials where $m_a=1$. In this case,
one can readily obtain the field $h_z$ for the circular current  
(\ref{ft_loop}) at the height $z_0$ above the film; its 2D Fourier transform
for $0<z< z_0$ is given by
\begin{eqnarray}
  h_z^r(k,z) &=& -k\varphi^re^{-kz} \nonumber\\&=&- 
\frac{4\pi^2Ia}{c}\, \frac{J_1(ka) }{1+
k\Lambda}\,e^{-k(z+z_0)}\, .
\label{hz}
\end{eqnarray}
This field can be measured by SSM.\cite{John} 

For completeness, we write down also the field on the opposite film side, i.e for
$z<0$, where the response potential is given by 
\begin{equation}
\varphi^r ({\bm r},z<0)=\int\frac{d^2{\bm k}}{(2\pi)^2}\,\varphi^r ({\bm k})\,
e^{i{\bm k}\cdot{\bm r}+kz}\,.
\label{Laplace2}
\end{equation}
The London boundary condition (\ref{L-film}) should now be written in terms of
the field components at $z=-0$:
\begin{equation} 
  h_z^s+h_z^r -\Lambda (m_{xx}\, h_{y,y}^r  +m_a\, h_{x,x}^r 
 )  =0\,,
\label{L-film-}
\end{equation}
which yields the 2D Fourier transform of the response potential
(\ref{potential-film}) with the minus sign. Proceeding as above, we obtain for
the isotropic case:
\begin{equation}
h_z^r(k,z<0) = -
\frac{4\pi^2Ia}{c}\,\frac{J_1(ka) }{1+
k\Lambda}\,e^{k(z-z_0)}\, .
\label{hz}
\end{equation}

The force acting upon a Pearl vortex situated in the film at a radial distance
$r$ from the current ring center is readily evaluated:
\begin{equation} 
F_r(r) =  
 -\frac{\phi_0Ia}{c}\,\int_0^{\infty} \frac{ 
k\,e^{-kz_0}}{1+k\Lambda}J_1(ka)J_1(kr) \,dk\,.
\label{Fr}
\end{equation}


\acknowledgements
These notes were prepared for modeling SSM experiments by John Kirtley and
Kathy Moler; I thank them for encouraging me to put the notes
in a publishable form. Discussions with Walter Hardy, Roman Mints, and Brian
Gardner helped me to focus on calculations which might be useful. In part, the
work was supported by Binational US-Israel Science Foundation. Ames Laboratory is
operated for US DOE by the Iowa State University under Contract No. W-7405-Eng-82.

\appendix

\section{}

The integral $u_0$ of Eq. (\ref{u12}) is known  (see
Ref.\onlinecite{Grad}, 6.612.3):
\begin{eqnarray}
I(r,z) &=& \int_0^{\infty} e^{-zt} J_1(t) J_1( rt)\,dt
%=\frac{1}{\pi\sqrt{r}}\,Q_{1/2}(u)
\nonumber\\
&=&\frac{1}{2^{5/2}\sqrt{r}u^{3/2}}\,_2F_1\Big({5\over 4},{3\over 4};2;{1\over
u^2}\Big)\,,
\label{I}
\end{eqnarray}
where the primes by $z$ and $r$ are omitted, and  $u=(z^2+r^2+1)/2r $.
Therefore: 
\begin{eqnarray}
&&  \int_0^{\infty}t e^{-zt} J_1(t) J_1( rt)\,dt =-\frac{\partial I}{\partial
z}\label{II}\\ 
&&=\frac{3z r^{-3/2}}{2^{7/2}u^{5/2}}\Big[\, _2F_1\Big({5\over
4},{3\over 4};2;{1\over u^2}\Big)+{5\over 8u^2}\,_2F_1\Big({9\over 4},{7\over
4};3;{1\over u^2}\Big)\Big].\nonumber 
\end{eqnarray}


\begin{references}

\bibitem{Kirtley} J.R. Kirtley,  C.C. Tsuei, M. Rupp,  J.Z. Sun, Lock See
Yu-Jahnes, A. Gupta, M.B. Ketchen, K.A. Moler, and M. Bhushan,
% {\it et al.} 
Phys. Rev. Lett. {\bf 76}, 1336 (1996).
 
\bibitem{KirtleyT} J.R. Kirtley, C.C. Tsuei,  K.A. Moler, V.G. Kogan, J.R.
Clem, A.J. Turberfield,
% {\it et al.} 
Appl. Phys. Lett. {\bf 74}, 4011 (1999).

\bibitem{Moler-Gardner} B.W. Gardner, J.C. Winn, D.A. Bonn, R. Liang, W.N. Hardy, J.R.
Kirtley, V.G. Kogan, and K.A. Moler,   
Appl. Phys. Lett. {\bf 80}, 1010 (2002).


\bibitem{Clem}J.R. Clem, M.W. Coffey, \prb {\bf 46}, 14662 (1992); M.W. Coffey, 
\prb {\bf 57}, 11648 (1998); M.W. Coffey, \prl {\bf 83}, 1648 (1999).

\bibitem{Milosevic}M.V. Milosevic, S.V. Yampolskii, and F.M. Peters, \prb, {\bf
66}, 174519 (2002).

\bibitem{Obukhov}Yu.V. Obukhov, J. Low Temp. Phys. {\bf 114}, 277 (1999).

\bibitem{K81}  V.G. Kogan, \prb, {\bf 24}, 1572 (1981).

\bibitem{KSL}  V.G. Kogan, A.Yu. Simonov, and M. Ledvij,
Phys. Rev. B {\bf 48}, 392 (1993).

\bibitem{Grad}I.S. Gradshteyn and I.M. Ryzhik, {\it Tables of
Integrals, Series, and Products}, Academic Press, 1980.

%\bibitem{LL}L.D. Landau and E.M. Lifshitz, {\it Electrodynamics of
%Continuous Media}, Pergamon, Oxford, New York, 1984.

%\bibitem{Field_Theory}L.D. Landau and E.M. Lifshitz, {\it The Classical Theory
%of Fields}, Pergamon, Oxford, New York, 1975.

\bibitem{John} J.R. Kirtley, unpublished.

\end{references}

\end{document}



 
\section{Meissner response to low frequency electromagnetic waves}

 Let us consider a plane electromagnetic wave falling on the 
superconducting surface $xy$ with two different penetration depths
$\lambda_a$ and $\lambda_b$. The incident wave vector  ${\bm
k}_0=(\omega/c)(\sin\theta,0,-\cos\theta)$ forms the angle $\theta$ with the $z$
axis and is situated in the plane $xz$; for the outgoing wave ${\bm
k}_1=(\omega/c)(\sin\theta,0,\cos\theta)$.  The problem of
reflection from  superconductors is relatively simple for low frequencies
$\omega$ for which one can use the concept of the surface impedance.\cite{LL} For
this, the London penetration depth should be small relative to the wave
length, $\lambda(T)\ll c/\omega$. Besides, the condition $\hbar\omega\ll
\Delta(T)$ ($\Delta(T)$ is the energy gap) should hold; otherwise the
superconductivity becomes irrelevant and the material behaves as normal metal. At
low temperatures, the latter restriction is dominant, which places the possible
frequencies in the deep infrared or lower. 

If the above conditions are met, instead of solving both the internal and 
external problems, one can use the 
concept of the surface impedance which suffices to find the outside
fields.\cite{LL} In our case, the surface impedance is a two-rank tensor
$\zeta_{\alpha\beta}$ with the eigenvalues
 \begin{equation}
 \zeta_{a} = i\,{\omega\over c}\,\lambda_a\,,\qquad  \zeta_{b} = i\,{\omega\over
c}\,\lambda_b\,. \label{zeta}
\end{equation}
Here, $a,b$ are the principal crystal directions of the interface. In the frame
$xy$, the components of this tensor are: 
\begin{eqnarray}
\zeta_{xx}&=& \zeta_a \cos^2\varphi+\zeta_b \sin^2\varphi\,,\nonumber\\
\zeta_{yy}&=& \zeta_a \sin^2\varphi+\zeta_b \cos^2\varphi\,, \label{zeta_ij} \\  
 \zeta_{xy}&=&  \zeta_{yx}=(\zeta_a -\zeta_b)\cos \varphi \sin \varphi\,,
\nonumber
\end{eqnarray}
where $\varphi$ is the angle between the axes $x$ and $a$. 

The surface impedance boundary condition
 $ E_{\alpha}=\zeta_{\alpha\beta}({\bm
H}\times {\hat {\bm z}})_{\beta} $ 
  relates the tangential components of the fields:
\begin{eqnarray}
E_x=\zeta_{xx}H_y- \zeta_{xy} H_x\,, \nonumber\\
E_y=\zeta_{yx}H_y- \zeta_{yy} H_x\,. \label{imp_bc}  
\end{eqnarray}
Since all components of the tensor $\zeta$ are imaginary, the averaged over time
energy flux $S_z=c\,{\rm Re} ({\bm E}_t\times {\bm H}^*_t)/8\pi$ through the
interface  is zero ($t$ stands for tangential components). In other words, the
reflection coefficient is unity. Clearly, this holds as long as the normal
currents of quasiparticles along with other possible dissipative mechanisms (as
e.g., the plasma resonance in layered compounds) can be neglected. 

Using the plane-wave relation ${\bm H}={\hat {\bm k}}\times {\bm E}$  for both
incident and outgoing waves one can express\cite{LL} the complex amplitudes of the
outgoing wave ${\bm E}_1,{\bm H}_1$ in terms of the incident fields ${\bm E}_0,
{\bm H}_0$. Consider first the $y$-polarized electric field: 
${\bm E}_0=E_0(0,1,0)$ and ${\bm H}_0=E_0(\cos\theta, 0,\sin\theta)$. 
 After straightforward algebra we obtain:
\begin{eqnarray}
H_{1x}&=& \cos\theta\,\frac{(\cos\theta-\zeta_{xx})(1+
\zeta_{yy}\cos\theta)+\zeta_{xy}^2\cos\theta}{(\cos\theta-\zeta_{xx})(1-
\zeta_{yy}\cos\theta)-\zeta_{xy}^2\cos\theta}\,E_0\,,\nonumber\\
H_{1y}&=& -\frac{2\cos\theta\,\zeta_{xy} }{(\cos\theta-\zeta_{xx})(1-
\zeta_{yy}\cos\theta)-\zeta_{xy}^2\cos\theta}\,E_0\,,\nonumber\\
H_{1z}&=&-\tan\theta\,H_{1x}\,. \label{H1} 
\end{eqnarray} 
The plane-wave formula ${\bm E}= {\bm H}\times{\hat {\bm k}}$ then gives:
\begin{eqnarray}
E_{1x}&=&- E_{1z}\cot\theta=H_{1y} \cos\theta\,,\nonumber\\  
E_{1y}&=&-H_{1x}\,/\cos\theta\,. \label{E1}
\end{eqnarray} 

These formulas   simplify if one takes into account that at frequencies
of our interest  the surface impedance is small:
\begin{equation}
|\zeta_{\alpha\beta}| \ll 1\,. \label{small_z}
\end{equation}
We then obtain instead of Eqs. (\ref{H1}) for the incidence angle not close to
$\pi/2$:
\begin{eqnarray}
H_{1x}&=&E_0 \cos\theta\, (1+2\zeta_{yy} \cos\theta) \, ,\nonumber\\ 
H_{1y}&=& - 2\,E_0 \,\zeta_{xy} \,, \label{H1_th}\\ 
H_{1z}&=&-E_0\sin\theta\, (1+2\zeta_{yy}\cos\theta) \,.\nonumber  
\end{eqnarray} 
and 
\begin{eqnarray}
  E_{1x}&=&-2E_0\zeta_{xy}\cos\theta\, ,\nonumber\\ 
  E_{1y}&=&-E_0(1+2\zeta_{yy} ),\label{E1_th}\\ 
  E_{1z}&=&2E_0\zeta_{xy}\sin\theta\,.\nonumber  
\end{eqnarray} 

For the grazing incidence, $\pi/2-\theta =\epsilon \ll 1$, we have:
\begin{eqnarray}
H_{1x}&=&E_0 \epsilon\,  ,\nonumber\\
H_{1y}&=& E_0 \,\frac{ 2\,\zeta_{xy}\,\epsilon}{\zeta_{xx} }
\,,\label{H1_eps}\\ 
H_{1z}&=&-E_0 \,  
\,.\nonumber  
\end{eqnarray} 

Equations (\ref{H1}) as well as (\ref{H1_th}) and (\ref{E1_th}) represent
elliptically polarized wave. Following the known procedure \cite{Field_Theory}
one can show that the ratio of the ellipse semi-axes is $2|\zeta_{xy}|$. Therefore,
 by recording the part of the reflected wave with the polarization
perpendicular to $y$, one can extract, in principle, information about the
superconductor anisotropy.  The intensity $I$ of this part depends on the
relative orientation of the plane $xz$ of the incident wave and the crystal,
$I\propto |\zeta_{xy}|^2\propto \sin^22\varphi$. 

Similar calculation gives for the complementary polarization ${\bm
H}_0=H_0(0,1,0)$ and ${\bm E}_0=H_0(-\cos\theta, 0,-\sin\theta)$: 
\begin{eqnarray}
H_{1x}&=& -\frac{2\cos^2\theta\,\zeta_{xy} }{(\cos\theta-\zeta_{xx})(1-
\zeta_{yy}\cos\theta)-\zeta_{xy}^2\cos\theta}\,H_0\,,\nonumber\\
H_{1y}&=&  \frac{(\cos\theta+\zeta_{xx})(1-
\zeta_{yy}\cos\theta)+\zeta_{xy}^2\cos\theta}{(\cos\theta-\zeta_{xx})(1-
\zeta_{yy}\cos\theta)-\zeta_{xy}^2\cos\theta}\,H_0\,,\nonumber\\
H_{1z}&=&-\tan\theta\,H_{1x}\,, \label{H1p} 
\end{eqnarray} 
and Eqs. (\ref{E1}) for ${\bm E}_1$ in terms of ${\bm H}_1$.  

 One easily checks that in the limit $\lambda\to 0$, cumbersome formulas
(\ref{H1}), (\ref{E1}), and (\ref{H1p}) give vanishing total fields ${\bm E}_t$ 
and $H_z$   as is the case for the normal metal with vanishing skin
depth.\cite{LL} Moreover, in this limit $E_z$ and ${\bm H}_t$ are the same in both
reflected and incident waves.  

We conclude this section by noting that the method of extracting the anisotropy
of the penetration depth from the data on reflection of polarized
low-frequency electromagnetic waves from anisotropic superconductors may provide
valuable information. Other known approaches to measure the $\lambda$ 
anisotropy are based on properties of vortices, and as such are bound to work
at high temperatures where the pinning is weak and the vortex system is in
equilibrium.  



\subsection{ Point sources}

Consider a point "magnetic charge" $\nu$ placed at $z=z_0$ above the
superconducting surface, formally the simplest possible source of the
external field. The potential $\varphi^s$ of this source is
\begin{equation}
\varphi^s=-\,\frac{\nu}{ \sqrt{r^2+(z-z_0)^2}}\,.
\label{phi_monopole}
\end{equation} 
Then, the 2D Fourier transform at $z=0$ reads:
\begin{eqnarray}
\varphi^s({\bm k})  &=& -\int \frac{\nu\, e^{i{\bm k}\cdot{\bm
r}}d{\bm r}}{\sqrt{r^2+z_0^2}}\nonumber\\
=-2\pi \nu && \int_0^{\infty}\,dr\,
\frac{rJ_0(kr)}{\sqrt{r^2+ z_0^2}}= -\frac{2\pi\nu}{ k}\,e^{-k z_0}
\end{eqnarray}
(see  Ref. \onlinecite{Grad}, 6.565). 
 
We can now study the response field $\varphi^r ({\bm r},z)$.  For
isotropic materials, $q_1=q_2=q=\sqrt{\lambda^{-2}+k^2}$, and we obtain
instead of (\ref{phi(k)an}):
 \begin{equation}
\varphi^r =\frac{q-k }{q+k}\,\varphi^s \,.
\label{phi(k)is}
\end{equation}
Therefore,
  \begin{eqnarray}
\varphi^r (r,z)  &=& -2\pi \nu\int \frac{d^2{\bm k}}{(2\pi)^2}\,\frac{q-k
}{k (q+k) }\,e^{- z_0k}\,e^{i{\bm k}\cdot{\bm r}-kz}\,\nonumber\\
 &=&-\nu\int_0^{\infty}dk\,\frac{q-k }{q+k
}\,e^{-k(z+z_0)}\,J_0(kr) \,.
\label{phi(r,z)is}
\end{eqnarray}
As in the general case, the relevant $k$'s here satisfy $k\leq
{\rm min}(1/r,1/z)$, and for  $r,z \gg\lambda$,  we can utilize 
$k\lambda \ll 1$: 
\begin{equation}
\frac{q-k }{q+k}\approx 1-2k\lambda\,.\label{kl<<1}
\end{equation}
Then, we obtain: 
 \begin{eqnarray}
\varphi^r (r,z) &=&  -\, \frac{\nu}{R} +2\lambda\,
\frac{\nu (z+z_0)}{ R^3}\,,\nonumber\\
 R &=& \sqrt{r^2+(z+z_0)^2}\,
\label{phi(r,z)is1}
\end{eqnarray}
(6.623 of Ref. \onlinecite{Grad} have been used). Clearly, the first term, 
which corresponds to $\lambda =0$, is the field of an image charge $\nu$ at
$z=-z_0$.  The first order correction due to a finite penetration depth
can be considered as the field of a magnetic dipole with a moment
$-2\lambda\nu {\hat z}$ placed at $z=-z_0$.  This result suggests that the
part of the response which carries the information about $\lambda$ is
quite weak. Indeed, the source image gives the part independent of
$\lambda$. To get the part of response  dependent on 
$\lambda$, one has to place a negative image displaced along $z$ with
respect to the image by 2$\lambda$.  

For the anisotropic case, the first term in Eq. (\ref{phi(small_k)an}) for
the response at distances large with respect to $\lambda$'s gives
$-\nu/R$, i.e.,  the field of the image at $z=-z_0$. In the second
term, the integration over the angle in the $k$ space is somewhat more
involved. The result is:
\begin{eqnarray}
\varphi^r (r,z) &=&  -\, \frac{\nu}{R} +(\lambda_c+\lambda_a)\,
\frac{\nu Z}{ R^3}\nonumber\\
&+&\nu (\lambda_c-\lambda_a)\cos 2\theta
\,\Big(\frac{2}{r^2}-\frac{2Z}{r^2R}+\frac{Z}{R^3}\Big)\,,
\end{eqnarray}
 where $Z=z+z_0$ and $\cos 2\theta =(x^2-y^2)/r^2$.



\subsection{$\lambda_x=\lambda_z\ne \lambda_y$}

This is the case of  the screening by the ``side surface" of a
uniaxial crystal like BSCCO (i.e., by the $ac$ plane). We can use the general
formulas (\ref{London_1})-(\ref{q's})  with the replacement
$\lambda_c\rightarrow\lambda_b$, and $\lambda_b\rightarrow\lambda_a$,
$\lambda_c\rightarrow\lambda_b$. After the replacement, we set
$\lambda_a= \lambda_b$ and use for both the common notation
$\lambda_{ab}$. We then obtain using Eq. (\ref{q's}):
\begin{equation}
q_1=\sqrt{\lambda_{ab}^{-2}+ k^2}\,,\quad 
q_2=\sqrt{1+\lambda_{ab}^2k_x^2+\lambda_c^2k_y^2}\Big/\lambda_c\,.
\label{q_1_2}
\end{equation}
These can be also obtained directly from 
Eqs. (\ref{FT1}) which now take the form:
\begin{eqnarray}
&&H_x\,(1+\lambda_{ab}^2k^2 -\lambda_{ab}^2q^2) = 0\,,\nonumber\\
&&H_x(\lambda_c^2-\lambda_{ab}^2)k_xk_y
+H_y(1+\lambda_{ab}^2k_x^2+\lambda_c^2k_y^2-\lambda_c^2q^2)=0\,,\nonumber
\\
&&i(k_xH_x +k_yH_y)+qH_z=0\,.   \label{FT2}
\end{eqnarray}
With $q=q_1$, this system gives
\begin{equation}
H_x^{(1)}=i\,\frac{q_1^2-k_y^2}{k_xq_1}\,H_z^{(1)}\,,\quad
H_y^{(1)}=i\,\frac{k_y}{q_1}\,H_z^{(1)}\,,\label{H1}
\end{equation}
and for $q=q_2$
\begin{equation}
 H_x^{(2)}=0\,,\qquad
H_y^{(2)}=i\,\frac{q_2}{k_y}\,H_z^{(2)}\,.\label{H1}
\end{equation}
Again, the two independent amplitudes $H_z^{(1)}$ and $H_z^{(2)}$  are found
along with $\varphi^r$ from the boundary conditions:
 \begin{eqnarray}
\varphi^r &=& \frac{k_x^2 q_1(q_2-q_1)   
+k(q_2-k)(k_y^2-q_1^2)}{(k_y^2-q_1^2)(k_y^2+kq_2)+k_x^2(k_y^2-q_1q_2)}\,
\varphi^s({\bm k})\,, \label{phi(k)an}\\
H_z^{(1)} &=&
\frac{q_1}{k_y^2-q_1q_2}\big[(k_y^2-kq_2)\varphi^s+
(k_y^2+kq_2)\varphi^r\big]\,,\label{Hz1}\\
H_z^{(2)} &=& \frac{k_y^2}{k_y^2-q_1q_2}\big[(k -q_1)\varphi^s-
(k+q_1)\varphi^r\big]\,.\label{Hz2}
\end{eqnarray}
Note that since $\varphi^s({\bm k})$ can be replaced with
$ h_z^s({\bm k})/k$, the response field can be expressed in terms of the
$z$ component of the source field at the interface.  







\bibitem{rem1} Take $F(z,\lambda)=$ const and differentiate it with respect to
$\lambda$: $F_{\lambda}+F_z\,z'=0$; this is Eq. (\ref{z'}). Differentiate it
again: $F_{\lambda\lambda}+2F_{\lambda z}\,z'+F_{zz}(z')^2+F_zz''=0$ and use 
$z'=-F_{\lambda}/F_z$ to obtain Eq. (\ref{z''}) for $z''$. The next
differentiation yields: $F_{\lambda\lambda\lambda}+3F_{\lambda\lambda z}\,z'
+3F_{zz\lambda}(z')^2+F_{zzz}(z')^3 +3F_{\lambda z}\,z'' +3 F_{zz}z'z''+F_zz'''=0$
in which again $z' $ and $z''$ are already known (in our case $z''=0$). Hence,
we obtain $z'''$.


For simplicity, consider the uniaxial material ($ m_a=m_b< m_c$)
with the $b$ axis perpendicular to the surface (i.e., $b=z$) and choose the
$x,y$ axes so that $x=c$ and $y=a$, see Fig. 1. In this situation, the mass
tensor is diagonal ($m_{xx}= m_c,\quad m_{yy}=m_{zz}= m_a$),
and Eqs. (\ref{general}) reduce to
\begin{eqnarray}
h_x &-&\lambda_a^2(h_{x,zz}+h_{x,yy}-h_{z,xz}-h_{y,yx})=0\,,\nonumber\\
h_y &-&\lambda_c^2(h_{y,zz}-h_{z,zy})-\lambda_a^2(h_{y,xx}-h_{x,yx})=0\,,
\label{system}\\
h_z &-&\lambda_c^2(h_{z,yy}-h_{y,yz})+\lambda_a^2(h_{x,xz}-h_{z,xx})=0\,,\nonumber
\end{eqnarray}
where   $\lambda_a^2=\lambda^2m_a$ and  $\lambda_c^2=\lambda^2m_c$. 

The London equations as written in Eq. (\ref{general}) were derived without
utilizing div${\bm h}=0$; doing this with the system (\ref{system}) we obtain:
\begin{eqnarray}
h_x &-&\lambda_a^2\nabla^2h_x =0\,,\label{system1}\\
h_y&-&\lambda_a^2(h_{y,xx}+h_{y,yy})-\lambda_c^2h_{y,zz}+
(\lambda_c^2-\lambda_a^2)h_{z,zy})=0\,,\nonumber\\
h_z&-&\lambda_c^2h_{z,yy}-\lambda_a^2(h_{z,xx}+h_{z,zz})+
(\lambda_c^2-\lambda_a^2)h_{y,yz})=0\,,\nonumber
\end{eqnarray}
This system has an advantage: the equation for $h_x$ separates from those for
$h_y$ and $h_z$.

After applying the $x,y$ Fourier transform to Eqs. (\ref{system}) (or to
(\ref{system1})), we are left with the linear system of ordinary 2nd order
differential equations with respect to the variable $z$ for $h_i({\bm k},z)$. The
solutions are linear combinations of simple exponentials:
\begin{equation}
h_i({\bm k},z) = \sum_n H_i^{(n)}({\bm k})e^{ q_nz} \,.
\label{exp}
\end{equation}
Substituting each term of this sum into the 2D Fourier transformed system
(\ref{system}) we obtain a linear homogeneous system for $H_i^{(n)}$
\begin{equation}
\Delta_{ij}^{(n)}H_i^{(n)}=0\, \label{hom_sys}
\end{equation}
with
\begin{eqnarray}
\Delta_{xx}^{(n)}&=&1+\lambda_a^2 k_y^2 -\lambda_a^2 q_n^2 \,,\nonumber\\ 
\Delta_{xy}^{(n)}&=&-\lambda_a^2 k_xk_y  \,,\nonumber\\ 
\Delta_{xz}^{(n)}&=&i\,\lambda_a^2 k_x  q_n  \,,\label{Delta}\\ 
\Delta_{yy}^{(n)}&=&1+\lambda_a^2 k_x^2 -\lambda_c^2 q_n^2 \,,\nonumber\\ 
\Delta_{yz}^{(n)}&=&i\,\lambda_c^2k_y q_n \,,\nonumber\\ 
\Delta_{zz}^{(n)}&=&1+\lambda_a^2 k_x^2 +\lambda_c^2k_y^2 \,.\nonumber 
  \end{eqnarray}
The determinant of the system (\ref{hom_sys}) must be zero; this provides
an equation for possible $ q_n$:\cite{KSL}
\begin{equation}
q_1=\frac{\sqrt{1 + \lambda_a^2k^2}}{\lambda_a}\,,\qquad
q_2=\frac{\sqrt{1 +
 \lambda_a^2k_x^2+\lambda_c^2k_y^2}}{\lambda_c} \,.\label{q's}
\end{equation}
 Note that the quantities $q$   must be positive
since the solution must vanish as $z\rightarrow -\infty$. 

For each of $q$'s, the system (\ref{hom_sys}) allows us to express two of
$H_i$'s in terms of the third. Following the standard procedure of solving
homogeneous linear systems we obtain for $n=1$:
 \begin{equation}
H_x^{(1)}=i\,\frac{q_1^2-k_y^2}{k_xq_1}\,H_z^{(1)},\qquad
H_y^{(1)}=i\,\frac{k_y}{q_1}H_z^{(1)}\,,
\label{H1}
\end{equation}
and for $n=2$
\begin{equation}
H_x^{(2)}=0\,,\qquad H_y^{(2)}=\,\frac{i\,q_2}{k_y}\,H_z^{(2)}\,. \label{H2}
\end{equation}
Therefore, to find the field inside the superconductor, we have to
determine only two independent quantities, $H_z^{(1,2)}$.  

These quantities along with the 
potential $\varphi^r$ for the response field outside can be  found from 
the boundary conditions of the field continuity at the interface $z=0$:
\begin{eqnarray}
h_x^s + i\,k_x\varphi^r &=& H_x^{(1)} + H_x^{(2)} = 
i\,\frac{q_1^2-k_y^2}{k_xq_1}\,H_z^{(1)}\,,\nonumber\\ 
h_y^s + i\,k_y\varphi^r &=& H_y^{(1)} + H_y^{(2)}
=i\,\frac{k_y}{q_1}H_z^{(1)} +i\,\frac{q_2}{k_y}\,H_z^{(2)}\,,
\label{bc}\\
 h_z^s - k\varphi^r &=& H_z^{(1)} +H_z^{(2)}\,.\nonumber
\end{eqnarray} 
We are interested only in the potential which is readily obtained:
 \begin{equation}
\varphi^r = \frac{i\,k_x h_x^s(k_y^2-q_1q_2)+i\,k_y h_y^s(k_y^2-q_1^2)
+q_2h^s_z(k_y^2-q_1^2)}{(k_y^2-q_1^2)(k_y^2+kq_2)+k_x^2(k_y^2-q_1q_2)}\,.
\label{phi(k)}
\end{equation}
This is the main result: knowing the unperturbed field at the
superconductor surface, we can evaluate the response field. 

As a consequence of div$\,{\bm h}^s=$ curl$\,{\bm h}^s=0$ div$\,{\bm
h}^s=0$ out of the source, the 2D Fourier components of ${\bm h}^s$ are
not independent. 
%To find the relation between the 2D Fourier components of
%${\bm h}^s$, we note that out of the sourse div$\,{\bm h}^s=$ curl$\,{\bm
%h}^s=0$ and, as with ${\bm h}^r$, 
As with the response field, we can look for ${\bm h}^s=$
grad$\,\varphi^s$ such that $\nabla^2\varphi^s=0$. In our case, 
 the source is situated in the upper half-space $z>0$, and we are
interested in the field ${\bm h}^s$ "under" the source. The general
solution of the Laplace equation which vanishes as $z\rightarrow -\infty$
is:
\begin{equation}
\varphi^s ({\bm r},z)=\int\frac{d^2{\bm k}}{(2\pi)^2}\,\varphi^s ({\bm k})\,
e^{i{\bm k}\cdot{\bm r}+kz}\,.
\label{Laplace_s}
\end{equation}
The 2D Fourier components of the field at the interface then read:
\begin{equation}
 h_x^s = i\,k_x \varphi^s ({\bm k})\,,\quad  h_y^s = i\,k_y \varphi^s
({\bm k})\,,\quad h_z^s =k\varphi^s ({\bm k})\,.\label{2dh}
\end{equation}
Therefore, we can express the response field $\varphi^r$ in terms of a
single function $\varphi^s$ instead of the three components of the source
field $h_i^s$:
%have 
%\begin{equation}
%i\,k_x h_x^s +i\,k_y h_y^s+k\,h_z^s =0\,.\label{divHs}
%\end{equation}

%We can, therefore, exclude one of $h^s_i$ from Eq. (\ref{phi(k)}):  
% \begin{equation}
%\varphi^r = \frac{i\,k_x q_1(q_1-q_2)h_x^s  
%+(q_2-k)(k_y^2-q_1^2)h^s_z}{(k_y^2-q_1^2)(k_y^2+kq_2)+
%k_x^2(k_y^2-q_1q_2)}\,.
%\label{phi(k)an}
%\end{equation}

 \begin{equation}
\varphi^r = \frac{k_x^2 q_1(q_2-q_1)   
+k(q_2-k)(k_y^2-q_1^2)}{(k_y^2-q_1^2)(k_y^2+kq_2)+k_x^2(k_y^2-q_1q_2)}\,
\varphi^s({\bm k}).
\label{phi(k)an}
\end{equation}
It is worth noting that since $\varphi^s({\bm k})$ can be replaced with
$k\,h_z^s({\bm k})$, the response field can be expressed in terms of the
$z$ component of the source field at the interface.  

 Equation (\ref{phi(k)an}) shows, in particular, that the total flux in
$z$ direction "reflected" by the superconductor is equal and opposite in
sign to the incident flux of the source crossing the interface. Indeed,
this flux is 
\begin{equation}
h_z^r|_{k=0}=-k\varphi^r|_{k=0} = -k \varphi^s|_{k=0}=-h_z^s|_{k=0} \,.
\end{equation}

The relevant $k$'s in the integral (\ref{Laplace_s}) obey $k\leq {\rm
min}(1/r,1/z)$. Therefore, for $r,z \gg\lambda$, in the relevant domain of
integration $k\lambda \ll 1$. Then, $q_{1,2}$ can be replaced with
$1/\lambda_{a,c}$, and the result (\ref{phi(k)an}) simplifies to
 \begin{equation}
\varphi^r({\bm k}) = \Big(1-2\frac{\lambda_ak_x^2 +
\lambda_ck_y^2}{ k }\Big)\, \varphi^s({\bm k})\,.
\label{phi(small_k)an}
\end{equation}
Note that the last expression cannot be used for $T\rightarrow T_c$ where 
$\lambda\rightarrow\infty$. We will consider this limit separately. 

\subsection{ Point sources}

Consider a point "magnetic charge" $\nu$ placed at $z=z_0$ above the
superconducting surface, formally the simplest possible source of the
external field. The potential $\varphi^s$ of this source is
\begin{equation}
\varphi^s=-\,\frac{\nu}{ \sqrt{r^2+(z-z_0)^2}}\,.
\label{phi_monopole}
\end{equation} 
Then, the 2D Fourier transform at $z=0$ reads:
\begin{eqnarray}
\varphi^s({\bm k})  &=& -\int \frac{\nu\, e^{i{\bm k}\cdot{\bm
r}}d{\bm r}}{\sqrt{r^2+z_0^2}}\nonumber\\
=-2\pi \nu && \int_0^{\infty}\,dr\,
\frac{rJ_0(kr)}{\sqrt{r^2+ z_0^2}}= -\frac{2\pi\nu}{ k}\,e^{-k z_0}
\end{eqnarray}
(see  Ref. \onlinecite{Grad}, 6.565). 
 
We can now study the response field $\varphi^r ({\bm r},z)$.  For
isotropic materials, $q_1=q_2=q=\sqrt{\lambda^{-2}+k^2}$, and we obtain
instead of (\ref{phi(k)an}):
 \begin{equation}
\varphi^r =\frac{q-k }{q+k}\,\varphi^s \,.
\label{phi(k)is}
\end{equation}
Therefore,
  \begin{eqnarray}
\varphi^r (r,z)  &=& -2\pi \nu\int \frac{d^2{\bm k}}{(2\pi)^2}\,\frac{q-k
}{k (q+k) }\,e^{- z_0k}\,e^{i{\bm k}\cdot{\bm r}-kz}\,\nonumber\\
 &=&-\nu\int_0^{\infty}dk\,\frac{q-k }{q+k
}\,e^{-k(z+z_0)}\,J_0(kr) \,.
\label{phi(r,z)is}
\end{eqnarray}
As in the general case, the relevant $k$'s here satisfy $k\leq
{\rm min}(1/r,1/z)$, and for  $r,z \gg\lambda$,  we can utilize 
$k\lambda \ll 1$: 
\begin{equation}
\frac{q-k }{q+k}\approx 1-2k\lambda\,.\label{kl<<1}
\end{equation}
Then, we obtain: 
 \begin{eqnarray}
\varphi^r (r,z) &=&  -\, \frac{\nu}{R} +2\lambda\,
\frac{\nu (z+z_0)}{ R^3}\,,\nonumber\\
 R &=& \sqrt{r^2+(z+z_0)^2}\,
\label{phi(r,z)is1}
\end{eqnarray}
(6.623 of Ref. \onlinecite{Grad} have been used). Clearly, the first term, 
which corresponds to $\lambda =0$, is the field of an image charge $\nu$ at
$z=-z_0$.  The first order correction due to a finite penetration depth
can be considered as the field of a magnetic dipole with a moment
$-2\lambda\nu {\hat z}$ placed at $z=-z_0$.  This result suggests that the
part of the response which carries the information about $\lambda$ is
quite weak. Indeed, the source image gives the part independent of
$\lambda$. To get the part of response  dependent on 
$\lambda$, one has to place a negative image displaced along $z$ with
respect to the image by 2$\lambda$.  

For the anisotropic case, the first term in Eq. (\ref{phi(small_k)an}) for
the response at distances large with respect to $\lambda$'s gives
$-\nu/R$, i.e.,  the field of the image at $z=-z_0$. In the second
term, the integration over the angle in the $k$ space is somewhat more
involved. The result is:
\begin{eqnarray}
\varphi^r (r,z) &=&  -\, \frac{\nu}{R} +(\lambda_c+\lambda_a)\,
\frac{\nu Z}{ R^3}\nonumber\\
&+&\nu (\lambda_c-\lambda_a)\cos 2\theta
\,\Big(\frac{2}{r^2}-\frac{2Z}{r^2R}+\frac{Z}{R^3}\Big)\,,
\end{eqnarray}
 where $Z=z+z_0$ and $\cos 2\theta =(x^2-y^2)/r^2$.

A similar problem can be solved for a more realistic source of a magnetic
moment $\vec{\mu}$ at the hight $z_0$ above a superconductor. If the
moment is directed along $z$, the potential of its field at $z=0$ is
\begin{equation}
\varphi^s=-\,\frac{\vec{\mu}\cdot \vec {R}}{R^3} = \frac{\mu\,z_0 }
 {(r^2+ z_0^2)^{3/2}}\,.
\label{phi_moment}
\end{equation} 
Then, the 2D Fourier transform is readily obtained:
\begin{equation}
\varphi^s({\bm k})  = 2\pi\mu \,e^{-k z_0}\,.
\label{FTphi}
\end{equation}
Using Eqs. (\ref{phi(k)is},\ref{kl<<1}) for the isotropic superconductors,
we find the potential of the response field:
 \begin{eqnarray}
\varphi^r (r,z) &=& \mu\Big( \frac{Z}{R^3} -2\lambda\,
\frac{2Z^2-r^2}{ R^5}\big)\,,\nonumber\\
 R &=& \sqrt{r^2+Z^2}\,,\qquad Z=z+z_0\,.
\label{phi(r,z)is2}
\end{eqnarray}
Here, the first term is the field of a moment $\vec{\mu}$ at
$z=-z_0$, i.e., the image source. The second term is  the field of a
proportional to $\lambda$ quadrupole at the same point. 

If the moment $\vec{\mu}$ is parallel to the interface, e.g., directed
along $x$, the source potential at $z=0$ reads
\begin{equation}
\varphi^s=-\,\frac{\mu\,x } {(r^2+ z_0^2)^{3/2}}\,.
\label{phi_moment_x}
\end{equation} 
It is straightforward to find the 2D Fourier transform
\begin{equation}
\varphi^s({\bm k})  = 2\pi\mu\, {k_x\over k}\,e^{-k z_0}\,
\end{equation}
and the response field:
 \begin{equation}
\varphi^r (r,z) = -\mu\,{x\over r}\,\Big( \frac{Z}{R^3} -2\lambda\,
\frac{2Z^2-r^2}{ R^5}\big)\,.
\label{phi(r,z)is3}
\end{equation}

\subsection{The source as a current loop}

According to (\ref{phi(k)an}), the response field everywhere can be
evaluated if the 2D Fourier transform of the unperturbed source field is
known. The latter is source specific. If the source is a loop of a current
$I$ of a size (or radius, if this a circle) $a$ situated in the plane
$z=z_0$, and if $a$ and $z_0$ are of the same order of magnitude
(practically, they are both of a  few-micron size),  modelling of the
source by a point-size magnetic moment does not suffice. The scalar
potential of the field created by a loop reads:\cite{LL}   
\begin{equation} 
\varphi^s ({\bm r}) = \frac{I}{c}\int\frac{d{\bm S}\cdot {\bm R}}{R^3}\,.
\label{LL}
\end{equation}
Here $d{\bm S}$ is an element of area, ${\bm R}$ is a radius-vector from
this element to the "observation point" ${\bm r}$ (in our case ${\bm
r}=(x,y,0)$), and the integral is over the area of the current contour.
The position of the element $d{\bm S}=d^2{\bm r}^{\prime}{\hat z}$ in our
situation is $({\bm r}^{\prime},z_0)$ so that
\begin{equation} 
 {\bm R} = {\bm r} - {\bm r}^{\prime} - z_0{\hat z}  \,,
\label{R}
\end{equation}
  and $R^2=({\bm r} - {\bm r}^{\prime})^2+z_0^2$. Then,
\begin{equation} 
\varphi^s ({\bm r}) = -\frac{Iz_0}{c}\int\frac{d^2{\bm r}^{\prime}  
 }{[({\bm r} - {\bm r}^{\prime})^2+z_0^2]^{3/2}}\,.
\label{LL1}
\end{equation}
For the circular loop of the radius $a$, the intergal here is over the
circle area. Comparing this with Eq. (\ref{phi_moment}) we see that the
source field can be considered as created by magnetic moments
distributed uniformily over the loop area with the density $-I/c$ so that
the total moment of the loop is $-\pi a^2I/c$. 

The 2D Fourier transform of $z_0[({\bm r} - {\bm
r}^{\prime})^2+z_0^2]^{-3/2}$ with respect to ${\bm r}$ is $2\pi
e^{-kz_0}e^{i{\bm k}\cdot{\bm r}^{\prime}}$ (see Eqs. (\ref{phi_moment},
\ref{FTphi})); therefore,
\begin{eqnarray} 
\varphi^s ({\bm k}) &=& -\frac{I}{c}\,2\pi e^{-kz_0}\int d^2{\bm
r}^{\prime} e^{i{\bm k}\cdot{\bm r}^{\prime}} \nonumber\\ 
&=& \frac{4\pi^2Ia}{ck} e^{-kz_0}J_1(ka) \,.
\label{ft_loop}
\end{eqnarray}
 
As an example, consider the response flux through the circular loop of the
source placed above the isotropic superconductor. Then, according to
(\ref{phi(k)is})
\begin{equation} 
\varphi^r ({\bm k}) =  
 \frac{4\pi^2Ia}{ck}\,\frac{q-k}{q+k}\, e^{-kz_0}J_1(ka) \,.
\label{ft_loop_r}
\end{equation}
and the $z$-component of the response field is
\begin{eqnarray} 
h_z^r ({\bm r},z) =  
 \frac{Ia}{c}\,\int  d^2{\bm k}\frac{q-k}{q+k}\, e^{-k(z+z_0)}J_1(ka)
e^{i{\bm k}\cdot{\bm r}}\nonumber\\
= \frac{2\pi Ia}{c}\,\int  dk\,k\,\frac{q-k}{q+k}\, e^{-k(z+z_0)}J_1(ka)
 J_0(kr)\,.
\label{ft_loop_hz}
\end{eqnarray}
Setting here $z=z_0$ and intergating over the area of the circle we obtain:
\begin{equation}
\Phi_z^r(\lambda,z_0;a)= \frac{4\pi^2 Ia^2}{c}\,\int_0^{\infty}  dk\,
\frac{q-k}{q+k}\, e^{-2k z_0}J_1^2(ka)\,,\label{38}
\end{equation}
the result of Clem obtained making use of the
cylindrical symmetry of the problem for the isotropic case. 

If the pick-up ring has a radius $a_p$ different from the radius $a_s$ of the
source, we obtain:
\FL\begin{eqnarray}
&&\Phi_z^r(\lambda,z_0;a_s,a_p)= \nonumber\\
&&\frac{4\pi^2Ia_sa_p}{c}\,\int_0^{\infty}  dk\,
\frac{q-k}{q+k}\, e^{-2k z_0}J_1(ka_s)J_1(ka_p)\,.\label{39}
\end{eqnarray}

\subsection{Discussion of the result (\ref{39})}

The flux $\Phi(\lambda,z_0;a)$ can be measured by SQUID to high precision.
Then, in principle, $\lambda$ can be evaluated for a given $z_0$ with the 
help of Eq.(\ref{39}). However, the hight $z_0$ is hard to  determine 
accurately. Hence, a high accuracy of $\Phi $ notwithstanding, one cannot
extract an accurate value of $\lambda$.

Physically, this difficulty is due to a simple fact that the same
change in the response can be  caused either by a change in $\lambda$
or, alternatively, by a variation of $z_0$.  Implicitly, Eq. (\ref{39})
determines only a certain function $z_0(\lambda)$. Numerical study suggests 
that this function is close to a simple linear dependence:
\begin{equation}
\lambda+z_0 \approx {\rm const}\,.\label{linear}
\end{equation}

We now show that this relation is valid indeed for $\lambda \ll z_0$,
i.e., in the practically relevant domain. To this end we consider the expansion
of the function $z_0(\lambda)$ near $\lambda=0$:
\begin{equation}
z_0 = L + A_1\lambda + A_2\lambda^2/2 + ... \,,\label{expansion}
\end{equation}
where $A$'s are corresponding derivatives evaluated at $\lambda=0$, while the
length $L$ is related to the response flux of an "ideal" Meissner surface with
$\lambda=0$. An equation for $L$ is obtained by setting $\lambda=0$ in Eq.
(\ref{39}):
\begin{eqnarray}
\frac{c\Phi_z^r}{4\pi^2Ia_sa_p}&=& \int_0^{\infty}dk\, e^{-2k
L }J_1(ka_p)J_1(ka_s )
\nonumber\\
&=&\frac{1}{\pi\sqrt{ a_p a_s}}\,Q_{1/2}(u)\,,\label{405}
\end{eqnarray}
where $Q$ is the Legendre function of the second kind and 
\begin{equation}
u=\frac{4L^2+a_p^2+a_s^2}{2 a_pa_s}\,.\label{u} 
\end{equation}
(formula  6.612.3 of Ref. \onlinecite{Grad} has been used). 
More convenient for the numerical work is the representation of $Q_{1/2}(u)$ in
terms of the hypergeometric function,\cite{Grad} which gives for Eq. (\ref{405}):
\begin{equation}
\frac{c\Phi_z^r\sqrt{2}}{\pi^2I\sqrt{a_sa_p}}=\frac{1}{u^{3/2}}\,\,
_2F_1\Big({3\over 4}, {5\over 4};2;{1\over u^2}\Big)\,.
\label{406}
\end{equation}

The coefficients $A$ can be obtained by using formulas for derivatives of
implicit functions.\cite{rem1} Denoting as $ F (z_0,\lambda)$ the integral
in the right-hand side of Eq. (\ref{39}) we have
\begin{equation}
A_1=\frac{dz_0}{d\lambda}\Big|_{\lambda=0}=-\frac{  F_{\lambda}}{F_z}
\Big|_{\lambda=0}\,,
\label{z'}
\end{equation}
where we use the notation $F_{\lambda}\equiv\partial
F/\partial\lambda$, $F_{zz}\equiv\partial^2 F/\partial z^2$ etc. and omit the
subscript 0 in $z_0$ for brevity. Further, we obtain:
\begin{eqnarray}
F_{\lambda}&=& -\int_0^{\infty}dk\, e^{-2kz_0}J_1(ka_p)J_1(ka_s)\frac{2
k}{\lambda^3q(q+k)^2}\Big|_{\lambda=0}\nonumber\\
&=& -2\int_0^{\infty}dk\,k\, e^{-2kL}J_1(ka_p)J_1(ka_s)\,,
\label{F_l}
\end{eqnarray}
and 
\begin{eqnarray}
F_z(0) &=& -2\int_0^{\infty}dk\,k\frac{q-k}{q+k} e^{-2kL}J_1(ka_p)J_1(ka_s)
\Big|_{\lambda=0}\nonumber\\
 &=& F_{\lambda}(0) \,.
\label{F_z}
\end{eqnarray}
Therefore, $A_1=-1$. 

Evaluation of $A_2$ is somewhat more involved, because the expression for the
second derivative of an implicit function reads:\cite{rem1}
\begin{equation}
\frac{d^2z_0}{d\lambda^2} 
=\frac{2F_{\lambda}F_zF_{z\lambda}-(F_{\lambda})^2F_{zz}-
 (F_z)^2F_{\lambda\lambda}}{(F_z)^3} \,.
\label{z''}
\end{equation}
Evaluating this at $\lambda=0$ we use $F_{\lambda}(0)=F_z(0)$:
\begin{equation}
\frac{d^2z_0}{d\lambda^2}\Big|_{\lambda=0} 
=\frac{2 F_{z\lambda}- F_{zz}-
  F_{\lambda\lambda}}{F_z}\Big|_{\lambda=0} \,. \label{z''}
\end{equation}
It is readily shown that the numerator of this expression vanishes. To see this,
it suffices to verify that
\begin{equation}
\Big[(\partial_z-\partial_{\lambda})^2\,\frac{q-k}{q+k}e^{-2kz}\Big]_{\lambda=0}=0\,.
\end{equation}
(with the substitution $\lambda k=\sinh u$, the expression in the square
parentheses reads
\begin{equation}
k^2e^{-2kz}\Big(2+\frac{1}{\cosh u}\,\partial_u\Big)^2e^{-2u}\,, 
\end{equation}
which makes the calculation easy). Thus, $A_2=0$. 

The procedure of evaluating higher order terms in the expansion (\ref{expansion})
, although increasingly cumbersome, can in principle be continued to any order. 
The third order term is still manageable because $z''(0)=0$ and
$F_z(0)=F_{\lambda} (0)$:\cite{rem1}
\begin{equation}
 z'''(0) 
=\frac{(\partial_z-\partial_{\lambda})^3
F }{F_z}\Big|_{\lambda=0} \,.
\label{z'''}
\end{equation}
After straightforward algebra, we obtain
\begin{equation}
z'''_0(0)={I_3\over I_1}\,,\quad I_n=\int_0^{\infty}dk\,k^n\,
e^{-2kL}J_1(ka_p)J_1(ka_s)\,,
\end{equation}
The integrals here are readily evaluated since $I_n=(-1/2)^nd^nI_0/dL^n$ and $I_0$
is given in Eq. (\ref{406}).

Thus we have
\begin{equation}
z_0 =  L\Big[1-{\lambda\over L} + 
\frac{i_3(L)}{6i_1(L)}\,\frac{\lambda^3}{L^3} +{\cal
O}\Big({\lambda^4\over L^4}\Big) \Big] \,,       \label{z(lam)}
\end{equation}
where
\begin{equation}
i_n=\int_0^{\infty}dt\,t^n\,
e^{-2t}J_1\Big(t{a_p\over L}\Big)J_1\Big(t{a_s\over L}\Big) \,.\label{i_n}
\end{equation}
The relevant part of the integration domain here is $t < 0.5$, which means that
$|i_3/i_1| < 1$. Therefore, the coefficient by $\lambda^3/L^3$ in is a small
number. 

For example, solving Eq. (\ref{406}) with $a_s=10\,\mu$, $a_p=5\,\mu$, and
$\Phi_z^rc/ I=2033\,\phi_0/$Amp, we obtain $L=1.82\,\mu$. Numerical evaluation
then yields $i_3(L)/6i_1(L) = -5.79\times 10^{-2}$.
 Since $\lambda/L\sim 10^{-1}$, the correction to linear dependence
$z_0=L-\lambda$ is of the relative order $10 ^{-4}$. 

\subsection{zzzzzz}

If a material has a sizable GL temperature domain where $\lambda (T)\propto
1/\sqrt{T_c-T}$, one can determine $z_0$ as follows:

Let $T_1$ and $T_2$ be two temperatures within this domain. Mesure the response
$\Phi_z^r$ at these $T$'s and evaluate corresponding $L(T)$'s with the help of
Eqs. (\ref{406}) and (\ref{u}). With the accuracy indicated above we have
\begin{eqnarray}
\lambda(T_1)=L(T_1)-z_0\,,\nonumber\\
\lambda(T_2)=L(T_2)-z_0\,.
\end{eqnarray}
Taking the ratio of these equations we obtain
\begin{equation}
z_0=\frac{L_1-L_2\alpha}{1-\alpha}\,,\qquad \alpha=\sqrt{{T_c-T_2\over
T_c-T_1}}\,.
\end{equation}



\bibitem{Abr} {\it Handbook of Mathematical Functions with
Formulas, Graphs, and Mathematical Tables}, edited by M. Abramowitz and
A. Stegun, Natl. Bur. Stand. Appl. Math. Ser. No. 55 (U.S. GPO,
Washington, D.C., 1965).


**********

% first rewrite 
% the integral (\ref{38}) using $t=k\lambda$ as the integration variable:
% \begin{equation}
% \Phi_z^r\frac{c}{4\pi^2Ia^2}={1\over\lambda}\int_0^{\infty}dt\,
% \frac{\sqrt{1+t^2}-t}{\sqrt{1+t^2}+t}\,
% e^{-2z_0t/\lambda}J_1^2\Big(t{a\over\lambda}\Big).
% \label{400}
%\end{equation}
 substitute the ansatz $z_0=L-\lambda$ in (\ref{39}):
\FL
\begin{eqnarray}
&&\frac{\Phi_z^rc}{4\pi^2Ia_sa_p}= \nonumber\\
&&\int_0^{\infty}dk\,
\frac{\sqrt{1+k^2\lambda^2}-k\lambda}{\sqrt{1+k^2\lambda^2}+k\lambda}\,
e^{2k\lambda}e^{-2 kL}J_1(k{a_s})J_1(k{a_p})\,,
\label{401}
\end{eqnarray}
  We now observe that in the relevant domain of integration
\begin{equation}
k\lambda \sim \lambda/z_0 \ll 1\,,
\end{equation}
 and one can expand the part of the integrand:
\begin{equation}
\frac{\sqrt{1+k^2\lambda^2}-k\lambda}{\sqrt{1+k^2\lambda^2}+k\lambda}\,
e^{2k\lambda}  =1 + {\cal O}(k^3\lambda^3)\,.
\label{402}
\end{equation}
The cancellation of linear and quadratic terms in this
expansion would not happen had the linear relation (\ref{linear})   been
taken as $z_0+C\lambda=L$ with $C\neq 1$. It is worth noting that this feature is
preserved for any shape of the source and the pick-up coil for isotropic
superconductors (as long as they are flat and parallel to the surface) because the
structure (\ref{402}) is inevitably present in expressions for the response field
at the pick-up coil. Therefore, the relation (\ref{linear}) is valid (within the
accuracy indicated) not only for concentric circular elements of the SQUID and the
field source. 




For $\lambda$ small relative to both $a$ and $z$, $(q-k)/(q+k)\approx
1-2k\lambda$. The integral (\ref{38}) then can be evaluated with the result
\begin{eqnarray}
\Phi_z^r(z_0,\lambda)&=& \frac{\pi^2
Ia}{4c}\,\frac{a^3}{z_0^3}\,\Big[F\Big({3\over 2}, {3\over 2};3;-{a^2\over
z_0^2}\Big)\nonumber\\ &-&3\,{\lambda\over z_0}\, F\Big({3\over
2}, {5\over 2};3;-{a^2\over z_0^2}\Big)\Big]\,.
\label{40}
\end{eqnarray}
For a given source size $a$, the response  depends on  $a/z_0$ and
$\lambda/z_0$. Using the known solution for the field of a circular source, one
can check that the first term is the field of the ring's image.\cite{LL}

 
It is not clear how this conclusion is changed for other possible
orientations of the source and SQUID relative to the surface. The simple linear
relation $z_0+\lambda\approx$const will not hold for anisotropic materials; it
can be seen from Eq. (\ref{phi(small_k)an}) that even the linear in
$k\lambda$ terms do not cancel. Hence, for the latter situation, the relationship
between $z_0$ and $\lambda$ which produce the same SQUID signal  should be
studied separately.

  Thus, provided $(\lambda/z_0)^3 \ll 1$, we have   
% \begin{equation}
% d= \int_0^{\infty}dk\, e^{-2 L k}J_1^2 ( ka
% )=\frac{a^2}{16L^3}F\Big({3\over 2},{3\over 2};3;-{a^2\over
% L^2}\Big)\, 
% \label{403}
% \end{equation}
One can say, therefore, that up to terms of the order $(\lambda/z_0)^3$, the
relation $z_0+\lambda=L$ satisfies indeed Eq. (\ref{39}). Moreover, for each data
point $\Phi_z^r$, one can solve Eq. (\ref{405}) numerically to find the length
$L(T)$.  

% For example, taking the low temperature value for $\Phi_z^r
% /I=2033\phi_0/{\rm Amp}$ and $a=5\,\mu$, we obtain $16da=3.411$; the equation
% \begin{equation}
% 3.411 =y^3F ( 3/2,3/2;3;-y^2 )\,,\qquad y=a/L\,
% \label{404}
% \end{equation}
% has a root $y=3.0388$ or $L= 1.645\,\mu$. 

The physical meaning of the lenght $L$ is quite simple. Consider the response
flux (\ref{39}) for $\lambda =0$ and compare it with Eq. (\ref{405}) for $L$. It
is seen that $L$ satisfies the same equation as the height $z_0$ for the case
$\lambda =0$. Hence, $L$ coincides with $z_0$ for a "pefectly reflecting"
superconductor. This statement follows, of course, from an approximate relation
$z_0+\lambda =L$. We then can interpret this relation as follows: the
outside magnetic response of a finite-$\lambda$ material at a distance
$z_0\gg\lambda$ is equivalent to the response of an ideal zero-$\lambda$
superconductor at a distance $z_0+\lambda$. In other words, a real material
reflects the field as if its surface were at the depth $\lambda$ and were
reflecting ideally. 

%The same argument as above gives for $L=z_0+\lambda$:\cite{Grad}

cond-mat0303381.extracted_bib
\bibitem{Kirtley} J.R. Kirtley,  C.C. Tsuei, M. Rupp,  J.Z. Sun, Lock See
Yu-Jahnes, A. Gupta, M.B. Ketchen, K.A. Moler, and M. Bhushan,
% {\it et al.} 
Phys. Rev. Lett. {\bf 76}, 1336 (1996).
 

\bibitem{KirtleyT} J.R. Kirtley, C.C. Tsuei,  K.A. Moler, V.G. Kogan, J.R.
Clem, A.J. Turberfield,
% {\it et al.} 
Appl. Phys. Lett. {\bf 74}, 4011 (1999).


\bibitem{Moler-Gardner} B.W. Gardner, J.C. Winn, D.A. Bonn, R. Liang, W.N. Hardy, J.R.
Kirtley, V.G. Kogan, and K.A. Moler,   
Appl. Phys. Lett. {\bf 80}, 1010 (2002).



\bibitem{Clem}J.R. Clem, M.W. Coffey, \prb {\bf 46}, 14662 (1992); M.W. Coffey, 
\prb {\bf 57}, 11648 (1998); M.W. Coffey, \prl {\bf 83}, 1648 (1999).


\bibitem{Milosevic}M.V. Milosevic, S.V. Yampolskii, and F.M. Peters, \prb, {\bf
66}, 174519 (2002).


\bibitem{Obukhov}Yu.V. Obukhov, J. Low Temp. Phys. {\bf 114}, 277 (1999).


\bibitem{K81}  V.G. Kogan, \prb, {\bf 24}, 1572 (1981).


\bibitem{KSL}  V.G. Kogan, A.Yu. Simonov, and M. Ledvij,
Phys. Rev. B {\bf 48}, 392 (1993).


\bibitem{Grad}I.S. Gradshteyn and I.M. Ryzhik, {\it Tables of
Integrals, Series, and Products}, Academic Press, 1980.

%
\bibitem{LL}L.D. Landau and E.M. Lifshitz, {\it Electrodynamics of
%Continuous Media}, Pergamon, Oxford, New York, 1984.

%
\bibitem{Field_Theory}L.D. Landau and E.M. Lifshitz, {\it The Classical Theory
%of Fields}, Pergamon, Oxford, New York, 1975.


\bibitem{John} J.R. Kirtley, unpublished.


\bibitem{rem1} Take $F(z,\lambda)=$ const and differentiate it with respect to
$\lambda$: $F_{\lambda}+F_z\,z'=0$; this is Eq. (\ref{z'}). Differentiate it
again: $F_{\lambda\lambda}+2F_{\lambda z}\,z'+F_{zz}(z')^2+F_zz''=0$ and use 
$z'=-F_{\lambda}/F_z$ to obtain Eq. (\ref{z''}) for $z''$. The next
differentiation yields: $F_{\lambda\lambda\lambda}+3F_{\lambda\lambda z}\,z'
+3F_{zz\lambda}(z')^2+F_{zzz}(z')^3 +3F_{\lambda z}\,z'' +3 F_{zz}z'z''+F_zz'''=0$
in which again $z' $ and $z''$ are already known (in our case $z''=0$). Hence,
we obtain $z'''$.


For simplicity, consider the uniaxial material ($ m_a=m_b< m_c$)
with the $b$ axis perpendicular to the surface (i.e., $b=z$) and choose the
$x,y$ axes so that $x=c$ and $y=a$, see Fig. 1. In this situation, the mass
tensor is diagonal ($m_{xx}= m_c,\quad m_{yy}=m_{zz}= m_a$),
and Eqs. (\ref{general}) reduce to
\begin{eqnarray}
h_x &-&\lambda_a^2(h_{x,zz}+h_{x,yy}-h_{z,xz}-h_{y,yx})=0\,,\nonumber\\
h_y &-&\lambda_c^2(h_{y,zz}-h_{z,zy})-\lambda_a^2(h_{y,xx}-h_{x,yx})=0\,,
\label{system}\\
h_z &-&\lambda_c^2(h_{z,yy}-h_{y,yz})+\lambda_a^2(h_{x,xz}-h_{z,xx})=0\,,\nonumber

\bibitem{Abr} {\it Handbook of Mathematical Functions with
Formulas, Graphs, and Mathematical Tables}, edited by M. Abramowitz and
A. Stegun, Natl. Bur. Stand. Appl. Math. Ser. No. 55 (U.S. GPO,
Washington, D.C., 1965).


**********

% first rewrite 
% the integral (\ref{38}) using $t=k\lambda$ as the integration variable:
% \begin{equation}
% \Phi_z^r\frac{c}{4\pi^2Ia^2}={1\over\lambda}\int_0^{\infty}dt\,
% \frac{\sqrt{1+t^2}-t}{\sqrt{1+t^2}+t}\,
% e^{-2z_0t/\lambda}J_1^2\Big(t{a\over\lambda}\Big).
% \label{400}
%