text:
# Theory of shell structure and of the “magic” effect in spherical nuclei 11footnote 1Zh. Eksp. Teor. Fiz. 61, 1303–1318 (1971) [Sov. Pys. JETP 34, No. 4, 694-701 (1972)]

## 1 Introduction. Macroscopic treatment of magic nuclei

The study of the nucleus considered as a macroscopic body originated with the well known work of von Weizsäcker . Recognition of the distinction between volume and surface contributions to the binding energy was an important step in understanding the saturation properties of nuclear forces and served as the basis for current ideas regarding uniform nuclear matter (see, for example, Landau and Lifshitz ). The subsequent inclusion of the electrostatic energy and addition of the symmetry term taking into account the difference in the numbers of protons and neutrons  gave these ideas a rather high degree of refinement which permitted the total binding energies of a very large number of nuclei to be described with high relative accuracy. This, however, only deepened interest in possible deviations from the Weizsäcker formula, and the question of these deviations has arisen repeatedly in one form or another during the past thirty-five years. The most easily perceived cause is the quantum-mechanical discreteness of elementary particles, i.e., simply the obvious fact that protons and neutrons can be added to a nucleus only in integral numbers. This, naturally has stimulated the consideration of other similar deviations as having a more or less clearly expressed quantum nature.

However, the question of which deviations from the Weizsäcker formula have a macroscopic origin lies essentially in a different plane. At the present time we have no serious basis for rejecting such a possibility. On the opposite, developments in recent years indicate that macroscopic quantum phenomena are rather typical for nuclei which are not too light. In the first place, the simplest model of a uniform body of this type with limited size—the independent-particle model with a constant general potential in the internal region—permits one to make rigorous calculation of several important characteristics. The instability demonstrated by Nosov  for a spherical configuration of this model has an explicitly macroscopic nature. It is true that determination of the true stable equilibrium shape of the body encounters serious difficulty that the problem has no small parameter. The attempts to ignore this fact, with resort to direct numerical calculations, are difficult to analyze and we will not set for ourselves this goal. However, as often happens in similar situations, experimental data come to our assistance. The “collective” macroscopic order of magnitude of quadrupole moments of non-spherical nuclei is well known; an individual quasiparticle could not produce a quadrupole moment of this order. Closely connected with this fact is the usually smooth behavior of the deformation, i.e., the fact that this macroscopic characteristic of an non-spherical nucleus is practically independent of the individual quantum state of the last nucleon (quasiparticle) of the Fermi distribution.

Also of great significance is the study, intensely carried out during the last fifteen years, of the rotation of the nucleus as a whole, and also of oscillations of the nuclear deformation. Satisfaction of the laws of quantum mechanics does not diminish, of course, the macroscopic nature of this type of motion. It is remarkable that, even from calculations based on a simple model , it can be clearly seen how the characteristics of the corresponding collective degree of freedom (for example, rigidity of the spherical configuration with respect to deformation) “develop” with increase of the number of particles in the nucleus. The macroscopic properties of a body should lose their meaning and physical content in the transition to too small number of particles.

The undoubted existence of spherical nuclei is not compatible with the primitive model of ref. , and the real situation in a nucleus must be more complex. The most important reason for that is probably the residual interaction between the quasiparticles, which is unavoidable in a body of finite size. The loss of equilibrium non-sphericity is difficult to imagine without some qualitative change in the nuclear structure. However, as was shown by the theory of the corresponding phase transition developed by Nosov , this complication of the physical picture only increases the role of macroscopic quantum effects. Actually, our attention is primarily attracted not only by the singularities of thermodynamic quantities (as occurs for ordinary temperature-dependent second-order phase transitions but also by the singularities of the macroscopic-body ground-state wave function (which depends parametrically on the deformation), observed at the Curie point (see, for example, Eq. (28) in ref. )222For ordinary condensed matter we can also give a physical example of a temperatureless second-order phase transition. When a single-domain ferromagnetic particle grows in size, beginning at some critical size, a smooth change occurs in the spin directions, which gradually destroys the uniformity of magnetization over the volume. The theory given by Frei et al.  permits wave functions of the phases to be determined near a transition point. It is easy to see that their properties are consistent with the semi-phenomenological reasoning developed by Nosov  for the case of a nucleus..

It is possible that the significance of macroscopic quantum effects in a nucleus is still somewhat underestimated. In any case, against the background of the progress described above, the widely held opinion that the so-called magic numbers are clearly of non-macroscopic origin  appears to be more a dogmatic assertion than a fact following directly from experiment. It is natural to relate a theoretical study of the macroscopic properties of a medium with the experimentally observed singularities of the thermodynamic variables. This facilitates observation and detailed investigation of characteristic phenomena such as phase transitions. A specific feature of a nucleus is the direct accessibility of absolute zero temperature (the ground state), and the most important thermodynamic variable turns out to be the corresponding energy, i.e., the nuclear mass. As far as can be judged from the experimental data shown schematically in Fig. 1, the singularities under discussion are apparently breaks, i.e., discontinuities in the derivative of the energy with respect to the number of particles. Where the curve becomes sluggish, almost horizontal, on one side of the singularity, we are entering the region of non-spherical nuclei, and the Curie point corresponds to the phase transition discussed by Nosov , a consequence of which, in particular, is a change in the equilibrium shape of the nucleus. In addition, in Fig. 1 we can clearly see the magic cusps, i.e., downward spikes at whose vertices are located the corresponding magic nuclei. We emphasize the difference of this singularity from the ordinary Curie point. A qualitative difference between the phases—different values of the chemical potential—arises only near an isolated transition point. Away from this point it is not possible to point to a quality which would be possessed, say, by one of the phases while the other would not have it. The actual state of the body apparently does not undergo a discontinuity (with macroscopic accuracy, of course; see the following section), so that in other respects we appear to be dealing with a second-order transition.

In opposition to the impression created by Fig. 1, the widely used explanation of magic numbers  appeals directly to the pattern of filling by fermions of $`(2j+1)`$-fold degenerate single-particle levels in the field of some spherically symmetric potential. After a shell is filled, the next nucleon approaches the lower edge of the region of continuum states and the nucleon binding energy $`\epsilon `$ (the chemical potential taken with a reversed sign) correspondingly drops. The weak point of this type of interpretation is that a similar situation should arise also after the filling of each level individually, and not only of an entire definite group of levels (a shell). In comparison with the gaps separating neighboring shells, the distances between the levels inside each shell cannot always be small; there is no way in which such a small parameter could be obtained in this case. Specific calculations of single-nucleon level schemes confirm the validity of this reasoning (see, for example, the neutron level scheme from the monograph by Hodgson ). As a rule, the gap between shells is only two or three times greater than the typical distance between levels, and we would have to find a whole series of closely spaced cusps corresponding to the filling of each of the levels. In reality, for not too light nuclei, beyond the magic number 28, only the cusps shown in Fig. 1 are observed, which correspond to the much less frequently encountered “true” magic numbers. It is symptomatic that even a small, far from qualitative, excess of the gaps between shells over the level spacing inside the shells is achieved by selection of a rather substantial number of parameters. The gaps following shells 2 and 8 turn out to be large for practically any reasonable potential , but for the magic numbers 28, 50, 82, and 126 it is already impossible to make this statement. The fact that even the average field acting on the nucleon is determined here by no less than four parameters (the potential-well radius and depth, the width of the transition layer at the surface, and the intensity of the spin-orbit interaction within its boundaries) naturally provides additional cause for skepticism. In final analysis, we can attempt to interpret almost any phenomenon as the result of appropriate selection of fitting parameters; however, we have a right to expect more from the theory. Therefore, the widespread explanation of magic numbers produces seems to be somewhat artificial333Here and subsequently we will refrain from making any far reaching parallels with the shell structure of an atom. The principal distinction from nuclear matter is due to the Coulomb field of the nucleus at the center, as a result of which the atom has a sharply expressed spatial non-uniformity..

The key to understanding the macroscopic nature of the magic cusps and the nature of the phenomenon itself, which does not depend on specific numerical values of the nuclear parameters (equilibrium density of nuclear matter, structure of the transition surface layer, value of spin-orbit coupling, and the like), lies in a well-known general property of Fermi systems. The thermodynamic variables of these systems (for example, the total magnetic moment, the energy, or number of particles $`N`$) usually reduce to sums over a large number of fermion-filled states characterized by definite sets of quantum numbers. The simple replacement of summation over quantum numbers by integration is a crude treatment which leads to a result depending smoothly on such variables as the limiting momentum $`k_f`$ of the Fermi distribution or the intensity of the applied magnetic field. However, a rigorous calculation of the sum by means of Poisson’s formula

$$\underset{n=0}{\overset{\mathrm{}}{}}\phi (n)=\frac{1}{2}\phi (0)+_0^{\mathrm{}}\phi (n)𝑑n+\underset{\nu =1}{\overset{\mathrm{}}{}}_0^{\mathrm{}}\left(e^{i2\pi \nu n}+e^{i2\pi \nu n}\right)\phi (n)𝑑n$$
(1)
shows that the thermodynamic variables, generally speaking, undergo oscillations as a result of the third term in (1). Here the main contribution is from the quasiparticles closest to the Fermi surface. Actually, the oscillating integral which occurs inside the summation over $`\nu `$ would be extremely small if $`\phi (n)`$ were a sufficiently smooth function. The latter, however, contains as a factor the statistical distribution of quasiparticles, which changes sharply in the vicinity of the Fermi surface. An increase in temperature smooths the $`\phi (n)`$ function, as a result of which the oscillating effect falls off exponentially. A well-known specific example of this phenomenon is the de Haas-van Alphen effect , i.e., quantum oscillations of the magnetic susceptibility of metals at low temperatures.

Since we have in mind first of all nuclear applications, it is particularly important to emphasize that, in addition to the sharpness of the Fermi distribution, there is another equally necessary condition. An oscillation clearly does not arise with only a single summation according to Eq. (1), i.e., if $`n`$ is the only quantum number of the particles being discussed (in the general case, of quasiparticles). An elementary calculation shows that in the one-dimensional case the last term in (1) describes only the aforementioned trivial fact of addition of particles one at a time; physically this corresponds to the absence of a macroscopic effect444A similar situation occurs also in a three-dimensional potential well of sufficiently irregular shape so that the particle energy is its only quantum number.. An important role is therefore played by the fact that the quasiparticle in the spherical nuclei discussed below has an orbital angular momentum $`l`$. From the calculations of the next section it is evident that there are enough quantum numbers in a spherical nucleus for occurrence of oscillations of a macroscopic scale. However, upon formal application of relation (1), there can arise also terms which are not macroscopic in nature, and which must be discarded. In making this distinction the decisive criterion is

$$\rho _f=k_fR1$$
(2)
($`R`$ is the nuclear radius); for brevity we will also refer to this dimensionless parameter as the limiting momentum of the Fermi distribution.

To avoid misunderstanding, we note that the uses of Eq. (1) in solid-state physics and in nuclear physics have somewhat different aims. In the study of metals by means of the de Haas van Alphen effect, problems of sharpness of the Fermi surface essentially do not arise, and the problem is to determine the shape of the Fermi surface in momentum space (more accurately, in quasi-momentum space) . On the other hand, in the spherical nuclei of interest to us, the obvious isotropy removes the question of the “shape” of the Fermi limit, but its detailed structure is nowhere near as clear, owing to the residual interaction (see above). Furthermore, it is apparently only because of the residual interaction that spherical nuclei in general exist—we noted above that in a model without the interaction  a sphere is unstable. However, the effect of this circumstance on the energy of the body as a whole and on the structure of its quantum state in the region of rapid fall-off of the quasi-particle distribution function, and also the limiting behavior in accordance with which the residual interaction may disappear with increasing nuclear radius, represent a field which is still little studied. Therefore, in the calculations which follow, we will take tentatively the ordinary Fermi distribution function unperturbed by the residual interaction. In the current state of the problem, only experiment can answer the question as to how closely we approach the ideal limit of a step-shaped distribution for the nuclear sizes achievable at the present day. Sections 3 and 4 of the present article are devoted to questions bearing on this point.

For what follows it is very important to comment on the various characteristics of the quasiparticle and the possible relation between them. The point is that in comparison with the primitive model of a Fermi gas in an external field, a real Fermi liquid (i.e., a bound system of fermions of finite size) has some qualitative differences in this respect. We have attempted to compare these situations schematically in Figs. 2a and 2b, respectively. In both cases we can speak of an energy, if we reckon it from the zero of the kinetic energy of an external free nucleon. At the limit of the Fermi distribution this quantity reduces to the chemical potential $`\epsilon _f=\epsilon `$, where $`\epsilon `$ is the reversible work function of the particle for removal from the nucleus (the nucleon binding energy). In addition, in the case corresponding to Fig. 2a, there exists also an energy $`\epsilon ^{}`$ relative to the bottom of the potential well. It determines the limiting momentum $`\rho _f\sqrt{\epsilon _f^{}}`$.

Turning to a nuclear Fermi liquid (see Fig. 2b), we encounter the, strictly speaking, inadequate nature of the concept of the “bottom” of the potential well. Actually, a situation with noninteracting quasiparticles can exist only in the immediate vicinity of the Fermi limit. Even if we have in some sense complete data on the energy spectrum of the nucleus (we assume that the energy $`E\{n_𝐤\}`$ of the body is known as a functional of the occupation numbers $`n_𝐤`$), it is not possible to give a reasonable definition of the quasiparticle energy $`\epsilon ^{}`$ relative to the “bottom.” Nevertheless, the idea of a limiting momentum still makes sense. In principle we could deduce the value of $`\rho _f`$, say, from the form of the wave function of the last quasiparticle.

However, if the “bottom,” strictly speaking, does not exist, and the quasiparticle energy $`\epsilon ^{}`$ calculated with respect to it is not an adequate quantitative notion, then what physically determines the limiting momentum $`\rho _f`$? The answer to this question is given by the major assumption on which the theory developed in the next section is based.

## 2 Theory of the shell structure of a spherical nucleus

As a starting point for what follows we will assume that the value of $`\rho _f`$ is completely determined by the total number of particles $`N`$. It is assumed that the function $`N(\rho _f)`$ does not undergo oscillations555For the model of an ideal gas located in an externally applied potential well, the oscillations would be given by Eq. (16).. In a purely formal sense, only differentiability of this function is required (see below, the transition from Eq. (22) to Eq. (24)). However, in connection with criterion (2), it is natural to be interested in its asymptotic representation for large $`\rho _f`$. In that case the corresponding expansion should be terminated with a finite number of terms:

$$N(\rho _f)=\frac{4}{9\pi }\rho _f^3s\rho _f^2+q\rho _f.$$
(3)
Actually the fifth term of the expansion, which is $`\rho _f^1`$, would give physically meaningless fractional additions to the number of particles. In a systematically macroscopic treatment it is obviously required to discard also the fourth term as describing the single-particle effect (see the discussion of this point in the Introduction). In what follows we will also drop all similar non-macroscopic terms of the type $`\rho _f^01`$666The correctness of the third term of Eq. (3) also could turn out to be questionable, since $`2l+1\rho _f`$ particles fill a single level in a spherically symmetric field. However, the unavoidable zero-point oscillations of the deformation have a scale $`\mathrm{\Delta }\alpha \rho _f^2`$, and the corresponding shift in particle energy is of the same order of magnitude. The degeneracy is as a result lifted to a sufficient degree, while, on the other hand, such small deformations do not yet violate the conservation of the integral of motion $`l`$ (see ref. ). Thus, the averaging of the function $`N(\rho _f)`$ implied in (3) is achieved automatically in a real situation..

The assumptions described above have a great similarity to the system of assumptions of the theory of an unbounded Fermi liquid . In particular, the coefficient in the first term of the right-hand side of (3) corresponds (if the additional spin doubling is taken into account) to the volume contribution to the number of cells of phase space, that is, it is equal to the corresponding expression for an ideal Fermi gas. However, we have no basis for assuming that the remaining terms are also independent, to the same degree, of the law of interaction between particles. On the contrary, the surface term $`s\rho _f^2`$, and also the curvature term $`q\rho _f`$, take into account the structure of the transition layer on the nuclear surface, spin-orbit coupling, and so forth. Numerical values of the coefficients $`s`$ and $`q`$ must be determined experimentally (see the following section).

We represent the ground-state energy of a spherical nucleus in the form of the sum

$$E=E_0+E_1$$
(4)
of a smooth part $`E_0`$ and an oscillating part $`E_1`$. In the calculation of the latter, it is convenient to express it initially in terms of the variable $`\rho _f`$. The coupling noted in the Introduction between the oscillations and quasi-particles approaching the Fermi limit permits the problem to be reduced to the model  of a gas in a well with a constant potential inside it (the last fact corresponds physically to uniformity of the spatial distribution of matter in the internal region of the nucleus). Subsequent use of Eq. (1) to carry out the summation over the quantum numbers $`l`$ and $`n`$ will lead to integrals of two substantially different types. Those of them which are due to practically the entire region $`0<\rho =kR<\rho _f`$ of values of the particle wave number do not permit generalization to a real Fermi liquid. However, they depend smoothly on $`\rho _f`$ and are not of interest here—the problem of explicit calculation of the function $`E_0(\rho _f)`$ does not arise here. The integrals which have an oscillating dependence on $`\rho _f`$ converge rapidly at $`\rho \rho _f`$, and the generalization of their contribution to the case of a Fermi liquid is obvious. This fact even allows us to limit ourselves to calculation of the number of particles

$$\stackrel{~}{N}=\stackrel{~}{N}_0+\stackrel{~}{N}_1$$
(5)
in this model. The point is that, in the terms which determine $`\stackrel{~}{N}_1`$, the energy of a single particle (subsequently a quasi-particle) would anyway be taken outside the integral sign as a slowly varying function. Equation (5) corresponds to the quantity

$$\phi (l,n)=(2l+1)w_f(l,n),$$
(6)
which must be summed over $`l`$ and $`n`$. The Fermi distribution

$$w_f(l,n)=\{\begin{array}{c}1\mathrm{for}\rho _{ln}<\rho _f,\hfill \\ 0\mathrm{for}\rho _{ln}>\rho _f,\hfill \end{array}$$
(7)
depends on these same quantum numbers777In view of the macroscopic nature of the effect being studied, the existence of spin in the nucleon will be taken into account easily by subsequent doubling of the result. We note that in problems of this type it is more natural to begin the enumeration of the principal quantum number $`n`$ not from unit but from zero. For example, the $`1s`$ state is assigned to $`n=0`$. In particular, this provides the possibility of formally symmetric application of Eq. (1) to summation over the two quantum numbers (see Eq. (11) below). through the eigenvalues $`\rho =\rho _{ln}`$. The Bohr-Sommerfeld quantization rule  serves to determine the eigenvalues,

$$_a^Rk_l(r)𝑑r=\pi (n+\gamma ),$$
(8)
where the internal turning point $`r=a`$ is due to the centrifugal barrier, and $`\gamma <1`$ defines an additional phase shift which depends on the nature of the boundary conditions. It is not necessary to calculate integral (8) from the beginning, since in this case both the wave functions themselves (spherical Bessel functions) and their quasi-classical asymptotic forms are well known. Using a notation close to that of ref. , we have

$$\begin{array}{c}\hfill \rho (\mathrm{sin}\beta \beta \mathrm{cos}\beta )=\pi (n+3/4),\\ \hfill \beta =\mathrm{arccos}\frac{l+1/2}{\rho }=\mathrm{arcsin}\sqrt{1\frac{(l+1/2)^2}{\rho ^2}}.\end{array}$$
(9)
The Jacobian of the transformation to the new variables also was determined in ref. :

$$dndl=\frac{1}{\pi }\rho d\rho \mathrm{sin}^2\beta d\beta .$$
(10)
Two-fold summation of the function (6) according to Eq. (1) gives

$$\begin{array}{c}\hfill \underset{n=0}{\overset{\mathrm{}}{}}\underset{l=0}{\overset{\mathrm{}}{}}\phi (l,n)=\frac{1}{4}\phi (0,0)+_0^{\mathrm{}}\phi (l,n)𝑑l𝑑n\\ \hfill +\underset{\nu =1}{\overset{\mathrm{}}{}}_0^{\mathrm{}}\left(e^{i2\pi \nu n}+e^{i2\pi \nu n}\right)\phi (l,n)𝑑l𝑑n+\underset{\lambda =1}{\overset{\mathrm{}}{}}_0^{\mathrm{}}\left(e^{i2\pi \lambda l}+e^{i2\pi \lambda l}\right)\phi (l,n)𝑑l𝑑n\\ \hfill +\frac{1}{2}\left\{_0^{\mathrm{}}\phi (0,n)𝑑n+\underset{\nu =1}{\overset{\mathrm{}}{}}_0^{\mathrm{}}\left(e^{i2\pi \nu n}+e^{i2\pi \nu n}\right)\phi (0,n)𝑑n\right\}\\ \hfill +\frac{1}{2}\left\{_0^{\mathrm{}}\phi (l,0)𝑑l+\underset{\lambda =1}{\overset{\mathrm{}}{}}_0^{\mathrm{}}\left(e^{i2\pi \lambda l}+e^{i2\pi \lambda l}\right)\phi (l,0)𝑑l\right\}\\ \hfill +\underset{\nu =1}{\overset{\mathrm{}}{}}\underset{\lambda =1}{\overset{\mathrm{}}{}}_0^{\mathrm{}}\left[e^{i2\pi (\lambda l+\nu n)}+e^{i2\pi (\lambda l+\nu n)}+e^{i2\pi (\lambda l\nu n)}+e^{i2\pi (\lambda l\nu n)}\right]\phi (l,n)𝑑l𝑑n.\end{array}$$
(11)
All terms in the right-hand side of (11) permit analysis and classification on the basis of the criteria mentioned above. The non-macroscopic nature of the first of these terms is obvious. The non-oscillatory behavior of the second term follows directly from its form. After carrying out the integration, we can verify that the third term also provides no macroscopic contribution. As the result of integration on the basis of the same formulas (9) and (10) and subsequent summation over $`\lambda `$, we find that the fourth term is $`\rho _f/6\pi `$. It depends smoothly on the limiting momentum. The same can be said of the two pairs of following terms, which are enclosed in curly brackets. In both cases the sum which adds to the integral actually describes its increase as the result of those levels which, with further advance of the Fermi limit, appear separately for $`\rho <\rho _f`$ 888This rather formal fact has already been mentioned in the Introduction and commented on in footnote 5..

The oscillations are described by the double sum over $`\lambda `$ and $`\nu `$ (the last term in Eq. (11)). Before we begin calculation of this sum, we will make one remark of a less formal nature. In addition to the Fermi limit $`\rho =\rho _f`$ (see Eq. (7)) the distribution of particles along the $`l`$ (or $`\beta `$) axis is also cut off—there are no negative orbital angular momenta. An intuitive interpretation of this quantum phenomenon must lie in the fact that the shell oscillations depend on the relative position of these two limits. As will be seen from Eq. (19) below, a magic nucleus arises when the two limits coincide. In view of this role of the lower limit of the angular momentum scale, where $`\beta \pi /2`$, it is convenient to introduce an additional angle

$$\stackrel{~}{\beta }=\pi /2\beta .$$
(12)
Then it is easy to see that the integral under the summation sign, generally speaking, has the same structure as in the second term of the right-hand part of Eq. (11), i.e., it is non-macroscopic. However, a saddle point arises for a definite ratio between $`\lambda `$ and $`\nu `$ at the lower limit of the angular momentum axis, and this provides a macroscopic contribution 999This pertains to the first two exponential terms under the integral sign. Some other terms of the double sum have saddle points at intermediate angular-momentum values, which do not affect the result. It is interesting to note that in the preceding work devoted to calculation of the rigidity of the spherical configuration of the simplest model  these same points played the role of poles of the integrand, and incidentally gave likewise contribution..

If we take into account Eqs. (9) and (12), we can write the expansion of the argument of the exponential in powers of $`\stackrel{~}{\beta }`$ up to and including quadratic terms as follows:

$$2\pi (\lambda l+\nu n)\pi \lambda \frac{3}{2}\pi \nu +2\nu \rho +\pi (2\lambda \nu )\rho \stackrel{~}{\beta }+\nu \rho \stackrel{~}{\beta }^2.$$
(13)
Vanishing of the linear term (the condition for existence of a saddle point at $`\stackrel{~}{\beta }=0`$) requires that

$$\nu =2\lambda .$$
(14)
Since the phase-space element is given by Eq. (10), we obtain from Eqs. (6), (7), (9), and (12)

$$e^{i2\pi (\lambda l+\nu n)}\phi (l,n)𝑑l𝑑n\frac{2}{\pi }^{\rho _f}𝑑\rho \rho ^2e^{i4\lambda \rho }_0e^{i2\lambda \rho \stackrel{~}{\beta }^2}\stackrel{~}{\beta }𝑑\stackrel{~}{\beta }=\frac{\rho _f}{8\pi }\frac{e^{i4\lambda \rho _f}}{\lambda ^2}.$$
(15)
Addition of the complex conjugate of the expression and summation over the single remaining free index leads to

$$\stackrel{~}{N}_1=\frac{\rho _f}{4\pi }\underset{\nu =1}{\overset{\mathrm{}}{}}\frac{\mathrm{cos}4\nu \rho _f}{\nu ^2}$$
(16)
according to Eqs. (5) and (6). For conversion to the oscillating part of the energy of the nucleus $`E_1`$ it is necessary to take into account the nucleon spin by doubling, and also to multiply Eq. (16) by the energy $`\epsilon _f=\epsilon `$ ($`\epsilon `$ is the nucleon binding energy) of a single quasiparticle:

$$E_1=\epsilon \frac{\rho _f}{2\pi }𝔐(\rho _f).$$
(17)
Thus, for an absolutely sharp step-shaped Fermi limit for the quasiparticles, the shell effects in a spherical nucleus are described by a certain universal periodic function

$$𝔐(\rho _f)=\underset{\nu =1}{\overset{\mathrm{}}{}}\frac{\mathrm{cos}4\nu \rho _f}{\nu ^2}$$
(18)
whose plot is shown in Fig. 3. Its derivative has discontinuities at values of the argument

$$k_fR=\frac{\pi }{2}p,p=2,3,4,\mathrm{},$$
(19)
which correspond to the magic cusps (see Eq. (17)). In other words, Eq. (19) is a unique quantization rule for magic values of the parameter $`\rho _f=k_fR`$.

We emphasize that the theory being developed is an expansion in inverse powers of the integer “number” of the magic nucleus $`p`$ (see also Eq. (2)). It cannot pretend to give a quantitative description of the case of extremely small values, especially since it is then evidently impossible in general to make any reasonable separation of “smooth” and “oscillating” effects. However, the number of nodes of the radial wave function of the nucleon is such a stable characteristic that we can hope to trace it, at least qualitatively, to the lightest nuclei. Since $`\rho _f=\pi `$ already corresponds to a $`1s`$ state (the doubly magic nucleus $`{}_{2}{}^{}\mathrm{He}_{2}^{4}`$ corresponds to this physically), the region of possible values of $`p`$ should begin with 2 101010The definition of the effective radius $`R`$ of the nucleus implied in (19) concerns the internal structure of the nucleus exclusively. Since we have reduced the problem to a model with an impenetrable wall, we must represent the wall as constructed where the wave function of the corresponding quasiparticle extrapolated from the internal region approaches zero. In other words, this effective boundary of the nucleus can always be disposed so that the additional phase shift from the Bohr-Sommerfeld rule (8), in accordance with (9), will take on the value $`\gamma =3/4`$ for the quasiparticles playing the principal role. In close connection with this fact, the equations expressing the shell oscillations in terms of $`\rho _f`$ are universal; they do not depend on the spin-orbit interaction or on the structure of the surface layer where. After transformation to the $`N`$-scale this universality is lost (see also Eq. (3) and the explanation for it)..

It is easy to verify that the right-hand part of Eq. (18) is a Fourier series for an elementary function which we will write out explicitly in a form valid for the two periods adjacent to the magic number $`p`$:

$$𝔐(\rho _f)=\frac{\pi ^2}{6}2\pi \left|\rho _f\frac{\pi }{2}p\right|+4\left(\rho _f\frac{\pi }{2}p\right)^2,\frac{\pi }{2}(p1)<\rho _f<\frac{\pi }{2}(p+1).$$
(20)
The absolute-value sign reflects the non-analyticity of the function at the magic cusp. Using the indices $`+`$ and $``$ to distinguish the values of the discontinuous function to the right and left respectively, we have

$$\frac{d𝔐}{d\rho _f}|_\pm =2\pi .$$
(21)
Now, returning to (17), it is not hard to obtain with the required accuracy an expression for the discontinuity in the derivative of the oscillating part of the energy of the nucleus:

$$\mathrm{\Delta }\left(\frac{dE_1}{d\rho _f}\right)=\rho _f(\epsilon _++\epsilon _{})=2\rho _f\overline{\epsilon }.$$
(22)
In order to convert from the limiting momentum to the true number of particles (3), it is necessary to multiply both sides by $`d\rho _f/dN`$. At the same time we recognize that the smooth component $`E_0`$ does not have a singularity at the cusp, so that Eq. (22) actually gives the discontinuity in the derivative of the entire energy $`E`$:

$$\frac{d\rho _f}{dN}\mathrm{\Delta }\left(\frac{dE}{d\rho _f}\right)=2\rho _f\frac{d\rho _f}{dN}\overline{\epsilon }.$$
(23)
This relation for the discontinuity in the chemical potential can be given a more compact form:

$$\mathrm{\Delta }\left(\frac{dE}{dN}\right)=\frac{\overline{\epsilon }}{dN/d(k_fR)^2}.$$
(24)
This can also be considered as a formula for the discontinuity in nucleon binding energy

$$\epsilon =\frac{dE}{dN}$$
(25)
in the vicinity of a magic nucleus. Then, obviously, we must assume $`\mathrm{\Delta }\epsilon \mathrm{\Delta }\left(dE/dN\right)=\epsilon _{}\epsilon _+`$.

## 3 Comparison with experiment

The most direct physical meaning is attached to the number-of-particles scale, as the most closely connected with experiment. We will discuss this after we have determined for the nuclear Fermi liquid, from experimental data, the numerical values of the coefficients $`s`$ and $`q`$ entering into Eq. (3). If we assume that, of the magic numbers found experimentally up to the present time, the most important in this connection are 82 and 126, then according to (19) and (3) these are sufficient to determine the two parameters. We finally obtain

$$s=1.1,q=6.8.$$
(26)
It is interesting to note that the surface term exceeds by only a factor of about two the value $`\stackrel{~}{s}=1/2`$ characteristic of an ideal Fermi gas adjacent to an impenetrable wall . With the values in (26), Eqs. (3) and (19) for the number of particles in a magic nucleus take the form

$$N=0.548p^32.79p^2+10.7p.$$
(27)
The results of a calculation on the basis of Eq. (27) are compared with the known magic numbers in Table I. Down to $`N=28`$ inclusive, the agreement must be considered satisfactory. For the lightest magic nuclei there is not even qualitative agreement 111111The frequently encountered separation of a single $`1f_{7/2}`$ level in an individual shell has always been problematical. The reasoning and formulas of the preceding section create the impression of rejecting the magic number 20 which arises here. It is necessary, however, to make the reservation that magic phenomenon can in any case not be regarded as macroscopic when it comes to the lightest nuclei..

If we take into account the two-component nature of nuclear matter, the oscillating term is given by the sum of expressions (17) for neutron and proton quasiparticles

$$E_1(N,Z)=\epsilon _N\frac{\rho _f^N}{2\pi }𝔐(\rho _f^N)\epsilon _Z\frac{\rho _f^Z}{2\pi }𝔐(\rho _f^Z).$$
(28)
We note in this connection that the agreement observed so far of the corresponding magic numbers evidently indicates practically identical values of the parameters $`s`$ and $`q`$ for the two components. We may suppose that this is due to the relative smallness of the effects distinguishing protons and neutrons in the nucleus, such as, for example, the Coulomb curve (the radial dependence of electrostatic potential in the internal region of the nucleus) 121212These remarks have pertained to the total number (3) of nucleons of a given kind in the nucleus. However, in regard to the more detailed properties of the residual interaction of quasiparticles close to the Fermi limit, the neutron and proton components may be far from identical; see below..

Clarification of the nature of the singularity in the energy $`E(N,Z)`$ at the point of occurrence of a doubly magic nucleus presents substantial theoretical interest. If we return to expression (20) and neglect in it the last quadratic term, we see that, according to (28), the surface $`E_1(N,Z)`$ has here the form of a pyramid of rhombic cross section with its axis directed vertically downward. As far as we can judge from the experimental data, the shape of the intersection of the mass surface with the plane $`E_1=\mathrm{const}`$ is actually close to a rhombus near a doubly magic nucleus; see, for example, ref. .

Incidentally, the separation of the oscillating part $`E_1`$ of the mass, from a purely practical point of view, suffers from a certain ambiguity. Therefore our attention is drawn to the relation (24), which does not depend on this procedure, for the discontinuity in binding energy of a nucleon of the corresponding type in a singly magic nucleus. Comparison with mass data (see, for example, ref. ) has shown that the theory based on the step-shaped distribution (7) does not agree with experiment for the magic numbers found up to the present time. Since the oscillations are due to the sharpness of the Fermi distribution, it is necessary to have in mind that the jump in quasiparticle occupation numbers from unity to zero is the maximum permissible by the Pauli principle. Consequently the oscillations calculated according to Eq. (7) are exaggerated, and the value $`(\mathrm{\Delta }\epsilon )_0`$ calculated from Eq. (24) actually exceeds everywhere the really observed binding-energy jump $`\mathrm{\Delta }\epsilon `$. It is natural to characterize this discrepancy by the quantity

$$\omega =\frac{(\mathrm{\Delta }\epsilon )_0}{\mathrm{\Delta }\epsilon }1,$$
(29)
which would approach zero in the ideal case (7) of absence of residual interaction between quasiparticles. Values of $`\omega `$ for 52 magic nuclei are listed in Table II 131313The well known even-odd oscillations of nuclear masses are not related to the phenomenon being discussed. However, the widely used methods of calculating the corresponding correction (see, for example, ref. ) are purely empirical and can hardly be considered reliable. We therefore selected only those cases in which there are data on the binding energy of two nucleons. On removal of two nucleons of the same kind, the type of nucleus does not change, and the even-odd correction to the mass drops out of the result..

Each magic number is associated with several data differing in the number of nucleons of different kinds. Although some systematic behavior is apparently observed in the dependence on the latter, we will assume tentatively in each case that these are deviations of an accidental and uncorrelated nature. The corresponding average values of $`\omega `$ are shown in Fig. 4, together with the corresponding variances. The rapid decrease of $`\omega `$ in the transition to heavy magic nuclei reflects the naturally expected drop in strength of the residual interaction with increasing nuclear size. The electrically charged protons are less affected by the residual interaction from the very start, and at the end of the periodic table it shows a sharper drop than in the case of neutrons. We will dwell briefly on the attempts to provide some simple empirical description of the behavior of the points plotted in Fig. 4. From the number of curves that asymptotically approach $`\omega =0`$, the possibility $`\omega =C/p^3`$ is suggested. For neutrons this gives $`C=336\pm 12`$; according to the criterion of likelihood  of the law $`\omega 1/p^3`$, the probability of $`\chi ^2`$ greater than the value observed in the present case is 20%. However, in the case of protons the corresponding probability of a similar law drops to 0.5%. Another possibility is that the proton curve $`\omega (p)`$ intersects the axis of abscissas, and we will discuss this in the next section.

## 4 Discussion

The problem of the energy spectrum of a quantum liquid consisting of fermions is rather complex in itself. In application to such a peculiar object as a nucleus it becomes particularly difficult and many-faceted. The accuracy of the description of its excitations in the model of an ideal gas of quasiparticles is limited by the residual interaction between them. It is natural to suppose that it is created only by the finite size of the system, and that in the transition to unbounded nuclear matter we would obtain the simplest variety of Fermi-liquid spectrum . Another variant of the energy spectrum takes into account the possibility that, even in an infinite quantum liquid, residual phenomena such as the Cooper phenomenon  are preserved. We must suppose that the deviation shown in Fig. 4 from the asymptotic value $`\omega =0`$ is not accidental. In spherical nuclei, which owe their very existence to the residual interaction, the latter must of course smear in some way the boundary of the Fermi distribution, and in so doing suppress the oscillations. As a result of this special role and, apparently, of the relatively large magnitude of the effect, it will probably be possible to clear up the question of the residual interaction, which is finally and with difficulty yielding to analysis.

The monotonic drop observed so far (Fig. 4) in the function $`\omega (p)`$ already permits some preliminary conclusions to be drawn. The point is that the shell oscillations calculated in the present article are an extremely fine instrument for analysis of the quasiparticle distribution near the limit $`\rho \rho _f`$, and with increasing nuclear size the accuracy of this method is increased. In fact, according to (1), (9), (11), and (15), the oscillations penetrate below the Fermi limit by a depth $`\mathrm{\Delta }\rho 1`$, i.e., the relative thickness of the layer of the Fermi distribution occupied by them is $`1/\rho _f1/R`$. On the other hand, in the Cooper phenomenon the width of the zone of smearing of the distribution of initial quasiparticles does not depend on the size of the system. It can be seen from this how greatly the oscillations would be weakened in the presence of such a phenomenon. The drop in the $`\omega (p)`$ curve would be delayed and, furthermore, it evidently would have the reverse behavior. If still heavier magic nuclei do not exhibit similar features and the monotonic nature of the $`\omega (p)`$ function is further confirmed, it will be necessary to give a negative answer to the question of the Cooper phenomenon (or any other similar instability of an absolutely sharp Fermi limit of quasiparticles in infinite nuclear matter).

In speaking of the residual interaction, it is necessary to recall that not all conceivable varieties of this interaction are, so to speak, physically real. At least some variants of the interaction may be excluded by a canonical transformation to new quasiparticles without loss of sharpness in the Fermi limit of their statistical distribution. The well-known case of repulsive forces of zero radius between fermions  serves as an example of this. This tendency of a repulsive interaction to be “renormalized” may turn out to be important in regard to the proton component of nuclear matter. Actually, the points referring to protons in Fig. 4 may reflect a competition between the nuclear and Coulomb interactions. If the latter become dominant with increasing nuclear size, this may be expressed in the intersection of the corresponding $`\omega (p)`$ curve with the abscissa. This intersection would be a termination point for the function $`\omega (p)`$, beyond which the quasiparticles would satisfy the ordinary Fermi distribution (7). On the other hand, according to the results of ref. , disappearance of the residual interaction between proton quasiparticles would not, apparently, favor the stable equilibrium nature of a spherical configuration. Even in the already studied region of nuclei beyond the doubly-magic lead, the data on the phase transition  appear rather symptomatic in this respect. Although the entire period from 82 to 126 contains 44 protons, addition to lead of only five of them turns out to be sufficient to remove the stability of the spherical nuclear shape. We may suppose that beyond the terminal point of the proton $`\omega (p)`$ curve the regions of existence of spherical nuclei narrow more rapidly or even will not exist at all.

We express our thanks to M.Ya. Amusya, A.I. Baz’, B.L. Birbrair, I.I. Gurevich, M.V. Kazarnovskii, L.P. Kudrin, I.M. Pavlichenkov, G.A. Pik-Pichak, V.P. Smilga, V.V. Sudakov, and O.B. Firsov for discussion of the results of this work. We also express our indebtedness to I.M. Lifshitz for pointing out to one of the authors (V.N.) the mathematical apparatus which apparently turns out to be most convenient for discussion of problems of this nature.

Figures captions

Fig. 1. The smooth part of the mass, described by the Weizsäcker formula, does not have singularities; it is taken as the reference zero of the “shell” effect pictured here. A more detailed and accurate plot would contain a whole set of curves corresponding to the different chemical elements (see, for example, the review of Myers and Swiatecki ). As a function of the number of protons $`Z`$ the behavior of the mass reveals no less sharply expressed features of a similar nature; see also Sections 3 and 4 of the present article.

Fig. 2.

Fig. 3.

Fig. 4. Values of the quantity $`\omega `$: $``$ — neutron magic numbers, $``$ — proton magic numbers.


format: text

text:
# Hysteresis of Finite Arrays of Magnetic Nano Dots.

## I Introduction

A ferromagnetic particle goes into a monodomain state if its size $`D`$ is below a critical value $`D_c=10100`$ nm. This is due to the competition between the exchange and dipolar energies. Therefore, a nanoparticle in a monodomain state may be viewed as a giant magnetic dipole with magnetic moment of thousands of Bohr magnetons. For an $`N\times N`$ array of well-separated nanoparticles the exchange energy is usually negligible in comparison with the dipolar and anisotropy and Zeeman energies. The study of such systems is of increasing importance because of their technological applications in data storage devices and magnetic field sensors. As the technology of these devices moves towards higher densities of stored information, it requires smaller particles of magnetic media SW , SSF , for which finite size effects become relevant. In finite arrays of such large dipole moment particles, the dipolar field of the array becomes comparable with the bulk anisotropy field. Dipolar effects in such systems affect the static and dynamics properties of the array; and thus must be taken into account.

Recently, Camley and Stamps in CS1 , CS2 , CS3  investigated the dynamics and magnetization processes of a finite planar array of $`N\times N`$ ferromagnetic nano dots, for $`N=3,4,5,6`$. The nano-dots were taken to interact only via the dipole-dipole interaction, and they were subject to an external field applied either along one side of the array or along its diagonal. They found rather complicated hysteresis loops with the magnetization reversal controlled by the shape anisotropy induced by the array itself. We have considered the same model, and have extended their results, for $`N=2\mathrm{}13`$. Our results for $`N=3`$ qualitatively agree with those of CS1 . We find that the behavior of these systems is surprisingly complex, both for small and for larger values of $`N`$, and we have studied the approach of these systems to $`N\mathrm{}`$ behavior.

In the present work each dot is taken to have a radius $`R_d`$, thickness $`d`$ and a single degree of freedom corresponding to the orientation of a magnet of saturation magnetization $`M_0`$. We consider only the case of zero temperature. The dots are arranged on a square lattice with lattice spacing $`a>2R_d`$, and the dots interact only via the dipole-dipole interaction. The equation of motion for the magnetic moment of each dot is governed by the Landau-Lifshitz-Gilbert equation (LLG), LL  which reads

$`{\displaystyle \frac{d𝐌}{dt}}=\gamma 𝐌\times 𝐇_{eff}\alpha {\displaystyle \frac{𝐌\times (𝐌\times 𝐇_{eff})}{M_s}}`$ (1)
where $`\gamma `$ is the gyromagnetic ratio, $`\alpha `$ is the damping coefficient, $`𝐌`$ is the magnetic moment of the dot and $`M_s=|𝐌|`$ is the saturation magnetization, $`𝐇_{eff}`$ is the average effective magnetic field acting on the dot. The average effective magnetic field acting on the $`i`$-th dot is due to the applied external field, the dipolar fields, and the anisotropy field

$`𝐇_{eff}^i=H_0\mathrm{cos}\theta \widehat{x}+H_0\mathrm{sin}\theta \widehat{y}𝐇_{dip}^i+2K_1{\displaystyle \frac{m_z^i}{M^2}}\widehat{z}.`$ (2)
Here the dipole field acting on the $`i`$-th dot due to all other dots in the array is given by

$`𝐇_{dip}^i={\displaystyle \underset{ji}{}}{\displaystyle \frac{𝐌_j}{r_{ij}^3}}3{\displaystyle \frac{(𝐌_j𝐫_{ij})𝐫_{ij}}{r_{ij}^5}}`$ (3)
The choice of anisotropy field is determined by the shape of the dot, which in our problem is directed along the symmetry axis of the (cylindrical) dots. We divide both sides of Eq.(1) by ($`\gamma M_s^2`$) and define a dimensionless time variable $`\tau =\gamma M_st`$. The LLG in these reduced units becomes

$`{\displaystyle \frac{d𝐦}{d\tau }}=𝐦\times 𝐡_{eff}{\displaystyle \frac{\alpha }{\gamma }}𝐦\times (𝐦\times 𝐡_{eff})`$ (4)
where $`𝐦=\frac{𝐌}{M_s}`$ and $`𝐡_{eff}=\frac{𝐇_{eff}}{M_s}`$. We employ a system of units where magnetic fields are given in units of $`M_s`$ and distances are given in units of the array’s lattice spacing $`a`$. The strength of the dipole field is characterized by $`h_{dip}=\frac{\pi R_d^2d}{a^3}`$, which is the ratio of the dot volume to the volume of a cubic cell with side $`a`$. For all arrays studied in this work, we take $`h_{dip}=0.5`$.

This article is organized as follows. Section II presents the numerical techniques and method of time integration. Section III presents an extensive discussion of magnetization processes and hysteresis for arrays of $`N\times N`$ nano dots ($`N=2,3,\mathrm{},13`$) when the external field is applied along one side of the array. Section IV considers the effects on the hysteresis loop of misalignment of the external field. Section V considers the relationship between the area of the hysteresis loop and $`N`$. A brief summary is given in section VI.

## II Numerical Techniques

We have employed two different approaches to study the magnetization processes of our $`N\times N`$ arrays of nano-dots. The first approach uses the second rank Runge-Kutta (RK) algorithm with fixed time step to integrate the LLG equation in reduced units. The second approach uses the “greedy algorithm“ to find the stable final state of the array. In both approaches the calculations were done using FORTRAN and Mathematica languages walf . The two approaches yielded similar results.

In the RK approach, the integration is done using a fixed time step $`\mathrm{\Delta }\tau =5\times 10^3`$ and a damping coefficient $`\frac{\alpha }{\gamma }=0.6`$, with an initial state in which the magnetic moment of the dot is randomly generated. The time integration proceeds until a stable final state is reached. The stability of the final state is checked by changes in the total energy of the system. Iterations are stopped when the difference between the total energy of the system from the ($`n+1`$)-th iteration and that of the $`n`$-th iteration is of the order of $`\mathrm{\Delta }E_n=10^5`$ in units of $`M_s^2a^3`$. Our solutions converged after almost $`10^3`$ iterations.

The greedy algorithm approach assumes the dots to be aligned along the direction of their total local field. In the initial state each dot is chosen to point randomly. Next, all components of the total local field are calculated for each dot. Finally each component of the field is normalized by the magnitude of the total local field to yield the new dot magnetization orientation. The final state of the $`n`$-th iteration is used as an initial one for the next iteration until the solution is converged. The convergence of the final state is checked in manner similar to that used in the first approach. Note that the greedy algorithm can be obtained from the LLG equation in the limit where the damping coefficient is very large. We find that the LLG approach converges faster than the greedy algorithm. This is probably related to a phenomenon known in RLC circuits, where a critically-damped circuit approaches equilibrium faster than a circuit with a larger amount of damping.

Calculation of the dipolar field at the $`i`$-th dot is the most computationally time-consuming aspect of both approaches since it requires summation of the dipole fields from all other dots in the array. However, due to the relatively small sizes of our dot arrays this calculation is performed rather quickly.

## III Hysteresis Loop and External Field Orientation Effect

Hysteresis loops $`M(H)`$ for $`N\times N`$ arrays of nano dots subject to an external magnetic field applied along one side have been calculated. Initially, a strong external field $`H_0`$ is applied to the array until saturation. Then the field is decreased to $`H_0`$, and then is increased again up to $`H_0`$. We have taken $`H_0=2M_s`$, and a fixed field-step of $`\mathrm{\Delta }H=2\times 10^3M_s`$ is used to simulate the sweeping process. For each value of the external field the system was iterated until a stable final state is reached. As seen in Fig. 1, where $`M(H)`$ is plotted (both in units of $`M_s`$), the odd and even $`N`$ arrays have somewhat different behaviors, especially for small $`N`$. One aspect of this is that the odd-$`N`$ systems display magnetization jumps as the field changes. The odd-$`N`$ and even-$`N`$ behavior becomes similar for larger values of $`N`$, something we study in a later section.

In Fig. 1, the angle of the field to one of the sides (the $`y`$-axis) is taken to be $`\theta =0`$. Experimentally, however, field misalignment is almost inevitable, so that we have also studied this phenomenon.

Fig. 2 shows results for $`N=3\mathrm{}6`$ and angles $`\theta =5^{}`$, $`30^{}`$, and $`45^{}`$. (We present only some of the more representative results; angles between 0 and $`45^{}`$ were studied in $`5^{}`$ increments.)

Comparison of Fig. 1 with Fig. 2 for $`\theta =5^{}`$ shows that a small misalignment of the applied field can change the hysteresis loop drastically.

For $`N=3`$, Fig.2 shows that the central part of the hysteresis loop shrinks as $`\theta `$ increases. For $`\theta =45^{}`$, the central part almost disappears completely, and new small loops start to develop away from the center of the hysteresis loop. Our results for $`N=3`$ agree with those given in CS1 , which studied the cases $`\theta =0^{}`$ and $`45^{}`$.

For $`N=4`$, at $`\theta =0^{}`$ there is no central loop, but there is a prominent loop at finite field. On the other hand, at $`\theta =5^{}`$, there is a central loop, and the finite field structure becomes rather complex. Further increase of $`\theta `$ leads to a filling out and connecting of various subloops. Also, note the appearance of jumps for non-zero $`\theta `$.

For $`N=5`$, a small field misalignment has an enormous effect, at $`\theta =5^{}`$ shrinking the loop to a relatively small central region. As $`\theta `$ increases, the central loop grows, but the loop for $`\theta =45^{}`$ pinches off to yield three subloops, as for $`N=3`$.

For $`N=6`$, again a small field misalignment has an enormous effect, at $`\theta =5^{}`$ shrinking the loop to a relatively small central region. As $`\theta `$ increases, the central loop grows, but in contrast to $`N=5`$, the loop for $`\theta =45^{}`$ does not pinch off, and closely resembles the loop for $`N=4`$.

These different types of behavior indicates that these are complex systems, for which it is difficult to make generalizations.

## IV Hysteresis and Even-Odd Signature in Finite Array of Nano-Dots

In the absence of an external magnetic field the array of $`N\times N`$ nano-dots favors antiferromagnetic ordering, thus minimizing its magnetostatic energy. A large external magnetic field applied to the array tends to orient the magnetic moments along the field, thus minimizing the Zeeman field energy. However, the spins at the array corners tip by a small angle, as shown in Fig. 3-Fig. 6, forming a two-dimensional “flower“ state aharoni , TLW . The flower state persists until the applied field falls to $`H_0=M_s`$. For lower values of the applied field, the competition between the dipole-dipole interaction and the Zeeman energy becomes significant, and changes the array ordering.

For $`N=2`$, Fig. 3 shows that below $`H_0=M_s`$ and for $`H_00`$ the spins form a snake-like domain structure winding clockwise or counterclockwise. At $`H_0=0`$, the array has zero net per-dot magnetization, due to a vortex-like structure that persists for $`3M_sH_00.3M_s`$.

The $`N=3`$ array was analysed by Camley and Stamps in CS1 . This array also shows snake-like arrangements below $`H_0=M_s`$. At zero field the final state of the array shows what we call a ”barrel” state in which the spins at the left and right columns are oriented opposite to the central column with a slight tipping of the corner spins, as shown in Fig. 4. This agrees with Ref. CS1,  except for the tipping of the corner spins. We have made a small angle expansion of the energy for the spins being nearly aligned, and indeed we find that the tipped corner spin state is more energy favorable. The tipping angle was determined by both iteration and analytical calculations, and agreement is found between the two values. The numerical value of the corner spin tipping angle, found both by iteration and small tipping angle analysis, is $`|\alpha |=9.1115^{}`$, with numerical error of order $`10^5`$. The spin snapshots in Fig. 4 and hysteresis loop analysis show that the barrel state switches to an inverted barrel state when the applied field changes sign.

Fig.5 shows that the $`N=4`$ array also features snake-like arrangements of the spins, for intermediate values of the applied field. However, in zero field the total per-dot magnetization is zero, which can be attributed to the formation of a vortex in the array’s central $`2\times 2`$ block. The magnetic moments of the rest of the dots in the array form a ring that surrounds the vortex with opposite circulation. This state is stable for applied fields $`H_0`$ satisfying $`0.2M_sH_00.2M_s`$.

For all $`N`$, the flower state appear at high fields (here, $`|H_0|=2.0M_s`$). The hysteresis loops shown in Fig. 1 show a subtle difference in shape between arrays with odd $`N`$ and arrays with even $`N`$. For odd $`N`$ the loops show well-defined jumps whereas for even $`N`$ this behavior is absent. This behavior is due to unpaired spins with uncompensated dipole fields. The jumps become less apparent as $`N`$ grows, and will eventually disappear for large $`N`$, where the distinction between even and odd $`N`$ becomes unimportant.

Note that, when an array was placed in zero external field and given random initial conditions, the solutions converged to the same states as obtained in the hysteresis-cycle calculations, up to the degeneracy of the system. Thus, for $`N=2,4`$ there are two degenerate metastable states of opposite chirality (winding) with zero net magnetization in zero field, each of which has a four-fold rotational symmetry. For $`N=6`$ there are two degenerate states of opposite chirality, with non-vanishing net magnetization, each of these state has no apparent rotational symmetry. For $`N=3`$ there are two degenerate barrel states, with no rotational symmetry, and $`N=5`$ is similar to $`N=6`$. The hysteresis loop area $`A_N`$ will produce further evidence that large system behavior commences with $`N=5`$ and $`N=6`$. Metastable states with vanishing net magnetization may appear for arrays with even $`N`$. However, our simulations showed that these states could appear only for $`N=2,4`$. For arrays with odd $`N`$ the unpaired dipoles prevent the formation of such states.

## V Hysteresis Loop Area $`A_N`$ vs. Particle Number $`N`$

Although the area of the hysteresis loop $`A_N`$ tends to zero for the $`N=2`$ array, it clearly is nonzero for all other arrays. We present the hysteresis loop areas in Fig. 8 as circles and squares. From Fig. 8, we see that the area of the hysteresis loop decreases with increasing $`N`$ except for $`N=3`$ for $`N`$ odd and $`N=2,4`$ for $`N`$ even. The $`N=5`$ and $`N=6`$ arrays, the first to show something like bulk behavior, have maximum $`A_N`$ for odd and even respectively; their spin arrangements are given in Fig. 6 and Fig. 7. We have fitted our data to the asymptotic form

$`A_N=A_{\mathrm{}}+{\displaystyle \frac{C}{N^p}}`$ (5)
where $`A_{\mathrm{}}`$, $`C`$ and $`p`$ are constants to be determined. If larger values of $`N`$ had been computationally feasible, we would have considered only large values of $`N`$ for the fit. The fit is not as good when $`N=5`$ is included, so we do not show this case. In practice, for $`N`$ odd the data are fit starting from $`N=7`$ and for $`N`$ even from $`N=6`$. Both fits are shown as solid lines in Fig. 8, where the values of $`A_{\mathrm{}}`$, $`C`$, and $`p`$ are given in the table.

The associated values of $`\chi ^2`$ are both less than $`10^5`$.

For odd $`N`$, $`A_N`$ varies approximately as $`A_NN^{\frac{3}{2}}`$ whereas for even $`N`$ it varies approximately as $`N^2`$. We attribute no fundamental significance to these values.

## VI Summary

We have studied the hysteresis and magnetization processes for $`N\times N`$ arrays (with $`N=2\mathrm{}13`$) of uniaxial ferromagnetic nano dots interacting via the dipole-dipole interaction. For an external magnetic field aligned or misaligned with one side of the array, the hysteresis loops are surprisingly complex. For arrays with odd $`N`$ the hysteresis loops possess jumps, whereas for even $`N`$ they do not. As $`N`$ increases, the area $`A_N`$ of the hysteresis loop begins to saturate, approaching a non-zero finite value determined from a data fit. The area of the hysteresis loop scales with $`N`$ approximately as $`N^{\frac{3}{2}}`$ for $`N`$ odd, and approximately as $`N^2`$ for even.

We would like to thank V. L. Pokrovsky, A. S. Kirakosyan and S. Erdin for fruitful discussions. This work was supported by NSF grants DMR 0103455 and DMR 0072115, DOE grant DE-FG03-96ER45598, and Telecommunication and Information Task Force at Texas A&M University.


text:
# A more detailed look at the Opacities for Enriched Carbon and Oxygen Mixtures.

## 1 Introduction.

The main source of a star’s energy is the nuclear fusion reactions occurring either in its core or in thin burning shells around the core. This energy is transported from the production site to the surface by radiative transfer or convection. The first transports energy in the form of photons, while the second is a cyclic macroscopic mass motion that carries the energy in bulk.

In regions where convection is stabilised, radiative transfer leads to the equation of stellar structure,

$$\frac{dT}{dr}=\frac{3\overline{\kappa _\mathrm{R}}\rho }{16\sigma T^3}\frac{L_r}{4\pi r^2},$$
(1)
where $`L_r`$ is the luminosity at radius $`r`$, $`\rho `$ the density, $`\sigma `$ the Stefan-Boltzmann constant, $`T`$ the temperature and $`\overline{\kappa _\mathrm{R}}`$ is the Rosseland mean opacity.

Opacity $`\kappa `$ is a measure of the degree to which matter absorbs photons. There are four main sources of opacity in stars

* Bound-bound transitions are the transitions of an electron in an atom, ion or molecule between energy levels which are accompanied by either the absorption or emission of a photon. Only a photon with the correct wavelength can cause a given transition, the process is wavelength dependent.
* Bound-free transitions, or photoionisation, occur when an incoming photon has enough energy to ionise an atom or ion and free an electron. The reverse process is the capture of an electron by an atom or ion. It will not occur until photons above a threshold energy are available and falls off as $`\nu ^3`$ where $`\nu `$ is the frequency.
* Free-free transitions are scattering processes which occur when an electron and photon interact near an atom or ion. The process is also known as bremsstrahlung. Again it is proportional to $`\nu ^3`$.
* Electron scattering is wavelength independent at low temperatures where it is Thomson scattering. The electron, with its low cross-section, is a small target and so only dominates at high temperatures when most atoms are ionised. At very high temperatures relativistic effects are important and Compton scattering dominates.

From this list it is possible to see that calculating the opacity in a stellar model is a difficult process and depends on the composition, temperature and density of the material. Deep in the star local thermodynamic equilibrium (LTE) is achieved and an average for radiative transfer over all wavelengths requires the Rosseland mean opacity $`\overline{\kappa _\mathrm{R}}`$ (Rosseland, 1924) expressed as

$$\frac{1}{\overline{\kappa _\mathrm{R}}}=\frac{{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{1}{\kappa _\nu }}{\displaystyle \frac{B_\nu }{T}}𝑑\nu }{{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{B_\nu }{T}}𝑑\nu },$$
(2)
where $`\kappa _\nu `$ is the opacity at frequency $`\nu `$, $`T`$ is the temperature and $`B_\nu `$ is the flux per unit area into unit solid angle per unit frequency in LTE.

The most comprehensive opacities available today are from the OPAL (Iglesias & Rogers, 1996) or OP (Cunto et al, 1993) groups who have made detailed models of the above processes. They provide tables of the Rosseland mean opacity variation with temperature, density and composition with the metal abundance usually scaled to solar compositions. Iglesias & Rogers (1993) took a step forward for the OPAL project team by providing tables that include mixtures enhanced in carbon and oxygen (C and O) relative to the base solar composition. We have incorporated their full range of composition tables in the Eggleton (1971) evolution code. In the last implementation of this code (Pols et al., 1995) there only 10 of the 265 tables available were used. This provides us with the opportunity to refine our models so as to accurately follow the changes of opacity which occur in the later stages of evolution. The effects are expected in many types of stars. While these may be small for the evolution of main-sequence stars and white dwarfs we expect larger differences to be found in AGB and Wolf-Rayet stars.

* Asymptotic Giant Branch (AGB) stars undergo third dredge-up which mixes helium burning products to the surface and forms carbon stars. Using the enhanced mixture tables we shall be able to model the thermal pulses and envelope evolution more accurately.
* Wolf-Rayet (WR) stars are massive and have lost their hydrogen envelopes exposing the helium cores. As time progresses helium burning products are slowly exposed at the surface and in some cases the stars are eventually mostly composed of carbon and oxygen.

We present our method of opacity table construction and detail its implementation in the Eggleton stellar evolution code. We discuss the effects on main-sequence stars, red giants and white dwarfs. We then present three tests of our work. The first is the effect of including extra carbon and oxygen on the structure and evolution of a low-mass population-III star, the second is a $`5\mathrm{M}_{}`$ thermally pulsing AGB star and the third a Wolf-Rayet star of $`40\mathrm{M}_{}`$ with mass loss.

## 2 Table Construction.

We start with the OPAL tables (Iglesias & Rogers, 1996) as the framework around which we construct our full tables in an approach similar to that employed by Pols et al. (1995). However we choose to use the variable $`=\rho /T_6^3`$ rather than the density. This has the advantage that our tables can be a third smaller in this dimension. This partly compensates for increased memory requirement of the full tables.

The OPAL tables extend in $`\mathrm{log}_{10}(T/\mathrm{K})`$ from $`3.75`$ to $`8.70`$ in steps of $`0.05`$ to $`0.2`$ and in $`\mathrm{log}_{10}(/\mathrm{g}\mathrm{cm}^3\mathrm{K}^3)`$ from $`7`$ to $`+1`$ in steps of $`0.5`$. We choose the range of the final tables to cover from $`3`$ to $`10`$ in $`\mathrm{log}_{10}(T/\mathrm{K})`$ with increments of $`0.05`$, and from $`8`$ to $`+7`$ in $`\mathrm{log}_{10}(/\mathrm{g}\mathrm{cm}^3\mathrm{K}^3)`$ with increments of $`0.5`$. This gives us tables of 141 by 31 elements.

To account for the effect of composition, the OPAL tables include typically 265 $`T`$ and $``$ grids for each metallicity. These tables are split into 5 groups with different hydrogen mass fractions X of $`0.0`$, $`0.03`$, $`0.1`$, $`0.35`$ and $`0.7`$. The tables in each group have different compositions in helium, carbon and oxygen.

The lowest temperature in the OPAL tables is $`\mathrm{log}_{10}(T/\mathrm{K})=3.75`$ and we use the tables of Alexander & Ferguson (1994), which extend from $`3`$ to $`4.1`$ in $`\mathrm{log}_{10}(T/\mathrm{K})`$, to complete them at low temperature. Where the tables overlap they match well (Iglesias & Rogers, 1996). It should be noted that the tables provided by Alexander & Ferguson (1994) do not include enhanced carbon and oxygen mixtures. Therefore at temperatures when $`\mathrm{log}_{10}(T/\mathrm{K})<4.0`$ our tables do not follow these mixtures. This is only important in the surface of AGB stars. We plan in future to deal with low temperature enhanced mixtures by including the work of Marigo (2002) before application to the structure of AGB star envelopes.

To complete each table we fill the region from $`8.70`$ to $`10`$ in $`\mathrm{log}_{10}(T/\mathrm{K})`$, where electron scattering dominates, according to Buchler & Yueh (1976). We then combine with a full table of the effective opacity owing to electron conduction for which we use the fits of Iben (1975) to the tables of Hubbard & Lamp (1969) and Canuto (1970). We combine the radiative opacity with the effective conductive opacity by reciprocal addition,

$$\frac{1}{\kappa _{\mathrm{eff}}}=\frac{1}{\kappa _{\mathrm{rad}}}+\frac{1}{\kappa _{\mathrm{cond}}}.$$
(3)
Before the complete tables are implemented in the stellar evolution code we perform one last process. We rescale the tables over the C and O plane to a more regular grid. In the OPAL grid the spacing used requires interpolation within irregular polygons, so to make the interpolation simpler and faster, we include extra tables in a similar but fully rectangular grid system. This means we have 61 tables for each X abundance or 305 tables for each metallicity.

## 3 Implementation into a Stellar Evolution Code.

We implement the tables in the stellar evolution code first described by Eggleton (1971) and most recently updated by Pols et al. (1995). Further details can be found in the references therein. The main drawback with implementation of the new tables is the increased memory requirement. Whereas before there were only 10 opacity tables there are now 305. The new tables are each a third of the size of the old tables so the overall increase is a factor of ten. This corresponds to a memory requirement increase from 17 Mb to 170 Mb but the average specification of today’s desktop computers comfortably accommodates this. One of the reasons why this memory requirement is so large is that we also store the spline coefficients for each of the tables. When we use the tables within the code we interpolate in the $``$ and $`T`$ plane with bicubic splines. This ensures the smoothness of our tables. We then use linear interpolation in the 3 composition dimensions to derive the final opacity value at each point.

To test the effect of these tables we use three different interpolation schemes,

* Method A uses 5 tables with $`X=0,\mathrm{\hspace{0.17em}0.03},\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}0.35}`$ and $`0.7`$ and $`Y=1XZ`$. Where $`X`$ is the hydrogen mass fraction, $`Y`$ the helium mass fraction and $`Z`$ the initial, solar mixture, metallicity. We add a further two tables, both with $`X=Y=0`$, the first with $`X_\mathrm{C}=1Z`$ and $`X_\mathrm{O}=0`$. The other with $`X_\mathrm{C}=0`$ and $`X_\mathrm{O}=1Z`$. Where $`X_\mathrm{C}`$ and $`X_\mathrm{O}`$ are the mass fractions for the enhanced amount of carbon and oxygen above that included in the metallicity mass fraction so the total carbon mass fraction is $`X_\mathrm{C}+Z(\mathrm{C})`$ and for oxygen $`X_\mathrm{O}+Z(\mathrm{O})`$. This gives a total of 7 tables allowing simple interpolation in composition similar to the old method.
* Method B is similar to method A but we include two extra tables for each hydrogen mass fraction. Again both have $`Y=0`$. One has $`X_\mathrm{C}=1XZ`$ and $`X_\mathrm{O}=0`$ and the other $`X_\mathrm{C}=0`$ and $`X_\mathrm{O}=1XZ`$. This leads to a total of 15 tables in all.
* Method C incorporates all 305 tables in the opacity calculations. We have 61 tables for each hydrogen abundance. When interpolating in composition we use three variables $`X`$, $`X_\mathrm{C}/(1XZ)`$ and $`X_\mathrm{O}/(1XZ)`$. The hydrogen abundances are as for method A. We take values of $`0,\mathrm{\hspace{0.17em}0.01},\mathrm{\hspace{0.17em}0.03},\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}0.2},\mathrm{\hspace{0.17em}0.4},\mathrm{\hspace{0.17em}0.6}`$ and $`1`$ for the carbon and oxygen grid planes. On this grid we perform a three-dimensional linear interpolation.

## 4 Testing & Results.

The new tables give almost identical results to the old for main-sequence stars and red giants. The timescale of core helium burning is altered. From method A to method C it increases by a factor of $`0.1`$ percent for a $`5\mathrm{M}_{}`$ star. This is reassuring. Current stellar evolution models describe well the evolution of the main sequence and the red giant branch. However we do find that the opacity tables enhance numerical stability in the later stages of evolution because they resolve changes with composition in greater detail. This leads to a smoother variation in opacity.

A similar result is found for white dwarfs. Between methods A and C for a carbon and oxygen white dwarf the radius increases by 0.2 percent. While conditions in the atmosphere do pass through regions of the greatest difference between the two methods (the dotted line in figure 3) there is not enough mass at these temperatures and densities to make a significant global difference.

### 4.1 Polluting a $`0.5\mathrm{M}_{}`$ Population-III Star.

For our first detailed test we take a zero metallicity $`0.5\mathrm{M}_{}`$ star and pollute it uniformly with carbon and oxygen. Such a star could have been formed just after the death of the first stars if the products of these stars were mainly carbon and oxygen. We choose such a low-mass star to limit the effect on the results of burning via the CNO cycle. Our initial model has the following properties, $`X=0.7`$, $`Y=0.3`$, $`R=0.45\mathrm{R}_{}`$, $`T_{\mathrm{eff}}=4780\mathrm{K}`$ and $`L=0.96\mathrm{L}_{}`$. The tables provided by OPAL do not include one with $`X=0.75`$, we could extrapolate to obtain the opacity but for this test we prefer to remain within the tables. From the initial model we perform three sub-tests. In the first we pollute the star at a constant C/O ratio of 1, then evolve each model and compare the time from zero-age main sequence to the helium flash for methods A, B and C. We also compare a sequence of models for stars of the same mass but with a metallicity of the same mass fraction as the carbon and oxygen. These results are presented in table 1. We see that with the pollution fraction increase, the timescale of hydrogen burning increases because the energy produced at the centre takes longer to reach the surface. However when comparing the difference between methods B & C we see only a comparatively small increase in timescale. When we compare with the scaled solar models we see a greater difference. This is due to the much larger opacity of heavier elements such as iron.

With the second set of models we look at the radius, surface temperature and luminosity of the star on the main sequence as we pollute it with either only carbon or only oxygen (tables 2 and 3). We use method A as the base to which we compare methods B and C. A similar trend is seen in the results in tables 2 and 3. We therefore see it is important to include the affect of enhanced C and O mixtures before hydrogen exhaustion. Similar mixtures are encountered in the envelopes of carbon stars of very low metallicity.

Finally we consider the effect of varying the C/O ratio of the pollution (table 4) at a constant pollutant mass fraction of $`0.001`$. Table 4 demonstrates the need to use method C to follow accurately the variation in opacity as the C/O ratio changes. This is important in AGB evolution during the thermal pulses when the C/O ratio changes gradually as carbon becomes dominant. The difference between methods B and C in this case comes about because of non-linear structure in the C and O opacities that can only be resolved with the full set of tables.

### 4.2 Thermal Pulses of an AGB Star.

Our thermally pulsing AGB model is of a $`5\mathrm{M}_{}`$ star. To follow a number of pulses we include convective overshooting (Schröder, Pols & Eggleton, 1997) during the pre-AGB evolution. This leaves the star with a larger core which reduces the rate at which the thermal pulses grow in strength. This means we do not find dredge-up in our models but we can evolve through the pulses. We obtain the average pulse period as the arithmetic mean of the times between the first ten helium shell flashes. Results are recorded in table 5.

First we find that the different methods affect the age of the star when it starts to undergo thermal pulses. Second we find that methods A and B have similar interpulse periods while method C’s is $`3.6`$ percent less. We attribute this to the change in opacity due to varying He/C/O ratios in the inter-shell region which alters the time evolution of the pulse.

### 4.3 Evolution of a Wolf-Rayet Star.

We have used the mass-loss rates of Nugis & Lamers (2000) to manufacture a Wolf-Rayet star of solar metallicity and an initial mass of $`40\mathrm{M}_{}`$. Some details are described by Dray & Tout (2003). Our tables are expected to make a difference because Wolf-Rayet stars become naked helium stars and helium burning products are mixed throughout and exposed at the surface of the star. We only use methods A and C and not method B because the carbon and oxygen mass fractions do not rise above those included in the base metallicity until after the hydrogen envelope has been removed and so our tables only begin to make a difference after this.

In figure 1 we present the variation of opacity through the star, just after core carbon ignition, with mesh point, the grid over which we solve the equations of stellar structure. The solid line is with method A. We then take this structure and use method C to calculate the opacities to give the dashed line. This model is not in hydrostatic equilibrium and the dot-dashed line shows the relaxed new structure in hydrostatic equilibrium from method C. Comparing the physical parameters of these calculations we find the main change to be an increase in radius of 3.0 percent. The surface temperature has decreased by 1.5 percent.

We then made two complete runs with similar initial conditions by methods A and C. In figure 2 we display the difference between Wolf-Rayet evolution as the hydrogen envelope is removed. The greatest changes are in the radius and surface temperature. The star is in the WN phase between $`14`$ to $`12\mathrm{M}_{}`$, after which it enters the WC phase. In the WN phase there are only small differences but upon entering the WC phase the differences build up so that the temperature is lower by around 2 percent and the radius about 1 percent larger. The luminosity also tends to be larger by around 0.5 percent. Figure 3 is a plot of where the difference between methods A and C are greatest in $`\rho `$ and $`T`$ when $`X=0`$. The plot is of the mean difference in the opacity calculated by methods A and C at all possible carbon and oxygen abundances divided by the opacity when $`Y=1Z`$. From the plot we see that the regions of greatest difference are small. The lines show where the Wolf-Rayet model calculated by method C lies at certain points of its evolution. The solid line is for the star at $`12.3\mathrm{M}_{}`$ and the dashes at $`6.1\mathrm{M}_{}`$. From figure 2 we see that the largest changes begin around the $`12\mathrm{M}_{}`$ point. This is when the model has moved into the region where the differences between A and C are greatest.

We find the opacity tables alter the timescales for the late burning stages. Method C increases the helium burning timescale by 0.1 percent while the carbon-burning time scale increases by 2.5 percent relative to method A. These changes come about because the opacity varies with composition and effects the temperature structure and burning rates. While there is no large difference in a single variable the overall effect on structure and evolution is more significant.

## 5 Conclusions.

Our major conclusion is that the changes induced by properly including opacities for varying C and O mixtures are small. However they are a good thing to include in stellar evolution codes because they help numerical stability by removing sharp changes in the variation of opacity with composition that are encountered when we use fewer tables and add another level of detail to models without much loss of computational speed.

The main computational cost is the extra memory to store the tables. Method C requires ten times the memory (about 200 Mb) compared to methods A and B. However there is little computational speed cost: relative to method A, methods B and C require about 4 percent more time for the evolution of a Wolf-Rayet star even though we must evaluate four times as many opacities for each $``$ and $`T`$ value than with method A.

The tables and stellar evolution code are freely available from http://www.ast.cam.ac.uk/$``$stars. Interested users are welcome to contact the authors for details on how to download and implement our tables.

## 6 Acknowledgements.

The authors would like to thank Onno Pols for informative suggestions. JJE would like to thank PPARC for his funding and his parents for helping to provide the computing power for this work. CAT thanks Churchill College for his fellowship.


text:
# A More Improved Lattice Action for Heavy Quarks

## 1 INTRODUCTION

One source of uncertainty in lattice QCD comes from the discretization effects in numerical simulations. Here we extend the Fermilab formalism to higher dimensional operators to obtain a more improved action for heavy quarks, aiming to reduce this uncertainty to 1–2%.

The method starts with Symanzik , who introduced a local effective Lagrangian

$$_{\mathrm{eff}}=_{\mathrm{QCD}}+_ja^{dim_j4}K_j_j$$
(1)
of a continuum field theory. The $`_j`$ are higher-dimensional operators. The lattice action is adjusted so that the short-distance $`K_j`$ coefficients vanish, at least approximately.

The Fermilab formalism  takes Wilson fermions  with the clover action  and introduces different couplings for space-like and time-like interactions. Ref.  considers operators up to dimension five. Ref.  adds dimension-six operators to this action. In this paper, we give a more complete action which has dimension six and some dimension seven operators and determine the couplings by matching at tree level.

## 2 QUARK BILINEARS

In this section, we consider the quark bilinears. The dimension-five (-six) interactions are listed in Ref.  (Ref. ). Redundant directions are identified by field transformations. They are used to eliminate the higher-order time derivatives, which would give large errors for large masses. We do not repeat this analysis here.

The corrections in Eq. (1) should be classified by a power-counting scheme. For light quark physics the only relevant ratio is $`\mathrm{\Lambda }a`$. For heavy quark physics, however, the power counting should follow HQET for heavy-light systems or NRQCD for quarkonium. In HQET, operators are assigned a power of $`\lambda =\mathrm{\Lambda }/m_Q`$; in NRQCD, a power of relative velocity $`v`$. In these power-counting schemes, some dimension-seven interactions are commensurate with dimension-six interactions. As seen in the list in Table 1, they scale as $`\lambda ^3`$ or $`v^6`$.

Further dimension-seven operators carry power $`\lambda ^4`$ (or higher) in HQET and $`v^8`$ (or higher) in NRQCD.

Including the needed dimension-seven operators, our lattice action takes the form

$`S_F`$ $`=`$ $`S_0+S_B+S_E`$ (2)
$`+`$ $`a^2c_1{\displaystyle \overline{\psi }\{𝜸𝑫,𝑫^2\}\psi }+a^2c_2{\displaystyle \overline{\psi }\gamma _iD_i^3\psi }`$
$`+`$ $`a^3r_3{\displaystyle \overline{\psi }(𝑫^2)^2\psi }+a^3r_4{\displaystyle \overline{\psi }D_i^4\psi }`$
$`+`$ $`a^2c_5{\displaystyle \overline{\psi }\{𝜸𝑫,i𝚺𝑩\}\psi }`$
$`+`$ $`a^2c_6{\displaystyle \overline{\psi }(𝑫𝑬𝑬𝑫)\psi }`$
$`+`$ $`a^2c_7{\displaystyle \overline{\psi }𝜸(𝑫\times 𝑩+𝑩\times 𝑫)\psi }`$
$`+`$ $`a^2c_8{\displaystyle \overline{\psi }\{𝜸𝑫,𝜶𝑬\}\psi }`$
$`+`$ $`a^3r_9{\displaystyle \overline{\psi }\{𝑫^2,i𝚺𝑩\}\psi }`$
$`+`$ $`a^3r_{10}{\displaystyle \overline{\psi }_{ij}\{i\mathrm{\Sigma }_iB_i,D_j^2\}\psi }`$
$`+`$ $`a^3r_{11}{\displaystyle \overline{\psi }_{ij}i\mathrm{\Sigma }_iD_jB_iD_j\psi }`$
$`+`$ $`a^3r_{12}{\displaystyle \overline{\psi }𝜸𝑫i𝚺𝑩𝜸𝑫\psi },`$
where conventional difference operators and clover fields are understood. We have included some terms that we expect to be redundant, to retain some flexibility to, say, speed up computer simulations.

For brevity we do not introduce terms, corresponding to the last six rows of Table 1. The operators $`\overline{Q}D_ii𝚺𝑩D_iQ`$ and $`\overline{Q}𝑫(𝑩\times 𝑫)Q`$ are redundant. Matching of the last four rows in Table 1 requires a calculation of Compton scattering, which is is still in progress.

## 3 MATCHING

We derive matching conditions for the improvement couplings in Eq. (2) by matching on-shell quantities for small momentum but any $`m_Qa`$. We do this by calculating the energy momentum relation and the scattering of a heavy quark off a background field. The energy momentum relation has been studied  and summarized  before. It gives us conditions on $`c_1`$, $`c_2`$, $`r_3`$ and $`r_4`$. Results for the chromoelectric terms $`c_E`$, $`c_6`$ and $`c_8`$ are given in Ref. , determined from the temporal component of the background field. Scattering from the spatial component of the background field yields results for the chromomagnetic terms. Our results yield the following matching conditions. The energy momentum relation gives:

$`16\zeta c_1={\displaystyle \frac{4\zeta ^4(1\zeta ^2)}{[m_0(2+m_0)]^2}}{\displaystyle \frac{3r_s^2\zeta ^4}{(1+m_0)^2}}`$ (3)
$`+`$ $`{\displaystyle \frac{\zeta ^3[2\zeta +4r_s(1+m_0)6r_s\zeta ^2/(1+m_0)]}{m_0(2+m_0)}}`$
$`+`$ $`{\displaystyle \frac{m_0(2+m_0)}{2(1+m_0)}}[{\displaystyle \frac{r_s^2\zeta ^2}{(1+m_0)}}{\displaystyle \frac{r_s^3\zeta ^3}{(1+m_0)^2}}8r_3].`$
where $`r_3`$ can be chosen to be zero. Also

$$0=\zeta ^2+6\zeta c_2+(r_s\zeta 24r_4)\frac{m_0(2+m_0)}{8(1+m_0)}.$$
(4)
We do not set $`r_4=0`$ or $`r_s\zeta /24`$ but choose it below so that $`r_{11}=0`$. For the chromoelectric interactions one finds

$`\zeta ^2c_E^2+c_8{\displaystyle \frac{2m_0(2+m_0)}{(1+m_0)}}`$ $`=`$ $`{\displaystyle \frac{\zeta ^2(\zeta ^21)}{m_0(2+m_0)}}+{\displaystyle \frac{r_s\zeta }{(1+m_0)}}`$ (5)
$`+`$ $`{\displaystyle \frac{r_s^2m_0(2+m_0)}{4(1+m_0)^2}},`$
$$c_6=0.$$
(6)
For the chromomagnetic interactions one finds

$`c_B`$ $`=`$ $`r_s,c_7=0,c_5=c_1,r_9=r_3,`$ (7)
$`r_{10}`$ $`=`$ $`{\displaystyle \frac{1}{48}}(r_s\zeta 24r_4)+{\displaystyle \frac{1}{12}}c_B\zeta ,`$ (8)
$`r_{11}`$ $`=`$ $`{\displaystyle \frac{1}{48}}(r_s\zeta 24r_4){\displaystyle \frac{1}{6}}c_B\zeta ,`$ (9)
$`r_{12}`$ $`=`$ $`0.`$ (10)
In summary, we may set $`c_6=c_7=c_8=r_3=r_9=r_{11}=r_{12}=0`$, leaving only five non-zero couplings: $`c_1`$, $`c_2`$, $`r_4`$, $`c_5`$ and $`r_{10}`$.

## 4 ERROR ESTIMATES

We now estimate the discretization effects. It is advantageous to use a heavy-quark description of cutoff effects , because then contributions from, say, $`\overline{Q}𝜶𝑬Q`$ and $`\overline{Q}\{𝜸𝑫,𝜶𝑬\}Q`$ are automatically combined. Thus, we express the uncertainties as

$$\mathrm{𝙴𝚛𝚛𝚘𝚛}_i=\left(C_i^{\text{lat}}C_i^{\text{cont}}\right)O_i$$
(11)
where $`C_i`$ are HQET/NRQCD short-distance coefficients for lattice gauge theory and continuum QCD . The mismatched coefficients can be obtained from the results of the previous section.

The matrix elements $`O_i`$ are estimated by using HQET (NRQCD) power counting for heavy-light systems (quarkonia). Since higher-order calculations are not yet available we multiply the tree-level mismatch with $`\alpha _s^l`$ to estimate the $`l`$loop mismatch. This is a conservative estimate because the most pessimistic asymptotic behavior for the coefficients of the improvement operators is the same at higher orders as in the tree level formulas of previous section. To obtain the numbers shown in Figure 1, we take $`\alpha _s=0.25`$.

Fig. 1 shows results for the chromomagnetic energy (labelled $`1/2m_B`$), chromoelectric energy (labelled $`1/4m_E^2`$), and the relativistic correction term (labelled $`1/8m_4^3`$). See Refs.  for details.

The vertical lines represent $`a=1/8`$ fm and $`a=1/11`$ fm, corresponding to the unquenched MILC lattices . The black (gray) curves represent the charmed (bottom) quark. Note that, in most cases, bottom quarks have smaller discretization effects than charmed quarks. Thus, Fig. 1 shows concretely that the Fermilab formalism avoids errors of order $`(m_Qa)^n`$.

## 5 CONCLUSIONS

We have presented a more improved action for heavy quarks with five new operators. Our error analysis shows that, to reduce errors below 1%, one-loop matching is essential for the coupling $`c_B`$ of the chromomagnetic clover term. Tree-level matching for the coupling of the chromoelectric term, $`c_E`$, is usually enough, but one-loop may be necessary for charmonium. Tree-level matching suffices for the other couplings in Eq. (2).

M.B.O. and A.X.K. are supported in part by the U.S. Department of Energy under contract DE-FG02-91ER40677. Fermilab is operated by the Universities Research Association Inc., under contract with the U.S. Department of Energy.


text:
# A search for LSB dwarf galaxies in various environments

## 1. Introduction

According to standard Cold Dark Matter (CDM) hierarchical clustering theory, there should be numerous low mass dark matter halos present in the Universe today. If these halos contain sufficient stars, they should be detectable as dwarf galaxies. Observationally this appears to be true for clusters of galaxies where the galactic density is high, but not so for the lower density environments. We conducted a search for these objects in the Millennium galaxy strip which runs along the celestial equator in the field, passing through filaments and voids. It is therefore an excellent data set for studies into the influence of the environment on dwarf galaxy populations. We compare these results with those from similar surveys carried out in the Virgo and Ursa Major (UMa) clusters. Our results are unique as the three surveys were conducted using the same instrument, same technique (exposure time, filter band) and same selection criteria, thus we can be sure that we are comparing ’like with like’.

Low surface brightness (LSB) galaxies are difficult to detect as their surface brightnesses are below that of the sky ($`23mag/arcsec^2`$). The detection algorithm that we developed for this project is optimised for the detection of faint, diffuse objects on CCD frames (see Sabatini et al. 2003). To ensure that the objects picked out by the algorithm are actually dwarf galaxies and not background contamination, a selection criteria based on morphology and magnitude is applied to the objects. This criteria was chosen following simulations of a cone of the universe randomly populated with galaxies, as detailed in Roberts et al (2003).

## 2. Results and Discussion

We have presented the results obtained for 3 surveys carried out in very different environments (table 1). We find a DGR in the field of 6:1, compared to a value of 20:1 in the Virgo cluster. This very large ratio of dwarf to giant galaxies found in the Virgo cluster indicates that this region is very different to lower density clusters such as UMa, and the field where we find relatively few dwarfs.

Our results for the DGR of the MGS are consistent with those derived from the recent redshift survey determinations of the field LF made by 2dF (Norberg et al. 2002) even though we sample to some two magnitudes fainter in central surface brightness and magnitude. There is no hidden population of dwarf galaxies that have been missed by the redshift surveys.

These observational results are in disagreement with most predictions made by CDM models, commonly referred to as the ’substructure’ problem. The models predict far more small dark matter halos than observations detect (Kauffman et al. 1993). A number of theories have been put forward to explain the apparent difference in the observed number of dwarf galaxies in different environments. These range from those which try and explain how more dwarfs may form in rich environments, such as galaxy harassment (Moore et al 1996) or the external pressure confining the ejected galaxy gas (Babul & Rees 1992), to those ideas which emphasise the suppression of dwarf galaxy formation in the field such as supernovae wind expulsion and galaxy squelching (Tully et al 2002). Detailed observations of dwarf galaxies provide a challenge to the ’concordance’ cosmological model to which CDM is central. Dwarf galaxies are found in large numbers in rich clusters, but not in less dense galactic environments. For the CDM model to remain viable it has to provide a satisfactory solution to this problem. At present, it is not clear which, if any, of the mechanisms described above would be the best to help provide this solution.

## References

Babul, A., & Rees, M.J., 1992, MNRAS, 255, 346B

Kauffmann, G., White, S.D.M., Guideroni, B., 1993, MNRAS, 264, 201

Moore, B., Lake, G., Quinn, T., Stadel,J., 1999, MNRAS, 304, 465M

Norberg et al. , 2002, MNRAS, 336, 907

Roberts, S., Davies, J., Sabatini, S., Baes, M., Linder, S., Smith, R., Evans, Rh., 2003, submitted

Sabatini, S., Davies, J., Scaramella, R., Smith, R., Baes, M., Linder, S.M., Roberts, S., Testa, V., 2003, MNRAS, 342, 981

Tully, R.B., Somerville, R.S., Trentham, N., Verheijen, M.A.W., 2002, ApJ, 569, 573


text:
# Pulsar Birthrate from Parkes Multi-beam Survey

## 1. Introduction

The Parkes multi-beam pulsar survey covers a $`10^{}`$-wide strip of the southern Galactic plane from $`l=260^{}`$ to $`l=50^{}`$ (Manchester et al. 2001). It utilises a 13-beam receiver operating in the 20-cm band on the Parkes 64-m radio telescope and is seven times more sensitive than any previous large-scale survey. The survey has discovered more than 600 new pulsars so far, including many that are young and energetic. Already published pulsars can be found in the on-line catalogue at http://www.atnf.csiro.au/research/pulsar/psrcat. In the following analysis we used 520 pulsars that were discovered by the Parkes survey, of which 415 are obtained from public catalogue and 105 are currently unpublished. These new discoveries, together with 210 previously known pulsars, result in a sample of $`N_{\mathrm{psr}}=730`$ pulsars, with a mean period of 0.76 s.

## 2. Pulsar current calculation

The observed distribution of pulsars in the Galaxy differs systematically from the true distribution due to various observational selection effects inherent in pulsar surveys (see the contribution by Lorimer in these proceedings for further details). We have accurately modeled the sensitivity threshold of the multibeam survey due to these effects and use this model to compute a weight or scale factor for each pulsar by: 1) placing it at large number of randomly selected locations in the model Galaxy; 2) calculating the effective dispersion measure and the interstellar scattering; and 3) recording the number of detections, i.e. those positions for which the predicted flux density exceeds the survey limit. For each pulsar, then, the scale factor is simply the ratio of the total number of locations in our model Galaxy (typically $`10^5`$) to the number of detections. This results in an estimate of the number of similar pulsars in the Galaxy.

The pulsar current, $`J`$, is the flow of objects from short to long spin periods. We can compute $`J`$ as a function of period, $`P`$, by binning pulsars into intervals of width $`\mathrm{\Delta }P`$. The current in a given bin is given by:

$$J(P)=\frac{1}{\mathrm{\Delta }P}\left(\underset{i=1}{\overset{n_{\mathrm{psr}}}{}}\frac{S_i\dot{P_i}}{f}\right),$$

where $`S_i`$ is the scale factor for each pulsar, $`n_{\mathrm{psr}}`$ is the number of pulsars in the bin and $`f`$ is the fraction of $`4\pi `$ sr covered by the emission beam. In this analysis, we assume for simplicity that $`f=0.2`$. Under steady-state conditions, the maximum value of $`J`$ is the birthrate of pulsars. The total number of active pulsars in the Galaxy is given by the sum of $`S_i/f`$. For further details, see Lorimer (these proceedings) and also Vivekanand & Narayan (1981).

## 3. Results

The results are summarized in Table 1 for both the TC93 distance model (Taylor & Cordes 1993) and the newer NE2001 model (Cordes & Lazio 2002). Our main conclusion is that pulsars with surface dipole magnetic fields greater than $`10^{12.5}`$ G account for half of the total birthrate. For both distance models, the birthrate of these high-field pulsars is 52% of the total population. Low-field ($`<10^{12}`$ G) pulsars contribute only 10% of the total birthrate using the TC93 model, or 17% for NE2001. Full details of this analysis are now being prepared for publication.

## References

Cordes, J. M. & Lazio, T. J. W. 2002, astro-ph/0207156

Manchester, R. N. et al. 2001, MNRAS, 328, 17

Taylor, J. H. & Cordes, J. M. 1993, ApJ, 411, 674

Vivekanand, M. & Narayan, R. 1981, J. Astrophys. Astr., 2, 315