\section{Introduction}
While there is a general consensus that strong correlation governs the basic physics of high $T_c$ cuprates~\cite{LeeRevModPhys}, the role of electron-phonon interaction in determining $T_c$ is still under debate. The isotope effect is often viewed as an useful tool  which can provide information on this important issue. There is extensive isotope effect data on hole doped cuprates and the following picture has emerged. There exists significant isotope effect on the transition temperature for underdoped cuprates, but none for overdoped ones. On the other hand, for all doping substantial isotope effect on the London penetration depth $\lambda_{ab}$ has been observed~\cite{Khasanov03,Khasanov04}. Recall that $\lambda^{-2}_{ab}\propto n_S/m^*$ is a direct measure of the ratio between the superfluid density $n_S$ and the carrier mass $m^*$. The unusual isotope effect on $T_c$ can be understood qualitatively by the following picture~\cite{LeeRevModPhys}. For underdoped samples the transition temperature is controlled by the phase stiffness  $K_S=\hbar^2 n_S/4m^*$~\cite{Emery}. Hence, the isotope effect on $T_c$ may be simply inherited from the isotope effect on $n_S/m^*$. On the other hand, in overdoped samples $T_c$ is controlled by the pairing gap which has no isotope effect if phonons do not contribute significantly to its origin. Thus the isotope effect of $T_c$ can be qualitatively understood provided we accept the isotope effect on $n_S/m^*$. Then the puzzle is shifted to the origin of the isotope effect on $n_S/m^*$. Up to now there has been very little discussion in the literature on the isotope effect on $n_S/m^*$ in superconductors. In this paper we take the first step to address this important issue.

In~\cite{Khasanov03,Khasanov04} the authors suggested that the isotope effect on   $n_S/m^*$ is due to an isotope effect on the effective mass~$m^*$. However, the effective mass is given by $m^*/m=1+\lambda$ where~$\lambda$ is the dimensionless electron-phonon coupling.  Usually~$\lambda$ is not considered to have an isotope effect~\cite{McMillan}. This is because $\lambda=g^2N(0)\omega_D$ and $g^2\approx\corr{I^2}/\omega_D$, where $I$ is the coupling of the electron density to lattice deformation, which has no isotope dependence. This is supported by the direct experimental observation of Iwasawa \emph{et al.}~\cite{Iwasawa08} Using high-resolution laser ARPES measurements of electronic dispersion in Bi$_2$Sr$_2$CaCu$_2$O$_{8+\delta}$ (Bi2212), these authors study the influence of isotope oxygen substitution on the boson coupling ``kink'' in the electronic dispersion. The experiment clearly reveals an isotopic dependence of the kink energy, while it observes no change in the effective electron mass. On the other hand this experiment provides evidence for the coupling of phonon mode~\footnote{Authors of Ref.~[\onlinecite{Iwasawa08}] suggest that this is so-called breathing mode with $\Omega\sim69$~meV coupled to the electrons with coupling parameter $\lambda\sim0.6$ that can be estimated from mass renormalization.} to electrons, which leads us to propose alternative explanation for the isotope effect on $n_S$.

Our model disregards any effects of strong correlation and considers a clean superconductor with $s$ or $d$-wave singlet pairing. We suppose that there exists one phonon mode weakly coupled to the electron system. However, the superconducting order parameter is not related to the phonon mode, but rather induced by some other mechanism. Therefore, we ignore the isotope effect on the gap $\Delta$ and $T_c$, even though these effects can easily be added to lowest order in $\lambda$. We assume the standard adiabatic approximation for electron-phonon coupling~\cite{[{We note that there is earlier works that ascribe the isotope effect to a change in hole density by appealing to the breakdown of the adiabatic approximation~[}]Kresin94}\phantom{\cite{*[{}][{]. We consider this scenario unlikely to be applicable to cuprates.}]Kresin98}} \!\!and neglect Coulomb interaction. The only quantity which depends on the isotope mass is $\omega_E$.

First, we present simple qualitative arguments that explain the existence of the isotope effect on $n_S$ due to phonons using the sum rule which relates total spectral weight (SW) to the electron density and (bare) electron mass,
\begin{equation}
  \frac{ne^2}{m}
  =
  \frac{2}{\pi}\int_0^\infty d\omega\,\sigma'(\omega),
  \label{EqSumRule0}
\end{equation}
where $\sigma'(\omega)$ is the dissipative part of the conductivity. The onset of superconductivity does not change the total SW, but is redistributed~\cite{Ferrell}. Some part of the SW now goes to the $\delta$-function response of the condensate, whose weight determines $n_S/m^*$, while the remaining goes to frequencies $\omega\geq2\Delta$. This is visualized in the FIG.~\ref{FigQual}. It displays the dissipative part of conductivity for normal metal, $\sigma'(\omega)$, that consists of a narrow Drude peak (since our sample is clean) along with a contribution from phonons at some characteristic frequency $\omega_E$. Qualitatively with the onset of superconductivity some SW from frequencies $\omega<\Delta$ goes to the $\delta$-peak and the SW from frequencies $\omega>\Delta$ goes to the optical conductivity of superconductor. Thus, if we have $\Delta\ll\omega_E$ (the case in weak coupling BCS theory), variation of $\omega_E$ would have no effect on SW that goes to the condensate response. This may be the reason why this problem has not received much attention in the literature up to now. If $\Delta$ and $\omega_E$ have the same order of magnitude, $n_S/m^*$ responds to the change of $\omega_E$ and isotope effect for $n_S/m^*$ is present.

In the following sections we present calculations for the isotope effect for $n_S/m^*$ based on the sum rule. In Section~\ref{Sec2} we introduce our toy model, then calculate conductivity and superfluid density. In Section~\ref{Sec3} we discuss results and various approximations that were made. Finally, in the Appendix we present detailed calculations of the response kernel and conductivity.


\section{Calculation for $s$ and $d$-wave superconductors \label{Sec2}}
\begin{figure}[t]
\hspace*{0.02\linewidth}%
\includegraphics[width=0.999\linewidth]{fig0.eps}
\caption{\label{FigQual} Qualitative form of the conductivity for a normal system. With the onset of superconductivity, all spectral weight to the left of the corresponding red line goes to the condensate response.}
\end{figure}

In this section we present our calculation of the isotope effect on the superfluid density $n_S$. We note, that both the experiment on the London penetration depth and the theory described in Eq.~(\ref{EqSumRule0}) give the combination $n_S/m^*$. From hereon we simply define $n_S$ as a product of $n_S/m^*$ and $m^*$ which is taken to be constant. We calculate the dissipative part of the conductivity due to phonons. After this, having the expression for the conductivity via the sum rule, we get access to the superfluid density.

The theoretical model is a superconductor with electron-phonon interaction, described by the standard  Fr\"olich Hamiltonian. The superconducting order parameter (singlet pairing with $s$ or $d$ wave gap symmetry) is assumed to be \emph{not related} to phonon mode under consideration. One may think about this pairing as induced by other phonon modes, strong correlation or some different mechanism. The resulting Hamiltonian can be written~as:
\begin{equation}
  H
  =
  H_\text{BCS}+H_\text{ph}+H_\text{e-ph}.
\end{equation}
$H_\text{BCS}$ is the BCS Hamiltonian, with the pairing in Cooper channel treated within mean-field theory (Coulomb interaction is neglected):
\begin{multline}
  H_\text{BCS}
  =
  \sum_\mathbf{p}\xi_\mathbf{p} \Psi^+_\mathbf{p}\tau_3 \Psi_\mathbf{p}
  \\
  -
  \frac12\sum_{\mathbf{p}}
  \Psi^+_\mathbf{p} (\Delta_\mathbf{p}\tau_+ + \Delta^*_\mathbf{p} \tau_-)\Psi_\mathbf{p},
\end{multline}
where $\xi_\mathbf{p}$ is the Bloch energy measured relative to the Fermi energy (we consider only one electron band that is coupled to one phonon mode), $\tau^3$ is the Pauli matrix. Electron field operators here are written as a vector in the Nambu space
\begin{equation}
  \Psi_\mathbf{p}
  =
  \left(
  \begin{matrix}
  c_{\mathbf{p}\uparrow}\\
  c^+_{-\mathbf{p}\downarrow}
  \end{matrix}
  \right),
\end{equation}
with $c_{\mathbf{p},\uparrow}$, $c^+_{\mathbf{p},\uparrow}$ being standard creation and annihilation operators for up (down) spin electrons.
The Hamiltonian that describes the phonon mode and electron-phonon interaction is written as
\begin{multline}
  H_\text{ph}+H_\text{e-ph}
  =
  \sum_{\mathbf{k}} \omega(\mathbf{k}) b^+_\mathbf{k}b_\mathbf{k}
  \\
  +
  \sum_{\mathbf{p},\mathbf{p}'} g(\mathbf{p},\mathbf{p}')
  (b_\mathbf{\mathbf{p}'-\mathbf{p}}+b^+_\mathbf{\mathbf{p}-\mathbf{p}'})
  \Psi^+_{\mathbf{p}'}\tau_3\Psi_{\mathbf{p}}.
  \label{EqHFrolich}
\end{multline}
Here $\omega(\mathbf{k})$ is  the (bare) phonon dispersion and  $g(\mathbf{p},\mathbf{p}')$ is the  electron-phonon coupling. In what follows we use the Einstein approximation for the dispersion of phonons
\begin{equation}
  \omega(\mathbf{k})=\text{const}=\omega_E,
\end{equation}
and assume that our characteristic frequency $\omega_E$ is on order of the superconducting gap, i.e. $\omega_E\sim\Delta$ (note that for BCS superconductivity $\omega_E\gg\Delta$). Moreover, we limit ourselves to leading order corrections in the electron-phonon coupling that is assumed to be weak. Thus we can neglect not only renormalization of the electron-phonon coupling~(that is small in the parameter $\omega_E/E_F$~\cite{ Migdal}), but renormalization of the phonon dispersion as well.


\subsection{Green's function and self-energy\label{Sec2A}}
\begin{figure}[t]
\hspace*{0.02\linewidth}%
\includegraphics[width=0.999\linewidth]{fig1.eps}
\caption{\label{FigDiag} Diagrams corresponding to the self-energy~(a) corrections and conductivity~(b).}
\end{figure}
Before we start our calculations of conductivity we need to consider the electron Green's function. Following standard procedure, we define the Green's function of electrons $G_\mathbf{p}(\tau)$~(which is a matrix in Nambu space) and phonons $D_\mathbf{q}(\tau)$ in imaginary time as
\begin{gather}
  G_\mathbf{p}(\tau)=-\corr{T_\tau(\Psi_\mathbf{p}(\tau)\Psi_\mathbf{p}^+(0))},
  \\
  D_\mathbf{q}(\tau)=-\corr{T_\tau(b^+_\mathbf{p}(\tau)b_\mathbf{p}(0))}.
\end{gather}
Fourier transform of bare Green's functions (which are denoted by index $(0)$) can be written in the form~\cite{AGD}:
\begin{gather}
  (G^0_\mathbf{q}(i\omega_m))^{-1}
  =
  i\omega_m - \xi_\mathbf{q} \tau_3 - \Delta \tau_1,
  \label{EqGSC0Def}
  \\
  D^0_\mathbf{q}(i\omega_\nu)
  =
  \frac{1}{i\omega_\nu-\omega(\mathbf{q})}
  -
  \frac{1}{i\omega_\nu+\omega(\mathbf{q})},
\end{gather}
where $i\omega_m=(2n+1)\pi T$ and $i\omega_\nu=2\nu \pi T$
are the Matsubara frequencies for fermions and bosons.

In order to find the Green's function of electrons in the presence of the electron-phonon interaction we note that the superconducting order parameter is not related to the phonon mode under consideration. Moreover, the electron-phonon coupling is small and we limit ourselves to the leading order correction. Therefore, we can use perturbation theory to get access to the electron Green's function and there is no need in doing the full Eliashberg strong coupling theory~\cite{Eliashberg-60,*Eliashberg-61}.  Renormalization of the Green's function is given by the self-energy $\Sigma_\text{ph}$ due to electron-phonon interaction:
\begin{multline}
  (G^R_\mathbf{p}(\omega))^{-1}
  =
  (G^{0R}_\mathbf{p}(\omega))^{-1}-\Sigma_\text{ph}^R(\mathbf{p},\omega)
  \\
  =
  {\tilde\omega}^R(\omega) - \xi_\mathbf{p} \tau_3 - \tilde\Delta^R(\omega) \tau_1,
  \label{EqGSCRenDef}
\end{multline}
where $\tilde\omega (\omega)$ and $\tilde \Delta (\omega)$ are the renormalized~$\omega$ and gap~$\Delta$.

It suffices to calculate the self-energy to the \emph{leading order} in electron-phonon coupling. The expression for the self-energy in the Matsubara diagram technique can be read from the diagram depicted in FIG.~\ref{FigDiag}~(a):
\begin{multline}
  \Sigma_\text{ph}(\mathbf{p},i\omega_n)
  =
  -T\sum_m
  \int(dp')
  |g(\mathbf{p},\mathbf{p}')|^2
  \\
   D^0_{\mathbf{p}-\mathbf{p}'}(i\omega_n-i\omega_m)\tau^3 G^0_{\mathbf{p}'}(i\omega_m)\tau^3.
   \label{EqSigma1}
\end{multline}
Here the Pauli matrix $\tau_3$ (note that the self-energy is a matrix in the Nambu space) comes from the interaction vertex. Since we work to leading order of perturbation theory, we use the bare electron and phonon Green's functions. We proceed with analytical continuation~(the detailed procedure is described in \cite{Eliashberg-61,Scalapino}). Sums over positive and negative $m$ are represented as contour integrals. Afterwards, employing the spectral representation for the electron and phonon Green's functions~(see book~\cite{Scalapino}), one has:
\begin{multline}
  \Sigma^{R,A}_\text{ph}(\mathbf{p},\omega)
  =
  -\int_{-\infty}^\infty (dx) \int_{0}^\infty d\Omega\, \alpha^2(\Omega)F(\Omega)
  \\
  \int d\xi_{\mathbf{p}'}
  \mathop{\rm Im}[\tau^3 G^{0R}_{\mathbf{p}'}(x)\tau^3]
  \\
  \times
  \left[
  \frac{\tanh\frac{x}{2T}-\coth\frac{\Omega}{2T}}{x-\omega-\Omega \mp i\delta}
  -
  \frac{\tanh\frac{x}{2T}+\coth\frac{\Omega}{2T}}{x-\omega+\Omega \mp i\delta}
  \right],
  \label{EqSigma2}
\end{multline}
where $\alpha^2(\Omega) F(\Omega)$ is the spectral function of electron-phonon interaction (Eliashberg function) that consists of the product of effective electron-phonon coupling and phonon density of states. It is expressed through the electron-phonon interaction vertex as
\begin{multline}
  \alpha^2(\Omega)F(\Omega)
  =
  \frac{\int_{FS} d^2p
  \int_{FS} \frac{d^2p'}{(2\pi)^3v'_F} g(\mathbf{p},\mathbf{p}')\delta(\Omega-\omega_{\mathbf{p}-\mathbf{p}'})}{\int_{FS}d^2p},
  \label{EqalphaF}
\end{multline}
where $v_F'$ is the Fermi velocity and integrations are over the Fermi surface. Finally, we integrate over the loop momentum~\footnote{Account for $\xi_\mathbf{p}$ in the numerator of Green's function in the leading order gives the correction to the chemical potential that is small by parameter $\Delta^2/\epsilon_F^2$, see Ref.~[\onlinecite{Scalapino}].}:
\begin{multline}
  \int d\xi_{\mathbf{p}'}\, \mathop{\rm Im}[\tau^3 G^{0R}_{\mathbf{p}'}(x)\tau^3]
  =
  -i\pi\mathop{\rm sign} x\frac{x+\tau^1 \Delta}{\sqrt{x^2-\Delta^2}},
\end{multline}
and put temperature $T=0$. Expanding $\Sigma^{R,A}_\text{ph}(\mathbf{p},\omega)$ in components we obtain for renormalized $\omega$ and $\Delta$, that were defined earlier in Eq.~(\ref{EqGSCRenDef}),
\begin{multline}
  {\tilde \omega}^{R,A} (\omega)
  =
  \omega
  +
  \int_{-\infty}^\infty
  dx\,
  d\Omega\,
  \alpha^2F(\Omega)
  \mathop{\rm Re}\left[
  \frac{x \mathop{\rm sign} x}{\sqrt{x^2-\Delta^2}}
  \right]
  \\
  \left[
  \frac{\theta(-x)}{x-\omega-\Omega \mp i\delta}
  +
  \frac{\theta(x)}{z-\omega+\Omega \mp i\delta}
  \right]
  ,
  \label{EqZR}
\end{multline}
\begin{multline}
  {\tilde \Delta}^{R,A} (\omega)
  =
  \Delta
  +
  \int_{-\infty}^\infty
  dx\,
  d\Omega\,
  \alpha^2F(\Omega)
  \mathop{\rm Re}\left[
  \frac{\Delta\mathop{\rm sign} x}{\sqrt{x^2-\Delta^2}}
  \right]
  \\
  \left[
  \frac{\theta(-x)}{x-\omega-\Omega \mp i\delta}
  +
  \frac{\theta(x)}{z-\omega+\Omega \mp i\delta}
  \right]
  .
  \label{EqDeltaR}
\end{multline}

It is convenient to introduce a dimensionless coupling constant $\lambda$ defined to be responsible for the electron mass renormalization due to phonons:
\begin{equation}
  m^*
  =
  m(1+\lambda+O(\lambda)),
  \label{EqLambdaDef}
\end{equation}
where $m$ is bare and $m^*$ is renormalized (measured in experiments, e.g. in ARPES~\cite{Iwasawa08}) mass. The dimensionless coupling is expressed through the Eliashberg function~(\ref{EqalphaF}) as:
\begin{equation}
  \lambda
  =
  2
  \int \frac{d\omega}{\omega} \alpha^2(\omega) F(\omega).
  \label{EqlambdaDef}
\end{equation}
We work in the Einstein approximation, where the dispersion of phonon mode does not depend on the momentum, $\omega_\mathbf{p}=\omega_E$. Using this, one can immediately infer from Eq.~(\ref{EqalphaF}) that Eliashberg function $ \alpha^2(\Omega) F(\Omega)\propto\delta(\Omega-\omega_E)$. The constant of proportionality can be read from  Eq.~(\ref{EqlambdaDef}), giving us the expression for the Eliashberg function that will be used in the remainder of this paper:
\begin{equation}
  \alpha^2(\Omega) F(\Omega)
  =
  \frac{\omega_E\lambda}{2}
  \,
  \delta(\Omega-\omega_E).
\end{equation}
In what follows we will make use of the imaginary part of the renormalized frequency and gap, that can be easily inferred from Eqs.~(\ref{EqZR})-(\ref{EqDeltaR}):
\begin{gather}
  {{\tilde \omega}^R}{}''(\omega)
  =
  \pi\omega_E\lambda  \frac{\theta(|\omega|-\omega_E-\Delta)(|\omega|-\omega_E)}{2\sqrt{(|\omega|-\omega_E)^2-\Delta^2}},
  \label{EqtildeomIm}
  \\
  {\tilde \Delta}^{R}{}''(\omega)
  =
  \pi\omega_E\lambda \frac{\theta(|\omega|-\omega_E-\Delta)\Delta \mathop{\rm sign}\omega}{2\sqrt{(|\omega|-\omega_D)^2-\Delta^2}},
  \label{EqtildeDIm}
\end{gather}
whereas, for advanced functions  we have
\begin{equation}
  {{\tilde \omega}^A}{}''(\omega)
  =
  -{{\tilde \omega}^R}{}''(\omega),
  \qquad
  {\tilde \Delta}^{A}{}''(\omega)
  =
  -{\tilde \Delta}^{R}{}''(\omega).
\end{equation}
Not surpisingly, the imaginary part of the self-energy is present only for $\omega\geq\omega_E+\Delta$, i.e. where the real exitation can be created. The square root singularity that is present at this threshold would be smeared for a more realistic phonon spectrum.
\subsection{Conductivity\label{Sec2B}}
Having the electron Green's functions at hand, we proceed to the calculations of the conductivity. Following the standard approach~\cite{AGD}, in order to calculate the conductivity we consider the response kernel $Q_{\alpha\beta}(\mathbf{k},\omega)$ that relates current response to the vector potential:
\begin{equation}
  j_\alpha(\mathbf{k},\omega)
  =
  -\frac{ne^2}{m}Q_{\alpha\beta}(\mathbf{k},\omega) A_\beta(\mathbf{k},\omega).
  \label{EqQDefine}
\end{equation}
The optical conductivity is expressed through $Q$ as:
\begin{equation}
  \sigma_{\alpha\beta}(\mathbf{k},\omega)
  =
  \frac{ne^2}{m}
  \frac{Q_{\alpha\beta}(\mathbf{k},\omega)}{i\omega}
  .
  \label{EqSigmaThroughQ}
\end{equation}
The response kernel $Q_{\alpha\beta}(\mathbf{k},\omega)$ can be calculated within the Kubo linear response method as a sum of the diamagnetic contribution and a current-current correlator~\cite{AGD}:
\begin{equation}
  \frac{ne^2}{m}
  Q_{\alpha\beta}(\mathbf{k},\omega)
  =
  \frac{ne^2}{m}\delta_{\alpha\beta}
   -
  {\cal P}^R_{\alpha\beta}(\mathbf{k},\omega)
  ,
\end{equation}
where ${\cal P}^R_{\alpha\beta}(\mathbf{k},\omega)$ is analytic continuation to the real frequency of the (Fourier transformed) current-current correlator
\begin{equation}
  {\cal P}_{\alpha\beta}(\mathbf{r}-\mathbf{r}',\tau-\tau')
  =
  \corr{T({\hat j}_{1\alpha}(\mathbf{r},\tau) {\hat j}_{1\beta}(\mathbf{r}',\tau'))}.
\end{equation}
The current operator $\hat {\mathbf{j}}_{1}(\mathbf{r})$ is defined as a paramagnetic part of the full current operator $ \hat {\mathbf{j}}(\mathbf{r})$:
\begin{multline}
 \hat {\mathbf{j}}(\mathbf{r},\tau)
  =
  \frac{ie}{2m}(\nabla_{\mathbf{r}'}-\nabla_{\mathbf{r}})_{\mathbf{r}'\rightarrow \mathbf{r} }
  c^+(\mathbf{r}')c(\mathbf{r})
  \\
  -
  \frac{e^2}{m}\mathbf{A}(r)  c^+(\mathbf{r})c(\mathbf{r})
  \equiv
  \hat {\mathbf{j}}_{1}(\mathbf{r},\tau)
  -
  \frac{e^2}{m} \mathbf{A}(\mathbf{r})  c^+(\mathbf{r})c(\mathbf{r})
  .
  \label{EqjDef}
\end{multline}
Since we consider the superconductor in the London limit, we can put $\mathbf{k}=0$. Denoting $Q_{xx}(0,\omega)\equiv Q(\omega)$ we proceed to the calculation of $Q(\omega)$.

The current-current correlator ${\cal P}^R_{\alpha\beta}(k,\omega)$ is given by the sum of all possible diagrams with 2 external vertices corresponding to current operators. If we neglect vertex corrections (see discussion in Section~\ref{Sec3}), in the leading order in the electron-phonon coupling we need to consider only the simplest diagram, FIG.~\ref{FigDiag}~(b). It is convenient to represent the diamagnetic contribution to $Q(\omega)$ through Green's functions of a normal metal~\cite{AGD,Nam-67}. Combining the diamagnetic term with the current-current correlator given by the diagram in FIG.~\ref{FigDiag}~(b) we have:
\begin{multline}
  Q(i\omega_n)
  =
  \frac12\mathop{\rm tr}\left[T\sum_{\omega_m}\int d\xi_\mathbf{p}\,
  (G_\mathbf{p}(i\omega_n+i\omega_m) G_\mathbf{p}(i\omega_m)
  \right.
  \\
  -
  \left.
  G^{\Delta=0}_\mathbf{p}(i\omega_n+i\omega_m) G^{\Delta=0}_\mathbf{p}(i\omega_m))
  \right].
  \label{EqQMatsubara}
\end{multline}
Starting from this expression we perform analytic continuation, and integrate over the loop momentum. Analytical continuation is done in the standard way and is presented in details in the literature~\cite{Nam-67,Scalapino} and in the Appendix. After lengthy but straightforward calculations we get the following expression for the imaginary part of the response kernel $Q(\omega)$ at zero temperature:
\begin{multline}
  Q''(\omega)
  =
  \frac12\mathop{\rm Im}
  \int_{\Delta-\omega}^{-\Delta} dz\,
  \left[
   -\frac{g^{AA}(z+\omega,z)-1}{\varepsilon^A(z+\omega)+\varepsilon^A(z)}
   \right.
   \\
   \left.
  +
    \frac{g^{RA}(z+\omega,z)+1}{\varepsilon^R(z+\omega)-\varepsilon^A(z)}
  \right].
  \label{EqQpp}
\end{multline}
Where $\varepsilon^{\alpha}(\omega)$ with $\alpha=R,A$ is defined as:
\begin{equation}
  \varepsilon^{\alpha}(z)=\mathop{\rm sign} z \sqrt{(\tilde \omega^{\alpha}(z))^2-(\tilde \Delta^{\alpha}(z))^2}
  ,
  \label{EqveRADef}
\end{equation}
and $g^{\alpha\beta}$ are the structure factors introduced in~\cite{Nam-67}:
\begin{equation}
  g^{\alpha\beta}(z_1,z_2)
  =
  \frac{{\tilde \omega}^\alpha(z_1) {\tilde\omega}^\beta(z_2) +\tilde \Delta^\alpha(z_1)\tilde \Delta^\beta(z_2)}{\varepsilon^\alpha(z_1)\varepsilon^\beta(z_2)}.
\end{equation}
Using the expressions for the imaginary part of $\tilde \omega^{R,A}$ and $\tilde \Delta^{R,A}$,  Eqs.~(\ref{EqtildeomIm})-(\ref{EqtildeDIm}), we expand to  leading order in the corrections to the Green's function. After some calculations shown in the Appendix we arrive at the final answer, that can be expressed via complete elliptic integrals of the second kind, $E(x)$ and $K(x)$:
\begin{multline}
 \sigma'(\omega)
  =
  \frac{ne^2}{m}
  \frac{\pi\lambda\omega_E\omega_-}{\omega^3}
  \\
  \left[
   E\left(\frac{\omega^2_-}{\omega^2_+}\right)
   \right.
   \left.
   -
   \frac{\Delta(\omega-\omega_E)}{\omega_+^2} K\left(\frac{\omega^2_-}{\omega^2_+}\right)
   \right]
  \label{EqSigma1Final}
  ,
\end{multline}
where we used shorthand notation
\begin{gather}
  \omega_\pm
  =
  \omega - \omega_E \pm 2\Delta.
\end{gather}
In case when there is no superconductivity in our system, putting $\Delta=0$ from~(\ref{EqSigma1Final}) we reproduce the well-known result for the leading order correction to the conductivity of normal metal due to phonons~\cite{Holstein}:
\begin{equation}
  \sigma'_\text{N,e-ph}(\omega)
  =
   \frac{n e^2}{m}
   \lambda
   \frac{\pi\omega_E (\omega-\omega_E)}{\omega^3}\theta(\omega-\omega_E).
  \label{EqSigmaD=0}
\end{equation}
$\sigma'(\omega)$ has been computed earlier by Allen~\cite{Allen-71}. However the second term in Eq.~(\ref{EqSigma1Final}) containing the elliptic function $K$ is missing in Allen's formula. This difference is negligible when $\Delta\ll\omega_E$ but gives a small noticeable correction when $\Delta\gtrsim\omega_E$.
\begin{figure}[t]
\hspace*{0.02\linewidth}%
\includegraphics[width=0.999\linewidth]{sigma.eps}
\caption{\label{FigSigma} Correction to the dissipative part of conductivity $\sigma'(\omega)$ due to phonons for normal metal (blue) and superconductor with $s$ or $d$ wave symmetry (red and green curves correspondingly) for $\omega_E=1.5\Delta$. Sigma normalized to $\sigma_0=\frac{ne^2}{m}\lambda$, Drude peak and $\delta$-function response of condensate are not shown.}
\end{figure}
\subsection{d-wave symmetry\label{Sec2C}}
The previous calculations can be easily generalized for the $d$-wave symmetry of pairing. We use the simplest possible model with the cylindrical Fermi surface with an axis parallel to the $c$-axis of a crystal. The order parameter is supposed to have $d_{x^2-y^2}$ structure. It can be written as
\begin{equation}
  \Delta(\mathbf{k})
  =
  \Delta f
  \quad
  \text{with}
  \quad
  f=\cos 2\theta,
\end{equation}
where $\theta$ is the angle between momentum (which lies in the $a-b$ plane) and $a$ axis.

First we recalculate our corrections to the electron Green's function. In contrary to the $s$-wave case, the gap $\Delta(\mathbf{k})$ is dependent on the direction of momentum. Therefore the previous answer for the $s$-wave case needs to be averaged over the Fermi surface. Averaging is defined as
\begin{equation}
  \corr{X}_{FS}
  =
  \int_0^{2\pi} (d\theta)\, X(|\cos2\theta|)
  \label{EqDefAveragingDwave}
  .
\end{equation}
Gap renormalization ${\tilde \Delta}^R(\omega)-\Delta = 0$ vanishes under averaging over angle, while the imaginary part of the renormalization of $\omega$ reads:
\begin{multline}
  {\tilde \omega}^{R}{}'' (\omega)
  =
  \frac{\pi\omega_E\lambda}{2} \mathop{\rm \theta}(|\omega|-\omega_E)
  \corr{\frac{|\omega|-\omega_E}{\sqrt{(|\omega|-\omega_E)^2-\Delta^2f^2}}}_{FS}
  \\
  =
  \frac{\pi\omega_E\lambda}{2} \mathop{\rm \theta}(|\omega|-\omega_E)\, \kappa\left(\frac{|\omega|-\omega_E}{\Delta}\right)
\end{multline}
where $\kappa(x)$ is defined as
\begin{equation}
  \kappa(x)
  =
  \frac2\pi
  \begin{cases}
  x K(x^2) \ \ \text{for} \ \ x\leq 1 \\
  K(x^{-2}) \ \ \text{for} \ \ x> 1
 \end{cases}
 .
\end{equation}

We repeat calculations of conductivity with $\tilde \Delta(\omega)=\Delta$ and the new expression for ${\tilde \omega}^{R}{}''(\omega)$. Correction to the conductivity to leading order in the electron-phonon coupling is
\begin{multline}
 \sigma'(\omega)
   =
 \frac{ne^2}{m}
  \frac{2}{\omega^3}
 \corr{ \int_{\Delta|f|-\omega}^{-\Delta|f|} dz\,
  \frac{z {\tilde \omega}^{R}{}''(z+\omega)}{\sqrt{z^2-\Delta^2}}
  }_{FS}.
\end{multline}
Interchanging averaging over FS and integration over frequency we have two different expressions for
$\omega_E\geq\Delta$ and $\omega_E<\Delta$. When $\omega_E\geq\Delta$ we obtain:
\begin{multline}
  \sigma'(\omega)
  =
  \frac{ne^2}{m}
  \frac{\pi\omega_E\lambda}{\omega^3}
  \mathop{\rm \theta}(\omega-\omega_E)
  \\
  \int_{\omega_E-\omega}^{0} dz\,
   \kappa\left(\frac{|z+\omega|-\omega_E}{\Delta}\right)
   \kappa\left(\frac{|z|}{\Delta}\right)
 .
\label{EqSigma1dwave}
\end{multline}
While for $\omega_E<\Delta$ we have a more complicated expression involving the elliptic integral of the first kind  $F(\sin\phi;k)$:
\begin{multline}
  \sigma'(\omega)
   =
  \frac{ne^2}{m}
  \frac{\pi\omega_E\lambda}{\omega^3}
  \theta(\omega-\omega_E)
  \\
  \left[
  \int_{\Delta-\omega}^{0} dz\,
   \kappa\left(\frac{|z+\omega|-\omega_E}{\Delta}\right)
   \kappa\left(\frac{|z|}{\Delta}\right)
   \right.
   \\
   \left.
   +
   \frac{2}{\pi}\int_{\omega_E-\omega}^{\Delta-\omega} dz\,
   \kappa\left(\frac{|z+\omega|-\omega_E}{\Delta}\right)
   F\left(\frac{z+\omega}{\Delta};\frac{\Delta^2}{z^2}\right)
 \right]
 .
 \label{EqSigma1dwave2}
\end{multline}

\subsection{Isotope effect for superfluid density $n_S$ \label{Sec2D}}
\begin{figure}[t]
\hspace*{0.02\linewidth}%
\includegraphics[width=0.999\linewidth]{ns.eps}
\caption{Correction to superfluid density for superconductor with $s$-wave~(blue) or $d$-wave~(red) symmetry.
\label{FignS}}
\end{figure}
Using our explicit results for the conductivity, by means of the sum rule, we get access to the superfluid density. Initially the sum rule is formulated as
\begin{equation}
  \frac{ne^2}{m}
  =
  \frac{2}{\pi}\int_0^\infty d\omega\,\sigma'(\omega),
  \label{EqSumRule}
\end{equation}
where $n$ and $m$ are mass and concentration of electrons, while $\sigma'(\omega)$ is the dissipative part of the conductivity.

Let us first apply the sum rule to our system when there is no superconductivity, i.e. $\Delta=0$. We assume that the mean free path is big, thus scattering time due to impurities  $\tau$ satisfies $\tau\omega_E\gg 1$. Then, in addition to the very narrow Drude peak at $\omega=0$, $\sigma_\text{D}(\omega)$, we have the small correction due to phonons, $\sigma_\text{N,e-ph}'(\omega)$, that starts at $\omega_E$~(see Eq.~(\ref{EqSigmaD=0})). The calculation of the spectral weight of the Drude peak from the sum rule goes as follows:
\begin{multline}
   \frac{2}{\pi}\int_0^\infty d\omega\,\sigma_\text{D}'(\omega)
   =
    \frac{ne^2}{m}
    -
    \frac{2}{\pi}\int_0^\infty d\omega\,\sigma_\text{N,e-ph}'(\omega)
    =
    \\
   = \frac{ne^2}{m} - \lambda \frac{ne^2}{m}
   =
  \frac{ne^2}{m^*}(1 + O(\lambda^2)),
  \label{EqSumRuleN}
\end{multline}
where $m^*=m(1+\lambda)$, see Eq.~(\ref{EqLambdaDef}). We see that the SW of $\sigma_\text{N,e-ph}'(\omega)$ does not depend on $\omega_E$ and the sum rule reproduces the correct answer for the spectral weight of the Drude peak with the renormalized mass.

With the onset of the superconductivity the picture does not change drastically. The dissipative part of the conductivity consists of the $\delta$-function with the weight given by $n_Se^2/m^*$~(response of the condensate) and the small contribution due to electron-phonon interaction $\sigma'(\omega)$~(calculated in Eq.~(\ref{EqSigma1Final}) and Eqs.~(\ref{EqSigma1dwave})-(\ref{EqSigma1dwave2}) for $s$ and $d$-wave pairing). The sum rule implies for $n_S$:
\begin{equation}
  n_S
  =
  n-\frac{2m}{\pi e^2}
  \int^\infty_{2\Delta} d\omega\, \sigma'(\omega).
  \label{EqSumRulenS}
\end{equation}
From Eqs.~(\ref{EqSigma1Final}) and (\ref{EqSigma1dwave}), the superfluid density $n_S$ can be easily computed. It depends on $\omega_E$, thus giving nonzero isotope effect on superconducting density. Results of the numerical computations of $(n_S-n)/\lambda$ as a function of $\omega_E/\Delta$ are shown in  FIG.~\ref{FignS}.

Based on the dependence of $n_S$ on $\omega_E$, an isotope coefficient  $\beta$  can be easily calculated. It is defined as~\cite{Khasanov03}
\begin{equation}
  \beta=-\frac{d\log \lambda_{ab}^{-2}}{d\log M_\text{O}}
  =
  \frac12\frac{d\log \lambda_{ab}^{-2}}{d\log \omega_E}
  \label{EqBetaDef},
\end{equation}
where the usual relation between the phonon characteristic frequency and the atomic mass $M_\text{O}$~(in experiment this is usually the atomic mass of oxygen) was assumed:
\begin{equation}
  \frac{\Delta \omega_E}{\omega_E}
  =
  -\frac12\frac{\Delta M_\text{O}}{M_\text{O}}.
\end{equation}
Results for $\beta$ are shown in the FIG.~\ref{FigBeta}. Notably, for both $s$ and $d$-wave symmetry, the isotope coefficient has maximum modulus when $\omega_E\sim\Delta$. While when $\omega_E\gg\Delta$, $\beta$ vanishes in either case. Such behavior naturally follows from our qualitative arguments, presented in the Introduction. However note that $\beta$ approaches zero from positive and negative values for the $s$-wave and $d$-wave pairing respectively.
\begin{figure}[t]
\hspace*{0.02\linewidth}%
\includegraphics[width=0.999\linewidth]{beta.eps}
\caption{\label{FigBeta} Isotope coefficient $\beta$ defined in main text as a function of $\omega_E$ in our model.}
\end{figure}

\section{Discussion \label{Sec3}}
\begin{figure*}[t]
\hspace*{0.02\linewidth}%
\includegraphics[width=0.999\linewidth]{fig5.eps}
\caption{\label{FigContour} Deformation of integration contour during analytical continuation. Branch cuts of Green's functions are shown in red and starts at $\omega_\Delta$ and $\omega_\Delta-i\omega_m$. When we do analytical continuation from $\omega_n \rightarrow -i\omega$, all four branch cuts lie on the real axis.}
\end{figure*}
Although our toy model proves that isotope effect on the superfluid density can exist due to phonons, it gives an effect of the opposite sign to the experimentally measured. For example in~\cite{Khasanov03} for YPrBaCuO  $\beta=0.38$ for $x=0.3$ and  $\beta=0.71$ for $x=0.4$, while our model predicts the isotope coefficient of order of $\beta \sim -0.1$ for $\lambda\sim1$.

The opposite sign of isotope coefficient leads us to a discussion of the limitations of our model.  We treated the electron-phonon interaction as being weak which is not really the case. As it can be inferred from the experiment~\cite{Iwasawa08,Note1}, the coupling constant is of order one. However, the qualitative features of our result should not change. The second approximation  was the neglect of vertex corrections. The justification of this has been given in~\cite{Allen-71}. Since for normal metal ladder diagrams give minor effects, we expect this to be the case for superconductors as well.

The most serious limitation of our model with respect to high $T_c$ superconductors is the neglect of strong correlation effects. The minimum requirement would be to reproduce the fact that $n_S$ is proportional to the doping hole concentration $x$ due to the proximity to the Mott transition, and not to the electron density $1-x$ as in our toy model.  Furthermore, from the point of view of the sum rule argument, it is known that strong correlation gives rise to an incoherent background at finite frequencies, the so-called mid-infrared peak, even before phonons are taken into account. This incoherent background is missing in our toy model. What is needed is clearly a model which includes both Coulomb repulsion and coupling to phonons, such as the Hubbard-Holstein model. There has been considerable progress in this difficult problem, as given in a recent review~\cite{Castellani}. It may be possible to extend the present model to include correlation effect in the future.
Even though our result does not agree with the experiment, the important message of the present paper is a proof of principle that isotope effect in $n_S/m^*$ is possible even if the pairing is not due to phonons.

Finally the isotope effect on the superfluid density discussed in this paper should be applicable to conventional superconductors in the clean limit with minor extensions. To linear order in $\lambda$, it is straightforward to include coupling to a distribution of phonon modes. Furthermore, we can easily include the isotope effect on the energy gap $\Delta$ in our consideration. In weak coupling BCS theory, $\Delta$ is proportional to $\omega_E$ and the ratio $\Delta/\omega_E$  has no isotope dependence. Then our theory trivially predicts no isotope effect on $n_S/m^*$.
Including $\mu^*=\mu/\log(\epsilon_F/\omega_E)$ in the $T_c$ formula will introduce some isotope dependence in $\Delta/\omega_E$. However, in most superconductors $\omega_E/\Delta\gg1$ and from FIG.~\ref{FigBeta} we see that the predicted effect is very small. Recently $s$-wave superconductors with relatively large $T_c$ and energy gap have been discovered. Examples are MgB$_2$ and doped fullerenes. While coupling to certain high frequency phonons may be responsible for the pairing, there exist in these materials lower frequency phonons with frequency $\omega_0$ which may bring us to the regime of intermediate $\omega_0/\Delta$. It will be interesting to search for the isotope effect on the penetration depth in these materials (provided they are in the clean limit). In addition, the new class of Fe-based superconductors have a rather large ratio of $\Delta/\omega_E$ and should exhibit the isotope effect on the penetration depth according to our theory.


