# Anharmonic collective excitation in a solvable model

## I Introduction

The existence of multiphonon states in the nuclear spectrum of excitation has been predicted since the introduction of collective models . Examples of low-lying nuclear vibrational states have been known for many years in nuclear spectra and are still being actively investigated with new generation of detectors; in particular, two-phonon multiplets and some three-phonon states based on the low-lying collective quadrupole and octupole modes have been found . Recently, it has also been beautifully demonstrated that multi-phonon excitations of these low-lying collective vibrations strongly influence heavy-ion fusion reactions at energies near and below the Coulomb barrier , through the so-called fusion barrier distribution analysis. It was pointed out that anharmonicities of vibrational excitations can significantly alter the shape of fusion barrier distribution and that thus sub-barrier fusion reactions offer an alternative method for extracting the static quadrupole moments of phonon states in spherical nuclei.

In the past 15 years, evidence has been collected for two-phonon giant resonances as well. This evidence stems from heavy-ion reactions at intermediate energy, pion-induced double charge exchange reactions, and relativistic heavy-ion reactions via electromagnetic excitations, in particular the excitation of the double giant dipole resonance (DGDR). In the experiments of the last type, a puzzling problem has been reported. Although the experimental data show that the centroid of the DGDR is about twice the energy of the single phonon resonance, theoretical calculations which assume the harmonic oscillator for giant resonances considerably underestimate the cross sections for double phonon states. In connection with this problem, anharmonic properties of giant resonances are attracting much interests.

Recently Bertsch and Feldmeier applied the time-dependent variational principle (TDVP)  to large amplitude collective motions and discussed anharmonicities of various giant resonances. One of the advantages of their approach is that solving the resultant equations and estimating the degree of anharmonicity are quite simple. They found that the relative frequency of the double phonon state scales as $`\mathrm{\Delta }\omega /\omega A^{4/3}`$, $`A`$ being the atomic number of a nucleus.

Earlier, Ponomarev et al.  noted that this quantity scales as $`A^{2/3}`$ in the nuclear field theory (NFT) , the same as the NFT result for the octupole mode . Reference  also remarks that a liquid-drop estimate gives an A<sup>-4/3</sup> dependences, implying that quantal effects are responsible for the difference. Both NFT and TDVP are quantal theories giving different results in refs.  and , so differences must be due either to inadeqate approximations or differences in the underlying Hamiltonians. We therefore undertook to try both methods on a solvable Hamiltonian. This will test the reliability of both methods, and if both give correct results, the disagreement is very likely attributable to the Hamiltonian assumptions.

The time-dependent variational approach was recently applied also to dipole plasmon resonances of metal clusters . In Ref. , it was shown that the anharmonicity of dipole plasmon resonances is very small and scales as $`\mathrm{\Delta }\omega /\omega A^{4/3}`$, which is identical to the result of the time-dependent variational approach for nuclear giant resonances. On marked contrast, Catara et al. claimed a large anharmonicity of plasmon resonances using a boson expansion method, which contradicts both to the result of Ref.  and to recent experimental observation . This fact also motivated us to compare these two methods, together with the nuclear field theory, on a solvable model in order to clarify the origin of the discrepancy.

As an aside, we mention also that there is a size-dependent anharmonicity in quantum electrodynamics, considering the photon spectrum in a small cavity. In QED the only dimensional quantity is the electron mass $`m_e`$, and the photon-photon interaction is proportional to $`m_e^4`$. Thus on dimensional grounds the relative shift of a two photon state in a cavity of side $`L`$ scales as $`\mathrm{\Delta }\omega /\omega 1/m_e^4L^4`$. Considering that sizes of a system scale as $`RA^{1/3}`$, the results of refs.  are $`\mathrm{\Delta }\omega /\omega R^4`$, identical to QED.

In this work we compare the nuclear field theory, the time-dependent variation approach, and the boson expansion method on the second Lipkin model, eq. (2.5) in ref. . The model is finite-dimensional and can be solved exactly. It has widely been used in literature to test a number of many-body theories .

The paper is organized as follows. In Sec. II, we first solve the model in the random phase approximation (RPA) and discuss the harmonic limit. In Sect. III we derive the collective Hamiltonian using the TDVP. We requantize the collective Hamiltonian and discuss the deviation from the harmonic limit. Numerical calculations are performed and compared with the exact solutions. In Sec. IV, we use the nuclear field theory as an alternative. There we see that it gives the same result as the TDVP, to leading order in the dependence on the number of particles. In Sec. V, we compare those results with the boson expansion approach. It will be shown that it leads to the identical result to the TDVP and the NFT. Finally we summarise the paper in Sec. VI.

## II HARMONIC LIMIT

Lipkin, Meshkov and Glick proposed two Hamiltonian models to describe $`N`$ particles, each of which can be in two states, making $`2^N`$ states in the basis. Using Pauli matrix representation for the operators in the two-state space, the second model has the Hamiltonian

$$H=\frac{1}{2}ϵ\sigma _z\frac{V}{2}\sigma _x\sigma _x.$$
(1)   
The first term is the single-particle Hamiltonian with an excitation energy $`ϵ`$, and the second term is a two-body interaction. The quasi-spin operators $`\sigma _z`$ and $`\sigma _x`$ are given by

$`\sigma _z`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}(a_{1i}^{}a_{1i}a_{0i}^{}a_{0i}),`$ (2)   
$`\sigma _x`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}(a_{1i}^{}a_{0i}+a_{0i}^{}a_{1i}),`$ (3)   
respectively. $`a_{1i}^{}`$ ($`a_{1i}`$) and $`a_{0i}^{}`$ ($`a_{0i}`$) are the creation (annihilation) operators of the $`i`$-th particle for the upper and the lower levels, respectively. For small $`V`$, the Hartree ground state $`|0>`$ is the fully spin-polarised state with matrix elements given by

$$<0|\sigma _i|0>=N\delta _{i,z}.$$
(4)   
A suitable basis for the exact diagonalization of $`H`$ is the set of eigenvectors of the angular momentum operators $`J^2`$ and $`J_z`$ with $`J=N/2`$. Then the dimension of the matrix diagonalization is reduced from $`2^N`$ to $`N+1`$, making the numerical problem very easy.

Before going to the anharmonicity, we note that the harmonic limit is obtained by solving the RPA equations. This was carried out in ref.  for the first Lipkin model. The RPA frequency for the second Lipkin model, Eq. (1), is obtained in the same manner. Setting the RPA excitation operator $`O^{}`$ as

$$O^{}=X\left(\underset{i=1}{\overset{N}{}}a_{1i}^{}a_{0i}\right)Y\left(\underset{i=1}{\overset{N}{}}a_{0i}^{}a_{1i}\right),$$
(5)   
the RPA equations read

$$\left(\begin{array}{cc}A& B\\ B& A\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right)=\omega \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right),$$
(6)   
where $`A`$ and $`B`$ are given by

$`A`$ $`=`$ $`{\displaystyle \frac{1}{4N}}<0|[\sigma _{},[H,\sigma _+]]|0>=ϵ(1\chi ),`$ (7)   
$`B`$ $`=`$ $`{\displaystyle \frac{1}{4N}}<0|[\sigma _{},[H,\sigma _{}]]|0>=ϵ\chi ,`$ (8)   
respectively. We have defined $`\sigma _{}`$, $`\sigma _+`$, and $`\chi `$ as

$`\sigma _{}`$ $`=`$ $`2{\displaystyle \underset{i=1}{\overset{N}{}}}a_{0i}^{}a_{1i},`$ (9)   
$`\sigma _+`$ $`=`$ $`2{\displaystyle \underset{i=1}{\overset{N}{}}}a_{1i}^{}a_{0i},`$ (10)   
$`\chi `$ $`=`$ $`V(N1)/ϵ,`$ (11)   
respectively. From these equations, the RPA frequency and the amplitude of the forward and the backward scatterings are found to be

$`\omega `$ $`=`$ $`\sqrt{(A+B)(AB)}=ϵ\sqrt{12\chi },`$ (12)   
$`X`$ $`=`$ $`{\displaystyle \frac{\omega +ϵ}{2\sqrt{N\omega ϵ}}},`$ (13)   
$`Y`$ $`=`$ $`{\displaystyle \frac{\omega ϵ}{2\sqrt{N\omega ϵ}}},`$ (14)   
respectively.

Fig. 1 compares the exact solution for the excitation energy of the first excited state with the RPA frequency given by eq. (12). As a typical example, we set the strength of the interaction $`V/ϵ`$ to be 0.01. The solid line shows the exact solutions for this particular choice of the Hamiltonian parameters, while the dashed line shows the RPA frequency. At $`N=51,\chi =0.5`$, the RPA frequency becomes zero and the system undergoes phase transition from spherical to deformed.

Figure 2 is the same as Fig. 1, but for a fixed $`\chi `$. We set $`\chi `$ to be 0.25, which corresponds to the isoscalar giant quadrupole resonance. We find significant deviation of the RPA frequency from the exact solution for small values of $`N`$, suggesting large anharmonicities. We discuss now the deviation from the harmonic limit.

## III TIME-DEPENDENT VARIATIONAL APPROACH

The time-dependent variational approach has been applied to the first Lipkin model by Kan et al. in ref. , but has never been applied to our knowledge to the second model. In keeping with the procedure of ref. , we postulate a time-dependent wave function of the form

$$|\alpha \beta =\mathrm{exp}(i\alpha (t)\sigma _x)\mathrm{exp}(i\beta (t)\sigma _y)|0.$$
(15)   
The motivation for this ansatz appears in ref. . The operator in the first term is the one that we wish to evaluate in the transition matrix elements. The operator in the second term is obtained by the commutator with the Hamiltonian,

$$[H,\sigma _x]=iϵ\sigma _y.$$
(16)   
The Lagrangian is given by

$$=\dot{\alpha }\beta |\sigma _x|\beta \alpha \beta |H|\alpha \beta .$$
(17)   
We reduce this with the help of the identity

$$e^{i\sigma _i\theta }\sigma _je^{i\sigma _i\theta }=\mathrm{cos}2\theta \sigma _j+\mathrm{sin}2\theta \sigma (\widehat{i}\times \widehat{j}),$$
(18)   
where $`ij`$ are Cartesian indices of the Pauli matrices. For example, the bracket in the first term of eq. (17) is reduced as

$$\beta |\sigma _x|\beta =0|e^{i\sigma _y\beta }\sigma _xe^{i\sigma _y\beta }|0=\mathrm{cos}2\beta 0|\sigma _x|0\mathrm{sin}2\beta 0|\sigma _z|0.$$
(19)   
The first term in the Lagrangian is

$$\dot{\alpha }\beta |\sigma _x|\beta =N\dot{\alpha }\mathrm{sin}2\beta .$$
(20)   
The second term is

$$\alpha \beta |H|\alpha \beta =ϵ\frac{N}{2}\mathrm{cos}2\alpha \mathrm{cos}2\beta +V\frac{N}{2}(\mathrm{cos}^22\beta +N\mathrm{sin}^22\beta ).$$
(21)   
The Lagrangian may then be expressed as

$$=N\dot{\alpha }\mathrm{sin}2\beta +ϵ\frac{N}{2}\mathrm{cos}2\alpha \mathrm{cos}2\beta +V\frac{N}{2}(\mathrm{cos}^22\beta +N\mathrm{sin}^22\beta ).$$
(22)   
The first Lagrangian equation is $`d/dt/\dot{\alpha }/\alpha =0`$. It reduces to

$$\dot{\beta }=\frac{ϵ}{2}\mathrm{sin}2\alpha .$$
(23)   
Similarly from the second Lagrange equation, $`d/dt/\dot{\beta }/\beta =0`$ we obtain

$$\dot{\alpha }\mathrm{cos}2\beta +\frac{ϵ}{2}\mathrm{cos}2\alpha \mathrm{sin}2\beta (N1)V\mathrm{cos}2\beta \mathrm{sin}2\beta =0.$$
(24)   
Next let us linearize and see what the RPA frequencies would be. The linearized equations are

$$\dot{\beta }=ϵ\alpha $$
(25)   
$$\dot{\alpha }+(ϵ2V(N1))\beta =0.$$

The equation for the frequency reads

$$\omega ^2=ϵ^22(N1)Vϵ=ϵ^2(12\chi )$$
(26)   
in agreement with the result of eq. (12).

A Hamiltonian corresponding to our Lagrangian can be seen by inspection, comparing eq. (22) to the form

$$=\dot{q}p(p,q).$$
(27)   
Equation (22) is already of this form with e.g., $`p=\frac{N}{2}\mathrm{sin}2\beta `$, $`q=2\alpha `$. The Hamiltonian is then given by

$$(p,q)=ϵ\frac{N}{2}\mathrm{cos}q\sqrt{1(2p/N)^2}V\frac{N}{2}\left((1(2p/N)^2)+N(2p/N)^2\right).$$
(28)   
We now expand $``$ in powers of $`q`$ and $`p`$ up to fourth order. Dropping the constant term, the expansion has the form

$$=\frac{p^2}{2m}+\frac{k}{2}q^2+aq^4+bp^4+cq^2p^2$$
(29)   
with coefficients

$`{\displaystyle \frac{1}{m}}`$ $`=`$ $`{\displaystyle \frac{2}{N}}\left(ϵ2V(N1)\right),`$ (30)   
$`k`$ $`=`$ $`ϵ{\displaystyle \frac{N}{2}},`$ (31)   
$`a`$ $`=`$ $`{\displaystyle \frac{ϵN}{48}},`$ (32)   
$`b`$ $`=`$ $`{\displaystyle \frac{ϵ}{N^3}},`$ (33)   
$`c`$ $`=`$ $`{\displaystyle \frac{ϵ}{2N}}.`$ (34)   
Note that we recover the linear frequency, eq. (12), immediately from $`\omega ^2=k/m`$.

In ref. , the anharmonicity was determined by requantizing the Hamiltonian with the Bohr-Sommerfeld condition,

$$_{q0}^{q1}p𝑑q=n\pi ,$$
(35)   
where $`p`$ and $`q`$ satisfy $`(p,q)=E`$ and $`q_0`$ and $`q_1`$ are the endpoints of the motion at energy $`E`$. However, here we find it more convenient to use the equivalent formula

$$_{<E}𝑑p𝑑q=2n\pi .$$
(36)   
In the same sense as the expansion of the Hamiltonian $`(p,q)`$ as done in eq. (29), we apply eq. (36) iteratively. We first consider only the harmonic part of the Hamiltonian and transform the integration region to a circle. We also use polar coordinates and write $`p^{}=p/\sqrt{m}=r\mathrm{sin}\theta ,xq=\sqrt{k}q=\mathrm{cos}\theta `$. The radius of the circle is then $`r_0=\sqrt{2E}`$ and the harmonic approximation gives

$$_{<E}𝑑p𝑑q=\sqrt{\frac{m}{k}}_{<E}𝑑p^{}𝑑q^{}=\frac{1}{\omega }_0^{2\pi }𝑑\theta _0^{r_0}r𝑑r=2\pi \frac{E}{\omega }.$$
(37)   
The nonlinearity can now be treated as a perturbation. To lowest order in the quartic terms, the radius to the boundary surface is given by

$$r_1^2r_0^22r_0^4\left(\frac{a}{k^2}\mathrm{cos}^4\theta +m^2b\mathrm{sin}^4\theta +\frac{mc}{k}\mathrm{sin}^2\theta \mathrm{cos}^2\theta \right).$$
(38)   
This integral is easily evaluated with the result

$$\frac{1}{2\pi }_{<E}𝑑p𝑑q\frac{E}{\omega }\frac{\omega }{2}\left(\frac{E}{\omega }\right)^2\left(3\frac{a}{k^2}+3bm^2+\frac{cm}{k}\right).$$
(39)   
Inserting the parameters from eqs. (29)-(34), the anharmonic term can be expressed

$$\frac{\omega }{8ϵN}\left(\frac{E}{\omega }\right)^2\left(1+3(ϵ/\omega )^42(ϵ/\omega )^2\right).$$
(40)   
Note that if there is no interaction, $`\omega =ϵ`$ and the anharmonicity vanishes. This is rather remarkable; the Hamiltonian in this case is the first term in eq. (28), which looks nonlinear. But the solution of the equations of motion are independent on excitation energy. It is not a harmonic oscillator spectrum, however, because the energy is bounded. These two properties correspond exactly to the quantum spectrum of the operator $`ϵJ_z`$.

We next quantize the above action to get

$$E_n=n\omega +n^2\frac{\omega ^2}{8ϵN}\left(1+3(ϵ/\omega )^42(ϵ/\omega )^2\right).$$
(41)   
Taking the second difference, this yields an anharmonicity of

$$\mathrm{\Delta }^{(2)}E=\frac{\omega ^2}{4ϵN}\left(1+3(ϵ/\omega )^42(ϵ/\omega )^2\right)=\frac{2\chi ϵ^3}{N\omega ^2}\left(1\frac{\chi }{2}\right).$$
(42)   
The exact value of anharmonicity $`\mathrm{\Delta }^{(2)}E`$ is compared with the value obtained from eq. (42) in Fig. 3. We can see that the time-dependent variational principle works very well.

## IV Nuclear Field Theory (NFT) Approach

The NFT is a formulation of many-body perturbation theory with vibrational modes summed to all orders in RPA. Its building blocks are RPA phonons and the single particle degrees of freedom which are described in the Hartree-Fock approximation. The coupling between them is treated diagramatically in the perturbation theory. For the Hamiltonian (1), the effective NFT Hamiltonian is given, to the lowest order, by

$$H_{NFT}=\frac{1}{2}ϵ\sigma _z+\omega O^{}O+H_{pv}.$$
(43)   
The first term in the $`H_{NFT}`$ describes single-particle spectrum. In writing down this term, we have used the fact that, for small value of $`V`$, the excitation energy in the HF is given by $`ϵ`$ and that the creation and the annihilation operators for the HF levels are the same as those for unperturbative levels. The second term describes the RPA phonons, with $`\omega `$ and $`O^{}`$ given by eqs. (12) and (5), respectively. The particle-vibration interaction $`H_{pv}`$ in eq. (43) is given as

$$H_{pv}=\mathrm{\Lambda }(O^{}+O)\sigma _x,$$
(44)   
where the coupling constant $`\mathrm{\Lambda }`$ is given by

$$\mathrm{\Lambda }=NV\sqrt{\frac{ϵ}{N\omega }}=ϵ\chi ^{}\sqrt{\frac{ϵ}{N\omega }},$$
(45)   
$`\chi ^{}`$ being $`NV/ϵ`$. This Hamiltonian is constructed by replacing $`\sigma _x`$ in the two-body interaction in the original Hamiltonian (1) as in Ref. . In general, there is also a residual interaction among particles and holes, but it does not contribute for the Hamiltonian (1) to the lowest order.

Each graph in the NFT contributes to a given order in $`1/N`$, but to all orders in $`\chi ^{}`$. Since the microscopic origin of the RPA phonon is a coherent sum of particle-hole excitations, bubble diagrams have to be excluded when one calculates physical quantities in the NFT . To the zero-order (in $`1/N`$), the phonon energy coincides with that in the RPA given by eq. (12). The anharmonicity begins with the leading $`1/N`$ diagrams, which are shown in Fig. 4. These diagrams are called “butterfly”graphs (see also refs. ). For each diagram shown in fig.4, there are 5 other diagrams which are obtained by changing the direction of the phonon lines. As already said, for the Hamiltonian (1) there is no diagram of order $`1/N`$ involving residual interaction among fermions.

The contribution from each diagram is most easily evaluated by using the Rayleigh-Schrödinger energy denominator, which are more suitable in the lowest order expansion . The four graphs in Fig. 4 have identical contributions, each given by

$$Graph(a)=\frac{N\mathrm{\Lambda }^4}{(2\omega \omega ϵ)^2(2\omega 2ϵ)}.$$
(46)   
In this equation, the minus sign appears because of the crossing of two fermion lines . By summing up the contributions from all diagrams, we obtain

$$\mathrm{\Delta }^{(2)}E=4N\mathrm{\Lambda }^4\frac{ϵ(\omega ^2+3ϵ^2)}{(\omega ^2ϵ^2)^3}=\frac{2\chi ^{}ϵ^3}{N\omega ^2}\left(1\frac{\chi ^{}}{2}\right).$$
(47)   
To compare with eq. (42), note that $`\chi ^{}=\chi `$ to the leading order of $`1/N`$. With this substitution, the two results are identical.

## V BOSON EXPANSION APPROACH

In the boson expansion method, each fermionic operators is replaced by corresponding operators which are written in terms of boson operators. There are several prescriptions to carry out the mapping from the fermionic to the bosonic spaces . Here we follow Refs.  which discussed the anharmonicities of the Lipkin model using the Holstein-Primakoff mapping. In this mapping, fermionic operators are mapped to bosonic operators so that the commutation relations among operators are preserved. The quasi-spin operators in the present two-level problem are then mapped as 

$`\sigma _+`$ $``$ $`2\sqrt{N}B^{}\sqrt{1{\displaystyle \frac{B^{}B}{N}}}`$ (48)   
$`\sigma _{}`$ $``$ $`2\sqrt{N}\sqrt{1{\displaystyle \frac{B^{}B}{N}}}B`$ (49)   
$`\sigma _z`$ $``$ $`N+2B^{}B,`$ (50)   
where the operators $`B`$ and $`B^{}`$ satisfy the boson commutation relation, i.e., $`[B,B^{}]=1`$. The Hamiltonian in the boson space which corresponds to the Lipkin Hamiltonian (1) is therefore obtained as

$`H_B`$ $``$ $`ϵ\left(1{\displaystyle \frac{NV}{ϵ}}\right)B^{}B{\displaystyle \frac{NV}{2}}\left(B^{}B^{}+BB\right)`$ (52)   
$`+VB^{}B+{\displaystyle \frac{V}{4}}\left(B^{}B^{}+BB\right)+VB^{}B^{}BB+{\displaystyle \frac{V}{2}}\left(B^{}B^{}B^{}B+B^{}BBB\right),`$   
to the first order in $`1/N`$.

A truncation of the expansion up to the leading order of the $`1/N`$ corresponds to the RPA which we discussed in Sec. II. To this order, the boson Hamiltonian (52) is given by

$$H_B^{(2)}=ϵ\left(1\frac{NV}{ϵ}\right)B^{}B\frac{NV}{2}\left(B^{}B^{}+BB\right).$$
(53)   
As is well known, this Hamiltonian can be diagonalised by introducing a transformation

$`B^{}`$ $`=`$ $`X_0O^{}+Y_0O`$ (54)   
$`B`$ $`=`$ $`X_0O+Y_0O^{},`$ (55)   
and imposing

$$[H_B^{(2)},O^{}]=\omega O^{},$$
(56)   
with a condition $`X_0^2Y_0^2=1`$. The frequency $`\omega `$ then reads

$$\omega =ϵ\sqrt{12\chi ^{}},$$
(57)   
$`\chi ^{}`$ being $`NV/ϵ`$, which was introduced in the previous section, together with

$`X_0^2`$ $`=`$ $`{\displaystyle \frac{\omega +ϵ(1\chi ^{})}{2\omega }},`$ (58)   
$`Y_0^2`$ $`=`$ $`{\displaystyle \frac{\omega +ϵ(1\chi ^{})}{2\omega }},`$ (59)   
$`X_0Y_0`$ $`=`$ $`{\displaystyle \frac{ϵ\chi ^{}}{2\omega }}.`$ (60)   
The frequency $`\omega `$ coincides with that in the RPA given by eq. (12) to the leading order of $`1/N`$.

As in the nuclear field theory approach, the anharmonicity begins with the next order of the $`1/N`$ expansion. In terms of RPA phonon creation and annihilation operators defined in eqs. (54) and (55), the boson Hamiltonian (52) can be rewritten as

$`H_B`$ $`=`$ $`\omega O^{}O+H_{11}O^{}O+H_{20}\left(O^2+O^2\right)`$ (62)   
$`+H_{40}\left(O^4+O^4\right)+H_{31}\left(O^3O+O^{}O^3\right)+H_{22}O^2O^2,`$   
where the first term is the leading order of the $`1/N`$ expansion given by eq. (53) and the rest are the higer order corrections. The coefficient $`H_{22}`$, for example, is given by

$$H_{22}=V(X_0^4+3X_0^3Y_0+4X_0^2Y_0^2+3X_0Y_0^3+Y_0^4)=\frac{Vϵ^2}{\omega ^2}\left(1\frac{\chi ^{}}{2}\right).$$
(63)   
In order to estimate the degree of anharmonicity, we use the perturbation theory. The first order perturbation gives the energy of the one and the two phonon states of

$`E_1`$ $`=`$ $`\omega +H_{11},`$ (64)   
$`E_2`$ $`=`$ $`2(\omega +H_{11})+2H_{22},`$ (65)   
respectively. Taking the second difference, the anharmonicity reads

$$\mathrm{\Delta }^{(2)}E=2H_{22}=\frac{2\chi ^{}ϵ^3}{N\omega ^2}\left(1\frac{\chi ^{}}{2}\right),$$
(66)   
which is identical to the result of the variational approach given by eq. (42) as well as that of the nuclear field theory, eq. (47).

## VI Conclusion

We have shown that the nuclear field theory, the time-dependent variational principle, and the boson expansion method give identical leading-order anharmonicities for the Lipkin model, and that the formulas agree well with the exact numerical solution. The anharmonicity is inversely proportional to the number of particles in the system, when the other parameters are fixed to keep the harmonic frequency the same. This clarifies the origin of the conflicting results for the $`A`$-dependence of the anharmonicity obtained in ref.  and . In ref.  the time-dependent method was applied to a Skyrme-like Hamiltonian involving all $`A`$ nucleons, and the result was $`\mathrm{\Delta }^{(2)}Ef(\omega )/A`$. In ref. , the Hamiltonian was restricted to a space of a single major shell for particle orbitals and similarly for the hole orbitals. Since the number of particles in the valence shell increases as $`A^{2/3}`$, the result was $`\mathrm{\Delta }^{(2)}Ef(\omega )/A^{2/3}`$. Finally, it should perhaps be emphasized that both methods predict that the anharmonicity is very small for giant resonances: both $`A^{2/3}`$ and $`A`$ are large numbers. This need not be the case for low-lying collective vibrations. The NFT can produce large effects when there are small energy denominators. Low-lying excitations are difficult to describe with a simple ansatz like eq. (15), so the time-dependent variational principle is not easily applied. Clearly this is an area that should be explored further.

As for the discrepancy between the time-dependent variational approach and the boson expansion method concerning anharmonicities of plasmon resonances of metal clusters, our study showed that the origin of the discrepancy should not be ascribed to the method used to solve the problem. The origin of the discrepancy, therefore, is not traceable at moment and further studies may be necessary in order to reconcile the discrepancy.

## Acknowledgments

We thank A. Bulgac for showing us the Hamiltonian construction and C. Volpe for discussions. K.H. acknowledges support from the Japan Society for the Promotion of Science for Young Scientists. G.F.B. acknowledges support from the U.S. Dept. of Energy under Grant DE-FG-06ER46561. Both G.F.B. and P.F.B. thank the IPN of Orsay for warm hospitality where this work was completed.

Figure Captions

Fig. 1: Energy of the one phonon state as a function of the number of particle. The solid line is the exact solution of the Lipkin model, whilist the dashed line is obtained in the RPA. The strength of the interaction $`V/ϵ`$ is set to be 0.01.

Fig. 2: Same as fig. 1, but for a fixed $`\chi `$ parameter, which is set to be 0.25.

Fig. 3: Anharmonicity of the double phonon state as a function of the number of particle for $`\chi =0.25`$. The dashed line is an approximate solution given by eq. (42), while the solid line is obtained by numerically diagonalising the Hamiltonian.

Fig. 4: The diagrams of order $`1/N`$ contributing to the anharmonicity of the double phonon state for the Hamiltonian (1). The wavy lines represent phonon propagations, while the particles and the holes are depicted by the allowed lines.