Quantum Pairing of Impurities in Quantum Crystals

We calculated the time of pairing by quantum diffusion of ortho-H2 impurities in solid para-H2. The important feature of the pairing process is a strong directional bias associated with the dependence of the hopping rates on energy mismatches caused by the interaction of the pairing particles. This bias at moderate temperatures is against a mutual approach of particles and creates a "kinetic barrier. " At lower temperatures, the corresponding diffusion mechanism freezes out, which leads to a rapid increase in pairing rates. This explains a well-developed, experimentally observed maximum in the pairing time as a function of temperature: a maximum that exists in spite of a monotonic temperature dependence of individual hopping rates. Our results are in good agreement with experimental data.


I. INTRODUCTION
Quantum diffusion is, basically, the tunneling of impurity particles in a host lattice at low temperatures. This immediately separates quantum diffusion from the classical case, which is thermally activated and depends exponentially on temperature.
The main quantitative and even qualitative difference between quantum diffusion and, for example, electron tunneling in metals and semiconductors is associated with extremely small values of tunneling frequencies and, therefore, narrow bandwidths for quantum diffusion of impurity particles. This, in turn, makes quantum diffusion very sensitive to any inhomogeneities of the potential relief. Even very small inhomogeneities may change the diffusion rates by several orders of magnitude.
Usually, inhomogeneities in the potential are random.
Then the global diffusion of quantum impurities is isotropic, though locally the instantaneous hopping probabilities are position and time dependent and are anisotropic. However, there are two notable exceptions when even the overall diffusion exhibits a strong directional bias. The first one is obvious and occurs in the presence of nonrandom external fields such as inhomogeneous deformation fields, electric fields, or gradients of magnetic fields for charged or magnetic impurities. The second exception is more subtle and is associated not with an external field, but rather with some particular internal processes. Here we have in mind, first of all, the quantum pairing of impurity particles in quantum crystalsthe formation of the nearest-neighbor pairs of impurities by means of quantum diffusion. In this case the strong bias in quantum diffusion is associated with the nonrandom distortion of a potential relief caused by the mutual interaction of the pairing particles.
The experimentally observable characteristic of the quantum pairing is the characteristic time of mutual approach of two impurity particles: the pairing time. The pairing time is an inherent characteristic of numerous processes such as clusterization of impurities in helium or hydrogen crystals, recombination of atomic hydrogen impurities in molecular hydrogen or metals, fusion of iso-topes of hydrogen in a metal matrix, etc. The understanding of probabilities and rates for all such processes requires a careful analysis of pairing times.
Usually, estimates of pairing times are done very roughly using some effective isotropic diffusion coefticient. However, these estimates may be off by several orders of magnitude because of the strong directional bias. This can affect the predictions for various reaction rates.
The cause of this strong bias is that the quantum pairing, at least in its latest stages, is not a random isotropic diffusion. When the distance r between the pairing particles is small, the major potential inhomogeneity for quantum tunneling is associated with the interaction U(r) of pairing particles. In many cases the tunneling frequency for quantum diffusion (i.e. , the bandwidth) is very small, coo &( U(r) (here and below, 6=k~= 1), at least for small distances. Then the final stages of pairing cannot be described by an isotropic random diffusion, but correspond to a highly anisotropic diffusion motion with a very strong bias induced by the interaction field U(r). This means that the characteristic diffusion time becomes very different from the usual unbiased ones if, as it is very often the case, the final stages of pairing account for a considerable fraction of the pairing time.
In this paper we report the results of a numerical analysis of the time of quantum pairing. The results, of course, are very sensitive to the type of a bias imposed by interaction of pairing particles, and the problem cannot be addressed in a general form. As an example, we have chosen a problem of quantum pairing of ortho-H2 impurities in the crystals of para-H2. This pairing process is best studied experimentally by the analysis of NMR signals from single impurities and nearest-neighbor pairs (see the review in Ref. l). An additional interest to this particular example of quantum pairing is caused by the following feature of the experimental results. The temperature dependence of the pairing time' shows a pronounced maximum at very low temperatures. This seems very surprising since all reasonable diffusion mechanisms lead to the diffusion rates and diffusion coefficients which are either monotonically de-47 2573 1993 The American Physical Society

II. PAIRING AND QUANTUM DIFFUSION OF ORTHO-H IMPURITIES
The tunneling frequency for ortho-H2 impurities is extremely small so that any imperfection in potential reliefenergy mismatches 5E between the neighboring lattice sitesprevents coherent tunneling of impurity particles. In the case of quantum pairing, the nonrandom deformation of the potential relief is caused by the electric quadrupole-quadrupole interaction U(r) = UQ(n)(a Ir ), where a is the lattice constant and U0(n) -1 K depends on the orientation of the pair of impurity molecules with respect to the crystalline axes. Therefore the energy mismatches hindering the tunneling of a molecule from the lattice site i to the neighboring site j are 5E; =U 5 5 creasing with decreasing temperatures or have already saturated. A standard response to such an experimental anomaly would be the suggestion of some new mechanism of quantum diffusion with a similar maximum. However, as far as we know, such attempts have been not very successful and have led to the pairing times with a saturation instead of a temperature maximum.
A different idea, not requiring any exotic diffusion mechanisms was put forward by one of the authors. The idea was that a strong directional bias due to anisotropy of the diffusion trajectories in the final stages of pairing may lead to a nonmonotonic temperature dependence of the pairing time even if the individual hopping rates depend monotonically on the temperature. In particular, it was pointed out that the temperature affects not only the hopping rates, but also the bias and relative probabilities of jumps in different directions, thus changing drastically the characteristic numbers of individual jumps necessary for pairing. The pairing was explained by an interplay of temperature-independent (tunneling with spontaneous emission of phonons ) and temperature-dependent (tunneling accompanied by scattering of thermal phonons ) diffusion mechanisms. Both types of hopping have strong directional biases: The temperature-independent one has a higher probability in the direction of a mutual approach, while the other one has a higher probability in the opposite direction. Therefore, freezing out of the temperature-dependent diffusion mechanism can actually speed up the pairing process, which explains the maximum on a temperature dependence of the pairing time.
The arguments suggested that this might be the source of the experimental anomaly, but the supporting analytic calculations corresponded to oversimplified models and could not be considered a reliable explanation of experimental results and a proof of the importance of a diffusion bias in quantum pairing.
In this paper we present the results of detailed numerical calculations which prove that the directional bias may cause the observed maximum even in the case of monotonic hopping rates.
If the energy (1) exceeds the tunneling frequency~0, then simple tunneling is impossible.
The impurity still can tunnel from i to j despite the large (on the scale of co0) energy mismatch if the excessive energy 5E, . is taken care of by phonons. There are two possible phonon processes: spontaneous emission of phonons with the frequency co=5E; (Ref. 4) and inelastic scattering of thermal phonons.
The former process is possible only for 5E; )0, while the latter is available at any sign of 5E, -. Another difference between these two types of phonon processes is the dependence of their probabilities on the energy mismatch. While the probability of the spontaneous emission of phonons increases with frequency (i.e. , 5E; ) proportionally to the cube of the frequency, the probability of the inelastic scattering drops as 1/5E, -. Therefore the mutual diffusion of two impurities is assisted by the phonon scattering at the moderate distances between impurities, while the final stage of pairing is dominated by the (irreversible) phonon emission.
The inverse time of the phonon-assisted tunneling of impurities from i to j is (2) =~o jJ (3) is by several orders of magnitude shorter than (2), which is quadratic in co0. The characteristic pairing time r(T) should be expressed via the individual hopping rates (2) and (3) assuming a random initial distribution of impurities.
III. AVERAGE PAIRING TIMES:

COMPUTATIONAL MODEL
We are going to solve the following diffusion problem.
Initially, at t =0, the impurity particles are distributed randomly through the host lattice. Then the impurities start to move with the diffusion rates (2) (the energy mismatches should be determined using the current distribution of impurities). With time, the number of nearest-neighbor pairs will increase, while the number of singles will go down. We will assume that once the nearest-neighbor bond has been formed it cannot be broken (a very strong bond, or some reaction takes place).
We will neglect the motion of pairs (which is very slow anyway) and the formation of larger clusters, which should be important only in the advanced stages of pairing. We will calculate the pairing times as a function of , where 0 is the Debye temperature and g, g, are some constants which are not known exactly and which should be used as fitting parameters (g) 12; see below). If the energy mismatch is much smaller than coo, then the tunneling time relative initial positions and extract the overall characteristic time. The calculation is tailored to describe the NMR (and, to some extent, thermal conductivity) experiments' in solid para-H2 which follow the time changes in the numbers of singles and nearest-neighbor pairs.
The experiments are performed at rather low impurity concentrations less than 1%. This strongly affects the accuracy of computations. The problem is twofold. First, at lower concentrations the calculations have to be done for larger lattices. Second, the ratio of two diffusion rates in Eq. (2) is proportional to oE, the fifth power of the energy mismatch (1). Therefore, when the distance between the pairing particles changes (in atomic units) from, say, 10 to 1, the ratio of the rates (2) changes by the factor 10 . As a result, the computation involves large matrices with elements that are different by many orders of magnitude. The severity of roundoff problems increases as the concentration decreases.
In this paper we calculate the pairing time for two impurities averaged over the initial random uniform distribution of impurities. In a two-particle setting, it does not matter which of the particles makes a step; what matters is a change in a relative positions of particles before and after each step. Therefore we will consider one of the impurities as immobile (the "origin"), while all the dynamics will be ascribed to the second one (the "impurity").
We will compute the average time of the approach of the impurity to the origin.
The main limitation of the present calculation concerns the properties of the rest of the impurities which we assume to remain in the fixed positions forming a periodic superlattice. The single impurity diffuses among all these fixed particles and can approach and pair with each one of them.
Technically, it is done by considering a single threedimensional (3D) cell containing one origin and one impurity and imposing periodic boundary conditions. The potential energy for the impurity within the cell, and therefore the mismatches, includes the electric quadrupole-quadrupole interaction of the impurity with all of the particles forming the superlattice, and not only with a local origin. In this way the impurity can diffuse from cell to cell and pair with any of the fixed origins. We track the motion of the impurity until it hits one of the origins and average the result over the initial positions of the impurity within one of the cells.
Of course, the above approximation of immobile origins is a rather rough one. We plan to lift this restriction in future work. However, as we will see below, even within such an approximation, the numerical results are in a very good agreement with experimental data and reproduce correctly the temperature profile of the pairing time.
Basically, the algorithm is very simple. Let t, be the average time for the impurity to reach the origin starting from site i within the cell. If a jump from site i to a nearest-neighbor site j takes a time r; [Eq. (2)], the average pairing times t; satisfy the equations t, = gp, , (r, +r, , ), where p; are the probabilities of jumps from i to j, 1/~, Pig j and the summation goes over all neighboring sites to i. The meaning of Eqs. (4) is quite simple: The average time to reach the origin from some site is equal to the time to reach the origin from a neighboring site plus the time to get to this neighboring site averaged over all neighbors. The additional constraint for the linear equations (4) is that the time t, -is equal to zero if i is an origin or one of its nearest neighbors.
The order of the set of linear equations (4) is proportional to the volume of our cell and is of the order of 1/x, where x is the impurity concentration. In our calculations we generated hcp lattices with 100 to more than 1000 lattice sites per cell; i.e. , we were dealing with impurity concentrations between 1% and less than 0.1%%uo. The set of linear equations (4) was solved using the conjugate gradient method. Special attention was paid to the accuracy of the solution in view of the presence of coefficients differing from each other by many orders of magnitude.
We would also like to point out an additional important detail. Occasionally, mostly because of a particular symmetry of a hcp lattice, an impurity can find itself on a lattice site for which one of the nearest neighbors has a very small, or even zero energy mismatch. Such a situation is illustrated in Fig. 1 (1)] and the origin (unlabeled solid circle) is the same for an impurity located at the site 1 or 2 because of obvious symmetry properties of the hcp lattice. The energy mismatch 6E&2 between sites 1 and 2 for the impurity is exactly zero. with the transition from 1 to 2 is exactly zero, 5E, z =0 (if we neglect the interaction with all other origins), while the mismatches with all other neighboring sites are still very large, on the scale of coo. As was mentioned above, the effective tunneling frequency for small mismatches 5E (coo is given by Eq. (3) rather than by Eq. (2). This means that the transitions between sites 1 and 2 take place for small mismatches with a frequency mo, which is much higher than the one for the phonon-assisted jumps given by Eq. (2). This provides a natural cutoff for Eq.
(2), which diverges when 5E~0. In fact, since the cutoff is very large, the results are not very sensitive to its value and correspond formally to the limit of the infinite cutoff. Therefore, the choice of cutoff for numerical calculations does not impose any real restrictions on the value of coo.

IV. RESULTS AND DISCUSSION
We calculated the number n [or the fraction y (t) of pairs; y (0) =0, y(~) =1] of paired impurities as a function of time assuming a random initial distribution of impurities over all sites except the origin and its nearest neighbors. If the pairing process is exponential in time with only one characteristic time~, then the time dependence y ( t) should look like t =s in[1y(t)] (6) and represent a straight line on a logarithmic plot. Figures 2 and 3 show the results for cells with 100 and 512 sites which correspond to impurity concentrations of 1% and 0.19%. As one can see, the real plot is more like a staircase than a straight line [the choice of parameters in Eq.
(2) was so that to ensure the best final fit to the experimental data; see below]. This means that the pairing process is characterized by a set of times rather than a single one. Therefore the precise definition of the pairing time cannot be chosen unambiguously. However, the same difficulty was encountered earlier in experiments where the time dynamics of the numbers of singles and pairs was exponential only in initial stages of pairing. (2) chosen so as to ensure the best fit with the experimental data. ' The fit was done by fixing the position and height of the maximum. As we see, the whole temperature dependence of the characteristic time is described amazingly well. Note that the experimental error bars, which are not given in the figure, are much larger than the distance between the curve and experimental points.
Generally, the well-developed temperature maximum is an inherent feature of the pairing time characterizing the initial stages of the pairing process. The shape of this maximum is rather stable and is not very sensitive to the values of parameters chosen for the calculation, although, of course, the position and height change with the parameters. Equation (2) contains three unknown parameters coo, g, g& of which only two are independent: the bare tunneling frequency coo and the ratio g&/g . Unfortunately, the third parameter still remains free, so that we cannot extract the exact values of the parameters. In Fig. 6 we plotted the values of coo as a function of g. The values of and g are not known theoretically. The value of g should be of the order of the coordination number, i.e. , larger than 12. ' The ratio g&/g determines the position of the maximum. The curves in Fig. 5 correspond to the value of this ratio equal to 2. 91 X 10, which was determined from the best fit to experimental data. However, as is obvious from the figure, the experimental data are not very accurate and are consistent with a position of the maximum anywhere between T=0.4 and 0.7 K. Unfortunately, the position of the maximum is not very sensitive to the value of the ratio g, /g, and such an experimental uncertainty in the position of the maximum can change the value of this ratio by a couple orders of magnitude. On the other hand, this insensitivity is a good sign: It demonstrates that the existence of the maximum is a very general effect which does not depend on some very particular choice of the parameters.
The main parameter of quantum diffusion is the bare tunneling frequency coo. The problem is that it is very difficult (or, in most of the cases, impossible) to extract the value of this parameter from direct measurements. On the other hand, the results of theoretical calculations of coo depend exponentially on some not very well known parameters, which makes such estimates not reliable either. All estimates of coo for known types of quantum diffusion are based on some indirect information and oversimplified assumptions about the diffusion process.    Flax. 6. coo as a function of g.
In this particular case of molecular hydrogen, the only estimate" of co0 provides for it the extremely small value, coo-10 s '. This would correspond (see Fig. 6) to a relatively large value of g-50. The values of g-12 would lead to the value of the tunneling frequency which is larger by two orders of magnitude, coo-10 s '. Probably, these two estimates give upper and lower bounds for co0, whose real value lies somewhere in between. For the sake of comparison, we note that the value of the tunneling frequency for impurities (isotopes) in helium crystals, coo -10 s ', is much higher than even the most optimistic estimate above.
Our results confirm the suggestion of Ref. 3 that the temperature maximum for the pairing time at low temperatures can be explained as a manifestation of a strongly directionally biased diffusion with the bias imposed by the interaction of the pairing particles. The effect is rather simple: One of the diffusion mechanisms, namely, the tunneling accompanied by inelastic phonon scattering, has a strong bias against the direction of increasing mismatches, i.e. , against the mutual approach of the pairing particles. This means that this diffusion mechanism can be, under certain conditions, counterproductive with respect to pairing, creating what was called a "kinetic barrier" in Ref.
3. At high temperatures, when this diffusion mechanism is the fastest one, the main result of a decrease in temperature should be a general slowdown of all motion, including the pairing rate. However, at lower temperatures, the other, temperature-independent, diffusion mechanism with spontaneous emission of phonons becomes equally or even more important. At this point, the freezing of the former mechanism, which is biased against pairing, will lead to a disappearance of the kinetic barrier and to an increase in pairing rates with decreasing temperatures. This effect is quite general and does not depend on details or the presence of other diffusion mechanisms.
Therefore one should expect that more detailed calculations, including the motion of the rest of the particles, will still preserve the well-developed temperature maximum while smoothing out the steps in the staircases in