text:
# 1 Introduction

## 1 Introduction

It is well known that in the processes of $`J/\psi `$ meson photoproduction on protons at high energies the photon-gluon fusion partonic subprocess dominates . In the framework of general factorization approach of QCD the $`J/\psi `$ photoproduction cross section depends on the gluon distribution function in a proton, the hard amplitude of $`c\overline{c}`$-pair production as well as the mechanism of a creation colorless final state with quantum numbers of the $`J/\psi `$ meson. Such a way, we suppose that the soft interactions in the initial state are described by introducing a gluon distribution function, the hard partonic amplitude is calculated using perturbative theory of QCD at order in $`\alpha _s(\mu ^2)`$, where $`\mu m_c`$, and the soft process of the $`c\overline{c}`$-pair transition into $`J/\psi `$ meson is described in nonrelativistic approximation using series in the small parameters $`\alpha _s`$ and $`v`$ (relative velocity of the quarks in $`J/\psi `$ meson). As it is statement in nonrelativistic QCD (NRQCD) , there are colour singlet mechanism, where $`c\overline{c}`$-pair is hardly produced in the colour singlet state, and colour octet mechanism, where $`c\overline{c}`$-pair is produced in the color octet state and at long distance it transforms into final colour singlet state in soft process. However, as it was shown in paper , the data from HERA collider  in the wide region of $`p_{}`$ and $`z`$ may be described well in the framework of the colour singlet model and not colour octet contribution is needed. Based on the above mentioned result we will take into account in our analysis only the colour singlet model contribution in the $`J/\psi `$ meson photoproduction . We consider a role of a proton gluon distribution function in the $`J/\psi `$ photoproduction in the framework of the conventional parton model as well as in the framework of the $`k_{}`$-factorization approach . Last one is based on BFKL evolution equation , which takes into account large terms proportional to $`\mathrm{log}1/x`$ and $`\mathrm{log}\mu ^2/\mathrm{\Lambda }_{QCD}^2`$ opposite to DGLAP evolution equation , where only large logarithmic terms $`\mathrm{log}\mu ^2/\mathrm{\Lambda }_{QCD}^2`$ are taken into account. In the process of the $`J/\psi `$ meson photoproduction one has $`x\mu ^2/s`$ and $`\mu ^2M^2`$, where $`M`$ is the $`J/\psi `$ meson mass, $`s`$ is the square of a total energy of colliding particles in center of mass frame.

## 2 The cross section for $`\gamma pJ/\psi X`$ in $`k_{}`$-factorization approach

Nowadays, there are two approaches in calculating of $`J/\psi `$ meson or heavy quark production cross sections at high energies. In the conventional parton model  it is suggested that hadronic cross section $`\sigma (\gamma pJ/\psi X,s)`$ and relevant partonic cross section $`\widehat{\sigma }(\gamma gJ/\psi g,\widehat{s})`$ are connected as follows:

$$\sigma ^{PM}(\gamma pJ/\psi X,s)=𝑑x\widehat{\sigma }(\gamma gJ/\psi g,\widehat{s})G(x,\mu ^2),$$
(1)
where $`\widehat{s}=xs`$, $`G(x,\mu ^2)`$ is the collinear gluon distribution function in a proton, $`x`$ is the fraction of a proton momentum, $`\mu ^2`$ is the typical scale of a hard process. The $`\mu ^2`$ dependence of the gluon distribution $`G(x,\mu ^2)`$ is described by DGLAP evolution equation . In the region of very high $`s`$ one has $`x<<1`$. This fact leads to BFKL evolution equation  for the unintegrated gluon distribution function $`\mathrm{\Phi }(x,𝐪_{}^\mathrm{𝟐},\mu ^\mathrm{𝟐})`$, where $`𝐪_{}^\mathrm{𝟐}`$ is the gluon virtuality. The unintegrated gluon distribution function can be related to the conventional gluon distribution by

$$xG(x,\mu ^2)=_0^{\mu ^2}\mathrm{\Phi }(x,𝐪_{}^2,\mu ^2)𝑑𝐪_{}^2.$$
(2)
The gluon 4-momentum is presented as follows

$$q=xp_N+q_{},$$

where $`q_{}=(0,𝐪_{},0)`$, $`p_N=(E_N,0,0,|𝐩_N|)`$ and $`q^2=q_{}^2=𝐪_{}^{}{}_{}{}^{2}`$. So called BFKL gluon is off mass-shell and it has polarization vector along its transverse momentum such as $`\epsilon ^\mu =q_{}^\mu /|𝐪_{}|`$. In $`k_{}`$-factorization approach hadronic and partonic cross sections are related by the following condition:

$$\sigma (\gamma pJ/\psi X)=\frac{dx}{x}_0^{\mu ^2}𝑑𝐪_{}^2\frac{d\varphi }{2\pi }\mathrm{\Phi }(x,𝐪_{}^2,\mu ^2)\widehat{\sigma }(\gamma g^{}J/\psi g,\widehat{s},𝐪_{}^2),$$
(3)
where $`\widehat{\sigma }(\gamma g^{}J/\psi g)`$ is $`J/\psi `$ meson photoproduction on BFKL gluon, $`\varphi `$ is azimuthal angle in the transverse $`XY`$ plane between vector $`𝐪_{}`$ and fixed $`OX`$ axis.

## 3 Unintegrated gluon distribution function

At the present time an exact form of the unintegrated gluon distribution $`\mathrm{\Phi }(x,𝐪_{}^2,\mu ^2)`$ is unknown yet because the relevant experimental analysis of the experimental data has never been carried out. There are several theoretical approximations for $`\mathrm{\Phi }(x,𝐪_{}^2,\mu ^2)`$, which are based on solving BFKL evolution equation . In the region of very small $`x0.01`$ and moderate $`𝐪_{}^2`$ ($`10`$ GeV2), which is relevant for $`J/\mathrm{\Psi }`$ photoproduction at HERA, all parameterizations  are approximately coincide (see, for example, discussions in  and ). For our purposes, we will use the unintegrated gluon distribution, which was obtained in the paper . The proposed method lies upon a straightforward perturbative solution of the BFKL equation where the collinear gluon density $`G(x,\mu ^2)`$ is used as the boundary condition in the integral form (2). Technically, the unintegrated gluon distribution is calculated as a convolution of collinear gluon distribution $`G(x,\mu ^2)`$ with universal weight factors:

$$\mathrm{\Phi }(x,𝐪_{}^2,\mu ^2)=_x^1𝒢(\eta ,𝐪_{}^2,\mu ^2)\frac{x}{\eta }G(\frac{x}{\eta },\mu ^2)𝑑\eta ,$$
(4)
where

$`𝒢(\eta ,𝐪_{}^2,\mu ^2)={\displaystyle \frac{\overline{\alpha }_s}{\eta 𝐪_{}^2}}J_0\left[2\sqrt{\overline{\alpha }_s\mathrm{ln}({\displaystyle \frac{1}{\eta }})\mathrm{ln}({\displaystyle \frac{\mu ^2}{𝐪_{}^2}})}\right],\text{ at }𝐪_{}^2\mu ^2`$ (5)
$`𝒢(\eta ,𝐪_{}^2,\mu ^2)={\displaystyle \frac{\overline{\alpha }_s}{\eta 𝐪_{}^2}}I_0\left[2\sqrt{\overline{\alpha }_s\mathrm{ln}({\displaystyle \frac{1}{\eta }})\mathrm{ln}({\displaystyle \frac{𝐪_{}^2}{\mu ^2}})}\right],\text{ at }𝐪_{}^2>\mu ^2,`$ (6)
$`J_0`$ and $`I_0`$ stand for Bessel functions (of real and imaginary arguments, respectively), and $`\overline{\alpha }_s={\displaystyle \frac{\alpha _s}{3\pi }}`$. As a input function $`G(x,\mu ^2)`$ in (5,6) we use the standard GRV parameterization . To test the method of calculation for $`\mathrm{\Phi }(x,𝐪_{}^2,\mu ^2)`$ we compare input collinear gluon distribution $`G(x,\mu ^2)`$ and collinear gluon distribution $`\stackrel{~}{G}(x,\mu ^2)`$, which is obtained using formula (2) from the unintegrated distribution $`\mathrm{\Phi }(x,𝐪_{}^2,\mu ^2)`$ after integration over $`𝐪_{}^2`$. In the Figure 1 the result of our calculation for the ratio $`R_1(x,\mu ^2)=\stackrel{~}{G}(x,\mu ^2)/G(x,\mu ^2)`$ is shown. It is visible, that at $`x<0.1`$ the ratio $`R_1`$ does not differ from 1 more than 1-2 %. Note, that in a process of the $`J/\psi `$ photoproduction on protons $`x{\displaystyle \frac{M}{\sqrt{s}}}`$ ($`M`$ is the mass of $`J/\psi `$ meson), and one has $`x0.03`$ at $`\sqrt{s}=100`$ GeV.

## 4 Amplitude for $`\gamma gJ/\psi g`$ process

There are six Feynman diagrams (Figure 2) which describere partonic process $`\gamma gJ/\psi g`$ at the leading order in $`\alpha _s`$ and $`\alpha `$. In the framework of the colour singlet model and nonrelativistic approximation the production of $`J/\psi `$ meson is considered as the production of a quark-antiquark system in the colour singlet state with orbital momentum $`L=0`$ and spin momentum $`S=1`$. The binding energy and relative momentum of quarks in the $`J/\psi `$ are neglected. Such a way $`M=2m_c`$ and $`p_c=p_{\overline{c}}={\displaystyle \frac{p}{2}}`$, where $`p`$ is 4-momentum of the $`J/\psi `$, $`p_c`$ and $`p_{\overline{c}}`$ are 4-momenta of quark and antiquark. Taking into account the formalism of the projection operator  the amplitude of the process $`\gamma gJ/\psi g`$ may be obtained from the amplitude of the process $`\gamma g\overline{c}cg`$ after replacement:

$$V^i(p_{\overline{c}})\overline{U}^j(p_c)\frac{\mathrm{\Psi }(0)}{2\sqrt{M}}\widehat{\epsilon }(p)(\widehat{p}+M)\frac{\delta ^{ij}}{\sqrt{3}},$$
(7)
where $`\widehat{\epsilon }(p)=\epsilon _\mu (p)\gamma ^\mu `$, $`\epsilon _\mu (p)`$ is 4-vector of the $`J/\psi `$ polarization, $`{\displaystyle \frac{\delta ^{ij}}{\sqrt{3}}}`$ is the color factor, $`\mathrm{\Psi }(0)`$ is the nonrelativistic meson wave function at the origin. The matrix elements of the process $`\gamma g^{}J/\psi g`$ may be presented as follows:

$$M_i=KC^{ab}\epsilon _\alpha (k_1)\epsilon _\mu ^a(q)\epsilon _\beta ^b(k_2)\epsilon _\nu (p)M_i^{\alpha \beta \mu \nu },$$
(8)
$`M_1^{\alpha \beta \mu \nu }=\text{Tr}\left[\gamma ^\nu (\widehat{p}+M)\gamma ^\alpha {\displaystyle \frac{\widehat{p}_c\widehat{k}_1+m_c}{(p_ck_1)^2m_c^2}}\gamma ^\mu {\displaystyle \frac{\widehat{p}_{\overline{c}}\widehat{k}_2+m_c}{(p_{\overline{c}}+k_2)^2m_c^2}}\gamma ^\beta \right]`$ (9)
$`M_2^{\alpha \beta \mu \nu }=\text{Tr}\left[\gamma ^\nu (\widehat{p}+M)\gamma ^\beta {\displaystyle \frac{\widehat{p}_c+\widehat{k}_2+m_c}{(p_c+k_2)^2m_c^2}}\gamma ^\alpha {\displaystyle \frac{\widehat{q}\widehat{p}_{\overline{c}}+m_c}{(qp_{\overline{c}})^2m_c^2}}\gamma ^\mu \right]`$ (10)
$`M_3^{\alpha \beta \mu \nu }=\text{Tr}\left[\gamma ^\nu (\widehat{p}+M)\gamma ^\alpha {\displaystyle \frac{\widehat{p}_c\widehat{k}_1+m_c}{(p_ck_1)^2m_c^2}}\gamma ^\beta {\displaystyle \frac{\widehat{q}\widehat{p}_{\overline{c}}+m_c}{(qp_{\overline{c}})^2m_c^2}}\gamma ^\mu \right]`$ (11)
$`M_4^{\alpha \beta \mu \nu }=\text{Tr}\left[\gamma ^\nu (\widehat{p}+M)\gamma ^\mu {\displaystyle \frac{\widehat{p}_c\widehat{q}+m_c}{(p_cq)^2m_c^2}}\gamma ^\alpha {\displaystyle \frac{\widehat{p}_{\overline{c}}\widehat{k}_2+m_c}{(k_2+p_{\overline{c}})^2m_c^2}}\gamma ^\beta \right]`$ (12)
$`M_5^{\alpha \beta \mu \nu }=\text{Tr}\left[\gamma ^\nu (\widehat{p}+M)\gamma ^\beta {\displaystyle \frac{\widehat{p}_c+\widehat{k}_2+m_c}{(p_c+k_2)^2m_c^2}}\gamma ^\mu {\displaystyle \frac{\widehat{k}_1\widehat{p}_{\overline{c}}+m_c}{(k_1p_{\overline{c}})^2m_c^2}}\gamma ^\alpha \right]`$ (13)
$`M_6^{\alpha \beta \mu \nu }=\text{Tr}\left[\gamma ^\nu (\widehat{p}+M)\gamma ^\mu {\displaystyle \frac{\widehat{p}_c\widehat{q}+m_c}{(p_cq)^2m_c^2}}\gamma ^\beta {\displaystyle \frac{\widehat{k}_1\widehat{p}_{\overline{c}}+m_c}{(k_1p_{\overline{c}})^2m_c^2}}\gamma ^\alpha \right]`$ (14)
where $`k_1`$ is the 4-momentum of the photon, $`q`$ is 4-momentum of the initial gluon, $`k_2`$ is the 4-momentum of the final gluon,

$$K=e_ceg_s^2\frac{\mathrm{\Psi }(0)}{2\sqrt{M}},C^{ab}=\frac{1}{\sqrt{3}}\text{Tr}[T^aT^b],e_c=\frac{2}{3},e=\sqrt{4\pi \alpha },g_s=\sqrt{4\pi \alpha _s}.$$

The summation on photon, $`J/\psi `$ meson and final gluon polarizations is carried out by covariant formulae:

$`{\displaystyle \underset{spin}{}}\epsilon _\alpha (k_1)\epsilon _\beta (k_1)=g_{\alpha \beta },`$ (15)
$`{\displaystyle \underset{spin}{}}\epsilon _\alpha (k_2)\epsilon _\beta (k_2)=g_{\alpha \beta },`$ (16)
$`{\displaystyle \underset{spin}{}}\epsilon _\mu (p)\epsilon _\nu (p)=g_{\mu \nu }+{\displaystyle \frac{p_\mu p_\nu }{M^2}}.`$ (17)
In case of the initial BFKL gluon we use the following prescription

$$\underset{spin}{}\epsilon _\mu (q)\epsilon _\nu (q)=\frac{q_\mu q_\nu }{𝐪_{}^2}.$$
(18)
For studing $`J/\psi `$ polarized photoproduction we introduce the 4-vector of the longitudinal polarization as follows

$$\epsilon _L^\mu (p)=\frac{p^\mu }{M}\frac{Mp_N^\mu }{(pp_N)}.$$
(19)
In the limit of $`s>>M^2`$ the polarization 4-vector satisfies usual conditions $`(\epsilon _L\epsilon _L)=1`$, $`(\epsilon _Lp)=0`$.

Traditionally for a description of quarkonium photoproduction processes the invariant variable $`z=(pp_N)/(k_1p_N)`$ is used. In the rest frame of the proton one has $`z=E_\psi /E_\gamma `$. In the $`k_{}`$-factorization approach the differential on $`p_{}`$ and $`z`$ cross section of the $`J/\psi `$ photoproduction may be written as follows

$$\frac{d\sigma }{dp_{}^2dz}=\frac{1}{z(1z)}_0^{2\pi }\frac{d\varphi }{2\pi }_0^{\mu ^2}𝑑𝐪_{}^2\mathrm{\Phi }(x,𝐪_{}^2,\mu ^2)\frac{\overline{|M|^2}}{16\pi (xs)^2}.$$
(20)
The numerical calculation is performed in the photon and proton center of mass frame where

$$p_N=\frac{\sqrt{s}}{2}(1,0,0,1),k_1=\frac{\sqrt{s}}{2}(1,0,0,1),$$

$$p=(E,𝐩_{},p_{}),q=(\frac{\sqrt{2}}{2}x,𝐪_{},\frac{\sqrt{s}}{2}x).$$

Here we take into account that the $`J/\psi `$ momentum $`𝐩`$ lies in $`XZ`$ plane and $`(𝐪_{}𝐩_{})=|𝐩_{}||𝐪_{}|\mathrm{cos}(\varphi )`$. The analytical calculation of the $`\overline{|M|^2}`$ is performed with help of REDUCE package and results are saved in the FORTRAN codes as a function of $`\widehat{s}=(k_1+q)^2`$, $`\widehat{t}=(pk_1)^2`$, $`\widehat{u}=(pq)^2`$, $`𝐩_{}^2`$, $`𝐪_{}^2`$ and $`\mathrm{cos}(\varphi )`$. We directly have tested that

$$\underset{𝐪_{}^20}{lim}_0^{2\pi }\frac{d\varphi }{2\pi }\overline{|M|^2}=\overline{|M_{PM}|^2},$$
(21)
where $`𝐩_{}^2={\displaystyle \frac{\widehat{t}\widehat{u}}{\widehat{s}}}`$ in the $`|M|^2`$ and $`\overline{|M_{PM}|^2}`$ is the square of the amplitude in the conventional parton model . In the limit of $`𝐪_{}^2=0`$ from formula (20) it is easy to find the differential cross section in the parton model too:

$$\frac{d\sigma ^{PM}}{dp_{}^2dz}=\frac{\overline{|M_{PM}|^2}xG(x,\mu ^2)}{16\pi (xs)^2z(1z)}.$$
(22)
However, making calculations in the parton model we use formula (20), where integration over $`𝐪_{}^2`$ and $`\varphi `$ is performed numerically, instead of (22). This method fixes common normalization factor for both approaches and gives direct opportunity to study effects connected with virtuality of the initial BFKL gluon in the partonic amplitude.

## 5 Results and discussion

After we fixed selection of the gluon distribution function $`G(x,\mu ^2)`$ there are two parameters only, which values determine the common normalization factor of the cross section under consideration: $`\mathrm{\Psi }(0)`$ and $`m_c`$. The value of the $`J/\psi `$ meson wave function at the origin may be calculated in a potential model or obtained from experimental well known decay width $`\mathrm{\Gamma }(J/\psi \mu ^+\mu ^{})`$. In our calculation we used the following choice $`|\mathrm{\Psi }(0)|^2=0.0876`$ GeV3 as the same as in Ref. . Concerning a charmed quark mass, the situation is not clear up to the end. From one hand, in the nonrelativistic approximation one has $`m_c={\displaystyle \frac{M}{2}}`$, but there are many evidences to take smaller value of a c-quark mass in the amplitude of a hard process, for example $`m_c=1.4`$ GeV. Taking into consideration above mentioned we perform calculations as at $`m_c=1.4`$ GeV as well as at $`m_c=1.55`$ GeV. The cinematic region under consideration is determined by the following conditions: $`0.4<z<0.9`$ and $`p_{}^2>1`$ GeV2, which correspond H1 and ZEUS Collaborations data . We assume that the contribution of the colour octet mechanism is large at the $`z>0.9`$ only. In the region of the small values of $`z<0.2`$ the contribution of the resolved photon processes  as well as the charm excitation processes  may be large too. All of these contributions are not in our consideration.

Figure 3 shows the ratio

$$R_2(\widehat{s},𝐪_{}^2)=\frac{\widehat{\sigma }(\widehat{s},𝐪_{}^2)}{\widehat{\sigma }(\widehat{s},0)},$$
(23)
as a function of $`𝐪_{}^\mathrm{𝟐}`$ at the different $`\widehat{s}`$. As well as it was necessary to expect ratio (23) decrease with growth of $`𝐪_{}^2`$ and the faster than less the value of $`\widehat{s}`$.

Figures 4–6 show our results which were obtained as in the conventional parton model as well as in the $`k_{}`$-factorization approach at two values of a charmed quark mass. The dependence of the results on selection of a hard scale parameter $`\mu `$ is much less than the dependence on selection of a c-quark mass. We put $`\mu ^2=M^2+𝐩_{}^\mathrm{𝟐}`$ in a gluon distribution function and in a running constant $`\alpha _s(\mu ^2)`$.

Figure 4 shows the dependence of the total $`J/\psi `$ photoproduction cross section on $`\sqrt{s}`$. It is visible that the difference between the parton model prediction and the result of $`k_{}`$-factorization approach is much less than between the results obtained at different values of c-quark mass in the both models. At the $`m_c=1.55`$ GeV the obtained cross sections in 1.5 –2 times are less than experimental data , but at the $`m_c=1.4`$ GeV the theoretical curves lie even a shade higher of the experimental points.

There are not contradictions between the theoretical predictions and data for the $`z`$-spectrum of the $`J/\psi `$ mesons. Figure 5 shows that the experimental points lie inside the theoretical corridor as in the parton model as in the $`k_{}`$-factorization approach.

The count of a transverse momentum of the BFKL gluons in the $`k_{}`$-factorization approach results in a flattening of the $`p_{}`$-spectrum of the $`J/\psi `$ as contrasted to by predictions of parton model. For the first time on this effect was indicated in the article , and later in the articles . Figure 6 shows the result of our calculation for the $`p_{}`$-spectrum of the $`J/\psi `$ mesons. Using the $`k_{}`$-factorization approach we have obtained the more hard $`p_{}`$-spectrum of the $`J/\psi `$ than it have been predicted in the LO parton model. It is visible that at large values of $`p_{}`$ only the $`k_{}`$-factorization approach gives correct description of the data . However, it is impossible to consider this visible effect as the direct indication on a nontrivial developments of the small-x physics. In the article  was shown that calculation in the NLO approximation gives more hard $`p_{}`$-spectrum of the $`J/\psi `$ meson too, which one will be agreed the data at the large $`p_{}`$.

As it was mentioned above, the main difference between the $`k_{}`$-factorization approach and the conventional parton model is nontrivial polarization of BFKL gluon. It is obviously, that such spin condition of the initial gluon should result in observed spin effects during birth of the polarized $`J/\psi `$ meson. We have performed calculations for the spin parameter $`\alpha `$, as a function $`z`$ or $`p_{}`$ in the conventional parton model and in the $`k_{}`$-factorization approach :

$$\alpha (z)=\frac{\frac{d\sigma _{tot}}{dz}3\frac{d\sigma _L}{dz}}{\frac{d\sigma _{tot}}{dz}+\frac{d\sigma _L}{dz}},$$
(24)
$$\alpha (p_{})=\frac{\frac{d\sigma _{tot}}{dp_{}}3\frac{d\sigma _L}{dp_{}}}{\frac{d\sigma _{tot}}{dp_{}}+\frac{d\sigma _L}{dp_{}}}$$
(25)
Here $`\sigma _{tot}=\sigma _L+\sigma _T`$ is the total $`J/\psi `$ production cross section, $`\sigma _L`$ is the production cross section for the longitudinal polarized $`J/\psi `$ mesons, $`\sigma _T`$ is the production cross section for the transverse polarized $`J/\psi `$ mesons. The parameter $`\alpha `$ controls the angle distribution for leptons in the decay $`J/\psi l^+l^{}`$ in the $`J/\psi `$ meson rest frame:

$$\frac{d\mathrm{\Gamma }}{d\mathrm{cos}(\theta )}1+\alpha \mathrm{cos}(\theta ).$$
(26)
Figure 7 shows the parameter $`\alpha (z)`$, which is calculated in the parton model (curve 2) and in the $`k_{}`$-factorization approach (curve 1). We see that both curves lie near zero at $`z<0.8`$ and increase at $`z>0.8`$. The large difference between predictions is visible only at $`z>0.9`$ where our consideration is not adequate. Let’s remark, that parameter is gentle depends on mass of a charmed quark and we demonstrate here only outcomes obtained at the $`m_c=1.55`$.

For the parameter $`\alpha (p_{})`$ we have found strongly opposite predictions in the parton model and in the $`k_{}`$-factorization approach, as it is visible in Figure 8. The parton model predicts that $`J/\psi `$ mesons should have transverse polarizations at the large $`p_{}`$ ($`\alpha (p_{})=0.6`$ at $`p_{}^2=20`$ GeV2), but $`k_{}`$-factorization approach predicts that $`J/\psi `$ mesons should be longitudinally polarized ( $`\alpha (p_{})=0.4`$ at $`p_{}^2=20`$ GeV2). Nowadays, a result of NLO parton model calculation in the case of the polarized $`J/\psi `$ meson photoproduction is unknown. It should be an interesting subject of future investigations. If the count of the NLO corrections will not change predictions of the LO parton model for $`\alpha (p_{})`$, the experimental measurement of this spin effect will be a direct signal about BFKL gluon dynamics.

Nowadays, the experimental data on $`J/\mathrm{\Psi }`$ polarization in photoproduction at large $`p_{}`$ are absent. However there are similar data from CDF Collaboration , where $`J/\psi `$ and $`\psi ^{}`$ $`p_{}`$-spectra and polarizations have been measured. Opposite the case of $`J/\psi `$ photoproduction, the hadroproduction dada needs to take into account the large colour-octet contribution in order to explain $`J/\psi `$ and $`\psi ^{}`$ production at Tevatron in the conventional collinear parton model. The relative weight of colour-octet contribution may be smaller if we use $`k_{}`$-factorization approach, as it was shown recently in . The predicted using collinear parton model transverse polarization of $`J/\psi `$ at large $`p_{}`$ is not supported by the CDF data, which can be roughly explain in $`k_{}`$-factorization approach . In conclusion, the number of theoretical uncertainties in the case of $`J/\psi `$ meson hadroproduction is much more than in the case of photoproduction and they need more complicated investigation, that is why the future experimental analysis of $`J/\psi `$ photoproduction at THERA will be clean check of the collinear parton model and $`k_{}`$-factorization approach.

This work has been supported in part by the programme ”Universities of Russia – Basic Researches” under Grant 02.01.003.


format: text

text:
# 1 Introduction

## 1 Introduction

Coarse-grained models are an important tool in theoretical studies of protein folding, for computational as well as conceptual reasons, and have been used to gain insights into the physical principles of folding (for a recent review, see Ref. ). These models are often lattice based. The main advantage of using a discrete conformational space is that exact calculations can be performed for short chains, by exhaustive enumeration of all possible conformations. One model that has been extensively studied this way is the minimalistic hydrophobic/polar HP model of Lau and Dill  on the square lattice. In previous work on this model, all sequences with unique ground states were determined for chains with up to $`N=18`$–20 monomers . Such sequences are called designing; they design their ground state conformations.

In this paper, we show how a few optimizations make it possible to extend these calculations to $`N=25`$, which corresponds to a 4000-fold increase in the number of possible sequence, conformation pairs.222The complete list of all designing $`N25`$ sequences and the corresponding conformations will be made electronically available at http://www.thep.lu.se/complex/wwwserver.html We then use this data set to address the question of how designing sequences differ from random ones statistically.

By analyzing the behavior of block variables, it has been found that designing $`N18`$ HP sequences  as well as real (globular) protein sequences  show negative hydrophobicity correlations. Therefore, one expects to find an increased number of hydrophobic and polar clumps along these chains. In this paper, we show that the average number of clumps is indeed larger for designing HP sequences than for random sequences. In particular, this implies that the earlier finding that designing sequences show negative hydrophobicity correlations remains unaffected when increasing $`N`$ to 25. This provides a non-trivial test of the robustness of this property.

In lattice model studies it is not uncommon to consider only maximally compact conformations, which for $`N=25`$ are confined to a $`5\times 5`$ square. This drastic reduction of conformational degrees of freedom leads to a sharp rise in the number of designing sequences. An interesting question is whether this reduction also affects the statistical properties of designing sequences. To address this question, we repeat the same statistical analysis for sequences that are designing when only maximally compact conformations are used. The results turn out to be qualitatively different in this case. In particular, this means that the agreement with the results obtained for real sequences gets lost when this reduced conformational space is used.

Finally, we also study the character of the folding transition for one of the designing $`N=25`$ sequences, which was selected by an optimization procedure. The thermodynamic behavior of this sequence is studied using Monte Carlo simulations.

The paper is organized as follows. In Section 2 we define the model and describe the algorithm and optimizations used for finding all designing sequences with $`N25`$. Our results are discussed in Section 3, which contains sequence and structure statistics, the statistical analysis of designing sequences, and the thermodynamic study of an optimized sequence. A summary is given in Section 4.

## 2 Enumerating designing sequences

In lattice models of proteins it is common to use a contact potential. This means that the energy that a sequence gets with a certain conformation, is given by what contacts exist in that conformation. That is,

$$E=\underset{i<j}{}C_{ij}U(\sigma _i,\sigma _j)$$
(1)
where the contact map $`C_{ij}`$ is defined as

$$C_{ij}=\{\begin{array}{cc}\hfill 1& \text{if monomers }i,j\text{ are neighbors on the lattice but }|ij|1\hfill \\ \hfill 0& \text{otherwise}\hfill \end{array}$$
(2)
and $`U(\sigma _i,\sigma _j)`$ is the interaction matrix. In the HP model there are two amino acids, hydrophobic (H) and polar (P). The interaction matrix is

$$U=\left(\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill 0\end{array}\right)$$
(3)
with H in the first row and column, so the energy is determined exclusively by the number of HH contacts.

Two conformations that have identical contact sets will have the same energy for all sequences, so they can be represented by a single contact set which is marked as degenerate. What this means is that a large number of conformations can be reduced to a smaller number of contact sets. A contact set that corresponds to a single conformation will be called unique.

### 2.1 Compressing conformational space

To find all designing sequences, we first determine all relevant conformations, which then are combined with sequences in Section 2.2. The conformation search is sometimes simplified by only considering maximally compact conformations, where the number of contacts is maximal. Looking only at those conformations corresponds to shifting the energies by adding a large negative term to all elements of the interaction matrix. An efficient method for enumerating compact conformations was recently proposed by Kloczkowski and Jernigan .

In this paper, we consider all possible self-avoiding walks, and not only those that are maximally compact. The space of all possible self-avoiding walks grows rapidly with $`N`$. Using contact sets rather than self-avoiding walks gives a significant reduction of the conformational space (see Table 2 below), but this reduction is not sufficient for our purposes; memory limitations prevented us from storing the complete list of all possible contact sets for $`N24`$. To be able to go to larger $`N`$, we therefore developed two procedures for reducing the number of contact sets to be stored. These procedures are based on the observation that as long as there are no repulsive forces (that is, as long as the elements of the interaction matrix are all non-positive), it is never energetically disadvantageous to add contacts as long as no existing contacts are broken.

Before discussing these two procedures, it is helpful to introduce some terminology. A contact set will be called maximal if it is not a proper subset of any other (realizable) contact set. A conformation that is designed by at least one sequence is designable. It can be readily seen that every designable self-avoiding walk corresponds to a unique and maximal contact set. Another important class are contact sets that can be safely ignored in the search for designing sequences. Such contact sets will be called redundant. A summary of our contact set terminology can be found in Table 1.

#### 2.1.1 Eliminating contact sets: Step 1

As mentioned above, it was impossible for us to store the complete list of all contact sets for $`N24`$. To circumvent this problem, we used a program that for each self-avoiding walk tries a carefully selected, predefined set of (local) moves. If any of these moves can be performed without destroying any existing contact (new contacts may form), the self-avoiding walk is discarded. All possible self-avoiding walks surviving this test are converted to contact sets. It is important to stress that the moves are chosen so that the resulting reduced list of contact sets has the property that any realizable contact set is a subset of some set in this list. In particular, the move set does not contain the inverse of any of its elements.

The list of contact sets obtained this way was indeed small enough to be stored up to $`N=25`$, but the procedure has the disadvantage that it may eliminate non-redundant contact sets. This is an unwanted property, but the problem is easy to solve. The solution is that once the sequences that have non-degenerate ground states with respect to the non-discarded conformations have been found, each sequence is combined with its conformation, and the chain is tested with the opposite of the tests used to discard conformations. That is, the program tests whether any of the opposite moves can be performed without breaking any of the existing HH contacts. By using the fact that the forces are repulsive, it can be seen that if no such move is possible, then the conformation has to be a unique ground state of this sequence. Hence, by performing this test, one can make sure that no sequence is falsely declared to be designing, even though there are non-redundant contact sets that are missing in the list used.

It is also important to note that all the self-avoiding walks with a given contact set are never discarded if the set is maximal. This means that all the contact sets that correspond to designable conformations are included in the list generated by this procedure. This is important because it implies that this reduced list of contact sets can be used without losing any designing sequence. What could go wrong when using this list is that some non-designing sequences were classified as designing, but this is avoided by performing the test discussed above.

#### 2.1.2 Eliminating contact sets: Step 2

Our second procedure removes contact sets from the list produced by the first procedure. This is done to speed up the calculations. All contact sets that are removed in this second step are redundant.

The procedure relies on the fact that all pair energies are non-positive. To see how that can be used, consider one set of contacts $`A`$ which is a subset of another contact set $`B`$. It is then clear that $`A`$ cannot represent a unique ground state, because for any sequence, $`B`$ has the same or lower energy. The set $`A`$ is nevertheless non-redundant in case $`B`$ is a unique contact set that would be falsely classified as the unique ground state of some sequence if $`A`$ were removed. If, on the other hand, the set $`B`$ is degenerate, then it follows that $`A`$ must be redundant.

Suppose now that $`A`$ is a subset of two other sets $`B`$ and $`C`$. Then, provided that $`C`$ is kept, there can be no sequence such that $`A`$ is needed in order to decide whether or not this sequence has $`B`$ as a unique ground state. This follows immediately from the fact that for any given sequence, both $`B`$ and $`C`$ have energies at least as low as that of $`A`$. Hence, if a contact set is a subset of two or more sets, then it has to be redundant. In particular, this is true if the set is a subset of a subset.

This reasoning gives us the following simple procedure for elimination of redundant contact sets.

* For each set $`A`$, find all sets of which $`A`$ is a subset.
* Keep $`A`$ if

1. no such sets are found, or
    2. one set is found, and this set corresponds to one conformation.
* Otherwise discard $`A`$.

Note that those of the surviving contact sets that meet condition 1 are maximal.

It should be stressed that, because all the contact sets removed by this procedure are redundant, one can still use the test in Section 2.1.1 to make sure that no sequence is falsely classified as designing.

Implementing the rules described above requires some care, since the number of sets grows rapidly with $`N`$. The solution we used is based on storing the sets in a tree, where each node in the tree represents the question of whether or not a certain contact is included. Before building the tree, it is useful to sort the sets in descending order by the number of contacts they contain. This sorting ensures that a set, once added to the tree, never has to be removed, since all sets of which it can be a subset are already added when it is examined. Having access to the tree, it is straightforward to search for supersets of a given set. One starts at the root of the tree and each time the tree branches one has to consider either one of the branches or both. The procedure is memory-consuming for large $`N`$ (see example below), but the CPU time required is relatively modest.

#### 2.1.3 Examples

To get an idea of the sizes of these different lists of contact sets, let us consider the case $`N=18`$. For this $`N`$, there are 5808335 self-avoiding walks and 170670 contact sets. After applying the redundancy test in Section 2.1.2 to the list of all possible contact sets, we are left with 33223 maximal contact sets (condition 1) and 6609 contacts sets with one superset corresponding to one conformation (condition 2). Among the 33223 maximal contact sets, there are 6181 sets representing more than one conformation. Subtracting the degenerate contact sets, we are left with $`332236181=27042`$ contact sets corresponding to possibly designable conformations. Of those 27042 sets, 5660 sets have a total of 6609 listed subsets. These subsets are needed because they may degenerate an otherwise unique ground state corresponding to one of the 5660 sets.

Here, the redundancy test was applied to the complete list of all possible contact sets. Alternatively, we may first use the program described in Section 2.1.1, which for $`N=18`$ generates a list of 51373 contact sets (which contains the 33223 maximal ones). When applying the redundancy test in Section 2.1.2 to this reduced list, we end up with $`33223+449`$ contact sets. Note that with this approach, we find only 449 of the 6609 non-redundant contact sets meeting condition 2. As a result, some conformations are erroneously found to be unique ground states. The test discussed in Section 2.1.1 removes these false unique ground states.

The CPU times required to generate these different lists on a Pentium III 800 MHz were as follows. Generating the list of all possible contact sets, by exhaustive enumeration, took 6 seconds, and reducing this list from 170670 to $`33223+6609`$ contact sets required 7 seconds. The time needed to run the program that generates the list described in Section 2.1.1 was 3 seconds, and reducing this list from 51373 to $`33223+449`$ contact sets took 1 second. The corresponding two numbers for $`N=25`$ were 40 and 60 minutes, respectively. In this case, building the tree used in the redundancy test required 220 megabytes of internal memory.

### 2.2 Searching sequence space

We now turn our attention to the sequence space. For each of the $`2^N`$ sequences we wish to determine what set of contacts gives the lowest energy. If this ground state energy can be accomplished only by a single contact set, and if that set corresponds to a single conformation, the sequence designs that conformation.

#### 2.2.1 Organizing the search

The most straightforward approach to finding all sequences with unique ground states is to go through all the sets of contacts for each sequence, and calculate the energy for each of the sets. By only storing the differences between consecutive sets of contacts, and by representing the sequences and contacts as numbers with one bit per position, the number of operations required for each combination can be kept small.

It is, however, also possible to use a very different approach. Represent each sequence by a binary number, and consider any given set of contacts. Between two consecutive sequences two bits are toggled on the average, which indicates that using information about the previous sequence and its energy will be a lot faster than recalculating the energy from scratch. This approach can be used if the whole sequence space is examined for one contact set at a time. The downside to doing this is that information needs to be stored for every sequence until all the contact sets have been considered.

Neither of the methods described above is bad, but they each contain an optimization that the other does not have. It is desirable to utilize both the similarity of consecutive sequences and the similarity of consecutive contact sets. The solution is to divide the sequence space into a number of blocks of fixed size, and apply the second method to each of those blocks. A block consists of $`2^M`$ sequences that have their first $`NM`$ residues in common. This part of the sequence will be referred to as the fixed part, and the remaining $`M`$ positions make up the variable part. We call a contact active if it connects two H monomers, and note that for each contact there are three possibilities, depending on the position and type of the monomers it connects:

* If both monomers are in the fixed part, the contact gives a contribution of $`1`$ to the energy for all the sequences of this block, if and only if both monomers are H.
* If at least one of the monomers is a P in the fixed part, the contact cannot be active for any sequence in the block.
* If both monomers are in the variable part, or if one of them is in the variable part and the other is an H in the fixed part, whether or not the contact is active depends on the variable part.

#### 2.2.2 The cutoff energy

This leads us to the next optimization, which has to do with the possible ground state energies, and can be seen as a sequence-dependent reduction of the conformational space. Clearly, a non-degenerate ground state with $`N3`$ cannot have an energy of 0. More generally, it is unreasonable that a small number of active contacts should be enough to give an arbitrarily long polymer a unique ground state. To see how this can be used to speed up the calculations, consider some contact set and sequence block, such that there are $`p`$ active contacts in the fixed part and $`q`$ contacts whose state depends on the variable part. If $`(p+q)`$ is larger than some cutoff value $`E_{\mathrm{max}}`$, none of the sequences in this block can have a unique ground state for this contact set. The major problem that arises with this optimization is to know what value to use for $`E_{\mathrm{max}}`$; the algorithm will find only those unique ground states that have energies $`EE_{\mathrm{max}}`$. For $`N20`$ all energies have been considered in the calculations, and it turns out that there are no unique ground states with $`E>4`$ for $`15N20`$ (see Table 3 below). For $`N>20`$ we have not proven that there can be no unique ground states with $`E>4`$, but it seems very reasonable that if there are any, most or all of them would have $`E=3`$. Therefore we have used a cutoff energy $`E_{\mathrm{max}}=3`$ for $`N>20`$. It turns out that for $`20<N25`$ there are no unique ground states with energy $`3`$, and this strongly indicates that $`4`$ is the highest possible ground state energy for any HP chain with $`N15`$.

To illustrate how the computer time required varies with $`E_{\mathrm{max}}`$, let us again take $`N=18`$ as an example. For this $`N`$, the sequence search is found to be about four times faster with $`E_{\mathrm{max}}=4`$ than with $`E_{\mathrm{max}}=1`$. With $`E_{\mathrm{max}}=4`$, the total time needed to find the 6349 designing $`N=18`$ sequences, including the time required to generate the conformations, is a few minutes on a Pentium III 800 MHz.

## 3 Results

Using the algorithm and optimizations described above, it was possible to determine all designing sequences for $`N25`$ within a reasonable amount of time. Previous work has covered $`N18`$ for the normal HP model  and $`N20`$ for shifted HP models . The increase in $`N`$ corresponds roughly to a 100-fold increase in the number of known designing sequences and conformations. This gives better confidence when doing statistics on the designing sequences, and it makes it possible to study how properties of the model depend on protein size.

### 3.1 Sequence and structure statistics

Sequence and structure statistics for $`4N25`$ are summarized in Tables 2 and 3. Column three of Table 2 shows the total number of contact sets, which has been studied before . It was estimated to grow as $`\mu ^N`$ with $`\mu =2.29\pm 0.02`$ for large $`N`$ . The number of maximal contact sets, column four of Table 2, appears to grow exponentially too, but slightly slower. A fit of our data for $`15N25`$ to the form $`\mu ^N`$ yields $`\mu 2.07`$. This growth with $`N`$ is considerably slower than that of the number of self-avoiding walks, for which the best available estimate is $`N^{\gamma 1}\mu ^N`$ with $`\gamma =43/32`$ and $`\mu =2.6381585`$ .

The fifth column of Table 2 shows the number of designing sequences. It turns out that the fraction of designing sequences varies between 2.27 and 2.57% for $`19N25`$, which is in line with previous results for smaller $`N`$ . The last column of Table 2 shows the number of designable conformations. The designability of a conformation is the number of sequences that designs it. From Table 2 it can be seen that the average designability of the conformations that are designable grows with $`N`$ and is $`765147/1073367.1`$ for $`N=25`$.

In Table 3 we show the distribution of ground state energies for different $`N`$, which is crucial for the optimization discussed in Section 2.2.2.

Figure 1 shows three designable $`N=25`$ conformations. The conformation (a) is the most designable one for $`N=25`$, with a designability of 326, whereas conformation (c) is designed by one sequence only.

For comparison, we also determined the sequences that are designing when only maximally compact conformations are used. In this case, it turns out that there are 6181800 designing $`N=25`$ sequences. The corresponding number is 765147 when the full conformational space is used (see Table 2), so the ratio between the number of designing sequences in the maximally compact ensemble and the number of truly designing sequences is $`6181800/7651478.1`$, for $`N=25`$. This ratio has previously been shown  to fluctuate between approximately 4 and 11 for $`N=11`$ through $`N=18`$.

It is worth noting that among the 765147 $`N=25`$ sequences that are designing when the full conformational space is used, there are only 605 sequences that design maximally compact conformations. Furthermore, it turns out that no maximally compact conformation is designed by more than 10 sequences, whereas the most designable conformation, as mentioned above, has a designability of 326. In fact, there are 19360 conformations that are more designable than the most designable one among the maximally compact conformations.

It is interesting to compare these results to those of Shahrezaei and Ejtehadi , who studied a shifted HP model where the interaction matrix is given by $`U(\mathrm{H},\mathrm{H})=2\gamma E_c`$, $`U(\mathrm{H},\mathrm{P})=1E_c`$ and $`U(\mathrm{P},\mathrm{P})=E_c`$. Using the full conformational ensemble, these authors found that the set of highly designable conformations was independent of the parameters $`\gamma `$ and $`E_c`$. In particular, this suggests that the set of highly designable conformations should remain the same when only maximally compact conformations are considered. Our results show that this conclusion does not hold in the original, unshifted HP model.

### 3.2 Statistical properties of designing sequences

In this section, we study the statistical properties of designing sequences by monitoring two different quantities. The first one is the total hydrophobicity $`M`$, defined as

$$M=\underset{i=1}{\overset{N}{}}\frac{1+\sigma _i}{2},$$
(4)
where $`\sigma _i=1`$ (H) or $`\sigma _i=1`$ (P). Our second quantity is the number $`j`$ of hydrophobic and polar clumps along the chain , which can be written as

$$j=\frac{N+1}{2}\frac{1}{2}\underset{i=1}{\overset{N1}{}}\sigma _i\sigma _{i+1}.$$
(5)
A similar analysis has previously been performed for $`N18`$ . By analyzing the fluctuations of block variables, evidence was found for negative $`\sigma _i,\sigma _j`$ correlations. Therefore, we expect the average $`j`$ to be larger for designing sequences than for random ones.

In Figure 2a we show the relative abundance of hydrophobic amino acids, $`M/N`$, as a function of $`N`$, where $``$ denotes an average for fixed $`N`$. For the sequences that are designing when the full conformational space is used, we see that the $`N`$ dependence of $`M/N`$ is quite weak if $`N`$ is not too small, which confirms a trend seen earlier . Furthermore, we see that these sequences, as expected, are more hydrophobic than those that are designing when only maximally compact conformations are used. Finally, it can also be seen that the sequences that design maximally compact conformations differ greatly from those that are designing when only such conformations are used. Figure 2b shows the frequency of different $`M`$ for $`N=25`$.

In Figure 3 we show the results of our clump analysis for $`N=25`$. The average number of clumps for fixed $`M`$ (and $`N`$), $`j_M`$, is indeed found to be larger for designing sequences than for random ones, which is in nice qualitative agreement with previous results for real protein sequences and model sequences with smaller $`N`$ . Sequences that are designing when the interaction energies are shifted so far that only the 1081 maximally compact conformations need to be considered, have, by contrast, a $`j_M`$ close to that of random sequences. Hence, the results obtained from studying only maximally compact conformations seem to be of less relevance, with respect to the applicability to real proteins, than the results obtained when considering all conformations.

Finally, we note that Buchler and Goldstein have compared various lattice models and found arguments against the use of only two letters , as in the HP model. These findings were all based on calculations for maximally compact conformations. However, it is not known to what extent the results of such calculations remain valid when the full set of conformations is used. Our results show that, in the HP model, both the set of highly designable conformations and the statistical properties of designing sequences depend strongly on which of the two conformational ensembles is used.

### 3.3 The character of the folding transition

The model studied in this paper has the important feature that there exists a significant number of designing sequences (this is not true on the triangular  and cubic  lattices), and that the corresponding conformations tend to show protein-like regularities . However, as in other two-dimensional models, the folding transition is not protein-like in character for the typical sequence; the folding process is not cooperative enough . On the other hand, the folding behavior is, at least to some extent, sequence dependent, and therefore we decided to look into the thermodynamic behavior of a carefully chosen sequence.

This sequence was obtained by applying a Monte Carlo-based sequence design algorithm  to the 326 sequences that design the most designable $`N=25`$ conformation (see Figure 1a). The design algorithm maximizes the stability of a given conformation with respect to sequence at a fixed non-zero temperature. The sequence we obtained by using this method is shown in Figure 1a. Subsequently, this sequence was subjected to Monte Carlo simulations at different temperatures. In Figure 4 we show the temperature-dependence of the specific heat, which is found to exhibit a pronounced peak. Also shown is the distribution of energy at $`kT=0.479`$, which is just above the specific-heat maximum. The energy distribution has one peak corresponding to the ground state, at $`E=13`$, and another, broader peak centered at $`E8`$. The coexistence of these states implies that the folding transition is much more cooperative for this sequence than for the typical sequence in this model.

## 4 Summary

By greatly reducing the conformational space, and by carefully optimizing the sequence space exploration, we were able to decrease the time needed to exhaustively search for designing sequences with $`N=18`$ roughly thousandfold compared to the most naïve methods. It became feasible to find all designing sequences for $`N`$ as large as 25 using only a small number of workstations. The results obtained by doing this were used to look at the statistical properties of designing sequences. We found that the average number of hydrophobic and polar clumps along the chains is larger for designing sequences than for random ones. In particular, this means that the finding that designing HP sequences, like real enzymes, show negative hydrophobicity correlations  remains unaffected when increasing $`N`$ from 18 to 25. By contrast, qualitatively different results were obtained when discarding conformations that are not maximally compact. This is of interest because restrictions to compact structures are common in both lattice model studies and determinations of statistical potentials from known protein structures. Finally, we saw an example of a folding behavior that is more cooperative than for the typical sequence in this model.

### Acknowledgements

We would like to thank Erik Sandelin for fruitful discussions, Björn Samuelsson for generously providing the program discussed in Section 2.1.1, and two anonymous referees for useful remarks. This work was in part supported by the Swedish Foundation for Strategic Research.


text:
# Minimum Bias Legacy of Search Results

## 1 INTRODUCTION AND RATIONALE

This contribution starts from an observation which I hope is shared by many other particle phycisists, and scientists in general. Science, and particle physics in the specific case, should be considered a global activity, with individuals cooperating at different levels. Certainly, those who have played relevant role in the research should be acknowledged and, in special exceptional cases, rewarded. But the outcome of the research should be considered, finally, intellectual property of the community which has allowed the research. This point of view implies that a result should be presented in such a way that the pieces of information provided by an experiment should be possibly exploited at best by other scientists. This does not mean that I agree with populistic “do it yourself analysis”, of which I have heard, as only those who know in depth detector performance and phenomenology can provide sensible results. Therefore, the problem is only that of making the result public in such a way that the data summaries can be used later in the most efficient way.

As far as searches are concerned, the experiment could report a spectacular effect which fits easily in the rest of the knowledge (“the net of beliefs”) and all members of the physics community believe to a discovery. In other, not controversial cases, there could be no hint at all, as it usually happens in most searches. But there could be the case in which there are some indications. The scientists are then in the doubt of making the claim, thus risking their reputation, or keeping quite, thus loosing the chance. “The experiment was inconclusive, and we had to use statistics”, said one once. But statistical methods are “inconclusive” by definition, if one interprets them rigorously: even 100 observed events over 1 expected background does not necesseraly implies that they are due to signal, though we all tend to believe so…. Therefore, it is important to distinguish between what we rationally believe from what the empirical facts teach us.

Besides personal preference for prudence or to risk, it is a matter of fact that the three situations sketched above (discovery, no effect and dubious hints) are not sharply separated. It is, therefore, important to understand what is the “optimal presentation”, following some desiderata which can be summarised saying that the information contained in the experimental data should be presented in the most powerful and unbiased way.

* The result should not depend on whether one believes that there is no effect (i.e. giving a limit), or that something has been found (a claim with an associated e.g. mass, cross-section, etc).
* The pieces of evidence coming from different experiments should be combined in the most efficient way:

+ If many independent data sets each provide a little evidence in favor of the searched-for signal, the combination of all data should enhance that hypothesis.
    + If the indications are incoherent, their combination should provide a stronger constraint against that hypothesis.

These desiderata are not satisfied by the usual way of reporting results of searches with “CL limits”. There are also other kinds of problems with “CL limits”, some of which will be illustrated in next section. In particular, there is a problem of interpretation, which is the main cause of false alarms in the past and of the spread of misleading information to the general public (consequently throwing bad shadows to the reputation of particle physicists).

## 2 PROBLEMS WITH LIMITS. EXAMPLES FROM THE CONTACT INTERACTION SEARCHES

The first time I was confronted with the problems of stating limits was in occasion of an overall analysis of results on electron compositeness. Analysing all results available at that time I found the situation not really satisfactory. In many cases the operative procedure used to evaluate limits was not described, and sometimes the numerical values of published limits seemed to disagree with differential cross-section to which they referred. I considered this problem particulary serious because the procedure which I understood had most consensus at that time had technical intrinsic problems in about 10% of the results, as it will be explained in a while. But out of the many douzens of results (many reactions $`\times `$ many coupling $`\times `$ many experiments), no technical anomaly was reported.

Other bad feature was the difficulty (or impossibility) to combine consistently the limits, with the consequent attitude to quote only the larger limit, that often was nothing but the largest statistical fluctuation, and not due to the experiment having the higher resolution power for that channel.

### 2.1 Costraining $`\mathrm{\Lambda }`$ to infinity

In a contact interaction analysis, no effect is obtained when the scale $`\mathrm{\Lambda }`$ is infinity. Therefore, one looks in the data for a compatibility of the measured $`\mathrm{\Lambda }`$ with infinity, thus providing a lower limit. Needless to say, a MINUIT minimization around infinity is not trivial at all, not to speak about the interpretation of results. As a consequence, I found that, at that times, there was much of kitchen to get a number to quote as lower limit out of the MINUIT printout. The worst case was an experiment which had got an upper limit of 1.3 TeV for a $`\mathrm{\Lambda }`$ in a certain coupling (this was the number I got reanalysing their data, due to a overfluctuation), but published exactly that number as lower limit.

This technical problem can be overcome working with the conjugate quantity $`ϵ=1/\mathrm{\Lambda }^2`$, which should come about zero in case of no effect. The choice of the second power of $`\mathrm{\Lambda }`$ is due to the observation that the new terms in the cross-section come with that weigth. As a consequence, any additive “noise” make $`ϵ`$ Gaussian distributed around the “true” value. Reporting the result on $`ϵ`$ is certainly a good empirical practice, also because the results can be, in most cases, easily combined with the standard weighted average formula. Moreover, the standard deviation of $`ϵ`$ is an intrinsic property of the experiment, a kind of “resolution power” depending on quality of the detector, luminosity, sensitivity to a particular reaction and level of background. And, in fact, my proposal was to use simply $`\sigma _ϵ`$ as measure of the resolution power of an individual experiment or of a combination of experiments.

The problem remains if we insist to report a CL limit. Calling $`ϵ_{}`$ the best fit value, and $`\sigma _ϵ`$ the standard deviation (the latter is related to the curvature - “width” - of the $`\chi ^2`$ parabola) the standard 95% lower limit for $`\mathrm{\Lambda }^\pm `$ is given by

$`\mathrm{\Lambda }^\pm `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1.64\sigma _ϵ\pm ϵ_{}}}}.`$ (1)
This implies that, if $`|ϵ_{}|`$ is approximately equal to $`1.64\sigma _ϵ`$, then either limit becomes very large. If $`|ϵ_{}|>1.64\sigma _ϵ`$ (10% of cases) there are problems with the standard procedure, including the fact that for either sign of coupling there should be an evidence, yielding an upper limit for $`\mathrm{\Lambda }`$. As a matter of fact, unwanted results are tamed using “prescriptions”, including “Bayesian prescriptions” – a contradiction in terms, in my view.

Besides the details of procedures, it is clear that this kind of approach can produce large limits just as statistical fluctuations, limits which have nothing to do with the effective resolution power of the experiment for a particular channel. It is also a matter of fact that the difficulty to combine limits obtained in such a way resulted in a general tendency to believe that a larger limit would make the lower ones obsolete, though the latter might result from higher resolution experiments.

## 3 MISINTERPRETATION OF “C.L. RESULTS”

A second problem with standard “C.L.” is their interpretation, as it resulted from a survey I made in 1998 . It came out, in fact, that most particle physicists believed that, e.g., a 95% C.L. lower bound on the Higgs mass meant that “the mass of the Higgs, provided it existed, has 95% chance to be above the limit, and 5% chance to be below”. So, the problem is: do we understand each other? Do we communicate the correct information to the general public of tax payer which financed our research?

Nowadays it is quite understood by many people that such probabilistic statements are erroneous and misleading, but there are still people which use such wrong statements based on frequentistic C.L.’s to report the results to the general public. So, for example, the 2000 hint of a 115 GeV Higgs was reported by a spokesperson of the LEP experiments saying that “It is a 2.6 sigma effect. So there’s still a 6 in 1000 chance that what we are seeing are background events, rather than the Higgs”. So, basically the problem persists, since also those who agree on what CL’s should not mean, but still stick on the frequentistic approach, have difficulty in explaining what those results do mean. In my opinion, the very reason of this matter of fact is that, from a genuine inferential point of view – which is what matters in Science  – frequentistic CL’s have no meaning,111Here is the most recent example I know of the series of misleading press releases in which frequentistic numbers (i.e. without a precise inferential meaning) are translated into probabilistic statements: “The experimenters reported a three-sigma discrepancy in sin2qW, which translates to a 99.75 percent probability that the neutrinos are not behaving like other particles.”  and only in some case, under some hidden hypotheses, they can be related to sensible probabilistic statements. This is the reason why I insist in my position that “the solution to the problem of confidence limits begins with removing the expression itself”.

## 4 BAYESIAN WAY OUT

In my opinion, the solution is not to propose a new prescription with the hope that it will be adopted by some influential friends, but rather to change radically the attitude. This means we should review critically how our beliefs about physics quantities or laws of nature are modified by empirical observation and, once we have understood this scheme, we should stick to it, using logic instead of tradition and/or authority. The most powerful tool to learn from data is – and this is not merely my opinion – a theorem, which states how our beliefs must be updated by new pieces of information. This is the undiscussed role of Bayes’ theorem, on which there is agreement also by those who disagree on the use of beliefs in physics (but this is a different story, on which there is little to discuss, since Science is nothing but a collection of beliefs based on empirical facts and intellectual constructions…).

Referring to a physics quantity of unknown value $`\mu `$, the Bayes’ theorem can be shortly expressed as

$$f(\mu |\text{data,\hspace{0.17em}I})f(\text{data}|\mu ,\text{I})\times f_{}(\mu |\text{I}),$$
(2)
where the three ingredients of the Bayesian inference are, from left to right, final pdf, likelihood and prior. Note that $`f()`$ stands, in this approach, for the the pdf expressing the relative beliefs.

It is self-evident that the likelihood has the role of re-shaping (re-weighting) the beliefs. This is a fundamental ingredient of inference, but not the only one. Usually it is considered “more objective” than priors, because it is easier to agree on the response of a detector than on purely speculative values of $`\mu `$. Nevertheless, $`f(\text{data}|\mu ,\text{I})`$ is a probability too, and, as such, tells how much we believe that some data could be observed, for every hypothesis on $`\mu `$.

At this point, the usual objection is that “there are the priors”. In my opinion, there is no real problem if we understand on the logically crucial, often practically irrelevant role of priors. Without them, it would be impossible to make the “probabilistic inversion” from $`f(\text{data}|\mu ,\text{I})`$ to $`f(\mu |\text{data,I})`$, which is the essence of the probabilistic inference. As far as its practical role, it depends on the different problems, and, more specifically, on the shape of the likelihood. This is the reason why I like to classify the inferential problems into closed and open likelihood.

### 4.1 Closed likelihood

The easy case is when, for a given set of data, the likelihood function, $`(\mu )=f(\text{data}|\mu ,I)`$ is closed in both sides, i.e. $`(\mu )0`$ when $`\mu `$ tends to the extremes of its physical region of definition (usually $`\mathrm{}<\mu <+\mathrm{}`$ or $`0<\mu <+\mathrm{}`$). The best understood example of this case is a Gaussian response of the apparatus. Under this condition and, in particular, when the width of $`(\mu )`$ is narrow with respect to any rational prior knowledge (the so Savage’s “precise experiment” situation), the conclusions do not depend practically (in the sense how the result is perceived) from any reasonable model of the prior. Obviously, there could be someone who has a pure mathematical approach to the problem and propose a fancy mathematical expression for the prior. But I do not think this is a problem for physicists (and, anyhow, I am very curious to meet a defender of such fancy prior to check how ready he/she is to defend it with a suitable combination of bets…).

### 4.2 Open likelihood

The question becomes really problematic when the likelihood is open in either side, i.e. $`(\mu )`$ goes to a constant at the edges of the physical region (if it is open in both sides the experiment has little relevance, unless it helps in enhancing a certain region selected by other experiment having closed likelihood). For example, in the Higgs search of the kind performed at LEP, an infinite mass produces a non zero pdf of observing what we do observe. This “bad” feature is due to background: whatever we observe, we are never certain that it is not due to background alone. If this is the case, we are never able to normalize the final pdf, unless we force the normalization with a prior (or other empirical evidence, like radiative correction in the specific problem of the Higgs , which forbids (or strongly inhibits) high masses. It follows that the sensitivity to the prior is usually so critical (for an extensive, introductory discussion of problem and proposed solution, see Ref. ) that one should refrain from publishing, and spread to the general public, probabilistic limits, or limit which are usually mis-interpreted as such (this is what happens practically always with CL’s limits).

In the case of open likelihood, only the likelihood itself should be reported. This information can be easily combined with that coming from other experiments, and satisfies the desiderata of Section 1. Alternatively, one could report the log-likelihood or the rescaled R function proposed in Refs. . This functions can be parametrized in a suitable way, and perhaps stored in web sites in the case of searches having the likelihood depending of several quantities. An example of published parametrized log-likelihood can be found in Ref. .

An alternative to probabilistic, or CL bounds, to quantify with a single number (or with contour lines in the case of multidimensional analysis) is to report the sensitivity bound. This quantity should give, though roughly, the edge after which the experiment looses sensitivity to the search .

## 5 CONCLUSION: A PLEA FROM A EU TAX PAYER

In conclusion, my plea to LEP and SLC colleagues, as well as to all other physicists involved in searches, is to report likelihoods for all searched channel, so that the effort of all community can be used at best for all future analyses.

I would like to remember that publishing likelihoods was, actually, a point about which there was unanimous agreement in the January 2000 workshop on confidence limits held at CERN, as explicitly asked by Massimo Corradi and recorded by Louis Lyons . It was, indeed, the only generally agreed conclusion of the workshop. 222In my point of view there were also another objective, though implicit, conclusion of that workshop. The fact itself that many “specialists”, all educated on the same books of statistics, meet to try to solve self-evident problems, absurdities and contradictions originated by the approach followed by those books, means that the basic concepts of those books should be critically reviewed, and that the gurus of that school of thought, responsible of such confusion, should simply keep silent for the rest of their life.

On the light of the well known properties of the likelihood and on the agreement during that meeting, attended by representative by all major experiments, it is surprising to realize that reporting likelihoods has not become the standard yet. Some say that I am too naïve, and that this will never happen, essentially for two reasons. First, publishing a likelihood is much more committing than just giving a single “95% CL limit” in the region where other experiments report similar limits. Second, there are people specialized in combination of results which do use likelihoods, but are afraid to loose this privileged position if likelihoods are available to any student. I hope it is not so, and that it is only due to some inertia. Therefore I am still optimistic (or naïve).


text:
# Superconductivity in High 𝑇_𝑐 Cuprates: The Cause is No Longer A Mystery

## Abstract

I discuss various direct calculations of the properties of the one-band Hubbard model on a square lattice and conclude that these properties sufficiently resemble those of the cuprate superconductors that no more complicated interactions are necessary to cause high Tc superconductivity.. In particular, I discuss phonon effects and conclude that these may be effective in reducing Tc and the gap in electron-doped materials.

The standard preamble for all kinds of papers on theory (or experiment, for that matter) in the field of high Tc usually contains the phrase “since there is no consensus on the cause of high Tc superconductivity” or words to that effect, and often proceeds to justify thereby yet another implausible conjecture as to some aspect of the phenomenon. Of course one man’s consensus is another man’s wild disagreement, and you can easily find those who do not agree that there is a consensus on relativity, the quantum theory, or the Standard Model, by doing a simple web search. Not, of course , to mention theories of evolution or of the big bang, which are unpopular in whole states and most legislatures. As a solid stater I thank my stars that the band theory of solids and the BCS model are such obscure targets that they do not garner this kind of disapprobation, although each had its powerful though somewhat irrational opponents in the past–names like Slater and Wigner, even, for the latter.

I would agree that there is or should be no real agreement as to the cause of much of the peculiar phenomenology of the cuprates. Each of us would-be experts has his favorite list of really puzzling questions about them. My own favorites are the peculiar insensitivity of Tc to disorder, and the strange transport properties in the normal state. Another mystery is the role of interlayer interaction effects in all kinds of ways, although we can now be reasonably certain that they do not cause superconductivity. Yet it is time , I feel, that it should no longer be legitimate to doubt the “first cause”, namely the minimal underlying model which produces, among a bewildering welter of other effects, superconductivity at still unprecedentedly high temperatures. I can speak with a certain lack of bias on this matter because the interlayer tunneling theory which I advocated for five years, and wrote a book about, turns out to be one of those which must be consigned to the dustbin. (Although before I was wrong, I was right, at least partially, as to the correct “first cause”.)

My confidence is based on two things. First is the rough agreement between a number of different simulations, extrapolations, and variational calculations, and the great similarity between these calculations and the actual experimental data, using a “bare-bones” model for the electrons and their interactions, a model for which there is much independent verification. The second is the experimental data themselves, which in many particulars seem to be trying to tell us the nature of the phenomena.

All of these successful calculations are based on the same simple model for the underlying physics. This model is the one-band Hubbard model, or, as can be shown to be essentially equivalent in the appropriate strong-coupling limit, the $`tJ`$ model. It is interesting that there seems to be no need to include electron-phonon coupling; I will discuss the reason and possible effects of the phonons later in this article. The appropriate Hamiltonian, then is

$$H=(i|t|j)c(i)c(j)+Un(i,+)n(i,)$$
$`(1)`$
$$H(tJ)=P(i|t|j)c(i)c(j)P+J(ij)S(i)S(j)$$
$`(2)`$
Here $`P`$ is the projection operator which removes double occupancy of any site. (2) is meaningless without the projection. $`+`$ and $``$ of course are the spin indices.

The most direct calculations with (2) are the quantum Monte-Carlo simulations of Scalapino and his group, and using a somewhat different method, those of Sorella and collaborators. There is considerable controversy between these two groups as to the details of the results; however, what strikes the disinterested observer such as myself is not these differences but the essential similarity of the two. The low-temperature phase for low doping is of course antiferromagnetic ((2) reduces to the Heisenberg model in that case) but with increasing doping one encounters first an inhomogeneous (probably striped) phase and then, at approximately the right temperature and hole density, d-wave superconductivity. The parameters of (1) or (2) are not arbitrary–experiments of many types give us $`J`$, and ARPES measurements give us a good estimate of the band structure at high doping, which is t–so that there is only a little freedom to manipulate the various components of t. There are some confirmatory simulations, particularly on ladder systems, which again give us a great deal of understanding of why d-wave superconductivity arises; there are also strong indications of a pseudogap phenomenon in the appropriate place. I am no real expert on this extremely tricky field, and in fact have been a pessimist about its ability to arrive at any decisive results, but at this point I have to confess myself impressed. Of course, the trouble with direct calculations is both that they do not really lead to deep understanding, nor can they say much about such things as excitation spectra and transport properties, but what they can say is: when you start out with this physics, you end up with the observed results. I am not at all an expert in the rather involved and subtle reasoning that goes into these calculations, so will not present them in detail, but there is no question that these are the most extensive direct calculations in existence.

A second approach has been taken by Randeria and Trivedi, by Masao Ogata, and to some extent by Sorella’s group. This is to take seriously the original ansatz of the RVB theory that the ground state wave function could be approximated by Gutzwiller projection of an appropriate product wave function of BCS type, and to determine variationally the energy gap and the chemical potential of the wave function to be projected (which, they point out, need not be the physical chemical potential) by calculating the energy and varying it. The Randeria group also improves the wave function somewhat by Monte Carlo methods. Again, as before the results resemble very strongly the experimental data insofar as these can be determined from the calculation; in particular, the prediction of d-wave superconductivity for an adequately doped sample is very robust. Ogata finds that for underdoped concentrations he can find mixtures of antiferromagnetism and superconductivity, as is observed in some cases.

There is still a third method of direct calculation which has been applied, starting with early work by Georges and Yedidia, and carried to an amazing degree of refinement by Bill Putikka in recent years. This is direct expansion of the power series in 1/T for the partition function, and extrapolation of the results with Pade approximants. This method has played a great role in the past in estimating critical behavior near phase transitions, but hasn’t previously been much used in quantum many-body systems. It has been very useful in disposing of red herrings which have appeared in various mean field approximations–showing that the $`tJ`$ Hamiltonian does not lead to phase separation until $`J`$ is unphysically large, for instance; and also identifying the very large $`U`$, low doping “Nagaoka” ferromagnetic region. But recently Puttika has carried his series out to ten to twelve orders and has been able to calculate n(k) with sufficient accuracy to spot the disappearance of the Fermi surface near the ($`\pi `$, 0) point which is characteristic of the pseudogap regime. The agreement of his results with those of ARPES at the same energy resolution is remarkable.

The model is of course not the real substance. Why does it (the model) work so well? In the first place, the question of layer interactions: these are remarkably weak, which we can take as an empirical fact on the basis that the Tl one-layer cuprate is experimentally seen to exhibit all the phenomena that the others do, though proven conclusively by Moler et. al to have very weak interlayer interaction. The same may be said about the gate-doped $`CaCuO2`$ sample of Kloc, Schon et. al, where there is only one doped layer, but perhaps total acceptance of these results should await confirmation on more samples and in other laboratories. The simplification to the simple one-band Hubbard model is often justified by the discussion of Zhang and Rice, which is more or less correct, but I prefer the use of projective canonical transformations a la O K Anderson; and actually the strongest argument was the very early exact calculations on small clusters by Michael Schluter et. al who showed that an effective one-band Hubbard model worked very well indeed.

There are two physical aspects which are not in the model and could be important. One is long-range Coulomb forces–one assumes effective screening as in a normal metal, and the only real justification is that it works. Some experimental data bear on this: we have good data on interlayer plasmons which show that the c-axis is a fairly insulating direction, actually, (except for supercurrents) but as far as we can see the intralayer plasmons are indeed high-frequency and may screen reasonably well.

The question which keeps coming up is phonons: why are there almost no relevant phonon phenomena? The answer to this, I would speculate, is two-fold. I believe that it is almost on the level of the calculational results above that the nature of the pseudogap is that it is a pairing phenomenon in the spin sector–i e, that the pseudogap region can be seen as an RVB with d-wave pairing. But the spin sector is not coupled to phonons in lowest order, because a displacement of the local potential does not break the spin degeneracy. In the second term of (2) the electron-phonon interaction couples only to J. It is worth confirming this experimentally, but it is known that in heavy-electron materials the Kondo spins do not couple to phonons. Therefore, the pseudogap is not affected one way or another by phonons (or ordinary impurities) In the point of view that I have called “RVB redux” the superconducting gap is caused by the kinetic energy cost of the opening of the pseudogap, which can only be restored by pair hopping via the anomalous “josephson” terms in the kinetic energy; thus at least for the higher doping regime the SC gap follows the pseudogap.

This is rather a speculative argument, and I am going immediately to say that in fact there are phonon effects on Tc. For most phonons these are reasonably small because, as has been exhaustively proven for ordinary superconductors, the electron-electron interaction due to phonons is very local and acts primarily between Wannier functions on the same site. In our Hubbard model, it acts simply as a modulation of U. But of course the reason why we have a d-wave is that it has zero amplitude at the origin and hence avoids the repulsive U, and neither U nor phonons are effective. But I think there is one phonon which can be expected to couple rather strongly to our d-electrons, and much more strongly to electrons than to holes (note that in our projective theory (2) there is no particle-hole symmetry.) This is the phonon which represents the Jahn-Teller displacement which breaks the $`d_{x^2y^2}`$ vs $`d_{z^2}`$ degeneracy (see Fig 1). An electron site has no ($`x^2y^2`$) hole and thus tends to relax the Jahn-Teller distortion, while a hole site simply adds to the distortion a little. As is easily seen from the figure, the phonon which couples best is at wave vector $`\pi ,\pi `$ and hence couples the peaks in the gap function $`\pi ,0>0,\pi `$ This phonon has been the subject of considerable discussion recently but rather off the point. What it will do is to cause a repulsive interaction for the d-wave gap and hence lower the superconducting Tc for electrons relative to holes., an effect which is observed but has seemed rather puzzling. It is noteworthy that in a couple of recent measurements (see, for instance, ) the pseudogap and T\* are about the same for electron and hole dopings at the same level, but $`T_c`$ is much reduced.

In fact, if the effect is rather localized in k-space we might expect the two gaps to behave as in Fig. 2, with the $`\pi ,0`$ peak depressed for the SC gap relative to the general point along the Fermi surface. This seems to explain the ARPES observations that the peak of the gap function is not at the $`0,\pi `$ points but rather there is a local minimum there, as we have sketched in Fig. 2 (not, of course, to scale)

One implication of this picture is that we would expect an isotope effect for the electron-doped materials, but a negative one. A crude calculation, assuming without any justification that the BCS gap equation is valid, would suggest

$$(\mathrm{}nT_c)/(\mathrm{}nfopt)=\mu /(J\mu )$$
$`(3)`$
where $`fopt`$ is the optical phonon frequency, $`\mu `$ the phonon coupling constant, and $`J`$ the effective antiferromagnetic exchange . (3) might be valid in the overdoped regime. For reasons mentioned but not emphasized in my book, (which include a great deal of experience) I am not very convinced of the accuracy or relevance of isotope effect measurements, but perhaps the gate-doping possibility, which avoids the necessity of making a new sample for every measurement, is a new opportunity.

In conclusion, my point here is that there is a great deal of consensus on the model which underlies high $`T_c`$ cuprate superconductivity, and there ought to be more: I think we have proved our point. But there is much more to be done and specifically, for instance, we are so far unable to give a closed-form gap equation.

FIGURE CAPTIONS

1. The phonon which is likely to interact strongly with electrons around the 0, points. Its displacements are those of the Jahn-Teller distortion caused by an $`x^2y^2`$ hole, so that it might be expected to be particularly soft for electron doping.
2. The reduction in the superconducting gap which might result from the repulsive interaction caused by the optical phonon of Fig. 1.


text:
# Curie Temperature Trends in (III,Mn)V Ferromagnetic Semiconductors

\[

## Abstract

We present a theoretical survey of ferromagnetic transition temperatures in cubic (III,Mn)V semiconductors based on a model with $`S=5/2`$ local moments exchange-coupled to itinerant holes in the host semiconductor valence band. Starting from the simplest mean-field theory of this model, we estimate the $`T_c`$ enhancement due to exchange and correlation in the itinerant-hole system, and the $`T_c`$ suppression due to collective fluctuations of the ordered moments. We show that high critical temperatures in these ferromagnetic semiconductors require both the large magnetic susceptibility contribution from the valence band’s heavy holes and the large spin stiffness that results from strong spin-orbit coupling in the valence band. Our calculations demonstrate that this model for the ferromagnetism of these systems is fully consistent with the room-temperature ferromagnetism reported for Mn doped nitrides.

\]

The 1992 discovery of hole-mediated ferromagnetic order in (In,Mn)As has motivated research on Mn doped GaAs and other III-V host materials. Ferromagnetic transition temperatures in excess of 100 K and long spin-coherence times in GaAs have fueled hopes that a new magnetic medium is emerging that could open new pathways for information processing and storage technologies. Recently, confirmation of the room-temperature ferromagnetism predicted in (Ga,Mn)N has added to interest in this class of materials.

When substituted on the cation site of a III-V semiconductor, it is generally accepted that Mn acts as an acceptor leaving a Mn2+ half-filled d-shell with angular momentum $`L=0`$ and spin $`S=5/2`$. There is also experimental evidence that ferromagnetism occurs in these materials because of interactions between Mn local moments that are mediated by holes in the semiconductor valence band. The local moments and the itinerant holes are coupled in this approach by an antiferromagnetic exchange interaction. The participation of itinerant holes in the ferromagnetism of these diluted magnetic semiconductors adds to their richness, leading to strong magnetotransport effects that might have important applications and, of particular interest to us here, to transition temperatures that are sensitive to the details of the host band structure.

The description of ordered states in (III,Mn)V systems is greatly simplified by a virtual-crystal type approximation in which the Mn ion distribution is replaced by a continuum with the same impurity density. Disorder effects associated with randomness in the Mn ion sites can then be treated perturbatively, when necessary. (For example randomness in the Mn distribution likely limits the conductivity at low-temperatures.) This strategy will fail in the limit of dilute Mn ions and also when the exchange interaction between band and Mn spins is too strong. It does, however, seem to be reliable in the limit of principal interest, that of high Mn densities and high critical temperatures, where the holes are metallic and their interaction with Mn acceptors will be effectively screened.

There is at present considerable activity directed toward the growth of many different (III,V) materials containing Mn. In this paper we present and discuss theoretical predictions for their ferromagnetic critical temperatures based on the continuum model. We go beyond the standard mean-field theory of this model by accounting for the role of Coulomb interactions holes in the valence band, which enhances the critical temperature, and for correlations in Mn ion orientations which reduce the energetic cost of unaligned spin configurations and therefore lower the critical temperature.

When disorder and hole correlations are neglected, the mean-field transition temperature of this model is given by

$$k_BT_c=\frac{N_{\mathrm{Mn}}S(S+1)}{3}\frac{J_{\mathrm{pd}}^2\chi _f}{(g\mu _B)^2},$$
(1)
where $`N_{\mathrm{Mn}}=4x/a_{\mathrm{lc}}^3`$ is the Mn density in MnxIII1-xV zinc-blende semiconductors with a lattice constant $`a_{\mathrm{lc}}`$, $`J_{\mathrm{pd}}`$ is the localized spin – itinerant spin exchange coupling constant, and $`\chi _f`$ is the itinerant hole magnetic susceptibility. To understand the qualitative implications of this $`T_c`$-equation (1), we first consider the case of a model itinerant electron system with a single spin-split band and an effective mass $`m^{}`$. The kinetic energy of the band holes gives a contribution to the susceptibility:

$$\frac{\chi _f^{\mathrm{kin}}}{(g\mu _B)^2}=\frac{m^{}k_F}{4\pi ^2\mathrm{}^2},$$
(2)
where $`k_F`$ is the Fermi wavevector. For weak interactions, the exchange energy of the spin-polarized parabolic-band model adds a contribution:

$$\frac{\chi _f^{\mathrm{ex}}}{(g\mu _B)^2}=\frac{e^2(m^{})^2}{4\pi ^3\epsilon \mathrm{}^4},$$
(3)
where $`\epsilon `$ is the dielectric constant of the host semiconductor. At high hole densities $`p`$, the kinetic-energy term dominates and the mean-field critical temperature, $`T_c^{\mathrm{MF}}`$, is proportional to the Fermi wavevector, i.e., to $`p^{1/3}`$. Equations (2) and (3) also show that the band contribution to the susceptibility increases linearly with $`m^{}`$ while the exchange correction is proportional to $`(m^{})^2`$. (Using this simple band model as a guide, hole correlation effects, which will not be discussed in detail here, suppress the critical temperature by several percent for typical experimental hole densities, $`p0.1`$ nm-3.)

To obtain quantitative predictions for the critical temperature, it is necessary to evaluate the kinetic and exchange contribution to the itinerant hole susceptibility using a realistic six-band Kohn-Luttinger model, instead of the parabolic-band model. The band Hamiltonian contains the spin-orbit splitting parameter $`\mathrm{\Delta }_{\mathrm{so}}`$ and three other phenomenological parameters, $`\gamma _1`$, $`\gamma _2`$, and $`\gamma _3`$, whose values for the specific III-V host can be found, e.g., in Refs. . Only zinc-blende crystals are considered here; in nitrides, which crystallize in both zinc-blende and wurtzite structures, the critical temperature is expected to be only weakly dependent on the lattice configuration. The non-interacting-hole mean-field, $`T_c^{\mathrm{MF}}`$, and the exchange-enhanced, $`T_c^{\mathrm{ex}}`$, critical temperatures for GaAs and GaN host semiconductors doped with 5% of Mn are plotted in Fig. 1 as a function of the hole density. Transition temperatures for other III-V hosts and for hole densities $`p=0.1`$ and 0.5 nm-3 are listed in Table I. We have assumed $`J_{\mathrm{pd}}=55`$ meV nm-3 for all III-V hosts. In the density range considered, only the two heavy-hole and two light-hole bands are occupied in arsenides and antimonides. However, the mixing between these four bands and the two spin-orbit split-off bands is strong and must be accounted for. In nitrides and phosphides, spin-orbit coupling is weaker and all six bands are occupied by holes. The numerical results for the transition temperature are consistent with the qualitative analysis based on the parabolic band model: $`T_c^{\mathrm{MF}}`$ follows roughly the $`p^{1/3}`$ dependence, the $`1050`$% exchange enhancement of the critical temperature has a weaker density dependence.

For the thoroughly studied (Ga,Mn)As material doped with 5% of Mn and with $`p=0.35`$ nm-3, we find $`T_c^{\mathrm{MF}}=102`$ K and $`T_c^{\mathrm{ex}}=112`$ K. This result is in remarkably good agreement with the experimental transition temperature of 110 K . Good agreement between experiment and mean-field theory is also found in (In,Mn)As compounds, that have lower critical temperatures compared to their (Ga,Mn)As counterparts, primarily because they have smaller band masses. An especially important system is (Ga,Mn)N, which appears to exhibit room-temperature ferromagnetism. Our calculations for this III-V semiconductor give $`T_c^{\mathrm{MF}}400`$-600 K and $`T_c^{\mathrm{ex}}600`$-900 K for Mn content $`x=5`$% over the range of hole densities studied. Although we cannot compare those numbers directly to the recently measured $`T_c=228370`$ K in (Ga,Mn)N samples, since they have inhomogeneous doping profiles (resulting from the chemical vapor deposition growth technique), the mean-field theory predictions have the correct order of magnitude but appear to overestimate the transition temperature.

We now establish the quantitative reliability of the mean-field theory in arsenides (and also in antimonides and lower-carrier-density phosphides) and explain why mean-field theory overestimates the transition temperature in the heavy effective mass nitrides. Our arguments are based on the estimates of the suppression of the transition temperature due to collective spin-wave fluctuations. We emphasize below that mean-field theory would fail badly in any of these materials if the appropriate band model were really a single parabolic band with a mass equal to the material’s heavy-hole mass. Strong spin-orbit coupling in the valence band, which yields band states whose orbital and spin character changes completely as the band wave-vector varies over the Brillouin zone, enhances the spin stiffness, increases the energy of correlated reduced magnetization configurations, and reduces the importance of corrections to the mean-field-theory result.

Isotropic ferromagnets have spin-wave Goldstone collective modes whose energies vanish at long wavelengths,

$$\mathrm{\Omega }_k=Dk^2+𝒪(k^4),$$
(4)
where $`k`$ is the wavevector of the mode. Spin-orbit coupling breaks rotational symmetry and leads to a finite gap. According to our numerical studies, this gap is small however, much smaller than $`T_c`$ for example, and plays a role in magnetic fluctuations only at very low temperatures. Spin-wave excitations reduce the total spin by one, at an energy cost that is, at least at long wavelengths, much smaller than the mean-field value, $`J_{\mathrm{pd}}s^{(0)}`$ where $`s^{(0)}`$ is the mean-field band spin-density. The importance of these correlated spin excitations, neglected by mean-field theory, can be judged by evaluating an approximate $`T_c`$ bound based on the following argument which uses a Debye-like model for the magnetic excitation spectrum. When spin-wave interactions are neglected, the magnetization vanishes at the temperature where the number of excited spin waves equals the total spin of the ground state:

$$N_{\mathrm{Mn}}S=\frac{1}{2\pi ^2}_0^{k_D}𝑑kk^2n(\mathrm{\Omega }_k),$$
(5)
where $`n(\mathrm{\Omega }_k)`$ is the Bose occupation number and the Debye cutoff, $`k_D=(6\pi ^2N_{\mathrm{Mn}})^{1/3}`$. It follows that the ferromagnetic transition temperature cannot exceed

$$k_BT_c^{\mathrm{coll}}=\frac{2S+1}{6}k_D^2D(T_c^{\mathrm{coll}}).$$
(6)
In applying this formula to estimate $`T_c`$ we have approximated the temperature dependence of the spin stiffness by

$$D(T)=D_0S(T)/S,$$
(7)
as we did in Ref. , where $`D_0`$ is the zero-temperature stiffness, and $`S(T)`$ is the mean-field Mn polarization at a temperature $`T`$. If the difference between $`T_c^{\mathrm{coll}}`$ and $`T_c^{\mathrm{MF}}`$ is large, the typical local valence-band carrier polarization will remain finite above the critical temperature and ferromagnetism will disappear only because of the loss of long-range spatial order, the usual circumstance for transition metal ferromagnetism for example.

In discussing corrections to mean-field-theory $`T_c`$ estimates, we compare spin-stiffness results obtained with the simple two-band and realistic six-band models. Details on the formalism used to calculate $`D_0`$ can be found in Refs. . We find that the zero-temperature spin stiffness is always much larger in the six-band model. For (Ga,Mn)As, for example, the two-band model underestimates $`D_0`$ by a factor of $``$10-30 over the range of hole densities considered. Furthermore, the trend is different: in the two-band model the stiffness decreases with increasing density, while for the six-band description the initial is followed by a saturation. Even in the limit of low carrier concentrations, it is not only the (heavy-hole) mass of the lowest band which is important for the spin stiffness. In the realistic band model, heavy-holes have their spin and orbital angular momenta aligned approximately along the direction of the Bloch wavevector. Exchange interactions with Mn spins mix the heavy holes with more dispersive light holes. Our calculations show that heavy-light mixing is responsible for the relatively large spin stiffnesses which are responsible for the general success of mean-field theory. Crudely, the large mass heavy-hole band dominates the spin susceptibility and enables local magnetic order at high temperatures, while the dispersive light-hole band dominates the spin stiffness and enables long-range magnetic order. The multi-band character of the semiconductor valence band plays an essential role in the ferromagnetism of these materials.

In the insets of Figs. 1(a) and (b), solid lines show the spin stiffness for GaAs and GaN, doped with 5% of Mn and with $`p=0.4`$ nm-3, as a function of temperature. The critical temperatures, $`T_c^{\mathrm{coll}}`$, were obtained from these curves and from Eqs. (6) and (7). $`T_c^{\mathrm{coll}}`$, calculated using a non-interacting hole approximation, and $`T_c^{\mathrm{ex}}`$, calculated including exchange interactions but based on a mean-field treatment of the Mn spin-orientation thermodynamics, represent approximate lower and upper bound for the ferromagnetic transition temperature in (III,Mn)V semiconductors. With the exception of nitrides and higher-carrier-density phosphides, exchange enhancement and correlated fluctuations change $`T_c`$ by less than 20%, and these two influences will tend to cancel. This property explains the remarkable success of the mean-field theory in these itinerant ferromagnets. In the heavy effective mass nitrides, mean-field theory looses its quantitative accuracy, but still predicts the correct order of magnitude for the ferromagnetic transition temperature.

The dot-dashed lines in Fig. 1, and Table I summarize critical temperature estimates $`T_c^{\mathrm{est}}`$ obtained from Eqs. (6) and (7) and by shifting the $`D(T)`$-curves in the insets of Fig. 1 along $`x`$-axis by $`T_c^{\mathrm{ex}}T_c^{\mathrm{MF}}`$. This approach accounts, in an approximate way, for both exchange interactions and correlated spin fluctuations. For arsenides, antimonides, and lower-carrier-density phosphides, $`T_c^{\mathrm{est}}`$ is very close to $`T_c^{\mathrm{MF}}`$ justifying the mean-field description of ferromagnetism in these semiconductor compounds. The estimated critical temperatures for (Ga,Mn)N, $`T_c^{\mathrm{est}}=200400`$ K, reflect a more significant suppression of ferromagnetism due to collective excitations. Our theoretical $`T_c^{\mathrm{est}}`$ are in good agreement with available experimental transition temperatures in studied (III,Mn)V diluted magnetic semiconductors, including (Ga,Mn)N.

The work was supported by the DARPA/ONR Award No. N00014-00-1-095, the Welch Foundation, NSF under grant DMR-0115947, the Deutsche Forschungsgemeinschaft under the Emmy-Noether program KO 1987/2-1, the EU COST program, the Grant Agency of the Czech Republic under grant 202/02/0912, and by the Ministry of Education of the Czech Republic under grant OC P5.10.