text:
# References

Deformed harmonic oscillators for metal clusters: Analytic properties and supershells

Dennis Bonatsos#111e-mail: bonat@mail.demokritos.gr, D. Lenis#, P. P. Raychev222e-mail: raychev@inrne.bas.bg, and P. A. Terziev333e-mail: terziev@inrne.bas.bg

# Institute of Nuclear Physics, N.C.S.R. “Demokritos”

GR-15310 Aghia Paraskevi, Attiki, Greece

Institute for Nuclear Research and Nuclear Energy, Bulgarian Academy of Sciences

72 Tzarigrad Road, BG-1784 Sofia, Bulgaria

Abstract

The analytic properties of Nilsson’s Modified Oscillator (MO), which was first introduced in nuclear structure, and of the recently introduced, based on quantum algebraic techniques, 3-dimensional $`q`$-deformed harmonic oscillator (3-dim $`q`$-HO) with uq(3) $``$ soq(3) symmetry, which is known to reproduce correctly in terms of only one parameter the magic numbers of alkali clusters up to 1500 (the expected limit of validity for theories based on the filling of electronic shells), are considered. Exact expressions for the total energy of closed shells are determined and compared among them. Furthermore, the systematics of the appearance of supershells in the spectra of the two oscillators is considered, showing that the 3-dim $`q`$-HO correctly predicts the first supershell closure in alkali clusters without use of any extra parameter.

PACS numbers: 36.40.-c, 36.40.Cg, 03.65.Fd

1. Introduction

Supershells, which are seen as beats of the deviation of the total energy of many particle systems from the function describing its average behavior vs. the number of particles $`N`$, are known to be a general property of the spectrum of potentials having sharp edges . Supershells in metal clusters have first been studied by Nishioka, Hansen and Mottelson  in terms of phenomenological mean field potentials.

In the other hand, using recently developed quantum algebraic techniques , it has been shown  that the magic numbers appearing in alkali clusters can be successfully reproduced up to 1500 (which is the expected limit of validity of theories based on the filling of electronic shells ) by the 3-dimensional $`q`$-deformed harmonic oscillator (3-dim $`q`$-HO), which possesses the uq(3) $``$ soq(3) symmetry . Furthermore, the magic numbers appearing in several divalent (Zn, Cd) and trivalent (Al, In) metal clusters have been satisfactorily reproduced  by the same model, in terms of only one free parameter, the deformation parameter $`\tau `$ (with $`q=e^\tau `$, where $`\tau `$ is a real number). It is therefore of interest to examine if the 3-dim $`q`$-HO can predict supershells and which these predictions are. It should be noticed that the calculation of supershells in the framework of the 3-dim $`q`$-HO will be parameter free, since the single parameter of the model has been fixed in reproducing the magic numbers for each kind of clusters .

In addition to the 3-dim $`q`$-HO, Nilsson’s Modified Oscillator (MO) , which has first been used in describing the structure of atomic nuclei, has also been early employed in describing atomic clusters  (after dropping the spin–orbit interaction, which plays an essential role in nuclear structure but is absent in the case of atomic clusters). It is therefore of interest to study the possible appearance of supershells in the framework of this model as well.

For the determination of supershells the method of Ref.  can be employed. Before doing so, one has however to examine the analytic properties of the spectra of the two oscillators, in order to be able to apply meaningful truncation schemes. Furthermore, the average behavior of the total energy of a system of many particles (an atomic cluster in the present case) as a function of the particle number $`N`$ should be determined, since it is needed in the procedure of the study of supershells. As a result, these tasks will be carried out for both models, before any attempt for the determination of supershells is made.

Supershells have been early predicted by Balian and Bloch  in the study of electrons moving in a spherical cavity, which by analogy can be used for the valence electrons in a metal cluster. The comparison of the stringent predictions of the theory of Balian and Bloch for various characteristics of the supershells  to the results of the present models turns out to be a fruitful testing procedure.

In Section 2 the analytic properties of Nilsson’s Modified Oscillator will be considered, while the corresponding properties of the 3-dim $`q`$-HO will be studied in Section 3. In Section 4 supershells in Nilsson’s Modified Oscillator will be studied, while supershells in the framework of the 3-dim $`q`$-HO will be considered in Section 5. Finally in Section 6 a discussion of the present results and plans for further work will be given.

2. Nilsson’s Modified Oscillator (MO)

The potential of the Modified Oscillator (MO) introduced in nuclear physics by Nilsson  (with the spin-orbit term omitted) is

$$V=\frac{1}{2}\mathrm{}\omega \rho ^2\mathrm{}\omega \mu ^{}(𝐋^2<𝐋^2>_n),\rho =r\sqrt{\frac{M\omega }{\mathrm{}}},$$
(1)
where

$$<𝐋^2>_n=\frac{n(n+3)}{2}.$$
(2)
The effect of the $`𝐋^2`$ term is to flatten the bottom of the potential. In addition it causes an overall compression of the spectrum, which is avoided through the addition of the $`<𝐋^2>_n`$ term, which will be further discussed later.

The energy eigenvalues of Nilsson’s modified harmonic oscillator are

$$E_{nl}=\mathrm{}\omega \left(n+\frac{3}{2}\right)\mathrm{}\omega \mu ^{}\left(l(l+1)\frac{1}{2}n(n+3)\right),$$
(3)
where $`n`$ is the number of vibrational quanta and $`l`$ is the eigenvalue of the angular momentum, obtaining the values $`l=n`$, $`n2`$, …, 0 or 1 (depending on $`n`$ being even or odd, respectively).

The number of states having energy $`\left(n+\frac{3}{2}\right)\mathrm{}\omega `$ in the case of even $`n`$ is

$$n_1=\underset{l=0,l=\mathrm{even}}{\overset{n}{}}(2l+1)=\frac{(n+1)(n+2)}{2},$$
(4)
where only the even values of $`l`$ are included in the summation and the sum

$$\underset{x=1}{\overset{n}{}}x=\frac{n(n+1)}{2}$$
(5)
has been used. The same result is obtained for odd $`n`$, in which case only the odd values of $`l`$ are included in the summation. Taking into account the spin of the particles the number of states having energy $`\left(n+\frac{3}{2}\right)\mathrm{}\omega `$ is

$$n_2=(n+1)(n+2).$$
(6)
The sum of the eigenvalues of $`𝐋^2`$ within each shell in the case of even $`n`$ is

$$L_1=\underset{l=0,l=\mathrm{even}}{\overset{n}{}}l(l+1)(2l+1)=\frac{n(n+1)(n+2)(n+3)}{4},$$
(7)
where only the even values of $`l`$ have been included in the summation and, in addition to Eq. (5), the sums

$$\underset{x=1}{\overset{n}{}}x^2=\frac{n(n+1)(2n+1)}{6}$$
(8)
and

$$\underset{x=1}{\overset{n}{}}x^3=\frac{n^2(n+1)^2}{4}$$
(9)
have been taken into account. The same result is obtained for $`n`$ being odd. Taking the spin of the particles into account one has

$$L_2=2L_1=\frac{n(n+1)(n+2)(n+3)}{2}.$$
(10)
Thus the average per particle of the square of the angular momentum within each shell is

$$<𝐋^2>_n=\frac{L_1}{n_1}=\frac{L_2}{n_2}=\frac{n(n+3)}{2},$$
(11)
a result which has already been used in Eqs. (1) and (3), in order to keep the “center of mass” of each shell constant, i.e. to counterbalance the overall compression of the spectrum caused by the $`𝐋^2`$ term alone.

The total number of particles which can be accommodated in the levels of the shells up to the $`n`$-th shell included, is

$$N=\underset{x=0}{\overset{n}{}}(x+1)(x+2)=\frac{(n+1)(n+2)(n+3)}{3},$$
(12)
where the spin of the particles has been taken into account and Eqs. (5), (6), (8) have been used. It should be remembered throughout the present work that $`N`$ stands for the total number of particles, while $`n`$ stands for the number of vibrational quanta.

The contribution of the first term of Eq. (3) to the total energy of the particles up to the $`n`$-th shell included (and taking the spin of the particles into account) is

$$E(n)=\mathrm{}\omega \underset{x=0}{\overset{n}{}}\left(x+\frac{3}{2}\right)(x+1)(x+2)=\mathrm{}\omega \frac{(n+1)(n+2)^2(n+3)}{4},$$
(13)
where Eq. (6) is used for the degeneracy within each shell and Eqs. (5), (8), (9) have been used for performing the summations. For later use we notice that omitting the ground state energy contribution in a similar manner one finds

$$E^{}(n)=\mathrm{}\omega \underset{x=0}{\overset{n}{}}x(x+1)(x+2)=\mathrm{}\omega \frac{n(n+1)(n+2)(n+3)}{4},$$
(14)
while a $`n^2`$ perturbation in the energy would have given an additional term

$$E_2(n)=\mathrm{}\omega \underset{x=0}{\overset{n}{}}x^2(x+1)(x+2)=\mathrm{}\omega \frac{n(n+1)(n+2)(n+3)(4n+1)}{20},$$
(15)
where the sum

$$\underset{x=1}{\overset{n}{}}x^4=\frac{n(n+1)(2n+1)(3n^2+3n1)}{30}$$
(16)
has been used in addition to Eqs. (5), (6), (8), (9).

The contribution of the second term of Eq. (3) to the total energy of the particles up to the $`n`$-th shell included is found by using Eq. (10)

$$E_3(n)=\mathrm{}\omega \mu ^{}\underset{x=0}{\overset{n}{}}\frac{x(x+1)(x+2)(x+3)}{2},$$
(17)
while the contribution of the third term of Eq. (3) will be

$$E_4(n)=+\frac{1}{2}\mathrm{}\omega \mu ^{}\underset{x=0}{\overset{n}{}}x(x+3)(x+1)(x+2),$$
(18)
where use of Eq. (6) has been made. We remark that $`E_3(n)`$ and $`E_4(n)`$ cancel. Thus we conclude that in Nilsson’s MO only the first term of Eq. (3) contributes to the total energy of the particles up to the $`n`$-th shell included.

The average energy per particle (up the $`n`$-th shell included) is then found using Eqs. (12) and (13) to be

$$<E>=\frac{E(n)}{N}=\mathrm{}\omega \left(\frac{3}{4}n+\frac{3}{2}\right),$$
(19)
i.e. the average energy per particle is increasing linearly with the shell number, which is the number of vibrational quanta, as it is expected for a harmonic oscillator, since the angular momentum terms make no contribution, as we have already seen. The same result is obtained from Eq. (14)

$$<E^{}>=\frac{E^{}(n)}{N}=\mathrm{}\omega \frac{3}{4}n,$$
(20)
where the $`3/2`$ term has been omitted already in Eq. (14), while a $`n^2`$ perturbation in energy, as seen from Eqs. (12) and (15), would have given an extra term

$$<E_2(n)>=\frac{E_2(N)}{N}=\mathrm{}\omega \frac{3}{20}n(4n+1),$$
(21)
which naturally contains a $`n^2`$ term.

The lower part of the spectrum of Nilsson’s Modified Oscillator, calculated from Eq. (3), is shown in Fig. 1(a) as a function of the parameter $`\mu ^{}`$, together with the magic numbers appearing at different parameter values. The following observations can be made:

a) The magic numbers at the left end of the figure are the ones of the 3-dimensional isotropic harmonic oscillator.

b) The magic numbers appearing around $`\mu ^{}=0.04`$ are in agreement with the magic numbers appearing in alkali clusters, up to 138 (see Ref.  for more details). The agreement is destroyed beyond this point.

c) Around the parameter value $`\mu ^{}=0.02`$ the magic numbers up to 138 are a mixture of magic numbers of the 3-dimensional isotropic harmonic oscillator and magic numbers appearing around $`\mu =0.04`$ (magic numbers of alkali clusters).

d) Eq. (3) can be rewritten as

$$E_{nl}=\frac{3}{2}\mathrm{}\omega +\mathrm{}\omega n\left(1+\frac{3\mu ^{}}{2}\right)+\mathrm{}\omega n^2\frac{\mu ^{}}{2}\mathrm{}\omega \mu ^{}l(l+1),$$
(22)
which clearly shows that in the case of $`\mu ^{}>0`$ the levels within each oscillator shell (characterized by a given value of $`n`$) are ordered according to the value of $`l`$, with the levels with higher values of $`l`$ lying lower in energy, because of the last term in Eq. (22). This is indeed the case in Fig. 1(a), although the levels have not been labelled by the quantum numbers $`n`$, $`l`$ because of lack of space.

e) The level with $`l=n`$ in particular lies lowest in energy within each shell and in general its energy is decreasing with increasing $`\mu ^{}`$, since in this case Eq. (3) takes the form

$$E_{nn}=\frac{3}{2}\mathrm{}\omega +\mathrm{}\omega n\left(1+\frac{\mu ^{}}{2}\right)\mathrm{}\omega n^2\frac{\mu ^{}}{2},$$
(23)
its derivative with respect to $`\mu ^{}`$ being

$$\frac{dE_{nn}}{d\mu ^{}}=\mathrm{}\omega \frac{1}{2}n(1n),$$
(24)
which is indeed negative for $`n>1`$. Indeed the levels which lie lowest within each oscillator shell in Fig. 1(a) are the levels with $`l=n`$, which also show negative slope with increasing $`\mu ^{}`$.

f) The fact that the $`n^2`$ term in Eq. (23) appears with a negative sign (for $`\mu ^{}>0`$) can cause difficulties if one tries to describe a system with a large number of particles in terms of this oscillator. The derivative of $`E_{nn}`$ with respect to $`n`$

$$\frac{dE_{nn}}{dn}=\mathrm{}\omega \left(1+\frac{\mu ^{}}{2}\mu ^{}n\right)$$
(25)
remains positive for

$$n<\frac{1}{\mu ^{}}+\frac{1}{2}.$$
(26)
Beyond this value of $`n`$ the derivative is negative, meaning that levels with higher values of $`n`$ will lie lower in energy, making it difficult to define a cut-off for the number of shells taken into account. For example, if $`\mu ^{}=0.04`$ (a value which has been found  relevant to the description of metal clusters), the derivative remains positive if $`n<25.5`$. It is then clear that a reasonable truncation of the spectrum is possible only if the number of shells to be taken into account is less than 26 (taking the $`n=0`$ shell into account), otherwise no truncation is possible. This drawback of Nilsson’s MO does not have any consequences in the case of nuclear structure, where the model has been first introduced, because of the small number of particles involved there (for which including the shells up to $`n=8`$, shown in Fig. 1(a), suffices), but it can cause difficulties if one tries to employ this model for the determination of supershells in metal clusters, as we shall see later in Sec. 4.

g) As an extention of f), one sees from Eq. (23) that $`E_{nn}`$ remains positive if $`n<\frac{2}{\mu ^{}}+1`$. Thus, in the case of $`\mu ^{}=0.04`$ one should have $`n<51`$, otherwise energies lower than the ground state energy will occur.

3. The 3-dimensional $`q`$-deformed harmonic oscillator (3-dim q-HO)

The space of the 3-dimensional $`q`$-deformed harmonic oscillator consists of the completely symmetric irreducible representations of the quantum algebra uq(3). In this space a deformed angular momentum algebra, soq(3), can be defined . The Hamiltonian of the 3-dimensional $`q`$-deformed harmonic oscillator is defined so that it satisfies the following requirements:

a) It is an soq(3) scalar, i.e. the energy is simultaneously measurable with the $`q`$-deformed angular momentum related to the algebra soq(3) and its $`z`$-projection.

b) It conserves the number of bosons, in terms of which the quantum algebras uq(3) and soq(3) are realized.

c) In the limit $`q1`$ it is in agreement with the Hamiltonian of the usual 3-dimensional harmonic oscillator.

It has been proved  that the Hamiltonian of the 3-dimensional $`q`$-deformed harmonic oscillator satisfying the above requirements takes the form

$$H_q=\mathrm{}\omega _0\left\{[𝒩]_qq^{𝒩+1}\frac{q(qq^1)}{[2]_q}C_q^{(2)}\right\},$$
(27)
where $`𝒩`$ is the number operator and $`C_q^{(2)}`$ is the second order Casimir operator of the algebra soq(3), while

$$[x]_q=\frac{q^xq^x}{qq^1}$$
(28)
is the definition of $`q`$-numbers and $`q`$-operators.

The energy eigenvalues of the 3-dimensional $`q`$-deformed harmonic oscillator are then

$$E_q(n,l)=\mathrm{}\omega _0\left\{[n]_qq^{n+1}\frac{q(qq^1)}{[2]_q}[l]_q[l+1]_q\right\},$$
(29)
where $`n`$ is the number of vibrational quanta and $`l`$ is the eigenvalue of the angular momentum, obtaining the values $`l=n,n2,\mathrm{},0`$ or 1.

In the limit of $`q1`$ one obtains $`\mathrm{lim}_{q1}E_q(n,l)=\mathrm{}\omega _0n`$, which coincides with the classical result.

For small values of the deformation parameter $`\tau `$ (where $`q=e^\tau `$) one can expand Eq. (29) in powers of $`\tau `$ obtaining

$$E_q(n,l)=\mathrm{}\omega _0n\mathrm{}\omega _0\tau \left(l(l+1)n(n+1)\right)$$

$$\mathrm{}\omega _0\tau ^2\left(l(l+1)\frac{1}{3}n(n+1)(2n+1)\right)+𝒪(\tau ^3).$$
(30)
The number of states characterized by a given value of $`n`$, i.e. the number of states in the $`n`$-th shell, is still given by Eq. (4) if the spin is not taken into account, and by Eq. (6) with spin taken into account.

The total number of particles which can be accommodated in the levels of all shells up to the $`n`$-th shell included is still given by Eq. (12), with spin taken into account.

The analogue of the sum of the eigenvalues of $`𝐋^2`$ within each shell in the case of even $`l`$ is, in analogy with Eq. (7), given by

$$L_{1,q}(n)=\underset{l=0,l=\mathrm{even}}{\overset{n}{}}[l]_q[l+1]_q(2l+1)=\frac{1}{(qq^1)^2(q^2q^2)^2}$$

$$((2n+1)(q^{2n+5}+q^{2n5})(2n+5)(q^{2n+1}+q^{2n1})+5(q+q^1)(q^5+q^5)$$

$$\frac{n(n+3)}{2}(q^5+q^5+q^3+q^32q2q^1)),$$
(31)
where $`[l]_q[l+1]_q`$ are the eigenvalues of the second order Casimir operator of soq(3), and use of $`(2l+1)`$ has been made for the degeneracy within a shell without taking spin into account. In performing the relevant summations one needs, in addition to Eq. (5), the sums

$$\underset{x=0}{\overset{n}{}}e^{\tau x}=\frac{e^{\tau (n+1)}1}{e^\tau 1},$$
(32)
$$\underset{x=0}{\overset{n}{}}xe^{\tau x}=\frac{1}{(e^\tau 1)^2}\left(ne^{\tau (n+2)}(n+1)e^{\tau (n+1)}+e^\tau \right),$$
(33)
of which the first is a simple geometric series, while the second can be derived from the first through differentiation with respect to the parameter $`\tau `$. One can easily see that for odd $`l`$ a result identical to the one given in Eq. (31) is obtained.

Using the definition of the $`q`$-numbers given in Eq. (28) the above result can be rewritten in the form

$$L_{1,q}(n)=\frac{q^2}{(q^21)^2}$$

$$\left((2n+1)\left[n+\frac{3}{2}\right]_{q^2}4\left[\frac{n+1}{2}\right]_{q^2}\left[\frac{n}{2}\right]_{q^2}\left(\frac{n}{2}+1\right)(n+1)\left(\left[\frac{3}{2}\right]_{q^2}+\left[\frac{1}{2}\right]_{q^2}\right)+\left[\frac{1}{2}\right]_{q^2}\right),$$
(34)
or equivalently

$$L_{1,q}(n)=\frac{1}{(qq^1)^2}\left((2n+1)\frac{[2n+3]_q}{[2]_q}4\frac{[n]_q}{[2]_q}\frac{[n+1]_q}{[2]_q}\left(\frac{n}{2}+1\right)(n+1)[2]_q+\frac{1}{[2]_q}\right),$$
(35)
where, for example, use of the identity

$$[n]_{q^2}=\frac{q^{2n}q^{2n}}{q^2q^2}=\frac{q^{2n}q^{2n}}{qq^1}\frac{qq^1}{q^2q^2}=\frac{[2n]_q}{[2]_q}$$
(36)
has been repeatedly made.

Eq. (31) can be rewritten in yet another form by using the definition for $`Q`$-numbers (see, for example, the review article in  for relevant details)

$$[n]_Q=\frac{Q^n1}{Q1},$$
(37)
where

$$Q=q^2.$$
(38)
Using this definition Eq. (31) takes the form

$$L_{1,q}(n)=\frac{Q^3}{(Q1)(Q^21)^2}$$

$$\left([n]_QQ^{1/2}\left((Q^25)+2n(Q^21)\right)+[n]_QQ^{1/2}\left((Q^25)+2n(Q^21)\right)\right)$$

$$\frac{Q^{1/2}}{2(Q1)(Q^21)}\left(n^2(Q+1)^2n(Q1)^2\right).$$
(39)
In the limit $`q1`$, keeping terms of order up to $`\tau ^2`$ one can see that Eq. (34) is reduced to Eq. (7), i.e. it is in agreement with the non-deformed case. For this calculation one finds helpful the Taylor expansion of $`q`$-numbers

$$[n]_q=n\pm \frac{\tau ^2}{6}(nn^3)+\frac{\tau ^4}{360}(7n10n^3+3n^5)\pm \frac{\tau ^6}{15120}(31n49n^3+21n^53n^7)+\mathrm{},$$
(40)
where the upper signs correspond to $`q`$ being a phase factor ($`q=e^{i\tau }`$ with $`\tau `$ being real), while the lower signs correspond to $`q`$ being real ($`q=e^\tau `$ with $`\tau `$ being real), as in the present case.

One can now proceed to the calculation of the total energy of the particles up to the $`n`$-th shell included. Using the identity

$$[n]_qq^{n+1}=Q[n]_Q,$$
(41)
where $`q`$-numbers of Eq. (28) ($`Q`$-numbers of Eq. (37) ) are used in the left (right) hand side and $`Q=q^2`$ (Eq. (38) ), one finds that the contribution of the first term of Eq. (29) to the total energy is

$$E_{1,q}(n)=\mathrm{}\omega _0\underset{x=0}{\overset{n}{}}[x]_qq^{x+1}(x+1)(x+2)=\mathrm{}\omega _0\underset{x=0}{\overset{n}{}}Q[x]_Q(x+1)(x+2)$$

$$=\mathrm{}\omega _0\frac{Q}{(Q1)^4}\left((n+1)(n+2)Q^{n+3}2(n+1)(n+3)Q^{n+2}+(n+2)(n+3)Q^{n+1}2\right)$$

$$\mathrm{}\omega _0\frac{Q}{3(Q1)}(n+1)(n+2)(n+3),$$
(42)
where, in addition to Eqs. (5), (6), (8), (32), (33) one also needs to use the sum

$$\underset{x=0}{\overset{n}{}}x^2e^{\tau x}=\frac{1}{(e^\tau 1)^3}\left(n^2e^{\tau (n+3)}(2n^2+2n1)e^{\tau (n+2)}+(n+1)^2e^{\tau (n+1)}e^{2\tau }e^\tau \right),$$
(43)
which is derived from Eq. (33) by differentiation with respect to the parameter $`\tau `$. Using Eq. (37) one can easily see that Eq. (42) can be put in the more symmetric form

$$E_{1,q}(n)=\mathrm{}\omega _0\frac{Q}{(Q1)^3}$$

$$\left((n+1)(n+2)[n+3]_Q2(n+1)[n+2]_Q(n+3)+[n+1]_Q(n+2)(n+3)\right)$$

$$\mathrm{}\omega _0\frac{Q}{3(Q1)}(n+1)(n+2)(n+3).$$
(44)
In the limit of $`Q1`$, keeping terms up to $`T^3`$ (where $`Q=e^T=q^2=e^{2\tau }`$ and thus $`T=2\tau `$) one finds that Eq. (44) agrees with Eq. (13) of the non-deformed case. In this calculation it is helpful to use the Taylor expansion of $`Q`$-numbers

$$[n]_Q=n+\frac{T}{2}(n^2n)+\frac{T^2}{12}(2n^33n^2+1)+\frac{T^3}{24}(n^42n^3+n^2)+\mathrm{}$$
(45)
The contribution of the second term of Eq. (29) to the total energy is found in a similar manner. One has

$$E_{2,q}(n)=\mathrm{}\omega _02\frac{q(qq^1)}{[2]_q}\underset{x=0}{\overset{n}{}}L_{1,q}(x)=\mathrm{}\omega _0\frac{2Q}{(Q^21)^3}$$

$$(Q^4[n]_Q(2(Q^21)n+(Q^22Q7))Q^2[n]_Q(2(Q^21)N+(Q^22Q^17))$$

$$\frac{1}{6}n^2(n+6)(Q+1)(Q^21)^2+\frac{1}{6}n(Q+1)(Q^4+6Q^3+34Q^2+6Q+1)),$$
(46)
where Eqs. (5), (8), (31), (32), (33) have been used and the spin of the particles has been taken into account.

The lower part of the spectrum of the 3-dimensional $`q`$-deformed harmonic oscillator, calculated from Eq. (29), is shown in Fig. 1(b) as a function of the parameter $`\tau `$, together with the magic numbers appearing at different parameter values, while in Fig. 1(c) the full spectrum up to about 1500 particles is exhibited. The following comments and comparisons to Nilsson’s MO are now in place:

a) The magic numbers at the left end of both figures are the ones of the 3-dimensional isotropic harmonic oscillator, as in Fig. 1(a).

b) The magic numbers appearing around $`\tau =0.04`$ are in agreement with the magic numbers appearing in alkali clusters, up to 1500 (see Ref.  for more details), which is the expected limit of validity for theories based on the filling of electronic shells , while in the case of Nilsson’s MO the agreement is limited to the magic numbers up to 138 (Fig. 1(a) ).

c) Around the parameter value $`\tau =0.02`$ the magic numbers up to 138 in Fig. 1(b) are a mixture of magic numbers of the 3-dimensional isotropic harmonic oscillator and magic numbers appearing around $`\tau =0.04`$ (magic numbers of alkali clusters). Beyond 138 other magic numbers appear.

d) Numerical calculations show that in the case of $`\tau >0`$ the levels within each oscillator shell (characterized by a given value of $`n`$) are ordered according to the value of $`l`$, with the levels with higher values of $`l`$ lying lower in energy, because of the last term in Eq. (29). This is indeed the case in Figs. 1(b) and 1(c), although the levels have not been labelled by the quantum numbers $`n`$, $`l`$ because of lack of space. The ordering of the levels within each shell is the same as the one appearing in the case of Nilsson’s MO (Fig. 1(a) ).

e) As in Nilsson’s MO, the level with $`l=n`$ lies lowest in energy within each shell. However, in the present case the energy of this level is not decreasing with increasing $`\tau `$. This is true for all levels of the 3-dim $`q`$-HO: Their energies increase with increasing $`\tau `$ (except in the cases of the levels with $`(n,l)=(0,0)`$ and (1,1), the energies of which remain constant with increasing $`\tau `$).

f) As a consequence of e), no difficulties related to truncation appear in the present case. Stopping the level scheme at the $`l=n`$ level of a given shell and taking into account all levels with lower $`n`$ (i.e. all levels of the shells lying below the given one), one makes sure that all levels up to the given level have been included. Therefore in the 3-dim $`q`$-HO reliable truncations can be made, allowing for the description of systems with many particles. This point is one of the main advantages of the 3-dim $`q`$-HO in comparison to Nilsson’s MO.

A few more comments on the comparison between the 3-dim $`q`$-HO and Nilsson’s MO are also in place:

a) The $`<𝐋^2>_n`$ term in Nilsson’s Hamiltonian (Eq. (1)) has been put in “by hand” in order to guarantee that the “center of mass” of each shell will remain constant, so that shells will not be compressed because of the presence of the $`𝐋^2`$ term. In the case of the 3-dim $`q`$-HO it is clear that the opposite effect is present: The shells are expanded, because of the extra terms added by the $`q`$-deformation. One way to see this is by comparing the first order corrections appearing in Eq. (30) to the last two terms in Eq. (3). In both cases the $`l(l+1)`$ term causes compression of the shells, while expansion of the shells is caused by the $`n(n+1)`$ term in the case of the 3-dim $`q`$-HO and by the $`n(n+3)/2`$ term in the case of Nilsson’s MO. Since the difference of these terms is

$$n(n+1)\frac{1}{2}n(n+3)=\frac{1}{2}n(n1),$$
(47)
which is positive for $`n>1`$, it is clear that the expansion in the case of the 3-dim $`q`$-HO will be stronger than the expansion in the case of Nilsson’s MO, which is exactly balanced by the compression caused by the $`l(l+1)`$ term. As a result, in the case of the 3-dim $`q`$-HO net expansion of the shells will occur, which will be of first order in the parameter $`\tau `$. Comparison of Figs. 1(a) and 1(b) clearly shows this effect.

b) In Nilsson’s MO the last two terms in the energy expression of Eq. (3) give no net contribution to the total energy up to the $`n`$-th shell included, and therefore to the average energy per particle as well, which turns out to be proportional to $`n`$ (see Eq. (19) ). In the 3-dim $`q`$-HO this dependence of the average energy per particle on $`n`$ is given by the lowest order contribution from Eq. (44), while the next order contribution from Eq. (44), as well as the lowest order contribution from Eq. (46), give terms with higher powers of $`n`$. This property will have to be taken into account when supershells will be considered.

c) In both models we have derived exact expressions for the total energy up to the $`n`$-th shell included, i.e. for systems of particles filling complete shells. In order to consider systems of particles for which the last shell is not full, one has to consider numerical methods, as the Strutinsky method , which are beyond the scope of the present study.

4. Supershells in Nilsson’s Modified Oscillator

For studying the existence and properties of supershells in Nilsson’s MO we are going to use the procedure employed by Nishioka et al. . For a given number of particles $`N`$ the single particle energies $`E_j(n,l)`$ of the $`N`$ occupied states are summed up

$$E(N)=\underset{j=1}{\overset{N}{}}E_j(n,l).$$
(48)
This sum is then divided into two parts: A smooth average part $`E_{av}`$ and a shell part $`E_{shell}`$, which will exhibit the supeshell structure

$$E(N)=E_{av}(N)+E_{shell}(N).$$
(49)
For the average part of the total energy a Liquid Drop Model expansion is used

$$E_{av}(N)=a_1N^{1/3}+a_2N^{2/3}+a_3N.$$
(50)
This expansion should be adequate in the case of Nilsson’s MO, for which the average energy per particle increases linearly with $`N`$ (see Eq. (19) ), but it should not suffice in the case of the 3-dim $`q`$-HO, for which the average energy per particle contains higher order terms, as we have seen in the previous section. For the latter case an expansion going up to a $`N^2`$ term

$$E_{av}(N)=a_1N^{1/3}+a_2N^{2/3}+a_3N+a_4N^{4/3}+a_5N^{5/3}+a_6N^2$$
(51)
should be more appropriate. In order to keep the calculations uniform and thus facilitate the comparisons between the two oscillators, we opted for using the expansion of Eq. (51) in all cases, although for Nilsson’s MO the first three terms would have been adequate. Therefore in the case of Nilsson’s MO very small values will be expected for $`a_6`$, the coefficient of the $`N^2`$ term.

The results for $`E_{shell}`$ obtained with Nilsson’s MO for 8 different values of the parameter $`\mu ^{}`$ (assuming $`\mathrm{}\omega =1`$) are shown in Fig. 2, while the relevant parameter values and rms deviations are exhibited in Table 1. As far as the way of calculation is concerned, the following comments apply:

a) In all cases the maximum shell used was $`n_{max}=30`$, implying that the last level included in the calculation was the one with $`(n,l)=(30,30)`$, in order to ensure that the complete spectrum up to this point has been taken into account, as discussed in Secs. 2 and 3.

b) From Eq. (26) one sees that for $`\mu ^{}=0.02`$ one should keep $`n<50.5`$, while for lower values of $`\mu ^{}`$ the limiting value of $`n`$ lies even higher. This means that including the shells up to $`n=30`$ is a reasonable truncation.

c) The procedure of the calculation was as follows: First the summations described by Eq. (48), resulting in the total energy $`E(N)`$ for each particle number $`N`$, have been performed. Subsequently, in order to reduce the size of the calculation approximately by a factor of 10, the average $`E(N)`$ was calculated every 11 points (i.e. for $`N=6`$, 17, 28, …), up to the point immediately below the cut-off of $`(n,l)=(30,30)`$, which is reported in Table 1 as $`N_{max}`$. These averaged values of $`E(N)`$ were subsequently fitted by the expansion of Eq. (51), resulting in the determination of $`E_{av}(N)`$ at these points. Finally $`E_{shell}(N)`$ has been obtained at these points as the difference $`E(N)E_{av}(N)`$ and plotted in Fig. 2.

On the contents of Table 1 the following comments can be made:

a) The behavior of the parameters as functions of $`\mu ^{}`$ is rather smooth, with the exception of $`a_6`$, the coefficient of $`N^2`$, which assumes very small values, as expected from the comments following Eq. (51).

b) The maximum number $`N_{max}`$ of particles below the cut-off is decreasing with increasing $`\mu ^{}`$, as expected from the fact that the level $`(n,l)=(30,30)`$ is getting lower with increasing $`\mu ^{}`$ (see Eq. (24) ).

c) The rms deviations are very small, given the fact that the relevant average energies range up to $`10^6`$.

On the contents of Fig. 2 the following comments apply:

a) The gradual development of supershells with increasing $`\mu ^{}`$ is clearly seen in Figs. 2(a)-2(d). While in Fig. 2(a) ($`\mu ^{}=0.0001`$) no supershell structure is seen up to $`N=10000`$, in Fig. 2(b) ($`\mu ^{}=0.001`$) the first supershell is seen to be completed around $`N=10000`$, in Fig. 2(c) ($`\mu ^{}=0.002`$) two supershells are completed up to $`N=10000`$, and in Fig. 2(d) ($`\mu ^{}=0.003`$) most of the third supershell is also completed by $`N=10000`$.

b) In Figs. 2(a)-2(g) ($`\mu ^{}=0.0001`$ \- 0.01) it is clear that $`N_{s1}`$, the value of $`N`$ at which the first supershell is completed, is decreasing as a function of $`\mu ^{}`$. The same is true for $`N_{s2}`$, the value of $`N`$ at which the second supershell is completed. These trends are not followed by Fig. 2(h) ($`\mu ^{}=0.02`$). An explanation of this effect can be gotten from Fig. 1(a), in which it is clear that for relatively small $`\mu ^{}`$ ($`\mu ^{}<0.01`$, for example) the order of the levels remains the same and only their individual energies change, while at $`\mu ^{}=0.02`$ and beyond the order of the levels changes drastically, especially at higher energies. This drastic mixing of the levels at $`\mu ^{}=0.02`$ and beyond can also explain the increasing difficulty in determining the closing of supershells for large $`\mu ^{}`$ especially at high energies, Fig. 2(h) ($`\mu ^{}=0.02`$) being the clearest example for this case.

c) According to the study by Balian and Bloch of electrons moving in a spherical cavity , which by analogy can be applied to the valence electrons in a metal cluster , the minima of the shell energy, $`E_{shell}`$, which correspond to shell closures, should appear at equidistant positions (i.e. they should exhibit a periodicity) as a function of $`N^{1/3}`$ within each supershell. From the data used for plotting Fig. 2 one can see that this condition is approximately fulfilled.

5. Supershells in the 3-dimensional $`q`$-deformed harmonic oscillator

The procedure described in the previous section has been used in exactly the same way for the determination of supershells in the case of the 3-dim $`q`$-HO.

The results for $`E_{shell}`$ obtained with the 3-dim $`q`$-HO for 10 different values of the parameter $`\tau `$ (assuming $`\mathrm{}\omega =1`$) are shown in Fig. 3, while the relevant parameter values and rms deviations are exhibited in Table 2. In order to facilitate comparisons between the two models, the first 8 values of $`\tau `$ used here are the same as the values of $`\mu ^{}`$ used in the previous section.

As far as the way of calculation is concerned, the maximum shell used in the first 8 cases was $`n_{max}=30`$, as in the previous section, in order to facilitate comparisons between the two models. Lower $`n_{max}`$ has been used only in the last two cases, as shown in Table 2, in order to keep the rms deviation small.

On the contents of Table 2 the following comments can be made:

a) The behavior of the parameters as functions of $`\tau `$ is smooth, even in the case of $`a_6`$, the coefficient of $`N^2`$, which in the case of Nilsson’s MO was not showing smooth behavior.

b) The maximum number $`N_{max}`$ of particles below the cut-off is decreasing with increasing $`\tau `$. This can be explained by looking at Fig. 1(c). The level used as the cut-off is the level with $`(n,l)=(n_{max},n_{max})`$, which, as explained in Sec. 3, is the lowest level within the $`n_{max}`$-th shell and, as seen in Fig. 1(c), increases very slowly with increasing $`\tau `$. It is then clear that the higher the value of $`\tau `$, the more will be the levels coming from below and crossing over the cut-off level, resulting in the decrease of $`N_{max}`$ with increasing $`\tau `$ seen in Table 2.

c) As in the case of Nilsson’s MO, the rms deviations are very small, given the fact that the relevant average energies range up to $`10^6`$.

On the contents of Fig. 3 the following comments apply:

a) The gradual development of supershells with increasing $`\tau `$ is clearly seen in Figs. 3(a)-3(d) ($`\tau =0.0001`$ \- 0.003), which look very similar to Figs. 2(a)-2(d) ($`\mu ^{}=0.0001`$ \- 0.003) of the previous section.

b) In Figs. 2(a)-2(g) ($`\tau =0.0001`$ \- 0.01) it is clear that $`N_{s1,q}`$, the value of $`N`$ at which the first supershell is completed, is decreasing as a function of $`\tau `$, almost in the same way as $`N_{s1}`$ is decreasing as a function of $`\mu ^{}`$ in the case of Nilsson’s MO, the only difference being that for numerically equal values of $`\tau `$ and $`\mu ^{}`$ the corresponding value of $`N_{s1,q}`$ is slightly higher that the relevant value of $`N_{1s}`$, a fact that can be explained by the general property of the spectrum of the 3-dim $`q`$-HO to expand more rapidly than the spectrum of Nilsson’s MO, as seen in Sec. 3. The same is true for $`N_{s2,q}`$, $`N_{s3,q}`$, $`N_{s4,q}`$, i.e. the values of $`N`$ at which the second, third, and fourth supershells are completed. These trends are not followed by Figs. 3(h) ($`\tau =0.02`$), 3(i) ($`\tau =0.038`$), 3(j) ($`\tau =0.05`$), which correspond to higher values of $`\tau `$, where the mixing of the levels is very strong in comparison to the picture existing at low values of $`\tau `$ (i.e. for $`\tau <0.01`$, as seen in Fig. 1(c) ).

c) Despite the fact that the systematics of supershell closures are modified beyond $`\tau =0.02`$ because of the strong mixing of levels, supershells are still seen beyond this point, while this was not possible in the case of Nilsson’s MO, because of truncation related problems, as we have seen in Sec. 3. It should be noticed that this is exactly the region of $`\tau `$ values which has been found relevant for the description of metal clusters .

d) In the case of alkali clusters, Na in particular, the first supershell is known to occur around $`N=1000`$ , which is quite in agreement with what is seen in Fig. 3(i), which corresponds to $`\tau =0.038`$, the parameter value for which the magic numbers of alkali clusters are correctly reproduced  up to 1500 particles, which is the expected limit of validity of theories based on the filling of electronic shells . It should be noticed that only one free parameter exists in the theory, namely $`\tau `$, which has been fixed in Ref.  in order to reproduce the magic numbers of alkali clusters. Therefore no free parameter has been left over in the calculation of supershells. The fact that the present calculation leads to a reasonable prediction of the position of the first supershell closure is a quite stringent test of the present theory.

e) As we have already mentioned, the theory by Balian and Bloch  for electrons moving in a spherical cavity, which by analogy can be applied to the valence electrons of metal clusters , predicts that the minima of the shell energy, $`E_{shell}`$, which correspond to shell closures, should appear within each supershell at equidistant positions (i.e. they should exhibit a periodicity) when plotted versus $`N^{1/3}`$. From the data used for plotting Fig. 3 one can see that this condition is approximately fulfilled.

f) A further prediction of the theory by Balian and Bloch  is that the plot of the cubic root $`N_i^{1/3}`$ of the magic numbers $`N_i`$, corresponding to shell closures, versus the index $`i`$, counting the shells, should be a straight line . In the case of alkali clusters, in particular, in which the triangular and squared closed orbits are supposed to dominate , the slope of the line should be 0.603 . In order to make a preliminary estimate of the degree to which this prediction is fulfilled, one can employ the data used in plotting Fig. 3(i), which corresponds to $`\tau =0.038`$, i.e. to the parameter value found appropriate  for reproducing the magic numbers of several alkali clusters. One can then see that this requirement is roughly fulfilled, although the slope appears to be gradually decreasing at large $`i`$, an effect which might be due to missing corrections mentioned in comment h).

g) In the case of Al clusters, magic numbers have been determined experimentally in detail up to 1200 electrons , while additional experimental results up to 2700 electrons exist . The slope of $`N_i^{1/3}`$ vs. $`i`$ in this case is considerably lower (around 0.32) , indicating that closed orbits other than the triangular and squared ones should be present in the Balian and Bloch approach . It will be interesting to examine if the 3-dim $`q`$-HO can provide any prediction for supershells in Al clusters, after choosing the value of the $`\tau `$ parameter in order to reproduce the right slope of $`N_i^{1/3}`$ vs. $`i`$.

h) In order to guarantee the reliability of such predictions, it is of interest to study in advance the influence of any modifications imposed by the Strutinsky method , which can be applied in cases in which the last shell is open, which are beyond the realm of the present analytic study of Section 3. The influence of the well known quantum mechanical effect that $`\mathrm{}\omega _0`$ should be decreasing with increasing number of particles in the cluster  should also be taken into account.

6. Discussion

The main results of the present study are in summary the following:

a) The 3-dimensional $`q`$-deformed harmonic oscillator, which is known  to describe very well the magic numbers of alkali clusters up to 1500 particles (the expected limit of validity for theories based on the filling of electronic shells), is found in the present study able to produce supershell structures, succesfully predicting the first supershell in alkali clusters without involving any free parameter in addition to the deformation parameter $`\tau `$, which has been fixed for reproducing the magic numbers.

b) It should be noticed that these successes of the 3-dim $`q`$-HO are largely due to terms in the Hamiltonian induced by the symmetry, which make the spectrum of the 3-dim $`q`$-HO to expand with increasing shell number $`n`$ more rapidly than the corresponding spectrum of Nilsson’s Modified Oscillator, allowing, among other things, for reliable truncations to be performed.

c) The successful prediction of the magic numbers can be considered as evidence that the 3-dimensional $`q`$-deformed harmonic oscillator owns a symmetry \[the uq(3) $``$ soq(3) symmetry, which is a nonlinear deformation of the u(3) symmetry of the spherical (3-dimensional isotropic) harmonic oscillator\] appropriate for the description of the physical systems under study. The use of this symmetry for predicting supershells in other kinds of metal clusters (Al clusters, for example) is an interesting open problem. The influence of using open shells (taken into account by the Strutinsky method ), as well as of the well known quantum mechanical fact that $`\mathrm{}\omega _0`$ should be decreasing with increasing number of particles in the cluster , should be taken into account before any final predictions can be made.

Acknowledgements

Support from the Bulgarian Ministry of Science and Education under Contract Nos. $`\mathrm{\Phi }`$-415 and $`\mathrm{\Phi }`$-547 is gratefully acknowledged.

Table 1

Parameters used for fitting the average part of the total energy (see Eq. (51)) in the case of Nilsson’s Modified Oscillator, for various values of the model parameter $`\mu ^{}`$ (see Eq. (3)), corresponding to the cases exhibited in Fig. 2. The parameters are dimensionless, since we have assumed $`\mathrm{}\omega =1`$ (see Eq. (3)) throughout. The highest shell $`n_{max}`$ and the number of particles $`N_{max}`$ included in each calculation, as well as the relevant rms deviation $`\sigma `$, are also shown. See Section 4 for further discussion.

| $`\mu ^{}`$ | $`a_1`$ | $`a_2`$ | $`a_3`$ | $`a_4`$ | $`10^3a_5`$ | $`10^5a_6`$ | $`n_{max}`$ | $`N_{Max}`$ | $`\sigma `$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0.0001 | 13.349 | -6.556 | -0.218 | 0.972 | 4.451 | 6.853 | 30 | 10027 | 22.06 |
| 0.001 | -0.025 | 0.281 | -1.482 | 1.080 | 0.080 | -0.132 | 30 | 10027 | 7.12 |
| 0.002 | -0.022 | 0.181 | -1.439 | 1.073 | 0.489 | 1.053 | 30 | 9609 | 4.19 |
| 0.003 | -1.281 | 1.080 | -1.651 | 1.094 | -0.445 | 0.282 | 30 | 9246 | 3.40 |
| 0.005 | -1.073 | 1.013 | -1.659 | 1.096 | -0.623 | 0.043 | 30 | 8641 | 2.76 |
| 0.007 | -0.830 | 0.910 | -1.654 | 1.097 | -0.724 | -0.626 | 30 | 7937 | 4.16 |
| 0.01 | -5.448 | 3.432 | -2.155 | 1.143 | -2.645 | 0.383 | 30 | 7189 | 4.88 |
| 0.02 | 5.692 | -3.891 | -0.376 | 0.938 | 9.135 | -38.541 | 30 | 4890 | 3.50 |

Table 2

Parameters used for fitting the average part of the total energy (see Eq. (51)) in the case of the 3-dimensional $`q`$-deformed harmonic oscillator for various values of the model parameter $`\tau `$ (see Eq. (29), with $`q=e^\tau `$), corresponding to the cases exhibited in Fig. 3. The parameters are dimensionless, since we have assumed $`\mathrm{}\omega _0=1`$ (see Eq. (29)) throughout. The highest shell $`n_{max}`$ and the number of particles $`N_{max}`$ included in each calculation, as well as the relevant rms deviation $`\sigma `$, are also shown. See Section 5 for further discussion.

| $`\tau `$ | $`a_1`$ | $`a_2`$ | $`a_3`$ | $`a_4`$ | $`10^3a_5`$ | $`10^5a_6`$ | $`n_{max}`$ | $`N_{Max}`$ | $`\sigma `$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0.0001 | 13.585 | -6.554 | -0.228 | 0.972 | 4.464 | -6.776 | 30 | 10027 | 22.12 |
| 0.001 | 0.753 | -0.061 | -1.424 | 1.073 | 0.867 | -0.272 | 30 | 10027 | 7.58 |
| 0.002 | 1.458 | -0.496 | -1.318 | 1.059 | 2.075 | -1.136 | 30 | 9686 | 4.68 |
| 0.003 | -3.358 | 2.186 | -1.855 | 1.106 | 0.615 | 2.501 | 30 | 9389 | 3.52 |
| 0.005 | 0.667 | 0.166 | -1.500 | 1.074 | 2.791 | 2.422 | 30 | 8795 | 3.49 |
| 0.007 | -2.631 | 1.965 | -1.867 | 1.106 | 2.204 | 7.580 | 30 | 8421 | 3.88 |
| 0.01 | -3.657 | 2.568 | -2.013 | 1.120 | 2.523 | 15.399 | 30 | 7893 | 6.75 |
| 0.02 | -17.409 | 11.661 | -4.267 | 1.379 | -10.052 | 88.547 | 30 | 7189 | 9.03 |
| 0.038 | -55.417 | 41.650 | -13.067 | 2.615 | -96.648 | 482.365 | 26 | 4648 | 17.98 |
| 0.05 | -64.123 | 55.364 | -19.160 | 3.777 | -202.994 | 999.866 | 22 | 3020 | 16.63 |

Figure captions

Fig. 1 (a) Energy spectrum of Nilsson’s Modified Oscillator (in units of $`\mathrm{}\omega `$, see Eq. (3)) as a function of the (dimensionless) model parameter $`\mu ^{}`$. Magic numbers are shown at the main gaps. (b) Energy spectrum of the 3-dimensional $`q`$-deformed harmonic oscillator (in units of $`\mathrm{}\omega _0`$, see Eq. (29)) as a function of the (dimensionless) model parameter $`\tau `$ (with $`q=e^\tau `$, where $`\tau `$ is real). (c) Same as (b), but extended to higher energy levels.

Fig. 2 Shell part ($`E_{shell}`$) of the total energy (in units of $`\mathrm{}\omega `$, see Eq. (3)) for Nilsson’s Modified Oscillator vs. the number of particles $`N`$. The values of the (dimensionless) parameter $`\mu ^{}`$ are the same as these listed in Table 1, together with the details of the calculation. See Section 4 for further discussion.

Fig. 3 Shell part ($`E_{shell}`$) of the total energy (in units of $`\mathrm{}\omega _0`$, see Eq. (29)) for the 3-dimensional $`q`$-deformed harmonic oscillator vs. the number of particles $`N`$. The values of the (dimensionless) parameter $`\tau `$ are the same as these listed in Table 2, together with the details of the calculation. See Section 5 for further discussion.


format: text

text:
# Possible evidence for the pulsed X-ray emission from outer gap in PSR B1937+21?

## 1 Introduction

Even more than 30 years after the discovering of high energy pulsars, the theoretical reproduction of X-ray and $`\gamma `$-ray emission from such pulsars is still a matter of debate. It is commonly agreed that there are two scenarios on modelling X-ray and $`\gamma `$-ray creation: the outer gap (OG) model (e.g., Cheng, Ho & Ruderman 1986a,b, Romani 1996, Cheng & Zhang 1999, Hirotani & Shibata 2001) and the polar cap (PC) model (e.g., Harding 1981, Sturner & Dermer 1994, Daugherty & Harding 1996, Luo et al. 2000, Zhang & Harding 2000). The fundamental difference between these two types of models is the location of the regions where particles are accelerated to relativistic energies and emit high energy photons. The early PC model (Harding 1981) assumed that the emission is produced just above the PC surface. In the present versions of PC model (e.g. Daugherty & Harding 1996, Harding & Muslimov 1998), it is proposed that the particle acceleration region may extend from the PC surface to several stellar radius above due to free charge flow and inertial frame dragging, so that wide double-peak X-ray and $`\gamma `$-ray light curves can be reproduced, given the inclination angle is not large. In contrast, the OG model (Cheng et al. 1986a, CHR) presumed that the gaps can exist in the outer magnetosphere between the null charge surface and the light cylinder. Later, CHR model was developed to the single OG models (e.g., Chiang & Romani 1994, Romani & Yadigaroglu 1995, Cheng et al. 2000), which claimed that the three dimension extents of OGs are constrained by the pair cascade processes, and a single OG can produce wide double-peak high energy light curves. Whereas theoretical considerations of more detailed physical processes for particle acceleration and photon emission are necessary, it is urgent and interesting to find new observational evidence for these models.

PSR B1937+21, with the period of 1.56 ms, is the fastest millisecond pulsar (MSP) known. At radio bands, it exhibits an interpulse emission that is roughly equal to the main pulse in intensity and separated from it by a phase of about $`180^\mathrm{o}`$. Recently, the ASCA observations detected pulsed X-ray emission from this pulsar (Takahashi et al. 2001). Only one nonthermal narrow pulse was observed, which is coincident with the radio interpulse in phase within the timing errors. The pulse width is about $`100\pm 61`$ $`\mu `$s ($`23^\mathrm{o}\pm 14^\mathrm{o}`$). Besides the narrow pulse, the light curve reveals two additional wide Gaussian-shaped bulges above the background level, with each phase interval being about 0.5 rotation period and each peak intensity $`1/4`$ of the nonthermal pulse peak.

Where does the X-ray emission of the fastest MSP come from, the PC or the OG? It has not been extensively studied as the other high energy pulsars such as the Crab pulsar and Vela. The previous work was done by Luo et al (2000), in which the PC model was modified for MSPs. They applied the theory to PSR B1937+21 and suggested that the X-ray emission probably originates from the location of one stellar radius above the PC.

In this paper, efforts on modelling the observational data at radio and X-ray bands are made for both the PC and OG models. The inclination angle between the rotation and magnetic axes is constrained in section 2, which is necessary for the modelling. The modelling is carried out in section 3. Conclusions and discussions are placed in section 4.

## 2 The inclination angle

The inclination angle ($`\alpha `$) between the rotation and magnetic axes is a necessary parameter for both the PC and OG models to give various high energy emission beams. Unfortunately, there is no agreement on the value of inclination angle of PSR B1937+21. In this section we reinvestigate the value of $`\alpha `$ under the double-pole model, viz., the radio interpulse and main pulse are considered to be from the opposite magnetic poles of a dipolar field, and the result is used in the calculation in section 3.

We assume that the radio emission beams from double poles are axisymmetric around the magnetic axis and have the same radius, namely, $`\rho _1=\rho _2`$. According to the geometry model (Gil et al. 1984, Lyne & Manchester 1988, LM88) one has

$$\mathrm{sin}^2\frac{\rho _1}{2}=\mathrm{sin}^2\frac{W_1}{4}\mathrm{sin}\alpha _1\mathrm{sin}(\alpha _1+\beta _1)+\mathrm{sin}^2\frac{\beta _1}{2}$$
(1)
$$\mathrm{sin}^2\frac{\rho _2}{2}=\mathrm{sin}^2\frac{W_2}{4}\mathrm{sin}\alpha _2\mathrm{sin}(\alpha _2+\beta _2)+\mathrm{sin}^2\frac{\beta _2}{2},$$
(2)
where $`\beta `$ is impact angle between the line of sight (LOS) and the magnetic axis, $`W`$ is the pulse width of the average profile, the subscripts ‘1’ and ‘2’ denote the main pulse and the interpulse, respectively. There are two simple geometrical relations between the inclination angles and the impact angles, namely, $`\alpha _2=\pi \alpha _1`$, and $`\beta _2=\alpha _1+\beta _1\alpha _2`$.

¿From the above relations $`\beta _1`$ can be derived as

$$\beta _1=\mathrm{tan}^1\left[\frac{A\mathrm{tan}^2\alpha _1}{(A+1)\mathrm{tan}\alpha _1}\right],$$
(3)
where

$$A=\frac{1}{\mathrm{sin}^2\frac{W_1}{4}\mathrm{sin}^2\frac{W_2}{4}}.$$
(4)
In the following we neglect the subscript ‘1’, so that all the $`\alpha `$, $`\beta `$ and $`\rho `$ below are referred to the main pulse except for special declaration.

To figure out the value of $`A`$, the profile at 1.5 GHz is chosen (downloaded from EPN) for its high time resolution and low dispersion smearing (Kramer et al. 1998). The pulse widths are measured at the level of 10% of their peak intensities, which are $`19^\mathrm{o}.4\pm 0^\mathrm{o}.4`$ and $`17^\mathrm{o}.8\pm 0^\mathrm{o}.4`$ for the main pulse and the interpulse respectively. So that $`\beta (\alpha )`$ and $`\rho (\alpha )`$ can be calculated, as shown by the solid curve and the curve dotted by circles respectively in Fig.1. The $`\beta (\alpha )`$ curve approximates to a linear function of $`\alpha +\beta 90^\mathrm{o}`$ when $`\alpha 85^\mathrm{o}`$, which is determined by the fact that the main pulse and interpulse are both narrow and differ only a little in pulse width.

For a dipolar field, the shape of the polar cap is found to change from a circle (for $`\alpha =0^\mathrm{o}`$) to an ellipse of which the longitudinal radius is about 1.6 times of the latitudinal radius (for $`\alpha =90^\mathrm{o}`$, Cheng et al. 2000). Since the deviation from circular shape is not essentially significant, we simply regard the polar cap as a circle in this section. So that the opening angle (between the magnetic axis and the tangent of the magnetic field line) of the polar cap $`\rho _{\mathrm{PC}}`$ can be determined by the last open field line on the plane containing the rotation and magnetic axes ($`𝛀\mu `$ plane). The radius $`\rho _{\mathrm{PC}}`$ is a function of the stellar radius R and the inclination angle $`\alpha `$, as shown by the dash curves in Fig.1 for $`R=3`$ km and 10 km.

We further assume that the boundary of the radio beam is defined by the last open field lines, then the beam radius should be greater than $`\rho _{\mathrm{PC}}`$. According to this criteria, the inclination angle is constrained to be $`\alpha 63^\mathrm{o}`$ for $`R=10`$ km and $`\alpha 77^\mathrm{o}`$ for $`R=3`$ km. As the radii of neutron stars (NSs) are currently believed to be about 10 km, we accept that $`\alpha 63^\mathrm{o}`$.

Can we determine the exact value of $`\alpha `$? It is sure if the maximum rate of position angle $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}`$ is exactly known, where $`\psi `$ is the position angle and $`\varphi `$ the pulse longitude. The maximum rate presents the second relationship between $`\alpha `$ and $`\beta `$, which reads

$$\left(\frac{\mathrm{d}\psi }{\mathrm{d}\varphi }\right)_{\mathrm{max}}=\frac{\mathrm{sin}\alpha }{\mathrm{sin}\beta }.$$
(5)
Combining with Eq.s 3 and 5, $`\alpha `$ and $`\beta `$ can be solved for a given $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}`$.

Although recent polarization observations present flat position angle sweeps (Thorsett & Stinebring 1990, Stairs et al. 1999), it can not be asserted that the real value of $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}`$ is small, because observations may give a less steep position angle gradient due to smearing of finite sampling time, to the frequency dispersion in pulse arrival time (Liu & Wu, 1999), and to the interstellar scattering (Gil 1985a). The $`\beta ^{}(\alpha )`$ curves derived from Eq. 5 for $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}=1,3`$ and 20 respectively are presented by the dotted curves in Fig.1. The intersections of $`\beta ^{}(\alpha )`$ and $`\beta (\alpha )`$ shows that a larger $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}`$ results in a larger $`\alpha `$ and smaller $`\beta `$, which means the LOS sweeps across the radio beam more closely to the beam center.

One may find that when the real value of $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}`$ is large enough, $`\alpha `$ would exceed $`63^\mathrm{o}`$ (for example, taken $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}=3`$ as proposed by Gil (1985a), $`\alpha `$ is $`71^\mathrm{o}`$) and hence contradict against $`\alpha 63^\mathrm{o}`$. However, if the radius is smaller, for example, $`R=3`$ km, this inconsistency would cancel. In fact, Xu et al. (2001) suggested that PSR B1937+21 is probably a strange star (SS) with low mass and small radius. The detailed discussion is placed in section 4.

The range of $`\alpha `$ presented above is different from the conventional consideration in double-pole model that $`\alpha `$ should be close to 90o (Stairs et al. 1999). Alternatively, there is another kind of so-called single pole model to interpret the interpulse, which suggests that the interpulse emission comes from the same pole as the main pulse. In the single pole model proposed by Gil (1985a) for PSR B1937+21, $`\alpha `$ only need to be 20o. The single pole model predicts that the separation between the main pulse and interpulse may be close to 180o and is frequency independent (Gil 1983, 1985b), which is coincident with the observation of PSR B1937+21 (Hankins & Fowler 1986). However, observations with high time resolution (e.g. Kramer et al. 1998, Stairs et al. 1999) failed to find the emission components between the main pulse and interpulse, which was reported by Stinebring et al. (1984) and was suggested to be a strong support to the single pole model (Gil 1985a). Therefore, in this paper the radio emission is considered to come from double poles. More confirmative estimations of $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}`$ are expected to determine $`\alpha `$ and $`\beta `$.

## 3 The origin of nonthermal X-ray emission from PSR B1937+21

In this section we calculate the X-ray beams in the frames of both the PC and OG models to find out whether they are able to reproduce the observational facts. The facts used here are:

(1) at 1.4 GHz the separation between the peaks of interpulse and main pulse is $`174^\mathrm{o}`$ (measured from the profile presented by Takahashi et al. (2001));

(2) at 1.5 GHz (EPN data) the 10% widths of the main pulse and interpulse are $`19.4^\mathrm{o}`$ and $`17.7^\mathrm{o}`$, respectively;

(3) the nonthermal X-ray pulse is nearly coincident with the interpulse, the X-ray pulse width is about $`23^\mathrm{o}`$ (Takahashi et al. 2001).

### 3.1 Origin from the PC accelerators?

Luo et al. (2000) discussed the viability of PC models for high energy emission from MSPs. They found that the maximum Lorenz factor of particles is limited by curvature radiation and not sensitive to the specific acceleration model, but the height where the Lorenz factor achieves the maximum is model dependent, which may be between $`0.01R`$ (for the inner vacuum gap) and above $`0.1R`$ (for the space-charge limited gap) from PC surface for pulsar period $`P=2`$ ms and a surface magnetic field $`B_\mathrm{s}=7.5\times 10^8`$ G. Assuming a space-charge limited flow, the pair cascades can occur at the typical distance (to the star center) of $`r(1.52.5)R`$, and high energy emission is radiated from this region. With respect to $`P=1.56`$ ms and $`B_\mathrm{s}=4.1\times 10^8`$ G for PSR B1937+21, their analysis applies to this pulsar.

Since the radio emission is also radiated from the region near the PC, homocentric radio and X-ray beams are produced (Fig.2a). Quantitatively, we assume the radio emission arises from the PC up to the distance $`r=2.5R`$, the X-rays may arise from the PC to a (a) relatively higher distance, e.g., $`r=3R`$, or (b) lower distance, e.g., $`r=2R`$. For simplicity we assume the emission regions are bounded by the last open field lines. In case (a) the X-ray beams are wider than the radio beams, therefore, when the line of sight sweeps across the both radio beams (to reproduce the main pulse and the interpulse) it would sweep across the both X-ray beams either, and hence gives double X-ray pulse, which is inconsistent with the observation. In case (b) the X-ray beams are narrower than the radio beams, and then, for a proper viewing angle (between LOS and the rotation axis), only one X-ray beam could be observed. In this case, could the fact (3) be explained as well? We analyze this issue by using the figure of the (phase-viewing angle) plane on which the emission beams are projected.

In Fig.3, the horizon axis is pulse phase, the vertical axis is is the viewing angle $`\zeta `$ ($`\zeta =\alpha +\beta `$). The dotted curves represent the boundaries of X-ray beams, the solid curves are the boundaries of radio beams, and the PCs are shown by the dashed curves. Retardation (due to distinct emission heights) and abberation effects (due to emission sources co-rotating with the pulsar) have been taken into account, both of which make the lower beam move towards the trailing edge of the higher beam. To obtain the figure, a moderate inclination angle is taken, $`\alpha =60^\mathrm{o}`$. Although there is $`\alpha 63^\mathrm{o}`$ as discussed in section 2, the angle should not be too small, or it would give too wide radio beams (Fig.1), which must be emitted from unbelievably high distances near the light cylinder.

According to Fig.3, when $`\zeta 80^\mathrm{o}`$ the X-ray peak is coincident with the radio pulse centered on phase about 0.9; when $`\zeta 100^\mathrm{o}`$ it is coincident with the radio pulse on phase about 0.4. However in both cases, the radio pulse associated with the X-ray pulse is the main pulse which is wider and more intense than the other one.

However, we should point out that the above analysis is based on a simple assumption that the radio (X-ray) beams from both poles have the same width and the emission pattern in the beams is symmetrical about the magnetic axis. Since the magnetosphere structure of pulsars and the detailed radio emission process are still unknown exactly, possibly significant deviation from the assumption can not be ruled out, and the validity of the PC model need to be considered further. We discuss such possibility in section 4.

### 3.2 Origin from the OG accelerators

In the original OG model CHR suggested that a global current flow through the magnetosphere can result in large regions of OGs between the null charge surface and the light cylinder along the last open field lines. Within the OGs particles with one kind of charge are accelerated outward from the star and give an outward emission beam, while those with the opposite charge are accelerated towards the star and give an inward beam. The high energy photons were proposed to be emitted from two OGs associated with the two poles so that double-peaked $`\gamma `$-ray pulse profile can be reproduced, of which one peak corresponds to the outward beam from one OG and the other peak to the inward beam from the opposite OG. CHR assumed that the OGs are active only near the $`𝛀\mu `$ plane. However, this assumption is merely valid for large inclination angles.

The currently prevalent OG models are the single OG models (e.g. Chiang & Romani 1992, Chiang & Romani 1994, Cheng et al. 2000). Generally the inward emission is not important in these models for the reason that the inward high energy photons can not pass freely through the inner magnetosphere due to magnetic pair production. The outward emission from the OG associated with a single pole can produce a broad, irregularly-shaped emission beam which is particularly dense near the edge. The OG regions can be supported along all the last open field lines, but the three dimension scales of OGs are limited by the pair production processes.

In the latest version of this type of model (Cheng et al. 2000), the fraction size ($`fh/R_\mathrm{L}`$) of the gap is $`f5.5P^{26/21}B_{12}^{4/7}\xi ^{1/7}`$, which can be estimated by the threshold of $`\gamma \gamma `$ pair production, $`E_\mathrm{X}(f)E_\gamma (f)(m_\mathrm{e}c^2)^2`$, where $`h`$ is the mean vertical extension perpendicular to the magnetic field, $`R_\mathrm{L}`$ is the radius of light cylinder, $`\xi =\mathrm{\Delta }\varphi /2\pi `$, $`\mathrm{\Delta }\varphi `$ is the transverse extension of the gap, $`E_\mathrm{x}`$ is the energy of the X-ray photons emitted from hot PCs, and $`E_\gamma `$ is the characteristic photon energy emitted by the relativistic particles. The radial scale of pair production is limited to a range between $`r_{\mathrm{in}}`$ and $`r_{\mathrm{lim}}6r_{\mathrm{in}}(\xi =0)`$, where $`r_{\mathrm{in}}`$ (the subscript ‘in’ means the inner boundary of the OG) is the distance of null charge surface, $`\xi =0`$ corresponds to the $`𝛀\mu `$ plane.

In the following our modelling is in the frame of single OG model. Only the outward emission beams from two OGs are considered, as illustrated by Fig.2b. For PSR B1937+21, we have a thin OG with $`f=0.16`$ $`\xi ^{1/7}`$, so that the X-rays can be simply regarded as being radiated from the last open field lines unless $`\xi `$ is too small. The radial scale $`r_{\mathrm{lim}}/r_{\mathrm{in}}`$ and the transverse scale $`\xi `$ are free parameters in calculating the X-ray beams.

First we consider a general situation of the OG scenario to test if the observation facts can be hopefully reproduced. From Fig.2b one can see that provided the observer’s viewing angle is not just $`90^\mathrm{o}`$, e.g., $`\zeta =83^\mathrm{o}`$, the LOS can sweep across both of the radio beams and only one X-ray beam. By assuming a group of reasonable values of parameters (see the first line of Table 1), namely, the inclination angle, the stellar radius, the distance of radio emission, and the radial and transverse extensions of the OGs, the X-ray beams are calculated, which are demonstrated by the line-shadowed areas in Fig.3. Retardation and abberation effects are also included. It shows clearly that the X-ray pulse could be associated with the radio interpulse (which has a smaller width than the other one), and the narrow X-ray pulse width may be obtained if proper extension of the OG is assumed.

Then we model the observational data. The parameters listed in the second line of Table 1 are found out to be able to reproduce the narrow X-ray pulse width, the radio pulse widths, and the coincidence between the X-ray pulse and the radio interpulse, which are in good agreement with the observational data. In order to simulate the observational profiles, we simply assume that the X-ray and radio pulses are Gaussian shapes, and assume additionally wide, weak, hot X-ray emission from both of the PCs. The theoretical profiles are plotted in Fig.4, together with the observational profiles for comparison.

It should be pointed out that modelling the X-ray pulse width is not sensitive to the value of $`r_{\mathrm{lim}}/r_{\mathrm{in}}`$ but to $`\xi `$, thereby the range of $`r_{\mathrm{lim}}/r_{\mathrm{in}}`$ is relaxed and a reasonable value is chosen. Other groups of values are also tried. It is found that for the moderate inclination angles $`40^\mathrm{o}\alpha 63^\mathrm{o}`$ the observational data can be reproduced with proper choice of the gap size ($`\xi =40^\mathrm{o}`$ in Table 1 is approximately the maximum transverse scale). Therefore, our calculation suggests that the nonthermal X-rays of PSR B1937+21 may be emitted from the OGs.

## 4 Conclusions and Discussions

The discovery of pulsed X-ray emission from the MSP PSR B1937+21 and the phase alignment between the X-ray pulse and the radio interpulse provide valuable information on the emission mechanism. In this paper we investigate the possible origin of the X-rays from both the PC and OG accelerators. In the frame of the prevailing OG model (e.g. Cheng et al. 2000), the X-rays from outer gap accelerators could account for the main observational facts by assuming proper size of the OGs for an oblique rotator. In the frame of PC model, by assuming symmetric geometry for radio and X-ray emission, the X-rays from polar cap accelerators would produce an X-ray peak aligned with the main radio pulse in phase, which is contradictive against the observation.

Some more discussions for both the OG and PC models are presented as follows. First we refer to the OG model.

There is a slight inconsistency as shown in Fig.4, i.e., in the OG model, the calculated separation between the interpulse and main pulse is $`180^\mathrm{o}`$, $`6^\mathrm{o}`$ greater than the observational value. This may be due to the retardation effect. The phase shift for different heights can be roughly estimated as $`\mathrm{\Delta }s=\mathrm{\Delta }r/(Pc)`$. A difference of $`\mathrm{\Delta }r=0.8R`$ is enough to produce the phase shift of $`6^\mathrm{o}`$.

A moderate inclination angle $`55^\mathrm{o}`$ is used in the modelling. But $`\alpha `$ could be larger if the real gradient of position angle is steeper than the present observations as suggested in section 2. Could the OG model be still valid for large value of $`\alpha `$?

As shown in Fig.1, when $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}4`$ it would result in a puzzling problem that the derived radio beam radius is considerably smaller than the PC radius, if $`R=10`$ km and a magnetic dipole are assumed. A much smaller stellar radius could cancel the problem, but this requires an SS scenario, because the smallest radius of NSs is generally believed to be $``$9 km while SSs can have much smaller radius due to their different equations of state. In fact, according to the observational limits on the radius and mass derived from the pulse width and $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}`$, Xu et al. (2001) suggested that PSR B1937+21 is probably a strange star with much low mass, small radius and weak magnetic moment. If PSR B1937+21 is an SS with small radius, for example, $`R=2`$ km, a group of parameters is found out to be able to model the observational data by OG model, which is listed in the third line in Table 1.

Then we turn to the PC model. When the emission geometry is not symmetrical, the PC model could be able to explain the observation as well. The possible asymmetry is discussed for the radio and the X-ray emission, respectively.

(1) It has been assumed that the radio emission pattern from the two poles is the same in our modelling; the different behaviors of the main pulse and the interpulse are thus geometrical origin. However, their difference may be caused intrinsically since the mechanism for the radio emission is still poorly understood. For example, a possible reason may be that the pulsar has a non-dipolar magnetic field, thus the distribution of emission intensity in the radio beams may be different from each other, which could lead to that for one beam only part of it is observed while for the other a larger part or the whole is observed.

(2) When PSR B1937+21 is an NS or an SS with crust and its inclination angle is near $`90^\mathrm{o}`$, the X-ray emission pattern may be significantly asymmetric about the magnetic axis. In this case, the space charge is negative on the side toward the rotation axis (where $`𝛀𝐁>0`$, hereafter side I) and positive on the side away from the axis (where $`𝛀𝐁<0`$, hereafter side II). On side II the ions can be pulled away from the surface by strong electric field if the binding energy is small enough. Therefore, the half beam from side II could be much less luminous than that from side I due to the much smaller Lorentz factors of the ions, and thus may be too weak to be observed. Notice that, if side I is above the equator on one pole, it should be below the equator on the other pole. Therefore, only single X-ray peak is observed which may coincide with the radio interpulse.

In recent years, it is suggested that pulsars may be bare strange stars (BSSs) (e.g., Xu 2002 and references therein). If PSR B1937+21 is also a BSS, the positive charge on side II is carried by positrons but not ions, then the emission pattern on both sides should be the same, and each beam may be symmetric around the magnetic axis. In this case, when symmetrical geometry is assumed for the radio emission, PC model could not account for the observational facts, otherwise, PC model may still be valid.

In general, further research on the pulsar’s magnetosphere structure and emission mechanisms will be helpful to understand the origin of its X-rays for PSR B1937+21. We are expecting that future polarization observations could provide more confirmative value of $`(\mathrm{d}\psi /\mathrm{d}\varphi )_{\mathrm{max}}`$, which is meaningful not only for constraining whether this pulsar is an NS or an SS but also for a better understanding of its X-ray emission.

Acknowledgments:

We are grateful to Prof. J.A. Gil and Dr. B. Zhang for their helpful comments and discussions. The valuable suggestions from an anonymous referee are also sincerely acknowledged. This work is supported by National Nature Science Foundation of China (10173002), and by the Special Funds for Major State Basic Research Projects of China.