text:
# Fermion zero modes in 𝑁=2 supervortices

## I Introduction

Stable magnetic flux tubes, or strings, occur in many high energy models. Fermionic zero modes on these strings can in some cases render them superconducting, and dramatically change their lifetime and interactions. In most of the examples that have been studied in a cosmological context the strings are topological, yet the effect of fermion zero modes on stable non–topological strings is equally powerful and interesting.

In this paper we investigate the fermionic zero modes on magnetic vortices in supersymmetric (SUSY) $`N=2`$ QED with two hypermultiplets of opposite charge. This model was proposed as a toy model for the low energy limit of type-II superstrings compactified on Calabi–Yau manifolds, where the possibility of magnetic vortices was first suggested in GMV96  (in the full low energy theory, there are sixteen hypermultiplets charged under fifteen U(1) groups, that is, fifteen copies of the model studied here).

In the absence of Fayet–Iliopoulos terms, it was proved in ARH98  that the strings appearing in that model are unstable. In ADPU01  the present authors studied the case where a (gauge) symmetry breaking Fayet-Iliopoulos (FI) term was added. In this case we found BPS string solutions. We also found that out of all the possible vacua, only certain choices can give rise to static strings, due to a vacuum selection effect PRTT96 . Even so, this effect was still not enough to obtain stable strings, unlike in the case of a comparable $`N=1`$ supersymmetric model with two chiral multiplets ADPU01 ; PU02 . The strings formed are semilocal strings VA91 , and in the BPS limit (equal scalar and gauge masses) they are only neutrally stable leese ; H92 .

In the present work, a study of the possible fermionic zero modes arising in this model has been carried out, in order to see whether they stabilise the string by further augmenting the vacuum selection effect. If this were so, the semilocal strings could possibly turn into chiral cosmic strings CP99 ; DKPS00 , long-lived vortices in which the fermions move in one direction only. However we have found that this is not the case. The fermions respect the same vacuum selection effect as the bosons, meaning that the underlying string remains semilocal. Moreover, the fermions move in both directions along the string, and therefore they could mix and leave the string. It is known that the fermion back-reaction can modify the stability of the background string N95 . We argue that this is not the case in the present context, because, unlike the systems studied in N95 , the Bogomol’nyi bound is protected by the remnant unbroken SUSY. The “topological charge” is identical for all the vortices in the family, and we do not expect any member of the family to be singled out by the fermion back-reaction.

Nevertheless, we have found that, due to the bosonic zero mode of the semilocal string, one of the fermions is coupled to a boson which is not zero at the core of the string, contrary to what might be expected from related assertions in the literature. To our knowledge this is the first time that the effect of charged, non–winding scalars has been considered on fermion zero modes, and it is interesting to point out that our model falls outside the scope of a number of index theorems for vortex backgrounds in scalar–fermion systems.

In section II we present a brief review of the bosonic sector of this model; we then obtain fermionic zero modes by SUSY transformations of the bosonic background in section III. We check that these in fact satisfy the fermionic equations of motion in section IV, and derive the mass matrix. The results are then discussed in section V.

## II Background

We will investigate an $`N=2`$ supersymmetric (SUSY) model consisting of two $`N=2`$ hypermultiplets $`(h_{ai},\psi _a,F_{ai})`$, labelled by $`a=1,2`$, with opposite charges $`q_a=\pm q`$, coupled to an $`N=2`$ abelian vector multiplet $`(A_\mu ,M,N,\lambda ^i,\stackrel{}{D})`$ (in Wess-Zumino gauge), $`i=1,2`$, plus a Fayet-Iliopoulos term $`\stackrel{}{k}\stackrel{}{D}`$, where $`\stackrel{}{k}`$ can be taken to be $`(0,0,k)`$ without loss of generality ADPU01 . This last term is responsible for breaking the $`U(1)`$ gauge symmetry, although it preserves the supersymmetry.

Using the conventions of sohnius , the Lagrangian we are interested in can be written as

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}D^\mu h_a^iD_\mu h_{ai}+i\overline{\psi }_a\gamma ^\mu D_\mu \psi _a+{\displaystyle \frac{i}{2}}\overline{\lambda }_i\gamma ^\mu _\mu \lambda ^i`$ (1)
$`{\displaystyle \frac{1}{4}}F^{\mu \nu }F_{\mu \nu }+i\widehat{q}_ah_a^i\overline{\lambda }_i\psi _ai\widehat{q}_a\overline{\psi }_a\lambda ^ih_{ai}`$
$`{\displaystyle \frac{1}{2}}\left(H_1^{\mathrm{\hspace{0.17em}1}}H_2^{\mathrm{\hspace{0.17em}2}}{\displaystyle \frac{1}{2}}\right)^2{\displaystyle \frac{1}{2}}\left(H_2^{\mathrm{\hspace{0.17em}1}}+H_1^{\mathrm{\hspace{0.17em}2}}\right)^2{\displaystyle \frac{1}{2}}\left(iH_2^{\mathrm{\hspace{0.17em}1}}iH_1^{\mathrm{\hspace{0.17em}2}}\right)^2,`$
where $`\widehat{q}_a=q_a/q`$, $`h_a^i=h_{aj}^{}`$, $`H_j^i=(\widehat{q}_a/2)h_a^ih_{aj}`$, $`D_\mu =_\mu +i\widehat{q}_aA_\mu `$ and $`A_\mu `$ is a U(1) gauge field. Auxiliary fields have been eliminated, all fields with zero vacuum expectation value (that is, $`M`$ and $`N`$), have been set to zero and suitable rescalings have been performed ADPU01 .

We can use a Bogomol’nyi argument to express the energy of the bosonic part of a static, straight vortex configuration as

$`E`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle }\mathrm{d}^2x\left[\right|(D_1+iD_2)h_{11}|^2+|(D_1iD_2)h_{12}|^2+|(D_1+iD_2)h_{21}|^2+|(D_1iD_2)h_{22}|^2+`$ (2)
$`[B+(H_1^{\mathrm{\hspace{0.17em}1}}H_2^{\mathrm{\hspace{0.17em}2}}+{\displaystyle \frac{1}{2}})]^2+(H_2^{\mathrm{\hspace{0.17em}1}}+H_1^{\mathrm{\hspace{0.17em}2}})^2+(iH_2^{\mathrm{\hspace{0.17em}1}}iH_1^{\mathrm{\hspace{0.17em}2}})^2]`$
$`{\displaystyle \frac{1}{2}}{\displaystyle \mathrm{d}^2xB},`$
where $`B=_1A_2_2A_1`$, and the integral in the last term is the magnetic flux, which is quantised in units of 2$`\pi `$. The Bogomol’nyi equations can then be immediately read off

$`(D_1+iD_2)h_{11}=0;(D_1iD_2)h_{12}=0;`$
$`(D_1+iD_2)h_{21}=0;(D_1iD_2)h_{22}=0;`$
$`H_1^{\mathrm{\hspace{0.17em}2}}=H_2^{\mathrm{\hspace{0.17em}1}}=0;B+H_1^{\mathrm{\hspace{0.17em}1}}H_2^{\mathrm{\hspace{0.17em}2}}+{\displaystyle \frac{1}{2}}=0`$ $`.`$
The only solutions to these equations have $`h_{12}=h_{21}=0`$ (this was dubbed a “vacuum selection effect” in ADPU01 ) ; the remaining equations are those of a semilocal string VA91 ; H92  in $`(h_{11},h_{22}^{})`$.

Up to global SU(2) transformations in the $`(h_{11},h_{22}^{})`$ space, the unit winding, cylindrically symmetric semilocal string solution can be expressed as

$`h_{11}`$ $`=`$ $`f(r)e^{i\theta };`$
$`h_{22}`$ $`=`$ $`g(r);`$
$`A_\theta `$ $`=`$ $`a(r).`$ (4)
with boundary conditions $`f(0)=a(0)=g^{}(0)=0`$ and $`f(\mathrm{})=1`$, $`g(\mathrm{})=0`$ and $`a(\mathrm{})=1`$. $`g(r)`$ is given by

$`g(r)=\kappa {\displaystyle \frac{f(r)}{r}}.`$ (5)
where the constant $`\kappa `$ essentially measures the width of the string, ranging from $`\kappa =0`$ –the Nielsen-Olesen NO73  string– to $`\kappa =\mathrm{}`$ –a $`CP^1`$ instanton H92 –. Note the lack of winding in $`h_{22}`$ and furthermore that $`h_{22}`$ does not necessarily have to be zero at $`r=0`$. Note also that $`B`$ is always negative with this choice of boundary condition, and its total flux does not depend on $`\kappa `$.

## III Fermionic zero modes

We are now interested in studying the fermionic zero mode solutions to this system to see whether they influence this vacuum selection effect. These solutions could be obtained by solving the explicit fermionic equations of motion, but we can also use the SUSY transformation to get the zero modes directly DDT97 . This will give us static configurations in the plane perpendicular to the string, and we will subsequently introduce $`t`$ and $`z`$ dependence on the solutions.

By zero mode solutions, we mean infinitesimal changes to the background configuration that preserve the action (for static configurations, that amounts to leaving the energy unchanged) and satisfy their equations of motion. We know that a SUSY transformation of a given configuration leaves the energy unchanged. Moreover, as we started with a static solution to the bosonic equations of motion, the fermions produced by this transformation must automatically satisfy their equations of motion.

The fermionic content of our system consists of two higgsinos (two Dirac fermions $`\psi _1`$ and $`\psi _2`$) coming from the hypermultiplets, and two gauginos (two symplectic Majorana fermions $`\lambda ^1`$ and $`\lambda ^2`$) coming from the gauge vector multiplet. Recall that the symplectic Majorana spinors are 4-component $`SU(2)`$ covariant spinors defined from 2-component spinors $`\lambda _{\alpha i}`$ and $`\overline{\lambda }_{\dot{\alpha }}^i(\lambda _{\alpha i})^{}`$ as

$`\lambda ^i\left(\begin{array}{c}i\epsilon ^{ij}\lambda _{\alpha j}\\ \overline{\lambda }^{\dot{\alpha }i}\end{array}\right);\overline{\lambda }_i=(\lambda _i^\alpha ,i\epsilon _{ij}\overline{\lambda }_{\dot{\alpha }}^j),`$ (6)
where $`\epsilon _{12}=\epsilon ^{12}=\epsilon _{\dot{1}\dot{2}}=\epsilon ^{\dot{1}\dot{2}}=1`$.

In $`N=1`$ SUSY language, the gauginos can be thought of as consisting of one gaugino coming from the $`N=1`$ gauge multiplet and one higgsino coming from a neutral (with respect to the gauge multiplet) chiral multiplet. The two SUSY generators ($`\xi _{\alpha (1)}`$, $`\xi _{\alpha (2)}`$) will be combined into two symplectic Majorana fermions $`ϵ^1`$, $`ϵ^2`$ (6).

Let us perform a SUSY transformation of the system in the background of the (bosonic) semilocal string obtained in the previous section. The bosonic fields do not transform. The higgsinos take the form

$`\delta \psi _{(a)}=i\gamma ^\mu D_\mu ϵ^{(i)}h_{ai},`$ (7)
while the gauginos may be written as

$`\delta \lambda ^{(i)}={\displaystyle \frac{i}{2}}\sigma ^{\mu \nu }ϵ^{(i)}F_{\mu \nu }iϵ^{(j)}\stackrel{}{\tau }_{(j)}^{(i)}\stackrel{}{D}.`$ (8)
Our conventions are sohnius

$`\gamma ^\mu =\left(\begin{array}{cc}0& \sigma ^\mu \\ \overline{\sigma }^\mu & 0\end{array}\right),{\displaystyle \frac{1}{2}}i[\gamma ^\mu ,\gamma ^\nu ]=\sigma ^{\mu \nu },`$ (9)
where, $`\sigma ^\mu =(1,\sigma )`$, $`\overline{\sigma }^\mu =(1,\sigma )`$ and $`\sigma ^i`$ are the Pauli matrices.

The higgsinos and the gauginos can be written more explicitly as

$`\delta \psi _{(1)}`$ $`=`$ $`{\displaystyle \frac{i}{2}}\left[\left(D_1iD_2\right)h_{11}\left(\gamma ^1+i\gamma ^2\right)ϵ^{(1)}+\left(D_1+iD_2\right)h_{12}\left(\gamma ^1i\gamma ^2\right)ϵ^{(2)}\right];`$
$`\delta \psi _{(2)}`$ $`=`$ $`{\displaystyle \frac{i}{2}}\left[\left(D_1+iD_2\right)h_{22}\left(\gamma ^1i\gamma ^2\right)ϵ^{(2)}+\left(D_1iD_2\right)h_{21}\left(\gamma ^1+i\gamma ^2\right)ϵ^{(1)}\right];`$
$`\delta \lambda ^{(1)}`$ $`=`$ $`\gamma ^1\gamma ^2ϵ^{(1)}Bi\left(H_1^{\mathrm{\hspace{0.17em}1}}H_2^{\mathrm{\hspace{0.17em}2}}+{\displaystyle \frac{1}{2}}\right)ϵ^{(1)};`$
$`\delta \lambda ^{(2)}`$ $`=`$ $`\gamma ^1\gamma ^2ϵ^{(2)}B+i\left(H_1^{\mathrm{\hspace{0.17em}1}}H_2^{\mathrm{\hspace{0.17em}2}}+{\displaystyle \frac{1}{2}}\right)ϵ^{(2)},`$ (10)
where we have used the Bogomol’nyi equations (II).

We expect the strings to be $`\frac{1}{2}`$BPS saturated WO78 , so let us try to obtain the broken and unbroken part of the SUSY transformation, i.e., let us obtain the fermionic zero modes. In order to do so, we can use the following projectors

$`P_\pm {\displaystyle \frac{1}{2}}\left(1\pm i\gamma ^1\gamma ^2\right),`$ (11)
which with our conventions, are given by

$`P_+\mathrm{diag}(1,0,1,0),P_{}\mathrm{diag}(0,1,0,1).`$ (12)
These projectors, besides $`P_\pm ^2=P_\pm ^{}=P_\pm `$ and $`P_\pm P_{}=0`$, have the following properties:

$`\gamma ^1P_\pm =P_{}\gamma ^1;`$
$`\gamma ^2P_\pm =P_{}\gamma ^2;`$
$`P_\pm \gamma ^1=\pm iP_\pm \gamma ^2.`$ (13)
$`i\gamma ^1\gamma ^2`$ is essentially a two-dimensional version of $`\gamma ^5`$, acting in the plane perpendicular to the string. Applying these projectors onto the fermions, we learn that

$`P_+\delta \psi _{(1)}`$ $``$ $`\delta \psi _{(1)+}=i\left(D_1iD_2\right)h_{11}\gamma ^1P_{}ϵ^{(1)}=2iD_1h_{11}\gamma ^1P_{}ϵ^{(1)};`$
$`P_{}\delta \psi _{(1)}`$ $``$ $`\delta \psi _{(1)}=i\left(D_1+iD_2\right)h_{12}\gamma ^1P_+ϵ^{(2)}=0;`$
$`P_+\delta \psi _{(2)}`$ $``$ $`\delta \psi _{(2)+}=i\left(D_1iD_2\right)h_{21}\gamma ^1P_{}ϵ^{(1)}=0;`$
$`P_{}\delta \psi _{(2)}`$ $``$ $`\delta \psi _{(2)}=i\left(D_1+iD_2\right)h_{22}\gamma ^1P_+ϵ^{(2)}=2iD_1h_{22}\gamma ^1P_+ϵ^{(2)},`$ (14)
and

$`P_+\delta \lambda ^{(1)}`$ $``$ $`\delta \lambda _+^{(1)}=i\left(B+H_1^{\mathrm{\hspace{0.17em}1}}H_2^{\mathrm{\hspace{0.17em}2}}+1\right)P_+ϵ^{(1)}=0;`$
$`P_{}\delta \lambda ^{(1)}`$ $``$ $`\delta \lambda _{}^{(1)}=i\left(BH_1^{\mathrm{\hspace{0.17em}1}}+H_2^{\mathrm{\hspace{0.17em}2}}1\right)P_{}ϵ^{(1)}=2iBP_{}ϵ^{(1)};`$
$`P_+\delta \lambda ^{(2)}`$ $``$ $`\delta \lambda _+^{(2)}=i\left(BH_1^{\mathrm{\hspace{0.17em}1}}+H_2^{\mathrm{\hspace{0.17em}2}}1\right)P_+ϵ^{(2)}=2iBP_+ϵ^{(2)};`$
$`P_{}\delta \lambda ^{(2)}`$ $``$ $`\delta \lambda _{}^{(2)}=i\left(B+H_1^{\mathrm{\hspace{0.17em}1}}H_2^{\mathrm{\hspace{0.17em}2}}+1\right)P_{}ϵ^{(2)}=0.`$ (15)
Note that $`\delta \psi _{(1)}`$ and $`\delta \psi _{(2)+}`$ vanish on BPS states due to the vacuum selection effect $`h_{12}=h_{21}=0`$.

It is clear that $`P_{}ϵ^{(1)}`$ and $`P_+ϵ^{(2)}`$ generate the fermionic zero modes, whereas $`P_+ϵ^{(1)}`$ and $`P_{}ϵ^{(2)}`$ are the generators of the unbroken SUSY. Note that there are fermionic zero modes of both chiralities, since N=2 SUSY is non–chiral.

## IV Fermionic equations of motion

In this section we investigate the structure of the fermionic zero modes by analysing the equations of motion directly, without reference to supersymmetry. We use two-spinor notation and define

$`\delta \psi _{(1)}\left(\begin{array}{c}\varphi _{\alpha (1)}\\ \overline{\chi }_{(1)}^{\dot{\alpha }}\end{array}\right)\delta \psi _{(2)}\left(\begin{array}{c}\varphi _{\alpha (2)}\\ \overline{\chi }_{(2)}^{\dot{\alpha }}\end{array}\right)\delta \lambda ^{(1)}\left(\begin{array}{c}i\mathrm{\Lambda }_{\alpha (2)}\\ \overline{\mathrm{\Lambda }}^{\dot{\alpha }(1)}\end{array}\right)\delta \lambda ^{(2)}\left(\begin{array}{c}i\mathrm{\Lambda }_{\alpha (1)}\\ \overline{\mathrm{\Lambda }}^{\dot{\alpha }(2)}\end{array}\right)`$ (16)
The projector operators can be defined in two spinor notation as

$`\sigma _+=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right);\sigma _{}=\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right).`$ (17)
In this notation the zero modes that we found previously are

$`\varphi _{\alpha (1)}`$ $`=`$ $`2iD_1h_{11}\sigma ^1\sigma _{}\overline{\xi }^{\dot{\alpha }(1)};`$
$`\overline{\chi }^{\dot{\alpha }(1)}`$ $`=`$ $`2D_1h_{11}\overline{\sigma }^1\sigma _{}\xi _{\alpha (2)};`$
$`\varphi _{\alpha (2)}`$ $`=`$ $`2iD_1h_{22}\sigma ^1\sigma _+\overline{\xi }^{\dot{\alpha }(2)};`$
$`\overline{\chi }^{\dot{\alpha }(2)}`$ $`=`$ $`2D_1h_{22}\overline{\sigma }^1\sigma _+\xi _{\alpha (1)};`$
$`\mathrm{\Lambda }_{\alpha (1)}`$ $`=`$ $`2iB\sigma _+\xi _{\alpha (1)};`$
$`\mathrm{\Lambda }_{\alpha (2)}`$ $`=`$ $`2iB\sigma _{}\overline{\xi }_{\alpha (2)},`$ (18)
where $`h_{11},h_{22}`$ and $`B`$ satisfy (II).

Consider the fermionic equations of motion derived from (1) in the bosonic background given by a solution to (II). Recall that $`h_{12}=h_{21}=0`$, and all other fields are independent of $`t`$ and $`z`$:

$`\left(\overline{\sigma }^\mu \right)D_\mu \varphi _{(1)}\overline{\mathrm{\Lambda }}^{(1)}h_{11}`$ $`=`$ $`0;`$
$`\left(\sigma ^\mu \right)D_\mu \overline{\chi }_{(1)}+i\mathrm{\Lambda }_{(2)}h_{11}`$ $`=`$ $`0;`$
$`\left(\overline{\sigma }^\mu \right)D_\mu \varphi _{(2)}+\overline{\mathrm{\Lambda }}^{(2)}h_{22}`$ $`=`$ $`0;`$
$`\left(\sigma ^\mu \right)D_\mu \overline{\chi }_{(2)}+i\mathrm{\Lambda }_{(1)}h_{22}`$ $`=`$ $`0;`$
$`\left(\sigma ^\mu \right)_\mu \overline{\mathrm{\Lambda }}^{(1)}+h_{11}^{}\varphi _{(1)}+i\chi _{(2)}h_{22}`$ $`=`$ $`0;`$
$`\left(\sigma ^\mu \right)_\mu \overline{\mathrm{\Lambda }}^{(2)}h_{22}^{}\varphi _{(2)}+i\chi _{(1)}h_{11}`$ $`=`$ $`0.`$ (19)
It can be checked that the static, $`z`$–independent configurations (18) satisfy these equations. Including the $`(t,z)`$ dependence into (19) we learn that three of the fermions are functions of $`(tz)`$ and therefore move in the positive $`z`$ direction, while the other three move in the negative $`z`$ direction:

$`\varphi _{\alpha (1)},\chi _{\alpha (2)},\mathrm{\Lambda }_{\alpha (1)}`$ $``$ $`t+z;`$
$`\varphi _{\alpha (2)},\chi _{\alpha (1)},\mathrm{\Lambda }_{\alpha (2)}`$ $``$ $`tz.`$ (20)
The fermions move in opposite directions, as expected, since $`N=2`$ SUSY is intrinsically non-chiral and it has not been broken in this model. We are interested in the $`(r,\theta )`$ dependence of the zero modes.

The usual Nielsen-Olesen string is one of the possible configurations in the family of semilocal strings, and it corresponds to the narrowest string. One can characterise this string as having $`h_{22}=0`$, annihilating two out of the six Weyl fermions $`\varphi _{(2)}=\overline{\chi }_{(2)}=0`$ hou . Moreover, if we remove one SUSY generator out of the two, we recover the Nielsen-Olesen string, but with chiral fermions DDT97  moving in the same direction. This is a good check of our results.

The situation is different when $`h_{22}0`$. The fermions $`\varphi _{(2)}`$ and $`\chi _{(2)}`$ are coupled to a field which is not zero at $`r=0`$, which might seem surprising. In terms of the general string ansatz (4), these fermions may be expressed in the form

$`\varphi _{(2)}`$ $`=`$ $`2i\kappa _r\left({\displaystyle \frac{f(r)}{r}}\right)e^{i\theta }\left(\begin{array}{c}0\\ \overline{\xi }^{\dot{1}(2)}\end{array}\right);`$
$`\chi _{(2)}`$ $`=`$ $`2\kappa _r\left({\displaystyle \frac{f(r)}{r}}\right)e^{i\theta }\left(\begin{array}{c}\overline{\xi }^{\dot{2}(1)}\\ 0\end{array}\right),`$ (21)
and, as can be seen from figure 1, they tend to $`0`$ at $`r=0`$. These two fermions are the only ones that wind.

In order to analyse the fermion zero modes, we convert the Dirac equations (19) into second order equations by acting with the operators $`\sigma D`$ and $`\overline{\sigma }D`$:

$`(\mathrm{}+\sigma ^3B+|h_{11}|^2)\varphi _{(1)}+ih_{11}h_{22}\chi _{(2)}(\sigma D)h_{11}\overline{\mathrm{\Lambda }}^{(1)}=0;`$
$`(\mathrm{}+\sigma ^3B+|h_{22}|^2)\chi _{(2)}ih_{11}^{}h_{22}^{}\varphi _{(1)}+i(\sigma D)h_{22}^{}\overline{\mathrm{\Lambda }}^{(1)}=0;`$
$`(\mathrm{}+|h_{11}|^2+|h_{22}|^2)\overline{\mathrm{\Lambda }}^{(1)}+(\overline{\sigma }D)h_{11}^{}\varphi _{(1)}+i(\overline{\sigma }D)h_{22}\chi _{(2)}=0;`$
$`(\mathrm{}\sigma ^3B+|h_{22}|^2)\varphi _{(2)}ih_{11}h_{22}\chi _{(1)}+(\sigma D)h_{22}\overline{\mathrm{\Lambda }}^{(2)}=0;`$
$`(\mathrm{}\sigma ^3B+|h_{11}|^2)\chi _{(1)}+ih_{11}^{}h_{22}^{}\varphi _{(2)}+i(\sigma D)h_{11}^{}\overline{\mathrm{\Lambda }}^{(2)}=0;`$
$`(\mathrm{}+|h_{11}|^2+|h_{22}|^2)\overline{\mathrm{\Lambda }}^{(2)}(\overline{\sigma }D)h_{22}^{}\varphi _{(2)}+i(\overline{\sigma }D)h_{11}\chi _{(1)}=0,`$ (22)
where $`\sigma D=\sigma ^1D_1+\sigma ^2D_2`$ and $`\overline{\sigma }^{1,2}=\sigma ^{1,2}`$.

We will use two dimensional projectors also in this case; and we define the projected two-spinors $`\mathrm{\Psi }_\pm `$ by $`\mathrm{\Psi }_\pm \sigma _\pm \mathrm{\Psi }`$, for any two-spinor $`\mathrm{\Psi }`$, where $`\sigma _\pm `$ are given in equation (17).

Half of the equations (22) correspond to the fermions coming from the unbroken SUSY symmetry:

$`(\mathrm{}B+|h_{11}|^2)\varphi _{(1)}+ih_{11}h_{22}\chi _{(2)}=0;`$
$`(\mathrm{}B+|h_{22}|^2)\chi _{(2)}ih_{11}^{}h_{22}^{}\varphi _{(1)}=0;`$
$`(\mathrm{}+|h_{11}|^2+|h_{22}|^2)\overline{\mathrm{\Lambda }}_+^{(1)}=0;`$
$`(\mathrm{}B+|h_{22}|^2)\varphi _{(2)+}ih_{11}h_{22}\chi _{(1)+}=0;`$
$`(\mathrm{}B+|h_{11}|^2)\chi _{(1)+}+ih_{11}^{}h_{22}^{}\varphi _{(2)+}=0;`$
$`(\mathrm{}+|h_{11}|^2+|h_{22}|^2)\overline{\mathrm{\Lambda }}_{}^{(2)}=0.`$ (23)
We can think of these equations as giving position–dependent masses (squared) for the fermions, in the sense that after diagonalisation, all six equations are of the form $`\mathrm{}\mathrm{\Psi }+M^\mathit{2}\mathrm{\Psi }=0`$, where

$`M^\mathit{2}=\mathrm{diag}(B+(|h_{11}|^2+|h_{22}|^2),B,|h_{11}|^2+|h_{22}|^2,B+(|h_{11}|^2+|h_{22}|^2),B,|h_{11}|^2+|h_{22}|^2).`$ (24)
Note that the matrix $`M^2`$ is not the same as the fermion mass matrix squared, due to the presence of derivative terms (in particular $`M^2`$ involves the magnetic field strength and is gauge invariant). Both matrices of course agree on constant bosonic backgrounds. As a simple check, far from the core we recover the fermion masses; indeed, $`B0`$, $`h_{11}1`$ and $`h_{22}0`$ as $`r\mathrm{}`$, and the diagonal terms become $`1,0,1`$, the correct masses for one higgsino, one goldstino and one gaugino of each chirality.

Recall that we are using signature $`(+,,,)`$. Since $`\mathrm{}=^2`$, the laplacian in the plane perpendicular to the string, (23) is a set of coupled time–independent Schrödinger equations in $`(r,\theta )`$ and we are looking for the zero energy states. Obviously, if the “potential” is everywhere positive there are no zero energy eigenstates. The functions in $`M^\mathit{2}`$ (24) are all greater than zero everywhere, and thus there are no non-zero normalizable solutions to (23), as can be shown by simply multiplying by $`\mathrm{\Psi }^{}`$ and integrating by parts.

Thus all of these fermionic components are automatically zero because of the equations of motion, and this agrees with the form of the solution obtained from SUSY transformations (19).

The other half of the equations reads

$`(\mathrm{}+B+|h_{11}|^2)\varphi _{(1)+}+ih_{11}h_{22}\chi _{(2)+}2\sigma ^2\overline{\mathrm{\Lambda }}_{}^{(1)}D_2h_{11}=0;`$
$`(\mathrm{}+B+|h_{22}|^2)\chi _{(2)+}ih_{11}^{}h_{22}^{}\varphi _{(1)+}+2i\sigma ^2\overline{\mathrm{\Lambda }}_{}^{(1)}D_2h_{22}^{}=0;`$
$`(\mathrm{}+|h_{11}|^2+|h_{22}|^2)\sigma ^2\overline{\mathrm{\Lambda }}_{}^{(1)}2\varphi _{(1)+}D_2h_{11}^{}2i\chi _{(2)+}D_2h_{22}=0;`$
$`(\mathrm{}+B+|h_{22}|^2)\varphi _{(2)}ih_{11}h_{22}\chi _{(1)}+2\sigma ^2\overline{\mathrm{\Lambda }}_+^{(2)}D_2h_{22}=0;`$
$`(\mathrm{}+B+|h_{11}|^2)\chi _{(1)}+ih_{11}^{}h_{22}^{}\varphi _{(2)}+2i\sigma ^2\overline{\mathrm{\Lambda }}_+^{(2)}D_2h_{11}^{}=0;`$
$`(\mathrm{}+|h_{11}|^2+|h_{22}|^2)\sigma ^2\overline{\mathrm{\Lambda }}_+^{(2)}+2\varphi _{(2)}D_2h_{22}^{}2i\chi _{(1)}D_2h_{11}=0,`$ (25)
where, again, the Bogomol’nyi equations (II) have been used extensively.

These may be written as two sets of three coupled equations<sup>1</sup><sup>1</sup>1Note that each set of equations correspond to fermions of one given chirality and coming from one of the supersymmetry generators and we again seek to diagonalize the resulting mass matrices.

$`\mathrm{}\left(\begin{array}{c}\varphi _{(1)+}\\ \chi _{(2)+}\\ \sigma ^2\overline{\mathrm{\Lambda }}_{}^{(1)}\end{array}\right)+\left(\begin{array}{ccc}B+|h_{11}|^2& ih_{11}h_{22}& 2D_2h_{11}\\ ih_{11}^{}h_{22}^{}& B+|h_{22}|^2& 2iD_2h_{22}^{}\\ 2D_2h_{11}^{}& 2iD_2h_{22}& |h_{11}|^2+|h_{22}|^2\end{array}\right)\left(\begin{array}{c}\varphi _{(1)+}\\ \chi _{(2)+}\\ \sigma ^2\overline{\mathrm{\Lambda }}_{}^{(1)}\end{array}\right)=0`$ (26)
$`\mathrm{}\left(\begin{array}{c}\varphi _{(2)}\\ \chi _{(1)}\\ \sigma ^2\overline{\mathrm{\Lambda }}^{(2)_+}\end{array}\right)+\left(\begin{array}{ccc}B+|h_{22}|^2& ih_{11}h_{22}& 2D_2h_{22}\\ ih_{11}^{}h_{22}^{}& B+|h_{11}|^2& 2iD_2h_{11}^{}\\ 2D_2h_{22}^{}& 2iD_2h_{11}& |h_{11}|^2+|h_{22}|^2\end{array}\right)\left(\begin{array}{c}\varphi _{(2)}\\ \chi _{(1)}\\ \sigma ^2\overline{\mathrm{\Lambda }}^{(2)_+}\end{array}\right)=0`$ (27)
For simplicity we just consider the first set of three equations. In the case of the Nielsen-Olesen string member of the semilocal family , for which $`\kappa =0`$, $`h_{22}=0`$, the spinors $`\chi _{(2)+}`$ and $`\varphi _{(2)}`$ do not couple to any spinors, and have mass squared B, which is negative in the core of the string, and zero at infinity. After diagonalisation, the other mass squared terms are

$`{\displaystyle \frac{1}{2}}(B+2|h_{11}|^2\pm \sqrt{B^2+16D_2h_{11}D_2h_{11}^{}})`$ (28)
which correspond to two (r-dependent) linear combinations of $`\varphi _{(1)+}`$ and $`\sigma ^2\overline{\mathrm{\Lambda }}_{}^{(1)}`$. At $`r=0`$, the signs of these masses are $`+,`$ respectively. At infinity these masses are 1, 1, and the mass states are simply the uncombined spinors. We can immediately see by the same reasoning as before that one combination of fermions is zero everywhere, since its mass squared is positive everywhere. The other one is a combination of the higgsino and the gaugino.

For the general case, $`\kappa 0`$, we have to diagonalize two sets of $`3\times 3`$ matrices. The signs of the “eigenvalues” of the $`M^2`$ matrices are $`+,,`$ at $`r=0`$ and $`+,0,+`$ at infinity, the “eigenvectors” being a combination of the three fermions. Once again, the “eigenvector” whose mass–squared function is always positive is zero, and thus we are left with two non-zero “eigenvectors”.

In all cases, the fermions at infinity have masses $`1,0,1,1,0,1`$, which agree with the masses of the fields $`h_{11},h_{12},A_\mu ,h_{22},h_{21},M+iN`$, as should be the case since supersymmetry is unbroken there.

## V Discussion

We have studied the fermion zero modes arising in SUSY $`N=2`$ QED with two hypermultiplets and a FI term. It is known that in such a model there is a vacuum selection effect in the bosonic sector, and only some of all the possible vacua form a family of neutrally stable vortices ADPU01 . We showed that the fermion zero modes do not improve the stability of the vortices. This is due to the fact that fermions obey the same vacuum selection effect as the vortices. On the other hand, as the SUSY is $`N=2`$, fermions with both chiralities are present in the vortex.

We do not expect the stability to be altered by the back-reaction of the fermions either. A recent detailed calculation of quantum corrections for SUSY kinks shows that the Bogomol’nyi bound is preserved G ), even if the mass receives corrections.

The model analysed in this paper contains a family of vortices that have the same topological charge, and all saturate the Bogomol’nyi bound. SUSY is half broken for all members of the family, and the unbroken SUSY will preserve the Bogomol’nyi bound in all cases, since the multiplets are shortened. Thus the BPS condition holds for all members of the family, and as they all have the same topological charge, the energy of all of them is also the same. Thus, no member of the family will be singled out by fermionic back-reaction.

The fermions in this system are also interesting because some are coupled to a scalar field which is not zero at the core of the vortex. It is sometimes stated in the literature that the reason why vortices support fermion zero modes in their core is because the fermion masses, which come from coupling to the scalars, are zero there. We have shown here that this heuristic argument is incorrect by calculating explicitly the zero modes that couple to the non–zero scalar at the core ($`h_{22}`$, if the string is in $`h_{11}`$).

It is interesting to try to relate the zero modes derived from supersymmetry to the zero modes we see in the equations of motion. The two fermionic zero modes obtained directly from supersymmetry are related to translational zero modes in the bosonic sector.

We can see from the second order equations of motion that in fact there are two more fermionic zero modes: In the bosonic sector, there are zero modes corresponding to changes in the parameter $`\kappa `$ of the semilocal string, expanding the string core. This bosonic zero mode corresponds to the fermionic zero modes that have negative mass at zero, and which are massless at infinity. One way to see this is to take the case where there is just one hypermultiplet $`(h_i,\psi ,F_i)`$. The string formed would then be an ordinary Nielsen-Olesen string, which would not possess the bosonic zero mode. Furthermore, the spinors $`\varphi _2`$ and $`\chi _2`$ would disappear, and the remaining eigenstates would have signs $`+,,+,`$ for their masses squared at zero and $`+,+,+,+`$ at infinity, so that only two of the eigenstates are non-zero. When we re-introduce the second hypermultiplet, the extra fermionic zero modes correspond to the extra bosonic zero mode associated with $`h_{22}`$.

It is also possible to consider supersymmetric transformations of the bosonic zero mode to give the fermionic zero mode. We perturb the field $`h_{22}`$, adjusting the other fields so that we retain the Bogomol’nyi equations (II), and hence keep the same energy. In the limit $`\kappa =0`$, the semilocal string has $`h_{22}=0`$, and an infinitesimal perturbation of $`h_{22}`$ of the form $`\delta h_{22}=\alpha h_{11(background)}/r`$ does not modify the other fields. Hence, a supersymmetry transformation of this zero mode gives fermionic zero modes in $`\varphi _{(2)}`$ and $`\chi _{(2)}`$ only.

At infinity the bosonic zero mode corresponding to changing the string width is a goldstone boson, and this agrees with what we have discovered from the fermionic equations of motion: that the eigenstates that have mass zero at infinity are pure $`\varphi _{(2)}`$ or $`\chi _{(2)}`$ in the $`\kappa =0`$ case.

For a general semilocal string $`(\kappa 0)`$, all the bosonic fields will be affected by a perturbation in the string width, and so the corresponding fermionic zero modes are combinations of each set of three spinors.

We note that the index theorem of Davis, Davis and Perkins DDP97  does not apply in the case of the semilocal string due to the fact that one of the fields doesn’t wind in the core – an assumption made in the derivation of the index theorem. Hence we can’t use it to ascertain the number of zero modes in this case. In the case where the index theorem does apply, namely the Nielsen–Olesen string (i.e. one hypermultiplet), our results agree with the index theorem. For the semilocal string there are two extra zero modes corresponding to the change of the string width; the fermion zero modes are those fermions corresponding to the bosonic zero mode. This has been constructed explicitly in the Nielsen–Olesen limit of the semilocal string. The model also falls outside the scope of the index theorem of Ganoulis and LazaridesGL88 , because $`h_{22}`$ is charged under $`U(1)`$ but goes to zero at infinity. Our results agree with Weinberg’s index theoremW81  (see also JR81 ), since the new zero modes have opposite chiralities and therefore do not affect the counting of net (left minus right) fermion zero modes. This is in agreement with the fact that the magnetic flux measured at infinity is the same for all semilocal strings, including the Nielsen–Olesen string.

## Acknowledgements

We thank F. Freire, T. Vachaspati and P. van Baal for very useful discussions. This work was partially supported by the ESF COSLAB programme (M.P.) and by grants AEN99-0315, FPA 2002-02037 and 9/UPV00172.310-14497/2002 (A.A., J.U.). J. Urrestilla and M. Pickles thank the University of Leiden for their hospitality.


format: text

text:
# Extension of the Random Phase Approximation including the self-consistent coupling to two-phonon contributions

## I Introduction

Recent interest in the spectroscopic factors for the removal of low-lying $`p_{1/2}`$ and $`p_{3/2}`$ strength from $`{}_{}{}^{16}\text{O}`$ points to a substantial discrepancy between experiment leus  and theory md94 ; rad94 ; pmd95 ; muetherO16 ; GeurtsO16 ; O16jast ; Fabr01 ; FaddRPAO16 ; muether02 ; udias ; GiustiG ; WimMarco . The contribution to the reduction of these spectroscopic factors due to short-range correlations is important and is uniformly calculated to be about 10% in $`{}_{}{}^{16}\text{O}`$ using different realistic interactions and theoretical approaches. The best agreement with the data is obtained when also long-range correlations are included in the calculations muetherO16 ; GeurtsO16 ; FaddRPAO16 ; muether02 . The results of Ref. FaddRPAO16  demonstrate a particularly strong correlation between the theoretically calculated excitation spectrum of $`{}_{}{}^{16}\text{O}`$ and the resulting one-hole spectral function, which describes the excitations of the residual $`A1`$ nucleus. The appearance of an additional $`p_{3/2}`$ fragment at small missing energy when a low-lying excited $`0^+`$ state is obtained is one such correlation FaddRPAO16 . Another one is the appearance of low-lying $`s_{1/2}`$ and $`d_{5/2}`$ strength when the first $`3^{}`$ state in $`{}_{}{}^{16}\text{O}`$ is correctly described GeurtsO16 ; FaddRPAO16 . Both features are in accordance with the experimental data leus .

The calculation of Ref. FaddRPAO16  employed a Random Phase Approximation (RPA) approach to describe the excitation spectrum of $`{}_{}{}^{16}\text{O}`$. The quality of the spectrum in this approach is, however, quite inadequate. Even with realistic $`G`$-matrices as residual particle-hole (ph) interaction one generates at most a few low-lying collective isoscalar states of negative parity ($`3^{}`$ and $`1^{}`$) and no low-lying isoscalar positive parity states are obtained czer . These considerations suggest that an adequate description of the sp fragmentation, at low energy, requires a better description of the excited states within the same framework of self-consistent Green’s function (SCGF).

This improvement of the SCGF method is also important because it can be applied to heavier nuclei. This is possible since the method describes the features of the low-energy spectrum in terms of interactions between a relatively small number of quasiparticles and collective excitations. These modes represent the experimentally observed excitations of the system and therefore they include both the effects of short- and long-range correlations. The fragmentation of the single-particle (sp) strength is thus self-consistently included. Early applications to $`{}_{}{}^{48}\text{Ca}`$ and $`{}_{}{}^{90}\text{Zr}`$ have been reported in Refs. brand88 ; brand90 . These calculations did not include self-consistency but nevertheless were able to yield a reasonable description of giant resonances and Gamow-Teller states as well as the low-lying collective states. Moreover a formalism based on the Bethe-Salpeter equation for ph excitations can be extended to include the effects of the continuum continuum . Although the SCGF method has only been applied in doubly closed-shell nuclei, it is feasible that increased computational power will allow applications to open-shell systems. An exploratory calculation for semi-magic isotopes of Ni and Sn was reported in JieTh . In the near future a huge amount of data will become available for unstable nuclei from radioactive beam facilities. For this reason it is important to develop techniques that are flexible enough to describe all the above effects.

In this paper, we report a first step toward an extension of the SCGF formalism in this direction. We chose to use $`{}_{}{}^{16}\text{O}`$ as a test system for the following reasons. First, a successful calculation for this nucleus will help resolving the aforementioned discrepancies for the spectroscopic factors. Second, this system has been studied by various approaches in the past. From the results of these calculations one can infer the relevant physical ingredients to be included in the SCGF method. Third, since the spectrum of $`{}_{}{}^{16}\text{O}`$ is rather complicated, meeting the challenge of describing this system will virtually guarantee good results when the same method is applied to heavier nuclei.

Shell-model calculations for this nucleus have been reported in Refs. zbg68 ; rw73 ; Haxton0+ ; wbm92  which indicate the particular importance of four-particle four-hole (4p4h) admixtures to the 0p0h ground state in generating the first excited $`0^+`$ state. Other positive parity states at low energy are dominated by 2p2h components, as indirectly confirmed by inelastic electron scattering O16ph86 . Some of this resulting physics was anticipated in terms of deformation effects which simulate these types of many-particle many-hole admixtures Brown0+ . A complementary point of view is given by the Interacting Boson Model of Refs. iach1 ; iach2 . There, the low-lying positive parity states are understood in terms of the coupling of $`3^{}`$ and $`1^{}`$ isoscalar states. One should note that, due to the predominant 1p1h nature of these excitations, the coupling of different phonons generates the npnh configurations relevant to this problem. In Refs. iach1 ; iach2  the $`3^{}`$ and $`1^{}`$ states were used as phenomenological boson degrees of freedom to generate the spectrum of the nucleus. The connection with the underlying fermionic description was indicated but not explored completely. This relation was considered in Ref. wim2ph . There it was argued that the microscopic ph interaction contains two-phonon exchange contributions which include the actually observed low-lying states of $`{}_{}{}^{16}\text{O}`$ themselves, thereby generating the correct number of low-lying states observed at low excitation energy.

In this work we begin to implement the physical ingredients proposed in Refs. iach1 ; iach2  by extending the RPA to include the coupling to two-phonon states. The basic idea of the present approach is illustrated in Fig. 1. This figure depicts the coupling of a ph state to two intermediate phonons that are described by the ph propagator itself. If the two intermediate phonons have been computed using the RPA equations, they will already provide a reasonable description of the low-lying collective isoscalar $`3^{}`$ and $`1^{}`$ excited states. These phonons can be sufficient to generate the quantum numbers of the most important positive parity states. Particularly relevant is the coupling of two $`3^{}`$ phonons, $`3^{}3^{}=0^+2^+4^+6^+`$, that represent some of the correlated 2p2h states of $`{}_{}{}^{16}\text{O}`$. In the framework of SCGF one employs dressed sp propagator in the construction of the microscopic phonons. This results in a Dressed RPA (DRPA) approach for the ph calculation and in its extension to the coupling to two-phonon states. In the present paper, we employ the self-consistent sp propagator of $`{}_{}{}^{16}\text{O}`$ computed in Refs. FaddRPA1 ; FaddRPAO16 . It should be noted that the incorporation of all the 4p4h effects in the present formalism requires a full four-phonon calculation. For this reason, a complete resolution in terms of a microscopic description of the spectrum may only be partially successful.

We note that there has been a tremendous progress in recent years in the microscopic description of p shell nuclei using Green’s Function Monte Carlo and no-core shell model methods GFMC ; NoCshell ; 3Nforce . A possible application of the no-core shell model to $`{}_{}{}^{16}\text{O}`$ would properly include such 4p4h effects. However the description of spectroscopic factors would still require the construction of effective operators to include the effects of short-range correlations on these quantities whereas these are automatically included in the SCGF method.

The paper is organized as follows: Sec. II introduces the formalism to account for two-phonon coupling in the calculation of the ph propagator. The approach employed here is based on a formalism first introduced by Baym and Kadanoff for the description of response functions in a many-body system at finite temperature bk61 ; baym62 ; BKBook . This framework provides a procedure to construct the effective interaction of the ph Bethe-Salpeter equation that generalizes the (D)RPA approach. This method is described in Sec. II.1 and the resulting equations in Sec. II.2. More technical details are left to the Appendix. Sec. III describes the results for the spectrum within the current approximation scheme. Sec. IV is devoted to a study of convergence properties related to the number of two-phonon configurations included, and the role of time-inversion diagrams. In Sec. IV.1 we discuss the possible appearance of instabilities, in particular for the $`0^+`$ state, that was also observed in Ref. FaddRPAO16 . Conclusions are drawn in Sec. V.

## II Extension of ph(D)RPA Formalism.

In this work the central quantity of interest is the two-time polarization propagator, whose Lehmann lehmann  representation is given by

$$\mathrm{\Pi }_{\alpha \beta ,\gamma \delta }(\omega )=\underset{n0}{}\frac{\left(𝒵_{\alpha \beta }^n\right)^{}𝒵_{\gamma \delta }^n}{\omega \epsilon _n^\pi +i\eta }\underset{n0}{}\frac{𝒵_{\beta \alpha }^n\left(𝒵_{\delta \gamma }^n\right)^{}}{\omega +\epsilon _n^\pi i\eta }.$$
(1)
In Eq. (1), the poles and residues contain the information on the response and excitation energies of the system with $`A`$ particles in terms of the quantities

$`𝒵_{\alpha \beta }^n`$ $`=`$ $`\mathrm{\Psi }_n^A|c_\alpha ^{}c_\beta |\mathrm{\Psi }_0^A,`$
$`\epsilon _n^\pi `$ $`=`$ $`E_n^AE_0^A,`$ (2)
where $`E_n^A`$ and $`|\mathrm{\Psi }_n^A`$ are the exact energies and eigenstates of the $`A`$-particle system, the subscript $`0`$ refers to the ground state and $`c_\alpha ^{}`$ ($`c_\alpha `$) is the creation (annihilation) operator of a particle in the state $`\alpha `$. For a clearer discussion of the formalism employed in this work, it is useful to first give a description of the relevant contributions in terms of Feynman diagrams.

### II.1 Diagrammatic contributions

The exact ph propagator (1) is a solution of the Bethe-Salpeter equation, depicted in Fig. 2 fetwa ; AAA . This equation can be written schematically as

$$\mathrm{\Pi }=\mathrm{\Pi }^f+\mathrm{\Pi }^fK^{(ph)}\mathrm{\Pi }$$
(3)
where $`\mathrm{\Pi }^f`$ represents the free propagation of a quasiparticle and a quasihole in the nuclear medium and the ph kernel $`K^{(ph)}`$ is, in general, a four-time quantity. According to the Baym-Kadanoff procedure, a solution for $`\mathrm{\Pi }`$ is obtained by first generating a self-consistent solution of the sp propagator using a given choice of the self-energy. From the functional derivative of this (self-consistent) self-energy with respect to the corresponding sp propagator, one then obtains the irreducible ph interaction $`K^{(ph)}`$ that generates the corresponding conserving approximation for the ph propagator (when used in the Bethe-Salpeter equation). The standard RPA approach is derived by applying this procedure to the the Hartree-Fock (HF) sp propagator and self-energy. This corresponds to approximating $`K^{(ph)}`$ with the bare interaction $`V`$ and employing bare (HF) sp propagators as external lines. Eq. (3) then generates the RPA series of ring diagrams shown in Fig. 3.

The extension of the RPA formalism, proposed and implemented by Brand et al. brand88 ; brand90 , was suggested by the observed fragmentation of the sp strength and the necessity to go beyond a lowest-order self-energy for a commensurate theoretical description. There, the Baym-Kadanoff procedure is applied to a second-order approximation for the self-energy. This yields contributions to the kernel $`K^{(ph)}`$ that automatically include all the terms which couple the ph states to the 2p2h ones, in accordance with the Pauli principle for the latter states. Thus, this formalism takes into account the mixing with 2p2h configurations in the construction of the ph propagator (the diagram of Fig. 4a gives an example). In the work of Brand et al. only a single-pole approximation to the self-consistent propagators was employed. It is the aim of the present work to take the effects of the sp fragmentation more completely into account and therefore a fully dressed sp propagator must be used. The one employed in the present work was obtained in Ref. FaddRPAO16  by means of a Faddeev expansion for the nuclear self-energy. When one applies the Baym-Kadanoff prescription to the latter self-energy, a large set of contributions to $`K^{(ph)}`$ is generated thesis . According to Ref wim2ph , one may expect that the most important of these terms involve the couplings to two ph phonons as depicted in Fig. 1. It is not difficult to see that different diagrams, similar to those in Fig. 1, can be obtained through Pauli exchange of the phonon’s external lines. In total there are sixteen such possible contributions, corresponding to all the possibilities of connecting two ph phonons to a ph state by means of a single interaction, both in the upper and lower part of the diagram. It is this approximation to the irreducible interaction that will be pursued in the present work.

An additional ingredient entering the 2p2h Extended RPA (ERPA) of Refs. brand88 ; brand90 ; EDRPAgeurts  requires further discussion. This involves diagrams similar to the one in Fig. 4b that are obtained from the 2p2h contribution by inverting the sense of propagation of either the incoming or the outgoing ph pair. These diagrams involve higher excitations (at least 3p3h, when combined in the expansion of Fig. 3) and are expected to give rather small contributions. Nevertheless, they represent corrections to the terms of the ph interactions that control the RPA correlations and they also add Pauli corrections to the RPA expansion of Fig. 3 at the 3p3h level. In Ref. brand90  it was found that they play a role in stabilizing some particular solutions.

The present calculation includes both the direct two-phonon contributions of the type depicted in Fig. 1 and the diagrams similar to the one in Fig. 4c. Obviously, all the 2p2h ERPA contributions are incorporated in this approach. Moreover, full two-phonon configurations are accouned for and a SCGF approach is applied. Thus, the present formalism is an extension of (and goes well beyond) the calculations of Refs. brand88 ; brand90 ; EDRPAgeurts ; norm1 . However, to avoid complications in the notation, we will still refer to it as “Extended RPA”, or as “two-phonon ERPA” whenever confusion may arise.

The Baym-Kadanoff procedure also generates other two-phonon contributions, for example those coupling a pp and a hh phonon. These could be mixed with the presently considered configurations by means of an all order expansion of the Faddeev-Yakubovsky type. Given our present knowledge, such a massive resummation of diagrams does not appear to be relevant for the understanding of the spectrum of $`{}_{}{}^{16}\text{O}`$. Such a study is in any case beyond the scope of the present paper.

### II.2 Two-phonon contributions to the ph propagator

The usual DRPA equations are obtained from Eq. (3) by choosing $`K_{\alpha \beta ,\gamma \delta }^{(ph)}=V_{\alpha \delta ,\beta \gamma }`$, as mentioned above, and by keeping dressed propagators as external lines. In this way, one is left with only two-time quantities and, after Fourier transformation, the DRPA equation becomes

$$\mathrm{\Pi }(\omega )_{\alpha \beta ,\gamma \delta }=\mathrm{\Pi }_{\alpha \beta ,\gamma \delta }^f(\omega )+\underset{\mu ,\rho ,\nu ,\sigma }{}\mathrm{\Pi }_{\alpha \beta ,\mu \rho }^f(\omega )V_{\mu \sigma ,\rho \nu }\mathrm{\Pi }_{\nu \sigma ,\gamma \delta }(\omega ),$$
(4)
where all the indices and summations are shown explicitly. In Eq. (4), the free polarization propagator $`\mathrm{\Pi }^f(\omega )`$ is also a two-time quantity, with the following Lehmann representation

$$\mathrm{\Pi }_{\alpha \beta ,\gamma \delta }^f(\omega )=\underset{n,k}{}\frac{\left(𝒳_\alpha ^n𝒴_\beta ^k\right)^{}𝒳_\gamma ^n𝒴_\delta ^k}{\omega \left(\epsilon _n^+\epsilon _k^{}\right)+i\eta }+\underset{k,n}{}\frac{𝒴_\alpha ^k𝒳_\beta ^n\left(𝒴_\gamma ^k𝒳_\delta ^n\right)^{}}{\omega +\left(\epsilon _n^+\epsilon _k^{}\right)i\eta },$$
(5)
where $`𝒳_\alpha ^n=\mathrm{\Psi }_n^{A+1}|c_\alpha ^{}|\mathrm{\Psi }_0^A`$ ($`𝒴_\alpha ^k=\mathrm{\Psi }_k^{A1}|c_\alpha |\mathrm{\Psi }_0^A`$) are the spectroscopic amplitudes for the excited states of a system with $`A+1`$ ($`A1`$) particles and the associated poles $`\epsilon _n^+=E_n^{A+1}E_0^A`$ ($`\epsilon _k^{}=E_0^AE_k^{A1}`$) correspond to the excitation energies with respect to the $`A`$-body ground state. The indices $`n`$ and $`k`$ label the eigenstates of the systems with $`A\pm 1`$ particles and enumerate the fragments associated with the one-particle and one-hole excitations, respectively. When a dressed propagator is used as input, its one-body overlap functions and quasiparticle energies already contain information about the coupling of sp motion to 2p1h, 2h1p and more complex configurations. As a consequence, contributions beyond the 1p1h case are already included in $`\mathrm{\Pi }^f(\omega )`$.

Methods based on an RPA-like expansion produce an infinite series of diagrams in which the direction of propagation can be reversed from backward to forward and vice versa. In the case of standard (D)RPA, the interaction kernel is simply given the potential $`V`$ and it is the same for every contribution to the diagrammatic expansion. Therefore the usual (D)RPA equation can be written in compact form, as in Eq. (4). This is no longer true when one aims to include additional contributions and at the same time insists on working with two-time quantities. The two-phonon diagrams of Figs. 1 and 4c have different analytical expressions. As a consequence one first needs to separate all the four possible time directions of the kernel $`K^{(ph)}`$ –forward to forward, backward to backward, and the two time-inversion cases– before including the relevant diagrams in the Bethe-Salpeter equation (3). This can be achieved by splitting the free ph propagator (5) into its forward- and backward-going parts, denoted by \> and \< respectively,

$$\mathrm{\Pi }^f(\omega )\mathrm{\Pi }^{f>}(\omega )+\mathrm{\Pi }^{f<}(\omega ).$$
(6)
By performing this substitution in Eq. (3) one obtains a similar separation for the complete ph propagator

$$\mathrm{\Pi }(\omega )\mathrm{\Pi }^>(\omega )+\mathrm{\Pi }^<(\omega ),$$
(7)
where \> and \< now refer to the sense of propagation of the final lines only. Suppressing the indices and summations one obtains

$`\mathrm{\Pi }^>(\omega )`$ $`=`$ $`\mathrm{\Pi }^{f>}(\omega )+\mathrm{\Pi }^{f>}(\omega )K^{(ph)}\mathrm{\Pi }(\omega ),`$
$`\mathrm{\Pi }^<(\omega )`$ $`=`$ $`\mathrm{\Pi }^{f<}(\omega )+\mathrm{\Pi }^{f<}(\omega )K^{(ph)}\mathrm{\Pi }(\omega ).`$ (8)
The last step consists in substituting Eq. (7) into (8) and approximating each component of the ph kernel $`K^{(ph)}`$ with the sum of the bare interaction and the corresponding two-phonon contributions. The result corresponds to the ERPA equations given by

$`\mathrm{\Pi }^>(\omega )`$ $`=`$ $`\mathrm{\Pi }^{f>}(\omega )+\mathrm{\Pi }^{f>}(\omega )\left\{\left(V+W^>(\omega )\right)\mathrm{\Pi }^>(\omega )+\left(V+H^{>,<}\right)\mathrm{\Pi }^<(\omega )\right\}`$
$`\mathrm{\Pi }^<(\omega )`$ $`=`$ $`\mathrm{\Pi }^{f<}(\omega )+\mathrm{\Pi }^{f<}(\omega )\left\{\left(V+H^{<,>}\right)\mathrm{\Pi }^>(\omega )+\left(V+W^<(\omega )\right)\mathrm{\Pi }^<(\omega )\right\}.`$ (9)
In Eq. (9), $`W^>(\omega )`$ represents the contribution of all the sixteen two-phonon diagrams, Fig. 1, in the forward direction. $`W^<(\omega )`$ corresponds to the contributions connecting the backward-going terms. Analogously $`H^{>,<}`$ represent the sum of the diagrams in which a backward-going hp state is inverted in a forward-going ph one (Fig. 4c illustrates such a case) while $`H^{>,<}`$ is the time-reversed contribution.

The practical implementation of Eqs. (9) requires additional manipulation to treat the freely propagating lines in Figs. 1 and 4c. This situation is completely analogous to the one already discussed in Ref. FaddRPA1 ; BaDLund  for the 2p1h expansion of the self-energy. Specific details for the ERPA equation (9) are given in the sections in the Appendix. The next section reports on the application of this formalism to $`{}_{}{}^{16}\text{O}`$.

## III Results

As in the work of Ref. FaddRPAO16 , the ERPA equations where solved within a model space consisting of a finite set of harmonic oscillator states. All the first four major shells (from $`1s`$ to $`2p1f`$) plus the $`1g_{9/2}`$ where included to account for the sp orbitals that are most relevant for low-lying excitations. The harmonic oscillator parameter was chosen to be $`b`$=1.76 fm. As a consequence of the truncation of the model space a Brueckner $`G`$-matrix was used as a microscopic effective interaction. This $`G`$-matrix was derived from a Bonn-C potential bonnc  and computed according to Ref. CALGM .

The contributions of two-phonon states were first studied by solving the ERPA equation (9) with an independent-particle model (IPM) propagator. This propagator was constructed from a Slater determinant composed of the lowest occupied harmonic oscillator (h.o.) orbitals in the model space. Where the effects of nuclear fragmentation were included, the calculations employed the fully dressed sp propagator computed in Ref. FaddRPAO16 . This propagator contains no more than two main quasihole (quasiparticle) fragments for each state in the $`p`$ ($`sd`$) shell. For these orbits, the sp propagator can therefore be well represented by means of one or two principal fragments and an effective pole that accounts for strength far from the Fermi energy. For quasiparticle states associated with the $`pf`$ shell, instead, the fragmentation is more substantial thesis  and the self-consistent propagator was approximated by including the most important poles (up to 4 main fragments for the $`p_{3/2}`$) plus two effective ones that gather all the background strength. We have checked that the results are not sensitive to the details of the fragmentation of the $`1f`$ orbitals. For the $`2p`$ orbitals, we found that the strength distribution is sufficiently represented by the adopted prescription and further discuss the sensitivity to the results in Sec. IV.1.

In both cases (IMP and dressed inputs), the standard (D)RPA equation was solved first, in order to generate the ph phonons that entered the final two-phonon ERPA calculation. When coupling different phonons, only the lowest few states of $`{}_{}{}^{16}\text{O}`$ have the right quantum numbers, and energies low enough, to generate two-phonon contributions with unperturbed energy below 20 MeV. These are the only states relevant for our purposes wim2ph . In practice, the calculations were performed by including, in each channel, all two-phonon configurations up to 30 MeV. The stability tests reported in Sec. IV.3 demonstrate that, for the case of a dressed input propagator, higher two-phonon contributions do not change the results for the low-lying states appreciably.

### III.1 Results for the particle-hole propagator

The results obtained from the h.o. input propagator are displayed in Fig. 5. The isoscalar eigenvalues obtained for energies below 15 MeV are displayed, for both the standard RPA and for the ERPA calculation. The ERPA calculation produces lower energies for some of the states that were already obtained in RPA. Both $`3^{}`$ and $`0^+`$ shifted down to about 0.5 MeV above the experimental energy. As discussed further in Sec. IV.3, this is to be considered rather fortuitous and one must remember that the dressing of the sp propagator, which has the effect to raise these eigenvalues, has not been taken into account yet. Still, this result indicates that correlations between ph and two-phonon states can be substantial and go in the right direction in explaining the experimental results.

It is also worth considering the total ph spectral strength of these states obtained by summing the corresponding amplitudes, Eq. (2), as follows

$$Z_{n_\pi }=\underset{\alpha \beta }{}\left|𝒵_{\alpha \beta }^{n_\pi }\right|^2.$$
(10)
Results for both the excitation energy and $`Z_{n_\pi }`$ of the principal levels in Fig. 5 are given in Table 1. For the ERPA case, also the fraction of quasiparticle-quasihole and two-phonon configuration that appear in the wave function are shown. Note that $`Z_{n_\pi }`$ is substantially bigger than one for the $`3^{}`$ and $`0^+`$ states and that these values increase further for the ERPA results. This signals an increase of the collective character of these solutions which may lead to an instability of the RPA equations for interactions that are even more attractive. The two-phonon ERPA approach generates a triplet of states at about 14 MeV, with quantum numbers $`0^+`$, $`2^+`$ and $`4^+`$. However, these levels are almost exclusively composed of two-phonon contributions and contain only small admixtures of ph states, resulting in a small ph spectral strength $`Z_{n_\pi }`$. The quantum numbers and energies of these states indeed correspond closely to those obtained by coupling two $`3^{}`$ RPA phonons, each at 7.14 MeV. We note that a similar triplet is found experimentally at 12.05, 11.52 and 11.10 MeV, which also corresponds to twice the experimental energy of the first $`3^{}`$ phonon. The first experimental $`2^+`$ is found at lower energy and its spectral strength is known to have relevant ph components O16ph86 . Thus it cannot be identified with any of the above two-phonon contribution. For all the lowest states that are not already reproduced by standard RPA a very small ph spectral strength has been found, due to a general lack of mixing between the ph and the two-phonon configurations. Of interest is also the $`2^+`$ that represents the giant quadrupole resonance at 20.7 MeV. In this case RPA and ERPA give 22.9 and 23.3 MeV for the main peak but with a lower $`Z_{n_\pi }`$ in the second case. For this state, part of the ph strength (about 20%) is shifted to two-phonon configurations representing the expected fragmentation of the giant resonance.

Fig. 6 and Table 2 show the analogous results when dressed sp propagators are used as input. In general, the main effect of fragmentation is to screen the nuclear interaction as a consequence of the quenching of spectroscopic factors for the input sp propagator. For the DRPA case, this results in increasing the lowest $`3^{}`$ and $`1^{}`$ solutions by $``$2 MeV. More substantial is the effect on the lowest $`0^+`$ state with a predominantly ph character which rises to about 17 MeV, confirming the sensitivity of this state to details of the fragmentation and the strength of the nuclear interaction. Unlike the IPM case, we have chosen to solve the two-phonon ERPA by first shifting the lowest solution for both the $`0^+`$, $`3^{}`$ and $`1^{}`$ states down to their relative experimental energies. This has the advantage of lowering the most important two-phonon configurations and allows to investigate their interplay with the ph ones.

As can be seen in Fig. 6, the ERPA equations still generate a triplet of $`0^+`$, $`2^+`$ and $`4^+`$ states at twice the energy of the first $`3^{}`$ phonon. Due to the screening of the ph interaction, these states mix very little with the ph configurations yielding an almost degenerate triplet. Table 2 also gives a comparison between the total two-phonon content of the states and the individual contributions of the most important configurations. This decomposition demonstrates that this triplet is formed by pure $`3^{}3^{}`$ states. This observation is in accordance with the 2p2h character of these states O16ph86 ; O162p2h . Since the overall energy of these states approximately agrees with experiment, one can expect that this calculation correctly represents the bulk properties of their wavefunctions. Nevertheless, it is clear that an additional interaction is needed in order to split this triplet as observed experimentally. It should be noted that no solution that can be identified with the $`2^+`$ and $`4^+`$ levels at 9.8 and 10.4 MeV. The $`2^+`$ level has been interpreted in terms of an $`\alpha `$ particle rotating around an excited $`{}_{}{}^{12}\text{C}`$ core O16ph86  and therefore involves correlations that may go well beyond the present calculation. With regard to the $`2^+`$ strength around 20 MeV, we note that the sp fragmentation included in the DRPA equation already generates a distribution of $`2^+`$ strength. Table 2 contains a few of these solutions and shows that also for the dressed case a sizable mixing of ph and two-phonon contributions is obtained at this energy. This mixing generates the spreading of ph strength over different solutions, consistent with the finite width of such a resonant state. It should be noted, however, that a discrete basis is used here and therefore it is not possible to properly describe the continuous strength distribution of a giant resonance, for which a continuum or a complex basis should be used berggren .

Also interesting are the results for the other isoscalar $`0^+`$ states. The lowest solution with a predominant ph character, that was found in DRPA at $``$17 MeV, is obtained at a similar energy but is now characterized by a partial mixing with two-phonon configurations. The ph and two-phonon contents of this state given in Table 2 show that this is the result of mixing with the the lowest solution, which ends up at $``$11 MeV. The latter is predominately a two-phonon state. It is also seen that in both cases the relevant configuration comes from the coupling of two $`0^+`$ phonons. Configurations involving two $`3^{}`$’s make smaller contributions while the other modes not reported in the table were negligible. We note that the wavefunctions for these states contain several relevant ph configurations, obtained from different quasiparticle fragments in the $`pf`$ ($`sd`$) shells combined with quasihole fragments of the $`p`$ ($`s`$) shells. Therefore, the situation is more complicated than the simple picture of only two levels interacting with each other. We observe that the shell-model calculations for the first excited state of $`{}_{}{}^{16}\text{O}`$ give very small contributions from both 0$`\mathrm{}\omega `$ and 2$`\mathrm{}\omega `$ configurations and a strong population of 4$`\mathrm{}\omega `$ states Haxton0+ ; wbm92 . On the other hand, inelastic electron scattering experiments O16ph86  clearly excite this state identifying its partial ph character even though the state is dominated by 4p4h components. From the point of view of the SCGF approach, the one-body response is completely described by the (dressed) ph propagator, Eq. (4), and therefore the total one-body strength must be represented by $`Z_{n_\pi }`$ (10). This observation points to the necessity of a stronger mixing between ph and two-phonon configurations than obtained in the present calculations. Still, the present results indicate that the two-phonon configuration obtained by coupling the two lowest $`0^+`$ phonons plays a role in the structure of the first $`0^+`$ state itself. This means that, starting from the dressed ph admixture contained in it, a fully self-consistent calculation would also generate some contributions of 2p2h, 4p4h, and beyond. The calculations of Ref. iach1  have shown that the bulk of the 4p4h contributions to the first excited state of $`{}_{}{}^{16}\text{O}`$ may come from the coupling of four different phonons with negative parity ($`3^{}`$ and $`1^{}`$). However, the self-consistent role of coupling positive parity states was not conidered in that work. Both this effect, the RPA correlations and the inclusion of the nuclear fragmentation allow for the –at least partial– inclusion of configurations beyond the 2p2h case even when no more than two-phonon coupling is considered, as in this paper. A study of the importance of three- and four-phonon configurations within this approach is beyond the scope of the present work.

The principal negative parity isoscalar states are still shifted down in energy by ERPA over DRPA, as expected, but the improvement is less than 0.5 MeV. The low-lying $`3^{}`$ and $`1^{}`$ levels remain substantially above the experimental energy at 9.23 and 10.90 MeV, respectively, and more correlations will be needed in order to lower their energy. We note that these wavefunctions contain a two-phonon admixture obtained by coupling the low-lying $`0^+`$ excitation to either the $`3^{}`$ or $`1^{}`$ phonons. According to the above interpretation of the first $`0^+`$ excited state this corresponds to the inclusion of 3p3h and beyond. At the same time two additional negative parity solutions, with two-phonon character, are found at higher energy in accordance with experiment.

The ERPA results for an h.o. input give appropriate corrections to several energy levels. However, such corrections tend to become negligible in the successive calculations employing dressed propagators. Presumably this happens because a sizable part of the correlations, which in the latter case are introduced by ERPA, are already accounted for by the dressing of the sp propagators. In general, the inclusion of fragmentation has also the effect of screening the interaction between ph and two-phonon configurations. However, the results of Table 2 show that a considerable mixing between the two can be obtained if all the relevant configurations are sufficiently low in energy, although no substantial mixing between different two-phonons states is seen. According to these considerations, the biggest deficiency of the present approach is probably the lack of an interaction within the multi-phonon space. Refs. iach1 ; iach2  have shown that such correlations are due to pp and hh interactions, neglected here, and that their inclusion has important effects on the final result for the spectrum. These terms may considerably lower most of the eigenstates of the ERPA equations. However, a second type of correlation between multi-phonon states, that was not accounted for in Refs. iach1 ; iach2 , is the Pauli exchange between the fermion lines included in different phonon propagators. These will probably generate the opposite effect, screening part of the pp and hh correlations.

Although the results obtained with a dressed input propagator for the lowest excited states are not completely successful, it must be recognized that the present ERPA calculation does represent a step forward. Indeed a proper description of two-phonon contributions, at the present level, had not been accomplished before. These states not only account for configurations at the 2p2h and higher level but appear to have some relevance for the description of states at higher energy. In addition, excitations are generated with quantum numbers that can not be obtained in the standard (D)RPA approach at these low energies. Given the results of the present work, it appears that a calculation including both pp and hh correlations, a more complete treatment of Pauli effects and the most relevant configurations up to the four-phonon states can be achieved with the presnt-day computers. Such a calculation is planned for the future.

## IV Stability of the results and truncation of the model space

### IV.1 RPA instability

The RPA approach may break down for strongly attractive interactions, by genereating an excess of collectivity and lowering the first excited state below the ground-state energy. The inclusion of fragmentation of the sp strength can affect this behavior in several ways. First, both the screening of the nuclear interaction and the splitting of the particle and hole spectral strengths over different fragments tends to stabilize the DRPA equations. Second, correlations lower most of the sp energies in comparison with the IPM. This results in smaller unperturbed ph energies and therefore has the opposite effect of pushing the DRPA approach toward instability. Third, the instability can be sensitive to the details of the spectroscopic amplitudes. This last point can be best illustrated by looking at the sp strength with $`j^\pi =3/2^{}`$, to which the present calculation is most sensitive. For the model space employed in this work a given particle overlap function is a superposition of two h.o. eigenstates,

$$\psi _{p_{3/2}}^n(𝐱)=𝒳_{1p_{3/2}}^n\varphi _{1p_{3/2}}^{h.o.}(𝐱)\pm 𝒳_{2p_{3/2}}^n\varphi _{2p_{3/2}}^{h.o.}(𝐱),$$
(11)
where we have stressed the fact the two components can sum either constructively or destructively. The fragmentation introduces a non-zero $`𝒳_{1p_{3/2}}^n`$ component even for particle states that in the IPM would be described only by the $`2p_{3/2}`$ subshell. By considering Eqs. (25) one can see that the ph interaction for the DRPA approach can be substantially changed by this new component by virtue of the strong matrix elements of the effective interaction between $`1p_{3/2}`$ states. In particular, this also applies to the off-diagonal terms that drive the RPA-like correlations. It is also worth nothing that the sign of this correction changes accordingly to the relative sign of the components in Eq. (11), therefore pushing the states either towards stability or instability, respectively.

We have observed that choosing a wrong sign for one or more of the $`p_{3/2}`$ particle fragments has drastic effects for the solutions of the isoscalar $`0^+`$ channel, with the choice of constructive interference in Eq. (11) leading toward instability of the DRPA approach. For the input propagator employed in this work, all of the background fragments have a similar wavefunction and a negative or positive interference sign according to whether their sp energy is lower or higher than 20 MeV, respectively thesis . Therefore it is natural to collect such background distributions in two different poles, one for each region of missing energy. A similar behavior applies also to the other quasiparticle states in the $`2p`$ and $`sd`$ shells. We have checked that the results of Sec. III are stable with respect to the number of main fragments and effective poles included in each subshell, therefore no instability affects the present work. However, in applications in which the number of poles has to be strongly reduced it may not be possible to chose effective poles of the sp propagator that represent the details of the spectral distribution sufficiently accurately. In such a case an instability in the DRPA approach may arise. This situation occurred with a previous DRPA calculation of the isoscalar $`0^+`$ channel, reported in Ref. FaddRPAO16 , in which a smaller number of fragments was included. Thanks to improvements in the computer code that solves the DRPA problem, it is now possible to consider a sufficiently large number of poles in the sp propagator and this problem has been overcome.

### IV.2 Time-inversion diagrams

The contribution of two-phonon terms to the time inversion diagrams of Fig. 4 can, in principle, generate other correction to the interaction between ph and hp states that drives the RPA correlations. Nevertheless, the effect is very small due to the large energy denominators that appear in Eqs. (27) and (28). We have tested their influence by neglecting the corresponding terms $`H^{>,<}`$ and $`H^{<,>}`$ in Eq. (9) and report the results in Table 3. No appreciable change is generated illustrated by differences with respect to the results of Sec. III of at most 0.1%. It is worth noting that the contribution of $`H^{>,<}`$ and $`H^{<,>}`$ can in principle carry information on the Pauli breaking at the level of 3p3h and beyond. Therefore, they may become more important for the case of three-phonon calculations. Nevertheless, Table 3 suggest that they are not likely to play an important role in the description of the spectrum of $`{}_{}{}^{16}\text{O}`$.

### IV.3 Stability vs. number of two-phonon configurations

Fig. 7 shows the results for selected solutions of the ERPA equations, obtained by employing different sets of two-phonon states. For any given point, only those configuration with energy $`\epsilon _{n_a}^\pi +\epsilon _{n_b}^\pi E_{cut}`$ have been included in the calculation. For $`E_{cut}`$=70 MeV, about 700 to 1000 two-phonon contributions have been included, depending on the channel. This is one order of magnitude larger than the number of ph states that enter the calculation. Since no two-phonon configuration has energy lower than $`10`$ MeV, the leftmost points in Fig. 7 correspond to the simple (D)RPA calculation. As the lowest few two-phonon contributions are included, the solutions for these levels show a sudden jump for all the lowest excited states. After this, the results obtained by using an IPM input still continue to exhibit a dependence on the value of the cut-off, roughly lowering states by $`1`$ MeV every time that $`E_{cut}`$ increases by $`20`$ MeV. In particular, the $`0^+`$ and $`3^{}`$ solutions decrease progressively, eventually heading to RPA instability. This confirms that the values reported in Table 1 for these states do not correspond to converged results. The situation is much better for the case of a dressed input propagator, for which the low-lying solutions are approximately stable for $`E_{cut}>`$ 20 MeV. This confirms that corresponding higher-energy two-phonon excitations do not have strong influence as their effects have already been included by the dressing of the sp motion.

In Fig. 7 we also show the results for the most relevant of the $`2^+`$ solutions that represent the resonance at about $`20`$ MeV. Here the variation is more appreciable, in particular for values of the cut-off $`E_{cut}`$ comparable with the energy of the state itself. Also here a more stable behavior is obtained for the case of a dressed input propagator. We conclude that when the nuclear fragmentation is accounted for the choice of $`E_{cut}=30`$ MeV adopted in Sec. III is a adequate for considering the low-lying spectrum.

### IV.4 Size of the model space

As a second possible source of uncertainty one may consider the size of the employed model space. The choice of the model space employed in Sec. III allows to take into account excitations up to the $`2p1f`$ shells. This means that ph configurations up to $`2\mathrm{}\omega 27`$ MeV are included, consistently with the choice for truncating two-phonon states discussed above. The calculations of Ref. brand88 ; EDRPAgeurts  also show that this is adequate and suggest that better convergence is to be expected for a calculation with a dressed input. To check this for the present calculation, we have augmented the model space by adding two more h.o. shells, up to $`\rho =2n+\mathrm{}`$=5, and recomputed the $`G`$-matrix elements accordingly. The input sp propagator used in the calculations was obtained by including an extra particle pole for every subshell added to both the IPM’s Slater determinant and the fragmented dressed propagator. The sp energies for the new shells were chosen by solving the Brueckner-Hartree-Fock equation within the new model space, in an approach analogous to the one of Ref. GeurtsO16 .

The results are compared to the ones obtained for the smaller space in Table 3. For the IPM, a sizable variation of the $`3^{}`$ excitation energy is found. Also, the high values of the ph strength $`Z_{n_\pi }`$ for both the $`0^+`$ and $`3^{}`$ solutions indicate that these states are approaching instability. As a consequence of the different $`3^{}`$ excitation energy, which also enters the two-phonon calculation, most of the ERPA solutions are shifted down in energy. The main conclusions of Sec. III, however, remain unchanged. When dressed sp propagators are employed, no dramatic change occurs for the low-lying solutions and only some states with dominant ph character are shifted by no more that 1 MeV. This suggests that the main conclusions of Sec. III are also not affected by the truncation of the model space, whence the bulk of the missing correlations has to be looked for in an extension of the diagrammatic expansion. We note, however, that a more attractive G-matrix is expected to be obtained if the continuum outside the model space were modified by accounting for self-energy corrections. This feature may lead to a further improvement of the present results. In addition, one may explore the inclusion of effective three-body force as due to propagation outside the model space. Such effective forces are critical in obtaining the correct binding energy in no-core shell model calculations reported in Ref. 3Nforce .

## V Conclusions

The dressed RPA equations have been extended to account for the coupling of two ph phonons in forming the excited states of a many-body system. The coupling among ph and multi-phonon configurations is conceptually similar to the Interacting Boson Model of Ref. iach1  but limited, in this application, to only two-phonon admixtures. Nevertheless, the present approach has the added advantage of taking into account both the effects of nuclear fragmentation and the RPA-like correlations, as generated by two-phonon fluctuations in the ground state.

The resulting Extended RPA formalism has been applied to study the excitation spectrum of $`{}_{}{}^{16}\text{O}`$. The results suggests a sizable mixing of ph and two-phonon configurations for the case of the $`2^+`$ quadrupole resonance at 20.7 MeV. Other solutions, carrying quantum numbers that cannot be generated by the simple DRPA, are obtained with this method, among which a few isoscalar negative parity states and a triplet with $`J^\pi `$=$`0^+`$, $`2^+`$ and $`4^+`$ near 12 MeV, in accord with experiment. In particular, the states in the triplet were seen to be almost pure two-phonon configurations obtained by coupling two $`3^{}`$ phonons. This feature confirms the 2p2h character of these states. The results are less satisfying for the low-lying positive parity states, which are known to require a proper description of 4p4h excitations. The present approach predicts the lowest solution for the isoscalar $`0^+`$ channel to be at 11.3 MeV, considerably above the first experimental excited state. A sizable component of this solution is seen to be generated by the coupling of two of the lowest $`0^+`$ phonon themselves. This suggests that important contribution of 4p4h excitation may be included already at the two-phonon level in a self-consistent fashion. Other correlations beyond the 2p2h level are partially included in the present calculation through the dressing of the propagators and the RPA approach. As expected, though, the present implementation is not sufficient to obtain a complete description of the low-lying excitation spectrum of $`{}_{}{}^{16}\text{O}`$. The results obtained in this work show that an interaction between multi-phonon configuration is still missing. This can be achieved by including pp and hh correlation between different ph phonons. Also the inclusion of up to four-phonon states is expected to be relevant for this system.

Finally, the stability of the present results has been tested with respect to the truncation of the model space and to the number of two-phonon configurations accounted for. It was found that the effects of nuclear fragmentation acts to ‘renormalize’ the sp propagator, making the solutions of the ERPA fairly independent of higher energy configurations. As discussed in Sec. IV, this feature generates stable solutions with respect to the size of the model space and the corresponding $`G`$-matrix interaction used in this work.

The four-phonon calculation of Refs. iach1 ; iach2  give a very good description of the excitation spectrum of $`{}_{}{}^{16}\text{O}`$. Those findings and the effects of fragmentation discussed in this paper suggest that the present ERPA formalism can be suitably extended to generate a satisfactory description of this nucleus, within the framework of SCGF. Such extensions are presently under consideration. Such a calculation may also significantly reduce the discrepancy between the measured and the theoretical sp spectral function for this nucleus FaddRPAO16 .

###### Acknowledgements.

This work was supported in part by the U.S. National Science Foundation under Grants No. PHY-9900713 and PHY-0140316 and in part by the Natural Sciences and Engineering Research Council of Canada (NSERC).

## Appendix A Extended DRPA Equation with Two-phonon Contributions

When reducing the kernel of the Bethe-Salpeter Eq. (3) to a two-time quantity, one has to deal with the fact that some lines in the diagrammatic expansion continue to propagate unperturbed while some interaction occurs between other particles. This is also the case for the diagrams of Figs. 1 and 4. This situation can be overcome by redefining the objects that appear in Eq. (9) and promoting the quantum numbers $`\{n,k\}`$, –that label particle and hole fragments– to external indices. The usual form of the ph propagator, Eq. (1), is obtained only as a last step of the calculation by performing the sum over the $`\{n,k\}`$ indices EDRPAgeurts ; FaddRPA1 .

The separation of the propagators $`\mathrm{\Pi }^f(\omega )`$ and $`\mathrm{\Pi }(\omega )`$ into forward and backward components, outlined in Sec. II.2, is a natural consequence of adopting this prescription. Indeed particle and hole external lines turn into each other by time inversion and become quantities that depend on different quantum numbers, i.e. the fragmentation indices $`\{n,k\}`$. For this reason, Eqs. (6) and (7) are only formal relations and the arrow can be substituted by an equal sign only before summing over all the particle and hole fragments. The free polarization propagator (5) naturally splits in two components that are purely forward and backward-going. The relevant Lehmann representations are

$`\mathrm{\Pi }_{\alpha n_\alpha \beta k_\beta ,\gamma n_\gamma \delta k_\delta }^{f>}(\omega )`$ $`=`$ $`\delta _{n_\alpha ,n_\gamma }\delta _{k_\beta ,k_\delta }{\displaystyle \frac{\left(𝒳_\alpha ^{n_\alpha }𝒴_\beta ^{k_\beta }\right)^{}𝒳_\gamma ^{n_\alpha }𝒴_\delta ^{k_\beta }}{\omega \left(\epsilon _{n_\alpha }^+\epsilon _{k_\beta }^{}\right)+i\eta }}`$ (12)
$``$ $`G^{}{\displaystyle \frac{1}{\omega D+i\eta }}G,`$
$`\mathrm{\Pi }_{\alpha k_\alpha \beta n_\beta ,\gamma k_\gamma \delta n_\delta }^{f<}(\omega )`$ $`=`$ $`\delta _{k_\alpha ,k_\gamma }\delta _{n_\beta ,n_\delta }{\displaystyle \frac{𝒴_\alpha ^{k_\alpha }𝒳_\beta ^{n_\beta }\left(𝒴_\gamma ^{k_\alpha }𝒳_\delta ^{n_\beta }\right)^{}}{\omega +\left(\epsilon _{n_\beta }^+\epsilon _{k_\alpha }^{}\right)i\eta }}`$ (13)
$``$ $`(G^{})^{}{\displaystyle \frac{1}{\omega +Di\eta }}G^{},`$
where no summation is implied and a shorter notation for the unperturbed poles and residues has been introduced. The quantities $`D`$ and $`G`$ can be thought as matrices whose elements contain all the unperturbed poles and residues, respectively. Note that in this case $`D`$ is diagonal and depends on the fragmentation indices $`\{n,k\}`$ only, while $`G`$ is rectangular because its column indices depend on the model space orbitals $`\{\alpha \}`$ as well.

The separation of the complete propagator is a little more complicated. The splitting of Eqs. (7) and (8) involves the time direction of the outgoing lines only. Since the RPA series contains contributions that can invert several times the sense of propagation of the ph diagrams, both forward and backward poles can appear in each component

$`\mathrm{\Pi }_{\alpha n_\alpha \beta k_\beta ,\gamma \delta }^>(\omega )={\displaystyle \underset{n0}{}}{\displaystyle \frac{\left(𝒵_{\alpha n_\alpha \beta k_\beta }^{>n}\right)^{}𝒵_{\gamma \delta }^n}{\omega \epsilon _n^\pi +i\eta }}{\displaystyle \underset{n0}{}}{\displaystyle \frac{𝒵_{\beta k_\beta \alpha n_\alpha }^{<n}\left(𝒵_{\delta \gamma }^n\right)^{}}{\omega +\epsilon _n^\pi i\eta }},`$ (14)
$`\mathrm{\Pi }_{\alpha k_\alpha \beta n_\beta ,\gamma \delta }^<(\omega )={\displaystyle \underset{n0}{}}{\displaystyle \frac{\left(𝒵_{\alpha k_\alpha \beta n_\beta }^{<n}\right)^{}𝒵_{\gamma \delta }^n}{\omega \epsilon _n^\pi +i\eta }}{\displaystyle \underset{n0}{}}{\displaystyle \frac{𝒵_{\beta n_\beta \alpha k_\alpha }^{>n}\left(𝒵_{\delta \gamma }^n\right)^{}}{\omega +\epsilon _n^\pi i\eta }}.`$ (15)
In Eqs. (14) and (15) the spectroscopic amplitude splits in two contributions $`𝒵^>`$ and $`𝒵^<`$. These appears unchanged in both equations due to the time-inversion symmetries obeyed by the ph states. In terms of these definitions, and applying the summation prescription, Eq. (7) can be exactly formulated as follow

$`𝒵_{\alpha \beta }^n`$ $`=`$ $`{\displaystyle \underset{n_1,k_2}{}}\left[𝒵_{\alpha n_1\beta k_2}^{>n}+𝒵_{\alpha k_2\beta n_1}^{<n}\right],`$ (16)
$`\mathrm{\Pi }_{\alpha \beta ,\gamma \delta }^>(\omega )`$ $`=`$ $`{\displaystyle \underset{n_1,k_2}{}}\left[\mathrm{\Pi }_{\alpha n_1\beta k_2,\gamma \delta }^>(\omega )+\mathrm{\Pi }_{\alpha k_2\beta n_1,\gamma \delta }^<(\omega )\right].`$ (17)
Finally, the contributions of two-phonon diagrams in the forward and backward direction can be expressed in Lehmann representation as well

$`W_{\alpha n_\alpha \beta k_\beta ,\gamma n_\gamma \delta k_\delta }^>(\omega )`$ $`=`$ $`{\displaystyle \underset{n_a,n_b}{}}{\displaystyle \frac{\left(K_{\alpha n_\alpha \beta k_\beta }^{>n_an_b}\right)^{}K_{\gamma n_\gamma \delta k_\delta }^{>n_an_b}}{\omega \left(\epsilon _{n_a}^\pi +\epsilon _{n_b}^\pi \right)+i\eta }}`$ (18)
$`=`$ $`K_{}^{>}{}_{}{}^{}{\displaystyle \frac{1}{\omega E+i\eta }}K^>,`$
$`W_{\alpha k_\alpha \beta n_\beta ,\gamma k_\gamma \delta n_\delta }^<(\omega )`$ $`=`$ $`{\displaystyle \underset{n_a,n_b}{}}{\displaystyle \frac{K_{\alpha k_\alpha \beta n_\beta }^{<n_an_b}\left(K_{\gamma k_\gamma \delta n_\delta }^{<n_an_b}\right)^{}}{\omega +\left(\epsilon _{n_a}^\pi +\epsilon _{n_b}^\pi \right)i\eta }}`$ (19)
$`=`$ $`K_{}^{<}{}_{}{}^{}{\displaystyle \frac{1}{\omega Ei\eta }}K^<,`$
where $`\epsilon _{n_a}^\pi `$ and $`\epsilon _{n_b}^\pi `$ are the energies of the intermediate ph phonons.

### A.1 ERPA Matrix

The eigenvalue equation for the ERPA is obtained in the usual way, by substituting Eqs. (12) to (15) into Eq. (9) and then extracting the residues of the ph poles $`\epsilon _n^\pi `$. The result is an eigenvalue equation in terms of the vectors $`𝒵^>`$ and $`𝒵^<`$.

To linearize the problem, we introduce the following components

$`X_{n_\alpha k_\beta }^{(1)}`$ $``$ $`{\displaystyle \frac{1}{\omega D+i\eta }}G\left\{\left(V+W^>(\omega )\right)\left(𝒵^>\right)^{}\left(V+H^{>,<}\right)\left(𝒵^<\right)^{}\right\},`$
$`Y_{k_\alpha n_\beta }^{(1)}`$ $``$ $`{\displaystyle \frac{1}{\omega +Di\eta }}G^{}\left\{\left(V+H^{<,>}\right)\left(𝒵^>\right)^{}\left(V+W^<(\omega )\right)\left(𝒵^<\right)^{}\right\},`$
$`X_{n_an_b}^{(2)}`$ $``$ $`{\displaystyle \frac{1}{\omega E+i\eta }}K^>\left(𝒵^>\right)^{},`$ (20)
$`Y_{n_an_b}^{(2)}`$ $``$ $`{\displaystyle \frac{1}{\omega +Ei\eta }}K^<\left(𝒵^<\right)^{},`$
where $`X^{(1)}`$ and $`Y^{(1)}`$ represent the ph amplitudes that appear in the standard (D)RPA equations schuckbook ; brand90  and $`X^{(2)}`$ and $`Y^{(2)}`$ are the analogous two-phonon amplitudes introduced by the ERPA approach. The components $`X^{(1)}`$ and $`Y^{(1)}`$ are related to $`𝒵^>`$ and $`𝒵^<`$ respectively by

$`\left(𝒵_{\alpha n_\alpha \beta k_\beta }^>\right)^{}`$ $`=`$ $`𝒳_\alpha ^{n_\alpha }𝒴_\beta ^{k_\beta }X_{n_\alpha k_\beta }^{(1)}=G^{}X^{(1)}`$
$`\left(𝒵_{\alpha k_\alpha \beta n_\beta }^<\right)^{}`$ $`=`$ $`𝒴_\alpha ^{k_\alpha }𝒳_\beta ^{n_\beta }Y_{k_\alpha n_\beta }^{(1)}=(G^{})^{}Y^{(1)}`$ (21)
Eqs. (9) can be put in the form of a linear eigenvalue equation

$$\omega \left(\begin{array}{c}X^{(1)}\\ X^{(2)}\\ Y^{(1)}\\ Y^{(2)}\end{array}\right)=𝐌\left(\begin{array}{c}X^{(1)}\\ X^{(2)}\\ Y^{(1)}\\ Y^{(2)}\end{array}\right),$$
(22)
where the matrix $`𝐌`$ is defined as

$$𝐌=\left[\begin{array}{cccc}GVG^{}+D& GK_{}^{>}{}_{}{}^{}& G[V+H^{>,<}](G^{})^{}& \\ K^>G^{}& E& & \\ G^{}[V+H^{<,>}]G^{}& & G^{}V(G^{})^{}D& G^{}K_{}^{<}{}_{}{}^{}\\ & & K^<(G^{})^{}& E\end{array}\right].$$
(23)
The off-diagonal 2$`\times `$2 blocks in Eq. (23) describe diagrams in which the time direction of propagation is inverted. In the present case the only non vanishing elements are the ones that involve the inversion of a single ph state into a hp one, or vice versa. These correspond to the sum of the first-order term $`V`$, which represents the kernel of the bare RPA, and the more complex diagrams of Fig. 4 . Blank spaces would in principle allow to include more complicated contributions, that involve time inversion of two-phonon diagrams into ph or hp configurations. These contributions are not expected to play a relevant role for the present problem. It must be noted that if the terms involving the matrix $`GK^{}`$ are discarded in Eq. (23), the components $`X^{(1)}`$ and $`Y^{(1)}`$ decouple from $`X^{(2)}`$ and $`Y^{(2)}`$. In this case, Eqs. (9) reduce to the ph-DRPA one (4) and the matrix (23) would take the form of the standard RPA matrix schuckbook . The normalization condition is derived in the usual way, by extracting the contribution of order zero of the expansion around a given pole and by employing the conjugate equation. One eventually obtains brand90

$$\underset{i=1,2}{}\left(X_{}^{(i)}{}_{}{}^{}X^{(i)}Y_{}^{(i)}{}_{}{}^{}Y^{(i)}\right)=1,$$
(24)
where the inner product of the vectors $`X^{(i)}`$ and $`Y^{(i)}`$ is implied.

### A.2 Matrix Elements for ph ERPA

In the following we give the explicit expression for the matrix elements of Eq. (23). The contributions originating from the standard (D)RPA equation are (here and below, summations over repeated greek indices are understood)

$`\left(GVG^{}\right)_{n_1k_2,n_3k_4}`$ $`=`$ $`𝒳_\alpha ^{n_1}𝒴_\beta ^{k_2}V_{\alpha \nu ,\beta \mu }\left(𝒳_\mu ^{n_3}𝒴_\nu ^{k_4}\right)^{},`$
$`\left(G^{}V(G^{})^{}\right)_{k_1n_2,k_3n_4}`$ $`=`$ $`\left(𝒴_\alpha ^{k_1}𝒳_\beta ^{n_2}\right)^{}V_{\alpha \nu ,\beta \mu }𝒴_\mu ^{k_3}𝒳_\nu ^{n_4},`$
$`\left(G^{}VG^{}\right)_{k_1n_2,n_3k_4}`$ $`=`$ $`\left(𝒴_\alpha ^{k_1}𝒳_\beta ^{n_2}\right)^{}V_{\alpha \nu ,\beta \mu }\left(𝒳_\mu ^{n_3}𝒴_\nu ^{k_4}\right)^{},`$ (25)
$`\left(GV(G^{})^{}\right)_{n_1k_2,k_3n_4}`$ $`=`$ $`𝒳_\alpha ^{n_1}𝒴_\beta ^{k_2}V_{\alpha \nu ,\beta \mu }𝒴_\mu ^{k_3}𝒳_\nu ^{n_4},`$
with the corresponding unperturbed ph energies

$$D_{n_1k_2}=diag\{\epsilon _{n_1}^+\epsilon _{k_2}^{}\}.$$

The interaction between two-phonon intermediate states and the ph ones is given by

$`\left(K^>G^{}\right)_{n_a^\pi n_b^\pi ,n_1k_2}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{[\left(X_{}^{(1)}{}_{n_\mu k_\rho }{}^{n_a^\pi }X_{}^{(1)}{}_{n_1k_ϵ}{}^{n_b^\pi }\right)^{}𝒳_\mu ^{n_\mu }𝒴_\rho ^{k_\rho }V_{\mu \nu ,\rho ϵ}𝒴_ϵ^{k_ϵ}\left(𝒴_\nu ^{k_2}\right)^{}`$ (26)
$`\left(X_{}^{(1)}{}_{n_\mu k_\rho }{}^{n_a^\pi }X_{}^{(1)}{}_{n_\nu k_2}{}^{n_b^\pi }\right)^{}𝒳_\mu ^{n_\mu }𝒴_\rho ^{k_\rho }V_{\mu \nu ,\rho ϵ}\left(𝒳_ϵ^{n_1}\right)^{}𝒳_\nu ^{n_\nu }]`$
$`+[n_a^\pi n_b^\pi ]\}`$
where $`n_a^\pi `$ and $`n_b^\pi `$ are the quantum numbers of the two phonons that form the intermediate state. The quantities $`X_{}^{(1)}{}_{n_1k_2}{}^{n^\pi }`$ are the forward-going amplitudes of the intermediate phonons (20) or (21) and are obtained from the previous solution of the polarization propagator. Note that $`K^>G^{}`$ is symmetric under the exchange of the two indices $`n_a^\pi `$ and $`n_b^\pi `$, as required by the boson character of the ph phonons. The factor $`\frac{1}{2}`$ assures that no double counting happens when only two free propagators $`\mathrm{\Pi }^f(\omega )`$ are coupled. Due to the time-inversion properties of npnh states, the corresponding contribution for backward-going propagation is simply related to above one:

$$\left(K^<(G^{})^{}\right)_{n_a^\pi n_b^\pi ,k_1n_2}=\left\{\left(K^>G^{}\right)_{n_a^\pi n_b^\pi ,n_2k_1}\right\}^{}.$$

The analytical expression for the contribution of the time-inversion diagrams is a little more complicated. First, we introduce the following quantity that corresponds to the generation of a two-phonon state by a fluctuation of the nuclear mean field

$$U_{n_a^\pi n_b^\pi }=\frac{1}{\left(\epsilon _{n_a}^\pi +\epsilon _{n_b}^\pi \right)}[\left(𝒴_\sigma ^{k_\sigma }𝒴_\lambda ^{k_\lambda }\right)^{}V_{\sigma \lambda ,\mu \nu }\left(𝒳_\mu ^{n_\mu }𝒳_\nu ^{n_\nu }\right)^{}X_{}^{(1)}{}_{n_\mu k_\sigma }{}^{n_a^\pi }X_{}^{(1)}{}_{n_\nu k_\lambda }{}^{n_b^\pi }.]$$
(27)
This is also symmetric under the exchange of the indices $`n_a^\pi `$ and $`n_b^\pi `$. Then, the matrix elements containing the time-inversion diagrams, $`H^{>,<}`$ and $`H^{<,>}`$, can be written as follows (a summation over $`n_a^\pi `$ and $`n_b^\pi `$ is also implied)

$`\left(GH^{>,<}(G^{})^{}\right)_{n_1k_2,k_3n_4}`$ $`=`$ $`\left\{\left((G^{})H^{<,>}G^{}\right)_{k_3n_4,n_1k_2}\right\}^{}=`$ (28)
$`=`$ $`{\displaystyle \frac{1}{2}}\{\left[𝒳_\alpha ^{n_1}𝒴_\beta ^{k_2}V_{\alpha ϵ,\beta \rho }\left(𝒳_\rho ^{n_\rho }𝒴_ϵ^{k_ϵ}\right)^{}X_{}^{(1)}{}_{n_\rho k_3}{}^{n_a^\pi }X_{}^{(1)}{}_{n_4k_ϵ}{}^{n_b^\pi }\right]`$
$`+[(n_1n_4)(k_2k_3)+(n_1n_4,k_2k_3)]`$
$`+\left[𝒳_\alpha ^{n_1}𝒴_\gamma ^{k_3}V_{\alpha ϵ,\gamma \rho }\left(𝒳_\rho ^{n_\rho }𝒴_ϵ^{k_ϵ}\right)^{}X_{}^{(1)}{}_{n_\rho k_ϵ}{}^{n_a^\pi }X_{}^{(1)}{}_{n_4k_2}{}^{n_b^\pi }\right]`$
$`+[n_1n_4,k_2k_3]`$
$`+𝒳_\alpha ^{n_1}𝒳_\delta ^{n_4}V_{\alpha \delta ,\rho \lambda }\left(𝒳_\rho ^{n_\rho }𝒳_\lambda ^{n_\lambda }\right)^{}X_{}^{(1)}{}_{n_\rho k_2}{}^{n_a^\pi }X_{}^{(1)}{}_{n_\lambda k_3}{}^{n_b^\pi }`$
$`+𝒴_\rho ^{k_\rho }𝒴_\lambda ^{k_\lambda }V_{\rho \lambda ,\beta \gamma }\left(𝒴_\beta ^{k_2}𝒴_\gamma ^{k_3}\right)^{}X_{}^{(1)}{}_{n_1k_\rho }{}^{n_a^\pi }X_{}^{(1)}{}_{n_4k_\lambda }{}^{n_b^\pi }\left\}\right(U_{n_a^\pi n_b^\pi })^{}.`$
Finally, the two-phonon unperturbed energies are given by

$$E_{n_a^\pi n_b^\pi }=diag\{\epsilon _{n_a}^\pi +\epsilon _{n_b}^\pi \}.$$
(29)