text:
# Exact solutions of generalized Hubbard Hamiltonian for diamond vacancies

## I introduction

Vacancy is one of the simplest and familiar intrinsic point defects in semiconductors and its effects have attracted high technological and theoretical interests.1-3 A vacancy in diamond in contrast to the situation in most semiconductor or insulator crystals is stable in room temperature.4 Since the formation energy of vacancies in diamond is high ($`67`$ eV) 5,6 they have been studied mainly in irradiated diamond. Optical absorption studies are used for determination of electronic excitation energies of vacancies7,8 and EPR, ENDOR experiments have been used for measuring the spin and spatial symmetry of electronic states.9-12 The most famous optical absorption lines in irradiated diamond are $`GR1`$ (with zero phonon line at $`1.673`$ eV)13 and $`ND1`$ (with zero phonon line at $`3.149`$ eV)14. The former is attributed to a neutral vacancy15 and the latter to a negatively charged one.16
In addition to these two famous lines many other absorption lines have been reported such as: N3 ($`2.985`$ eV)8, H3 ($`2.463`$ eV)17, NV ($`1.945`$ eV)18 which have been attributed to vacancy-nitrogen19,20, also GR2 to GR8 lines7 which have been attributed to the natural-charged vacancy complexes.21 Despite the simplicity of the system, different models have been suggested to explain the properties of these optical absorption lines. Generally there have been two main approaches to this problem, The first group of approaches are localized models22-26 and the second ones are extended models.27-31 In the first group, the electrons of the broken bonds which are in the tetrahedrally symmetric local potential of the vacant site are considered as an isolated molecule (many body approaches). In the second group i.e. extended models, the effect of lattice is more important on the vacancy and it is considered as a missing atom in a supercell with few hundreds of atoms (single particle approaches). Similarity between these two categories is their group theoretical aspects for accounting symmetries of the system while their main difference is in the importance of e-e correlation effects. Lannoo and Bourgoin 32have developed a model with four parameters which briliantly explaines the validity range of the single and many particle approaches for this problem. According to Messmer and Watkins27 a successful model for describing vacancy energy states should explain the following points: The position of the defect energy levels relative to the valence band edge, the electronic wavefunction of defect electrons for comparing with experimental data, lattice relaxation and distortion in neighborhood of vacant site and it must be the basis of a practical computational scheme. The localized approaches have relatively successful answers to these questions.
The first suggested model in this regard is the famous model of the Coulson and Kearsly33. Their model suggests a computational scheme that considers e-e interaction in the system accurately. They used group theoretical methods to manually construct symmetry and spin adapted wavefunctions for including spin and symmetry considerations of vacancies. By applying configuration interaction technique to the calculation results, they predicted qualitatively good results for electronic states of the vacancies, their model used two semi-empirical parameters to obtain satisfactory result for energy of the $`GR1`$ transition. By using these semi-empirical parameters their predictions for $`ND1`$ transiton and also the posision of the low lying $`{}_{}{}^{3}T_{1}^{}`$ state disagrees with avaiable experimental data. Up to now this model has been the most successful computational scheme for describing physical properties of the optical absorption lines of vacancies in diamond,3,21 and the next suggested localized models21-26 are based on the framework of this model.
Localized models, unlike the nowadays commonly used density functional theory (DFT) methods which pay more attention on the ground state information, can also give the excitation energies of the vacancies beyond the one electron approximation.
For improving the quantitative agreement between theory and experiment we introduce a new approach for localized models. This approach is based on a generalized Hubbard Hamiltonian for electrons of the vacancy system with atomic orbital bases and uses computer adaptable many body techniques instead of usually used molecular orbital method (fourth item of messmer and watkins). In this paper at first we review the computational scheme of the present model and discuss the advantages of the new notation which is used for this problem. Then we apply this new scheme to solve generalized form of Hubbard Hamiltonian for neutral and charged vacancy systems. For evaluation of solutions of this Hamiltonian we start with parameters which were already reported. Our using value for the two semi-empirical parameters slightly differs from already rpoeted values but other six parameters are as the same as the theory calculation. We report the effect of the hopping parameter on the energy spectrum of this Hamiltonian by varying this parameter.
The model predictions for $`GR1`$ and $`ND1`$ transitions with justified parameters and comparison the results with experimental data will be discussed. The location of the low lying $`{}_{}{}^{3}T_{1}^{}`$ state is predicted and we disscuss its comparision with other theoretical reports and experimental evidences. The results of our model can justify the semiempirical values that already were reported.33,41 We also will show that, this model is independent of familiar configuration interaction technique and the electronic configurations mixing of each state comes out in a natural manner from the formalism. The model predicts whole of the exact eigenstates and enery levels of the systems with the justified semi-empirical values. Finally the new information concerning electronic configurations of the ground and excited states of $`GR1`$ transition are reported which can be used by other approaches to this problem.
In our calculation the effects of lattice relaxation were ignored, which according to recent theoretical works are small.34-36

## II Computational model

The present computational model takes into account electron correlation effects accurately and gives the correct spin and spatial symmetry of eigenstates. It is also suitable for small atomic cluster systems or molecules provided the hopping, Coulombic and exchange integrals are known. In the localized model framework, it is assumed that physical properties of neutral (charged) vacancy depend on four (five) electrons of adjacent dangling bounds. The exact Hamiltonian of such a system can be written as Eq. (1). A more simplified form of this Hamiltonian recently is applied to carbon nanotube systems 36,37 .

$`H={\displaystyle \underset{ij\sigma }{}}t_{ij}c_{i\sigma }^{}c_{j\sigma }+{\displaystyle \frac{1}{2}}{\displaystyle \underset{ijlm\sigma \sigma ^\mathrm{`}}{}}V_{ijlm}c_{i\sigma }^{}c_{j\sigma ^\mathrm{`}}^{}c_{m\sigma ^\mathrm{`}}c_{l\sigma }`$ (1)
Where $`i,j,l,m`$ are atomic sites indices which go from $`1`$ to $`4`$ , $`\sigma `$ is spin index which is $`1/2`$ or $`+1/2`$ , $`c_{i\sigma }^{}`$ and $`c_{i\sigma }`$ are operators which create and annihilate electrons with specified spin $`\sigma `$ on site $`i`$. $`t_{ij}`$ and $`V_{ijlm}`$ parameters are hopping and $`ee`$ interaction overlap integrals which are defined as follows:

$`t_{ij}={\displaystyle 𝑑\stackrel{}{r}\psi _i^{}(\stackrel{}{r})(\frac{\mathrm{}^2}{2m}^2+V(r))\psi _j(\stackrel{}{r})}`$ (2)
$`V_{ijlm}={\displaystyle 𝑑\stackrel{}{r}𝑑\stackrel{}{r^{}}\psi _i^{}(\stackrel{}{r})\psi _j^{}(\stackrel{}{r})\frac{1}{|rr^{}|}\psi _l(\stackrel{}{r^{}})\psi _m(\stackrel{}{r^{}})}`$ (3)
in the above integrals $`\psi _i`$ is single electronic wave function in a specified site. The advantage of using second quantized form of Hamiltonian for this problem is its capability to include the spatial symmetry and spin information of the system in the form of Hamiltonian, hence the eigenstates have automatically the desired spin and symmetry properties. Also the effect of the changing each parameter directly can be observed on the electronic energy levels and also wave functions. This new approach uses atomic orbital bases to solve the Hamiltonian. In the previous localized models it was common as a start point, to manually construct symmetry and spin adapted molecular wave functions as correct basis for accounting symmetry considerations of the system then these basis should be used in evaluating Hamiltonian matrix elements. This process is a laborious task and can not be incorporated in a computational software easily.
Unlike previous methods, present computational scheme needs significantly less effort and is straightforward in converting to computer software language. The starting point of the model is constructing appropriate many-body basis needed for solving the Hamiltonian.
All of the possible distribution of electrons in the four adjacent atomic orbital make a set of complete basis in the configuration space. In contrast to previous models that have used molecular orbital for such a system we have developed atomic orbital bases by taking into account $`z`$ component of the spin degree of freedom with occupation condition of each orbital. Following notation can be used to represent each configuration state. Ignoring the spin degree of freedom simply reduces the states to the occupation number or the Fock space representation. Fock space representation is commonly used in modern molecular calculations.39

$$|\mathrm{\Psi }_i>=|a_{i0},a_{i1},a_{i2},\mathrm{},a_{i8}>$$
(4)
Parameters $`a_{i1}`$ to $`a_{i8}`$ are -1 or +1 (for spin down and up on the site $`i`$ respectively) which show occupation condition of each state. In fact these coefficients not only show extension of electronic wave function on each neighboring site of vacancy but also include spin information of the system. $`a_{i0}`$ is the sum of all $`a_i`$’s i.e. the total spin of each $`\mathrm{\Psi }_i`$ along $`z`$ axis. We will show the importance of this quantum number in our computation scheme later. Other capability of this notation is estimation of the maximum dimension of configuration space which one needs for calculating the Hamiltonian for $`V^0`$ $`(V^{})`$. In the $`V^0`$ $`(V^{})`$ system there are four (five) electrons respectively in four orbital, so the possible configurations of these electrons in the $`8`$ accessible states ($`2`$ spin and $`4`$ space or site degree of freedom) according to statistical combination rules are:

$$C_8^4=70(V^0)$$
(5)
$$C_8^5=56(V^{})$$
(6)
The fractions of $`70`$ $`(56)`$ configuration states in $`V^0`$ $`(V^{})`$ system which belong to each $`S_z`$ block are summarized in Table I. This counting method is in contrast to the old molecular orbital approaches which one needs to allow construction of symmetry and spin adapted basis by hand and then counting them.
After constructing the set of complete basis, we will calculate two form of Hamiltonian in this space. The first one is approximate form of Eq. (1), i.e. a simple extension of Hubbard model and the second is the exact form of the Hamiltonian (Eq. (1)) which we call it generalized Hubbard Hamiltonian.

## III Generalized Hubbard Hamiltonian

The main goal of the Hubbard Hamiltonian40 is resuming atomistic nature of the solid besides the free electron gas theory. This model assumes that the most important part of $`ee`$ interaction terms is on site term on the Columbic parts of Hamiltonian. The conventional form of this Hamiltonian consists of two parts, hopping term or $`t`$ term and on site term or $`U`$ term which is $`V_{iiii}`$ term of Eq. (3).

$$H=\underset{ij\sigma }{}t_{ij}c_{i\sigma }^{}c_{j\sigma }+\underset{i}{}U_in_in_i$$
(7)
Famous Hubbard model with two parameters $`t`$ and $`U`$ is not sufficient for solving the vacancy problem, the results of our model with these two parameters are not satisfactory. This result is physically expecting since assuming only two parameters of the Hubbard model for this system means that the electrons interact with each other only when they are at the same site of at the same atom of the vacancy and this eliminates the interaction of the electrons when they are in the different atomic site of the vacancy. Therefore we started with extended Hubbard Hamiltonian with one extra parameter $`V`$ which is $`V_{ijij}`$ term of Eq. (3). This Hamiltonian is equivalent to the Lannoo model32 where the $`U`$ parameter of their model is called $`V`$ in our model. This means that the on site and two site direct overlap integrals are the most dominant ones and the rest parameters are negligible (Eq.(8)).

$$H=\underset{ij\sigma }{}t_{ij}c_{i\sigma }^{}c_{j\sigma }+\underset{i}{}U_in_in_i+\frac{1}{2}\underset{ij\sigma \sigma ^\mathrm{`}}{}V_{ij}n_{i\sigma }n_{j\sigma ^\mathrm{`}}$$
(8)
Since the atomic sites in vacancy are from the vertex of a symmetrical tetrahedral, the distance between each two of them is the same, and hopping $`(t_{ij})`$, on site $`(U_i)`$ and two site $`(V_{ij})`$ parameters of Hamiltonian are all independent of $`i,j`$ indices and can be put outside of the sum Eq. (8).

$$H=t\underset{ij\sigma }{}c_{i\sigma }^{}c_{j\sigma }+U\underset{i}{}n_in_i+\frac{1}{2}V\underset{ij\sigma \sigma ^\mathrm{`}}{}n_{i\sigma }n_{j\sigma ^\mathrm{`}}$$
(9)
As discussed in sec. II, for $`V^0`$ $`(V^{})`$ system the size of the set which includes complete basis for solving this Hamiltonian is $`70`$ $`(56)`$. As a result solving Hamiltonian in this set results in a $`70\times 70`$ ($`56\times 56`$) Hamiltonian matrix for $`V^0`$ ($`V^{}`$) respectively. As the starting point we used numerical values of Hamiltonian parameters $`t`$, $`U`$ and $`V`$ which were reported by semi-empirical and theoretical works on diamond vacancies.33,41
This formalism can investigate the effect of variation of each parameter on the energy spectrum directly. By this means we found that variation of parameter $`t`$ has not any effect on the $`GR1`$ transition energy. But variation of parameters $`U`$ and $`V`$ changes the transition energy of $`GR1`$. These variations do not affect the sequence of the levels. The results of calculation for $`V^0`$ show that for a wide variation range of variables, $`t`$ greater than $`12.5`$ eV, $`U`$ greater than $`8.2`$ eV and $`V`$ greater than $`4.0`$ eV, the ground state of the system has $`{}_{}{}^{1}E`$ symmetry (double degenerate spin-less) and we can obtain appropriate transition to an excited state $`{}_{}{}^{1}T_{2}^{}`$ (triple degenerate spin-less). With our model we were able to obtain the experimental value of $`1.673`$ eV for the GR1 transition.
For $`V^{}`$ system, variation of $`t`$ has a dramatic effect on $`ND1`$ transition energy, but does not alter the ground and excited states sequence. Unlike this, $`U`$ and $`V`$ parameters which have only small effect on the values of the transition energies of the system. For the same range of variation of parameters we had a transition from $`{}_{}{}^{4}A_{2}^{}`$ ground state to a $`{}_{}{}^{4}T_{1}^{}`$ excited state ($`ND1`$ transition) but By the same parameters which we obtained the transition energy of $`GR1`$ ($`1.673`$ eV), there is a wide gap between the ground and excited states of $`ND1`$ transition (about $`30`$ eV).
Only with parameter $`t`$ close to $`1`$ was the result close to the experimental value. This high value for $`t`$ seems unreasonable for the vacancy problem and is very far from reported range.33,41
These results of the extended Hubbard Hamiltonian are in agreement with the results of the Lannoo and Bourgoin 32 model which has simplified the full Hamiltonian in the first quantization form by four parameters. They have used symmetrical molecular orbital base and also configuration interaction similar to Coulson and Kearsly model. The benefith of the reduction of the number of Hamiltonian parameter in their model was in this point that the validity range of single and many particle approches could be investigated. They investigated the results of their model in a full range of parameters and also used the semi-empirical parameters of the Coulson and Kearsly model. They successfully explained where the many body model works and where the single particle approaches can be valid. Lowther21 applied a similar model to the Lannoo model to explain some related optical phenomenon of the vacancies in diamond. Mainwood and Stonham26 applied Lannoo model directly to the vacancies electronic states in diamond. It seems that lack of the quantitative agreement with observed $`ND1`$ transition energy ($`3.15`$ eV) from $`{}_{}{}^{4}A_{2}^{}`$ ground state to a $`{}_{}{}^{4}T_{1}^{}`$ excited state in $`V^{}`$ system is due to elimination of exchange terms in Hamiltonian of vacancy system. Since calculation results of $`V^{}`$ system show weaker agreement with the experiment than those for $`V^0`$, we can deduce that the role of the exchange terms in charged electron vacancy which has higher electron density is more important than that the four electron (neutral) one. Quantitative failure of the results of this simple extension of Hubbard Hamiltonian reveals that for vacancy system $`ee`$ correlation effects are important and can not be neglected in contrast to the other cases which we are able to apply the Hubbard model. This point also was emphasized by previous models which have attempted diamond vacancy problem.21-26,33,41 Comparing numerical values of overlap integrals of Hubbard model40 and diamond vacancy problem33,41 (Table II) we find that the values of exchange integrals in this problem are comparable to the values of direct integrals so their elimination are not reasonable.

This leads us to keep all of the, $`t`$, $`U`$ and $`V`$ terms besides exchange terms in four (neutral) and five (charged) electronic Hamiltonian of the system similar to Coulson and Kearsly model.33,41 We will call it generalized Hubbard Hamiltonian (Eq. (1)). Although this significantly increases the volume and complexity of the calculation, but since we take into account exact $`ee`$ correlation effects, it makes our model more realistic.
We return to Hamiltonian form of Eq. (1). Similar symmetrical argument is valid for a system with $`T_d`$ symmetry and we can put hopping parameter $`(t)`$ out of the sum (Eq. (10)).

$$H=t\underset{ij\sigma }{}c_{i\sigma }^{}c_{j\sigma }+\frac{1}{2}\underset{ijlm\sigma \sigma ^\mathrm{`}}{}V_{ijlm}c_{i\sigma }^{}c_{j\sigma ^\mathrm{`}}^{}c_{m\sigma ^\mathrm{`}}c_{l\sigma }$$
(10)
Symmetry of defect molecule, reduces the number of the independent parameters $`V_{ijlm}`$, from $`4^4`$ $`(256)`$ components to only $`7`$ independent ones, hence symmetry considerations can be included in the generalized form of Hamiltonian with a total of $`8`$ parameters. These consist of the hopping $`t`$, two direct Coulombic integrals $`(U,V)`$ and five exchange $`(X1,\mathrm{},X5)`$ integrals as follows:

$`U=ii|{\displaystyle \frac{1}{r}}|ii`$
$`V=ij|{\displaystyle \frac{1}{r}}|ij`$
$`X1=ij|{\displaystyle \frac{1}{r}}|ji`$
$`X2=ii|{\displaystyle \frac{1}{r}}|ij`$
$`X3=ij|{\displaystyle \frac{1}{r}}|ik`$
$`X4=ii|{\displaystyle \frac{1}{r}}|jk`$
$`X5=ij|{\displaystyle \frac{1}{r}}|kl`$
Since the Hamiltonian is spin independent, useful conserving quantum numbers

$$[H,S^2]=[H,S_z]=0$$
(11)
are total spin of four electron system $`(S^2)`$ and the $`z`$ component of total spin $`(S_z)`$. The important computational point about Eq. (11) is that it should hold for each arbitrary set of parameters of Eq. (10).
By these conservation rules we can easily decrease dimension of the nonzero block of Hamiltonian matrix by using conservation of total spin ($`S^2`$) and $`z`$ component of the total spin ($`S_z`$). Therefore the volume of calculation significantly reduces and matrix element evaluation only confines in subspaces where wave functions have the same values of $`S_z`$. As a result, the maximum dimension of nonzero matrix block in the case of $`V^0`$ ($`V^{}`$) systems becomes $`36`$ $`(24)`$ respectively. The final step in solving the Hamiltonian of the $`V^0`$ ($`V^{}`$) for obtaining the resultant eigenstates with definite $`S_z`$ and $`S^2`$ values is the transformation of Hamiltonian matrix from $`S_z`$ to ($`S^2,S_z`$) basis.
We have used the parameters which were calculated by Coulson and Larkins41 to investigate energy spectrum and eigenstates of the system. These parameters have been calculated using Slater and Clementi type functions as atomic orbitals.
The parameter calculation is independent from the Hamiltonian solution formalism. Since the parameters are calculated two times by the Coulson 33,41 during a long period, other models similar to us have not attempted to re-estimate their values and all of them have used the results of the Coulson estimation for the parameters. 21,22,26,32The calculation of parameters were based on symmetric and antisymmetric molecular orbital basis ($`a,t_x,t_y,t_z`$)which were combinations of $`a,b,c,d`$ atomic orbitals. These atomic orbitals belong to the four nearest neighbor atoms of vacant site and are oriented inward to its center, however in the generalized Hubbrad Hamiltonian Eq. (1), parameters are calculated directly by using overlap integrals of these $`a,b,c,d`$ atomic orbitals. These overlap integrals can be extracted from the older symmetric and anti-symmetric ones by changing the calculation basis. Related calculations have been performed and the resultant parameters are used in the computational scheme.

## IV Results and discussion

The two key parameters of this model are $`t`$ and $`U`$. These two have been evaluated semi-empirically in the previous models and were starting point in calculation of the other parameters.33,41 The reported values for hopping parameter $`t`$ cover a wide range of variations, for example both $`7.13`$ eV and $`16.34`$ eV values at the same time are reported for this parameter.33,41 Concerning parameter $`U`$, the semi-empirical value is estimated to be $`13.29`$ eV and after revising $`12.85`$ eV while the theoretically calculated value is $`19`$ eV .33,41

In evaluation of the $`t`$ parameter despite of other $`7`$ Coulombic parameters, one not only need the information of the spatial distribution of the wave functions but also the form of the psedupotential of carbon atoms is needed, so the uncertainty in the value of this is more than other parameters. In the present model the variation of the energy spectrum with $`t`$ has been investigated and results are shown in Fig. 1. The values of the other parameters are according to table III. According to this figure, for obtaining correct results for ground state of $`V^0`$ i.e. $`{}_{}{}^{1}E`$, parameter $`t`$ should be confined in the following range.

$$6.25eVt10.22eV$$
(12)
Choosing $`t`$ out of this range converts the sequence of two lowest states i.e. $`{}_{}{}^{1}E`$ and $`{}_{}{}^{3}T_{1}^{}`$, and gives wrong ground state ($`{}_{}{}^{3}T_{1}^{}`$) for the system. Increasing value of the $`t`$ will lower the experimentally observed spin quintet $`{}_{}{}^{5}A_{2}^{}`$ state.11 It is interesting to note that above the upper limit of $`t`$, $`{}_{}{}^{5}A_{2}^{}`$ competes with $`{}_{}{}^{3}T_{1}^{}`$ to be the ground state and at $`t`$ value equals to $`5.96`$ eV, $`{}_{}{}^{5}A_{2}^{}`$ becomes ground state of $`V^0`$ system. It seems that by increasing $`t`$ we reach an area where the Hund’s rule can be applied (according to this rule $`{}_{}{}^{5}A_{2}^{}`$ must be the ground state). This might be due to the point that $`t`$ is the indicator of the looseness of the electrons from the nuclei so by increasing $`t`$ we can consider the electrons as more delocalized electrons.
For charged vacancy case, the calculation results are shown in Fig.2. From this

figure we conclude that for obtaining $`{}_{}{}^{4}A_{2}^{}`$ as ground state for the system, parameter $`t`$ should be restricted as:

$$6.95eVt16.34eV$$
(13)
The lower limit is the negative of the ionization energy of electron in C-C bound of diamond. From this figure we can find that increasing $`t`$ over upper limit converts the ground state to $`{}_{}{}^{4}T_{1}^{}`$ which disagrees with the experiment.
By comparing the results of the calculation for both $`V^0`$ and $`V^{}`$ we can put new boundaries to the single particle hopping parameter $`t`$ as follows:

$$6.95eVt10.22eV$$
(14)
This restriction on $`t`$ parameter, which is the most uncertain parameter in the vacancy electronic states problem, is tighter than already reported range ( $`7.13`$ eV to $`16.14`$ eV).33,41 Variation of $`t`$ in this allowed range does not have significant effect on $`GR1`$ transition energy. This is in agreement with other theoretical works.33,41 However this variation has a dramatic effect on $`ND1`$ transition energy.
Numerical optimization of parameter $`t`$ in this range and also parameter $`U`$ gives $`7.747`$ eV and $`13.58`$ eV for $`t`$ and $`U`$ respectively. These values are very close to the semi-empirically derived values by Coulson and Larkins (the difference between our optimized value by semi-empirical values in this two parameters are less than 8 and 5 percent respectively). 41 Other six parameters are as the same as the theoretically calculated values by them using Slater type functions 41. Calculation methods of these parameters are completely independent from the Hamiltonian solution formalism so the previous models 21,22,26,32 similar to us have used the parameters which are reported semi-empirically and theoretically by Coulson and Kearsly 33and Coulson and Larkins.41 These parameters are calculated by them during a 17 year period and with the Slater type functions so it seems that the parameters are well justified. The difference between our suggested parameters with two semi-empirical parameters of their model is very less than the wide range of variation of these parameters in their model.( parameter t from $`7.13`$ eV to $`16.34`$ eV and parameter U which smi-empirically is $`12.855`$ eV and theoretically is $`19`$ eV)
With these justified values the present model can simultaneously obtain the exact transition energies of the $`GR1`$ ($`1.67`$ eV) and $`ND1`$ ($`3.15`$ eV) for neutral and charged vacancy system.
The set of our parameters and parameters that were previously computed based on Slater and Clementi type functions33,41,42 are listed in Table (III). The model predicts the wrong ground state with parameters which are based on Clementi type functions, i.e. $`{}_{}{}^{3}T_{1}^{}`$ instead of the experimentally observed $`{}_{}{}^{1}E`$ state for $`V^0`$. This is similar to the Coulson and Larkins’s report. 41 However parameters which are based on the Slater type functions predict correct spin and symmetry for ground and excited states of $`GR1`$ and $`ND1`$ absorption lines but the values of the transition energies are not satisfactory. The results of the present model calculation are summarized in Fig. 3,4.
For neutral vacancy as it is shown in Fig. 3, there is a transition from double degenerate spinless ground state to a triple

degenerate spinless excited state $`(^1E^1T_2)`$ with a transition energy equal to $`1.67`$ eV.
We also see that there is another level $`{}_{}{}^{3}T_{1}^{}`$, very close to $`{}_{}{}^{1}E`$ ground state with a very little energy difference i.e. about $`113`$ meV above it. This value agrees with the recent experimental expectations which predict it will be more than $`100`$ meV .43 However the results of Coulson models gives $`40`$ meV which disagree with the experiment.33,41 This very close energy level to ground state also has been used to explain the behavior of $`R2`$ center absorption line in diamond where it was attributed to the strongly perturbed vacancy.42 By this assumption the energy difference of this level with the ground state has been predicted to be between $`40`$ meV and $`200`$ meV.42 Mainwood and Stonham 26 have used the simple model of Lannoo 31 and successfully explain many physical features of the vacancy. Lannoo for simplifying the problem do not enters the exchange parameters into energy levels and use only two parameters for describing energy levels. This is unlike to Coulson 32,40 and our model which reserve whole of the e-e terms. Besides explaining many physical feature of the vacancy they gave some estimation for the position of this level. Using the model of Lannoo they positioned $`{}_{}{}^{3}T_{1}^{}`$ at $`200`$ meV above the $`{}_{}{}^{1}E`$, Although this value is as the same magnitude as the error in the positions of the levels in their calculations.26
Our reported location for $`{}_{}{}^{3}T_{1}^{}`$ arises in a natural manner similar to $`GR1`$ and $`ND1`$ after evaluation of the exact Hamiltonian of the system with the justified parameters. The error which can be attributed to this value in our calculation is the difference between our parameters and semi-empirical parameters of the Coulson. This error is less than $`8`$ percent which is significantly less than already reported error for the position of this level in similar calculations.26 In obtaining this value we have not assumed that the vacancy is highly perturbed as it already was reported.42 This value can help EPR experimental investigations that expect to observe this state in temperatures higher than $`300`$ K.43

In the charged vacancy system, as it is shown in Fig. 4, the model predicts a transition from a non-degenerate ground state to a triple degenerate excited state with a spin equal to $`3/2`$ $`(^4A_2^4T_1)`$. The energy difference is $`3.15`$ eV which is the famous observed $`ND1`$ absorption line. Previous models that evaluate e-e interaction exactly can not obtain this value by any set of parameters. 33,41
In summary as it has been mentioned by other authors 32,26 the previous models that evaluate e-e interaction exactly 33,41 can only explain GR1 transition energy quantitatively. They could not predict $`ND1`$ transition energy and also they have predicted the experimentally wrong energy for the low-lying $`{}_{}{}^{3}T_{1}^{}`$ level. These unsatisfactory quantitative results are obtained not only by the semi-empirical parameters but also in a wide variation range of these two semi-empirical parameters.33,41 This variation ranges is very wider than the difference of our optimized values with semi-empirically reported ones for $`t`$ and $`U`$. Other six Coulombic parameters in our model and the previous ones are evaluated by theory.33,41 Our model can simultaneously give accurate transition energies for both of the $`GR1`$ and $`ND1`$ transitions besides the experimentally acceptable energy for the low-lying $`{}_{}{}^{3}T_{1}^{}`$ level with almost the same parameters ( only two semi-empirical values differs less than 8 percent).
The results of our model have justified the Coulson and Kearsly semi-empirical parameters and also their guess33 which expected that the $`t=8`$ eV value be more realistic than $`t=16.34`$ eV . Both values of this parameter have been used in their models.33,41
All of these qualitative and quantitatively satisfactory results in our model have obtained without configuration interaction. The good quantitative results for $`GR1`$ and even other reasonable qualitative results in the Coulson and Kearsly model and all of the previous molecular approaches 22-26,32 can only be achieved with assistance of the Configuration interaction technique. This is for the first time that a molecular approach put aside the configuration interaction technique and could predict more experimental data. Configuration interaction has been widely used in all of the previous theoretical works to overcome the disagreement between theoretical results and experiment even at the qualitative levels.22-26,33,41,42 Although precise review of the previous exact models 33,41 results on the charged vacancy with parameter $`t=16.34`$ eV shows that applying configuration interaction gives wrong ground state for $`V^{}`$.
Similar ambiguity also exists in the electronic configuration of the ground and excited state of the GR1. While the models assume that the ground and excited states of $`GR1`$ arise from the single configuration $`a^2t^2`$, at the same time they try to find or use some semi-empirically coefficients for the contribution of the other possible configurations in the ground and excited states of $`GR1`$.21,26 Here we refer to Davies 4 remark in this regard:
”The importance of configuration interaction at the vacancy in diamond was controversial, and alternative approaches were investigated in which emphasize was on the benefit of using large cluster”.
The main advantages of the presented formalism which is based on molecular model is removing the configuration interaction besides it’s satisfactory quantitative results. By this means our model have solved one important shortcoming of the molecular approaches with respect to cluster models or density functional approaches.
In the present model, the resultant eigenstates of the Hamiltonian have unique expansion on the starting constructed basis in configuration space, so configuration interaction comes out in a natural manner from the formalism. From this point, it is possible to find the contribution of each electronic configuration in ground and excited states of the system.
For $`GR1`$ transition, the results of calculation for the ground state $`{}_{}{}^{1}E`$ and the excited state $`{}_{}{}^{1}T_{2}^{}`$ with the justified parameters are illustrated in Fig. 5. The possible electronic configuration for states with value of $`S=0`$ are: $`(1,1,1,1),(2,1,1,0),(2,2,0,0)`$ where the numbers in parenthesis show the occupation number of each tetrahedral orbital of vacancy system. As it is shown in Fig. 5,

all possible electronic configurations are present in the ground state $`{}_{}{}^{1}E`$. The most probable configuration in ground state arise from $`(1,1,1,1)`$ configuration ($`61`$ $`percent`$) which has only singly occupied orbitals. Two doubly occupied configurations have very low contribution (3 percent) in the ground state. In $`(1,1,1,1)`$ configuration the electrons have maximum separation distance, which results in lowering the Coulombic repulsion energy of the system. In contrast to this, the excited state of $`GR1`$, i.e. $`{}_{}{}^{1}T_{2}^{}`$ is orthogonal to $`(1,1,1,1)`$ configuration and is almost a purely doubly occupied configuration state. This configuration has the highest Coulombic repulsion energy due to its lowest separation distance between electrons of the system.
Roughly speaking, we can state that $`GR1`$ transition is resultant of a transition from higher probability configuration $`(1,1,1,1)`$ and $`(2,1,1,0)`$ configurations of the system to $`(2,1,1,0)`$ configuration .
Our results modifies the common belief about ground and excited states of the $`GR1`$ transition that assumes they both belong to $`a^2t^2`$ configuration.21,26,33,41,42,44,45
As the results of our calculation show, the absence of $`(1,1,1,1)`$ configuration in the ground state of the system is unexpected, since it has minimum Coulombic repulsion energy among the other configurations of the system. Also it is the only compatible configuration of the unique $`{}_{}{}^{5}A_{2}^{}`$ excited state, which in the delocalization limit ($`t>5.96`$ eV) is the ground state of the system.
For $`ND1`$ transition, since the ground and exited states have the spin equal to $`3/2`$ the only configuration which can be encountered is $`(2,1,1,1)`$. Therefore the ground and excited states are only a combination of these configurations with different $`S_z`$ values and the results are same as other models .
The new information about the ground and excited state of the $`GR1`$ is very important, It can open a new opportunity for ab-initio density functional approaches which have not been attempted to estimate the value of this transition.35 This is due to the results of the previous molecular orbital models which have assumed the ground and excited state belong to the same configuration $`a^2t^2`$. However by the new result of present model, they can try to obtain the transition energy of the $`GR1`$ transition similar to $`ND1`$. These information are also useful for theoretical approaches which start from a guess wave function for ground and excited states of such a molecular system.

## V Conclusion

Theoretical studies of the vacancy electronic states in diamond have been a challenging problem for the five decades. Localized models which originate from Coulson and Kearsly pioneering work, predict correctly the spatial symmetry and spin of electronic states but can not give quantitatively good results for electronic transitions in diamond vacancies. In their model despite of excellent theoretical framework and strong physical intuition, the computational method are very old and can not be adapted with nowadays powerful computer assisted methods.
We have developed a new approach based on generalized Hubbard Hamiltonian and many body techniques. A complete set of atomic orbital basis for describing this system is constructed according to $`S_z`$ representation of each electronic configuration in the four vacancy orbital. This confines the calculation only in definite $`S_z`$ subspace. This representation can simply estimate number of states of the system and also each definite $`S_z`$ block size of Hamiltonian. This is in contrast to previous models that uses molecular orbital bases and manually construct each symmetry and spin adapted wave functions.
The generalized Hubbard Hamiltonian which contains whole of the $`ee`$ correlation terms is included with spatial symmetry information of the vacancy therefore the number of the independent parameters of the Hamiltonian is reduced to $`8`$. By this means the effect of the parameters on the energy states directly can be examined. A detailed computer software is developed to solve this Hamiltonian . The resultant Hamiltonian matrix is transformed to change basis from $`S_z`$ to $`S^2`$ representation. As the result, the eigenstates have definite and correct spin and symmetry properties.
The effect of changing each Hamiltonian parameters on the electronic states of the system directly can be examined. This calculation shows that for obtaining experimentally correct results for the ground state of $`V^0`$ and $`V^{}`$, parameter $`t`$ should be in the $`10.22`$ eV to $`6.95`$ eV range, although increasing $`t`$ in this range has low effect on $`GR1`$ but significantly increases $`ND1`$ transition energy. According to this model the optimum value for two key semi-empirical parameters $`t`$ (in allowed range) and $`U`$ are $`7.747`$ eV, $`13.58`$ eV respectively. Two famous absorption lines of $`V^0`$ and $`V^{}`$ namely $`GR1`$ ($`1.67`$ eV) and $`ND1`$ ($`3.15`$ eV) can be obtained simultaneously and exactly with these parameters in addition to the location of the low-lying $`{}_{}{}^{3}T_{1}^{}`$ state at $`113`$ meV above the ground state of the $`V^0`$. These results have been obtained for the time by such justifiable parameters. These parameters also justify already reported semi-empirical values that only could estimate $`GR1`$ transition and gave experimentally wrong results for $`{}_{}{}^{3}T_{1}^{}`$.
For the first time the configuration interaction procedure, has not used in such a calculations and it arises in a natural from from the formalism. Despite of its ambiguous physical nature, the configuration interaction techniques plays a fundamental and vital role in reaching to an agreement with experimental data in previous localized models. This vital role even exist at the level of qualitatively acceptable results. This point enables us to estimate the contribution of each electronic configuration in ground and excited states of the system. The results for $`GR1`$ absorption line show that it is almost a transition from configurations that consist of a pure singly occupied orbital and a singly double occupied orbital to a singly occupied orbital configuration. This new information can be used by the density functional approache to this problem. The results of our model justifies the first suggested molecular model by the Coulson and Kearsly is the five decades ago and shows that their model has enough physical depth to accompany with a new computational methods to obtain very satisfactory quantitative results. The results also show the benefit of the localized models with respect to cluster or density functional approaches by no need of configuration interaction .

## VI Acknowledgment

We would like to thanks Dr.A. Shafikhani and Mr.M. Saiari for their valuable comments about this work. We acknowledge the help of Mr.A. Rezakhani(Phd student) and Mr.M. Ghasempour for computer programming and paper preparation.
e-mail: Heydaris@mehr.sharif.edu
e-mail: Vesaghi@sharif.edu
e-mail: K1@sharif.edu


format: text

text:
# 1 Introduction

## 1 Introduction

Spurred on by the famous EMC experiment  at CERN in 1987, there has been a huge growth of interest in polarized DIS experiments which yield more refined information about the partonic structure of the nucleon, i.e., how the nucleon spin is divided up among its constituents, quarks and gluons. Many experiments have been carried out at SLAC, CERN and DESY to measure the longitudinal ($`A_{}`$) and transverse ($`A_{}`$) asymmetries and to extract from them the photon-nucleon asymmetries $`A_1(x,Q^2)`$ and $`A_2(x,Q^2)`$ as well as the nucleon spin-dependent structure functions $`g_1(x,Q^2)`$ and $`g_2(x,Q^2)`$. Many theoretical analyses of the world data on $`A_1`$ and $`g_1`$ based on leading (LO) and next-to-leading order (NLO) calculations in perturbative QCD have been performed in order to test the spin properties of QCD and extract from the data the polarized parton densities. Note that the theoretical analyses have been mainly concentrated on the $`A_1`$($`g_1`$) data because the measurements of the quantities $`A_2`$($`g_2`$) are much less accurate with the exception of the very recent data of E155 Collaboration at SLAC . Another reason is that the theoretical treatment of $`g_2`$ is much more complicated. It was demonstrated that the polarized DIS data are in an excellent agreement with the pQCD predictions for $`A_1^N(x,Q^2)`$ and $`g_1^N(x,Q^2)`$. What also follows from these analyses is that the limited kinematic range and the precision of the present generation of inclusive DIS experiments are enough to determine with a good accuracy only the polarized parton densities $`(\mathrm{\Delta }u+\mathrm{\Delta }\overline{u})(x,Q^2)`$ and $`(\mathrm{\Delta }d+\mathrm{\Delta }\overline{d})(x,Q^2)`$. The polarized strange sea density $`(\mathrm{\Delta }s+\mathrm{\Delta }\overline{s})(x,Q^2)`$ as well as the polarized gluon density $`\mathrm{\Delta }G(x,Q^2)`$ are still weakly constrained, especially $`\mathrm{\Delta }G`$. The non-strange polarized sea-quark densities $`\mathrm{\Delta }\overline{u}`$ and $`\mathrm{\Delta }\overline{d}`$ cannot be determined, even in principle, from the inclusive DIS experiments alone without additional assumptions.

There is, however, an important difference between the kinematic regions of the unpolarized and polarized data sets. While in the unpolarized case we can cut the low $`Q^2`$ and $`W^2`$ data in order to eliminate the less known non-perturbative higher twist effects, it is impossible to perform such a procedure for the present data on the spin-dependent structure functions without loosing too much information. This is especially the case for the HERMES, SLAC and Jefferson Lab experiments. So, to extract the polarized parton densities from the experimental data on $`g_1^N(x,Q^2)`$ the higher twist corrections have to be included in the data fits. Note that the polarized parton densities in QCD are related to the leading-twist expression of $`g_1`$.

It was shown  that to avoid this problem and to determine polarized parton densities less sensitive to higher twist effects it is better to analyze data on $`A_1(g_1/F_1)`$ using for the $`g_1`$ and $`F_1`$ structure functions their leading twist (LT) expressions. It is found that if for $`(g_1)_{\mathrm{LT}}`$ an NLO approximation is used, the ”effective higher twist” corrections to $`A_1`$, extracted from the data, are negligible and consistent with zero within the errors, which means that the higher twists corrections (HT) to $`g_1`$ and $`F_1`$ approximately cancel in the ratio $`g_1/F_1`$, or more precisely, $`(g_1)_{\mathrm{HT}}/(g_1)_{\mathrm{LT}}(F_1)_{\mathrm{HT}}/(F_1)_{\mathrm{LT}}`$.

In this paper we present a detailed study of the higher twist contributions $`h^N(x)/Q^2`$ to the nucleon structure function $`g_1^N(x,Q^2)`$. The quantities $`h^N(x)`$ have been extracted from the data in a model independent way. The role of higher twists in the determination of the polarized parton densities is discussed.

## 2 Connection between Theory and Experiment

The nucleon spin-dependent structure function $`g_1^N(x,Q^2)`$ is a linear combination of the asymmetries $`A_{}^N`$ and $`A_{}^N`$ (or the related virtual photon-nucleon asymmetries $`A_1^N`$ and $`A_2^N`$) measured with the target polarized longitudinally or perpendicular to the lepton beam, respectively. The most direct way to confront the QCD predictions to the data is a fit to data on the ratio of the structure functions, $`g_1^N/F_1^N`$. Such data have been directly presented by SLAC/E143 and SLAC/E155 experiments . Most of the Collaborations, however, have presented data on the asymmetry $`A_1^N`$ which, in practice, are data on $`A_{}^N/D`$. The photon-nucleon asymmetry $`A_1^N`$ and the ratio $`g_1^N/F_1^N`$ are related to the measured longitudinal asymmetry $`A_{}^N`$ by

$`A_1^N`$ $`=`$ $`{\displaystyle \frac{A_{}^N}{D}}\eta A_2^N,`$ (1)
$`(1+\gamma ^2){\displaystyle \frac{g_1^N}{F_1^N}}`$ $`=`$ $`{\displaystyle \frac{A_{}^N}{D}}+(\gamma \eta )A_2^N,`$ (2)
where D denotes the photon depolarization factor, $`\eta `$ and $`\gamma `$ are kinematic factors. $`\eta `$ is proportional to $`\gamma `$ and $`\gamma `$ is given by

$$\gamma ^2=\frac{4M_N^2x^2}{Q^2}.$$
(3)
In (3) $`\mathrm{M}_\mathrm{N}`$ is the nucleon mass. It should be noted that in the SLAC and HERMES kinematic regions $`\gamma ^2`$ cannot be neglected on LHS of (2).

The magnitude of $`A_2^N`$ has been measured by SMC, SLAC/E143 and SLAC/E155 and found to be small . Then to a good approximation its contribution to the RHS of Eqs. (1) and (2) can be neglected and $`A_1^N`$ and $`g_1^N/F_1^N`$ can be expressed as

$`A_1^N`$ $``$ $`{\displaystyle \frac{A_{}^N}{D}},`$ (4)
$`(1+\gamma ^2){\displaystyle \frac{g_1^N}{F_1^N}}`$ $``$ $`{\displaystyle \frac{A_{}^N}{D}}.`$ (5)
It is important to note that due to the additional small factor $`(\gamma \eta )`$ in (2) the ratio $`g_1^N/F_1^N`$ is better approximated by the measured asymmetry $`A_{}^N`$ (Eq. (5)) than the virtual photon-nucleon asymmetry $`A_1^N`$ (Eq. (4)).

Using (4) and (5) we reach the well known relation

$$A_1^N(x,Q^2)(1+\gamma ^2)\frac{g_1^N(x,Q^2)}{F_1^N(x,Q^2)}$$
(6)
usually used in the literature. However, as was already mentioned, we have to keep in mind that the presented experimental values on $`A_1^N`$  neglecting $`A_2`$, $`{}_{}{}^{\prime \prime }A_{1}^{N}(x,Q^2)_{exp}^{\prime \prime }`$, are really the experimental values of $`A_{}^N(x,Q^2)/D`$ and that the latter quantity is very well approximated by $`(1+\gamma ^2)g_1^N/F_1^N`$.

Using the relation between the unpolarized structure function $`F_1(x,Q^2)`$ and the usually extracted from unpolarized DIS experiments $`F_2(x,Q^2)`$ and $`R(x,Q^2)`$

$$2xF_1^N=F_2^N(1+\gamma ^2)/(1+R^N)(N=p,n,d)$$
(7)
Eq. (6) can be rewritten as

$$A_1^N(x,Q^2)\frac{g_1^N(x,Q^2)}{F_2^N(x,Q^2)}2x[1+R^N(x,Q^2)].$$
(8)
Up to now, two approaches have been mainly used to extract the polarized parton densities (PPD) from the world polarized DIS data. According to the first  the leading twist LO/NLO QCD expressions for the structure functions $`g_1^N`$ and $`F_1^N`$ have been used in (6) in order to confront the data

$`\left[{\displaystyle \frac{g_1(x,Q^2)}{F_1(x,Q^2)}}\right]_{exp}`$ $``$ $`{\displaystyle \frac{g_1(x,Q^2)_{\mathrm{LT}}}{F_1(x,Q^2)_{\mathrm{LT}}}},`$
$`A_1(x,Q^2)_{exp}`$ $``$ $`(1+\gamma ^2){\displaystyle \frac{g_1(x,Q^2)_{\mathrm{LT}}}{F_1(x,Q^2)_{\mathrm{LT}}}}.`$ (9)
In (9) we have dropped the nucleon target label $`N`$. It was shown  that in this case the extracted from the data “effective” HT corrections $`h^{A_1}(x)`$ to $`A_1`$

$$A_1(x,Q^2)=(1+\gamma ^2)\frac{g_1(x,Q^2)_{\mathrm{LT}}}{F_1(x,Q^2)_{\mathrm{LT}}}+\frac{h^{A_1}(x)}{Q^2}$$
(10)
are negligible and consistent with zero within the errors, $`h^{A_1}(x)0`$ (see Fig. 1). What follows from this result is that the higher twist corrections to $`g_1`$ and $`F_1`$ compensate each other in the ratio $`g_1/F_1`$ and the PPD extracted this way are less sensitive to higher twist effects. We stress again that the polarized parton densities in QCD are related only to the leading-twist part of $`g_1`$.

According to the second approach , $`g_1/F_1`$ and $`A_1`$ data have been fitted using phenomenological parametrizations of the experimental data for $`F_2(x,Q^2)`$ and $`R(x,Q^2)`$

$`\left[{\displaystyle \frac{g_1(x,Q^2)}{F_1(x,Q^2)}}\right]_{exp}`$ $``$ $`{\displaystyle \frac{g_1(x,Q^2)_{\mathrm{LT}}}{F_2(x,Q^2)_{exp}}}2x{\displaystyle \frac{[1+R(x,Q^2)_{exp}]}{(1+\gamma ^2)}},`$
$`A_1(x,Q^2)_{exp}`$ $``$ $`{\displaystyle \frac{g_1(x,Q^2)_{\mathrm{LT}}}{F_2(x,Q^2)_{exp}}}2x[1+R(x,Q^2)_{exp}].`$ (11)
Note that such a procedure is equivalent to a fit to $`(g_1)_{exp}`$, but it is in principle better than the fit to the $`g_1`$ data themselves actually presented by the experimental groups. The point is that most of the experimental data on $`g_1`$ have been extracted from the $`A_1`$ and $`g_1/F_1`$ data using the additional assumption that the ratio $`g_1/F_1`$ does not depend on $`Q^2`$. Also, different experimental groups have used different parametrizations for $`F_2`$ and $`R`$.

In the analyses  a procedure has been used which is somehow a mixture between the two methods above, but bearing in mind the sensitivity of the results to higher twist effects it is analogous to the second one. In these fits the leading twist expression for $`F_2`$ instead of its experimental values has been used

$`\left[{\displaystyle \frac{g_1(x,Q^2)}{F_1(x,Q^2)}}\right]_{exp}`$ $``$ $`{\displaystyle \frac{g_1(x,Q^2)_{\mathrm{LT}}}{F_2(x,Q^2)_{\mathrm{LT}}}}2x{\displaystyle \frac{[1+R(x,Q^2)_{exp}]}{(1+\gamma ^2)}},`$
$`A_1(x,Q^2)_{exp}`$ $``$ $`{\displaystyle \frac{g_1(x,Q^2)_{\mathrm{LT}}}{F_2(x,Q^2)_{\mathrm{LT}}}}2x[1+R(x,Q^2)_{exp}].`$ (12)
It was shown by GRSV  that if the second approach (11) is applied to the data ($`F_2`$ and $`R`$ are taken from experiment) allowing at the same time “effective higher twist” contribution $`h^{A_1}(x)/Q^2`$ to the RHS of (11), $`h^{A_1}(x)`$ is found to be sizeable and important in the fit. In other words, bearing in mind that a lot of data on $`A_1`$ and $`g_1/F_1`$ are at small $`Q^2`$ special attention must be paid to higher twist corrections to the structure function $`g_1`$. To extract correctly the polarized parton densities from the $`g_1`$ data these corrections have to be included into data fits. Note that a QCD fit to the data in this case, keeping in $`g_1(x,Q^2)_{QCD}`$ only the leading-twist expression, leads to some ”effective” parton densities which involve in themselves the HT effects and therefore, are not quite correct.

## 3 Higher twist effects in $`g_1(x,Q^2)`$

The usual pQCD expression for the nucleon structure function $`g_1^p(x,Q^2)`$, in terms of polarized quark and gluon densities, arises from the contribution of the leading twist ($`\tau =2`$) QCD operators and in NLO has the form (a similar formula holds for $`g_1^n`$):

$$g_1^p(x,Q^2)_{\mathrm{pQCD}}=\frac{1}{2}\underset{q}{\overset{N_f}{}}e_q^2[(\mathrm{\Delta }q+\mathrm{\Delta }\overline{q})(1+\frac{\alpha _s(Q^2)}{2\pi }\delta C_q)+\frac{\alpha _s(Q^2)}{2\pi }\mathrm{\Delta }G\frac{\delta C_G}{N_f}],$$
(13)
where $`\mathrm{\Delta }q(x,Q^2),\mathrm{\Delta }\overline{q}(x,Q^2)`$ and $`\mathrm{\Delta }G(x,Q^2)`$ are quark, anti-quark and gluon polarized densities in the proton, which evolve in $`Q^2`$ according to the spin-dependent NLO DGLAP equations. $`\delta C(x)_{q,G}`$ are the NLO spin-dependent Wilson coefficient functions and the symbol $``$ denotes the usual convolution in Bjorken $`x`$ space. $`\mathrm{N}_\mathrm{f}`$ is the number of active flavors. In LO QCD the coefficients $`\delta C(x)_q`$ and $`\delta C(x)_G`$ vanish and the polarized parton densities in (13) evolve in $`Q^2`$ according to the spin-dependent LO DGLAP equations.

It is well known that at NLO and beyond, the parton densities as well as the Wilson coefficient functions become dependent on the renormalization (or factorization) scheme employed.Of course, physical quantities such as the virtual photon-nucleon asymmetry $`A_1(x,Q^2)`$ and the polarized structure function $`g_1(x,Q^2)`$ are independent of choice of the factorization convention. Two often used schemes are the $`\overline{\mathrm{MS}}`$ and the JET schemes. Both the NLO polarized coefficient functions  and the NLO polarized splitting functions (anomalous dimensions)  needed for the calculation of $`g_1(x,Q^2)`$ in the $`\overline{\mathrm{MS}}`$ scheme are well known at present. The corresponding expressions for these quantities in the JET scheme can be found in .

However, there are other contributions to $`g_1`$, arising from QCD operators of higher twist (HT), namely $`\tau 3`$, which are related to multi-parton correlations in the nucleon. It can be shown that these give rise to contributions to $`g_1^N(x,Q^2)`$ that decrease like inverse powers of $`Q^2`$. The leading term has the form $`h^N(x,Q^2)/Q^2`$, where $`h^N(x,Q^2)`$ could have a slow, logarithmic dependence on $`Q^2`$.

In the kinematic regime where such terms might be relevant it is important for consistency to realise that the QCD expression (13) is derived under the assumption that $`Q^2>>M_N^2`$. There will thus be purely kinematic corrections to (13), which involve a power series in $`M_N^2/Q^2`$ with small coefficients. The leading term of these so-called target mass corrections (TMC) therefore has a $`Q^2`$ behavior similar to the genuine HT terms, but it is not a dynamical HT effect.

In view of this we shall write

$$g_1(x,Q^2)=g_1(x,Q^2)_{\mathrm{LT}}+g_1(x,Q^2)_{\mathrm{HT}},$$
(14)
where we have dropped the nucleon target label N. In (14)

$$g_1(x,Q^2)_{\mathrm{LT}}=g_1(x,Q^2)_{\mathrm{pQCD}}+h^{\mathrm{TMC}}(x,Q^2)/Q^2,$$
(15)
where $`h^{\mathrm{TMC}}(x,Q^2)`$ is exactly calculable  and

$$g_1(x,Q^2)_{\mathrm{HT}}=h(x,Q^2)/Q^2.$$
(16)
As mentioned above, $`h(x,Q^2)`$ denotes the dynamical higher twist power corrections to $`g_1`$ which represent the multi-parton correlations in the target. The latter are non-perturbative effects and their calculation is model dependent (see, e.g.,  and references therein). That is why a model independent extraction of the dynamical higher twists $`h(x)`$ from the experimental data is important not only for a better determination of the polarized parton densities but also because it would lead to interesting tests of the non-perturbative QCD regime.

## 4 Method of Analysis

In this Section we will briefly describe the method of our analysis of the data on inclusive polarized DIS taking into account the higher twist corrections to the nucleon structure function $`g_1^N(x,Q^2)`$. In our fit to the data we have used the following expressions for $`g_1/F_1`$ and $`A_1`$:

$`\left[{\displaystyle \frac{g_1^N(x,Q^2)}{F_1^N(x,Q^2)}}\right]_{exp}`$ $``$ $`{\displaystyle \frac{g_1^N(x,Q^2)_{\mathrm{LT}}+h^N(x)/Q^2}{F_2^N(x,Q^2)_{exp}}}2x{\displaystyle \frac{[1+R(x,Q^2)_{exp}]}{(1+\gamma ^2)}},`$
$`A_1^N(x,Q^2)_{exp}`$ $``$ $`{\displaystyle \frac{g_1^N(x,Q^2)_{\mathrm{LT}}+h^N(x)/Q^2}{F_2^N(x,Q^2)_{exp}}}2x[1+R(x,Q^2)_{exp}].`$ (17)
where $`g_1^N(x,Q^2)_{\mathrm{LT}}`$ is given by the leading twist expression (15). In (17) $`h^N(x)`$ are a measure of the dynamical higher twists. In our analysis their $`Q^2`$ dependence is neglected. It is small and the accuracy of the present data does not allow to determine it. For the unpolarized structure functions $`F_2^N(x,Q^2)_{exp}`$ and $`R(x,Q^2)_{exp}`$ we have used the NMC parametrization  and the SLAC parametrization $`\mathrm{R}_{1998}`$ , respectively.

As in our previous analysis , for the input LO and NLO polarized parton densities at $`Q_0^2=1GeV^2`$ we have adopted a simple parametrization

$`x\mathrm{\Delta }u_v(x,Q_0^2)`$ $`=`$ $`\eta _uA_ux^{a_u}xu_v(x,Q_0^2),`$
$`x\mathrm{\Delta }d_v(x,Q_0^2)`$ $`=`$ $`\eta _dA_dx^{a_d}xd_v(x,Q_0^2),`$
$`x\mathrm{\Delta }s(x,Q_0^2)`$ $`=`$ $`\eta _sA_sx^{a_s}xs(x,Q_0^2),`$
$`x\mathrm{\Delta }G(x,Q_0^2)`$ $`=`$ $`\eta _gA_gx^{a_g}xG(x,Q_0^2),`$ (18)
where on RHS of (18) we have used the MRST98 (central gluon)  and MRST99 (central gluon)  parametrizations for the LO and NLO($`\overline{\mathrm{MS}}`$) unpolarized densities, respectively. The number of active flavors is $`\mathrm{N}_\mathrm{f}=3`$. The normalization factors $`A_i`$ in (18) are fixed such that $`\eta _i`$ are the first moments of the polarized densities. To fit better the data in LO QCD, an additional factor $`(1+\gamma _vx)`$ on RHS is used for the valence quarks. Bearing in mind that the light quark sea densities $`\mathrm{\Delta }\overline{u}`$ and $`\mathrm{\Delta }\overline{d}`$ cannot be, in principle, determined from the present inclusive data (in the absence of polararized charge current neutrino experiments) we have adopted the convention of a flavor symmetric sea

$$\mathrm{\Delta }u_{sea}=\mathrm{\Delta }\overline{u}=\mathrm{\Delta }d_{sea}=\mathrm{\Delta }\overline{d}=\mathrm{\Delta }s=\mathrm{\Delta }\overline{s}.$$
(19)
The first moments of the valence quark densities $`\eta _u`$ and $`\eta _d`$ are constrained by the sum rules

$$a_3=g_A=\mathrm{F}+\mathrm{D}=1.2670\pm 0.0035\text{[23]},$$
(20)
$$a_8=3\mathrm{F}\mathrm{D}=0.585\pm 0.025\text{[11]},$$
(21)
where $`a_3`$ and $`a_8`$ are non-singlet combinations of the first moments of the polarized parton densities corresponding to $`3^{\mathrm{rd}}`$ and $`8^{\mathrm{th}}`$ components of the axial vector Cabibbo current

$$a_3=(\mathrm{\Delta }u+\mathrm{\Delta }\overline{u})(Q^2)(\mathrm{\Delta }d+\mathrm{\Delta }\overline{d})(Q^2),$$
(22)
$$a_8=(\mathrm{\Delta }u+\mathrm{\Delta }\overline{u})(Q^2)+(\mathrm{\Delta }d+\mathrm{\Delta }\overline{d})(Q^2)2(\mathrm{\Delta }s+\mathrm{\Delta }\overline{s})(Q^2).$$
(23)
The sum rule (20) reflects isospin SU(2) symmetry, whereas (21) is a consequence of the $`SU(3)_f`$ flavor symmetry treatment of the hyperon $`\beta `$-decays. While the isospin symmetry is not in doubt, there is some question about the accuracy of assuming $`SU(3)_f`$ in analyzing hyperon $`\beta `$-decays. We have previously studied the sensitivity of the polarized parton densities to the deviation of $`a_8`$ from its SU(3) flavor symmetric value (0.58). The results are given in . In this analysis we will use for $`a_8`$ its SU(3) symmetric value (21).

In our past papers we have used the Jacobi polynomial method to yield the structure functions $`g_1^N(x,Q^2)_{\mathrm{LT}}`$ from their Mellin moments in n space. The details of this procedure are given in . But the accuracy of this method is limited in the low x region, $`x<0.01`$, so we have now used the inverse Mellin-transformation method (see, e.g., ) which reconstructs very precisely $`g_1^N(x,Q^2)_{\mathrm{LT}}`$ from its moments in the whole $`x`$ region. We have repeated our fits without including HT corrections in order to compare to our previous results  obtained by the Jacobi polynomial method. We have found very good agreement between the results obtained by both methods. The reason is that the present kinematic $`x`$ region of the polarized DIS data coincides with the domain where the Jocobi polynomial method works well. Also, the difference between the structure functions calculated by Jacobi and inverse Mellin-transformation methods is much smaller than the accuracy achieved in the present polarized DIS experiments. Nevertheless, bearing in mind its universality and, in particular, its applicability to the semi-inclusive DIS processes, we have decided to use the inverse Mellin-transformation approach in this analysis.

The unknown higher twists $`h^N(x)`$ in (17) have been extracted from the data following the method used in  and  for the higher twist corrections to the unpolarized structure functions $`F_2`$ and $`xF_3`$, respectivelyNote that the moments of the $`g_1`$ higher twists have been studied in the SLAC/E143 paper  as well as in . The measured $`x`$ region has been split into 5 bins and to any $`x`$-bin two parameters $`h_i^{(p)}`$ and $`h_i^{(n)}`$ have been attached. We have found that for a deutron the relation $`h_i^{(d)}=0.925(h_i^{(p)}+h_i^{(n)})/2`$ is a good approximation. So, to the parameters connected with the input PPD (18) we add the parameters $`h_i^{(p)}`$ and $`h_i^{(n)},(i=1,2,..,5`$).

All free parameters

$$\{a_u,a_d,a_s,a_g,\eta _s,\eta _g(\gamma _u,\gamma _d);h_i^{(p)},h_i^{(n)}\}$$
(24)
have been determined from the best fit to $`g_1/F_1`$ and $`A_1`$ data using (17), i.e., effectively by fitting $`(g_1)_{exp}`$. Note that in the calculations of $`g_1(x,Q^2)_{\mathrm{LT}}`$ we have used for the strong coupling constant $`\alpha _s(Q^2)`$ the same procedure as in our previous paper  (see the details in LSS2001 FORTRAN code at http://durpdg.dur.ac.uk/HEPDATA/PDF).

## 5 Results

In this section we present the numerical results of our fits to the world data on $`g_1/F_1`$  and $`A_1`$ . The data used (185 experimental points) cover the following kinematic region:

$$0.005x0.75,1<Q^258GeV^2.$$
(25)
The total (statistical and systematic) errors are taken into account. The systematic errors are added quadratically.

We prefer to discuss the results of the NLO analysis in the JET (or so-called chirally invariant) factorization scheme . In this scheme the first moment of singlet $`\mathrm{\Delta }\mathrm{\Sigma }(Q^2)`$, as well as the strange sea polarization $`(\mathrm{\Delta }s+\mathrm{\Delta }\overline{s})(Q^2)`$, are $`Q^\mathit{2}`$ independent quantities. Then, it is meaningful to directly interpret $`\mathrm{\Delta }\mathrm{\Sigma }`$ as the contribution of the quark spins to the nucleon spin and to compare its value obtained from the DIS region with the predictions of the different (constituent, chiral, etc.) quark models at low $`Q^2`$. Later we will briefly comment on the scheme dependence effects on the results of the analysis.

### 5.1 Higher twist effects

The numerical results of our fits to the data are summarized in Tables 1 and 2. As seen from the values of $`\chi ^2`$ per degree of freedom $`\chi _{\mathrm{DF}}^2`$ in Table 1, a very good description of the data is achieved. The best LO and NLO(JET) fits correspond to $`\chi _{\mathrm{DF},\mathrm{LO}}^2=0.892`$ and to $`\chi _{\mathrm{DF},\mathrm{NLO}}^2=0.858`$. We have found that the fit to the data is significantly improved, especially in the LO case, when the higher twist corrections to $`g_1^N`$ are included in the analysis (see Table 2). In contrast to the case when the HT corrections to $`g_1`$ are not taken into account in the fits, the value of $`\chi _{\mathrm{LO}}^2(\mathrm{HT})`$ is very close to that of $`\chi _{\mathrm{NLO}}^2(\mathrm{HT})`$, which is an indication that the tail of the neglected higher order logarithmic corrections to $`g_1`$ resemble a power behavior of order $`𝒪(1/Q^2)`$ . A similar behavior of $`\chi ^2`$ has been observed in the QCD analysis of the unpolarized structure function $`xF_3(x,Q^2)`$ in .

The extracted higher twist corrections to the proton and neutron spin structure functions, $`h^p(x)`$ and $`h^n(x)`$, are shown in Fig. 2. As seen from Fig. 2 the corrections for the proton and neutron have a different shape. While $`h^p(x)`$ changes sign in the LO as well in the NLO case, $`h^n(x)`$ is non-negative in the measured $`x`$ region in both cases. One can see also that the HT corrections to the proton structure function $`g_1^p`$ appear to be smaller when for $`(g_1)_{\mathrm{LT}}`$ the NLO approximation is used. In Fig. 3 we demonstrate how the choice of the factorization scheme for the perturbative calculation of $`(g_1)_{\mathrm{LT}}`$ influences the higher twists results. The results are presented for the JET and $`\overline{\mathrm{MS}}`$ schemes. It is seen that the HT corrections to $`g_1`$ in both cases coincide within the errors. The small deference between the central values could be considered as an estimation of the NNLO effects in $`(g_1)_{\mathrm{LT}}`$ In  the HT terms have been discussed using for them two phenomenological parametrizations. The authors conclude that they do not find a significant higher twist contribution to $`g_1`$ in an NLO treatment of $`(g_1)_{\mathrm{LT}}`$. On other hand, studying the higher twist effects in the moments of $`g_1`$, it was shown in  that while the first moment of higher twist is quite small, the higher order moments are relevant at $`Q^2`$few $`GeV^2`$..

### 5.2 NLO polarized parton densities

Let us discuss now the polarized parton densities extracted from the data in the presence of the HT corrections to $`g_1`$. We will call this set of parton densities PD($`g_1^{\mathrm{LT}}+\mathrm{HT}`$). In Fig. 4 we compare the NLO(JET) polarized PD($`g_1^{\mathrm{NLO}}+\mathrm{HT}`$) with those obtained in our analysis  where we performed fits to the data according to (9). We will call the latter PD($`g_1^{\mathrm{NLO}}/F_1^{\mathrm{NLO}}`$). As seen from Fig. 4 the two sets of polarized parton densities are very close to each other. This is a good illustration of the fact that a fit to the data on $`A_1(g_1/F_1)`$ using for the $`g_1`$ and $`F_1`$ structure functions their NLO leading twist expressions ($`\chi _{\mathrm{DF},\mathrm{NLO}}^2=0.859`$) is equivalent to a fit to the $`g_1`$ data taking into account the higher twist corrections to $`g_1`$ ($`\chi _{\mathrm{DF},\mathrm{NLO}}^2=0.858`$). In other words, this analysis confirms once more that the higher twist corrections to $`g_1`$ and $`F_1`$ approximately cancel in the ratio $`g_1/F_1`$.

### 5.3 LO polarized parton densities

Let us turn now to the LO polarized parton densities. In LO QCD $`\mathrm{\Delta }G(x,Q^2)`$ does not contribute directly to $`g_1`$ and the gluons cannot be determined from DIS data alone. For this reason the LO fit to the data was performed using for the input polarized gluon density $`\mathrm{\Delta }G(x,Q_0^2)`$ the one extracted in the NLO fit to the data:

$$\mathrm{\Delta }G(x,Q_0^2)_{\mathrm{LO}}=\mathrm{\Delta }G(x,Q_0^2)_{\mathrm{NLO}(\mathrm{JET})}.$$
(26)
It is important to note that in the polarized case the LO approximation has some peculiarities compared to the unpolarized one. As a consequence of the gluon axial anomaly, the difference between NLO anti-quark polarizations $`\mathrm{\Delta }\overline{q}_i`$ in different factorization schemes can be quite large, comparable in magnitude to the $`\mathrm{\Delta }\overline{q}_i`$ themselves (see, e.g., ). In this case the leading order will be a bad approximation, at least for the polarized sea-quark densities. Also, bearing in mind that in polarized DIS most of the data points are at low $`Q^2`$, lower than the usual cuts in the analyses of unpolarized data ($`Q^245GeV^2`$), the NLO corrections to all polarized parton densities are large in this region and it is better to take them into account. Nevertheless, the LO polarized parton densities may be useful for some practical purposes; e.g., for preliminary estimations of the cross sections in future polarized experiments, etc. They are also needed for comparison with those extracted from semi-inclusive DIS data , where the NLO QCD analysis is still very complicated. The extracted LO polarized parton densities PD($`g_1^{\mathrm{LO}}+HT`$) are shown in Fig. 5. Also shown in Fig. 5 are the LO polarized PD($`g_1^{\mathrm{LO}}/F_1^{\mathrm{LO}}`$) obtained in our analysis . In contrast to the NLO case, the two sets of LO polarized densities are significantly different. As a result we obtain different theoretical curves for $`g_1`$ (see Fig. 6). To illustrate how these curves fit the data, the SLAC/E143 experimental proton data at $`Q^2=5GeV^2`$ are also shown. As seen from Fig. 6, the proton structure function $`g_1^p(x,Q^2)_{LO}`$ calculated using the LO polarized PD($`g_1^{\mathrm{LO}}/F_1^{\mathrm{LO}}`$) does not agree with the data for $`x<0.25`$. Note that at the same time the ratio $`g_1^{\mathrm{LO}}/F_1^{\mathrm{LO}}`$ fits the world data on $`g_1/F_1`$ and $`A_1`$ quite well in the measured $`(x,Q^2)`$ region ($`\chi _{\mathrm{DF}}^2=0.921`$) . The main reason for this peculiarity is that the LO approximations for the unpolarized structure functions $`F_2(x,Q^2)`$ and $`F_1(x,Q^2)`$ (or $`R(x,Q^2))`$ presented in the literature are not self consistent. The unpolarized parton densities in leading order QCD (more correctly in leading logarithmic approximation LLA) are usually extracted including in the data set of the analysis only the experimental data on $`F_2`$ and ignoring the data on R (or F1). Remember that in LLA of QCD the structure functions satisfy the Callan-Gross relation $`2xF_1(x,Q^2)=F_2(x,Q^2)`$ which leads to $`R=4M^2x/Q^2`$. However, the experimental data on $`F_2`$ and $`R`$ do not satisfy these relations in a large kinematic region. They are approximately satisfied only at large $`x`$ and/or large $`Q^2`$. At small $`x`$ and $`Q^2`$ the experimental values of $`F_2`$ are larger than those of $`2xF_1`$ by up to 30%. That is why the extracted sets of LO unpolarized parton densities (MRST, CTEQ, etc.) are not quite consistent. While they fit well the data on $`F_2`$, they badly fail to describe the $`R(F_1)`$ data in the region of small $`x`$ and $`Q^2`$. One way to improve the situation would be to perform a LO QCD fit including in the data set the $`R(F_1)`$ data too. Also, if the data at low $`Q^2`$, lower than the usual cuts ($`Q^245GeV^2`$), are included in the analysis, the higher twist corrections to $`F_2`$ and $`F_1`$ should be taken into account.

## 6 Conclusion

We have analyzed the world data on inclusive polarized deep inelastic lepton-nucleon scattering in leading and next-to-leading order of QCD including in the analysis the higher twist $`h^N(x)/Q^2`$ and the target mass corrections to the nucleon spin structure function $`g_1^N(x,Q^2)`$. We have found that the fit to the data on $`g_1`$ is essentially improved, especially in the LO case, when the higher twist terms are included in the analysis. The $`x`$-dependence of the higher twists $`h^N(x)`$ have been extracted from the data in a model independent way. It is shown that the size of their contribution to $`g_1`$ is not negligible and their shape depends on the target:$`h^p(x)`$ changes sign while $`h^n(x)`$ is a non-negative function in the measured $`x`$ region.

We have found that the polarized parton densities depend on whether the higher twist terms are or are not included in the analysis of $`g_1`$. Moreover, the NLO polarized parton densities extracted from the $`g_1`$ data in the presence of higher twist terms are in good agreement with those determined by our previous fits  to the data on $`g_1/F_1`$ and $`A_1`$ using for the structure functions $`g_1`$ and $`F_1`$ only their leading twist expressions in NLO QCD. This observation confirms once more that the higher twist corrections to $`g_1/F_1`$ and $`A_1`$ are negligible so that in the analysis of $`g_1/F_1`$ and $`A_1`$ data it is enough to account only for the leading twist of the structure functions $`g_1`$ and $`F_1`$. On the other hand, in fits to the $`g_1`$ data themselves the higher twist contribution to $`g_1`$ must be taken into account. The latter is especially important for the LO QCD analysis of the inclusive and semi-inclusive DIS data.

Acknowledgments

This research was supported by a UK Royal Society Collaborative Grant, by the JINR-Bulgaria Collaborative Grant, by the RFBR (No 00-02-16696), INTAS 2000 (No 587) and by the Bulgarian National Science Foundation under Contract Ph-1010.

|  | Table 1. Parameters of the LO and NLO(JET) input parton densities at |
| --- | --- |
|  | $`Q^2=1GeV^2`$ as obtained from the best fits to the world $`g_1/F_1`$ and $`A_1^N`$ data |
|  | including the HT corrections to $`g_1`$. The errors shown are total (statistical and |
|  | systematic). The parameters marked by (\*) are fixed. Note that the TMC are |
|  | included in $`(g_1)_{\mathrm{LT}}`$. |

| Fit | $`(g_1)_{LO}+h(x)/Q^2`$ | $`(g_1)_{NLO}+h(x)/Q^2`$ |
| --- | --- | --- |
| $`\mathrm{DF}`$ | 185 - 16 | 185 - 16 |
| $`\chi ^2`$ | 150.7 | 145.0 |
| $`\chi ^2/\mathrm{DF}`$ | 0.892 | 0.858 |
| $`\eta _u`$ | 0.926 | $`0.926^{}`$ |
| $`a_u`$ | 0.000 $`\pm `$ 0.002 | 0.312 $`\pm `$ 0.048 |
| $`\gamma _u`$ | 1.556 $`\pm `$ 0.261 | $`0^{}`$ |
| $`\eta _d`$ | \- 0.341 | $`0.341^{}`$ |
| $`a_d`$ | 0.000 $`\pm `$ 0.005 | 0.000 $`\pm `$ 0.049 |
| $`\gamma _d`$ | 2.808 $`\pm `$ 1.249 | $`0^{}`$ |
| $`\eta _s`$ | \- 0.072 $`\pm `$ 0.008 | \- 0.045 $`\pm `$ 0.007 |
| $`a_s`$ | 0.601 $`\pm `$ 0.064 | 1.583 $`\pm `$ 0.434 |
| $`\eta _g`$ | $`0.803^{}`$ | 0.803 $`\pm `$ 0.244 |
| $`a_g`$ | $`0.376^{}`$ | 0.376 $`\pm `$ 0.503 |
| $`x_i`$ | $`h^p(x_i)[GeV^2]`$ | |
| 0.028 | 0.013 $`\pm `$ 0.036 | 0.064 $`\pm `$ 0.044 |
| 0.100 | \- 0.076 $`\pm `$ 0.032 | \- 0.007 $`\pm `$ 0.034 |
| 0.200 | \- 0.145 $`\pm `$ 0.032 | \- 0.060 $`\pm `$ 0.035 |
| 0.350 | \- 0.030 $`\pm `$ 0.035 | \- 0.008 $`\pm `$ 0.038 |
| 0.600 | 0.035 $`\pm `$ 0.019 | 0.026 $`\pm `$ 0.021 |
| $`x_i`$ | $`h^n(x_i)[GeV^2]`$ | |
| 0.028 | 0.234 $`\pm `$ 0.073 | 0.178 $`\pm `$ 0.078 |
| 0.100 | 0.192 $`\pm `$ 0.048 | 0.199 $`\pm `$ 0.050 |
| 0.200 | 0.035 $`\pm `$ 0.056 | 0.079 $`\pm `$ 0.059 |
| 0.325 | 0.072 $`\pm `$ 0.071 | 0.055 $`\pm `$ 0.073 |
| 0.500 | 0.023 $`\pm `$ 0.043 | \- 0.020 $`\pm `$ 0.040 |

|  | Table 2. The values of $`\chi ^2`$ for the LO and NLO QCD fits without HT included |
| --- | --- |
|  | compared to those when the HT corrections to $`g_1`$ are taken into account. The |
|  | TMC are included in $`(g_1)_{\mathrm{LT}}`$ . |

| Fit | LO(HT=0) | NLO(HT=0) | LO+HT | NLO+HT |
| --- | --- | --- | --- | --- |
| $`\chi ^2`$ | 244.5 | 218.8 | 150.9 | 145.0 |
| $`\mathrm{DF}`$ | 185 - 6 | 185 - 6 | 185 - 16 | 185 - 16 |
| $`\chi ^2/\mathrm{DF}`$ | 1.36 | 1.22 | 0.893 | 0.858 |

Figure Captions

Fig.1. Effective higher twist contribution $`h^{A_1}(x)`$ to the spin asymmetry $`A_1^N(x,Q^2)`$ extracted from the data. Compared to our 1999 result (Fig. 4  and Fig. 1 ) we present here the results of a new analysis including in the data set the SLAC/155 proton data not available at that time.

Fig. 2. Higher twist corrections to the proton and neutron $`g_1`$ structure functions extracted from the data on $`g_1`$ in the case of LO and NLO QCD approximation for $`g_1(x,Q^2)_{\mathrm{LT}}`$.

Fig. 3. Higher twist corrections to the proton and neutron $`g_1`$ structure functions extracted from the data when the leading twist calculations of $`g_1(x,Q^2)_{\mathrm{NLO}}`$ are performed in different factorization schemes.

Fig. 4. NLO(JET) polarized parton densities PD($`g_1^{\mathrm{NLO}}+\mathrm{HT}`$) (solid curves) compared to PD($`g_1^{\mathrm{NLO}}/F_1^{\mathrm{NLO}}`$) (dashed curves) at $`Q^2=1GeV^2`$ (see the text).

Fig. 5. LO polarized parton densities PD($`g_1^{\mathrm{LO}}+\mathrm{HT}`$) (solid curves) compared to PD($`g_1^{\mathrm{LO}}/F_1^{\mathrm{LO}}`$) (dashed curves) at $`Q^2=1GeV^2`$ (see the text).

Fig. 6. Comparison of the proton structure function $`g_1^{\mathrm{LO}}`$ calculated using the polarized parton densities PD($`g_1^{\mathrm{LO}}+\mathrm{HT}`$) (solid curve) and PD($`g_1^{\mathrm{LO}}/F_1^{\mathrm{LO}}`$) (dashed curve) with SLAC/E143 proton data.


text:
# remote information concentration by GHZ state and by bound entangled state

## I Introduction

As a peculiar feature of quantum formalism without classical counterpart, quantum entanglement has already been at the heart of the emerging quantum information theory. The two-qubit maximally entangled state has produced highly nonintuitive effects such as quantum teleportation, quantum dense coding, quantum cryptography. In practice, one usually deals with noisy entanglement represented by mixed state of a composite system. All the noisy entanglement of two-qubit system can be distilled to the singlet form. But, for the mixed entangled states of more-than-two-qubit systems, there are two qualitatively different kinds of entanglement : free entanglement which can always be distilled and bound one which cannot be brought to the singlet form only by local quantum operations and classical communication (LOCC). In the last few years, the bound entangled states of the multipartite systems is under extensive research because it arouses a deeper understanding of the entanglement and the nonlocality of quantum states.

On the other hand, the bound entanglement seems to be useless for quantum information work such as reliable transmission of quantum data via teleportation. However, it is not the case. Resent research has shown that in a sense many copies of a bound entangled state can be pumped into a free entangled state. And it has been found that entanglement can be distilled from two copies of a bound entangled state distributed to different parties. Namely, even a small amount of bound entanglement can be activated and used to process quantum information. More recently, it has been demonstrated that a single copy of a bound entangled state can perform information work better than classically correlated states. That is, a single bound entangled state alone is useful for quantum information processing. Here we give an example to demonstrate that in a sense a single bound entangled state performs as good as the maximally entangled GHZ state.

In this paper we consider that quantum information originally from a single qubit, but now asymmetrically distributed into three spatially separated qubits, is remote concentrated back to a single qubit by an initially shared four-partite GHZ state without performing any global operations. This quantum information work can also be achieved by a initially shared four-partite unlockable bound entangled state. We discuss remote information concentration in detail by GHZ state in Sec.II and by the bound entangled state in Sec.III, then give our discussion and conclusion in Sec.IV.

## II Remote information concentration by GHZ state

For simplicity, we focus on the protocol of three parties.

Supposing that three separate parties Alice, Bob and Charlie hold three qubits A, B and C, respectively. The three qubits are in an asymmetric telecloning state

$$|\psi _{ABC}=\alpha |\varphi _0_{ABC}+\beta |\varphi _1_{ABC},$$
(1)
where $`\alpha `$ and $`\beta `$ are real numbers and satisfies $`\alpha ^2+\beta ^2=1`$. The states $`|\varphi _0_{ABC}`$ and $`|\varphi _1_{ABC}`$ is defined as

$`|\varphi _0_{ABC}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N}}}(|0_A|0_B|0_C+p|1_A|0_B|1_C+q|1_A|1_B|0_C),`$ (2)
$`|\varphi _1_{ABC}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N}}}(|1_A|1_B|1_C+p|0_A|1_B|0_C+q|0_A|0_B|1_C),`$ (3)
where $`q=1p,p>q,N`$ is a normalization factor given by $`N=1+p^2+q^2`$. That is, the information of an unknown state $`|\chi =\alpha |0+\beta |1`$ is diluted to a composite system consisting of qubits A, B and C by an asymmetric telecloning process. Now the task is to concentrate the information back to a single qubit without the collective operations between the three qubits. The schematic picture of the remote information concentration is illustrated in Fig.1.

The first protocol we proposed is to share in advance a four-partite GHZ state

$$|\phi _{DEFG}=\frac{1}{\sqrt{2}}(|0_D|0_E|0_F|0_G+|1_D|1_E|1_F|1_G)$$
(4)
between three parties Alice, Bob and Charlie, which hold the qubits E, F and G, respectively. Let David hold the qubit D. Then the three parties Alice, Bob and Charlie perform the Bell-state measurements (BSMs) on the respective pairs of qubits, and inform David of the results of the measurements. Each of the three parties gets one of four possible outputs $`\left\{|\mathrm{\Phi }^i_{(i=0,1,2,3)}\right\}`$ of the BSM, where $`|\mathrm{\Phi }^i_{(i=0,1,2,3)}`$ represent the four Bell states $`|\mathrm{\Phi }^{0,1}=\left(|00\pm |11\right)/\sqrt{2}`$, $`|\mathrm{\Phi }^{2,3}=\left(|01\pm |01\right)/\sqrt{2}`$. Each of the four possible outputs is associated with the corresponding Pauli operators $`\left\{\sigma _{(i=0,1,2,3)}^i\right\}`$, where $`\sigma _0=I,\sigma _1=\sigma _z,\sigma _2=\sigma _x,\sigma _3=\sigma _y`$. According to their results of the measurements David rotates his qubit D, finally obtains the state $`|\chi `$ on the qubit D.

The information concentration can be understood in a formular way. We write the states $`|\psi _{ABC}`$ and $`|\phi _{DEFG}`$ together:

$`|\psi _{ABC}|\phi _{DEFG}`$ $`=`$ $`\left(\alpha |\varphi _0_{ABC}+\beta |\varphi _1_{ABC}\right){\displaystyle \frac{1}{\sqrt{2}}}(|0_D|0_E|0_F|0_G+|1_D|1_E|1_F|1_G)`$ (5)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2N}}}\{\alpha (|00_{AE}|00_{BF}|00_{CG}|0_D+p|10_{AE}|00_{BF}|10_{CG}|0_D+q|10_{AE}|10_{BF}|00_{CG}|0_D`$ (9)
$`+|01_{AE}|01_{BF}|01_{CG}|1_D+p|11_{AE}|01_{BF}|11_{CG}|1_D+q|11_{AE}|11_{BF}|01_{CG}|1_D)`$
$`+\beta (|10_{AE}|10_{BF}|10_{CG}|0_D+q|00_{AE}|00_{BF}|10_{CG}|0_D+p|00_{AE}|10_{BF}|00_{CG}|0_D`$
$`+|11_{AE}|11_{BF}|11_{CG}|1_D+q|01_{AE}|01_{BF}|11_{CG}|1_D+p|01_{AE}|11_{BF}|01_{CG}|1_D)\}.`$
A lengthy but straightward calculation gives that if the combining result of the three BSMs is one of the set $`\{\mathrm{\Phi }_{AE}^0\mathrm{\Phi }_{BF}^0\mathrm{\Phi }_{CG}^0,\mathrm{\Phi }_{AE}^0\mathrm{\Phi }_{BF}^1\mathrm{\Phi }_{CG}^1,\mathrm{\Phi }_{AE}^1\mathrm{\Phi }_{BF}^0\mathrm{\Phi }_{CG}^1,\mathrm{\Phi }_{AE}^1\mathrm{\Phi }_{BF}^1\mathrm{\Phi }_{CG}^0,\mathrm{\Phi }_{AE}^2\mathrm{\Phi }_{BF}^0\mathrm{\Phi }_{CG}^2,\mathrm{\Phi }_{AE}^3\mathrm{\Phi }_{BF}^0\mathrm{\Phi }_{CG}^3,\mathrm{\Phi }_{AE}^3\mathrm{\Phi }_{BF}^1\mathrm{\Phi }_{CG}^2,`$

$`\mathrm{\Phi }_{AE}^2\mathrm{\Phi }_{BF}^1\mathrm{\Phi }_{CG}^3,\mathrm{\Phi }_{AE}^2\mathrm{\Phi }_{BF}^2\mathrm{\Phi }_{CG}^0,\mathrm{\Phi }_{AE}^3\mathrm{\Phi }_{BF}^3\mathrm{\Phi }_{CG}^0,\mathrm{\Phi }_{AE}^2\mathrm{\Phi }_{BF}^3\mathrm{\Phi }_{CG}^1,\mathrm{\Phi }_{AE}^3\mathrm{\Phi }_{BF}^2\mathrm{\Phi }_{CG}^1\}`$, which are associated with the Pauli operator $`\sigma _0=I`$, the state of the qubit D is projected to the state $`|\chi =\alpha |0+\beta |1`$. Analogously, there is a Pauli operator $`\sigma ^i`$ pertaining to the combining result of the three BSMs from each of other sets. The Pauli operator $`\sigma ^i`$ is the product (up to a global phase) of the respective Pauli operators pertaining to the three BSMs. That is, up a global phase factor, $`\sigma ^i`$ is equal to $`\sigma _{AE}^l\sigma _{BF}^j\sigma _{CG}^k`$, where $`l,j,k=0,1,2,3`$ and the subscripts denote the BSMs on the corresponding pairs of qubits. Then David performs the operator $`\sigma ^i`$ on his qubit D to retrieve the state $`|\chi `$. So the protocol achieves the task: $`|\psi _{ABC}|\chi `$ with certainty by pre-sharing a four-partite GHZ state. The task of remote information concentration in  is the special case when $`p=q=\frac{1}{2}`$.

It is worth noting that the distribution probability of the combining results of the three BSMs reveals some information of the state $`|\psi _{ABC}`$. For instance, while we get $`\mathrm{\Phi }_{AE}^0\mathrm{\Phi }_{BF}^0\mathrm{\Phi }_{CG}^0`$ with the probability of $`\frac{1}{16N}`$, we obtain $`\mathrm{\Phi }_{AE}^2\mathrm{\Phi }_{BF}^0\mathrm{\Phi }_{CG}^2`$ with the probability of $`\frac{p^2}{16N}`$ and $`\mathrm{\Phi }_{AE}^2\mathrm{\Phi }_{BF}^2\mathrm{\Phi }_{CG}^0`$ with the probability of $`\frac{q^2}{16N}`$.

Consequently, by utilizing a pre-sharing four-partite GHZ state a remote information concentration by LOCC can be achieved with certainty. And the distribution of the combing results of the three BSMs reveals some information of the state $`|\psi _{ABC}`$.

## III Remote information concentration by bound entangled state

In the above discussion a remote information concentration is achieved by sharing a GHZ state and LOCC. In the following we demonstrate that the information work can also be implemented by an unlockable bound entangled state and LOCC.

In order to concentrate the information diluted in the state $`|\psi _{ABC}`$ back to the state $`|\chi `$ on a single qubit, a four-partite unlockable bound entangled state

$$\rho _{DEFG}^{ub}=\frac{1}{4}\underset{i=0}{\overset{3}{}}|\mathrm{\Phi }_i_{DE}\mathrm{\Phi }_i||\mathrm{\Phi }_i_{FG}\mathrm{\Phi }_i|,$$
(10)
where $`|\mathrm{\Phi }_i`$ is defined as before, is pre-shared between Alice, Bob, Charlie and David. In the same way as above, the qubit E is sent to Alice, the qubit F to Bob, the qubit G to Charlie and the qubit D to David. No joint operations between qubits belonging to different parties is allowed. The three parties Alice, Bob, and Charlie perform the Bell-state measurements on their respective pairs of qubits in hand. Likewise, each of them obtains one of the possible outputs $`\left\{|\mathrm{\Phi }^i_{(i=0,1,2,3)}\right\}`$ of the BSM, which is associated with a corresponding Pauli operator in the set $`\left\{\sigma _{(i=0,1,2,3)}^i\right\}`$, and communicates the result with David, respectively. David determines a Pauli operator $`\sigma ^i`$ on his qubit D for retrieving the state $`|\chi `$ on the qubit D, according the product of the three Pauli operators pertaining to each one of the three BSMs. Finally, the pure state $`|\chi `$ comes out on the qubit D.

The process also can be expressed in a formular way. We joint the state $`|\psi _{ABC}`$ with $`\rho _{DEFG}^{ub}`$:

$`|\psi _{ABC}\rho _{DEFG}^{ub}_{ABC}\psi |`$ $`=`$ $`(\alpha |\varphi _0_{ABC}+\beta |\varphi _1_{ABC})\rho _{DEFG}^{ub}(\alpha _{ABC}^{}\varphi _0|+\beta _{ABC}^{}\varphi _1|)`$ (11)
$``$ $`{\displaystyle \frac{1}{64}}\{(|II|)_{AE,BF,CG}\sigma _0(\alpha |0_D+\beta |1_D)(\alpha _D^{}0|+\beta _D^{}1|)\sigma _0`$ (15)
$`+(|IIII|)_{AE,BF,CG}\sigma _1(\alpha |0_D+\beta |1_D)(\alpha _D^{}0|+\beta _D^{}1|)\sigma _1`$
$`+(|IIIIII|)_{AE,BF,CG}\sigma _2(\alpha |0_D+\beta |1_D)(\alpha _D^{}0|+\beta _D^{}1|)\sigma _2`$
$`+(|IVIV|)_{AE,BF,CG}\sigma _3(\alpha |0_D+\beta |1_D)(\alpha _D^{}0|+\beta _D^{}1|)\sigma _3\},`$
where$`|II|`$ marks a set of the combining results of the three BSMs as

$`|II|`$ $`=`$ $`|\mathrm{\Phi }^0_{AE}\mathrm{\Phi }^0||\mathrm{\Phi }^0_{BF}\mathrm{\Phi }^0||\mathrm{\Phi }^0_{CG}\mathrm{\Phi }^0|+|\mathrm{\Phi }^0_{AE}\mathrm{\Phi }^0||\mathrm{\Phi }^1_{BF}\mathrm{\Phi }^1||\mathrm{\Phi }^1_{CG}\mathrm{\Phi }^1|`$ (23)
$`+|\mathrm{\Phi }^1_{AE}\mathrm{\Phi }^1||\mathrm{\Phi }^0_{BF}\mathrm{\Phi }^0||\mathrm{\Phi }^1_{CG}\mathrm{\Phi }^1|+|\mathrm{\Phi }^1_{AE}\mathrm{\Phi }^1||\mathrm{\Phi }^1_{BF}\mathrm{\Phi }^1||\mathrm{\Phi }^0_{CG}\mathrm{\Phi }^0|`$
$`+|\mathrm{\Phi }^2_{AE}\mathrm{\Phi }^2||\mathrm{\Phi }^2_{BF}\mathrm{\Phi }^2||\mathrm{\Phi }^0_{CG}\mathrm{\Phi }^0|+|\mathrm{\Phi }^2_{AE}\mathrm{\Phi }^2||\mathrm{\Phi }^3_{BF}\mathrm{\Phi }^3||\mathrm{\Phi }^1_{CG}\mathrm{\Phi }^1|`$
$`+|\mathrm{\Phi }^3_{AE}\mathrm{\Phi }^3||\mathrm{\Phi }^2_{BF}\mathrm{\Phi }^2||\mathrm{\Phi }^1_{CG}\mathrm{\Phi }^1|+|\mathrm{\Phi }^3_{AE}\mathrm{\Phi }^3||\mathrm{\Phi }^3_{BF}\mathrm{\Phi }^3||\mathrm{\Phi }^0_{CG}\mathrm{\Phi }^0|`$
$`+|\mathrm{\Phi }^0_{AE}\mathrm{\Phi }^0||\mathrm{\Phi }^2_{BF}\mathrm{\Phi }^2||\mathrm{\Phi }^2_{CG}\mathrm{\Phi }^2|+|\mathrm{\Phi }^0_{AE}\mathrm{\Phi }^0||\mathrm{\Phi }^3_{BF}\mathrm{\Phi }^3||\mathrm{\Phi }^3_{CG}\mathrm{\Phi }^3|`$
$`+|\mathrm{\Phi }^1_{AE}\mathrm{\Phi }^1||\mathrm{\Phi }^2_{BF}\mathrm{\Phi }^2||\mathrm{\Phi }^3_{CG}\mathrm{\Phi }^3|+|\mathrm{\Phi }^1_{AE}\mathrm{\Phi }^1||\mathrm{\Phi }^3_{BF}\mathrm{\Phi }^3||\mathrm{\Phi }^2_{CG}\mathrm{\Phi }^2|`$
$`+|\mathrm{\Phi }^2_{AE}\mathrm{\Phi }^2||\mathrm{\Phi }^0_{BF}\mathrm{\Phi }^0||\mathrm{\Phi }^2_{CG}\mathrm{\Phi }^2|+|\mathrm{\Phi }^2_{AE}\mathrm{\Phi }^2||\mathrm{\Phi }^1_{BF}\mathrm{\Phi }^1||\mathrm{\Phi }^3_{CG}\mathrm{\Phi }^3|`$
$`+|\mathrm{\Phi }^3_{AE}\mathrm{\Phi }^3||\mathrm{\Phi }^0_{BF}\mathrm{\Phi }^0||\mathrm{\Phi }^3_{CG}\mathrm{\Phi }^3|+|\mathrm{\Phi }^3_{AE}\mathrm{\Phi }^3||\mathrm{\Phi }^1_{BF}\mathrm{\Phi }^1||\mathrm{\Phi }^2_{CG}\mathrm{\Phi }^2|.`$
It is clear that the combing results of the three BSMs $`|II|`$ correspond to the case that the product $`\sigma _{AE}^l\sigma _{BF}^j\sigma _{CG(l,j,k=0,1,2,3)}^k`$ of the three Pauli operators pertaining to each one of the three BSMs amounts to $`\sigma ^0`$. So the state of the qubit D is projected to the state $`|\chi `$. The rest may be deduced by analogy that $`|IIII|`$ corresponds to the product of $`\sigma ^1`$, $`|IIIIII|`$ to the product of $`\sigma ^2`$, $`|IVIV|`$ to the product of $`\sigma ^3`$. Correspondingly, David does a Pauli operator $`\sigma ^1`$, $`\sigma ^2`$ and $`\sigma ^3`$ on his qubit D in order to get the correct state $`|\chi `$, respectively. In the end, David gets a qubit D in the state $`|\chi `$ with certainty.

Hence, with exploiting an unlockable bound entangled state the same remote information concentration is achieved only by LOCC as with a maximally entangled GHZ state.

However, in the protocol with the bound entangled state, the distribution of the combining results donot reveal any information of the state $`|\psi _{ABC}`$. A complicated calculation displays that each of all combining results of the three BSMs comes with the equal probability of $`\frac{1}{64}`$. That is, taken for example, the combining result $`\mathrm{\Phi }_{AE}^0\mathrm{\Phi }_{BF}^0\mathrm{\Phi }_{CG}^0`$ comes out with the probability of $`\frac{1}{64}`$, the other combining results $`\mathrm{\Phi }_{AE}^2\mathrm{\Phi }_{BF}^0\mathrm{\Phi }_{CG}^2`$ and $`\mathrm{\Phi }_{AE}^2\mathrm{\Phi }_{BF}^2\mathrm{\Phi }_{CG}^0`$ do so. Allowing for communication security, the unlockable bound entangled state is more suitable for the remote information concentration than the maximally entangled GHZ state.

## IV Conclusion

In the paper, we consider a quantum information work, where the quantum information initially in the state $`|\chi `$ of a single qubit, but now distributed into three spatially separated qubits, is remotely concentrated back in the state $`|\chi `$ of a single qubit. It is found that a maximally entangled GHZ state and an unlockable bound entangled state all can do the work with LOCC. It gives a demonstration that in a sense, a single copy of a bound entangled state is still as useful in quantum information processing as a maximally entangled GHZ state. Even in the view of communication security, the bound entangled state does the work better than the GHZ state, because for the former the distribution of the combining results of the three local BSMs does not reveal any information about the input state $`|\psi _{ABC}`$.

As analyzed in , it is the entanglement existing in the state $`|\psi _{ABC}`$ that liberates the bound entangled state for transmitting the quantum information. The entanglement in the error correction state $`|\psi _{ABC}`$ is crucial for the remote information concentration with the bound entangled state and LOCC. It evokes another question about to what extent the bound entangled state cannot work for this information work. The consideration is helpful for deep understanding the relation between the free entanglement and the bound one, and quantum entanglement itself. We hope that our work will stimulate more research into the nature of entanglement and its applications in the quantum communication.

Acknowledgments

This work has been financially supported by the National Natural Science Foundation of China under the Grant No.10074072.