text:
# Protein secondary structure prediction based on quintuplets

## 1 Introduction

Methods for predicting the secondary structure of a protein from its amino acid sequence have been developed for 3 decades. Besides neural network models and nearest-neighbor methods, the Chou-Fasman/GOR statistical method is well-established and commonly used. In 1974, assuming an oversimplified independency to cope with the large size 20 of the amino acid alphabets at a small size of database, Chou and Fasman (1974) derived a table of propensity for a particular residue to be in a given secondary structure state. By combining with a set of rules, the protein secondary structure was predicted using this propensity. Later, in the first version of the GOR program (Garnier, Osguthorpe, and Robson, 1978), the state of a single residue $`a_i`$ was predicted according to a window from $`i8`$ to $`i+8`$ surrounding the residue. Unlike Chou-Fasman which assumes that each amino acid individually influences its own secondary structure state, GOR takes into account the influence of the amino acids flanking the central residue on the central residue state by deriving an information score from the weight matrix describing 17 individual amino acid frequencies at sites $`i+k`$ with $`8k+8`$. By using a single weight matrix, the correlation among amino acids within the window was still ignored. In the later version GOR III (Gibrat, Garnier, and Robson, 1987), instead of single weight matrix for every structure state, 20 weight matrices, each of which corresponds to a specific type of the central residue, were used. These conditional weight matrices take the pair correlation between the central residue and a flanking one into account. In the most recent version of GOR (GOR IV, Garnier, Gibrat, and Robson, 1996), all pairwise combinations of amino acids in the flanking region were included.

Hidden Markov models (HMMs) (Rabiner, 1989) have been applied to molecular biology, in particular in gene-finding. There is a constant tendency in developing HMMs for protein structure prediction (Asai et al.,1993; Stultz et al., 1993; Karplus et al., 1997; Camproux et al. 1999). A probabilistic approach similar to the gene finder Genscan has been developped for protein secondary structure prediction without using sliding windows (Schmidler, Liu and Brutlag, 2000). In terms of Bayesian segmentation, the mothod integrated explicit models for secondary structure classes helices, sheets and loops with other observed structure aspects such as segment capping signals and length distributions, and reached an accuracy comparable to GOR.

Compared with DNA sequences, protein sequences are generally short, and their amino acid alphabet is of a large size 20. The range of lengths of secondary structure segments is rather small. The effect of duration might play a less important role. Here we develop a simple hidden Markov model with higher order correlations included for the secondary structure prediction. We propose several schemes for the amino acid alphabet reduction in order to incorporate residue correlation in the model. While the model is much simpler than the Bayesian segmentation model, its performance is still competitive.

## 2 Methods

A simplified version of the model can be constructed in the frame of the Chou-Fasman propensity scheme. We shall start with this model to explain several key points, and then discuss more realistic models.

In the Chou-Fasman approach, discriminant thresholds and post-prediction filtering rules are required. They can be avoided in a full probabilistic model.

As in the most methods, we consider 3 states $`\{h,e,c\}`$ generated from the 8 states of Kabsch and Sander (1983) by the coarse-graining $`H,G,Ih`$, $`Ee`$ and $`X,T,S,Bc`$. Let $`RR_{1:n}=R_1R_2\mathrm{}R_n`$ be a sequence of $`n`$ amino acid residues, and its corresponding secondary structure sequence be $`S=S_1S_2\mathrm{}S_n`$. The structure prediction is the mapping from $`R`$ to $`S`$. The main restriction to a structure sequence is that the shortest length of the consecutive state $`h`$ must be 3, and that of $`e`$ be 2. To cope with this restriction, we use triplet states instead of the 3 single states $`c`$, $`e`$ and $`h`$. In the total 27 triplets, only 19 are legitimate. The forbidden 8 are of the type $`\overline{e}e\overline{e}`$ or $`\overline{h}h\overline{h}`$, where $`\overline{e}`$ indicates a non-$`e`$, i.e. either $`c`$ or $`h`$, and the meaning of $`\overline{h}`$ is analogous. The first order Markov model for the triplets is the third order Markov model for the original mono-states. Any of the 19 triplets can only transit either to 3 or to 1 state. That is, the transition matrix is rather sparse. The 19 triplets and their transited states are listed in Table 1.

Denoting by $`\sigma _i`$ the triplet states, we may translate $`S_{1:n}`$ to $`\mathrm{\Sigma }_{2:n1}=\sigma _2\sigma _3\mathrm{}\sigma _{n1}`$ of length $`n2`$. (Note that the subscripts $`i`$ for $`\sigma _i`$ are from 2 to $`n1`$ for convenience.) The Markov process for $`\sigma `$ is characterized by the set of probabilities for initial states

$$\pi (l)=\mathrm{Prob}(\sigma _2=l),$$
(1)
and the transition rates

$$T_{kl}=\mathrm{Prob}(\sigma _i=l|\sigma _{i1}=k)=T(\sigma _{i1},\sigma _i).$$
(2)
Sequence $`R`$ is then related to $`\mathrm{\Sigma }`$ or $`S`$ by the emission probabilities

$$P(x|\delta )=\mathrm{Prob}(R_i=x|\sigma _i=\delta ),$$
(3)
which generate $`S_{2:n1}`$. Extra probabilities

$$Q(x|\delta )=\mathrm{Prob}(R_j=x|\sigma _{j\pm 1}=\delta ),j\{1,n\}$$
(4)
are needed for generating $`R_1`$ and $`R_n`$.

In this model the probability for the state sequence $`S`$ or $`\mathrm{\Sigma }`$ is

$$P(\mathrm{\Sigma })=\pi (\sigma _2)\underset{i=2}{\overset{n2}{}}T(\sigma _i,\sigma _{i+1}),$$
(5)
and the likelihood for $`R`$ to be at $`S`$ is

$$P(R|\mathrm{\Sigma })=Q(S_1|\sigma _2)Q(S_n|\sigma _{n1})\underset{i=2}{\overset{n1}{}}P(S_i|\sigma _i).$$
(6)
The joint probability is then

$$P(R,\mathrm{\Sigma })=P(R,S)=P(R|\mathrm{\Sigma })P(\mathrm{\Sigma }).$$
(7)
The predicted structure is infered as

$$\mathrm{\Sigma }^{}=\mathrm{arg}\mathrm{max}_\mathrm{\Sigma }P(\mathrm{\Sigma }|R)\mathrm{arg}\mathrm{max}_\mathrm{\Sigma }P(R,\mathrm{\Sigma }).$$
(8)
By means of the recursion relation

$`\mathrm{\Gamma }_2(\sigma )`$ $`=`$ $`Q(S_1|\sigma )\pi (\sigma )P(R_2|\sigma ),`$ (9)
$`\mathrm{\Gamma }_i(\sigma )`$ $`=`$ $`\mathrm{max}_\delta \mathrm{\Gamma }_{i1}(\delta )T(\delta ,\sigma )P(R_i|\sigma ),2<i<n,`$ (10)
$`\mathrm{\Gamma }`$ $``$ $`\mathrm{\Gamma }_n=\mathrm{max}_\delta \mathrm{\Gamma }_{n1}(\delta )Q(S_n|\delta ),`$ (11)
and by recording the pre-state leading to the maximal $`\mathrm{\Gamma }_i(\sigma )`$, the ‘best’ path $`\mathrm{\Sigma }^{}`$ can be traced back from the last $`\sigma _{n1}^{}`$. This is the so-called Viterbi algorithm for dynamic programming.

According to the standard forward-backward algorithm for HMMs, the forward and backward variables $`A_i`$ and $`B_i`$ may be defined as follows.

$`A_2(\sigma )`$ $`=`$ $`Q(S_1|\sigma )\pi (\sigma )P(R_2|\sigma ),`$ (12)
$`A_i(\sigma )`$ $`=`$ $`{\displaystyle \underset{\delta }{}}A_{i1}(\delta )T(\delta ,\sigma )P(R_i|\sigma ),2<i<n.`$ (13)
Similarly,

$`B_{n1}(\sigma )`$ $`=`$ $`Q(S_n|\sigma ),`$ (14)
$`B_i(\sigma )`$ $`=`$ $`{\displaystyle \underset{\delta }{}}B_{i+1}(\delta )T(\sigma ,\delta )P(R_{i+1}|\delta ),2i<n1.`$ (15)
It can be seen that for non-ending $`i`$, $`A_i(\sigma )=\mathrm{Prob}(R_{1:i},\sigma _i=\sigma )`$, $`B_i(\sigma )=\mathrm{Prob}(R_{i+1:n}|\sigma _i=\sigma )`$, and the partition function

$`Z`$ $``$ $`{\displaystyle \underset{\mathrm{\Sigma }}{}}P(R,\mathrm{\Sigma })={\displaystyle \underset{\sigma }{}}A_i(\sigma )B_i(\sigma )`$ (17)
$`=`$ $`{\displaystyle \underset{\sigma }{}}A_i(\sigma )B_i(\sigma ),2in1.`$ (18)
Denoting $`\sigma _i^{(0)}`$ the center $`S_i`$ of the triplet $`\sigma _i=S_{i1}S_iS_{i+1}`$, and introducing the characteristic function for $`z\{c,e,h\}`$

$$\delta (\sigma _i,z)=\{\begin{array}{cc}1,\hfill & \mathrm{if}\sigma _i^{(0)}=z,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}$$
(19)
we may infer single residue state from the marginal posterior

$$\mathrm{Prob}(S_i=z|R)\underset{\sigma }{}A_i(\sigma )B_i(\sigma )\delta (\sigma ,z).$$
(20)
This is the Baum-Welch algorithm for single residue states.

So far, only the correlation of conformation states has been considered. Residue triple will involve $`20^3=8000`$ parameters for each state $`\sigma `$. To avoid large training sets and model overfitting, a reduced amino acid alphabet is desired. For example, the reduction of 20 amino acids to 3 classes leads to only 27 combinations. However, there are as many as

$$\frac{1}{3!}\underset{j=0}{\overset{3}{}}C_3^j(1)^j(3j)^{20}5.8\times 10^8$$

ways of clustering 20 amino acids into 3 classes (Duran and Odell, 1974). For a given clustering $`\{a_i\}_{i=0}^{19}\{b_j\},b_j\{0,1,2\}`$, denoting by $`\rho _i=r_{i1}r_ir_{i+1}`$ the reduced residue triple corresponding to state $`\sigma _i`$, we may calculate the mutual information between the reduced residue triple $`\rho `$ and the triple state $`\sigma `$:

$$I(\rho ,\sigma )=H[\rho ]+H[\sigma ]H[\rho ,\sigma ],$$
(21)
where $`H[x]`$ is the Shannon entropy of $`x`$. (If the clustering independent $`H[\sigma ]`$ is ignored, $`I`$ would become the conditional entropy $`H[\sigma |\rho ]`$.) The best clustering may be determined by maximizing the objective function $`I(\rho ,\sigma )`$. The replacement of the above $`P(R_i|\sigma _i)`$ with $`P(\rho _i|\sigma _i)`$ leads to a version which takes some residue correlation into account.

A more realistic model uses quintuplets. Among the total $`3^5=243`$ quintuplets of the conformation states, only 75 are legitimate. Exclusion of 7 rare ones ($`eceeh`$, $`hceeh`$, $`heece`$, $`heech`$, $`heeeh`$, $`heehh`$, and $`hheeh`$) further reduces the total number of states into 68, which are listed in Table 2. We shall still use the same notation $`\sigma _i`$ for these 68 conformation states. To take residue quintuplet correlation into account, we substitute the central residue score $`P(R_i|\sigma _i)`$ with $`P(R_i|\sigma _i,r_{i2}r_{i1}r_{i+1}r_{i+2})`$, where $`r_i`$ stands for reduced residue classes. More words need to be said about the amino acid clustering. We have observed the fact that counts of $`ccccc`$, $`eeeee`$ and $`hhhhh`$ in databases are dominant over those of remaining 65. The various propensities of amino acid residues to different conformations imply that amino acid clustering should depend on conformations. We want to cluster amino acids separately for each conformation. For this purpose, for example, to find the best clustering at conformation $`c`$, we collect a subset of residue quintuplets whose conformation is $`ccccc`$. Denote by $`R_0`$ the central residue of a residue quintuplet, and by $`r_2r_1r_1r_2`$ the reduced classes of the other 4 residues. Taking the mutual information $`I(R_0,r_2r_1r_1r_2)`$ as the objective function, we determine the best clustering at $`c`$. To find at which position the residue depends most strongly on others, we calculate mutual information for nonreduced residue placed at different positions of a residue quintuplet. While the largest $`I`$ is found when the nonreduced residue is at the center for conformation $`e`$, the position of the nonreduced residue which gives the largest $`I`$ is the second position for conformation $`c`$, and is the fifth for $`h`$. However, for either $`c`$ or $`h`$, the largest $`I`$ is still very close to $`I(R_0,r_2r_1r_1r_2)`$. The mutual information excesses with respect to $`I(R_0,r_2r_1r_1r_2)`$ for different positions at conformations $`c`$, $`e`$ and $`h`$ are listed in Table 3. Thus, for simplicity, we always place the nonreduced residue at the center for all conformations. We then calculate refined residue scores from the conformation-dependent clustering. For example,

$$P(R_i|\sigma _i=ccchh,r_{i2}r_{i1}r_{i+1}r_{i+2})P(R_i|\sigma _i,r_{i2}^cr_{i1}^cr_{i+1}^hr_{i+2}^h),$$
(22)
where the superscript of $`r_i`$ indicates its conformation. The whole procedure for the dynamic programming remains almost the same, except for some care when dealing with two more end sites.

## 3 Result

We create a nonredundant set of 1612 non-membrane proteins for training parameters from PDB\_SELECT (Hobohm and Sander, 1994) with amino acid identity less than 25% issued on 25 September of 2001. The secondary structure for these sequences are taken from DSSP database (Kabsch and Sander, 1983). As mentioned above, the eight states of DSSP are coarse-grained into 3 states: $`h`$, $`e`$ and $`c`$. This learning set contains 268031 residues with known conformations, among which 94415 are $`h`$, 56510 are $`e`$, and 117106 are $`c`$. The size of the learning set is reasonable for training our parameters. There are 296 unknown residues. We add an extra ‘unknown’ amino acid category called $`X`$ to the 20 known ones.

In order to assess the accuracy of our approach, we use the following 2 test sets: Sets 1 and 2. A set of 124 nonhomologous proteins is created from the representative database of Rost and Sander (1993) by removing subunits A and B of hemagglutinin 3hmg, which are designated as membrane protein by SCOP (Murzin et al, 1995). The 124 sequences and the learning set are not independent of each other according to HSSP database (Dodge, Schneider and Sander, 1998). That is, some proteins of the 124 sequences are homologous with certain proteins in the learning set. Removing the homologous proteins from the 124 sequences and 5 seuqences with unknown amino acid segments longer than 6, we construct Set 1 of 76 proteins. Nonredundant 34 proteins with known structures of the CASP4 database issued in December of 2000 are taken as Set 2 (CASP4, 2000).

### 3.1 Amino acid clustering

The first method for clustering uses the mutual information $`I(\rho ,\sigma )`$ between the reduced oligo-peptide $`\rho `$ and its corresponding conformation $`\sigma `$. Setting the number of reduced classes at 3, 4 and 5, and fixing $`\sigma `$ to be triplets, we find the conformation-independent clustering of amino acids as shown in Table 4. Roughly speaking, class 0 is the hydrophobic, and is the same for the all 3 clusterings (except for the ‘unknown’ $`X`$ in the clustering into 5). Increasing the number of classes from 3 to 4, 5 results in new classes forming by single special amino acids P(Pro) and G(Gly). The similar clustering may be conducted in terms of quintuplets, but the results do not coincide with those from triplets.

The second method for clustering is conformation-dependent. We collect residue quintuplets of conformations $`ccccc`$, $`eeeee`$ and $`hhhhh`$ separately. The objective function for clustering now is the mutual information $`I(R_0,r_2r_1r_1r_2|\sigma )`$ with conformation $`\sigma `$ fixed. The results of clustering into 3 and 4 classes are listed in Table 5. Indeed, the cluster patterns for different conformations are quite dissimilar.

### 3.2 Secondary structure prediction

We shall index different models by $`N_s`$-$`N_a`$, where $`N_s`$ and $`N_a`$ are the numbers of conformation states and residue combinations, respectively. For example, the simplest model is model 19-21, which uses $`P(R_i|S_{i1}S_iS_{i+1})`$ to score residues. We may also use the propensity scores $`P(R_i|S_{i1}S_iS_{i+1})/P(R_i)`$, which result in an extra factor $`_{i=0}^nP(R_i)`$. Since the factor is independent of conformation sequences, it brings no new effect. When the scores are replaced by $`P(r_{i1}r_ir_{i+1}|S_{i1}S_iS_{i+1})`$ with $`r_i`$ being the reduced 3-class residues, we have model 19-$`3\times 3\times 3`$. The first and last sites of the predicted conformation of any sequences are always set at conformation $`c`$. For model 19-21 we determine the best conformation sequence as whole by the Viterbi algorithm and single residue conformation states by the Baum-Welch algorithm. To assess prediction methods, we calculate for each conformation the sensitivity $`s_n`$ and specificity $`s_p`$

$$s_n=\frac{TP}{TP+FN},s_p=\frac{TP}{TP+FP},$$
(23)
where $`TP`$, $`FP`$ and $`FN`$ are site counts of the ‘true positive’, ‘false positive’ and ‘true negative’ with respect to the observed real conformation. The results of model 19-21 are listed in Table 6, where the total sensitivity $`Q_3`$ for all conformations is also given. It is clearly seen that the inference from the Baum-Welch marginal posterior is significantly superior to that from the Viterbi algorithm in $`Q_3`$ value.

The next examined model uses reduced residue triplets. Reducing amino acids into 3 and 4 classes, we have models 19-27 and 19-64, respectively. The prediction accuracies of these two models are also listed in Tabel 6. We see that the inclusion of residue correlation dramatically improves the prediction accuracy.

The remaining part of Table 6 shows the prediction accuracies of quintuplet models. We examine the models on both test Sets 1 and 2. For all the models we take $`P(R_i|r_{i2}r_{i1}r_{i+1}r_{i+2},\sigma _i)`$ as the residue scores. We first compare conformation-independent with conformation-dependent clustering of amino acids. We find that the conformation-dependent clustering gains about 1 percent in the prediction accuracy for model 68-81$`\times `$21. The conformation-independent clustering is then not considered later on. Model 68-256$`\times `$21 contains more information about correlated residues, and has a better performance. The accuracies obtained on Set 1 are generally higher than those on Set 2. Besides the popular predictor GOR IV, there is another secondary structure predictor SSP (Solovyev and Salamov, 1991, 1994) based on discriminant analysis using single sequence. To compare with them, their accuracies on the same test sets are also listed in the table.

## 4 Discussions

We have presented simple hidden Markov models to predict secondary structure using single protein sequence. The hidden sequence is generated by a Markov process of multi-site conformation states. Considering that structure segments of proteins are generally short, we have ignored the duration effect, and focused on short range correlations. We proposed several schemes for coarse-graining the amino acid alphabet in order to include multi-residue correlation. Such reduction has been used in the Bayesian segmentation of protein secondary structure (Schmidler, et. al., 2000). However, here we derived the coarse-graining schemes specially for scoring residues to fit conformations. We have discussed only the principle of taking proper mutual information as an objective function for clustering, but did not exhaust all possibilities for clustering. For example, to diminish parameters, one may consider residue triplets at conformation quintuplet states. Another possibility is to use quartuples for both residues and conformations. One can also cluster quintuplet conformation states to less states.

There are rooms for further improvement of our approach. Simple weights (3 values) may be introduced to adjust residue scores according to its single-site conformation. Moreover, we may divide a training set into several, say 2, subsets according to residue statistics. Again, for this purpose the coarse-graining schemes help. The two subsets are then used separately for training to get refined models. We first classify a query sequence into one of the two categories, and then apply to it the corresponding refined model. We have tested this on the simple triplet models. The primitive result of up to 2% in accuracy improvement is encouraging. This, and the protein family recognition using multiple reduced amino acid residues, are under study.

> This work was supported in part by the Special Funds for Major National Basic Research Projects and the National Natural Science Foundation of China.

Table 1. Triplet states and their transited states.

| State | | Transited | State | | Transited |
| --- | --- | --- | --- | --- | --- |
| 0 | CCC | 0, 1, 2 | 10 | EEH | 11 |
| 1 | CCE | 3 | 11 | EHH | 18 |
| 2 | CCH | 4 | 12 | HCC | 0, 1, 2 |
| 3 | CEE | 8, 9,10 | 13 | HCE | 3 |
| 4 | CHH | 18 | 14 | HCH | 4 |
| 5 | ECC | 0, 1, 2 | 15 | HEE | 8, 9,10 |
| 6 | ECE | 3 | 16 | HHC | 12,13,14 |
| 7 | ECH | 4 | 17 | HHE | 15 |
| 8 | EEC | 5, 6, 7 | 18 | HHH | 16,17,18 |
| 9 | EEE | 8, 9,10 |  |  |  |

Table 2. 68 quintuplet states.

| 0 | CCCCC | 17 | CHHHE | 34 | EEECH | 51 | HEECC |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | CCCCE | 18 | CHHHH | 35 | EEEEC | 52 | HEEEC |
| 2 | CCCCH | 19 | ECCCC | 36 | EEEEE | 53 | HEEEE |
| 3 | CCCEE | 20 | ECCCE | 37 | EEEEH | 54 | HHCCC |
| 4 | CCCHH | 21 | ECCCH | 38 | EEEHH | 55 | HHCCE |
| 5 | CCEEC | 22 | ECCEE | 39 | EEHHH | 56 | HHCCH |
| 6 | CCEEE | 23 | ECCHH | 40 | EHHHC | 57 | HHCEE |
| 7 | CCEEH | 24 | ECEEC | 41 | EHHHE | 58 | HHCHH |
| 8 | CCHHH | 25 | ECEEE | 42 | EHHHH | 59 | HHEEC |
| 9 | CEECC | 26 | ECHHH | 43 | HCCCC | 60 | HHEEE |
| 10 | CEECE | 27 | EECCC | 44 | HCCCE | 61 | HHHCC |
| 11 | CEECH | 28 | EECCE | 45 | HCCCH | 62 | HHHCE |
| 12 | CEEEC | 29 | EECCH | 46 | HCCEE | 63 | HHHCH |
| 13 | CEEEE | 30 | EECEE | 47 | HCCHH | 64 | HHHEE |
| 14 | CEEEH | 31 | EECHH | 48 | HCEEC | 65 | HHHHC |
| 15 | CEEHH | 32 | EEECC | 49 | HCEEE | 66 | HHHHE |
| 16 | CHHHC | 33 | EEECE | 50 | HCHHH | 67 | HHHHH |

Table 3. Mutual information excesses ($`\times 10^3`$) for different positions of the nonreduced residue with respect to the mutual information for the nonreduced residue at the quintuplet center. The center is referred to as position 3.

|  | Position of the nonreduced | | | |
| --- | --- | --- | --- | --- |
| Conformation | 1 | 2 | 4 | 5 |
| c | -19.5 | 2.7 | -3.0 | -7.4 |
| e | -46.1 | -10.2 | -15.5 | -82.3 |
| h | -49.6 | -33.8 | -0.5 | 1.8 |

Table 4. Conformation-independent clustering of amino acids into 3, 4 and 5 classes.

| Amino Acid | A | V | C | D | E | F | G | H | I | W | K | L | M | N | Y | P | Q | R | S | T | X |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 3-class | 1 | 0 | 0 | 2 | 1 | 0 | 2 | 1 | 0 | 0 | 1 | 0 | 0 | 2 | 0 | 2 | 1 | 1 | 1 | 1 | 1 |
| 4-class | 1 | 0 | 0 | 2 | 1 | 0 | 2 | 2 | 0 | 0 | 1 | 0 | 0 | 2 | 0 | 3 | 1 | 1 | 2 | 2 | 1 |
| 5-class | 1 | 0 | 0 | 2 | 1 | 0 | 4 | 2 | 0 | 0 | 1 | 0 | 0 | 2 | 0 | 3 | 1 | 1 | 2 | 2 | 0 |

Table 5. Conformation-dependent clustering of amino acids into 3 and 4 classes.

| Amino Acid | A | V | C | D | E | F | G | H | I | W | K | L | M | N | Y | P | Q | R | S | T | X |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| c, 3-class | 1 | 2 | 0 | 0 | 1 | 2 | 0 | 0 | 2 | 2 | 1 | 2 | 2 | 0 | 2 | 0 | 1 | 1 | 1 | 1 | 2 |
| e, 3-class | 0 | 2 | 0 | 1 | 1 | 2 | 0 | 1 | 2 | 2 | 1 | 0 | 0 | 1 | 2 | 0 | 1 | 0 | 1 | 1 | 1 |
| h, 3-class | 0 | 2 | 0 | 1 | 1 | 2 | 1 | 1 | 2 | 2 | 1 | 2 | 2 | 1 | 2 | 1 | 1 | 1 | 1 | 1 | 0 |
| c, 4-class | 1 | 2 | 2 | 0 | 1 | 2 | 3 | 1 | 2 | 2 | 1 | 2 | 2 | 0 | 2 | 0 | 1 | 1 | 1 | 1 | 2 |
| e, 4-class | 0 | 1 | 2 | 3 | 3 | 1 | 2 | 0 | 1 | 0 | 3 | 2 | 2 | 3 | 1 | 2 | 3 | 0 | 0 | 3 | 2 |
| h, 4-class | 1 | 2 | 1 | 3 | 0 | 2 | 3 | 3 | 2 | 2 | 0 | 2 | 2 | 3 | 1 | 3 | 3 | 0 | 3 | 0 | 1 |

Table 6. Accuracy of secondary structure predictions for different models. Here, VI and BW stand for ‘Viterbi’ and ‘Baum-Welch’ algorithms, respectively. For all models 68-k, conformation-dependent reductions are used except for model 68-81$`\times `$21\* where the single conformation-independent reduction is used.

| Model | Test set | $`S_n^c`$ | $`S_p^c`$ | $`S_n^e`$ | $`S_p^e`$ | $`S_n^h`$ | $`S_p^h`$ | $`Q_3`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 19-21, VI | 2 | 48.69 | 62.92 | 7.94 | 64.09 | 84.17 | 48.07 | 53.24 |
| 19-21, BW | 2 | 70.67 | 62.27 | 30.74 | 54.55 | 66.12 | 60.11 | 60.45 |
| 19-27, VI | 2 | 62.50 | 63.87 | 56.02 | 46.58 | 55.82 | 61.48 | 58.63 |
| 19-64, VI | 2 | 64.10 | 65.31 | 50.11 | 47.59 | 62.45 | 63.06 | 60.50 |
| 68-81$`\times `$21\*, VI | 2 | 60.12 | 69.16 | 44.54 | 54.89 | 75.48 | 60.26 | 62.53 |
| 68-81$`\times `$21, VI | 2 | 62.01 | 68.21 | 44.26 | 57.79 | 76.51 | 61.98 | 63.64 |
| 68-81$`\times `$21, BW | 2 | 74.51 | 61.71 | 47.18 | 59.99 | 61.52 | 68.93 | 63.83 |
| 68-81$`\times `$21, VI | 1 | 70.92 | 66.51 | 54.27 | 53.45 | 62.00 | 69.07 | 64.46 |
| 68-81$`\times `$21, BW | 1 | 71.01 | 68.98 | 52.97 | 57.58 | 67.37 | 66.69 | 66.00 |
| 68-256$`\times `$21, VI | 2 | 65.35 | 68.41 | 54.73 | 58.41 | 70.11 | 64.58 | 64.86 |
| 68-256$`\times `$21, BW | 2 | 70.90 | 68.15 | 57.60 | 60.86 | 68.85 | 69.84 | 67.30 |
| 68-256$`\times `$21, VI | 1 | 68.34 | 69.97 | 55.34 | 57.68 | 70.12 | 66.23 | 66.18 |
| 68-256$`\times `$21, BW | 1 | 73.22 | 69.77 | 57.40 | 59.61 | 67.20 | 70.38 | 67.90 |
| GOR4 | 1 | 79.3 | 66.1 | 54.7 | 55.3 | 63.3 | 68.5 | 66.2 |
| SSP | 1 | 59.2 | 52.8 | 69.0 | 55.3 | 67.0 | 68.1 | 60.0 |
| GOR4 | 2 | 81.9 | 62.0 | 43.0 | 54.6 | 67.1 | 64.3 | 63.4 |
| SSP | 2 | 74.7 | 58.8 | 45.7 | 55.6 | 66.3 | 63.3 | 61.4 |


format: text

text:
# Theory of Two-Step Magnetic Ordering Phenomena in a Geometrically Frustrated Heisenberg Pyrochlore Antiferromagnet with Long Range Dipolar Interactions

## Abstract

A model for a Heisenberg antiferromagnet on a pyrochlore lattice with exchange and dipole-dipole interactions is studied via mean-field theory. In treating the dipoles by use of the Ewald method, a soft (critical) mode with a unique ordering wave vector is selected for all strengths of the dipole-dipole coupling. For weak dipoles a partially ordered, three sublattice spin structure (P state), with $`𝒒_{\mathrm{ord}}=\frac{1}{2}\frac{1}{2}\frac{1}{2}`$, is selected. A fully ordered, four sublattice spin structure (F state), with $`𝒒=000`$, competes with the P state and becomes the stable structure as the temperature is reduced. Our results are compared against other theoretical calculations and connection to recent experimental results for Gd2Ti2O7 are discussed.

Frustration is ubiquitous in condensed matter and arises whenever a system cannot minimize its total classical ground state energy by minimizing the energy of its pairwise interactions, pair by pair. In the context of magnetism, frustration commonly occurs in triangular or tetrahedral unit based lattice structures with antiferromagnetic (AFM) interactions, geometric frustration, and prevents the magnetic moments (spins) from pointing antiparallel to each other. Two and three dimensional lattices of corner sharing triangles or tetrahedra with nearest neighbor AFM Heisenberg exchange are highly frustrated and are particularly interesting. There, mean-field theory (MFT) finds a massive degeneracy of soft (critical) modes with no unique preferred state chosen at some critical temperature Reimers et al. (1991); Gingras and den Hertog (2001). As a consequence, frustrated systems exhibit a propensity for thermal and quantum fluctuations Moessner and Sondhi (2001); Hermele et al. . The possibility that exotic novel quantum states arise in frustrated quantum spin systems is currently attracting much attention.

It is generally expected that perturbations to the nearest neighbor spin Hamiltonian select a unique ordered state. In the insulating pyrochlore systems R2M2O7, where R3+ is a magnetic trivalent rare earth ion (R=Gd, Tb, Ho, Dy), and M=Ti4+ or Sn4+, non-magnetic (Fig. 1a), the leading perturbation to the nearest-neighbor Hamiltonian is the long-range magnetic dipole-dipole interaction. Due to its (i) long-range $`1/r^3`$ and (ii) anisotropic $`𝒓𝑺`$ spin$``$space coupling nature, one would naively expect dipolar interactions to select a unique and robust soft magnetization density (mode) at an ordering wave vector $`𝒒_{\mathrm{ord}}`$. However, explicit calculations in two specific cases have found dramatically different behavior.

In the Ho2(Ti,Sn)2O7 and Dy2Ti2O7 pyrochlores with local $`111`$ Ising spins, it has been found that the $`1/r^3`$ dipolar interactions are to a very large extend self-screened Gingras and den Hertog (2001). As a result, the system does not develop long-range order at a critical temperature $`T_cD_{\mathrm{nn}}`$, where $`D_{\mathrm{nn}}2`$ K is the dipolar energy scale at nearest-neighbors. Instead, the system develops an extensive quasi-degenerate “spin ice” state at $`TD_{\mathrm{nn}}`$ analogous to the proton disordered state in common ice water with the same residual extensive Pauling entropy Ramirez et al. (1999); Bramwell and Gingras (2001). Interestingly, detailed calculations (mean-field Gingras and den Hertog (2001) and Monte Carlo Melko et al. (2001)) find that the dipolar interactions do not produce an exact symmetry and thus allows for the weak selection of an ordered state at ordering wave vector $`𝒒_{ord}=001`$. Monte Carlo simulations employing nonlocal loop dynamics find a first order transition, removing all residual entropy, to the $`𝒒_{ord}=001`$ long-range ordered state at the critical temperature $`T_c/D_{\mathrm{nn}}0.071`$ Melko et al. (2001).

In the Gd2(Ti,Sn)2O7 pyrochlores, the Gd3+ ions have $`S=7/2`$ and $`L=0`$, hence negligible single-ion anisotropy. They are, therefore, good realizations of antiferromagnetically coupled Heisenberg spins with weak long-range magnetic dipole-dipole interactions ($`20\%`$ of nearest-neighbor AFM exchange) Raju et al. (1999). Here too, the effect of dipolar interactions are far from obvious. Two independent mean-field (Gaussian approximation) calculations have found that long-range dipole-dipole interactions truncated at large distance (i.e., 100 nearest-neighbors) give rise to a line of degenerate soft modes along the cubic $`111`$ direction Raju et al. (1999); Palmer and Chalker (2000). Palmer and Chalker argued that including quartic terms in the theory leads to the selection of a four sublattice Néel ordered state at $`𝒒_{ord}=000`$, see Fig. 1b. Experiments provide a different picture: (i) Neutron diffraction measurements at $`T=50`$mK find a partially ordered phase at $`𝒒=\frac{1}{2}\frac{1}{2}\frac{1}{2}`$, see Fig. 1c, with one disordered sublattice Champion et al. (2001). (ii) Mössbauer experiments Bertin et al. (2002) find strong evidence for a paramagnetic (PM) fluctuating spin down to $`35`$mK, too. (iii) Specific heat measurements Ramirez et al. (2002) find two transitions at $`T=0.7`$K and $`T=1.0`$K. All these experimental results are seemingly incompatible with the prediction of Ref. Palmer and Chalker, 2000. In view of these perplexing results from theory and experiments, we have revisited the mean-field calculations of Refs. Raju et al., 1999 and Palmer and Chalker, 2000. We have found that the degeneracy of soft modes along the $`𝒒=111`$ direction reported in Raju et al., 1999 and Palmer and Chalker, 2000 is only approximate. Instead, we find that long-range dipolar interactions, when properly re-summed to infinite distance, give rise to a weakly selected unique soft mode at $`𝒒_{ord}=\frac{1}{2}\frac{1}{2}\frac{1}{2}`$ with one disordered sublattice, similar to what is found in experiments (Fig. 1c) Champion et al. (2001). Just $`0.7`$% below that first transition, we find a second transition to the $`𝒒=000`$ state of Fig. 1b. We also find that weak second and third nearest neighbor exchange only acts to suppress or enhance the transition temperature of the $`𝒒=000`$ mode relative to that of the $`𝒒_{ord}=\frac{1}{2}\frac{1}{2}\frac{1}{2}`$ state.

We consider a model for large $`𝑺`$-spins on a pyrochlore lattice with exchange and dipolar interactions but negligible single-ion anisotropy. The appropriate Heisenberg Hamiltonian, $`H`$, is given by $`H=\frac{1}{2}_{i,j}_{a,b}_{u,v}𝒥_{uv}^{ab}(i,j)S_i^{a,u}S_j^{b,v}`$, where indices $`(i,j)`$ denote fcc Bravais lattice points, $`(a,b)`$ describe the tetrahedral basis, and $`(u,v)`$ define the components of classical Heisenberg spins, $`S_i^{a,u}`$. The $`𝒥_{uv}^{ab}(i,j)`$ interaction matrix contains exchange and dipolar contributions and is given by $`𝒥_{uv}^{ab}(i,j)=J_N\delta _{ij,N\mathrm{nn}}\delta ^{uv}D_{\mathrm{dd}}\left[\delta ^{uv}|𝑹_{ij}^{ab}|^33𝑹_{ij}^{ab,u}𝑹_{ij}^{ab,v}|𝑹_{ij}^{ab}|^5\right]`$, where $`J_N\delta _{ij,N\mathrm{nn}}`$ restricts the exchange to the $`N`$-th nearest neighbor distance, $`D_{\mathrm{dd}}DR_{\mathrm{nn}}^3`$ with $`D=\mu _\mathrm{o}\mu ^2/4\pi R_{\mathrm{nn}}^3`$, $`\mu `$ is the magnetic moment occupying the pyrochlore lattice site, $`R_{\mathrm{nn}}`$ is the nearest neighbor distance, and $`𝑹_{ij}^{ab}`$ is the vector separation between spins $`𝑺_i^a`$ and $`𝑺_j^b`$. The nearest neighbor exchange, $`J_1`$, sets the energy scale ($`J_1<0`$ is AFM). The effects of weak long-range exchange ($`J_2`$ and $`J_3`$) are also considered.

We first investigate the soft-mode spectrum of $`H`$ by studying the momentum dependence of the static susceptibility, $`\chi (𝒒)`$ within MFT Enjalran and Gingras ; Melko et al. . The mean-field free energy is $`=\mathrm{Tr}\{\rho H\}+T\mathrm{Tr}\{\rho \mathrm{ln}\rho \}`$ Reimers et al. (1991); Gingras and den Hertog (2001), where $`\mathrm{Tr}`$ is a trace over spin configurations and $`\rho `$ is the mean-field (one-particle) density matrix. The expression for $``$ is expanded to quadratic order in $`𝒎_i^a`$, the vector order parameter, Fourier transformed, and then diagonalized by a normal mode transformation. The diagonal form of $``$ is used to calculate the correlation function $`𝑺_i^a𝑺_j^b`$ in the definition of $`\chi (𝒒)`$,

$$\chi (𝒒)=(\beta /N_{\mathrm{cell}})\underset{(i,a),(j,b)}{}𝑺_i^a𝑺_j^b\mathrm{e}^{ı𝒒𝑹_{ij}^{ab}},$$
(1)
where $`\beta =1/T`$ in units of $`k_\mathrm{B}=1`$ and $`N_{\mathrm{cell}}`$ is the number of Bravais lattice points, and $`𝒒`$ is wave vector in the first Brillouin zone. The final form for the $`𝒒`$-dependent susceptibility is $`\chi (𝒒)/N_{\mathrm{cell}}=(\beta /3)_{\alpha ,\mu }|𝑭_\mu ^\alpha (𝒒)|^2(1\lambda _\mu ^\alpha (𝒒)/3T)^1`$, where the Greek indices label the normal modes ($`\alpha =1,2,3,4`$ and $`\mu =1,2,3`$) at each value of $`𝒒`$, and $`𝑭_\mu ^\alpha (𝒒)=_{a,u}U_{u,\mu }^{a,\alpha }(𝒒)=_a𝑼_\mu ^{a,\alpha }(𝒒)`$ is a 3-component vector. $`U(𝒒)`$ is the unitary matrix that diagonalizes $`𝒥(𝒒)`$ in the sublattice space with eigenvalues $`\lambda _\mu ^\alpha (𝒒)`$. The static susceptibility is a maximum at the ordering wave vector, $`𝒒_{\mathrm{ord}}`$, which also defines the maximum eigenvalue, i.e., $`\lambda ^{\mathrm{max}}(𝒒_{\mathrm{ord}})=\mathrm{max}_𝒒\{\lambda ^{\mathrm{max}}(𝒒)\}`$, where $`\lambda ^{\mathrm{max}}(𝒒)`$ is the maximum eigenvalue at $`𝒒`$ and $`\mathrm{max}_𝒒`$ selects the global maximum for all $`𝒒`$. The mean-field critical temperature is given by, $`T_c^{\mathrm{MF}}=\lambda ^{\mathrm{max}}(𝒒_{\mathrm{ord}})/3`$.

To find $`𝒒_{\mathrm{ord}}`$ for our model, we search for the maximum in the soft mode spectrum ($`\lambda ^{\mathrm{max}}(𝒒)`$). We consider AFM first neighbor exchange, $`J_1=1.0`$, and a variable dipole-dipole interaction strength, $`D`$. The dipolar interactions in $`𝒥(𝒒)`$ are summed directly to a cut off separation distance, $`R_c`$, measured in units of the nearest neighbor distance and to the infinite range limit through the use of the Ewald method Enjalran and Gingras . The regime of weak dipole-dipole interactions is of greatest interest to us because the value $`D/|J_1|=0.20`$ faithfully approximates Gd2Ti2O7Raju et al. (1999). At this point, we find a spectrum in the $`(hhl)`$ plane that appears degenerate along the $`(hhh)`$ diagonal, in agreement with earlier works Raju et al. (1999); Palmer and Chalker (2000) (Fig. 6 in Ref. Raju et al., 1999). However, when we examine the spectrum along the $`(hhh)`$ diagonal in greater detail, we observe the selection of a global maximum as one increases the cut off distance of the dipolar lattice sum. Figure 2 displays results for several values of $`R_c`$ and the Ewald limit. An ordering wave vector occurs at $`𝒒_{\mathrm{ord}}=\frac{1}{2}\frac{1}{2}\frac{1}{2}`$ (in units of $`2\pi /a`$), and emerges from a smoothly varying spectrum along $`(hhh)`$ in the Ewald limit. We note that for dipole-dipole interactions summed to $`R_c<500`$, $`\lambda ^{\mathrm{max}}(𝒒)`$ displays substantial variations making it difficult to ascertain the existence of a critical ordering wave vector. These fluctuations are damped with increasing $`R_c`$ and approach the Ewald limit for $`R_c>1000`$. A similar sensitivity of the soft-mode spectrum to the dipole-dipole cut off has been observed in the spin ice materialsEnjalran and Gingras . In a broader expanse of phase space, we find three regions for an ordering wave vector: (i) $`𝒒_{\mathrm{ord}}=\frac{1}{2}\frac{1}{2}\frac{1}{2}`$ for $`D/|J_1|5.60`$, (ii) $`𝒒_{\mathrm{ord}}=001`$ for $`5.60<D/|J_1|<5.85`$, (iii) $`𝒒_{\mathrm{ord}}=000`$ for $`D/|J_1|5.85`$. Earlier work that imposed a cut off on the dipole-dipole lattice sum found a degenerate spectrum along the $`(hhh)`$ diagonal for $`D/|J_1|5.70`$, with the selection of an ordering wave vector near $`(000)`$ for larger values of $`D`$, (see Fig. 2 in Ref. Palmer and Chalker, 2000).

The spin structure for $`𝒒_{\mathrm{ord}}=\frac{1}{2}\frac{1}{2}\frac{1}{2}`$ can be obtained from the corresponding eigenvectors, $`U_{u,\mu }^{a,\alpha }(\frac{1}{2}\frac{1}{2}\frac{1}{2})`$. About a single tetrahedron one infers that three sublattices have equal moments with all spins confined to a plane and a fourth sublattice that is PM, i.e., zero moment, see Fig. 1c. When this configuration is extended over the pyrochlore lattice and modulated by $`𝒒_{\mathrm{ord}}`$, one has ordered kagomé planes along the $`[111]`$ direction with a reversal of spin direction in alternating planes, all interstitial sites are PM. This partially ordered three sublattice structure (P state) has been observed in neutron scattering experimentsChampion et al. (2001) and is supported by MössbauerBertin et al. (2002) experiments on Gd2Ti2O7. The mean-field transition temperature to this state is $`T_c^\mathrm{P}=\lambda ^{\mathrm{max}}(\frac{1}{2}\frac{1}{2}\frac{1}{2})/3=0.94008|J_1|`$ ($`D/|J_1|=0.2`$).

The mode $`𝒒=000`$ corresponds to a fully ordered (four sublattice) equal moment structure (F state) initially predicted to be the ground state of a $`D/|J_1|=0.20`$ model Palmer and Chalker (2000). We find the limit of stability for the paramagnetic state against an ordering to the F state occurs at $`T_c^\mathrm{F}=\lambda ^{\mathrm{max}}(000)/3=0.93923|J_1|`$, which is very close to $`T_c^\mathrm{P}=0.94008|J_1|`$ num . Given the close proximity of $`T_c^\mathrm{P}`$ and $`T_c^\mathrm{F}`$ at the Gaussian level, we study the evolution of the respective free energies as a function of temperature below both values of $`T_c`$.

The real-space mean-field theory (RSMFT) begins with the mean-field decomposition $`(S_i^{a,u}m_i^{a,u})(S_j^{b,v}m_j^{b,v})=0`$ applied to $`H`$. The free energy is determined from $`f=(1/\beta )\mathrm{ln}Z_{\mathrm{MF}}`$, where $`Z_{\mathrm{MF}}`$ is the RSMFT partition function. For our model, one obtains

$$f=\frac{1}{2}\underset{i,j}{}\underset{a,b}{}\underset{u,v}{}𝒥_{uv}^{ab}(i,j)m_i^{a,u}m_j^{b,v}\frac{1}{\beta }\underset{i,a}{}\mathrm{ln}Z_i^a,$$
(2)
where $`Z_i^a=(4\pi /\beta |𝑩_i^a|)\mathrm{sinh}(\beta |𝑩_i^a|)`$ is the single-site partition function and $`𝑩_i^a`$ is the local field. The sublattice magnetization $`𝒎_i^a`$ (order parameter) is obtained from differentiation of $`f`$ with respect to local $`𝐁_i^a`$ field. One finds $`𝒎_i^a=\frac{𝑩_i^a}{|𝑩_i^a|}\left[\mathrm{coth}(\beta |𝑩_i^a|)\frac{1}{\beta |𝑩_i^a|}\right]`$, where the field at site $`(i,a)`$ is due to all the other moments in the system acting at site $`(i,a)`$, i.e., $`𝑩_i^a=_uB_i^{a,u}=_u_{j,b,v}𝒥_{uv}^{ab}(i,j)m_j^{b,v}`$. It is convenient to factor out the magnitude of the order parameter and write, $`𝑩_i^a=m𝒃_i^a`$, where $`𝒎_i^a=m\widehat{n}_i^a`$ and $`\widehat{n}_i^a`$ is a unit vector that defines the orientation of the order parameter at site $`(i,a)`$. Using this form for $`𝑩_i^a`$, one expands the right-hand-side of the above expression for $`𝒎_i^a`$ for small $`m`$ to obtain $`T_c`$ for RSMFT, $`T_c^{\mathrm{MF}}=|𝒃_i^a|/3`$.

We work on a lattice of $`L\times L\times L`$ cubic cells, or $`16L^3`$ spins, with periodic boundary conditions and $`D/|J_1|=0.20`$. The lattice is tiled with one of the above order parameters. We consider systems of size $`L=4,6`$ and calculate the dipole-dipole interactions via the Ewald method in real-space Melko et al. . Our procedure is to solve self-consistently for the sublattice magnetization, $`𝒎_i^a`$, at temperature $`T/|J_1|`$ and then calculate the corresponding free energy from Eq. 2, i.e., $`f/|J_1|`$. The real-space transition temperatures are $`T_c^\mathrm{P}=|𝒃_i^{a(\mathrm{P})}|/3=0.94002|J_1|`$ and $`T_c^\mathrm{F}=|𝒃_i^{a(\mathrm{F})}|/3=0.93917|J_1|`$ for $`L=6`$, in good agreement with the above $`𝒒`$-space results. For $`L=8`$, the real-space and $`𝒒`$-space results for $`T_c^\mathrm{P}`$ and $`T_c^\mathrm{F}`$ agree within $`10^7`$. Free energy versus temperature curves for the F and P ordered states on a $`L=4`$ lattice are displayed in Fig. 3. At high temperatures, near either $`T_c^\mathrm{P}`$ or $`T_c^\mathrm{F}`$, the free energy curves appear degenerate. However, the insets to Fig. 3 demonstrate that the respective curves actually cross. We performed linear fits to points from the two curves about the crossing temperature. Solving these equations we find the temperature $`T_{\mathrm{cross}}0.93358|J_1|`$. Above $`T_{\mathrm{cross}}`$, the P state has a lower free energy and is thus the preferred state. Below $`T_{\mathrm{cross}}`$, the free energy of the F state drops below that of the P state and eventually becomes the stable ground state of the system. We find $`T=0`$ energies per spin of $`E_o^\mathrm{F}=1.40856|J_1|`$ and $`E_o^\mathrm{P}=1.05737|J_1|`$. Hence, in this very simple mean-field picture, the three sublattice P state is the initially preferred structure at high temperatures, but the nearby four sublattice F state competes with it and very rapidly wins out as the temperature is lowered.

The proximity of the P and F states in MFT suggests it would be difficult to realize their experimental differentiation (or in Monte Carlo simulations). However, the inclusion of weak Greedan et al. (1998) long-range exchange in our model separates the transition temperatures of these two states and effectively reduces their competition. We note that a $`J_2`$ or $`J_3`$ has no effect on the transition temperature of the P state gdt . In the Fourier transform of the $`J_2`$ and $`J_3`$ matrices, the particular combination of the lattice symmetry and the $`q=\pi \pi \pi `$ phase yields zero matrix elements for $`J_2(\pi \pi \pi )`$ and $`J_3(\pi \pi \pi )`$ in the above calculation of $`\chi (𝒒)`$ at $`𝒒=\pi \pi \pi `$. Hence, we find that $`J_3=0.02J_1`$ (AFM) suppresses the F state to $`T_c^\mathrm{F}=0.859237|J_1|`$ ($`L=4`$). A small FM $`J_2`$ has a similar effect but of reduced magnitude. Conversely, a small AFM $`J_2`$ or FM $`J_3`$ enhances the F state by raising its transition temperature above $`T_c^\mathrm{P}`$.

Hence, experimentally reasonable values for long-range exchange Greedan et al. (1998) could allow for the detection of two well-separated transitions in Gd2Ti2O7. Although two transitions have been observed in specific heat measurements at $`T1`$K and $`T0.7`$K Ramirez et al. (2002), neutron scattering results only find the three sublattice P state down to $`T50`$mK Bramwell . We note that a $`1`$-$`𝒌`$ structure is purposed experimentally for the P state, but a $`4`$-$`𝒌`$ spin state is a possibility. In our RSMFT, a $`4`$-$`𝒌`$ spin structure corresponds to a unit cubic cell in which one of the tetrahedra is fully disordered bra . It is unclear at this time whether it is a material (sample) specific effect or something else that causes the P state to freeze out in the neutron experiments as the material is cooled from $`T1`$ K down to $`50`$mK Champion et al. (2001). Recent ESR results have also found peculiar behavior in Gd2Ti2O7, where a large $`[111]`$ anisotropy develops spontaneously below $`80`$K Hassan et al. . Clearly, more experimental and theoretical work is needed to elucidate the peculiar thermodynamic properties of Gd2Ti2O7.

In conclusion, our calculations above have identified another system, apart from spin ice, where the long-range $`1/r^3`$ nature of dipolar interactions cause a large degree of symmetry restoration on the $`𝒒`$dependent susceptibility. However, this symmetry restoration is incomplete, and gives rise ultimately to a well defined ordering state upon cooling from the fully disordered paramagnetic state.

###### Acknowledgements.

We would like to acknowledge Steve Bramwell, Ying-Jer Kao, Adrian del Maestro, Jean-Yves Delannoy, and Hamid Molavian for many useful discussions. This work is supported by NSERC of Canada, Research Corporation and the Province of Ontario.


text:
# Proton strange form factors and SAMPLE experiments

## I Introduction

Recently, the SAMPLE Collaboration reported the experimental data of the proton strange magnetic form factor through parity violating electron scattering at a small momentum transfer $`Q_S^2=0.1(\mathrm{GeV}/\mathrm{c})^2`$ sample01

$$G_M^s(Q_S^2)=+0.14\pm 0.29(\mathrm{stat})\pm 0.31(\mathrm{sys})\mathrm{n}.\mathrm{m}..$$
(1.1)
On the other hand, baryons were described by topological solitons skyrme61 ; witten83npb ; brown83prl ; hong98 ; hong02pr  and the MIT bag model chodos74  was later unified with the Skyrmion model to yield the chiral bag model (CBM) gerry791 , which then includes the pion cloud degrees of freedom and the chiral invariance consistently. Moreover, the soliton was exploited to yield superqualiton hong01  in color flavor locking phase rajagopal00 .

The QCD is the basic underlying theory of strong interaction, from which low energy hadron physics should be attainable. Moreover, for hadron structure calculations, the coupling constant $`g`$ is not a relevant expansion parameter of QCD. Long ago, ’t Hooft noted that $`1/N_c`$ could be regarded as expansion parameter of QCD thooft74  where $`N_c`$ is the number of colors and $`gN_c^2`$ is kept constant. The properties of large $`N_c`$ limit of the QCD can be satisfied by the meson sector of the nonlinear sigma model such as the Skyrmion model.

Now, in order to study the hadron physics phenomenology, we start with the SU(3) Skyrmion lagrangian of the form witten83npb

$$L=\mathrm{d}^3x\left[\frac{f_\pi ^2}{4}\mathrm{tr}(l_\mu l^\mu )+\frac{1}{32e^2}\mathrm{tr}[l_\mu ,l_\nu ]^2\right]+L_{WZW}$$
(1.2)
where $`l_\mu =U^{}_\mu U`$ and $`U`$ SU(3) is described by pseudoscalar meson fields $`\pi _a`$ $`(a=1,2,\mathrm{},8)`$ and the topological aspects can be included via the WZW action witten83npb . Assuming maximal symmetry, we introduce the hedgehog ansatz $`U_0`$ embedded in the SU(2) isospin subgroup of SU(3) to yield the topological charge

$$Q=\frac{1}{2\pi }\chi _E(\theta \mathrm{sin}\theta \mathrm{cos}\theta )=1$$
(1.3)
where $`\theta `$ is the chiral angle and $`\chi _E`$ is the Euler characteristic being an inter two in the spherical bag surface.

## II BRST symmetries in Skyrmion

In order to define the spin and isospin one can quantize, in the SU(2) Skyrmion for instance, the zero modes via $`U_0AU_0A^{}`$ and $`A(t)=a^0+i\stackrel{}{a}\stackrel{}{\tau }`$ with $`a^\mu `$ being the collective coordinates. One can then obtain the Lagrangian $`L=m_0+2i_1\dot{a}^\mu \dot{a}^\mu `$ where the static mass $`m_0`$ and the moment of inertia $`i_1`$ are calculable in the Skyrmion model. Introducing the canonical momenta $`\pi ^\mu `$ one can obtain the canonical Hamiltonian $`H=m_0+\frac{1}{8i_1}\pi ^\mu \pi ^\mu `$. Note that the second-class geometrical constraints

$$\mathrm{\Omega }_1=a^\mu a^\mu 10,\mathrm{\Omega }_2=a^\mu \pi ^\mu 0$$
(2.4)
should be treated via the Dirac brackets dirac64 . However, in the Dirac quantization, one has difficulties in finding the canonically conjugate pair, which were later overcome bft  by introducing pair of auxiliary Stuckelberg fields $`\theta `$ and $`\pi _\theta `$ with $`\{\theta ,\pi _\theta \}=1`$. In the Skyrmion the first-class constraints $`\stackrel{~}{\mathrm{\Omega }}_1=a^\mu a^\mu 1+2\theta `$ and $`\stackrel{~}{\mathrm{\Omega }}_2=a^\mu \pi ^\mu a^\mu a^\mu \pi _\theta `$ were constructed hong99prd  to satisfy the strongly involutive Lie algebra $`\{\stackrel{~}{\mathrm{\Omega }}_1,\stackrel{~}{\mathrm{\Omega }}_2\}=0`$. Similarly, the first-class Hamiltonian was formulated to yield the baryon mass spectrum

$$m_B=m_0+\frac{1}{2i_1}\left[J(J+1)+\frac{1}{4}\right]$$
(2.5)
with the isospin quantum number $`J`$. Here note that an additional global shift is due to the Weyl ordering correction. Following the BRST quantization scheme brst  with (anti)ghost and their Lagrangian multiplier fields, we obtain the BRST symmetric Lagrangian hong99prd ,

$$L_{eff}=m_0+\frac{2i_1\dot{a}^\mu \dot{a}^\mu \dot{\theta }\dot{b}}{12\theta }\frac{2i_1\dot{\theta }^2}{(12\theta )^2}2i_1(12\theta )^2(b+2\overline{c}c)^2+\dot{\overline{c}}\dot{c}$$
(2.6)
invariant under the transformations, $`\delta _\lambda a^\mu =\lambda a^\mu c`$, $`\delta _\lambda \theta =\lambda (12\theta )c`$, $`\delta _\lambda \overline{c}=\lambda b`$ and $`\delta _\lambda c=\delta _\lambda b=0`$. (For more details of the BRST quantization of the SU(2) and SU(3) Skyrmions, see hong99prd  and hong01prd , respectively.)

## III Baryon strange form factors

Next, we consider the CBM which is a hybrid of two different models: the MIT bag model at infinite bag radius on one hand and Skyrmion model at vanishing radius on the other hand. (The explicit CBM lagrangian is given in Ref. hong02pr  for instance.) In the CBM the total topological charge $`Q`$ in Eq. (1.3) is now splitted into the meson and quark pieces to satisfy the Cheshire cat principle ccp . Moreover, the quark fractional charge is given by sum of integer one (from valence quarks) and the quark vacuum contribution, which is also rewritten in terms of the eta invariant atiyah .

In the collective quantization of the CBM, we explicitly obtain the proton magnetic moment hong93 ; hong93npa

$$\mu _p=\frac{1}{90}(9I_1+24I_2+12I_3+16I_44I_5)+\frac{2I_6}{1125}\left(9I_1+4I_28I_3\right)$$
(3.7)
with the inertia parameters $`I_n`$ $`(n=1,\mathrm{},6)`$ calculable in the CBM. Similarly we construct the baryon octet magnetic moments to reproduce the Coleman-Glashow sum rules coleman ; hong93npa  such as $`U`$-spin symmetries, $`\mu _{\mathrm{\Sigma }^+}=\mu _p`$, $`\mu _{\mathrm{\Xi }^0}=\mu _n`$ and $`\mu _\mathrm{\Xi }^{}=\mu _\mathrm{\Sigma }^{}`$.

Now we define the Dirac and Pauli EM form factors via

$$p+q|\widehat{V}^\mu |p=\overline{u}(p+q)\left[F_{1B}(q^2)\gamma ^\mu +\frac{i}{2m_B}F_{2B}(q^2)\sigma ^{\mu \nu }q_\nu \right]u(p)$$
(3.8)
where $`q`$ is momentum transfer and $`\sigma ^{\mu \nu }=\frac{i}{2}(\gamma ^\mu \gamma ^\nu \gamma ^\nu \gamma ^\mu )`$ and $`m_B`$ is baryon mass. The Sachs form factors are then given by $`G_M=F_{1B}+F_{2B}`$ and $`G_E=F_{1B}+\frac{q^2}{4m_B^2}F_{2B}`$ so that, at zero momentum transfer, the Pauli strange form factor is identical to the Sachs strange form factor: $`F_{2B}^s(0)=G_M^s(0)`$. In the SAMPLE experiment, they measured the neutral weak form factor

$$G_M^{Z,p}=\left(\frac{1}{4}\mathrm{sin}^2\theta _W\right)G_M^p\frac{1}{4}G_M^n\frac{1}{4}G_M^s$$
(3.9)
with $`G_M^p`$ and $`G_M^n`$ being the proton and neutron Sachs form factors, to predict the proton strange form factor (1.1) which is positive value contrary to the negative values from most of the model calculations except the predictions hong93 ; hong97plb  of the SU(3) CBM and the recent predictions of the chiral quark soliton model kim98  and the chiral perturbation theory meissner00 ; vankolck . (See hong02pr  for more details.) In the CBM the proton strange form factor is given by hong93

$$F_{2N}^s(0)=\frac{1}{60}(21I_14I_22I_34I_42I_5)+\frac{I_6}{2250}(129I_1+76I_252I_3)$$
(3.10)
which, after some algebra with the other baryon octet strange form factors, yields the sum rule for the proton strange form factor in terms of the baryon octet magnetic moments only (for the other baryon sum rules see hong01hep )

$$F_{2N}^s(0)=\mu _p\mu _\mathrm{\Xi }^{}(\mu _p+\mu _n)\frac{1}{3}(\mu _{\mathrm{\Sigma }^+}\mu _{\mathrm{\Xi }^0})+\frac{4}{3}(\mu _n\mu _\mathrm{\Sigma }^{}).$$
(3.11)
Explicitly calculating the inertia parameters $`I_n`$ numerically in (3.10), we predict the proton strange form factor, 0.30 n.m. and exploiting the experimental data for the baryon octet magnetic moments in (3.11) we obtain 0.32 n.m.. These predictions are comparable to the SAMPLE experimental data (1.1) and are shown in Table 1, together with those of the other baryon strange form factors.

## IV Conclusions

In conclusion, we discussed the SAMPLE experiments in the topological solitons such as the Skymion and chiral models to predict baryon strange form factors. We also exploited the Dirac quantization associated with ghost field degrees of freedom to construct the BRST invariant chiral Lagrangian. In future, it will be interesting to consider in these chiral models anapole form factors discussed in the literatures vankolck .

This work is financially supported in part by the Korea Science and Engineering Foundation Grant (R01-2000-00015).


text:
# Untitled Document

Abstract

We have performed a search for the leptonic $`B`$ decay $`B^+\mu ^+\nu _\mu `$ with data collected at the $`\mathrm{{\rm Y}}(4S)`$ resonance by the B A B AR experiment at the PEP-II storage ring. With an integrated luminosity of approximately 81.4 $`\text{ fb}^1`$ (88.4 million $`B\overline{B}`$ pairs), we find no convincing evidence for a signal and set a preliminary upper limit on the branching fraction of $`(B^+\mu ^+\nu _\mu )<6.6\times 10^6`$ at the 90% confidence level.


text:
# Energy Dependence of the Cronin Effect from Non-Linear QCD Evolution1footnote 11footnote 1We dedicate this work to the memory of Ian Kogan.

## Abstract

The non-linear evolution of dense partonic systems has been suggested as one of the novel physics mechanisms relevant to the dynamics of hadron–nucleus and nucleus–nucleus collisions at collider energies. Here we study to what extent the description of Cronin enhancement in the framework of this non-linear evolution is consistent with the recent observation in $`\sqrt{s}=200`$ GeV d–Au collisions at the Relativistic Heavy Ion Collider. We solve the Balitsky-Kovchegov (BK) evolution equation numerically for several initial conditions encoding Cronin enhancement. We find that the properly normalized nuclear gluon distribution is suppressed at all momenta relative to that of a single nucleon. Calculating the resulting spectrum of produced gluons in p–A and A–A collisions, we establish that the nonlinear QCD evolution is unable to generate a Cronin type enhancement, and that it quickly erases any such enhancement which may be present at lower energies.

The observation that the ratio of particle yields in p–A and A–A, scaled by the number of collisions, exceeds unity in an intermediate transverse momentum range of a few GeV, is commonly referred to as the Cronin effect. This was first seen at lower fixed target energies Cronin:zm  and was recently confirmed in $`\sqrt{s}=200`$ GeV d–Au collisions at RHIC dA . The current interest focuses mainly on comparing this Cronin enhancement in d–Au to the relative suppression of produced hadrons in Au–Au collisions at the same center of mass energy and in the same transverse momentum range gold . The opposite trend of the two effects and their centrality dependence suggests that d–Au data may serve as an efficient benchmark measurement to distinguish between the two different physical mechanisms suggested for the relative suppression of hadron spectra in Au–Au collisions: initial state parton saturationklm  and final state jet quenching Wang:2003aw .

The physics of dense partonic systems and their non-linear perturbative evolution to higher energy has motivated several attempts at understanding bulk properties of ultra-relativistic heavy ion collisions such as the multiplicity, rapidity distribution and centrality dependence of particle production kln ; klm . In particular it has been suggested that saturation effects can account for the suppression of the high-$`p_T`$ hadronic spectra in Au–Au collisions at RHIC. On the other hand, it is known that saturation models based on multiple scattering (the so called Glauber-Mueller glaubermuller  or McLerran-Venugopalan mv  models) exhibit Cronin enhancement in p–A Gelis:2002nn ; bkw ; kkt  and A–A bkw ; Jalilian-Marian:2003mf . In these models, the saturation of low $`p_T`$ gluons is the result of a redistribution of gluons in transverse phase space kovchegov ; jkmw  which does not change the total number of gluons, thus resulting in a compensating enhancement at momenta just above the saturation momentum $`Q_s`$. What is not fully understood is i) whether such Cronin enhancement encoded in the initial condition of a nuclear wave function persists in the non-linear perturbative QCD evolution to higher energy and ii) whether such Cronin enhancement can be generated by the non-linear evolution itself. This paper goes beyond earlier discussions bkw ; kkt ; Jalilian-Marian:2003mf  by providing the first complete (numerical) answer to these questions. We do not address other approaches to Cronin enhancement othercronins .

We start from the Balitsky-Kovchegov (BK) evolution equation balitsky ; k , which describes the evolution of the forward scattering amplitude $`N(𝐫,y)`$ of a QCD dipole of transverse size $`|𝐫|`$ with rapidity $`Y`$ and $`y=(\alpha _sN_c/\pi )Y`$,

$`{\displaystyle \frac{dN(|𝐫|,y)}{dy}}={\displaystyle \frac{1}{2\pi }}{\displaystyle d^2𝐳\frac{(𝐫𝐳)𝐳}{(𝐫𝐳)^2𝐳^2}}`$ (1)
$`\times [N(|𝐫𝐳|)+N(|𝐳|)N(|𝐫|)N(|𝐫𝐳|)N(|𝐳|)].`$
The unintegrated gluon distribution is related to the inclusive gluon distribution $`\varphi (k)\frac{d(xG(x,k^2))}{d^2kd^2b}`$ and is given in terms of the dipole amplitude

$$\varphi (k)=\frac{d^2r}{2\pi r^2}\mathrm{exp}\{i𝐫𝐤\}N(r).$$
(2)
In the following, we also use the modified gluon distribution

$$h(k)=k^2_k^2\varphi (k).$$
(3)
The two definitions coincide for the leading order perturbative distribution $`\varphi (k)\frac{1}{k^2}`$, but are different in general, and especially at low momenta.

Using the second order Runge-Kutta algorithm Braun:2000wr , we solve the BK equation (1) numerically with 8000 equally spaced intervals in $`\mathrm{ln}k`$-space between $`15`$ and $`35`$ and a step $`\mathrm{\Delta }y=0.0025`$. The accuracy of this algorithm is better than 2 % in the entire range of $`k`$ discussed below. We evolve two initial conditions given by the McLerran-Venugopalan mv  (MV) and Golec-Biernat–Wüsthoff gbw  (GBW) model respectively:

$`N_{MV}^{Q_s}`$ $`=`$ $`1\mathrm{exp}\left[{\displaystyle \frac{Q_s^2r^2}{4}}\mathrm{ln}\left({\displaystyle \frac{1}{r^2\mathrm{\Lambda }_{QCD}^2}}+e\right)\right],`$ (4)
$`N_{GBW}^{Q_s}`$ $`=`$ $`1\mathrm{exp}\left[{\displaystyle \frac{Q_s^2r^2}{4}}\right],`$ (5)
where $`\mathrm{\Lambda }_{\mathrm{QCD}}=0.2\mathrm{GeV}`$. For momenta $`kO(1\mathrm{GeV})`$, the sensitivity on the infrared cut off $`e`$ is negligible. The amplitudes $`N_{MV}`$ and $`N_{GBW}`$ are similar for momenta of order $`Q_s`$, but differ strongly in their high $`k`$ behavior; $`\varphi _{GBW}(k)`$ decays exponentially while $`\varphi _{MV}`$ has a power-like tail $`1/k^2`$.

Fig.1 shows the evolution of $`h(k,y)`$ and $`\varphi (k,y)`$ for different initial conditions. The solutions for $`h(k,y)`$ quickly approach a universal soliton-like shape and do not change further except uniformly moving in $`k`$ on the logarithmic plot. The position of the maximum is the evolved value of the saturation momentum $`Q_s(y)`$. The solutions for different initial conditions and different rapidities scale as a function of the scaling variable $`\rho =k/Q_s(y)`$. The shape of the initial condition affects only the value of the saturation momentum $`Q_s(y)`$, but not the shape of the evolved function $`h(\rho ,y)`$. The $`y`$-dependence of $`h(\rho ,y)`$ is very weak: the function evolves fast towards a scaling form $`h(\rho )`$. As the rapidity changes between $`y=4`$ and $`y=10`$, the ratio $`h(\rho ,y_1)/h(\rho ,y_2)`$ varies by at most $`40\%`$ over three orders of magnitude of the scaling variable $`\rho `$. Similar behavior was found for $`\varphi `$ (results not shown). This is consistent with previous numerical works ab ; gbs .

To get a quantitative idea of the behavior of the scaling solution, we fitted the numerical solution of $`\varphi (\rho )`$ to two analytical expressions: $`s_1(\rho )=a\rho ^{2(1\lambda )}`$ and $`s_2(\rho )=a\mathrm{ln}(b\rho )\rho ^{2(1\lambda )}`$ for $`\rho >5`$. The functional form $`s_1`$ with $`\lambda =0.37`$ and $`\mathrm{ln}Q_sy`$ describes the scaling behavior of solutions of the linear BFKL equation i^2m . It was argued in Ref. mt  that $`s_2`$ with the same value of $`\lambda `$ and $`\mathrm{ln}Q_s\frac{2\chi (\lambda )}{1\lambda }y\frac{3}{2(1\lambda )}\mathrm{ln}y`$ accounts for the effects of nonlinearities in Eq. (1). We find that $`s_1`$ does not give an acceptable fit to $`\varphi (\rho )`$ in any extended range of $`\rho `$. For values of $`\rho `$ between $`1`$ and $`10^3`$ the value of $`\lambda `$ varies between 0.39 and 0.46. This is in contrast to the BFKL equation, where we find numerically that $`s_1`$ with $`\lambda =0.37`$ does indeed approximate the solution over several orders of magnitude with very good accuracy (results not shown). On the other hand, for $`5<\rho <1000`$, $`s_2`$ gives a good fit with $`\lambda =0.32`$. If, following mt  we fix $`\lambda =0.37`$, the fit is still good.

To study the effect of the evolution on the Cronin enhancement, we consider two initial conditions, $`N_{MV}^q(r)`$ and $`N_{MV}^Q(r)`$ with $`q^2=0.1\mathrm{GeV}^2`$ and $`Q^2=2\mathrm{GeV}^2`$. Since $`q\mathrm{\Lambda }_{QCD}`$ and $`Q`$ is of order of the estimated saturation momentum for a gold nucleus kln , this choice mimics the gluon distributions of a proton and of a nucleus respectively. At large transverse momenta the ratio of the corresponding Fourier transforms is given by the ratio of the saturation momenta,

$$\frac{h^Q(k,y=0)}{h^q(k,y=0)}=\frac{Q^2}{q^2}=A^{1/3}.$$
(6)
This relation also holds for $`\varphi `$. As discussed in bkw ; kkt , these initial conditions exhibit Cronin enhancement, namely $`\frac{h^Q(k,y=0)}{A^{1/3}h^q(k,y=0)}>1`$ for $`kQ`$. We solve the BK equation with these two initial conditions and construct the ratio $`R(k,y)=h^Q(k,y)/A^{1/3}h^q(k,y)`$ and the corresponding ratio for $`\varphi `$, see Fig. 2. The initial Cronin enhancement at rapidity $`y=0`$ is seen to be wiped out very quickly by the evolution. Within less than half a unit of rapidity $`y`$ the ratios show uniform suppression for all values of transverse momentum. The observed behaviour persists if different amounts of Cronin enhancement are included in the initial condition.

As seen in the lower panel of Fig. 2, the Cronin enhancement also disappears rapidly with rapidity when gluon distributions are evolved according to the linear BFKL equation. Qualitative differences between the BFKL and BK dynamics are only visible at momenta $`k<Q_s`$, where saturation effects are important. For larger momenta $`k`$, the ratios are very similar for linear and non-linear QCD evolution. We thus conclude that the wiping out of the initial enhancement is primarily driven by the linear BFKL dynamics which is contained in the BK equation as well.

For the evolved gluon distributions determined above, we have calculated the yield of produced gluons in p–A and A–A collisions at central rapidity according to the factorized expressions glr

$`{\displaystyle \frac{dN_{pA}}{dyd^2pd^2b}}`$ $``$ $`{\displaystyle \frac{1}{p^2}}{\displaystyle d^2kh^q(y,k)h^Q(y,pk)},`$ (7)
$`{\displaystyle \frac{dN_{AA}}{dyd^2pd^2b}}`$ $``$ $`{\displaystyle \frac{A^{2/3}}{p^2}}{\displaystyle d^2kh^Q(y,k)h^Q(y,pk)}.`$ (8)
From these spectra, we calculate the $`p`$\- and $`y`$-dependent ratios

$`R_{pA}`$ $`=`$ $`{\displaystyle \frac{\frac{dN_{pA}}{dyd^2pd^2b}}{A^{1/3}\frac{dN_{pp}}{dyd^2pd^2b}}},R_{AA}={\displaystyle \frac{\frac{dN_{AA}}{dyd^2pd^2b}}{A^{4/3}\frac{dN_{pp}}{dyd^2pd^2b}}}.`$
As seen in Fig. 3, the non-linear BK evolution quickly wipes out any initial Cronin enhancement not only on the level of single parton distribution functions but also on the level of particle spectra. We have checked that this behaviour is generic by evolving different initial conditions corresponding to different initial amounts of enhancement. We note that in calculations of the gluon production in p–A in the eikonal approximation kovchegov ; km ; braun , the gluon distribution $`h`$ rather than $`\varphi `$ enters the right hand side of (7), but no similar statement exists of nucleus-nucleus collisions. We have checked that the results using $`\varphi `$ are very close to those shown in Fig. 3. More generally, the expressions (7) and (8) are based on rather strong approximations discussed in Ref. bkw . However, our conclusion about the disappearance of Cronin enhancement during QCD evolution is likely to persist in more refined ways of calculating particle spectra, since it is rooted directly in the rapidity dependence of gluon distributions.

We now comment on a recent formal argument kkt  which \- in contrast to our numerical findings - suggests that Cronin enhancement survives the non-linear evolution. It is based on the observation that at very short distances $`r0`$, the dipole amplitude $`N(r)`$ is not affected by evolution. Thus, the integral of the gluon distribution function $`\varphi `$ over the transverse momentum is expected to be rapidity independent,

$$d^2k\varphi (k)=\frac{1}{r^2}N(r)|_{r=0}.$$
(9)
One thus obtains the sum rule

$$d^2k\varphi _A(k,y)=A^{1/3}d^2k\varphi _p(k,y)$$
(10)
valid for any rapidity, since it is satisfied by the initial condition $`\varphi _{MV}`$. Since the nonlinear evolution leads to the depletion of the gluon distribution $`\varphi _A(k)`$ relative to $`A^{1/3}\varphi _p(k)`$ at low momenta, it must follow that in some range of momenta this effect is compensated by enhancement of $`\varphi _A`$. Similar arguments can be made about the distribution $`h`$ and also about the gluon yield in p–A and A–A collisions.

However, this argument breaks down since the quantity defined in Eq. (9) is infinite. As such Eq. (10) relates only the divergent parts of the integrals which are dominated by ultraviolet, and carries no information about possible behaviour at finite momentum. To be more specific, we use the scaling property $`\varphi (k,y)=\varphi (k/Q_s(y))`$ of the solution of the BK equation established above. It is known that the ratio of the saturation momenta for any two solutions is preserved by the BK evolution ab ; gbs ; i^2m . For our two solutions representing a nucleus and a nucleon, this implies $`\frac{Q_s^A(y)}{Q_s^p(y)}=A^{1/6}`$. We now rewrite the sum rule (10) by regulating the divergent integrals with a large but finite UV cutoff $`aQ_s^A(y)`$,

$`{\displaystyle _0^{a^2(Q_s^A)^2}}d^2k\varphi _A(k,y)=(Q_s^A)^2{\displaystyle _0^{a^2}}d^2\rho \varphi (\rho )`$ (11)
$`=A^{1/3}(Q_s^p)^2{\displaystyle _0^{a^2}}d^2\rho \varphi (\rho )=A^{1/3}{\displaystyle _0^{a^2(Q_s^p)^2}}d^2k\varphi _p(k,y).`$
In the formal limit $`a\mathrm{}`$, we recover Eq. (10). However, since $`Q_s^AQ_s^p`$, the regularized sum rule (11) is easily satisfied even if the nuclear distribution is suppressed relative to that of a single nucleon uniformly at all momenta. Thus, the sum rule (10) carries no information about either presence or absence of the Cronin enhancement.

In summary, we have found that the nonlinear QCD evolution to high energy is very efficient in erasing any Cronin type enhancement which may be present in the initial conditions. For realistic initial conditions this disappearance occurs within half a unit of rapidity. We note that in our units the evolution from $`130`$ GeV to $`200`$ GeV corresponds to $`\delta y0.1`$ for $`\alpha _s=0.2`$, and thus is not sufficient to completely eliminate an initial enhancement at central rapidity. For forward rapidity, $`\delta y`$ is greater. The evolution to the LHC energy corresponds to $`\delta y1`$. Thus the BK evolution suggests the reduction of the Cronin effect in d–Au for forward rapidities at RHIC and it predicts the disappearance of the Cronin effect for p–A collisions at LHC.

We thank the CERN TH Division (JLA and AK) and the INT Seattle for its hospitality and the DOE for partial support (NA, AK, CAS and UAW). JLA is supported by MECD of Spain, grant no. AP2001-3333, and CAS by a Marie Curie Fellowship no. HPMF-CT-2000-01025 of the European Community TMR program. Useful discussions with R. Baier, E. Iancu, D. Kharzeev, Y. Kovchegov, P. Jacobs and L. McLerran are gratefully acknowledged. We thank M. A. Braun for discussions and a numerical cross-check.


text:
# On Heavy Element Enrichment in Classical Novae

## 1. Introduction

Classical novae are a manifestation of thermonuclear runaways in accreted hydrogen shells on the surfaces of white dwarf stars in close binary systems (see, e.g., the review by Gehrz et al., 1998). Compelling observational data indicate that the material ejected by some classical novae can be significantly enriched in C, N, O, and Ne, by $`30\%`$ by mass (Livio & Truran 1994; Gehrz et al. 1998). Since the abundance of CNO catalysts constrains the rate of energy release, such high levels of CNO enrichment are required in the fastest novae (Truran, 1982), for which the hydrogen burning reactions (via the CNO cycle) increase the entropy of the envelope faster than it can adjust. It was recognized early that if some of the underlying white dwarf matter111By “white dwarf,” we mean the predominantly C/O substrate; we shall refer to the accumulated H/He layer as the “atmosphere” or “envelope.” could be mixed into the accreted envelope prior to the final stages of the thermonuclear runaway, then the explosion would be more energetic and the ejecta more enriched in CNONe elements (Starrfield et al., 1972). One-dimensional models that best reproduce observations typically accrete material “seeded” with a super-solar composition (e.g., Hernanz et al., 1996; Starrfield et al., 1998).

The question of how the accreted envelope is enriched has challenged theory for several decades (see, e.g. the review by Livio & Truran, 1994). For very slow accretion ($`10^{10}M_{}\mathrm{yr}^1`$), the downward diffusion of H into the underlying C/O or O/Ne layers (Prialnik & Kovetz, 1984; Kovetz & Prialnik, 1985) could trigger ignition in the H diffusive tail and drive convective mixing of heavy elements into the envelope during the early stages of runaway. It is unclear, however, whether there is time at higher accretion rates for this process to be relevant. The possibility of mixing via convective overshoot was first considered by Woosley (1986). Numerical simulations of convective penetration and mixing at the core-envelope interface have been carried out in both two dimensions (Glasner, Livne, & Truran, 1997; Kercek, Hillebrandt, & Truran, 1998) and three (Kercek, Hillebrandt, & Truran, 1999). From comparisons of two- and three-dimensional simulations and a careful resolution study, Kercek et al. (1998, 1999) concluded that convective penetration would not significantly enrich the accreted envelope.

The general problem of shear mixing was considered by many (Kippenhahn & Thomas, 1978; MacDonald, 1983; Fujimoto, 1988), but to date no self-consistent calculations of these effects have been performed. Kippenhahn & Thomas (1978) considered shearing from differential rotation in the envelope and estimated the amount of mass mixed by assuming that the envelope was marginally stable according to the Richardson criterion. MacDonald (1983) considered non-axisymmetric perturbations and the resulting redistribution of matter and angular momentum. Fujimoto (1988) emphasized the role of barotropic and baroclinic instabilities in transporting angular momentum through the envelope. A variation on this was envisaged by Sparks & Kutter (1987) and Kutter & Sparks (1989), who suggested that convection just prior to the runaway would transport angular momentum inward and lead to a large horizontal shear above the C/O white dwarf.

While these studies examined the generation of shear in the envelope, the details of how C/O is mixed into the envelope and the quantity of mass mixed, as functions of the shear profile, were not resolved. In this paper, we examine a promising mechanism for effecting the entrainment: a resonant interaction between large-scale shear flows in the accreted envelope and interfacial gravity waves (Rosner et al., 2001). The greater compositional buoyancy in the C/O white dwarf means that the interface sustains gravity waves. Miles (1957) showed that in the presence of a shear flow (i.e., a “wind”), gravity waves with a group velocity matching a velocity in the shear flow are resonantly amplified. These waves eventually form cusps and break. When the waves break, they inject, analogously to ocean waves, a spray of C/O into the H/He atmosphere.

The source of the shear could arise from a number of mechanisms, such as those described above. A self-consistent calculation of the accretion boundary layer and the wind profile in the atmosphere is far beyond the scope of current simulations. Indeed, as we explain in § 2, the resolution needed to simulate the interaction is so fine as to prohibit the direct simulation of this effect within a macroscopically large (i.e., the entire atmosphere of the white dwarf) computational domain. What this paper provides is a scaling relation (eq. ) for the mixing rate per unit area, for a given wind velocity. This relation may be used to construct a sub-grid model for large-scale multi-dimensional simulations, which can in turn address how matter and angular momentum are transported over the surface.

In § 2 we present the two-dimensional models and demonstrate how the mixed mass depends on the velocity of the flow. We then, in § 3, incorporate this relation into a one-dimensional simulation and discuss two scenarios: that mixing only occurs if the convective zone reaches the C/O interface; and that the mixing occurs prior to the onset of convection. We implement these results as a sub-grid model in one-dimensional novae simulations performed with a stellar evolution code and discuss the implications, in § 4, for the ejected mass of C/O. Because the shear in the envelope prior to the onset of convection is unknown, we compute the accretion and runaway for several different masses of the mixed layers.

## 2. Breaking Gravity Waves

To calculate the amount of mass mixed by wave breaking as a function of time, we performed a suite of numerical simulations of wind-driven gravity waves. The initial configuration consists of two stably stratified layers of compressible fluids in hydrostatic equilibrium in a uniform gravitational field $`g`$. The upper fluid is composed of H/He and the lower fluid is composed of C/O. The equation of state is an ideal gamma-law of internal energy $`\rho ϵ=P/(\gamma 1)`$. We choose $`\gamma =5/3`$ and for the initial condition we construct an isentropic envelope in hydrostatic equilibrium. The density is then given by $`\rho (z)=\rho _i[1(\gamma 1)z/H]^{1/(\gamma 1)}`$, where $`\rho _{i=1}`$ is the density of the upper fluid just above the interface, $`\rho _{i=2}`$ is the density of the lower fluid just below the interface, $`H=\gamma ^1c_\mathrm{s}^2/g`$ is the pressure scale height, and $`c_\mathrm{s}`$ is the adiabatic sound speed. This choice makes the Brunt-Väisälä frequency zero, i.e., there are no internal gravity waves within each of the two layers, which isolates the interface modes for study.

We impose a wind in the upper fluid with velocity

$$U(z)=U_{\mathrm{max}}(1e^{z/\delta }),$$
(1)
where $`z`$ is the vertical direction and $`z=0`$ corresponds to the location of the initial interface; $`U_{\mathrm{max}}`$ is the asymptotic maximum velocity at $`z\delta `$, and $`\delta `$ is the length scale of the shear boundary layer. The form of the wind profile (eq. ) is chosen to model the boundary layer formed at the surface of the white dwarf, and is motivated by studies of the wind-wave resonant instability in oceanography (Miles, 1957). This instability is the principal mechanism that drives gravity waves and air-water mixing in oceans. In the limit $`\delta 0`$, one recovers the Kelvin-Helmholtz wind profile. Without a study of angular momentum transport in the accreting envelope, we do not know the form of the wind profile and in this work investigate the mixing for a range of $`\delta `$, with $`\delta /H1`$.

Alexakis, Young, & Rosner (2002) performed a linear stability analysis of the wave amplification for the novae problem, and in a forthcoming paper (Alexakis et al., in preparation) we will present extensive numerical studies of wave-wind interactions. Here we apply those results to the conditions at the base of the white dwarf envelope. The simulations are performed using the Flash code (Fryxell et al., 2000; Calder et al., 2000, 2002), a parallel, multi-dimensional, adaptive-mesh code for simulating compressible, inviscid flow. The computational grid was a uniform mesh of resolution $`1024\times 1024`$ with periodic boundary conditions along the sides of the box and hydrostatic boundaries (Zingale et al., 2002) at the top and bottom. The interface is perturbed with specified modes and a velocity perturbation added to the surface of the WD representing the motion of the gravity waves. The wind profile is specified as an initial condition and in particular is not forced. For this study, we set the Mach number $`\mathrm{𝖬𝖺}=U_{\mathrm{max}}/c_\mathrm{s}=0.5`$ and consider four values of $`\delta /H`$: 0.005, 0.01, 0.02, and 0.04. These values vary the reciprocal Froude number $`(\delta g/U_{\mathrm{max}}^2)^{1/2}=(\delta /H)/(\gamma \mathrm{𝖬𝖺}^2)`$, which measures the available kinetic energy for mixing.

The unstable modes have wavenumbers (Alexakis et al., 2002) $`kg/U_{\mathrm{max}}^2(1\rho _1/\rho _2)`$. We chose the size of our computational domain to be roughly equal to this maximum wavelength, so that we have 1024 grid points along the largest possible wave. In practice, smaller waves become dominant, with the wavelength of the dominant mode depending on the wind profile. In all cases this dominant mode is well-resolved. Figure 1 shows fully developed waves breaking, generation of the “spray”, and mixing of the white dwarf substrate up into the atmosphere. Mixing is not only effected by the overturning of the wave but also by the wind breaking off the tip of the wave. Figure 2 shows contours (*black lines*) of the carbon mass fraction at 0.49, 0.20, and 0.02, respectively, for the simulation depicted in Figure 1, but at a later time. The contour at 0.49 corresponds to the carbon mass fraction of the underlying white dwarf, while the smaller contour values indicate how far outward into the accreted material the white dwarf substrate is mixed. The contour at 0.2 delimits the region that contains most of the initial enrichment.

Figure 3 shows the surface mass density (the mass of C/O in the mixed layer per unit area), $`dM_{\mathrm{CO}}/dS`$, averaged over the horizontal direction. The mixed layer is defined here as the region in which the carbon mass fraction is between 0.49 and 0.01. The amount of mixed mass depends on these delimiters; reducing the upper limit to 0.4 decreases the mixed mass by less than a factor of 2. The C/O is mixed rapidly until it saturates; further mixing occurs on diffusive timescales. Our two-dimensional simulations show, for the range of parameters examined, that the total mass of white dwarf material that becomes mixed is independent of the lengthscale $`\delta `$. The rate of mixing, i.e., the initial slopes of the curves in Figure 3, does show, however, some dependence on $`\delta `$. Using dimensional analysis and the numerical results, we find that for a fixed density ratio $`\rho _1/\rho _2=0.6`$ (appropriate for the degenerate plasma prior to ignition) the total mass per unit area, $`dM_{\mathrm{CO}}/dS`$, mixed into the accreted H/He saturates at

$$\frac{dM_{\mathrm{CO}}}{dS}=\alpha \frac{U_{\mathrm{max}}^2}{g}\rho _{\mathrm{WD}},$$
(2)
where $`\alpha `$ is a non-dimensional constant222In general, $`\alpha `$ depends on $`\rho _2/\rho _1`$ and possibly $`\mathrm{𝖬𝖺}`$. A parameter study, outside the scope of this paper, is necessary to determine this dependence. that we determine from the simulations to be $`\alpha 0.6`$ (see Figure 3). The timescale to reach saturation is far shorter than the timescale of either the accretion phase ($`>10^4\mathrm{yr}`$) or the pre-peak convective phase ($`100\mathrm{yr}`$) of a typical classical nova. For example, if $`\delta /H=0.01`$ and $`\mathrm{𝖬𝖺}=0.5`$, then the saturation timescale is of order seconds for our nova setup.

## 3. One-Dimensional Nova Models

We now incorporate our simulations of the wave breaking and mixing into simulations of the thermonuclear burning of a nova. In order to explore the global properties of this local mixing mechanism, we compute several one-dimensional models of novae with a modified version of the KEPLER stellar evolution code (Weaver, Zimmerman, & Woosley, 1978) that includes a large network of $`100`$ light isotopes (Rauscher et al., 2002) implicitly coupled into the stellar structure solver. This allows us to follow throughout the entire convective envelope the radioactive decays of nuclei that are formed at the base of the burning zone. The underlying white dwarf has a mass of $`1.0M_{}`$, a radius of $`5000\mathrm{km}`$, a luminosity of $`10^{31}\mathrm{ergs}\mathrm{s}^1`$, and is composed of a 50%/50% $`{}_{}{}^{12}\mathrm{C}`$/$`{}_{}{}^{16}\mathrm{O}`$ mixture. We retain the outer $`0.005M_{}`$ on the computational grid to follow the thermal inertia of the outer WD layers and allow the model to relax until it has a constant luminosity before we start the accretion. The accretion rate is $`10^9M_{}\mathrm{yr}^1`$, and in all but one model (see item 3 below), the accreted material consists of 70 % $`{}_{}{}^{1}\mathrm{H}`$, 28.7 % $`{}_{}{}^{4}\mathrm{He}`$, 0.3 % $`{}_{}{}^{12}\mathrm{C}`$, 0.1 % $`{}_{}{}^{14}\mathrm{N}`$, and 0.9 % $`{}_{}{}^{16}\mathrm{O}`$, by mass. Convection is modeled using the Ledoux criterion for stability and mixing length theory. We assume no convective overshooting and that semiconvection is too slow to cause mixing.

We investigate three scenarios, as summarized in Table 1, for generating the gravity wave induced mixed layer.

1. In the first scenario, the only shearing considered is that from the convective cells driving a wind at the interface between the H-rich atmosphere and the C/O substrate. This is the scenario envisaged by Rosner et al. (2001), in which the C/O, after being entrained, is then distributed throughout the convective zone. Figure 4 summarizes the evolution of the accreted layer. After accretion for about $`6.9\times 10^4\mathrm{yr}`$ a convective zone forms about $`2.5\times 10^5M_{}`$ above the WD interface. The total envelope mass accreted at ignition is $`7\times 10^5M_{}`$, in good agreement with the estimates of Fujimoto (1982a, b). Within $`100\mathrm{yr}`$ from the onset of convection, the convective zone extends upward to the surface and downward to about $`5\times 10^6M_{}`$ above the interface, at which time the runaway reaches its peak rate of energy generation. This peak phase evolves on a timescale of $`1\mathrm{h}`$, but the convection does not reach the WD interface. Only days later, when the nova envelope is already significantly expanding, does the burning layer reach the WD. This downward movement, in mass, of the burning layer is mediated by heat conduction into these deeper H-rich layers. The convective zone also moves deeper, but never quite reaches the WD interface333At the very bottom of the H shell the nuclear luminosity from H burning can never become large enough to drive convection, as it goes to zero at the interface; the only possibility for mixing to occur is if a steep temperature gradient arises from the thermal inertia of the heated WD core below an expanding H envelope.. As a result, there is no injection of C/O into the H-rich envelope, by construction.

The envelope expands to large radii ($`>10^{12}\mathrm{cm}`$, large enough to engulf the secondary) by the end of the simulation, but the result will only be a slow nova, in rough agreement with the semi-analytical calculations of Fujimoto (1982a, b). We do not find any contact of the convective region with the WD substrate; as a result, we do not expect that any convection-driven wave mixing would occur, unless there is a large redistribution of angular momentum in the accreted envelope (Sparks & Kutter, 1987; Kutter & Sparks, 1989).
2. In the second scenario, the shearing originates from a wind at the base of the accreted envelope blowing across the underlying C/O substrate. We assume that the wind persists throughout the H/He layer with velocity sufficient to drive mixing on a timescale much less than that to accrete a critical mass of fuel. In this case the mixed layer is generated prior to ignition. We also assume a linear (in Lagrangian mass coordinate) gradient in the mass fraction of C/O between the WD and the accreted envelope. Because the shear profile, and hence the amount of mass mixed, is unknown, we consider a range of mixed masses (see Table 1). Our ignorance of the fluid motions in the envelope and substrate prevents us from saying where and how the accreted material spreads over the surface, and so we cannot determine the actual shear at the base of the envelope. In all cases that we consider, the shear velocity is subsonic, which is the range of validity of the 2-d simulations. For such velocities, the mixed layer is always thinner than a pressure scale height (eq. ).

Figure 5 shows the case with the largest pre-enrichment: the mixed material comprises $`4.6\times 10^{28}\mathrm{g}`$ (i.e., a total of $`2.3\times 10^{28}\mathrm{g}`$ of WD material is being mixed with H-rich material). This is about 25% of the envelope mass, which is less than that generated by the calculation described in section 2. From equation (2), a 25% enrichment corresponds to $`\mathrm{𝖬𝖺}=0.4`$ if the velocity were uniform over the surface. We note that this velocity is much less than Keplerian, $`U_{\mathrm{max}}0.05(GM/R)^{1/2}`$, as expected in the envelope well below the accretion disk boundary layer (see, e.g., Popham & Narayan, 1995). If the velocity were not uniform over the surface, then a higher maximum velocity would be required to inject the same percentage of the total envelope mass.

The mixed layers are added to the surface using the same accretion rate as the rest of the envelope ($`10^9M_{}\mathrm{yr}^1`$); to prevent a runaway while these layers are being added, we suspend nuclear energy generation during the accretion of the first $`5\times 10^{28}\mathrm{g}`$. About 50 $`\mathrm{yr}`$ prior to peak energy generation a semiconvective region forms in the outer 60 % (by mass) of the enriched layer; about 8 months before the peak a convection zone starts at the base of this semiconvective layer (at $`7\times 10^6M_{}`$ above the WD interface) that extends upward to the surface and downward to about $`5\times 10^7M_{}`$ above the interface. In contrast to the first case, the pre-mixed envelope with the largest enriched mass ignites at a smaller accumulated mass, $`5\times 10^5M_{}`$.

The energy generation in the enriched layer is dominated by $`{}_{}{}^{12}\mathrm{C}(\mathrm{p},\gamma )^{13}\mathrm{N}(\beta ^+)^{13}\mathrm{C}`$. The peak energy generation timescale is about $`3\mathrm{min}`$—20 times faster than in the first case—with a specific energy generation rate about 20 times higher. Figure 6 shows the specific nuclear energy generation in the accreted envelope, with the time centered about the peak energy generation rate. The larger rate of energy generation roughly corresponds to the $`20`$ times greater metal enrichment. In this scenario the entire 25 % of pre-enriched C/O material is spread throughout the hydrogen envelope during the runaway. This leads to a a corresponding enrichment in the ejecta. The ejecta are unbound, i.e., they have a positive velocity at very large radii.

As in the first case, the convection never reaches the WD interface (although the base of the convective zone is much deeper), so there is no additional convection-induced wave mixing of the C/O substrate with the hydrogen-rich envelope. We did experiment with pre-runaway, pre-enriched layers as small as 1/10,000 the mass used in the case described above (Table 1). As before, we disregarded nuclear burning for the first $`10^{28}\mathrm{g}`$ of accreted material, with the exception of the case with $`3.77\times 10^6M_{}`$ of enrichment, where we took twice that mass to accommodate the entire enriched layer. The difference in masses for which nuclear burning is suspended is the reason for the case with the largest pre-enriched mass, $`1.15\times 10^5M_{}`$ having a slightly larger envelope mass at runaway than that of the case with $`3.77\times 10^6M_{}`$ of enriched mass.

In none of these cases did the convection reach the WD interface. Eventually, for the case of an enriched mass of $`10^9M_{}`$, the runaway behaves as in the non-enriched case: the energy release rate in the mixed layers is too small to drive convection above and heat conduction cools it efficiently so that the runaway does not occur in its vicinity. The slightly higher envelope masses obtained in the limit of small enrichment is likely to be an artifact of our suspension of nuclear energy generation for the first fraction of the accretion phase.
3. In a third scenario we assume accretion of material enriched to a mass fraction of 50 % in C/O. Allowing for nuclear energy generation instantaneously in this material leads to a rapid runaway after only a thin layer has been accumulated at the surface of the star. If this magnitude of enrichment was due to instantaneous mixing with a colder WD substrate, then a lower temperature in the enriched material and thus a significantly later runaway could result. We include this run to demonstrate how qualitatively different the accretion of pre-enriched material behaves as compared to first accreting unenriched material and then, only later after it has settled, mixing the accreted layers with the C/O substrate.

## 4. Implications and Summary

Using a constant wind profile that blows across the surface of a white dwarf, we have performed a two-dimensional parameter study of the mechanism proposed by Rosner et al. (2001). Our primary results of the mixing rate and the maximum mixed mass, equation (2), suggest that this process can mix about $`10^6`$$`10^5M_{}`$ of the underlying C/O into the hydrogen-rich envelope. From this, we investigated two scenarios: that the wind is the result of convection, and that the wind originates during the accretion phase.

We find that if no enrichment occurs prior to the onset of convection, then the convective zone does not reach the C/O interface, and no additional mixing occurs (in the one-dimensional model) in the absence of convective overshoot. The result in this case will be a slow nova, with little enrichment of the ejecta. In contrast, an envelope with a mixed layer at the C/O interface, consistent with the scalings from high-resolution numerical studies, provides a more violent runaway and the ejecta are enriched by $`25\%`$ in C/O, consistent the enrichment observed in some nova ejecta. Such an event does require a strong shear velocity, but the saturation amount of mass mixed is roughly independent of the shape of the shear profile, so long as its thickness is much less than a pressure scale height. Our runs, outlined in Table 1, show that when the metals are concentrated at the base of the accreted envelope, the amount of mass accreted prior to runaway is *larger* than if the C/O were uniformly distributed over the envelope, i.e., if the white dwarf were to accrete material with a supersolar metallicity. The reason is that the opacity of the envelope is larger than when the C/O are concentrated at the base.

There are a number of issues that we have not yet addressed or investigated. First, our results might depend on the dimensionality of the system. We would therefore expect that although the generation of the gravity waves is captured by two-dimensional dynamics, the energy cascade of the waves and the advection of spray and vorticity will be different in a three dimensions. Second, this work investigates only the density ratio for an accreted envelope in thermal equilibrium. During the early phase of accretion and during the runaway the density ratio will likely be different. Third, although our parameter study covers more than one order of magnitude in $`\delta `$, it is important to know if our results still hold at even larger values. As $`\delta `$ becomes a respectable fraction of $`H`$, we expect that the stratification of the envelope, and in particular the effect of a non-zero Brunt-Väisälä frequency, will become important. All of these differences can affect the mass and thickness of the mixed layer, and the investigation of these issues is the subject of ongoing work.

In our one-dimensional nova simulations, we assumed that a convective roll will only entrain the mixed layer if the base of the convective zone (as computed from a mixing-length formalism) reaches the interface. Future multidimensional studies can inform us as to how convection interacts with the interfacial gravity waves, and to whether this assumption is in fact correct. Moreover, our sub-grid model, when incorporated into a one-dimensional calculation, does not account for spatial variations of the wind. Such a variation would occur, for example, if the fuel accretes non-uniformly over the surface of the star and then spreads. In reality the amount of mass mixed will vary over the surface of the star. Our one-dimensional calculations must of necessity take the mixed mass as a free parameter. Because the amount of C/O entrained depends quadratically on the velocity in the wind profile, our calculation underscores the need for detailed simulations of the shear profile of an accreting white dwarf. These issues are clearly important and are the next steps to pursue.

We thank Lars Bildsten and Stan Woosley for helpful discussions. This work is supported in part by the U.S. Department of Energy under Grant No. B341495 to the Center for Astrophysical Thermonuclear Flashes at the University of Chicago. A. H. is supported by the Department of Energy under contract W-7405-ENG-36. L. J. D. acknowledges support by the Krell Institute CSGF. M. Z. is supported by DOE grant number DE-FC02-01ER41176 to the Supernova Science Center/UCSC. P. M. R. acknowledges support from the University of Illinois at Urbana-Champaign and the National Center for Supercomputing Applications (NCSA). K. Olson acknowledges partial support from NASA grant NAS5-28524.