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

In the holographic approach to condensed matter physics there is great freedom in how a bulk theory is constructed. This has the advantage of allowing the construction of holographic models which capture a wide variety of physical phenomena. The drawback, of course, is that the detailed phenomenology of such constructions typically depends on the particular model employed. For this reason a desirable feature of any phenomenological holographic programme is the existence of some `universal result', at least for a particular class of bulk theories. A celebrated example is the ratio of shear viscosity to entropy density for strongly coupled, isotropic holographic plasmas arising from two derivative bulk Lagrangians \cite{Kovtun:2003wp, Buchel:2003tz}.

Holographic superconductivity on the other hand has no simple universal sector; in the minimal example of a single $U(1)$-charged bulk scalar field there is freedom to specify the potential and gauge field coupling. In a bottom-up approach, whilst the small-field expansion of these functions can be easily parameterised, the behaviour away from the critical temperature depends strongly on their nonlinear characteristics. Turning to a top-down approach one may hope to overcome such difficulties. However, as observed in \cite{Donos:2011ut}, different consistent sub-truncations of a particular top-down model exhibit radically different results for the behaviour of the superconducting order parameter with temperature. To understand if there are any generic features for the landscape of consistent holographic superconductors, one must delve deeper into larger families of consistent truncations. See \cite{Denef:2009tp, Donos:2011ut, Aprile:2011uq, DallAgata:2011aa} for some interesting work along these lines.

In this paper we take a step in this direction by considering a truncation of $\mathcal{N}=8$ $SO(6)$ gauged supergravity in 5D containing 20 scalars and 15 gauge fields. The truncation contains a three-charge family of black branes \cite{Behrndt:1998jd}, dual to a normal phase, about which we look for the emergence of new branches of superconductors. Such branches are indicated by the existence of marginal modes for the various charged scalar fields at finite temperature. This calculation can be viewed as a first step towards understanding the sequence of preferred phases of the full supergravity theory.

When cooled sufficiently the black branes we employ to model the normal phase become thermodynamically unstable.\footnote{At high temperatures in the regime of local thermodynamic stability there are indications that this system is only metastable \cite{Yamada:2006rx,Yamada:2008em,Hollowood:2008gp}. With the interpretation of the three R-charges as angular momenta of D3-branes spinning in the $S^5$ of the full 10D theory, this instability corresponds to the one-by-one fragmentation of D3-branes from the rotating stack. The superconducting branches we seek in this work are local instabilities of the black brane solutions, and whilst they may be subject to arguments vis-\`a-vis metastability, a study of the critical temperatures and chemical potential dependence is not directly affected by it. The family of black branes used here represents a convenient local minimum for our purposes.} Note that in many cases the stable black brane simply ceases to exist at low temperatures, and so the system is largely immune to any $T=0$ test of stability, such as Breitenlohner-Freedman (BF) bound violation for any near horizon $AdS_2\times \mathbb{R}^3$. We focus on superconducting instabilities that occur in the normal phase at sufficiently high temperatures where it is locally thermodynamically stable. We find that a superconducting black brane branch always exists in this region of the normal phase.

In addition, the black brane solutions we employ for the normal phase have three R-charges, which allows us to explore the dependence of the various superconducting phase transitions on the conjugate chemical potentials.  Some modes in the truncation are charged under a single $U(1)$ and do not couple to the other gauge fields at the linearised level. As such, turning on an additional chemical potential for another $U(1)$ provides an effective holographic description of a state with a density imbalance. An example in condensed matter physics is an imbalanced Fermi mixture in which there is a mismatch of chemical potentials for spin-up and spin-down constituents \cite{citeulike:2350785}. Phenomenological holographic models of such mixtures were investigated in bottom-up bulk and probe brane constructions in \cite{Erdmenger:2011hp, Bigazzi:2011ak}, where the additional $U(1)$ plays the role of up-down spin imbalance. One feature of these models is that the superconducting phase ceases to exist if the imbalance is too great.

Here we have an imbalance as a natural consequence of the three chemical potentials.
As the imbalance for a given superconducting mode is increased, its critical temperature decreases until it reaches the threshold of thermodynamic instability, at finite temperature. However, increasing the imbalance for one mode decreases the imbalance for another, and is a consequence of the $SO(6)$ symmetry. Surprisingly there is no set of chemical potentials for which the imbalance is simultaneously too great for all modes.

The paper proceeds as follows. In section \ref{sec:trunc} we present the consistent truncation and the three-chemical potential black brane solution we will be employing as the background normal phase. In section \ref{sec:fluctuations} we calculate the field equations for homogeneous fluctuations of the other fields in the truncation about this background that do not source gravitational fluctuations. The structure of the equations for the fluctuations and the instabilities found naturally leads to two different types of mode: modes that are continuously connected to charged modes are discussed in \ref{sec:minusmodes}, while those that are continuously connected to neutral modes are discussed in section \ref{sec:plusmodes}. In both cases we map out the regions in the space of chemical potentials which are unstable to each mode and detail their critical temperatures. We conclude in section \ref{sec:discussion}.

\section{Truncation and the normal phase}\label{sec:trunc}
In addition to gravity, the bosonic sector of $\mathcal{N}=8$ $SO(6)$ gauged supergravity in 5D \cite{Pernici:1985ju, Gunaydin:1984qu,Gunaydin:1985cu} contains 42 scalars, 15 $SO(6)$ gauge fields as well as twelve 2-form gauge potentials. The scalars are organised into irreps of $SO(6)$ which have masses $m^2L^2 = -4,-3,0$ in an $AdS_5$ vacuum, the lightest of which saturates the BF-bound. We will employ a further consistent  truncation of this theory down to 15 gauge fields and 20 of the BF-bound saturating scalars. In 10D ingredients this truncation corresponds to retaining only the metric and 5-form \cite{Cvetic:2000nc}.

This truncation provides a convenient starting point in our search for superconducting black brane branches. First, the truncation contains a three-parameter family of black brane backgrounds which we may employ as the normal phase. Second, only the lightest scalars remain and these should display a greater susceptibility to forming condensates in the black brane backgrounds.\footnote{Superconducting branches for special cases of the other scalar field sectors are investigated in \cite{Aprile:2011uq}.} Finally, we are able to make contact with special cases studied elsewhere at the nonlinear level.

The 20 scalars and 15 gauge fields can be packaged as a symmetric unimodular $SO(6)$ tensor $T_{ij}$ and an antisymmetric $SO(6)$ tensor $A_{ij}$, respectively.  The equations of motion for this consistent truncation may be derived from the following action, with the adjoint $SO(6)$ indices written in matrix form:
\begin{equation} \label{eq:action}
S = \frac{1}{16\pi G_5} \int d^5x \sqrt{-g} \left(R - \frac{1}{4} \text{tr}\, T^{-1} \left(D_\mu T\right) \,T^{-1} \left(D^\mu T\right) +\frac{1}{8} \text{tr}\, T^{-1} F_{\mu\nu} T^{-1} F^{\mu\nu} - V \right)
\end{equation}
We have omitted a Chern-Simons term that will not be important for the electric, homogeneous superconductors sought here.  The scalar potential is given by
\begin{equation}
V = \frac{g_c^2}{2}\left(2\text{tr} T^2 - (\text{tr} T)^2\right),
\end{equation}
with the $SO(6)$ gauge-covariant derivative $D_\mu T = \nabla_\mu T + g_c \left[A_\mu,T\right]$ and field strength $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu + g_c\left[A_\mu,A_\nu\right]$. We set the dimensionful gauge coupling $g_c=1$ throughout, which fixes the AdS length.

When written in 2$\times$2 blocks, the background three-charge asymptotically AdS black brane solution \cite{Behrndt:1998jd} takes the form
\begin{eqnarray}
T = T_0 &\equiv& \text{diag}(X_1 \mathds{1},X_2 \mathds{1} ,X_3 \mathds{1}) \label{eq:T0}
\\
A = A_0 dt &\equiv& \text{diag}(A_1 \sigma, A_2 \sigma, A_3 \sigma)dt \label{eq:A0}
\end{eqnarray}
where $\mathds{1}$ is the 2$\times$2 unit matrix and $\sigma \equiv i \sigma_2$.
The unimodularity condition on $T$ enforces $X_1X_2X_3=1$.  The black branes are parametrised by three charge parameters $q_I$ and take the form
\begin{equation}
\begin{gathered}\label{gathered}
ds^2 = -H^{-2/3} f dt^2 + \frac{H^{1/3}}{f} dr^2 + r^2 H^{1/3} d \vec{x}_3^2 \\
X_I = \frac{H^{1/3}}{H_I}, \quad A_I =  \sqrt{\frac{r_+^4 H(r_+)}{q_I}} \left( \frac{1}{H_I(r)} - \frac{1}{H_I(r_+)} \right)
\end{gathered}
\end{equation}
We have defined
\begin{equation}
H_I(r) = 1 + \frac{q_I}{r^2}, \quad H(r) = H_1 H_2 H_3, \quad f(r) = - \frac{r_+^4 H(r_+)}{r^2} + r^2 H(r)
\end{equation}
and $I=1,2,3$ labels each 2$\times$2 block. The horizon is at $r=r_+$; we use the coordinate freedom in $r$ to fix $r_+=1$ and thus $\kappa_I = q_I$ in the conventions of \cite{Son:2006em}.

The thermodynamics of these backgrounds  was originally  studied in \cite{Chamblin:1999tk,Cvetic:1999ne} in the context of holography, with the study of spinning D3-branes, their counterparts in  the full 10D theory, initiated in \cite{Gubser:1998jb,Cai:1998ji} and extended in \cite{Cvetic:1999rb,Harmark:1999xt}.

\subsection{Multiple black brane branches}\label{sec:branches}
We work in the grand canonical ensemble throughout.
The chemical potentials for this system can be expressed in terms of the charge parameters
\begin{equation}
\mu_I = \frac{\sqrt{q_I (1+q_1)(1+q_2)(1+q_3)}}{1+q_I},
\end{equation}
and are conjugate to the three R-charges.   In this ensemble the relevant thermodynamic potential is the Gibbs potential, $\Omega$, with density $\omega =\frac{ \Omega}{\mathrm{vol}_3}$. We employ a dimensionless expression $\hat{\omega}$:
\begin{equation}
\omega = -\frac{ (1+q_1)(1+q_2)(1+q_3)}{16 \pi G_5}, \qquad \hat{\omega} \equiv \frac{16\pi G_5\, \omega}{(\mu_1^2+\mu_2^2+\mu_3^2)^2}.
\end{equation}
Similarly we employ a dimensionless version of the temperature, $\hat{T}$:
\begin{equation}
T = \frac{2+q_1+q_2+q_3 -q_1 q_2 q_3 }{2\pi\sqrt{(1+q_1)(1+q_2)(1+q_3)}}, \qquad  \hat{T} = \frac{T}{\sqrt{\mu_1^2+\mu_2^2+\mu_3^2}}.
\end{equation}

It is well known that this black brane family is not thermodynamically stable everywhere. For a stable thermodynamic equilibrium, the theory should lie at a local minimum of the thermodynamic potential, at fixed temperature and chemical potential. Constructing the Hessian matrix for variations of the energy density $\varepsilon$,
\begin{equation}
\mathcal{H} = \frac{\partial^2 \varepsilon}{\partial s\partial \rho_I}
\end{equation}
with respect to the entropy density $s$ and the R-charge densities $\rho_I$ following \cite{Son:2006em}, a necessary condition for local thermodynamic stability is given by
\begin{equation}
\det\,\mathcal{H} \propto 2-q_1-q_2-q_3 + q_1 q_2 q_3 >0. \label{detHess}
\end{equation}
However note that this is not a sufficient condition for local thermodynamic stability as there may be an even number of negative eigenvalues of $\mathcal{H}$.\footnote{We wish to thank Mukund Rangamani for bringing this point to our attention.}
At a fixed set of chemical potentials, $(\mu_1,\mu_2,\mu_3)$, we find that there is a minimum temperature for the existence of a thermodynamically stable dominant branch of solutions.

Without loss of generality we focus on  $\mu_I\geq0$ at which we only consider those solutions with $q_I\geq0$. This ensures that the solutions are regular. For this ensemble we fix two dimensionless ratios of the chemical potentials, $\zeta_1 = \frac{\mu_1}{\mu_3}$ and $\zeta_2 = \frac{\mu_2}{\mu_3}$, or any function thereof, in addition to $\hat{T}$. In general there are four black brane solutions at any grand canonical coordinate, $(\zeta_1,\zeta_2, \hat{T})$. To explore the phases of the system we adopt the coordinates
\begin{equation}
p_1=\frac{\mu_1+\mu_2-2\mu_3}{\sqrt{6}\sqrt{\mu_1^2+\mu_2^2+\mu_3^2}},\qquad p_2= \frac{\mu_1-\mu_2}{\sqrt{2}\sqrt{\mu_1^2+\mu_2^2+\mu_3^2}}\label{pcoords}
\end{equation}
which are dimensionless quantities and depend only on $\zeta_{1,2}$.

At a representative fixed $(p_1,p_2)$ we show the free energy as a function of the temperature in the left panel of figure \ref{fig:freeenergy}, illustrating that the dominant line of black branes are only stable above a threshold temperature. This is the threshold for thermodynamic instability as described by the condition \eqref{detHess}, as well as the positivity of the eigenvalues of $\mathcal{H}$. In appendix  \ref{app:branch} we show that there is always one dominant branch of black branes with a minimum temperature given by the threshold for thermodynamic instability.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.45\textwidth]{freeenergygeneric}
\includegraphics[width=0.45\textwidth]{freeenergyequalq}
\setlength{\unitlength}{0.1\columnwidth}
\begin{picture}(0.1,0.25)(0,0)
\put(-6.,2.85){\makebox(0,0){$\hat{T}$}}
\put(-1.,2.85){\makebox(0,0){$\hat{T}$}}
\put(-9.5,1.8){\makebox(0,0){$\hat{\omega}$}}
\end{picture}
\caption{\emph{Left:} Gibbs free energy density for the various black brane branches that exist at fixed $(p_1,p_2)=(-0.2,-0.1)$, defined in \eqref{pcoords}. Solutions that are thermodynamically unstable are shown in blue. A locally thermodynamically stable branch has a minimum temperature marked by the dashed line. \emph{Right:} Fixed $p_1=p_2=0$. The equal-charge black brane is only thermodynamically stable above the temperature indicated by the dashed line, whilst below it is both thermodynamically unstable and subdominant in the grand canonical ensemble.
}
\label{fig:freeenergy}
\end{center}
\end{figure}

It is worth noting that in some cases there is a branch of black branes which satisfy the condition \eqref{detHess}, but which are not thermodynamically stable due to an even number of negative eigenvalues of $\mathcal{H}$. For example,  the three equal chemical potential case,  where $p_1=p_2=0$, is shown in the right panel of figure \ref{fig:freeenergy}.  The dominant black brane above the critical temperature $\hat{T}\simeq 0.09$ is locally thermodynamically stable and has equal charges $q_1=q_2=q_3$ and is nothing but the Reissner-Nordstr\"om-$AdS_5$ (RN) black brane. Below this critical temperature $\mathcal{H}$ picks up two equal negative eigenvalues with eigenvectors $\delta(s,\rho_1,\rho_2,\rho_3) = (0,-1,0,1),(0,-1,1,0)$. These eigenvectors for this thermodynamic instability are naturally associated with a Gubser-Mitra dynamical instability \cite{Gubser:2000ec, Gubser:2000mm}, as well as the generation of neutral zero modes at the threshold of stability which we discuss further in section \ref{sec:plusmodes}.

On a different tack, the free energy analysis presented here reveals a worrisome issue associated with the RN solution, as considered in this top-down approach. Working in the further consistent truncation of our theory corresponding to three equal charge parameters, RN becomes the only black brane solution and therefore appears to be thermodynamically dominant. However, below the critical temperature for the thermodynamic instability discussed above, $\hat{T}\simeq 0.09$, RN is not thermodynamically dominant in our more complete theory.

To understand the coordinates \eqref{pcoords} better, we project the condition $\mu_I\geq0$ into $(p_1,p_2)$ space. This region is formed by the intersection of three ellipses, with equations
\begin{equation}
3 p_1^2 + p_2^2 = 1, \quad 3 p_1^2 + 5 p_2^2 \pm  2\sqrt{3} p_1 p_2
 = 2
\end{equation}
and is shown in figure \ref{fig:background}. Within this region we find there is always a stable branch that dominates the ensemble above some critical temperature, as in the explicit examples shown in figure \ref{fig:freeenergy}. The minimum temperature for this dominant stable branch as a function of $p_1$ and $p_2$ is illustrated in the right panel of figure \ref{fig:background}.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.48\textwidth]{picknellipseRaster}
\includegraphics[width=0.48\textwidth]{detDensityRaster}
\setlength{\unitlength}{0.1\columnwidth}
\begin{picture}(0.1,0.25)(0,0)
\put(-7.3,-0.2){\makebox(0,0){$p_1$}}
\put(-10.1,2.4){\makebox(0,0){$p_2$}}
\put(-2.45,-0.2){\makebox(0,0){$p_1$}}
\end{picture}
\caption{\emph{Left:} Region where black brane backgrounds with $\mu_I\geq0$ exist.  \emph{Right:} Critical temperatures for the thermodynamic instability --- \emph{i.e.}\ the minimum temperature, $\hat{T}$, for which the dominant branch is thermodynamically stable. \label{fig:background}}
\end{center}
\end{figure}

\section{Fluctuations about the normal phase}\label{sec:fluctuations}
We now detail the fluctuations of the 20 scalars and 15 gauge fields about the background normal phase solution given in section \ref{sec:trunc}.
We parameterise deviations from the scalar field matrix $T_0$ \eqref{eq:T0} by the fluctuation matrix $\chi(r)$, and deviations about the gauge field matrix $A_0$ \eqref{eq:A0} with radial electric fluctuations, $\alpha(r)$:
\begin{equation} \label{eq:Tflucs}
T = T_0 e^{\chi}\, , \qquad \text{tr} \chi = 0\, ; \qquad A = (A_0 + \alpha)dt
\end{equation}
where the trace condition enforces unimodularity of $T$ at all orders in the fluctuation. The action for the linearised fluctuations can be found in appendix \ref{app:details}. Since the background is block diagonal, the fluctuations in different 2$\times $2 blocks do not couple.

In general the fluctuations $\chi$ and $\alpha$ source metric fluctuations. For general background $\{q_I\}$ these may be switched off by enforcing the conditions
\begin{equation}
\text{tr}\, \chi_{II} =  \text{tr}\, \sigma\alpha_{II} = 0 \label{diagonalconstraints}
\end{equation}
where the trace is performed within the 2$\times$2 block. In detail these conditions are so restrictive that we have $\alpha_{II} = 0$ and $\sigma\chi_{II} + \chi_{II}\sigma = 0$, leaving the only non-zero diagonal block fluctuation as $ \sigma\chi_{II} - \chi_{II}\sigma$, which as we shall see shortly is a charged mode. For special cases of the background the conditions relax: for two or three equal charges we find
\begin{equation}
\sum_I \text{tr}\, \chi_{II} = \sum_I \text{tr}\, \sigma\alpha_{II} = 0
\end{equation}
where the sum is performed over equal $q_I$ blocks, with any remaining diagonal blocks constrained by \eqref{diagonalconstraints}.

\subsection{Diagonal fluctuations}
Without loss of generality we focus on fluctuations living in the $(1,1)$ block,
\begin{equation}
\chi = \left(
\begin{array}{ccc}
\chi_{11} & 0 & 0 \\
0 & 0 & 0\\
0 & 0 & 0
\end{array}
\right),\qquad
\alpha = \left(
\begin{array}{ccc}
\sigma\hat{\alpha}_{11} & 0 & 0 \\
0 & 0 & 0\\
0 & 0 & 0
\end{array}
\right)\label{diagonal}
\end{equation}
where $\chi_{11}$ is a symmetric 2$\times$2 matrix and $\hat{\alpha}_{11}$ is proportional to the 2$\times$2 identity matrix. To satisfy the backreaction and unimodularity constraints away from equal charges we must impose $\text{tr} \chi_{11} = 0$.

The equations of motion for these fluctuations may be partially diagonalised by moving to the new variables
\begin{eqnarray}
\textbf{Neutral:}\quad\phi_{11}^{(+)} = \sigma \chi_{11} + \chi_{11} \sigma,&&\quad \alpha_{11}^{(+)} = \sigma \hat{\alpha}_{11} +\hat{\alpha}_{11} \sigma, \label{diagonalNewVars}\\
\textbf{Charged:}\quad\phi_{11}^{(-)} = \sigma \chi_{11} - \chi_{11} \sigma, &&\quad \alpha_{11}^{(-)} = \sigma \hat{\alpha}_{11} - \hat{\alpha}_{11} \sigma = 0.
\end{eqnarray}
The gauge potential for the $(-)$ mode vanishes due to the antisymmetry of $\sigma\hat{\alpha}_{11}$. Note that in what follows we have specialised to radial dependence in deriving these equations, despite their being presented in partially covariantised form.

\subsubsection{Charged equations of motion}
The equation of motion for the charged scalar $\phi_{11}^{(-)}$ is
\begin{equation}
\Box \phi_{11}^{(-)} -4 \left(X_1 (X_1-X_2-X_3)+ A_1^2\right)\phi_{11}^{(-)} = 0.\label{Weom}
\end{equation}
Various special cases of this mode were considered nonlinearly in \cite{Aprile:2011uq} for the one, two and three equal-charge background; for these cases we have numerically verified the critical temperatures quoted. The unimodularity and backreaction constraint $\text{tr} \chi_{11} = 0$ can be imposed without trivialising this mode.

\subsubsection{Neutral equations of motion}
The equations of motion for the neutral scalar $\phi_{11}^{(+)}$ together with the gauge field fluctuation are
\begin{eqnarray}
\Box \phi_{11}^{(+)} - 4X_1(X_1-X_2-X_3)\phi_{11}^{(+)}-\frac{F_1^2}{X_1^2}\phi_{11}^{(+)}+\frac{F_1\cdot h}{X_1^2}  = 0 \label{plusdiag1}\\
\frac{1}{\sqrt{-g}}\partial_r \left(\sqrt{-g}\;\frac{h^{rt}-2F_1^{rt}\phi_{11}^{(+)}}{X_1^2}\right)=0\label{plusdiag2}
\end{eqnarray}
where we have defined $h = d \alpha_{11}^{(+)}$.

If the mode is to be considered alone as in \eqref{diagonal}, we must impose $\text{tr} \chi_{11} = 0$ for which we have $\phi_{11}^{(+)}=0$ by \eqref{diagonalNewVars}. In order to satisfy the backreaction and unimodularity constraints without trivialising these modes, we must therefore include other diagonal blocks, but this can only achieved if two or more of the background charges are equal. Such a sum of modes at two equal charges is closely related to the Gubser-Mitra instability \cite{Gubser:2000ec, Gubser:2000mm} --- a connection that we explore further in section \ref{sec:plusmodes}.

\subsection{Off-diagonal fluctuations}
Without loss of generality we focus on fluctuations living in the $(1,2)$ block (and $(2,1)$ by symmetry),
\begin{equation}
\chi = \left(
\begin{array}{ccc}
0 & \chi_{12} & 0 \\
\chi_{12}^T & 0 & 0\\
0 & 0 & 0
\end{array}
\right),\qquad
\alpha = \left(
\begin{array}{ccc}
0 & \sigma\hat{\alpha}_{12} & 0 \\
\hat{\alpha}_{12}^T\sigma & 0 & 0\\
0 & 0 & 0
\end{array}
\right)\label{offdiagonal}
\end{equation}
for which the condition $\text{tr} \chi = 0$ and backreaction constraints are automatically satisfied.

As for the diagonal fluctuations, the equations of motion for these fluctuations may be partially diagonalised by moving to the new variables
\begin{eqnarray}
\textbf{($+$) mode:}\quad\phi_{12}^{(+)} = \sigma \chi_{12} + \chi_{12} \sigma,&&\quad \alpha_{12}^{(+)} = \sigma \hat{\alpha}_{12} + \hat{\alpha}_{12} \sigma,\\
\textbf{($-$) mode:}\quad\phi_{12}^{(-)} = \sigma \chi_{12} - \chi_{12} \sigma, &&\quad \alpha_{12}^{(-)} = \sigma \hat{\alpha}_{12} - \hat{\alpha}_{12} \sigma.
\end{eqnarray}
where in contrast to the diagonal modes, here both the $(+)$ mode and $(-)$ mode are charged and couple to gauge field fluctuations in general. Dropping the indices and defining
\begin{equation}
A_\pm \equiv \frac{1}{2}(A_1 \mp A_2),\qquad F_\pm \equiv \frac{F_1}{X_1}\pm\frac{F_2}{X_2},\qquad h_\pm \equiv d \alpha_\pm
\end{equation}
we obtain the coupled equations of motion
\begin{eqnarray}
\Box \phi_\pm  - \left(4 A_\pm^2 - 2 (X_1+X_2)X_3 +\frac{1}{4}F_\pm^2\right)\phi_\pm && \nonumber\\
+\frac{X_1+X_2}{ X_1 X_2}\left((X_1-X_2) A_\pm\cdot \alpha_
\pm  + \frac{1}{4}F_\pm\cdot h_\pm \right) & = & 0, \\
\frac{1}{2\sqrt{-g}}   \partial_r  \left[  \frac{\sqrt{-g}}{X_1X_2} \left(2h^{rt}_\pm -(X_1+X_2)  F^{rt}_\pm \phi_\pm \right) \right] && \nonumber\\
 - \frac{(X_1-X_2)^2}{X_1 X_2} g^{tt} \alpha_\pm + \frac{(X_1+X_2)(X_1-X_2)}{ X_1 X_2} A^t_\pm \phi_\pm & = & 0.
\end{eqnarray}

When $q_1 = q_2$  we have that $\phi_{12}^{(-)}$ satisfies the same equation of motion as the charged mode $\phi_{11}^{(-)}$ \eqref{Weom}, and $\phi_{12}^{(+)}$ satisfies the same equation of motion as the neutral mode $\phi_{11}^{(+)}$ \eqref{plusdiag1}.

\subsection{Marginal modes}
Collecting all the modes above we have
\begin{center}
\begin{tabular}{c|c|c}
Type& Scalar field & See section\\
\hline
Charged diagonal & $\phi_{11}^{(-)},\phi_{22}^{(-)},\phi_{33}^{(-)}$  & \ref{sec:chargeddiag} \\
Off diagonal ($-$) & $\phi^{(-)}_{12},\phi^{(-)}_{13},\phi^{(-)}_{23}$  & \ref{sec:offdiagminus}\\
Neutral diagonal & $\phi_{11}^{(+)},\phi_{22}^{(+)},\phi_{33}^{(+)}$   & \ref{sec:plusmodes} \\
Off diagonal ($+$) & $\phi^{(+)}_{12},\phi^{(+)}_{13},\phi^{(+)}_{23}$ &  \ref{sec:plusmodes}\\
\end{tabular}
\end{center}
Note that the equation of motion for the $I\neq J$ mode $\phi_{IJ}^{(\pm)}$  reduces to the equation of motion for the $\phi_{II}^{(\pm)}$ when the background charges on the legs coincide: $q_I = q_J$.

We seek time-independent solutions that are regular at the horizon and normalisable at the AdS boundary for each of the fluctuations detailed in the table above. Fluctuations of the gauge field at the boundary are required to vanish so that the chemical potentials are not perturbed. Such marginal modes indicate new branches of superconducting and normal black brane solutions. Developing series solutions we find the near-horizon behaviour (here the horizon position has been fixed to $r_+=1$)
\begin{align}
\phi &= a_1 + \ldots  +  \log(r-1) (a_2 + \ldots ) \\
\alpha &=a_3 (r-1)+ \ldots
\end{align}
where the parameters $a_n$ are undetermined in the expansion and the ellipses denote determined higher order terms.  For regularity we set $a_2=0$, and  since the system of equations is  linear we are free to scale $a_1=1$. Additionally, for the diagonal modes we have fixed a shift symmetry in $\alpha$ such that $\alpha(r_+) = 0$. The scalar fluctuations in the UV vacuum all have mass $m^2 = -4$, which saturates the Breitenlohner-Freedman bound there.  The boundary series takes the form
\begin{align}
\phi &= r^{-2} (b_1 + \ldots ) +  r^{-2} \log r (b_2 + \ldots ) \\
\alpha &= b_3  + \ldots  +  r^{-2} (b_4 + \ldots )
\end{align}
with data $b_n$. Normalisability demands that we set $b_2=0$ and we must also impose $b_3=0$ to preserve our choice of chemical potentials.

\section{Sector $(-)$}\label{sec:minusmodes}
In this section we consider fluctuations of the fields labelled by $(-)$.
Before considering general ratios of the chemical potentials and general $p_1,p_2$ coordinates it is instructive to study a slice of the allowed region (figure \ref{fig:background}) at $\mu_1=\mu_2$, for which $p_2=0$. We find marginal modes with critical temperatures illustrated in figure \ref{fig:slice}.

On this slice we can consider the effect of an imbalance on the mode $\phi_{33}^{(-)}$, which is charged under $A_3$ and does not couple directly to any other gauge field. At $p_2=0$, $p_1$ becomes a function of the imbalance for this mode, $\zeta_{1,2} = \frac{\mu_{1,2}}{\mu_3}$. The imbalance $\zeta_{1,2} = 0$ at $p_1 = -\sqrt{2/3}$, $\zeta_{1,2} = 1$ at $p_1=0$ and $\zeta_{1,2} \to \infty$ as $p_1 \to \sqrt{1/3}$. As shown in figure \ref{fig:slice} the critical temperature as a function of this imbalance meets the temperature of thermodynamic instability at sufficiently high $\zeta_{1,2}$. Note however that before this value of $\zeta_{1,2}$ is reached we see superconducting branches for the modes $\phi_{11}^{(-)},\phi_{22}^{(-)}$ and $\phi_{12}^{(-)}$, which are charged under $A_{1,2}$.

Extending the analysis to the full allowed region of $(p_1,p_2)$ space, we find the modes with the highest critical temperature as labelled in figure \ref{fig:highestTc}.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.9\textwidth]{slice}
\setlength{\unitlength}{0.1\columnwidth}
\begin{picture}(0.1,0.25)(0,0)
\put(-1.9,-0.1){\makebox(0,0){$p_1$}}
\put(-9.3,4.0){\makebox(0,0){$\hat{T}$}}
\end{picture}
\caption{A $\mu_1=\mu_2$ slice ($\mu_3\neq 0$) as a warm-up to the full $(p_1,p_2)$ space of normal phase solutions. The light blue region is thermodynamically unstable. The lines indicate critical temperatures of \emph{green:} $\phi_{33}^{(-)}$, \emph{red:} $\phi_{13}^{(-)},\phi_{23}^{(-)}$, \emph{blue:} $\phi_{11}^{(-)},\phi_{22}^{(-)},\phi_{12}^{(-)}$. Finally we show the line of neutral scalar zero modes (see section \ref{sec:plusmodes}) in \emph{black:} $\phi_{11}^{(+)},\phi_{22}^{(+)},\phi_{12}^{(+)}$ . \label{fig:slice}}
\end{center}
\end{figure}

\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.6\textwidth]{highestTcRaster}
\setlength{\unitlength}{0.1\columnwidth}
\begin{picture}(0.1,0.25)(0,0)
\put(-3.,-0.2){\makebox(0,0){$p_1$}}
\put(-6.6,3.0){\makebox(0,0){$p_2$}}
\put(-2.5,4.2){\makebox(0,0){$\phi^{(-)}_{11}$}}
\put(-2.5,2.0){\makebox(0,0){$\phi^{(-)}_{22}$}}
\put(-4.5,3.1){\makebox(0,0){$\phi^{(-)}_{33}$}}
\put(-2.2,3.1){\makebox(0,0){$\phi^{(-)}_{12}$}}
\put(-3.5,3.8){\makebox(0,0){$\phi^{(-)}_{13}$}}
\put(-3.5,2.2){\makebox(0,0){$\phi^{(-)}_{23}$}}
\put(-0.9,3.05){\makebox(0,0){$(\mu,\mu,0)$}}
\put(-4,4.7){\makebox(0,0){$(\mu,0,\mu)$}}
\put(-4,1.3){\makebox(0,0){$(0,\mu,\mu)$}}
\put(-5.8,3.2){\makebox(0,0){$(0,0,\mu)$}}
\put(-1.8,5.2){\makebox(0,0){$(\mu,0,0)$}}
\put(-1.8,0.8){\makebox(0,0){$(0,\mu,0)$}}
\end{picture}
\caption{Modes which give the highest critical temperature. The diagonal modes $\phi_{11}^{(-)},\phi_{22}^{(-)},\phi_{33}^{(-)}$ give the highest critical temperature in the blue regions as shown. In addition, along the red lines the modes  $\phi_{12}^{(-)},\phi_{13}^{(-)},\phi_{23}^{(-)}$ linearly coincide with the diagonal modes and have the same critical temperature. Dots indicate special cases of the chemical potentials, $(\mu_1,\mu_2,\mu_3)$, with all equal at the origin.\label{fig:highestTc}}
\end{center}
\end{figure}

\subsection{$\phi^{(-)}_{11},\phi^{(-)}_{22},\phi^{(-)}_{33}$}\label{sec:chargeddiag}
The critical temperatures for $\phi^{(-)}_{33}$ are plotted in figure \ref{fig:chargeddiag}. It is important to emphasise that the region where  $\phi^{(-)}_{33}$ is unstable is not restricted to the segment indicated in figure \ref{fig:highestTc} --- there is a slight overlap of $\phi^{(-)}_{11},\phi^{(-)}_{22}$ and $\phi^{(-)}_{33}$ near lines of equal chemical potential.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.6\textwidth]{phi33Raster}
\setlength{\unitlength}{0.1\columnwidth}
\begin{picture}(0.1,0.25)(0,0)
\put(-3.,-0.2){\makebox(0,0){$p_1$}}
\put(-6.6,3.0){\makebox(0,0){$p_2$}}
\end{picture}
\caption{Critical temperatures for the mode $\phi^{(-)}_{33}$. The modes $\phi^{(-)}_{11},\phi^{(-)}_{22}$ can be obtained by permuting the labels and appear predominantly in the regions as labelled in figure \ref{fig:highestTc}. \label{fig:chargeddiag}}
\end{center}
\end{figure}

\subsection{$\phi^{(-)}_{12},\phi^{(-)}_{13},\phi^{(-)}_{23}$}\label{sec:offdiagminus}
We find that marginal modes exist in a window of allowed $p_1$ and $p_2$ values.  This region and their critical temperatures are plotted in figure \ref{fig:offdiagminus}. Note that the critical temperatures all lie below those found in section \ref{sec:chargeddiag}, except along the locus of equal chemical potentials on the two legs; for instance along $p_2=0$ the critical temperatures for $\phi^{(-)}_{12}$ and $\phi^{(-)}_{11}$ and $\phi^{(-)}_{22}$ coincide.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.9\textwidth]{offDiagMinusRaster}
\setlength{\unitlength}{0.1\columnwidth}
\begin{picture}(0.1,0.25)(0,0)
\put(-1.6,-0.05){\makebox(0,0){$p_1$}}
\put(-9.4,2.25){\makebox(0,0){$p_2$}}
\end{picture}
\caption{Region where marginal modes for $\phi^{(-)}_{12}$ exist, with critical temperatures indicated by the contours. Regions for the other modes $\phi^{(-)}_{23}$ and  $\phi^{(-)}_{13}$ can be found by permuting the indices. The blue shaded region is the region shown in figure \ref{fig:background}.  \label{fig:offdiagminus}}
\end{center}
\end{figure}

\section{Sector $(+)$ and Gubser-Mitra instabilities}\label{sec:plusmodes}
In this section we consider fluctuations of the fields labelled by $(+)$. The fields $\phi_{11}^{(+)},\phi_{22}^{(+)},\phi_{33}^{(+)}$ are neutral, while $\phi^{(+)}_{12},\phi^{(+)}_{13},\phi^{(+)}_{23}$ become neutral when two chemical potentials coincide. We find marginal modes along the loci
\begin{eqnarray}
\phi_{11}^{(+)},\phi_{22}^{(+)},\phi^{(+)}_{12} &\quad \mu_1=\mu_2 >\mu_3\nonumber\\
\phi_{11}^{(+)},\phi_{33}^{(+)},\phi^{(+)}_{13} &\quad \mu_1=\mu_3 >\mu_2\\
\phi_{22}^{(+)},\phi_{33}^{(+)},\phi^{(+)}_{23} &\quad \mu_2=\mu_3 >\mu_1\nonumber
\end{eqnarray}
where all fluctuations are neutral. The backreaction and unimodularity constraints for the diagonal modes can be satisfied by turning on two modes, e.g. $\phi_{22}^{(+)} = -\phi_{11}^{(+)}$.

It is perhaps surprising that unlike the $(-)$ modes of section \ref{sec:minusmodes}, these marginal modes exist only for one dimensional regions of the $(p_1,p_2)$ plane.  However, they do not indicate a true instability of the black brane background.  Instead, as we shall demonstrate, they have a natural interpretation as zero modes which move in the family of background black branes. The existence of these zero modes is closely related to the threshold of thermodynamic instability; indeed, the critical temperatures are given precisely by the lowest temperature at which the normal phase is locally thermodynamically stable.

We can analytically construct this class of marginal modes in the grand canonical ensemble by perturbing the charges and the horizon position whilst ensuring that the chemical potential and temperature remain fixed. A 4$\times$4 matrix for this variation is
\begin{equation}
\mathcal{M} = \frac{\delta \{\mu_I, T\}}{\delta \{q_I,r_+\}} \label{muThess}
\end{equation}
and, if the fluctuations are to be marginal modes in the grand canonical ensemble, we require $\det \mathcal{M} = 0$, which is achieved for background values of $q_I$ at the threshold of the thermodynamic instability given by $\det \mathcal{H}=0$. We emphasise that this applies to modes which include gravitational fluctuations at any values of the chemical potentials and is easily derived in this way. For example, when $q_3 = \frac{2-q_1-q_2}{1-q_1q_2}>0$ we have a surface of zero modes generated by perturbations to \eqref{gathered}, $H_I \to H_I + \delta H_I$, $f\to f+\delta f$, where
\begin{eqnarray}
\delta H_i &=& \frac{q_i}{(q_i-1)}\frac{1}{r^2},\qquad\qquad i=1,2\\
\delta H_3 &=& -\frac{q_1+q_2-2}{(q_1-1)(q_2-1)}\frac{1}{r^2}\\
\delta f &=& 2 + \frac{ q_1 q_2 (q_1+q_2-2)}{(q_1q_2-1)}\left(\frac{1}{r^2}-\frac{1}{r^4}\right)
\end{eqnarray}
Specialising to the non-gravitational modes of interest here which exist only along loci of the chemical potentials, we obtain\footnote{The no-backreaction condition imposes that any two of the background charges are equal; let us here choose $q_1 = q_2$.  Subsequent conditions on the fluctuations themselves become $\delta q_1+\delta q_2 = 0$ and $\delta q_3=0$. Further imposing that the chemical potential is not perturbed we have $q_2 = 1$. }
\begin{eqnarray}
\frac{\delta X_1}{X_1} = -\frac{\delta X_2}{X_2} &=& \frac{1}{1+r^2}, \qquad \delta X_3 = 0\\
\delta A_1 = -\delta A_2 &=&  \frac{\sqrt{(1+q_3)}(r^2-1)}{\left(1+r^2\right)^2}, \qquad \delta A_3 = 0
\end{eqnarray}
at temperature $\hat{T} = \frac{1}{\pi \sqrt{2(1+4q_3+q_3^2)}}$. This zero mode corresponds to a redistribution of the black brane charges when two of them are equal. These modes satisfy  equations \eqref{plusdiag1} and \eqref{plusdiag2} and consequently it is precisely this mode which gives rise to the black line in figure \ref{fig:slice}. It is also the zero mode generated by the thermodynamically unstable directions discussed in section \ref{sec:branches}.  We expect that there is also a true instability associated with these fluctuations \cite{Gubser:2000ec, Gubser:2000mm} but it is not captured by a zero mode analysis. We also anticipate a true instability away from the lines of equal chemical potentials, provided gravitational fluctuations are included.

\section{Discussion}\label{sec:discussion}
We have found new branches of superconducting black branes in $\mathcal{N}=8$ $SO(6)$ gauged supergravity, coming from a truncation involving 20 scalar fields and 15 gauge fields. Superconductivity persists to any values of the chemical potentials conjugate to the three R-charge densities of the background black brane.  The critical temperatures at which these black brane solutions arise is always higher than the temperature at which the three-charge background becomes thermodynamically unstable.

It is known that in some special cases --- such as two and three non-zero equal chemical potentials --- one has `retrograde condensation', \emph{i.e.}\ the branch of superconducting black branes is thermodynamically subdominant and the condensate exists only above a critical temperature \cite{Aprile:2011uq} (see also \cite{Buchel:2009ge, Donos:2011ut} in 4D).\footnote{For systems exhibiting retrograde condensation, the expectation is that superconducting ground states do not exist in the planar theory, but can exist when the boundary has spherical topology \cite{Gentle:2011kv}.} The other previously studied case of this system is that of one non-zero chemical potential, where the superconducting black brane is thermodynamically dominant \cite{Aprile:2011uq}. Here, in all known cases of retrograde condensation, we find that there is an additional marginal mode with the same critical temperature (see figure \ref{fig:highestTc}). We leave a nonlinear investigation of this branch, and whether it is dominant, to future work.

For a given mode we may interpret dialling the chemical potentials as introducing an imbalanced mixture as initiated in \cite{Erdmenger:2011hp, Bigazzi:2011ak}. The critical temperature surface cannot be extended to arbitrarily low temperatures as the imbalance is dialled, due the presence of a thermodynamic instability in the normal phase at finite temperature. At a given set of chemical potentials, different modes in the consistent truncation experience a different effective imbalance. Interestingly, we find that before a superconducting instability vanishes as a function of its imbalance, another appears to take its place.

The difference of the chemical potentials sets a scale for the system. It would be interesting to see if this scale sets a value of spatial momentum in this truncation resulting in an inhomogeneous, FFLO-like phase \cite{loff}.\footnote{For examples of spatially modulated phases in a truncation of IIB supergravity on $S^5$ see \cite{Donos:2011ff}.} It would also be interesting to construct all of the superconducting branches found nonlinearly, away from the critical temperature. This would allow the direct comparison of the free energy with the normal phase, and consequently the identification of the true phase boundary for the spontaneous breaking of a $U(1)$-bulk symmetry.

For many of the linear modes studied in this paper, labelled by $(+)$, we have found only neutral marginal modes, whose role is to move one infinitesimally within the family of background black branes. These modes occur precisely at the threshold temperature for thermodynamic instability. As such the marginal modes do not generate an actual instability of the black brane. We have shown that it is possible to populate the entire minimum temperature surface (as a function of chemical potentials) with zero modes of the background black brane, which can be found analytically and involve gravitational fluctuations in general. One expects to find a true instability associated with the thermodynamic instability in these cases.

Our work was motivated in part by the search for universal features in top-down holographic superconductivity. We wish to highlight a basic problem we have encountered when using consistent truncations: a subdominant solution in a given theory can appear as the dominant solution in a consistent truncation of that theory. This is exemplified by the normal phase solutions employed in this paper.

In the case of three equal chemical potentials the normal phase solution is simply RN at high temperatures. Whilst the RN branch exists at low temperatures too, it is no longer the dominant contribution to the ensemble. This has important consequences for the equal-$q_I$ consistent truncation where RN is the only normal phase solution remaining. Evidently in the equal-$q_I$ truncation RN appears as the dominant black brane branch at all temperatures, when in fact, at low temperatures in the full theory it is subdominant to a thermodynamically unstable black brane.

\subsection*{Acknowledgements}
It is a pleasure to thank Aristomenis Donos, Jerome Gauntlett, Mukund Rangamani and Julian Sonner for helpful discussions. SAG is supported by a STFC studentship. BW is supported by the Royal Commission for the Exhibition of 1851.


\section{Introduction}

Millisecond pulsars (MSPs) are neutron stars (NSs)
characterized by short
spin periods ($P_{\rm spin} < 20 $ ms) and weak surface
magnetic fields ($B < 10^{9}$ G), which are often found in
binaries with a white dwarf (WD) companion. It is generally
agreed that MSPs are old NSs,
recycled by accretion of mass and angular momentum from
the donor stars through Roche-lobe overflow (RLOF) during the previous
low-mass X-ray binary (LMXB) evolution. The NS was spun
up, due to mass accretion, into a MSP, while the donor evolved
to be a He or CO WD \citep[see e.g., ][for reviews]{bv91,tv06}.

The stability of mass transfer in an X-ray binary
depends on the ratio of the masses of the donor star and
the NS, and the initial orbital period of the system.
Traditionally, it was thought that, if the
mass ratio is large enough ($\gtrsim 1.5$), the mass
transfer is likely to be unstable, resulting in a common envelope (CE)
evolution \citep{p76,w84,il93}, where the NS
spirals into the envelope of the donor on a
very short timescale ($<10^{3}$ yr).
More recent investigations \citep[e.g.,][]{tvs00,p00,kdkr00,prp02,prp03}
show that  X-ray binary systems with
intermediate-mass (up to $\sim 5 M_{\odot}$) donor stars
(i.e., IMXBs) can avoid the spiral-in phase and experience
rapid mass transfer on a thermal timescale,
successfully evolve to become
LMXBs.

For LMXBs, a CE evolution is even unavoidable if the separation of
binary components is large enough, so that the donor reaches the asymptotic
giant branch (AGB) phase with a deep convective envelope
before RLOF (i.e. Case C RLOF).
Stable mass transfer in LMXBs usually occurs on a timescale
$\sim 10^{8}-10^{10}$ yr when the donor star is
on the main sequence or (sub)giant branch at the onset
of RLOF (i.e., Case A or B RLOF),
and forming MSPs seems to be feasible
in this process.
It was found by \citet{ps88,ps89} that there exists a critical
bifurcation orbital period ($P_{\rm bif}$),
the initial orbital period that
separates the formation
of converging LMXBs (which evolve with
decreasing orbital periods until the donor star
 becomes degenerate or an ultra-compact
binary is formed) from diverging LMXBs. The
value of $P_{\rm bif}$ is found to be $\sim 1$
d, but depends heavily on the processes of tidal
interactions and the mechanisms and efficiency of
orbital angular momentum loss \citep{e98,prp02,vs05,ml09}.

The final products of LMXB evolution are binary pulsars.
The distributions of the spin periods of the pulsars, the orbital
periods and the WD masses can be used to testify
the models of I/LMXB evolution.
\citet{d08} compared the theoretical expectations of
I/LMXB evolution to the populations of Galactic binary pulsars.
He showed that a significant population
of binary pulsars with 1 d $\lesssim P_{\rm orb} \lesssim 100$ d
are generally consistent with being the
descendants of long-period LMXBs or IMXBs.
However, there remain quite a few unresolved puzzles.
For example, binary pulsars with $P_{\rm orb}\gtrsim$ 60 d seem
to have WD companions less massive
than predicted by theory,
as pointed out previously by \citet{ts99}, and those
with 0.1 d $\lesssim P_{\rm orb}\lesssim$ 1 d are inconsistent
with any I/LMXB evolution.

In the previous studies on I/LMXB evolution, a canonical NS
(of mass $\sim 1.3-1.4\,M_{\odot}$)
was usually adopted. This seems to be supported
by the finding that the NS mass distribution
is consistent with a narrow Gaussian at
$1.35\pm 0.04\,M_{\odot}$ \citep{tc99}.
However, both
observations \citep[see][and references therein]{z11,k10,spr10,o12} and
theories \citep{n84,tww96,h03,wj05} suggest that the initial masses
of NSs may occupy a large range, probably
originating from two different mechanisms of
forming NS: iron-core collapse supernovae and
electron-capture supernovae.
NSs in  high-mass X-ray binaries (HMXBs)  have experienced
very little accretion because of their young ages, so their masses should be
very close to those at birth.
The measured masses of NSs in HMXBs range from $1.06^{+0.11}_{-0.10}
\,M_{\odot}$ for SMC X$-$1 \citep{v07} to
$1.86\pm 0.16\,M_{\odot}$ for Vela X$-$1 \citep{b01,q03}.
Recently \citet{r11} present an improved method for determining the mass of
NSs in eclipsing X-ray pulsar binaries and apply it
to six systems. They find that the NS masses  range from
$0.87\pm 0.07\,M_{\odot}$
(eccentric orbit) or $1.00\pm0.10\,M_{\odot}$ (circular orbit)
for 4U1538$-$52 to $1.77\pm 0.08\,M_{\odot}$ for Vela X$-$1.
So in a proper investigation, the influence of the NS masses
should be included. Recent evolutionary calculations by
\citet{db10} have provided evidence
that the evolution of I/LMXBs
depends upon the NS mass. In this work we perform
systematic calculations of I/LMXB evolution and discuss
the properties of the produced binary and millisecond pulsars (BMSPs), taking into
account different initial NS masses. We adopt the initial donor
masses to be $1.0 - 6.0 M_{\odot}$, and two
different initial masses (1.0 and 1.8 $M_{\odot}$) for the NS.

During the mass transfer processes, part of the
transferred mass from the donor star may escape
from the binary system, carrying away the orbital
angular momentum. The formation of BMSPs is closely
related to the mechanisms of mass and angular momentum
loss, and the stability of mass transfer is also
dependent on the angular mementum loss rate
\citep{spv97}.
Generally, it is assumed that all the mass transferred
from the donor will accrete onto the NS unless in a
super-Eddington mass transfer phase, during which the NS
will accrete at the Eddington rate
($\dot{M}_{\rm Edd}\sim 1.5\times 10^{-8}\,M_{\odot}$\,yr$^{-1}$
for a $1.4\,M_{\odot}$ NS), and the residual mass
escapes from the binary system, carrying the
NS's specific orbital angular momentum. This is
the so-called ``isotropic reemission model''.
Actually, there may also be mass loss even during the
sub-Eddington mass transfer phase. For example,
the accretion disk can become thermally
and viscously unstable
when the orbital period is larger than a critical value
\citep{v96,dlhc99}, leading to limit cycle behavior of the mass
transfer rate.
The NS can accrete mass only
during outbursts, while most matter may be ejected
out of the binary systems during quiescence by the
radiation and magnetic pressure of the rapidly rotating
NS  \citep{rst89,bdb02}.
In our calculations, mass loss due to super-Eddington
mass transfer, radio ejection caused by an unstable disk
\citep{bdb02} or bump-induced detachment
\citep{dvl06}, and in some cases outflow from the $L_{2}$
Lagrangian point  are included.

This paper is organized as follows. In Section 2
 we describe the stellar evolution code and the
binary model used in this paper. We present the calculated results
in Section 3 and discuss their possible applications
in the formation of BMSPs in Section 4,
and summarize in Section 5.

\section{Binary Evolutionary Calculations}

\subsection{The Stellar Evolution Code}

All calculations were carried out with an updated
version of the stellar evolution code developed
by \citet{e71,e72} \citep[see also][]{h94,p95}.
We set initial solar chemical compositions
(i.e., $X = 0.7$, $Y = 0.28$, and $Z = 0.02$) for the
donor star, and take the ratio of the mixing length to
the pressure scale height to be 2.0, and the general
convective overshooting parameter to be 0.12. In our
calculations, we have considered a number of
binary interactions in order to
follow the details of the mass transfer process,
including orbital angular momentum
loss due to gravitational radiation (GR),
magnetic braking (MB), mass loss, and the effect of disk
instability.

\subsection{ The Input Physics}

Each of the binaries initially consists of an NS
of mass $M_{1}$ ($=1.0$ or
1.8 $M_{\odot}$) and a zero-age main-sequence (ZAMS) donor of
mass $M_{2}$ ($\sim 1.0-6.0 M_{\odot}$). The efective
radius of the Roche lobe $R_{L}$ of the donor is given by the
following formula  \citep{e83},
\begin{equation}
R_{L}/a = \frac{0.49q^{-2/3}}{0.6q^{-2/3} + \ln(1 + q^{-1/3})},
\end{equation}
where $a$ is the orbital separation of the binary,
and $q = M_{2}/M_{1}$ is the mass ratio of the
binary components.
As usual we assume that tides keep the binary orbit
circular \citep{k88}, so that the orbital angular momentum is
\begin{equation}
J_{\rm orb} = \frac{M_{1}M_{2}}{M}\Omega a^{2},
\end{equation}
where $M =M_{1}+M_{2}$, and $\Omega=\sqrt{GM/a^{3}}$ is
the orbital angular velocity. Logarithmic
differentiation of Eq.~(2) with time gives the rate of
change in the orbital separation
\begin{equation}
\frac{\dot{a}}{a} = 2\frac{\dot{J}_{\rm orb}}{J_{\rm orb}}-
2\frac{\dot{M_{1}}}{M_{1}}-2\frac{\dot{M_{2}}}{M_{2}}+\frac{\dot{M}}{M},
\end{equation}
where the total changing rate in the orbital angular momentum is
\begin{equation}
\frac{\dot{J}_{\rm orb}}{J_{\rm orb}} =\frac{\dot{J}_{\rm GR}}{J_{\rm orb}}+
\frac{\dot{J}_{\rm MB}}{J_{\rm orb}}+
\frac{\dot{J}_{\rm ml}}{J_{\rm orb}}.
\end{equation}
The three terms on the right-hand-side of Eq.~(4) represent
angular momentum loss due to GR, MB,
and mass loss, respectively.
The rate of angular momentum loss due to GR
is calculated according to the standard
formula \citep{ll59,f71}
\begin{equation}
\frac{\dot{J}_{\rm GR}}{J_{\rm orb}}=
-\frac{32G^{3}}{5c^{5}}\frac{M_{1}M_{2}M}{a^{4}},
\end{equation}
where $G$ and $c$ are the gravitational constant and
the speed of light, respectively.
The prescription of \citet{vz81} is adopted to
calculate the angular momentum loss due to
MB,
\begin{equation}
\frac{\dot{J}_{\rm MB}}{J_{\rm orb}}=
-3.8\times10^{-30}\frac{GR_2^{4}M^{2}}{a^{5}M_{1}} \,{\rm s}^{-1},
\end{equation}
where $R_2$ is the radius of the donor.

For IMXBs and wide LMXBs, the mass transfer rate
$|\dot{M_{2}}|$ (note that $\dot{M}_2<0$) can be larger than the
Eddington accretion rate
$\dot{M}_{\rm Edd}$. Thus the
mass loss rate from the binary system is $\dot{M} =\dot{M_{2}}+
\dot{M}_{\rm Edd}$. Here, we adopt the isotropic reemission
model, assuming that the extra material leaves the binary in the form of
isotropic  wind from the NS, carrying  off
the NS's specific  orbital angular momentum
$j_{1}=(M_{2}/M_{1}M)J_{\rm orb}$.  The related angular momentum loss
rate can be derived to be
\begin{equation}
\frac{\dot{J}_{\rm ml,1}}{J_{\rm orb}}=
\dfrac{M_{2}}{M_{1}M}(\dot{M}_{2}+\dot{M}_{\rm Edd}).
\end{equation}

In some cases part of the material lost from the donor star
may escape the system through the $L_{2}$ Lagrangian point
rather accrete onto the NS.
Assuming that a fraction $\delta$ of the matter flow escapes
from the binary with the specific orbital angular momentum
$j_{2}=a_{L_{2}}\Omega$,
the angular momentum loss rate due to the $L_{2}$
point outflow is given by \citep{sl12}
\begin{equation}
\frac{\dot{J}_{\rm ml,2}}{J_{\rm orb}}=\delta\frac{a_{L_{2}}
^{2}}{a^{2}}\frac{M\dot{M}_{2}}{M_1M_2},
\end{equation}
where $a_{L_{2}}$ is the distance between the mass
center of binary and the $L_{2}$ point.

The transferred material from the donor will form
an accretion disk surrounding the NS. If the effective temperature
in the accretion disk is below $\sim6500$ K (the
hydrogen ionization temperature), the  accretion disk is likely to be
thermally and viscously unstable \citep{l01}. In LMXBs
irradiation from the
NS may help stabilize the disk to some extent
\citep{v96,k97,l01,r08}. The critical mass transfer
rate for the disk instability is given by \citep{dlhc99}
\begin{equation}
\dot{M}_{\rm cr}\simeq 3.2\times10^{-9}(\frac{M_{1}}{1.4M_{\odot}})^{0.5}(\frac{M_{2}}
{1.0M_{\odot}})^{-0.2}(\frac{P_{\rm orb}}{1.0\,{\rm d}})^{1.4} M_{\odot}\, {\rm yr}^{-1}.
\end{equation}
If the mass transfer rate is less than $\dot{M}_{\rm cr}$,
the X-ray binary is assumed to become transient,
experiencing short outbursts separated by
long quiescent intervals.
Material accumulates in the disk during the quiescent
phase while the NS
accretes mostly during outbursts. Define the duty
cycle $d$ \citep[$\sim 0.001-0.1$, see][]{k03}
to be the ratio of the outburst timescale to the recurrence
time, the accretion rate during outbursts can be estimated
as $\sim - \dot{M}_{2} /d$ . Of course, the mass
accretion rate in this case is also limited by the Eddington
accretion rate.

\section{ Results of I/LMXB Evolution Calculations}

\subsection{The Initial $M_{2}-P_{\rm orb}$ Parameter Space for
Successful Evolution}

In Fig.~1 we outline the results of our calculations
by showing the fate of the I/LMXB evolution in the initial
$M_{2}-P_{\rm orb}$ diagram, similar as in \citet{tvs00}.
The left and right panels correspond to $1.0\,M_{\odot}$ and
$1.8\,M_{\odot}$ NS, respectively. It is seen that the allowed
space for successful evolution into binary pulsars (i.e.,
without CE evolution) is larger
for $M_1=1.8\,M_{\odot}$ than for $M_1=1.0\,M_{\odot}$,
since the lower mass ratio in the former case can
stabilize mass transfer during the IMXB evolution.
When $M_1=1.0\,M_{\odot}$,
systems with a low mass companion ($M_{2}\lesssim 1.3 M_{\odot}$)
can form binary radio pulsars with a He/CO WD companion
if the initial orbital period is $<1000$ d,
while the lower limit of $M_{2}$ increases to
$\sim 2.2\,M_{\odot}$ when $M_1=1.8\,M_{\odot}$.
IMXBs can avoid spiral-in and CE evolution
when $1.5M_{\odot}\lesssim M_{2}\lesssim 3.3M_{\odot}$
and $2\, {\rm d}\lesssim  P_{\rm orb} \lesssim 10$ d for $M_1=1.0\,M_{\odot}$,
and $2.4M_{\odot} \lesssim M_{2}\lesssim 5.5M_{\odot}$ and
$2\,{\rm d}\lesssim P_{\rm orb} \lesssim 40$ d
for $M_1=1.8\,M_{\odot}$, respectively.
In both cases, if the initial $P_{\rm orb}$
is too short, the binary
systems will experience either a CE phase, or become X-ray binaries
with degenerate hydrogen stars \citep{e98}, possibly
forming ultra-compact binaries. On the other hand, if the initial
$P_{\rm orb}$ is too long, the donor will develop
a deep convective envelope at the onset of RLOF, and a runaway mass transfer
is initiated, leading to CE evolution. In \citet{tvs00}
systems with $M_1=1.3\,M_{\odot}$,
$2M_{\odot} \lesssim M_{2}\lesssim 5M_{\odot}$
and $1\,{\rm d}\lesssim P_{\rm orb} \lesssim 20$ d can survive
unstable mass transfer. These systems are right
 between the limiting cases
presented in Fig.~1.

To see in more detail how the NS mass affects the binary evolution,
we show in Fig.~2 the evolutionary paths of two IMXBs
with the same donor mass $M_{2} = 4\,M_{\odot}$ and orbital period
$P_{\rm orb}$ = 10 d, but different NS mass
(1.0 $M_{\odot}$ and 1.8 $M_{\odot}$).
At the age of $\sim 175.1$ Myr, the donor star starts to
transfer mass via RLOF. In the upper two panels with a 1.0
$M_{\odot}$ NS, mass transfer occurs rapidly, rising up to
$>10^{-4}M_{\odot}$\,yr$^{-1}$ within $\sim 10^4$ yr.
Meanwhile, the orbital period reduces to $\lesssim 4$ d.
This dynamically unstable mass transfer will result in
CE evolution and probably merging of the NS and the
He core of the donor.
In the lower two
panels with a 1.8 $M_{\odot}$ NS, mass transfer initially
proceeds at a rate of $\sim3\times10^{-5}M_{\odot}$\,yr$^{-1}$
for $\sim 10^{5}$ yr, during which most of the hydrogen-rich
envelope ($\sim 3
\,M_{\odot}$) is removed from the donor star.
The mass ratio inverts during this phase,
and subsequently the mass transfer rate decreases
to be several $10^{-7}M_{\odot}$\,yr$^{-1}$ for
$\sim 10^{6}$ yr.  During this time the NS  accretes about $0.012
M_{\odot}$ mass, and this is the primary process to spin up the
NS. The final product is a binary consisting of a $1.815
\,M_{\odot}$ (mildly) recycled NS and a $\sim 0.6M_{\odot}$
CO WD in an orbital period of $\sim36$ d.

\subsection{The $P_{\rm orb}^{final}-M_{\rm WD}$ Diagram}

Figure 3 shows the calculated correlation between the final orbital
period $P_{\rm orb}^{final}$ and the mass of the WD (the remnant of the
donor) $M_{\rm WD}$, for two different initial NS masses, $1.0
M_{\odot}$ (left panel) and $1.8 M_{\odot}$ (right panel).
As mentioned above, low-mass
donor stars ($M_{2}\la 1.3M_{\odot}$ or $\la 2.2 M_{\odot}$ for
$1.0 M_{\odot}$ or $1.8 M_{\odot}$ NS, respectively)
may evolve to be He WDs.
However, when the mass of He is accumulated to exceed
$\sim 0.4-0.5M_{\odot}$ in the core of the donor, He flash
will occur, giving rise to the formation of a CO core. The orbital period
will reach $\ga 10^{3}$ days.
Intermediate-mass donors can avoid
He flash due to their higher temperature, forming CO WDs
with the mass $M_{\rm WD}\gtrsim 0.33\,M_{\odot}$.
Their $P_{\rm orb}^{final}-M_{\rm WD}$ distribution
deviates from that of the low-mass branch obviously.

\subsection{Radio Ejection during the ``Bump-related" Detachment?}

When its hydrogen shell reaches a discontinuity
in the hydrogen content at the time of the first
dredge-up, a low-mass star will suffer
a temporary contraction, thus producing a
``bump'' in the luminosity function of the red giants.
\citet{dvl06} suggested that the orbital period gap
($\sim 20 -60$ days) of BMSPs may be
related to the bump-induced
detachment of the donor star from its RL.
If the NS has already accreted sufficient mass,
it may turn on as a MSP during the detachment phase.
Material transferred from the
donor, once it expands
again to refill its Roche lobe, may be inhibited by the pulsar's
radiation pressure and ejected at the inner lagrangian point $L_{1}$
\citep{rst89,bdb02},  so that no further mass accretion
would occur. Similarly, this so-called ``radio ejection" process
may also occur in LMXBs with an unstable accretion
disk \citep{bdb02}. In both cases the mass loss will influence
the orbital evolution of the binary and the spin evolution
of the NS. \citet{dvl06} considered the radio ejection only
during the bump-related detachment. Our calculations show that the disk
instability generally appears earlier than the bump phase.
Two examples are presented in Fig.~4. In the upper two panels, the
binary system initially consists of
a 1.0 $M_{\odot}$ NS and a 1.0 $M_{\odot}$
ZAMS donor star in an orbital period of 5.0 d. RLOF
initiates  at the age $\sim 1.186\times10^{10}$ yr, and the bump
occurs $\sim7\times10^{7}$ yr later, at which the donor mass has
decreased to be $\sim 0.52
M_{\odot}$. However, the disk becomes thermally unstable when the
donor mass is $\sim 0.9M_{\odot}$, much earlier than the
bump phase. If the NS's spin was accelerated during the prior mass
transfer, the accumulated mass in the disk during
quiescence will be ejected out of
the binary by the radiation pressure of the rapidly spinning
NS. We find a similar situation in the case of
a binary containing a 1.8
$M_{\odot}$ NS and a 1.5 $M_{\odot}$ donor star shown
in the lower two panels.
Compared with the effect of unstable disk accretion,
it seems that the bump-related detachment
might play a less important role in the evolution of LMXBs.

We finally note that the process of pulsar-driven mass ejection is
highly uncertain, since its condition and efficiency depend on several
unknown parameters \citep{fl11}, such as the value of the equilibrium period
that a NS will reach during the mass transfer. However, as
seen below, there is controversy on the estimate of the
equilibrium period for NSs in LMXBs, and one should be
cautious when considering the effect of radio ejection
on the binary evolution.

\subsection{The effect of outflow from the $L_{2}$ point}

In the binary evolution the stability of mass transfer
strongly depends on mass loss and related angular
momentum loss. Although we have taken into account
 various ways for mass loss, including super-Eddington
 mass transfer, unstable mass
transfer due to disk instability, and mass ejection due to
the propeller effect and the pulsar's radiation,
the detail processes are  complicated
and uncertain, and  simplified treatment
might not reveal the realistic situation.
In the following, we show an example of how outflows from the
binary systems through the $L_{2}$ point influence the evolution of
mass transfer, and change the conditions of forming BMSPs.

MSP J1614$-$2230 is a 3.15 ms pulsar of mass
1.97 $\pm$ 0.04 $M_{\odot}$ \citep{d10}. The mass of
its WD companion ($0.500 \pm 0.006M_{\odot}$) and the
orbital period (8.7 d) were also measured accurately.
\citet{l11} and \citet{tlk11} systematically investigated
the formation channels of this pulsar from an IMXB. They showed
that NS this massive are not easy to produce
in spite of the initially high mass of the donor star,
unless they were already born as a relatively massive NS.
However, \citet{tlk11} found that, for the system with a 1.8 $M_{\odot}$ NS,
the final orbital period will be always larger than 10
d, inconsistent with the observation of the pulsar.
This conclusion will not hold if we let a small fraction $\delta$
of the transferred mass leave the system from the $L_{2}$ point
\citep[e.g.,][]{b89}.
As an illustration,
in Fig.~5 we compare the evolution of an IMXB containing
a 1.8 $M_{\odot}$ NS and a 4.5 $M_{\odot}$ companion
star in an orbit of 2.4 d, with $\delta = 0$ (upper two panels)
and $\delta = 0.04$ (lower two panels), respectively.
In the former case the primary angular momentum
loss mechanism is the isotropic reemission around the NS
during the super-Eddington accretion phase.
The final  system consists of a 2 $M_{\odot}$ NS and 0.5 $M_{\odot}$
CO WD, with a orbital period 15 d. In the latter case
angular momentum loss due
to the $L_{2}$ point outflow plays a role when the system
evolves into an LMXB, and reduces the orbital period to $\sim 8.7$ days.

\section{Comparison with BMSPs}

\subsection{The $P_{\rm orb}-M_{\rm WD}$ relation for low-mass BMSPs}

In wide LMXBs, the donor will climb to
the red-giant branch (RGB) in the HR diagram before RLOF.
For low-mass stars ($< 2.3 M_{\odot}$) on the RGB,
there is a well known relationship between the mass of the degenerate
He core and the radius of the giant star, which is almost
entirely independent of the mass of the hydrogen-rich envelope
\citep{rw71,w83}. Based on this relationship, a specific
correlation between the orbital period $P_{\rm orb}$
and the WD mass $M_{\rm WD}$ is obtained
\citep{j87,r95,ts99}. Comparison with the observations
shows that a significant population of BMSPs with He WD
companion is generally
consistent with this $P_{\rm orb}-M_{\rm WD}$ relation.
However, there seems to be a systematic deviation
from the correlation for
pulsars with $P_{\rm orb} \gtrsim 60$ d, which
seem to have WD companions lighter than expected
\citep{t96}.

Both systematic small values of the orbital inclination
$i$ and large NS mass can increase $M_{\rm WD}$
for the given observed mass functions.
Since there does not seem to be any observational selection
effect favoring small inclination angle $i$ \citep{t96}, we first examine
whether the $P_{\rm orb}-M_{\rm WD}$ correlation can be accounted
for if the long-period BMSPs have CO WD
companions. For example,
\citet{s05} noticed that the \citet{ts99} $P_{\rm orb}-M_{\rm WD}$
relation is incompatible at the $99.5\%$ level  with
a uniform distribution of $\cos i$ if the pulsar masses are drawn
from a Gaussian distribution centered on $1.35 M_{\odot}$
with width $0.04 M_{\odot}$, and better agreement with
uniformity in $\cos i$ can be reached if the pulsar masses are
large on average (e.g., $1.75\pm 0.04M_{\odot}$).
An extreme example is PSR B0820$+$02, which has
a $0.6\,M_{\odot}$ CO WD
companion in a very wide orbit with $P_{\rm orb}\simeq 1232$ d
\citep{kr00}. However, in such  wide binaries
it is difficult for the NS to accrete  enough matter
\citep{lw98,ts99}, and the NS must be born
heavy if this interpretation is correct. In Fig.~6 we compare the
relations between
$P_{\rm orb}$ and $M_{\rm WD}$
in the cases that the NS has an initial mass of 1.0 $M_{\odot}$
(left panel) and 1.8 $M_{\odot}$ (right panel).
Also plotted are binary pulsars with measured $P_{\rm orb}$
and $M_{\rm WD}$ (90\% probability
mass range for randomly oriented orbits) for a fixed NS mass
of 1.2 and
2.0 $M_{\odot}$, respectively (i.e., in each case the NS is
assumed to have accreted  $0.2 M_{\odot}$
during the mass transfer).
It is seen that some binary pulsars with $P_{\rm orb} >$ 100 d can fairly
match the relation if they have massive NSs
($\sim2M_{\odot}$) and heavy
CO WDs (although in some cases small orbital inclination angles
may be required),
while for those with $P_{\rm orb} <$ 20 d, statistically lighter NSs
($\sim1.2M_{\odot}$)
seem to follow the relation better.

It is expected that the NS has to accrete at least  a few $\sim
0.01\,M_{\odot}$ mass to reach a millisecond period, but
this is difficult to achieve for NSs in wide binaries \citep[see also][]{lc11},
because the mass transfer rate (which increases with increasing
orbital period) is likely to be super-Eddington,
and the accretion disk is likely to be unstable.
In Fig.~7 we show the mass transfer rate
$|\dot{M}_{2}|$, and the accreted mass $\Delta M_{1}$
of the NS as a function of the final
orbital period $P_{\rm orb}$.
Note that here we plot the mass
transfer rate only for the stable mass transfer phase -
when the accretion disk becomes unstable,
we take it to be the critical value $\dot{M}_{\rm cr}$
at the onset of disk instability.
Hence both the mass accretion
rate and $\Delta M_{1}$ are limited by the value
of $\dot{M}_{\rm cr}$: when
$|\dot{M}_{2}|$ is less than the $\dot{M}_{\rm cr}$,
unstable disk accretion occurs, and the
NS is assumed to accrete mass only during outbursts,
and part of the matter will be ejected out of the binary
systems if $|\dot{M}_{2}/d|$ is super-Eddington.

The left panel of  Fig.~7 shows the results for a
$1.0\,M_{\odot}$ NS
with a companion star of initial mass $M_2=1.0$, 2.0 and 3.0
$M_{\odot}$, respectively. Generally more massive
donor stars result in higher mass transfer rate,
which also increases with longer orbital
period $P_{\rm orb}$.
It is noted that (1) $|\dot{M}_{2}|$
is always $\gtrsim 0.1\dot{M}_{\rm Edd}$,
and (2) $\Delta M_{1}\lesssim 0.3 M_{\odot}$,
which decreases with $P_{\rm orb}$, because in wider
systems the mass accretion rates (with both  stable and
unstable accretion disks) are more likely to be super-Eddington.
The right panel is for the systems with a $1.8 M_{\odot}$
NS. In the case of  light companion star
($\sim 1.0 M_{\odot}$), the initial orbital period should be larger
than $\sim1$ d, so that the critical mass transfer
rate $\dot{M}_{\rm cr}$ is always larger than the
mass transfer rate, since $\dot{M}_{\rm cr}$
increases during the mass transfer process
(with increasing orbit period),
while the mass transfer rate $|\dot{M}_{2}|$ decreases
all the way \citep{w83}.
The accreted mass $\Delta M_{1}$ by the NS is thus
significantly lower than
those with more massive donors.
So only results with  2.0 and 3.0
$M_{\odot}$ donor stars are presented,
in which $\Delta M_{1}\lesssim
0.35 M_{\odot}$, and $\dot{M}_{2}\gtrsim
0.3 \dot{M}_{\rm Edd}$.

Figure~7 shows that in general
$\dot{M}_{2}\gtrsim 0.1 \dot{M}_{\rm Edd}$,
and $\Delta M_{1} \lesssim 0.3 M_{\odot}$.
The  distribution of the mass transfer rate is considerably
higher than that ($\dot{M}_{2} \sim10^{-11}-10^{-8} M_{\odot}$\, yr$^{-1}$)
obtained by \citet{prp02}, but more compatible with
the observations of persistent LMXBs. The main reason is that
we have taken into account
the effect of disk instability.  The $\Delta M_{1}$ distribution
is roughly in line with \citet{lc11},
who found that $\Delta M_{1}$ is generally
less than 0.6 $M_{\odot}$ in their calculations for
systems  with a
1.4 $M_{\odot}$ NS and a $1.0-2.0 M_{\odot}$
donor star. Moreover, we find that systems with initially massive
NSs ($1.8\,M_{\odot}$)
may accrete enough mass to evolve into BMSPs
 with 60 d $\la P_{\rm orb}\la 200 $ d, , while light
NSs in wide binaries are more likely to be partially
recycled. This seems to
support our conjecture that some of the wide BMSPs
might be born massive. Hopefully accurate measurements of both the
pulsar and the WD masses will help settle this issue.

\citet{l05} performed timing observations of the
BMSP J1640$+$2224 (with an orbital period of
175 d), and constrained the WD mass to be
$0.15^{+0.08}_{-0.05}\,M_{\odot}$
($1 \sigma$ uncertainties), which indicates that
the companion is very likely to be a low-mass He WD.
In this case the massive NS + CO WD model obviously
does not work, and one has to explore other possibilities.
Evaporation of the companion star from a
wind of relativistic particles after the pulsar turns on may
decrease the companion mass significantly, but it is unlikely
for wide BMSPs, since the evaporation timescale
would be longer than the Hubble time
\citep{t96}.

Here we suggest another possible solution.
It is interesting to note that solar-type stars  are
usually found to be surrounded by sub-stellar companions
(usually planets and/or brown dwarfs) \citep{c12}.
One may expect  that in some relatively wide
LMXBs the companion star had
possessed substellar companion(s) in close orbits
like ``hot Jupiters".
When the star evolved on the giant branch it would become
big enough to capture its planet/brown dwarf.
The planet/brown dwarf spiraled into the envelope of the
giant to initiate a CE phase. The
frictional drag arising from its motion through the
CE would lead to loss of its orbital angular
momentum and deposit of orbital energy in the envelope.
If there was enough orbital energy, the spiral-in process would
expel the envelope of the giant, leaving a WD remnant
\citep{nt98}\footnote{In addition, \citet{bs10} suggest that the
binary systems may reach
stable synchronized orbits before the
onset of the CE phase. Such stable synchronized orbits allow the
RGB star to lose mass prior to the onset of
the CE phase. Even after the secondary enters the giant envelope,
the rotational velocity is high enough to
cause an enhanced mass-loss rate.}.
If the initial separation between the star and the substellar
object(s) is less than tens of Solar radii,
the final outcome would be an under-massive WD
with or without the surrounding planet/brown dwarf, depending on
whether it evaporated, filled its own Roche-lobe, or survived.
The discovery of a $0.053 M_{\odot}$ brown dwarf
in a short (0.08 d) period orbit
around a $0.39 M_{\odot}$ WD 0137$-$349 \citep{m06} presents
strong observational evidence for this interaction.

During the planet/brown dwarf-involved CE phase, the companion's
envelope expanded rapidly and filled its Roche-lobe,
leading to mass transfer onto the NS. Since this phase
was very short ($<10^3$ yr), the mass transfer rate would be
much higher than the Eddington limit rate for the NS,
so that the NS accreted very small mass $\sim 10^{-5} M_{\odot}$
(unless the CE phase lasted much longer time).
If this is the case, we have to require that
{\em some MSPs were born this way, rather recycled during the
LMXB evolution}. In the literature, this idea has already been discussed
by \citet{mh01}, who showed that, the existence
of the innermost, moon-sized planet in the PSR 1257$+$12 system
suggests that the pulsar was born with
approximately its current spin frequency and magnetic field.
A schematic view of the formation of BMSPs
with planet/brown dwarf-involved CE evolution is shown in Fig. 8.

Not only having an impact on the formation of MSPs,
the substellar objects,
if really exist in low-mass binaries, may also
play a role in the CE evolution, especially influence the
estimate of the CE efficiency parameter $\alpha_{\rm CE}$.
Recently \citet{dkk12} reconstructed
the CE phase for the current sample
of post-CE binaries (PCEBs) with observationally
determined component masses and orbital periods.
Searching for correlations between $\alpha_{\rm CE}$ and the
binary parameters, they found that,  when the
internal energy of the progenitor primary envelope is taken into
account,  $\alpha_{\rm CE}$
decreases with increasing mass $M_{\rm p}$ of the primary
(i.e., the progenitor of the WD) \citep[see however][]{dm11}, and
$\alpha_{\rm CE}\gtrsim 1$ for $1\lesssim M_{\rm p}/M_{\odot} \lesssim 2$,
which seems to be in contrast with $\alpha_{\rm CE}\sim 0.25$
derived by \citet{rt12} from numerical simulations.
If there are planets/brown dwarfs around these
low-mass primaries, they can contribute extra orbital
energy to
help expel the primary's envelope during the RGB/AGB
phase, making the CE efficiency parameter within the canonical
range $0 < \alpha_{\rm CE}< 1$.

\subsection{The Rebirth Periods of BMSPs}

Accretion onto the NS in I/LMXBs changes both the mass and spin
of the NS. The spin evolution of an accreting NS depends on
the interaction between the magnetosphere and the accretion
disk, as described as follows \citep[e.g.,][]{bv91}: the accretion disk is truncated
at the magnetospheric radius $R_{\rm m}$, which is
close to the Alfv\'en radius of the disk
\[
R_{\rm m}\simeq R_{\rm A}=(\frac{B_{\rm s}^{2}R_{\rm s}^{6}}
{\dot{M}_{1}\sqrt{2GM_{1}}})^{2/7},
\]
 where $B_{\rm s}$ and $R_{\rm s}$ are
the surface magnetic field and the radius of the NS,
respectively. Material is channeled to the NS at the inner radius,
producing a spin-up torque. If the NS is spinning very rapidly,
the disk may be truncated outside the corotation radius
$R_{\rm co}=(GM_1P^{2}/4\pi^{2})^{1/3}$, and the NS experiences
a centrifugal barrier that can inhibit
accretion \citep[i.e., the so-called ``propeller" effect;][]{is75}.
So the NS will eventually reach the
equilibrium spin period such
that $ R_{\rm m}=R_{\rm co}$, or
\begin{equation}
P_{\rm eq} \simeq 2.4 B_{9}^{6/7}R_{6}^{18/7}M_{1}^{-5/7}
(\dot{M}_{1}/\dot{M}_{\rm Edd})^{-3/7}\,{\rm ms},
\end{equation}
where $B_{9}=B_{\rm s}/10^{9}$ G, and
$R_{6}=R_{\rm s}/10^{6}$ cm. This period can be regarded
as the beginning or rebirth period of recycled
pulsars when accretion terminates.

The rebirth period of a MSP can be
derived only when its actual spin-down
time is known. \citet{bk11} recently reported the optical discovery of the
companion to the $\sim 2 M_{\odot}$  MSP J1614$-$2230.
The optical colors show that the $0.5 M_{\odot}$
companion is a 2.2 Gyr old CO WD.
From the age of the WD, \citet{bk11} calculated the period of the
pulsar at birth, and found that pulsar should be born with a spin close
to its current value, implying
that the final accretion rate was $<10^{-2}\dot{M}_{\rm Edd}$.
This value is two orders of magnitude smaller than the estimate from
theoretical calculations by \citet{l11} and \citet{tlk11}.
These authors
suggested that the system
began as an IMXB consisting of a NS and a $\sim 4 M_{\odot}$
main-sequence secondary, which evolved to
be a CO WD with He envelope. The NS gained the most
mass during the final LMXB phase lasting $\sim 5-10$ Myr
at near-Eddington rates.

We present a systematic view on the relation
between $P_{\rm eq}/B_{9}^{6/7}$
and $P_{\rm orb}$ derived from theory and observations in Fig.~9.
The solid lines outline the theoretically
expected distribution of $P_{\rm eq}/B_{9}^{6/7}$
($\varpropto\dot{M}_{1}^{-3/7}M_{1}^{-5/7}$,  here $R_{6}$ is
taken to be 1) from binary evolution calculations, and
the symbols represent the observations of recycled pulsars
\citep[data are taken from the ATNF Pulsar Catalogue;][]{m05}.
Comparison between observations and theory shows
that, except in a few cases,
the observational $P_{\rm spin}/B_{9}^{6/7}$ is generally
larger than $P_{\rm eq}/B_{9}^{6/7}$ by $\sim 1-2$ orders of
magnitude, challenging the simple recycling theory\footnote{Here
we assume that the current spin periods of MSPs
are not far from their initial ones. This holds if their characteristic
ages are longer than the cooling
ages of the WDs \citep[e.g.][]{l95} or the Hubble time.}.
Thus, some other mechanisms must work to
produce efficient spin down torque(s) or reduce the spin-up
torque due to mass accretion, if the stellar evolution models
are correct. The proposed explanations can be
summarized as follows.

(1) There is no definite spin equilibrium since GR
may remove the angular momentum from
the NS. This idea was first suggested by \citet{pp78} to account
for the cutoff in the spin distribution of NSs in LMXBs.
The main emission mechanisms involve crustal mountains
\citep{b98}, magnetic deformations \citep{c02}, and
unstable $r$-modes \citep{a98} in the NSs. All these processes can produce
a substantial mass quadrupole moment and thus a
spin-down torque due to GR. The problem with
this interpretation is that,
recent observational results on some millisecond
X-ray pulsars have shown that the efficiency of GR
induced spin-down might be too low to be responsible
for balancing the spin-up process during outbursts
\citep[see][and references therein]{hp11,phd12}. So the
following models focus on modifications of the equilibrium
period (Eq.~[10]) that the NS will finally reach.

(2) If the NS magnetic field lines can thread the
accretion disk, an extra magnetic torque will be
exerted on the NS \citep{gl79a}, and the equilibrium period will
be $\sim 1-3$ times that in Eq.~(10) \citep{gl79b,w95,lw96}.
\citet{a05} further suggested that the inner disk region
may be geometrically thick
and sub-Keplerian, and dominated by radiation pressure,
if the mass transfer
rate is above a few percent Eddington accretion rate.
The coupling between the disk and magnetic field can reduce
the amount of angular momentum deposited
onto the NS from accretion by a factor $A$ (the
ratio of orbital angular velocity to Keplerian rotational
speed in the disk, $A<1$), compared with the
case of thin disks ($A = 1$), thus a decreased spin-up
torque results.
\citet{phd12} show that the existence of spin equilibrium
as set by the disk$-$magnetosphere interaction is able to
explain the observations of millisecond X-ray pulsars, if
the spin-down torque coming from the interaction between
the disk and field outside corotation is sufficiently large.
However, the radial extent of the coupling between the
star and the disk is highly controversial \citep[e.g.][]{mp05,g07}.
Stellar and disk winds may also take away the angular
momentum of the NS \citep{g97,g99,u06,r09,zl10}. In
addition, the abrupt torque reversals observed in
disk-fed X-ray pulsars Her X-1, Cen X-3, and 4U 1626$-$67
\citep{b97} suggest that current models of
disk$-$magnetosphere interaction
may have severe limitations.

(3) If Eq.~(10) does apply, a decaying mass transfer rate
at the end of the mass transfer process will result in relatively
large rebirth period of BMSPs, provided that the $\dot{M}_2$ evolution
time is longer than the NS spin evolution timescale \citep{j86,rst89}.
More recently \citet{t12} suggested that
this equilibrium will even be broken when
$|\dot{M}_2|$ decreases substantially,
so that the magnetospheric radius $R_{\rm m}$ becomes
larger than the corotation radius $R_{\rm co}$, causing the NS
to enter the propeller phase.
A centrifugal barrier arises to expel matter entering
the magnetosphere with a braking torque
\begin{equation}
N\simeq \dot{M}_2R_{\rm m}^{2}\Omega_{\rm K}(R_{\rm m})
\end{equation}
(where $\Omega_{\rm K}(R_{\rm m})$ is the Keplerian
angular velocity at $R_{\rm m}$) to act to slow down the
NS.
An alternatively
appropriate form for the accretion torque is \citep[e.g.][]{m99}
\begin{equation}
N\simeq -\dot{M}_2R_{\rm m}^{2}[\Omega_{\rm K}(R_{\rm m})
-\Omega_{\rm s}],
\end{equation}
where $\Omega_{\rm s}$ is the angular velocity of the NS.
With Eqs.~(11) and (12) we recalculate the spin evolution of the NS when the
donor star gradually decouples from its RL as in
\citet{t12}. In Fig.~10 the blue and red
lines denote the results calculated with Eqs.~(11) and (12),
respectively\footnote{The discontinuity at the point where
the two lines meets is due to the fact that we use
Eq.~(12) to calculate the spin evolution even before the
propeller phase in both cases.}. It is seen that the decaying mass transfer
rate indeed causes an increase in the spin period, which
is larger than the equilibrium period corresponding to
the average mass transfer rate before decoupling. However,
The calculated discrepancy seems not to be large enough to explain the
deviation of the rebirth periods of BMSPs from that given by
Eq.~(10), especially if we adopt the torque form of Eq.~(12).

(4) A decaying magnetic field may also change the rebirth
period of the BMSPs, if the field decay
timescale is shorter than the spin evolution timescale,
so that the final spin period deviates from the equilibrium
period. During the decay of the NS magnetic field induced by accretion,
the spin-up timescale
increases as $R_{\rm m}$ decreases because the accreting
matter carries less specific angular momentum, while the
field decay timescale decreases  as $R_{\rm m}$ decreases
(field decay stops when $R_{\rm m}=R_{\rm s}$), resulting in a
significant departure from the spin-up line \citep{b96,kb99}.

(5) The rebirth periods of BMSPs may be related to the long-term
spin equilibrium with varying mass accretion rate. A
significant fraction of LMXBs are likely to be transient systems subject to
disk instability \citep{cfd12}. For example, an NS LMXB in a 10 hr
orbit would be transient if the average mass transfer rate is lower
than $\sim 10^{17}$ gs$^{-1}$.
In addition, in narrow LMXBs, X-ray irradiation of the
donor star could destabilize the mass transfer, and lead to mass transfer
cycles \citep{h93} - mass transfer is spasmodic with phases
of high mass transfer driven by the thermal expansion of the
convective envelope of the irradiated donor alternating with
phases with low or no mass transfer, during  which the donor
readjusts towards thermal equilibrium of the un-irradiated star.
The final spin period is thus determined by the
spin-up during outbursts/high mass transfer phase
balanced by the spin-down during
quiescence/low mass transfer phase, which may be considerably
longer than that
attained via stable mass accretion \citep{l98,phd12}.

Currently our knowledge about the evolution of LMXBs
and accretion disks is not sufficient to tell which one is the
dominant factor in determining the rebirth periods of BMSPs.
Perhaps they are influenced by a combination of
(at least some of) the above-mentioned mechanisms.

\subsection{Origin of Intermediate-Mass BMSPs with $P_{\rm orb} \lesssim$ 1 d}

It is well known that a large fraction of IMXBs will produce
binary pulsars with a CO WD companion.
As shown in \textbf{Fig.~3}, the orbital periods of these intermediate-mass
binary pulsars (IMBPs) are generally $>\sim 3-6$ d. This leaves the puzzle
that the binary pulsars with $P_{\rm orb} \lesssim$ 1 d and
a CO WD companion
are hard to explain reasonably in any binary population \citep{d08}.
In Table 1 we list the parameters of IMBPs
with $P_{\rm orb} \lesssim$ 1 d. Except
PSR B0655$+$64, other pulsars share the characteristics of
rapid spin and low magnetic field.

The small orbital periods combined with massive WDs
suggest that this type of systems may have evolved
through a CE and spiral-in phase \citep{vt84}.
Similar to that of the Double NSs, their direct progenitors could be
binaries consisting of a He star and a NS star in a close and circular orbit,
as a result of spiral-in of a wide binary in late Case B or Case C
mass transfer. The core of the companion was a He star which could
already have gone through quite some He burning. After the spiral-in one
then had a He star with already CO in its core.
During He-shell burning, the envelopes of the He stars, if their masses are
$ < 3.5\,M_{\sun}$, will expand and slowly transfer mass to the NS,
which will have caused the spin-up and recycling of the NS.
This interpretation seem to be
responsible for PSR B0655$+$64 \citep{vt84}, which has the longest
spin period (195.6 ms) and highest magnetic field
($B\sim 10^{10}$ G) among the pulsars in Table 1,
implying that it has been partially
recycled with little mass accreted.
For other pulsars, the
short periods and low magnetic fields ($B\sim 10^{8}$ G)
require that at least some $10^{-2}\, M_{\odot}$ material has been
accreted by the NSs. Recently
\citet{clx11} presented the detailed binary evolution calculation
for a binary consisting of a NS (of mass 1.3 $M_{\odot}$) and a
low-mass He star
(of mass 1.0 $M_{\odot}$) with an initial orbital period of 0.5 d.
They showed that the mass transfer seems to be able to
spin up the NS's spin to milliseconds, producing BMSPs like PSR J1802$-$2124,
which has a short orbital period $P_{\rm orb} \lesssim$ 1 d.

Another possible way to form compact IMBPs
invokes evolution of IMXBs with anomalous MB.
Traditionally MB is thought to work in low-mass,
main-sequence stars.
However, There are many intermediate-mass stars
which have anomalously high magnetic fields ($\gtrsim$ 100 G), i.e.,
Ap/Bp stars. In IMXBs, the irradiation-driven stellar wind from the
surface of the Ap/Bp donor star by the NSs  may couple with the
high magnetic fields, resulting in angular momentum
loss at a rate \citep{jrp06},
\begin{equation}
\dot{J}_{\rm MB} = -B_2(\dfrac{\psi\dot{M}_{2}M_1}{a^{3}})^{1/2}
(\frac{R_2^{15}}{GM_{2}^{3}})^{1/4},
\end{equation}
where $B_2$ is the magnetic field strength at the surface
of the donor star, and $\psi$ is a parameter with typical value
$\psi/c^{2}\sim 10^{-6}$,
combining the wind-driving energy conversion efficiency and
irradiation geometry. \citet{jrp06} explored the effect
of anomalous MB in the formation
of compact BH LMXBs, in which the orbital energy
of the companion star, if of low-mass initially, would be
unable to expel the envelope of the BH progenitor.
For the black hole LMXB XTE J1118$+$480, \citet{g12} measured a
decay rate of the orbital period
$\dot{P}_{\rm orb}=-1.83\pm0.66 $ ms\,yr$^{-1}$.
This value is much larger than predicted by
most conventional MB and mass loss models
for the binary parameters of XTE J1118$+$480, unless
the companion star has a surface magnetic field
$B_2 \geq 10-20$ kG to enhance MB.

In Fig.~11 we show the evolution of a binary initially
consisting of a 1.8 $M_{\odot}$ NS and a 3.5 $M_{\odot}$ donor
with $P_{\rm orb}\simeq 2.5$ d. We assume that the
donor possesses a surface magnetic field $B_2= 500$ G.
RLOF starts at the age $\sim2.4\times10^{8}$ yr
with a mass transfer rate $\sim10^{-5} M_{\odot}$\,yr$^{-1}$.
Because the donor is initially more massive than the NS,
the orbital period decreases rapidly. Most of the transferred
mass is assumed to be ejected out of the binary
in the vicinity of the NS, taking away both mass and orbital
angular momentum. The orbit starts to widen when
$M_2\sim 2\,M_{\odot}$.
When the donor mass decreases to be $\sim 0.8 M_{\odot}$,
the mass transfer rate reduces to be
$< 10^{-7}M_{\odot}$\,yr$^{-1}$. In the subsequent
$\sim8\times10^{7}$ yr, MB begins to dominate the evolution.
The orbital period keeps decreasing,  the NS accretes
most mass in this phase and is spun up to milliseconds.
The final outcome is a tight binary system
consisting  of a $\sim 2.1\,M_{\odot}$ NS and
a $\sim 0.5 M_{\odot}$ WD with $P_{\rm orb}\sim$ 0.4 d.

The anomalous MB model will not work if the
NS is initially of low-mass ($\sim 1\,M_{\odot}$), since the
allowed parameter space for forming binary pulsars will be
very small (see Fig.~1).
Additionally we note that models involving
accretion-induced collapse of an ONeMg WD
cannot account for the $P_{\rm orb} $ distribution
of these BMSPs \citep{sl00}.

\section{Summary}

This work is motivated by the fact the evolution of I/LMXBs with
canonical NSs seems to meet difficulties in explaining some of the
observational characteristics of BMSPs, and the
measurements of the NS masses indicate a
wide distribution $\sim 1-1.8\,M_{\odot}$.
We have preformed numerical
calculations of the evolution of I/LMXBs
consisting of a 1.0 or $1.8 M_{\odot}$ NS and a
$1.0-6.0 M_{\odot}$ donor star, to investigate its dependence
on the initial NS mass, and the properties of the descendent
BMSPs.
The main results can be summarized as follows.

1. The allowed parameter space in the initial $P_{\rm orb}
-M_2$ diagram for forming recycled pulsars
increases with increasing NS mass. This
may help explain the formation of BMSPs with
$P_{\rm orb}\gtrsim 60$ d  and their distribution in the
 $P_{\rm orb}-M_{\rm WD}$ diagram.
Alternatively, some of these wide binary pulsars may be formed
through mass transfer driven by planet/brown dwarf-involved
CE evolution.

2. The equilibrium spin periods $P_{\rm eq}$ of accreting
NSs in LMXBs derived from the standard
magnetosphere-accretion disk interaction model are in general
shorter than the observed spin periods  of BMSPs
by more than one order of magnitude. This implies that either the simple
equilibrium spin model does not apply for the spin evolution
in accreting NSs, or there are
other mechanisms/processes to spin down the NSs when forming
BMSPs.

3. Some of the compact IMBPs
might have evolved from IMXBs in which the companion
star were strongly magnetized Ap/Bp stars with enhanced MB.

4. Our calculations doubt the suggestion that the orbital
period gap ($\sim20 -60$ d) of BMSPs is related to the
occurrence of the bump-related detached phase in the
LMXB evolution, since  the accretion disk would become
unstable at earlier time.

\begin{acknowledgements}

We are grateful to Ed van den Heuvel for helpful discussions and
the referee for constructive comments.
This work was supported by the Natural Science Foundation of China
under grant number 11133001, the
National Basic Research Program of China (973 Program 2009CB824800),
and the Qinglan project of Jiangsu Province.

\end{acknowledgements}