arXiv:0704.0152v2  [nucl-th]  18 Dec 2007
Kinetic equation for finite systems of fermions
with pairing
V. I. Abrosimov a, D. M. Brink b, A. Dellafiore c,∗, F. Matera c,d
aInstitute for Nuclear Research, 03028 Kiev, Ukraine
bOxford University, Oxford, U.K.
cIstituto Nazionale di Fisica Nucleare, Sezione di Firenze
dDipartimento di Fisica, Università degli Studi di Firenze, via Sansone 1, I 50019
Sesto F.no (Firenze), Italy
Abstract
The solutions of the Wigner-transformed time-dependent Hartree–Fock–Bogoliubov
equations are studied in the constant-∆ approximation. This approximation is
known to violate particle-number conservation. As a consequence, the density fluc-
tuation and the longitudinal response function given by this approximation contain
spurious contributions. A simple prescription for restoring both local and global
particle-number conservation is proposed. Explicit expressions for the eigenfrequen-
cies of the correlated systems and for the density response function are derived and
it is shown that the semiclassical analogous of the quantum single–particle spectrum
has an excitation gap of 2∆, in agreement with the quantum result. The collective
response is studied for a simplified form of the residual interaction.
Key words: Pairing, Vlasov equation
PACS: 21.10.Pc, 03.65.Sq
1 Introduction
The problem of extending the Vlasov equation to systems in which pairing
correlations play an important role has been tackled some time ago by Di
Toro and Kolomietz [1] in a nuclear physics context and, more recently, by
Urban and Schuck [2] for trapped fermion droplets. These last authors derived
∗ Corresponding author
Email address: della@fi.infn.it (A. Dellafiore).
Preprint submitted to Elsevier 16 November 2021
http://arxiv.org/abs/0704.0152v2
the TDHFB equations for the Wigner transform of the normal density matrix
ρ and of the pair correlation function κ (plus their time-reversal conjugates)
and used them to study the dynamics of a spin-saturated trapped Fermi gas.
In the time-dependent theory one obtains a system of four coupled differential
equations for ρ, κ, and their conjugates [2] and, if one wants an analytical so-
lution, some approximation must be introduced. Here we try to find a solution
of the equations of motion derived by Urban and Schuck in the approximation
in which the pairing field ∆(r,p, t) is treated as a constant. It is well known
that such an approximation violates both particle-number-conservation and
gauge invariance (see e.g. sect. 8-5 of [3] and [4]), nonetheless we study it
because of its simplicity, with the aim of correcting the final results for its
shortcomings. Moreover, the constant-∆ approximation is not satisfactory for
describing long wavelength pairing modes in a large system. Such modes have
frequencies which are much less than the pairing frequency ∆/~ and for their
study it is essential to use a self consistent theory where the gap ∆ is related
to the pair density κ through the pairing interaction. The phases of ∆ and
κ are particularly important because they describe the collective superfluid
currents. On the other hand nuclei are small systems. Shell gaps are large
compared with ∆, or equivalently giant resonance frequencies are large com-
pared with the pairing frequency. The constant-∆ approximation is much more
reasonable in such systems.
In Sect. 2, the basic equations are recalled and reformulated in terms of the
even and odd components of the normal density ρ. In Sect. 3, the static limit
is studied by following the approach of [5] and the constant-∆ approximation
is introduced. In Sect. 4, the simplified dynamic equations resulting from the
constant-∆ approxmation are derived and their solutions are determined in lin-
ear approximation. In Sect. 5, these solutions are studied in a one-dimensional
model and the problem of particle-number conservation is examined in detail.
By studying the energy-weighted sum rule (in the Appendix), we find that the
constant-∆ approximation introduces some spurious strength into the density
response of the system. A simple prescription, based on the continuity equa-
tion, is proposed in order to eliminate the spurious strength. The resulting
strength function gives the same energy-weighted sum rule as for the uncor-
related systems. In Sect. 6, the general solution found in Sect. 4 is re-written
for spherical systems, where the angular integrations can be performed ex-
plicitly, leading to expressions containing only radial integrations. In Sect.
7, the collective response function of spherical nuclei is derived for a simple
multipole-multipole residual interaction. In Sect. 8, the quadrupole and oc-
tupole channels, that are the ones most affected by the pairing correlations
are shown explicitly. Finally, in Sect. 9 conclusions are drawn.
2 Basic equations
We assume that our system is saturated both in spin and isospin space and
do not distinguish between neutrons and protons, so we can use directly the
equations of motion of Urban and Schuck.
We start from the equations of motion derived in Ref. [2] for the Wigner-
transformed density matrices ρ = ρ(r,p, t) and κ = κ(r,p, t), with the warning
that the sign of κ that we are using agrees with that of Ref. [1], hence it is
opposite to that of [2]. Moreover we find convenient to use the odd and even
combinations of the normal density introduced in [2]:
ρev =
[ρ(r,p, t) + ρ(r,−p, t)] , (1)
ρod =
[ρ(r,p, t)− ρ(r,−p, t)] . (2)
Thus, the equations of motion given by Eqs.(15a...d) of Ref. [2] read
i~∂tρev = i~{h, ρod} − 2iIm[∆∗(r,p, t)κ] (3)
i~∂tρod = i~{h, ρev}+ i~Re{∆∗(r,p, t), κ} (4)
i~∂tκ=2(h− µ)κ−∆(r,p, t)(2ρev − 1) + i~{∆(r,p, t), ρod} . (5)
Here h is the Wigner-transformed Hartree–Fock hamiltonian h(r,p, t), while
∆(r,p, t) is the Wigner-transformed pairing field. Since the time-dependent
part of κ is complex, κ = κr+iκi, the last equation gives two separate equations
for the real and imaginary parts of κ.
Moreover, from the supplementary normalization condition ([6], p. 252)
R2 = R (6)
satisfied by the generalized density matrix R, the two following independent
equations are obtained:
ρodκ+ i
{ρev, κ}=0 , (7)
ρev(ρev − 1) + ρ2od + κκ∗=0 . (8)
We shall use the equations of motion (3–5), together with these equations, as
our starting point, but first we notice that, in the limit of no pairing, both ∆
and κ vanish, the third equation of motion reduces to a trivial identity, while
the first two give the Vlasov equation for normal systems, expressed in terms
of the even and odd components of ρ:
∂tρev = {h, ρod} , (9)
∂tρod = {h, ρev} . (10)
A solution of the linearized Vlasov equation for normal systems (i. e. without
pairing) has been obtained in Ref. [7] and our aim here is to study the changes
introduced by the pairing interaction in the solution of [7].
Moreover, before studying the time-dependent problem, it is useful to look at
the static limit.
3 Static limit
In this section we follow the approach of Ref. [5]. At equilibrium we have
ρev = ρ0(r,p) , (11)
ρod =0 , (12)
h=h0(r,p) , (13)
κ=κ0(r,p) (14)
∆0 =∆0(r,p) (15)
and equations (3–5) give
0=−2iIm(∆∗0κ0) , (16)
0= i~{h0, ρ0}+ i~Re{∆∗0, κ0} , (17)
0=2(h0 − µ)κ0 −∆0(2ρ0 − 1) , (18)
while Eqs. (7, 8) give
{ρ0, κ0}=0 , (19)
ρ0(ρ0 − 1) + |κ0|2=0 . (20)
Equation (16) is satisfied if we assume that ∆0 and κ0 are real quantities,
while Eqs. (18) and (20), taken as a system, have the solution:[5]
ρ0(r,p) =
h0(r,p)− µ
E(r,p)
κ0(r,p) = −
∆0(r,p)
2E(r,p)
, (22)
with the quasiparticle energy
E(r,p) =
∆20(r,p) + (h0(r,p)− µ)2 . (23)
It can be easily checked that Eqs. (21, 22) satisfy also Eqs. (17) and (19), that
{h0, ρ0}+ {∆0, κ0} = 0 (24)
{ρ0, κ0} = 0 . (25)
The (semi)classical equilibrium phase-space distribution is closely related to
ρ0(r,p):
f0(r,p) =
(2π~)3
ρ0(r,p) (26)
and the statistical factor 4 takes into account the fact that there are two kinds
of fermions.
The parametrer µ is determined by the condition
drdpf0(r,p) , (27)
where A is the number of particles. This integral should keep the same value
also out of equilibrium (global particle-number conservation).
3.1 Constant-∆ approximation
In a fully self-consistent approach, the pairing field ∆(r,p, t) is related to
κ(r,p, t), however here we introduce an approximation and replace the pairing
field of the HFB theory with the phenomenological pairing gap of nuclei, hence
in all our equations we put
∆(r,p, t) ≈ ∆0(r,p) ≈ ∆ = const. , (28)
with ∆ ≈ 1MeV.
In the constant-∆ approximation the equilibrium distributions become
ρ0(ǫ) =
1− ǫ− µ
, (29)
κ0(ǫ) =−
2E(ǫ)
and the quasiparticle energy
E(ǫ) =
∆2 + (ǫ− µ)2 , (31)
ǫ = h0(r,p) =
+ V0(r)
the particle energy in the equilibrium mean field.
In the following we shall use the relation:
κ0(ǫ) =
E2(ǫ)
dρ0(ǫ)
. (32)
4 Dynamic equations
Always in the approximation where ∆ is constant and real, the time-dependent
equations (3–5) become
i~∂tρev = i~{h, ρod} − 2i∆Im(κ) (33)
i~∂tρod = i~{h, ρev} (34)
i~∂tκ=2(h− µ)κ−∆(2ρev − 1) . (35)
This is the simplified set of equations that we want to study here. The sum of
the first two equations gives an equation that is similar to the Vlasov equation
of normal systems, only with the extra term −2i∆Im(κ). This extra term
couples the equation of motion of ρ with that of κ, thus, instead of a single
differential equation (Vlasov equation), now we have a system of two coupled
differential equations (for ρ and κi).
Our aim here is that of determining the effects of pairing on the linear response
of nuclei, thus we assume that our system is initially at equilibrium, with
densities given by Eqs. (29,30), and that at time t = 0 a weak external field
of the kind
δV ext(r, t) = βδ(t)Q(r) (36)
is applied to it. This simple time-dependence is sufficient to determine the
linear response of the system. In a self-consistent approach, we should take
into account also the changes of the mean field surrounding each particle
induced by the external force and consider a perturbing hamiltonian of the
δh = δV ext + δV int , (37)
however we start with the zero-order approximation
δh = δV ext (38)
and will consider collective effects in a second stage.
Since we want to solve Eqs. (33–35) in linear approximation, we consider small
fluctuations of the time-dependent quantities about their equilibrium values
and neglect terms that are of second order in the fluctuations. Hence, in Eqs.
(33–35) we put:
h=h0 + δh , (39)
ρev = ρ0 + δρev , (40)
ρod = δρod , (41)
κ=κ0 + δκ = κ0 + δκr + iδκi . (42)
Then, the linearized form of Eqs. (33–35) is
i~∂tδρev = i~{h0, δρod} − 2i∆δκi (43)
i~∂tδρod= i~{h0, δρev}+ i~{δh, ρ0} (44)
−~∂tδκi=2(ǫ− µ)δκr + 2κ0δh− 2∆δρev , (45)
~∂tδκr =2(ǫ− µ)δκi . (46)
Taking the sum of the first two equations gives
i~∂tδρ(r,p, t) = (47)
i~{h0, δρ(r,p, t)}+ i~{δh(r,p, t), ρ0} − 2i∆δκi(r,p, t) ,
which can be regarded as an extension of the linearized Vlasov equation stud-
ied in [7].
In order to make the comparison with [7] easier, from now on we change the
normalization of the phase-space densities and define
f(r,p, t)=
(2π~)3
ρ(r,p, t) (48)
χ(r,p, t)=
(2π~)3
κ(r,p, t) , (49)
moreover, we put
F (ǫ) =
(2π~)3
ρ0(ǫ) . (50)
In terms of the new functions Eqs. (47) and (45) read
i~∂tδf(r,p, t) = i~{h0, δf(r,p, t)}+ i~{δh(r,p, t), f0}
− 2i∆δχi(r,p, t) , (51)
−~∂tδχi(r,p, t) = 2(ǫ− µ)δχr(r,p, t) + 2χ0δh(r,p, t)
− 2∆δfev(r,p, t) . (52)
The function fev is given by the obvious extension of Eq. (1). In order to get
a closed system of equations, we still need an extra equation for δχr(r,p, t).
This can be obtained from the linearized form of the supplementary condition
(8) that reads
δρev(2ρ0 − 1) = −2κ0δκr , (53)
δκr(r,p, t) =
1− 2ρ0(ǫ)
δρev(r,p, t) = −
δρev(r,p, t) . (54)
The last expression has been obtained with the help of Eq. (18). In terms of
the new functions f and χ, the last equation reads
δχr(r,p, t) = −
δfev(r,p, t) . (55)
Equations (51, 52) and (55) are the set of coupled equations for the phase-
space densities that we have to solve.
Replacing Eq. (55) into Eq. (52), and using Eq. (30), gives the following system
of coupled differential equations:
∂tδf(r,p, t)= {h0, δf}+ {δh, f0} − 2
δχi(r,p, t) , (56)
∂tδχi(r,p, t)=
E2(ǫ)
[δf(r,p, t) + δf(r,−p, t)]− 2
δh(r,p, t) . (57)
Taking the Fourier transform in time, gives
−iωδf(r,p, ω)= {h0, δf}+ {δh, f0} − 2
δχi(r,p, ω) , (58)
−iωδχi(r,p, ω))=
E2(ǫ)
[δf(r,p, ω) + δf(r,−p, ω)]
− 2χ0
δh(r,p, ω) , (59)
or (for ω 6= 0)
−iωδf(r,p, ω) + {δf, h0}=−iωd2
[δf(r,p, ω) + δf(r,−p, ω)]
+F ′(ǫ)[{δh, h0}+ iωd2δh] , (60)
(Ω(ǫ)
Ω(ǫ) = 2
. (62)
This frequency plays a crucial role in our approach, its minimum value is
2∆/~.
In Eq. (60) we have used the relation {f0, δh} = F ′(ǫ){ho, δh} as well as Eq.
(32).
By comparing Eq. (60) with the analogous equation for normal systems
−iωδf(r,p, ω) + {δf, h0} = F ′(h0){δh, h0} , (63)
we can see that the only effect of pairing in the constant-∆ approximation is
that of adding the terms proportional to d2.
The normal Vlasov equation (63) can be solved in a very compact way by using
the method of action-angle variables [8], [7]. In that approach one expands
δhn(I)e
in·Φ , (64)
where I and Φ are the action and angle variables, respectively. Moreover
δf(r,p, ω) =
δfn(I, ω)e
in·Φ (65)
{δf, h0} =
i(n · ~ω)δfn(I, ω)ein·Φ , (66)
where the vector ~ω has components
. (67)
Then Eq.(63) gives
δfn(I, ω) =
(n · ~ω)
(n · ~ω)− (ω + iε)
F ′(ǫ)δhn(I) . (68)
The (zero-order) eigenfrequencies of the (normal) physical system are
ωn = n · ~ω . (69)
Here want to use the same method to solve the more complicated equation
(60). Since that equation contains also the function δf(r,−p, ω), we need also
the analogous equation for this other quantity:
−iωδf(r,−p, ω)− {δf, h0}r,p=−iωd2
[δf(r,p, ω) + δf(r,−p, ω)]
+F ′(ǫ)[−{δh, h0}r,p + iωd2δh] . (70)
By expanding δf(r,p, ω) and δf(r,−p, ω) as
δf(r,±p, ω) =
δf±n (I, ω)e
in·Φ , (71)
Eqs. (60) and (70) give
[−ω(1− d
) + ωn]δf
n + ω
δf−n = F
′(ǫ)[ωn + ωd
2]δhn , (72)
]δf+n + [−ω(1−
)− ωn]δf−n = F ′(ǫ)[−ωn + ωd2]δhn , (73)
which is a system of two coupled algebraic equations for the coefficients δf+n
and δf−n . Its solution is
δf+n =
ω̄2n + ωωn
ω̄2n − ω2
F ′(ǫ)δhn , (74)
δf−n =
ω̄2n − ωωn
ω̄2n − ω2
F ′(ǫ)δhn , (75)
where
ω̄2n = ω
n + Ω
2(ǫ) (76)
are the (squared) eigenfrequencies of the correlated system. These eigenfre-
quencies are in agreement with the enegy spectrum of a superfluid infinite ho-
mogeneous Fermi gas (see e. g. Sect. 39 of [10]) and they lead to a low-energy
gap of 2∆ in the excitation spectrum of the correlated systems. However,
as anticipated, we expect problems with particle-number conservation. These
problems are better discussed in one dimension, where formulae are simpler.
5 One-dimensional systems and particle-number conservation
In one dimension, Eq.(60) reads
−iωδf(x, p, ω) + ẋ∂xδf −
dV0(x)
∂pδf = (77)
−iωd2 1
[δf(x, p, ω) + δf(x,−p, ω)] + F ′(ǫ)( p
∂xδh + iωd
2δh) .
In zero-order approximation δh(x, p, ω) = βQ(x), moreover, in one dimension
F ′(ǫ) =
. (78)
The vectors ~ω and Φ have only one component:
ω0(ǫ) =
T (ǫ)
Φ(x) = ω0τ(x) , (80)
τ(x) =
v(ǫ, x′)
, (81)
v(ǫ, x) =
[ǫ− V0(x)] . (82)
The time T (ǫ) is the period of the bound motion of particles with enegy ǫ
in the equilibrium potential well V0(x): T = 2τ(x2). The points x1,2 are the
classical turning points for the same particles. Instead of the action variable
I(ǫ) = 1
dxp(ǫ, x), it is more convenient to use the particle energy ǫ as
constant of motion. As pointed out in [7], the range of values of τ can be
extended to the whole interval (0, T ), by defining
τ(x) =
v(ǫ, x′)
when τ > T
. With this extension, the angle variable Φ(x) takes values between
0 and 2π, as it should.
In one dimension, Eqs. (74, 75) give
δf±n (ǫ, ω) =
ω̄2n ± ω ωn
ω̄2n − ω2
F ′(ǫ)δhn , (84)
δhn = βQn
v(ǫ, x)
e−iωnτ(x)
v(ǫ, x)
cos[ωnτ(x)] , (85)
The frequencies ωn are the eigenfrequencies of the uncorrelated system:
ωn = nω0 , (86)
while ω̄n are the new eigenfrequencies modified by the pairing correlations:
ω̄n = ±
ω2n + Ω
2(ǫ) . (87)
Note that, since δh−n = δhn, then δf
n = δf
By using the solutions (84), we can also obtain an expansion for the even and
odd parts of δf :
δfev(x, ǫ, ω) =
An(ω) cosnω0τ(x) , (88)
δfod(x, ǫ, ω) =
Bn(ω) sinnω0τ(x) , (89)
An(ω) =
ω̄2n − ω2
F ′(ǫ)δh′n , (90)
Bn(ω) = iω
ω̄2n − ω2
F ′(ǫ)δh′n (91)
δh′n =2δhn , n 6= 0 , (92)
= δhn , n = 0 . (93)
Note that, while Bn=0(ω) = 0, we have An=0(ω) 6= 0, and this fact leads to
an unphysical fluctuation of the number of particles, induced by the applied
external field. These fluctuations are given by
δA(ω) =
dxδ̺(x, ω) , (94)
where δ̺(x, ω) is the density fluctuation at point x:
δ̺(x, ω) =
dp δf(x, p, ω) = 2
v(ǫ, x)
δfev(x, ǫ, ω) . (95)
Equation (88) gives
δA(ω) = 2
An(ω)
dτ cos nω0τ . (96)
Since the integrals
∫ T/2
0 dτ cosnω0τ vanish when n 6= 0, the term with n = 0
is the only one contributing to this sum, thus givig an unphysical fluctuation
of the number of particles. This problem could be solved simply by excluding
the mode n = 0 from the sum in Eq. (88), however this would not be sufficient
to solve all problems with particle-number conservation, since we can easily
check that the solutions (88, 89) do not satisfy the continuity equation
iω̺(x, ω) = ∂xj(x, ω) . (97)
The density fluctuation involves only the even part of δf , while the current
density j(x, ω) involves only the odd part:
j(x, ω) =
δf(x, p, ω) = 2
dǫδfod(x, ǫ, ω) . (98)
The fact that the continuity equation is violated is a very serious shortcom-
ing of the constant-∆ approximation. However, since we have seen that this
approximation leads to very simple equations and to rather satisfactory ex-
pressions for the eigenfrequencies of the correlated systems, we still use it,
but with the following prescription: when calculating the longitudinal response
function, the density fluctuations should be evaluated by using Eq. (97), instead
of Eq, (95). Then, the density fluctuations (95) should be replaced by
δ ¯̺(x, ω) =
dǫ∂xδfod(x, ǫ, ω) . (99)
In practice we are proposing to evaluate the longitudinal response function
in terms of the transverse response function. It is well known that also the
more familiar BCS approximation gives a more accurate description of the
transverse response (see e.g. sect. 8-5 of [3]). In the Appendix we show that the
longitudinal response function resulting from the present prescription satisfies
the same energy-weighted sum rule as the uncorrelated response function. This
would not necessarily happen if, instead of changing only the even part of δf ,
we had modified also its odd part.
It is interesting to see how the solutions (84) are changed by our prescription.
By using Eq. (89) forfod, Eq. (99) gives
δ ¯̺(x, ω) = 2
v(ǫ, x)
δf̄ev(x, ǫ, ω) , (100)
δf̄ev(x, ǫ, ω) =
Ān(ω) cosnω0τ(x) , (101)
Ān(ω) =
Bn(ω) , (102)
(note that Ān=0(ω) = 0). Then
δf̄(x,±p, ω) = δf̄ev(x, ǫ, ω)± δfod(x, ǫ, ω) (103)
and Eq. (84) is replaced by
δf̄±n (ǫ, ω) =
ω2n ± ω ωn
ω̄2n − ω2
F ′(ǫ)δhn . (104)
By comparing this expression to Eq. (84), we can see that the fluctuations
of the phase-space density given by the constant-∆ approximation contain an
extra contribution that we identify as spurious:
δf±n (ǫ, ω) = δf̄
n (ǫ, ω) + δf
n (ǫ, ω) , (105)
δf spurn (ǫ, ω) =
Ω2(ǫ)
ω̄2n − ω2
F ′(ǫ)δhn . (106)
The spurious character of δf spurn is suggested also by sum-rule arguments
(see Appendix). The term f spurn (ǫ, ω) contributes to all modes of the density
strength function: the contribution to the mode n = 0 gives a fluctuation of
the particle-number integral (global paticle-number violation), while the other
modes give a spurious contribution to the density strength function, increas-
ing the sum rule and violating the continuity equation (local particle-number
violation). Note that the spurious contribution (106) affects only the even part
of the pase-space density, not the odd part.
6 Spherical Systems
The method of action-angle variables gives a very compact solution of the lin-
earized Vlasov equation both in the uncorrelated and correlated cases, however
it may be useful to make a connection between the results given by this method
and the more explicit treatment of spherical nuclei given in [7]. For uncorre-
lated system this has been done in [11]. Here we follow that approach in order
to derive useful expressions for correlated spherical systems. The components
of the vecor n are (n1, n2, n3), the first point to notice is that, because of the
degeneracy associated with any central-force field, the vector ~ω has only two
non vanishing components:
~ω = (0, ωϕ(ǫ, λ), ω0(ǫ, λ)) . (107)
With λ we denote the magnitude of the particle angular momentum. According
to Eq. (69), the eigenfrequencies of the uncorrelated system are [7]
ωn = ωn3,n2(ǫ, λ) = n3ω0 + n2ωϕ , (108)
while Eq. (76) gives the correlated eigenfrequencies
ω̄n = ω̄n3,n2 = ±
ω2n3,n2 + Ω
2(ǫ) . (109)
In three dimensions, the Fourier coefficients analogous to (85) are
Qn(I) =
(2π)3
dΦe−in·ΦQ(r) . (110)
The external field Q(r) can be expanded in partial waves as
Q(r) =
QL(r)YLM(r̂) , (111)
giving
Qn(I) =
Q(LM)n , (112)
with [11]
Q(LM)n =
dLMN(β
′)δM,n1δN,n2Q
. (113)
By using this last equation (and changing n3 → n), the expansion (71) be-
comes
δf(r,±p, ω) =
(114)
δfL±nN (ǫ, λ, ω)e
iφnN (r)
DLMN(α, β ′, γ)
δfL±nN (ǫ, λ, ω) =
ω̄2nN ± ωωnN
ω̄2nN − ω2
βF ′(ǫ)QLnN (115)
and QLnN the semiclassical limit of the radial matrix elements:
QLnN =
vr(r)
e−iφnN (r)QL(r) ,
vr(r)
cos[φnN(r)]QL(r) . (116)
Here T is the period of radial motion, vr(r) the radial velocity
vr(r) =
ǫ− V0(r)−
(117)
and the phases φnN(r) are given by
φnN(r) = ωnNτ(r)−Nγ(r) , (118)
where
τ(r) =
vr(r′)
(119)
γ(r) =
vr(r′)
. (120)
The frequencies ω0 and ωϕ are given by
τ(r2)
, (121)
γ(r2)
τ(r2)
. (122)
The Wigner rotation matrix elements in Eq. (114) are given by [9]
DLMN(α, β ′, γ) = e−iMαdLMN(β ′)e−iNγ , (123)
where (α, β ′, γ) are the Euler angles introduced in [7].
On the basis of the discussion in Sect. (5), we expect that the solution (114)
will contain some spurious strength introduced by the constant-∆ approxima-
tion. In order to eliminate the spurious contributions, we should replace the
coefficients (115) with
δf̄L±nN (ǫ, λ, ω) =
ω2nN ± ωωnN
ω̄2nN − ω2
βF ′(ǫ)QLnN . (124)
These modified coefficients allow us to obtain the modified zero-order propa-
gator
D̄0L(r, r
′, ω) =
dǫF ′(ǫ)
d̄LnN(r, r
ω − ω̄nN + iε
, (125)
d̄LnN(r, r
′) = (126)
(4π)2
2L+ 1
|YLN(
(−2ωnN
)(ωnN
)cosφnN(r)
r2vr(r)
cosφnN(r
r′2vr(r′)
and the corresponding response and strength functions:
R̄0L(ω)=
drdr′r2QL(r)D̄
L(r, r
′, ω)r′2QL(r
′) , (127)
S̄0L(ω)=−
ImR̄0L(ω) . (128)
For multipole response: QL(r) = r
For normal systems, the zero-order propagator D0L(r, r
′, ω) is given by Eqs.
(125) and(126) where ω̄nN is replaced by ωnN and F
′(ǫ) is proportional to a
δ-function[7].
7 Collective response
Up to now, we have been concerned only with the zero-order approximation,
which corresponds to the single-particle approximation of the quantum ap-
proach. In this approximation, the perturbing part of the hamiltonian is given
only by the external field, while a more consistent approach would require
taking into account also the mean-field fluctuation induced by the external
force, so that
δh = δV ext(r, ω) + δV int(r, ω) . (129)
In the Hartree approximation,
δV int(r, ω) =
dr′v(r− r′)δ̺(r′, ω) . (130)
where v(r− r′)is the (long-range) interaction between constituents.
For consistency, we take
δV int(r, ω) =
dr′v(r− r′)δ ¯̺(r′, ω) , (131)
then the collective propagator for correlated systems satisfies the same kind
of integral equation as for normal systems [7]:
D̄L(r, r
′, ω) = (132)
D̄0L(r, r
′, ω) +
dyy2D̄0L(r, x, ω)vL(x, y)D̄L(y, r
′, ω) .
Here vL(x, y) is the partial-wave component of the interaction between par-
ticles. We assume that this interaction can be approximated by a separable
form of the kind
vL(x, y) = κLx
LyL , (133)
where κL is a parameter that determines the strength of the interaction. Then,
the integral equation (132) gives an algebraic equation for the collective cor-
related response function
R̄L(ω) =
drr2rL
dr′r′2r′LD̄L(r, r
′, ω)
leading to the expression
R̄L(ω) =
R̄0L(ω)
1− κLR̄0L(ω)
. (134)
8 Results
Here we compare the multipole strength functions given by our simplified
model of pairing correlations with that of the corresponding uncorrelated sys-
tem. This comparison is made for the quadrupole and octupole strength func-
tions, since these channels are the ones that are most affected.
The static nuclear mean field is approximated with a spherical cavity of radius
R = 1.2A
3 fm and the A nucleons are treated on the same footing, i. e., we do
not distinguish between neutrons and protons. Moreover, we chose A = 208
for ease of comparison with previous calculations of uncorrelated response
functions [12,13]. Shell effects are not included in our semiclassical picture and
the results shown below should be considered as an indication of the qualitative
effects to be expected in heavy nuclei. For the uncorrelated calculations, the
Fermi energy is determined by the parametrization chosen for the radius as
ǫF ≈ 33.33 MeV, while for the correlated case, the parameter µ is determined
by the condition (27); with the value of ∆ = 1 MeV used here, the value of µ
is practically coincident with that of ǫF , so we have used µ = ǫF = 33.33 MeV
in the calculations below. Moreover, the small parameter ε appearing in Eq.
(125) has been given the value ε = 0.1 MeV. This value is chosen to simplify
the evaluation of the response function by smoothing out discontinuities in its
dependence on ω.
In the evaluation of the collective response, the value of parameters κL is
the same as in [12,13], that is: κ2 = −1 × 10−3 MeV/fm4 and κ3 = −2 ×
10−5 MeV/fm6.
8.1 Quadrupole response
Figure 1 show the longitudinal quadrupole strength function evaluated in the
zero-order approximation (corresponding to the quantum single-particle ap-
proximation). The dashed curve shows the uncorrelated response evaluated
according to the theory of [7], while the full curve shows the result of the
present correlated calculation. As we can see the effect of pairing correlations
on this zero-order strength function is rather small, however, since pairing af-
fects also the real part of the zero-order response function,in Fig. 2 we plot
0 5 10 15 20 25 30
ÑΩ @MeVD
Fig. 1. Quadrupole stength function in zero-order approximation. The dashed curve
gives the response of a normal system of A = 208 nucleons contained in a spherical
cavity, while the solid curve includes the effects of pairing correlations in constant-∆
approximation.
0 5 10 15 20 25 30
ÑΩ @MeVD
Fig. 2. Collective quadrupole strength function showing the giant quadrupole res-
onance. The solid curve involves also pair correlations, the dashed curve has no
pairing.
also the collective strength function given by Eq. (134). Again, the effect of
pairing correlations is very small, in agreement with the results of [1].
0 1 2 3 4 5
ÑΩ @MeVD
0.0025
0.005
0.0075
0.0125
0.015
0.0175
Fig. 3. Same as Fig.1, at low excitation energy.
The main difference between the uncorrelated and correlated responses occurs
at small excitation energy, Fig. 3 shows a detail of Fig. 1 at low excitation
energy. The correlated strength function displays a gap of about 2 ∆ , the
very small strength extending below 2 MeV is entirely due to the finite value
of the small parameter ε used in the numerical evaluation of the propagator
(125).
8.2 Octupole response
Figures 4 and 5 show the zero-order and collective octupole strength functions,
both correlated and uncorrelated. As we can see, in this case too the effect is
rather small.
0 5 10 15 20 25 30
ÑΩ @MeVD
Fig. 4. Zero-order octupole strength function. The solid curve involves pair correla-
tions, the dashed curve has no pairing.
0 5 10 15 20 25 30
ÑΩ @MeVD
Fig. 5. Collective octupole strength function. The solid curve involves pair correla-
tions, the dashed curve has no pairing.
9 Conclusions
The solutions of the semiclassical time-dependent Hartree–Fock–Bogoliubov
equations have been studied in a simplified model in which the pairing field
∆(r,p, t) is treated as a constant phenomenological parameter. Such an ap-
proximation is known to violate some important constraints, like global (particle-
number integral) and local (continuity equation) particle-number conserva-
tion. In a linearized approach , we have shown that the global particle-number
violation is related only to one particular mode of the density fluctuations,
while the violation of the continuity equation gives a spurious contribution
to all modes of the density response. Both global and local particle-number
conservation can be restored by introducing a new density fluctuation that is
related to the current density by the continuity equation. This prescription
changes the strength associated with the various eigenmodes of the density
fluctuations, but not the eigenfrequencies of the system. We have shown in a
one-dimensional model that the energy-weighted sum rule calculated accord-
ing to this prescription has exactly the same value as for normal, uncorrelated
systems, thus we conclude that our prescription eliminates all the spurious
strength introduced by the constant-∆ approximation.
In a simplified model of nuclei, the effects of pairing correlations on the
isoscalar strength functions has been studied in detail for the quadrupole and
octupole channels, in the region of giant resonances. In both cases the effects
of pairing are rather small. More sizable effects are found at lower excitation
energy, in the region of surface modes which have not been included in the
present model, but will certainly be more affected by pairing correlations.
Appendix
In this Appendix we show that the correlated zero-order response function
given by the modified density fluctuation (99) satisfies the same energy-weighted
sum rule (EWSR) as the uncorrelated response function. We assume that par-
ticles move in one-dimensional square-well potential, so that formulae become
simpler because the particle velocity does not depend on position: v(ǫ, x) =
v(ǫ) =
Uncorrelated sum rule
The uncorrelated propagator [7]
D0(x, x′, ω) = (135)
dǫF ′(ǫ)
−2nω0
cos[nω0τ(x)]
v(ǫ, x)
ω − nω0(ǫ) + iε
cos[nω0τ(x
v(ǫ, x′)
gives the uncorrelated strength function S0(ω) = − 1
dxdx′Q(x)D0(x, x′, ω)Q(x′)
and the first moment
dωωS0(ω)
dǫF ′(ǫ)
(−2nω0(ǫ)
T (ǫ)
)(T (ǫ)
Q2nnω0(ǫ). (136)
Correlated sum rule
The modified density fluctuation (99) allows us to evaluate the correlated
propagator D̄0 through the relation
δ ¯̺(x, ω) = β
dx′D̄0(x, x′, ω)Q(x′) , (137)
giving
D̄0(x, x′, ω) = (138)
dǫF ′(ǫ)
−2nω0
)cosnω0τ(x)
v(ǫ, x)
ω − ω̄n + iε
cosnω0τ(x
v(ǫ, x′)
and the correlated first moment
M̄1 = 2
dǫF ′(ǫ)
(−2nω0(ǫ)
T (ǫ)
)(T (ǫ)
(nω0(ǫ)
ω̄n(ǫ)
ω̄n(ǫ) . (139)
The only difference between this expression and Eq. (136) is in the form of
F ′(ǫ), which is proportional to a δ-function in (136), while it is smoother in
the correlated case, however, if the parameter µ is determined by the one-
dimensional version of Eq. (27), then it can be easly found that, for a square-
well mean field,
M̄1 = M1 . (140)
The detailed argument goes as follows: both for correlated and uncorrelated
systems, the number of particles is given by
dxdpF (ǫ) =
dǫT (ǫ)F (ǫ) , (141)
with F (ǫ) = 4
θ(ǫF − ǫ) for uncorrelated fermions and F (ǫ) = 42π~ρ0(ǫ) in the
correlated case, while the moments (136, 139) are given by
dǫF ′(ǫ)G(ǫ) , (142)
G(ǫ) = −
(2π)2
T (ǫ)
Q2n (143)
in both cases, F ′(ǫ) obviously differs in the two cases.
Integrating by parts the last expression in (141), gives
A = J(ǫ)F (ǫ)
F ′(ǫ)J(ǫ) , (144)
J(ǫ) =
dǫT (ǫ) . (145)
For a square-well potential of size L:
T (ǫ) =
, (146)
J(ǫ) =
2mL 2
ǫ . (147)
Since
J(ǫ)F (ǫ) = lim
J(ǫ)F (ǫ) = 0 , (148)
both for the correlated and uncorrelated distributions, we have
A = −
dǫF ′(ǫ)J(ǫ) . (149)
The explicit expressions of J(ǫ) and T (ǫ) give
A=−2L
dǫF ′(ǫ)
ǫ , (150)
(2π)2√
dǫF ′(ǫ)
ǫ (151)
for both distributions. From these relations follows that, for a square-well mean
field, the relation (140) is exact. If we had used the density fluctuation (95),
instead of (99), to evaluate the correlated propagator, we would have obtained
a different value of the first moment because the additional term (106) in
the phase-space density gives an extra contribution to the density response
function and hence to the EWSR. Because of the fundamental character of
the EWSR, as well as of the continuity equation, we conclude that this term
is a spurious contribution generated by the constant-∆ approximation.
References
[1] M. Di Toro and V.M. Kolomietz, Zeit. Phys. A-Atomic Nuclei 328 (1987) 285
[2] M. Urban and P. Schuck, Phys. Rev. A 73 (2006) 013621
[3] J.R. Schrieffer, Theory of superconductivity, (W.A. Benjamin, Inc., New York,
1964)
[4] R. Combescot, M. Yu. Kagan, S. Stringari, Phys. Rev. A 74, (2006) 042717
[5] R. Bengtsson and P. Schuck, Phys. Lett. 89B (1980) 321
[6] P. Ring and P. Schuck, The Nuclear Many-Body Problem, (Springer, New
York, 1980)
[7] D.M. Brink, A. Dellafiore, M. Di Toro, Nucl. Phys. A456 (1986) 205
[8] V.L. Polyachenko, I.G. Schuckman, Sov. Astron. 25 (1981) 533
[9] D.M. Brink and G.R. Satchler, Angular Momentum (Oxford University Press,
Oxford,U.K., 1968), p. 147
[10] E.M. Lifshitz and L.P. Pitaevskii, Statistical Physics, Part 2 (Pergamon Press,
Oxford, 1980)
[11] A. Dellafiore, F. Matera, D.M. Brink, Phys. Rev. A 51 (1995) 914
[12] V.I. Abrosimov, A. Dellafiore, F. Matera, Nucl. Phys. A697 (2002) 748
[13] V.I. Abrosimov, O.I. Davidovskaya, A. Dellafiore, F. Matera, Nucl. Phys. A727
(2003) 748220

  The solutions of the Wigner-transformed time-dependent
Hartree--Fock--Bogoliubov equations are studied in the constant-$\Delta$
approximation. This approximation is known to violate particle-number
conservation. As a consequence, the density fluctuation and the longitudinal
response function given by this approximation contain spurious contributions. A
simple prescription for restoring both local and global particle-number
conservation is proposed. Explicit expressions for the eigenfrequencies of the
correlated systems and for the density response function are derived and it is
shown that the semiclassical analogous of the quantum single--particle spectrum
has an excitation gap of $2\Delta$, in agreement with the quantum result. The
collective response is studied for a simplified form of the residual
interaction.

Introduction
The problem of extending the Vlasov equation to systems in which pairing
correlations play an important role has been tackled some time ago by Di
Toro and Kolomietz [1] in a nuclear physics context and, more recently, by
Urban and Schuck [2] for trapped fermion droplets. These last authors derived
∗ Corresponding author
Email address: della@fi.infn.it (A. Dellafiore).
Preprint submitted to Elsevier 16 November 2021
http://arxiv.org/abs/0704.0152v2
the TDHFB equations for the Wigner transform of the normal density matrix
ρ and of the pair correlation function κ (plus their time-reversal conjugates)
and used them to study the dynamics of a spin-saturated trapped Fermi gas.
In the time-dependent theory one obtains a system of four coupled differential
equations for ρ, κ, and their conjugates [2] and, if one wants an analytical so-
lution, some approximation must be introduced. Here we try to find a solution
of the equations of motion derived by Urban and Schuck in the approximation
in which the pairing field ∆(r,p, t) is treated as a constant. It is well known
that such an approximation violates both particle-number-conservation and
gauge invariance (see e.g. sect. 8-5 of [3] and [4]), nonetheless we study it
because of its simplicity, with the aim of correcting the final results for its
shortcomings. Moreover, the constant-∆ approximation is not satisfactory for
describing long wavelength pairing modes in a large system. Such modes have
frequencies which are much less than the pairing frequency ∆/~ and for their
study it is essential to use a self consistent theory where the gap ∆ is related
to the pair density κ through the pairing interaction. The phases of ∆ and
κ are particularly important because they describe the collective superfluid
currents. On the other hand nuclei are small systems. Shell gaps are large
compared with ∆, or equivalently giant resonance frequencies are large com-
pared with the pairing frequency. The constant-∆ approximation is much more
reasonable in such systems.
In Sect. 2, the basic equations are recalled and reformulated in terms of the
even and odd components of the normal density ρ. In Sect. 3, the static limit
is studied by following the approach of [5] and the constant-∆ approximation
is introduced. In Sect. 4, the simplified dynamic equations resulting from the
constant-∆ approxmation are derived and their solutions are determined in lin-
ear approximation. In Sect. 5, these solutions are studied in a one-dimensional
model and the problem of particle-number conservation is examined in detail.
By studying the energy-weighted sum rule (in the Appendix), we find that the
constant-∆ approximation introduces some spurious strength into the density
response of the system. A simple prescription, based on the continuity equa-
tion, is proposed in order to eliminate the spurious strength. The resulting
strength function gives the same energy-weighted sum rule as for the uncor-
related systems. In Sect. 6, the general solution found in Sect. 4 is re-written
for spherical systems, where the angular integrations can be performed ex-
plicitly, leading to expressions containing only radial integrations. In Sect.
7, the collective response function of spherical nuclei is derived for a simple
multipole-multipole residual interaction. In Sect. 8, the quadrupole and oc-
tupole channels, that are the ones most affected by the pairing correlations
are shown explicitly. Finally, in Sect. 9 conclusions are drawn.
2 Basic equations
We assume that our system is saturated both in spin and isospin space and
do not distinguish between neutrons and protons, so we can use directly the
equations of motion of Urban and Schuck.
We start from the equations of motion derived in Ref. [2] for the Wigner-
transformed density matrices ρ = ρ(r,p, t) and κ = κ(r,p, t), with the warning
that the sign of κ that we are using agrees with that of Ref. [1], hence it is
opposite to that of [2]. Moreover we find convenient to use the odd and even
combinations of the normal density introduced in [2]:
ρev =
[ρ(r,p, t) + ρ(r,−p, t)] , (1)
ρod =
[ρ(r,p, t)− ρ(r,−p, t)] . (2)
Thus, the equations of motion given by Eqs.(15a...d) of Ref. [2] read
i~∂tρev = i~{h, ρod} − 2iIm[∆∗(r,p, t)κ] (3)
i~∂tρod = i~{h, ρev}+ i~Re{∆∗(r,p, t), κ} (4)
i~∂tκ=2(h− µ)κ−∆(r,p, t)(2ρev − 1) + i~{∆(r,p, t), ρod} . (5)
Here h is the Wigner-transformed Hartree–Fock hamiltonian h(r,p, t), while
∆(r,p, t) is the Wigner-transformed pairing field. Since the time-dependent
part of κ is complex, κ = κr+iκi, the last equation gives two separate equations
for the real and imaginary parts of κ.
Moreover, from the supplementary normalization condition ([6], p. 252)
R2 = R (6)
satisfied by the generalized density matrix R, the two following independent
equations are obtained:
ρodκ+ i
{ρev, κ}=0 , (7)
ρev(ρev − 1) + ρ2od + κκ∗=0 . (8)
We shall use the equations of motion (3–5), together with these equations, as
our starting point, but first we notice that, in the limit of no pairing, both ∆
and κ vanish, the third equation of motion reduces to a trivial identity, while
the first two give the Vlasov equation for normal systems, expressed in terms
of the even and odd components of ρ:
∂tρev = {h, ρod} , (9)
∂tρod = {h, ρev} . (10)
A solution of the linearized Vlasov equation for normal systems (i. e. without
pairing) has been obtained in Ref. [7] and our aim here is to study the changes
introduced by the pairing interaction in the solution of [7].
Moreover, before studying the time-dependent problem, it is useful to look at
the static limit.
3 Static limit
In this section we follow the approach of Ref. [5]. At equilibrium we have
ρev = ρ0(r,p) , (11)
ρod =0 , (12)
h=h0(r,p) , (13)
κ=κ0(r,p) (14)
∆0 =∆0(r,p) (15)
and equations (3–5) give
0=−2iIm(∆∗0κ0) , (16)
0= i~{h0, ρ0}+ i~Re{∆∗0, κ0} , (17)
0=2(h0 − µ)κ0 −∆0(2ρ0 − 1) , (18)
while Eqs. (7, 8) give
{ρ0, κ0}=0 , (19)
ρ0(ρ0 − 1) + |κ0|2=0 . (20)
Equation (16) is satisfied if we assume that ∆0 and κ0 are real quantities,
while Eqs. (18) and (20), taken as a system, have the solution:[5]
ρ0(r,p) =
h0(r,p)− µ
E(r,p)
κ0(r,p) = −
∆0(r,p)
2E(r,p)
, (22)
with the quasiparticle energy
E(r,p) =
∆20(r,p) + (h0(r,p)− µ)2 . (23)
It can be easily checked that Eqs. (21, 22) satisfy also Eqs. (17) and (19), that
{h0, ρ0}+ {∆0, κ0} = 0 (24)
{ρ0, κ0} = 0 . (25)
The (semi)classical equilibrium phase-space distribution is closely related to
ρ0(r,p):
f0(r,p) =
(2π~)3
ρ0(r,p) (26)
and the statistical factor 4 takes into account the fact that there are two kinds
of fermions.
The parametrer µ is determined by the condition
drdpf0(r,p) , (27)
where A is the number of particles. This integral should keep the same value
also out of equilibrium (global particle-number conservation).
3.1 Constant-∆ approximation
In a fully self-consistent approach, the pairing field ∆(r,p, t) is related to
κ(r,p, t), however here we introduce an approximation and replace the pairing
field of the HFB theory with the phenomenological pairing gap of nuclei, hence
in all our equations we put
∆(r,p, t) ≈ ∆0(r,p) ≈ ∆ = const. , (28)
with ∆ ≈ 1MeV.
In the constant-∆ approximation the equilibrium distributions become
ρ0(ǫ) =
1− ǫ− µ
, (29)
κ0(ǫ) =−
2E(ǫ)
and the quasiparticle energy
E(ǫ) =
∆2 + (ǫ− µ)2 , (31)
ǫ = h0(r,p) =
+ V0(r)
the particle energy in the equilibrium mean field.
In the following we shall use the relation:
κ0(ǫ) =
E2(ǫ)
dρ0(ǫ)
. (32)
4 Dynamic equations
Always in the approximation where ∆ is constant and real, the time-dependent
equations (3–5) become
i~∂tρev = i~{h, ρod} − 2i∆Im(κ) (33)
i~∂tρod = i~{h, ρev} (34)
i~∂tκ=2(h− µ)κ−∆(2ρev − 1) . (35)
This is the simplified set of equations that we want to study here. The sum of
the first two equations gives an equation that is similar to the Vlasov equation
of normal systems, only with the extra term −2i∆Im(κ). This extra term
couples the equation of motion of ρ with that of κ, thus, instead of a single
differential equation (Vlasov equation), now we have a system of two coupled
differential equations (for ρ and κi).
Our aim here is that of determining the effects of pairing on the linear response
of nuclei, thus we assume that our system is initially at equilibrium, with
densities given by Eqs. (29,30), and that at time t = 0 a weak external field
of the kind
δV ext(r, t) = βδ(t)Q(r) (36)
is applied to it. This simple time-dependence is sufficient to determine the
linear response of the system. In a self-consistent approach, we should take
into account also the changes of the mean field surrounding each particle
induced by the external force and consider a perturbing hamiltonian of the
δh = δV ext + δV int , (37)
however we start with the zero-order approximation
δh = δV ext (38)
and will consider collective effects in a second stage.
Since we want to solve Eqs. (33–35) in linear approximation, we consider small
fluctuations of the time-dependent quantities about their equilibrium values
and neglect terms that are of second order in the fluctuations. Hence, in Eqs.
(33–35) we put:
h=h0 + δh , (39)
ρev = ρ0 + δρev , (40)
ρod = δρod , (41)
κ=κ0 + δκ = κ0 + δκr + iδκi . (42)
Then, the linearized form of Eqs. (33–35) is
i~∂tδρev = i~{h0, δρod} − 2i∆δκi (43)
i~∂tδρod= i~{h0, δρev}+ i~{δh, ρ0} (44)
−~∂tδκi=2(ǫ− µ)δκr + 2κ0δh− 2∆δρev , (45)
~∂tδκr =2(ǫ− µ)δκi . (46)
Taking the sum of the first two equations gives
i~∂tδρ(r,p, t) = (47)
i~{h0, δρ(r,p, t)}+ i~{δh(r,p, t), ρ0} − 2i∆δκi(r,p, t) ,
which can be regarded as an extension of the linearized Vlasov equation stud-
ied in [7].
In order to make the comparison with [7] easier, from now on we change the
normalization of the phase-space densities and define
f(r,p, t)=
(2π~)3
ρ(r,p, t) (48)
χ(r,p, t)=
(2π~)3
κ(r,p, t) , (49)
moreover, we put
F (ǫ) =
(2π~)3
ρ0(ǫ) . (50)
In terms of the new functions Eqs. (47) and (45) read
i~∂tδf(r,p, t) = i~{h0, δf(r,p, t)}+ i~{δh(r,p, t), f0}
− 2i∆δχi(r,p, t) , (51)
−~∂tδχi(r,p, t) = 2(ǫ− µ)δχr(r,p, t) + 2χ0δh(r,p, t)
− 2∆δfev(r,p, t) . (52)
The function fev is given by the obvious extension of Eq. (1). In order to get
a closed system of equations, we still need an extra equation for δχr(r,p, t).
This can be obtained from the linearized form of the supplementary condition
(8) that reads
δρev(2ρ0 − 1) = −2κ0δκr , (53)
δκr(r,p, t) =
1− 2ρ0(ǫ)
δρev(r,p, t) = −
δρev(r,p, t) . (54)
The last expression has been obtained with the help of Eq. (18). In terms of
the new functions f and χ, the last equation reads
δχr(r,p, t) = −
δfev(r,p, t) . (55)
Equations (51, 52) and (55) are the set of coupled equations for the phase-
space densities that we have to solve.
Replacing Eq. (55) into Eq. (52), and using Eq. (30), gives the following system
of coupled differential equations:
∂tδf(r,p, t)= {h0, δf}+ {δh, f0} − 2
δχi(r,p, t) , (56)
∂tδχi(r,p, t)=
E2(ǫ)
[δf(r,p, t) + δf(r,−p, t)]− 2
δh(r,p, t) . (57)
Taking the Fourier transform in time, gives
−iωδf(r,p, ω)= {h0, δf}+ {δh, f0} − 2
δχi(r,p, ω) , (58)
−iωδχi(r,p, ω))=
E2(ǫ)
[δf(r,p, ω) + δf(r,−p, ω)]
− 2χ0
δh(r,p, ω) , (59)
or (for ω 6= 0)
−iωδf(r,p, ω) + {δf, h0}=−iωd2
[δf(r,p, ω) + δf(r,−p, ω)]
+F ′(ǫ)[{δh, h0}+ iωd2δh] , (60)
(Ω(ǫ)
Ω(ǫ) = 2
. (62)
This frequency plays a crucial role in our approach, its minimum value is
2∆/~.
In Eq. (60) we have used the relation {f0, δh} = F ′(ǫ){ho, δh} as well as Eq.
(32).
By comparing Eq. (60) with the analogous equation for normal systems
−iωδf(r,p, ω) + {δf, h0} = F ′(h0){δh, h0} , (63)
we can see that the only effect of pairing in the constant-∆ approximation is
that of adding the terms proportional to d2.
The normal Vlasov equation (63) can be solved in a very compact way by using
the method of action-angle variables [8], [7]. In that approach one expands
δhn(I)e
in·Φ , (64)
where I and Φ are the action and angle variables, respectively. Moreover
δf(r,p, ω) =
δfn(I, ω)e
in·Φ (65)
{δf, h0} =
i(n · ~ω)δfn(I, ω)ein·Φ , (66)
where the vector ~ω has components
. (67)
Then Eq.(63) gives
δfn(I, ω) =
(n · ~ω)
(n · ~ω)− (ω + iε)
F ′(ǫ)δhn(I) . (68)
The (zero-order) eigenfrequencies of the (normal) physical system are
ωn = n · ~ω . (69)
Here want to use the same method to solve the more complicated equation
(60). Since that equation contains also the function δf(r,−p, ω), we need also
the analogous equation for this other quantity:
−iωδf(r,−p, ω)− {δf, h0}r,p=−iωd2
[δf(r,p, ω) + δf(r,−p, ω)]
+F ′(ǫ)[−{δh, h0}r,p + iωd2δh] . (70)
By expanding δf(r,p, ω) and δf(r,−p, ω) as
δf(r,±p, ω) =
δf±n (I, ω)e
in·Φ , (71)
Eqs. (60) and (70) give
[−ω(1− d
) + ωn]δf
n + ω
δf−n = F
′(ǫ)[ωn + ωd
2]δhn , (72)
]δf+n + [−ω(1−
)− ωn]δf−n = F ′(ǫ)[−ωn + ωd2]δhn , (73)
which is a system of two coupled algebraic equations for the coefficients δf+n
and δf−n . Its solution is
δf+n =
ω̄2n + ωωn
ω̄2n − ω2
F ′(ǫ)δhn , (74)
δf−n =
ω̄2n − ωωn
ω̄2n − ω2
F ′(ǫ)δhn , (75)
where
ω̄2n = ω
n + Ω
2(ǫ) (76)
are the (squared) eigenfrequencies of the correlated system. These eigenfre-
quencies are in agreement with the enegy spectrum of a superfluid infinite ho-
mogeneous Fermi gas (see e. g. Sect. 39 of [10]) and they lead to a low-energy
gap of 2∆ in the excitation spectrum of the correlated systems. However,
as anticipated, we expect problems with particle-number conservation. These
problems are better discussed in one dimension, where formulae are simpler.
5 One-dimensional systems and particle-number conservation
In one dimension, Eq.(60) reads
−iωδf(x, p, ω) + ẋ∂xδf −
dV0(x)
∂pδf = (77)
−iωd2 1
[δf(x, p, ω) + δf(x,−p, ω)] + F ′(ǫ)( p
∂xδh + iωd
2δh) .
In zero-order approximation δh(x, p, ω) = βQ(x), moreover, in one dimension
F ′(ǫ) =
. (78)
The vectors ~ω and Φ have only one component:
ω0(ǫ) =
T (ǫ)
Φ(x) = ω0τ(x) , (80)
τ(x) =
v(ǫ, x′)
, (81)
v(ǫ, x) =
[ǫ− V0(x)] . (82)
The time T (ǫ) is the period of the bound motion of particles with enegy ǫ
in the equilibrium potential well V0(x): T = 2τ(x2). The points x1,2 are the
classical turning points for the same particles. Instead of the action variable
I(ǫ) = 1
dxp(ǫ, x), it is more convenient to use the particle energy ǫ as
constant of motion. As pointed out in [7], the range of values of τ can be
extended to the whole interval (0, T ), by defining
τ(x) =
v(ǫ, x′)
when τ > T
. With this extension, the angle variable Φ(x) takes values between
0 and 2π, as it should.
In one dimension, Eqs. (74, 75) give
δf±n (ǫ, ω) =
ω̄2n ± ω ωn
ω̄2n − ω2
F ′(ǫ)δhn , (84)
δhn = βQn
v(ǫ, x)
e−iωnτ(x)
v(ǫ, x)
cos[ωnτ(x)] , (85)
The frequencies ωn are the eigenfrequencies of the uncorrelated system:
ωn = nω0 , (86)
while ω̄n are the new eigenfrequencies modified by the pairing correlations:
ω̄n = ±
ω2n + Ω
2(ǫ) . (87)
Note that, since δh−n = δhn, then δf
n = δf
By using the solutions (84), we can also obtain an expansion for the even and
odd parts of δf :
δfev(x, ǫ, ω) =
An(ω) cosnω0τ(x) , (88)
δfod(x, ǫ, ω) =
Bn(ω) sinnω0τ(x) , (89)
An(ω) =
ω̄2n − ω2
F ′(ǫ)δh′n , (90)
Bn(ω) = iω
ω̄2n − ω2
F ′(ǫ)δh′n (91)
δh′n =2δhn , n 6= 0 , (92)
= δhn , n = 0 . (93)
Note that, while Bn=0(ω) = 0, we have An=0(ω) 6= 0, and this fact leads to
an unphysical fluctuation of the number of particles, induced by the applied
external field. These fluctuations are given by
δA(ω) =
dxδ̺(x, ω) , (94)
where δ̺(x, ω) is the density fluctuation at point x:
δ̺(x, ω) =
dp δf(x, p, ω) = 2
v(ǫ, x)
δfev(x, ǫ, ω) . (95)
Equation (88) gives
δA(ω) = 2
An(ω)
dτ cos nω0τ . (96)
Since the integrals
∫ T/2
0 dτ cosnω0τ vanish when n 6= 0, the term with n = 0
is the only one contributing to this sum, thus givig an unphysical fluctuation
of the number of particles. This problem could be solved simply by excluding
the mode n = 0 from the sum in Eq. (88), however this would not be sufficient
to solve all problems with particle-number conservation, since we can easily
check that the solutions (88, 89) do not satisfy the continuity equation
iω̺(x, ω) = ∂xj(x, ω) . (97)
The density fluctuation involves only the even part of δf , while the current
density j(x, ω) involves only the odd part:
j(x, ω) =
δf(x, p, ω) = 2
dǫδfod(x, ǫ, ω) . (98)
The fact that the continuity equation is violated is a very serious shortcom-
ing of the constant-∆ approximation. However, since we have seen that this
approximation leads to very simple equations and to rather satisfactory ex-
pressions for the eigenfrequencies of the correlated systems, we still use it,
but with the following prescription: when calculating the longitudinal response
function, the density fluctuations should be evaluated by using Eq. (97), instead
of Eq, (95). Then, the density fluctuations (95) should be replaced by
δ ¯̺(x, ω) =
dǫ∂xδfod(x, ǫ, ω) . (99)
In practice we are proposing to evaluate the longitudinal response function
in terms of the transverse response function. It is well known that also the
more familiar BCS approximation gives a more accurate description of the
transverse response (see e.g. sect. 8-5 of [3]). In the Appendix we show that the
longitudinal response function resulting from the present prescription satisfies
the same energy-weighted sum rule as the uncorrelated response function. This
would not necessarily happen if, instead of changing only the even part of δf ,
we had modified also its odd part.
It is interesting to see how the solutions (84) are changed by our prescription.
By using Eq. (89) forfod, Eq. (99) gives
δ ¯̺(x, ω) = 2
v(ǫ, x)
δf̄ev(x, ǫ, ω) , (100)
δf̄ev(x, ǫ, ω) =
Ān(ω) cosnω0τ(x) , (101)
Ān(ω) =
Bn(ω) , (102)
(note that Ān=0(ω) = 0). Then
δf̄(x,±p, ω) = δf̄ev(x, ǫ, ω)± δfod(x, ǫ, ω) (103)
and Eq. (84) is replaced by
δf̄±n (ǫ, ω) =
ω2n ± ω ωn
ω̄2n − ω2
F ′(ǫ)δhn . (104)
By comparing this expression to Eq. (84), we can see that the fluctuations
of the phase-space density given by the constant-∆ approximation contain an
extra contribution that we identify as spurious:
δf±n (ǫ, ω) = δf̄
n (ǫ, ω) + δf
n (ǫ, ω) , (105)
δf spurn (ǫ, ω) =
Ω2(ǫ)
ω̄2n − ω2
F ′(ǫ)δhn . (106)
The spurious character of δf spurn is suggested also by sum-rule arguments
(see Appendix). The term f spurn (ǫ, ω) contributes to all modes of the density
strength function: the contribution to the mode n = 0 gives a fluctuation of
the particle-number integral (global paticle-number violation), while the other
modes give a spurious contribution to the density strength function, increas-
ing the sum rule and violating the continuity equation (local particle-number
violation). Note that the spurious contribution (106) affects only the even part
of the pase-space density, not the odd part.
6 Spherical Systems
The method of action-angle variables gives a very compact solution of the lin-
earized Vlasov equation both in the uncorrelated and correlated cases, however
it may be useful to make a connection between the results given by this method
and the more explicit treatment of spherical nuclei given in [7]. For uncorre-
lated system this has been done in [11]. Here we follow that approach in order
to derive useful expressions for correlated spherical systems. The components
of the vecor n are (n1, n2, n3), the first point to notice is that, because of the
degeneracy associated with any central-force field, the vector ~ω has only two
non vanishing components:
~ω = (0, ωϕ(ǫ, λ), ω0(ǫ, λ)) . (107)
With λ we denote the magnitude of the particle angular momentum. According
to Eq. (69), the eigenfrequencies of the uncorrelated system are [7]
ωn = ωn3,n2(ǫ, λ) = n3ω0 + n2ωϕ , (108)
while Eq. (76) gives the correlated eigenfrequencies
ω̄n = ω̄n3,n2 = ±
ω2n3,n2 + Ω
2(ǫ) . (109)
In three dimensions, the Fourier coefficients analogous to (85) are
Qn(I) =
(2π)3
dΦe−in·ΦQ(r) . (110)
The external field Q(r) can be expanded in partial waves as
Q(r) =
QL(r)YLM(r̂) , (111)
giving
Qn(I) =
Q(LM)n , (112)
with [11]
Q(LM)n =
dLMN(β
′)δM,n1δN,n2Q
. (113)
By using this last equation (and changing n3 → n), the expansion (71) be-
comes
δf(r,±p, ω) =
(114)
δfL±nN (ǫ, λ, ω)e
iφnN (r)
DLMN(α, β ′, γ)
δfL±nN (ǫ, λ, ω) =
ω̄2nN ± ωωnN
ω̄2nN − ω2
βF ′(ǫ)QLnN (115)
and QLnN the semiclassical limit of the radial matrix elements:
QLnN =
vr(r)
e−iφnN (r)QL(r) ,
vr(r)
cos[φnN(r)]QL(r) . (116)
Here T is the period of radial motion, vr(r) the radial velocity
vr(r) =
ǫ− V0(r)−
(117)
and the phases φnN(r) are given by
φnN(r) = ωnNτ(r)−Nγ(r) , (118)
where
τ(r) =
vr(r′)
(119)
γ(r) =
vr(r′)
. (120)
The frequencies ω0 and ωϕ are given by
τ(r2)
, (121)
γ(r2)
τ(r2)
. (122)
The Wigner rotation matrix elements in Eq. (114) are given by [9]
DLMN(α, β ′, γ) = e−iMαdLMN(β ′)e−iNγ , (123)
where (α, β ′, γ) are the Euler angles introduced in [7].
On the basis of the discussion in Sect. (5), we expect that the solution (114)
will contain some spurious strength introduced by the constant-∆ approxima-
tion. In order to eliminate the spurious contributions, we should replace the
coefficients (115) with
δf̄L±nN (ǫ, λ, ω) =
ω2nN ± ωωnN
ω̄2nN − ω2
βF ′(ǫ)QLnN . (124)
These modified coefficients allow us to obtain the modified zero-order propa-
gator
D̄0L(r, r
′, ω) =
dǫF ′(ǫ)
d̄LnN(r, r
ω − ω̄nN + iε
, (125)
d̄LnN(r, r
′) = (126)
(4π)2
2L+ 1
|YLN(
(−2ωnN
)(ωnN
)cosφnN(r)
r2vr(r)
cosφnN(r
r′2vr(r′)
and the corresponding response and strength functions:
R̄0L(ω)=
drdr′r2QL(r)D̄
L(r, r
′, ω)r′2QL(r
′) , (127)
S̄0L(ω)=−
ImR̄0L(ω) . (128)
For multipole response: QL(r) = r
For normal systems, the zero-order propagator D0L(r, r
′, ω) is given by Eqs.
(125) and(126) where ω̄nN is replaced by ωnN and F
′(ǫ) is proportional to a
δ-function[7].
7 Collective response
Up to now, we have been concerned only with the zero-order approximation,
which corresponds to the single-particle approximation of the quantum ap-
proach. In this approximation, the perturbing part of the hamiltonian is given
only by the external field, while a more consistent approach would require
taking into account also the mean-field fluctuation induced by the external
force, so that
δh = δV ext(r, ω) + δV int(r, ω) . (129)
In the Hartree approximation,
δV int(r, ω) =
dr′v(r− r′)δ̺(r′, ω) . (130)
where v(r− r′)is the (long-range) interaction between constituents.
For consistency, we take
δV int(r, ω) =
dr′v(r− r′)δ ¯̺(r′, ω) , (131)
then the collective propagator for correlated systems satisfies the same kind
of integral equation as for normal systems [7]:
D̄L(r, r
′, ω) = (132)
D̄0L(r, r
′, ω) +
dyy2D̄0L(r, x, ω)vL(x, y)D̄L(y, r
′, ω) .
Here vL(x, y) is the partial-wave component of the interaction between par-
ticles. We assume that this interaction can be approximated by a separable
form of the kind
vL(x, y) = κLx
LyL , (133)
where κL is a parameter that determines the strength of the interaction. Then,
the integral equation (132) gives an algebraic equation for the collective cor-
related response function
R̄L(ω) =
drr2rL
dr′r′2r′LD̄L(r, r
′, ω)
leading to the expression
R̄L(ω) =
R̄0L(ω)
1− κLR̄0L(ω)
. (134)
8 Results
Here we compare the multipole strength functions given by our simplified
model of pairing correlations with that of the corresponding uncorrelated sys-
tem. This comparison is made for the quadrupole and octupole strength func-
tions, since these channels are the ones that are most affected.
The static nuclear mean field is approximated with a spherical cavity of radius
R = 1.2A
3 fm and the A nucleons are treated on the same footing, i. e., we do
not distinguish between neutrons and protons. Moreover, we chose A = 208
for ease of comparison with previous calculations of uncorrelated response
functions [12,13]. Shell effects are not included in our semiclassical picture and
the results shown below should be considered as an indication of the qualitative
effects to be expected in heavy nuclei. For the uncorrelated calculations, the
Fermi energy is determined by the parametrization chosen for the radius as
ǫF ≈ 33.33 MeV, while for the correlated case, the parameter µ is determined
by the condition (27); with the value of ∆ = 1 MeV used here, the value of µ
is practically coincident with that of ǫF , so we have used µ = ǫF = 33.33 MeV
in the calculations below. Moreover, the small parameter ε appearing in Eq.
(125) has been given the value ε = 0.1 MeV. This value is chosen to simplify
the evaluation of the response function by smoothing out discontinuities in its
dependence on ω.
In the evaluation of the collective response, the value of parameters κL is
the same as in [12,13], that is: κ2 = −1 × 10−3 MeV/fm4 and κ3 = −2 ×
10−5 MeV/fm6.
8.1 Quadrupole response
Figure 1 show the longitudinal quadrupole strength function evaluated in the
zero-order approximation (corresponding to the quantum single-particle ap-
proximation). The dashed curve shows the uncorrelated response evaluated
according to the theory of [7], while the full curve shows the result of the
present correlated calculation. As we can see the effect of pairing correlations
on this zero-order strength function is rather small, however, since pairing af-
fects also the real part of the zero-order response function,in Fig. 2 we plot
0 5 10 15 20 25 30
ÑΩ @MeVD
Fig. 1. Quadrupole stength function in zero-order approximation. The dashed curve
gives the response of a normal system of A = 208 nucleons contained in a spherical
cavity, while the solid curve includes the effects of pairing correlations in constant-∆
approximation.
0 5 10 15 20 25 30
ÑΩ @MeVD
Fig. 2. Collective quadrupole strength function showing the giant quadrupole res-
onance. The solid curve involves also pair correlations, the dashed curve has no
pairing.
also the collective strength function given by Eq. (134). Again, the effect of
pairing correlations is very small, in agreement with the results of [1].
0 1 2 3 4 5
ÑΩ @MeVD
0.0025
0.005
0.0075
0.0125
0.015
0.0175
Fig. 3. Same as Fig.1, at low excitation energy.
The main difference between the uncorrelated and correlated responses occurs
at small excitation energy, Fig. 3 shows a detail of Fig. 1 at low excitation
energy. The correlated strength function displays a gap of about 2 ∆ , the
very small strength extending below 2 MeV is entirely due to the finite value
of the small parameter ε used in the numerical evaluation of the propagator
(125).
8.2 Octupole response
Figures 4 and 5 show the zero-order and collective octupole strength functions,
both correlated and uncorrelated. As we can see, in this case too the effect is
rather small.
0 5 10 15 20 25 30
ÑΩ @MeVD
Fig. 4. Zero-order octupole strength function. The solid curve involves pair correla-
tions, the dashed curve has no pairing.
0 5 10 15 20 25 30
ÑΩ @MeVD
Fig. 5. Collective octupole strength function. The solid curve involves pair correla-
tions, the dashed curve has no pairing.
9 Conclusions
The solutions of the semiclassical time-dependent Hartree–Fock–Bogoliubov
equations have been studied in a simplified model in which the pairing field
∆(r,p, t) is treated as a constant phenomenological parameter. Such an ap-
proximation is known to violate some important constraints, like global (particle-
number integral) and local (continuity equation) particle-number conserva-
tion. In a linearized approach , we have shown that the global particle-number
violation is related only to one particular mode of the density fluctuations,
while the violation of the continuity equation gives a spurious contribution
to all modes of the density response. Both global and local particle-number
conservation can be restored by introducing a new density fluctuation that is
related to the current density by the continuity equation. This prescription
changes the strength associated with the various eigenmodes of the density
fluctuations, but not the eigenfrequencies of the system. We have shown in a
one-dimensional model that the energy-weighted sum rule calculated accord-
ing to this prescription has exactly the same value as for normal, uncorrelated
systems, thus we conclude that our prescription eliminates all the spurious
strength introduced by the constant-∆ approximation.
In a simplified model of nuclei, the effects of pairing correlations on the
isoscalar strength functions has been studied in detail for the quadrupole and
octupole channels, in the region of giant resonances. In both cases the effects
of pairing are rather small. More sizable effects are found at lower excitation
energy, in the region of surface modes which have not been included in the
present model, but will certainly be more affected by pairing correlations.
Appendix
In this Appendix we show that the correlated zero-order response function
given by the modified density fluctuation (99) satisfies the same energy-weighted
sum rule (EWSR) as the uncorrelated response function. We assume that par-
ticles move in one-dimensional square-well potential, so that formulae become
simpler because the particle velocity does not depend on position: v(ǫ, x) =
v(ǫ) =
Uncorrelated sum rule
The uncorrelated propagator [7]
D0(x, x′, ω) = (135)
dǫF ′(ǫ)
−2nω0
cos[nω0τ(x)]
v(ǫ, x)
ω − nω0(ǫ) + iε
cos[nω0τ(x
v(ǫ, x′)
gives the uncorrelated strength function S0(ω) = − 1
dxdx′Q(x)D0(x, x′, ω)Q(x′)
and the first moment
dωωS0(ω)
dǫF ′(ǫ)
(−2nω0(ǫ)
T (ǫ)
)(T (ǫ)
Q2nnω0(ǫ). (136)
Correlated sum rule
The modified density fluctuation (99) allows us to evaluate the correlated
propagator D̄0 through the relation
δ ¯̺(x, ω) = β
dx′D̄0(x, x′, ω)Q(x′) , (137)
giving
D̄0(x, x′, ω) = (138)
dǫF ′(ǫ)
−2nω0
)cosnω0τ(x)
v(ǫ, x)
ω − ω̄n + iε
cosnω0τ(x
v(ǫ, x′)
and the correlated first moment
M̄1 = 2
dǫF ′(ǫ)
(−2nω0(ǫ)
T (ǫ)
)(T (ǫ)
(nω0(ǫ)
ω̄n(ǫ)
ω̄n(ǫ) . (139)
The only difference between this expression and Eq. (136) is in the form of
F ′(ǫ), which is proportional to a δ-function in (136), while it is smoother in
the correlated case, however, if the parameter µ is determined by the one-
dimensional version of Eq. (27), then it can be easly found that, for a square-
well mean field,
M̄1 = M1 . (140)
The detailed argument goes as follows: both for correlated and uncorrelated
systems, the number of particles is given by
dxdpF (ǫ) =
dǫT (ǫ)F (ǫ) , (141)
with F (ǫ) = 4
θ(ǫF − ǫ) for uncorrelated fermions and F (ǫ) = 42π~ρ0(ǫ) in the
correlated case, while the moments (136, 139) are given by
dǫF ′(ǫ)G(ǫ) , (142)
G(ǫ) = −
(2π)2
T (ǫ)
Q2n (143)
in both cases, F ′(ǫ) obviously differs in the two cases.
Integrating by parts the last expression in (141), gives
A = J(ǫ)F (ǫ)
F ′(ǫ)J(ǫ) , (144)
J(ǫ) =
dǫT (ǫ) . (145)
For a square-well potential of size L:
T (ǫ) =
, (146)
J(ǫ) =
2mL 2
ǫ . (147)
Since
J(ǫ)F (ǫ) = lim
J(ǫ)F (ǫ) = 0 , (148)
both for the correlated and uncorrelated distributions, we have
A = −
dǫF ′(ǫ)J(ǫ) . (149)
The explicit expressions of J(ǫ) and T (ǫ) give
A=−2L
dǫF ′(ǫ)
ǫ , (150)
(2π)2√
dǫF ′(ǫ)
ǫ (151)
for both distributions. From these relations follows that, for a square-well mean
field, the relation (140) is exact. If we had used the density fluctuation (95),
instead of (99), to evaluate the correlated propagator, we would have obtained
a different value of the first moment because the additional term (106) in
the phase-space density gives an extra contribution to the density response
function and hence to the EWSR. Because of the fundamental character of
the EWSR, as well as of the continuity equation, we conclude that this term
is a spurious contribution generated by the constant-∆ approximation.
References
[1] M. Di Toro and V.M. Kolomietz, Zeit. Phys. A-Atomic Nuclei 328 (1987) 285
[2] M. Urban and P. Schuck, Phys. Rev. A 73 (2006) 013621
[3] J.R. Schrieffer, Theory of superconductivity, (W.A. Benjamin, Inc., New York,
1964)
[4] R. Combescot, M. Yu. Kagan, S. Stringari, Phys. Rev. A 74, (2006) 042717
[5] R. Bengtsson and P. Schuck, Phys. Lett. 89B (1980) 321
[6] P. Ring and P. Schuck, The Nuclear Many-Body Problem, (Springer, New
York, 1980)
[7] D.M. Brink, A. Dellafiore, M. Di Toro, Nucl. Phys. A456 (1986) 205
[8] V.L. Polyachenko, I.G. Schuckman, Sov. Astron. 25 (1981) 533
[9] D.M. Brink and G.R. Satchler, Angular Momentum (Oxford University Press,
Oxford,U.K., 1968), p. 147
[10] E.M. Lifshitz and L.P. Pitaevskii, Statistical Physics, Part 2 (Pergamon Press,
Oxford, 1980)
[11] A. Dellafiore, F. Matera, D.M. Brink, Phys. Rev. A 51 (1995) 914
[12] V.I. Abrosimov, A. Dellafiore, F. Matera, Nucl. Phys. A697 (2002) 748
[13] V.I. Abrosimov, O.I. Davidovskaya, A. Dellafiore, F. Matera, Nucl. Phys. A727
(2003) 748220
