text:
# Analytic Approaches to Kaon Physics

## 1 INTRODUCTION

In the Standard Model, the electroweak interactions of hadrons at very low energies are conveniently described by an effective chiral Lagrangian, which has as active degrees of freedom the low lying $`SU(3)`$ octet of pseudoscalar particles, plus leptons and photons. The underlying theory is a $`SU(3)_C\times SU(2)_L\times U(1)_W`$ gauge theory which is formulated in terms of quarks, gluons and leptons, together with the massive gauge fields of the electroweak interactions and the hitherto unobserved Higgs particle. Going from the underlying Lagrangian to the effective chiral Lagrangian is a typical renormalization group problem. It has been possible to integrate the heavy degrees of freedom of the underlying theory, in the presence of the strong interactions, perturbatively, thanks to the asymptotic freedom property of the $`SU(3)`$–QCD sector of the theory. This brings us down to an effective field theory which consists of the QCD Lagrangian with the $`u`$, $`d`$, $`s`$ quarks still active, plus a string of four quark operators and mixed quark–lepton operators, modulated by coefficients which are functions of the masses of the fields which have been integrated out and the scale $`\mu `$ of whatever renormalization scheme has been used to carry out this integration. We are still left with the evolution from this effective field theory, appropriate at intermediate scales of the order of a few GeV, down to an effective Lagrangian description in terms of the low–lying pseudoscalar particles which are the Goldstone modes associated with the spontaneous symmetry breaking of chiral–$`SU(3)`$ in the light quark sector. The dynamical description of this evolution involves genuinely non–perturbative phenomena and it is mostly studied using the techniques of Lattice QCD. In this talk, I shall review the progress which has been made in approaching this last step using analytic methods; in particular when the problem is formulated within the context of QCD in the limit of a large number of colours $`N_c`$.

The suggestion to keep $`N_c`$ in QCD as a free parameter was first made by G. ’t Hooft  as a possible way to approach the study of non–perturbative aspects of QCD. The limit $`N_c\mathrm{}`$ is taken with the product $`\alpha _\text{s}N_c`$ kept fixed. In spite of the efforts of many illustrious theorists who have worked on the subject, QCD in the large–$`N_c`$ limit ($`\text{QCD}_{\mathrm{}}`$) still remains unsolved; but many interesting properties have been proved, which suggest that, indeed, the theory in this limit has the bulk of the non–perturbative behaviour of full QCD. In particular, it has been shown that, if confinement persists in this limit, there is spontaneous chiral symmetry breaking . It should be stressed that $`\text{QCD}_{\mathrm{}}`$ is not just a “wild extrapolation” from $`N_c=3`$ to $`N_c=\mathrm{}`$. In fact, $`N_c`$ is really used as a label to select specific topologies among Feynman diagrams. The topology which corresponds to the highest power in the $`N_c`$–label is the one which selects planar diagrams only; and the claim is that this class provides already a good approximation to the full theory.

The spectrum of $`\text{QCD}_{\mathrm{}}`$ consists of an infinite number of narrow stable meson states which are flavour nonets . This spectrum looks a priori rather different to the one in the real world: the examples of the vector and axial–vector spectral functions measured in $`e^+e^{}`$ hadrons and in the hadronic $`\tau `$–decay which are shown below, have indeed a richer structure than just a sum of narrow states. There are, however, many instances where one can show that the observables one is interested in, are given by smooth integrals of specific hadronic spectral functions, though unfortunately in most cases, the specific hadronic spectral functions in question are not accessible from experiments. Typical examples of such observables, are the coupling constants of the effective chiral Lagrangian of QCD at low energies, as well as the coupling constants of the effective chiral Lagrangian of the electroweak interactions of pseudoscalar particles in the Standard Model. What is needed for the evaluation of these coupling constants is not so much the detailed point–by–point knowledge of the relevant hadronic spectrum, but rather a good approximation consistent with the asymptotic properties of QCD, both at short– and long–distances. It is in this sense that the simple $`\text{QCD}_{\mathrm{}}`$–spectrum becomes useful. It provides a simple parameterization of the physical hadronic spectrum, based on first principles.

## 2 THE CHIRAL LAGRANGIAN

The strong and electroweak interactions of the Goldstone modes at very low energies are described by an effective Lagrangian which has terms with an increasing number of derivatives (and quark masses if explicit chiral symmetry breaking is taken into account.) Typical terms of the chiral Lagrangian are

$`_{\text{eff}}`$ $`=`$ $`\underset{\pi \pi \pi \pi ,K\pi e\nu }{\underset{}{{\displaystyle \frac{1}{4}}F_0^2\text{tr}\left(D_\mu UD^\mu U^{}\right)}}`$ (3)
$`+\underset{\pi e\nu \gamma }{\underset{}{L_{10}\text{tr}\left(U^{}F_{R\mu \nu }UF_L^{\mu \nu }\right)}}+\mathrm{}`$
$`+\underset{e^2C\frac{2}{F_0^2}\left(\pi ^+\pi ^{}+K^+K^{}\right)}{\underset{}{e^2C\text{tr}\left(Q_RUQ_LU^{}\right)}}+\mathrm{}`$
$$\underset{K\pi \pi ,K\pi \pi \pi }{\underset{}{\frac{G_\text{F}}{\sqrt{2}}V_{ud}V_{us}^{}g_{\underset{¯}{8}}F_0^4\left(D_\mu UD^\mu U^{}\right)_{23}}}+\mathrm{},$$
(4)
where $`U`$ is a $`3\times 3`$ unitary matrix in flavour space which collects the Goldstone fields and which under chiral rotations transforms as $`UV_RUV_L^{}`$; $`D_\mu U`$ denotes the covariant derivative in the presence of external vector and axial–vector sources. The first line is the lowest order term in the sector of the strong interactions , $`F_0`$ is the pion–decay coupling constant in the chiral limit where the light quark masses $`u`$, $`d`$, $`s`$ are neglected ($`F_090\text{MeV}`$); the second line shows one of the couplings at $`𝒪(p^4)`$ ; the third line shows the lowest order term which appears when photons and $`Z^{}s`$ are integrated out ($`Q_L=Q_R=\text{diag.}[(2/3,1/3,1/3]`$), in the presence of the strong interactions ; the fourth line shows one of the lowest order terms in the sector of the weak interactions. The typical physical processes to which each term contributes are indicated under the braces. Each term is modulated by a coupling constant: $`F_0^2`$, $`L_{10}`$,… $`C`$$`g_{\underset{¯}{8}}`$…, which encodes the underlying dynamics responsible for the appearance of the corresponding effective term. The evaluation of these couplings from the underlying theory is the question we are interested in. The coupling $`g_{\underset{¯}{8}}`$ for example, governs the strength of the dominant $`\mathrm{\Delta }I=1/2`$ transitions for $`K`$–decays to leading order in the chiral expansion.

### 2.1 Two crucial observations

There are two crucial observations to be made concerning the relation of these low energy constants to the underlying theory.

i) The low–energy constants of the Strong Lagrangian, like $`F_0^2`$ and $`L_{10}`$, are the coefficients of the Taylor expansion of appropriate QCD Green’s Functions. For example, with $`\mathrm{\Pi }_{LR}(Q^2)`$ the correlation function of a left–current with a right–current in the chiral limit:

$`{\displaystyle d^4xe^{iqx}0|T\left(\overline{u}_L\gamma ^\mu d_L(x)\overline{u}_R\gamma ^\nu d_R(0)^{}\right)|0}`$
$`={\displaystyle \frac{1}{2i}}(q^\mu q^\nu g^{\mu \nu }q^2)\mathrm{\Pi }_{LR}(Q^2),`$ (5)
the Taylor expansion

$$Q^2\mathrm{\Pi }_{LR}(Q^2)|_{Q^20}=F_0^24L_{10}Q^2+\mathrm{},$$
(6)
defines the constants $`F_0^2`$ and $`L_{10}`$.

ii) By contrast, the low–energy constants of the Electroweak Lagrangian, like e.g. $`C`$ and $`g_8`$, are integrals of appropriate QCD Green’s Functions. For example ,

$$C=\frac{3}{32\pi ^2}_0^{\mathrm{}}𝑑Q^2\frac{M_Z^2}{Q^2+M_Z^2}\left(Q^2\mathrm{\Pi }_{LR}(Q^2)\right).$$
(7)
Their evaluation appears to be, a priori, quite a formidable task because they require the knowledge of Green’s functions at all values of the euclidean momenta; i.e. they require a precise matching of the short–distance and the long–distance contributions to the underlying Green’s functions.

These two observations are generic in the case of the Standard Model, independently of the $`1/N_c`$–expansion. The large–$`N_c`$ approximation helps, however, because it restricts the analytic structure of the Green’s functions in general, and $`\mathrm{\Pi }_{LR}(Q^2)`$ in particular, to be meromorphic functions: they only have poles as singularities; e.g., in $`\text{QCD}_{\mathrm{}}`$,

$$\mathrm{\Pi }_{LR}(Q^2)=\underset{V}{}\frac{f_V^2M_V^2}{Q^2+M_V^2}\underset{A}{}\frac{f_A^2M_A^2}{Q^2+M_A^2}\frac{F_0^2}{Q^2},$$
(8)
where the sums are extended to an infinite number of states. In practice, however, in the case of Green’s functions which are order parameters of spontaneous chiral symmetry breaking (like $`\mathrm{\Pi }_{LR}`$ in the chiral limit), these sums will be restricted to a finite number of states.

### 2.2 Sum Rules

There are two types of important restrictions on Green’s functions like $`\mathrm{\Pi }_{LR}(Q^2)`$. One type follows from the fact that, as already stated above, the Taylor expansion at low euclidean momenta must match the low energy constants of the strong chiral Lagrangian. This results in a series of long–distance sum rules like e.g.

$$\underset{V}{}f_V^2\underset{A}{}f_A^2=4L_{10}.$$
(9)
Another type of constraints follows from the short–distance properties of the underlying Green’s functions. The behaviour at large euclidean momenta of the Green’s functions which govern the low energy constants of the chiral Lagrangian can be obtained from the operator product expansion (OPE) of local currents in QCD. In the large–$`N_c`$ limit, this results in a series of algebraic sum rules  which restrict the coupling constants and masses of the hadronic poles. In the case of the $`LR`$–correlation function in Eq. (8) one has e.g.,

$$f_V^2M_V^2f_A^2M_A^2F_0^2=0,$$
(10)
$$f_V^2M_V^4f_A^2M_A^4=0,$$
(11)
$$f_V^2M_V^6f_A^2M_A^64\pi \alpha _s\overline{\psi }\psi ^2.$$
(12)
The sum rules in Eqs. (10) and (11) are the celebrated Weinberg sum rules, which follow from the fact that in the chiral limit, there are no $`𝒪(1/Q^2)`$ terms and no $`𝒪(1/Q^4)`$ terms in the OPE of $`\mathrm{\Pi }_{LR}(Q^2)`$. The third sum rule in Eq. (12, follows from the matching between the $`𝒪(1/Q^6)`$–terms in the $`\text{QCD}_{\mathrm{}}`$ expression in Eq. (8) and the corresponding one in the OPE (evaluated at leading $`𝒪(\alpha _\text{s}N_c^2)`$ ). In principle there are an infinite number of sum rules in $`\text{QCD}_{\mathrm{}}`$ which relate the masses and couplings of the $`\text{QCD}_{\mathrm{}}`$ narrow states to the local order parameters which appear in the OPE and to the non–local order parameters which govern the chiral expansion in the strong sector.

## 3 THE $`\pi ^+\pi ^0`$ MASS DIFFERENCE AS A THEORETICAL LABORATORY

The physical effect of the coupling $`C`$ in Eq. (3) is to give a mass, mostly of electromagnetic origin, to the charged pions:

$$m_{\pi ^\pm }^2|_{\text{EM}}=\frac{\alpha }{\pi }\frac{3}{4F_0^2}_0^{\mathrm{}}𝑑Q^2\left[Q^2\mathrm{\Pi }_{LR}(Q^2)\right].$$
(13)
The first evaluation of this integral was done by F. Low and collaborators in 1967 . What I am going to do next is, simply, to put their phenomenological calculation within the context of $`\text{QCD}_{\mathrm{}}`$.

As already stated, the function $`\mathrm{\Pi }_{LR}(Q^2)`$ in $`\text{QCD}_{\mathrm{}}`$ is a meromorphic function. That means that within a finite radius in the complex $`Q^2`$–plane it only has a finite number of poles. The natural question which arises is: what is the minimal number of poles required to satisfy the OPE constraints? The answer to that follows from a well known theorem in analysis , which when applied to our case, states that the function

$$Q^2\mathrm{\Pi }_{LR}(Q^2)\mathrm{\Delta }[z]\text{with}z=\frac{Q^2}{M_V^2},$$
(14)
and $`M_V`$ the mass of the lowest narrow state, has the property that

$$𝒩𝒫=\frac{1}{2\pi i}\frac{\mathrm{\Delta }^{}[z]}{\mathrm{\Delta }[z]}𝑑z,$$
(15)
where $`𝒩`$ is the number of zeros and $`𝒫`$ is the number of poles inside the integration contour in the complex $`z`$–plane (a zero and/or a pole of order m is counted m times.) For a circle of radius sufficiently large, we simply have that $`𝒩𝒫=p,`$ where $`p`$ denotes the asymptotic power fall–off in $`1/z`$ of the $`\mathrm{\Delta }[z]`$ function. Identifying the power $`p`$ with the leading power predicted by the OPE, and from the fact that $`𝒩0`$, there follows that $`𝒫p_{\text{OPE}}`$. In our case $`p_{\text{OPE}}=2𝒫2`$ and the minimal hadronic approximation (MHA) compatible with the OPE requires two poles 111Notice that with the definition of $`\mathrm{\Delta }[z]`$ in Eq. (14) the pion pole is removed.: one vector state and one axial–vector state. The MHA to the $`\text{QCD}_{\mathrm{}}`$ expression in Eq. (8) is then the simple function

$$Q^2\mathrm{\Pi }_{LR}(Q^2)=F_0^2\frac{M_V^2M_A^2}{(Q^2+M_V^2)(Q^2+M_A^2)}.$$
(16)
Inserting this function in Eq. (13) gives a prediction for the electromagnetic $`\pi ^+\pi ^0\mathrm{\Delta }m_\pi `$ mass difference, with the result 222This is the result for $`F_0=(87\pm 3.5)\text{MeV}`$, $`M_V=(748\pm 29)\text{MeV}`$ and $`g_A=\frac{M_V^2}{M_A^2}=0.50\pm 0.06`$. These values follow from an overall fit to predictions of the low energy constants in the MHA to $`\text{QCD}_{\mathrm{}}`$ .

$$\mathrm{\Delta }m_\pi =(4.9\pm 0.4)\text{MeV},$$
(17)
to be compared with the experimental value

$$\mathrm{\Delta }m_\pi =(4.5936\pm 0.0005)\text{MeV}.$$
(18)
### 3.1 The real world and the MHA to $`\text{QCD}_{\mathrm{}}`$

The spectral function of the correlation function defined in Eq. (2.1) can be obtained from measurements of the hadronic $`\tau `$–decay spectrum (vector–like decays minus axial–vector like decays). We plot in Figure 1 the experimental determination of $`\frac{1}{\pi }\text{Im}\mathrm{\Pi }_{LR}(t)`$, obtained from the ALEPH collaboration data , versus the invariant hadronic mass squared $`t`$ in the accessible region $`0tm_\tau ^2`$.

In this case, the corresponding spectrum of the MHA to $`\text{QCD}_{\mathrm{}}`$ consists of the pion pole (not shown in the figure), a vector narrow state(the $`\rho `$) and an axial–vector narrow state (the $`A_1`$). At this level of approximation, and in the chiral limit, the rest of the vector and axial–vector states are taken to be degenerate and cancel in the difference, a fact which is simulated by the horizontal continuum line in the figure. Looking at this plot, one can hardly claim that this approximation reproduces the details of the experimental data. However, with $`\mathrm{\Pi }_{LR}(Q^2)`$ determined from the spectral function by the unsubtracted dispersion relation

$$\mathrm{\Pi }_{LR}(Q^2)=_0^{\mathrm{}}𝑑t\frac{1}{t+Q^2}\frac{1}{\pi }\text{Im}\mathrm{\Pi }_{LR}(t),$$
(19)
the corresponding plot of the function $`\frac{Q^2}{F_0^2}\mathrm{\Pi }_{LR}(Q^2)`$ versus the euclidean variable $`Q^2=t`$ is shown in Fig 2. The solid curve is the one corresponding to the simple MHA in Eq. (16) and the dotted curve the one from the experimental data in Fig 1. One can see that, by contrast to what happens in the Minkowski region shown in Fig. 1, the corresponding curves in the euclidean region are both very smooth and in fact the MHA already provides a rather good interpolation between the asymptotic regimes where, by construction, it has been constrained to satisfy the lowest order chiral behaviour in Eq. (6), and the two Weinberg sum rules: the OPE constraints in Eqs. (10) and (11). This good interpolation is the reason why the integral in Eq. (7), evaluated at the MHA, already reproduces the experimental result rather well .

### 3.2 Other Analytic Approaches

At this stage, it is illustrative to compare the large–$`N_c`$ approach we have discussed so far with other analytic approaches in the literature. The $`\mathrm{\Pi }_{LR}`$ correlation function provides us with an excellent theoretical laboratory to do this comparison. The different shapes of $`\frac{Q^2}{F_0^2}\mathrm{\Pi }_{LR}(Q^2)`$ predicted by other analytic approaches are collected in Fig. 2. Here, several comments are in order.

a) The suggestion to use a large–$`N_c`$ QCD framework combined with $`\chi `$PT cut–off loops, was first proposed by Bardeen, Buras and Gérard in a series of seminal papers . The same method has been applied by the Dortmund group , in particular to the evaluation of $`ϵ^{}/ϵ`$. In this approach the hadronic ansatz to the Green’s functions consists of Goldstone poles alone and, therefore, integrals like the one in Eq. (7) become UV–divergent (often quadratically divergent) since the correct QCD short–distance behaviour is not implemented. In practice the integrals are cut at a physical cut off: $`\mathrm{\Lambda }1\text{GeV}`$. In the case we are considering, the predicted shape of the LR–correlation function normalized to its value at $`Q^2=0`$ is a constant all the way up to the chosen cut–off value, as shown by the BBG, HKPSB line (the horizontal dotted line) in Fig. 2.

b) The Trieste group evaluate the relevant Green’s functions using the constituent chiral quark model (C$`\chi `$QM) proposed in refs.  and . They have obtained a long list of predictions , in particular $`ϵ^{}/ϵ`$. The model gives an educated first guess of the low–$`Q^2`$ behaviour of the Green’s functions, as one can judge from the C$`\chi `$QM–curve (dashed curve) in Fig. 2, but it fails to reproduce the short–distance QCD–behaviour. Another objection to this approach is that the “natural matching scale” to the short–distance behaviour in this model should be $`4M_Q^2`$ ($`M_Q`$ the constituent quark mass), which is too low to be trusted.

c) The extended Nambu–Jona-Lasinio (ENJL) model was developed as an improvement on the C$`\chi `$QM, since in a certain way it incorporates the vector–like fluctuations of the underlying QCD theory which are known to be phenomenologically important (see e.g. refs.  where other references can also be found). The model is rather successful in predicting the low–energy $`O(p^4)`$ constants of the chiral Lagrangian . It has, indeed, a better low–energy behaviour than the C$`\chi `$QM, as the ENJL–curve (dot–dashed) in Fig. 2 shows; but it fails to reproduce the short–distance behaviour of the OPE in QCD. Arguments, however, to do the matching to short–distance QCD have been forcefully elaborated in refs. , which also give a lot of numerical predictions; a large value for $`ϵ^{}/ϵ`$ in particular.

The problem with the ENJL model as a plausible model of large–$`N_c`$ QCD, is that the on–shell production of unconfined constituent quark $`Q\overline{Q}`$ pairs that it predicts, violates the large–$`N_c`$ QCD counting rules. In fact, as shown in ref. , when the unconfining pieces in the ENJL spectral functions are removed by adding an appropriate series of local counterterms, the resulting theory is entirely equivalent to an effective chiral meson theory with three narrow states V, A and S; very similar to the phenomenological Resonance Chiral Lagrangians proposed in refs. . These Lagrangians can be viewed, therefore, as particular models of large–$`N_c`$ QCD. They predict the same Green’s functions as the MHA to $`\text{QCD}_{\mathrm{}}`$ discussed above, in some particular cases but not in general (see e.g. the three–point functions discussed in ref. ).

In view of the difficulties which the analytic approaches discussed in a), b) and c) above, have in reproducing the shape of the simplest Green’s function one can think of, it is difficult to attribute more than a qualitative significance to their “predictions”; $`ϵ^{}/ϵ`$ in particular, which requires the interplay of several other Green’s functions much more complex than $`\mathrm{\Pi }_{LR}(Q^2)`$.

## 4  METHODOLOGY, APPLICATIONS

The $`\text{QCD}_{\mathrm{}}`$ approach that we are proposing in order to compute a specific coupling of the chiral electroweak Lagrangian consists of the following steps:

1. Identify the relevant Green’s functions.

In most cases of interest, the underlying QCD Green’s functions in question are two–point functions with additional zero momentum insertions of vector, axial vector, scalar and pseudoscalar currents. The higher the power in the chiral expansion, the higher will be the number of zero momentum insertions. This step is totally general and does not invoke any large–$`N_c`$ approximation.
2. Work out the short–distance behaviour and the long–distance behaviour of the relevant Green’s functions.

The long–distance behaviour is governed by the Goldstone singularities and can be obtained from $`\chi `$PT. The short–distance behaviour is governed by the OPE of the currents through which the hard momenta flows. Again, this step is well defined independently of the large–$`N_c`$ expansion; in practice, however, the calculations simplify when restricted to the appropriate order in the $`1/N_c`$–expansion one is interested in. In particular, chiral loops are subleading in the $`1/N_c`$–expansion.
3. Introduce a large–$`N_c`$ approximation for the underlying Green’s functions.

As already mentioned, Green’s functions in $`\text{QCD}_{\mathrm{}}`$ are meromorphic functions. The minimal hadronic approximation (MHA) that we are proposing consists in limiting the number of poles to the minimum required to satisfy the leading power fall–off at short–distances, as well as the appropriate $`\chi `$PT long–distance constraints.

The three steps above can be done analytically, which helps to unravel the intricacies of the underlying dynamics. The method is, in principle, improvable (unlike the other models discussed above). It can be improved by adding more constraints from the next–to–leading short–distance inverse power behaviour in the OPE and/or higher orders in the chiral expansion. This has been tested within a toy model of $`QCD_{\mathrm{}}`$ in ref.

We have checked this approach with the calculation of a few low–energy observables:

i) The electroweak $`\mathrm{\Delta }m_\pi `$ mass difference which we have already discussed .

ii) The hadronic vacuum polarization contribution to the anomalous magnetic moment of the muon $`a_\mu `$. The MHA in this case requires one vector–state pole and a perturbative QCD continuum. The absence of dimension two operators in QCD in the chiral limit, constrains the threshold of the continuum. The result thus obtained  is compatible with the more precise phenomenological determinations which use experimental input.

iii) The $`\pi ^0e^+e^{}`$ and $`\eta \mu ^+\mu ^{}`$ decay rates. These processes are governed by a $`PVV`$ three–point function, with the $`Q^2`$–momentum flowing through the two $`V`$–currents. The MHA in this case requires a vector–pole and a double vector–pole. The predictions  of the branching ratios $`\frac{\mathrm{\Gamma }(Pl^+l^{})}{\mathrm{\Gamma }(P\gamma \gamma )}`$ are in good agreement with the present experimental determinations.

These successful predictions have encouraged some of us to pursue a systematic analysis of $`K`$–physics observables within the same large–$`N_c`$ framework. So far, the following calculations have been made:

* The $`B_K`$–Factor in the Chiral Limit .
* Weak Matrix Elements of the Electroweak Penguin Operators $`Q_7`$ and $`Q_8`$ .
* Hadronic Light–by–Light Scattering Contribution to the Muon $`g2`$ .
* Electroweak Hadronic Contributions to the Muon $`g2`$ .

We also have preliminary results on:

* The Long–Distance Contribution to the $`K_L\mu ^+\mu ^{}`$ Decay Rate .
* Weak Matrix Elements of the Strong Penguin Operators $`Q_4`$ and $`Q_6`$ .

## 5 FACTOR $`B_K`$ OF $`K^0\overline{K}^0`$ MIXING

The factor in question is conventionally defined by the matrix element of the four–quark operator $`Q_{\mathrm{\Delta }S=2}(x)=(\overline{s}_L\gamma ^\mu d_L)(\overline{s}_L\gamma _\mu d_L)(x)`$:

$$\overline{K}^0|Q_{\mathrm{\Delta }S=2}(0)|K^0=\frac{4}{3}f_K^2M_K^2B_K(\mu ).$$
(20)
To lowest order in the chiral expansion the operator $`Q_{\mathrm{\Delta }S=2}(x)`$ bosonizes into an $`O(p^2)`$ term:

$$\frac{F_0^4}{4}g_{\mathrm{\Delta }S=2}(\mu )\underset{\frac{2}{F_0^2}^\mu K^0_\mu K^0(r^\mu )_{23}(r_\mu )_{23}+\mathrm{}}{\underset{}{\left[(D^\mu U^{})U\right]_{23}\left[(D_\mu U^{})U\right]_{23}}},$$
(21)
with $`g_{\mathrm{\Delta }S=2}(\mu )`$ the low energy constant (which depends on the renormalization scale $`\mu `$) to be determined. A convenient choice of the underlying Green’s function here is the four–point function $`W_{LRLR}^{\mu \alpha \nu \beta }(q,l)`$ of two left–currents $`(\overline{s}_L\gamma ^\mu d_L)`$ and $`(\overline{s}_L\gamma ^\nu d_L)`$, which carry the $`q`$–momentum one has to integrate over, and two right–currents $`(\overline{d}_R\gamma ^\alpha s_R)`$ and $`(\overline{d}_R\gamma ^\beta s_R)`$ with soft $`l`$–momentum insertions, which couple to the external $`(r_\alpha )_{32}`$– and $`(r_\beta )_{32}`$–sources. The coupling constant $`g_{\mathrm{\Delta }S=2}(\mu )`$, which has to be evaluated in the same renormalization scheme as the Wilson coefficient $`C_{\mathrm{\Delta }S=2}(\mu )`$ has been evaluated, is then given by the integral

$$g_{\mathrm{\Delta }S=2}(\mu )=1\frac{1}{32\pi ^2F_0^2}[\frac{\left(4\pi \mu ^2\right)^{ϵ/2}}{\mathrm{\Gamma }(2ϵ/2)}_0^{\mathrm{}}\frac{dQ^2}{(Q^2)^{ϵ/2}}\times $$

$$\underset{W_{LL}(Q^2)}{\underset{}{(\frac{g_{\alpha \beta }}{3})𝑑\mathrm{\Omega }_q\underset{l0}{lim}g_{\mu \nu }W_{LRLR}^{\mu \alpha \nu \beta }(q,l)}}]_{\overline{MS}},$$
(22)
conceptually similar to the one which determines the electroweak constant $`C`$ in Eq. (7). The renormalization–scale invariant $`\widehat{B}_K`$–factor is then given by the product

$$\widehat{B}_K=\frac{3}{4}C_{\mathrm{\Delta }S=2}(\mu )\times g_{\mathrm{\Delta }S=2}(\mu ).$$
(23)
The large–$`N_c`$ hadronic approximation of the Green’s function $`W_{LL}(Q^2)`$ in Eq. (22), which fulfills the leading OPE short–distance constraint and the long–distance constraints which fix $`W_{LL}(0)`$ and $`W_{LL}^\mathrm{`}(0)`$ in $`\chi `$PT, requires one vector–pole, a double vector–pole and a triple vector–pole 333In fact, this goes beyond the strict MHA which here only requires a vector–pole. It is the extra input of $`W_{LL}(0)`$ and $`W_{LL}^{}(0)`$, as known from $`\chi `$PT, which allows us to improve on the strict MHA. The integral in Eq. (22) can then be evaluated, with the result

$$\widehat{B}_K=0.38\pm 0.11.$$
(24)
The $`\mu `$–scale and renormalization scheme dependence in $`g_{\mathrm{\Delta }S=2}(\mu )`$ cancels with the one in $`C_{\mathrm{\Delta }S=2}(\mu )`$, when evaluated at the same approximation in the $`I/N_c`$–expansion.

When comparing this result to other determinations, specially in Lattice QCD, it should be realized that the unfactorized contribution in Eq. (22) is the one in the chiral limit. It is possible, in principle, to calculate chiral corrections within the same large–$`N_c`$ approach, but this has not yet been done.

The result in Eq. (24) is compatible with the old current algebra prediction  which, to lowest order in $`\chi `$PT, relates the $`B_K`$–factor to the $`K^+\pi ^+\pi ^0`$ decay rate. In fact, our calculation can be viewed as a successful prediction of the $`K^+\pi ^+\pi ^0`$ decay rate !

As discussed in ref. , the bosonization of the four–quark operator $`Q_{\mathrm{\Delta }S=2}`$ and the bosonization of the operator $`Q_2Q_1`$ which generates $`\mathrm{\Delta }I=1/2`$ transitions, are related to each other in the combined chiral limit and next–to–leading order in the $`1/N_c`$–expansion, when mixing with the penguin operators is neglected. Lowering the value of $`\widehat{B}_K`$ from the factorized $`𝒪(N_c^2)`$ prediction: $`\widehat{B}_K=3/4`$, to the result in Eq. (24), is correlated with an increase of the coupling constant $`g_{\underset{¯}{8}}`$ in the lowest order effective chiral Lagrangian (see Eq. (4)) which generates $`\mathrm{\Delta }I=1/2`$ transitions, and provides a first step towards a quantitative understanding of underlying dynamics of the $`\mathrm{\Delta }I=1/2`$ rule.

## 6 ELECTROWEAK PENGUINS

We shall be next concerned with the four–quark operators generated by the so called electroweak penguin diagrams of the Electroweak Theory

$$\mathrm{}C_7(\mu )Q_7+C_8(\mu )Q_8,$$
(25)
with $`C_7(\mu )`$, $`C_8(\mu )`$ the Wilson coefficients of the operators

$$Q_7=6(\overline{s}_L\gamma ^\mu d_L)\underset{q=u,d,s}{}e_q(\overline{q}_R\gamma _\mu q_R),$$
(26)
$$Q_8=12\underset{q=u,d,s}{}e_q(\overline{s}_Lq_R)(\overline{q}_Rd_L).$$
(27)
These operators generate terms of $`O(p^0)`$ in the effective chiral Lagrangian; therefore, their matrix elements, although suppressed by an $`e^2`$ factor, are chirally enhanced. Furthermore, the Wilson coefficient $`C_8`$ has a large imaginary part induced by the top–quark integration, which makes the matrix elements of the $`Q_8`$ operator to be particularly important in the evaluation of $`ϵ^{}/ϵ`$.

Within the large–$`N_c`$ framework, the bosonization of these operators produce matrix elements with the following counting

$$Q_7|_{O(p^0)}=\underset{¯}{O(N_c)}+O(N_c^0),$$
(28)
$$Q_8|_{O(p^0)}=\underset{¯}{O(N_c^2)}+\underset{\text{Zweig suppressed}}{\underset{}{O(N_c^0)}}.$$
(29)
### 6.1 Bosonization of $`Q_7`$

The bosonization of the $`Q_7`$ operator to $`O(p^0)`$ in the chiral expansion and to $`O(N_c)`$ is very similar to the calculation of the $`Z`$–contribution to the coupling constant $`C`$ in Eq. (7). An evaluation which also takes into account the renormalization scheme dependence has been recently made in ref.  with the result:

$$Q_7|_{𝒪(p^0)}=6\underset{O_1(\mu )}{\underset{}{0|(\overline{s}_L\gamma ^\mu d_L)(\overline{d}_R\gamma _\mu s_R)|0}}\times $$

$$\text{tr}\left(U\lambda _L^{(23)}U^{}Q_R\right)^{},\text{and}\lambda _L^{(23)}=\delta _{i2}\delta _{j3}.$$
(30)
Here, the vev $`O_1(\mu )`$ is given by the integral

$$O_1(\mu )=\{\frac{3(ϵ3)}{16\pi ^2}\frac{(4\pi \mu ^2)^{ϵ/2}}{\mathrm{\Gamma }(2ϵ/2)}\times $$

$$_0^{\mathrm{}}dQ^2(Q^2)^{1ϵ/2}(Q^2\mathrm{\Pi }_{LR}(Q^2))\}_{\overline{MS}},$$
(31)
with $`\mathrm{\Pi }_{LR}(Q^2)`$ the same correlation function as in Eq. (2.1). In $`\text{QCD}_{\mathrm{}}`$ this integral can be, formally, evaluated exactly. When restricted to the same MHA which has been used to calculate $`\mathrm{\Delta }m_\pi `$, one gets the simple result

$$O_1=\frac{3}{32\pi ^2}\left[f_V^2M_V^6\mathrm{ln}\frac{\mathrm{\Lambda }^2}{M_V^2}f_A^2M_A^6\mathrm{ln}\frac{\mathrm{\Lambda }^2}{M_A^2}\right],$$
(32)
where $`\mathrm{\Lambda }^2=\mu ^2\mathrm{exp}(1/3+\kappa )`$, with $`\kappa =1/2`$ in the naive dimensional renormalization scheme (NDR) and $`\kappa =+3/2`$ in the ’t Hooft–Veltman scheme (HV).

### 6.2 Bosonization of $`Q_8`$

To lowest $`𝒪(p^0)`$ in the chiral expansion, the four–quark operator $`Q_8`$ bosonizes as follows

$$Q_8|_{𝒪(p^0)}=12\underset{0|(\overline{s}_Ls_R)(\overline{d}_Rd_L)|0}{\underset{}{O_2(\mu )}}\text{tr}\left(U\lambda _L^{(23)}U^{}Q_R\right)^{}.$$

As noted in refs. , the vev $`O_2(\mu )`$ also appears in the Wilson coefficient of the $`1/Q^6`$ term in the OPE of the same $`\mathrm{\Pi }_{LR}(Q^2)`$ correlation function as in Eq. (2.1), for which the MHA to $`\text{QCD}_{\mathrm{}}`$ gives a rather good approximation, as we have already discussed. This offers the possibility of obtaining an estimate of the vev $`O_2(\mu )`$ beyond the large–$`N_c`$ approximation, where $`O_2\frac{1}{4}\overline{\psi }\psi ^2`$, and, therefore, without having to fix a value for the $`\overline{\psi }\psi `$–condensate, which is poorly known. Unfortunately, as one can judge from the results in Table 1, different methods to extract the value of $`O_2(\mu )`$ lead, at present, to rather different values for the weak matrix elements:

$$M_{7,8}(\pi \pi )_{I=2}|Q_{7,8}|0,\text{at}\mu =2\text{GeV}.$$
(33)
In view of the diversity of results obtained with the dispersive approach, depending on which type of sum rules are used; and the unknown systematic errors of the quenched approximation and the chiral limit extrapolations in the lattice results, it looks like it will take some time before one can claim that these matrix elements are known reliably.

### 6.3 Test of the Zweig Rule

Besides the important issue of getting an accurate determination of the $`Q_8`$ matrix elements, a reliable determination of the vev $`0|(\overline{s}_Ls_R)(\overline{d}_Rd_L)|0`$ would be most welcome as a test of the Zweig rule in the scalar and pseudoscalar sectors. More precisely,

$$O_2(\mu )=\frac{1}{4}\overline{\psi }\psi ^2+0|(\overline{s}_Ls_R)(\overline{d}_Rd_L)|0_c,$$
(34)
where the unfactorized contribution (the second term) involves Feynman diagrams which require gluon exchanges between at least two quark loops. These are the so called Zweig–suppressed contributions, which are indeed $`𝒪(N_c^0)`$ in the $`1/N_c`$ expansion. There are reasons to suspect that subleading terms in the $`1/N_c`$ expansion involving Green’s functions of scalar (and pseudoscalar) density operators might be important, unlike those which only involve vector (and axial-vector) currents. There are several phenomenological examples of this: the $`\eta ^{}`$ mass, the possible existence of a broad $`\sigma `$ meson, large final state interactions in states with $`J=0`$ and $`I=0`$, etc. We are considering the possibility that the appropriate expansion for these exceptional Green’s functions could be a $`1/N_c`$ expansion in which $`n_f/N_c`$ is held fixed, where $`n_f`$ denotes the number of light flavours; the kind of expansion originally advocated by G. Veneziano .

Formally, the connected part of the vev $`0|(\overline{s}_Ls_R)(\overline{d}_Rd_L)|0_c`$ is defined by the integral

$$\frac{1}{i}\left(\frac{d^Dq}{(2\pi )^D}\mathrm{\Psi }_{ds}(Q^2)\right)_{\overline{\text{MS}}}^{\text{ren.}},$$
(35)
with $`\mathrm{\Psi }_{ij}(Q^2)`$ the two–point function

$$\frac{i}{4}d^4xe^{iq.x}\{0|T\left[\overline{d}d(x)\overline{s}s(0)\right]|0$$

$$0|T\left[\overline{d}\gamma _5d(x)\overline{s}\gamma _5s(0)\right]|0\}.$$
(36)
It would be very helpful to get some information on the vev $`0|(\overline{s}_Ls_R)(\overline{d}_Rd_L)|0_c`$ from lattice QCD. Ultimately, this could provide a way to focus on the origin of the discrepancies in Table 1.

## 7 GLUONIC PENGUINS

Finally, I shall make a few analytic remarks concerning the sector of the four–quark operators generated by the strong Penguins of the Electroweak Theory

$$\mathrm{}C_4(\mu )Q_4+C_6(\mu )Q_6,$$
(37)
with $`C_4(\mu )`$, $`C_6(\mu )`$ the Wilson coefficients of the operators

$$Q_4=4\underset{q=u,d,s}{}(\overline{s}_L\gamma ^\mu q_L)(\overline{q}_L\gamma ^\mu d_L),$$
(38)
$$Q_6=8\underset{q=u,d,s}{}(\overline{s}_Lq_R)(\overline{q}_Rd_L).$$
(39)
These operators generate contributions to the coupling constant $`g_{\underset{¯}{8}}`$ of the $`𝒪(p^2)`$ chiral Lagrangian in Eq. (4). The Wilson coefficient $`C_6(\mu )`$ has also an important imaginary part generated by the integration of the heavy flavour degrees of freedom, which makes the operator $`Q_6`$ particularly relevant for the evaluation of $`ϵ^{}/ϵ`$ in the Standard Model.

Seen from the large–$`N_c`$ framework point of view, the only terms which have been retained (so far) in analytic evaluations of weak matrix elements of four–quark operators, is best illustrated by showing the terms which, equivalently, would be retained in the lowest order anomalous dimension matrix of these operators. Non–zero entries in this matrix appear then in three blocks, in the sector of the $`Q_1`$, $`Q_2`$, $`Q_4`$, $`Q_6`$ and $`Q_8`$ operators, as follows:

$`\left(\begin{array}{ccccc}0& \frac{3}{N_c}& & & \\ \frac{3}{N_c}& 0& & & \\ & & \frac{1}{3}\frac{n_f}{N_c}& \frac{1}{3}\frac{n_f}{N_c}& \\ & & \frac{1}{3}\frac{n_f}{N_c}& 3+\frac{1}{3}\frac{n_f}{N_c}& \\ & & & & 3+\frac{3}{N_c^2}\end{array}\right)`$ (45)
More precisely, the mixing between matrix elements of the $`(Q_2,Q_1)`$–operators on the one hand and those of the penguin–like operators on the other, is neglected. Furthermore, in the $`(Q_2,Q_1)`$–sector, only the leading and next–to–leading terms in the $`1/N_c`$–expansion have been retained . Mixing between the strong penguin sector and the electroweak penguin sector is also neglected, but terms $`𝒪(\frac{n_f}{N_c})`$ in the mixing of $`Q_4`$ and $`Q_6`$ are retained. Exceptionally, for reasons already discussed above, the Zweig suppressed $`𝒪(1/N_c^2)`$–corrections in $`Q_8`$ are also retained. I can now discuss some recent analytic observations which have been made.

### 7.1 Final State Interactions and $`ϵ^{}/ϵ`$

The importance of final state interactions in the evaluation of $`K\pi \pi `$ matrix elements has been reexamined in a series of recent papers . The main observation is that a simple numerical estimate of $`ϵ^{}/ϵ`$ can be obtained if one proceeds as follows:

i) The evolution of the relevant Wilson coefficients from the $`M_W`$–scale to the $`m_c`$–scale is retained, exactly, to next to leading order in pQCD.

ii) The bosonization of the relevant $`Q_6`$ and $`Q_8`$ operators is done at the leading $`𝒪(N_c^2)`$ only. Notice that this corresponds to the further restriction where the matrix in Eq. (45) becomes diagonal, with only non–zero entries in the ($`Q_6`$, $`Q_8`$) sector: $`\left(\begin{array}{cc}3& 0\\ 0& 3\end{array}\right)`$. Starting at this level, the only $`1/N_c`$ contributions which are then retained, are those generated by an approximate estimate of the the $`\pi \pi `$ final state interactions.

iii) The $`\pi \pi `$ final state interactions are approximated by an Omnès–like resummation of the leading $`\pi \pi `$ chiral loops with non–zero discontinuities. The authors of ref.  justify this procedure by the fact that it reproduces rather well some of the phenomenologically known $`𝒪(p^4)`$ local terms of the electroweak chiral Lagrangian.

The prediction thus obtained

$$\text{Re}(\frac{ϵ}{ϵ})=(1.7\pm \underset{\mu \text{scale}}{\underset{}{0.2}}\underset{\overline{\psi }\psi }{\underset{}{\begin{array}{c}+0.8\\ 0.5\end{array}}}\pm \underset{1/N_c}{\underset{}{0.5}})\times 10^3,$$
(46)
where their estimates of the various sources of errors are indicated in underbracing, agrees within errors with the present experimental world average  from NA31, NA48 and KTeV:

$$\text{Re}(ϵ/ϵ)|_{\text{Exp. }}=(1.66\pm 0.16)\times 10^3.$$
(47)
The obvious question here is to know how stable is this prediction when terms of $`𝒪(N_c)`$ are also retained in the bosonization of $`Q_6`$ and $`Q_8`$. The scenario with large deviations from the naive $`𝒪(N_c^2)`$ bosonization of $`Q_6`$ and $`Q_8`$, which cancel in their overall contribution to $`ϵ^{}/ϵ`$, cannot be excluded. Let us not forget that the restriction to $`𝒪(N_c^2)`$ terms in the bosonization of the four–quark operators, fails dramatically to reproduce the $`\mathrm{\Delta }I=1/2`$ rule; a fact which, perhaps, should not be so surprising since, after all, it is only when the $`𝒪(N_c)`$ terms are incorporated as well, that the $`\text{QCD}_{\mathrm{}}`$–properties of planarity apply.

### 7.2 Bosonization of $`Q_6`$

Related to the previous discussion is the question of the bosonization of the $`Q_6`$–operator. It has been known for sometime that this operator, to $`𝒪(N_c^2)`$, gives a contribution to the coupling constant $`g_{\underset{¯}{8}}`$ in Eq. (4) which is modulated by the product of the ratio $`\overline{\psi }\psi ^2/F_0^6`$ times the $`𝒪(p^4)`$ $`L_5`$–coupling of the strong chiral Lagrangian (known from the $`f_\pi /f_K`$ ratio). It can be shown that, in the presence of unfactorized terms, this contribution is modified as follows :

$$g_{\underset{¯}{8}}|_{Q_6}=C_6(\mu )\{16L_5\frac{\overline{\psi }\psi ^2}{F_0^6}+[\frac{1}{2\pi ^2F_0^4}\frac{\left(4\pi \mu ^2\right)^{ϵ/2}}{\mathrm{\Gamma }(2ϵ/2)}$$

$$\times _0^{\mathrm{}}dQ^2(Q^2)^{1ϵ/2}𝒲_{DG}(Q^2)]_{\overline{MS}}\},$$
(48)
where the function $`𝒲_{DG}(Q^2)`$ is the analog of the function $`W_{LL}(Q^2)`$ in Eq. (22) and the function $`\mathrm{\Pi }_{LR}(Q^2)`$ in Eqs. (7) and (2.1). Here, $`𝒲_{DG}(Q^2)`$ is the invariant function associated with a four–point function $`𝒲_{DGRR}^{\alpha \beta }(q,l)`$ of two density–currents $`(\overline{s}_Lq_R)`$ and $`(\overline{q}_Rd_L)`$, which carry the $`q`$–momentum one has to integrate over, and two right–currents $`(\overline{d}_R\gamma ^\alpha u_R)`$ and $`(\overline{u}_R\gamma ^\beta s_R)`$ with soft $`l`$–momentum insertions, which couple to the external $`(r_\alpha )_{12}`$– and $`(r_\beta )_{31}`$–sources. We have found that the effect of the unfactorized piece in Eq. (48) is very important. The large effect can already be seen in the low $`Q^2`$ behaviour of the function $`𝒲_{DG}(Q^2)`$ in $`\chi `$PT ($`B\frac{|\overline{\psi }\psi |}{F_0^2})`$:

$$\underset{Q^20}{lim}𝒲_{DG}(Q^2)=\frac{n_f}{8}\left(\frac{BF_0}{Q^2}\right)^2$$

$$n_f\left(L_5\frac{5}{2}L_3\right)\frac{B^2}{Q^2}+\mathrm{},$$
(49)
where the factor $`L_5\frac{5}{2}L_311\times 10^3`$, which governs the behaviour of $`Q^2𝒲_{DG}(Q^2)`$ at the origin (the integral of the first term vanishes in dimensional regularization), is very large as compared to the individual values of $`L_5`$ and $`L_3`$. This large value results mostly from the effect of the lowest vector state which the first $`𝒪(N_c^2)`$ term in Eq. (49) does not see. Clearly, this has serious implications both for the $`\mathrm{\Delta }I=1/2`$ rule and $`ϵ^{}/ϵ`$. An evaluation, along these lines, of the coupling $`g_{\underset{¯}{8}}|_{Q_6}`$ is in progress.

## 8 CONCLUSIONS and OUTLOOK

We think that large–$`N_c`$ QCD provides a very useful framework to formulate approximate calculations of the low–energy constants of the effective chiral Lagrangian, both in the strong and electroweak sector.

The approach that we propose has been tested successfully with the calculation of a few low–energy observables discussed in section 4. Here, we have described various applications to the evaluation of weak matrix elements which, for the first time in analytic calculations of this type, explicitly show the cancellation of the renormalization scale between the short–distance contributions and the long–distance contributions. The analytic implementation of the approach is sufficiently simple so as to provide an explanation of the dominant underlying physical effects. It may also help, as a guidance, to unravel important physical effects in Lattice QCD numerical simulations.

Acknowledgements

My knowledge on the subjects reported here owes much to work and discussions with my collaborators Santi Peris and Marc Knecht; and also, at different stages, with M. Perrottet, A. Pich, J. Bijnens, M. Golterman, T. Hambye, A. Nyffeler, and my students: B. Phily, S. Friot and D. Greynat. I take this opportunity to thank them all. Special thanks to Laurent Lellouch and Santi Peris for a careful reading of the manuscript.


format: text

text:
# Self-Consistent Hartree-Fock Based Random Phase Approximation And The Spurious State Mixing

## I Introduction

Hartree-Fock (HF) based random phase approximation (RPA) has been a very successful theory in providing microscopic description of phenomena associated with collective motion in nuclei . Accurate information for important physical quantities can be obtained by comparing the experimentally deduced strength function distribution, $`S(E)`$, with the results obtained from HF-RPA theory. In particular, the strength function distributions of the isoscalar giant monopole resonance (ISGMR) and the isoscalar giant dipole resonance (ISGDR) are quite sensitive to the value of the nuclear matter incompressibility coefficient, $`K_{nm}`$, , a very important physical quantity since it is directly related to the curvature of the equation of state.

Over the last two decades, a significant amount of experimental work was carried out to identify the strength distributions of the isoscalar giant resonances in nuclei, particularly the ISGMR  and ISGDR . The main development in the area of experimental investigation of the isoscalar giant resonances is the high accuracy data, of excitation cross section, by $`\alpha `$particle scattering, obtained at Texas A&M University using a beam analysis system (BAS), a multipole-dipole-multipole (MDM) spectrometer and broad range multiwire proportional counter. The new system improved the signal to background ratio by more than a factor of 15. This led to the discovery of a high lying structure in the strength function of the ISGMR and the location of the ISGMR in light nuclei. Also, accurate data on the ISGDR has been obtained for a wide range of nuclei . This has led to renewed interest in the nuclear response function and the need to carry out detailed and accurate calculations of $`S(E)`$ and the transition density, $`\rho _t`$, within the HF-RPA theory. In particular there have been quite a few recent HF-RPA and relativistic mean field (RMF) based RPA calculations of the ISGDR , considering the issues of (i) spurious state mixing (SSM), (ii) the strength of the lower component (at 1$`\mathrm{}\omega `$) and (iii) the value of $`K_{nm}`$ deduced from the centroid energy $`E_1`$ of the ISGDR compression mode (at 3$`\mathrm{}\omega `$).

Comparison between the recent data on the ISGMR and the results of HF based RPA calculations confirms the value of $`K_{nm}=210\pm 20`$ MeV, determined earlier in . It was first pointed out in Ref.  that the HF-RPA results for $`E_1`$, obtained with interactions adjusted to reproduce the ISGMR data, are higher than the experimental values  by more than 3 MeV and thus this discrepancy between theory and experiment raises doubts concerning the unambiguous extraction of $`K_{nm}`$ from energies of compression modes. This discrepancy between theory and experiment was also reported in more recent experiments . Recently, Shlomo and Sanzhur have addressed this discrepancy by carrying out accurate microscopic calculations for $`S(E)`$ and the excitation cross section $`\sigma (E)`$ of the ISGDR, within the folding model (FM) distorted-wave-Born-approximation (DWBA), with $`\rho _t`$ obtained from HF-RPA calculations and corrected for the SSM. They demonstrated that the calculated $`\sigma (E)`$ drops below the experimental sensitivity in the region of high excitation energy containing 30-40% of the ISGDR energy weighted sum rule (EWSR). This missing strength leads to a reduction of more than 3.0 MeV in the value of $`E_1`$ and thus explains the discrepancy between theory and experiment.

Clearly accurate calculations of $`S(E)`$ and $`\sigma (E)`$ are needed. In fully self-consistent HF-RPA calculations, the spurious isoscalar dipole ($`T=0`$, $`L=1`$) state (associated with the center of mass motion) appears at energy $`E=0`$ and no SSM in the ISGDR occurs. It was pointed out in  that none of the calculations carried out for $`S(E)`$ and $`\rho _t`$ and published in the literature are fully self-consistent. In some RPA calculations the mean field and the particle-hole interaction $`V_{ph}`$ are chosen independently. Although this approach can provide physical insight on the structure of giant resonances, it can not be used to accurately determine important physical quantity such as $`K_{nm}`$. In self-consistent HF-RPA calculation one starts by adopting specific effective nucleon-nucleon interaction, $`V_{12}`$, such as the Skyrme interaction, and carries out HF calculations. The parameters of the interaction are determined by a fit to properties of nuclei (binding energies, radii, etc.). Then one solves the RPA equation using the particle-hole (p-h) interaction $`V_{ph}`$ which corresponds to $`V_{12}`$. However, although not always stated in the literature, self-consistency is violated in actual implementations of the RPA (and relativistic RPA) calculations. One usually makes the following approximations: (i) use a $`V_{ph}`$ which is not consistent with $`V_{12}`$. It is common to neglect the two-body Coulomb and spin-orbit interactions in $`V_{ph}`$ and approximate the momentum dependent parts in $`V_{ph}`$, (ii) limiting the p-h space in a discretized calculation by a cut-off energy $`E_{ph}^{max}`$, and (iii) introducing a smearing parameter (i.e., a Lorentzian with $`\mathrm{\Gamma }/2`$). The consequences of these violations of self-consistency on $`S(E)`$ and $`\rho _t`$ and of numerical inaccuracy are usually ignored in the literature.

In this work we present results of detailed investigations of the consequences of common violations of self-consistency in actual implementations of HF based RPA, for determining the response functions $`S(E)`$ and $`\rho _t`$ of isoscalar multipole (L=0,1 and 2) giant resonances. In particular, we consider the ISGDR and concentrate on the effects of the SSM. We determine the effects of a violation of self-consistency by comparing the calculated results for $`S(E)`$and $`\rho _t`$ with those obtained from highly accurate fully self-consistent HF- continuum RPA (HF$``$CRPA) calculations. We also extend the projection method for eliminating the effects of SSM, described in Ref. , to properly normalize $`S(E)`$ and $`\rho _t`$ and determine the mixing amplitude of the spurious state in the ISGDR.

In Section II we present an extension of the Green’s function based derivation of the projection operator method for eliminating the effects of the SSM, described in , to also account for the proper normalization of the $`S(E)`$ and $`\rho _t(𝐫)`$ of the ISGDR and determine the mixing amplitude of the spurious state, obtained in HF-RPA calculations. We emphasize here that the method is quite general and applicable for any scattering operator $`F`$ and for any numerical method used in carrying out the RPA calculation, such as configuration space RPA, coordinate space (continuum and discretized) RPA and with and without the addition of smearing. We also provide in this section the basic expressions used in the calculations and the presentation of our results.

In Section III we present and discuss our results. We first present the results of a highly accurate and fully self-consistent HF$``$CRPA calculation of $`S(E)`$ and $`\rho _t(r)`$ in $`{}_{}{}^{80}Zr`$, which we use as a basis for a comparison with results obtained with common violations of self-consistency. These accurate fully self-consistent HF$``$CRPA results were obtained using $`\mathrm{\Gamma }=0`$ (i.e., no smearing) and very small mesh sizes of $`dr_{HF}=0.04`$ fm and $`dr_{RPA}=0.04`$ fm with corresponding number of mesh points $`N_{HF}=900`$ and $`N_{RPA}=300`$, used in the HF and the CRPA calculations, respectively. We note that the values of $`S(E)`$ and $`\rho _t(r)`$ associated with a bound RPA state were deduced from the residue of the RPA Green’s function. Next, we present our results of fully self-consistent HF$``$CRPA calculations (with $`\mathrm{\Gamma }=0`$) carried out using various mesh sizes $`dr_{HF}`$ and $`dr_{RPA}`$ and discuss the issue of numerical accuracy. We then present and discuss the results obtained with certain violations of self-consistency in CRPA and discretized RPA (DRPA) calculations and assess the effects on $`S(E)`$ and $`\rho _t(E)`$ by comparing with the highly accurate fully self-consistent results over the whole range of excitation energies. We point out that comparing the total energy weighted transition strength with the EWSR may lead to incorrect conclusions. Very recently the accuracy of the projection operator method in eliminating the effects of the SSM on $`S(E)`$ and $`\rho _t`$ of the ISGDR was investigated in Refs. . However, in these works, the calculations carried out using mesh sizes $`dr0.1`$ fm, were not fully self-consistent. We emphasize that in the present work we have carried out highly accurate self-consistent calculations, established the accuracy of the projection operator method and provide assessments on the effects on $`S(E)`$, $`E_L`$ and $`\rho _t`$ of the isoscalar resonances with $`L=0`$, 1 and 2, which are due to common violation of self-consistency in actual implementation of HF$``$RPA often encountered in the literature. We note that preliminary results of the present work were presented earlier . In section IV we state our conclusion.

## II Formalism

The RPA Green’s function $`G`$  is given by,

$$G=G_0(1+V_{ph}G_0)^1,$$
(1)
where $`G_0`$ is the free p-h Green’s function given by,

$$G_0(𝒓,𝒓^{},E)=\underset{h}{}\varphi _h(𝒓)\left[\frac{1}{H_0ϵ_h\omega }+\frac{1}{H_0ϵ_h+\omega }\right]\varphi _h(𝒓^{}).$$
(2)
Here $`H_0`$ is the HF hamiltonian and $`ϵ_h`$ and $`\psi _h`$ are the single particle energy and the wave function of the occupied state, respectively. The continuum effects (particle escape width) are included by using

$$G_{lj}(𝒓,𝒓^{},E)=\frac{1}{H_0E}=\frac{2m}{\mathrm{}^2}u_{lj}(r_<)v_{lj}(r_>)/w,$$
(3)
where $`r_<`$ and $`r_>`$ are the the lesser and the greater of $`r`$ and $`r^{}`$, respectively, $`u`$ and $`v`$ are the regular and irregular solution of $`H_0`$, with the appropriate boundary conditions, respectively and $`w`$ is the Wronskian. The strength function $`S(E)`$ and transition density $`\rho _t`$, associated with the scattering operator,

$$F=\underset{i=1}{\overset{A}{}}f(𝒓_i),$$
(4)
are given by,

$$S(E)=\underset{n}{}\left|0|F|n\right|^2\delta (EE_n)=\frac{1}{\pi }\text{Im}\left[\text{Tr}(fGf)\right],$$
(5)
$$\rho _t(𝒓,E)=\frac{\mathrm{\Delta }E}{\sqrt{S(E)\mathrm{\Delta }E}}f(𝒓^{})\left[\frac{1}{\pi }\text{Im}G(𝒓^{},𝒓,E)\right]𝑑𝒓^{}.$$
(6)
Note that $`\rho _t(𝒓,E)`$, as defined in (6), is associated with the strength in the region of $`E\pm \mathrm{\Delta }E/2`$ and is consistent with

$$S(E)=\left|\rho _t(𝒓,E)f(𝒓)𝑑𝒓\right|^2/\mathrm{\Delta }E.$$
(7)
It is important to note that $`S(E)`$ and $`\rho _t`$ of a state at energy $`E_n`$ below the particle escape threshold (or having a very small width) can be obtained from Eqs. (5) and (6), respectively, by replacing $`\frac{1}{\pi }\text{Im}G(𝒓^{},𝒓,E)`$ with

$$\underset{EE_n}{lim}\text{Re}G(𝒓^{},𝒓,E)(EE_n),$$
(8)
The energy weighted sum rule (EWSR) associated with the operator $`f_{LM}=f(r)Y_{LM}`$ is given by,

$$EWSR(fY_{LM})=ES_{LM}(E)𝑑E=\frac{\mathrm{}^2}{2m}\frac{A}{4\pi }\left[0|\left(\frac{df}{dr}\right)^2+L(L+1)\left(\frac{f}{r}\right)^2|0\right].$$
(9)
Using the equation of continuity and assuming that there is only one collective state  with energy $`E_{coll}`$, exhausting 100% of the EWSR associated with the scattering operator $`f_{LM}=f(r)Y_{LM}`$, one obtains the form for the corresponding transition density,

$$\rho _t^{coll}(r)=\frac{\mathrm{}^2}{2m}\sqrt{\frac{2L+1}{EWSR(f_{LM})E_{coll}}}\left[\left(\frac{1}{r}\frac{d^2}{dr^2}(rf)\frac{L(L+1)}{r^2}f\right)\rho _0+\frac{df}{dr}\frac{d\rho _0}{dr}\right].$$
(10)
Let us consider scattering operators, Eq. (4), with

$$f(𝒓)=f(r)Y_{1M}(\mathrm{\Omega }),f_1(𝒓)=rY_{1M}(\mathrm{\Omega }),$$
(11)
and write $`{\displaystyle \frac{1}{\pi }}`$Im$`G`$ as the sum of separable terms

$$R(𝒓^{},𝒓,E)=\frac{1}{\pi }\text{Im}G(𝒓^{},𝒓,E)=\underset{n}{}d_n(E)\rho _n(𝒓)\rho _n(𝒓^{}).$$
(12)
Note that $`d_n(E)`$ accounts for the energy dependence of $`R(𝒓^{},𝒓,E)`$. In the case of a well defined resonance, or in a discretized continuum calculation, the sum in Eq. (12) has only one term. In this case $`\rho _n`$ is proportional to the transition density associated with the resonance and may contain a spurious state contribution. In general, due to the smearing with $`\mathrm{\Gamma }/2`$, the sum in Eq. (12) may contain quite a few terms. We now write $`\rho _n`$ as

$$\rho _n(𝒓)=a_n\rho _{n3}(𝒓)+b_n\rho _{n1}(𝒓),$$
(13)
with

$$a_{n}^{}{}_{}{}^{2}+b_{n}^{}{}_{}{}^{2}=1.0.$$
(14)
Note that $`\rho _{n1}(𝒓)`$ is due to SSM and $`\rho _{n3}`$, associated with the ISGDR, fulfills the center of mass condition (for all $`n`$ )

$$f_1\rho _{n3}=f_1(𝒓)\rho _{n3}(𝒓)𝑑𝒓=0.$$
(15)
We point out that in the projection method for eliminating the effects of SSM, described in Ref. , it was assumed that $`a_n=1.0`$ (in (13)).

Following the derivation described in Ref. , we first note that all $`\rho _{n1}`$ coincide with the coherent spurious state transition density $`\rho _{ss}(𝒓)`$

$$\rho _{n1}(𝒓)=\rho _{ss}(𝒓)=\sqrt{\frac{\mathrm{}^2}{2m}\frac{4\pi }{AE_{ss}}}\frac{\rho _0}{r}Y_{1M}(\mathrm{\Omega }),$$
(16)
where $`E_{ss}`$ is the spurious state energy and $`\rho _0`$ is the ground state density of the nucleus. Note that $`\rho _{ss}`$ in (16) is normalized to 100% of the energy weighted sum rule (see (9) and (10)),

$$EWSR(rY_{1M})=\frac{\mathrm{}^2}{2m}\frac{3}{4\pi }A.$$
(17)
Looking for a projection operator that projects out $`\rho _{n1}(𝒓)`$,

$$F_\eta =\underset{i=1}{\overset{A}{}}f_\eta (𝒓_i)=F\eta F_1,$$
(18)
with $`f_\eta =f\eta f_1`$, we find that the value of $`\eta `$ associated with $`\rho _{ss}`$ is given by

$$\eta =f\rho _{ss}/f_1\rho _{ss}.$$
(19)
Using (15) and (19) we have

$$S_\eta (E)=f_\eta Rf_\eta =fR_{33}f,$$
(20)
where,

$$R_{33}=d_n(E)a_n^2\rho _{n3}(𝒓)\rho _{n3}(𝒓^{}).$$
(21)
To determine $`\rho _t`$ for the ISGDR we first use (6), (12), (13), (15) and (19) with $`F_\eta `$ and obtain

$$\rho _\eta (𝒓)=\frac{\mathrm{\Delta }E}{\sqrt{S_\eta (E)\mathrm{\Delta }E}}c_na_n[a_n\rho _{a3}(𝒓)+b_n\rho _{ss}(𝒓)],$$
(22)
with $`c_n=d_n(E)f_\eta \rho _{n3}`$. To project out the spurious term from (22) we make use of (15) and obtain

$$\rho _t(𝒓)=\rho _\eta (𝒓)b\rho _{ss},b=f_1\rho _\eta /f_1\rho _{ss}.$$
(23)
To properly normalize $`S_\eta (E)`$ and $`\rho _t`$, we have to determine the mixing amplitudes $`b_n`$ of the spurious state in the ISGDR. These amplitude can be obtained from the response function to the scattering operator $`f_1`$. Using (13), (15) and (16) we obtained from (12),

$$S_1(E)=f_1Rf_1=f_1R_{11}f_1=d_n(E)b_n^2f_1\rho _{ss}^2.$$
(24)
Note that $`f_1\rho _{ss}`$ can be obtained from the EWSR, Eq. (17),

$$f_1\rho _{ss}^2=\mathrm{}^22m\frac{3}{4\pi }A/E_{ss},$$
(25)
and the SSM probabilities from

$$b_n^2=\frac{S_1(E_n)}{f_1\rho _{ss}^2}.$$
(26)
In the present work we limit our discussion to the operator $`F_3=\underset{i=1}{\overset{A}{}}f_3(𝒓_i)`$, where $`f(𝒓)=f_3(𝒓)=r^3Y_{1M}(\mathrm{\Omega })`$. For this operator, the value of $`\eta `$ associated with the spurious transition density (16) is

$$\eta =\frac{5}{3}r^2,$$
(27)
and

$$S_\eta (E)=S_3(E)2\eta S_{13}(E)+\eta ^2S_1(E),$$
(28)
where $`S_3(E)=f_3Rf_3`$ is the strength function associated with $`f_3`$ and $`S_{13}=f_1Rf_3`$ is the non-diagonal strength function.

## III Results and discussions

In the following, we present our results for isoscalar giant resonances ($`L=0,1`$ and 2) obtained within the HF based RPA framework as briefly outlined in the previous section. Calculations are performed for $`{}_{}{}^{80}Zr`$ ($`N=Z=40`$). The two-body interaction $`V_{12}`$ is taken to be of a simplified Skyrme type,

$$V_{12}=\delta (\stackrel{}{r}_1\stackrel{}{r}_2)\left[t_0+\frac{1}{6}t_3\rho ^\alpha (\frac{\stackrel{}{r}_1+\stackrel{}{r}_2}{2})\right],$$
(29)
where $`\alpha =1/3`$, $`t_0=1800`$ MeVfm3 and $`t_3=12871`$ MeVfm3(α+1). For these values of the interaction parameters the nuclear matter equation of state has a minimum at $`E/A=15.99`$ MeV, $`\rho _0=0.157`$ fm-3 with $`K_{nm}=226`$ MeV, where $`E/A`$, $`\rho _0`$ and $`K_{nm}`$ being the binding energy per nucleon, matter saturation density and incompressibility coefficient for symmetric nuclear matter, respectively. This choice of the two-body interaction enables us to use the continuum RPA method to carry out a fully self-consistent calculation for giant resonances. Following Ref.  one can write the mean field potential $`V_{mf}`$ as,

$$V_{mf}=\frac{3}{4}t_0\rho (r)+\frac{\alpha +2}{16}t_3\rho ^{\alpha +1}(r)$$
(30)
and the particle-hole interaction $`V_{ph}`$ contributing to the isoscalar channel is given by

$$V_{ph}=\delta (\stackrel{}{r}_1\stackrel{}{r}_2)\left[\frac{3}{4}t_0+\frac{(\alpha +1)(\alpha +2)}{16}t_3\rho ^\alpha \right].$$
(31)
To begin with, we consider our results for isoscalar giant monopole, dipole and quadrupole resonances which are fully self-consistent and numerically accurate. Then, we shall analyze the influence of various numerical approximations on the centroid energies and transition densities for these resonances. Finally, we shall illustrate the possible effects of the violation of self-consistency on the properties of these isoscalar giant resonances (ISGR).

### A Self-consistent continuum RPA results

We now present our results of fully self-consistent HF$``$CRPA calculations for $`{}_{}{}^{80}Zr`$, using the Skyrme interaction of Eq. (29) with spin-orbit and Coulomb interactions switched off. It was pointed out in  that in order to have cancellations of the hole-hole transitions occurring in $`G_0`$ (Eq. (2)) and obtain numerically accurate results, it is important to employ the same mean-field and the same integration algorithm for the bound states and the single-particle Green’s function, using a small mesh size in double precision calculations. In the following we first present our results of highly accurate calculations obtained using $`dr_{HF}=0.04`$ fm and $`dr_{RPA}=0.04`$ fm, and with no smearing ($`\mathrm{\Gamma }=0`$ MeV), which we use in the following as a basis for comparison with other calculations. We note that in common implementations of HF-RPA one usually adopts the values of $`(dr_{HF},dr_{RPA})`$=(0.1 fm, 0.3 fm) and a smearing parameter of $`\mathrm{\Gamma }/21.0`$ MeV. In the following we use the notation $`dr=(dr_{HF},dr_{RPA})`$, with the values given in units of fm.

To facilitate our discussions we have displayed in Table I the HF single-particle energies for $`{}_{}{}^{80}Zr`$ obtained for $`dr_{HF}=0.04`$ fm. In Table II we give the values for the density radial moments $`r^2`$, $`r^4`$ and EWSRs (Eq. (9)) for various multipoles evaluated for different values of mesh size in the HF calculation. In Table III we present the values of energy weighted transition strengths (EWTS) for free and CRPA responses obtained using the operators $`f_3`$, $`f_1`$ and $`f_\eta `$ with $`dr=(0.04,0.04)`$ and $`\mathrm{\Gamma }=0`$ MeV. The quantities $`S_1^{EW}`$, $`S_3^{EW}`$, $`S_{13}^{EW}`$ and $`S_\eta ^{EW}`$ in Table III denote the EWTS for the corresponding strength functions $`S_1`$, $`S_3`$, $`S_{13}`$ and $`S_\eta `$, respectively, see Eq. (28). The transition strengths associated with sharp transitions were determined from the residues of the Green’s function, using its real part (see Eq. (8)). For the free response we get sharp peaks at the bound state single particle-hole transitions associated with $`L=1`$. These transitions can be easily identified from Table I as $`0g0f`$ (10.83), $`1d1p`$ (11.35), $`2s1p`$ (12.70), $`1d0f`$ (17.43), $`1d0p`$ (35.16), and $`2s0p`$ (36.52), with corresponding transition energies given in brackets in MeV. We checked that the values of the EWTS for these sharp transitions agree with the corresponding values obtained directly from the particle and hole wave functions.

For CRPA response, the sharp peaks occur below the particle threshold at 15.33 MeV. In addition to these sharp transitions, we have contributions from the continuum starting at the particle threshold. We obtained the contributions from the continuum by integrating the energy weighted strength function using small enough energy steps of $`dE=0.01`$ MeV. It is seen from Table III that the spurious state mixing is significantly larger for the free response (see 3rd column). Once the spurious state mixing is eliminated using the projection operator $`f_\eta `$ we find from the 2nd and 5th columns of this table that most of the strengths of the free response in the $`1\mathrm{}\omega `$ region of excitation energy ($`E<20`$ MeV) is spurious in nature. Only $`6.8\%`$ of the EWTS for the operator $`f_3`$ contributes to the intrinsic excitations for $`E<20`$ MeV. On the other hand, in the case of CRPA, since the calculation is fully self-consistent and numerically very accurate, the resonance occurring at 0.079 MeV is fully spurious and it exhausts $`99.99\%`$ of the EWSR associated with the operator $`f_1`$. For $`E>0.08`$ MeV, the values of $`S_1`$ and $`S_{13}`$ are so small that SSM is negligible. We see from the Table III that the values of $`b_n^2`$ is $`10^8`$ (see Eq.( 26)) indicating that the SSM is so small that one need not renormalize the strength $`S_\eta `$. For $`E>0.08`$ MeV, the values of the CRPA EWTS for the operators $`f_3`$ and $`f_\eta `$ are the same within $`1\%`$. We would like to emphasize that though the spurious state mixing is significantly large for the free response, it is fully eliminated by using the projection operator $`f_\eta `$ giving rise to $`99.95\%`$ of the expected EWSR which is quite close to the CRPA results. It may also be added that the fraction energy weighted sum rule, FEWSR = EWTS/EWSR, for the operator $`f_\eta `$ is $`8.4\%`$ and $`7.4\%`$ for $`E<20`$ MeV in the case of free and CRPA responses, respectively.

A proper test of a fully self-consistent calculation is to check how close $`\rho _t(r,E_{ss})`$ is to $`\rho _{ss}`$, where $`\rho _t(r,E_{ss})`$ is obtained from Eqs. (6) and (8) at the spurious state energy $`E_{ss}`$ using $`f_1`$. In Fig. 13 we compare the CRPA result for the $`\rho _t(r,E_{ss})`$ with the coherent state transition density calculated using Eq. (16). It is seen in Fig. 13 that in this highly accurate HF$``$CRPA calculation $`\rho _t(r,E_{ss})`$ coincides with $`\rho _{ss}`$, indicating a very negligible SSM.

We shall now present some plots for the strength functions for various multipoles obtained from our most accurate calculations. For plotting purpose we used a very small smearing width $`\mathrm{\Gamma }/2=0.025`$ MeV. In Figs. 13 and 13 we have shown the free and RPA response for the ISGDR, respectively. We see from Fig. 13 that most of the spurious components lie in the low energy region ($`E<20`$ MeV). As mentioned before, we see that the response for the operators $`r^3`$ and $`(r^3\eta r)`$ are indistinguishable in the case of a fully self-consistent HF based CRPA calculation. It also appears from these figures that particle-hole correlations do not alter the ISGDR strength distribution (shown in Fig. 13) very much which suggests that the isoscalar dipole state is not a very collective one. In Figs. 13 and 13 we have shown the plots for the ISGMR and ISGQR response functions, respectively. We have also carried out these calculations for $`dr=(0.24,0.24)`$ (not shown here) and find that they can not be distinguished from our most accurate calculations. We also note that at the surface the transition density for ISGMR looks like $`3\rho _0+rd\rho _0/dr`$ as given by Eq. (10), whereas, the ISGQR transition density looks more like $`d\rho _0/dr`$ rather than $`rd\rho _0/dr`$ as given by Eq. (10). We point out that Eq. (10) was derived under the assumption that one collective state exhausts the EWSR.

We have repeated the fully self-consistent calculations for $`\mathrm{\Gamma }=0`$ MeV, using various values of $`dr_{HF}`$ and $`dr_{RPA}`$. In Table IV we present CRPA results for the EWTS only for $`dr=(0.24,0.24)`$ and $`(0.04,0.24)`$ with $`N_{RPA}=50`$. We see from Table IV that the results for the operator $`f_3`$ for the different combinations of the mesh size differ by about $`2.5\%`$. The spurious state for $`dr=(0.24,0.24)`$ occurs at 0.7 MeV and its excitation energy becomes imaginary for $`dr=(0.04,0.24)`$. By multiplying the particle-hole interaction by a constant factor $`V_{sc}=0.9916`$ we push the spurious state to 0.1 MeV for $`dr=(0.04,0.24)`$ calculations. Nevertheless, we find that the SSM, or equivalently $`b_n^2`$ is very small ($`10^6`$). Once the spurious components are eliminated using the projection operator $`f_\eta `$, we get $`99.40\%`$ and $`99.76\%`$ of the expected EWSR for $`dr=(0.24,0.24)`$ and $`(0.04,0.24)`$, respectively. So far we have demonstrated that (i) as long as the calculation is fully self-consistent and numerically highly accurate, there is practically no spurious state mixing and (ii) the spurious state mixing introduced due to the use of a large mesh size (0.24 fm) in a CRPA calculation can be projected out using the operator $`f_\eta `$.

In Tables V and VI we have collected the centroid energies and the FEWSR, respectively, for the isoscalar resonances with $`L=0`$, 1 and 2 calculated using different combinations of the mesh size and a fixed value of $`\mathrm{\Gamma }/2=0.025`$ MeV. We notice that as long as the particle-hole interaction is not renormalized (i.e., $`V_{sc}=1.0`$) the centroid energies of the resonances do not deviate by more than 0.5$`\%`$ compared with the most accurate values. Though the energy of spurious state is sensitive to the values of the mesh sizes and increases from $`0.08`$ MeV to 0.71 MeV with the increase of radial mesh size from 0.04 fm to 0.24 fm, the centroid energy for ISGDR changes only by about 0.08 MeV. Even if $`V_{sc}`$ is used to shift the spurious peak to 2.0 MeV, the centroid energy for $`L=0`$ and $`L=1`$ resonance do not change appreciably. However, the centroid energy for the $`L=2`$ resonance goes up by about $`2\%`$ (0.3 MeV) . It is also clear from Table V and Figs. 13, 13 and 13 that the peak energy for ISGMR and ISGDR is higher than their centroid energies by about 0.5 and 0.15 MeV, respectively. From the Table VI we find that when $`dr=(0.04,0.24)`$, the values of the total EWTS for ISGMR and ISGQR are overestimated by $`12\%`$.

### B Influence of the smearing parameter $`\mathrm{\Gamma }`$

One of the requirements to avoid any SSM is that one must not use any smearing parameter (i.e., $`\mathrm{\Gamma }=0`$) and the calculations should be performed using a very fine mesh in the co-ordinate space while solving HF and RPA equations. However, one typically uses $`\mathrm{\Gamma }/21.0`$ MeV and the mesh $`dr0.1`$ fm. If the smearing width is finite, the spurious state would have a long energy tail which can give rise to large SSM. Because, $`\rho _{ss}d\rho _0/dr`$, which is a surface peaked functions, and has a large matrix element for the operator $`f_3`$. One must project out the SSM by making use of the projection operator $`f_\eta `$.

In Fig. 13 we plot CRPA results for the spurious state and ISGDR strength functions calculated using radial mesh size of 0.04 fm and smearing parameter $`\mathrm{\Gamma }/2=1`$ MeV. We clearly see from the figure that the strength function for the spurious state is extended up to a very high energy. The SSM caused due to the energy tail of the spurious state is eliminated using the operator $`f_\eta `$. In Table VII we give the values of FEWSR, associated with the scattering operator $`f_\eta `$, for the ISGDR for various energy ranges up to 150 MeV obtained using different values for the mesh size and the smearing parameter in the HF$``$CRPA calculation. Considering the values of the FEWSR in each energy range $`\omega _1\omega _2`$ of Table VII it can be easily seen that these values are practically the same as those obtained with $`\mathrm{\Gamma }=0`$, i.e., the SSM due to non-zero $`\mathrm{\Gamma }`$ is completely projected out. For $`\mathrm{\Gamma }/2=1.0`$ MeV the values for FEWSR for $`E=018`$ MeV is lower by about $`1\%`$ as compared to that for $`\mathrm{\Gamma }=0`$. This is because of the resonance at $`17.0`$ MeV (see Fig. 13). If we integrate the energy weighted strength for $`E=020`$ MeV, this difference reduces from $`1\%`$ to about $`0.5\%`$. We also note that for $`\mathrm{\Gamma }/2=1.0`$ MeV the total FEWSR obtained by integrating up to $`E=150`$ MeV is about $`1\%`$ lower than the one obtained for $`\mathrm{\Gamma }=0`$. Of course, this is because of the remaining strength beyond 150 MeV. For instance, in the case of $`dr=(0.24,0.24)`$ and $`\mathrm{\Gamma }/2=1.0`$ MeV we get FEWSR = 0.48$`\%`$ for the region for $`E=150300`$ MeV.

We point out that due to $`\mathrm{\Gamma }0`$, the transition density $`\rho _t`$ calculated using Eq. (6) depends on the scattering operator $`f`$. The consequences of this on the $`S(E)`$ and $`\rho _t`$ of the ISGDR was investigated and discussed in detail in Ref.  and we do not repeat it here. We thus demonstrated that using the projection scattering operator $`f_\eta `$ one can accurately eliminate the SSM effects occurring due to the use of a finite smearing parameter $`\mathrm{\Gamma }/2`$.

### C HF$``$DRPA results

The continuum can be discretized by confining the nucleus in a box of finite size. One can satisfactorily reproduce the continuum results, provided the calculations are carried out using a box of very large size (i.e., dense discretization) and the cut-off for the particle-hole excitation energy ($`E_{ph}^{max}`$) set to be reasonably high. We now consider our results obtained by discretizing the continuum using boxes of different sizes. The length of the box is given by $`N_{HF}`$ times $`dr_{HF}`$, where $`N_{HF}`$ is the number of radial mesh point used in a HF calculation. In the following, we present the results for discretized RPA calculations obtained using $`dr=(0.08,0.24)`$ with $`N_{HF}=150`$ and 900 (box sizes of 12 and 72 fm, respectively). In Figs. 13a, and 13b we show the ISGDR response for box sizes of 12 and 72 fm and smearing parameter $`\mathrm{\Gamma }/2=0.25`$ and 1.0 MeV, keeping $`E_{ph}^{max}=200`$ MeV, together with the corresponding results obtained in HF$``$CRPA. We see that the DRPA results obtained for the larger box coincide with the results obtained within CRPA. The transition strength gets fragmented if the discretization is carried out using a small box. To avoid misleading interpretation of the fragmentation and obtain agreement with the CRPA results, one needs to use a larger value of the smearing parameter, consistent with the size of the box. To examine more closely the effects of discretization on the response function we present in Table VIII our DRPA results for the FEWSR over various energy ranges up to 150 MeV. It is evident from this table that the total FEWSR increases significantly when $`E_{ph}^{max}`$ is increased from 50 MeV to 200 MeV. This increase is about $`56\%`$ and $`910\%`$ for $`\mathrm{\Gamma }=0.5`$ and 2 MeV, respectively. With a further increase in $`E_{ph}^{max}`$ there is no noticeable change in the value of the total FEWSR. It can be easily verified from this table that the FEWSR associated with the low-lying ISGDR component($`E<20`$ MeV) increases from $`6.4\%`$ to $`6.9\%`$ when $`E_{ph}^{max}`$ is increased from 50 to 200 MeV for the case of $`N_{HF}=150`$ and $`\mathrm{\Gamma }=0.5`$ MeV and it further increases to $`7.2\%`$ for $`N_{HF}=900`$ (see also Table VII). Comparing the Tables VII and VIII we can conclude (see also Figs. 13a and 13b) that with the proper choice of discretization and $`E_{ph}^{max}`$ one can mimic the continuum even for smaller values of $`\mathrm{\Gamma }0.51.0`$ MeV. Comparing the values of FEWSR, in each of the energy range $`\omega _1\omega _2`$, of Table VIII with Table VII we conclude that using $`f_\eta `$ one accurately eliminates SSM occurring due to the use of a low value for $`E_{ph}^{max}`$.

In Table IX we have displayed the values of $`E_{ss}`$ and the centroid energies for the $`L=0`$, 1 and 2 isoscalar giant resonances. These results are obtained using $`N_{HF}=900`$, $`\mathrm{\Gamma }/2=0.25`$ MeV with different values of $`E_{ph}^{max}`$ in HF-DRPA calculations. The corresponding HF$``$CRPA results are given in the last row of the table. We clearly see that as $`E_{ph}^{max}`$ increases, the centroid energies $`E_0`$, $`E_1`$ and $`E_2`$ converge to their corresponding exact values obtained using HF-CRPA. However, this convergence is slower for the spurious state energy $`E_{ss}`$. For low values of $`E_{ph}^{max}`$ we observe that the centroid energy for ISGMR is overestimated by 0.5 MeV, which can significantly effect the value of nuclear incompressibility. We also notice that $`E_1=35.3`$ MeV is little low for $`E_{ph}^{max}=50`$ MeV, because of the fact that the resonance energy for the ISGDR compressional mode is about 35.5 MeV (see also Fig. 13).

We saw in the previous subsection that the spurious transition density $`\rho _t(r,E_{ss})`$ obtained using a fully self-consistent CRPA calculation is indistinguishable from the corresponding collective model form for $`\rho _{ss}`$ which is proportional to $`d\rho _0/dr`$. In Fig. 13 we show some of the DRPA results for $`\rho _t(r,E_{ss})`$ and compare them with the $`\rho _{ss}`$. We see that for $`E_{ph}^{max}=50`$ MeV $`\rho _t(r,E_{ss})`$ deviates from $`d\rho _0/dr`$ even for the case of $`N_{HF}=900`$. However, for $`E_{ph}^{max}=200`$ MeV, the $`\rho _t(r,E_{ss})`$ from the DRPA is almost identical to the collective model results. Thus, one must use a reasonably large value for the cut-off energy, $`E_{ph}^{max}`$, in order to fully eliminate from the intrinsic excitations the contribution due to SSM.

### D Effects of violation of self-consistency

So far we have examined the various effects of numerical approximation on the properties of the isoscalar giant resonances of multipolarity $`L=02`$ and established the validity of the projection operator method in eliminating the SSM effects from the ISGDR. Here we report our investigations of the influence of certain violations of self-consistency on the strength function for isoscalar giant monopole ($`L=0`$), dipole ($`L=1`$) and quadrupole ($`L=2`$) resonances. These investigations are quite important in view of the fact that one often performs non self-consistent calculations for giant resonances such as the use of a phenomenological nuclear mean field (e.g., Woods-Saxon potential) and Landau-Migdal particle-hole interaction. Moreover, one often come across HF$``$RPA calculations carried out using particle-hole interaction not consistent with the mean field used in HF. We present below the results for HF based CRPA calculations carried out with the two-body interaction given in Eq. (29). We use the parameter $`V_{sc}`$ to renormalize the particle-hole interaction (i.e., $`t_0t_0V_{sc}`$ and $`t_3t_3V_{sc}`$ in Eq. (31)) so that the position of the spurious state can be adjusted close to zero. To study the consequences of the violation of self-consistency we vary $`t_0`$ and $`t_3`$ only in the particle-hole interaction (only in Eq. (31)).

In Table X we summarize our results for the centroid energies for isoscalar giant resonances for $`L=02`$. The quantity $`K_{nm}^{}`$ is the nuclear matter incompressibility coefficient associated with the renormalized parameters $`t_0V_{sc}`$ and $`t_3V_{sc}`$. Here, $`t_0`$ and $`t_3`$ are the values used in Eq. (31). Let us first consider the results obtained by varying $`t_0`$ by $`\pm 5\%`$ and $`\pm 10\%`$ and keeping $`t_3=12871`$ MeVfm4. It can be clearly seen from the table that the centroid energies for ISGMR and ISGDR significantly differ from their corresponding self-consistent values even if $`V_{sc}`$ is adjusted to give $`E_{ss}=0.1`$ MeV. On the other hand, the centroid energy for ISGQR reattains its self-consistent value when $`V_{sc}`$ is adjusted to yield $`E_{ss}=0.1`$ MeV. One may understand this discrepancy in terms of the incompressibility coefficient. With the renormalization of $`V_{ph}`$, though, $`E_{ss}`$ becomes close to zero, but values of $`K_{nm}^{}`$ in the RPA calculation remain quite different then the HF value of 226 MeV. In Fig. 13 we plot the values of $`E_0`$ and $`E_1`$ versus $`\sqrt{K_{nm}^{}}`$ for the cases with $`E_{ss}=0.1`$ MeV. This plot clearly depicts the systematic increase in $`E_0`$ and $`E_1`$ with increase in $`K_{nm}^{}`$. One may be tempted to infer at this point that as long as the nuclear matter incompressibility associated with the particle-hole interaction and the mean field is the same, centroid energies for the resonances considered here may come out to be reliable. In order to verify this, we adjust $`t_3`$ in particle-hole interaction in such a way that $`K_{nm}^{}`$ becomes 226 MeV when $`t_0`$ is varied by $`\pm 10\%`$. We see from Table X that even if $`K_{nm}^{}`$ is adjusted to 226 MeV, the values of $`E_0`$ and $`E_1`$ are off by about $`10\%`$ and $`3.5\%`$, respectively. This is due to the fact that the shape of the particle-hole interaction is not the same, though, the $`K_{nm}^{}`$ is kept constant. We note that if the ISGMR centroid energy is determined within $`10\%`$ accuracy, the value of nuclear matter incompressibility will be correct only within $`20\%`$.

Apart from the centroid energies for the giant resonances, it is also important to investigate the effects on the strength function itself when the self-consistency is not maintained. We looked into the plots for the strength functions $`S`$ and $`S_\eta `$ for the operators $`f_3`$ and $`f_\eta `$, respectively, for the different cases listed in Table X. We find that $`S_3S_\eta `$ or $`S_3<S_\eta `$, depending on the sign of interference between the spurious state and the intrinsic state (i.e., sign of the non-diagonal strength $`S_{13}`$). As an illustrative example, we show in Fig. 13 our results for the case in which $`t_0`$ is varied by $`10\%`$ and $`V_{sc}=1.7118`$. Similar is the case when $`t_0`$ is varied by $`5\%`$ and $`V_{sc}=1.2938`$. These values of $`V_{sc}`$ were chosen so that $`E_{ss}=0.1`$ MeV. In Figs. 13a, 13b and 13c we compare the fully self-consistent results for isoscalar giant resonances with those obtained by varying $`t_0`$ by $`\pm 5\%`$ in Eq. (31) and $`V_{sc}`$ is adjusted to yield $`E_{ss}=0.1`$ MeV. We see that the strength function for ISGMR and ISGDR are significantly different compared with their corresponding self-consistent results. Whereas, in case of ISGQR not only their centroid energies, but also the strength function seem to agree well with the corresponding self-consistent results. It is very important to point out that the violation of self-consistency causes redistribution of the strength in such a way that the total EWTS remains unaltered. This redistribution may be crucial in determining the energy weighted strengths associated with the low energy and the high-lying energy components of the ISGDR. For example, the fraction of the EWSR (in percent) for the energy range $`E=020`$ MeV is 6.94, 9.33 and 12.42 for $`t_0=1710`$, $`1800`$ and $`1890`$ MeV fm3, respectively, and for $`E=0150`$ MeV we have for the FEWSR = $`99.76\%`$ in these three cases.

We now focus on the influence of self-consistency violation when the continuum is discretized. As seen above, the discretization introduces two additional constraints, namely, the box size used in HF calculations and the maximum allowed particle-hole energy $`E_{ph}^{max}`$. We present here only the results for box size of 12 fm with $`E_{ph}^{max}=50`$ and $`200`$ MeV. In Fig. 13 we compare the ISGDR response function obtained for $`t_0=1620`$, -1800 and -1980 MeVfm3, keeping $`E_{ph}^{max}=50`$ MeV. Similar results are shown in Fig. 13 but obtained by raising $`E_{ph}^{max}`$ to 200 MeV. From $`S_3(E)`$ we see clearly that when the particle-hole interaction is in accordance with the mean-field potential, the SSM is only due to $`\mathrm{\Gamma }0`$. For the cases with $`t_01800`$ MeV fm3 one can immediately see a marked enhancement in spuriocity at $`E=1012`$ MeV. Furthermore, it is startling to see that the total FEWSR associated with operator $`f_\eta `$ for $`t_0=1620`$ and $`1980`$ MeVfm3 is 94.97$`\%`$ and $`58.97\%`$, respectively, and it is 95.13 $`\%`$ for $`t_0=1800`$ MeVfm3. For $`t_0=1710`$ MeVfm3 we get a total FEWSR = $`88.39\%`$ which is once again too much off compared to its expected value. We repeated the same analysis for box size 72 fm keeping $`E_{ph}^{max}=50`$ MeV but did not find any appreciable change in the values of the total FEWSR. However, when we raised the $`E_{ph}^{max}`$ from 50 to 200 MeV, we get the total FEWSR $`99.63\%`$, $`100.52\%`$ and $`99.94\%`$ for $`t_0=1620`$, $`1800`$ and $`1980`$ MeVfm3, respectively.

We also calculate the SSM probabilities (i.e., $`b_n^2`$) when self-consistency is not maintained. The values of $`b_n^2`$ are extracted using an extremely small smearing parameter. In case of $`t_0=1620`$ MeVfm3 and $`E_{ph}^{max}=50`$ MeV used in DRPA calculation, we find that $`E_{ss}=9.84`$ MeV. We get from Eq. (26), $`b_n^2=2.4\%`$ for the state occurring at $`13`$ MeV. When $`V_{sc}`$ is adjusted to push the spurious state energy $`E_{ss}`$ to about $`0.1`$ MeV, the EWTS of the 13 MeV state, associated with SSM, remains unchanged. Consequently $`b_n^2`$ reduces by two orders of magnitude. We thus conclude that since the values of $`b_n^2`$ are less than a few percent even with large violation of self-consistency, the renormalization of the strength function is not needed.

We have considered the effects on the ISGDR strength function when Coulomb/spin-orbit interaction is switched on in the HF calculation, but, ignoring it in the particle-hole interaction. We find that when spin-orbit interaction is included, the strength function obtained using $`\mathrm{\Gamma }/2=1`$ MeV hardly gets affected at any energy and the differences can not be seen on the plots (not shown here). This is due to the fact that the nucleus in question, $`{}_{}{}^{80}Zr`$, is spin saturated, i.e., the single-particle states with $`j=l\pm 1/2`$ are occupied. However, this may not be the case for non spin-saturated heavy nuclei. When we carried out similar exercise with the Coulomb interaction, the mean field changes significantly and we find that the strength functions gets shifted towards lower energy by about 2.0 MeV. We note that, with the inclusion of Coulomb interaction, the particle threshold for protons reduces from 15.33 MeV to 3.5 MeV.

## IV Conclusions

We have carried out self-consistent HF based CRPA calculations for isoscalar giant resonances with multipolarities $`L=0`$, 1 and 2 for $`{}_{}{}^{80}Zr`$ nucleus as an example. We demonstrate that if a self-consistent calculation is performed using zero smearing width and a very fine radial mesh size ($`dr=0.04`$ fm), the spurious state occurs at $`E_{ss}=0.08`$ MeV and the ISGDR response for the operators $`f_3`$ and $`f_\eta `$ are essentially the same for the energy $`E>E_{ss}`$ which indicates no SSM and the corresponding EWSR is reproduced remarkably well. When we use $`dr=0.24`$ fm in HF and CRPA calculations, $`E_{ss}`$ becomes about 0.7 MeV and there exists a small SSM. The amplitude of this SSM (i.e., $`b_n^2`$) is $`10^6`$, which is negligibly small and one need not renormalize the projected strength function. Although the position of the spurious state is quite sensitive to the radial mesh size and smearing parameter $`\mathrm{\Gamma }`$, the centroid energy for the isoscalar resonances for $`L=0`$, 1 and 2 do not change by more than $`0.5\%`$.

We have also performed the calculation for $`L=0`$, 1 and 2 isoscalar resonances by discretizing the continuum using boxes of different sizes (12 and 72 fm) with $`E_{ph}^{max}`$ ranging from $`50600`$ MeV. We found that the strength distribution is fragmented over a wide energy range for the case of the smaller box irrespective of $`E_{ph}^{max}`$. For the case of discretization in a large box (72 fm) with $`E_{ph}^{max}=200`$ MeV we find that the strength distribution agrees reasonably well with the corresponding one obtained from CRPA, if a moderate value of the smearing parameter ($`\mathrm{\Gamma }/21`$ MeV) is used. The spurious state occurs at about $`4.5`$ MeV for $`E_{ph}^{max}=50`$ MeV for both the small as well as large box discretization considered. With the increase of $`E_{ph}^{max}`$ to 600 MeV, we find that $`E_{ss}`$ approaches the corresponding value obtained within the CRPA. The centroid energies for $`L=0`$, 1 and 2 resonances converge to their corresponding exact values obtained from HF$``$CRPA. This convergence is little slow in case of spurious state energy. For $`E_{ph}^{max}=50`$ MeV, the transition density $`\rho _t(r,E_{ss})`$ at the spurious state energy obtained using discretized RPA differs from the corresponding CRPA results (which reproduce $`\rho _{ss}`$). However, with increase of $`E_{ph}^{max}`$ to 200 MeV, DRPA results for the spurious state transition density $`\rho _t(r,E_{ss})`$ become quite close to the CRPA results. We also point out that one should use $`E_{ph}^{max}200`$ MeV in order to calculate the centroid energies of the isoscalar $`L=0`$, 1 and 2 resonances with the accuracy of 0.1 MeV, comparable to the experimental uncertainties.

We have demonstrated that the spurious state mixing due to the non-zero smearing width and a choice of a coarse sized radial mesh can be accurately eliminated using the projection operator $`f_\eta `$. Furthermore, we show that the SSM due to a small value of $`E_{ph}^{max}`$ used in the DRPA calculation can be fully eliminated by applying the projection method.

We have investigated the consequences of violation of self-consistency on the $`S(E)`$ and $`\rho _t`$ of the isoscalar $`L=0`$,1 and 2 giant resonances by varying the parameter $`t_0`$ by $`\pm 5\%`$ and $`\pm 10\%`$ in the patrticle-hole interaction. We find that if the self-consistency is not maintained then the values of $`E_{ss}`$ and centroid energies for the $`L=0`$ , 1 and 2 isoscalar resonances are significantly different compared with their self-consistent values. Even if the particle-hole interaction is renormalized to shift $`E_{ss}`$ close to its self-consistent value, the centroid energies for $`L=0`$ and 1 resonances could not be corrected. This is due to the fact that though the renormalization corrects the value of $`E_{ss}`$, the nuclear matter incompressibility coefficient, $`K_{nm}^{}`$ associated with the particle-hole interaction is quite different than the one with the mean-field. However, the $`L=2`$ resonance is not very sensitive to the self-consistency violation as long as the particle-hole interaction is renormalized to shift $`E_{ss}`$ close to its self-consistent value. It is also important to point out that the violation of self-consistency causes a significant redistribution of the transition strength. In particular, energy weighted transition strength of the lower energy component ($`E<20`$ MeV) of the ISGDR response function may differ by $`50\%`$. The values of the SSM probabilities $`b_n^2`$ were found to be less than $`12\%`$. Therefore, one can neglect the renormalization of the ISGDR strength function obtained using projection operator $`f_\eta `$. Further, we found that the total energy weighted transition strength for the operator $`f_\eta `$ remains unaltered even with the violation of self-consistency.

Calculations were also carried out by changing the parameters appearing in the particle-hole interaction in such a way that the nuclear matter incompressibility coefficient associated with it remains unaltered. We find that though the incompressibility coefficient associated with the particle-hole and the mean field are kept the same, due to the lack of self-consistency, the centroid energy for $`L=0`$ and 1 isoscalar resonances are off by $`10\%`$ and $`3.5\%`$, respectively, compared to their self-consistent values. We may remark that if the ISGMR centroid energy is determined with accuracy of $`10\%`$, the value of $`K_{nm}`$ deduced from a comparison with experimental data is then accurate within only $`20\%`$.

This work was supported in part by the US Department of Energy under grant no. DOE-FG03-93ER40773.

$`{}_{}{}^{}V_{sc}^{}=0.9916`$ and $`E_{ss}=0.1`$ MeV.