STATISTICAL THEORY OF ATOMIC NUCLEI*)

P. Gombás

§ 1. Introduction

Statistical methods of calculation have long been used to explain the properties and regularities of atomic nuclei.

The first attempt to apply statistical theory to the atomic nucleus was made by Majorana¹, who obtained only a few qualitative, but very important, results. The deep foundations for the further development of the statistical theory of the nucleus were laid, on the one hand, by Heisenberg², who calculated explicitly the exchange energy for various attraction potentials between nucleons, and, on the other hand, by Weizsäcker³, who introduced a correction to the kinetic energy of the nucleus, the so-called “inhomogeneity correction.” Subsequent calculations⁴–⁸ proceeded from the basic assumptions developed in these works; moreover, as the attraction energy between two nucleons, expressions of the form \(-\gamma e^{-r_{ij}/r_0}\) and \(-\gamma e^{-(r_{ij}/r_0)^2}\) were for the most part chosen, where \(r_{ij}\) is the distance between nucleons, and \(\gamma\) and \(r_0\) are certain constants**). With the exception of work⁷, all these calculations were carried out under very crude assumptions; for example, a very simplified expression was used for the density of the distribution of nucleons, or it was even assumed that, when the “inhomogeneity correction” was neglected, the density was altogether constant. In this connection, the already-mentioned work⁷ represented a definite step forward, since no simplifying assumptions were made in it concerning the distribution density; on the contrary, the variation of the density and energy of nuclei was established by means of the Ritz method for light nuclei up to \(Si^{28}\). Common to all these calculations is the fact that with their aid the energy of nuclei can be represented

*) P. Gombás, Acta Physica Hungarica, 1, No. 4 (1952).
**) For further literature references, see, for example,⁹˒¹⁰ and also¹¹.

only in the form of such a formula which, after a certain choice of some free parameter, gives noticeable deviations from the experimental data either in the region of light nuclei or in the region of heavy nuclei. Thus, for example, in paper\(^7\), where the choice of the free parameter is made in the region of the lightest nuclei, already for \(\mathrm{Si}^{28}\)—the heaviest of the nuclei considered—the deviation of the energy from the experimental data is about 12%. The exception, of course, is those works whose aim is to find a formula, best agreeing with experiment, containing sometimes three, four, or even more empirical parameters. With the aid of such a formula it is, naturally, possible to achieve excellent agreement with the experimental data, beginning with the lightest and ending with the heaviest nuclei\(^*)\). However, by this means one obtains, in the main, only a formal description of empirical results, but not their explanation with the aid of general laws. From this point of view these works appear to be of little interest, although, of course, their heuristic significance cannot be denied.

The aim of the present work is the calculation of the binding energy of nuclei, and also the location of stable isobars, on the basis of statistical theory and under the assumption of a scalar Yukawa potential of attraction between nucleons. The calculations are based on the general principle of requiring a minimum of energy and are carried out by the Ritz method without introducing any arbitrary assumptions. In so doing, those parts of the interaction energy between neutrons and protons, neutrons and neutrons, and protons and protons which depend only on the mutual distance between the particles are taken to be equal. In the expression for the interaction energy of two nucleons only one constant is regarded as a free parameter, which at the same time is the only free parameter in the whole theory. It will be shown below that, with the appropriate choice of this parameter, results are obtained for the binding energy of nuclei which are in very good agreement with experiment both for the lightest and for the heaviest nuclei. The deviations from the experimental values, excluding the very light nuclei, lie within 7%. For the very light nuclei, whose binding energy is not a monotonic function of the number of particles, the statistical theory, naturally, can give, and in fact gives, only a certain average value in comparison with the strongly fluctuating empirical results. The position of the stable isobars is also determined in a very satisfactory manner.

The statistical theory developed in the present work is, by its very nature, applicable above all in the region of nuclei with large

\(^*)\) See, for example,\(^{10}\).

number of particles. The fact that it gives good results also for light nuclei is a consequence of our introduction of two corrections, namely: we introduced a correction to the exchange energy of the nucleon interaction, by means of which the energy associated with the self-exchange of particles was excluded; and, further, we introduced a correction also to the kinetic energy, by means of which the statistical expression for the kinetic energy of the He\(^4\) nucleus and of lighter nuclei goes over into the exact wave-mechanical expression. Owing to these corrections it is achieved that the statistical expression for the total energy of light nuclei well approximates the wave-mechanical expression and gives a good mean value of the empirical results.

The statistical theory of the atomic nucleus, in its essential features, is based on the behavior of a nucleon gas at absolute zero temperature and on the interaction between nucleons. Therefore, below we shall deal first of all with the general principles of the statistical method of considering a free nucleon gas and, in connection with this, shall examine the kinetic energy of a free nucleon gas. Then the general principles of the theory of nucleon interaction and the statistical calculation of the interaction energy will follow. Starting from these foundations, with the aid of variational methods, a statistical model of the nucleus will then be developed and the calculation of the energy and of the density distribution of nucleons will be carried out, beginning with the lightest and ending with the heaviest nuclei. The energy of some isobars will also be calculated.

§ 2. Fundamentals of the statistical method for considering a free nucleon gas

In the statistical model of the atomic nucleus, the nucleons are considered in the form of a nucleon gas at absolute zero temperature and are described by means of statistical methods. The statistical theory of neutron or proton gases is based on Fermi–Dirac statistics, which in turn is based on the Pauli principle and on the fact of the indistinguishability of particles. The Pauli principle, as is well known, means that in a completely quantized system each quantum state can be occupied by no more than one particle. The quantum state is here determined by specifying the state of spatial motion and the sign of the spin. The connection between quantum theory and statistics is given by the well-known proposition that to a volume \(h^3\) of phase space there correspond, taking spin into account, two neutron quantum states with antiparallel spins and two proton quantum states likewise with antiparallel spins; here \(h\) denotes Planck’s constant. A cell of phase space of volume \(h^3\) can thus be occupied by at most two neutrons and two protons.

In what follows we shall deal with a free nucleon gas, consisting of \(Q_n\) neutrons and \(Q_p\) protons, which are contained in a volume \(\Omega\). With respect to the walls of this volume we shall assume that they are impenetrable to nucleons. The nucleon gas, as has already been said, we shall regard as free, i.e., we shall assume that inside \(\Omega\) there is a constant potential, which may be taken equal to zero. Here we shall restrict ourselves to the case of very low temperatures, i.e., to the case when the particles occupy the lowest energy levels. At very low temperatures considerable simplifications arise, since in this case the Fermi–Dirac statistics is practically reduced to the condition that in an elementary cell of volume \(h^3\) there can be at most two neutrons and two protons.

The lowest energy levels can be described in the following way. According to our assumption, the space in which the nucleons are located contains no potentials, as a result of which the energy of any nucleon is purely kinetic energy, and for it the relation holds

\[ u = p^2/2M, \tag{1} \]

where \(p\) denotes the magnitude of the momentum, and \(M\) the mass of the nucleon; here \(M\) may be assumed the same for the proton and the neutron and equal to the mean value of the masses of the proton and neutron. The quantity \(u\) is a function of \(p\) and does not depend on the coordinates. Consequently one may restrict oneself to consideration in momentum space.

Let us first consider the case of neutrons; for protons everything will be entirely analogous. In view of the fact that the space is free from forces, all directions of motion of neutrons are equivalent. Since, moreover, the energy of neutrons depends only on the magnitude of the momentum \(p\), but not on the direction of the momentum, the lowest energy quantum states are located in momentum space inside a sphere whose center coincides with the origin of momentum space and whose radius \(p_{\mu n}\) is equal to the maximum momentum of the neutrons. Each of these lowest energy quantum states at absolute zero of temperature is occupied by at most one neutron; all quantum states outside the indicated sphere are free.

The quantity \(p_{\mu n}\) can be determined in the following way. The volume of the sphere in momentum space is equal to \(\dfrac{4}{3}\pi p_{\mu n}^3\). To the neutrons located in the volume \(\Omega\) there corresponds a volume of phase space equal to

\[ \Omega\,\frac{4}{3}\pi p_{\mu n}^3. \]

The number of quantum states is obtained by dividing this quantity by \(\dfrac{1}{2}h^3\). Since each of the quantum states under consideration contains one neutron, we have:

\[ 2\frac{4\pi p_{\mu n}^3}{3h^3}\,\Omega = Q_n. \tag{2} \]

Hence, for \(p_{\mu n}\) one obtains

\[ p_{\mu n}=\frac{1}{2}\left(\frac{3}{\pi}\right)^{1/3} h\rho_n^{1/3}, \tag{3} \]

where

\[ \rho_n=Q_n/\Omega \tag{4} \]

denotes the density of the neutron gas.

For protons everything proceeds in an entirely analogous manner. If the magnitude of the maximum momentum of the protons is denoted by \(p_{\mu p}\), then the relation will hold

\[ p_{\mu p}=\frac{1}{2}\left(\frac{3}{\pi}\right)^{1/3} h\rho_p^{1/3}, \tag{5} \]

where

\[ \rho_p=Q_p/\Omega \tag{6} \]

denotes the density of the proton gas.

§ 3. Kinetic Energy of the Nucleon Gas

We now proceed to the calculation of the kinetic energy of the nucleon gas and first of all take up the free nucleon gas already described in the preceding paragraph, which is contained in the volume \(\Omega\) and consists of \(Q_n\) neutrons and \(Q_p\) protons. First we shall consider a neutron gas alone. Since the neutrons may be regarded as completely free, the kinetic energy of the neutron gas \(U_K^n\) is nothing other than the so-called Fermi kinetic energy, which can be represented in the form\(^{1}\)

\[ U_K^n=\int u\,dn, \tag{7} \]

where \(u\) is a function of \(p\), and \(dn\) denotes the number of neutron quantum states to which the momentum magnitude between \(p\) and \(p+dp\) corresponds; the integration extends over all occupied neutron quantum states. Taking into account (1) and the relation

\[ dn=\frac{8\pi\Omega}{h^3}p^2\,dp, \tag{8} \]

for the Fermi kinetic energy of the neutrons we obtain:

\[ U_K^n=\frac{4\pi\Omega}{Mh^3}\int_0^{p_{\mu n}} p^4\,dp =\frac{4\pi\Omega}{5Mh^3}p_{\mu n}^{5}. \tag{9} \]

Substituting here expression (3) for \(p_{\mu n}\), we find:

\[ U_K^n=\varkappa\rho_n^{5/3}\Omega, \tag{10} \]

where $\varkappa_K$ denotes the constant

\[ \varkappa_K=\frac{3}{40}\left(\frac{3}{\pi}\right)^{2/3}\frac{h^2}{M}. \tag{11} \]

For a proton gas everything proceeds in exactly the same way; accordingly, for the Fermi kinetic energy of the protons we obtain the expression

\[ U_K^p=\varkappa_K \rho_p^{5/3}\Omega . \tag{12} \]

The Fermi kinetic energy of the entire nucleon gas is thus equal to

\[ U=U_K^n+U_K^p=\varkappa_K\left(\rho_n^{5/3}+\rho_p^{5/3}\right)\Omega . \tag{13} \]

Everything set forth above referred to a completely free nucleon gas, i.e., to a nucleon gas with a constant density distribution. The question immediately arises: what form does the Fermi kinetic energy take for a non-free nucleon gas? To derive this expression, let us divide the volume occupied by the nucleon gas, by means of a system of partitions, into partial volumes in such a way that in each partial volume element $dv$ there are still many neutrons and protons, and so that in this volume element the potential is practically constant. Then the nucleons located in individual volume elements may be regarded as a free nucleon gas. If the density of neutrons or protons in $dv$ is denoted by $\rho_n$ or $\rho_p$, then for the Fermi kinetic energy of the neutron (or proton) gas contained in $dv$, according to (10) and (12), one obtains the expression $\varkappa_K \rho_n^{5/3}dv$ (or $\varkappa_K \rho_p^{5/3}dv$). The Fermi kinetic energy of all neutrons (or all protons) is obtained by integrating this expression over the whole volume. Thus we have:

\[ U_K^n=\varkappa_K\int \rho_n^{5/3}dv;\qquad U_K^p=\varkappa_K\int \rho_p^{5/3}dv. \tag{14} \]

For the Fermi kinetic energy of the entire nucleon gas there follows the expression

\[ U_K=U_K^n+U_K^p=\varkappa_K\int\left(\rho_n^{5/3}+\rho_p^{5/3}\right)dv. \tag{15} \]

For a completely free nucleon gas, i.e., in the case of a constant density distribution, these expressions, of course, reduce to the corresponding expressions (10), (12), and (13).

Expression (15) gives the full kinetic energy of the nucleon gas only when the nucleons in the volume element $dv$ are completely free. If, however, the nucleons in the volume element $dv$ cannot be regarded as entirely free, then it is necessary to supplement this expression with the so-called “inhomogeneity correction”

Weizsäcker. This correction was justified in detail by Weizsäcker*), and we shall not go into details here. Instead, we shall choose a very simplified way of deriving this correction, which, precisely because of its simplicity, makes it possible to clarify the essence of the matter very well.

Let us again first consider neutrons. Without taking account of the Pauli principle, i.e., in the case of Bose statistics, all neutrons in equilibrium would occupy the lowest energy quantum level. If the eigenfunction of this state is denoted by \(\psi\), then for the density \(\rho_n\) and the Schrödinger kinetic energy \(U_j^n\) of the neutrons we shall have:

\[ \rho_n = Q_n |\psi|^2, \tag{16} \]

\[ U_j^n = Q_n \frac{h^2}{8\pi^2 M}\int |\operatorname{grad}\psi|^2\,dv. \tag{17} \]

Eliminating \(\psi\), we obtain:

\[ U_j^n = \chi_j \int \frac{(\operatorname{grad}\rho_n)^2}{\rho_n}\,dv, \tag{18} \]

where \(\chi_j\) denotes the constant

\[ \chi_j=\frac{h^2}{32\pi^2 M}. \tag{19} \]

For \(M\) here one may also take the mean value of the proton and neutron masses.

\(U_j^n\) is the Weizsäcker inhomogeneity correction for neutrons.

In an entirely analogous way we obtain the inhomogeneity correction for protons

\[ U_j^p=\chi_j \int \frac{(\operatorname{grad}\rho_p)^2}{\rho_p}\,dv. \tag{20} \]

The inhomogeneity correction for the whole nucleon gas is therefore equal to

\[ U_j = U_j^n + U_j^p. \tag{21} \]

It is obvious that \(U_j\) vanishes for a constant density of the nucleon distribution.

As follows from the derivation, \(U_j\) is the kinetic energy of a nucleon gas for whose particles the Pauli principle does not apply. Taking account of the Pauli principle, according to which not all particles can be in the energetically lowest quantum state, proceeds in such a way that, along with the quantity \(U_j\), in all calculations it is also necessary to take the Fermi kinetic energy \(U_k\), which is a consequence of the Pauli principle. With the aid of the energy \(U_k\), the particles will occupy, in accordance with the Pauli principle,

*) See, for example, 3, and also 12, where a list of the relevant literature is given.

higher-energy quantum states. In doing so it is assumed that \(\operatorname{grad}\rho_n\) and \(\operatorname{grad}\rho_p\) do not depend on the momentum vector, i.e. do not depend on the quantum state. All this demonstrates quite clearly the role of both constituent parts of the kinetic energy.

When both parts \(U_J\) and \(U_K\) of the kinetic energy are taken into account simultaneously, a certain error is nevertheless made, which is easiest to detect for the case of two neutrons and two protons. In this case both neutrons, as well as both protons, are in the lowest-energy quantum state, so that the energy \(U_J\) represents the total kinetic energy of the particles, while the energy \(U_K\) should vanish. This also follows from the fact that in this case the Pauli principle plays no role, since none of the particles is obliged to rise to a higher energy state, and therefore the part of the energy \(U_K\) due to the Pauli principle must vanish. In reality, however, \(U_K\) does not vanish. This circumstance is a consequence of the statistical method of investigation chosen by us, according to which in phase space, or, respectively, in momentum space, a finite volume already corresponds to a single particle; as a result, the kinetic energy of two particles in the lowest quantum state is contained not only in \(U_K\), but also in \(U_J\). Hence it is immediately clear that this error occurs not only for the case of two neutrons and two protons, but also for any other numbers of particles. Of course, it is significant only for small numbers of particles, whereas for large numbers of particles the error is small in comparison with the total kinetic energy.

The question arises: in what simplest way can this error be corrected? One way to do this is to set the Fermi kinetic energy of a particle in the lowest energy quantum state equal to the average Fermi kinetic energy per particle, and to subtract this quantity for both neutrons and protons in the lowest energy state. The average Fermi kinetic energy per neutron or per proton is \(U_K^n/Q_n\) or \(U_K^p/Q_p\). According to our proposal, the correction then consists in subtracting from the energy \(U_K\) the quantity \(2U_K^n/Q_n + 2U_K^p/Q_p\), or, equivalently, in (15) the quantity \(U_K^n\) is multiplied by the factor \(k_n = 1 - \dfrac{2}{Q_n}\), and \(U_K^p\) by the factor \(k_p = 1 - \dfrac{2}{Q_p}\). The total kinetic energy of our nucleon gas will therefore take the form

\[ U_K = k_n U_K^n + k_p U_K^p + U_J . \tag{22} \]

The factors \(k_n\) and \(k_p\) for \(Q_n=2\) and \(Q_p=2\) vanish, as a result of which in this case \(U_K\) disappears and the correction here thus leads to the exact value of the kinetic energy. Further, the quantity \(k_n U_K^n\) must also disappear for \(Q_n=1\), and the quantity \(k_p U_K^p\) for \(Q_p=1\), which, however, is not ensured by the correction factors. Consequently we must make the additional stipulation that the correction factor \(k_n\) also for \(Q_n<2\), and the correction factor \(k_p\) also for \(Q_p<2\), are to be set equal to zero. The kinetic energy then, in these cases as well, is given by the expression \(U\), which in these cases is in fact the exact expression for the kinetic energy. For large numbers of particles this correction becomes inexact, which, however, is immaterial, since for large numbers of particles it is small in comparison with the total kinetic energy.

§ 4. General foundations for calculating the interaction of nucleons

First of all we shall deal with the calculation, in the most general form, of the interaction of nucleons on the basis of wave mechanics, which we shall need for the statistical calculation of the energy of interaction of nucleons. We shall again assume that our nucleon gas consists of \(Q_n\) neutrons and \(Q_p\) protons. For the wave-mechanical description of the complete system we shall start from the Hartree–Fock approximation, i.e. we shall describe the system by means of one-particle eigenfunctions; namely, we shall suppose that \(Q_n\) functions \(\varphi_1, \varphi_2, \ldots, \varphi_{Q_n}\) are eigenfunctions of \(Q_n\) neutron states, and \(Q_p\) functions \(\psi_1, \psi_2, \ldots, \psi_{Q_p}\) are eigenfunctions of \(Q_p\) proton states. Taking account of the Pauli principle, the eigenfunction of the complete system can then be represented in the following form\({}^{13}\):

\[ \Psi=\frac{1}{\sqrt{Q_n!Q_p!}} \left| \begin{array}{cccc} \varphi_1(q_n^1) & \varphi_1(q_n^2) & \ldots & \varphi_1(q_n^{Q_n})\\ \varphi_2(q_n^1) & \varphi_2(q_n^2) & \ldots & \varphi_2(q_n^{Q_n})\\ \cdot & \cdot & \cdot & \cdot\\ \varphi_{Q_n}(q_n^1) & \varphi_{Q_n}(q_n^2) & \ldots & \varphi_{Q_n}(q_n^{Q_n}) \end{array} \right| \times \left| \begin{array}{cccc} \psi_1(q_p^1) & \psi_1(q_p^2) & \ldots & \psi_1(q_p^{Q_p})\\ \psi_2(q_p^1) & \psi_2(q_p^2) & \ldots & \psi_2(q_p^{Q_p})\\ \cdot & \cdot & \cdot & \cdot\\ \psi_{Q_p}(q_p^1) & \psi_{Q_p}(q_p^2) & \ldots & \psi_{Q_p}(q_p^{Q_p}) \end{array} \right|, \tag{23} \]

where \(q_n^i\) denotes the spatial \(\mathbf r_n^i(x_n^i, y_n^i, z_n^i)\) and spin \(s_n^i\) coordinates of the \(i\)-th neutron, while \(q_p^i\) denotes the spatial \(\mathbf r_p^i(x_p^i, y_p^i, z_p^i)\) and spin \(s_p^i\) coordinates of the \(i\)-th proton; \(\mathbf r_n^i\) and \(\mathbf r_p^i\) denote the coordinate vectors of the \(i\)-th neutron or the \(i\)-th proton. With respect to the one-particle eigenfunctions \(\varphi_i\) and \(\psi_i\), we may assume, without loss of generality, that both the neutron eigenfunctions \(\varphi_i\) among themselves and the proton eigenfunctions \(\psi_i\) among themselves are orthogonal. In addition, we shall assume that they are also normalized; then \(\Psi\) will also be normalized.

In the following derivation of the general expression for the interaction energy of nucleons, we shall first abstract from the Coulomb interaction of protons and restrict ourselves to considering only the nuclear interaction of nucleons. To derive the general expression for the interaction energy of nucleons, it is necessary to make a certain assumption concerning the nucleon interaction. This interaction must be chosen with account of the experimental data, which testify, first, to the saturation property of nuclear forces and, second, to their very small range of action. The first fact may be taken into account by an appropriate choice of the interaction mechanism, namely by assuming a certain exchange interaction between nucleons (cf. § 6). The second fact may be taken into account by an appropriate choice of the nucleon interaction potential. With respect to the nucleon interaction energy, we shall first make the completely general assumption that it is composed of pair interactions of particles in the following way:

\[ W = \sum_{i=1}^{Q_n} \sum_{k=1}^{Q_p} G_{ik}\bigl(q_n^i, q_p^k\bigr) + \frac{1}{2} \sum_{i=1}^{Q_n} \sum_{k \ne i}^{Q_n} N_{ik}\bigl(q_n^i, q_n^k\bigr) + \]

\[ + \frac{1}{2} \sum_{i=1}^{Q_p} \sum_{k \ne i}^{Q_p} P_{ik}\bigl(q_p^i, q_p^k\bigr), \tag{24} \]

where \(G_{ik}\) denotes the interaction energy of the \(i\)-th neutron and the \(k\)-th proton, \(N_{ik}\) the interaction energy of the \(i\)-th and \(k\)-th neutrons, and \(P_{ik}\) the interaction energy of the \(i\)-th and \(k\)-th protons. The factors \(1/2\) before the second and third sums on the right-hand side of (24) are included in order to avoid double counting of neutron and proton pairs.

With the aid of expression (24) and the eigenfunction (23), the wave-mechanical calculation of the mathematical expectation of the nucleon interaction energy \(W_A\) proceeds in complete analogy with the case of the interaction energy of electrons in many-electron atoms \(^{13}\).

We obtain:

\[ \begin{aligned} W_\Lambda^{*}={}& \sum_{i=1}^{Q_n}\sum_{k=1}^{Q_p} \iint \varphi_i^{*}(q_n)\psi_k^{*}(q_p)\, G_{ik}(q_n,q_p)\, \varphi_i(q_n)\psi_k(q_p)\,d\tau_n\,d\tau_p \\ &+\frac{1}{2}\sum_{i=1}^{Q_n}\sum_{k\ne i}^{Q_n} \iint \varphi_i^{*}(q_n)\varphi_k^{*}(q_n')\, N_{ik}(q_n,q_n')\, \varphi_i(q_n)\varphi_k(q_n')\,d\tau_n\,d\tau_n' \\ &-\frac{1}{2}\sum_{i=1}^{Q_n}\sum_{k\ne i}^{Q_n} \iint \varphi_k^{*}(q_n)\varphi_i^{*}(q_n')\, N_{ik}(q_n,q_n')\, \varphi_i(q_n)\varphi_k(q_n')\,d\tau_n\,d\tau_n' \\ &+\frac{1}{2}\sum_{i=1}^{Q_p}\sum_{k\ne i}^{Q_p} \iint \psi_i^{*}(q_p)\psi_k^{*}(q_p')\, P_{ik}(q_p,q_p')\, \psi_i(q_p)\psi_k(q_p')\,d\tau_p\,d\tau_p' \\ &-\frac{1}{2}\sum_{i=1}^{Q_p}\sum_{k\ne i}^{Q_p} \iint \psi_k^{*}(q_p)\psi_i^{*}(q_p')\, P_{ik}(q_p,q_p')\, \psi_i(q_p)\psi_k(q_p')\,d\tau_p\,d\tau_p' . \tag{25} \end{aligned} \]

Here, in integration over \(\tau_n\), \(\tau_p\), etc., in addition to integration over the spatial coordinates of neutrons and protons, summation over the spin coordinates is also understood. The first of the five double sums gives the interaction energy of neutrons and protons with one another; the second and third give the interaction energy of neutrons; and, finally, the fourth and fifth give the interaction energy of protons. It should be noted that the third and fifth double sums represent the so-called exchange energy of neutrons and protons, respectively. It occurs only between identical particles, and the corresponding terms arise because of the determinant representation of the proper function of these particles.

As can be shown (see § 6), the saturation property of the forces acting between nucleons can be taken into account by assuming a purely exchange interaction as the interaction between nucleons. Accordingly, we shall assign to the terms \(G_{ik}\) the form of a certain exchange interaction; namely, we choose for them an interaction of the Majorana type*), i.e., we shall assume that when the positions of a proton and a neutron are interchanged the spin turns out to be coupled with the charge. With the aid of the operator \(G_{ik}\), only the spatial coordinates of the \(i\)-th neutron and the \(k\)-th proton can then be exchanged. Thus we have:

\[ G_{ik}(q_n,q_p)=-J_{np}(r_{np})\,\Pi(\mathbf r_n,\mathbf r_p), \tag{26} \]

where \(\Pi(\mathbf r_n,\mathbf r_p)\) denotes the permutation operator, which

*) For more details on various types of nuclear forces see, for example, \(^{10}\) or \(^{11}\).

permutes only the positions of the neutron and proton, but leaves the spin coordinates unchanged; \(J_{np}(r_{np})\) is an ordinary function depending only on the mutual distance \(r_{np}=|\mathbf r_n-\mathbf r_p|\) of the neutron and proton; we shall discuss it in more detail below. With the aid of this expression for \(G_{ik}\) we have:

\[ G_{ik}(q_n,q_p)\,\varphi_i(\mathbf r_n,\sigma_n)\psi_k(\mathbf r_p,\sigma_p) = \varphi_i(\mathbf r_p,\sigma_n)\psi_k(\mathbf r_n,\sigma_p)J_{np}(r_{np}). \tag{27} \]

If, in the first approximation, the state of the spatial motion is regarded as independent of the spin state, then the eigenfunctions \(\varphi_i\) and \(\psi_i\) may be represented as the product of two functions: one depending only on the spatial coordinates, and the other depending only on the spin coordinate:

\[ \varphi_i=u_i(\mathbf r)\eta(\sigma),\qquad \psi_i=w_i(\mathbf r)\eta(\sigma). \tag{28} \]

Here \(u_i(\mathbf r)\) and \(w_i(\mathbf r)\) are the parts of the neutron and, respectively, proton eigenfunctions that depend only on the spatial coordinates; \(\eta(\sigma)\) denotes the spin function, which for the two values of the spin coordinate \(\sigma=+\tfrac12\) and \(\sigma=-\tfrac12\) we denote by \(\eta_+\) and \(\eta_-\), respectively. With the aid of (27) and (28) one can easily carry out the summation over the spin coordinates in the first integral of expression (25) and obtain, for the interaction energy of the \(i\)-th neutron and the \(k\)-th proton, the result

\[ \varepsilon_{ik}^{np} = \iint u_i^*(\mathbf r_n)w_k^*(\mathbf r_p) u_i(\mathbf r_p)w_k(\mathbf r_n) J_{np}(r_{np})\,dv_n\,dv_p, \tag{29} \]

where \(dv_n\) and \(dv_p\) denote the volume elements of the corresponding coordinate spaces. As is easy to see, \(\varepsilon_{ik}^{np}\) has the form of an exchange energy arising from the exchange of a neutron and a proton.

We now turn to the discussion of the neutron—neutron and proton—proton interactions, i.e. to the last four double sums in expression (25). This interaction between identical particles may be assumed either to depend on the spins or not to depend on them; moreover, let us recall once more that here we are for the time being disregarding the spin-independent Coulomb interaction. In the case of an interaction independent of the spins, together with the exchange energy of the interaction there would appear the ordinary energy of interaction of the particles, which would lead to a completely incorrect dependence of the binding energy of atomic nuclei on the number of particles. It is therefore necessary to assume that the interaction of identical particles depends on the spins. It proves expedient to make the following assumption:

\[ N_{ik}(q_n,q_n')=-\frac13\,J_{nn'}(r_{nn'})(S,S'), \tag{30} \]

\[ P_{ik}(q_p,q_p')=-\frac13\,J_{pp'}(r_{pp'})(S,S'). \tag{31} \]

Here \(q_n\) and \(q'_n\) denote the coordinates of the \(i\)-th and \(k\)-th neutrons, respectively, with the spin coordinates of the particles also included in \(q_n, q'_n\); \(r_{nn'}\) denotes the mutual distance of two neutrons, i.e. \(r_{nn'} = |\mathbf r_n - \mathbf r_{n'}|\); finally, \(S\) and \(S'\) denote the Pauli spin operators referring to the two particles. The quantities appearing in (31) with the indices \(p, p'\) refer to the \(i\)-th and \(k\)-th protons and have completely analogous meanings for them. The factor \(-\frac13\) has been introduced so that the interaction between two neutrons with oppositely directed spins, as well as between two protons with oppositely directed spins, should be equal to \(J_{nn}\) or, respectively, \(J_{pp}\). It is easy to see that this is indeed satisfied, since in this case \(S + S'\) is equal to zero, as a result of which the relation holds*)

\[ 2(S,S')=(S+S')^2-S^2-S'^2=-6. \tag{32} \]

To calculate the interaction energy of identical particles, we again represent the eigenfunction of a particle in the form of a product of a function depending only on the spatial coordinates and a spin function (cf. (28)). For later purposes we note that the spin functions satisfy the relations:

\[ S_x\eta_+ = \eta_-, \qquad S_x\eta_- = \eta_+, \tag{33} \]

\[ S_y\eta_+ = i\eta_-, \qquad S_y\eta_- = -i\eta_+, \tag{34} \]

\[ S_z\eta_+ = \eta_+, \qquad S_z\eta_- = -\eta_-. \tag{35} \]

Hence it follows, for example,

\[ (S,S')\,\eta_+(\sigma)\eta_+(\sigma')=\eta_+(\sigma)\eta_+(\sigma') \tag{36} \]

and

\[ (S,S')\,\eta_+(\sigma)\eta_-(\sigma') = -\eta_+(\sigma)\eta_-(\sigma')+2\eta_-(\sigma)\eta_+(\sigma'). \tag{37} \]

With the aid of these relations one can very easily calculate the interaction energy of identical particles. We shall carry out the calculations here only for the case of neutrons; for protons everything will be completely analogous. First of all, let us calculate the interaction energy of a neutron with positive spin, in the state \(\varphi_i = u_i\eta_+\), with two neutrons in the states \(\varphi_{k_1}=u_{k_1}\eta_+\) and \(\varphi_{k_2}=u_{k_2}\eta_-\), which differ from one another only in the directions of spin.

The ordinary interaction energy of both neutron pairs, according to the second line of formula (25), is represented in the form of the sum of the following—

*) \(S\) is twice the spin angular momentum of the nucleon, measured in units \(h/2\pi\). Taking into account that the eigenvalue of the square of the spin angular momentum in units \(h/2\pi\) is \(l(l+1)=3\), we obtain \(S^2=3\).

following double integrals:

\[ -\frac{1}{3}\iint |u_i(\mathbf r)|^2 |u_k(\mathbf r')|^2 \eta_+(\sigma)\eta_+(\sigma')\times \]
\[ \times (S,S')\eta_+(\sigma)\eta_+(\sigma')J_{nn}(r_{nn'})\,d\tau_n\,d\tau'_{n}, \tag{38} \]

\[ -\frac{1}{3}\iint |u_i(\mathbf r)|^2 |u_k(\mathbf r')|^2 \eta_+(\sigma)\eta_-(\sigma')\times \]
\[ \times (S,S')\eta_+(\sigma)\eta_-(\sigma')J_{nn}(r_{nn'})\,d\tau_n\,d\tau'_{n}, \tag{39} \]

where it is necessary to remember that the indicated integration also denotes summation over the spin coordinates. Taking into account (36) and (37), as well as the fact that the spin functions \(\eta_+\) and \(\eta_-\) are orthogonal to each other, we find that summation over the spin coordinates gives \(+1\) in the first integral and \(-1\) in the second integral. Since in all other respects both integrals are identical, the double integrals in the second double sum of expression (25) cancel one another. At the same time, the entire second double sum in (25) also vanishes, i.e., the ordinary interaction energy of the neutrons.

The exchange energy of both neutron pairs is obtained according to the third line of (25) in the form of the sum of the following double integrals:

\[ \frac{1}{3}\iint u_k^*(\mathbf r)u_i^*(\mathbf r')u_i(\mathbf r)u_k(\mathbf r') \eta_+(\sigma)\eta_+(\sigma')\times \]
\[ \times (S,S')\eta_+(\sigma)\eta_+(\sigma')J_{nn}(r_{nn'})\,d\tau_n\,d\tau'_{n}, \tag{40} \]

\[ \frac{1}{3}\iint u_k^*(\mathbf r)u_i^*(\mathbf r')u_i(\mathbf r)u_k(\mathbf r') \eta_-(\sigma)\eta_+(\sigma')\times \]
\[ \times (S,S')\eta_+(\sigma)\eta_-(\sigma')J_{nn}(r_{nn'})\,d\tau_n\,d\tau'_{n}. \tag{41} \]

In the same way as in the preceding case, after summation over the spin coordinates one obtains \(+1\) in the first integral and \(+2\) in the second. Since in all other respects the integrals here too are identical, for the sum of the exchange energies of both neutron pairs, i.e., for the sum of (40) and (41), we obtain:

\[ \iint u_k^*(\mathbf r)u_i^*(\mathbf r')u_i(\mathbf r)u_k(\mathbf r')J_{nn}(r_{nn'})\,d\tau\,d\tau'. \tag{42} \]

This expression is the exchange energy of a neutron in the state \(\varphi_i=u_i\eta_+\), with two neutrons in the states \(\varphi_{k_1}=u_k\eta_+\) and \(\varphi_{k_2}=u_k\eta_-\), i.e., with two neutrons which are in one and the same state of spatial motion but possess antiparallel spins. The average exchange energy between neutrons \(i\) and \(k\) is therefore equal to

\[ \varepsilon_{ik}^{m}=\frac{1}{2}\iint u_k^*(\mathbf r)u_i^*(\mathbf r')u_i(\mathbf r)u_k(\mathbf r')J_{nn}(r_{nn'})\,d\tau\,d\tau'. \tag{43} \]

STATISTICAL THEORY OF ATOMIC NUCLEI

With the aid of assumption (31), all calculations for protons are carried out in a completely analogous way. The ordinary interaction energy also vanishes in the case of protons, and for the mean exchange energy of two protons \(i\) and \(k\) one obtains the expression

\[ \varepsilon^{pp}_{ik} = \frac{1}{2} \iint w_k^*(\mathbf r)\, w_i^*(\mathbf r')\, w_i(\mathbf r)\, w_k(\mathbf r')\, J_{pp}(r_{pp'})\, dv\, dv' . \tag{44} \]

For the total energy \(W_A^{np}\) arising from the neutron—proton interaction, as well as for the total energies \(W_A^{nn}\) or \(W_A^{pp}\) arising from the neutron—neutron or, respectively, proton—proton interaction, we obtain:

\[ W_A^{np} = \sum_{i=1}^{Q_n}\sum_{k=1}^{Q_p} \varepsilon^{np}_{ik}, \qquad W_A^{nn} = \frac{1}{2} \sum_{i=1}^{Q_n}\sum_{k\ne i}^{Q_n} \varepsilon^{nn}_{ik}, \]

\[ W_A^{pp} = \frac{1}{2} \sum_{i=1}^{Q_p}\sum_{k\ne i}^{Q_p} \varepsilon^{pp}_{ik}. \tag{45} \]

Hence, for the total energy arising from the interaction of all nucleons, we obtain:

\[ W_A = W_A^{np} + W_A^{nn} + W_A^{pp}. \tag{46} \]

We must now still specify in greater detail those parts of the nuclear interaction which depend only on the mutual distance between nucleons, i.e. the functions \(J_{np}(r_{np})\), \(J_{nn}(r_{nn})\), and \(J_{pp}(r_{pp})\). Since, according to experimental data, nuclear forces have a very small range of action, namely, the region of their action is of the order of \(10^{-13}\) cm, it is necessary to choose for \(J_{np}\), \(J_{nn}\), and \(J_{pp}\) functions which decrease very rapidly as the mutual distance of the nucleons increases. We shall choose here for the interaction of nucleons the scalar central Yukawa potential, i.e. we shall assume that the exchange interaction of neutrons with protons is effected by means of the exchange of electrically charged mesons, namely \(\pi\)-mesons, and that the exchange interaction between identical nucleons is effected by means of the exchange of electrically neutral mesons of the same mass. Further, we shall assume that the part of the interaction between neutrons—protons, neutrons—neutrons, and protons—protons which depends only on the mutual distance of the nucleons is the same for all three cases. In other words, we shall put:

\[ J_{np} = J_{nn} = J_{pp} = J = \gamma \frac{1}{|\mathbf r-\mathbf r'|} e^{-|\mathbf r-\mathbf r'|/r_0}, \tag{47} \]

where \(|\mathbf r-\mathbf r'|\) denotes the mutual distance between the corresponding nucleons, and \(\gamma\) is a certain constant of the dimension of energy,

multiplied by a length; \(r_0\) is likewise a certain constant, equal, according to Yukawa’s theory, to the Compton wavelength of the \(\pi\)-meson divided by \(2\pi\). If for the mass of the \(\pi\)-meson, \(M_\pi\), we take the value of 285 electron masses,\(^{14}\) then we have:

\[ r_0=\frac{h}{2\pi M_\pi c}=1.355\cdot 10^{-13}\ \text{cm}, \tag{48} \]

where \(c\) denotes the velocity of light. We shall make use below of the constant \(\gamma\), which has not yet been defined more precisely.

Here it may be noted that, in applications of statistical theory to atomic nuclei, in most works the expression usually chosen for the interaction energy of two nucleons was

\[ J=-\varepsilon e^{-|\mathbf r-\mathbf r'|/r_0} \tag{49} \]

or

\[ J=-\varepsilon e^{-(\mathbf r-\mathbf r'/r_0)^2}, \tag{50} \]

where \(\varepsilon\) is a certain constant of the dimension of energy, and \(r_0\) is again a certain length of order \(10^{-13}\ \text{cm}\).

In addition to the interaction between nucleons discussed up to now, there is also the Coulomb electrostatic interaction between protons, which is a consequence of the presence of an electric charge on protons. This interaction, in turn, is formed in the known way from pair interactions of particles, quite analogously to what we assumed for the interaction of nucleons considered up to now, and the consequence of which was the formal writing of expression (24). While in this respect both interactions behave in the same way, in another respect there is an essential difference between them. Indeed, the interaction considered above between identical nucleons had to be assumed to depend on spins, whereas the electrostatic Coulomb interaction of protons has no spin dependence. Correspondingly, in this case there appears not only the exchange energy, but also the ordinary energy of interaction of the particles. To compute the Coulomb energy of interaction it is necessary to substitute in expression (25), instead of \(I_{ik}\), the Coulomb interaction energy of two protons \(e^2/|\mathbf r-\mathbf r'|\), where \(e\) denotes the elementary positive charge. For the ordinary electrostatic Coulomb energy of interaction of the protons we then obtain:

\[ W_C=\frac{1}{2}e^2\sum_{i=1}^{Q_p}\sum_{\substack{k=1\\k\ne i}}^{Q_p}\iint |\psi_i(q_p)|^2|\psi_k(q'_p)|^2 \frac{1}{|\mathbf r_p-\mathbf r'_p|}\,d\tau_p\,d\tau'_p \tag{51} \]

and for the exchange energy arising from the electrostatic interaction ...

for the energy we have:

\[ W_R=-\frac{1}{2}e^2 \sum_{i=1}^{Q_p}\sum_{k\ne i}^{Q_p} \iint \psi_k^*(q_p)\psi_i^*(q'_p)\psi_i(q_p)\psi_k(q'_p) \times \frac{1}{\left|\mathbf r_p-\mathbf r'_p\right|} \,d\tau_p\,d\tau'_p . \tag{52} \]

Assuming (28), we obtain for \(W_C\) the expression

\[ W_C=\frac{1}{2}e^2 \sum_{i=1}^{Q_p}\sum_{k\ne i}^{Q_p} \iint |w_i(\mathbf r_p)|^2\,|w_k(\mathbf r'_p)|^2 \frac{1}{\left|\mathbf r_p-\mathbf r'_p\right|} \,d\tau_p\,d\tau'_p . \tag{53} \]

For \(W_R\), under assumption (28), it follows that only those integrals vanish in which the integration extends over states with parallel spins. We obtain:

\[ W_R=-\frac{1}{2}e^2 \sum_{\substack{i=1\\ \sigma_k=\sigma_i}}^{Q_p} \sum_{k\ne i}^{Q_p} \iint w_k^*(\mathbf r_p)w_i^*(\mathbf r'_p)w_i(\mathbf r_p)w_k(\mathbf r'_p) \times \frac{1}{\left|\mathbf r_p-\mathbf r'_p\right|} \,d\tau_p\,d\tau'_p , \tag{54} \]

where the double sum must be extended only over two states that have parallel spins. The fact that in \(W_R\), in contrast to \(W_A^{pp}\) (or \(W_A^{nn}\)), the exchange interaction arises from the interaction of particles only with parallel spins follows from the fact that the Coulomb interaction is independent of the spins, unlike the nucleon interaction \(P_{ik}\) (or \(N_{ik}\)).

Thus, we have at our disposal a wave-mechanical expression for the interaction energy between nucleons. Regarding the ordinary Coulomb interaction energy \(W_C\), we shall say nothing further. However, we must still examine in detail the calculation of the exchange energy, in particular the exchange energy of a free nucleon gas.

Let us recall here that we have not yet shown in what way the exchange interaction leads to the saturation property of nuclear forces. We shall discuss this in one of the following paragraphs; moreover, it is simplest to show this on the basis of the explicit expression for the energy of the exchange interaction.

§ 5. Interaction energy of a nucleon gas

In this paragraph we shall undertake the concrete calculation of the interaction energy of particles of the previously considered free nucleon gas, which consists of \(Q_n\) neutrons and \(Q_p\) protons and is contained in a volume \(\Omega\). We shall again assume that the gas is under

at absolute zero temperature, i.e. in the very lowest energy state. This means that each neutron state of motion is occupied by two neutrons with antiparallel spins, and each proton state of motion is occupied by two protons with antiparallel spins.*) Any \(k\)-th state of motion of a free neutron or proton may be described, for example, by means of plane waves

\[ u_k(\mathbf r)=\frac{1}{\sqrt{\Omega}}\,e^{\frac{2\pi i}{h}(\mathbf p_k\mathbf r)},\qquad w_k(\mathbf r)=\frac{1}{\sqrt{\Omega}}\,e^{\frac{2\pi i}{h}(\mathbf p_k\mathbf r)}, \tag{55} \]

where \(\mathbf p_k\) denotes the momentum vector of the neutron or proton. These wave functions do not satisfy the boundary conditions that the eigenfunctions vanish on the planes bounding the volume \(\Omega\). If, however, \(\Omega\) is assumed to be very large, then these boundary conditions can already be satisfied by means of a small change in the wave functions, namely by integrating the wave functions over an infinitely small region of variation of the propagation vector. As a result, the disappearance of the wave functions at infinity is ensured. Further, the eigenfunctions of different neutron states, as well as of different proton states, must be mutually orthogonal. This is achieved by the fact that, for example, for neutrons we shall regard the momentum vectors \(\mathbf p_i\) and \(\mathbf p_k\) in phase space as belonging to different cells of momentum space.

To calculate the exchange energy \(W_A\) of our free nucleon gas, we introduce, for simplicity, the following notation:

\[ \begin{aligned} \rho_{ik}^{np}(\mathbf r)&=u_i^*(\mathbf r)w_k(\mathbf r),\\ \rho_{ik}^{nn}(\mathbf r)&=u_i^*(\mathbf r)u_k(\mathbf r),\\ \rho_{ik}^{pp}(\mathbf r)&=w_i^*(\mathbf r)w_k(\mathbf r). \end{aligned} \left. \right\} \tag{56} \]

With this notation one may write \(\varepsilon_{ik}^{np}\), \(\varepsilon_{ik}^{nn}\), and \(\varepsilon_{ik}^{pp}\) in the following form:

\[ \varepsilon_{ik}^{np}= \iint \rho_{ik}^{np}(\mathbf r)\rho_{ik}^{np*}(\mathbf r')J(|\mathbf r-\mathbf r'|)\,dv\,dv', \tag{57} \]

\[ \varepsilon_{ik}^{nn}= \iint \rho_{ik}^{nn}(\mathbf r)\rho_{ik}^{nn*}(\mathbf r')J(|\mathbf r-\mathbf r'|)\,dv\,dv', \tag{58} \]

\[ \varepsilon_{ik}^{pp}= \iint \rho_{ik}^{pp}(\mathbf r)\rho_{ik}^{pp*}(\mathbf r')J(|\mathbf r-\mathbf r'|)\,dv\,dv'. \tag{59} \]

To find \(W_A\), we must compute these three integrals,

*) In the case when the number of neutrons or protons is odd, the relatively highest energy state is occupied by only one particle; this, however, below plays no role in view of the assumption of a large number of particles.

where in the transition densities \(\rho^{np}_{ik}\), \(\rho^{pn}_{ik}\), and \(\rho^{pp}_{ik}\) we shall substitute plane waves of the form (55).

We shall begin with the calculation of \(\varepsilon^{np}_{ik}\), i.e., the exchange energy of the \(i\)-th neutron and the \(k\)-th proton, which, taking (47) into account, can be written in the following form:

\[ \varepsilon^{np}_{ik} = -\gamma \int \rho^{np}_{ik}(\mathbf r)\, \rho^{np*}_{ik}(\mathbf r')\, \frac{e^{-|\mathbf r-\mathbf r'|/r_0}}{|\mathbf r-\mathbf r'|} \,d\tau\,d\tau', \tag{60} \]

where

\[ \rho^{np}_{ik}(\mathbf r) = \frac{1}{\Omega} e^{\frac{2\pi i}{h}(\mathbf p_k-\mathbf p_i,\mathbf r)} . \tag{61} \]

As the distribution potential of \(\rho^{np*}_{ik}(\mathbf r)\) one may introduce the quantity

\[ V_{ik}(\mathbf r) = \int \rho^{np*}_{ik}(\mathbf r')\, \frac{e^{-|\mathbf r-\mathbf r'|/r_0}}{|\mathbf r-\mathbf r'|} \,d\tau' \tag{62} \]

and with its aid represent \(\varepsilon^{np}_{ik}\) in the form

\[ \varepsilon^{np}_{ik} = -\gamma \int \rho_{ik}(\mathbf r)V_{ik}(\mathbf r)\,d\tau . \tag{63} \]

The calculation of the potential \(V_{ik}\) presents no difficulty. Quite analogously to the fact that the potential of a certain spatial charge distribution with Coulomb forces satisfies Poisson’s equation, the potential \(V_{ik}\) of the distribution \(\rho^{np*}_{ik}\) with Yukawa forces satisfies the generalized Poisson equation of the form

\[ \Delta V_{ik} - \frac{1}{r_0^2}V_{ik} = -4\pi \rho^{np*}_{ik}. \tag{64} \]

The solution of this equation is

\[ V_{ik}(\mathbf r) = \frac{h^2}{\Omega\pi}\, \frac{1}{p_0^2+|\mathbf p_i-\mathbf p_k|^2}\, \rho^{np*}_{ik}, \tag{65} \]

where

\[ p_0=\frac{h}{2\pi}\cdot \frac{1}{r_0}. \tag{66} \]

Substituting (65) into (63), we obtain:

\[ \varepsilon^{np}_{ik} = -\frac{\gamma h^2}{\Omega\pi}\, \frac{1}{p_0^2+|\mathbf p_i-\mathbf p_k|^2}. \tag{67} \]

The total exchange energy \(W^{np}_A\) arising from the exchange interaction of neutrons and protons is obtained, according to (45),

by summing \(\varepsilon^{np}_{ik}\) over all completely occupied neutron and proton states, i.e., by summing over all neutron momenta \(\mathbf p_i\) and over all proton momenta \(\mathbf p_k\).

We carried out the derivation of the expression for \(\varepsilon^{np}_{ik}\) in the preceding paragraphs on the basis of wave mechanics; the calculation of \(W^{np}_A\), i.e., the summation over neutron and proton momenta, will be performed on the basis of statistical considerations. This means that we shall replace the summation over \(\mathbf p_i\) and \(\mathbf p_k\) by integration, which can be done in the following way. If the absolute values of \(\mathbf p_i\) and \(\mathbf p_k\) are denoted by \(p_i\) and \(p_k\), and the angle between \(\mathbf p_i\) and \(\mathbf p_k\) by \(\vartheta\), then we shall have

\[ |\mathbf p_i-\mathbf p_k|^2=p_i^2+p_k^2-2p_ip_k\cos\vartheta . \tag{68} \]

Let us introduce a polar coordinate system with polar axis directed along \(\mathbf p_i\). The number of neutrons whose momentum directions lie between \(\vartheta\) and \(\vartheta+d\vartheta\), and whose momentum magnitudes lie between \(p_i\) and \(p_i+dp_i\), is equal to

\[ dn=\frac{4\pi\Omega}{h^3}\,p_i^2\sin\vartheta\,dp_i\,d\vartheta . \tag{69} \]

With the aid of (69) one may replace the summation over \(p_i\) by the following integration:

\[ \sum_{i=1}^{Q_n}\varepsilon^{np}_{ik} = -\frac{4\gamma}{h} \int_0^{p_{\mu n}} dp_i\,p_i^2 \int_0^\pi \frac{\sin\vartheta\,d\vartheta} {p_0^2+p_i^2+p_k^2-2p_ip_k\cos\vartheta} = \]

\[ = -\frac{2\gamma}{h} \left[ \frac{1}{2p_k}(p_{\mu n}^2-p_k^2-p_0^2) \ln \frac{(p_{\mu n}+p_k)^2+p_0^2} {(p_{\mu n}-p_k)^2+p_0^2} + \right. \]

\[ \left. {}+2p_{\mu n} -2p_0\operatorname{arctg}\frac{p_{\mu n}+p_k}{p_0} -2p_0\operatorname{arctg}\frac{p_{\mu n}-p_k}{p_0} \right], \tag{70} \]

where \(p_{\mu n}\) denotes the magnitude of the maximum momentum of the neutrons.

This expression gives the exchange energy of all \(Q_n\) neutrons with a proton having momentum \(\mathbf p_k\). The total energy arising from the exchange interaction of all \(Q_n\) neutrons with all \(Q_p\) protons will be obtained if this expression is summed over the proton momenta \(\mathbf p_k\). We shall again replace this summation by integration. For this purpose we must multiply the expression standing on the right-hand side of (70) by \(8\pi\Omega p_k^2/h^3\) and integrate with respect to \(p_k\) from zero to \(p_{\mu p}\), where \(p_{\mu p}\) denotes the magnitude of the maximum momentum of the protons.

After rather tedious integration one obtains

\[ W_A^{np}=\sum_{i=1}^{Q_n}\sum_{k=1}^{Q_p}\varepsilon_{ik}^{np}= \]

\[ =-\frac{16\pi\gamma\Omega}{h^4}\int_{0}^{p_{\mu p}} \left[ 2p_{\mu n}+\frac{1}{2p_k}\left(p_{\mu n}^{\,2}-p_k^{\,2}+p_0^{\,2}\right) \ln\frac{p_0^{\,2}+(p_{\mu n}+p_k)^2}{p_0^{\,2}+(p_{\mu n}-p_k)^2} \right. \]

\[ \left. -2p_0\operatorname{arctg}\frac{p_{\mu n}+p_k}{p_0} -2p_0\operatorname{arctg}\frac{p_{\mu n}-p_k}{p_0} \right]p_k^{\,2}\,dp_k= \]

\[ =-\frac{16\pi\gamma\Omega}{h^4} \left\{ \frac{1}{2}\left(p_{\mu n}^{\,3}p_{\mu p}+p_{\mu n}p_{\mu p}^{\,3}\right) -\frac{1}{6}p_0^{\,2}p_{\mu n}p_{\mu p} +\right. \]

\[ \left. +\frac{1}{24}\left[p_0^{\,4}+6p_0^{\,2}\left(p_{\mu n}^{\,2}+p_{\mu p}^{\,2}\right) -3\left(p_{\mu n}^{\,2}-p_{\mu p}^{\,2}\right)^2\right] \ln\frac{p_0^{\,2}+(p_{\mu n}+p_{\mu p})^2}{p_0^{\,2}+(p_{\mu n}-p_{\mu p})^2} \right. \]

\[ \left. -\frac{2}{3}p_0\left(p_{\mu n}^{\,3}+p_{\mu p}^{\,3}\right) \operatorname{arctg}\frac{p_{\mu n}+p_{\mu p}}{p_0} + \right. \]

\[ \left. +\frac{2}{3}p_0\left(p_{\mu n}^{\,3}-p_{\mu p}^{\,3}\right) \operatorname{arctg}\frac{p_{\mu n}-p_{\mu p}}{p_0} \right\}. \tag{71} \]

If we introduce the dimensionless functions

\[ \omega_n=\frac{p_{\mu n}}{p_0}=(3\pi^2)^{1/3}r_0\rho_n^{1/3} \quad\text{and}\quad \omega_p=\frac{p_{\mu p}}{p_0}=(3\pi^2)^{1/3}r_0\rho_p^{1/3}, \tag{72} \]

then \(W_A^{np}\) can be written in the following form:

\[ W_A^{np}=-\frac{\gamma\Omega}{r_0^{\,4}}\frac{1}{6\pi^3} \left\{ 3\left(\omega_n^{\,3}\omega_p+\omega_n\omega_p^{\,3}\right)-\omega_n\omega_p+ \right. \]

\[ \left. +\frac{1}{4}\left[1+6\left(\omega_n^{\,2}+\omega_p^{\,2}\right) -3\left(\omega_n^{\,2}-\omega_p^{\,2}\right)^2\right] \left\{\ln\left[1+(\omega_n+\omega_p)^2\right]\right. \right. \]

\[ \left. \left. -\ln\left[1+(\omega_n-\omega_p)^2\right]\right\} -4\left(\omega_n^{\,3}+\omega_p^{\,3}\right)\operatorname{arctg}(\omega_n+\omega_p)+ \right. \]

\[ \left. +4\left(\omega_n^{\,3}-\omega_p^{\,3}\right)\operatorname{arctg}(\omega_n-\omega_p) \right\}. \tag{73} \]

The calculation of \(W_A^{nn}\) and \(W_A^{pp}\) is carried out in an entirely analogous manner. First of all, in these cases as well, \(\varepsilon_{ik}^{nn}\) and \(\varepsilon_{ik}^{pp}\) are calculated, for which one obtains formally the same expressions as for \(\frac{1}{2}\varepsilon_{ik}^{np}\) [see (67)], with, however, the essential difference that in the expression for \(\varepsilon_{ik}^{nn}\) both momenta \(\mathbf p_i\) and \(\mathbf p_k\) denote neutron momenta, while in the expression for \(\varepsilon_{ik}^{pp}\) they denote proton momenta. Thereafter everything proceeds entirely analogously to the calculation of \(W_A^{np}\), except that in the case \(W_A^{nn}\) one has to integrate with respect to \(p_i\) and \(p_k\) both times from zero to \(p_{\mu n}\), and in the case \(W_A^{pp}\) both times from zero

to \(p_{\mu p}\). Taking (45) into account, we obtain:

\[ W_A^{nn}=-\frac{4\pi\gamma\Omega}{h^4}\left[ p_{\mu n}^4-\frac{1}{6}p_0^2p_{\mu n}^2+ \frac{1}{24}\left(p_0^4+12p_0^2p_{\mu n}^2\right) \ln\frac{p_0^2+4p_{\mu n}^2}{p_0^2} -\frac{4}{3}p_0p_{\mu n}^3\operatorname{arctg}\frac{2p_{\mu n}}{p_0} \right], \tag{74} \]

\[ W_A^{pp}=-\frac{4\pi\gamma\Omega}{h^4}\left[ p_{\mu p}^4-\frac{1}{6}p_0^2p_{\mu p}^2+ \frac{1}{24}\left(p_0^4+12p_0^2p_{\mu p}^2\right) \ln\frac{p_0^2+4p_{\mu p}^2}{p_0^2} -\frac{4}{3}p_0p_{\mu p}^3\operatorname{arctg}\frac{2p_{\mu p}}{p_0} \right]. \tag{75} \]

If we introduce here as well the dimensionless functions \(\omega_n\) and \(\omega_p\), then we obtain:

\[ W_A^{nn}=-\frac{\gamma\Omega}{r_0^4}\frac{1}{24\pi^3} \left[ 6\omega_n^4-\omega_n^2+ \frac{1}{4}\left(1+12\omega_n^2\right)\ln\left(1+4\omega_n^2\right) -8\omega_n^3\operatorname{arctg}2\omega_n \right], \tag{76} \]

\[ W_A^{pp}=-\frac{\gamma\Omega}{r_0^4}\frac{1}{24\pi^3} \left[ 6\omega_p^4-\omega_p^2+ \frac{1}{4}\left(1+12\omega_p^2\right)\ln\left(1+4\omega_p^2\right) -8\omega_p^3\operatorname{arctg}2\omega_p \right]. \tag{77} \]

Everything set forth above is valid only for the case of a free nucleon gas, i.e. for the case when the density of the distribution of nucleons is constant. For a non-free nucleon gas the exchange energy can be calculated analogously to the way we proceeded in the case of the Fermi kinetic energy of a non-free gas. As before, we divide the volume occupied by the gas, by means of some system of partitions, into partial volumes \(dv\), containing still many neutrons and protons, which can already be regarded as a free nucleon gas, i.e. as a nucleon gas with constant density. The exchange energy of the nucleons in the separate partial volumes can be represented in the form (73), (76), and (77), with the difference that instead of \(\Omega\) there stands everywhere \(dv\). The total exchange energy is then obtained by integration over the entire volume occupied by the nucleon gas. If, for brevity, we introduce the notation

\[ f(\omega_n,\omega_p)=\frac{1}{24\pi^3}\left\{ 3\left(\omega_n^3\omega_p+\omega_n\omega_p^3\right)-\omega_n\omega_p+ \right. \]

\[ +\frac{1}{4}\left[1+6(\omega_n+\omega_p)^2-3(\omega_n^2-\omega_p^2)^2\right] \left\{\ln\left[1+(\omega_n+\omega_p)^2\right] -\ln\left[1+(\omega_n-\omega_p)^2\right]\right\} \]

\[ \left. -4(\omega_n^3+\omega_p^3)\operatorname{arctg}(\omega_n+\omega_p) +4(\omega_n^3-\omega_p^3)\operatorname{arctg}(\omega_n-\omega_p) \right\}, \tag{78} \]

then we shall have:

\[ W_A^{np}=-\frac{4\gamma}{r_0^4}\int f(\omega_n,\omega_p)\,dv, \tag{79} \]

\[ W_A^{nn}=-\frac{\gamma}{r_0^4}\int f(\omega_n,\omega_n)\,dv, \tag{80} \]

\[ W_A^{pp}=-\frac{\gamma}{r_0^4}\int f(\omega_p,\omega_p)\,dv. \tag{81} \]

Whereas in (45), in the expressions for \(W_A^{nn}\) and \(W_A^{pp}\), the terms arising from the self-exchange of neutrons or protons drop out owing to the cancellation of the terms with \(k=i\) in the double sum, expressions (80) and (81), like the corresponding expressions preceding them, also contain the energy arising from self-exchange, since in both cases one has to integrate twice over one and the same interval of momenta. In view of the absence of an exact method for introducing the correction for self-exchange, one may proceed as follows: equate the energy arising from the self-exchange of one neutron or proton to the mean exchange energy per one neutron or proton, i.e. to the quantity \(W_A^{nn}/Q_n^2\) or \(W_A^{pp}/Q_p^2\).

Thus the energy arising from the self-exchange of all neutrons or protons is \(W_A^{nn}/Q_n\) or \(W_A^{pp}/Q_p\). To correct for self-exchange we must subtract these energies from (80) and (81), and therefore the correction for self-exchange consists in simply multiplying expressions (80) and (81) by the quantities \(s_n=1-\frac{1}{Q_n}\) or, respectively, by \(s_p=1-\frac{1}{Q_p}\). Thus, for the total exchange energy one obtains

\[ W_A=W_A^{np}+s_n W_A^{nn}+s_p W_A^{pp}. \tag{82} \]

Of course, the correction introduced in this way is very crude, although it should nevertheless be emphasized that for the case \(Q_n=1\) and \(Q_p=1\), when the correction plays a relatively most important role, the exact value is obtained, since in this case the correction factors \(s_n\) and \(s_p\) become zero, and with them the corrected expressions \(s_n W_A^{nn}\) and \(s_p W_A^{pp}\) vanish, as they should. For large numbers of particles this correction gives only a rough approximate value, which, however, is of little significance, since, obviously, the larger the number of particles, the smaller the significance of this correction.

For completeness we shall also give here the functions \(f(\omega_n,\omega_p)\) for the case of an interaction energy of the form (49) and (50). In the case of an energy

for the interaction (49) we have:

\[ f(\omega_n,\omega_p)=-\frac{r_0}{24\pi^3}\left\{4\omega_n\omega_p-\left[1+3(\omega_n^2+\omega_p^2)\right]\left\{\ln\left[1+(\omega_n+\omega_p)^2\right]\right.\right. \]
\[ \left.\left.-\ln\left[1+(\omega_n-\omega_p)^2\right]\right\}+4(\omega_n^3+\omega_p^3)\operatorname{arctg}(\omega_n+\omega_p)\right. \]
\[ \left.-4(\omega_n^3-\omega_p^3)\operatorname{arctg}(\omega_n-\omega_p)\right\} \tag{83} \]

and in the case of the interaction energy (50) we find:

\[ f(\omega_n,\omega_p)=\frac{r_0}{6\pi^{1/2}}\left[ \frac{\pi^{1/2}}{2}(\omega_n^3+\omega_p^3)\Phi\left(\frac{\omega_n+\omega_p}{2}\right)+ \right. \]
\[ +(\omega_n^2+\omega_p^2-\omega_n\omega_p-2)e^{-\frac{\omega_n^2+\omega_p^2}{2}}- \]
\[ -\frac{\pi^{1/2}}{2}(\omega_n^3-\omega_p^3)\Phi\left(\frac{\omega_n-\omega_p}{2}\right)- \]
\[ \left. -(\omega_n^2+\omega_p^2+\omega_n\omega_p-2)e^{-\left(\frac{\omega_n-\omega_p}{2}\right)^2} \right], \tag{84} \]

where \(\Phi(x)\) denotes the Gaussian transcendental function

\[ \Phi(x)=\frac{2}{\pi^{1/2}}\int_0^x e^{-u^2}\,du. \tag{85} \]

The derivation of these expressions, i.e. of the exchange energies for the cases of the interactions (49) and (50), can be carried out in a somewhat different manner than in the case of the Yukawa interaction considered by us. The calculation of \(\varepsilon_{ik}^{np}\) in these cases can be performed by direct integration of expression (60)\(^9\). For the interaction energy (49), for example, one can find \(W_A^{pp}\), \(W_A^{nn}\), and \(W_A^{np}\) also by differentiating, with respect to the parameter \(-1/r_0\), the expressions (71), (76), and (77), corresponding to the Yukawa interaction (47), i.e. to \(-2\pi p_0/h\), as is immediately evident from a comparison of (47) and (49).

The exchange energy \(W_R\) arising from the electrostatic Coulomb interaction of protons can be readily calculated for the case of free protons by means of the expression for \(\varepsilon_{ik}^{pp}\)*). The double integral in (54), which for brevity we denote by \(\varepsilon_{ik}^{R}\), for protons with parallel spins is then equal to \(2\varepsilon_{ik}^{pp}\), if in the expression for \(\varepsilon_{ik}^{pp}\) one substitutes, instead of the constants \(\gamma\) and \(1/r_0\), the values \(\gamma=e^2\) and \(1/r_0=0\), i.e. \(p_0=0\). In this case \(V(|\mathbf r-\mathbf r'|)\) is precisely transformed into the Coulomb interaction

*) Concerning the expression for \(\varepsilon_{ik}^{pp}\), compare, for example, § 9.

\(-\dfrac{e^2}{|\mathbf r-\mathbf r'|}\). Thus, for protons with parallel spins one obtains

\[ \varepsilon^R_{ik}=-\frac{e^2h^2}{\Omega\pi}\,\frac{1}{|\mathbf p_i-\mathbf p_k|^2}. \tag{86} \]

For protons with antiparallel spins, \(\varepsilon^R_{ik}=0\).

To calculate \(W_R\), it is necessary to sum \(\varepsilon^R_{ik}\) over all proton states with parallel spins. The summation may be replaced by integration, quite analogously to the calculation of \(W^p_A\); however, it must be noted that now the summation, i.e. the integration, extends only over proton states with parallel spins. We obtain\({}^{12}\)

\[ W_R=-\frac{1}{2}\,e^2\sum_{i=1}^{Q_p}\sum_{\substack{k=1\\ \sigma_k=\sigma_i}}^{Q_p}\varepsilon^R_{ik} = -\frac{4\pi e^2\Omega}{h^4}\,p_{\mu p}^4 . \tag{87} \]

If here as well, instead of \(p_{\mu p}\), one introduces the dimensionless function \(\omega_p\), then for a free proton gas it follows that

\[ W_R=-\frac{e^2\Omega}{r_0^4}\,\frac{1}{4\pi^3}\,\omega_p^4 . \tag{88} \]

For a nonuniform proton gas, i.e. for a proton gas with an arbitrary distribution density, we obtain, quite analogously to (79)—(81), the expression

\[ W_R=-\frac{e^2}{r_0^4}\,\frac{1}{4\pi^3}\int \omega_p^4\,dv . \tag{89} \]

Exactly as in (80) and (81), this expression contains the energy arising from the self-exchange of the protons, for the same reasons as there. Nevertheless, there is no need to correct the error thereby made in any special way; rather, one may proceed as follows. In the usual expression for the Coulomb energy \(W_C\), restore the terms corresponding to electrostatic self-interaction, which in any case cannot be avoided in the statistical method of treatment. This means that in expression (53) we shall also take into account the terms with \(k=i\), as a result of which, for an arbitrary, i.e. for both a free and a nonuniform, proton gas, we obtain the expression

\[ W_C=\frac{1}{2}e^2\sum_{i=1}^{Q_p}\sum_{k=1}^{Q_p} \iint |w_i(\mathbf r)|^2 |w_k(\mathbf r)|^2 \frac{1}{|\mathbf r-\mathbf r'|}\,dv\,dv' = \]

\[ =\frac{1}{2}e^2 \iint \frac{\rho_p(\mathbf r)\rho_p(\mathbf r')}{|\mathbf r-\mathbf r'|} \,dv\,dv', \tag{90} \]

where, according to the wave-mechanical definition of the particle density, there has been substituted

\[ \sum_{i=1}^{Q_p} |w_i(\mathbf r)|^2=\rho_p(\mathbf r). \tag{91} \]

The terms with \(k=i\) in \(W_C\), apart from the sign, are exactly the same as the terms with \(k=i\) in \(W_R\), as can easily be established by comparing (90) with (54). Since, moreover, these terms in \(W_C\) and \(W_R\) have different signs, it follows that in the sum \(W_C+W_R\)—and only such a sum will occur for us in what follows—they cancel completely. True, this cancellation takes place completely only in the wave-mechanical expressions, while in the statistical expressions it is only partial. This nevertheless need not trouble us, since these self-energies are significant only for small numbers of protons, and in that case, as can be shown, both \(W_C\) and \(W_R\) are small in comparison with \(W_A\).

Concerning expression (90) it should also be noted that \(W_C\) can also be written in the form

\[ W_C=\frac{1}{2}\,e\int V(\mathbf r)\rho_p(\mathbf r)\,dv, \tag{92} \]

where

\[ V(\mathbf r)=e\int \frac{\rho_p(\mathbf r')}{|\mathbf r-\mathbf r'|}\,dv' \tag{93} \]

denotes the electrostatic potential of the protons at the point \(\mathbf r\).

§ 6. The Saturation Property of Exchange Forces

As is known, the experimental data indicate that the energy of a nucleus is approximately proportional to the number of nucleons. This indicates that the attractive forces acting between nucleons, just like chemical valence forces, must possess a certain saturation, i.e., by means of these forces one nucleon can bind only a limited number of nucleons. Here we wish to show that the exchange forces described above, in contrast to ordinary forces, possess such a saturation property; moreover, we shall demonstrate this indirectly, by investigating the dependence of the exchange energy on the number of particles. In doing so we may put the correction factor \(s\) equal to unity.

Let us again consider the already discussed case of \(Q_n\) neutrons and \(Q_p\) protons, located in a volume \(\Omega\) and having a constant distribution density. In addition, we shall make one more simplification, assum—

assuming that the numbers of neutrons and protons are equal, i.e. \(Q_n=Q_p=Q/2\), where \(Q\) is the total number of nucleons. Then we shall have:

\[ \rho_n=\rho_p=\frac{1}{2}\rho=\frac{Q}{2\Omega} \quad\text{and}\quad \omega_n=\omega_p=\omega=\left(\frac{3\pi^2 r_0^3}{2}\right)^{1/3} \left(\frac{Q}{\Omega}\right)^{1/3}, \tag{94} \]

where \(\rho\) denotes the total density of nucleons.

For the exchange energy of the nucleons, from (73), (76), and (77), taking (94) into account, one obtains the expression

\[ W_A=W_A^{np}+W_A^{nn}+W_A^{pp} =-\frac{6\gamma}{r_0^4} f(\omega,\omega)\Omega = \]

\[ = -\frac{3}{8\pi}\frac{\gamma}{r_0}Q \left[ 6\omega-\frac{1}{\omega} +\frac{1}{4}\left(\frac{1}{\omega^3}+\frac{12}{\omega}\right) \times \right. \]

\[ \left. \times \ln(1+4\omega^2)-8\arctg 2\omega \right]. \tag{95} \]

As is seen from this expression, for large \(\omega\) we have:

\[ W_A \simeq -\frac{9}{4}\frac{\gamma}{r_0}Q\omega. \tag{96} \]

Thus, for large \(\omega\), \(W_A\) is proportional to \(Q\omega\). Since, furthermore, in view of the constant density, \(\omega\) is proportional to \(Q^{1/3}\), it follows that \(W_A\) is proportional to \(Q^{4/3}\).

Thus, whereas the ordinary energy of interaction of nucleons is proportional to \(r^2\), i.e. to \(Q^2\), \(W_A\) is proportional to \(Q\) only to the power \(4/3\). In this the saturation property of the exchange forces finds its expression. True, in this case the saturation is not complete, since with complete saturation of the interaction forces the interaction energy should have been proportional to the number of particles, i.e. proportional to \(Q\), whereas in our case proportionality to \(Q^{4/3}\) is obtained. This is a consequence of the fact that we represented the interaction energy between two nucleons in the form (47). If, however, we had taken the interaction energy in the form (49) or (50), then the exchange energy for \(\omega\gg 1\) would prove to be \(^{14}\) proportional to \(Q\), and the saturation in these cases for \(\omega\gg 1\) would thus be complete.

All this can easily be shown. Indeed, we need to obtain the expression for the total exchange energy \(W_A\) for the interactions (49) and (50), analogous to expression (95). With the aid of (83) and (84), and taking (94) into account, for the case of the interaction energy (49) one obtains

\[ W_A=-\frac{3}{8\pi}\varepsilon Q \left[ \frac{4}{\omega} -\left(\frac{1}{\omega^3}+\frac{6}{\omega}\right) \ln(1+4\omega^2)+8\arctg 2\omega \right] \tag{97} \]

and in the case of the interaction energy (50) one obtains

\[ W_A=-\frac{3}{2\pi^{1/2}}\varepsilon Q \left[ \pi^{1/2}\Phi(\omega) +\left(\frac{1}{\omega}-\frac{2}{\omega^3}\right)e^{-\omega^2} -\frac{3}{\omega}+\frac{2}{\omega^3} \right]. \tag{98} \]

For large \(\omega\), these expressions pass asymptotically, respectively, into

\[ W_A \simeq -\frac{3}{8\pi}\,\varepsilon Q\, 8\arctg 2\omega \simeq -\frac{3}{2}\,\varepsilon Q \tag{99} \]

and

\[ W_A \simeq -\frac{3}{2\pi^{1/2}}\,\varepsilon Q\,\tau^{1/2}\Phi(\omega) \simeq -\frac{3}{2}\,\varepsilon Q . \tag{100} \]

In both cases, therefore, \(W_A\) is proportional to \(Q\).

For the cases (95), (99), and (100), the quantity \(W_A/Q\) as a function of \(\omega\) is presented in Fig. 1, which makes the dependences quite clear.

Fig. 1. The exchange energy \(W_A/Q\) per one particle of a nucleon gas with constant density as a function of \(\omega\): curve 1 — \(W_A/Q\), calculated from (95) in units of \(\gamma/r_0\); curve 2 — \(W_A/Q\), calculated from (97) in units of \(\varepsilon\); curve 3 — \(W_A/Q\), calculated from (98) in units of \(\varepsilon\).

Thus, for large \(\omega\), there is a substantial difference between the interaction energies (49) and (50), on the one hand, and the Yukawa interaction energy (47) adopted by us as the basis, on the other hand. Partly on the basis of this difference, one may hope that the calculations carried out later in the present work will lead to results more satisfactory than those which could be obtained with the aid of the other interaction energies mentioned.

The reason for this different asymptotic behavior of the indicated interaction energies evidently lies in the fact that, for \(|\mathbf r-\mathbf r'|=0\), the interaction energies (49) and (50) pass into a certain constant, whereas the interaction energy (47) tends to infinity as \(1/|\mathbf r-\mathbf r'|\).

§ 7. Statistical model of the atomic nucleus

The statistical theory of the atomic nucleus is based on the assumption that the nucleons of the nucleus may be regarded as a degenerate neutron and proton gas at absolute zero temperature. It is assumed that in this nucleon gas the particles are distributed continuously, i.e., it is assumed that the particles are spread over the whole volume. The nucleon gas is held together compactly by means of the exchange energies of the kind discussed in the two preceding sections. The repulsive energies that balance the attraction are nothing other than the kinetic energies of the nucleon gas discussed in § 3 and the Coulomb interaction energy of the protons considered in § 4 and at the end of § 5. Accordingly, individual properties of nucleons are completely smoothed out in the statistical model of the nucleus. It follows further from the basic premises of the theory that it can be applied only to such nuclei in which the numbers of neutrons and protons are large and, consequently, the statistical method of treatment is justified.

In what follows we shall consider a nucleus consisting of \(N\) neutrons and \(Z\) protons. The total number of nucleons in the nucleus, i.e. the mass number, will be denoted by \(A\); thus we have:

\[ A=N+Z. \tag{101} \]

In consequence of the corrections to the Fermi kinetic energy and the exchange energy of identical particles already discussed above, we can consider not only heavy nuclei, but also take light nuclei into account.

The main problem consists in determining the energy of the nucleus and finding the density of the distribution of nucleons in the nucleus. The solution of this problem can be found with the aid of a certain variational principle. We shall arrive at this variational principle after we form expressions for the energy of the nucleus. To derive this expression we introduce a system of partitions, by means of which we divide the nucleon gas, in exactly the same way as in §§ 3 and 5, into partial volumes \(d\tau\), in such a way that each spatial volume element would still contain many neutrons and protons, which we shall regard as a free neutron or proton gas. We shall for the time being disregard the difficulty that this condition is not fulfilled in regions considerably removed from the center of the nucleus, in view of the small number of nucleons, since these regions

when calculating the energy of the nucleus, which is of primary interest to us, are of little importance. We can indicate the separate components of the nuclear energy at once with the aid of the results obtained in §§ 3 and 5. Before we do this, let us state several more general propositions that partly simplify the matter.

In calculating the energy it is convenient to express all the constituent parts of the energy in the same energy units. We shall put \(g=4e\) and choose as the unit of energy the quantity

\[ \varepsilon_0=\frac{g^2}{r_0}=2.724\cdot 10^{-5}\ \text{erg}=17.01\ \text{MeV}. \tag{102} \]

This is practically twice the average binding energy per nucleon of the nucleus.

Next, we shall simplify the expression for the Fermi kinetic energy and the expression for the exchange energy of identical particles by introducing the following simplifications for the correction factors: \(k_n, k_p, s_n\) and \(s_p\). These correction factors are significant only for light nuclei. In the most stable light nuclei \(N=Z=A/2\); therefore, in calculating the energies of the most stable light nuclei we may put

\[ k_n=k_p=k=1-\frac{4}{A}\quad \text{for } A\geq 4, \tag{103} \]

\[ k_n=k_p=k=0\quad \text{for } A<4, \tag{104} \]

and also

\[ s_n=s_p=s=1-\frac{2}{A}\quad \text{for } A\geq 2, \tag{105} \]

\[ s_n=s_p=s=0\quad \text{for } A<2. \tag{106} \]

These simplified correction factors \(k\) and \(s\) may also be retained for heavy nuclei, in which \(N\ne Z\), since for them the whole correction is insignificant. The cases \(A<4\) or \(A<2\) are indicated only for completeness; in practice they have no significance, since with such small numbers of particles the statistical theory becomes wholly inadmissible, and in our further calculations we shall consider only, as the most limiting case, the value \(A=4\). Therefore, in what follows we shall simply put

\[ k=1-\frac{4}{A}\quad \text{and}\quad s=1-\frac{2}{A}. \]

Finally, in the constituent parts of the energy, instead of the density functions \(\rho_n\) and \(\rho_p\), we shall introduce everywhere the dimensionless functions \(\omega_n\) and \(\omega_p\), which are defined by means of (72).

Taking into account all that has been said above, we shall now give the separate constituent parts of the energy. If we denote the Fermi kinetic energy of the nucleons by \(E_K\), the Weizsäcker correction to the kinetic

energy due to inhomogeneity by \(E_J\), the exchange energy arising from the neutron–proton interaction by \(E_A^{np}\), the exchange energy arising from the neutron–neutron interaction by \(E_A^{nn}\), the exchange energy arising from the proton–proton interaction by \(E_A^{pp}\), the ordinary electrostatic Coulomb energy of interaction of the protons by \(E_C\), and, finally, the exchange energy arising from this interaction by \(E_R\), then we shall have:

\[ E_K = k\varkappa \int \left(\rho_n^{5/3}+\rho_p^{5/3}\right)\,dv = \frac{k\mu_K}{(3\pi^2)^{5/3}r_0^5} \int \left(\omega_n^5+\omega_p^5\right)\,dv = \frac{k\mu_K}{(3\pi^2)^{5/3}r_0^5}\,\varepsilon_0 \int \left(\omega_n^5+\omega_p^5\right)\,dv, \tag{107} \]

\[ E_J=\chi_J\int\left[ \frac{1}{\rho_n}\left(\frac{d\rho_n}{dr}\right)^2 + \frac{1}{\rho_p}\left(\frac{d\rho_p}{dr}\right)^2 \right]\,dv = \frac{3\chi_J}{\pi^2 r_0^3} \int\left[ \omega_n\left(\frac{d\omega_n}{dr}\right)^2 + \omega_p\left(\frac{d\omega_p}{dr}\right)^2 \right]\,dv = \frac{3\mu_J}{\pi^2 r_0}\,\varepsilon_0 \int\left[ \omega_n\left(\frac{d\omega_n}{dr}\right)^2 + \omega_p\left(\frac{d\omega_p}{dr}\right)^2 \right]\,dv, \tag{108} \]

\[ E_A^{np} = -\frac{4\gamma}{r_0^4}\int f(\omega_n,\omega_p)\,dv = -\frac{4\lambda}{r_0^3}\,\varepsilon_0 \int f(\omega_n,\omega_p)\,dv, \tag{109} \]

\[ E_A^{nn} = -\frac{\gamma s}{r_0^4}\int f(\omega_n,\omega_n)\,dv = -\frac{\lambda s}{r_0^3}\,\varepsilon_0 \int f(\omega_n,\omega_n)\,dv, \tag{110} \]

\[ E_A^{pp} = -\frac{\gamma s}{r_0^4}\int f(\omega_p,\omega_p)\,dv = -\frac{\lambda s}{r_0^3}\,\varepsilon_0 \int f(\omega_p,\omega_p)\,dv, \tag{111} \]

\[ E_C = \frac{1}{2}e^2 \iint \frac{\rho_p(r)\rho_p(r')}{|r-r'|}\,dv\,dv' = \frac{e^2}{18\pi^4 r_0^6} \iint \frac{\omega_p^3(r)\omega_p^3(r')}{|r-r'|}\,dv\,dv' = \frac{1}{2}\frac{e}{3\pi^2 r_0^3} \int V(r)\omega_p^3(r)\,dv = \frac{1}{96\pi^2 e r_0^3}\,\varepsilon_0 \int V(r)\omega_p^3(r)\,dv, \tag{112} \]

\[ E_R = -\frac{e^2}{4\pi^3 r_0^4} \int \omega_p^4\,dv = -\frac{1}{64\pi^3 r_0^3}\,\varepsilon_0 \int \omega_p^4\,dv. \tag{113} \]

Here \(\mu_K\), \(\mu_J\), and \(\lambda\) denote the following dimensionless quantities:

\[ \mu_K=\frac{\varkappa}{r_0^2\varepsilon_0}=3.814;\qquad \mu_J=\frac{\chi_J}{r_0^2\varepsilon_0}=0.1660;\qquad \lambda=\frac{\gamma}{r_0\varepsilon_0}=\frac{1}{g^2}. \tag{114} \]

\(V(r)\) is the electrostatic potential of the protons, determined by means of (93); the integration is extended over the whole region where \(\rho_n\) and \(\rho_p\) do not vanish.

Instead of \(r_0\), a quantity \(\lambda\) has thus been introduced, which is for the time being an undetermined dimensionless parameter; moreover, it is the only arbitrary parameter in the theory, the value of which we shall determine from experimental data.

The total energy \(E\) of the nucleus is

\[ E=E_K+E_J+E_A^{np}+E_A^{nn}+E_A^{pp}+E_C+E_R . \tag{115} \]

The problem consists in determining such distribution functions \(\rho_n\) and \(\rho_p\), or \(\omega_n\) and \(\omega_p\), which give a minimum of the energy of the nucleus; moreover it is necessary to remember that for \(\rho_n\) and \(\rho_p\), or \(\omega_n\) and \(\omega_p\), there are additional conditions

\[ \int \rho_n\,d\tau=\frac{1}{3\pi^2 r_0^3}\int \omega_n^3\,d\tau=N, \tag{116} \]

\[ \int \rho_p\,d\tau=\frac{1}{3\pi^2 r_0^3}\int \omega_p^3\,d\tau=Z, \tag{117} \]

by means of which the total number of neutrons or protons in the nucleus is established. The problem, therefore, may be formulated in the form of the following variational principle:

\[ \delta\left(E+b_nN+b_pZ\right)=0, \tag{118} \]

where \(b_n\) and \(b_p\) denote Lagrange multipliers, and the variation must be carried out either with respect to \(\rho_n\) and \(\rho_p\), or, correspondingly, with respect to \(\omega_n\) and \(\omega_p\). Further, it should be noted that the electrostatic potential of the protons must satisfy the Poisson equation

\[ \Delta V=-4\pi e\rho_p=-\frac{4e}{3\pi r_0^3}\,\omega_p^3 . \tag{119} \]

It would not be difficult to derive, for finding \(\rho_n\) and \(\rho_p\), or \(\omega_n\) and \(\omega_p\), the corresponding system of differential equations with certain boundary conditions. However, solving this system of equations would lead to very lengthy and difficult numerical computations, which can be avoided by directly solving the variational problem by means of the Ritz method, as we shall do in the following paragraphs. In doing so, alongside conditions (116) and (117), it is sufficient to take into account also the following boundary conditions. First, from considerations of symmetry the following conditions must be satisfied:

\[ \left(\frac{\partial \rho_n}{\partial r}\right)_{r=0}=0 \quad \text{and} \quad \left(\frac{\partial \rho_p}{\partial r}\right)_{r=0}=0, \tag{120} \]

or the conditions

\[ \left(\frac{\partial \omega_n}{\partial r}\right)_{r=0}=0 \quad \text{and} \quad \left(\frac{\partial \omega_p}{\partial r}\right)_{r=0}=0, \tag{121} \]

where \(r\) denotes the distance from the center of the nucleus. Secondly, \(\rho_n\) and \(\rho_p\), or \(\omega_n\) and \(\omega_p\), must vanish at an infinite distance from the nucleus.

§ 8. A statistical nucleus with constant nucleon density

Before proceeding to the solution of the variational problem by the Ritz method, we shall first discuss the zero approximation, whose exact solution is an everywhere constant distribution density. This zero approximation is obtained if in the expression for the total energy (115) only the most important constituent parts of the energy are taken into account, namely \(E_K\) and \(E_A=E_A^{np}+E_A^{nn}+E_A^{pp}\). If very light nuclei are not considered, then, in comparison with these quantities, all the remaining components of the energy are only corrections, and in the zero approximation they may be neglected. Correspondingly, one may also set the correction factors \(k\) and \(s\) equal to unity. In addition, we shall make a further simplifying assumption, taking the numbers of neutrons and protons to be equal to one another, i.e., \(N=Z=A/2\). Then

\[ \rho_n=\rho_p=\frac{1}{2}\rho \quad\text{and}\quad \omega_n=\omega_p=\omega, \tag{122} \]

where \(\rho\) is the total density of nucleons.

For the energy of the nucleus, according to this simplified model, one obtains

\[ E= \frac{2\mu_K\varepsilon_0}{(3\pi^2)^{5/3}r_0^3} \int \omega^5\,dv - \frac{6\lambda\varepsilon_0}{r_0^3} \int f(\omega,\omega)\,dv, \tag{123} \]

where \(\omega\) must satisfy the additional condition

\[ \int \rho\,dv = \frac{2}{3\pi^2 r_0^3} \int \omega^3\,dv = A. \tag{124} \]

The variational principle now states

\[ \delta(E+bA)=0, \tag{125} \]

where \(b\) denotes the Lagrange multiplier.

From (125) one obtains the relation

\[ \frac{10\mu_K\varepsilon_0}{(3\pi^2)^{5/3}}\omega^4 - 6\lambda\varepsilon_0\frac{\partial f}{\partial \omega} + b\frac{2}{3\pi^2}\omega^2 = 0. \tag{126} \]

Taking into account that \(\frac{\partial f}{\partial\omega}\) is a function of only \(\omega\), while the Lagrange multiplier \(b\) is a constant, we find that \(\omega\), i.e. \(\rho\), is constant.

This result, which we have derived from the variational principle, could also have been obtained without any calculation. From expression (123) for the energy in the zero approximation it follows that the energy density, i.e. the energy per unit volume at each point, depends only on the density of particles at that very point and is completely independent of the distribution of density in the neighborhood.

mentioned point. Thus, in the nucleus there is no energetically distinguished point, whence it follows, in accordance with the result obtained above, that throughout the whole nucleus there exists everywhere one and the same density, namely that which leads to the least energy.

For a constant density of nucleons it follows from the normalization condition (124) that the density exists only inside a certain sphere, outside of which it vanishes everywhere. If the volume of this sphere is denoted by \(\Omega\), then from (124) one obtains the relation

\[ \rho=\frac{2}{3\pi^{2}r_0^3}\,\omega^3=\frac{A}{\Omega}. \tag{127} \]

Taking into account that \(\omega\) is constant, we obtain, by means of (127), from (123) for the energy of the nucleus

\[ E=\frac{2\mu_K\varepsilon_0}{(3\pi^2)^{5/3}r_0^3}\,\omega^5\Omega -\frac{6\lambda\varepsilon_0}{r_0^3}\,f(\omega,\omega)\,\Omega = \]

\[ = A\varepsilon_0\left\{ \frac{\mu_K}{(3\pi^2)^{2/3}}\omega^2 -\frac{3\lambda}{8\pi}\frac{1}{\omega^3} \left[ 6\omega^4-\omega^2 +\frac{1}{4}(1+12\omega^2)\ln(1+4\omega^2) -8\omega^3\operatorname{arctg}2\omega \right]\right\}. \tag{128} \]

The only arbitrary parameter in this expression is \(\lambda\). If one puts

\[ \lambda=3.31, \tag{129} \]

then for the energy minimum one obtains the relation

\[ E_0=-0.50\,A\varepsilon_0=-8.5A\ \text{MeV}, \tag{130} \]

which, except for the case of the very lightest nuclei, is on average well confirmed by experimental data.

The minimum of the energy lies at the following value of \(\omega\):

\[ \omega=2.39, \tag{131} \]

to which there corresponds the density

\[ \rho=0.922\,\frac{1}{r_0^3}. \tag{132} \]

For the limiting radius up to which the density extends, one obtains from (124), for \(\omega=\omega_0\),

\[ R=\frac{(9\pi)^{1/3}}{2\omega_0}\,A^{1/3}r_0 =0.637\,A^{1/3}r_0. \tag{133} \]

Thus, the limiting radius in this approximation, in accordance with the experimentally found nuclear radii, turns out to be proportional to \(A^{1/3}\), and also has the correct order of magnitude. To what extent it is permissible to compare this limiting radius directly with the experimental nuclear radii, we shall discuss at the end of § 9.

§ 9. Calculation of the distribution of the density and energy of a nucleus by Ritz’s method

We shall now proceed to the solution of the full variational problem (64) by Ritz’s method. Ritz’s method consists in assuming that the unknown functions—in our case the distribution densities \(\rho_n\) and \(\rho_p\), or \(\omega_n\) and \(\omega_p\)—have some form satisfying the boundary conditions and containing certain initially undetermined parameters. With the aid of this assumption the energy is calculated as a function of these parameters, while the parameters themselves are determined from the requirement that the energy be a minimum. Thus the whole method reduces to an ordinary minimum problem. For an exact determination of the functions sought it would be necessary to introduce infinitely many initially undetermined parameters, for example by expanding the function sought in a series in suitably chosen functions \(f_i\) with initially undetermined coefficients, i.e. by representing the function sought in the form \(\sum c_i f_i\). The functions \(f_i\) must satisfy the same boundary conditions as the function sought; otherwise they may be arbitrary. Since usually the behavior of the function sought is approximately known, in practice it is most expedient to choose such a formula of interpolating functions as approximates the exact function most closely. In this way it is possible to obtain a good approximation to the unknown function for the most part already with the aid of a small number of parameters. The integral whose minimum is sought—in our case the energy—will be the more rapidly approximated to the exact value, the closer the approximating functions are to the exact solution. The energy values obtained in approximate solutions by Ritz’s method naturally lie above the exact value of the energy.

Our problem possesses spherical symmetry, so that we have only one independent variable, namely the distance from the center of the nucleus \(r\). It proves expedient to introduce, instead of \(r\), the dimensionless variable

\[ x = \frac{1}{3^{1/2}}\,a\,\frac{r}{r_0}, \tag{134} \]

where \(a\) denotes a certain variational parameter*).

For \(\omega_n\) and \(\omega_p\) we shall make the following assumptions:

\[ \omega_n = -\omega_{n0} e^{-x^2}\sum_{i=0}^{m}\gamma_{ni}x^{2i}, \qquad \omega_p = -\omega_{p0} e^{-x^2}\sum_{i=0}^{m}\gamma_{pi}x^{2i}. \tag{135} \]

*) The factor \(\dfrac{1}{3^{1/2}}\) has been introduced in order to simplify the numerical coefficients.

where \(\omega_{n0}\) and \(\omega_{p0}\) are normalization factors, while the coefficients \(\gamma_{ni}\) and \(\gamma_{pi}\) denote further variational parameters; the coefficients \(\gamma_{n0}\) and \(\gamma_{p0}\) we shall set equal to unity. The normalization factors are determined by the normalization conditions (116) and (117). The variational parameters \(a\), \(\gamma_{ni}\), and \(\gamma_{pi}\) are determined from the requirement that the energy be a minimum. The indicated assumptions satisfy the boundary conditions, since, on the one hand, the functions (135), as is easy to see, satisfy the conditions (121), and, on the other hand, the factor \(e^{-x^2}\) in \(\omega_n\) and \(\omega_p\) guarantees that these functions vanish at infinity. For \(m=1, 2, 3, \ldots\), i.e., taking into account \(1, 2, 3, \ldots\) terms of the sums in the expressions (135), we shall obtain approximations of progressively increasing degree. We shall restrict ourselves here to the simplest case, when all the coefficients \(\gamma_{ni}\) and \(\gamma_{pi}\), with the exception of \(\gamma_{n0}\) and \(\gamma_{p0}\), are equal to zero. Then we have

\[ \omega_n=\omega_{n0}e^{-x^2},\qquad \omega_p=\omega_{p0}e^{-x^2}. \tag{136} \]

This means that we assume identical distributions of the densities of neutrons and protons and allow only a certain difference in the density functions, since we normalize them differently, in accordance with the different numbers of neutrons and protons.

The assumptions (136) contain only one variational parameter, namely the parameter \(a\), which, according to (134), is explicitly contained in \(x\) and, naturally, also enters into the normalization factors \(\omega_{n0}\) and \(\omega_{p0}\), for which, in view of (116) and (117), the relations hold

\[ \omega_{n0}=(9\pi)^{1/6}N^{1/3}a,\qquad \omega_{p0}=(9\pi)^{1/6}Z^{1/3}a. \tag{137} \]

In what follows, for shortening and simplifying certain formulas, it proves useful to introduce the following notation:

\[ c_n=2\omega_{n0}=2(9\pi)^{1/6}N^{1/3}a, \tag{138} \]

\[ c_p=2\omega_{p0}=2(9\pi)^{1/6}Z^{1/3}a, \tag{139} \]

\[ c=2(9\pi)^{1/6}\left(\frac{N+Z}{2}\right)^{1/3}a =2(9\pi)^{1/6}\left(\frac{A}{2}\right)^{1/3}a, \tag{140} \]

\[ c_1=\frac12(c_n+c_p), \tag{141} \]

\[ c_d=\frac12(c_n-c_p). \tag{142} \]

Next, it is expedient to express \(\omega_n\) and \(\omega_p\), as well as \(\rho_n\) and \(\rho_p\), through the parameters \(c_n\) and/or \(c_p\), and also to represent \(\rho_n\) and \(\rho_p\) as functions of \(r\). We have:

\[ \omega_n=\frac12 c_n e^{-x^2},\qquad \omega_p=\frac12 c_p e^{-x^2}, \tag{143} \]

\[ \rho_n=\frac{1}{3\pi^2r_0^3}\omega_n^3 =\frac{1}{24\pi^3r_0^3}c_n^3e^{-3x^2} =\frac{N}{\pi^{5/2}r_0^3}a^3e^{-a^2r^2/r_0^2}, \tag{144} \]

\[ \rho_p=\frac{1}{3\pi^2r_0^3}\omega_p^3 =\frac{1}{24\pi^3r_0^3}c_p^3e^{-3x^2} =\frac{Z}{\pi^{5/2}r_0^3}a^3e^{-a^2r^2/r_0^2}. \tag{145} \]

After all these preparations we shall proceed to the discussion of our original problem of calculating the energy as a function of the variational parameters; for this purpose it is convenient, instead of \(a\) as the variational parameter, to introduce the quantity \(c\), proportional to \(a\).

We start from formulas (107)—(113). After carrying out the elementary integration, the component parts of the energy \(E_K\), \(E_J\), \(E_C\), and \(E_R\) can be very simply expressed as functions of \(c\).

For \(E_K\) one obtains

\[ E_K = \frac{4\cdot 3^{1/2}\pi x_K k}{(3\pi^2)^{1/3} r_0^5}\, \frac{1}{a^3} \int_0^\infty \left(\omega_n^5+\omega_p^5\right)x^2\,dx = \]

\[ = \frac{3}{20}\left(\frac{3}{5}\right)^{1/2} \frac{\mu_K}{(3\pi^2)^{2/3}}\, 2^{2/3}\, \frac{N^{5/3}+Z^{5/3}}{A^{5/3}}\, kAc^2\varepsilon_0 . \tag{146} \]

The expression \(2^{2/3}(N^{5/3}+Z^{5/3})/A^{5/3}\), which for \(N=Z\) is equal to unity, can be simplified. If the neutron excess is denoted by \(n\), i.e. if one puts

\[ n=N-Z \tag{147} \]

and the expression just mentioned is expanded in a series in powers of \(n/A\), then one obtains

\[ 2^{2/3}\,\frac{N^{5/3}+Z^{5/3}}{A^{5/3}} = 1+ \frac{5}{9}\left(\frac{n}{A}\right)^2 + \frac{5}{243}\left(\frac{n}{A}\right)^4 +\cdots . \tag{148} \]

If this series is cut off after the second term and the expression thus obtained is substituted into (146), then we shall have:

\[ E_K = \frac{3}{20}\left(\frac{3}{5}\right)^{1/2} \frac{\mu_K}{(3\pi^2)^{2/3}} \left[ 1+\frac{5}{9}\left(\frac{n}{A}\right)^2 \right]kAc^2\varepsilon_0 = \]

\[ = 0.04630 \left[ 1+\frac{5}{9}\left(\frac{n}{A}\right)^2 \right] kAc^2\varepsilon_0 . \tag{149} \]

For \(E_J\) one obtains

\[ E_J = \frac{12\cdot 3^{1/2}x_J}{\pi r_0^2}\, \frac{1}{a} \int_0^\infty \left[ \omega_n\left(\frac{d\omega_n}{dx}\right)^2 + \omega_p\left(\frac{d\omega_p}{dx}\right)^2 \right]x^2dx = \]

\[ = \frac{3\mu_J}{(18\pi)^{1/3}}\, A^{1/3}c^2\varepsilon_0 = 0.1297\,A^{1/3}c^2\varepsilon_0 . \tag{150} \]

The component energy \(E_C\), after simple transformations, can be represented in the form

\[ E_C=-\frac{2e}{3\pi r_0^3}\int_0^\infty V(x)\omega_p^3(x)x^2\,dx =-\frac{e^2}{(2\pi)^{1/2}r_0}Z^2a \]
\[ =-\frac{1}{(32\pi)^{1/3}(72\pi^4)^{1/6}} \left(\frac{2Z}{A}\right)^2 A^{4/3}c\varepsilon_0, \tag{151} \]

where \(V(x)\) is the electrostatic potential, determined with the aid of (93), of the electric charge density
\[ e\rho_p=\frac{e}{3\pi^2r_0^3}\omega_p^3 \]
at the point \(x\). If the factor \(2Z/A\) is expressed, with the aid of (147) and the dependence \(A=N+Z\), in terms of \(n\) and \(A\), then one obtains

\[ \frac{2Z}{A}=1-\frac{n}{A}. \tag{152} \]

After substitution into (151) it follows that

\[ E_C=-\frac{1}{(32\pi)^{1/3}(72\pi^4)^{1/6}} \left[1-2\frac{n}{A}+\left(\frac{n}{A}\right)^2\right]A^{4/3}c\varepsilon_0 \]
\[ =-0.002250\left[1-2\frac{n}{A}+\left(\frac{n}{A}\right)^2\right]A^{4/3}c\varepsilon_0. \tag{153} \]

For \(E_R\) one obtains

\[ E_R=-\frac{3^{1/2}e^2}{\pi^2r_0}\frac{1}{a^3} \int_0^\infty \omega_p^4x^3\,dx =-\frac{27}{32\pi}\left(\frac{\pi}{3}\right)^{1/6} \frac{e^2}{r_0}Z^{4/3}a \]
\[ =-\frac{27}{2048\cdot 3^{1/2}\pi} \left(\frac{2Z}{A}\right)^{1/3}Ac\varepsilon_0. \tag{154} \]

The factor \((2Z/A)^{1/3}\) can be expanded in a series in powers of \(n/A\). Taking (152) into account, we obtain:

\[ \left(\frac{2Z}{A}\right)^{1/3} =1-\frac{4}{3}\frac{n}{A} +\frac{2}{9}\left(\frac{n}{A}\right)^2+\cdots \tag{155} \]

If this series is cut off after the third term and the expression thus obtained is substituted into (154), then it follows that

\[ E_R=-\frac{27}{2048\cdot 3^{1/2}\pi} \left[1-\frac{4}{3}\frac{n}{A} +\frac{2}{9}\left(\frac{n}{A}\right)^2\right]Ac\varepsilon_0 \]
\[ =-0.002423\left[1-\frac{4}{3}\frac{n}{A} +\frac{2}{9}\left(\frac{n}{A}\right)^2\right]Ac\varepsilon_0. \tag{156} \]

We must still calculate the exchange energies \(E_A^{np}\), \(E_A^{nn}\), and \(E_A^{pp}\) as functions of \(c\), which is somewhat more complicated. From (109), (110)

and (111), taking (78) into account, one obtains, in the notation (138)—(142),

\[ E_A^{np}=-16\cdot 3^{3/2}\pi\,\frac{\gamma}{r_0}\,\frac{1}{a^3} \int_0^\infty f(\omega_n,\omega_p)x^2\,dx = \]

\[ =-\frac{24}{\pi}\left(\frac{3}{\pi}\right)^{1/2}A\lambda\varepsilon_0\, \frac{1}{c^3}\left\{ \frac{3\pi^{1/2}}{512}(c_n^3c_p+c_nc_p^3) -\frac{1}{32}\left(\frac{\pi}{2}\right)^{1/2}c_nc_p \right. \]

\[ \left. +\frac{1}{4}\,[L(c_s)-L(c_d)] +\frac{3}{8}(c_n^2+c_p^2)[K(c_s)-K(c_d)] \right. \]

\[ \left. -\frac{3}{64}(c_n^4-2c_n^2c_p^2+c_p^4)[H(c_s)-H(c_d)] \right. \]

\[ \left. -\frac{1}{2}(c_n^3+c_p^3)M(c_s) +\frac{1}{2}(c_n^3-c_p^3)M(c_d) \right\}, \tag{157} \]

\[ E_A^{nn}=-4\cdot 3^{3/2}\pi\,\frac{\gamma}{r_0}\,\frac{s}{a^3} \int_0^\infty f(\omega_n,\omega_n)x^2\,dx = \]

\[ =-\frac{6}{\pi}\left(\frac{3}{\pi}\right)^{1/2}sA\lambda\varepsilon_0\, \frac{1}{c^3}\left[ \frac{6\pi^{1/2}}{512}c_n^4 -\frac{1}{32}\left(\frac{\pi}{2}\right)^{1/2}c_n^2 +\frac{1}{4}L(c_n) +\frac{3}{4}c_n^2K(c_n) -c_n^3M(c_n) \right], \tag{158} \]

\[ E_A^{pp}=-4\cdot 3^{3/2}\pi\,\frac{\gamma}{r_0}\,\frac{s}{a^3} \int_0^\infty f(\omega_p,\omega_p)x^2\,dx = \]

\[ =-\frac{6}{\pi}\left(\frac{3}{\pi}\right)^{1/2}sA\lambda\varepsilon_0\, \frac{1}{c^3}\left[ \frac{6\pi^{1/2}}{512}c_p^4 -\frac{1}{32}\left(\frac{\pi}{2}\right)^{1/2}c_p^2 +\frac{1}{4}L(c_p) +\frac{3}{8}c_p^2K(c_p) -c_p^3M(c_p) \right], \tag{159} \]

where \(L\), \(K\), \(H\), and \(M\) denote the following integrals:

\[ L(q)=\int_0^\infty x^2\ln(1+q^3e^{-2x^2})\,dx, \tag{160} \]

\[ K(q)=\int_0^\infty x^2e^{-2x^2}\ln(1+q^2e^{-2x^2})\,dx, \tag{161} \]

\[ H(q)=\int_0^\infty x^2e^{-4x^2}\ln(1+q^3e^{-2x^2})\,dx, \tag{162} \]

\[ M(q)=\int_0^\infty x^3e^{-3x^2}\operatorname{arctg}(1+qe^{-x^2})\,dx, \tag{163} \]

and the parameter \(q\) stands for \(c_s, c_d, c_n\), and \(c_p\).

For \(q>1\) these integrals cannot be evaluated in closed form, and for finding them one has to resort to numerical or graphical methods.

We have found, by a numerical method with high accuracy, the values of the integrals \(L(q)\), \(K(q)\), and \(M(q)\) for 18 different values of \(q\) in the interval \(1.3<q<11\), which is essential for our calculations. The values of the integral \(H(q)\), which, as will be seen, plays only a subordinate role in expression (157), were found by us for 8 different values of \(q\). The results are presented in Table IV of the Mathematical Appendix. For the integrals \(K(q)\), \(H(q)\), and \(M(q)\) one can give simple approximate formulas which, in the interval of values of \(q\) most important for the calculations, give a maximum error of \(0.7\%\). The integral \(L(q)\) can be expressed with the aid of a certain convenient approximate expression. We shall also discuss this in the appendix.

For \(q\leqslant 1\) these integrals can be represented in the form of the following series:

\[ L(q)=\frac{\pi^{1/2}}{2}\sum_{k=1}^{\infty}(-1)^{k-1}\frac{1}{(2k)^{5/2}}\,q^{2k}, \tag{164} \]

\[ K(q)=\frac{\pi^{1/2}}{4}\sum_{k=1}^{\infty}(-1)^{k-1}\frac{1}{k(2k+2)^{3/2}}\,q^{2k}, \tag{165} \]

\[ H(q)=\frac{\pi^{1/2}}{4}\sum_{k=1}^{\infty}(-1)^{k-1}\frac{1}{k(2k+4)^{1/2}}\,q^{2k}, \tag{166} \]

\[ M(q)=\frac{\pi^{1/2}}{4}\sum_{k=1}^{\infty}(-1)^{k-1}\frac{1}{(2k-1)(2k+2)^{1/2}}\,q^{2k+1}. \tag{167} \]

Now we shall return again to expressions (157)—(159). Our aim is to transform these expressions and to represent \(E_A^{np}\), \(E_A^{nn}\), and \(E_A^{pp}\) as functions of \(c\). This can be done step by step in the following way.

First of all we shall consider the quantity \(E_A^{np}\) and introduce into the expressions for \(E_A^{np}\), instead of \(c_n\) and \(c_p\), first the parameters \(c_s\) and \(c_d\). In view of (138), (139), (141), and (142), the relations

\[ c_n=c_s+c_d \quad \text{and} \quad c_p=c_s-c_d, \tag{168} \]

hold, with the aid of which the quantities occurring in (157) and indicated below...

the expressions can be transformed as follows:

\[ \begin{gathered} c_n^3 c_p + c_n c_p^3 = 2c_s^4 - 2c_d^4,\\ c_n c_p = c_s^2 - c_d^2,\\ c_n^2 + c_p^2 = 2c_s^2 + 2c_d^2,\\ c_n^4 + c_p^4 - 2c_n^2 c_p^2 = 16c_s^2 c_d^2,\\ c_n^3 + c_p^3 = 2c_s^3 + 6c_s c_d,\\ c_n^3 - c_p^3 = 6c_s^2 c_d + 2c_d^3. \end{gathered} \tag{169} \]

If, for brevity, we introduce the notation

\[ F(q)=-\frac{36}{\pi}\left(\frac{3}{\pi}\right)^{1/2}\frac{1}{q^3} \left[ \left(\frac{6\pi}{512}\right)^{1/2}q^4 -\frac{1}{32}\left(\frac{\pi}{2}\right)^{1/2}q^2 +\frac{1}{4}L(q)+\frac{3}{4}q^2K(q)-q^3M(q) \right], \tag{170} \]

\[ S(q)=\frac{72}{\pi}\left(\frac{3}{\pi}\right)^{1/2}\frac{1}{q^3} \left[ \frac{1}{4}q^4H(q) +q^3M(q)-\frac{1}{4}q^2K(q) \right], \tag{171} \]

then, taking (169) into account, one can write \(E_A^{np}\) in the form

\[ E_A^{np}=-\left[ \frac{2}{3}\left(\frac{c_s}{c}\right)^3F(c_s) -\frac{2}{3}\left(\frac{c_d}{c}\right)^3F(c_d) +\frac{c_s}{c}\left(\frac{c_d}{c}\right)^2S(c_s) +\frac{c_d}{c}\left(\frac{c_s}{c}\right)^2S(c_d) \right]A\lambda\varepsilon_0. \tag{172} \]

A substantial simplification of this expression can be made by taking into account the following two circumstances. First, among the parameters \(c\), \(c_s\), and \(c_d\) entering this expression, the parameter \(c_d\) is of a smaller order of magnitude than the parameters \(c_s\) and \(c\), which are of the same order of magnitude*); namely, everywhere \(c_d/c<1/10\); for \(N=Z\), \(c_d\) vanishes. Secondly, in the region \(c\simeq 6.5\), which is essential for the subsequent calculations, the parameter \(c_d\) is constantly less than unity. Taking all this into account, one may expand \(E_A^{np}\) in powers of \(c_d/c\) into a rapidly convergent series which, as it turns out, can be truncated without any appreciable error already after the term with \((c_d/c)^3\).

As is easy to verify, the principal term in (172) is the first term; the following three terms play the role of corrections**). We shall first expand precisely these correction terms in \(c_d/c\).

*) In this connection see equations (140)—(142), and also equations (182) and (183).

**) Concerning \(F\) and \(S\), see the expressions (187), (188) given below, taking (189) into account.

In the second term in (172) the expansion of \(F(c_d)\) reduces to the expansion of \(L(c_d)\), \(K(c_d)\), and \(M(c_d)\), which can be carried out with the aid of formulas (164), (165), and (167). It turns out that in \(F(c_d)\) the first non-vanishing term is of order \(c_d^2\). Taking this into account, we have that the second term in (172) is of order \(c_d^3(c_d/c)^3\) and, in comparison with the first term, is entirely negligible, in view of which it may be discarded.

The third term in (172), strictly speaking, need not be expanded in a series in \(c_d/c\), since \(c_d\) enters it only in the form of the factor \((c_d/c)^3\). It is possible, however—since we are dealing only with correction terms—to carry out certain simplifications, so that in this term \(c_s\) is replaced everywhere by \(c\). The difference between \(c_s\) and \(c\) is essentially proportional to \(c_d^2\)*; thus the distinction between \(c_s\) and \(c\) in this correction term will manifest itself only in a quantity proportional to \((c_d/c)^4\), which may be neglected. Thus we have:

\[ \frac{c_s}{c}\left(\frac{c_d}{c}\right)^2 S(c_s) = \left(\frac{c_d}{c}\right)^2 S(c). \tag{173} \]

In the fourth term in (172), one may again expand the integrals \(K(c_d)\), \(H(c_d)\), and \(M(c_d)\) with respect to \(c_d\) by means of formulas (165)—(170). Taking into account that here as well, as in the third term in (172), one may put \(c_s=c\), after simple calculations one obtains:

\[ \frac{c_d}{c}\left(\frac{c_s}{c}\right)^2 S(c_d) = -\tau\left(\frac{c_d}{c}\right)^2, \tag{174} \]

where

\[ \tau = \frac{72}{\pi}\left(\frac{3}{\pi}\right)^{1/2} - \frac{3\pi^{1/2}}{128} = \frac{27\cdot 3^{1/2}}{16\pi} = 0.93037. \tag{175} \]

The next non-vanishing term of the expansion of the left-hand side of (174) would be a quantity proportional to \((c_d/c)^4\), and consequently negligibly small in comparison with the term standing on the right-hand side of (174).

Taking all these results into account, we thus obtain from (172):

\[ E_A^{np} = -\frac{2}{3}\left(\frac{c_s}{c}\right)^3 F(c_s)A\lambda\varepsilon_0 + [S(c)-\tau c]\left(\frac{c_d}{c}\right)^2 A\lambda\varepsilon_0 . \tag{176} \]

Before transforming in this expression \(F(c_s)\) from the variable \(c_s\) to the variable \(c\), it is expedient also to take into account the components of the energy \(E_A^{nn}\) and \(E_A^{pp}\), and only then carry out the indicated transformation in the general expression for the total exchange energy.

The expressions for the energies \(E_A^{nn}\) and \(E_A^{pp}\) defined in (158) and (157) can be written in the notation (170) in the following form:

\[ E_A^{nn} = -\frac{s}{6}\left(\frac{c_n}{c}\right)^3 F(c_n)A\lambda\varepsilon_0, \]

\[ E_A^{pp} = -\frac{s}{6}\left(\frac{c_p}{c}\right)^3 F(c_p)A\lambda\varepsilon_0 . \tag{177} \]

\[ \text{*) We shall speak of this below, cf. (182) and (183).} \]

Using relations (168), one can transform these expressions from the variables \(c_n\) and \(c_p\) to the variable \(c_d\), and then expand \(E_A^{nn}\) and \(E_A^{pp}\) in a series in powers of \(c_d/c\), which we again cut off after the terms with \((c_d/c)^3\). After simple calculations it follows that

\[ E_A^{nn}+E_A^{pp} = -\frac{s}{3}\left(\frac{c_s}{c}\right)^3 F(c_s)A\lambda\varepsilon_0 - \left[ F(c)+cF'(c)+\frac{1}{6}c^3F''(c) \right] \left(\frac{c_d}{c}\right)^2 A\lambda\varepsilon_0 . \tag{178} \]

Here again the first term is the most essential, while the second plays the role of a correction; in it, just as in (173) and (174), \(c_s=c\) is set; \(F'\) and \(F''\) denote the first and second derivatives of \(F(c)\) with respect to \(c\). In the correction term the correction factor \(s\) is set equal to unity.

With the aid of (176) and (178), we can now represent the total exchange energy in the following form:

\[ E_A = E_A^{np}+E_A^{nn}+E_A^{pp} = -\left(\frac{2}{3}+\frac{s}{3}\right) \left(\frac{c_s}{c}\right)^3 F(c_s)A\lambda\varepsilon_0 + \]

\[ +\left\{ S(c)-\tau c - \left[ F(c)+cF'(c)+\frac{1}{6}c^3F''(c) \right] \right\} \left(\frac{c_d}{c}\right)^2 A\lambda\varepsilon_0 . \tag{179} \]

It remains for us to transform the first term from the variable \(c_s\) to the variable \(c\).

Before taking this last step and, with the aid of the identity \(c_s=c+(c_s-c)\), expanding the first term on the right-hand side of (179) in a series in the small quantity \((c_s-c)/c\) and simultaneously introducing the variable \(c\) instead of \(c_s\), it is expedient to express \(c_d/c\) and \((c_s-c)/c\) in terms of the neutron excess \(n=N-Z\) and the mass number \(A=N+Z\). This can be done very simply.

For this purpose we express \(N\) and \(Z\) in terms of \(n\) and \(A\); we have:

\[ N=\frac{1}{2}(A+n), \qquad Z=\frac{1}{2}(A-n). \tag{180} \]

If these expressions are substituted into (138) or (139), \(N^{1/3}\) or \(Z^{1/3}\) is expanded in a series in the quantity \(n/A\), and the series is cut off after the term linear in \(n/A\), then for \(c_d/c\) one obtains

\[ \frac{c_d}{c} = \frac{c_n-c_p}{2c} = \frac{1}{3}\frac{n}{A}. \tag{181} \]

The next nonvanishing term of the expansion would be of order \((n/A)^3\), and it may be neglected.

To calculate \((c_s-c)/c\), we analogously expand \(c_s\) in a series in the quantity \(n/A\). If the series is cut off after the term with \((n/A)^2\),

then we obtain:

\[ c_s=-\frac{1}{2}(c_n+c_p)=c-c\frac{1}{9}\left(\frac{n}{A}\right)^2, \tag{182} \]

\[ \frac{c_s-c}{c}=-\frac{1}{9}\left(\frac{n}{A}\right)^2=-\left(\frac{c_d}{c}\right)^2. \tag{183} \]

The next nonvanishing term of the expansion is proportional to \((n/A)^4\) and may be discarded. Thus, whereas \(c_d/c\) is proportional to \(n/A\), the quantity \((c_s-c)/c\) proves to be proportional to \((n/A)^2\).

Taking this into account, and also our assumption that in the expression for the energy we neglect terms of order smaller than \((c_d/c)^2 \sim \frac{1}{10}(n/A)^2\), one may cut off the expansion of the expression for the energy (179) in \(c_s-c\) after the term linear with respect to \(c_s-c\). Thus, the part of the first term in (179) that depends on \(c_s\) may be represented in the following form:

\[ \left(\frac{c_s}{c}\right)^3 F(c_s) = F(c)-\frac{1}{3}F(c)\left(\frac{n}{A}\right)^2 -\frac{1}{9}cF'(c)\left(\frac{n}{A}\right)^2. \tag{184} \]

In this expression the most essential term is the first summand, while the terms containing \((n/A)^2\) play the role of a correction.

When this expression is substituted into (179), in the correction term one may again set the correction factor \(s\) equal to unity. If, further, in the relation thus obtained one expresses \(c_d/c\) in terms of \(n/A\) by means of (181), then for \(E_A\), as a function of \(c\), one finally obtains

\[ E_A = -\left(\frac{2}{3}+\frac{s}{3}\right)F(c)A\lambda\varepsilon_0 +\frac{1}{9}\left[ S(c)-\pi c+2F(c) -\frac{1}{6}c^2F''(c) \right] \left(\frac{n}{A}\right)^2 A\lambda\varepsilon_0. \tag{185} \]

The function \(F(c)\) required for knowing this energy is given in the Mathematical Appendix, with accuracy to four figures, for 18 different values of \(c\) (see Table IV).

With the aid of some orienting calculations it is easy to establish that all values of the quantity \(c\), determined from the requirement that the energy be a minimum, lie in the interval

\[ 5.4<c<7.5 \tag{186} \]

for nuclei beginning with the lightest and ending with the heaviest; thus, in what follows we may restrict ourselves to values of \(c\) from this interval. To simplify the calculations, one may proceed so as to approximate in this interval the functions \(F(c)\) and \(S(c)\), which are contained in the integrals (160)—(163), by means of simpler functions. The approximation of \(F(c)\) should be carried out with particular care, for \(F(c)\) determines the principal part \(E_A\). It turns out that precisely this approximation is made

most easily, in view of the fact that in the interval (186) \(F(c)\) varies approximately linearly. Concerning this see Fig. 2, in which \(F(c)\) is represented graphically. We find that both functions can be represented in the following form:

\[ F(c)=a_1c^2+a_2c+a_3, \tag{187} \]

\[ S(c)=b_1c^3+b_2c, \tag{188} \]

and thereby we do not introduce any noticeable error into the results.

Fig. 2. Behavior of the function \(F(c)\).

For the coefficients, by means of the method of least squares, the following values are obtained:

\[ \begin{aligned} a_1&=0.006238, & b_1&=0.05375,\\ a_2&=0.09723, & b_2&=0.93456,\\ a_3&=-0.1858, \end{aligned} \tag{189} \]

With such coefficients both polynomials in the interval (186) approximate the exact functions excellently*), the maximum error

*) One could also proceed as follows: for the integrals \(L, K, M\) in \(F(c)\), substitute the approximate expressions given in the Appendix. However, the approximate expression thereby obtained for \(F(c)\) is less accurate than expression (187).

in both cases less than \(0.2\%\). Let us also mention that the correction term in \(E_A\) contains the second derivative of \(F(c)\), and for it, naturally, the approximation is less accurate than for \(F(c)\). However, this is of no significance, since we are dealing only with a correction term.

With the aid of functions (187) and (188), after substituting the value \(s=1-\dfrac{2}{A}\), for the correction factor one obtains \(E_A\) in interval (186) in the following final form:

\[ \left. \begin{aligned} E_A &=-\left\{\left(1-\frac{2}{3A}\right)(a_1c^2+a_2c+a_3) -\frac{1}{9}\left[\left(\frac{5}{3}a_1+b_1\right)c^3\right.\right.\\ &\qquad\left.\left.+(2a_2+b_2-\tau)c+2a_3\right]\left(\frac{n}{A}\right)^2\right\}A\lambda\varepsilon_0\\ &=-\left[\left(1-\frac{2}{3A}\right)(0.006238c^2+0.09723c-0.1858)\right.\\ &\qquad\left.-(0.007127c^2+0.02208c-0.04129)\left(\frac{n}{A}\right)^2\right]A\lambda\varepsilon_0 . \end{aligned} \right\} \tag{190} \]

After summing the expressions for the energies (149), (150), (153), (156), and (190), the total energy of the nucleus \(E\) can be represented as a function of \(c\) in the following form:

\[ E=E_K+E_J+E_A+E_C+E_R=(P_1c^2-P_2c+P_3)A\varepsilon_0, \tag{191} \]

where \(P_1\), \(P_2\), and \(P_3\) depend only on the mass number \(A\) and the ratio \(n/A\), and have the following form:

\[ P_1=\alpha_0+\alpha_2\left(\frac{n}{A}\right)^2;\qquad P_2=\beta_0+\beta_1\frac{n}{A}+\beta_3\left(\frac{n}{A}\right)^3; \]

\[ P_3=\gamma_0+\gamma_2\left(\frac{n}{A}\right)^3, \tag{192} \]

\[ \left. \begin{aligned} \alpha_0&=0.04630-0.006238\lambda+0.1297\frac{1}{A^{1/3}} -(0.1852-0.004159\lambda)\frac{1}{A},\\ \alpha_2&=0.02572+0.007127\lambda-0.1029\frac{1}{A},\\ \beta_0&=0.09723\lambda+0.002423 -0.002250A^{2/3}-0.06482\lambda\frac{1}{A},\\ \beta_1&=-0.003231+0.004500A^{2/3},\\ \beta_2&=-0.02208\lambda+0.00054-0.002250A^{2/3},\\ \gamma_0&=0.1858\lambda-0.1239\lambda\frac{1}{A},\\ \gamma_2&=0.04129\lambda . \end{aligned} \right\} \tag{193} \]

where, let us note, for the correction factor \(k\) the value \(k=1-\dfrac{4}{A}\) has been substituted. The quantities \(a_0, a_2, \beta_0, \beta_1, \beta_2, \gamma_0\), and \(\gamma_2\) are thus functions only of \(A\). In the expressions for \(P_1, P_2\), and \(P_3\), the terms containing \(n/A\) or \((n/A)^2\) again represent corrections that vanish when the neutron excess is set equal to zero.

The determination of the variational parameter \(c\) is carried out by means of the requirement that the energy \(E\) be a minimum, i.e., from the equation

\[ \frac{dE}{dc}=2P_1c-P_2=0, \tag{194} \]

from which it follows that

\[ c_0=\frac{P_2}{2P_1}. \tag{195} \]

With the aid of this value of \(c\), for the minimum of \(E\) one obtains the expression

\[ E_0=-\left(\frac{P_2^2}{4P_1}-P_3\right)A\varepsilon_0, \tag{196} \]

which, like \(c_0\), is a function of \(A\) and \(n/A\).

§ 10. Results relating to the energy and the density distribution of the nucleus

The energy of the nucleus is given by the expression for \(E_0\). After substituting (192) and (193) into (196), \(E_0\) is obtained in the form of a function of the mass number \(A=N+Z\) and the neutron excess \(n=N-Z\). For a given \(A\), therefore, the energy of the nucleus is also a function of \(n\), which gives the energy of isobars with mass number \(A\).

The sequence of the most stable nuclei, i.e., those values of \(n\) or \(Z\) which give the minimum energy for a given value of \(A\), is obtained from the condition

\[ \left(\frac{\partial E_0}{\partial n}\right)_A=0. \tag{197} \]

To compute from this \(Z\) as a function of \(A\), we expand \(E_0\) in a series in \(n/A\) and cut it off after the term with \((n/A)^3\). We shall have

\[ E_0=-\left[ \frac{1}{4}\frac{\beta_0^2}{a_0}-\gamma_0 +\frac{1}{2}\frac{\beta_0\beta_1}{a_0}\left(\frac{n}{A}\right) + \frac{1}{4}\frac{1}{a_0} \left( \beta_1^2+2\beta_0\beta_2-\frac{a_2}{a_0}\beta_0^2-4a_0\gamma_2 \right) \left(\frac{n}{A}\right)^2 \right]A\varepsilon_0. \tag{198} \]

Table I

Comparison of the calculated values of nuclear charge and energy with experimental data for nuclei possessing the lowest energy. Energies are expressed in MeV*)

\(A\)	Calculated \(c_{0n}\)	Calculated \(Z_{\mathrm{eff}}\)	Calculated \(Z\)	Calculated \(-E_{00}/A\)	Experimental \(Z\)	Experimental \(-E_{00}/A\)
4	5,546	1,97	2	4,769	2	7,05
5	5,254	2,45	2	4,098	(2)	5,48
6	5,556	2,94	3	5,030	3	5,32
7	5,470	3,42	3	4,891	3	5,59
8	5,680	3,91	4	5,516	(4)	7,039
9	5,660\(_5\)	4,39	4	5,542	4	6,442
10	5,814	4,87	5	5,972	5	6,443
11	5,822	5,35	5	6,072	5	6,709
12	5,938	5,83	6	6,372	6	7,515
13	5,959	6,30	6	6,506	6	7,435
14	6,050	6,78	7	6,718	7	7,477
15	6,077	7,26	7	6,871	7	7,671
16	6,149	7,73	8	7,015	8	7,948
17	6,180	8,21	8	7,177	8	7,724
20	6,316	9,62	10	7,498	10	7,999
21	6,350	10,09	10	7,663	10	7,933
28	6,540	13,37	13	8,203	14	8,426
40	6,803	18,89	19	8,788	18	8,556
41	6,796	19,35	19	8,804	19	8,551
52	6,928	24,33	24	9,040	24	8,830
60	7,010	27,90	28	9,126	28	8,752
64	7,042	29,67	30	9,141	28	8,760
80	7,113	36,67	37	9,130	34; 36	—
96	7,150	43,52	44	9,119	40; 42; 44	—
100	7,134	45,21	45	8,984	42; 44	—
110	7,132	49,40	49	8,873	46; 48	—
120	7,157	53,33	54	8,759	50; 52	—
125	7,153	55,58	56	8,697	52	—
140	7,135	61,64	62	8,494	58	—
160	7,105	69,52	70	8,209	64; 66	8,202
170	7,059	73,39	73	8,051	68; 70	—
180	7,042	77,19	77	7,907	72; 74	—
200	7,005	84,63	85	7,610	80	—
209	6,972	87,90	88	7,473	83	7,799
215	6,951	90,05	90	7,381	(85)	7,719
220	6,941	91,83	92	7,308	(86)	7,669
240	6,879	98,78	99	7,014	(96)	7,496
242	6,854	99,46	99	6,973	(96)	7,488

*) In those cases where, owing to the absence of experimental data for the given \(A\), the nucleus with the lowest energy has not been established, the charge values are given for all stable isotopes; the charge values of the nuclei of unstable isotopes are placed in parentheses.

Substituting (198) into (197), differentiating with respect to \(n\) and then replacing, according to the definition, \(n\) by \(N-Z\), we obtain for the atomic number of the most stable nuclei the expression

\[ Z_{\mathrm{eff}}=\left(1+\frac{a_0\beta_0\beta_1}{a_0\beta_1^2+2a_0\beta_0\beta_2-a_2\beta_0^2-4a_0^2\gamma_2}\right)\frac{A}{2}. \tag{199} \]

Since the atomic number calculated in this way is not an integer, as the true atomic number of the most stable nucleus with mass number \(A\) we shall take that integer which lies closest to \(Z_{\mathrm{eff}}\), and as the neutron excess we shall take the value obtained from this integer \(Z\) by means of the relation \(n=A-2Z\).

Taking this into account, we can find the atomic number \(Z\), or the neutron excess \(n=A-2Z\), for the most stable nuclei for any value of \(A\). The energy \(E_{00}\) of the most stable nuclei is obtained if, for each prescribed \(A\), one calculates from (196) the energy by using \(n\), rounded in the indicated manner to an integer.

Of course, these calculations can be carried out only when the constant \(\lambda\), which has so far not yet been determined, is established. We shall choose it so that the calculated energies of the most stable nuclei over the entire interval from the lightest to the heaviest nuclei deviate as little as possible from the corresponding experimental values. We set

\[ \lambda=4.140; \tag{200} \]

for this value our requirement is satisfied quite satisfactorily. The fact that the value of \(\lambda\) obtained here is larger than the corresponding value obtained in § 8 for the model of a nucleus with constant density is a consequence of the fact that there we neglected the repulsion energies \(E_J\) and \(E_C\). In the exact calculation carried out here, however, they were taken into account, as a result of which they had to be compensated by increasing the energy of attraction, i.e. by increasing the value of \(\lambda\).

Table I gives the mean calculated energies \(E_{00}/A\), per particle of the nucleus that has the lowest energy for the given \(A\) (the most stable isobars), as well as the corresponding experimental values*; next, the values of \(c_0\) for these nuclei, denoted by us as \(c_{00}\), are given. In addition, \(E_{00}/A\) as a function of \(A\) is shown graphically in Fig. 3, where the experimental values are also given. As the mean experimental energy of the nucleus for a prescribed \(A\) we have taken

\[ \text{* Relative to the experimental values, cf. }{}^{11}. \]

the mean energy per particle of the stable nucleus with the lowest energy. For those values of \(A\) for which stable nuclei do not exist, we indicated the mean energy per particle of the unstable nucleus with the lowest energy. As is seen from the comparison, the calculated energies agree very well with the experimental values. For medium and heavy nuclei the maximum deviation of the calculated energies from the experimental ones is everywhere less than \(7\%\). For light nuclei, statistical theory can, of course, give only a certain

Fig. 3. Mean energy per particle \(E_{00}/A\) as a function of \(A\) for nuclei having, for a given \(A\), the lowest energy: line — calculated values, ● — experimental values for stable nuclei, ○ — experimental values for unstable nuclei.

average of strongly fluctuating empirical values, which in fact is quite satisfactory for the theory. The above-mentioned rounding of the values of \(Z_{\mathrm{eff}}\), calculated from (199), leads to the fact that for very light nuclei, where the difference in energy of individual isobars is relatively large, jumps are observed in the behavior of the energy. In fact, the mean energies per particle \(E_{00}/A\) of the most stable isobars with successive values of \(A\) do not lie on any smooth curve; on the contrary, for even mass numbers these energies lie somewhat lower, and for odd ones somewhat higher, than would correspond to a smooth behavior.

The reason why the theory also leads to such satisfactory results for light nuclei lies, basically, in the fact that, on the one hand, we have corrected the expression for the kinetic energy, by means of which, for \(A \leq 4\), it passes into the exact wave-mechanical expression, and, on the other hand, we

they introduced a correction also to the exchange energy of identical particles, thereby eliminating the self-exchange of particles.

In discussing the results obtained it is necessary to emphasize that they were achieved with the aid of only one single empirical parameter—the parameter \(\lambda\).

The formula for the energy (196), in view of (193), is valid only for \(A \geq 4\). It is now possible, mainly for orientation, to extend the calculations also to the case of nuclei with \(A < 4\), and it must be taken into account that for \(A < 4\) the correction factor \(k\) determined by formulas (103) and (104) must be set equal to zero. It is noteworthy that these calculations in no way lead to meaningless results, which may be attributed to the two corrections mentioned above.

In Table I we have also given the values of \(Z_{\mathrm{eff}}\), calculated with the aid of (199), for the most stable nuclei, which, as was already mentioned, are in general not integers. The theoretical values of \(Z\), i.e. the integers lying closest to \(Z_{\mathrm{eff}}\), are indicated in the fourth column of the table. For comparison the table also gives the experimental atomic numbers\({}^{11}\) of the most stable isobars. In those cases where, owing to the absence of experimental data, the most stable isobar has not been established, the atomic numbers of all stable isobars are indicated. Comparison of the theoretical values of \(Z\) with the experimental ones leads us, in this case as well, to a very satisfactory result—the maximum deviation of the calculated atomic numbers from the empirical ones is less than \(8\%\).

In addition to the results indicated here, the magnitudes of the separate component parts of which the energy of the nucleus is composed are also of interest. They are given in Table II (see p. 436), again for the case of the most stable isobars. As can be seen from these data, with the exception of the lightest nuclei, the most important components of the energy are \(E_A\) and \(E_K\), while the components \(E_J\), \(E_C\) and especially \(E_R\), in comparison with them, have considerably smaller weight. The exception is only the very lightest nuclei, for which, along with \(E_A\), not \(E_K\) but \(E_J\) is the most important term. It should also be noted that, of the energy components \(E_J\) and \(E_C\), the first has a relatively large value for light nuclei, and the second—for heavy nuclei, as is also directly evident from formulas (150) and (151), according to which, in view of the practical constancy of the values of \(c\), the energy \(E_J\) is proportional to \(A^{1/3}\), while the energy \(E_C\) is proportional to \(A^{5/3}\).

The calculations that we have presented here for the case of the most stable isobars can, of course, be extended to any other nuclei as well. For example, for a sequence of isobars with constant \(A\) one can calculate the energy for various prescribed values of \(Z\). We have done this for the cases of small, medium, and large

Table II

Total energy and component parts of the energy for various nuclei. All energies are expressed in MeV

\(A\)	\(Z\)	\(-E_A\)	\(E_K\)	\(E_J\)	\(E_C\)	\(-E_R\)	\(-E_0\)
4	2	128,012	0	107,726	2,140	0,914	19,060
5	2	147,935	22,227	104,155	1,882	0,804	20,475
6	3	205,457	48,629	123,770	4,213	1,374	30,219
7	3	234,623	71,490	126,266	3,940	1,285	34,212
8	4	293,192	101,638	142,359	6,937	1,873	44,111
9	4	328,799	127,038	147,039	6,666	1,794	49,850
10	5	388,030	159,736	160,673	10,328	2,396	59,689
12	6	488,073	222,170	178,107	14,295	2,937	76,438
16	8	699,63	357,33	210,20	23,91	4,05	112,24
20	10	922,27	502,90	239,00	35,63	5,21	149,95
40	19	2113,87	1313,80	319,14	109,93	10,47	351,47
60	28	3313,73	2172,66	424,43	214,92	15,81	547,53
80	37	4573,39	3037,24	480,87	345,95	21,14	730,47
100	45	5739,13	3868,94	521,12	476,45	25,78	898,40
120	54	6930,5	4705,1	557,4	647,7	30,8	1051,1
140	62	8038,4	5492,4	583,2	808,6	35,0	1189,2
160	70	9116,0	6255,7	604,6	981,7	39,2	1313,2
180	77	10083,0	6953,0	617,7	1132,0	42,4	1422,7
200	85	11102,7	7668,8	633,1	1324,9	46,5	1522,4
220	92	12006,3	8316,8	641,8	1489,8	49,6	1607,7
240	99	12883,1	8943,8	648,7	1660,8	52,6	1682,7

of the mass number, namely for \(A=16\), \(A=80\), and \(A=200\); the corresponding energy curves are presented in Figs. 4, 5, and 6. As is seen from a comparison of these figures, the energy curve for heavy isobars has a considerably flatter form than for light ones. The energy minimum for the isobar \(A=16\) lies at \(Z=8\), for the isobar \(A=80\) at \(Z=37\), and for the isobar \(A=200\) at \(Z=86\). The fact that the last value of \(Z\) deviates by approximately \(1\%\) from the value of \(Z\) computed from formula (199) (see Table I) is a consequence of the use, in deriving formula (199), of the expansion (198) of the energy expression (196), whereas the calculations carried out here were based on the general energy expression (196).

We now turn to the discussion of the density of distribution of nucleons in the nucleus. With the aid of equations (144) and (145), as well as (138)—(140), the distribution density of all nucleons (neutrons + protons) can be represented in the following form:

\[ \rho=\rho_n+\rho_p=\rho_0 e^{-3x^2}=\rho_0 e^{-a^2(r/r_0)^2}, \tag{201} \]

where

\[ \rho_0=\frac{c_n^3+c_p^3}{24\pi^2 r_0^3}=\frac{c^3}{12\pi^2 r_0^3}. \tag{202} \]

Fig. 4. Average energy per particle \(E_0/A\) as a function of \(Z\) for the isobar \(A = 16\).

Fig. 5. Average energy per particle \(E_0/A\) as a function of \(Z\) for the isobar \(A = 80\).

Fig. 6. Average energy per particle \(E_0/A\) as a function of \(Z\) for the isobar \(A = 200\).

denotes the density of nucleons at the center of the nucleus. For the parameters \(c\) and \(a\) we shall take the corresponding values \(c_0\) and \(a_0\) that minimize the energy; moreover, between \(c_0\) and \(a_0\), according to (140), there is the relation

\[ a_0=\frac{2^{1/3}}{2(9\pi)^{1/6}}\,\frac{c_0}{A^{1/3}} . \tag{203} \]

According to (202), \(\rho_0\) is proportional to \(c^3\), with a universal proportionality factor independent of \(A\).

The values of \(c_{00}\) for several of the most stable isobars are given in Table I. As may be seen, \(c_{00}\) is practically independent of \(A\) for \(A>50\), i.e., for medium and heavy nuclei. For these nuclei \(c_{00}\) has a maximum relative deviation from the value \(c_{00}\simeq 7.0\) of order \(\pm 2\%\). As a consequence of this behavior of \(c_{00}\) for \(A>50\), it follows that for these nuclei the density of nucleons

Fig. 7. Nucleon density \(\rho\) as a function of \(r\) for nuclei \(A=16,\ Z=8;\ A=80,\ Z=37;\ A=200,\ Z=83\). The numbers standing near the curves are the mass numbers of the corresponding nuclei.

at the center of the nucleus, \(\rho_0\), is practically constant, having a maximum relative deviation from the value

\[ \rho_0=\frac{73}{12\pi^2 r_0^3}=2.90\,\frac{1}{r_0^3}=1.16\cdot 10^{39}\ {\rm cm}^{-3} \tag{204} \]

of order \(\pm 6\%\). In the region \(A>50\), \(\rho_0\) has a very flat maximum at \(A\simeq 120\). From \(A=50\) to \(A=4\), the dependence of the nucleon density at the center of the nucleus on the mass number appears much more strongly; namely, \(\rho_0\) decreases monotonically in the direction of smaller mass numbers and at \(A=4\) reaches approximately half the value (204). For all this, see Table III, where the values of \(\rho_0\) are given. For the lightest nuclei, jumps are observed in the behavior of \(\rho_0\) as a function of \(A\), with the values of \(\rho_0\) for even \(A\) lying somewhat higher, and for odd \(A\) somewhat lower, than would correspond to a smooth behavior. The reason for this is the same as in the case of the jumps of the mean energies mentioned above.

The constancy of the nucleon density at the center of the nucleus for medium and heavy nuclei finds its expression in the well-known empirical

Table III

Nuclear radii and nucleon densities at the center of the nucleus. \(R\) in units of \(r_0\); \(\rho_0\) in units of \(1/r_0^3\)

\(A\)	\(Z\)	\(\rho_0\)	\(R\), calculated from (208)	\(R\), calculated from (210)	Empirical \(R\)
4	2	1.44	1.79	1.34	1.67
5	2	1.22	2.04	1.59	1.80
6	3	1.45	2.04	1.65	1.91
8	4	1.55	2.20	1.87	2.10
10	5	1.66	2.32	2.03	2.26
12	6	1.77	2.41	2.16	2.40
16	8	1.96	2.56	2.39	2.65
20	10	2.13	2.69	2.56	2.85
40	19	2.66	3.14	3.26	3.59
60	29	2.91	3.49	3.74	4.11
80	37	3.04	3.79	4.17	4.52
100	45	3.07	4.07	4.55	4.87
120	54	3.10	4.31	4.88	5.18
140	62	3.07	4.55	5.21	5.45
160	70	3.03	4.78	5.54	5.70
180	77	2.95	5.01	5.85	5.93
200	85	2.90	5.22	6.13	6.14
220	92	2.82	5.44	6.46	6.34
240	99	2.75	5.65	6.76	6.53

...ical law, according to which, for such nuclei, the density of nucleons inside the nucleus does not depend on the mass number.

After establishing the facts indicated, we shall now consider the behavior of the nucleon density in various nuclei and, in connection with this, the calculation of nuclear radii. The density \(\rho\) as a function of \(r\) is presented in Fig. 7 for the cases of a light, a medium, and a heavy nucleus, namely for nuclei with mass numbers 16, 80, and 200.

In addition to the density \(\rho\), the radial density of the distribution is also of importance,

\[ D = 4\pi r^2 \rho, \tag{205} \]

which we have presented in Fig. 8, likewise for nuclei with mass numbers 16, 80, and 200. Whereas \(\rho\) is the number of nucleons per unit volume, \(Ddr\) gives the number of nucleons lying between the surfaces of spheres with radii \(r\) and \(r + dr\).

After having become acquainted with the behavior of the density and the radial distribution of nucleons in light, medium, and heavy nuclei, we shall now speak of nuclear radii, whose determination is by no means as simple as is assumed in some works.

In the case of constant nuclear density, one usually defines as the radius of the nucleus that radius up to which the constant nuclear density extends. According to (133), this radius proves to be proportional to \(A^{1/3}\). Thus, in this case the radius of the nucleus coincides with the radius of a sphere uniformly filled with nucleons. If one starts from this definition, then in the general case one could take as the radius of the nucleus that

Fig. 8. Radial nucleon density \(D\) as a function of \(r\) for nuclei \(A=16,\ Z=8;\ A=80,\ Z=37;\ A=200,\ Z=85\). The numbers standing near the curves indicate the mass numbers of the corresponding nuclei.

value \(r\) for which, on approaching the nucleus, the density begins to become appreciable in magnitude. However, now, when the nucleon density (in contrast to the above-mentioned case of constant density) is not constant, but decreases continuously to zero, there are already several possibilities for defining the nuclear radius; of these we shall discuss the two that seem to us the best justified.

The first possibility consists in the assumption that the nuclear radius \(R\) is proportional to the mean statistical value \(r\). For the mean statistical value of \(r\), taking into account (201), (202), and (140), one obtains

\[ \bar r = \frac{1}{A}\int r\rho\,dv = \frac{4\pi}{A}\int_0^\infty \rho r^3\,dr = \frac{4\pi\rho_0}{A}\int_0^\infty r^3 e^{-a^2 r^2/r_0^2}\,dr = 4\left(\frac{3}{2\pi}\right)^{1/3}\frac{1}{c_0}A^{1/3}r_0 . \tag{206} \]

If one restricts oneself to the medium and heavy, most stable nuclei, then for \(c_0\) one may substitute here the value \(c_{00}=7\), and moreover

the error thereby made does not exceed 2%, and obtain

\[ \bar r = 0.447 A^{1/3} r_0 . \tag{207} \]

For \(A=80\) and \(A=200\), this gives \(\bar r = 1.92 r_0\) and, respectively, \(\bar r = 2.61 r_0\). As is seen from Fig. 7, at these distances the radial nucleon densities by no means vanish yet, but have a considerable value; therefore \(\bar r\) cannot be taken as the radius of the nucleus. However, from Fig. 7 it is easy to see that a length approximately twice as large as \(\bar r\) could, with much greater success, be taken as the radius of the nucleus. Then we would have

\[ R = 2 \bar r = 8 \left( \frac{3}{2\pi} \right)^{1/3} \frac{1}{c_0} A^{1/3} r_0, \tag{208} \]

whence, for medium and heavy nuclei, at \(c_0 = c_{00} = 7.0\), there follows the expression

\[ R = 0.894 A^{1/3} r_0 . \tag{209} \]

Thus, with the exception of light nuclei, the nuclear radii would be proportional to \(A^{1/3}\). For light nuclei, however, this proportionality does not hold, since for them the quantity \(c_{00}\) cannot be regarded as constant; this is indicated, in particular, by the values of \(R\) presented in the second column of Table III, which are calculated from (208) with the aid of the values of \(c_{00}\) indicated in Table I.

Another possibility for defining the radius of a nucleus in the above-mentioned sense consists in taking it to be the radius of such a sphere outside of which there is the same, for all nuclei, number of nucleons \(\eta\) (of the order of one nucleon). Thus \(R\) is determined by means of the equation

\[ \int_R^\infty D\,dr = \eta . \tag{210} \]

The quantity \(R\), in view of the steep fall of \(D\) for large \(r\), is not very sensitive to changes in \(\eta\). For example, for the nucleus indicated in Table I, \(A=120\), \(Z=54\) (\(c_{00}=7.157\)), for \(\eta=1\) one obtains the value \(R=4.63 r_0\), while for \(\eta=\frac{1}{2}\) one obtains the value \(R=4.88 r_0\).

We put \(\eta=\frac{1}{2}\) and then obtained the results given in the third column of Table III, where for \(Z\) and \(c_{00}\) the values indicated in the fourth and, respectively, fifth column of Table I were used.

When comparing the nuclear radii obtained in this way with the empirical ones, we encounter a difficulty which exists

not only here, but also everywhere where, as the radius of the nucleus, one conventionally takes that value \(r\) for which, on approaching the nucleus, the density of nucleons begins to become appreciable. The difficulty is that the empirical nuclear radii were obtained from experiments on the scattering of \(\alpha\)-particles and neutrons, and also from data on \(\alpha\)-decay; moreover, as the radius of the nucleus one considered such a value \(r\) at which the Coulomb forces acting between two nucleons became appreciable. Such a definition agrees, of course, in order of magnitude with the theoretical definition given above, in view of the very small range of action of the nuclear forces; however, exact agreement cannot be expected. Nevertheless, one may hope that the theoretical nuclear radii will show the same dependence on \(A\) as the empirical ones. In order to arrive at nuclear radii defined in the same way as the empirical radii, it would be necessary to calculate the scattering power with respect to \(\alpha\)-particles, protons, or neutrons of a sphere filled with a nucleon gas having a definite density distribution. Such a calculation is proposed for a subsequent paper.

Here we shall content ourselves with comparing the nuclear radii calculated above with empirical ones, which may be represented by the formula

\[ R = 1.42 A^{1/3} 10^{-13}\ \text{cm} = 1.05 A^{1/3} r_0, \tag{211} \]

which indicates the proportionality of the empirical nuclear radii to the quantity \(A^{1/3}\). The nuclear radii obtained from this formula for several values of \(A\) are placed in the sixth column of Table III. As can be seen, the theoretical values obtained above agree well with the experimental ones. In view of what was said above, what appears important is not so much the good numerical agreement as the circumstance that the theoretical nuclear radii show practically the same dependence on \(A\) as the empirical ones. In this respect the nuclear radii calculated by means of formula (210) and placed in the fifth column of Table III behave especially satisfactorily. This is all the more satisfactory since the definition of nuclear radii underlying formula (210), owing to its comparatively slight dependence on \(\eta\), appears physically quite probable.

§ 11. Comparison with the Results of the Wave-Mechanical Calculation of the Deuteron

It is also of interest to compare the parameter \(\lambda\) obtained above with the corresponding parameter obtained on the basis of wave mechanics for some very light nucleus, for example for the deuteron. For this purpose we shall calculate, under the assumption

(47) the Schrödinger energy of the deuteron. We obtain

\[ E=4\pi\,\frac{h^2}{8\pi^2 M_r}\int_0^\infty P'^2\,dr +4\pi\int_0^\infty V(r)P^2\,dr, \tag{212} \]

where \(M_r=\frac{1}{2}M\); \(P\) is the eigenfunction multiplied by \(r\), normalized to unity; \(P'\) is the derivative of \(P\) with respect to \(r\).

We shall determine here also \(E\) and the parameter \(\lambda\) by the Ritz method; for this purpose we make the assumption

\[ P=\left(\frac{a}{2r_0}\right)^{3/2}\frac{1}{\pi^{1/2}}\,r e^{-a\frac{r}{2r_0}}, \tag{213} \]

where \(a\) denotes the variational parameter.

Taking into account (213), as well as (47), (19), (114), and (102), one can represent the energy in the following form:

\[ E=2\chi_J\frac{1}{r_0^2}a^2-\frac{\gamma}{2r_0}\frac{a^3}{(1+a)^2} = \left[ 2\mu_K a^2-\frac{\lambda}{2}\frac{a^3}{(1+a)^2} \right]\varepsilon_0. \tag{214} \]

It now remains to find such a value of \(a\) for which \(E\) has a minimum, and to determine the free parameter so that this minimum coincides with the experimental value of the deuteron binding energy \(E_0=-2.19\) MeV. We obtain

\[ a=1.67 \quad \text{and} \quad \lambda=3.228. \tag{215} \]

This value of \(\lambda\) is only about \(3\%\) smaller than the one we obtained in § 8 for the statistical model of the nucleus with constant density. In comparison with the value \(\lambda\) (200), obtained with the aid of the exact expression for the energy, the difference is larger; however, even in this case it is only about \(22\%\). This difference will decrease in the next approximation of the statistical model, when the neutron and proton densities are varied independently of one another and the correlation of nucleons is taken into account (cf. § 10), since in this case an additional attraction arises, leading to a decrease in the value of \(\lambda\) obtained in § 10.

§ 12. Discussion of the results and further possibilities

The results for the energies and radii of nuclei obtained with the aid of the first approximation of the variational method developed here may be regarded as quite satisfactory, especially if one takes into account that only one empirical parameter was used and that no arbitrary assumptions were made anywhere. Also worthy of attention is the circumstance that the values of the empirical parameter \(\lambda\) in the case of the statistical model and in the case of the wave-mechanical model of the deuteron differ little from one another.

one from another, and this difference becomes still smaller in the following approximations.

The question arises as to what these following approximations consist in and what may be expected from them. The next step would be the independent variation of the neutron and proton densities, as a result of which there arises some decrease in the Coulomb energy of repulsion. This decrease in energy becomes noticeable chiefly for heavy nuclei, for which the Coulomb energy is relatively large. Thus, in this second approximation one may expect a small decrease in the energy, increasing in magnitude with increasing \(A\). At the same time the rise of the energy curve from medium to heavy nuclei will be somewhat reduced. If the energy curve obtained in this way is shifted in the direction of increasing energy, which can be achieved by a certain decrease of the parameter \(\lambda\), then the new energy curve will approximate the empirical values still better.

A further improvement of the behavior of the density in higher approximations of the variational method is achieved by refining the variational assumption (136), and in expression (135) the following terms of the sum must be taken into account; consequently, the density function will have better possibilities for approximating the true function.

A somewhat different refinement of the assumptions made in the present work consists in introducing a correction to the attraction energy by supplementing it with an energy arising from the correlation of nucleons\(^{15}\) (cf. also \(^{13}\)); this leads to the appearance of a certain additional attraction.

It would also be possible to improve the expression for the kinetic energy by subdividing the nucleons according to the secondary quantum number*).

Furthermore, one may try to imagine that some of the nucleons are present in the nucleus in the form of ready-made \(\alpha\)-particles, and determine the percentage of these nucleons again from the principle of minimality.

Finally, it should also be mentioned that, according to Heisenberg\(^{14}\), it is necessary to assume that nuclear forces are very probably caused not by one single kind of particle, but by an entire spectrum of heavy particles (mesons) of different masses and properties. Among these particles the \(\pi\)-mesons are distinguished only by the fact that they are the lightest, and therefore the forces due to them have the greatest range. It follows from this that assumption (47) is valid only for not too small distances between particles. Allowance for these additional forces with very small ranges may be carried out

*) This can be done in the same way as in the case of atoms. See in this connection \(^{12}\) and further \(^{16}\).

only when the properties of these forces have been studied in greater detail.

Those higher approximations or refinements of the theory that can be carried out without further discussion are to be dealt with in subsequent papers.

MATHEMATICAL APPENDIX

Calculation of the integrals \(L(q)\), \(K(q)\), \(H(q)\), and \(M(q)\)

The integrals \(L(q)\), \(K(q)\), \(H(q)\), and \(M(q)\) are defined by (160)—(163). For \(q>1\) they cannot be calculated analytically. We determined the values of the integrals \(L(q)\), \(K(q)\), and \(M(q)\) by a numerical method, with an accuracy of up to four digits, for 18 different values of \(q\) \((q>1)\), and the integral \(H(q)\) for 8 different values of \(q\) \((q>1)\). The values of the integrals, as well as of the function \(F(q)\) computed with their aid, are presented in Table IV. The results given in § 10 were obtained by means of these values of the integrals.

Table IV

Values of the integrals \(L(q)\), \(K(q)\), \(H(q)\), and \(M(q)\), as well as the values of the function \(F(q)\), found numerically with an accuracy of up to 4 digits.

\(q\)	\(L(q)\)	\(K(q)\)	\(H(q)\)	\(M(q)\)	\(F(q)\)
1.3355	0.2228	0.07044	0.03619	0.05908	0.02530
1.8951	0.3889	0.1149	0.05735	0.07302	0.06126
2.8427	0.6930	0.1857	0.08901	0.08818	0.1526
3.6000	0.9392	0.2357	0.1111*	0.09602	0.2464
4.0998	1.099	0.2656	0.1227	0.09991	0.3164
4.7378	1.300	0.3006	0.1369	0.1039	0.4131
5.4000	1.503	0.3336	0.1510*	0.1071	0.5205
6.0000	1.682	0.3611	0.1619*	0.1096	0.6229
6.6329	1.865	0.3879	0.1714	0.1117	0.7349
7.1000	1.997	0.4064	0.1794*	0.1130	0.8203
7.5000	2.107	0.4215	0.1852*	0.1141	0.8945
7.7500	2.176	0.4306	0.1888*	0.1147	0.9416
8.0000	2.242	0.4394	0.1921*	0.1152	0.9891
8.2500	2.310	0.4481	0.1954*	0.1158	1.037
8.5280	2.384	0.4574	0.1990*	0.1163	1.091
9.0000	2.505	0.4727	0.2048*	0.1172	1.183
9.3636	2.601	0.4840	0.2079	0.1178	1.255
11.9173	3.216	0.5540	0.2340	0.1211	1.774

The asterisk marks those values of \(H(q)\) that were calculated with the aid of the approximate formula (222); these approximate values are approximately \(0.7\%\) greater than the exact ones.

For the integrals \(K(q)\), \(H(q)\), and \(M(q)\) one can derive approximate formulas which closely approximate the exact integrals in the range most important for computations, \(q(5<q<8)\). If we confine ourselves to this interval, then in all three integrals the integrands will be the product of a certain function having a narrow and steep maximum of the type \(f(x)=x^{2}e^{-n x^{4}}\) \((n=2,3,4)\), and vanishing practically everywhere except in the neighborhood of this maximum, and also a function \(\ln(1+q^{2}e^{-2x^{2}})\) or \(\operatorname{arctg}(1+qe^{-x^{2}})\), which varies comparatively slowly in the neighborhood of the maximum of \(f(x)\) and for which we shall introduce the abbreviated notation \(g(x)\). By the mean-value theorem of integral calculus, in the integral

\[ \int_{0}^{\infty} f(x)g(x)\,dx \]

one may replace one of the functions, for example the slowly varying function \(g(x)\), by a suitably chosen constant value. In view of the behavior of the function \(f(x)\) described above, a good approximate value for the integral will be obtained if, as the constant value of the function \(g(x)\), one chooses the value assumed by \(g(x)\) at the point \(x_m\) of the sharp maximum of \(f(x)\). Thus one obtains

\[ \int_{0}^{\infty} f(x)g(x)\,dx = g(x_m)\int_{0}^{\infty} f(x)\,dx = g(x_m)\frac{\pi^{1/2}}{4n^{3/2}}, \tag{216} \]

\[ x_m=\frac{1}{n^{1/2}}. \tag{217} \]

The accuracy of this approximation is still insufficient for our purposes. However, we immediately obtain an approximate formula possessing sufficient accuracy if in formula (216) we replace \(g(x_m)\) by \(\frac{1}{2}[g(x_1)+g(x_2)]\), where \(x_1\) and \(x_2\) are those values of \(x\) at which \(f(x)\) is equal to half its value at the maximum, i.e., is equal to \(\frac{1}{2}f(x_m)\). Then we have:

\[ \int_{0}^{\infty} f(x)g(x)\,dx = \frac{1}{2}\{g(x_1)+g(x_2)\}\frac{\pi^{1/2}}{4n^{3/2}}, \tag{218} \]

where \(x_1\) and \(x_2\) denote the two roots of the transcendental equation

\[ f(x)=x^{3}e^{-n x^{2}}=\frac{1}{2}f(x_m)=\frac{1}{2n}e^{-1}, \tag{219} \]

whence it follows that

\[ x_1=\left(\frac{0.2320}{n}\right)^{1/2},\quad x_2=\left(\frac{2.678}{n}\right)^{1/2}. \tag{220} \]

STATISTICAL THEORY OF ATOMIC NUCLEI

If in (216) the original functions are substituted for \(f(x)\) and \(g(x)\), then for \(K(q)\), \(H(q)\), and \(M(q)\) the following very good approximate formulas are obtained:

\[ \begin{aligned} K(q) &= \int_{0}^{\infty} x^3 e^{-2x^2}\ln\left(1+q^2 e^{-2x^2}\right)\,dx = \\ &= \frac{1}{16}\left(\frac{\pi}{2}\right)^{1/2} \left[ \ln\left(1+q^2 e^{-2x_1^2}\right) + \ln\left(1+q^2 e^{-2x_2^2}\right) \right], \end{aligned} \tag{221} \]

where

\[ x_1=\left(\frac{0,2320}{2}\right)^{1/2},\qquad x_2=\left(\frac{2,678}{2}\right)^{1/2}; \]

\[ \begin{aligned} H(q) &= \int_{0}^{\infty} x^3 e^{-4x^2}\ln\left(1+q^2 e^{-2x^2}\right)\,dx = \\ &= \frac{\pi^{1/2}}{64} \left[ \ln\left(1+q^2 e^{-2x_1^2}\right) + \ln\left(1+q^2 e^{-2x_2^2}\right) \right], \end{aligned} \tag{222} \]

where

\[ x_1=\left(\frac{0,2320}{4}\right)^{1/2},\qquad x_2=\left(\frac{2,678}{4}\right)^{1/2}; \]

\[ \begin{aligned} M(q) &= \int_{0}^{\infty} x^2 e^{-3x^2}\operatorname{arctg}\left(1+q e^{-x^2}\right)\,dx = \\ &= \frac{1}{24}\left(\frac{\pi}{3}\right)^{1/2} \left[ \operatorname{arctg}\left(1+q e^{-x_1^2}\right) + \operatorname{arctg}\left(1+q e^{-x_2^2}\right) \right], \end{aligned} \tag{223} \]

where

\[ x_1=\left(\frac{0,2320}{3}\right)^{1/2},\qquad x_2=\left(\frac{2,678}{3}\right)^{1/2}. \]

The approximate formula for \(K(q)\) gives a maximum deviation from the exact value in the interval \(5.4<q<10\) of less than \(0.6\%\), and in the interval \(4<q<5.4\) of less than \(1.2\%\). The approximate formula for \(H(q)\) gives a maximum error in the interval \(4<q<10\) of less than \(0.7\%\). The best approximation is achieved by the approximate formula for \(M(q)\), which gives a maximum error in the interval \(1.5<q<10\) of less than \(0.4\%\).

Since the integral \(H(q)\) enters the expression for the energy of the nucleus only in a term of subordinate importance, for \(H(q)\) one may use the approximate formula (222) instead of the exact value. Some values calculated by means of this approximate formula are given in Table IV.

For the integral \(L(q)\), of course, the approximation derived here cannot be applied, since in the integrand in this

in this case there is no factor \(e^{-nx^2}\), which ensures the rapid disappearance of \(f(x)\), essential for our approximation, as \(x\) increases. Since \(L(q)\) in the indicated interval of \(q\) varies approximately linearly, one can approximate \(L(q)\) well in this interval by a simple formula whose coefficients are determined by the method of least squares. We find that in the interval \(1.5<q<10\) the approximate formula

\[ L(q)=-0.0051755q^2+0.35220q-0.24814 \tag{224} \]

represents the exact integral with an error of less than \(1\%\). Since, in calculating \(F(q)\), this integral has a subordinate significance in comparison with the other terms in \(F(q)\), this accuracy proves to be quite sufficient.

We have given here approximate formulas for \(K(q)\), \(H(q)\), \(M(q)\), and \(L(q)\), on the one hand because of their simplicity, and on the other in order to show that all calculations by means of the variational method can be carried out with an insignificant expenditure of labor, which is essential for the further development of this method.

The numerical calculations were carried out by my assistants O. Kunvari, E. Magori, V. Molnar, and E. Szabo; the drawings were made by assistants G. Knapek and L. Zelenka. I would also like here to express to all of them my gratitude for their work.

References

1. E. Majorana, Zeits. f. Phys., 82, 137 (1933).
2. W. Heisenberg, Rapport du VII-ième Congrès Solvay, Paris, 1934.
3. C. F. v. Weizsäcker, Zeits. f. Phys., 96, 431 (1935).
4. C. F. v. Weizsäcker, Zeits. f. Phys., 95, 431 (1935).
5. G. C. Wick, Nuovo Cimento Nr. 4 (1934); Rend. Accad. Lincei, 19, 319 (1934); 21, 170 (1935).
6. K. Nakabayasi, Zeits. f. Phys., 97, 211 (1935).
7. S. Flügge, Zeits. f. Phys., 96, 459 (1935).
8. F. S. Wang, Zeits. f. Phys., 100, 736 (1936).
9. H. A. Bethe und R. F. Bacher, Rev. Mod. Phys., 8, 82 (1936).
10. C. F. v. Weizsäcker, Die Atomkerns, Physik und Chemie und ihre Anwendungen in Einzeldarstellungen, T. II, Akad. Verlagsges. Leipzig, 1937.
11. L. Rosenfeld, Nuclear forces, Monographs on Theoretic and Applied Physics, T. I. North-Holland Publishing Co., Amsterdam, 1948.
12. P. Gombás, Statistical Theory of the Atom and Its Applications. IL, 1951.
13. P. Gombás, The Many-Particle Problem in Quantum Mechanics. IL, 1952.
14. W. Heisenberg, Theorie des Atomkerns, Max-Planck-Institut für Physik, Göttingen, 1951.
15. E. Wigner, Phys. Rev. (2) 46, 1002 (1934), Trans. Faraday Soc. 34, 678 (1938).
16. P. Gombás, Acta Phys. Hung. 1, No. 3 (1952).