Higgs/amplitude mode dynamics from holography

Second order phase transitions are universally driven by an order parameter which becomes trivial at the critical point. At the same time, collective excitations which involve the amplitude of the order parameter develop a gap which smoothly closes to zero at criticality. We develop analytical techniques to study this “Higgs” mode in holographic systems which undergo a continuous phase transition at finite temperature and chemical potential. This allows us to study the linear response of the system at energy scales of the order of the gap. We express the Green’s functions of scalar operators in terms of thermodynamic quantities and a single transport coefficient which we fix in terms of black hole horizon data.


Introduction
Holography provides a powerful framework to study strongly coupled quantum field theories [1] at finite temperature and chemical potential [2]. In particular, in the large N limit, it allows us to perform exact computations which are impossible to carry out with standard quantum field theory techniques.
Second order phase transitions are driven by an order parameter which becomes non trivial in the broken phase. Depending on the nature of the system, the phase transition can be accompanied by the breaking of a global symmetry. In such a case, new gapless modes appear in the spectrum of the system in the form of Goldstone modes. The universal feature is that the phase transition is always accompanied by the emergence of a gapped mode which is exactly gapless at the critical point, the amplitude/Higgs mode (from now on we will call it the Higgs mode). At infinite wavelengths, the corresponding mode 1 is homogeneous and it is related to fluctuating degrees of freedom which are intimately JHEP08(2022)246 connected to the amplitude of the order parameter and which decouple from the entropy and charge density.
The above highlights that the Higgs mode dominates critical phenomena as its decay rate is parametrically small close to the critical point [3,4]. Despite its fundamental importance, the Higgs mode was observed in condensed matter systems only very recently. For the case of superfluids see e.g. [5,6]. In view of these results, it is very important to obtain theoretical predictions on potentially universal properties of the critical dynamics. In this paper, we will study aspects of the dynamics close to a second order phase transition by exploiting the holographic principle.
In holography, continuous phase transitions can be realised through black hole instabilities. With fixed deformation parameters and chemical potential, a unique black hole solution dominates the phase diagram at high temperatures. Lowering the Hawking temperature, certain bulk perturbations can become unstable below a critical temperature T c . Following the instability leads to new branches of black hole solutions corresponding to the broken phase. In the context of applied holography, some of the most well studied examples are the holographic superfluids [7][8][9] and phases which spontaneously break translations [10][11][12].
It is natural to ask how much we can learn about the Higgs mode by considering holographic systems near their critical points. The Higgs mode has been identified in several holographic studies in the past by using numerical techniques [13,14] or analytically when full solutions in the bulk can be obtained [15]. Given the power of holography in obtaining exact results, the ultimate goal would be to construct an enlarged theory of hydrodynamics which incorporates this universal nearly gapless mode in its description.
In this paper we will make a significant step in this direction by examining the dynamics of holographic systems close to their phase transitions. More specifically, we will obtain analytic expressions for the low frequency retarded Green's functions of scalar operators. As we show in section 5, these are entirely determined by the thermodynamic properties of the system and a single transport coefficient which is fixed by the black hole horizon data. The linear response of the system is dominated by a single pole which corresponds to the Higgs mode which decays slowly near the critical point. Our main tool will be the techniques which have been recently exploited in [16] to obtain a hydrodynamic description of holographic superfluids at zero charge density. The main contribution was to write the transport coefficients of [17] in terms of black hole horizon data.
In that construction, the basic ingredient is the symplectic current density of the bulk theory [18]. This can be considered for any pair of perturbative solutions in the bulk. Being a generalisation of Liouville's theorem for classical mechanics, the symplectic current is divergence free. This feature will be very important to us since it allows us to extract useful information for an unknown perturbative solution when it is parametrically close to a known one.
In our context, the role of the known solution will be played by the static solutions we can construct by varying the thermodynamic parameters of the thermal state close to the critical point. At exactly the critical point, this variation reduces to the static bulk perturbative solution that emerges at the critical point and which leads to the broken JHEP08(2022)246 phase black holes. This, is also the limit of the time dependent Higgs mode we are after when the system is exactly at the transition point. This suggests that the bulk dual of the Higgs mode is parametrically close to the solution generated by varying the thermodynamic parameters and will therefore be the unknown solution in the symplectic current.
The paper is organised in five sections as follows. In section 2 we introduce the holographic framework in which we plan to study phase transitions driven by a scalar order parameter. In this case, the phase transition is not accompanied by the breaking of a continuous symmetry. In section 3 we analyse the time dependent perturbation which captures the Higgs mode and we use our technique involving the symplectic current to extract its gap. In section 4 we generalise our holographic model to describe the breaking of a global symmetry on the boundary at finite chemical potential. This is particularly interest since, as we will see, the lack of particle-hole symmetry mixes the Goldstone mode. In section 5 we compute the Green's functions of the scalar operators in our holographic theory highlighting their domination by the pole relevant to the Higgs mode. Finally, in section 6 we perform numerical checks of our analytic results for the gap and the Green's functions we compute in section 5. We conclude the paper with a summary of our main results and outlook in section 7.

Setup
We will consider a bulk theory which can model the phase transition of a conformal field theory at finite chemical potential and which has been deformed by a scalar operator O φ of dimension ∆ φ with bulk dual φ. In addition to that, our bulk theory will contain the metric g µν , gauge field A µ and scalar field ρ which will spontaneously give a VEV to the field theory dual operator O ρ .
Without loss of generality, we will consider the bulk action, where F = dA is the field strength of A µ . In general, the gauge coupling Z and potential V can depend on the scalar fields φ and ρ. In order for the bulk equations of motion to admit solutions which asymptote to unit radius AdS 4 , we will impose that for small values of the scalar fields, Moreover, without any loss of generality, for small values of ρ we can assume that,

JHEP08(2022)246
The background black holes that can capture the different phases of our theory are described by the background geometries with, The scalar ρ is identically zero in the normal phase, above the critical temperature T c . We will choose coordinates so that the event horizon with Hawking temperature T is at r = 0. Near the horizon, regularity dictates the Taylor expansions, At the other end of the geometry, where r → ∞, we want to impose AdS 4 asymptotics with unit radius, The conformal dimensions ∆ φ and ∆ ρ of the dual operators O φ and O ρ are determined in terms of the masses according to The constants of integration φ (s) and ρ (s) correspond to sources for the boundary operator O φ and the "modulus" of the operator O ρ . We will be mostly setting the latter to zero except for section 5 where will turn it on perturbatively in order to extract a few of the retarded Green's functions of the system. At this point, it is useful to note that the horizon charge density (0) and entropy density s satisfy, where Z (0) denotes the value of the coupling function Z evaluated on the horizon values of the background scalars φ and ρ. For normal fluids, the horizon charge density ρ (0) coincides with the field theory charge density = (v) which can be extracted from the asymptotics (2.6). For the superfluid phase that we will consider in section 4, this is no longer true. In contrast, the entropy density can always be computed as a horizon quantity from (2.7).

The Higgs mode from holography
In this section we will develop our tools that will lead to the main results of our paper. The time dependent perturbation capturing the Higgs mode we are after will be constructed JHEP08(2022)246 in an expansion of the thermodynamic and deformation parameters around the transition point. In section 3.1 we will discuss the static perturbations we can obtain by varying the thermodynamic and deformation parameters of the black hole solutions which are dual to the thermal states of our system. These will prove particularly useful in constructing the next to leading order terms in the expansion of the time dependent perturbation in section 3.2. Finally, in section 3.3 we will make use of the symplectic current to fix the frequency of the quasinormal mode as a function of our thermodynamics and black hole horizon data.

Expansions around the critical point
The gapped mode that we will study in the following sections corresponds to a time dependent perturbation around the background geometry (2.4). We will be interested in studying this perturbation very close to the critical point (T c (µ, φ (s) ), µ, φ (s) ). For this reason, we will introduce a parameter ε which parametrises a curve in the (T, µ, φ (s) ) plane that originates from the point (T c (µ, φ (s) ), µ, φ (s) ). As we vary the parameter ε we move in the space of thermodynamic parameters according to (T c (µ, φ (s) ) + δT (ε), µ + δµ(ε), φ (s) + δφ (s) (ε)). To make our notation clear we consider the curve for small values of the expanision parameter with, For a general function Φ of the thermodynamic parameters, along the curve we can write the variation δΦ = ε 2 2 δΦ (2) + · · · with, Right at the critical point, there exists a static perturbative solution δρ * (0) to the equation of motion for the scalar ρ with source free boundary conditions, just like in (2.6) for the non-linear problem. Our goal is to track this perturbative mode right below the critical temperature where it acquires a gap δω g . As we will see later more explicitly the gap is itself an expandable function of ε and we will have δω g ∝ ε 2 .
In studying our gapped mode, we will find useful to consider the expansion of the background geometry itself in ε. Constructing this perturbative expansion is a rather nontrivial task. However, here we will only need the fact that it exists and that it is analytic in ε,

JHEP08(2022)246
parameter ε * according to, Apart from the broken phase black holes, it is useful to consider the expansion of the normal phase black hole solutions below T c , Similarly to the case of the broken phase, we will define the expansion parameter ε # according to, It is also useful to write thermodynamic quantities in terms of the perturbative expansions of the horizon data. In particular, we can simply expand the charge and the entropy densities of equation (2.7) in ε to obtain, (3.7) The above horizon constants come from e.g. the near horizon expansions δU * (n) (r) = 4π δT * (n) r + · · · , δg * (n) (r) = δg (0) * (n) + · · · , δa * (n) (r) = δa (0) * (n) r + · · · , δφ * (n) = δφ (0) * (n) + · · · , (3.8) and similarly for the normal phase ε expansion. Close to the conformal boundary we will have the asymptotic expansions, (3.9) We will mostly keep δρ (s) * (n) = 0, except in section 5 where we will study the retarded Green's functions of our scalar operators. Note that similar expansions hold for the functions appearing in the expansion around the normal phase black holes.

JHEP08(2022)246
For the VEVs of the scalar operators we can write the expansions, (3.10)

Time dependent perturbation
In this section we will consider perturbative solutions of the equations of motion around the backgrounds (2.4). Since our background geometries possess spacetime translational invariance, it makes sense to consider the Fourier decomposition of our perturbations. Since we are only interested in extracting the gap, we will only consider Fourier modes with zero wavenumber according to, where F repsesents perturbations of the scalars as well as the metric and gauge field components. The function S(r) is chosen so that it drops faster that O(1/r 3 ) close to the conformal boundary and it therefore doesn't interfere with holographic renormalisation. However, close to the horizon, it is chosen so that it approaches S(r) → 1 4πT ln r and the combination t + S(r) is regular and ingoing.
For the system we are interested in, it is consistent to set the component perturbations δg ti , δg ri , δg xy and δa i to zero. Moreover, due to the unbroken rotational symmetry in the x-y plane we can set δg xx = δg yy . Imposing regular ingoing boundary conditions close to the horizon, leads to the near horizon expansions, δg tt (r) = 4πT r δg δg ii (r) = δg (0) + r δg (1) + · · · , δg tr (r) = δg which are compatible with the equations of motion. In order to achieve regularity, the above conditions need to be supplemented by, The time dependent perturbation which is capturing the Higgs mode can be expanded JHEP08(2022)246 (3.14) The above expansion makes clear that our perturbation has to go to the zero mode δρ * (0) at the critical point with ε * = 0. Notice that the choice of δg rr in the ansatz (3.14) fixes the choice of the radial coordinate at order O(ε * ). The functions in equation (3.14) satisfy the boundary conditions (3.12) and (3.13). For ε * = 0 the mode becomes exactly the static mode which drives the transition. This is precisely the way the above expansion is singling out the specific mode we wish to study from the rest of the tower of the quasi-normal modes that we expect to exist in the lower half-plane. As we will see later, the leading term in the frequency expansion ω [1] is equal to zero. An important point in understanding the first non-trivial ε * correction in our perturbation (3.14) is that the equations that the functions δŨ (2) , δg (2) , δã (2) , δφ (2) and δρ (2) satisfy are a superset of the equations saitisfied by the functions δU * (2) , δg * (2) , δa * (2) , δφ * (2) and δρ * (2) . The two extra equations are constraints and impose charge and energy conservation from the field theory point of view. As such these constraints can be imposed on any constant r surface. However, as we will see in section 4, this is no longer true for the equation of motion of the 1-form field due to the fact that it becomes massive in the bulk. This is directly related to the spontaneous symmetry breaking and the involvement of the field theory Goldstone mode.
Returning to the case at hand, apart from the two constraints, this is an inhomogeneous set of ordinary differential equations which are sourced by quadratic terms in δρ * (0) . An important observation to make here is that the homogeneous part of these equations is also solved by the perturbations δU #(2) , δg #(2) , δa #(2) , δφ #(2) defined by normal phase expansions of (3.5) defined by any choice of expansion parameter ε # . For such a solution of the homogeneous part we will have a trivial contribution to δρ (2) .
The final goal is the construction of a linear superposition of perturbations coming from the first corrections of the broken phase branch (3.3) and normal branch (3.5) such that it doesn't carry any net charge or entropy density. More precisely, after setting the value of the perturbative parametres equal ε * = ε # , we have,

JHEP08(2022)246
with the parameters δT #(2) and δµ # (2) such that the total perturbation does not carry any net charge or entropy density. Moreover, from the form of the solution (3.15), we can see that the source for scalar field φ in the perturbation is, Therefore, when we look for quasi-normal modes when trying to compute the gap, we will need to set δφ (s)(2) = 0 in the perturbation. Notice that in order to obtain a source free perturbation for the gauge field, we also performed a gauge transformation with parameter Λ = (−δµ * (2) + δµ #(2) ) (t + S(r)) with S(r) as in (3.11). This step will not be possible in section 4 where we will be dealing with superfluids and the bulk 1-form field will become massive due to the spontaneous symmetry breaking. As we will see, the non-trivial asymptotics of the bulk 1-form field will be rather related to the field theory Goldstone which we expect to be involved in the case of a superfluid at finite chemical potential, with no particle-hole symmetry.
In general we can write variations of the entropy, the charge densities and the scalar VEV in terms of varations of the temperature, the chemical potential and the scalar source according to, where Ξ * and Ξ # are the thermodynamic susceptibility matrices along the broken and normal phases correspondingly evaluated at the critical point. Using the above we can express temperature and chemical potential variations in terms of scalar deformations and entropy and charge densities according to, In the above we have defined the vectors, which according to our earlier discussion, have to be equal to each other when constructing the time dependent perturbation. Finally, using equations (3.17) and (3.18), the variations of the scalar VEVs can be written as,

JHEP08(2022)246
The above allows us to relate the thermodynamic susceptibilities between the two different ensembles through, From the horizon point of view, the zero net charge and entropy density translates to the boundary conditions, 3 This will be useful in the next section where we employ the symplectic current to fix the term ω [2] in the expansion (3.14) and also show that ω [1] = 0. At this point, we understand the gapped mode in the bulk up to order ε * in the perturbative expansion. The ingredients we need are the static solutions of the broken and the normal phase bulk geometries. The final piece of information is the characteristic frequency ω [2] which will be the task of the next subsection.

The symplectic current
In this section, we will employ the Crnkovic-Witten sympectic current to exract the final piece of information we are after, the frequency of the gapped mode close to the phase transition. To define it, we consider a generic classical Lagrangian field theory of a collection ϕ I of fields and two perturbative solutions δ 1 φ I and δ 2 φ I around a background φ I b . If the Lagrangian density L(ϕ I , ∂ µ ϕ I ) can be written in terms of the fields and their first derivatives then the vector density, is divergence free, The symplectic current (3.23) is antisymmetric in field space and as such, when the second perturbative solution is infinitesimally close to the first one, the symplectic form is expandable around zero. As we will see, this observation is going to be particularly useful for us. In particular, the role of the background will be played by the black hole solution in the broken phase. We will take the first of the two perturbations in (3.23) to be the static perturbation which JHEP08 (2022)246 is simply a derivative of the expansion (3.3) with respect to ε * giving, The second perturbative solution will be the time dependent perturbation (3.14) and the expansion parameter will naturally be given by ε * . In order to obtain the symplectic current for the theory of equation (2.1), we first need to write it in a form where all the fields will appear with their first derivative at most. This will yield an equivalent actionS b such that, where we have introduced the Gibbons-Hawking term, and ϕ I include the metric along with the rest of the matter fields of our bulk theory. As usual, K = ∇ µ n µ is the trace of the extrinsic curvature of the conformal boundary ∂M and n = dr/ √ N is its normal one-form of unit norm. The integration measure is with respect to the induced metric h µν = g µν − n µ n ν .
For our bulk gravitational theory originating from (2.1), the sympectic current will read, In appendix A we evaluate the derivatives of the bulk action density with respect to the partial derivatives of our fields. The next step is to evaluate the symplectic current (3.28) for the pair of perturbations (3.25) and (3.14) around the expanded broken phase background (3.3). Before we do that, it is good to understand how we can benefit from considering the condition (3.24) in our setup. Since we are Fourier expanding our time dependent modes, we will write the non-trivial components as, P t = e −iω(t+S(r)) P t (r), P r = e −iω(t+S(r)) P r (r) . (3.29)

JHEP08(2022)246
The continuity equation (3.24) then becomes, Considering now the expansion in ε * , we can write, from where we see that at leading order (3.30) implies, We therefore see that the term involving the time derivative in (3.24), is subleading in the ε * expansion. Keeping up to order ε * in the radial component of the symplectic current we obtain, After integrating (3.32) from the horizon up to infinity, we can show that ω [1] = 0, as promised. This implies that after moving on to higher order in the perturbative expansion we will have that, which will allow us to determine ω [2] , as we will see shortly. Expanding the sympectic current at next to leading order we obtain, Integrating equation (3.34) from the horizon to infinity yields a term on the horizon and a term evaluated at the conformal boundary which have to be equal to each other. The expression we obtain reads,

JHEP08(2022)246
where we imposed the absence of scalar sources on the boundary by setting δφ (s) * (2) = δφ (s)#(2) in equation (3.16). After recognising the variations of the horizon charge and entropy densities (2.7), we can write, − i ω [2] s c 4π δρ with the variations of the scalar VEVs defined according to (3.10). After introducing appropriate factors of ε * , setting, and using equations (3.18) and (3.20), the frequency of our Higgs takes the remarkably simple form, In the above, we have defined the horizon quantity, as well as As we show in appendix B, ∆E is the energy density difference between the normal and broken phase thermal states at entropy density s + δs * , charge density + δ * and scalar deformation parameter φ (s) + δφ (s) * . As we would expect, the frequency of the Higgs mode (3.41) is in the lower half of the complex plane as long as the energy of the broken phase black holes is lower than then normal phase in the microcanonical ensemble. As we will see later, apart from the energy, there will be other quantities which will have to compare between the two phases. For this reason, if a thermodynamic quantity O can be expressed in terms of the entropy density, the charge density and the scalar deformation, ∆O will denote the difference, with O * and O # the values of O in the broken and normal phase respectively. Moreover, s c is the entropy at the critical point and is set by ρ and φ (s) . In our work we will only need the leading order approximation in ε * where as we know δs * ∼ δ * ∼ δφ (s) * ∼ O(ε 2 * ).

JHEP08(2022)246 4 Charged superfluids
In this section, we will examine the case in which the broken phase describes a superfluid. To achieve this, the bulk theory will need to contain a complex scalar field ψ with a U(1) gauged symmetry which will drive the phase transition. For concreteness, we will consider the bulk action, In the superfluid phase, the operator O ψ which is dual to the bulk field ψ takes a nontrivial VEV yielding a corresponding non-trivial bulk field ψ. In order to proceed, we will perform a field transformation to write the complex scalar field in a polar decomposition according to, with ρ > 0 and 0 ≤ θ < 2π/q. The above transformation brings the bulk action to the form, where we have set B = A + ∂θ and therefore F = dB. This form is very similar to the bulk action (2.1) that we considered in the previous sections for neutral scalar fields which are dual to real order parameters. The crucial difference is that in the broken phase, the 1-form field does not enjoy gauge invariance due to the Stueckelberg mechanism.
We would now like to consider the asymptotic expansion of 1-form field B α fluctuations. For the components along the conformal boundary directions we obtain, with the radial component B r completely determined by this. The asymptotic expansion of the 1-form field B is therefore fixed by a scalar function θ (s) and the two 1-forms A α +∂ α θ (v) and j α . The function θ (v) and the parameter ρ (v) that appears in the expansion (2.6) of the modulus ρ of the complex field ψ combine together to parametrise the VEV, from which we obtain the perturbative epxression VEV δθ (v) . In the absence of external sources, the constants of integration j α are directly related to VEV of the boundary charge density J α satisfying the continuity equation ∂ α j α = 0. The treatment of the gapped mode is identical to that of the neutral scalar field with the gauge field A replaced by the massive 1-form B. The whole argument would go through apart from one technical point of particular physical significance. This is related to the gauge transformation that we dicuss below equation (3.15) and which removes the sources from the asymptotic of the gauge field A. Since we cannot perform gauge transformations to B, the same treatment is not possible in this case. However, as we can see from the expansion (4.4), a non zero value for B on the conformal boundary can be treated as the Goldstone mode δθ (v) as long as it is an exact form i.e. it can be expressed as the derivative of the function δθ (v) . This is certainly possible for the perturbations we are studying in this paper since only the time and the radial components of the 1-form B will be involved. The constant on the boundary coming from the time component can then always be expressed as a partial derivative of a phase δθ (v) with respect to time.
To proceed, we introduce the Fourier modes of equation (3.11) according to, Another important ingredient in order to construct the perturbation of order ε * for the 1-form field is the discussion below equation (3.14). We can indeed form zero total charge linear combinations of the static perturbations δa * (2) and δa #(2) as we did in the case of the neutral order parameter in order to obtain solutions of the time component of the equations of motion. The important difference now is that the radial component is no longer a constraint which can be solved everywhere in the bulk by imposing charge continuity on the boundary. This is a signal that, in addition to our previous considerations, the phase δθ in the bulk is going to get involved in the perturbation at order ε. More concretly, we can expand the bulk phase in ε * according to, where δθ (0) is constant everywhere in the bulk and δθ (2) an analytic function of the radial coordinate. The above implies that at order ε * the bulk perturbation for the one form field takes the form, [2] δθ (0) , δb r(2) = −i S ω [2] δθ (0) + δθ (2) . (4.9) A crucial point to note here is that the terms related to δθ (0) and δθ (2) enter the time component of the equation of motion of δB at order ε 3 * . This shows that the above solves the time dependent equation at order ε * even after the addition of the terms which are an exact form.
On the other hand, as one can easily see, both δθ (0) and δθ (2) will enter the radial component of the equation of motion of δB at order ε 3 * . Expanding this equation close to JHEP08(2022)246 the horizon with the boundary conditions given in (2.5), (3.12) and (3.13) we obtain, To obtain the first equality in the second line of (4.10), we used equation (2.7) which relates the charge density perturbation of the normal phase δ #(2) to the flux density of the gauge field at the horizon. At the same time, we have defined δ (0) * (2) , the perturbation of a horizon charge density related to the 1-form field B in the broken phase. In the case of superfluids, an important point is that this is not equal to the charge density of the dual field theory which needs to be read off from the bulk fields after expanding close to the conformal boundary. Finally, to obtain the last equality, we have used that the two perturbations that we obtain from the static backgrounds have to be chosen so that δ #(2) = δ * (2) . We see that the bulk constant δθ (0) is directly proportional to the difference between the horizon and the conformal boundary charge densities.
We now turn our attention to the symplectic current built from the background solution and the time dependent perturbation of our gapped mode. We follow the same logic as in section 3.3 for the neutral case. In fact, the bulk actions (2.1) and (4.3) have identical kinetic terms and they will therefore yield the same form for the symplectic current. Keeping terms up to order ε 2 * yields the same form (3.35) for the radial component. Following the argument of section 3.3, we equate the value of the radial component P r (2) on the horizon to its value at the conformal boundary. In this case however, there is a non-zero contribution from the conformal boundary due to the non-trivial asymptotics of the vector field component δb t (2) . This allows us to write, − i ω [2] e 2g (0) c δρ (0) 2 * (0) − iω [2] (4.11) After using equation (4.10) for the zero order perturbation of the phase δθ (0) , we obtain the simple result, where we have defined the transport coefficient, 13) and ∆E remains the energy difference of equation (3.41). This transport coefficient is intimately connected to the dynamics of the Higgs mode. A certain way to understand it is through the decay rate of equation (4.12). Finally, by using the equation of motion for the bulk one form field, it is easy to show that close to the phase transitions (0) − ∼ δρ 2 ∼ ε 2 giving that ∼ ε 2 .

JHEP08(2022)246 5 Green's functions of scalar operators
In this section we wish to study the effect of the Higgs mode we constructed in the previous section on the linear response of the system near the critical point. For this reason, we will need to introduce boundary sources for the scalar fields in the time dependent perturbation that we constructed in section 3.2. As we will see shortly, the scalar sources will have to be turned on at orders δρ (s) ∼ O(ε 2 * ) and δφ (s) ∼ O(ε * ) in the ε * expansion of equation (3.14). Including sources does not change the arguments around the construction of sections 3.2 and 4. As we have seen there, at next to leading order in ε * , the perturbation can is simply a linear combination of a variation of the broken and the unbroken phase black hole backgrounds. The former variation is fixed by how we move on the phase diagram and does not allow us to introduce new sources in the perturbation.
This shows that the only way to introduce perturbative source for the scalar φ, which are not present in the background thermal state, is through the variation of the normal phase branch of the black holes δφ #(2) in equation (3.15). However, for the scalar which is relevant to the amplitude of the order parameter O ρ , the source will need to be introduced in the O(ε 3 ) term of equation (3.14). This can be achieved through the asymptotics of the bulk function δρ (2) .
We will start our treatment by considering a source δρ (s) for the operator O ρ . As described above, this appears as a constant of integration in the asymptotic expansion, for the function δρ (2) defined through equation (3.14). This proves our earlier statement that the source for O ρ should be of order ε 2 * with δρ (s) = ε 2 * 2 δρ (s) (2) . By equating the value of the radial component of the symplectic current (3.35) on the horizon and the conformal boundary we obtain, − i ω [2] e 2g (0) c δρ (0) 2 * (0) − iω [2] with δθ (0) given once again by equation (4.10). The above allows us to relate the frequency ω = ε 2 * ω [2] to the source δρ (s) . On the other hand, the expansion (3.10) for the VEVs of the operators O ρ and O φ is enough to fix the retarded Green's function, where the symbol ∆ is as defined in equation (3.42). Once again, our ε * of equation (3.14) has singled out the pole which is relevant to the dominant Higgs mode. We expect a whole

JHEP08(2022)246
tower of quasi-normal modes with much larger decay rates and which will not influence the low energy dynamics of our system. In the above, we have defined the transport coefficient and gap frequency ω gap according to, This is the first Green's function we determine and as we anticipated it is dominated by a single simple pole at ω = −i ω gap agreeing with the gap of equation (4.12). Notice, that equation (5.3) allows us to write the expressions, where χ OρOρ is the susceptibility of the operator O ρ . The above allows us to write the Kubo formula, noting that the left hand side remains finite at the critical point, given the fact that ∼ ε 2 . This ratio is very closely related to the coefficient Γ of [4] which the authors assume to remain finite at the critical point. The next Green's function we can compute after having turned on the source for O ρ is, where once again we used the relation (5.2) for the source δρ (s) . In the current situation, we are not turning on the source δφ (s) for the scalar operator. Therefore, in this case we simply have, allowing us to write, Once again, we see that the retarded Green's functions are determined by the transport coefficient and thermodynamic properties of the broken and normal phases. We now turn our attention to the source of the scalar operator O φ . Introducing a source δφ (s)(2) as in equation (3.16) allows us to express the scalar source variation for the normal phase as, δφ (s)# = δφ (s) * − δφ (s) . (5.10)

JHEP08(2022)246
The final ingredient left in order to compute the retarded Green's function is a relation that we can obtain from the symplectic current (3.35). By evaluating it on the horizon and the conformal boundary we obtain the relation, − i ω [2] e 2g (0) c δρ (0) 2 * (0) − iω [2] (5.11) After using equations (3.18) and (3.20) as well as the definitions (3.42) and (5.4), the above takes the form, The combination of this and relation (3.10) allows us to compute the retarded Green's functions, (5.14) From the form of the Green's functions in equations (5.9) and (5.14), it is trivial to check the validity of the Onsager relation G OρO φ (ω) = G * O φ Oρ (−ω). It is interesting to examine the behaviour of the boundary theory Goldstone mode which is represented by the angle δθ (v) of equation (4.6). This can be read off from the asymptotics of the time component of the 1-form field δb t and interpreting the constant term as the time derivative ∂ t δθ (v) . This gives, with the variation of the normal phase chemical potential δµ #(2) being fixed by equation (3.18). The VEV of the operator O ρ is related to the amplitude of the order parameter, according to the decomposition of fluctuations (4.6). More specifically, in the broken phase, where the polar coordinate decomposition (4.2) makes sense, we have,

JHEP08(2022)246
It is interesting to also consider the operator O Y which is related to the phase of the order parameter and therefore the Goldstone mode of the theory. For the fluctuations of this operator we can write, Given the expression (5.15) for the phase, we can easily compute the retarded Green's functions measuring the response of the operator O Y against sources for the operators O ρ and O φ , In the above we have set, An obvious remaining question regards turning on a source δs Y for the operator O Y . As explained in [16], in the spontaneous case this is achieved by considering a non-trivial source term δθ (s) in the UV expansion (4.4) for the bulk vector field. However, as shown in [16], turning a source λ for the complex operator O ψ would imply the non-conservation of electric charge through the Ward identity, Applying this for the case of our perturbative setup gives, Having it mind that the frequency is of order O(ε 2 * ) and the leading correction to the charge density is of order O(ε * ), we see that the source δs Y should be of order O(ε 2 * ) and we will write δs Y = ε 2 * 2 δs Y (2) + · · · . Combining the above ingredients leads us to the charge imbalance, when linearly superposing the static perturbations obtained from the broken and the normal phase black holes. We see that this source implies a slight modification of our previous derivations to account for the non-zero difference between the charge densities. One of the crucial steps concerns equation (4.10) which now gets modified to, 2 ω [2] δs Y (2) .
We are now in the position to evaluate the retarded Green's functions, The last equality can be used to show that our expression for G O Y O Y leads to positive spectral weight. Moreover, the double pole in the same Green's function is due to the presence of the Goldstone mode, directly related to the phase of the condensate. After noting that O Y transforms as a pseudoscalar operator under time reversal, our expressions for the retarded Green's functions satisfy the Onsager relations. 4 Finally, our Green's function G O Y O Y agrees 5 with our previous result in [16] at zero chemical potential and when we take the near critical point and zero pinning limit of the result presented in [16].

Numerical checks
The aim of this section is to provide numerical evidence in support of the analytic expression for the gapped mode given in equation (3.39) and (4.12). We consider the model (2.1) and (4.3) respectively, together with, In order to check this, it is useful to note the Maxwell type of relation coming from the enlarged first law (B.2) 5 In order to see the agreement, one has to write the transport coefficient Ξ of [16] as Ξ ≈ γ χ 2 QQ /2 in the near critical point limit. Moreover, at zero chemical potential we have that ∂µ * ∂ = χ −1 QQ . Note also that there is a mismatch of a factor of 2 due to the different normalisation of the kinetic term of the complex scalar in our bulk action (4.1).

JHEP08(2022)246
and m 2 ρ = −2 , m 2 φ = −2. As described in section 2, we choose to deform our boundary theory by a chemical potential µ = 1 and a relevant operator, O φ , with scaling dimension ∆ φ = 2; in particular we pick the source for φ to be given by φ s = 1. The corresponding backreacted solution will then be given by black holes supported by a non-trivial profile for the scalar field φ, while the scalar field ρ remains trivial. These solutions correspond to the normal phase of the system.

Ungauged model
For the ungauged model with ζ ρ = 1, ζ φ = 0, λ φ = 1, λ ρ = 1, λ = −1, the normal phase black holes exhibit a second order phase transition at temperature T c = 0.115 to a configuration supported by a non-trivial scalar field ρ. The latter describe the broken phase that we are interested in.
With these broken phase black holes at hand, we now turn our attention to studying perturbations around them. In order to check numerically the validity of the analytic expressions for the gap, we need to construct the quasinormal modes. In particular, we consider the following perturbations, where v is the infalling Eddington-Finkelstein coordinate defined as, These perturbations satisfy three first order and two second order differential equations, which we solve using a shooting method subject to appropriate boundary conditions. In particular, in the horizon (located at r = 0) we impose, where c 1 , c 2 , c 3 are fixed, while close to the boundary we impose,

JHEP08(2022)246
For the computation of quasinormal modes, we need to ensure that we remove all the sources from the UV expansion up to a combination of coordinate reparametrisations and gauge transformations, where the gauge transformations are of the form, for ζ, λ constants. This requirement boils down to the sources appearing in (6.5) taking the form, Overall, in addition to the frequency ω, we have two constants in the IR (δρ h , δφ h ) and five constants in the UV (ζ 1 , ζ 2 , λ, δρ (v) , δφ (v) ). This gives a total of eight free constants, out of which one can be set to unity due to the linearity of the equations. The remaining seven parameters match exactly the integration constants of the problem, giving rise to a single (discrete set) solution. In figure 1, we compare the results for gap coming from the direct calculation of the quasinormal modes and from the analytic expression (3.39) for ζ ρ = 1, ζ φ = 0, λ φ = 1, λ ρ = 1, and λ = −1. In particular, we plot the ratio ω analytic /ω numerics as a function of T /T c . We see good quantitative agreement as T → T c .

Superfluid model
For superfluids with ζ ρ = 1, ζ φ = 0, λ φ = 1, λ ρ = 1, λ = −1, the critical temperature is T c = 0.128 and it also corresponds to a second order phase transition. The broken phase is supported by a non-trivial charged condensate, ρ. Given these backreacted black holes, we now turn our attention to studying perturbations around them.
Similarly to above, to check the validity of the analytic expressions for the gap, we construct the quasinormal modes. In particular, we consider the following perturbations, where v is the infalling Eddington-Finkelstein coordinate. These perturbations satisfy three first order and three second order differential equations, which we solve using a shooting method subject to the following boundary conditions. Close to the horizon (located at r = 0) we impose,  Demanding the absence of sources, the sources appearing in (6.11) take the form, Overall, in addition to the frequency ω, we have three constants in the IR (δρ h , δφ h , δb h t ) and six constants in the UV (ζ 1 , ζ 2 , λ, δρ (v) , δφ (v) , δθ (v) ). Due to the linearity of the equations, one out of these ten free constants can be set to unit. One is then left with nine parameters which match exactly the integration constants of the problem, giving rise to a single (discrete set) solution. Just like above, in figure 2, we plot the ratio ω analytic /ω numerics coming from the analytic expression (4.12) and from the direct calculation of the quasinormal modes for ζ ρ = 1, ζ φ = 0, λ φ = 1, λ ρ = 1, and λ = −1. We see good quantitative agreement as T → T c .
For this model with ζ ρ = 0, ζ φ = 0, λ φ = 1, λ ρ = 1/4, λ = −3/4, and T c = 0.0484121, we also compute the Green's functions for the scalar fields for T /T c = 0.999957, corresponding to ω gap = 4 · 10 −5 . In this case, the IR and UV expansions stay unchanged. For the metric and gauge field sources we impose absence of sources just like above. For the scalar field sources, we impose, and we set either (s 1 , s 2 ) = (1, 0) or (s 1 , s 2 ) = (0, 1) depending on which correlator we want to compute. Note that under the reparametrisation (6.6), the VEVs for the scalars transform as Overall, for fixed s 1 , s 2 , we have nine constants in the expansion δρ h , δφ h , δb h t , ζ 1 , ζ 2 , λ, δρ (v) , δφ (v) , δθ (v) in addition to the frequency ω. Given that the problem is fixed in terms of nine integration constants, we conclude that we expect to find a one parameter family of solutions labeled by ω. In figure 3 and 4, we plot the real and imaginary parts of the four Green's functions in terms of the frequency: the red dots correspond to the numerical results, while the solid blue lines indicate the analytic expressions of the previous section. We see excellent quantitative agreement. Note that the agreement is expected to improve as one approaches the critical temperature. The location of the peak in the imaginary part of the Green's functions corresponds to the location of the gap.

Discussion
We considered the perturbative dynamics of holographic systems which are parametrically close to a thermal phase transition. In particular, we focused on gravitational fluctuations capturing the Higgs mode which universally emerges in the broken phase of all second order phase transitions.
One of the crucial ingredients in our construction is the utilisation of the Crnkovic-Witten symplectic current. In the case of the source free quasinormal mode, dual to the Higgs mode, this allowed us to obtain the gap of equations (3.39) and (4.12). In the case of section 5, where the frequency was fixed by an external source, the Crnkovic-Witten symplectic current allowed us to relate the external source to the VEVs of our scalar operators, fixing the form of the Green's functions. The second important ingredient was the construction of the terms which are next to leading order in the ε expansion of equation (3.14). As we explained in section 3.2, these are simply linear combinations of the static perturbations that can be obtained by varying the thermodynamic parameters of the broken and normal phase black holes.
The results of section 5 confirm that the Higgs/amplitude mode dominates the linear response of the scalar operators in the infinite wavelength limit. At finite wavelengths, the Higgs mode will mix with the charge and current densities of theory, affecting the hydrodynamics. It is natural to consider holographic techniques in order to derive an enlarged hydrodynamic description which includes the dynamics of the Higgs mode. This, would increase the range of validity of hydrodynamics close to the critical point.
A different physical situation in which a parametrically small gap exists is when a symmetry is broken pseudospontanesouly. This is the case when apart from the sponta-JHEP08(2022)246 neous breaking, a small explicit parameter that breaks the symmetry is introduced in the system. Our recent work [16] makes clear that our techniques are applicable in the case of superfluids with the internal broken pseudosponaneously. It would be interesting to apply our techniques in setups where the weakly broken symmetry is translations [19][20][21].
In this paper, we have considered phase transitions which preserve the spacetime symmetries of the boundary theory. However, over the past few years a plethora of gravitational instabilities which spontaneously break translations have been discovered in the context of holography [10,11,22,23]. These have been shown to lead to second order phase transitions in general [24][25][26][27][28]. It is interesting to examine the dynamics of the Higgs/amplitude mode for the order parameter of density waves.
From a slightly different angle, in [16] we used the Crnkovic-Witten current to obtain the dissipative terms in the constitutive relations for the electric current in the superfluid phase at zero charge density. Even though that setup is simple due to the non-mixing of the normal and the superfluid components, it can be seen as a proof of principle in extracting the transport coefficients of holographic theories.We expect that the techniques we have developed here and in [16] will be central in the directions outlined above and we wish to report more on these in the near future.
Finally, it is necessary to make the connection between our results and the field theory approach on critical phenomena e.g. [4,29,30]. In particular, it would be very interesting to understand the role of our transport coefficient in the framework of field theory.

JHEP08(2022)246
density is, satisfying the First Law of thermodynamics, given that the gravitation free energy F is the appropriate potential for the grand canonical ensemble, where in the last line we used equation (3.18) to compute the partial derivatives of T and µ at fixed s and . Putting everything together we have that, The above shows that after expanding along the broken and the normal phases and taking the difference, the terms linear in the variations will cancel out and the quadratic ones will give the result of equation (