text:
# X-ray Absorption Near-Edge Structure calculations with pseudopotentials. Application to 𝐾-edge in diamond and 𝛼-quartz

## 1 Introduction

XANES spectroscopy is a powerful technique to probe the empty states in solids. Moreover XANES gives information about the atomic arrangement around the probed atom up to the medium range order (i.e., about 8 Å). Nevertheless the interpretation of XANES spectra is not straightforward and often requires sophisticated simulation tools. To calculate a XANES spectrum, one needs to compute the absorption cross section. This is given by the Fermi golden rule as a sum of probabilities per unit of time of making a transition from an initial state to an unoccupied final state through an interaction Hamiltonian . Electronic transitions involved in X-ray absorption spectroscopy are mainly governed by the electric dipole operator. The main difficulty of any absorption cross section calculation lies in the solution of the Schrödinger equation for the final and initial states. The choice of the method used to solve the Schrödinger equation depends on the localized or delocalized character of final states. Crystal field multiplet theory  is adapted to the calculation of localized final states, i.e., with strong electron-electron interactions. It is typically the case of $`L_{2,3}`$ edges of transition elements  and $`M_{4,5}`$ edges of rare earth elements . If the final state is delocalized (the case of $`K`$ and $`L_1`$ edges) multielectronic interactions are weak and monoelectronic approaches based on the density functional theory (DFT) are usually employed. Among monoelectronic methods, one distinguishes the real space (cluster) approach from the reciprocal space (band structure) approach. Real space multiple scattering theory  has been extensively used in the past twenty years. However, multiple scattering theory suffers from an unavoidable drawback, i.e., the “muffin-tin” approximation of the potential. Recently, the “muffin-tin” approximation has been overtaken in a real space approach by using the finite difference method (MDF) to solve the Schrödinger equation . Although it yields good results , the method presented in Ref.  requires significant computing power. At present this is what restricts the applications to small clusters (around 50 atoms without any symmetry, at maximum).

With regards to band structure calculations, local projected densities of empty states (commonly called LDOS) have been often used to interpret XANES spectra (see for instance Refs. ). Recently several band structure codes have incorporated the radial matrix elements  to properly calculate XANES or ELNES (electron energy-loss near-edge structure). It is notably the case of `Wien2k` , which is based on the full-potential linearized augmented plane wave method , the orthogonalized-linear-combination-of-atomic-orbital program developed by Mo and Ching  and `CASTEP`, whose ELNES part has been developed by Pickard . This last approach uses pseudopotentials and reconstructs all electron wave functions within the projector augmented wave (PAW) method of Blöchl . Within these ab initio schemes the treatment of core hole-electron interactions is usually achieved through the inclusion of an excited atom into a supercell. The absorbing atom is then considered as an impurity. Beside these methods, Shirley  and more recently Soininen and Shirley  have proposed alternative formalisms that treat core-hole interactions in a different way. Indeed their methods incorporate core-hole interactions in the two-particle Bethe-Salpeter (BS) equation. Such methods are of great importance for optical absorption calculations . However, in the case of $`K`$-edge XANES spectra, the core-hole is frozen at one atomic site and the two-particle BS approach can be reduced to a single particle calculation.

By using supercell approaches, large systems (several hundreds of atoms) can be a priori treated. The charge density for such large cells is easy to calculate self-consistently by using pseudopotentials. On the other hand, the computation of the cross section is at present limited by the diagonalisation of the Hamiltonian for several empty states at many $`𝒌`$-points in the Brillouin zone. Indeed in Refs. pseudopotential calculations of ELNES spectra in diamond and $`c`$-BN have been performed with a small 16-atom supercell. Here we propose to use the recursion method of Haydock, Heine and Kelly  to make the calculation of the cross section computationally tractable for larger supercells. Thanks to this iterative technique, the calculation time of XANES becomes negligible compared with the calculation of the self-consistent charge density in large systems.

We incorporated the recursion method into a pseudopotential scheme. We first show in section 2 that the absorption cross section can be easily computed within a pseudopotential formalism using the PAW method. Then, by using the recursion method, we express the absorption cross section as a continued fraction. We present applications to $`K`$-edges in section 3. Two materials have been chosen : diamond which is an interesting example for method testing, and $`\alpha `$-quartz for which we have measured polarized XANES spectra.

## 2 Method

### 2.1 X-ray absorption cross section in the impurity model

In a monoelectronic approach, the X-ray absorption cross section $`\sigma (\omega )`$ can be written as

$$\sigma (\omega )=4\pi ^2\alpha _0\mathrm{}\omega \underset{f}{}|M_{if}|^2\delta (E_fE_i\mathrm{}\omega ),$$
(1)
where $`\mathrm{}\omega `$ is the photon energy, $`\alpha _0`$ is the fine structure constant and $`M_{if}`$ are the transition amplitudes between an initial core state $`|\psi _i`$ with energy $`E_i`$, localized on the absorbing atom site $`𝑹_0`$, and an all electron final state $`|\psi _f`$ with energy $`E_f`$

$$M_{if}=\psi _f|𝔒|\psi _i.$$
(2)
$`𝔒`$ is a transition operator coupling initial and final states. In the electric quadrupole approximation, $`𝔒`$ is given by $`\widehat{𝜺}𝒓(1+\frac{i}{2}𝒌𝒓)`$, where $`\widehat{𝜺}`$ and $`𝒌`$ are the polarisation vector and the wave vector of the photon beam, respectively. Within the frozen core approximation, $`|\psi _i`$ is a core state that can be taken from an all electron ground state atomic calculation. In the impurity model, $`|\psi _f`$ is an excited empty state that is solution of the Schrödinger equation for a potential that includes a core-hole on the absorbing atom.

### 2.2 PAW formalism

In the following we want to express the transition amplitude (Eq.2) within the PAW formalism, as originally described by Blöchl . We only recall the main aspects of the method that are needed to give a simple expression for the $`M_{if}`$ terms.

Within the PAW formalism, the final state all electron wave functions $`|\psi _f`$ are related to the corresponding final pseudo wave functions $`|\stackrel{~}{\psi }_f`$ through a linear operator $`𝒯`$:

$$|\psi _f=𝒯|\stackrel{~}{\psi }_f.$$
(3)
$`𝒯`$ differs from identity by a sum of local atom-centered contributions, that act only within spherical core regions centered on each atomic site $`𝑹`$, called augmentation regions or $`\mathrm{\Omega }_𝑹`$:

$$𝒯=\mathrm{𝟏}+\underset{𝑹,n}{}\left[|\varphi _{𝑹,n}|\stackrel{~}{\varphi }_{𝑹,n}\right]\stackrel{~}{p}_{𝑹,n}|.$$
(4)
Here $`|\varphi _{𝑹,n}`$ and $`|\stackrel{~}{\varphi }_{𝑹,n}`$ are the all electron and pseudo partial waves, respectively, which coincide outside $`\mathrm{\Omega }_𝑹`$. The vectors $`\stackrel{~}{p}_{𝑹,n}|`$, called projector functions , are equal to zero outside $`\mathrm{\Omega }_𝑹`$ and satisfy the condition $`\stackrel{~}{p}_{𝑹,n}|\stackrel{~}{\varphi }_{𝑹^{},n^{}}=\delta _{𝑹𝑹^{}}\delta _{nn^{}}`$. The index $`n`$ refers to the angular momentum quantum numbers $`(\mathrm{},m)`$ and to an additional number, used if there is more than one projector per angular momentum channel. The $`|\varphi _{𝑹,n}`$ form a complete basis for any physical non-core all electron wave function within $`\mathrm{\Omega }_𝑹`$. Therefore the $`|\stackrel{~}{\varphi }_{𝑹,n}`$ are also a complete basis for any physical pseudo wave function $`|\stackrel{~}{\psi }`$ within $`\mathrm{\Omega }_𝑹`$, i.e. for any $`𝒓|\chi _𝑹`$ function centered on an atomic site $`𝑹`$ and equal to zero outside $`\mathrm{\Omega }_𝑹`$,

$$\underset{n}{}\stackrel{~}{\psi }|\stackrel{~}{p}_{𝑹,n}\stackrel{~}{\varphi }_{𝑹,n}|\chi _𝑹=\stackrel{~}{\psi }|\chi _𝑹.$$
(5)
Substituting Eq. (4) in Eq. (3) and then Eq. (3) in Eq. (2), the transition amplitude $`M_{if}`$ becomes

$`M_{if}`$ $`=`$ $`\stackrel{~}{\psi }_f|𝔒|\psi _i`$ (6)
$`+{\displaystyle \underset{𝑹,n}{}}\stackrel{~}{\psi }_f|\stackrel{~}{p}_{𝑹,n}\varphi _{𝑹,n}|𝔒|\psi _i`$
$`{\displaystyle \underset{𝑹,n}{}}\stackrel{~}{\psi }_f|\stackrel{~}{p}_{𝑹,n}\stackrel{~}{\varphi }_{𝑹,n}|𝔒|\psi _i.`$
In Eq.6, the initial wave function $`𝒓|\psi _i`$ is localized on the site of the absorbing atom, $`𝑹_\mathrm{𝟎}`$, then only the $`𝑹_\mathrm{𝟎}`$ term has to be considered in each sum. Furthemore it should be noticed that $`𝒓|𝔒|\psi _i`$ is zero outside the $`\mathrm{\Omega }_{𝑹_0}`$ region. Therefore we can make use of Eq.(5) for the third term of Eq.(6), which thus vanishes with the first term. The transition amplitude $`M_{if}`$ is then reduced to one term. Now, introducing

$$|\stackrel{~}{\phi }_{𝑹_0}=\underset{n}{}|\stackrel{~}{p}_{𝑹_0,n}\varphi _{𝑹_0,n}|𝔒|\psi _i,$$
(7)
we obtain the following simple expression for the X-ray absorption cross section

$$\sigma (\omega )=4\pi ^2\alpha _0\mathrm{}\omega \underset{f}{}|\stackrel{~}{\psi }_f|\stackrel{~}{\phi }_{𝑹_0}|^2\delta (E_fE_i\mathrm{}\omega ).$$
(8)
The calculation of XANES spectra from Eq.8 presents us with the problem of determining many empty states $`|\stackrel{~}{\psi }_f`$, as has been mentioned in section 1. Indeed, the addition of unoccupied bands significantly increases computing time and then limits the size of the supercells. In the following, we show how the recursion method permits the cross section (Eq.8) to be rewritten as a continued fraction, so that only occupied bands have to be calculated.

### 2.3 Recursion method

The recursion method comprises of a powerful recursive algorithm that can be applied to a Hermitian matrix in order to transform it into a tridiagonal form. This method has been widely used in solid states physics. The pioneering works in this field were mainly by Haydock, Heine and Kelly . The first application to absorption cross section calculation was carried out by Filipponi within the multiple scattering formalism . In a close field Benedict and Shirley have implemented it for $`\epsilon _2(\omega )`$ calculation including core-hole interaction  within the BS approach. For simplicity, in the following, we assume that the norm of all electron partial waves coincide with that of the pseudo-partial waves.

In order to use the recursion method in our scheme, we shall introduce the Green operator into Eq.(8) by making the substitution

$$\underset{f}{}|\stackrel{~}{\psi }_f\delta (E_fE_i\mathrm{}\omega )\stackrel{~}{\psi }_f|=\frac{1}{\pi }𝔪[\stackrel{~}{G}(E)],$$
(9)
with

$$\stackrel{~}{G}(E)=(E\stackrel{~}{H}+i\gamma )^1.$$
(10)
In Eq.(10) $`\stackrel{~}{G}(E)`$ is the Green operator associated with the pseudo-Hamiltonian $`\stackrel{~}{H}=𝒯^{}H𝒯`$, which is Hermitian. The energy $`E`$ is given by $`E=E_i+\mathrm{}\omega `$ and $`\gamma `$ is an infinitesimal positive number. The cross section (Eq.8) can be rewritten as

$$\sigma (\omega )=4\pi \alpha _0\mathrm{}\omega 𝔪\left[\stackrel{~}{\phi }_{𝑹_\mathrm{𝟎}}|(E\stackrel{~}{H}+i\gamma )^1|\stackrel{~}{\phi }_{𝑹_\mathrm{𝟎}}\right].$$
(11)
Following the original work of Lanczos , the recursion method sets up a new basis in which the pseudo-Hamiltonian $`\stackrel{~}{H}`$, has a tridiagonal representation, from which the matrix elements $`\phi _{𝑹_\mathrm{𝟎}}|(\stackrel{~}{H}Ei\gamma )^1|\phi _{𝑹_\mathrm{𝟎}}`$ are very simply derived. This new basis (called Lanczos basis in the following) is obtained by the repeated action of $`\stackrel{~}{H}`$ onto the normalized initial vector $`|u_0=|\stackrel{~}{\phi }_{𝑹_\mathrm{𝟎}}/\sqrt{\stackrel{~}{\phi }_{𝑹_\mathrm{𝟎}}|\stackrel{~}{\phi }_{𝑹_\mathrm{𝟎}}}`$ through the symmetric three-term recurrence relation

$$\stackrel{~}{H}|u_i=a_i|u_i+b_{i+1}|u_{i+1}+b_i|u_{i1}.$$
(12)
where $`\{a_i\}`$ and $`\{b_i\}`$ are two sets of real parameters given by $`a_i=u_i|\stackrel{~}{H}|u_i`$ and $`b_i=u_i|\stackrel{~}{H}|u_{i1}=u_{i1}|\stackrel{~}{H}|u_i`$. The tridiagonal matrix representation of $`\stackrel{~}{H}`$ in the $`\{|u_i\}`$ basis then leads to the following form of the matrix elements of Eq. (11)

$$\stackrel{~}{\phi }_{𝑹_\mathrm{𝟎}}|(\stackrel{~}{H}Ei\gamma )^1|\stackrel{~}{\phi }_{𝑹_\mathrm{𝟎}}=\frac{\stackrel{~}{\phi }_{𝑹_\mathrm{𝟎}}|\stackrel{~}{\phi }_{𝑹_\mathrm{𝟎}}}{a_0Ei\gamma \frac{b_1^2}{a_1Ei\gamma \frac{b_2^2}{\mathrm{}}}}.$$
(13)
A simple terminator has been used to finish the continued fraction . In particular, if $`N`$ is the number of iterations required to converge the calculation, we consider that the coefficients $`(a_i,b_i)`$ are equal to $`(a_N,b_N)`$ for $`i>N`$, this leads to an analytical form of the terminator. It should be noted that the number of iterations $`N`$ strongly depends on the broadening parameter $`\gamma `$. With the iterative technique of Haydock the main part of XANES calculation involves the computation of the Hamiltonian acting on a single vector. This means that the computing time is considerably reduced compared to that can be required by Eq. 8 with an explicit diagonalisation.

### 2.4 Comparison with Bethe-Salpeter approach

Recently, the BS approach has been successfully applied to the calculation of XANES . Since the formalism of such two particles excited states calculations has been recently detailed in a review , here we only recall the main aspects. The BS procedure  consists of a ground state DFT calculation in order to obtain the Kohn-Sham orbitals, a $`GW`$ calculation to correct the eigenvalues , and a solution of BS equation using $`GW`$ eigenvalues and a statically screened potential.

Our method differs from the BS approach of Soininen and Shirley  by three points. The first two differences are substantial, the last one is just methodological. In particular, (i) we use the DFT eigenvalues, whereas Ref.  uses $`GW`$ corrected eigenvalues . If the $`GW`$ correction on the empty states can be described by a rigid shift, as it is often the case, this difference will not affect the XANES spectra. (ii) We screen the core-hole by the valence electron response at all orders computed self-consistently within DFT. In Ref. , the core-hole is screened at the linear order using the random phase approximation . (iii) We develop the final state wave functions $`|\stackrel{~}{\psi _f}`$ in a plane wave basis set of a supercell. In Ref. , the BS calculation uses as basis the eigenstates of the ground state Hamiltonian.

In the next section, our results will be compared with BS calculations as far as possible.

## 3 Applications

Our method was applied to the carbon $`K`$-edge in the diamond phase and to the silicon and oxygen $`K`$-edges in $`\alpha `$-quartz. XANES calculations were performed in two steps: (i) the evaluation of the self-consistent charge density for supercells including one $`1s`$ core-hole, (ii) the construction of a converged $`\{|u_i\}`$ basis leading to the tridiagonal representation of the pseudo-Hamiltonian. The transition matrix elements were determined in the electric dipole approximation. Calculations were carried out using an ab initio total energy code based on density functional theory within the local density approximation (LDA). Norm-conserving Troullier-Martins pseudopotentials  with a single component for each ($`\mathrm{},m`$) component were used. For the excited atom a pseudopotential with only one $`1s`$ electron was generated. In both diamond and $`\alpha `$-quartz cases, the size of the supercell was increased until the neighboring excited potentials did not interact with each other. The wave functions were expanded in plane waves with an energy cutoff of 50 Ry and 70 Ry for diamond and $`\alpha `$-quartz respectively. Reciprocal-space integrations were performed using a Monkhorst-Pack $`𝒌`$-point grid .

### 3.1 Diamond

Figure 1 shows calculations performed with four different sizes of supercell, all including a $`1s`$ core-hole on one of the carbon atoms. The size of the primitive cell – trigonal cell ($`a`$=2.524 Å) containing two carbon atoms – was successively multiplied in the three directions by a factor two (2$`\times `$2$`\times `$2 supercell, 16 atoms), three (3$`\times `$3$`\times `$3 supercell, 54 atoms), four (4$`\times `$4$`\times `$4 supercell, 128 atoms) and five (5$`\times `$5$`\times `$5 supercell, 250 atoms). Further increasing the supercell size does not provide significant modifications in the calculated spectrum. $`𝒌`$-point convergence of the charge density was obtained with a 2$`\times `$2$`\times `$2 $`𝒌`$-point grid for both 2$`\times `$2$`\times `$2 and 3$`\times `$3$`\times `$3 supercells, and with only one $`𝒌`$-point for larger supercells. $`𝒌`$-point convergence of the Lanczos basis was reached with a 10$`\times `$10$`\times `$10 $`𝒌`$-point grid for the 2$`\times `$2$`\times `$2 supercell, a 8$`\times `$8$`\times `$8 $`𝒌`$-point grid for the 3$`\times `$3$`\times `$3 supercell, and 6$`\times `$6$`\times `$6 $`𝒌`$-point grids for both 4$`\times `$4$`\times `$4 and 5$`\times `$5$`\times `$5 supercells. The number of iterations of the Lanczos basis was checked with $`\gamma =`$0.3 eV. The required number of $`|u_i`$ vectors varied from around 800 to 1400, while increasing the supercell size from 16 atoms to 250 atoms.

Calculations are compared with the carbon $`K`$-XANES spectrum of diamond measured by Ma et al. . The overall agreement between the experimental spectrum and the 5$`\times `$5$`\times `$5 supercell calculated one is quite good. In particular, the features from $`d`$ to $`i`$ are correctly reproduced as well in relative intensity as in energy position. Nevertheless in the first four eV of the spectrum the agreement between experiment and calculation remains unsatisfactory. The intensity of peak $`b`$ is exagerated and the exciton peak $`a`$ is not reproduced. Looking at the various supercell calculations, one notices that the first 25 eV of the spectrum are affected by the size of the cell, i.e. the interaction between core-holes belonging to neighboring cells. This gives information about the volume seen by the photoelectron during the absorption process: the photoelectron interacts with atoms located within a 6 Å radius sphere centered on the absorbing atom.

Several calculations of the carbon $`K`$-edge in diamond have been reported in the literature. Weng et al.  present pseudo-atomic-orbital density of $`p`$ empty states together with multiple scattering calculations. On the one hand, “muffin-tin” potentials as required in multiple scattering theory appear to be inappropriate in that case. This result was confirmed by a 300 atom cluster calculation performed by the FEFF8 code . On the other hand, although it neglects core-hole effects, LDOS calculation by Weng et al.  seems to correctly reproduce features $`h`$ and $`i`$. Pickard shows that including core-hole effects into a 2$`\times `$2$`\times `$2 supercell significantly improves the agreement with experiment. However, as the supercell size increases, peak $`b`$ becomes much too intense (figure 1). This may suggest an over estimation of the core-hole effects in that case. This result has also been observed by Soininen and Shirley while core-hole interactions is treated within BS approach. Their carbon $`K`$ edge calculation in diamond, is however very similar to our 128 and 250 atom supercell spectra. It seems that the corrections brought to excited states by $`GW`$ calculations has only a negligible influence on the energy positions of calculated features. This result is not surprising since previous studies have shown that the $`GW`$ approximation used for band structure calculation in diamond, amounts to applying a scissors-operator shift which is valid to within a 0.5 eV error. This error is comparable to the energy resolution of XANES features. The exageration of peak $`b`$ in both calculations requires more understanding and the problem of the treatment of core-hole interaction in diamond is at this time still unclear.

### 3.2 $`\alpha `$-quartz

In the case of $`\alpha `$-quartz, calculations are compared with polarized Si and O $`K`$-XANES experiments. $`\alpha `$-quartz ($`\alpha `$-SiO<sub>2</sub>) has a trigonal space group with a hexagonal unit cell. Since the point group is $`32`$ (or $`D_3`$), $`\alpha `$-quartz is a dichroic compound . Therefore any absorption spectrum can be expressed as a linear combination of the two spectra $`\sigma _{}`$ and $`\sigma _{}`$ corresponding to the polarisation vector $`𝜺`$ of the photon beam parallel and perpendicular to the ternary axis of the hexagonal cell (i.e. the $`c`$ axis), respectively. While calculations of isotropic XANES spectra of $`\alpha `$-quartz have been the subject of several studies in the past , little attention has been paid to the angular dependence of absorption. Lagarde et al.  show that the differences observed between $`\sigma _{}`$ and $`\sigma _{}`$ spectra are similar to that observed between the spectra of silica and densified silica. In Ref.  polarization effects are interpreted in terms of geometrical considerations. However, to our knowledge no convincing calculations of polarized data have been reported yet.

Si and O $`K`$-edge measurements were carried out respectively on the SA32 and SA72 beam lines at the Super-ACO facility of LURE (Orsay, France). The storage ring was operating at 800 MeV ($`\lambda _c=18.6`$ Å) and 200 mA electron current. The sample was a synthetic ($`10\overline{1}0`$) oriented single crystal of $`\alpha `$-quartz. The SA32 beam line was equipped with a double InSb (111) crystal monochromator. The instrumental resolution was estimated to be around 0.7 eV. The incident beam was monitored by an ionisation chamber while the absorption of the sample was measured in the electron drain current mode. Energy calibration was checked after each scan by measuring the $`c`$-Si spectrum for which the inflection point of the absorption edge was set to 1839 eV. The spectra were collected over a photon energy range of 1830-1910 eV with 0.2 eV steps and 1 s integration time. The sample was placed on a rotating holder, the rotation axis of which was parallel to the photon beam. The ($`10\overline{1}0`$) plane of the sample was set normal to the photon beam and the spectra were measured for the polarisation vector $`𝜺`$ either parallel to the $`c`$ axis ($`\sigma _{}`$) or perpendicular to the $`c`$ axis ($`\sigma _{}`$). We verified that the combination $`\frac{2}{3}\sigma _{}+\frac{1}{3}\sigma _{}`$ reproduced the powder sample well. The X-ray natural linear dichroism (XNLD) is obtained by the difference $`\sigma _{}\sigma _{}`$. The SA72 beam line was equipped with a toroidal grating monochromator. The estimated resolution was around 0.3 eV. Polarized O $`K`$-XANES were measured in the fluorescence mode over a photon energy range of 525-580 eV with 0.11 eV steps. The first 12 eV were slightly affected by self-absorption and the signal to noise was not good enough to permit an unambiguous extraction of the XNLD.

Calculations were performed with various sizes of cell built from the cristallographic parameters given in Ref. . Convergence in terms of supercell size was obtained for a 2$`\times `$2$`\times `$2 supercell (72 atoms) at both silicon and oxygen $`K`$-edges. In order to show the influence of core-hole effects we also present results obtained from the primitive cell (9 atoms) with or without core-hole (ground state calculation). Charge density was calculated with 2$`\times `$2$`\times `$2 $`𝒌`$-point grid for primitive cells (with or without core-hole) and only one $`𝒌`$-point for the 2$`\times `$2$`\times `$2 supercell. The Lanczos basis was determined with 4$`\times `$4$`\times `$4 $`𝒌`$-point grid for the primitive cell and 3$`\times `$3$`\times `$3 $`𝒌`$-point grid for the supercell. The broadening parameter $`\gamma `$, was set to 1 eV, this led to the construction of around 500 and 600 $`|u_i`$ vectors for the primitive cell and the supercell calculations respectively.

Experimental Si $`K`$-edge X-ray absorption spectra are compared with 2$`\times `$2$`\times `$2 supercell calculations in figure 2. From top to bottom one can see $`\sigma _{}`$ contribution (figure 2a), $`\sigma _{}`$ contribution (figure 2b), and the XNLD signal (figure 2c). Figure 2c shows that $`\alpha `$-quartz presents a strong angular dependence at the silicon $`K`$-edge, close to 9% of the white line’s intensity. Experimental and calculated spectra – polarized and dichroic– are in excellent agreement. The strength and the energy position of all the features are correctly reproduced. At the oxygen $`K`$-edge (figures 3a and 3b) the agreement with experiment is also satisfactory. $`\sigma _{}`$ and $`\sigma _{}`$ spectra differ essentially in the 12-19 eV energy range. These differences are well reproduced by the calculations.

Figure 4 compares three calculated $`\sigma _{}`$ spectra at both silicon and oxygen $`K`$-edges: the 2$`\times `$2$`\times `$2 supercell spectrum (solid line), the primitive cell calculation including core-hole effects (dashed line) and the ground state calculation (dotted line). At the silicon $`K`$-edge (upper panel), the ground state calculation does not reproduce the intense white line, which is clearly a signature of core hole-electron interaction. Including a $`1s`$ core-hole on one of the three silicon atoms of the unit cell drastically modifies the spectrum and results in the appearance of the white line. Increasing the size of the cell yields changes that are essentially in the 8-18 eV energy range, showing that the resonances occuring in that region are related to medium range organisation around the absorbing atom. At the oxygen $`K`$-edge (lower panel), although core-hole effects are less important than at the silicon $`K`$-edge, they definitely have to be taken into account in the calculation. Furthemore our simulations show the sensitivity to medium range order of the 12-20 eV region just above the main peak.

The importance of core-hole effects in $`\alpha `$-quartz was first noticed by Jollet and Noguera . This point was emphasized in multiple scattering studies . Our calculations confirm these previous considerations and demonstrate that reciprocal space calculations can also evidence the influence of medium range order, as is usually shown in cluster size analysis performed within the multiple scattering formalism .

Figures 2 and 3 display a good agreement between calculated and experimental spectra. Nevertheless the agreement in energy position may be discussed to a certain extent. Indeed calculated features appear slightly too contracted compared with experiment. For example, 1 eV shift is observed at 25 eV in figure 2. This contraction becomes more pronounced as energy increases. This small discrepancy could be due to the use of DFT eigenvalues. Unfortunately, in $`\alpha `$-quartz, the $`GW`$ corrections have been computed only for empty states that are close to the conduction band edge , thus we can not verify the accuracy of a scissor-operator approximation far from the edge.

Recently Mo and Ching  have calculated isotropic Si and O $`K`$-XANES and the corresponding LDOS using a supercell approach within an ab initio orthogonalized linear combination of atomic orbitals (OLCAO) method. A good agreement has been obtained between experiments and XANES calculations. With regards to LDOS calculations, the situation is quite different. They have shown that the Si $`K`$-edge can be reasonably represented by Si $`p`$ LDOS, while O $`p`$ LDOS surprisingly fails to reproduce any feature of the experimental spectrum. This result would mean that the radial dipole matrix element are strongly energy dependent in that case. This last consideration may be taken cautiously since the energy dependence of the dipole matrix element as well as the LDOS are widely dependent on the method used for their calculation. In particular, the LDOS depend on the size of the integration region used for the angular momentum projection. Only if this region is sufficiently small, the LDOS are independent of the size of the region, appart from a multiplicative factor. We calculated Si and O $`p`$ LDOS, using an integration sphere of 0.4 Å radius. The resulting LDOS are similar to the corresponding total cross sections computed with our approach, in contradiction with the conclusions of Mo and Ching . To explain this discrepancy, it should be noticed that Mo and Ching used the Mulliken population to compute the LDOS, which corresponds to a much larger effective integration region.

## 4 conclusion

We have presented an ab initio framework based on pseudopotentials for calculating XANES spectra. The scheme uses self-consistent full-potentials. It was applied to the calculation of $`K`$-edge spectra in diamond and $`\alpha `$-quartz. The results were compared with isotropic (diamond) and dichroic ($`\alpha `$-quartz) experimental data and a good agreement was obtained. In the case of diamond, our 250 atom supercell calculated spectrum is comparable with the BS calculated one of Soininen and Shirley. These first applications have shown that the introduction of the recursion method permits the treatment of large supercells. This opens new fields of applications like surfaces, small agregates or amorphous materials, for which large supercells are required. The continued fraction may also permit the investigation of more-extended x-ray absorption regions, without adding significant computation time. The results suggest that our scheme is quite robust and could be further successfully applied in XANES and ELNES spectroscopies.

We would like to thank G. Tourillon and C. Laffon for having provided polarized O $`K`$-edge experimental spectra. We also thank J.J. Rehr for encouraging and fruitful discussions. The english was polished by S. Joyce. We are grateful to Ch. Brouder for his thorough reading of the manuscript. This is an IPGP contribution #XXXX.


format: text

text:
# 1 Introduction

## 1 Introduction

In its weak–field and slow–motion approximation General Relativity predicts that, among other things, the orbit of a test particle freely falling in the gravitational field of a central rotating body is affected by the so called gravitomagnetic dragging of the inertial frames or Lense–Thirring effect. More precisely, the longitude of the ascending node $`\mathrm{\Omega }`$ and the argument of the perigee $`\omega `$ of the orbit undergo tiny precessions according to \[Lense and Thirring, 1918\]

$`\dot{\mathrm{\Omega }}_{\mathrm{LT}}`$ $`=`$ $`{\displaystyle \frac{2GJ}{c^2a^3(1e^2)^{\frac{3}{2}}}},`$ ($`1`$)
$`\dot{\omega }_{\mathrm{LT}}`$ $`=`$ $`{\displaystyle \frac{6GJ\mathrm{cos}i}{c^2a^3(1e^2)^{\frac{3}{2}}}},`$ ($`2`$)
in which $`G`$ is the Newtonian gravitational constant, $`J`$ is the proper angular momentum of the central body, $`c`$ is the speed of light $`invacuum`$, $`a,e`$ and $`i`$ are the semimajor axis, the eccentricity and the inclination, respectively, of the orbit of the test particle. The Lense–Thirring precessions for the LAGEOS satellites amount to

$`\dot{\mathrm{\Omega }}_{\mathrm{LT}}^{\mathrm{LAGEOS}}`$ $`=`$ $`31\text{mas/y},`$ ($`3`$)
$`\dot{\mathrm{\Omega }}_{\mathrm{LT}}^{\mathrm{LAGEOS}\mathrm{II}}`$ $`=`$ $`31.5\text{mas/y},`$ ($`4`$)
$`\dot{\omega }_{\mathrm{LT}}^{\mathrm{LAGEOS}}`$ $`=`$ $`31.6\text{mas/y},`$ ($`5`$)
$`\dot{\omega }_{\mathrm{LT}}^{\mathrm{LAGEOS}\mathrm{II}}`$ $`=`$ $`57\text{mas/y}.`$ ($`6`$)
The first measurement of this effect in the gravitational field of the Earth has been obtained by analyzing a suitable combination of the laser-ranged data to the existing passive geodetic satellites LAGEOS and LAGEOS II \[Ciufolini et al., 1998\]. The observable \[Ciufolini, 1996\] is a linear trend with a slope of 60.2 milliarcseconds per year (mas/y in the following) and includes the residuals of the nodes of LAGEOS and LAGEOS II and the perigee of LAGEOS II111The perigee of LAGEOS was not used because it introduces large observational errors due to the smallness of the LAGEOS eccentricity \[Ciufolini, 1996\] which amounts to 0.0045.. The total relative accuracy of the measurement of the solve-for parameter $`\mu _{\mathrm{LT}}`$, introduced in order to account for this general relativistic effect, is of the order of $`2\times 10^1`$ \[Ciufolini et al., 1998\].

In this kind of satellite–based space experiments the major source of systematic errors is represented by the aliasing trends due to the classical secular precessions \[Kaula, 1966\] of the node and the perigee induced by the mismodelled even zonal harmonics of the geopotential $`J_2,J_4,J_6,\mathrm{}`$ Indeed, according to the present knowledge of the Earth’s gravity field based on the EGM96 model \[Lemoine et al., 1998\], they amount to a large part of the gravitomagnetic precessions of interest, especially for the first two even zonal harmonics. In the performed LAGEOS–LAGEOS II Lense–Thirring experiment the adopted observable allows for the cancellation of the static and dynamical effects of $`J_2`$ and $`J_4`$. The remaining higher degree even zonal harmonics affects the measurement at a $`13\%`$ level.

In order to achieve a few percent accuracy, in \[Ciufolini, 1986\] it was proposed to launch a passive geodetic laser-ranged satellite- the former LAGEOS III - with the same orbital parameters of LAGEOS apart from its inclination which should be supplementary to that of LAGEOS.

This orbital configuration would be able to cancel out exactly the classical nodal precessions, which are proportional to $`\mathrm{cos}i`$, provided that the observable to be adopted is the sum of the residuals of the nodal precessions of LAGEOS III and LAGEOS

$$\delta \dot{\mathrm{\Omega }}^{\mathrm{III}}+\delta \dot{\mathrm{\Omega }}^\mathrm{I}=62\mu _{\mathrm{LT}}.$$
($`7`$)
Later on the concept of the mission slightly changed. The area-to-mass ratio of LAGEOS III was reduced in order to make less relevant the impact of the non-gravitational perturbations and the eccentricity was enhanced in order to be able to perform other general relativistic tests: the LARES was born \[Ciufolini, 1998\]. The orbital parameters of LAGEOS, LAGEOS II and LARES are in Table 1.

Recent developments of the concept of the twin satellites in supplementary orbits have led to the discover of new, possible gravitomagnetic observables based on the use of the perigees as well \[Iorio, 2002a\] and of unexpected connections with the gravitomagnetic clock effect \[Iorio and Lichtenegger, 2002\].

Unfortunately, at present we do not know if the LARES mission will be approved by any space agency. Although much cheaper than other proposed and/or approved space–based missions, funding is the major obstacle in implementing the LARES project. The most expensive part is the launching segment.

Very recently the possibility of launching the LARES satellite into an orbit with $`a=8,378`$ km, $`i=90`$ deg, $`e=0.04`$ and using as observable its node has been considered \[Lucchesi and Paolozzi, 2001\]. In the following we will name POLARES the LARES satellite in such proposed polar orbit. The Lense–Thirring secular rate of the POLARES node would amount to 96.9 mas/y. Then, the observable would be

$$\delta \dot{\mathrm{\Omega }}^{\mathrm{PL}}=96.9\mu _{\mathrm{LT}},$$
($`8`$)
i.e. a linear secular trend with a slope of 96.9 mas/y. The choice of a so low altitude is motivated by the need of using a cheap rocket launcher; so, it must certainly be considered as an admirable and important further effort towards the practical realization of the LARES project. The polar orbit would allow to prevent the aliasing effects of the mismodelled classical secular precessions of the node induced by the even $`l=2n`$ zonal $`m=0`$ coefficients of the multipolar expansion of the static part of the geopotential. The non–gravitational perturbations \[Lucchesi, 2001; 2002; Lucchesi and Paolozzi,  2001\] would represent a minor problem.

In this paper we wish to critically analyze such important and interesting evolution of the LARES concept.

## 2 The POLARES

### 2.1 The node–only scenario

Here we will show that the proposed POLARES configuration and the use of the node only as gravitomagnetic observable would present some drawbacks due to the gravitational static and time–varying perturbations.

In regard to the classical secular precessions induced by the centrifugal oblateness of the Earth, it should be considered that the choice of the inclination would be of crucial importance. Such orbital element is affected neither by the secular perturbations induced by the even zonal harmonics of the geopotential \[Kaula, 1966\] nor by the semisecular 18.6–year and 9.3–year tidal perturbations \[Iorio, 2001\]. On the other hand, due to possible orbital injection errors which are closely related to the quality of the rocket launcher to be used, the POLARES inclination would be far from being exactly 90 deg. The relatively low altitude of POLARES would enhance the impact of even small departures of $`i_{\mathrm{PL}}`$ from its projected nominal value on the systematic error due to the classical static even zonal harmonics of geopotential $`\delta \mu _{\mathrm{LT}}^{\mathrm{even}\mathrm{zonals}}`$. This is clearly shown in Figure 1 in which it has been calculated by means of the covariance matrix of the EGM96 Earth gravity model up to degree $`l=20`$.

It should be considered that, contrary to the LAGEOS satellites which are almost insensitive to the even zonal harmonics of degree higher than $`l=20`$, this would not be the case for the POLARES, due to its projected altitude of 2,000 km. Then, it might turn out that the estimates of Figure 1 are optimistic.

If we consider that the impact of the mismodelled even zonal harmonics of the geopotential on the current LAGEOS–LAGEOS II Lense–Thirring experiment is of the order of 13$`\%`$, according to EGM96, a certain weakness of the POLARES node–only scenario becomes apparent. In regard to the possibilities offered by the forthcoming new gravity models from the CHAMP and GRACE missions, it should be considered that the major improvements in the accuracy of the even zonal coefficients of the geopotential are expected for the first degrees. Now, in the LAGEOS–LAGEOS II Lense–Thirring experiment the first two even zonal harmonics are cancelled out and the terms of degree higher than $`l=20`$ are not relevant because the LAGEOS and LAGEOS II orbits are not sensitive to them. So, the systematic error due to the geopotential, in this case, is due to the even zonal harmonics from $`l=6`$ to $`l=20`$. Then, the CHAMP and GRACE missions would certainly improve such error. On the other hand, the POLARES configuration would be sensitive both to the first two even zonal harmonics, which should be greatly improved by the CHAMP and GRACE results, and to the terms of degree higher than $`l=20`$, for which the improvements should be less relevant. So, even in this perspective, maybe that a simple reanalysis of the LAGEOS and LAGEOS II data should be more fruitful than the node–only POLARES choice.

Such conclusion is enforced also by a simple evaluation of the impact of some tidal perturbations222The amplitude of the tidal perturbations on the node is proportional to $`(1/\mathrm{sin}i)dF_{lmp}/di`$ where $`F_{lmp}(i)`$ are the inclination functions \[Kaula, 1966\]. For the even zonal ($`m=0`$), tesseral ($`m=1`$) and sectorial ($`m=2`$) terms they are $`F_{201}=(3/4)\mathrm{sin}^2i(1/2)`$, $`F_{211}=(3/2)\mathrm{sin}i\mathrm{cos}i`$ and $`F_{221}=(3/2)\mathrm{sin}^2i`$, respectively. \[Iorio, 2001\]. Putting aside the fact that the 18.6-year lunar tide, which has a period of 18.6 years due to the lunar node only motion and a large amplitude, is a $`l=2,m=0`$ constituent which would affect the POLARES node for $`i_{\mathrm{PL}}90`$ deg, let us draw our attention to the $`l=2,m=1`$ $`K_1`$ tesseral tide. It is an important constituent which exerts a relevant perturbing action on the node of a satellite333For LAGEOS it has a nominal amplitude of the order of 103 mas \[Iorio, 2001\] and has the same period of just the node of the satellite. It turns out that for, say, $`i_{\mathrm{PL}}=89.8`$ deg the period of the POLARES node would amount to 26,801 days, i.e. 73.3 years. This means that, even by not considering at all the even zonal perturbations, the mismodelled part of the $`K_1`$ tide itself would resemble a superimposed aliasing trend over a reasonable observational time span of a few years and would completely mask the Lense–Thirring trend.

### 2.2 The combined residuals scenario

Let us try to see what could happen by inserting the node and the perigee of POLARES in a suitable combination of orbital residuals with the nodes of LAGEOS and LAGEOS II and the perigee of LAGEOS II along the lines sketched in \[Ciufolini, 1996; Iorio, 2002b\]. Recall that if $`N`$ orbital elements are present in such combinations, the effects of the first $`N1`$ even zonal harmonics of the geopotential are cancelled out, irrespectively of the orbital geometry of the employed satellites.

It turns out that by combining the nodes of LAGEOS, LAGEOS II and POLARES and the perigees of LAGEOS II and POLARES in

$$\delta \dot{\mathrm{\Omega }}^\mathrm{I}+c_1\delta \dot{\mathrm{\Omega }}^{\mathrm{II}}+c_2\delta \dot{\mathrm{\Omega }}^{\mathrm{PL}}+c_3\delta \dot{\omega }^{\mathrm{II}}+c_4\delta \dot{\omega }^{\mathrm{PL}}=X_{\mathrm{LT}}\mu _{\mathrm{LT}}$$
($`9`$)
the inclination $`i_{\mathrm{PL}}=90`$ deg is singular in the sense that the slope of the relativistic effect diverges because the coefficient $`c_2`$ with which the node of POLARES enters the combination diverges. For values of $`i_{\mathrm{PL}}`$ close to 90 deg Figure 2 and Figure 3 show that the systematic error due to the remaining $`\delta J_{10},\delta J_{12},\mathrm{}`$, according to the covariance matrix of EGM96 up to degree $`l=20`$, would seriously affect the measurement of the Lense–Thirring linear trend. Moreover, the impact of all the gravitational and non–gravitational time–varying perturbations would be greatly enhanced by the large values of $`c_2`$.

The situation would not be better even if we would include only the perigee of POLARES in the combined residuals

$$\delta \dot{\mathrm{\Omega }}^\mathrm{I}+c_1\delta \dot{\mathrm{\Omega }}^{\mathrm{II}}+c_2\delta \dot{\omega }^{\mathrm{II}}+c_3\delta \dot{\omega }^{\mathrm{PL}}=X_{\mathrm{LT}}\mu _{\mathrm{LT}},$$
($`10`$)
as shown by Figure 4.

## 3 Conclusions

In this paper the proposal of putting the LARES satellite into a polar, elliptical orbit with an altitude of 2,000 km in order to look at its gravitomagnetic secular node shift has been critically analyzed.

The key point is that the mismodelled classical nodal secular precessions induced by the even zonal coefficients of the multipolar expansion of the terrestrial gravitational field vanished if and only if the inclination of the satellite would be exactly $`90`$ deg. Of course, mainly due to possible orbital injection errors, this could never happen. It turns out that the low altitude of the proposed orbital configuration, and the consequent high sensitivity to the higher even degree zonal terms of the geopotential, would greatly enhance the impact of even small departures of the real values of the inclination from the nominal value of 90 deg. For example, for just $`i_{\mathrm{PL}}=90\pm 0.2`$ deg the systematic gravitational error would be of the order of 5$`\%`$–10$`\%`$, according to the EGM96 Earth gravity model up to $`l=20`$, which is of the same order of magnitude of the present LAGEOS–LAGEOS II experiment (Its total error, including various systematic gravitational and non–gravitational perturbations, is of the order of 20$`\%`$–30$`\%`$). Of course, such a situation would be further made critical by the possible use of a low–cost launcher which, unavoidably, would induce not negligible orbital injection errors. Moreover, the tesseral $`K_1`$ tidal perturbation, which has the same period of the satellite’s node, would induce a secular aliasing trend over an observational time span of a few years because its period would amount to several tens of years for near polar orbits.

Another important point is that the proposed POLARES would not yield substantial improvements also in the context of the combined residuals approach which allows to cancel out the contribution of the first even zonal coefficients of the geopotential irrespectively of the inclination of the satellites. Indeed, it turns out that a combination including the nodes of LAGEOS, LAGEOS II, POLARES and the perigees of LAGEOS II and POLARES would be not defined for $`i_{\mathrm{PL}}=90`$ deg because the coefficient weighing the node of the LARES would go to infinity. For very small deviations of $`i_{\mathrm{PL}}`$ from such a critical value the systematic error induced by the remaining even zonal harmonics would amount to 10$`\%`$–20$`\%`$. Moreover, the coefficient with which the node of POLARES would enter the combination would be much more larger than unity and would greatly enhance the impact of the gravitational and non–gravitational time–dependent perturbations. Even the use of the nodes of LAGEOS and LAGEOS II and the perigees of LAGEOS II and POLARES, in which case $`i_{\mathrm{PL}}=90`$ deg would not create problems, would not yield benefits because the systematic gravitational error would be of the order of 20$`\%`$–30$`\%`$.

Moreover, the accuracy of the practical data reduction from such version of LARES would be affected by the atmospheric drag.

Finally, we can conclude that the proposed low–cost version of the LARES mission would not yield any significant improvements in the measurement of the elusive Lense–Thirring effect, according to the present–day level of knowledge of the Earth’s gravitational field. Perhaps, the situation could improve to a certain extent with the new, more accurate gravity models from CHAMP and GRACE missions which should become available in the next few years. On the contrary, the original concept of the couple of supplementary satellites would deserve grater attention thanks to its much richer spectrum of high accuracy relativistic observables.


text:
# Interesting thermomagnetic history effects in the antiferromagnetic state of SmMn2Ge2

## Abstract

We present results of magnetization measurements showing that the magnetic response of the antiferromagnetic state of SmMn2Ge2 depends on the path used in the field(H)-temperature(T) phase space to reach this state. Distinct signature of metastablity is observed in this antiferromagnetic state when obtained via field-cooling/field-warming paths. The isothermal M-H loops show lack of end-point memory, reminiscent of that seen in metastable vortex states near the field-induced first order phase transition in various type-II superconductors.

The intermetallic compound SmMn2Ge2 with its interesting magnetic properties has been a subject of intensive study during last two decades. In low applied magnetic fields it shows at least three magnetic transitions as a function of temperature. First it undergoes a paramagnetic (PM) to ferromagnetic (FM1) transition at around 350K, followed by an FM1 to antiferromagnetic (AFM) transition at around 160 K (TN1). On reduction of the temperature further this AFM state transforms again into another ferromagnetic (FM2) state around 100 K(TN2). There is a large spread in the reported values of the transition temperatures from FM1 to AFM and AFM to FM2 states. Quality of the samples may be a possible source for the reported differences in the transition temperatures, especially when it is known that the microscopic magnetic properties of RE(rare-earth)1-2-2 compomuds are quite sensitive to their underlying crystal lattice structure. On the other hand, there exist now enough evidences from various studies that both of these transitions are probably first order in nature . The first order nature of these magnetic phase transitions can also provide a natural explanation to the reported spread in the transition temperatures. Supercooling(superheating) can take place down(up) to a temperature T (T∗∗) while cooling (heating) across a first order transition point (TN) . The extent of supercooled/superheated phases will depend on the path followed in the field(H)-temperature(T) phase space . In addition in the samples with defect structures the lower(higher) temperature phase will start nucleating around these defect structures once the sample is cooled(heated) across TN. This nucleation of the lower(higher) T phase will be completed at T(T∗∗), and in the temperature regime TN-T (T∗∗-TN) there will be co-existence of two phases. All these properties will give rise to thermal hysteresis, and such thermal hysteresis is actually observed across the FM2-AFM and AFM-FM1 phase transitions in SmMn2Ge2 .

Confined between two FM phase at low and high temperature and reached via first order phase transitions, the AFM phase in SmMn2Ge2 is something special. In this paper, based on our careful dc magnetization studies we shall show that the magnetic response of this AFM state actually depends on the path used in the (H,T) phase space to reach this state. Distinct trace of the high(low) temperature FM1(FM2) state remains well inside the AFM state, when this state is reached via a field cooling(warming) path. We seek a possible explanation of the observed behaviour in terms of supercooling/superheating and phase coexistence across a first order phase transition.

The SmMn2Ge2 sample used in our present study was prepared by argon arc melting and characterized by X-ray diffraction(XRD) measurements . Dc magnetization measurements were performed with a commercial SQUID-magnetometer (MPMS-5, Quantum Design) using a 4 cm scan length.

A low field ( H=50 Oe) magnetization(M) versus temperature(T) measurement reveal the FM1-AFM and AFM-FM2 transition temperatures in the present SmMn2Ge2 sample to be approximately 150 K and 105 K respectively (see Fig.1). We thus choose a temperature T=120 K which is well inside the AFM state away from both the upper and lower temperature phase boundaries. Moreover, the M-H measurements at 120K reveal a field-induced ferromagnetic transition around H$``$4 kOe (see the inset of Fig.1). To keep the sample well inside the AFM state in the (H,T) space, the applied field in the present measurements is limited to a maximum value of Hmax=1 kOe . We present in Fig.2(a) the M-H plot at T=120 K, measured after reaching this temperature in zero field cooled (ZFC) condition. There is a small non-linearity as well as a small but distinct hysteresis in this M-H curve, which indicates that the magnetic state is not pure antiferromagnetic in nature and there probably exists a finite amount of spin-canting. The same feature may also arise from very small amount (1-2%) of ferromagnetic phase which can go undetected in standard XRD measurements. In fact similar non-linearity in the M-H curves in the antiferromagnetic state of various CeFe2 based pseudobinary alloys, has earlier been attributed to small amount of ferromagnetic impurity phase which could not be detected in the XRD studies. However, subsequent studies have indicated that this feature can be of intrinsic origin. In our present system also various studies (to be narrated below) indicate that the unconventional properties of the antiferromagnetic state in SmMn2Ge2 cannot be explained in terms of a small amount of ferromagnetic impurity phase.

We shall now study the magnetic response of the AFM state at 120K after reaching this state following the three different experimental protocols :

1. the temperature 120 K is reached from temperature well above TN1 in absence of any applied field, and then a field of 1 kOe is switched on.
2. a field of 1 kOe is switched on within the FM2 state at 4.5K and the AFM state is reached subsequently by warming up the sample unidirectionally across TN2 to 120K.
3. a field of 1 kOe is switched on within the FM1 state at 300 K and the AFM state is reached subsequently by cooling the sample unidirectionally across TN1 to 120K.

In all these experimental protocols sufficient wait time is given after reaching 120K to ensure complete temperature stability. Furthermore, in experimental protocol no. 1 and 3 temperature is reduced slowly in small steps to avoid any temperature oscillation.

The values of magnetization measured at 120 K with H=1 kOe in both the protocols 2 and 3 are higher than that measured after switching on the field of 1 kOe at 120K in the ZFC condition i.e.with the protocol 1. This observation can be rationalized in terms of supercooling (superheating) of the FM1(FM2) state. While cooling(warming) acrosss the FM1(FM2)-AFM transition temperatures TN1(TN2) some amount of the FM1(FM2) state will supercool (superheat) into the temperature regime well beyond the transition temperature. The extent of temperature regime of supercooling/superheating actually widens in presence of the applied magnetic field .

To support the above conjecture we have studied the detailed nature of the magnetic state at T=120K obtained by three different experimental protocols mentioned above. We use a ”$`\mathrm{𝑚𝑖𝑛𝑜𝑟}`$ $`\mathrm{ℎ𝑦𝑠𝑡𝑒𝑟𝑒𝑠𝑖𝑠}`$ $`\mathrm{𝑙𝑜𝑜𝑝}`$ $`\mathrm{𝑡𝑒𝑐ℎ𝑛𝑖𝑞𝑢𝑒}`$” to show that the AFM state at 120K obtained in ZFC condition is a stable magnetic state, while the magnetic states obtained at 120K by cooling/heating in presence of 1 kOe field are metastable in nature. Taking the M-H curve drawn in the ZFC AFM state at 120K by field cycling between 1 kOe and 0 kOe as the envelope curve, we draw minor hysteresis loops (MHLs), terminating the M-H curve at various field points ($`0<H<1kOe`$) on the the descending field leg of the envelope curve (see Fig. 2(a) and 3(a)). A distinct ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$” is observed for all the MHls, namely on completion of the cycle the MHLs show the same end point magnetization value on the envelope M-H curve at 1 kOe. This kind of ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$” is common with various kinds of hysteretic systems including hard ferromagnets and superconductors . To emphasise this point we show in Fig. 2 (b) and 2 (c) similar M-H curves obtained by field cycling between 1 k Oe and 0 Oe at T=20K and 220 K which are well inside the FM2 and FM1 phase respectively . Clear ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$” is observed in both the cases.

We shall now draw the envelope curve and the MHLs in the AFM state at 120K obtained by heating from 4.5 K in the presence of a field of 1 kOe i.e. following the protocol no.2. In contrast to the ZFC case, the envelope curve and the MHLs are of very different nature (see Fig.3(b)).The most prominent difference is the absence of ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$”. Starting from the H=1 kOe point and returning back to this point by drawing MHLs of increasing field amplitude, a steadily decreasing value of magnetization is observed at the end point i.e. H=1 kOe. The end point memory is in fact lost in the process of drawing the first MHL by cycling between 1 kOe and 0.8 kOe (see the paths marked by 1 and 2 in the inset of fig. 3(b)). The envelope curve drawn by lowering the field from 1 kOe after this first cycle (see path marked by 3 in the inset of Fig.3) is quite different from the initial envelope curve. The differences in both the end point magnetization and the envelope curves rise steadily with further cycling of field with succesively larger field amplitude (see path marked by 5,6,7,8,9 and 10 in the inset of Fig.3(b)). Same field cycling process does not cause any effect on the end point magnetization and envelope curve in the magnetic state obtained with the experimental protocol no.1 (see Fig. 3(a)). This clearly shows that the initial magnetic state obtained with protocol no. 2 is a metastable state, and the energy fluctuation introduced while drawing the MHLs steadily push this state towards the stable ZFC state. This kind of the lack of ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$” has earlier been observed across the vortex matter phase transition from one kind of vortex solid to another in a type-II superconductor CeRu2 (see Fig. 7 of Ref.22 ). This was attributed to the existence of metastable states across a disorder broadened first order transition . These metastable states were shattered through the energy fluctuation generated while drawing the MHLs in the concerned (H,T) regime. The same behaviour has been observed subsequently across a disorder induced first order transition in the vortex matter phase space of another type-II superconductor NbSe2 .

The same lack of ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$” effect is observed at 120 K on cooling form 300 K well inside the FM1 state i.e. following the protocol no.3 (see Fig. 3(c)). On the other hand same kind of experiments after preparing the magnetic states well inside the ferromagnetic regions FM1 and FM2 by crossing the transition temperatures TN1 and TN2 in presence of field, show clear ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$” effect (data not shown here for the sake of clarity and conciseness). These results also negate any possible contribution from an impurity ferromagnetic phase in the observed metastable behaviour of the AFM state. Metastability (if any) related to the hindrance of domain motions in a ferromagnet is known to be more while warming up from the low temperature ZFC state than while cooling down from the high temperature region in the presence of an applied field.

This lack of ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$” is now accepted as a signature of metastability associated with a first order transition , which is taken as a support for establishing the first order nature of certain vortex matter phase transitions in type-II superconductors . In this paper we have looked at the AFM state of SmMn2Ge2 which can be reached from both high and low temperature FM state through magnetic phase transitions. The first order nature of these two transitions is already considered . We have now shown the lack of ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$” in the minor hysteresis loops and associated metastable behaviour in this AFM state of SmMn2Ge2. The observed behaviour highlights the interesting status of the AFM state in SmMn2Ge2 sandwiched between two FM states in the (H,T) phase space and reached via first order phase transitions.

Summarising our results, we find interesting thermomagnetic history effects well inside the AFM state of SmMn2Ge2. When this AFM state is reached either from the FM1(FM2) state by cooling (heating) in presence of an applied field, the traces of the FM1(FM2) state remain as supercooled(superheated) state. This is an example of phase coexistence of the AFM-FM1(FM2) state, and the resulting magnetic state is metastable in nature. On drawing MHLs in this metastable state one introduces energy fluctuations, which drive the domains of metastable FM1(FM2) state to the stable AFM state. Such metastability in the form of the lack of ”$`\mathrm{𝑒𝑛𝑑}`$ $`\mathrm{𝑝𝑜𝑖𝑛𝑡}`$ $`\mathrm{𝑚𝑒𝑚𝑜𝑟𝑦}`$” (which is also observed across solid to solid vortex matter phase transitions in various type-II superconductors (ref.22 and 23)) may serve as characteristic signatures of a first order transition in samples with substantial defect structures where the detection of latent heat as the canonical signature of a first order transition is relatively difficult .


text:
# COMPARISON OF DIFFERENT GOODNESS-OF-FIT TESTS

## 1 INTRODUCTION

Goodness-of-fit (gof) tests are designed to measure the compatibility of a random sample with a theoretical probability distribution function (pdf). The null hypothesis $`H_0`$ is that the sample follows the pdf. Under the assumption that H0 applies, the fraction of wrongly rejected experiments - the probability of committing an error of the first kind - is fixed to typically a few percent. A test is considered powerful if the probability of accepting $`H_0`$ when $`H_0`$ is wrong - the probability of committing an error of the second kind - is low. Of course, without specifying the alternatives, the power cannot be quantified.

A discrepancy between a data sample and the theoretical description can be of different origin. The problem may be in the theory which is wrong or the sample may be biased by measurement errors or by background contamination. In natural sciences we mainly have the latter situation. Even though the statistical description is the same in both cases the choice of the specific test may be different. In our applications we are mainly confronted with “slowly varying” deviations between data and theoretical description whereas in other fields where for example time series are investigated, “high frequency” distortions are more likely.

Goodness-of-fit tests are based on classical statistical methods and are closely related to classical interval estimation, but they contain also Bayesian elements. Those, however, are only related with some prejudice on the alternative hypothesis which affects the purity of the accepted decisions and not the error of the first kind.

The power of one dimensional tests is not always invariant against transformations of the variates. In more than one dimension (number of variates), an invariant description is not possible.

Tests are classified in distribution dependent and distribution free tests. The former are adapted to special pdfs like Gaussian, exponential or uniform distributions. We will restrict our discussion mainly to distribution free tests and tests which can be adapted to arbitrary distributions. Here we distinguish tests applied to binned data and binning free tests. The latter are in principle preferable but so far they are almost exclusively limited to one dimensional distributions. A further distinction concerns the alternative hypothesis. Usually, it is not restricted but there exist also tests where it is parametrized.

Physicists tend to be content with $`\chi ^2`$ tests which are not necessarily optimum in all cases. A very useful and comprehensive survey of goodness-of-fit tests can be found in Ref.  from 1986. Since then, some new developments have occurred and the increase in computing power has opened the possibility to apply more elaborate tests.

In Section 2 we summarize the most important tests. To keep this article short we do not discuss tests based on the order statistic and spacing tests. In Section 3 we introduce two new tests, the three region test and the energy test. To compare the tests we apply them in Section 4 to some specific alternative hypotheses. We do not consider explicitly composite hypotheses.

## 2 SOME RELEVANT TESTS

### 2.1 Chi-squared test

The $`\chi ^2`$ test is closely connected to least square fits with the difference that the hypothesis is fixed. The test statistic is

$$\chi ^2=\underset{i=1}{\overset{B}{}}\frac{(Y_it_i)^2}{\delta _i^2}$$

with the random variable $`Y_i`$, $`t_i`$ the expectation $`E(Y_i)`$ and $`\delta _i^2`$ the expectation $`E((Y_it_i)^2)`$. Obviously the expectation value of $`\chi ^2`$ is

$$E(\chi ^2)=B$$

In the Gaussian approximation $`Y_i`$ follows a Gaussian with mean $`t_i`$ and variance $`\delta _i^2`$ and the test statistic follows a $`\chi ^2`$ distribution function $`F_B(\chi ^2)`$ with $`B`$ degrees of freedom. The probability of an error of the first kind $`\alpha `$ (significance level, p-value) defines $`\chi _0^2`$ with $`F_B(\chi _0^2)=1\alpha `$. The null hypothesis is rejected if in an actual experiment we find $`\chi ^2>\chi _0^2`$.

We obtain the Pearson test when the random variables $`N_i`$ are Poisson distributed.

$$\chi ^2=\underset{i=1}{\overset{B}{}}\frac{(N_it_i)^2}{t_i}$$

In the large sample limit the test statistic $`\chi ^2`$ has approximately a $`\chi ^2`$ distribution with $`B`$ degrees of freedom. When we have a histogram, where a total of $`N`$ events are distributed according to a multinomial distribution among the $`B`$ bins with probabilities $`p_i`$ we get

$$\chi ^2=\underset{i=1}{\overset{B}{}}\frac{(N_iNp_i)^2}{Np_i}$$

which again follows asymptotically a $`\chi ^2`$ distribution, this time with $`B1`$ degrees of freedom. The reduced number of degrees of freedom is due to the constraint $`\mathrm{\Sigma }N_i=N`$.

Nowadays, the distribution function of the test statistic can be computed numerically without much effort. The $`\chi ^2`$ test then can also be applied to small samples. The Gaussian approximation is no longer required.

The $`\chi ^2`$ test is very simple and needs only limited computational power. A big advantage compared to most of the other methods is that is can be applied to multidimensional histograms. There are however also serious drawbacks:

* Its power in detecting slowly varying deviations of a histogram from predictions is rather poor due to the neglect of possible correlations between adjacent bins.
* Binning is required and the choice of the binning is arbitrary.
* When the statistics is low or the number of dimensions is high, the event numbers per bin may be low. Then the asymptotic properties are no longer valid and systematic deviations are hidden by statistical fluctuations.

There are proposals to fix the bin widths by the requirement of equal number of expected entries per bin. This is not necessarily the optimum choice . Often there are outliers in regions where no events are expected which would be hidden in wide bins.

For the number of bins a dependence on the sample size $`n`$

$$B=2n^{2/5}$$

is proposed in Ref. . Our experience is that in most experiments the number of bins is chosen too high. The sensitivity to slowly varying deviations roughly goes with $`B^{1/4}`$ . In multidimensional cases the power of the test often can be increased by applying it to the marginal distributions.

There is a whole class of $`\chi ^2`$ like tests. Many studies can be found in the literature. The reader is referred to Ref. .

### 2.2 Binning-free empirical distribution function tests

The tests described in this section have been taken from the article by Stephens in Ref. .

Supposing that a random sample of size $`n`$ is given, we form the order statistic $`X_1<X_2<\mathrm{}<X_n`$. We consider the empirical distribution function (EDF)

$$F_n(x)=\frac{\text{\# of observations }x}{n}$$

or

$`\begin{array}{cc}F_n(x)=0\hfill & x<X_1\hfill \\ F_n(x)=\frac{i}{n}\hfill & X_ix<X_{i+1}\hfill \\ F_n(x)=1\hfill & X_nx\hfill \end{array}`$

$`F_n(x)`$ is a step function which is to be compared to the distribution $`F(x)`$ corresponding to $`H_0`$.

The EDF is consistent and unbiased. The tests discussed in this section are invariant under transformation of the random variable. Because of this feature, we can transform the distribution to the uniform distribution and restrict our discussion to the latter.

#### 2.2.1 Probability integral transformation

The probability integral transformation (PIT)

$$Z=F(X)$$

transforms a general pdf $`f(X)`$ of $`X`$ into the uniform distribution $`f^{}(Z)`$ of $`Z`$ .

$`f^{}(Z)`$ $`=`$ $`1;\text{ }0Z1`$
$`F^{}(Z)`$ $`=`$ $`Z`$
The underlying idea of this transformation is that the new EDF of $`Z`$, $`F^{}(Z)`$ is extremely simple and that it conserves the distribution of the test quantities discussed in this section. It is easily seen that

$$F_n(x)F(x)=F_n^{}(z)z$$

Note, however, that the PIT does not necessarily conserve all interesting features of the gof problem. Resolution effects are washed out and for example in a lifetime distribution, an excess of events at small and large lifetimes may be judged differently but are treated similarly after a PIT. It is not logical to select specific gofs for specific applications but to transform all kinds of pdfs to the same uniform distribution. The PIT is very useful because it permits standardization but one has to be aware of its limitations.

#### 2.2.2 Supremum statistics

The maximum positive (negative) deviation of $`F_n(x)`$ from $`F(x)`$ $`D_+`$ ($`D_{}`$) (see Fig. 1) are used as tests statistics. Kolmogorov (Kolmogorov-Smirnov test) has proposed to use the maximum absolute difference. Kuiper uses the sum $`V=D_++D_{}`$. This test statistic is useful for observations “on the circle” for example for azimuthal distributions where the zero angle is a matter of definition.

$`D_+`$ $`=`$ $`\underset{x}{sup}\left\{F_n(x)F(x)\right\}`$
$`D_{}`$ $`=`$ $`\underset{x}{sup}\left\{F(x)F_n(x)\right\}`$
$`D`$ $`=`$ $`\underset{x}{sup}\left\{\left|F_n(x)F(x)\right|\right\}\text{(Kolmogorov)}`$
$`V`$ $`=`$ $`D_++D_{}\text{(Kuiper)}`$
The supremum statistics are invariant under the PIT.

#### 2.2.3 Integrated deviations - quadratic statistics

The Cramer-von Mises family of tests measures the integrated quadratic deviation of $`F_n(x)`$ from $`F(x)`$ suitably weighted by a weighting function $`\psi `$:

$$Q=n_{\mathrm{}}^{\mathrm{}}\left[F(x)F_n(x)\right]^2\psi (x)𝑑F$$

In the standardized form we have

$$Q=n_0^1\left[zF_n^{}(z)\right]^2\psi (z)𝑑z$$
(1)
Since the construction of $`F_n^{}(Z)`$ includes already an integration, $`F_n^{}(z_i)`$ and $`F_n^{}(z_k)`$ are not independent and the additional integration in Equation 1 is not obvious.

With $`\psi _{CvM}=1`$ we get the Cramer-von Mises statistic $`W^2`$ and $`\psi _{AD}=\left[z(1z)\right]^1`$ leads to the Anderson-Darling statistic $`A^2`$.

$`W^2`$ $`=`$ $`n{\displaystyle _0^1}\left[zF_n^{}(z)\right]^2𝑑z\text{(Cramer-von Mises)}`$
$`A^2`$ $`=`$ $`n{\displaystyle _0^1}{\displaystyle \frac{\left[zF_n^{}(z)\right]^2}{z(1z)}}𝑑z\text{(Anderson - Darling)}`$
The Anderson-Darling statistic $`A^2`$ weights strongly deviations near $`z=0`$ and $`z=1`$. This is justified because there the experimental deviations are small due to the constraints $`\left[zF_n^{}(z)\right]=0`$ at $`z=0`$ and $`z=1`$.

Watson has proposed a quadratic statistic on the circle:

$$U^2=n_0^1\left\{F_n^{}(z)z_0^1\left[F_n^{}(z)z\right]𝑑z\right\}^2𝑑z\text{(Watson)}$$

### 2.3 The Neyman statistic test

This test is different from all previously discussed tests. It parametrizes the alternative hypothesis and applies the likelihood ratio test. The alternative hypothesis corresponds to a pdf of the exponential family:

$$g_k(z)=C(\theta _1,\theta _2\mathrm{}\theta _k)\mathrm{exp}\left[\underset{i=1}{\overset{k}{}}\theta _i\pi _i(z)\right]$$

$`g_k(z)`$ are smooth alternatives to uniformity. The functions $`\pi _i`$ are Legendre polynomials of order $`i`$, $`\theta _i`$ are free parameters and $`C`$ is a normalization function. The number $`k`$ of parameters is selected by the user.

The likelihood ratio leads to the test statistic

$$N_k=\frac{1}{n}\underset{i=1}{\overset{k}{}}\left(\underset{j=1}{\overset{n}{}}\pi _i(z_j)\right)^2$$

Asymptotically, for large values of $`N_k`$, $`N_k`$ is distributed according to the $`\chi ^2`$ distribution with $`k`$ degrees of freedom.

## 3 NEW TESTS

### 3.1 Three region test

Often experimental distributions are biased by an excess or lack of events in certain regions of the random variable. We have designed a test which subdivides the variable space into three pieces, containing $`n_1,n_2,n_3=nn_1n_2`$ events, such that the deviation between data and prediction from $`H_0`$ is maximum. The test quantity is

$$O=\underset{n_1,n_2}{sup}\left\{w_1(n_1np_1)^2+w_2(n_2np_2)^2+w_3(n_3np_3)^2\right\}$$

where $`np_k`$ are the expectation values and $`w_k`$ weights depending on $`np_k`$. The specific choice

$`w_k`$ $`=`$ $`{\displaystyle \frac{1}{np_k}}`$
$`O_\chi `$ $`=`$ $`\underset{n_1,n_2}{sup}\left\{{\displaystyle \frac{(n_1np_1)^2}{np_1}}+{\displaystyle \frac{(n_2np_2)^2}{np_2}}+{\displaystyle \frac{(n_3np_3)^2}{np_3}}\right\}`$
maximizes $`\chi ^2`$ of the three bins. In the comparison below we have chosen weights equal to one.

Of course the test can be generalized to a higher number of subregions.

### 3.2 Minimum energy test

#### 3.2.1 The idea

Let us assume that we have a continuous charge distribution $`\rho (\stackrel{}{r})`$ of positive electric charges and a sample of negative point charges with total charge equal to minus the integrated positive charge. The potential energy is minimum when the negative point charges follow $`\rho `$. Then, up to effects due to the discrete nature of the point charges, the charge density is zero everywhere. Any displacement of charges would increase the energy. We use this property to construct a binning free test procedure.

We simulate the theoretical distribution by $`m`$ charges of charge $`1/m`$ each. Usually, these charges are distributed using a Monte Carlo simulation. To the $`n`$ experimental sample points (data points) we associate charges $`1/n`$. The test quantity $`\varphi `$ corresponds to the potential energy. It contains two terms $`\varphi _1,\varphi _2`$ corresponding to the interaction of the experimental charges with each other and to the interaction of the experimental charges with the positive simulation charges.

$`\varphi _1`$ $`={\displaystyle \frac{1}{n^2}}{\displaystyle \underset{i<j}{}}R(d_{ij})`$ (2)
$`\varphi _2`$ $`={\displaystyle \frac{1}{nm}}{\displaystyle \underset{i,j}{}}R(t_{ij})`$ (3)
$`\varphi `$ $`=\varphi _1+\varphi _2`$ (4)
Here $`d_{ij}`$ is the distance between two data points and $`t_{ij}`$ is the distance between a data point and a simulation point and $`R`$ is a correlation function defined below. The sums run over all combinations.

Remark: The minimum energy requirement for the equality of experimental and theoretical distribution is strictly correct only when the number $`m`$ of simulation charges is equal to the number $`n`$ of experimental charges. For the general case with a continuous theoretical distribution or simulation sample and experimental sample of different size, the optimum agreement of the two distributions is not well defined and there is a slight dependence of the minimum energy configuration on the correlation function. This is however a purely academic problem, the test statistic $`\varphi `$ remains a powerful indicator for an incompatibility of the experimental sample with $`H_0`$.

#### 3.2.2 The correlation function

We note that the minimum energy configuration does not depend on the application of the one-over-distance power law of electrostatics. We may apply a wide class of correlation functions $`R(r)`$ with the only requirement that $`R`$ has to decrease monotonically with the distance $`r`$.

We have investigated three different types of correlation functions, power laws, a logarithmic dependence and Gaussians.

$`R_\kappa (r)`$ $`={\displaystyle \frac{1}{r^\kappa }}`$ (5)
$`R_l(r)`$ $`=\mathrm{ln}r`$ (6)
$`R_s(r)`$ $`=e^{r^2/(2s^2)}`$ (7)
The first type is motivated by the analogy to electrostatics, the second is long range and the third emphasizes a limited range for the correlation between different points. The power $`\kappa `$ of the denominator in Equ. 5 and the parameter $`s`$ in Equ. 7 may be chosen differently for different dimensions of the sample space and different applications. For long range distortions a small value of $`\kappa `$ around $`0.1`$ is recommended. For short range deviations the test quantity with larger values around $`0.3`$ is more sensitive.

The inverse power law and the logarithm have a singularity at $`r`$ equal to zero. Very small distances, however, should not be weighted too strongly since distortions with sharp peaks are not expected and usually inhibited by the finite experimental resolution. We eliminate the singularity by introducing a lower cutoff $`d_{\mathrm{min}}`$ for the distances $`d`$ and $`t`$. Distances less than $`d_{\mathrm{min}}`$ are replaced by $`d_{\mathrm{min}}`$. The value of this parameter is not critical, it should be of the order of the average distance $`d`$ in the regions where the $`f_0`$ is maximum and not less than the experimental resolution.

The energy test with Gaussian correlation function is closely related to the Pearson $`\chi ^2`$ test. A more detailed description of the energy test is presented in Ref. .

### 3.3 Comparison of uni-variate tests

We have tested the null hypothesis of a uniform distribution in the interval $`[0,1]`$ using a uniform distribution contaminated by the background distributions displayed in Figure 2.

Background hypothesis A modifies the mean, hypotheses B, C change the variance of $`H_0.`$

The power of various tests described above is presented in Figure 3.

As expected, non of the tests is optimum for all kind of distortions. Several tests perform better than the $`\chi ^2`$ test. The Neyman test, the Anderson-Darling test and the Kolmogorov-Smirnov test are sensitive to a shift in the mean. Anderson’s test detects especially deviations at the borders of the interval. Watson’s and Kuiper’s tests are useful for the detection of distortions of the variance. The two new tests compare favorably with the standard ones.

## 4 MULTIVARIATE TESTS

The Mardia test  and the Neyman smooth test  can be used to investigate two-dimensional Gaussian distributions. The only distribution free test known to us which is independent of the dimensions of the variate space is the $`\chi ^2`$ test. A generalized Kolmogorov-Smirnov test  depends on the ordering scheme of the variates. The binning free energy test developed by us is also independent of the number of variates, however the distribution of test statistic has to be computed for the specific sample distribution under study.

### 4.1 Comparison of multivariate tests

We have used a two-dimensional Gaussian null hypothesis and contaminated the sample with the background distributions shown in Figure 4.

The power of the two Mardia tests, the Neyman smooth test and the energy test with logarithmic and Gaussian correlation function is presented in Figure 5.

In most cases the two energy tests perform better than the alternatives even though those have been designed for a Gaussian null hypothesis.

### 4.2 Example: Comparing experimental data to a Monte Carlo prediction

In Figure 6, left hand side, we compare the position and momentum of a few $`J/\psi `$ decay tracks to a Monte Carlo simulation.The right hand plot compares the energy computed from the distribution on the left hand side to a Monte Carlo simulation of the null hypothesis. The experimental point, indicated by the arrow, is larger than all Monte Carlo values. Apparently, the data do not follow the prediction.

## 5  CONCLUSIONS

The $`\chi ^2`$ test suffers from the requirement to choose a binning. In one dimension it should be replaced by the well established binning free tests like the Kolmogorov test. The choice of a specific test has to depend on the expected kind of possible distortion of the theoretical distribution. For a localized background we advise to use the new three region test. For multivariate applications the new energy test is an attractive alternative to the $`\chi ^2`$ test.


text:
# 1 Introduction

## 1  Introduction

The deep inelastic scattering (DIS) leptons on hadrons is the basical process to study the values of the parton distribution functions (PDF) which are universal (after choosing of factorization and renormalization schemes) and can be used in other processes. The accuracy of the present data for deep inelastic structure functions (SF) reached the level at which the $`Q^2`$-dependence of logarithmic QCD-motivated terms and power-like ones may be studied separately (for a review, see the recent papers  and references therein).

In the present letter we sketch the results of our analysis  at the next-to-leading order (NLO) of perturbative QCD for the most known DIS SF $`F_2(x,Q^2)`$ 111Here $`Q^2=q^2`$ and $`x=Q^2/(2pq)`$ are standard DIS variables, where $`q`$ and $`p`$ are photon and hadron momentums, respectively. taking into account experimental data - of SLAC, NM, BCDMS and BFP Collaborations. We stress the power-like effects, so-called twist-4 (i.e. $`1/Q^2`$) contributions. To our purposes we represent the SF $`F_2(x,Q^2)`$ as the contribution of the leading twist part $`F_2^{pQCD}(x,Q^2)`$ described by perturbative QCD, when the target mass corrections are taken into account (and coincides with $`F_2^{tw2}(x,Q^2)`$ when the target mass corrections are withdrawn), and the nonperturbative part (“dynamical” twist-four terms):

$$F_2(x,Q^2)F_2^{full}(x,Q^2)=F_2^{pQCD}(x,Q^2)\left(1+\frac{\stackrel{~}{h}_4(x)}{Q^2}\right),$$
(1)
where $`\stackrel{~}{h}_4(x)`$ is magnitude of twist-four terms.

Contrary to standard fits (see, for example, \- ) when the direct numerical calculations based on Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation  are used to evaluate structure functions, we use the exact solution of DGLAP equation for the Mellin moments $`M_n^{tw2}(Q^2)`$ of SF $`F_2^{tw2}(x,Q^2)`$:

$$M_n^k(Q^2)=_0^1x^{n2}F_2^k(x,Q^2)𝑑x(k=full,pQCD,tw2,\mathrm{})$$
(2)
and the subsequent reproduction of $`F_2^k(x,Q^2)`$ at every needed $`Q^2`$-value with help of the Jacobi Polynomial expansion method  (see similar analyses at the NLO level  and at the next-next-to-leading order (NNLO) level and above .

In this letter we do not present exact formulae of $`Q^2`$-dependence of SF $`F_2`$ which are given in . We note only that the moments $`M_n^{tw2}(Q^2)`$ at some $`Q_0^2`$ is theoretical input of our analysis and the twist-four term $`\stackrel{~}{h}_4(x)`$ is considered as a set of free parameters (one constant $`\stackrel{~}{h}_4(x_i)`$ per $`x_i`$-bin): $`\stackrel{~}{h}_4^{free}(x)=_{i=1}^I\stackrel{~}{h}_4(x_i)`$, where $`I`$ is the number of bins.

## 2  Fits of $`F_2`$: procedure

Having the QCD expressions for the Mellin moments $`M_n^k`$ we can reconstruct the SF $`F_2^k(x)`$ as

$$F_2^{k,N_{max}}(x,Q^2)=x^a(1x)^b\underset{n=0}{\overset{N_{max}}{}}\mathrm{\Theta }_n^{a,b}(x)\underset{j=0}{\overset{n}{}}c_j^{(n)}(\alpha ,\beta )M_{j+2}^k\left(Q^2\right),$$
(3)
where $`\mathrm{\Theta }_n^{a,b}`$ are the Jacobi polynomials 222We note here that there is similar method , based on Bernstein polynomials. The method has been used in the analyses at the NLO level in  and at the NNLO level in . and $`a,b`$ are fitted parameters.

First of all, we choose the cut $`Q^21`$ GeV2 in all our studies. For $`Q^2<1`$ GeV2, the applicability of twist expansion is very questionable. Secondly, we choose quite large values of the normalization point $`Q_0^2`$: our perturbative formulae should be applicable at the value of $`Q_0^2`$. Moreover, the higher order corrections $`\alpha _s^k(Q_0^2)`$ and $`(\alpha _s(Q^2)\alpha _s(Q_0^2))^k`$ ($`k2`$) should be less important at higher $`Q_0^2`$ values.

We use MINUIT program  for minimization of $`\chi ^2(F_2)=\left|(F_2^{exp}F_2^{teor})/\mathrm{\Delta }F_2^{exp}\right|^2`$. We consider free normalizations of data for different experiments. For the reference, we use the most stable deuterium BCDMS data at the value of energy $`E_0=200`$ GeV ($`E_0`$ is the initial energy lepton beam). Using other types of data as reference gives negligible changes in our results. The usage of fixed normalization for all data leads to fits with a bit worser $`\chi ^2`$.

## 3  Results of fits

Hereafter we choose $`Q_0^2`$ = 90 GeV2 ($`Q_0^2`$ = 20 GeV2) for the nonsinglet (combine nonsinglet and singlet) evolution, that is in good agreement with above conditions. We use also $`N_{max}=8`$.

### 3.1  BCDMS $`{}_{}{}^{12}C+H_2+D_2`$ data

We start our analysis with the most precise experimental data  obtained by BCDMS muon scattering experiment at the high $`Q^2`$ values. The full set of data is 762 (607) points (for the bounded $`x`$ range: $`x0.25`$).

It is well known that the original analyses given by BCDMS Collaboration itself (see also Ref. ) lead to quite small values $`\alpha _s(M_Z^2)=0.113`$. Although in some recent papers (see, for example, ) more higher values of the coupling constant $`\alpha _s(M_Z^2)`$ have been observed, we think that an additional reanalysis of BCDMS data should be very useful.

Based on study  we proposed in  that the reason for small values of $`\alpha _s(M_Z^2)`$ coming from BCDMS data was the existence of the subset of the data having large systematic errors. We studied this subject by introducing several so-called $`Y`$-cuts 333Hereafter we use the kinematical variable $`Y=(E_0E)/E_0`$, where $`E`$ is scattering energies of lepton. (see ). Excluding this set of data with large systematic errors leads to essentially larger values of $`\alpha _s(M_Z^2)`$ and very slow dependence of the values on the concrete choice of the $`Y`$-cut (see below).

We use the following $`x`$-dependent $`Y`$-cuts:

$`y0.14\text{ when }0.3<x0.4,y0.16\text{ when }0.4<x0.5`$
$`yY_{cut3}\text{ when }0.5<x0.6,yY_{cut4}\text{ when }0.6<x0.7`$
$`yY_{cut5}\text{ when }0.7<x0.8`$ (4)
and several $`N`$ sets for the cuts at $`0.5<x0.8`$:

The systematic errors for BCDMS data were given  as multiplicative factors to be applied to $`F_2(x,Q^2)`$: $`f_r,f_b,f_s,f_d`$ and $`f_h`$ are the uncertainties due to spectrometer resolution, beam momentum, calibration, spectrometer magnetic field calibration, detector inefficiencies and energy normalization, respectively. For this study each experimental point of the undistorted set was multiplied by a factor characterizing a given type of uncertainties and a new (distorted) data set was fitted again in agreement with our procedure considered in the previous section. The factors ($`f_r,f_b,f_s,f_d,f_h`$) were taken from papers  (see CERN preprint versions in ). The $`\alpha _s`$ values for the distorted and undistorted sets of data are given in the Figs. 1 and 2 (for the cases of nonsinglet and complete evolutions, respectively) together with the total systematic error estimated in quadratures.

From the Figs. 1 and 2 we can see that the $`\alpha _s`$ values are obtained for $`N=1÷6`$ of $`Y_{cut3}`$, $`Y_{cut4}`$ and $`Y_{cut5}`$ are very stable and statistically consistent. The case $`N=6`$ of the Table 1 reduces the systematic error in $`\alpha _s`$ by factor $`1.8`$ and increases the value of $`\alpha _s`$, while increasing the statistical error on the 30%.

After the cuts have been implemented (we use the set $`N=6`$ of the Table 1), we have 590 (452) points (for the bounded $`x`$ range: $`x0.25`$). Fitting them in agreement with the same procedure considered in the previous Section, we obtain the following results:

from fits, based on nonsinglet evolution (i.e. when $`x0.25`$):

$`\alpha _s(M_Z^2)`$ $`=`$ $`0.1153\pm 0.0013\text{(stat)}\pm 0.0022\text{(syst)}\pm 0.0012\text{(norm)},`$
from fits, based on combined singlet and nonsinglet evolution:

$`\alpha _s(M_Z^2)`$ $`=`$ $`0.1175\pm 0.0014\text{(stat)}\pm 0.0020\text{(syst)}\pm 0.0011\text{(norm)},`$ (5)
where hereafter the symbol “norm” marks the error of normalization of experimental data.

The results are agree each other within considered errors. In Ref.  we have also analyzed the combine SLAC, NM and BFP data and found good agreement with (5). So, we have a possibility to fit together all the data. It is the subject of the following subsection.

### 3.2  SLAC, BCDMS, NM and BFP data

After these $`Y`$-cuts have been incorporated (with $`N=6`$) for BCDMS data, the full set of combine data is 1309 (797) points (for the bounded $`x`$ range: $`x0.25`$).

To verify the range of applicability of perturbative QCD, we analyze firstly the data without a contribution of twist-four terms, i.e. when $`F_2=F_2^{pQCD}`$. We do several fits using the cut $`Q^2Q_{cut}^2`$ and increase the value $`Q_{cut}^2`$ step by step. We observe good agreement of the fits with the data when $`Q_{cut}^210÷15`$ GeV2 (see the Figs. 3 and 4). Later we add the twist-four corrections and fit the data with the standard cut $`Q^21`$ GeV2. We have find very good agreement with the data. Moreover the predictions for $`\alpha _s(M_Z^2)`$ in both above procedures are very similar (see the Figs. 3 and 4).

The results of the fits are compiled in Summary (see Eqs. (6)-(9)).

## 4  Summary

We have demonstrated several steps of our study  of the $`Q^2`$-evolution of DIS structure function $`F_2`$ fitting all modern fixed target experimental data.

From the fits we have obtained the value of the normalization $`\alpha _s(M_Z^2)`$ of QCD coupling constant. First of all, we have reanalyzed the BCDMS data cutting the range with large systematic errors. As it is possible to see in the Fig. 1, the value of $`\alpha _s(M_Z^2)`$ rises strongly when the cuts of systematics were incorporated. In another side, the value of $`\alpha _s(M_Z^2)`$ does not dependent on the concrete type of the cut within modern statistical errors.

We have found that at $`Q^210÷15`$ GeV2 the formulae of pure perturbative QCD (i.e. twist-two approximation together with target mass corrections) are in good agreement with all data. 444We note that at small $`x`$ values, the perturbative QCD works well starting with $`Q^2=1.5÷2`$ GeV2 and higher twist corrections are important only at very low $`Q^2`$: $`Q^20.5`$ GeV2 (see  and references therein). As it is was observed in  (see also discussions in ) the good agreement between perturbative QCD and experiment seems connect with large effective argument of coupling constant at low $`x`$ range. The results for $`\alpha _s(M_Z^2)`$ are very similar (see ) for the both types of analyses: ones, based on nonsinglet evolution, and ones, based on combined singlet and nonsinglet evolution. They have the following form:

* from fits, based on nonsinglet evolution:

$`\alpha _s(M_Z^2)`$ $`=`$ $`0.1170\pm 0.0009\text{(stat)}\pm 0.0019\text{(syst)}\pm 0.0010\text{(norm)},`$ (6)
* from fits, based on combined singlet and nonsinglet evolution:

$`\alpha _s(M_Z^2)`$ $`=`$ $`0.1180\pm 0.0013\text{(stat)}\pm 0.0021\text{(syst)}\pm 0.0009\text{(norm)},`$ (7)

When we have added twist-four corrections, we have very good agreement between QCD (i.e. first two coefficients of Wilson expansion) and data starting already with $`Q^2=1`$ GeV2, where the Wilson expansion should begin to be applicable. The results for $`\alpha _s(M_Z^2)`$ coincide for the both types of analyses: ones, based on nonsinglet evolution, and ones, based on combined singlet and nonsinglet evolution. They have the following form:

* from fits, based on nonsinglet evolution:

$`\alpha _s(M_Z^2)`$ $`=`$ $`0.1174\pm 0.0007\text{(stat)}\pm 0.0019\text{(syst)}\pm 0.0010\text{(norm)},`$ (8)
* from fits, based on combined singlet and nonsinglet evolution:

$`\alpha _s(M_Z^2)`$ $`=`$ $`0.1177\pm 0.0007\text{(stat)}\pm 0.0021\text{(syst)}\pm 0.0009\text{(norm)},`$ (9)

Thus, there is very good agreement (see Eqs. (6), (7), (8) and (9)) between results based on pure perturbative QCD at quite large $`Q^2`$ values (i.e. at $`Q^210÷15`$ GeV2) and the results based on first two twist terms of Wilson expansion (at $`Q^21`$ GeV2, where the Wilson expansion should be applicable).

We would like to note that we have good agreement also with the analysis  of combined H1 and BCDMS data, which has been given by H1 Collaboration very recently. Our results for $`\alpha _s(M_Z^2)`$ are in good agreement also with the average value for coupling constant, presented in the recent studies (see  and references therein) and in famous Altarelli and Bethke reviews .

The last result (9) based on all data with $`Q^21`$ GeV2 can be considered as “best value” for the coupling constant $`\alpha _s(M_Z^2)`$ coming in our analysis.

Acknowledgments. The study is supported in part by the Heisenberg-Landau program. Authors would like to express their sincerely thanks to the Organizing Committee of Vth Int. Conference “Renormalization group 2002” for the kind invitation, the financial support at such remarkable Conferences, and for fruitful discussions. A.V.K. was supported in part by Alexander von Humboldt fellowship and INTAS grant N366.