text:
# The structural phase transition and loss of magnetic moments in NpO2: ab initio approach to the crystal and mean field

## I Introduction

Since the discovery of the low temperature phase transition in NpO2 at $`T_c=25.5`$ K almost fifty years ago, Wes  the magnetic properties of this compound represent a challenge for theoretical interpretation. Erd ; San1  NpO2 and its neighbor in the actinide series, UO2, Fra ; Cac  are both crystallized in the cubic fluorite (CaF2) structure in the high temperature phase. The temperature dependence of the magnetic susceptibilities of these compounds is suggestive for a Néel transition to an antiferromagnetic state. Erd2 ; Fra  However, while it was found that UO2 orders antiferromagnetically ($`T_N=30.8`$ K) with 1.74$`\mu _B`$/U atom, Fra  no magnetic ordering Cac2  was confirmed for NpO2 (see for a review Ref. San1, ). As follows from Mössbauer experiments Fri  the upper limit on the magnetic moments of Np in the ordered phase is only 0.01$`\mu _B`$/Np. This implies almost a complete quenching of the 3.0$`\mu _B`$ magnetic moment of Np, which it exhibits in the paramagnetic high temperature phase. The demagnetization of Np magnetic moments seems to be incompatible with the widely accepted viewpoint that the Np cations are in the tetravalent valence state, Np4+, with three localized 5f electrons left at each neptunium site. With the $`5f^3`$ configuration, Np becomes a Kramers ion having a magnetic doublet ($`\mathrm{\Gamma }_6`$, $`\mathrm{\Gamma }_7`$) or quadruplet ($`\mathrm{\Gamma }_8`$) ground state. Because of this, a divergent susceptibility at $`T0`$ will be present in any model of the phase transition. San1 ; San2

To describe the disappearance of magnetic moments of Np at $`T<T_c`$ Santini and Amoretti put forward an idea of a magnetic octupole order parameter, San2  which on one hand, is not invariant under time reversal symmetry and, on the other hand, is different from the magnetic dipolar order parameter which brings about the ordinary magnetic ordering. The idea of breaking down of Kramers’ degeneracy by the magnetic octupole seems to be the only solution in the framework of the $`5f^3`$ model. San2 ; San1

An alternative explanation was offered by Friedt et al. Fri . They suggested that the magnetic anomaly in NpO2 could be explained by a splitting of the ground state $`\mathrm{\Gamma }_8`$ quadruplet of cubic symmetry into two doublets in the crystal field of lower symmetry of the ordered phase. The splitting was ascribed to an internal distortion of the oxygen sublattice at $`T_c`$. This viewpoint was inspired by the experimental observation that the phase transition in UO2 is accompanied by an internal distortion of the oxygen cube that surrounds the U cation, Fab  while the external cubic structure of uranium dioxides survives the transition. However, the mechanism mainly weakens the magnetic response from the Np sites at 20 K, while the $`T0`$ divergence of magnetic susceptibility still persists in the model in contradistinction with experiment. Erd2 ; Cac2  On the other hand, no evidence for an internal or external crystallographic distortion in NpO2 has been found by synchrotron experiments. Man  Thus, the phase transition in NpO2 appears to be isostructural like the $`\gamma \alpha `$ phase transition in elemental cerium Kos ; Mal  ($`T_c100`$ K) or the isostructural expansion in YbInCu4 at $`T=42`$ K. Sar  It is worth mentioning that in all these compounds the phase transition is accompanied by a loss of magnetic moments in the ordered phase.

However, very recent resonant X-ray scattering (RXS) experiments at the Np $`M_{IV}`$ and $`M_V`$ edges in NpO2 indicated an unexpected result: the phase transition is not isostructural. Pai  In the low temperature phase a long range order of Np electric quadrupoles was revealed by the growth of superlattice Bragg peaks. Man ; Pai  The space symmetry of the ordered phase was identified as $`Pn\overline{3}m`$. Pai  The symmetry lowering is a special one. The centers of mass positions of neptunium and oxygen remain in the cubic CaF2 structure as in the high temperature phase. However, electronic quadrupoles of the Np sites have four different orientations, which allows us to distinguish four different sublattices of Np cations. Pai

This experimental finding gives rise to a question of correlations between structural and magnetic properties in solids, which we have investigated theoretically in our model of the $`\gamma \alpha `$ phase transition in Ce. NM1 ; NM2 ; NM3  Unlike NpO2, pristine cerium is a metal and the phase transition there at normal pressure is accompanied by a huge volume change. Kos ; Mal  The volume anomaly in NpO2 (contraction) is only 0.018%. Man  However, the disappearance of the magnetic moments and the “isostructural” character of the phase transition makes them similar. In this respect it is interesting to notice that in our model for Ce we have predicted the $`Pa\overline{3}`$ space symmetry for the ordered $`\alpha `$ phase NM1  which is very close to the $`Pn\overline{3}m`$ structure reported for NpO2 in Ref. Pai, . The active irreducible representation also belongs to the $`X`$ point of the Brillouin zone (BZ). While the existence of superstructure reflections for $`\alpha `$-Ce remains an open question which has to be investigated experimentally, NM1 ; NM2 ; NM3  here we want to apply our theoretical concepts for the study of the crystal- and mean-field in neptunium dioxide.

We will use the technique of multipole expansion of the Coulomb interaction. NM1 ; NM2 ; NM3 ; NM4  The multipole expansion represents a unified description of the crystal field effects and atomic term splitting. The intrasite Coulomb repulsion which is responsible for Hund’s rules, the spin-orbit coupling and the crystal field effects are included on equal footing. (For a single site the technique is equivalent to the classical description of atomic terms. CS ) In order to calculate the corresponding electronic spectra we use many Slater determinants, which indicates that our scheme is a genuine many electron approach NM3 ; NM4  corresponding to the configuration interaction (CI). (Notice that ordinarily used single-determinant Hartree-Fock approach is not sufficient for these purposes).

The paper is organized as follows. In Sec. II we discuss the triple-$`\stackrel{}{q}^X`$ structures ($`Pn\overline{3}m`$ or $`Pa\overline{3}`$) of the ordered phases and derive the selection rules for resonant X-ray scattering (RXS) experiments. In Sec. III we describe our method for treating local electron configurations at neptunium sites. Next (Sec. IV) we examine thoroughly the case of three localized $`5f`$ electrons. In Sec. V we generalize the $`5f^3`$ model by including the multipole interactions with a valence electron which is present instantaneously on a Np site. We show that in the framework of this extended model the disappearance of the magnetic moments can be understood and ascribed to the trigonal symmetry lowering at 25 K. Our conclusions are summarized in Sec. VI.

## II Structural phase transitions and RXS experiments

In this section we discuss the allowed space groups of the low temperature phase of NpO2. Our considerations here are based on a general group theoretical approach which is independent of the model assumptions.

In the disordered phase ($`T>25`$ K) the crystal structure is cubic ($`Fm\overline{3}m`$) and the site group is $`O_h`$. The electron density of neptunium is given by the spherically symmetric component ($`Y_0^0`$) and the cubic harmonics $`K_l(\mathrm{\Omega })`$ with $`l=4,6`$. There is no contribution to the density from the quadrupole functions $`Y_{l=2}^\tau `$, where $`\tau =0`$, $`(m,c)`$ or $`(m,s)`$, $`m=1,2`$. (Here and below we work with the real spherical harmonics with the phase definition of Ref. Bra, .) At $`T_c=25`$ K the transition to a new phase sets in. The new phase was characterized as a triple-$`\stackrel{}{q}`$ antiferro- quadrupolar ordering of $`T_{2g}`$ ($`\mathrm{\Gamma }_{5g}`$) electric quadrupoles at the Np sites. Pai

In real space the ordering is characterized by four different sublattices of the simple cubic structure. We label these sublattices which contain the sites (0,0,0) $`(a/2)`$(0,1,1), $`(a/2)`$(1,0,1) and $`(a/2)`$(1,1,0) by $`\{\stackrel{}{n}_p\}`$, $`p=14`$, respectively. The most significant feature of the ordered phase is the existence of only one three-fold axis of symmetry $`C_3`$ at each Np site which is also a cube diagonal. The only quadrupole function compatible with the symmetry lowering is $`Y_2^0(\mathrm{\Omega }^{})`$ in the coordinate system where the $`z^{}`$-axis coincides with one of the three-fold axes (cube diagonals): $`[111]`$, $`[1,1,1]`$, $`[1,1,1]`$, and $`[1,1,1]`$. Consequently, there are four such functions which are given by

$`𝒮_a(\mathrm{\Omega })={\displaystyle \frac{1}{\sqrt{3}}}(Y_2^{1s}(\mathrm{\Omega })+Y_2^{1c}(\mathrm{\Omega })+Y_2^{2s}(\mathrm{\Omega })),`$ (1a)
$`𝒮_b(\mathrm{\Omega })={\displaystyle \frac{1}{\sqrt{3}}}(Y_2^{1s}(\mathrm{\Omega })Y_2^{1c}(\mathrm{\Omega })+Y_2^{2s}(\mathrm{\Omega })),`$ (1b)
$`𝒮_c(\mathrm{\Omega })={\displaystyle \frac{1}{\sqrt{3}}}(Y_2^{1s}(\mathrm{\Omega })Y_2^{1c}(\mathrm{\Omega })Y_2^{2s}(\mathrm{\Omega })),`$ (1c)
$`𝒮_d(\mathrm{\Omega })={\displaystyle \frac{1}{\sqrt{3}}}(Y_2^{1s}(\mathrm{\Omega })+Y_2^{1c}(\mathrm{\Omega })Y_2^{2s}(\mathrm{\Omega })).`$ (1d)
The spherical harmonics $`Y_2^{1s}`$, $`Y_2^{1c}`$, $`Y_2^{2s}`$ belong to a three-dimensional irreducible representation $`T_{2g}`$ of $`O_h`$. They are proportional to the Cartesian components $`yz`$, $`zx`$, and $`xy`$.

In Ref. Pai,  the space group of the ordered phase was identified as $`Pn\overline{3}m`$ (No. 224, Ref. Tables, ). The corresponding ordering of $`𝒮_a𝒮_d`$ functions is shown in Fig. 1, left panel.

However, a close examination reveals that there is another possibility which was not considered by the authors of Ref. Pai, . This is the $`Pa\overline{3}`$ structure (No. 205 of Ref. Tables, ) depicted in the right panel of Fig. 1. The $`Pa\overline{3}`$ structure also realizes a triple-$`\stackrel{}{q}`$ quadrupolar order, the site symmetry of the Np sites is $`S_6=C_3\times i`$.

The problem of finding the appropriate space group can be simplified to the task of distributing four orientational functions $`𝒮_a𝒮_d`$ among four distinct sublattices $`\{\stackrel{}{n}_p\}`$, $`p=14`$. Keeping the three-fold rotation symmetry about the $`[111]`$ axis and $`𝒮_a`$ for the first sublattice, we distinguish only three different choices here. The first is when $`(i)`$ $`𝒮_b`$ corresponds to $`\{\stackrel{}{n}_4\}`$, $`𝒮_c\{\stackrel{}{n}_2\}`$, and $`𝒮_d\{\stackrel{}{n}_3\}`$. This gives the $`Pn\overline{3}m`$ structure, Fig. 1, left panel. Another assignment is $`(ii)`$ $`𝒮_b\{\stackrel{}{n}_2\}`$, $`𝒮_c\{\stackrel{}{n}_3\}`$, and $`𝒮_d\{\stackrel{}{n}_4\}`$. This scheme corresponds to the $`Pa\overline{3}`$ structure as shown in Fig. 1, right panel. The third choice is $`(iii)`$ $`𝒮_b\{\stackrel{}{n}_3\}`$, $`𝒮_c\{\stackrel{}{n}_4\}`$, and $`𝒮_d\{\stackrel{}{n}_2\}`$. This is simply another variant (or domain) of the $`Pa\overline{3}`$ structure.

From the mathematical point of view we are analyzing the symmetry lowering and the condensation schemes at the $`X`$ point of the BZ of the fcc lattice, which involves the density components of the $`T_{2g}`$ symmetry. In both cases ($`Pn\overline{3}m`$ and $`Pa\overline{3}`$) the triple-$`\stackrel{}{q}`$ mode which drives the structural phase transition belongs to the $`X`$ point of BZ $`{}_{}{}^{}\stackrel{}{q}_{}^{X}`$, and involves the three arms of the latter, $`\stackrel{}{q}_x^X=(2\pi /a)(1,0,0)`$, $`\stackrel{}{q}_y^X=(2\pi /a)(0,1,0)`$, and $`\stackrel{}{q}_z^X=(2\pi /a)(0,0,1)`$. To describe the three quadrupolar components $`S_1=Y_2^{1s}`$, $`S_2=Y_2^{1c}`$, $`S_3=Y_2^{2s}`$ at site $`\stackrel{}{n}`$ we introduce the functions $`S_i(\stackrel{}{n})`$ ($`i=13`$) and consider their Fourier transforms,

$`S_i(\stackrel{}{q})={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{\stackrel{}{n}}{}}e^{i\stackrel{}{q}\stackrel{}{X}(\stackrel{}{n})}S_i(\stackrel{}{n}),`$ (2)
where the position vector $`\stackrel{}{X}`$ runs over $`N`$ sites of the face centered cubic (fcc) neptunium lattice. The little group of $`{}_{}{}^{}\stackrel{}{q}_{}^{X}`$ is $`D_{4h}`$ ($`4/mmm`$). At $`\stackrel{}{q}_x^X`$ the functions $`S_2`$ and $`S_3`$ belong to the $`E_g`$ small representation, while the function $`S_1`$ belong to the $`B_{2g}`$ small representation. Consequently, we distinguish two different irreducible representations of the space group $`Fm\overline{3}m`$. The first one, $`X_5^+`$ in notations of Stokes and Hatch, Sto  comprises six orientational functions, with two functions from every arm of $`{}_{}{}^{}\stackrel{}{q}_{}^{X}`$. The functions are $`S_2(\stackrel{}{q}_x^X)`$, $`S_3(\stackrel{}{q}_x^X)`$; $`S_1(\stackrel{}{q}_y^X)`$, $`S_3(\stackrel{}{q}_y^X)`$; $`S_1(\stackrel{}{q}_z^X)`$, and $`S_2(\stackrel{}{q}_z^X)`$. The second irreducible representation, $`X_4^+`$, has only three components, with one function from each arm of $`{}_{}{}^{}\stackrel{}{q}_{}^{X}`$, namely $`S_1(\stackrel{}{q}_x^X)`$, $`S_2(\stackrel{}{q}_y^X)`$ and $`S_3(\stackrel{}{q}_z^X)`$. For the symmetry lowering to the $`Pn\overline{3}m`$ structure we obtain the following condensation scheme:

$`Fm\overline{3}m:X_4^+`$
$`(S_1(\stackrel{}{q}_x^X)=S_2(\stackrel{}{q}_y^X)=S_3(\stackrel{}{q}_z^X)=\rho _1\sqrt{N}0)`$
$`Pn\overline{3}m(Z=4).`$ (3)
Here $`\mathrm{}`$ stand for a quantum and thermal average and $`\rho _1`$ is the order parameter amplitude. The $`X_4^+`$ active representation corresponds to the irreducible star $`\{\stackrel{}{k}_{10}\}`$ and the loaded one-dimensional representation $`\widehat{\tau }^7`$ in Kovalev’s notation. Kov  There are four domains for the $`Pn\overline{3}m`$ structure.

The condensation scheme for the $`Pa\overline{3}`$ structure is given by

$`Fm\overline{3}m:X_5^+`$
$`(S_3(\stackrel{}{q}_x^X)=S_1(\stackrel{}{q}_y^X)=S_2(\stackrel{}{q}_z^X)=\rho _2\sqrt{N}0,`$
$`S_2(\stackrel{}{q}_x^X)=S_3(\stackrel{}{q}_y^X)=S_1(\stackrel{}{q}_z^X)=0)`$
$`Pa\overline{3}(Z=4)\text{domain I}.`$ (4)
Here again $`\rho _2`$ is the order parameter amplitude. This condensation scheme corresponds to a domain shown in Fig. 1, and there are eight possible domains of $`Pa\overline{3}`$. In particular, there is a domain $`(iii)`$ which we have already considered before ($`𝒮_a\{\stackrel{}{n}_1\}`$, $`𝒮_b\{\stackrel{}{n}_3\}`$, $`𝒮_c\{\stackrel{}{n}_4\}`$, and $`𝒮_d\{\stackrel{}{n}_2\}`$) when we analyzed the variants of distribution of $`𝒮_a𝒮_d`$ over four sublattices. This domain is obtained as a result of the following symmetry breaking

$`Fm\overline{3}m:X_5^+`$
$`(S_3(\stackrel{}{q}_x^X)=S_1(\stackrel{}{q}_y^X)=S_2(\stackrel{}{q}_z^X)=0,`$
$`S_2(\stackrel{}{q}_x^X)=S_3(\stackrel{}{q}_y^X)=S_1(\stackrel{}{q}_z^X)=\rho _2\sqrt{N}0)`$
$`Pa\overline{3}(Z=4)\text{domain II}.`$ (5)
Previously the symmetry change $`Fm\overline{3}mPa\overline{3}`$ has been considered in a number of publications, Refs. Zie, ; Mic1, ; NM1, ; NM2, . The $`X_5^+`$ active representation is defined by the star $`\{\stackrel{}{k}_{10}\}`$ and the loaded two-dimensional representation $`\widehat{\tau }^9`$ in Kovalev’s notation. Kov

Although the two structures look similar there is one very important difference between them. The $`Pa\overline{3}`$ ordering corresponds to an effective attraction of the quadrupoles, while the $`Pn\overline{3}m`$ ordering results in a repulsive interaction between them, see Appendix A for details. This is the main reason why the $`Pa\overline{3}`$ structure (but not $`Pn\overline{3}m`$) is found in many molecular solids such as H2, N2,Scott  NaO2,Zie  and C60.Sach1 ; Dav ; Sach2 ; Hei2 ; Mic1  Here we speak about the direct bilinear electronic quadrupole-quadrupole interaction. However, in case of NpO2 there are indirect (superexchange super ) interactions via oxygen which may lead to an effective attraction and therefore, the $`Pn\overline{3}m`$ structure can not be ruled out. Nevertheless, the $`Pa\overline{3}`$ space group is a good candidate for the symmetry of the ordered phase of NpO2. Provided that the superexchange interaction remains the same for two structures, the $`Pa\overline{3}`$ symmetry becomes preferable. It is interesting to notice that at first the symmetry of the ordered phase of the pristine C60 molecular crystal was identified as $`Pn\overline{3}`$.Hei2  Only later it was corrected and proved to be $`Pa\overline{3}`$. Sach1 ; Dav ; Sach2

Now we will study how the $`Pn\overline{3}m`$ and $`Pa\overline{3}`$ structures manifest themselves in resonant X-ray scattering experiments. Tem1 ; Pai  The tensor of scattering on the quadrupolar density of a neptunium atom in a cubic lattice is given by Dmi1 ; Dmi2 ; Tem2

$`\widehat{f}_n=\stackrel{~}{f}\left(\begin{array}{ccc}0& \rho _1& \rho _2\\ \rho _1& 0& \rho _3\\ \rho _2& \rho _3& 0\end{array}\right),`$ (6)
where $`\rho _i=\pm 1`$. We have the following correspondence Dmi2 ; Tem2  between the functions $`𝒮_n`$, Eq. (1a-d), and the scattering tensors: the function $`𝒮_a`$ corresponds to $`\widehat{f}_a`$, where $`\rho _1=\rho _2=\rho _3=1`$, the function $`𝒮_b`$ to $`\widehat{f}_b`$ with $`\rho _1=1`$, $`\rho _2=\rho _3=1`$; the function $`𝒮_c`$ to $`\widehat{f}_c`$ with $`\rho _1=\rho _2=1`$, $`\rho _3=1`$; and $`𝒮_d`$ to $`\widehat{f}_d`$ with $`\rho _1=\rho _3=1`$, $`\rho _2=1`$. The scattering tensor structure amplitude is found as

$`(h,k,l)={\displaystyle \underset{n=1}{\overset{N}{}}}\widehat{f}_ne^{i\stackrel{}{K}\stackrel{}{X}(n)},`$ (7)
where $`\stackrel{}{K}=(2\pi /a)(h,k,l)`$ stands for the vectors of the reciprocal lattice and $`a`$ is the cubic lattice constant. In performing the summation in (7) we distinguish four contributions from four sublattices $`\{n_p\}`$, $`p=14`$. We are mainly interested in superstructure Bragg reflections. We obtain then that in general

$`(h,k,l)=\stackrel{~}{F},`$ (8)
where $`\stackrel{~}{F}=\stackrel{~}{f}N`$ and the matrix $``$ is either $`A`$, $`B`$, or $`C`$:

$`A=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right),B=\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right),C=\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right).`$
(9)
In Table 1 we quote which of the matrices ($`A,B,C`$) occurs for every particular case of $`(h,k,l)`$ for the $`Pn\overline{3}m`$ and $`Pa\overline{3}`$ structures.

We recall that the conditions Tables  for the isotropic scattering on the fcc lattice of Np (from the spherically symmetric densities) are $`h+k,k+l=2n`$, that corresponds to two first lines of Table 1.

Finally, we consider the polarization dependencies for resonant X-ray scattering experiments. We assume that a crystal has a flat (001) surface and the azimuthal angle $`\psi `$ is counted from the $`x`$axis defined by $`\stackrel{}{e}_x=(1,0,0)`$. Pai  The incident beam direction is given by ($`\mathrm{cos}\mathrm{\Theta }\mathrm{cos}\psi `$, $`\mathrm{cos}\mathrm{\Theta }\mathrm{sin}\psi `$, $`\mathrm{sin}\mathrm{\Theta }`$), the scattered beam direction by ($`\mathrm{cos}\mathrm{\Theta }\mathrm{cos}\psi `$, $`\mathrm{cos}\mathrm{\Theta }\mathrm{sin}\psi `$, $`\mathrm{sin}\mathrm{\Theta }`$). We introduce the standard polarization vectors $`\stackrel{}{e}_i`$ ($`i=13`$) parallel ($`\stackrel{}{\pi }_i`$, $`\stackrel{}{\pi }_f`$) or perpendicular ($`\stackrel{}{\sigma }`$) to the scattered plane,

$`\stackrel{}{\sigma }=(\mathrm{sin}\psi ,\mathrm{cos}\psi ,\mathrm{\hspace{0.33em}0}),`$ (10a)
$`\stackrel{}{\pi }_i=(\mathrm{sin}\mathrm{\Theta }\mathrm{cos}\psi ,\mathrm{sin}\mathrm{\Theta }\mathrm{sin}\psi ,\mathrm{cos}\mathrm{\Theta }),`$ (10b)
$`\stackrel{}{\pi }_f=(\mathrm{sin}\mathrm{\Theta }\mathrm{cos}\psi ,\mathrm{sin}\mathrm{\Theta }\mathrm{sin}\psi ,\mathrm{cos}\mathrm{\Theta }).`$ (10c)
A structure factor amplitude $`F(h,k,l)`$ is found as Dmi1 ; Dmi2 ; Tem2

$`F_{pp^{}}(h,k,l)=\stackrel{~}{F}\stackrel{}{e}_p^{}^T\stackrel{}{e}_p,`$ (11)
where $`T`$ stands for transpose and $`pp^{}`$ denotes one of the four polarization channels, $`\sigma \sigma `$, $`\sigma \pi _f`$, $`\pi _i\sigma `$, and $`\pi _i\pi _f`$. The structure factor amplitudes for the three cases ($`=A,B,C`$) are quoted in Table 2.

The corresponding intensities of the superstructure Bragg reflections are found as $`I_{pp^{}}(h,k,l)=|F_{pp^{}}(h,k,l)|^2`$.

From Tables 1 and 2 one can easily obtain all necessary dependencies for intensities of different polarizations. For example, the intensity of the $`(003)`$ reflection of $`Pn\overline{3}m`$ is exactly the same as the intensity of the same polarization of the $`(300)`$ reflection of domain I $`(ii)`$ of $`Pa\overline{3}`$, and the intensity of the $`(030)`$ reflection of domain II $`(iii)`$ of $`Pa\overline{3}`$. Furthermore, the $`(013)`$ reflection of domain I $`(ii)`$, and the $`(103)`$ reflection of domain II $`(iii)`$ of $`Pa\overline{3}`$ also have the same intensity, and so on. Since both symmetries are very similar, we believe that special care should be taken to discriminate between the $`Pn\overline{3}m`$ and $`Pa\overline{3}`$ structures.

Tables 1 and 2 are also useful in considering the contributions from the domains of the same group. We start with the $`Pn\overline{3}m`$ structure. The four domains differ by the orientation of the $`Y_2^0(\mathrm{\Omega }^{})`$ quadrupolar function at $`\stackrel{}{n}=0`$. This function may be chosen to align along four main cube diagonals, which are connected with each other through rotations by $`\pi /2`$ about the $`z`$axis. Hence, all four domains are obtained from the first, $`(i)`$, by applying three consecutive rotations by $`\pi /2`$ about the $`z`$axis. The effect of the domains in RXS experiments is equivalent to a superposition of four structure amplitudes, Table 2, with angles $`\psi `$, $`\psi +\pi /2`$, $`\psi +\pi `$, and $`\psi +3\pi /2`$. (Here we suggest that the domains produce the coherent scattering.) We label these domains by indices $`d1d4`$, and introduce their populations, $`P_{d1}P_{d4}`$. The population of a domain $`d`$ is defined as $`P_d=N_d/N`$, where $`N_d`$ is the total number of neptunium atoms in the domain. Then we obtain from Eq. (11) that the structure amplitude for the $`\sigma \sigma `$ and $`\sigma \pi _f`$ channels of the $`(003)`$ reflection is given by

$`F_{\sigma \sigma }(003)`$ $`=`$ $`\stackrel{~}{F}\stackrel{~}{P}\mathrm{sin}(2\psi ),`$ (12a)
$`F_{\sigma \pi _f}(003)`$ $`=`$ $`\stackrel{~}{F}\stackrel{~}{P}\mathrm{cos}(2\psi )\mathrm{sin}\mathrm{\Theta },`$ (12b)
where
$`\stackrel{~}{P}=\left(P_{d1}P_{d2}+P_{d3}P_{d4}\right).`$ (12c)
If $`P_{d1}+P_{d3}=P_{d2}+P_{d4}`$, $`\stackrel{~}{P}=0`$ and the $`(003)`$ Bragg reflection is suppressed. Otherwise, the $`\psi `$ and $`\mathrm{\Theta }`$dependencies in (12a,b) are exactly the same as for a single domain, Tables 1, 2. For the $`(300)`$ reflection of $`Pn\overline{3}m`$ one obtains that

$`F_{\sigma \pi _f}(300)=\stackrel{~}{F}P^{}\mathrm{cos}\mathrm{\Theta }\mathrm{cos}(\psi +\psi _0),`$ (13a)
where
$`\psi _0=\mathrm{arccos}([P_{d1}P_{d3}]/P^{}),`$ (13b)
$`P^{}=\sqrt{[P_{d1}P_{d3}]^2+[P_{d2}P_{d4}]^2}.`$ (13c)
Comparing this result with the expression for a single domain, we observe that the main effect is the phase shift $`\psi _0`$ given by Eq. (13b). The condition for suppression of the reflection is $`P_{d1}=P_{d3}`$ and $`P_{d2}=P_{d4}`$. The polarization dependencies of the other reflections of $`Pn\overline{3}m`$ can be figured out analogously.

There are eight domains of $`Pa\overline{3}`$ structure. Earlier we have already considered the two basic variants of $`Pa\overline{3}`$: I $`(ii)`$, and II $`(iii)`$. The others are obtained by rotating these two variants by the angles $`\pi /2`$, $`\pi `$, and $`3\pi /2`$ about the $`z`$axis. Applying the rotations to I and II, we arrive at two series. We label the population of the series of domains by the indices $`P_{d1}^IP_{d4}^I`$, and $`P_{d1}^{II}P_{d4}^{II}`$, respectively. The two series result in two distinct scattering matrices, for example, for $`(003)`$ they are $`B`$ and $`A`$, for $`(300)`$ $`C`$ and $`B`$, etc, see Table 1. The domain pattern produces a superposition of 8 terms (amplitudes). Each term corresponds to one of the two matrices and to one of the four azimuthal angles: $`\psi `$, $`\psi +\pi /2`$, $`\psi +\pi `$, and $`\psi +3\pi /2`$. The relevant expressions are obtained in the same way which we have used to derive Eqs. (12a-c) and (13a-c). For example,

$`F_{\sigma \pi _f}(003)`$ $`=`$ $`0,`$ (14a)
$`F_{\sigma \pi _f}(003)`$ $`=`$ $`\stackrel{~}{F}\mathrm{cos}\mathrm{\Theta }(P^{}{}_{}{}^{I}\mathrm{sin}(\psi +\psi _0^I)`$ (14b)
$`+P^{}{}_{}{}^{II}\mathrm{cos}(\psi +\psi _0^{II})),`$
The latter expression can be transformed to a shifted $`sin`$ or $`cos`$ like functions, i.e. $`F_{\sigma \pi _f}(003)\mathrm{cos}\mathrm{\Theta }\mathrm{sin}(\psi +\psi _0^{})`$, where $`\psi _0^{}`$ is a phase shift depending on the domain pattern of $`Pa\overline{3}`$. $`P^{}^I`$, $`P^{}^{II}`$, and $`\psi _0^I`$, $`\psi _0^{II}`$ are given by Eq. (13c) and (13b), where $`P_{di}`$ ($`i=14`$) are replaced by $`P_{di}^I`$ or $`P_{di}^{II}`$. For the $`(300)`$ reflection of $`Pa\overline{3}`$ we find

$`F_{\sigma \sigma }(300)=\stackrel{~}{F}\stackrel{~}{P}^I\mathrm{sin}(2\psi ),`$ (15a)
$`F_{\sigma \pi _f}(300)=\stackrel{~}{F}\stackrel{~}{P}^I\mathrm{cos}(2\psi )\mathrm{sin}\mathrm{\Theta }`$
$`\stackrel{~}{F}P^{}{}_{}{}^{II}\mathrm{sin}(\psi +\psi _0^{II})\mathrm{cos}\mathrm{\Theta }.`$ (15b)
Here again $`\stackrel{~}{P}^I`$ is given by (12c) for $`P_{di}^I`$, $`i=14`$. As before, the superstructure Bragg reflections are suppressed in the case of equal population of 8 domains in the crystal.

## III Configuration interaction calculation

In this section we describe in detail how we construct the basis consisting of many determinantal wave functions for a many electron system and how we calculate the relevant matrix elements. First, we perform a band structure calculation in order to determine the conduction electron charge distribution in the muffin-tin (MT) sphere around Np. Subsequently, we apply the configuration interaction method to treat the many electron system consisting of the localized 5$`f^3`$ electrons and the conduction electron in the MT sphere.

### III.1 Electron band structure calculation

We have started by performing an electron band structure calculation of NpO2 using our linear augmented plane wave (LAPW) code. Ion  The calculation has been done assuming the muffin-tin (MT) shape of the one-electron potential and the Barth-Hedin expression for exchange, Bar  which is a variant of the local density approximation (LDA). The equal MT radii $`R_{MT}^{Np}=`$2.2206 and $`R_{MT}^O=`$2.2206 in atomic units (a.u.) were chosen for Np and O, with the cubic lattice constant $`a=10.2567`$ a.u. (or 5.4276 Å). Man  The MT potential and density of Np and O have been obtained self-consistently using a LAPW basis of $``$300 functions on a 20-point mesh of the irreducible part of the Brillouin zone. The three $`5f`$ electrons of Np were treated as core states, which adjust self-consistently to the conduction electron density.

As a result of the calculation we obtained that NpO2 is an insulator, with the energy gap $`\mathrm{}E=0.789`$ eV. The width of the occupied electron band below the Fermi level is $`E_TE_B=5.953`$ eV. The spin-orbit splitting between $`5f_{7/2}`$ and $`5f_{5/2}`$ one-electron states is $`\mathrm{}_f=`$0.983 eV. The main goal of the calculation however was the electron charge density distribution inside the neptunium MT-sphere. The calculated partial charges of different angular symmetry ($`l=03`$) are quoted in Table 3.

An important result is that on average there is approximately one conduction electron present inside the Np MT-sphere. Therefore, the localized $`5f^3`$ configuration of Np can not be considered separately from this valence electron, which can be of $`7s`$, $`7p`$ or $`6d`$ type. The instantaneous configuration at the neptunium site becomes $`7s5f^3`$, $`7p5f^3`$ or $`6d5f^3`$.

In all cases this additional electron experiences a strong Coulomb repulsion with the three localized partners. This interaction is not fully accounted for by the band structure calculations because it requires a multideterminant treatment or configuration interaction (CI). NM3  Therefore, we have to follow a different route and below we study the electron spectrum using the multipole expansion of electronic densities.

### III.2 Many electron basis states

Our method of multipole expansion of the Coulomb interaction has been used before. NM1 ; NM2 ; NM3 ; NM4  Here we formulate it in detail for the $`sf^3`$ configuration following Refs. NM3, ; NM4, . In Sec. IV we deal with the $`5f^3`$ configuration which is easily obtained from $`sf^3`$ by omitting one $`s`$ electron. In Sec. V we consider $`7s5f^3`$, $`7p5f^3`$ and $`6d5f^3`$ configurations. For the latter two cases we will describe the important differences with $`sf^3`$.

We start by considering a face centered cubic (fcc) crystal of $`N`$ Np atoms. Each atomic site possesses one 7$`s`$ and three 5$`f`$ electrons. The position vector of an electron near a crystal lattice site $`\stackrel{}{n}`$ is given by

$`\stackrel{}{R}(\stackrel{}{n})=\stackrel{}{X}(\stackrel{}{n})+\stackrel{}{r}(\stackrel{}{n}).`$ (16)
Here $`\stackrel{}{X}(\stackrel{}{n})`$ is the lattice vector which specifies the centers of the atoms (Np-nuclei) on a rigid fcc lattice. The radius vector $`\stackrel{}{r}(\stackrel{}{n})`$ is given in polar coordinates by $`(r(\stackrel{}{n}),\mathrm{\Omega }(\stackrel{}{n}))`$, where $`r`$ is the length and $`\mathrm{\Omega }=(\mathrm{\Theta },\varphi )`$ stands for the polar angles. We label the basis ket-vectors at the lattice site $`\stackrel{}{n}`$ by a single index $`I`$ or, alternatively, by four one electron indices $`(i_1^f,i_2^f,i_3^f;i^s)`$:

$`|I_\stackrel{}{n}=|i_1^f,i_2^f,i_3^f;i^s_\stackrel{}{n}.`$ (17)
The index $`i^f=(m^f,s_z^f)`$ stands for the orbital ($`m^f=17`$) and spin projection ($`s_z=\pm 1/2`$) quantum numbers of one $`f`$ electron. Therefore, there are 14 states which we label by $`i^f=114`$. Two states of the $`s`$ electron are labeled by $`i^s=1,2`$. The many electron basis wave functions are

$`\stackrel{}{r}_1,\stackrel{}{r}_2,\stackrel{}{r}_3,\stackrel{}{r}_4|I_\stackrel{}{n}={\displaystyle \frac{1}{\sqrt{N_a}}}{\displaystyle \underset{a}{}}(1)^{P(a)}`$
$`\times {\displaystyle \underset{t=1}{\overset{3}{}}}\stackrel{}{r}_t|i_t^f_\stackrel{}{n}\stackrel{}{r}_4|i^s_\stackrel{}{n},`$ (18)
where $`a`$ stands for a permutation of four electrons, the factor $`(1)^P`$ takes into account the parity of the permutation, $`N_a`$ is the number of the permutations, and

$`\stackrel{}{r}|i^f_\stackrel{}{n}=_f(r(\stackrel{}{n}))\widehat{n}|i^f,`$ (19a)
$`\stackrel{}{r}^{}|i^s_\stackrel{}{n}=_s(r^{}(\stackrel{}{n}))\widehat{n}^{}|i^s.`$ (19b)
Here $`_f`$ and $`_s`$ are the radial components of the 5$`f`$ and 7$`s`$ electrons, respectively; $`\widehat{n}`$ stands for $`\mathrm{\Omega }(\stackrel{}{n})`$. The $`5f`$ spin-orbitals can be written as

$`\widehat{n}|i^f=\widehat{n}|m_fu_s(s_z(f)),`$ (20)
Here $`u_s`$ is the spin function ($`s=\pm `$). The $`f`$orbital parts, $`\widehat{n}|m_f`$, are expressed in terms of spherical harmonics $`Y_l^m(\mathrm{\Omega })=\widehat{n}|l,m`$. We find it convenient to work with real spherical harmonics Bra  $`Y_l^\tau `$, where $`\tau =0`$, $`(m,c)`$ or $`(m,s)`$.

The order of indices in (17) is important. For example, as follows from the dynamical equivalence of the electrons the state $`|i_1^f,i_2^f,i^s,i_3^f`$ can be reduced to $`|i_1^f,i_2^f,i_3^f;i^s`$ by permuting the third and the fourth electrons, i.e.

$`|i_1^f,i_2^f,i^s,i_3^f=|i_1^f,i_2^f,i_3^f;i^s,`$ (21)
and so on. To describe the same quantum state we will use the basis vectors (17) and apply the permutation law (21) when needed. Alternatively, one can use the corresponding Slater determinants for the four electron wave functions, Eq. (18). However, the permutation relations of the type of Eq. (21) are more efficient for our purposes. Excluding equivalent states, we find only $`(141312/3!)\times 2=728`$ independent functions, or determinants, for $`7s5f^3`$. (These are 364, 2184 and 3640 for $`5f^3`$, $`7p5f^3`$, $`6d5f^3`$, respectively.) Notice, that every basis wave function is in fact a Slater determinant, Eq. (18).

### III.3 Multipole repulsion between electrons

Now we take into account the Coulomb intrasite and intersite repulsion by expanding the interaction in multipole series. As was discussed in Ref. NM3, , these interactions are treated exactly in the chosen quantum space ($`7s5f^3`$).

The Coulomb interaction between two electrons at sites $`\stackrel{}{n}`$ and $`\stackrel{}{n}^{}`$ is given by

$`V(\stackrel{}{R}(\stackrel{}{n}),\stackrel{}{R}^{}(\stackrel{}{n}^{}))={\displaystyle \frac{1}{|\stackrel{}{R}(\stackrel{}{n})\stackrel{}{R}^{}(\stackrel{}{n}^{})|}}.`$ (22)
The multipole expansion in terms of site symmetry adapted functions (SAF’s) Bra  is

$`V(\stackrel{}{R}(\stackrel{}{n}),\stackrel{}{R}^{}(\stackrel{}{n}^{}))={\displaystyle \underset{\mathrm{\Lambda }\mathrm{\Lambda }^{}}{}}v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})S_\mathrm{\Lambda }(\widehat{n})S_\mathrm{\Lambda }^{}(\widehat{n}^{}),`$
(23)
where

$`v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})={\displaystyle 𝑑\mathrm{\Omega }(\stackrel{}{n})𝑑\mathrm{\Omega }^{}(\stackrel{}{n}^{})\frac{S_\mathrm{\Lambda }(\widehat{n})S_\mathrm{\Lambda }^{}(\widehat{n}^{})}{|\stackrel{}{R}(\stackrel{}{n})\stackrel{}{R}^{}(\stackrel{}{n}^{})|}}.`$
(24)
The SAF’s are linear combinations of spherical harmonics and transform as irreducible representations of a site point group, Ref. Bra, . The index $`\mathrm{\Lambda }`$ stands for $`(l,\tau )`$, with $`\tau =(\mathrm{\Gamma },\mu ,k)`$. Here $`l`$ accounts for the angular dependence of the multipolar expansion, $`\mathrm{\Gamma }`$ denotes an irreducible representation (in the present case the group is $`O_h`$), $`\mu `$ labels the representations that occur more than once and $`k`$ denotes the rows of a given representation.

The intrasite case corresponds to $`\stackrel{}{n}=\stackrel{}{n}^{}`$. The interaction function $`v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(r,r^{})v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n}=\stackrel{}{n}^{};r,r^{})`$ then becomes particular simple,

$`v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(r,r^{})=\left({\displaystyle \frac{r_<^l}{r_>^{(l+1)}}}\right){\displaystyle \frac{4\pi }{2l+1}}\delta _{\mathrm{\Lambda }\mathrm{\Lambda }^{}},`$ (25)
where $`r_>=max(r,r^{})`$, $`r_<=min(r,r^{})`$ and $`\delta _{\mathrm{\Lambda }\mathrm{\Lambda }^{}}=\delta _{\tau \tau ^{}}\delta _{ll^{}}`$. The last expression is also site independent.

There is no simple analytical expression for the intersite case, $`\stackrel{}{n}\stackrel{}{n}^{}`$. Hei  The intersite multipole interactions are anisotropic and for practical purposes it is important to use the following dependence Hei

$`v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{}){\displaystyle \frac{(r)^l(r^{})^l^{}}{|\stackrel{}{X}(\stackrel{}{n})\stackrel{}{X}(\stackrel{}{n}^{})|^{l+l^{}+1}}}.`$ (26)
### III.4 Intrasite matrix elements

For the Coulomb interaction between four electrons on a same site $`\stackrel{}{n}`$ we have a sum of six two-body terms,

$`V^{(4)}={\displaystyle \frac{1}{2}}{\displaystyle \underset{t=1}{\overset{4}{}}}{\displaystyle \underset{p(t)=1}{\overset{4}{}}}V(\stackrel{}{r}_t,\stackrel{}{r}_p),`$ (27)
where each term is given by the multipole expansion (23). In order to calculate the matrix elements of $`V^{(4)}`$, $`i_1^f,i_2^f,i_3^f;i^s|V^{(4)}|j_1^f,j_2^f,j_3^f;j^s`$, one has to classify the electronic transitions. Following Ref. NM4,  where the energy terms of molecular ions C$`{}_{}{}^{m\pm }{}_{60}{}^{}`$, $`m=25`$, were calculated, we consider four possibilities for the fourth $`s`$electron: (1) $`i^sj^s`$, (2) $`i^sj_3^f`$, (3) $`i^sj_2^f`$, and (4) $`i^sj_1^f`$, which we label by the index $`a_4=14`$. The $`a_4=2`$ and $`a_4=4`$ transitions involve odd number of transpositions among $`j_1^f,j_2^f,j_3^f;j^s`$, and the parity is $`P(a_4=2)=P(a_4=4)=1`$. For two other transitions the number of transpositions is even and $`P(a_4=1)=P(a_4=3)=1`$. After this we are left with only three $`j`$-states, which we label as $`j_1^{}`$, $`j_2^{}`$, and $`j_3^{}`$. For the next electron, $`i_3^f`$, we can consider three possibilities ($`i_3^fj_3^{}`$, $`i_3^fj_2^{}`$, $`i_3^fj_1^{}`$) which we label by the index $`a_3=13`$. In this way we continue until we exhaust all four electrons. As a result, each subcase (or electron transition) is classified by the three index label $`a(a_4,a_3,a_2)`$, and its parity is $`P(a)=P(a_4)P(a_3)P(a_2)`$. Mathematically, we reduce a permutation of four electrons to a product of transpositions. The matrix element $`I|V^{(4)}|J`$ is found as QC49

$`I|V^{(4)}|J={\displaystyle \underset{a}{}}P(a)I|V^{(4)}|J^{(a)},`$ (28)
where $`_a=_{a_4=1}^4_{a_3=1}^3_{a_2=1}^2`$, and

$`I|V^{(4)}|J^{(a)}`$ $`=`$ $`{\displaystyle \underset{l,\tau }{}}v_l^{s,j_{a4}f,j_{a3}}c_{l,\tau }(i^sj_{a4})c_{l,\tau }(i_3^fj_{a3})`$ (29)
$`\times \delta (i_2^f,j_{a2})\delta (i_1^f,j_{a1})+p.i.`$
Here $`p.i.`$ stands for the other pair Coulomb interactions, Eq. (27). (The explicitly written term in Eq. (29) corresponds to the interaction between fourth and third electron.) The elements $`c_\mathrm{\Lambda }(ij)c_{l,\tau }(ij)`$ are defined by

$`c_\mathrm{\Lambda }(ij)={\displaystyle 𝑑\mathrm{\Omega }i|\widehat{n}S_\mathrm{\Lambda }(\widehat{n})\widehat{n}|j}.`$ (30)
For the $`7s5f^3`$ configuration there are three types of these coefficients. For the $`ss`$ transition it is only one integral $`s|Y_0^0|s=1/\sqrt{4\pi }`$, which is not zero. For the $`ff`$ transitions and real spherical harmonics, Bra  the coefficients $`c_{l,\tau }(i^fj^f)`$ were tabulated in Ref. NM1, . Finally, there are $`fs`$ and $`sf`$ transitions which require the evaluation of $`c_{l,\tau }(i^fj^s)`$. From the orthogonality of spherical harmonics we find that

$`0,0|Y_3^\tau |3,\tau ={\displaystyle \frac{1}{\sqrt{4\pi }}},`$ (31)
where $`\tau =0`$, $`(m,c)`$, $`(m,s)`$, $`m=13`$, and zero otherwise. The matrix quantities (30) were first introduced by Condon and Shortley for the description of atomic spectra, CS  but they are also at the center of the calculation of the crystal electric field effects. NM1

In Eq. (29) $`v_l^{s,j_{a4}f,j_{a3}}`$ stands for a radial average. The general expression is

$`v_l^{a,bc,e}={\displaystyle 𝑑rr^2𝑑r^{}r_{}^{}{}_{}{}^{2}_a(r)_b(r)}`$
$`_c(r^{})_e(r^{})v_l(r,r^{}),`$ (32)
where $`_a`$, $`_b`$, $`_c`$ and $`_e`$ are radial components, and $`v_l(r,r^{})=v_{\mathrm{\Lambda }\mathrm{\Lambda }}(r,r^{})`$ is given by Eq. (25). For the $`7s5f^3`$ configuration the indices $`a,b,c,e`$ refer either to $`7s`$ or to $`5f`$ electron radial components. There are only two types of radial integrals, corresponding to nontrivial multipolar terms ($`l0`$) of $`7s5f^3`$, which are

$`v_l^{ffff}`$ $`=`$ $`{\displaystyle 𝑑rr^2𝑑r^{}r_{}^{}{}_{}{}^{2}_f^2(r)_f^2(r^{})v_l(r,r^{})},`$ (33a)
$`v_l^{}^{fssf}`$ $`=`$ $`{\displaystyle 𝑑rr^2𝑑r^{}r_{}^{}{}_{}{}^{2}_s(r)_f(r)}`$ (33b)
$`_f(r^{})_s(r^{})v_l(r,r^{}).`$
Here $`l=2,4,6`$ and $`l^{}=3`$, as follows from the selection rules for $`c_\mathrm{\Lambda }(ff)`$ and $`c_\mathrm{\Lambda }^{}(fs)`$. The radial integrals $`v_l^{ffff}`$ and $`v_l^{}^{fssf}`$ are proportional to the quantities $`F_l`$ and $`G_l^{}`$ introduced by Condon and Shortley in Ref. CS, . It is important to notice that even the $`7s`$ electron with the trivial dependence of its angular part is strongly coupled to the three $`5f`$ electrons through $`fs`$ and $`sf`$ transitions.

The classification scheme for electronic transitions which we have introduced here is very useful for handling the single particle interactions as well. The main difference is that now the interaction occurs to a single electron while the rest of them produce Kronecker factors, Ref. NM4, . This group of interactions includes the spin-orbit coupling $`H_{so}`$, the crystal electric field $`V_{CF}`$ and the mean field $`V_{MF}`$. The latter two interactions are dealt with in Secs. IV and V. The spin-orbit coupling is

$`H_{so}={\displaystyle \underset{i}{}}V_{so}(i),`$ (34)
where the sum runs over all electrons, $`V_{so}`$ being the corresponding one-electron spin-orbit operator. The $`s`$electron does not experience the spin-orbit coupling and in the $`7s5f^3`$ case the summation includes only three $`5f`$ terms $`V_{so}(i^f)`$, $`i^f=13`$, where

$`V_{so}(i^f)=\zeta _f\stackrel{}{L}(i^f)\stackrel{}{S}(i^f).`$ (35)
Here $`\stackrel{}{L}(i^f)`$ and $`\stackrel{}{S}(i^f)`$ are the one-electron operators of orbital and spin momentum, $`\zeta _f`$ is the constant of the spin-orbit coupling. The full intrasite Hamiltonian is given by $`H_{intra}=V^{(4)}+H_{so}`$. It describes the $`7s5f^3`$ configuration of a free neptunium ion. Since the present method is not based on perturbation theory it extends the classical calculations of Condon and Shortley. CS

### III.5 Intersite matrix elements

We start with expression (23) and write it in the space of many electron basis vectors $`|I`$, Eq. (17). Carrying out the angular integrations $`d\mathrm{\Omega }(\stackrel{}{n})`$, $`d\mathrm{\Omega }^{}(\stackrel{}{n})`$, $`d\mathrm{\Omega }(\stackrel{}{n}^{})`$, $`d\mathrm{\Omega }^{}(\stackrel{}{n}^{})`$, we obtain

$`I|_\stackrel{}{n}I^{}|_\stackrel{}{n}^{}V(\stackrel{}{R}(\stackrel{}{n}),\stackrel{}{R}^{}(\stackrel{}{n}^{}))|J^{}_\stackrel{}{n}^{}|J_\stackrel{}{n}={\displaystyle \frac{1}{N_a}}`$
$`\times `$ $`{\displaystyle \underset{a(\stackrel{}{n})}{}}{\displaystyle \underset{a^{}(\stackrel{}{n}^{})}{}}P(a_\stackrel{}{n})P(a_\stackrel{}{n}^{}^{}){\displaystyle \underset{\alpha =1}{\overset{4}{}}}{\displaystyle \underset{\alpha ^{}=1}{\overset{4}{}}}{\displaystyle \underset{\mathrm{\Lambda }\mathrm{\Lambda }^{}}{}}v_{\mathrm{\Lambda }}^{\alpha \alpha }{}_{\mathrm{\Lambda }^{}}{}^{\alpha ^{}\alpha ^{}}(\stackrel{}{n}\stackrel{}{n}^{})`$
$`\times `$ $`\left\{c_\mathrm{\Lambda }(i_\alpha j_\alpha ^a){\displaystyle \underset{\beta =1}{\overset{3}{}}}\delta (i_\beta j_\beta ^a)\right\}\left\{c_\mathrm{\Lambda }^{}(i_\alpha ^{}^{}j^{}{}_{\alpha ^{}}{}^{a^{}}){\displaystyle \underset{\beta ^{}=1}{\overset{3}{}}}\delta (i_\beta ^{}^{}j^{}{}_{\beta ^{}}{}^{a^{}})\right\}.`$
Here the sum over $`a`$ means the summation over all permutations of indices $`j_1^f,j_2^f,j_3^f,j^s`$ at site $`\stackrel{}{n}`$ transforming them to indices $`j_\alpha ^a`$ ($`\alpha =14`$). Analogously, the sum over $`a^{}`$ implies the summation over all permutations of $`j_{}^{}{}_{1}{}^{f},j_{}^{}{}_{2}{}^{f},j_{}^{}{}_{3}{}^{f},j_{}^{}{}_{}{}^{s}`$ at site $`\stackrel{}{n}^{}`$ transforming them to $`j^{}_\alpha ^a^{}`$ ($`\alpha =14`$). $`P(a)`$ and $`P(a^{})`$ stand for the parities of the permutations. Indices $`\alpha `$ and $`\alpha ^{}`$ indicate which electron at site $`\stackrel{}{n}`$ interacts with which electron at site $`\stackrel{}{n}^{}`$. The other electrons labeled by $`\beta =13`$ at site $`\stackrel{}{n}`$ and by $`\beta ^{}=13`$ at site $`\stackrel{}{n}^{}`$ do not contribute to the interaction and produce the Kronecker delta symbols. The coefficients $`c_\mathrm{\Lambda }`$ are defined by Eq. (30), and the intersite ($`\stackrel{}{n}\stackrel{}{n}^{}`$) interaction element $`v_{\mathrm{\Lambda }}^{\alpha \alpha }{}_{\mathrm{\Lambda }^{}}{}^{\alpha ^{}\alpha ^{}}`$ is given by

$`v_{\mathrm{\Lambda }}^{\alpha \alpha }{}_{\mathrm{\Lambda }^{}}{}^{\alpha ^{}\alpha ^{}}(\stackrel{}{n}\stackrel{}{n}^{})={\displaystyle 𝑑rr^2𝑑r^{}r_{}^{}{}_{}{}^{2}}`$
$`\times _\alpha ^2(r)_\alpha ^{}^2(r^{})v_{\mathrm{\Lambda }\mathrm{\Lambda }^{}}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{}).`$ (37)
For the four electron space of $`7s5f^3`$ only even $`l`$ and $`l^{}`$ (in $`\mathrm{\Lambda },\mathrm{\Lambda }^{}`$) are retained in Eqs. (LABEL:2.19) and (37), and $`l,l^{}=0,2,4,6`$. Two very important examples of intersite Coulomb interactions, namely, the crystal electric field and the mean field will be considered in sections IV and V.

## IV Crystal and mean field of the $`5f^3`$ neptunium configuration

In this section we study the model where we assume that there are only three localized $`5f`$ electrons at each neptunium site. Although we believe that the model is not adequate for a realistic description of NpO2, especially in the part concerning the loss of the magnetic moments in the ordered phase, it is nevertheless very instructive to consider it in detail. The $`5f^3`$ configuration being relatively simple offers an opportunity to be studied thoroughly and to understand the interplay between the disordered and ordered phases, or between the crystal and quadrupolar mean field. On the other hand, the configurations $`7p5f^3`$ and $`6d5f^3`$ involve too many basis states and consume to much time to be processed self-consistently for any temperature.

### IV.1 Free ion electron energy spectrum

The basis states for the $`5f^3`$ configuration are given by

$`|I=|i_1^f,i_2^f,i_3^f,`$ (38)
where as before, $`i^f=114`$. The total number of basis vectors is $`14\times 13\times 12/3!=364`$. We treat the $`5f^3`$ configuration in the way which was specified in Sec. III.

There are only $`ff`$ transitions described by four radial integrals (33a): $`v_l^{ffff}`$, $`l=0,2,4,6`$. The others are zero due to the selection rules imposed by the coefficients $`c_\mathrm{\Lambda }(ij)`$, Eq. (30). The radial integral $`v_0^{ffff}`$ (Hubbard $`U`$) is not important here since it does not result in term splittings. In the following we will use the condense notation $`F`$ for $`ff`$ and thus $`v_l^{ffff}v_l^{FF}`$. These quantities are connected with the Slater (Condon-Shortley) parameters CS  $`F^l(5f,5f)`$ through the following relation:

$`v_l^{FF}={\displaystyle \frac{4\pi }{2l+1}}F^l.`$ (39)
In particular, the Slater parameters $`F^2`$, $`F^4`$, and $`F^6`$ of Amoretti et al. Amo  correspond to $`v_2^{FF}=14.007`$ eV, $`v_4^{FF}=7.091`$ eV and $`v_6^{FF}=3.168`$ eV. (In order to obtain the exact term splitting quoted in Table 4 of Ref. Amo,  we had to scale their Slater parameters by a factor 0.9755.) Alternatively, the quantities $`v_l^{FF}`$ can be calculated by using the radial dependence of the $`5f`$ electrons $`_f`$, Eq. (33a). We have done such a calculation and then corrected the parameters by comparing the splittings of the $`f^3`$ configuration with experimental data for Pr3+ and Nd4+ (details are given in Appendix B). We arrived at

$`v_2^{FF}=18.164,v_4^{FF}=8.578,v_6^{FF}=3.362,`$
$`\zeta _f=0.2547,\text{in eV}.`$ (40)
After calculating the matrix elements of the Coulomb repulsion and the spin-orbit coupling, we diagonalize the matrix

$`H_{intra}=V^{(3)}+H_{so}`$ (41)
and obtain the electronic spectrum of $`5f^3`$. The 9 lowest and two highest eigenvalues are shown in Table 4, where for comparison we also quote the spectrum of free Np ion of Ref. Amo, .

Notice, that in comparison with the spectrum Amo  of Amoretti at al., $`{}_{}{}^{4}F_{3/2}^{}`$ and $`{}_{}{}^{4}F_{5/2}^{}`$ are higher than $`{}_{}{}^{4}I_{13/2}^{}`$ and $`{}_{}{}^{2}H_{9/2}^{}`$, which is in better agreement with the sequence of terms for the $`4f^3`$ configurations of Pr3+ and Nd4+ known from experiment. La\_sp

### IV.2 CEF excitations in the disordered phase ($`T>25`$ K)

In the disordered phase ($`T>25`$ K) the electron density of the $`5f^3`$ configuration adopts the cubic ($`O_h`$) site symmetry. This density modulation is induced by the cubic crystal electric field experienced by three $`5f`$ electrons. In terms of the multipole intersite expansion (23) it implies that for a given Np site $`\stackrel{}{n}`$ we treat 12 Np neighbors ($`\stackrel{}{n}_1^{}=112`$) and 8 oxygen neighbors ($`\stackrel{}{n}_2^{}=18`$) in spherical approximation, i.e. $`l^{}=0`$ and $`S_\mathrm{\Lambda }^{}(\stackrel{}{n}^{})`$ reduces to $`Y_0^0=1/\sqrt{4\pi }`$. The coefficients $`c_\mathrm{\Lambda }^{}`$, Eq. (30), become simple,

$`c_0(i_\alpha ^{}j_\alpha ^{})={\displaystyle \frac{1}{\sqrt{4\pi }}}\delta (i_\alpha ^{},j_\alpha ^{}).`$ (42)
Here we write $`0`$ for $`\mathrm{\Lambda }^{}(l^{}=0,A_{1g})`$. At the central site $`\stackrel{}{n}`$ we expand the CEF in terms of SAF’s $`S_{\mathrm{\Lambda }_1}(\stackrel{}{n})`$, $`\mathrm{\Lambda }_1(l,A_{1g})`$, where $`A_{1g}`$ stands for the unit representation of the cubic site group $`O_h`$. The selection rules for the coefficients $`c_\mathrm{\Lambda }(i_\alpha ,j_\alpha )`$ of the $`ff`$ transitions imply that there remain only two nontrivial functions $`S_{\mathrm{\Lambda }_1}`$ with $`l=4`$ and $`l=6`$, which correspond to the cubic harmonics $`K_4(\mathrm{\Omega })`$ and $`K_6(\mathrm{\Omega })`$. The multipole two-center expansion (23) becomes

$`V(\stackrel{}{R}(\stackrel{}{n}),\stackrel{}{R}^{}(\stackrel{}{n}^{}))={\displaystyle \frac{1}{\sqrt{4\pi }}}{\displaystyle \underset{\mathrm{\Lambda }_1}{}}v_{\mathrm{\Lambda }_10}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})S_{\mathrm{\Lambda }_1}(\stackrel{}{n}),`$
(43a)
where
$`v_{\mathrm{\Lambda }_10}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})=`$
$`{\displaystyle \frac{1}{\sqrt{4\pi }}}{\displaystyle 𝑑\mathrm{\Omega }(\stackrel{}{n})𝑑\mathrm{\Omega }^{}(\stackrel{}{n}^{})\frac{S_{\mathrm{\Lambda }_1}(\widehat{n})}{|\stackrel{}{R}(\stackrel{}{n})\stackrel{}{R}^{}(\stackrel{}{n}^{})|}}.`$ (43b)
Here $`v_{\mathrm{\Lambda }_10}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})`$ has a same value for all 12 Np neighbors ($`\stackrel{}{n}_1^{}=112`$), and a same value for all 8 oxygen neighbors ($`\stackrel{}{n}_2^{}=18`$). As follows from Eq. (26) $`v_{\mathrm{\Lambda }_10}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})`$ is independent of $`r^{}`$. Eq. (37) then can be written in the following form:

$`v_{\mathrm{\Lambda }_1}^{\alpha \alpha }{}_{0}{}^{\alpha ^{}\alpha ^{}}(\stackrel{}{n}\stackrel{}{n}^{})=v_{\mathrm{\Lambda }_1}^{\alpha \alpha }{}_{0}{}^{}Q_\alpha ^{},`$ (44)
where

$`v_{\mathrm{\Lambda }_1}^{\alpha \alpha }{}_{0}{}^{}={\displaystyle }drr^2_\alpha ^2(r)v_{\mathrm{\Lambda }_1\mathrm{\hspace{0.17em}0}}(\stackrel{}{n},\stackrel{}{n}^{};r,\mathit{}^{})`$ (45a)
and
$`Q_\alpha ^{}={\displaystyle 𝑑r^{}r_{}^{}{}_{}{}^{2}_\alpha ^{}^2(r^{})}.`$ (45b)
Here the integrations are taken over $`0<r^{}<R_{MT}`$, where $`R_{MT}`$ is the radius of the muffin-tin sphere of neptunium or oxygen. (The influence of the interstitial region will be discussed later.) $`Q_\alpha ^{}`$ refers to an electron at site $`\stackrel{}{n}^{}`$ which interacts with one of the three $`5f`$ electrons at $`\stackrel{}{n}`$. We then can perform a summation over all electrons at $`\stackrel{}{n}^{}`$ and include also to this term the interaction with the nucleus. This results in replacing $`Q_\alpha ^{}`$ by $`eQ_{MT}`$ in Eq. (44), $`Q_{MT}`$ and $`e`$ being the total charge inside the MT sphere and the electron charge ($`e=1`$), respectively. From Eq. (26) it follows that

$`v_{\mathrm{\Lambda }_1}^{\alpha \alpha }{}_{0}{}^{}=v_{\mathrm{\Lambda }_1\mathrm{\hspace{0.17em}0}}(\stackrel{}{n},\stackrel{}{n}^{};R_{MT},\mathit{}^{}){\displaystyle \frac{q_l^\alpha }{R_{MT}^l}},`$ (46a)
where $`l`$ in the index $`\mathrm{\Lambda }_1`$ is 4 or 6, and
$`q_l^\alpha ={\displaystyle 𝑑r^{}r_{}^{}{}_{}{}^{(l+2)}_\alpha ^2(r^{})}.`$ (46b)
Therefore, the CEF operator for any neptunium site ($`r<R_{MT}^{Np}`$) can be written explicitly as

$`V_{CF}(\stackrel{}{R}(\stackrel{}{n}))={\displaystyle \underset{l=4,6}{}}B_lS_{(l,A_{1g})}(\widehat{n})\left({\displaystyle \frac{r}{R_{MT}^{Np}}}\right)^l,`$ (47a)
where
$`B_l=B_l^{Np}+B_l^O,`$ (47b)
and
$`B_l^{Np}={\displaystyle \frac{12}{\sqrt{4\pi }}}Q_{eff}^{Np}ev_{\mathrm{\Lambda }_1\mathrm{\hspace{0.17em}0}}(\stackrel{}{n},\stackrel{}{n}_1^{};R_{MT}^{Np},\mathit{}^{}),`$ (47c)
$`B_l^O={\displaystyle \frac{8}{\sqrt{4\pi }}}Q_{eff}^Oev_{\mathrm{\Lambda }_1\mathrm{\hspace{0.17em}0}}(\stackrel{}{n},\stackrel{}{n}_2^{};R_{MT}^{Np},\mathit{}^{}).`$ (47d)
We quote all relevant parameters of CEF in Table 5.

As given by Eq. (47a), the CEF operator $`V_{CF}`$ is a one-electron quantity. New1 ; New2  CEF acts along with the Coulomb intrasite repulsion, Eq. (41). Therefore, the total Hamiltonian for the disordered phase becomes

$`H^{dis}(\stackrel{}{n})=H_{intra}+V_{CF}(\stackrel{}{n}).`$ (48)
Although we have considered CEF from first principles there is still an ambiguity related to the charge distribution in the interstitial region. A more rigorous treatment of the problem is given in Refs. Wei, ; ND, . A careful consideration of the problem based on the solution of a periodic Poisson’s equation leads to a renormalization of the charges inside the MT spheres, Ref. ND, . In other words, in Eqs. (47a-d) the effective charges for Np and O are given by

$`Q_{eff}=Q_{MT}{\displaystyle \frac{4\pi R_{MT}^3}{3}}\rho _I(\stackrel{}{K}=0)`$
$`4\pi R_{MT}^2{\displaystyle \underset{\stackrel{}{K}0}{}^{}}{\displaystyle \frac{j_1(KR_{MT})}{K}}\rho _I(\stackrel{}{K}),`$ (49)
where $`\rho _I(\stackrel{}{K})`$ is the Fourier series expansion of the electron density in the interstitial region, $`j_{l=1}`$ is the spherical Bessel function. $`\rho _I(\stackrel{}{K}=0)`$ is the average density in the interstitial region, $`Q_{out}/V_{out}`$, where $`Q_{out}`$ and $`V_{out}`$ is the charge and volume of the interstitial region.

The calculation of $`Q_{eff}`$ according to Eq. (49) is quite laborious since it requires the evaluation of the Fourier coefficients $`\rho _I(\stackrel{}{K})`$. Instead, below we consider two approximation to (49). In the first approximation we assume that $`Q_{eff}^{Np}(I)=Q_{MT}^{Np}`$ and $`Q_{eff}^O(I)=Q_{MT}^O`$, where $`Q_{MT}^{Np}`$ and $`Q_{MT}^O`$ are the total charges inside the MT spheres of neptunium and oxygen. Here the electron charge in the interstices is completely ignored. In the second approximation we take

$`Q_{eff}(II)=Q_{MT}{\displaystyle \frac{4\pi R_{MT}^3}{3}}\rho _I(\stackrel{}{K}=0).`$ (50)
This expression corresponds to the homogeneous electron density distribution in the interstitial region. However, the modification of effective charges in this approximation is too strong: $`Q_{eff}^{Np}(I)=+4.108|e|`$ and $`Q_{eff}^{Np}(II)=+5.337|e|`$ for neptunium, $`Q_{eff}^O(I)=0.284|e|`$ and $`Q_{eff}^O(II)=+0.945|e|`$ for oxygen. As was discussed in Refs. NM5, ; NM3,  the CEF splitting is overestimated in the second approximation. In reality the charge density in interstices is highly inhomogeneous, concentrated mainly in the proximity to oxygen and neptunium. This leads to $`\rho _I(\stackrel{}{K}0)0`$, and the last term in Eq. (49) acts in the opposite direction, decreasing $`Q_{eff}`$ backward to $`Q_{MT}`$ values, which correspond to the first approximation.

The exact calculation of $`Q_{eff}`$ according to Eq. (49) is beyond the scope of the present study. Instead, we have studied the crystal field effects as a function of $`Q_{eff}`$ by introducing

$`Q_{eff}^{Np}(x_{eff})=Q_{MT}^{Np}+x_{eff}(Q_{eff}^{Np}(II)Q_{MT}^{Np}),`$ (51a)
$`Q_{eff}^O(x_{eff})=Q_{MT}^O+x_{eff}(Q_{eff}^O(II)Q_{MT}^O),`$ (51b)
where $`0<x_{eff}<1`$. Diagonalizing $`H^{dis}`$, Eq. (48), we have found that 41 terms of $`5f^3`$ are split into 120 distinct sublevels of $`\mathrm{\Gamma }_6`$, $`\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_8`$ symmetry of the cubic double group $`O_h^{}`$. In particular two lowest atomic-like levels are split according to the following scheme,

$`{}_{}{}^{4}I_{9/2}^{}\mathrm{\Gamma }_8+\mathrm{\Gamma }_8+\mathrm{\Gamma }_6,`$ (52a)
$`{}_{}{}^{4}I_{11/2}^{}\mathrm{\Gamma }_8+\mathrm{\Gamma }_7+\mathrm{\Gamma }_6+\mathrm{\Gamma }_8.`$ (52b)
The resulting splittings and the dependence of CEF on $`x_{eff}`$ is shown in Fig. 2. The splittings of two lowest terms of the $`5f^3`$ configuration, Eq. (52a,b), is also given in Tables 6 and 7 for $`x_{eff}=0`$ and $`x_{eff}=0.5`$, respectively.

The most comprehensive study of the crystal field of the $`5f^3`$ configuration was performed by Amoretti et al., Ref. Amo, . Comparing their results with ours, we obtain the following relations connecting $`B_4`$ and $`B_6`$ with $`V_4`$ and $`V_6`$ used there,

$`B_4`$ $`=`$ $`8\sqrt{{\displaystyle \frac{12}{7}}}V_4,`$ (53a)
$`B_6`$ $`=`$ $`16\sqrt{8}V_6.`$ (53b)
We observe that for a realistic choice of $`Q_{eff}`$, which corresponds to $`x_{eff}00.5`$, the CEF splitting is a few times smaller than the value 55 meV considered for CEF excitations in Ref. Amo, . Correspondingly, the calculated parameters $`B_4`$ and $`B_6`$ (Tables 6 and 7) are smaller. We discuss a possible explanation to this fact in the Conclusions, Sec. VI. Notice that it is not possible to relate the feature at 55 meV with the $`{}_{}{}^{4}I_{11/2}^{}`$ splittings because it is situated at much higher energy $`650`$ meV, Table 4. Most likely, the experimental excitations at 55 meV refer to the valence electrons delocalized on the Np-O bonds, while the lowest CEF excitations of Np lie at a smaller energy $`6`$ meV, Tables 6, 7.

Finally, we would like to notice that our CEF calculation is based on first principles, and in that respect it differs from the others, Amo ; Lea ; Ste  which use fitting from experiment to extract the CEF parameters. There are also other technical differences. In contrast to Refs. Lea, ; Ste,  we do not assume that the full momentum $`J`$ is a good quantum number, that allows for a mixing of components belonging to different $`J`$. In Ref. Amo,  the basis was truncated to the first 11 low-lying levels. In our approach we do not have these limitations.

### IV.3 Mean field and the structural phase transition at 25 K

Now we consider the intersite quadrupole interactions $`V^{QQ}(\stackrel{}{n},\stackrel{}{n}^{})`$ between a central Np site (sublattice $`\{n_1\}`$) and its 12 nearest Np neighbors belonging to sublattices $`\{n_p^{}\}`$, $`p^{}=2,3,4`$. (Here we will not take into account the interactions involving higher spherical harmonics because they are considerably smaller, Eq. (26).) We will handle the $`Pn\overline{3}m`$ and $`Pa\overline{3}`$ symmetries together since both of them result in the same local quadrupolar function (1a) for neptunium.

The problem for the $`Pa\overline{3}`$ structure has been considered in Refs. NM1, ; NM3,  where it is shown that this spatial order of quadrupoles gives an effective attraction between them. This direct quadrupole-quadrupole coupling can be calculated from first principles. However, we will see that its strength is not sufficient and we have to assume a substantial reinforcement via oxygen mediated interaction. For the $`Pn\overline{3}m`$ structure the direct interaction is repulsive and we have to resort to the indirect coupling from the very beginning. The important fact which we exploit in this section is that the rhombohedral (trigonal) mean field for both structures can be described by the same expression (59) below.

We start by deriving an effective mean field for direct Coulomb coupling between quadrupoles in the $`Pa\overline{3}`$ structure. Following Ref. NM3, , one obtains from Eq. (LABEL:2.19) the following expression for the quadrupolar interaction operator between neptunium sites at $`\stackrel{}{n}_1`$ and $`\stackrel{}{n}_p^{}`$

$`V^{QQ}(\stackrel{}{n}_1,\stackrel{}{n}_p^{})={\displaystyle \frac{\gamma ^{ff}}{3}}\rho _f^Q(\stackrel{}{n}_1)\rho _f^Q(\stackrel{}{n}_p^{}).`$ (54)
Here $`\stackrel{}{n}_1\{n_1\}`$, $`\stackrel{}{n}_p^{}\{n_p^{}\}`$, and the quadrupolar density operator $`\rho _f^Q`$ at site $`\stackrel{}{n}_p`$ is given by

$`\rho _f^Q(\stackrel{}{n})={\displaystyle \underset{I,J}{}}|I{\displaystyle \underset{a}{}}P(a){\displaystyle \underset{\alpha =1}{\overset{3}{}}}c_p(i_\alpha j_\alpha ^a){\displaystyle \underset{\beta =1}{\overset{2}{}}}\delta _{i_\beta j_\beta ^a}J|,`$
(55)
where again $`a`$ is a permutation of $`j_1^f,j_2^f,j_3^f`$ transforming them to a new order given by $`j_\kappa ^a`$, $`\kappa =13`$. $`P(a)`$ is the parity of the permutation, the index $`\alpha `$ stands for the interacting electron at the site. (The second interacting electron belongs to a neighboring Np site.) The other (noninteracting) electrons ($`\beta =1,2`$) at $`\stackrel{}{n}`$ produce the product of the Kronecker delta symbols. The quadrupolar $`ff`$ coefficients are defined as

$`c_p(i_\alpha j_\alpha )=i_\alpha |𝒮_p|j_\alpha .`$ (56)
There are four types of such coefficients (i.e. $`p=14`$) as follows from Eqs. (1a-d). Finally, the interaction constant $`\gamma ^{ff}`$ in Eq. (54) is given by

$`\gamma ^{ff}={\displaystyle 𝑑rr^2𝑑r^{}r_{}^{}{}_{}{}^{2}_f^2(r)_f^2(r^{})v_{\mathrm{\Lambda }\mathrm{\Lambda }}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})},`$
(57)
with $`v_{\mathrm{\Lambda }\mathrm{\Lambda }}(\stackrel{}{n},\stackrel{}{n}^{};r,r^{})`$ defined by Eq. (24), where $`\stackrel{}{n}=(0,0,0)`$, $`\stackrel{}{n}^{}=(a/2)(0,1,1)`$ and $`\mathrm{\Lambda }=(l=2,T_{2g},k=1)`$. The corresponding SAF is $`S_\mathrm{\Lambda }=Y_2^{1s}`$. Using the property (26) we rewrite $`\gamma ^{ff}`$ as

$`\gamma ^{ff}={\displaystyle \frac{q_2^f}{(R_{MT}^{Np})^2}}v^{QQ}(R_{MT}^{Np}){\displaystyle \frac{q_2^f}{(R_{MT}^{Np})^2}},`$ (58)
where the short notation $`v^{QQ}(R_{MT}^{Np})`$ stands for $`v_{\mathrm{\Lambda }\mathrm{\Lambda }}(\stackrel{}{n},\stackrel{}{n}^{};R_{MT}^{Np},R_{MT}^{Np})`$, and the “quadrupole charge” $`q_{l=2}^f`$ of $`5f`$ electron is given by Eq. (46b).

In the mean-field approximation after summing over 12 nearest neighbors belonging to 3 sublattices we arrive at the effective bilinear quadrupole-quadrupole operator

$`U^{QQ}(\stackrel{}{n}_p)=\lambda ^{ff}\rho _f^Q\rho _f^Q(\stackrel{}{n}_p),`$ (59)
where $`\lambda ^{ff}=4\gamma ^{ff}>0`$ and the calculated parameters are quoted in Table 8.

The same expression holds for $`Pn\overline{3}m`$, but in that case for the direct quadrupole interaction $`\lambda ^{ff}<0`$ (repulsion), Appendix A. $`\rho _f^Q`$ stands for an expectation value of the quadrupole operator (55). At zero temperature it is the quantum average

$`\rho _f^Q=I_{gs}|\rho _f^Q(\stackrel{}{n}_p)|I_{gs},`$ (60)
where $`|I_{gs}`$ refers to the ground state of the full mean-field Hamiltonian

$`H^{MF}(\stackrel{}{n})=U^{QQ}(\stackrel{}{n})+V_{CF}(\stackrel{}{n})+H_{intra}(\stackrel{}{n}).`$ (61)
The intrasite part of the interactions $`H_{intra}`$ is given by Eq. (41). For CEF we used the values $`x_{eff}=0`$ and $`Q_{eff}^{Np}=Q_{MT}^{Np}`$, $`Q_{eff}^O=Q_{MT}^O`$. (As we discussed in Sec. IV.B this gives the most realistic estimate for CEF.) For both $`Pa\overline{3}`$ and $`Pm\overline{3}n`$ structures the mean-field Hamiltonian has the $`S_6=C_3\times i`$ trigonal common site symmetry, and lifts the 4-fold degeneracy of the quartet states,

$`\mathrm{\Gamma }_8E+E,`$ (62)
where $`E`$ stands for the two-fold degenerate irreducible representation of $`C_3`$. The states $`\mathrm{\Gamma }_6`$ remain unsplit as a consequence of the Kramers theorem, $`\mathrm{\Gamma }_6E`$.

We now obtain a system of equations which can be solved self-consistently. At first we introduce an average $`\rho _f^Q`$, which defines the interactions (59), (61). In the space of the $`5f^3`$ configuration (364 state vectors), we diagonalize the total Hamiltonian $`H^{MF}`$, Eq. (61), and obtain the eigenvectors $`|K`$, $`K=1364`$:

$`H^{MF}|K=ϵ_K|K,`$ (63)
where the lowest value of $`ϵ_K`$ corresponds to $`K=1,2`$. This is the Kramers doublet of the ground state. We then calculate the quantities

$`\rho _f^Q(K)=K|\rho _f^Q(\stackrel{}{n}_p)|K,`$ (64)
which evaluate the quadrupolar moments of the states $`K`$. Next, we find an improved value for $`\rho _f^Q`$ which is $`\rho _f^Q(K=1)=\rho _f^Q(K=2)`$, Eq. (60). The procedure continues until the input and output values for $`\rho _f^Q`$ converge. As a consequence of symmetry the expectation value $`\rho _f^Q`$ is independent of the sublattice $`\{n_p\}`$ chosen for calculations, but the Hamiltonian and eigenvectors do depend on the choice. This is because the quadrupoles have different orientations for different sublattices.

For nonzero temperature $`T`$ the mean-field equation for $`\rho _f^Q`$ becomes

$`\rho _f^Q={\displaystyle \frac{1}{Z}}{\displaystyle \underset{K}{}}\rho _f^Q(K)e^{ϵ_K/T},`$ (65a)
where
$`Z={\displaystyle \underset{K}{}}e^{ϵ_K/T}.`$ (65b)
The results of the calculations are quoted in Tables 9, 10, and in Figs. 3, 4, 5. With $`\lambda ^{ff}=`$1750.5 K, Table 8, which corresponds to the direct quadrupole-quadrupole attraction of the $`Pa\overline{3}`$ structure, we found that the transition temperature $`T_Q`$ is only 0.44 K. A typical dependence of $`\rho _f^Q`$ is shown in Fig. 3. The phase transition is of first order, with a discontinuity of the order parameter amplitude $`\rho _f^Q=0.0075`$ at $`T_Q`$. Comparing the present calculation with that for cerium, NM3  we observe that the low value of $`T_Q`$ is due to a small quadrupolar susceptibility of the ground state quartet $`\mathrm{\Gamma }_8`$, since the difference in $`\lambda `$ is not that much (for cerium $`\lambda ^{ff}=2241`$ K.)

The calculated transition temperature is very small in comparison with the experimental value of 25 K. We then conclude that in the framework of the $`5f^3`$ model the structural phase transition can not be explained by means of the direct bilinear quadrupole coupling. Most likely, there is indirect coupling between $`5f`$ electron densities on neptunium via the oxygen atoms (superexchange interaction super ). The microscopic consideration of the superexchange interaction is beyond the scope of the present work. Instead, we model it by increasing the value of $`\lambda ^{ff}`$. In such a case $`\lambda ^{ff}`$ becomes a phenomenological parameter which we denote as $`\lambda _{eff}`$. By changing $`\lambda _{eff}`$ we change the transition temperature as shown in Fig. 4. We have found that the experimental value of 25 K is achieved for $`\lambda _{eff}/\lambda ^{ff}10`$, which indicates a substantial increase of the effective bilinear coupling, Eq. (59). The relevant parameters of such strong mean-field are given in Table 10.

Finally, we would like to mention that the mean-field calculations have been done assuming that the CEF is weak, i.e. $`x_{eff}=0`$. Increasing $`x_{eff}`$ leads to an increase of the CEF splittings, which results in a strong suppression of the transition temperature $`T_Q`$, Fig. 5. Notice that at $`x_{eff}(\mathrm{\Gamma }_8\mathrm{\Gamma }_6)=0.241`$ the ground state changes to the $`\mathrm{\Gamma }_6`$ doublet, Fig. 2, which is apparently unfavorable for the quadrupolar order. With further increase of $`x_{eff}`$ beyond the $`x_{eff}(\mathrm{\Gamma }_8\mathrm{\Gamma }_6)`$ point the transition temperature goes fast to zero and the ordered phase disappears.

## V Four electron configurations at neptunium site

As follows from the previous section the model which takes into account only the $`5f^3`$ configuration at each neptunium site is not capable to explain the disappearance of the magnetic moments in the ordered phase. On the other hand, the charge distribution inside the MT sphere centered at the neptunium nucleus indicates that there is always approximately one valence electron instantaneously present along with the three localized $`5f`$ electrons, Table 3. Even if the electron is in the $`s`$-orbital state, it experiences strong coupling with the $`5f`$ electrons via intrasite $`sf`$ transitions. Therefore, the excitation spectrum at each neptunium site differs from that for $`5f^3`$ considered in Sec. IV. In the following we model the couplings with the valence electron by considering $`7s5f^3`$, $`7p5f^3`$, and $`6d5f^3`$ instantaneous configurations. Here we will not study the mean-field in such detail as for the $`5f^3`$ configuration. Our main objective is to show that the ground state can be nonmagnetic and separated from the magnetic excitations by an energy gap larger than 25 K.

### V.1 $`7s5f^3`$ configuration

The $`7s5f^3`$ configuration has been considered in detail in Sec. III. We are working here in the space of 728 state vectors $`|I=|i_1^f,i_2^f,i_3^f;i^s`$. The $`5f`$-states are coupled to two $`s`$ states through the $`fs`$ transitions accompanied by the multipole Coulomb interactions with $`l=3`$. The strength of this interaction was estimated from the LDA calculation of a Np atom, Eq. (33b),

$`v_3^{fssf}=1.313\text{eV}.`$ (66)
The electron energy spectrum of the $`7s5f^3`$ configuration consists of 82 distinct levels, see Table 11.

In the cubic crystal field two lowest levels ($`{}_{}{}^{5}I_{4}^{}`$ and $`{}_{}{}^{5}I_{5}^{}`$) are split as quoted in Table 12.

Notice that the ground state is nonmagnetic. However, at higher temperatures ($`T>25`$ K) two low lying excitations of $`T_2`$ and $`T_1`$ symmetry contribute to the Curie law for the magnetic susceptibility. Erd2

In the ordered phase the local symmetry of the Np site changes to $`S_6`$ ($`Pa\overline{3}`$) or $`D_{3d}`$ ($`Pn\overline{3}m`$). As we discussed in Sec. IV.C the mean-field can be expanded in a multipole series and from symmetry it follows that the strongest interaction is of quadrupolar type,

$`U^{QQ}(\stackrel{}{n}_p)=\mathrm{\Lambda }^f\rho _f^Q(\stackrel{}{n}_p),`$ (67)
where $`\mathrm{\Lambda }^f=\lambda ^{ff}\rho _f^Q`$ \[compare with Eq. (59)\]. The quadrupolar operator $`\rho _f^Q(\stackrel{}{n}_p)`$ belonging to the sublattice $`\{n_p\}`$ is given again by Eq. (55) with the corresponding orientational function $`𝒮_p`$, Eqs (1a-d). The only difference is that the $`s`$electron produces an additional Kronecker symbol, i.e. in Eq. (55) $`\beta =13`$. Notice that both CEF, Eq. (47a), and the mean field, Eq. (67), act only on the $`5f`$ electrons. Taking $`\mathrm{\Lambda }^f=5252`$ K we diagonalized the full Hamiltonian

$`H^{MF}(\stackrel{}{n})=U^{QQ}(\stackrel{}{n})+V_{CF}(\stackrel{}{n})+H_{intra}(\stackrel{}{n}).`$ (68)
The resultant electronic spectrum is shown in Table 13.

We observe that the ground state is nonmagnetic, while the first magnetic excitation ($`E`$) does not contribute to the magnetic susceptibility if $`T<25`$ K.

### V.2 $`7p5f^3`$ configuration

In case of the $`7p5f^3`$ configuration, we construct $`(14\times 13\times 12/3!)\times 6`$=2184 basis vectors

$`|I=|i_1^f,i_2^f,i_3^f;i^p.`$ (69)
As before, $`i^f`$ refers to a $`5f`$ electron, $`i^p=(k,s_z)`$ to the $`7p`$ electron ($`k`$ is its orbital part, $`k=13`$, and $`s_z`$ is the spin part).

The interaction between $`5f`$ electrons was described in detail earlier. In addition, $`pp`$ and $`ff`$ transitions lead to intrasite multipole interactions with the $`l=0`$ and $`l=2`$ components, $`fp`$ transitions give the Coulomb multipole couplings with $`l=2`$ and $`l=4`$. The relevant Slater parameters were extracted from the radial dependences of a Np ion,

$`v_2^{fp}=1.442,v_4^{fp}=0.739,`$
$`v_2^{ffpp}=5.237\text{(in eV)}.`$ (70)
Also, the $`p`$-electron experiences the spin-orbit interaction with $`\zeta _p=1.2795`$ eV.

First, we calculated the electron spectrum of the free ion, $`H_{intra}`$ and obtained 242 distinct levels, Table 14.

In the cubic crystal the $`{}_{}{}^{5}K_{4}^{}`$ and $`{}_{}{}^{5}I_{4}^{}`$ levels are split as quoted in Table 15. (The procedure of treating CEF effects is outlined in Sec. IV.B, CEF does not act on the $`p`$electron.)

Notice, that the ground state now is a magnetic triplet of the $`T_2`$ symmetry, which together with two other magnetic $`T_1`$ excitations, Table 15, gives the Curie law for the magnetic susceptibility.

At $`T_c=25`$ K the structural phase transition occurs and the symmetry of the Np sites is reduced. This symmetry change is accompanied by lifting degeneracies of some cubic levels. In particular, the ground state $`T_2`$ triplet is split in a doublet and a single level as demonstrated in Table 16. The mean-field is approximated by its quadrupolar electric part,

$`U^{QQ}(\stackrel{}{n}_i)=\mathrm{\Lambda }^f\rho _f^Q(\stackrel{}{n}_i)\mathrm{\Lambda }^p\rho _p^Q(\stackrel{}{n}_i),`$ (71)
where

$`\mathrm{\Lambda }^f=\lambda ^{ff}\rho _f^Q+\lambda ^{fp}\rho _p^Q,`$ (72a)
$`\mathrm{\Lambda }^p=\lambda ^{pf}\rho _f^Q+\lambda ^{pp}\rho _p^Q.`$ (72b)
The following parameters of the interaction were assumed,

$`{\displaystyle \frac{q_2^p}{(R_{MT}^{Np})^2}}=0.1604,\mathrm{\Lambda }^f=6612\text{K},\mathrm{\Lambda }^p=3426\text{K}.`$
(73)
Since the charge density expansion of the $`p`$-electron has quadrupolar components, $`c_k(i^pj^p)=i^p|𝒮_k|j^p0`$ ($`𝒮_k`$ are given by Eqs (1a-d)), the mean-field expansion (71) includes the quadrupolar projection on $`p`$-states,

$`\rho _p^Q(\stackrel{}{n}_k)={\displaystyle \underset{I,J}{}}|I{\displaystyle \underset{a}{}}P(a)c_k(i^pj^p){\displaystyle \underset{\beta =1}{\overset{3}{}}}\delta _{i_\beta ^fj^a_\beta ^f}J|.`$
(74)
Here $`a`$ is a permutation of $`j_1^f,j_2^f,j_3^f`$ to $`j^a_\beta ^f`$, $`\beta =13`$; $`P(a)`$ is the parity of the permutation. The quadrupolar operator for the $`5f`$ electrons is given again by Eq. (55) where the index $`\beta `$ comprises the additional $`p`$electron, i.e. $`\beta =13`$.

Notice that in the ordered phase the ground state level is single and nonmagnetic, Table 16. This mechanism can explain the loss of magnetic moments because the magnetic excitations of the $`7p5f^3`$ configuration lie too high in energy.

### V.3 $`6d5f^3`$ configuration

The basis vectors here are

$`|I=|i_1^f,i_2^f,i_3^f;i^d,`$ (75)
where index $`i^f`$ stands for $`5f`$ states ($`i^f=114`$), while the index $`i^d=(k,s_z)`$ refers to five $`d`$-orbitals and the spin projection $`s_z`$. Thus, $`i^d=110`$, and in total there are 3640 nonequivalent basis vectors $`|I`$.

We start by considering the intrasite interactions $`H_{intra}`$. Here, in addition to $`ff`$ interactions we distinguish two groups. The first group arises between $`dd`$ and $`ff`$ transitions. It is described by the multipole Coulomb repulsion with the $`l=0`$, 2 and 4 angular components (SAF’s). The second group is due to the $`fd`$ and $`df`$ transitions. The corresponding multipole interactions are with $`l=1`$, 3 and 5. The relevant parameters were extracted from the LDA calculation of the Np ion in the $`6d5f^3`$ configuration,

$`v_1^{fddf}=11.322,v_3^{fddf}=3.701,v_5^{fddf}=1.794,`$
$`v_2^{ffdd}=11.289,v_4^{ffdd}=3.482\text{(in eV)},`$
$`\zeta _d=0.3497\text{eV}.`$ (76)
The parameters for the $`ff`$ interactions were kept unchanged. We then diagonalized the $`3640\times 3640`$ matrix of $`I|H_{intra}|J`$ and obtained 383 distinct levels. The 5 lowest and the highest levels are quoted in Table 17.

The CEF splittings of the lowest $`{}_{}{}^{5}L_{6}^{}`$ and $`{}_{}{}^{5}K_{5}^{}`$ levels are given in Table 18.

It should be noted that unlike before, the CEF operator acts not only on the $`5f`$ electrons but also on the $`6d`$ one, NM3

$`U^{CEF}(\stackrel{}{n})=B_4^f\rho _f^4(\stackrel{}{n})+B_6^f\rho _f^6(\stackrel{}{n})+B_4^d\rho _d^4(\stackrel{}{n}).`$ (77)
Here $`\rho _f^l(\stackrel{}{n})`$, $`l=4,6`$, and $`\rho _d^4(\stackrel{}{n})`$ are cubic projectors on the $`f`$ and $`d`$ states, respectively, given by

$`\rho _f^l(\stackrel{}{n})={\displaystyle \underset{I,J}{}}|I{\displaystyle }_aP(a){\displaystyle }_{\alpha =1}^3i_\alpha ^f|K_l|j^a{}_{\alpha }{}^{f}{\displaystyle }_{\beta =1}^3\delta _{i_\beta j_\beta ^a}J|,`$
(78a)
$`\rho _d^4(\stackrel{}{n})={\displaystyle \underset{I,J}{}}|I{\displaystyle }_aP(a)i^d|K_4|j^a{}_{}{}^{d}{\displaystyle }_{\beta =1}^3\delta _{i_\beta ^fj^a_\beta ^f}J|,`$
(78b)
where $`K_l(\mathrm{\Omega })`$ refers to the cubic harmonics with $`l=4`$ and 6. Here we keep the same notations as before \[Eqs (55) and (74)\], i.e. the permutation $`a`$ transforms the indices $`j_1^f,j_2^f,j_3^f`$ to $`j^a_\kappa ^f`$, $`\kappa =13`$. The permutations which interchange the $`d`$ and $`f`$ indices are excluded because they give zero contribution to (78a,b). The parameter $`B_4^d`$ was calculated by the method described in Sec. IV.B. For $`x_{eff}=0`$ we have found that

$`{\displaystyle \frac{q_4^d}{(R_{MT}^{Np})^4}}=0.1287,B_4^d=232.7\text{K}.`$ (79)
Notice that for the ground state the CEF gives a nonmagnetic single level, Table 18, but there are two low lying magnetic levels ($`T_2`$ and $`T_1`$) at 9.7 and 15.6 K, which contribute to the Curie law of the magnetic susceptibility at $`T>25`$ K.

Below 25 K the local symmetry of Np is lowered. The mean field is given by

$`U^{QQ}(\stackrel{}{n}_p)=\mathrm{\Lambda }^f\rho _f^Q(\stackrel{}{n}_p)\mathrm{\Lambda }^d\rho _d^Q(\stackrel{}{n}_p),`$ (80)
where

$`\mathrm{\Lambda }^f=\lambda ^{ff}\rho _f^Q+\lambda ^{fd}\rho _d^Q,`$ (81a)
$`\mathrm{\Lambda }^d=\lambda ^{df}\rho _f^Q+\lambda ^{dd}\rho _d^Q.`$ (81b)
Here again $`\rho _f^Q(\stackrel{}{n})`$ and $`\rho _d^Q(\stackrel{}{n})`$ are quadrupolar projection on the $`f`$ and $`d`$ states, respectively. They are given by expressions similar to (78a-b), where we replace $`K_l(\widehat{n}_p)`$ by $`𝒮_p`$, Eqs (1a-d), for 4 sublattices $`\{n_p\}=14`$ of $`Pa\overline{3}`$ or $`Pn\overline{3}m`$, Sec. II. \[Compare also with Eqs (55) and (74).\] Below we approximated the parameters of this interaction by

$`{\displaystyle \frac{q_2^d}{(R_{MT}^{Np})^2}}=0.1713,\mathrm{\Lambda }^f=6220\text{K},\mathrm{\Lambda }^d=3442\text{K}.`$
(82)
We then diagonalized the whole Hamiltonian $`H^{MF}`$ ($`H^{MF}=U^{QQ}(\stackrel{}{n})+V_{CF}+H_{intra}`$) and obtained the lowest energy levels quoted in Table 19.

Notice that now the first magnetic excitation of $`E`$ symmetry is separated from the nonmagnetic ground state by an energy gap of $`33`$ K, which implies again disappearance of the Curie law for the magnetic susceptibility of the ordered phase ($`T<25`$ K).

## VI Discussion and conclusions

In resonant X-ray scattering (RXS) experiments the transition to the ordered phase at $`T_c=25.5`$ K in NpO2 manifests itself by the appearance of superstructure Bragg reflections \[like (003) and others\], which are not compatible with the $`Fm\overline{3}m`$ structure of the disordered phase. Man ; Pai  Paixão et al., Ref. Pai, , ascribed the symmetry of the ordered phase to the $`Pn\overline{3}m`$ space group (No. 224, Ref. Tables, ), but did not elaborate on the theory. In this work (Sec. II), besides $`Pn\overline{3}m`$ we discuss also the space group $`Pa\overline{3}`$ (No. 205), which is the other possible candidate for the triple-$`\stackrel{}{q}^X`$ quadrupole ordering. $`Pn\overline{3}m`$ and $`Pa\overline{3}`$ are close symmetries. They imply condensations of different modes ($`X_4^+`$ and $`X_5^+`$, Ref. Sto, ) at the same $`X`$point of the Brillouin zone. We have considered the dependence of the scattering amplitude for different polarizations on the azimuthal angle $`\psi `$ and the Bragg angle $`\mathrm{\Theta }`$ taking into account the domain pattern of both symmetries. $`Pn\overline{3}m`$ and $`Pa\overline{3}`$ produce superstructure Bragg reflections at the same sites $`(h,k,l)`$ of the reciprocal lattice (Tables 1 and 2). The relation between the space groups is such that in RXS experiments the $`(00\mathrm{})`$ reflection of $`Pn\overline{3}m`$ behaves like the $`(\mathrm{}00)`$ and $`(0\mathrm{}0)`$ reflections of $`Pa\overline{3}`$ ($`\mathrm{}=2n+1`$). Therefore, we believe that a special care should be taken to distinguish between these two symmetries. We have shown that the two structures are completely different with respect to the direct bilinear quadrupole-quadrupole interactions. The $`Pn\overline{3}m`$ symmetry leads to a repulsion between the quadrupoles while the $`Pa\overline{3}`$ structure implies an attraction. (This is the main reason why the $`Pa\overline{3}`$ symmetry occurs in many molecular solids: NaO2 Zie , N2 Scott , C60. Hei2 ; Sach1 ; Dav ; Sach2 ; Mic1 )

We have presented an ab initio approach to crystal- and mean field for the structural phase transition at 25 K. The method is based on the technique of expanding the Coulomb repulsion between electrons in a multipolar series, Sec. III.

In the disordered phase we considered the Hamiltonian which includes the crystal electric field (CEF) effects and the intrasite Coulomb repulsion responsible for Hund’s rules on equal footing, Sec. IV. The crystal electric field levels above the ground state are in fact the lowest excitations of the neptunium electron complex. The typical splittings there are of the order of ten kelvins and thus the electron spectrum is very sensitive to the crystal symmetry. In this paper we have pursued the simple idea that the symmetry lowering produces a splitting of the many electron spectrum, which can explain the difference in the behavior of the magnetic susceptibility and the loss of the magnetic moments. Our band structure calculation (Sec. III.A) indicates that besides the three localized $`5f`$ electrons there is always approximately one conduction electron at each Np site. This changes the effective instantaneous configuration from the three electron $`5f^3`$ to a four electron one ($`7s5f^3`$, $`7p5f^3`$ or $`6d5f^3`$) and opens a possibility to obtain a nonmagnetic ground state without invoking the concept of the octupole order parameter.San2  From this point of view it represents an alternative to the latter, and we believe that both approaches deserves a thorough experimental consideration and verification.

Our finding that the CEF splitting in NpO2 is rather small ($`50150`$ K) contrasts the commonly accepted value of 55 meV. Amo  We think that the excitation at 55 meV observed experimentally by neutron scattering Amo  is a higher excitation of the neptunium electronic complex possibly involving the valence electrons on the Np-O bonds. The scale of CEF excitations calculated in the present work is a better match for the width of peaks in magnetic neutron scattering cross sections. Amo  We believe that the crystal and mean field should be of the same order of magnitude and the transition temperature of 25 K gives a natural estimation for it. An implicit support of our viewpoint is the fact that in the ordered phase there appears an inelastic peak centered at about 6.4 meV. Amo

The phase transition to the ordered phase sets in at 25 K and reduces the basic local symmetry of the neptunium sites to $`S_6=C_3\times i`$ (or $`\overline{3}`$). This symmetry holds for both $`Pn\overline{3}m`$ and $`Pa\overline{3}`$ space groups. \[The addition of three two-fold axes or three mirror planes for $`Pn\overline{3}m`$ changes the point group to $`D_{3d}`$ (or $`\overline{3}2`$) but this is not of principal importance.\] We have found that the structural change can not be ascribed to the bilinear quadrupole-quadrupole Coulomb repulsion, which is too weak to drive the transition. Therefore, one has to resort to the superexchange interactions super  via oxygen as a driving force for the transformation. In this work we have introduced an effective (enhanced) quadrupole interaction and studied the interplay between it and the crystal electric field. The influence of both interactions on the transition temperature has been investigated in detail for the $`5f^3`$ configuration, Sec. IV and figures 3, 4, and 5.

The most intriguing question is the disappearance of the neptunium magnetic moments below 25 K. This however can be explained if we consider a four electron complex at the neptunium site. We have demonstrated that all relevant configurations ($`7s5f^3`$, $`7p5f^3`$, $`6d5f^3`$) can lead to a nonmagnetic ground state separated from the magnetic excitations by an energy gap larger than $`25`$ K, Sec. V. Perhaps, the most clear example is the $`7p5f^3`$ configuration (Tables 15, 16). In the disordered phase the ground state is a triplet (Table 15), while in the ordered phase it becomes a singlet (Table 16). The general idea for the loss of magnetic moments is similar to the one used by Kondo and Anderson and often referred to as the Kondo effect. Ful  Notice however, that here we are dealing with the intrasite interactions treated on ab initio level. In particular, we replace the Anderson hybridization And2  which is linear in terms of creation/annihilation operators for valence and localized electrons, by the Coulomb intrasite repulsion, which being a density-density coupling is bilinear in terms of these operators. Another important theoretical ingredient of our model is the symmetry lowering which modifies the excitation spectrum of the electron system at low temperatures. This part is absent in the Kondo mechanism.

CEF and mean field have been objects of theoretical interest for many years Ste ; Hut ; New1 ; New2  and we would like to mention here some important relations between our model and other approaches. We have shown that CEF effects can be perceived as a first meaningful term of the intersite multipole expansion, when all neighbors of a neptunium site are considered in the spherical approximation ($`l^{}=0`$). It is then reduced to a single particle potential. New1 ; New2  The intersite nonspherical terms are also included in the full potential (FP) electron band structure calculations like FP-LMTO (linear muffin-tin orbital method) and FP-LAPW (linear augmented plane wave method). Sin  Therefore, in principle one can say that the CEF effects are equivalent to the full potential treatment. Wei ; ND  However, there are two very important caveats here. First, in the band structure calculations the nonspherical terms of the potential apply to itinerant electrons in the ground state, while CEF effects are considered usually for localized electrons in the ground and excited states. The second more important objection is that practically all band structure calculations are based on the single determinant approximation. This intrinsic feature does not allow to describe the intrasite interactions fully. In particular, the atomic term structure and Hund’s rules are excluded from the consideration. This shortcoming does not apply to our treatment which is based on a many determinant approach. For the intrasite part of interactions our model is very close to the scheme described by Condon and Shortley for the electron spectra of atoms and ions although there are some unimportant differences. The main approximation and limitation of our approach is the choice of basis. Once the many electron basis states are fixed, the work of calculating matrix elements is reduced mainly to classifications of all possible permutations, which is done numerically. Thus, the fundamental group of electron permutations Ell  is explicitly taken into account. No additional approximations used sometimes for crystal field calculations are made. This distinguishes our approach from Stevens’ Ste  where the CEF is expressed in terms of equivalent operators $`J_x`$, $`J_y`$, and $`J_z`$. The latter approach as well as the work of Lea, Leask and Wolf for cubic CEF Lea  based on it are justified only if $`J`$ is a good quantum number. Notice also that the approach of Stephens starts with the symmetry arguments while the interactions are introduced later in a phenomenological manner.

However, the present calculation scheme does not take into account chemical bonding in an intrinsic way. Therefore, further development of the method should be focused on this problem.

###### Acknowledgements.

We thank C. Detlefs, P. Santini and other authors of Ref. Pai,  for helpful discussions and communications on this problem. This work has been financially supported by the Fonds voor Wetenschappelijk Onderzoek, Vlaanderen.

## Appendix A Bilinear quadrupole interactions for $`Pn\overline{3}m`$ and $`Pa\overline{3}`$

Here we calculate the direct bilinear electronic quadrupole-quadrupole interactions for $`Pn\overline{3}m`$ and $`Pa\overline{3}`$ structures, Fig. 1. We consider the quadrupolar components $`S_i`$ of $`T_{2g}`$ symmetry ($`i=13`$) at site $`\stackrel{}{n}=0(0,0,0)`$ ($`\stackrel{}{n}\{n_1\}`$). There are 12 nearest neighbors of $`\stackrel{}{n}`$ belonging to three different sublattices $`\{n_p\}`$, $`=24`$. The interactions between three components $`S_i`$ centered at $`\stackrel{}{n}=0`$ and those ($`S_j`$) located at four nearest neighbors ($`\stackrel{}{n}^{}=14`$) of the second sublattice ($`\stackrel{}{n}^{}\{n_2\}`$) are given in Table 20.

The matrix $`S_i^{}(\stackrel{}{n})S_j^{}(\stackrel{}{n}^{})`$ for the fourth sublattice $`\{n_4\}`$ is given by the same Table provided that the cyclic permutation $`S_1S_3^{}`$, $`S_2S_1^{}`$, $`S_3S_2^{}`$ is performed. Here we label four nearest neighbors $`\stackrel{}{n}^{}=58`$ of $`\stackrel{}{n}`$ as $`5a(1/2,1/2,0)`$, $`6a(1/2,1/2,0)`$, $`7a(1/2,1/2,0)`$, and $`8a(1/2,1/2,0)`$. The matrix $`S_i^{\prime \prime }(\stackrel{}{n})S_j^{\prime \prime }(\stackrel{}{n}^{})`$ for the third sublattice $`\{n_3\}`$ is obtained from Table 20 by replacing $`S_1S_2^{\prime \prime }`$, $`S_2S_3^{\prime \prime }`$, $`S_3S_1^{\prime \prime }`$. Here $`\stackrel{}{n}^{}`$ runs over the sites $`9a(1/2,0,1/2)`$, $`10a(1/2,0,1/2)`$, $`11a(1/2,0,1/2)`$, and $`12a(1/2,0,1/2)`$.

For the $`Pn\overline{3}m`$ order of quadrupoles we have $`𝒮_a`$, Eq. (1a), at $`(0,0,0)`$ and $`𝒮_d`$, Eq. (1d), for $`\stackrel{}{n}^{}=14`$. By means of Table 20 we find that this interaction is given by

$`E_0^{QQ}(\{n_2\})={\displaystyle \frac{4}{3}}(\gamma 2\alpha ).`$ (83)
The same expression is obtained for the interactions with four neighbors of third and forth sublattices, $`E_0^{QQ}(\{n_3\})`$ and $`E_0^{QQ}(\{n_4\})`$. The total interaction then is

$`E^{QQ}(Pn\overline{3}m)`$ $`=`$ $`E_0^{QQ}(\{n_2\})+E_0^{QQ}(\{n_3\})+E_0^{QQ}(\{n_4\})`$ (84)
$`=`$ $`4(\gamma 2\alpha )>0,`$
because $`\gamma >0`$ and $`\alpha <0`$ independently of the lattice constant $`a`$.

For the $`Pa\overline{3}`$ structure we have the function $`𝒮_b`$, Eq. (1b), at $`\stackrel{}{n}^{}=14`$, and $`E_0^{QQ}(\{n_2\})=4\gamma /3`$. The same value is obtained for $`E_0^{QQ}(\{n_3\})`$ and $`E_0^{QQ}(\{n_4\})`$. As a result we arrive at

$`E^{QQ}(Pa\overline{3})=4\gamma <0.`$ (85)
We conclude that the $`Pn\overline{3}m`$ structure always leads to a repulsion while the $`Pa\overline{3}`$ to an attraction.

## Appendix B Correction of Slater integrals for Np

Since the experimental data on the energy splittings of the $`4f^3`$ configuration of Pr3+ and Nd4+ are available from the atomic database of NIST, Ref. La\_sp, , while there is no such information for $`5f^3`$, we have performed calculations of $`v_l^{FF}`$ ($`l=2,4,6`$), Eq. (33a), by using the radial dependence of $`_f`$ obtained from LDA calculations of atoms. After this we diagonalize the Hamiltonian of the free ion ($`V^{(3)}+H_{so}`$) and compared our calculated spectra with the experimental ones. We have noticed that the comparison is improved (the sequence of terms corresponds to the experimental one) if we reduce $`v_2^{FF}`$ and $`v_6^{FF}`$ by a factor of 0.75 while keep $`v_4^{FF}`$ almost the same (factor 0.975). Therefore, we have used the same scale factors for Np in NpO2 and obtained parameters given by Eq. (40).