Higgs response and pair condensation energy in superfluid nuclei

The pairing correlation in nuclei causes a characteristic excitation, known as the pair vibration, which is populated by the pair transfer reactions. Here we introduce a new method of characterizing the pair vibration by employing an analogy to the Higgs mode, which emerges in infinite superconducting/superfluid systems as a collective vibrational mode associated with the amplitude oscillation of the Cooper pair condensate. The idea is formulated by defining a pair-transfer probe, the Higgs operator, and then describing the linear response and the strength function to this probe. We will show that the pair condensation energy in nuclei can be extracted with use of the strength sum and the static polarizability of the Higgs response. In order to demonstrate and validate the method, we perform for Sn isotopes numerical analysis based the quasi-particle random phase approximation to the Skyrme-Hartree-Fock-Bogoliubov model. We discuss a possibility to apply this new scheme to pair transfer experiment.

with respect to the phase and the amplitude, respectively. The Anderson-Bogoliubov mode is a mass-less phonon with a linear dispersion relation whereas the Higgs mode is massive, i.e. the lowest excitation energy is 2∆ [10,[13][14][15][16].
An analogy to the superfluidity/superconductivity has been brought into finite nuclei by applying the BCS theory to the nuclear Hamiltonian [4,5,17,18]. In this description, nuclei may be in superconducting/superfluid or normal phase depending on whether they are closed or open shell nuclei, and also on other conditions such as rotational frequency [1,2].
The odd-even mass difference is a typical indicator of the nuclear superfluid phase [4,5] since it corresponds to the pair gap in the single-particle excitation spectrum. Two-nucleon transfer (pair transfer) reaction has been considered as a probe for the pair correlation in nuclei [1,4,[19][20][21][22][23]. A typical example is the (p, t) and (t, p) reaction on the even-even Sn isotopes, which shows enhanced cross sections for transitions between 0 + ground states of adjacent isotopes [1,[19][20][21]. In parallel to the rotational motion of deformed nuclei, the concept of the pair rotation is introduced [1,20,21], as a nuclear counter part of the Anderson-Bogoliubov mode [10,11]. It has been argued that the pair rotation is a phase rotation of the gap, and is characterized with a rotational band, a series of 0 + ground states along the isotope chain. The transition strength of the pair transfer within the rotational band is related to the order parameter, the pair gap.
Another collective mode associated with the pair correlation is the pair vibration [1,20,21], which is introduced in an analogy to the shape vibration modes in nuclei. It is a collective vibrational state with spin parity 0 + and excitation energy around ∼ 2∆. It is regarded as oscillation of the amplitude of the pair gap ∆. The lowest-lying excited 0 + states populated by the (p, t) and (t, p) reactions or other two-neutron transfer are identified as the pair vibrational states. Although not stated explicitly in the preceding works, the low-lying pair vibrational mode might be regarded as a counter part of the Higgs mode in superfluids/superconductors.
We remark here that microscopic theories such as the quasiparticle random phase approximation predict also another pair vibrational mode with high excitation energy, called often the giant pair vibration (GPV) [24][25][26][27], which arises from the shell structure inherent to finite nuclei. It is therefore a non-trivial issue how the pair vibrations in nuclei, either low-lying or high-lying, or both, are related to the macroscopic picture of the Higgs mode. Recently an observation of the giant pair vibrations in a two-neutron transfer reaction populating 14 C and 15 C has been reported [25,27]. Identification in medium and heavy nuclei is an open question [26] and further experimental study is awaited.
In the present study, we intend to describe the pair vibrations, including both low-lying and high-lying ones, from a view point of the Higgs mode. We shall discuss then that this viewpoint provides a new perspective to the nature of the pair correlation in nuclei.
Our study proceeds in two steps. Firstly, we introduce a new scheme of exploring a response that probes the Higgs mode in a pair correlated nucleus. One of commonly adopted methods to describe the pairing collective modes, including the pair vibrations, is to consider addition or removal of a two-nucleon pair, or describe transition strength for pair-addition and pair-removal operators. In the present work, we propose a variant of these operators, which is more suitable for probing the Higgs mode in finite nuclei. We formulate then the strength function for this new two-nucleon transfer operator, named Higgs strength function, and we then characterize the pair vibrations, including both low-lying and high-lying, in terms of the Higgs strength function. The formulation is based on Skyrme-Hartree-Fock-Bogoliubov mean-field model and the quasiparticle random phase approximation, which are widely used to describe the ground state and collective excitations, including pair vibrations, in medium and heavy nuclei [5,28,29].
Secondly, we shall argue that the Higgs strength function carries an important information on the pair correlation in nuclei, in particular, the pair condensation energy or the effective potential energy, which is the energy gain obtained by the condensation of Cooper pairs.
Similarly to the Ginzburg-Landau theory of superconducting systems, an effective potential energy as a function of the order parameter, i.e, the pair condensate, can be considered also in finite nuclei, and it is straightforwardly evaluated in the framework of the Hartree-Fock-Bogoliubov mean-field model. On another hand, the Higgs strength function probes the small amplitude oscillation of the order parameter and one can expect that it can be related to the curvature of the effective potential. Motivated by this idea, we examine in detail the effective potential energy obtained with the present Skyrme-Hartree-Fock-Bogoliubov model. As shown later, the effective potential has a rather simple structure, parametrized with a fourth order polynomials of the pair condensate. Combining these observations, we propose a novel method with which one can evaluate the pair condensation energy using the Higgs strength function. We shall discuss the validity of this method with numerical calculations performed for neutron pairing in Sn isotopes.

II. THEORETICAL FRAMEWORK
As a theoretical framework for our discussion, we adopt the Hartree-Fock-Bogoliubov (HFB) theory in which the pair correlation is described by means of the Bogoliubov's quasiparticles and the selfconsistent pair field [5]. It has the same structure as that of the Kohn-Sham density functional theory for superconducting electron systems with an extension based on the Bogoliubov-de Genne scheme [28]. The adopted model is essentially the same as the one employed in our previous studies [30][31][32][33]. We here describe an outline of the model with emphasis on some aspects which are necessary in the following discussion.
In the HFB framework, the superfluid/superconducting state is described with a generalized Slater determinant |Ψ , which is a determinantal state of Bogoliubov's quasiparticles. The total energy of the system for |Ψ is a functional of one-body densities of various types, and we adopt the Skyrme functional E Skyrme combined with the pairing energy E pair . Assuming that the pairing energy originates from a contact two-body interaction, E pair [ρ(r),ρ(r),ρ * (r)] is an functional of the local pair condensate (also called pair density in the literature) and local one-body density where ψ † (rσ) and ψ(rσ) are nucleon creation and annihilation operators with the coordinate and the spin variables. Note that the local pair condensateρ(r) andρ * (r) is finite if the condensation of Cooper pairs occurs. Here and hereafter we do not write the isospin index for simplicity. Notation follows Ref. [34].
The one-body density and the pair density are evaluated as a sum of the quasiparticle wave functions: The HFB ground state |Ψ 0 is obtained by solving the HFB equation with an iterative procedure.

U(1) gauge symmetry and its violation
All the states e iθN |Ψ 0 transformed from one realization of the HFB ground state |Ψ 0 are also degenerate HFB ground state. This can be seen in the fact that the above equations are invariant for ϕ 1,i (rσ) → e iθ ϕ 1,i (rσ), ϕ 2,i (rσ) → e −iθ ϕ 2,i (rσ) together with Eq.(7). In the following we choose one of the HFB ground states in which the pair condensateρ(r) is real.
In fact, the pair density and the pair potential are assumed to be real in our numerical code for the HFB ground state.
The situation is illustrated schematically in Fig.1. Here the energy gain or the effective potential U (p) = E(p) − E(0) is shown as a function of an "order parameter" p, which could be an average value of the pair condensateρ(r) or the pair potential ∆(r) (Details will be specified later). The order parameter p is a complex variable and it transforms as p → e 2iθ p under the U(1) gauge transformation. The origin p = 0 corresponds to the Hartree-Fock ground state where the pair condensate is absent, and p = p 0 corresponds to the HFB ground state. The difference U cond = E(p 0 ) − E(0) is the energy gain associated with the condensation of the Cooper pairs.

Excitation modes
We describe excitation modes build on top of the HFB ground state |Ψ 0 by means of the quasi-particle random-phase approximation (QRPA). The QRPA is equivalent to describing a linear response under an external perturbation in the frame work of the time-dependent HFB theory [28,34].
We consider a time-dependent perturbation µe −iωtV ext where a perturbing fieldV ext belongs to a class of generalized one-body operators expressed in terms of the local densityρ(r) and the pair densitiesP (r) andP † (r): The perturbation causes time-dependent fluctuations δρ(r, ω), δρ(r, ω) and δρ (r, ω) (expressed in the frequency domain) in the three densities ρ(r),ρ(r) andρ * (r). It induces also fluctuations in the self-consistent field δV sc (ω) = dr δΓ(r, ω)ρ(r) + δ∆ * (r, ω)P (r) + δ∆(r, ω)P † (r) (9) through the densities. Here we assume that fluctuations in the Hartree-Fock potential and the pair potential, δΓ and δ∆, are proportional to the density fluctuations δρ(r, ω), δρ(r, ω) and δρ (r, ω). Then the linear response of the system is governed by Here R αβ 0 (r, r , ω) is unperturbed density response function with E i in the denominator being the excitation energy of the quasiparticle state i. The numerators in the r.h.s. are matrix elements of the density operatorsρ α (r) =ρ(r),P (r),P † (r) between the HFB ground state |0 = |Ψ 0 and two-quasiparticle configurations |ij = Here we put the spectral representation of the unperturbed density response function expressed with discretized spectrum even for unbound quasiparticle states. The continuum nature of the unbound quasiparticle states can be taken into account explicitly using the method of the continuum QRPA formalism [34]. See Appendix for the two-quasiparticle matrix elements and the continuum QRPA.
The excited states |ν of the system appear as poles of the density fluctuations as a function of the excitation energyhω. The transition matrix elements for the perturbing field is obtained through the strength function which can be calculated in terms of the solution of the linear response equation:

Model parameters and numerical details
We discuss the pairing correlation of neutrons in semi-magic Sn isotopes, for which the ground states are expected to be spherical. We adopt the SLy4 parameter set [35] for the Skyrme energy functional. The pairing interaction defining the pairing energy functional is the density dependent delta interaction (DDDI), the contact force whose interaction strength depends on the position through the nucleon density: Correspondingly the pair potential is given by ∆ q (r) = V q (r)ρ q (r). In the actual calculation, we assume the spherical HFB mean-field and solve all the equations using the polar coordinate system and the partial wave expansion. The radial coordinate is truncated at r = R max with R max = 20fm. The partial waves are taken into account up to the maximum angular quantum number l max , j max = 12, 25/2. We truncate the quasiparticle states by maximum quasiparticle energy E max = 60 MeV. The radial coordinate is discretized with an interval ∆r = 0.2fm. The

III. PAIR VIBRATIONS AND HIGGS RESPONSE
A. Response to pair transfer operators Sensitive probes to the pairing correlation may be expressed in terms of the pair field We define a pair-addition and a pair-removal operator These operators bring about a transition changing particle number by ∆N = ±2. Here we introduce a form factor f (r) which is effective in a spatial region where the nucleon density is finite, but diminishes far outside the nucleus as we are interested in a process where a nucleon pair is added to or removed from a nucleus. Specifically we choose a Woods-Saxon function, with R = 1.27 × A 1/3 fm and a = 0.67 fm, but as shown later the main conclusion does not depend on detailed form of f (r). In the present study we describe the pair vibration with the lowest multipolarity J π = 0 + . Here Y 00 is the spherical harmonics with rank 0 and bothP ad andP rm carry the angular and parity quantum numbers 0 + . We describe the pair vibration with spin parity J π = 0 + .
Response of a nucleus against these operators are represented by the strength functions for the pair-addition process, and for the pair-removal process. Here |ν is the QRPA excited states with excitation energy E ν whereas ν|P ad |0 and ν|P rm |0 is the matrix elements of the pair-add and pair-removal operators between the HFB ground state |0 = |Ψ 0 and the QRPA excited states. Transition densities representing amplitudes of pair-addition and -removal are also useful measures to look into structure of the excited states. It is seen that the pair-addition transition density and the pair-removal one has the same shape as functions of r, but they have the opposite phases. This is what is expected for the pair rotation, since the phase change in the pair condensateρ(r) gives δρ(r) = δe 2iδθρ (r) ∼ 2iδθρ(r) and δρ * (r) = δe −2iδθρ (r) ∼ −2iδθρ(r), i.e., opposite sign in these quantities with a common shape proportional to the ground state pair condensateρ(r).
Other peaks in Fig.2

B. Higgs response
We shall here address a question how the pair vibrations discussed above are related to the macroscopic picture as an amplitude oscillation of the order parameter, the Higgs mode.
However, the existence of multiple states, including both the low-lying and high-lying pair vibrations, makes it difficult to relate this macroscopic picture to one of the pair vibrational states. We need a new way of characterizing the pair vibrations, which is more suitable for obtaining the macroscopic information.
As is illustrated by Fig.1, response of  has two different directions with respect to the order parameter of the pair condensation: the phase mode changing the phase of the pair gap (or the pair condensate), and the other changing the amplitude of the pair gap/condensate. The pair-add operators,P ad ∼ ψ † ψ † , and the pair-removal operatorsP rm ∼ ψψ, do not probe separately these two different modes since both the pair rotation (the phase mode) and the pair vibrations (possible amplitude mode) are seen in both of the pair-addition and pair-removal strengths functions.
Instead we introduce a pair field operator that is a linear combination ofP ad andP rm : Note that fluctuation of the amplitude |ρ(r)| of the pair condensate reads δ|ρ(r)| = (ρ(r)δρ * (r, t) +ρ * (r)δρ(r, t)) /2|ρ(r)| = (δρ * (r, t) + δρ(r, t))/2 for realρ(r), and therefore the Higgs operatorP H defined by Eq.(23) probes the amplitude fluctuation. In the following, we callP H the Higgs operator after the nomenclature that the amplitude mode of the pair field is often called Higgs mode [15,16]. The strength function for the Higgs operator is defined by Secondly we observe that the Higgs strength function is close to the sum of the pairaddition and pair-removal strengths although there are some differences. A significant and important difference is seen for the pair rotation mode (the peak at E = 0.60 MeV in this example), for which the Higgs strength almost invisible. Indeed, it should vanish in principle due to the very nature of the phase mode associated with the U(1) gauge symmetry. As discussed above, the pair-add and pair-removal amplitudes of the phase mode have the same shape but with opposite sign (cf. Fig.3), and consequently these two amplitudes cancel each other in the matrix element of the Higgs operator ν|P H |0 = ν|P ad |0 + ν|P rm |0 = 0.
For the low-lying pair vibration (the peak at E = 2.44 MeV), contrastingly, the Higgs strength is enhanced. Namely it is more than the sum of the pair-addition and pair-removal strengths (It is approximately two times the sum). This is because the both pair-addition and pair-removal transition densities of this mode have the same sign, as seen in where P H is a change of the expectation value under the static perturbationV ext = µP H .
The static polarizability is related to the strength function through the inversely energy weighted sum is obtained. We note that the inversely energy-weighted sum and the static polarizability has been discussed intensively for the case of the E1 strength function, i.e. the electric dipole excitations in nuclei, as a probe of neutron skin and nuclear matter properties [36][37][38][39][40].
In the present study the inversely weighted sum is calculated by taking an integral in an energy interval from E min = 1.1 MeV to the upper limit of the excitation energy E max = 60 MeV. The lower bound is set to exclude the contribution from the pair rotation, which should have vanished if the complete selfconsistency is fulfilled in the numerical calculation.  As the order parameter we adopt the expectation value of the pair removal operator which also represents a spatial integral (or average) of the pair condensateρ(r). (The factor of two is put here just for later convenience). The order parameter p is a complex variable as it follows p → e 2iθ p under the U(1) gauge transformation. We then consider the total energy calculation is plotted with the diamond symbol. The dashed curve is a quartic function, Eq.(29), fitted to the CHFB results. Figure 6 shows the effective potential calculated for 120 Sn, 132 Sn and 140 Sn. Let us focus on the representative example of 120 Sn. The effective potential in this example has a so-called Mexican hat shape, which has a minimum at finite p = p 0 where p 0 = P H corresponds to the HFB ground state |Ψ 0 . Since the shape is smooth as a function of p, it may be approximated by a polynomial. Guided by the Ginzburg-Landau theory of the superconductivity [9], we here assume a fourth-order polynomial and we fit this function to the calculated effective potential E(p). As seen in Fig.6, the fourth-order polynomial is a good representation of the CHFB results. We have checked and confirmed good accuracy of Eq.(29) for other even-even Sn isotopes, see two other examples in Fig.6.

B. Pair condensation energy
In the case where the effective potential is approximated well by the fourth-order polynomial there holds a simple relation between the potential parameters, a and b, and the order parameter p 0 and the curvature C of the effective potential at the potential minimum: As a consequence, the potential depth D can be also related to the two parameters characterizing the minimum: We remark here that the parameters p 0 and C can be derived from the Higgs response.
The order parameter p 0 at the minimum is the expectation value of the Higgs operator for the ground state The potential curvature C at the minimum p = p 0 is related to the Higgs polarizability α H , Eq. (25), through the Hellmann-Feynman theorem [5,41]. Consequently the potential depth D, which can be interpreted as the pair condensation energy U cond = E(p 0 ) − E(0), can be evaluated as is U Higgs cond = 1.51 MeV. It is compared the pair condensation energy obtained directly from the CHFB calculation, U CHFB cond = 1.58 MeV. It is found that the evaluation using the Higgs response reproduces the actual value ( obtained directly from the CHFB calculation) within an error less than 10 %. It suggests that the pair condensation energy may be evaluated via the Higgs response with good accuracy.

V. SN ISOTOPES
We shall apply the above method of evaluating the pair condensation energy to the neutron pair correlation in even-even Sn isotopes covering from doubly-magic neutron-deficient 100 Sn to neutron-rich 150 Sn. The ground state value of the order parameter p 0 = P H as well as the average pairing gap ∆ = drρ(r)∆(r)/ drρ(r) are shown in Fig.7 (a). It is seen that these two quantities are essentially proportional to each others. They are zero in neutron magic nuclei 100 Sn and 132 Sn, where the pair condensate vanishes due to a large shell gap at the neutron Fermi energy, corresponding to the "normal" phase. closed-shell nuclei 100 Sn and 132 Sn is larger than neighboring isotopes. It is also remarkable that the largest value is in 140 Sn where the pairing gap takes an intermediate value. Figure 8 is the pair condensation energy U Higgs cond of neutrons evaluated from the Higgs response, i.e. using Eq. (35). It is compared with the pair condensation energy U CHFB cond which is calculated directly using the constrained Hartree-Fock-Bogoliubov method for several isotopes. We confirm that the proposed method works well not only for 120 Sn discussed above, but also for the other isotopes.
Another point is that the isotopic trend of the pair condensation energy shows some The dotted line is an analytic estimate 1 2 g∆ 2 with the oscillator estimate of the single-particle level density g = 0.015A MeV −1 .
resemblance and difference from that of the pair gap. For comparison we plot also an estimate of the pair condensation energy 1 2 g∆ 2 based on the equidistant single-particle model [5] where g is the single-particle level density. This estimate reflects the pair correlation through the pair gap alone. It is seen that the microscopic condensation energy differs from this estimate with respect to both the magnitude and the isotopic dependence. This points to that the pair condensation energy is a physical quantity which carries information not necessarily expected from the pairing gap.
A noticeable example is 140 Sn, where the pair condensation energy is very small (about 0.1 MeV). The small condensation energy is due to the large value of the Higgs polarizability (see curvature ( Fig.6(a)), manifests itself as the presence of the low-lying pair vibration which has a very large Higgs strength, several times larger than that in 120 Sn.
We also note 132 Sn, in which the HFB ground state has no pair condensate. Correspondingly the effective potential shown in Fig. 6 (b) has the minimum at p 0 = 0, and there is no pair rotation in this case. The Higgs polarizability (or the potential curvature at the minimum) in this case is slightly larger (smaller) than that for 120 Sn (see Fig.7 (b) and Fig.6 ). This feature can be traced to the presence of the pair vibrations with relatively large pair-transfer strengths as seen in Fig.9.

VI. DISCUSSION
The method to evaluate the pair condensation energy from the Higgs response should not depend on details of the definitions of the pair transfer operatorsP ad ,P rm andP H . We have checked other choices of the form factor f (r) appearing in Eqs. (16) and (17)  We envisage application of the present method to experimental studies using two-neutron transfer reactions, such as (p, t) and (t, p) reactions. As discussed in Ref. [1,19,21], the cross sections of a two-neutron stripping and pick-up reactions may be related to matrix elements of the pair-addition and the pair-removal operators if one assumes a one-step process. In this simple picture, the transition amplitude for the ground-to-ground transition N → N ± 2 may correspond to matrix elements 0 + gs , N + 2|P ad |0 + gs , N and 0 + gs , N − 2|P rm |0 + gs , N , which is approximated by the expectation value Ψ 0 |P ad |Ψ 0 = Ψ 0 |P rm |Ψ 0 = P H /2 for the HFB ground state |Ψ 0 . Thus the ground-to-ground pair transfer is expected to provide the order parameter p 0 = P H . Similarly the two-neutron transfer cross sections populating excited 0 + states may be related to matrix elements 0 + ν , N + 2|P ad |0 + gs , N and 0 + ν , N − 2|P rm |0 + gs , N , which correspond to the matrix elements ν|P ad |0 and ν|P rm |0 for the QRPA excited states |ν . We thus expect that two-neutron stripping and pick-up reactions populating excited 0 + states provides the pair-addition and pair-removal strength functions, S ad (E) and S rm (E). The Higgs strength function S H (E) is then evaluated as a sum of S ad (E) and S rm (E) (cf. Fig.4). Here transition to the low-lying pair vibration needs to be treated separately since for this case the Higgs transition amplitude is a coherent sum of the amplitudes ν|P ad |0 and ν|P rm |0 . Once the Higgs strength function is obtained in this way, the Higgs polarizability α H can be evaluated. Consequently, combining p 0 and α H thus obtained, we may estimate the neutron pair condensation energy U Higgs cond using Eq. (35). Note that this argument based on the one-step picture may be simplistic from a viewpoint of the reaction mechanisms, such as the finite range effect and the two-step processes [23]. More quantitative analysis of the two-nucleon transfer reactions needs to be explored in detail. However it is beyond the scope of the present work, and we leave it for a future study.

VII. CONCLUSIONS
We have discussed a new idea of the Higgs response, which probes the Cooper pair con-
The continuous nature of the unbound quasiparticle states can be taken into account in the unperturbed response function by using the Green's function for the quasiparticle states: where G 0 (r σ rσ, E) is the Green's function for the HFB equation (4), and C dE is a contor integral in the complex E plane. Details are given in Ref. [34].