Pairing induced superconductivity in holography

We study pairing induced superconductivity in large $N$ strongly coupled systems at finite density using holography. In the weakly coupled dual gravitational theory the mechanism is conventional BCS theory. An IR hard wall cut-off is included to ensure that we can controllably address the dynamics of a single confined Fermi surface. We address in detail the interplay between the scalar order parameter field and fermion pairing. Adding an explicitly dynamical scalar operator with the same quantum numbers as the fermion-pair, the theory experiences a BCS/BEC crossover controlled by the relative scaling dimensions. We find the novel result that this BCS/BEC crossover exposes resonances in the canonical expectation value of the scalar operator. This occurs not only when the scaling dimension is degenerate with the Cooper pair, but also with that of higher derivative paired operators. We speculate that a proper definition of the order parameter which takes mixing with these operators into account stays finite nevertheless.


Introduction
The puzzles posed by strongly correlated electron systems have been considerably illuminated in recent years by the application of gauge-gravity duality. This "holography", which translates the challenging strongly coupled dynamics to an equivalent weakly coupled gravitational theory in one dimension higher, has given qualitative new insights into quantum critical transport [1,2], superconductivity beyond the weak coupling Bardeen-Cooper-Schrieffer (BCS) paradigm [3][4][5], and non-Fermi liquids [6,7].
A simple way to pose the challenge of strongly coupled systems is that the familiar weakly coupled particles no longer exist as controlled excitations in this regime of the theory. Our microscopic understanding of the observed macroscopics in condensed matter usually rests on the notion of an electron(ic quasi)-particle -a charged spin 1/2 fermionas the fundamental degree of freedom. The theory of Fermi-liquids and the BCS description JHEP09(2014)106 of superconductivity are good examples of such weakly coupled systems. Even in strongly correlated phases, parts of this electron quasi-particle picture survive. The transition from such a strongly correlated phase to a superconducting phase is still thought to arise from fundamental electron pairing at the microscopic level. After all, these are the only relevant charge carriers in the system. The open puzzle in strongly correlated electron systems such as high Tc superconductors is the nature of the "glue": the interaction that allows pairs to form.
In this article we take this suggestion that simple pairing mechanisms should survive in strongly coupled systems to heart. While staying ignorant on the glue, it is a very natural step to incorporate the BCS theory in the holographic framework. A straightforward reason to do so is to use this very well understood standard theory of superconductivity as a benchmark and inroad into a deeper understanding of holographic fermions. Although AdS/CF T models of superconductivity that have been constructed up to now are quite successful in capturing the main universal properties of real superconductors, they describe physics on the Landau-Ginzburg level of a scalar order parameter. In doing so it manifestly cannot reveal details of the underlying microscopic mechanisms that drive the superconducting instability, but it also ignores the Cooper pair origin of the order parameter. Our specific question here is whether holographic BCS can fill in the latter gap while being agnostic on the former, and serve as a good foothold for further research on this topic.
The most straightforward implementation of Cooper pairing in holography is to incorporate an attractive four-fermi interaction in the gravitational dual theory. In essence one now has a weakly coupled BCS interaction in the dual description of the strongly coupled theory. Pairing instabilities in this set-up were studied in [8], and the formation of a gap in the fermion spectral functions in a fixed Landau-Ginzburg holographic superconductor background, charactertistic of the broken groundstate, was shown in [9]; see also [11].
Both these studies consider the fermions as probes. Since then our understanding of holographic fermions has increased and we now understand that some of the peculiar holographic effects, in particular the non-Fermi-liquid behavior, arise from a coupling to an interacting critical IR [12]. We shall use that improved understanding to go beyond the probe limit and study the full condensation of any paired state, its subsequent groundstate and the self-consistent gap in the fluctuations around it. One way to fully treat the fermion physics is to approximate the fermions in the gravitational dual in a macroscopic fluid limit [13,14]. In this electron star approximation it is possible to understand the full macroscopic features of the system as it includes gravitational backreaction. A companion article takes this approach [15]. The drawback of the fluid limit is that it essentially describes a system with infinitely many Fermi surfaces -one for each mode in the extra radial AdS direction. This is very unusual from a condensed matter point of view.
Here we pursue an approach that allows us to concentrate on the dynamics of a single Fermi surface. This requires us to consider the fermions quantum-mechanically. In the straightforward holographic set-up this "quantum electron star" is fraught with subtle issues due to zero-point energy renormalization and its effect on the gravitational background [16,17]. From the perspective of the field theory side this difficulty is the interaction with the large number of surviving IR degrees in addition to the Fermi-surface quasiparti-JHEP09(2014)106 cle. As our first goal is to simply recover the physics of regular BCS in the dual description, the straightforward solution is to lift these extra IR degrees of freedom, and start with a regular confined Fermi-liquid. This can be done by the addition of a hard-wall [12,16]. This also discretizes the infinite number of Fermi surfaces dual to each radial mode that the AdS theory describes. We then tune the chemical potential such that only a single Fermi surface is occupied. This has the added virtue that the gravitational backreaction will be small, and we are allowed to neglect it. In this straightforward set-up the bulk AdS computation reduces to a standard Hartree BCS calculation but with relativistic fermions in an "effective box" that is spatially curved. This has several technical consequences: working in d = 3 + 1 bulk dimensions, there is an effective spin-splitting in that the up and down spin fermions have different Fermi-momenta [18,19]. Furthermore the non-trivial wavefunctions of the fermions enter into the gap equation. Accounting for this, we shall show that in this hard wall model conventional BCS maps cleanly between the dual gravitational theory and the strongly interacting field theory on the boundary.
To connect this closer to previous study [9] including the standard Landau-Ginzburg holographic superconductor, we next allow the gap-operator to become dynamical: i.e. we introduce a kinetic term for the scalar field in the gravitational bulk. The interpretation of this in the dual field theory is that we have explicitly added an additional charged scalar operator in the theory, that can independently condense. The characteristic quantum number of this new scalar operator in the strongly coupled critical theory is its scaling dimension. Following the well-known AdS/CFT dictionary, this translates into the mass of dual scalar field in the gravitational bulk. For very high mass/dimension the field/operator decouples and we have the conventional BCS scenario constructed earlier. For low masses, the field/operator starts to mix with the Cooper pair operator, and we observe a BCS/BEC crossover. Here we find a novel result. When the operator dimension is strictly degenerate with the that of the Cooper pair, the expectation values of each diverge. Nevertheless their sum -equal to the order parameter -and the gap stay finite. In effect the extra scalar and the Cooper pair act as a π-Josephson pair in that the relative phase of the condensates is opposite. 1 However, when the operator dimension is degenerate with that of a higher derivative cousin of the Cooper-pair -higher conformal partial wave -there is another resonance where the naive expectation values of each diverge. Arguably the gap should stay finite for any value of the scaling dimension. A direct application of AdS/CFT rules does not extract the gap cleanly and indicates that a clearer definition of the order parameter vev is needed in the AdS/CFT dictionary. We will address this in future work. Here we conclude by shownig that one can easily construct an expression that has the right order parameter property in that it stays finite. This postulated gap shows a clean BCS/BEC crossover.
2 Review of Fermion spectra in the AdS dual: spin splitting To start we shall recall a lesser known point of spectra of holographic fermions: the spectra depend on the spin [18,19]. The spectra follow from the simplest AdS model of fermions,

JHEP09(2014)106
Einstein Dirac-Maxwell theory -we shall add the BCS interaction in later. The action is Here the covariant derivative equals D µ = ∂ µ + 1 4 ω ab µ Γ ab − iqA µ , and Ψ = iΨ † Γ 0 . For the background we choose a pure AdS 4 spacetime with AdS radius L equal to one, and cut-off by a hard wall at a finite value of the holographic direction z = z w .
We shall consider a large charge q κ where it is consistent to ignore gravitational backreation. The cut off at z w plays a double role. Together with the AdS potential well, it renders the interval along the holographic coordinate 0 < z < z w effectively finite. This leads to quantization of fermionic energy bands ω n (k) (where n is the discrete band number). Therefore, on the one hand, we have well-defined sharp long living quasiparticles, and on the other hand the removal of the geometry beyond z w corresponds to a gapping out of normally present low energy deconfined degrees of freedom. This fundamental gap is also present in the fermion spectra itself. See figure 1(a). In this set-up we can arrive at the dual description of a single Fermi liquid by tuning the chemical potential such that exactly one band is partially occupied [12]. The charge density produced by the occupied fermions backreacts on the gauge field and its profile and the subsequent adjustment in the fermion spectra can be determined in a self-consistent Hartree manner [12]. Changing z w changes the size of the gap and the level spacing (larger values of z w correspond to smaller gap), but does not affect the qualitative picture. Only for strictly infinite z w do we enter a new critical regime which requires a completely different analysis [16,17]. We will keep z w finite throughout and therefore set z w = 1 for most of the remainder without loss of generality. Since all our computations will only depend on the combination qA 0 , we also set q = 1 in every numerical calculation from hereon.
As we shall review now, due to the spin carried by the relativistic fermions there are actually two Fermi liquids. Moreover, the (background or self-generated) electric field provides a spin-orbit coupling that renders them slightly non-degenerate in the curved background geometry. In addition the lowest energy state is at a non-zero momentum value; this is known as the plasmino mode [18,19]. This non-degeneracy of the different spin Fermi surfaces will be important in that it leads to a more complex pairing of the fermions.
The spectrum of the fermions is given by normalizable solutions to the Dirac equation. Eliminating the spin connection by rescaling Fourier transforming along the boundary directions, and making the assumption that the only non-vanishing component of the vector potential is A 0 , the Dirac equation reduces to the eigenvalue problem

JHEP09(2014)106
Hereinafter we use tangent-space gamma-matrices, and i = 1, 2 refers to the boundary spatial indices. Due to the impenetrability of the hard wall we choose the canonical momenta to vanish at z = z w : At the boundary z = 0 we demand that the fermion and scalar fields are normalizable (i.e. vanish sufficiently fast), and the boundary value of the gauge field sets the chemical potential in dual field theory: A 0 (0) = µ. The fermion spectra are determined together with the gauge field profile selfconsistently by (numerical) iteration [12]: solve the Dirac equation for a given gauge field profile (for the initial profile A 0 (z) = µ). Then solve Maxwell equations ∇ µ F µν = −iq ΨΓ ν Ψ with the source determined from the normalizable wave-functions. This gives a new gauge field profile for A 0 , etc. the result converges to a self-consistent solution after a few iterations (figure 2).
The interesting feature of the spectrum is that each band has a fine structure. To understand the origin of this splitting we examine profiles of the two spinor modes corresponding to the first band. Fermion spectra are frequently analyzed using rotational invariance to rotate the momentum k i parallel to the x-axis and choosing an appropriate basis of the gamma matrices one can simplify the problem [6]. It will, however, be useful for us to keep the rotational symmetry manifest. Our objective is to separate the radial evolution of the fermion from its spinorial structure as much as possible. We can solve the Dirac equation (2.4) with the ansatz where A ± z, | k| and B ± z, | k| are functions of the radial coordinate and u ± k i are spinors (with unit norm) independent of z. The latter are defined by the following properties wherek i is a unit (boundary) vector pointing to the direction of the momentum k i . In the basis (3.4) (which we will use later in this paper) and with a momentum parallel to the x-axis u + (u − ) is the spinor with only fourth (first) nontrivial component. The Dirac equation implies that Provided the electrostatic potential is regular near the AdS boundary at z = 0, the asymptotic behavior of the solution is  . We can observe that degeneracy of the two spin states is resolved, and state of a minimal energy is at non-zero momentum. The red and blue curves correspond to positive u + (k) and negative u − (k) modes respectively. (When the electric field is self-generated by the fermions the effect is smaller, see figure 2(a)) Normalizable solutions are those with a = 0. Note that the scaling dimension of the original fermion is ∆ Ψ = m Ψ + 3 2 and we obtained the powers of z above as a result of the rescaling (2.3). In the IR, the boundary condition (2.5) implies that A ± z w , | k| = 0.
In the absence of an electric field (i.e. A 0 (z) is constant), the positive and negative modes have the same energy. In this case we can actually solve our problem exactly in terms of Bessel functions [12]   A ±,n z, | k| with the dispersion relation ω n = (j n /z w ) 2 + k 2 − qµ. Here j n is the n-th zero of the Bessel function J m Ψ −1/2 , and N ± is the normalization constant.
However, in the presence of an electric field in the bulk (A 0 (z) = 0) the positive and negative modes no longer have the same energy anymore. The reason is that the densities of the two modes (2.10) have different radial profiles. The "effective chemical potential" A 0 (z) felt by each mode is therefore different, if the gauge field has a non-trivial z dependence, and this results in a different energy shift for the two modes ( figure 1(b)).

Majorana interaction
To study pairing driven superconductivity we now add a quartic contact fermionic interaction in the bulk of AdS: ψ C here is a charge conjugated spinor, and the z 6 factor is due to the rescaling (2.3). One can also consider the naive relativistic generalization of the Cooper pair ψ C ψ. However to boil down to standard BCS in non-relativistic limit, where the coupling occurs in s-wave channel between states time-reversed to each other, the unique Lorentz invariant term is actually the Majorana coupling ψ C Γ 5 ψ (see e.g. [20] for details). We therefore focus only on this term. As was shown in [18] the direction of the spin of each of the slightly offset modes is perpendicular to the momenta and the two modes have opposite spin. The zero-momentum pairing therefore occurs between opposite spin, without any mixing of the two fermion modes, see figure 3.
To analyze the interacting theory, we perform the standard Hubbard-Stratonovich transformation with the introduction of an auxiliary the scalar field φ(z) with charge q φ = 2q dual to the superconducting condensate. The scalar part of the action thus takes the form This is the theory studied in [8,9] with the kinetic term for the scalar turned off. We shall reintroduce this kinetic term in section 4.3.

Nambu-Gorkov formalism
The resulting system differs from standard BCS in that, as before, we are including the backreaction of the finite density fermions on the gauge field. Assuming translational invariance in the boundary directions, and restrict the scalar and the gauge field to depend only on z-coordinate, the holographic BCS system is formed by The fermionic expectation values are assumed to only depend on z as well; they are averaged over all other directions. To compute them, it is convenient to rewrite the action in a quadratic form in terms of the Nambu-Gorkov spinors. We choose the following basis of gamma-matrices (3.4) and rewrite the fermionic part of the action as

JHEP09(2014)106
where the Nambu-Gorkov spinor χ equals Taking the pure AdS metric (2.2) explicitly, and using rotational invariance of the problem to set k y = 0, the kinetic matrix K equals The fermionic expectation values can be written in terms of the Nambu-Gorkov Green's function, which satisfies the equation Note the additional factor of iΓ 0 in our definition. We determine the Green's function by spectral decomposition. For this we solve the Dirac eigenvalue problem in presence of both the (backreacted) scalar and gauge field H(i k, z)χ k,n (z) = ω k,n χ k,n (z). (3.11) Note, that the Nambu-Gorkov formalism flips the signs of some pieces of the spectrum. Figure 4(a) shows how the two low-lying energy bands in figure 2(a) look like in the Nambu-Gorkov formalism. It is convenient to write (3.11) in terms of (α 1 , α 2 , α 3 , α 4 ) = (χ 1 , iχ 2 , χ 3 , iχ 4 ). In this way the redefined "Hamiltonian" H is real (but we will still denote it with H).
We will construct the spectrum numerically, but it is instructive to first consider a toy example. We wish to show that the fermion spectrum becomes gapped in the presence of a condensate for φ. Consider the special case when the gauge field is constant A 0 = µ, and the scalar field profile is linear φ(z) = z. Then it is possible to solve the Dirac equation exactly, and the dispersion relation (corresponding to the first band) takes the form (figure 4(b)): where j 1 is the first zero of the Bessel-function J m Ψ −1/2 . We visibly see the eigenvalue repulsion responsible for the opening of a gap.

Perturbative calculation of the scalar source
In the Nambu-Gorkov formalism it is straightforward to compute the form of fermionic bilinears sourcing the electric and scalar fields (see appendix A for details).
where the sum is over the various bands (i.e. radial modes). The sum in the Cooper pair condensate needs to be cut-off at a momentum scale Λ in order to be well-defined. This momentum cut-off corresponds to an energy cut-off ω D . 2 From now on we will be using real coupling constant η * 5 = η 5 . A direct discretization of the momentum integral in (3.14) is not the most reliable way to numerically computing the fermionic source for the scalar field because contributions from different momenta are sharply peaked around the Fermi surfaces. For higher numerical accuracy and analytical control we solve (3.11) perturbatively in the scalar field. For this we split the Hamiltonian into an unperturbed piece and an interaction piece H = H 0 + V , H 0 = H| η 5 =0 . The typical spectrum for the unperturbed operator looks like the one in figure 4(a). With our choice of Gamma-matrices, the eigenspinor with the unperturbed energy ω (0) k and momentum parallel to the x-axis takes the form (we omit the band index)

JHEP09(2014)106
where ξ k is a two component spinor. There is also a mode k , for which only the lower two components are non-zero. Using nearly degenerate perturbation theory we find the matrix-element controlling the effect of the scalar field: The new energy levels are so the size of the gap is V k F . We show in the appendix B that the scalar source has the following form in terms of the unperturbed wave-functions (considering only one fermion mode): 3.4 Analytical study of the non-dynamical scalar: double gap equation The key difference is the way the spatial profiles ξ k of the fermion wavefunctions modify both the gap V k and the spatially varying profile of the pairing vev ψC Γ 5 ψ . Since the AdS geometry together with the hard wall confine the wavefunction, what we have essentially done is solve a relativistic BCS in a non-trivial potential. There is one additional subtlety, in that the Fermi surfaces corresponding to the updown spin are slightly split. Assuming, as is conventional, that the cut-off frequency is small enough, we are allowed to approximate V k and ξ k by their values at the Fermi surfaces. Doing so we can solve the gap equation A brief inspection reveals that the gap equation only depends on the dimensionless combinations η 5 m φ and η 5 ω D (these would be η 5 m φ L and η 5 ω D L if we had not set the AdS radius to unity). In appendix B.2 we show that the solution of the gap equation can be found in a form of linear combination of the two fermionic wave functions (up to an additional z 3 factor) For C 1 C 2 (C 2 C 1 ) the condensate profile is more similar to the wave-function at the first (second) Fermi surface. We obtain the coefficients where x is the ratio of the two gaps x = V 1 /V 2 , satisfying the following equation x log x. In figure 5(b) we show the perturbative solutions to the gap equation for µ = 4.5, q = 1, m Ψ = 1 and for two different couplings. (In principle there are two solutions but one of these contains a node and is presumably energitically unfavored). We can see a cross-over when we tune the coupling η 5 /m φ (see also figure 6). For small (large) coupling the profile of the condensate is dominated by ρ 2 (ρ 1 ). Note that the gap at the first Fermi-surface (with fermion wave-function ρ 1 ) is always smaller than the gap at the second Fermi-surface.
The analysis above is all from the perspective of the bulk AdS physics. All the data of the dual strongly coupled field theory is directly inferred from it. The spectral condition for a normalizable mode is the same [12], hence a gap in the bulk spectra equals a gap in the boundary fermion spectrum. The CFT order parameter is by construction the leading non-zero component of the fermion bilinear vev O U(1) = lim z→0 z −2∆ Ψ Ψ C Γ 5 Ψ , where ∆ Ψ is the scaling dimension of the single trace fermionic operator O Ψ dual to the AdS Dirac field (each normalizable fermion wavefunction behaves as z ∆ Ψ ) [21,22]. We thus neatly see how a bulk BCS coupling holographically encodes standard BCS in the dual CFT.

Fermionic ordering in holography
To establish a closer connection to previous works [9,10] on fermionic aspects in holographically ordered ground states, we now introduce by hand a kinetic term for the scalar field φ. With this additional bosonic degree of freedom we can describe a BCS/BEC crossover.
Physically it corresponds to the situation where the coherence length (which is inversely proportional to the energy gap of the fermion spectrum) of the Cooper pair is comparable to the relevant cut-off (1/ω D in our case). In flat space the relativistic BCS/BEC crossover has been studied in [36,37]. From the dual boundary field theory perspective this new ingredient corresponds to an explicit scalar operator of scaling dimension We reserve the symbol ∆ for the scaling dimensions of operators. It is not to be confused with the value of the gap. The mass (scaling dimension) will be our main tuning parameter. This (UV) quantity determines the IR mass of the scalar (together with the hard-wall position z w ), but also as we will see it mixes with the fermion bilinear operator in the UV. Again assuming translational invariance in the boundary directions, the bosonic equations now take the form where q φ = 2q. In addition one has the Dirac equation

JHEP09(2014)106
through which one defines the bulk expectation values on the right hand side. Here K(φ, A 0 ) is the kinetic matrix in (3.7), The distinction between the model with a dynamical and non-dynamical scalar field is two-fold: (1) Although physically the order parameter in the broken state cannot distinguish between a fermionic Cooper pair origin and a condensed scalar, in this holographic model they mathematically arise at different orders in the 1/N expansion. Recall that the coupling constant expansion in AdS/CFT maps to the 1/N matrix expansion of the dual field theory, whereas each AdS field is dual to a single trace composite operator. A Cooper pair is thus dual to double trace operator in the dual field theory which are always 1/N suppressed. This distinction is the same distinction between classical spontaneous symmetry breaking in a scalar field theory, and "quantum pairing" in BCS.
(2) Physically, strictly put the scalar is an additional degree of freedom (it will show up in the free energy). If the energy gap of the fermions is comparable to the relevant cutoff, one should indeed introduce this operator separately. In this "strong coupling" (equal to small coherence length) limit, the dynamical scalar field can condense by itself. In the formulation here this is controlled by its mass. For high mass the field should decouple. This is dual to the statement that in the dual field theory the corresponding operator will have a very high scaling dimension and become extremely irrelevant. All the IR dynamics is then controlled by the fermions and we recover the standard BCS of the previous section. For low mass, however, the boson dynamics will start to compete with the fermion pairing and rapidly take over the symmetry breaking dynamics in the IR.
For concreteness recall that in the previous section we derived (among other quantities) the size of the gap in the BCS limit. We can see from (B.12) and (B.17) that the gap to cut-off ratio (V /ω D ) sensitively depends on γ, which is a function of the dimensionless coupling η 5 , the Fermi momentum, and the mass of the scalar. If the latter is big enough the gap is exponentially small as compared with ω D , and the result remains valid even after introducing the boson kinetic term. Tuning the scalar mass therefore controls a crossover between pure BCS theory and a classic BEC spontaneous symmetry breaking. Qualitatively one can thus consider the mass/scaling dimension of the scalar operator as a proxy for the coherence length of the Cooper pair. When it is large, the dynamics is pure BCS; as it becomes comparable to and smaller than the relevant cut-off, one should introduce the paired operator independently. Writing out the spin components explicitly the full system of equations that we are attempting to solve is dk|k|Θ (ω k,n ) (α k,n,1 α k,n,4 − α k,n,2 α k,n,3 ) , Here all fields depend on only on the radial direction z. For completeness we recall boundary conditions for each of the fields. At the impenetrable hard wall all canonical momenta should vanish. For the bosons this means for the fermions this can be achieved by the choice At the AdS boundary, all field should be normalizable: they should vanish as a positive power of z. (For two of the fermion components this is automatic, see eq. (2.9)). We will approach the fully interacting scalar-fermion system in three steps: we first set all fermions to vanish and construct the purely scalar holographic superconductor. Next we include fermions, but hold the BCS coupling η 5 = 0; this exhibits bose-fermi competition in the system. Finally we will analyze fully interacting system at η 5 = 0. Details of the numerical calculations are discussed in appendix C.

Purely scalar holographic superconductor
First, as the scalar field in our system is a fully dynamical degree of freedom, it should condense for small enough mass even in absence of fermions [4,23,24]. This hardwall superconductor will be useful for later comparison.
Since we consider a pure hardwall AdS 4 spacetime without a black hole horizon, we study a T = 0 groundstate as a function of the mass/conformal dimension of the scalar field/dual scalar operator. Any phase transition is therefore of quantum origin. Note that the hard gap due to the hardwall directly implies that the physics is the same for any temperature T < 1/z w . Only when T > 1/z w will the the black hole horizon become relevant to the geometry, see e.g. [25].
The numerics of the pure scalar system is particularly simple as there is no need to solve the integro-differential equations iteratively. Varying the scalar conformal dimension we indeed find a condensate value below a critical value ( figure 7). We see a sharp second order phase transition as expected for spontaneous symmetry breaking. Scalar operators with smaller conformal dimensions (dual to lighter bulk scalar fields) are more likely to condense and yield an order parameter with higher density.  dimension of the fermionic operator ∆ Ψ = m Ψ + 3/2 = 5/2 we obtain a scalar condensate shown on figure 8. Comparing, the two condensate values become identical with the pure hardwall superconductor without fermions for low enough ∆ φ . For these values the bulk scalar field is so light that it consumes all the energy in the system. Ceteris paribus we would need a higher chemical potential to make fermions occupy the first band and backreact on A 0 .

Bose-Fermi competition
At larger values of ∆ φ there is still a scalar condensate, but it is suppressed compared to the pure hardwall superconductor ( figure 8(a)). This can be easily understood in terms of canonical ensemble. For fixed total electromagnetic charge of the system, adding new constituents (fermions) would redistribute the available charge ( figure 8(b)) and the con-

JHEP09(2014)106
densate of the original degrees of freedom would be suppressed. This effect has also been observed in a holographic set-up where the fermions are approximated in the fluid [26,27]

A dynamical BCS scalar and a BCS/BEC crossover
Now we analyze the most interesting case and include the full dynamics for the scalar field φ. Let us give another reason why this is quite natural from the field theory perspective. The evolution in the radial direction in AdS captures the (leading matrix large N contribution to the) RG flow of the corresponding operator in the field theory. The BCS gap, proportional to the vev of the scalar field is certainly sensitive to the RG scale. Hence one expects it to change dynamically as a function of the radial direction. Strictly speaking the double trace pairing operator which sets the value of the gap is a subleading operator in large N and any running that deviates from its semiclassical scaling is therefore a quantum effect in the AdS gravity theory. This is the situation we studied in section 3.4. The natural value of the mass (conformal dimension) of the Hubbard-Stratonovich scalar field in the BCS limit is of order N . When the scalar mass become small enough (of order N 0 ), the pair operator will become dynamical and qualitatively it ought to be given by the dynamical scalar we study here.
We will see a very interesting effect occurs in doing so. Because the scalar is sourced by the Cooper pair condensate, this changes near-boundary fall off of φ, and the standard holographic prescription for boundary field theory condensates has to be modified. Without the presence of a Cooper pair condensate, the zero momentum scalar mode equation in AdS 4 is a homogeneous (linear) differential equation Near the AdS boundary its solutions have the following form (4.10) and in the standard quantization scheme the coefficient A of the non-normalizable solution corresponds to the source J O φ for the operator O φ dual to φ, and the coefficient B of the normalizable solution to the vev O φ . Spontaneous symmetry breaking due to a condensation of the operator occurs for a solution in the absence of a source, i.e. with A = 0 as a boundary condition. For the interacting scalar-fermion system this simple one-to-one correspondence between bulk asymptotics and boundary condensates needs modification. We must now consider the inhomogeneous differential equation the particular solution will behave in the same way (assuming 2∆ Ψ = ∆ φ ): This particular solution will control the dominant normalizable near boundary behavior for ∆ φ > 2∆ Ψ . This raises the question what we should use as the vev for the corresponding operator. The canonical AdS/CFT prescription no longer gives a viable answer. Let us exhibit this in detail. As an aside, note that the near-boundary behavior of the fermions does not change provided the solution for φ(z) is normalizable.
Denoting the coefficient B of the normalizable homogeneous solution with B = H 1 we extract these coefficients from numerical solutions to the scalar and fermionic equations. (see figure 10,11). Immediately noticable are the singularities at ∆ φ = 2∆ Ψ and ∆ φ = 2∆ Ψ + 2. Strictly speaking when ∆ φ = 2∆ Ψ + n the expansion (4.13) breaks down and the solution has an extra logarithmic term The singular divergence of coefficients is a precursor of this logarithm. There is no singularity at 2∆ Ψ + 1 because P 2 happens to vanish in our case. The indisputable presence of these singularities or resonances can be readily seen by considering a simplified version of the scalar equation. Computing the series solution to the equation one directly finds the "resonances" , .
Note that they are Feschbach-like resonances in that the singularity is a single rather than a double pole. We can see that at finite chemical potential P 3 is singular not only around ∆ φ = 2∆ Ψ + 2 but also around ∆ φ = 2∆ Ψ + 2. The question is how to extract the information of the strongly coupled dual field theory from this asymptotic behavior of the AdS scalar wavefunction. Despite these singularities in the coefficients, by construction the bulk scalar wavefunction is regular at all points JHEP09(2014)106 ( figure 9). It is therefore physically natural to have regular observables in the boundary field theory as well. There are two obvious points to make here.
(1) When the system is far from the resonance points, the mixing between fermion bilinears and scalar is small, and approximately we can regard them as two different physically distinguishable primary operators (thought in a case of a real condensed matter system it might be hard to suggest an experiment which could differ between them). However, if we look at the vicinity of a resonance point, the mixing effect is substantial, and the "physically real" operators are linear combinations of those two, so one can not experimentally distinguish between O φ and O ψ 2 .
(2) Mathematically, the regularity of the bulk solution directly implies that the homogeneous component H 1 must have a similar resonance but with an opposite sign. Therefore it can not be the correct value of the order parameter even in the case ∆ φ < 2∆ψ.
An obvious and physically motivated choice is to postulate that the actual order parameter is the simply the sum of the naive condensates, with the Cooper pair condensate S 1 renormalized to P 1 : i.e.
Taking this linear combination does in fact lead to a cancelation of "resonances" and a smooth function at ∆ φ = 2∆ Ψ (see figure 12). However, the reflection of the next resonance ∆ φ = 2∆ Ψ + 2 in the homogenous solution H 1 remains. Likewise, a similar partial resolution occurs for the linear combination H 1 + P 3 . These "resonances" in P i and their cancellation by (part of) the homogeneous solution H 1 will in fact occur at every order of the expansion from the AdS boundary z = 0. It hints that the proper definition of the superconducting condensate should be given by a ∆ φ -dependent linear combination of the homogeneous H 1 and particular coefficients P 1 , P 3 , P 5 , ... that is regular for all ∆ φ . One can readily construct such a combination, e.g.
see figure 14(a). as a demonstration of the existence of a non-singular combination; though there is no proof at all that this constitutes the actual physical observable. 4 Normally the strict application of the AdS/CFT dictionary does not assign any role to such higher order coefficients in the bulk wavefunction. It is clear, however, that the singularities arise solely from the extraction of the coefficients, whereas the full AdS wavefunctions at any finite z are regular for ∆ φ = 2∆ Ψ + N. Let us now give an argument why the coefficient rule can receive modification.The right way to interpret the linear combination H 1 + P 1 is as a mixing of the two independent operators dual to the fundamental scalar operator and the bilinear (double trace) Cooper pair operator. This suggests that we should think in a similar way about the resonance at ∆ φ = 2∆ Ψ + 2. There should be another Cooper-pair like operator in the theory which mixes with the fundamental scalar, such that the linear combination that constitutes the order parameter is finite.
In AdS/CFT this connection between mixing and resonances is in fact cleanly seen in correlation functions of bilinear operators [28,29]. These bilinear operators are also known as double trace operators, since in the models where we know the dual CFT, each operator dual to an AdS field is a single trace over an N × N matrix valued combination of fields. Bilinear operators are thus the normal-ordered product of two single trace operators. Each

JHEP09(2014)106
pair of single trace CFT operators O Ψ , however, gives rise to an infinite tower of independent primary double trace operators: (4.20) These conformal partial waves are all the higher derivative bilinear operators that cannot be written as a descendant (a derivative) of the a lower order primary. All these operators have the same global quantum numbers as the simple pair operator with scaling dimension 2∆ Ψ , but increase their dimension by two integer units each time. The correlation function study [28,29] in particular shows that in the case of an interacting purely scalar bulk theory, all these linearly independent double trace primaries mix in as well and cause single-pole Feschbach resonances in s-wave scattering of single trace operators. The correspondence between the 2n difference in scaling dimension 5 between each successive primary and the location of the resonance in the leading part of the bulk scalar wavefunction supports that this mixing is the right interpretation of the resonance. 6 We do not yet have a controlled method to extract the quantative expectation value of these higher order double trace primaries from the constituent single trace fields. As usually, one should expect that the proper value of the order parameter results from the introduction of higher order boundary counterterms of the type where Ψ ± are eigenspinors of Γ 5 . To construct this correct set of counterterms and deduce the appropriate extraction of the vev in the boundary field theory is an interesting question to pursue. In order to make it possible we have to figure our what the set of primary operators is by diagonalizing the infinite tower of mixed operators [28,29]. It will be a subject of one of our subsequent papers. The conclusion is that the resulting condensate ought to be of the form in figure 14. Qualitatively this result shows the BCS/BEC crossover as a function of the scalar scaling dimension ∆ φ . Though our set up is rather abstract in that scalar field here is an additional degree of freedom introduced by hand, instead of emerging from microscopic dynamics, it captures the BCS/BEC physics. For small scaling dimension the scalar operator 5 As we mentioned one also expects a resonance at 2∆Ψ + 3 for high enough chemical potential. This is due to the effect of the electric field on the fermion wave functions. From the boundary perspective this could be a result of mixing with O ψ Jµ( ← ∂ µ − → ∂ µ )O ψ type operator which has the right scaling dimension (∆J = 2). 6 The conformal partial wave operators share a resemblance with operators relevant for Fulde-Ferrel-Larkin-Ovchinnikov pairing [30,31]. In the original FFLO set-up one considers the Zeeman splitting of spin-up/spin-down electrons and this causes an offset in their Fermi surfaces of the same form seen here.
The discussion about the mixing in of these higher order partial waves does not rely on the split degeneracy of Fermi surfaces. The mixing is therefore not correlated with an FFLO-like phenomenon.
that has a remarkable overlap with the η 5 = 0 solution at low ∆ φ as desired. All parameters are as in figure 10.
O φ dominates the Bose-Fermi competition, whereas at large scalar conformal dimensions corresponding to weak coupling regime, η 5 /m φ 1, the dynamics of the boson field are suppressed, and its order parameter expectation value is dominated by fermions as shown on figure 12. The most interesting region is just to the right of the red curve. Here there is no bosonic contribution to the order parameter, but there is an enhanced Cooper pair contribution (due to the proximity effect). This is the most notable region where we have pairing induced superconductivity in holography. At larger scalar conformal dimension the order parameter exponentially decreases with increasing of ∆ φ , although it never vanishes. In the strict m φ → ∞ limit we have the standard BCS scenario of section 3.4.
Let us finally briefly comment on the dependence on the UV cut off ω D . In the previous section we discussed that at very large bulk scalar mass all dynamics depends only on two parameters, η 5 /ω D and η 5 /m φ . For a dynamical scalar the dependence is more complicated, but we can still qualitatively infer what will happen. We know that most of the contribution to the pairing operator is located near the Fermi surfaces. Increasing ω D means taking into account states lying far away from k F 's. The physical picture will therefore only change minimally; to first approximation it can be compensated by adjusting η 5 such that η 5 /ω D stays constant. A non-trivial effect does happen when ω D becomes so large that the integral becomes sensitive to fermions in the second band (for instance, see figure 1), but this is beyond the scope of this paper.

Conclusions
We have constructed a holographic model of superconductivity which explicitly takes into account fermionic pairing driving the phase transition. In the simplest holographic models, the microscopic mechanism of superconductivity is not addressed. Specific top-down mod-JHEP09(2014)106 els may shed light on the strong coupling dynamics and a possible pairing mechanism [32][33][34][35], but generic holographic models operate at a Landau-Ginzburg order parameter level.
Even so, the physics of fermionic pairing and condensation should also be explicitly representable in holographic systems. The most straightforward way to do so is to mimic the classic BCS mechanism. This is what we have done here. By introducing an attractive four-fermion interaction in the AdS bulk, we show that this directly reduces to a pairing induced superconducting groundstate both in the bulk and the dual boundary. To cleanly separate the fermion physics, we introduced a hard wall cut-off. This essentially guaranteed this results as the low energy theory in both sides is just a Fermi liquid in the absence of the four-fermion interaction. The one technical difference with textbook BCS is the relativistic nature of the underlying fermion theory.
Next we introduced separately a kinetic term for the AdS dual of the order parameter. Physically the paired operator should become dynamical if the coherence length is much shorter than the scales of interest. One should find a BCS/BEC crossover as one tunes between these regimes. Here that control parameter is the scaling dimension of the order parameter field (relative to the scaling dimension of the Cooper pair operator). For large scaling dimension the kinetics of the dual AdS field is suppressed and we have the BCS physics found earlier. For low scaling dimension the scalar dynamics should be energetically favored compared to pairing condensation, and one should find a regular BEC (holographic) superconductor.
In observing this BCS/BEC crossover we encountered a surprise. At specific values ∆ φ = 2∆ Ψ and ∆ φ = 2∆ Ψ + 2 of the control parameter the independent scalar O φ and pairing O Ψ C O Ψ vevs diverge. In fact the naive order parameter O φ + O Ψ C O Ψ remains divergent at ∆ φ = 2∆ Ψ + 2. The mathematics is clear and suggests that these divergences can also occur at higher value of the scaling dimension. Physically, a plausible explanation is that higher order primaries that arise in the OPE of the product of two single fermion operators, mix in with the scalar vev and the lowest order To establish this concretely requires a more detailed study of single and double trace operator mixing in AdS/CF T . We aim to address this in a future publication. We can nevertheless readily construct an extraction rule for a finite order parameter that interpolates between the BCS and BEC regimes.
In both aspects the physics that holographic system describes is very conventional. It is again an excellent proving ground for AdS/CF T that it does so, but by construction it does not uncover any unconventional or exotic physics. The main reason it does not do so is the presence of the hard wall. It ensures that the groundstate dynamics closely follows standard Fermi liquid and Landau-Ginzburg theory. It would be very interesting, but technically challenging [16,17], to try to remove the hard wall. This would reintroduce the low energy dynamics that could yield exotic and novel behaviour. In particular, it might be an important step towards a holographic fermionic theory of unconventional superconductivity.

JHEP09(2014)106
where χ k,n (z) solves the Dirac equation (3.11) and form an orthonormal basis We can immediately perform the ω integral to get (supposing that t > t ) Substituting this into (A.5) we obtain (3.13) and (3.14).

B.1 Perturbative fermion spectrum, AdS-gap equation
We will solve the fermionic equation of motion (3.11) perturbatively in the scalar interaction and determine the gap equation. It is convenient for this to write the eigenvalue problem in terms of (α 1 , α 2 , α 3 , α 4 ) = (χ 1 , iχ 2 , χ 3 , iχ 4 ). The redefined Hamiltonian is real (but we will still denote it with H).
Our Hamiltonian can be split as H = H 0 +V , where H 0 = H(η 5 = 0). The perturbation is coming from the Majorana coupling where is the 2x2 matrix = iσ 2 The solution of the unperturbed problem for a given momentum takes the form We will focus on n = 1 and will omit this index. When doing the perturbation theory we should be careful because near the Fermisurface different bands are crossing each other. Therefore we start with two modes with unperturbed energy ω

JHEP09(2014)106
The new energy levels are so the size of the gap is V k F . The normalized wave-functions are Using this perturbative result we can express the scalar source with the unperturbed fermion wave functions: Here Λ(ω D ) is a momentum cut-off corresponding to the energy scale ω D . In our numerics we sample discrete number of momenta and sum over it. In order to capture the contribution around k F accurately we can use the following discretization

JHEP09(2014)106
In this case the perturbation matrix element is V 1 = 2η 5 (C 1 I 11 + C 2 I 12 ) , V 2 = 2η 5 (C 2 I 22 + C 1 I 12 ) , (B.12) where I 11 = In the limit of ω D η 5 our gap-equations take the following form For the ratio x = V 1 /V 2 we obtain We can now solve our equations easily to obtain C Numerical methods

C.1 General strategy
To solve the equations (4.5) numerically, we resort to an iterative Hartree resummation: • At a constant A 0 = µ and zero scalar field, we find the unperturbed spectrum of fermions. As a result we get a set of fermionic wavefunctions for a discrete array of energies and momenta (k i , ω n,i ).
• With these wavefunctions we construct the source terms on the right hand side of the first two equations in (4.5) and solve for A 0 (z) and φ(z). Both UV cut offs in both k and ω should be imposed to render the sums in the source terms finite.
• Substitute the new A 0 (z) and φ(z) into the Dirac equation and find the new spectrum.
Once the system converges sufficiently, we can extract the information of the dual theory by a fit to the near boundary behavior of the resulting wavefunctions.
We have optimized our numerics in several ways: The most time-consuming part of the algorithm is the repeated calculation of the Dirac fermion spectrum. A significant improvement is obtained using the perturbative prescription described in a previous section.

JHEP09(2014)106
We exclude the φ(z) field from the Dirac equation, and instead of four coupled ODE we get for fermions two identical decoupled systems of a second order. Then we construct the corrected wavefunctions. In addition, we do not need to take equally dense sampling in k, because most of the fermionic spectral weight is concentrated around k F (remember that we have two slighly different Fermi momenta in the theory), and we may take sparser k-sampling away from these points without loss in accuracy.
Empirically we found that different numerical schemes to fermionic and bosonic subsystems was the most efficient. For the fermionic spectrum we use the shooting method: we impose boundary conditions dependent on a free parameter at the boundary cut off z = , and scan over this parameter to make the resulting solution satisfy physical boundary conditions at the hard wall.
However, the shooting method in the gauge field and scalar sector often leads the system to converge to some higher harmonics instead of the groundstate. The Newton method is much more stable in that case: we impose both AdS-infinity and hardwall boundary conditions at the same time, approximate differential equations by finite differences, and solve the resulting system of linear algebraic equation with a relaxation algorithm. For our purposes a grid of N p = 3000 points in z-direction (for z w = 1) was chosen, in which case the relaxation algorithm converges after 5 − 6 iterations.
Once the bulk wave functions are obtained, it is still not a trivial question how to extract the leading boundary behavior from this data. This is what contains the information of the dual field theory. The analytical puzzles related to this problem were discussed in section 4.3. Here we focus on corresponding numerical issues.
We are interested in coefficients H 1 , P 1 , P 3 defined in (4.17). The function φ(z) is known in a form of discrete list of values {z i , φ(z i )} of the length N p = 3000, therefore our accuracy is limited and naive use of the standard fitting schemes of Mathematica leads to large errors.
Instead we first determine the expansion coefficients of the fermionic bilinear sourcing the scalar field − iη 5 z 3 ψ c Γ 5 ψ = S 1 z 5 + S 3 z 7 + ... (C.1) These can be easily found, as contra to the scalar field profile the fermionic bulk wave functions are derived with a great accuracy due to the use of the shooting method. Then we use the algebraic relations (4.17) to obtain the "particular" coefficients on the base of S i .
Knowing P 1 and P 3 we can subtract these from the scalar wave function and run the Newton relaxation scheme one more time for φ(z) = φ(z) − P 1 z 5 − P 3 z 7 . (C.2) We now need to fit only for the single coefficient H 1 . This can be easily done even for moderate number of discretized points N p .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.