A Simple Holographic Superconductor with Momentum Relaxation

We study a holographic superconductor model with momentum relaxation due to massless scalar fields linear to spatial coordinates($\psi_I = \beta \delta_{Ii} x^i$), where $\beta$ is the strength of momentum relaxation. In addition to the original superconductor induced by the chemical potential($\mu$) at $\beta=0$, there exists a new type of superconductor induced by $\beta$ even at $\mu=0$. It may imply a new `pairing' mechanism of particles and antiparticles interacting with $\beta$, which may be interpreted as `impurity'. Two parameters $\mu$ and $\beta$ compete in forming a superconducting phase. As a result, the critical temperature behaves differently depending on $\beta/\mu$. It decreases when $\beta/\mu$ is small and increases when $\beta/\mu$ is large, which is a novel feature compared to other models. After analysing ground states and phase diagrams for various $\beta/\mu$, we study optical electric($\sigma$), thermoelectric($\alpha$), and thermal($\bar{\kappa}$) conductivities. When the system undergoes a phase transition from a normal to a superconducting phase, $1/\omega$ pole appears in the imaginary part of the electric conductivity, implying infinite DC conductivity. If $\beta/\mu<1$, at small $\omega$, a two-fluid model with an imaginary $1/\omega$ pole and the Drude peak works for $\sigma$, $\alpha$, and $\bar{\kappa}$, but if $\beta/\mu>1$ a non-Drude peak replaces the Drude peak. It is consistent with the coherent/incoherent metal transition in its metal phase. The Ferrell-Glover-Tinkham (FGT) sum rule is satisfied for all cases even when $\mu=0$.


Introduction
Holographic methods (gauge/gravity duality) has provided novel tools to study diverse strongly correlated systems [1][2][3][4]. In particular, after a pioneering model of superconductor in holographic methods by Hartnoll, Herzog, and Horowitz(HHH) [5,6], there have been extensive development of the model. We refer to [2,3,7] for reviews and references.
The HHH model is translationally invariant. Because a translationally invariant system with finite charge density cannot relax momentum, the HHH model will exhibit an infinite electric DC conductivity even in the normal metal phase. Therefore, to construct more realistic superconductor models, it is important to incorporate momentum relaxation in the framework of holography.
One way to include momentum relaxation is to break translational invariance by imposing explicit inhomogeneous boundary conditions such as a spatially modulated scalar field or temporal U (1) gauge field A t , which mimicks an ionic lattice [8,9]. These models successfully yield a finite DC conductivity as well as interesting features in optical conductivity such as a Drude-like peak at small ω and some scaling laws at intermediate ω. We refer to [10][11][12] for further development. This idea (in particular, explicit optical lattice) was applied to the HHH model in [13] and, interestingly, many properties of the bismuth-based cuprates were observed.
In this method, however, because of inhomogeneity of dynamic fields, the equations of motion become a complicated coupled partial differential equations(PDE). It is technically involved and less flexible than ordinary differential equations(ODE), though conceptually clear. Therefore, it will be efficient and complementary if we can analyse the system with ODEs. In this line a few ideas have been proposed and developed.
Massive gravity models [14][15][16][17] introduces mass terms for some gravitons. It breaks bulk diffeomorphism invariance and consequently breaks translation invariance in the boundary field theory. Holographic Q-lattice models [18,19] exploit a continuous global symmetry of the bulk theory, where, for example, the global phase of a complex scalar field breaks translational invariance. Models with massless scalar fields linear in spatial coordinate [20][21][22][23] take advantage of the shift symmetry 1 . This model is related to Q-lattice models. For example, a massless complex scalar with constant ϕ in (2.6) of [18] gives a massless axion linear in a spatial direction. Also there are models utilising a Bianchi VII 0 symmetry to construct black holes dual to helical lattices [27][28][29]. All these models give us ODEs and yields a finite DC conductivity as expected. Furthermore, for a large class of models, the analytic DC conductivity formulas are available in terms of the data of the black hole horizon [30].
Building on this development based on ODE, it is natural to revisit the holographic superconductor models. The superconductor model combined to the massive gravity models and Q-lattice models have been studied in [31] and [32,33] respectively. For massless scalar models, an anisotropic background case with one scalar field was addressed in [34,35]. As in non-superconducting cases, the properties of theses ODE-based superconductor models qualitatively agree to the PDE-based model with ionic lattice [13].
In this paper, we study a holographic superconductor model based on a massless scalar model for isotropic background. The model consists of two parts: The HHH action [13] and two massless real scalar fields [20]. The HHH action is a class of Einstein-Maxwell-complex scalar action with negative cosmological constant. An on shell massless real scalar field(ψ) is linear to spatial coordinate with proportionality constant β, for example, ψ = βx. To have isotropic bulk fields, the identical scalar field is introduced for every spatial direction in field theory so there is only one parameter, β, which controls the strength of momentum relaxation.
The HHH model without massless real scalar sector has been studied extensively. See [2,3] for review. It has two phases, normal metallic phase and superconducting phase. Without the complex scalar field, the black hole is Reissner-Nordstrom type and the system is in normal metal phase. With a finite complex scalar hair, the system is in superconductor phase. With massless real scalar fields, a normal metal phase still exist [20,24] and its thermodynamic and transport coefficients were studied: the DC electric conductivity [20], DC thermoelectric and thermal conductivity [30], optical electric conductivity [22], optical electric, thermoelectric and thermal conductivities [23].
Having studied the metal phase of the model, we want to investigate the superconducting phase. First, we will examine the condition in which a superconductor phase may exist at finite β. Second, we will study the effect of β on superconducting phase transition and the properties of superconductor.
Interestingly, we find that there exists a new type of superconducting phase even when µ = 0 at finite β in addition to the original superconductivity at β = 0 and finite µ. i.e. large β > T induces superconductivity as large µ > T does. Because two parameters µ and β compete each other in forming superconducting phase the critical temperature may behave differently for β/µ > 1 and for β/µ < 1. At zero temperature the β suppresses the superconductivity. At finite temperature for small β/µ the critical temperature decreases, while for large β/µ, the critical temperature increases. It is a novel feature compared to other models.
Our main tool to analyse the superconducting phase is optical conductivities: electric(σ), thermoelectric(α), and thermal(κ) conductivities. In holographic framework, the computation of these conductivities are related to the classical dynamics of three coupled bulk field fluctuations(metric, gauge, scalar fields). By computing the on-shell quadratic action for these fluctuations we can read off the retarded Green's functions relevant to three conductivities. Most papers deal with only electric optical conductivity. However, for better understanding, it will be good to have a complete set of three conductivities: σ, α, andκ. For a class of models, analytic formulas for DC conductivities are available [30,36], but not for optical conductivities. In [23], a systematic numerical method to compute all three conductivities in a system with a constraint were developed based on [37,38]. The method was applied to the normal metallic phase of our model, producing numerical conductivities reliably [23], and we will use the same method for the superconducting phase in this paper 2 .
One of the main results of [23] is numerical demonstrations of coherent/incoherent metal transition 3 . In [42] metal without quasi-particle at strong coupling was classified by two classes: coherent metal with a Drude peak and incoherent metal without a Drude peak. At small ω, when momentum dissipation is weak (β < µ), all three optical conductivities fit well to the Drude form modified by K 0 : where K 0 is the contribution from pair production affected by net charge density. For β µ, K 0 can be ignored and a modified Drude form is reduced to the standard Drude from. If β > µ optical conductivities do not fit to (1.1) and goes to the incoherent metal phase 4 , which agrees to [42]. In superconducting phase, we find a similar result. If β < µ all three optical conductivities fit to where K s is supposed to be proportional to superfluid density. In superconducting phase, it turns out that K s = 0 for σ and κ, but K s = 0 for α. We have confirmed numerically the Ferrell-Glover-Tinkham (FGT) sum rule is satisfied in all cases we considered.
where σ s (σ n ) is the electric conductivity at T < T c (T > T c ). This paper is organised as follows. In section 2, we introduce our holographic superconductor model (Einstein-Maxwell-complex scalar action with negative cosmological constant) incorporating momentum relaxation by massless real scalar fields. Background bulk solutions corresponding to superconducting phase and normal phase are obtained. By comparing on-shell actions of both solutions, we identify the phase transition temperature as a function of chemical potential and momentum relaxation parameter, which yields 3dimensional phase diagrams. In superconducting phase, we also compute condensates as a function of temperature for given chemical potential and momentum relaxation parameter. In section 3, we compute optical electric, thermoelectric, and thermal conductivities in superconducting phase and normal phase. In particular, in superconducting phase, we discuss the effect of momentum relaxation on conductivity in several aspects such as the appearance of infinite DC conductivity, Drude-nature of optical conductivity in small frequency range, two-fluid model, and Ferrell-Glover-Tinkham(FGT) sum rule. We also present a general numerical method to compute retarded Green's functions when many fields are coupled. In section 4 we conclude.
Note added: After this paper was completed, we became aware of [33] which has overlap with ours.

Metal/superconductor phase transition
We start with the original holographic superconductor model proposed by Hartnoll, Herzog, and Horowitz(HHH) [5] where F = dA is the field strength for a U (1) gauge field A and Φ is a complex scalar field. We have chosen units such that the gravitational constant 16πG = 1. The second action, S GH , is the Gibbons-Hawking term, which is required for a well defined variational problem with Dirichlet boundary conditions. γ is the determinant of the induced metric γ µν at the boundary, and K denotes the trace of the extrinsic curvature. To impose a momentum relaxation effect, we add the action of free massless scalars proposed in [20] 3) The action S HHH + S ψ yields the equations of motion for matters 5 where the covariant derivative is defined by D M Φ = (∇ M − iqA M ) Φ, and the Einstein's equation Since we are mainly interested in 2 + 1 dimensional systems, we will set d = 3 from now on. In order to construct a plane(x, y)-symmetric superconducting background, we take the following ansatz, where non-zero A t (r) is introduced for a finite chemical potential or charge density and non-zero Φ(r) will be a hair of black hole, yielding a finite superconducting order parameter. A special form of ψ I is included for momentum relaxation, where β may be interpreted as a strength of impurity. Plugging the ansatz (2.8) into the equations of motion (2.4)-(2.7), we have four equations (2.9)-(2.12). Maxwell's equation (2.4) yields where Φ can be taken to be real, since r component of Maxwell's equation implies that the phase of Φ is constant. The complex scalar field equation (2.5) becomes Massless real scalar equations (2.6) is satisfied by the ansatz (2.8). The tt and rr compo-5 Index convention: M, N, · · · = 0, 1, 2, r, and µ, ν, · · · = 0, 1, 2, and i, j, · · · = 1, 2.
nents of Einstein's equations (2.7) give (2.12) We will numerically solve a set of coupled equations of two second order differential equations(2.9-2.10) and two first order differential equations(2.11-2.12) by integrating out from the horizon(r h ), which is defined by G(r h ) = 0, to infinity. It turns out that, by the regularity condition at the horizon, six initial conditions are determined by three initial values: with A t (r h ) = 0, for a given r h , β, m 2 , and L. From regularity of the Euclidean on-shell action the Hawking temperature(T H ) is given by Near boundary(r → ∞) the equations (2.9)-(2.12) are solved by the following series expansion.
where we considered a case with m 2 = −2/L 2 to be concrete. The mass m 2 is related to the conformal dimension ∆ of an operator in the dual field theory, m 2 = ∆(∆−d) , and the scalar field falls off as r −∆ . For m 2 = −2/L 2 and d = 3, ∆ = 1 or 2 as shown in (2.19).
The coefficients of (2.16)- (2.19) are identified with the field theory quantities as follows.
where µ, ρ, , O (i) and J (i) are chemical potential, charge density, energy density, ∆ = i operator and its source, respectively. So J (i) should vanish for condensation of O (i) . See appendix A for details. χ (0) should be set to be zero to identify the Hawking temperature of the black hole with the temperature of boundary field theory. Equivalently, we may rescale the time where a = e −χ (0) /2 , to set χ (0) = 0. In practical computation, this rescaling method is much easier since we have to shoot out from horizon. Then the field theory temperature is computed as where T H is defined in (2.15). However, analytic formulas from here are presented by assuming χ (0) = 0. There are two scaling symmetries of the equations of motion. The first symmetry is 23) and the second one is By taking a 1 = 1/L and a 2 = 1/r h we may define new scaled variables with tilde.
(2.25) While performing numerics we work in terms of these tilde-variables. In practice, it is equivlalent to set L = r H = 1 in numerical computation. However, when we interpret final results, we need to scale them back carefully.

Normal (metallic) phase
Here we minimally summarize the properties of the metal phase to set up the stage for our paper, referring to [20] for more details.
The normal phase of the system corresponds to the solution without condensate 26) and the analytic solution is given by where m 0 is determined by the condition G(r h ) = 0: The normal phase of the system is described by a charged black brane solution with nonvanishing β, which describes a metal state with a finite DC conductivity. The temperature is given by the Hawking temperature (2.14) : There is an important property of the geometry which we will rely on when discussing the instability of the metal phase. In the zero temperature limit, the near horizon geometry of an extremal black brane becomes AdS 2 × R d−1 with the effective radius of AdS 2 given by Notice that even at µ = 0 case, the near horizon geometry remains to be AdS 2 × R d−1 due to a finite β. In this sense, µ and β will play a similar role as far as instability is concerned.

Criterion for instability and quantum phase transitions
Before performing a full numerical analysis to search superconducting phase, we first want to address a simpler question: When can the normal phase( (2.26)-(2.30)) be unstable by small scalar field perturbations, which may develop into a hairy black hole with nonzero Φ? To address this question we perform the analysis presented in [2,44]. We look for an unstable mode of Φ = φ(r)e −iωt of the equation (2.5) in the background (2.26)-(2.30): The normal phase is unstable if there is a normalisable solution with incoming boundary conditions at the horizon such that ω has a positive imaginary part. Because we are interested in determining the critical temperature it is enough to search a static(ω = 0) normalisable mode. Our numerical procedure is as follows. The equation (2.35) depends only on four dimensionless quantities: ∆, q, β/µ, and T /µ. For fixed β/µ, ∆, and q, we shoot from horizon to boundary keeping dialing T /µ. Near boundary, φ falls off as r ∆−3 and r −∆ . We search for the largest value of T /µ for which the coefficient of r ∆−3 vanishes. We repeat this procedure for different pairs of values of ∆, and q. Figure 1(a) shows the line of constant critical temperature (T /µ) in ∆-q space when β/µ = 1.
In Figure 1(a) there is a special curve(solid one) representing quantum phase transition at T = 0. We may understand the quantum phase transition by considering the Breitenlohner-Freedman(BF) bound. First, the effective mass of scalar field, which is read from (2.35), near horizon at zero temperature is Second, the near horizon geometry of an extremal black brane is given by AdS 2 × R 2 with the effective radius of AdS 2 given by (2.33) Third, recall that, real scalar field in the AdS d+1 space with the radius L d+1 is unstable with mass(M ) below the BF bound From these three data (2.36)-(2.38), we conclude that the scalar field is unstable near horizon ifm where m 2 L 2 must be greater than − 9 4 for scalar to be stable at boundary AdS 4 space. Alternatively, this result can be obtained from (2.35). For the scalar field near horizon of the extremal black hole, (2.35) becomes The equation (2.40) is an equation for a real scalar field in AdS 2 spacetime with the curvature radius of unity. Thus, we reach the  same conclusion as (2.39) by (2.38). From the inequality (2.39), we can infer that qualitatively a larger mass (or ∆) 6 , smaller q, or larger β/µ suppresses instability(superconductivity) at zero temperature. Quantitative phase boundaries for several β/µ are plotted in Figure 1(b), where the region above (below) the curve is stable or normal phase (unstable or superconducting) phase. It is interesting that one can tune β/µ to trigger a qunatum phase transition. For example, let us consider the system with q = 3 and ∆ = 4, which is (c) in Figure 1(b). When β/µ is small the system is in normal phase, but when β/µ is large the system is in superconducting phase. The transition occurs between β/µ = 1 and β/µ = 2. The system at (d) in Figure  1(b) must be always in the normal phase at zero temperature regardless of β/µ. However the system with q = 2, ∆ = 2 and q = 3, ∆ = 2 ((a) and (b) in Figure 1(b) respectively) must be always in the superconducting phase at zero temperature. Notice that, in this case, the superconducting transition occurs even at µ = 0, i.e. β/µ = ∞. It also can be seen in Figure 2.

Superconducting phase
In the previous subsection we investigated the possibility of superconducting phase at finite temperature and zero temperature. Based on our result on quantum phase transition, Figure 1(b), we may anticipate which values of parameters (q, ∆, β, µ) allow thermal phase transition. In this subsection we want to confirm our anticipation by constructing explicit superconducting background at finite temperature.
Our numerical analysis is performed as follows. By shooting out from horizon we find a numerical solutions satisfying (2.9)-(2.12) and boundary condition Φ (1) = 0 and 6 The mass m 2 is related to the conformal dimension ∆ of an operator in the dual field theory, m 2 =  consider the condensate of the operator of dimension two, O (2) . See (2.20). We may choose the different boundary condition, Φ (2) = 0, but we will not deal with the case in this paper. At high temperature we obtain only one solution, which agrees to an analytic solution of normal state (2.26)-(2.30). At low temperature we find another solution with Φ = 0(superconducting phase) in addition to a normal state solution (2.26)-(2.30). In this case it turns out that the superconducting solution always has a lower grand potential and becomes a ground state. The phase transition is continuous at a critical temperature(T c ). Figure 2 shows typical examples of phase diagrams for three points (a),(b), and (c) in Figure 1(b). Let us start with the point (a) and (b) in Figure 1(b). They are always in superconducting phase for all β and µ at zero temperature. As temperature increases we expect that the system undergoes a phase transition from superconducting phase to normal phase. Our numerical analysis confirms it and the phase diagram is shown in Figure 2(a)(b), where the meshed surface is the phase boundary at the critical temperature. Dark region below the surface is superconducting phase while region above the surface is normal phase. Figure  2(a)(b) focuses on the phase structure for small β. In Figure 2 The red mesh is above the black mesh, which means that a large q enhances superconductivity, as at the zero temperature in Figure 1. However, the phase transition line coincides at µ = 0, because the effect of q enters in the combination of qµ as shown in (2.35).
Notice that the superconducting phase exists even when µ = 0. In the HHH model [6] the phase transition is understood as a competition between µ and T and in this case it is  Figure 3. Condensate for three values of µ/β and ∆ = 2, q = 3, which is the case of Figure 2(a). The color here matches the color of the lines in Figure 2(a). In other words, we compute condensate along the vertical line(temperature) standing on the colored-lines in Figure 2(a) a competition between β and T . The 'pairing mechanism' of two cases must be different, because when β = 0 it would be due to a particle-particle pair, while when µ = 0 it would be due to a particle-anti-particle pair also interacting with β, which may be interpreted as 'impurity' [23]. In general, at finite µ and β, two mechanisms will compete.
This competition is reflected on the phase boundary surfaces in Figure 2. In words, the dependence of the critical temperature on β/µ is not monotonic. the critical temperature decreases when β/µ is small and increases when β/µ is large. In graphics, see the line at µ = 2 in Figure 2(d) or the lines at µ = 0 and µ = 2 in Figure 2(a). It is different from the previous studies. In Q-lattice model [32,33] and single scalar model [34] the critical temperature decreases as momentum relaxation effect increases while in ionic lattice model [13] the critical temperature increases monotonically as lattice effect increase.
The point (c) is different from (a) and (b) in Figure 1(b), in that the system at (c) could be in superconducting phase or normal phase depending on the value of β/µ. Large β/µ suppresses superconductivity and Figure 2(c) shows both quantum and thermal phase transitions. If we include negative values of β, we obtain the phase diagram Figure 2(e). It looks similar to the superconducting dome in cuprate superconductor phase diagram when we identify β with a doping parameter. Note that β is a tunable continuous free parameter in the solution, while m and ∆ is fixed in the action. We also see large ∆ suppresses superconductivity by comparing Figure 2(a) and (c).
In superconducting phase there is finite condensate, Φ (2) , which is an order parameter. For example, in Figure 3, we show the condensate as a function of temperature. At the critical temperature condensate starts forming and increases continuously as temperature goes down. At very small temperature our numerics becomes not reliable and we did not plot in that range. The plot is for ∆ = 2, q = 3 which is the case of Figure 2(a). Three values of µ/β is chosen: µ/β = 0(blue), µ/β = 1(red), and µ/β = 10(green). The color here matches the color of the lines in Figure 2(a). In other words, we compute condensate along the vertical line(temperature) standing any point on the colored-lines in Figure 2(a). The condensate increases as β increases, which agrees to an anisotropic case [34]. When β → ∞ the condensate approaches to the finite upper bound, while when β → 0 it approaches to the lower bound.
We finish this subsection by discussing on the on-shell action of the ground state. To calculate a thermodynamic potential for the black hole solutions we calculate the on-shell Euclidean action(S E ) by analytically continuing to Euclidean time(τ ) where S ren consists of four-dimensional S bulk and three dimensional S bdy : (2.42) The first three terms are defined in (2.1), (2.3) and (2.2) respectively, and the last term is the counter term for holographic renormalisation [6]: 43) which cancels the divergence of the bulk action. To fix Φ (1) on the boundary we choose (η 1 , η 2 ) = (0, −1) while to fix Φ (2) we choose (η 1 , η 2 ) = (1, 1). n M = (0, 0, 0, G(r)) is an outgoing normal vector. See appendix A for more details.
Let us first consider the Euclidean bulk action which defines L bulk . It can be computed following [6]. The xx-component of the Einstein equation gives a useful relation: where G µν is the Einstein tensor. The trace of the Einstein equation yields Thus the bulk Lagrangian is In superconducting phase After adding S E bdy the total on-shell action is finite and reads where we consider a homogeneous and equilibrium system. Thus V 2 is the volume of the system and the length of the Euclidean time circle is identified with inverse temperature(T ). W is a thermodynamic potential per unit volume: (2.50) In our case, one of Φ (1) or Φ (2) is zero and W is simplified In normal sate, where G (1) = −m 0 and χ = 0, (2.51) becomes With the relations for energy density( ), charge density (ρ), and entropy density(s) which agrees to [20].
In superconducting phase, where G (1) and A (1) are numerical values. Lack of the analytic relation such as (2.55) make it difficult to check if W = Ω. Therefore, we numerically checked and found that Ω = − T s − µρ = G (1) − β 2 r h = W. However, as far as the phase diagram is concerned, this difference does not matter. It turns out W < Ω in superconducting phase so we may use Ω as our criteria for phase transition. To study thermodynamical quantities it will be important to understand physics of the difference between W and Ω. We leave it for future study.

Optical conductivity
In this section we study electric(σ), thermoelectric(α), thermal(κ) conductivity by considering small fluctuations of relevant gauge, metric, scalar fields around the normal and superconducting background we obtained in the previous section. From here on, we set L = 1 and use the scaled variables (2.25) without tilde.

Fluctuations for optical conductivity: equations and on-shell action
Electric conductivity is related to a small bulk gauge field fluctuation δA x (t, r) of which boundary dual operator is electric current. The fluctuation is chosen to be independent of x and y, which is allowed since all the background fields affecting the equations of motion are independent of x and y. Because of rotational symmetry in x-y plane, it is enough to consider δA x . The gauge field fluctuation(δA x (t, r)) sources to metric(δg tx (t, r)) and scalar field(δψ 1 (t, r)) fluctuation and all the other fluctuations can be decoupled. Notice that h tx (ω, r) is defined to approach constant as r goes to infinity. In momentum space, the linearized equations of motion around the background (2.8) are derived from (2.4)-(2.7): where G, χ, A t , Φ are background field obtained in the previous section. For normal phase we have the analytic solutions ((2.26)-(2.30)) but for superconducting phase we have it numerically. We solve these equations with two boundary conditions: incoming boundary conditions at the black hole horizon and the Dirichlet boundary conditions at the boundary. First, near the black hole horizon (r → 1) the solutions are expanded as where ν ± = ±iω e χ(1)/2 −G (1) and the incoming boundary condition corresponds to ν = ν + . Next, near the boundary (r → ∞) the asymptotic solutions read and we fix the values of the leading terms as boundary conditions. Plugging the solutions into the renormalized action (2.42), we have a quadratic order on-shell action where we discarded the contribution from the horizon as the prescription for the retarded Green's function [45]. In particular, with the spatially homogeneous ansatz (3.1)-(3.3), the quadratic action in momentum space yields where V 2 is the two dimensional spatial volume dxdy and we omit the term proportional tx since we are studying the case with Φ (1) = 0. The range of ω is chosen to be positive following the prescription in [45].
The on shell action (3.14) plays a role of the generating functional for two-point Green's functions sourced by a (0) x , h (0) tx , and ξ (0) . We may simply read off part of the two point functions from the first two terms in (3.14). The other three terms are nontrivial and we need to know the dependence of {a tx , ξ (0) }. However, thanks to linearity of equations (3.4)-(3.6), we can always find out the linear relation between tx , ξ (0) }. We will first explain our numerical method to find such a relationship in a more general setup in the following subsection and continue the computation in that setup.

Numerical method
A systematic numerical method for a system with multi fields and constraints were developed in [23] based on [37,38]. We summarise it briefly and refer to [23] for more details.
To develop a systematic method in a general setup, let us start with N fields Φ a (x, r), (a = 1, 2, · · · , N ), which satisfy a set of coupled N independent second order diffrential equations: where the index a may include components of higher spin fields. For convenience, r p is multiplied such that the solution Φ a k (r) goes to constant at boundary. For example, p = 2 in (3.2).
Near horizon(r = 1), solutions can be expanded as where we omitted the subscript k for simplicity and ν a± correspond to incoming/outgoing boundary conditions. To compute the retarded Green's function we choose the incoming boundary condition [45], fixing N initial conditions. The other N initial conditions, denoted by ϕ a i , (i = 1, 2, · · · , N ), can be chosen, for example, as Every column vector ϕ a i yields a solution, denoted by Φ a i (r), which is expanded as where S a i are the sources(leading terms) of i-th solution and O a i are the operator expectation values corresponding to sources(δ a ≥ 1). Notice that S and O can be written as regular matrices of order N , where the superscript a runs for row index and the subscript i runs for column index.
Since all N solutions {Φ a i } becomes a basis set, a general solution yields ≡ J a + · · · + R a r δa + · · · , (3.20) with real constants c i 's. For any given J a we can always find c i (3.21) and the corresponding response R a is expressed as A general on-shell quadratic action in momentum space has the form of where A and B are regular matrices of order N . In matrix notation, J a −k can be understood as a row matrix. For example, the action (3.14) can be written in the form (3.23) with with where the index ω is suppressed. With (3.22) the action (3.23) becomes where the range of ω is chosen to be positive following the prescription in [45]. See also [29] for a careful derivation of the gauge invariant Green's function matrix. Notice that for one field case without mass term, this is the well known structure of the retarded Green's function: A = 0 and G R ∼ O/S. In summary, to compute the retarded Green's function. We need four square matrices of order N (the number of fields): A, B, S, O. A and B can be read off from the action (3.23), after taking care of all divergences by counter terms. To construct regular matrices, S and O, we solve a set of differential equations N times with independent initial conditions. The retarded Green's function is schematically A + B · O · S −1 . There is one subtlety in our procedure.
For our equations there is one subtlety caused by a symmetry of the system. Solving the equations near horizon with the expansions (3.7)-(3.9) we find that only two of a (I) x , χ (I) , and h (I) tx are free. Therefore, we cannot make a complete basis to construct a general J a . It is due to the gauge fixing g rx = 0. However, it turns out that there is a residual gauge transformation keeping g rx = 0, which is generated by the vector field ξ µ of which non-vanishing component is ξ x = e −iωt . So we may add a constant vector (S a 0 ) along the residual gauge orbit.
S a 0 = (0, 1, iβ/ω) T , (3.26) since L ξ g tx = −iωr 2 ξ x and L ξ ϕ = βξ x . Notice that S a 0 satisfies the equations of motion (2.4)-(2.7), since the residual gauge transformation leaves the linearised equation of motion invariant. Therefore, our procedure is equivalent to formally adding a constant 'solutions' of the equations to the solution set {S a i }. With the matrices S and O, which is numerically computed, we may construct a 3 × 3 matrix of the retarded Green's function. We will focus on the 2×2 submatrix corresponding to a (0) From the linear response theory, we have the following relation between the response functions and the sources: (3.28) We want to relate these Green's functions to the electric (σ), thermal (κ), thermoelectric (α,ᾱ) conductivities defined as where Q x is the heat current, E x is an electric field and ∇ x T is a temperature gradient. As shown in [2,3,23], by taking into account a diffeomorphism invariance, (3.29) can be expressed as (3.30) The comparison of (3.28) and (3.30) yields (3.31) Figure 4 shows examples of electric optical conductivities(σ(ω)) for three cases of β/µ: β/µ = 0.1, 1 and, ∞(µ = 0). This choice of parameters also corresponds to the cases with the green(β/µ = 0.1), red(β/µ = 1), and blue(β/µ = ∞) lines in Figure 2(a) and 3. The color of curves represents temperature ratio, T /T c , where T c is the critical temperature.

Electric/thermal/thermoelectric conductivites
The numerical values of temperature ratio are shown in the caption. In particular the dotted black curve 7 is for the temperature above T c , which is in metal phase and the red  curve corresponds to the critical temperature (in practice, it is slightly higher than the transition temperature.). The first row shows the real part of electric conductivity (Re [σ]) and the second row shows the imaginary part of electric conductivity (Im[σ]).
One common important feature in Figure 4 is the appearance of 1/ω pole in Im[σ] below the critical temperature, while the disappearance of 1/ω pole above the critical temperature. By the Kramers-Kronig relation, 1/ω pole in Im[σ] implies the delta function at ω = 0 in Re [σ]. Therefore, in normal phase the DC conductivity is finite due to momentum relaxation and in superconducting phase the DC conductivity is infinite, which is one of the hallmarks of superconductor.
Roughly at T /T c > 0.5, in addition to the delta function at ω = 0, there is still a finite value of DC Re[σ] in superconducting phase. It may be interpreted as a contribution from normal components in superconducting phase, implying a two-fluid model. For small ω there is a Drude-like peak in some cases in (a) and (b) of Figure 4. For smaller β/µ or at higher temperature, the peak becomes sharper. For the sake of comparison we used a similar scales in (a),(b) and (c) of Figure 4, which hides the structure of (a) near ω = 0. Therefore we zoom in Figure 4(a) in Figure 5(a). The data points well fit to the formula (solid lines) whereω ≡ ω/µ and K s and K n are supposed to be proportional to the superfluid density and normal fluid density. For β/µ = 1 the formula (3.32) does not work and the data(red, orange, green) better fit to which is shown in Figure 5(b). K 0 is related to pair creation and it was necessary also in metal phase. The existence of K 0 is most apparent in Figure 4(c), where µ = K n = 0. Indeed (3.32) is understood as an approximation of (3.33) when K 0 is negligible compared to K n τ . If β > µ, it was shown that the coherent metal phase becomes incoherent, where the Drude peak becomes a non-Drude peak [23]. Therefore it is expected that, If β > µ, (3.33) does not work in superconductor phase either. Indeed we see this is the case. When β = µ: the fit of Figure 5 (b) is not as good as (a) and starts deviating from (3.33) As temperature is lowered K n and K 0 is reduced while K s is enhanced. K 0 becomes zero at low temperature (green and blue line in Figure 4(c)), but it is possible that K n is finite even at zero T [13,31,32]. As temperature goes down (T < T c ) the spectral weight of Re[σ] is reduced while K s of Im[σ] is enhanced. This transfer of the spectral weight to K s may be quantified by the Ferrell-Glover-Tinkham (FGT) sum rule: where σ s (σ n ) is the conductivity at T < T c (T > T c ). σ n can be taken for any temperature for T > T c since the spectral weight is constant in metal phase [23]. We computed (3.34) numerically for all cases in Figure 4 and showed that the FGT sum rule is satisfied up to 10 −3 in Figure 6. At µ = 0 and β = 0 (Figure 4 (c)) there is no net charge and no Drude peak. The plots are very similar to the case β = 0 at finite µ. For example, see Figure 6 in [2] or Figure 1 in [23], where the infinite DC conductivity (1/ω pole in the imaginary part) is due to translation invariance with finite µ. Here, there is no µ and no translation invariance. So the delta function must have a different origin, which may be a new type of superconductivity. Interestingly, even in this case, the FGT the sum rule works. The deficit of spectral function may be interpreted as a deficit of 'particle-anti-particle pairs', which will condense. It may imply a new 'pairing mechanism' of particles and anti particles interacting with β, which may be interpreted as kind of 'impurity' [23].  Figure 9. The data points fit to (3.32). At µ = 0, the DC values of Re[κ] decreases quickly as temperature goes down as shown in Figure 10. The thermoelectric conductivity α vanishes, if µ = 0, which is due to particle-hole symmetry. As ω → ∞, α andκ approaches to which agrees to the Ward identity [46]. There is no 1/ω pole in Im[α] for all cases. Forκ, at µ = 0 there is no 1/ω pole in Im[κ] while in (a) and (b) there is 1/ω pole only in the superconducting phase. The appearance of this 1/ω pole inκ may be understood from (3.31) as follows. In superconducting phase there is a pole in Im[σ] so Re[G 11 (0)] = 0. From Figure 7 we see that Re[G 12 (0)] = Re[G 11 (0)]µ and Re[G 12 (0)] is finite if µ is finite. Then, Im[κ(0)] = µReG 12 (0) ω ∼ µ 2 Ks ω . Therefore, if we subtract µReG 12 (0) ω from Im[κ], 1/ω pole is expected to be disappeared, which we have confirmed numerically.

Conclusions
In this paper, we studied a simple holographic superconductor model incorporating momentum relaxation. The model consists of two parts: The HHH model [5,6] and massless scalar fields sector for momentum relaxation [20], where the strength of momentum relaxation is parameterised by β.
One of the interesting features of the model is that the existence of a new type of superconductor induced by β at µ = 0. While in the HHH model the superconducting phase transition is understood as a competition between µ and T , in this new case, it is a competition between β and T . The 'pairing mechanism' of two cases must be different. In the new type(β = 0 and µ = 0), there is no net charge and the mechanism will be due to particle-anti-particle pairs also interacting with β, which may be interpreted as kind of 'impurity' [23]. The electric optical conductivity of this new superconductor satisfies the FGT sum rule too. The deficit of the spectral weight may be interpreted as a deficit of particle-anti-particle pairs which are condensed.
With finite µ and β together, two superconducting mechanisms will compete. As a result, the dependence of the critical temperature on β/µ is not monotonic: the critical temperature decreases when β/µ is small and increases when β/µ is large. It is different from the previous studies. As momentum relaxation effect increases, in a Q-lattice model [32,33] and a single scalar model [34] the critical temperature decreases while in the ionic lattice model [13] the critical temperature increases. The condensate has the upper bound when β/µ → ∞ and the lower bound when β/µ → 0.
We studied optical electric(σ), thermoelectric(α), and thermal(κ) conductivities. For all three conductivities, at small ω, a two-fluid model with a modified Drude peak works if β/µ < 1: where K s and K n are supposed to be proportional to the superfluid density and normal fluid density and K 0 is related to pair creation. For β/µ 1, K 0 becomes negligible. The restriction β/µ < 1 for (4.1) is consistent with the result in metal phase, where coherent metal becomes incoherent metal when β/µ > 1 and the Drude peak does not work [23]. However, the Ferrell-Glover-Tinkham (FGT) sum rule is satisfied for all cases regardless of β/µ.
We have not fully analyzed the parameters, K s , K n , K 0 , τ , their physical meanings and relations. The temperature dependence of the parameters are of physical importance. For example, K n may be related to the energy gap as studied in [13,[31][32][33]. The β dependence of τ is relevant to the nature of dissipation. The correct identification of the superfluid density, proportional to K s , will be essential to investigate Homes' law [47] holographically [48]. It will be also useful to obtain analytic formula for DC conductivities from the horizon data in superconducting phase as in metallic phase [30]. While the model we considered shows many interesting features as metal and superconductor, it also has shortcomings. The electric DC conductivity in normal phase is temperature independent and the insulator phase is lacking. It would be interesting to consider superconducting phase without those shortcomings. Indeed, there is a simple generalization of the model that provides a temperature dependent DC conductivity and an insulating phase at small temperature [49], so it would be interesting to construct a superconductor model based on this background. We leave these for future study.

A One point functions
We briefly summarize how to compute one point functions holographically. Let us consider an ADM decomposition as follows: where, in our case, the shift vectors vanish and N is the lapse function given by 1 √ G(r) .
The outward pointing normal vector is n M = (0, 0, 0, G(r)) and the extrinsic curvature is given as K µν = − 1 N γ µν . Let us consider a renormalised action(S ren ) consisting of a bare action(S 0 ) and a counter action(S ct ) by taking into account the holographic renormalization: and where the variations from the counter action are The expectation values of the energy momentum tensor, the current and the scalar operators in the dual field theory can be computed as (A.10)