text:
# High-Order Coupled Cluster Method (CCM) Calculations Via Parallel Processing: An Application To The Kagomé Antiferromagnet

## I Introduction

The coupled cluster method (CCM) refc1 ; refc2 ; refc3 ; refc4 ; refc5 ; refc6 ; refc7 ; refc8 ; refc9  is a highly successful and widely applied technique of modern-day quantum many-body theory. In particular, the CCM has been applied to quantum magnetic systems at zero temperature ccm1 ; ccm2 ; ccm999 ; ccm3 ; ccm4 ; ccm5 ; ccm6 ; ccm7 ; ccm8 ; ccm9 ; ccm10 ; ccm11 ; ccm12 ; ccm13 ; ccm14 ; ccm15 ; ccm16 ; ccm17 ; ccm18 ; ccm19 ; ccm20 ; ccm21  over the last ten years or so. Indeed, the use of computational approaches of the CCM for quantum systems has been found to be particularly efficacious. We note that the subsequent application of these highly accurate computational CCM techniques to other types of quantum system is ripe.

In this article, we present a simple parallelisation procedure for high-order coupled cluster method (CCM) calculations. This approach is investigated and illustrated by the application to the Heisenberg antiferromagnet on the Kagomé lattice with nearest- and next-nearest-neighbour bonds. In particular, we note that we may consider up to approximately 10000 fundamental clusters, although we believe that further increases (of at least another order of magnitude) in the number of fundamental clusters might be possible in future by using such techniques.

## II Formalism

The exact ket and bra ground-state energy eigenvectors, $`|\mathrm{\Psi }`$ and $`\stackrel{~}{\mathrm{\Psi }}|`$, of a general many-body system described by a Hamiltonian $`H`$,

$$H|\mathrm{\Psi }=E_g|\mathrm{\Psi };\stackrel{~}{\mathrm{\Psi }}|H=E_g\stackrel{~}{\mathrm{\Psi }}|,$$
(1)
are parametrised within the single-reference CCM as follows:

$`|\mathrm{\Psi }=\mathrm{e}^S|\mathrm{\Phi }`$ ; $`S={\displaystyle \underset{I0}{}}𝒮_IC_I^+,`$
$`\stackrel{~}{\mathrm{\Psi }}|=\mathrm{\Phi }|\stackrel{~}{S}\mathrm{e}^S`$ ; $`\stackrel{~}{S}=1+{\displaystyle \underset{I0}{}}\stackrel{~}{𝒮}_IC_I^{}.`$ (2)
The single model or reference state $`|\mathrm{\Phi }`$ is required to have the property of being a cyclic vector with respect to two well-defined Abelian subalgebras of multi-configurational creation operators $`\{C_I^+\}`$ and their Hermitian-adjoint destruction counterparts $`\{C_I^{}(C_I^+)^{}\}`$. Thus, $`|\mathrm{\Phi }`$ plays the role of a vacuum state with respect to a suitable set of (mutually commuting) many-body creation operators $`\{C_I^+\}`$. Note that $`C_I^{}|\mathrm{\Phi }=0`$, $`I0`$, and that $`C_0^{}1`$, the identity operator. These operators are furthermore complete in the many-body Hilbert (or Fock) space. Also, the correlation operator $`S`$ is decomposed entirely in terms of these creation operators $`\{C_I^+\}`$, which, when acting on the model state ($`\{C_I^+|\mathrm{\Phi }\}`$), create excitations from it. We note that although the manifest Hermiticity, ($`\stackrel{~}{\mathrm{\Psi }}|^{}=|\mathrm{\Psi }/\mathrm{\Psi }|\mathrm{\Psi }`$), is lost, the normalisation conditions $`\stackrel{~}{\mathrm{\Psi }}|\mathrm{\Psi }=\mathrm{\Phi }|\mathrm{\Psi }=\mathrm{\Phi }|\mathrm{\Phi }1`$ are explicitly imposed. The correlation coefficients $`\{𝒮_I,\stackrel{~}{𝒮}_I\}`$ are regarded as being independent variables, and the full set $`\{𝒮_I,\stackrel{~}{𝒮}_I\}`$ thus provides a complete description of the ground state. For instance, an arbitrary operator $`A`$ will have a ground-state expectation value given as,

$$\overline{A}\stackrel{~}{\mathrm{\Psi }}|A|\mathrm{\Psi }=\mathrm{\Phi }|\stackrel{~}{S}\mathrm{e}^SA\mathrm{e}^S|\mathrm{\Phi }=\overline{A}\left(\{𝒮_I,\stackrel{~}{𝒮}_I\}\right).$$
(3)
We note that the exponentiated form of the ground-state CCM parametrisation of Eq. (2) ensures the correct counting of the independent and excited correlated many-body clusters with respect to $`|\mathrm{\Phi }`$ which are present in the exact ground state $`|\mathrm{\Psi }`$. It also ensures the exact incorporation of the Goldstone linked-cluster theorem, which itself guarantees the size-extensivity of all relevant extensive physical quantities.

The determination of the correlation coefficients $`\{𝒮_I,\stackrel{~}{𝒮}_I\}`$ is achieved by taking appropriate projections onto the ground-state Schrödinger equations of Eq. (1). Equivalently, they may be determined variationally by requiring the ground-state energy expectation functional $`\overline{H}(\{𝒮_I,\stackrel{~}{𝒮}_I\})`$, defined as in Eq. (3), to be stationary with respect to variations in each of the (independent) variables of the full set. We thereby easily derive the following coupled set of equations,

$`\delta \overline{H}/\delta \stackrel{~}{𝒮}_I=0`$ $``$ $`\mathrm{\Phi }|C_I^{}\mathrm{e}^SH\mathrm{e}^S|\mathrm{\Phi }=0,I0;`$ (4)
$`\delta \overline{H}/\delta 𝒮_I=0`$ $``$ $`\mathrm{\Phi }|\stackrel{~}{S}\mathrm{e}^S[H,C_I^+]\mathrm{e}^S|\mathrm{\Phi }=0,I0.`$ (5)
Equation (4) also shows that the ground-state energy at the stationary point has the simple form

$$E_g=E_g(\{𝒮_I\})=\mathrm{\Phi }|\mathrm{e}^SH\mathrm{e}^S|\mathrm{\Phi }.$$
(6)
It is important to realize that this (bi-)variational formulation does not lead to an upper bound for $`E_g`$ when the summations for $`S`$ and $`\stackrel{~}{S}`$ in Eq. (2) are truncated, due to the lack of exact Hermiticity when such approximations are made. However, one can prove that the important Hellmann-Feynman theorem is preserved in all such approximations.

In the case of spin-lattice problems of the type considered here, the operators $`C_I^+`$ become products of spin-raising operators $`s_k^+`$ over a set of sites $`\{k\}`$, with respect to a model state $`|\mathrm{\Phi }`$ in which all spins points “downward” in some suitably chosen local spin axes. The CCM formalism is exact in the limit of inclusion of all possible such multi-spin cluster correlations for $`S`$ and $`\stackrel{~}{S}`$, although in any real application this is usually impossible to achieve. It is therefore necessary to utilise various approximation schemes within $`S`$ and $`\stackrel{~}{S}`$. The three most commonly employed schemes previously utilised have been: (1) the SUB$`n`$ scheme, in which all correlations involving only $`n`$ or fewer spins are retained, but no further restriction is made concerning their spatial separation on the lattice; (2) the SUB$`n`$-$`m`$ sub-approximation, in which all SUB$`n`$ correlations spanning a range of no more than $`m`$ adjacent lattice sites are retained; and (3) the localised LSUB$`m`$ scheme, in which all multi-spin correlations over all distinct locales on the lattice defined by $`m`$ or fewer contiguous sites are retained.

## III A “Brute-Force” Parallelisation Procedure

A simple “brute-force” approach in order to solve the CCM equations in parallel is firstly to rearrange Eqs. (4) and (5) for the ket and bra states, respectively, where

$`𝒮_I`$ $`=`$ $`f_I(𝒮_1,\mathrm{},𝒮_{I1},𝒮_{I+1},\mathrm{},𝒮_{N_f});`$ (7)
$`\stackrel{~}{𝒮}_I`$ $`=`$ $`\stackrel{~}{f}_I(𝒮_1,\mathrm{},𝒮_{N_f};\stackrel{~}{𝒮}_1,\mathrm{},\stackrel{~}{𝒮}_{I1},\stackrel{~}{𝒮}_{I+1},\mathrm{},\stackrel{~}{𝒮}_{N_f});`$ (8)
and where $`N_f`$ is the number of fundamental configurations. Note that there are therefore $`N_f`$ equations for both $`f_I`$ and $`\stackrel{~}{f}_I`$, which refer to the ket and bra states, respectively. We note that each equation contains a finite-number of terms and that, indeed, this is always the case for the CCM for Hamiltonians that contain a finite number of creation and destruction operators.

The “brute-force” approach is to iterate the set of equations in Eq. (7) to convergence first in order to obtain the ket-state correlation coefficients. We then iterate the set of equations Eq. (8) to convergence. This approach appears to function adequately as long as one is in the region where the model state is a good starting point. Clearly, this approach will become more difficult as we approach any critical points in the CCM equations. A simple parallelisation technique is now to split the problem of determining and solving the $`N_f`$ equations for both the bra and ket states of Eqs. (7) and (8) between each processor in the parallel cluster. Thus, each equation for $`I=\{1,\mathrm{},N_f\}`$ is ever dealt with by one node in the cluster only and the results of each node are collected together at each iteration in the “brute-force” algorithm. (Again, we remember that we iterate the ket-state equations first and then the bra-state equations.) Thus, the set of equations for the ket and bra states are solved independently and in parallel.

This approach has a number of advantages. Firstly, we do not need to define a Jacobian (e.g., as for the Newton-Raphson technique) and so we save on RAM by using this procedure. Secondly, the task of evaluating and saving the CCM equations is reduced by a factor of the number of machines used. Indeed, in our implementation we saved all of the data for a particular set of equations to disk locally on each machine, although we note that one could even re-derive the CCM equations as and when necessary in order to further reduce disk usage. Finally, this a simple approach and so is easy to amend existing code and so be parallelised.

This very simple approach was found to be surprisingly successful, although we note that the number of iterations for such a direct iteration method is at least an order of magnitude greater than that for the Newton-Raphson technique, for example. Indeed, we were able to reproduce those results determined on a single machine up to the LSUB7 level of approximation by using our parallelised implementation. This was found to be an excellent test of the validity of our results. However, we were able to solve the LSUB8 approximation for the ground-state energy only by using parallel processing. We note that this approximation contained 10707 fundamental configurations and that this is twice as big as the previous “largest” CCM calculation for spin systems. We believe that further increases in the number of fundamental configurations by at least another order of magnitude might be possible even by using moderate levels of parallelisation (e.g., 10 or 20 nodes) and for our “brute-force” approach.

Finally, we would wish to implement a more sophisticated solution of this problem using parallel processing in the future. A possible way of doing this would be to use PLAPACK, for example, in order to implement the Newton-Raphson technique “in parallel.” Other strategies might include numerical approximations to the Jacobian for the Newton-Raphson technique or other less memory-intensive ways of determining and solving the CCM equations of Eqs. (4) and (5).

## IV The $`J_1`$$`J_2`$ Kagomé Model

The spin-half $`J_1`$$`J_2`$ Kagomé Model model is is given by

$$H=J_1\underset{i,j}{}𝐬_i𝐬_jJ_2\underset{\{i,k\}}{}𝐬_i𝐬_k,$$
(9)
where $`i,j`$ runs over all nearest-neighbour (n.n.) bonds and $`\{i,k\}`$ runs over all next-nearest-neighbour (n.n.n.) bonds on the Kagomé lattice. Note that each bond is counted once and once only, and that we now set $`J_1=1.0`$.

We choose a model state $`|\mathrm{\Phi }`$ in which the lattice is divided into three sublattices, denoted $`\{`$A,B,C$`\}`$. The spins on sublattice A are oriented along the negative z-axis, and spins on sublattices B and C are oriented at $`+120^{}`$ and $`120^{}`$, respectively, with respect to the spins on sublattice A. Our local axes are chosen by rotating about the $`y`$-axis the spin axes on sublattices B and C by $`120^{}`$ and $`+120^{}`$ respectively, and by leaving the spin axes on sublattice A unchanged. The n.n.n. bonds for $`J_2>0`$ stabilise this type of Néel ordering in the quantum ground state.

Careful CCM and finite-sized calculations ccm18 ; kaf1 ; kaf2 ; kaf3  have been performed for the spin-half Kagomé antiferromagnet (KAF) ($`J_2=0`$), and these results indicate that the classical Néel-like ordering observed, for example, in the triangular-lattice antiferromagnet (TAF) is not seen for the quantum KAF model. The best estimate for the ground-state energy of the KAF via finite-sized calculations kaf3  stands at $`E_g/N`$=$``$0.43. Furthermore, CCM and series expansion results ccm18 ; kaf2  indicate that the ground-state of the KAF is disordered. Indeed, a variational calculation kaf4  which utilised a dimerised basis also found that the ground state of the KAF is some sort of spin liquid. Thus, it is highly likely that a quantum phase transition will occur at or near to $`J_2=0`$.

Our results for the ground-state energy per spin are given in Fig. 1 using the LSUB$`m`$ approximation with $`m=\{2,3,4,5,6,7\}`$. We see that our results converge rapidly with $`m`$ across a wide range of $`J_2`$ up to an including the point $`J_2=0`$. Results for the special case $`J_2=0`$ are given in Table 1. (Note that these results are slightly different than previous CCM calculations in Ref. ccm18  due to a different underlying “interpolating” lattice in this case.) We are able to consider 10707 fundamental configurations using our parallel approach and a simple extrapolation in the limit $`m\mathrm{}`$ gives $`E_g/N=0.43104`$, which should be considered to previous results of the best of other calculations ccm18 ; kaf3 , namely, $`E_g/N0.43`$.

The sublattice magnetisation is defined by,

$$M=\frac{2}{N}\underset{k=1}{\overset{N}{}}s_k^z.$$
(10)
where $`k`$ runs over all lattice sites on the Kagomé lattice. Results for the sublattice magnetisation are given in Fig. 2 using the LSUB$`m`$ approximation with $`m=\{2,3,4,5,6,7\}`$. Results for the special case $`J_2=0`$ are given in Table 1. Note that CCM equations break down in the region $`J<0`$ at their critical points. In the limit $`m\mathrm{}`$ we expect that this would occur at $`J_2=0`$. The “upturn” in $`M`$ for LSUB6 that occurs in the region $`J_2<0`$ is an artifact and such effects occur when the model state becomes an increasingly bad starting point. The CCM equations for LSUB7 break down at a critical point near to $`J_2=0`$, and we note that our results are consistent with the assertion that this model does not demonstrate Néel ordering at $`J_2=0`$.

## V Conclusion

We have shown that it is possible to parallelise the high-order CCM problem implemented using computational techniques. A simple “brute-force” method was used, although more sophisticated approaches were briefly considered.

An application to the spin-half $`J_1`$$`J_2`$ model on the Kagomé lattice was considered. Results for the ground-state energy and sublattice magnetisation (with respect to three-sublattice Néel ordering, as for the triangular lattice antiferromagnet) was considered. It was seen that our results were highly converged over a large range of $`J_2`$. Our results were consistent with the assertion that this model demonstrates no Néel ordering at $`J_2=0`$. All results were determined using the parallelised approach and with a single processor, except for a single result for the ground-state energy at $`J_2=0`$ which contained 10707 fundamental configurations and which was only able to be carried out “in parallel.” Our best estimate for the ground-state energy per spin of the spin-half Kagomé lattice antiferromagnet with nearest-neighbour terms only ($`J_2=0`$) was $`E_g/N=0.43104`$.

The increasing accuracy of the high-order CCM technique and its ability to be applied to even strongly frustrated systems (such as the model considered here) are strong advantages. However, there are many possible future extensions of high-order CCM which would make a parallel implementation even more useful. The CCM has been used previously ccm10  to determine the excitation energies of the XXZ model on one-, two-, and three-dimensional lattices with great accuracy using “localised” approximation schemes. We note that it is, in principle, straightforward to extend these treatments in order to determine the excitation spectra of even highly frustrated spin systems with equal accuracy and for lattices of complex crystallographic symmetries. Indeed, the use of existing codes would make this process simpler.

Such determination of the the excitation spectra of quantum spin systems using high-order CCM techniques might be made even more useful by the utilisation of “non-Néel” model states. The model considered here uses just such a “Néel” model state and typically we perform some kind of notational rotation of the spin axis such that all of the spin “point downwards”. This is then the starting point for our calculations. We note that non-Néel model states (such as the dimer solid state) have already been used at low orders of approximation to consider spin systems which demonstrate novel ordering. It would again be a straightforward process to develop our high-order CCM approach for the “localised” approximation schemes in order to consider such non-Néel model states.

Another extension of the current calculations is to bosonic or fermionic systems. Indeed, the formalism has already been developed for lattice bosonic systems, which we note have much in common (from a “CCM notational” point-of-view) which large quantum spin number, $`s`$, spin systems previously considered using high-order CCM techniques ccm20 .

The above calculations were presented for the infinite-lattice case, $`N\mathrm{}`$. We note that it is also possible to consider the case where $`N`$ is finite, and, again, high-order techniques would be invaluable here.

The application of the CCM to quantum spin systems at non-zero temperature might also be accomplished in future. Indeed, the application of the CCM at both zero- and non-zero temperatures might help to explain the subtle interplay of quantum and thermal fluctuations in driving phase transitions over a wide range of physical parameters.


format: text

text:
# Enhanced production of HD and D2 molecules on small dust grains in diffuse clouds

## 1 Introduction

The formation of molecular hydrogen on the surfaces of dust grains in the interstellar medium is a process of fundamental importance because the formation rate in the gas phase cannot account for the observed abundance (Gould & Salpeter, 1963). Once produced on grain surfaces and ejected or desorbed into the gas, hydrogen molecules are a necessary component in the initiation of gas-phase chemical reaction networks that give rise to the chemical complexity observed in interstellar clouds. The formation of molecular hydrogen has been studied both theoretically (Hollenbach & Salpeter, 1970, 1971; Hollenbach et al., 1971; Williams, 1968; Smoluchowski, 1981, 1983; , 1985; Duley & Williams, 1986; Pirronello & Avera, 1988; Sandford & Allamandola, 1993; Takahashi et al., 1999) and experimentally (Brackmann, 1961; Schutte et al., 1976; Pirronello et al., 1997a, b, 1999; Manico et al., 2001). Although formation of H2 is clearly the most studied surface reaction pertaining to interstellar clouds, there are many other reactions that are thought to take place, for the most part involving H atoms and other species. Reactions are typically considered to occur via the Langmuir-Hinshelwood mechanism, in which weakly-bound surface species diffuse towards one another. Some authors have considered different mechanisms, where chemisorption plays a role (Cazaux & Tielens, 2002).

The computational modeling of diffusive chemical reactions on dust grains in the interstellar medium is typically done using rate equation models (Pickles & Williams, 1977; d’Hendecourt et al., 1985; Brown, 1990; Brown & Charnley, 1990; Hasegawa et al., 1992; Hasegawa & Herbst, 1993a, b; Caselli et al., 1993; Willacy & Williams, 1993; Shalabiea & Greenberg, 1994). As long as a sizable number of reactive species exists on a grain surface, the rate equation results are valid. However, in the limit of small average surface populations ($`1`$) of reactive species, rate equations are not always suitable, because they take into account only average densities and ignore the fluctuations as well as the discrete nature of the populations of the atomic and molecular species (Tielens, 1995; Tielens & Hagen, 1982; Charnley et al., 1997; Caselli et al., 1998; Shalabiea et al., 1998; Stantcheva et al., 2001, 2002). Depending on the parameters chosen, this limit can pertain mainly to small grains or can even to “classical” grains, namely grains of 0.1$`\mu `$ radius. For example, as the number of H atoms, $`N_\mathrm{H}`$, on a grain fluctuates in the range of 0, 1 or 2, the H2 formation rate cannot be obtained from the average number alone. This can be easily understood, since recombination requires at least two H atoms simultaneously on the surface.

Recently, a master equation approach was proposed that is suitable for the simulation of diffusive chemical reactions on interstellar grains in the accretion limit (, 2001, 2001; Biham & Lipshtat, 2002). The master equation takes into account both the discrete nature of the reactive species as well as the fluctuations in their populations. (2001) considered a simple system in which H and O atoms land on classical grains to produce H2, OH, and O2. Biham and co-workers focused on the formation of H2. With a slow diffusion rate for H based on the analysis by Katz et al. (1999) of experimental work by Pirronello, Vidali, & co-workers (Pirronello et al., 1997a, b), they found the master equation approach to be needed mainly for small grains. For the formation of H2, the dynamical variables are the probabilities $`P_\mathrm{H}(N_\mathrm{H})`$ that there are $`N_\mathrm{H}`$ atoms on the grain at time $`t`$. The time derivatives $`\dot{P}_\mathrm{H}(N_\mathrm{H})`$, $`N_\mathrm{H}=0,1,2,\mathrm{},`$ are expressed in terms of the adsorption, reaction and desorption terms. The master equation provides the time evolution of $`P_\mathrm{H}(N_\mathrm{H})`$, $`N_\mathrm{H}=0,1,2,\mathrm{}`$, from which the recombination rate can be calculated. The master equation has also been used for the study of reaction networks involving multiple species (Stantcheva et al., 2002; Stantcheva & Herbst, 2003).

In addition to the master equation method, a Monte Carlo method (Tielens & Hagen, 1982) was proposed (Charnley, 2001), and applied to surface reaction networks involving multiple species (Caselli et al., 2002; Stantcheva & Herbst, 2003). The weakness of the Monte Carlo approach to surface chemistry is that it cannot be used easily for models in which gas-phase reactions take place simultaneously. Finally, a semi-empirical approach, known as the modified rate method, was proposed by Caselli et al. (1998). Although artificial, the modified rate method is easy to employ and has been used with mixed success (Stantcheva et al., 2001; Caselli et al., 2002).

In addition to H, which is the most abundant atom in the universe, there is a density of D atoms, which is smaller by a factor of about $`10^5`$ compared with the density of H atoms in regions where the gas is primarily atomic. If deuterated hydrogen molecules, namely HD and D2, are formed on grain surfaces in the same way as H2, these molecules are expected to be rare due to the low abundance of D atoms. Of course, some differences between the formation rates of HD and D2 and that of H2 are expected because of the possibly different diffusion and desorption rates of H and D atoms on the surface. It was realized early in the study of interstellar chemistry, however, that the formation of HD in diffuse interstellar clouds occurs more efficiently in the gas phase via the sequence of reactions

$`\mathrm{H}^+`$ $`+`$ $`\mathrm{D}\mathrm{D}^++\mathrm{H}`$
$`\mathrm{D}^+`$ $`+`$ $`\mathrm{H}_2\mathrm{HD}+\mathrm{H}^+,`$ (1)
which can lead to enhancements of $`10^2`$ over surface formation (Watson, 1973; Dalgarno et al., 1973) and explain the observed HD/H2 abundance ratio in diffuse clouds. By the stage in which dense interstellar clouds have been formed, HD has become the major repository of deuterium, and is the starting point for a rich chemistry of deuterium fractionation that occurs in the gas phase via ion-molecule reactions (Watson, 1973; Millar et al., 1989; Charnley et al., 1997; Roberts & Millar, 2000a, b; Rodgers & Charnley, 2001; Roberts et al., 2002). Since gas-phase fractionation produces a large D/H atomic abundance ratio in dense clouds, it then becomes possible that surface processes can produce multiply deuterated species in high abundance following the adsorption of D atoms (Charnley et al., 1997; Rodgers & Charnley, 2001; Caselli et al., 2002; Stantcheva & Herbst, 2003).

Although it has been thought for some time that the formation of HD in diffuse clouds occurs primarily in the gas, the surface processes for the formation of HD and D2 have not been looked at carefully. In this paper we propose a mechanism taking place on the surfaces of dust grains that, under certain conditions, can provide a significantly enhanced production of deuterated molecules. We demonstrate the mechanism for the simple case that involves only H and D atoms in diffuse clouds, and examine the formation rates of HD and D2 vs. that of H2. The mechanism is based on the assumption that D atoms stick to the surface more strongly than H atoms so that their desorption rate is lower. An isotope effect of this type has been observed in various experimental situations (Koehler et al., 1988; Hoogers et al., 1995).

Due to the very low cosmic abundance ratio of D vs. H ($`10^5`$) most of the atoms that hit the grain surfaces in diffuse clouds are H atoms. However, the residence time of D atoms on grains is expected to be longer than that of H atoms, due to their lower desorption rate. Thus, the D/H atomic abundance ratio on the grains is enhanced compared with the gas-phase ratio. As a result, newly adsorbed H (or D) atoms are more likely to find D atoms already residing on the grains, giving rise to an enhanced production of HD molecules. This enhancement turns out to be even larger on very small grains on which only a handful of H and D atoms typically reside. This effect is expected to be significant in diffuse clouds, where the density of the gas phase is low and the grain size distribution is broad and dominated by very small grains (Li & Draine, 2001).

The paper is organized as follows. The interaction between hydrogen atoms and dust grains in interstellar clouds in described in Sec. 2. The formation of molecular hydrogen on large grains is described in Sec. 3 using the rate equation model. The formation of molecules on small grains is presented in Sec. 4 using the master equation approach, followed by an asymptotic analysis for very small grains in Sec. 5. The simulations and results are presented in Sec. 6. The results are discussed in Sec. 7 and summarised in Sec. 8. In Appendix A we consider the relevance of the isotope effect studied here to related experiments.

## 2 The interaction between hydrogen atoms and dust grains

Consider a diffuse interstellar cloud dominated by a density $`n_\mathrm{H}`$ (cm-3) of H atoms and also including a density $`n_\mathrm{D}=fn_\mathrm{H}`$ of D atoms, where $`f`$ is the D$`/`$H abundance ratio (the cosmic ratio is approximately $`f=10^5`$). The typical velocity of the hydrogen isotope I (H or D) in the gas phase is denoted by $`v_\mathrm{I}`$ (cm s-1). It is given by

$$v_\mathrm{I}=\sqrt{\frac{8}{\pi }\frac{k_BT_{\mathrm{gas}}}{m_\mathrm{I}}},$$
(2)
where the mass $`m_\mathrm{I}`$ is $`m_\mathrm{H}=1.67\times 10^{24}`$ (grams) for H atoms and $`m_\mathrm{D}=3.34\times 10^{24}`$ (grams) for D atoms and the temperature of the gas is denoted by $`T_{\mathrm{gas}}`$. The cloud includes a density $`n_{\mathrm{gr}}`$ (cm-3) of dust grains. To evaluate the flux of atoms onto grain surfaces we will assume, for simplicity, that the grains are spherical with a radius $`r`$ (cm). The cross section of a grain is $`\sigma =\pi r^2`$ and its surface area is $`4\sigma `$. The number of adsorption sites on a grain is denoted by $`S`$, and their density $`s`$ (sites cm-2) on the surface is given by $`s=S/(4\sigma )`$. The fluxes $`F_\mathrm{H}`$ (atoms s-1) and $`F_\mathrm{D}`$ of H and D atoms onto the surface of a single grain are given by $`F_\mathrm{I}=n_\mathrm{I}v_\mathrm{I}\sigma `$. The atoms stick to the surface and hop as random walkers between adjacent sites until they either desorb as atoms or recombine into molecules. The desorption rates of H and D atoms on the surface are given by

$$W_\mathrm{I}=\nu \mathrm{exp}[E_1(\mathrm{I})/k_BT],$$
(3)
where $`\nu `$ is the attempt rate (typically assumed to be $`10^{12}`$ s-1), $`E_1(\mathrm{I})`$ is the energy needed for desorption of an atom of isotope I and $`T`$ is the surface temperature. The hopping rates of the atoms from site to site are

$$a_\mathrm{I}=\nu \mathrm{exp}[E_0(\mathrm{I})/k_BT],$$
(4)
where $`E_0(\mathrm{I})`$ is the activation energy barrier for diffusion of the isotope I. Here we assume that diffusion occurs only by thermal hopping, in agreement with recent experimental results (Katz et al., 1999). In the faster diffusion rates often adopted by astrochemists, tunneling is allowed to occur (Hasegawa et al., 1992) and results in very different rates for H and D. The rate $`A_\mathrm{I}=a_\mathrm{I}/S`$ is approximately the inverse of the time $`t_s(\mathrm{I})`$ required for an atom of isotope I to visit nearly all the adsorption sites on the grain surface, because in two dimensions the number of distinct sites visited by a random walker is linearly proportional to the number of steps, up to a logarithmic correction (Montroll & Weiss, 1965).

In this paper, we take the density of adsorption sites on the surface as $`s=5\times 10^{13}`$ (sites cm-2) (, 2001). The density of H atoms is $`n_\mathrm{H}=50`$ (cm-3) and the density of D atoms is $`n_\mathrm{D}=fn_\mathrm{H}`$, where $`f=10^5`$ is the cosmic abundance ratio, which according to Eq. (2) gives rise to to the ratio $`F_\mathrm{D}/F_\mathrm{H}=10^5/\sqrt{2}`$. The gas temperature is taken as $`T_{\mathrm{gas}}=90`$ K, thus $`v_H=1.37\times 10^5`$ and $`v_D=9.7\times 10^4`$ (cm s-1). The surface temperature is $`T=18`$ K, although variations from this value are considered. These values are within the range of typical parameters for diffuse interstellar clouds.

The parameters used for the interaction between the adsorbed H and D atoms and the grain are crucial to the determination of the isotope effect. Since D is heavier than H, it sits lower than H in any potential surface when zero-point energy is included. The smaller zero-point energy of D leads to larger values for the desorption energy and the barrier against diffusion, although the differences between H and D are uncertain. In the analysis of their experimental results, for example, Katz et al. (1999) assumed no difference between H and D atoms at all. In a previous paper, Caselli et al. (2002) argued for a difference of 2 meV between both the desorption energies and barriers against diffusion. Here, the adopted energy barriers for diffusion are $`E_0(H)=35`$ meV for H atoms and $`E_0(D)=37`$ meV for D atoms, while the desorption energies are $`E_1(H)=50`$ meV for H atoms and $`E_1(D)=60`$ meV for D atoms. The values for H atoms are intermediate between those obtained experimentally on olivine and amorphous carbon surfaces (Katz et al., 1999). The difference between the desorption energies is assumed to be 10 meV; an argument in favour of this larger value is that the barrier against diffusion balances the zero-point energy of a potential well with that of a saddle point while desorption does not possess a saddle point. In any case, the difference in desorption energies between H and D is crucial to the isotope effect discussed here.

## 3 Molecular hydrogen formation on large grains

After being adsorbed onto the surface, the H and D atoms hop between adsorption sites until they either desorb or meet one another and form molecules. The number of atoms of isotope I on the surface of a grain is denoted by $`N_\mathrm{I}`$. The rate equation analysis deals with the expectation values $`N_\mathrm{I}`$, I $`=`$ H, D obtained by averaging over a large ensemble of identical grains under given physical conditions. The rate equations for the expectation values of the H and D populations take the form

$`{\displaystyle \frac{dN_\mathrm{H}}{dt}}`$ $`=`$ $`F_\mathrm{H}W_\mathrm{H}N_\mathrm{H}2A_\mathrm{H}N_\mathrm{H}^2(A_\mathrm{H}+A_\mathrm{D})N_\mathrm{H}N_\mathrm{D}`$
$`{\displaystyle \frac{dN_\mathrm{D}}{dt}}`$ $`=`$ $`F_\mathrm{D}W_\mathrm{D}N_\mathrm{D}2A_\mathrm{D}N_\mathrm{D}^2(A_\mathrm{H}+A_\mathrm{D})N_\mathrm{H}N_\mathrm{D}.`$ (5)
These equations provide the time derivatives of the sizes of the populations of reactive species. Each equation includes a flux term, a desorption term and two reaction terms. The first reaction term accounts for reactions between atoms of the same specie, while the second reaction term accounts for the reaction that produces HD molecules. The rate of the latter reaction is proportional to the rate at which H and D atoms encounter each other, which is proportional to the sum of the hopping rates of the two species. Therefore, the pre-factor of this term is the sum of $`A_\mathrm{H}`$ and $`A_\mathrm{D}`$. By integration of the rate equations for a given initial condition one obtains the time dependence of the population sizes. The rate equations provide a good description of the surface reactions with our chosen parameters as long as the grain is not too small and the flux is not too low. The rate equations are strictly valid when both the D and H populations per grain are greater than unity. However, we have found that, despite a very small D population in even our largest grains, as long as the H population per grain is not too much smaller than unity, the rate equations provide a reasonable approximation to the master equation results.

We will now solve the rate equations for the case of steady state conditions, in which $`dN_\mathrm{I}/dt=0`$ for I $`=`$ H, D. Due to the small ratio between the fluxes of $`D`$ and $`H`$ atoms, the effect of the D atoms on the density $`N_\mathrm{H}`$ can be neglected. Thus, neglecting the HD formation term in the first equation in (5) we obtain

$$N_\mathrm{H}=\frac{1}{4}\left(\frac{W_\mathrm{H}}{A_\mathrm{H}}\right)\left[1+\sqrt{1+8\frac{(F_\mathrm{H}/W_\mathrm{H})}{(W_\mathrm{H}/A_\mathrm{H})}}\right].$$
(6)
Inserting this solution into the second equation in (5) we obtain

$$N_\mathrm{D}=\frac{1}{4}\left(\frac{W_{\mathrm{eff}}}{A_\mathrm{D}}\right)\left[1+\sqrt{1+8\frac{(F_\mathrm{D}/W_{\mathrm{eff}})}{(W_{\mathrm{eff}}/A_\mathrm{D})}}\right],$$
(7)
where

$$W_{\mathrm{eff}}=W_\mathrm{D}+(A_\mathrm{H}+A_\mathrm{D})N_\mathrm{H}$$
(8)
is the effective desorption coefficient for the D atoms. The ratio $`N_\mathrm{D}/N_\mathrm{H}`$ between the population sizes of D and H on the grain surface is obtained from Eqs. (6) and (7). This ratio turns out to be higher than the gas phase ratio $`D/H`$ due to the lower desorption rate of D atoms. The production rate of HD molecules per grain is

$$R_{\mathrm{HD}}=(A_\mathrm{H}+A_\mathrm{D})N_\mathrm{H}N_\mathrm{D}$$
(9)
and of $`H_2`$ molecules is

$$R_{\mathrm{H}_2}=A_\mathrm{H}N_\mathrm{H}^2.$$
(10)
The ratio between these production rates is thus

$$\frac{R_{\mathrm{HD}}}{R_{\mathrm{H}_2}}=\left(1+\frac{A_\mathrm{D}}{A_\mathrm{H}}\right)\frac{N_\mathrm{D}}{N_\mathrm{H}}.$$
(11)
Due to the assumption that $`E_0(H)<E_0(D)`$, namely that the diffusion of H atoms is faster than of D atoms, we obtain that the additional enhancement factor multiplying $`N_\mathrm{D}/N_\mathrm{H}`$ is in the range $`1<1+A_\mathrm{D}/A_\mathrm{H}<2`$. The production rate of D2 molecules is

$$R_{\mathrm{D}_2}=A_\mathrm{D}N_\mathrm{D}^2,$$
(12)
and the ratio between the D2 and H2 production rates is

$$\frac{R_{\mathrm{D}_2}}{R_{\mathrm{H}_2}}=\frac{A_\mathrm{D}}{A_\mathrm{H}}\frac{N_\mathrm{D}^2}{N_\mathrm{H}^2}.$$
(13)
## 4 Molecular hydrogen formation on small grains

The master equation that describes the hydrogen-deuterium system consists of a two-dimensional matrix of equations. These equations describe the time derivatives of the probabilities $`P(N_\mathrm{H},N_\mathrm{D})`$ of having $`N_\mathrm{H}`$ hydrogen atoms and $`N_\mathrm{D}`$ deuterium atoms simultaneously on the surface. The master equation takes the form

$`\dot{P}(N_\mathrm{H},N_\mathrm{D})`$ $`=`$ $`F_\mathrm{H}[P(N_\mathrm{H}1,N_\mathrm{D})P(N_\mathrm{H},N_\mathrm{D})]`$
$`+`$ $`F_\mathrm{D}[P(N_\mathrm{H},N_\mathrm{D}1)P(N_\mathrm{H},N_\mathrm{D})]`$
$`+`$ $`W_\mathrm{H}[(N_\mathrm{H}+1)P(N_\mathrm{H}+1,N_\mathrm{D})N_\mathrm{H}P(N_\mathrm{H},N_\mathrm{D})]`$
$`+`$ $`W_\mathrm{D}[(N_\mathrm{D}+1)P(N_\mathrm{H},N_\mathrm{D}+1)N_\mathrm{D}P(N_\mathrm{H},N_\mathrm{D})]`$
$`+`$ $`A_\mathrm{H}[(N_\mathrm{H}+2)(N_\mathrm{H}+1)P(N_\mathrm{H}+2,N_\mathrm{D})N_\mathrm{H}(N_\mathrm{H}1)P(N_\mathrm{H},N_\mathrm{O})]`$
$`+`$ $`A_\mathrm{D}[(N_\mathrm{D}+2)(N_\mathrm{D}+1)P(N_\mathrm{H},N_\mathrm{D}+2)N_\mathrm{D}(N_\mathrm{D}1)P(N_\mathrm{H},N_\mathrm{D})]`$
$`+`$ $`(A_\mathrm{H}+A_\mathrm{D})[(N_\mathrm{H}+1)(N_\mathrm{D}+1)P(N_\mathrm{H}+1,N_\mathrm{D}+1)N_\mathrm{H}N_\mathrm{D}P(N_\mathrm{H},N_\mathrm{D})],`$
where $`N_\mathrm{H},N_\mathrm{D}=0,1,2,\mathrm{}`$. Each of these equations includes three sets of terms. The first set describes the effect of the incoming fluxes $`F_H`$ and $`F_D`$ on the probabilities. The probability $`P(N_\mathrm{H},N_\mathrm{D})`$ increases when an H atom sticks to the surface of a grain that already has $`N_\mathrm{H}1`$ adsorbed H atoms and $`N_\mathrm{D}`$ adsorbed D atoms. There is a similar term for the adsorption of a D atom. The probability $`P(N_\mathrm{H},N_\mathrm{D})`$ decreases when a new atom is adsorbed on a grain with $`N_\mathrm{H}`$ H atoms and $`N_\mathrm{D}`$ D atoms that are already adsorbed. The second set of terms describes the effect of desorption. An H \[D\] atom that is desorbed from a grain on which there are $`N_\mathrm{H}`$ adsorbed H atoms and $`N_\mathrm{D}`$ adsorbed D atoms decreases the probability $`P(N_\mathrm{H},N_\mathrm{D})`$ and increases the probability $`P(N_\mathrm{H}1,N_\mathrm{D})`$ \[$`P(N_\mathrm{H},N_\mathrm{D}1)`$\]. The third set includes three terms describing the reactions that lead to the formation of H2, D2 and HD, respectively. The rate of each of these reactions is proportional to the sum of the sweeping rates of the atoms involved ($`A_\mathrm{H}`$ or $`A_\mathrm{D}`$) times the number of pairs of atoms that may participate in the reaction. In the equations in which $`N_\mathrm{H}=0,1`$ or $`N_\mathrm{D}=0,1`$, some of the terms vanish. These terms can be easily identified since their initial or final states include a negative number of H or D atoms on the grain.

By suitable summation over the probabilities, the master equation provides the moments of the population size distribution, given by

$$N_\mathrm{H}^kN_\mathrm{D}^{\mathrm{}}=\underset{N_\mathrm{H}=0}{\overset{\mathrm{}}{}}\underset{N_\mathrm{D}=0}{\overset{\mathrm{}}{}}N_\mathrm{H}^kN_\mathrm{D}^{\mathrm{}}P(N_\mathrm{H},N_\mathrm{D}),$$
(15)
where $`k,\mathrm{}=0,1,2,\mathrm{}`$ \[note that in case that $`k=0`$ ($`\mathrm{}=0`$) the corresponding power, $`N_\mathrm{H}^k`$ ($`N_\mathrm{D}^{\mathrm{}}`$), is removed from the equation\]. The order of each moment is specified by the sum of these powers, namely by $`m=k+\mathrm{}`$. Clearly, there are two first-order moments, $`N_\mathrm{H}`$ and $`N_\mathrm{D}`$, three second-order moments, $`N_\mathrm{H}^2`$, $`N_\mathrm{D}^2`$ and $`N_\mathrm{H}N_\mathrm{D}`$, as well as $`m+1`$ moments of any higher order $`m>2`$. The average rates of production of the H2, HD and D2 molecules on the surface of a single grain can be obtained by a proper summation over the third set of terms in the master equation (4). These production rates can be expressed in terms of the first and second moments of the population size distribution (, 2001; Stantcheva et al., 2001, 2002):

$`R_{\mathrm{H}_2}`$ $`=`$ $`A_\mathrm{H}\left(N_\mathrm{H}^2N_\mathrm{H}\right)`$
$`R_{\mathrm{HD}}`$ $`=`$ $`\left(A_\mathrm{H}+A_\mathrm{D}\right)N_\mathrm{H}N_\mathrm{D}`$
$`R_{\mathrm{D}_2}`$ $`=`$ $`A_\mathrm{D}\left(N_\mathrm{D}^2N_\mathrm{D}\right).`$ (16)
The fact that the production rates depend on the second moments explains the failure of the rate equations in the limit in which the average surface populations of reactive species become $`1`$. In the rate equations, which account only for the first moments, second moments are approximated by products of the corresponding first moments. The equations fail to account for the production rates when this approximation is no longer valid.

The production rates given by Eqs. (16) can be obtained by direct integration of the master equation (4) as a set of time dependent coupled differential equations, using a standard Runge-Kutta stepper. In this paper we focus on the production rates under steady state conditions. Therefore, the data collection is done only after the populations of the reactive species have reached their asymptotic sizes under the simulated conditions. In order to limit the number of equations, suitable cutoffs are imposed, either on the population size of each reactive specie or on the entire population of adsorbed atoms on the grain surface (Stantcheva et al., 2001). The production rates can also be obtained from Monte Carlo simulations, based on the master equation (Charnley, 2001).

## 5 The limit of very small grains

In the limit of very small grains, one can reach the situation in which both $`N_\mathrm{H}1`$ and $`N_\mathrm{D}1`$. In this limit the formation rate of molecules is reduced and the dominant mechanism for removing atoms from the surface is thermal desorption. The formation of molecules can then be considered as a perturbation on top of the adsorption and desorption processes. This limit can be described by the master equation using cutoffs that do not allow more than a single atom on the grain. The master equation then consists of only three coupled equations, for $`\dot{P}(0,0)`$, $`\dot{P}(1,0)`$ and $`\dot{P}(0,1)`$. With these cutoffs, the average population sizes are related to the probabilities by

$`N_\mathrm{H}`$ $`=`$ $`P(1,0)`$
$`N_\mathrm{D}`$ $`=`$ $`P(0,1).`$ (17)
Solving the three coupled equations for steady state conditions we obtain that

$`N_\mathrm{H}`$ $`=`$ $`F_\mathrm{H}/W_\mathrm{H}`$
$`N_\mathrm{D}`$ $`=`$ $`F_\mathrm{D}/W_\mathrm{D}.`$ (18)
Consider the situation in which there is one atom on the surface of a grain. Since the grain is very small an atom can easily scan the whole surface before it desorbs, namely $`A_\mathrm{I}W_\mathrm{I}`$. Thus, when a second atom is adsorbed, it is most likely to form a molecule with the atom that is already on the surface. Therefore, the production rate does not depend on the diffusion rates, a point noted by Tielens & Hagen (1982) and subsequently used to define the term “accretion” limit rigorously, as distinguished from the “reaction” limit. In the former limit, the production rates are given by the flux terms of H and D atoms in the master equation, multiplied by the probabilities of having a single atom (either H or D) on the surface, namely, $`F_\mathrm{I}P(0,1)`$ and $`F_\mathrm{I}P(1,0)`$, where I $`=`$ H, D. Using Eqs. (17), the production rates of molecules take the form

$`R_{\mathrm{H}_2}`$ $`=`$ $`N_\mathrm{H}F_\mathrm{H}`$
$`R_{\mathrm{HD}}`$ $`=`$ $`N_\mathrm{H}F_\mathrm{D}+N_\mathrm{D}F_\mathrm{H}`$
$`R_{\mathrm{D}_2}`$ $`=`$ $`N_\mathrm{D}F_\mathrm{D}.`$ (19)
The HD molecule can form in two channels: an H atom on the surface reacts with a newly adsorbed D atom, or a D atom on the surface reacts with a newly adsorbed H atom. Due to the longer residence time of a D atom on the surface, the second channel is dominant. Therefore, the production ratio $`R_{\mathrm{HD}}/R_{\mathrm{H}_2}`$ is approximately equal to the ratio between the population sizes $`N_\mathrm{D}/N_\mathrm{H}`$ on the surface. This amounts to an enhancement of the HD production by a factor of $`W_\mathrm{H}/W_\mathrm{D}`$ compared to the flux ratio, namely

$$\frac{R_{\mathrm{HD}}}{R_{\mathrm{H}_2}}=\frac{W_\mathrm{H}}{W_\mathrm{D}}\frac{F_\mathrm{D}}{F_\mathrm{H}}.$$
(20)
The production rate of HD is thus enhanced with respect to that of H2 by a factor given by the ratio of the desorption rates of H and D atoms. Since the desorption rates depend exponentially on the activation energy barriers for desorption $`E_1(H)`$ and $`E_1(D)`$, a small difference between these barriers can produce a large enhancement in HD production. The enhancement in D2 production is given by

$$\frac{R_{\mathrm{D}_2}}{R_{\mathrm{H}_2}}=\frac{W_\mathrm{H}}{W_\mathrm{D}}\left(\frac{F_\mathrm{D}}{F_\mathrm{H}}\right)^2.$$
(21)
giving rise to a similar enhancement compared to the expected production rate based on the ratio between the fluxes of D and H.

## 6 Simulations and results

To evaluate the production rates of H2, HD and D2 molecules we have used the master equation approach given by Eq. (4). The results presented in this paper were obtained by direct numerical integration of the set of differential equations until they reached steady state conditions. For a range of very small grains, namely $`S<5000`$, we also performed kinetic Monte Carlo simulations based on the master equation, and found perfect agreement between the two methods.

In Fig. 1 we show the average population sizes of the atomic species on the grain surface, $`N_\mathrm{H}`$ (solid line) and $`N_\mathrm{D}`$ (dashed line), as a function of the number of adsorption sites, $`S`$, and grain radius, $`r`$, under steady state conditions. The results are presented for a range of small grains on which the average number of H atoms is smaller than one while the average number of D atoms is of the order of $`10^4`$ and less. Under these conditions the master equation is required in order to evaluate the production rates correctly, although the rate equations (dotted lines) are progressively better as the size of the grain increases. Specifically, one can see that the rate equation approach appears to handle the average hydrogen abundance well for all grain sizes, while becoming much worse for the average deuterium abundance at small grain sizes. One should note that the actual production rate of molecules does not depend on $`N_\mathrm{H}`$ and $`N_\mathrm{D}`$ only \[see Eq. (16)\], so one cannot argue about the suitability of the rate equation method for the determination of production rates based on the results shown in Fig. 1.

The production rates of H2, HD and D2, molecules per grain vs. $`S`$ and $`r`$ are shown in Fig. 2. Clearly, the production rates decrease as $`S`$ is reduced. For large grains that support sizable populations of H and D, and can be described by rate equations, the production rates per grain are proportional to the surface area, namely linear in $`S`$. In the limit of diminishing grain sizes, the production rates become quadratic in $`S`$, as can be seen by substitution of $`N_\mathrm{H}`$ and $`N_\mathrm{D}`$ into Eq. (19), which leads to a quadratic dependence on flux. The range of grain sizes presented in the simulation results reaches both the limit of large grains and the limit of small grains, covering the intermediate range between these two limits.

In order to quantify the relative enhancement in the formation rate of deuterated molecules vs. H2 we present in Fig. 3 the ratios $`R_{\mathrm{HD}}/R_{\mathrm{H}_2}`$ and $`R_{\mathrm{D}_2}/R_{\mathrm{H}_2}`$ vs. $`S`$ and $`r`$. It is observed that as $`S`$ decreases these production ratios increase by over a factor of 20. In the limit of large grains, the production ratio converges towards the results of the rate equations (dotted line on right-hand side). In the limit of very small grains it converges towards the results of the asymptotic analysis presented in Sec. 5 (dotted line on left-hand side).

To examine the dependence of the efficiency of the molecular formation on the grain size we show in Fig. 4 the production rates per site of H2, HD and D2 on the grain as a function of $`S`$ and $`r`$. In the limit of large grains this quantity saturates and does not depend on the grain size. However, as the grain size decreases, the production rate per site goes down. The production rate of H2 is reduced more quickly than that of the deuterated species. This gives rise to an enhancement in the relative production rate of HD and D2 vs. H2.

To account for the production rate of molecular hydrogen in interstellar clouds one should integrate over the entire distribution of grain sizes in the cloud (Mathis et al., 1977; Mathis, 1990, 1996; Draine & Lee, 1984; Draine & Anderson, 1985). Consider a cloud in which the grain size distribution is given by $`n_{\mathrm{gr}}(r)`$ (cm-4) in the range $`r_{\mathrm{min}}<r<r_{\mathrm{max}}`$, where $`r_{\mathrm{min}}`$ and $`r_{\mathrm{max}}`$ are the lower and upper cutoffs of the distribution \[here $`n_{\mathrm{gr}}(r)dr`$ (cm-3) is the number density of grains of sizes in the range (r,r+dr)\]. The formation rates $`_{\mathrm{H}_2}`$, $`_{\mathrm{HD}}`$ and $`_{\mathrm{D}_2}`$, (cm-3 s-1 ) in the cloud are given by

$$_{\mathrm{H}_2}=_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}n_{\mathrm{gr}}(r)R_{\mathrm{H}_2}(r)𝑑r,$$
(22)
$$_{\mathrm{HD}}=_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}n_{\mathrm{gr}}(r)R_{\mathrm{HD}}(r)𝑑r,$$
(23)
and

$$_{\mathrm{D}_2}=_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}n_{\mathrm{gr}}(r)R_{\mathrm{D}_2}(r)𝑑r,$$
(24)
respectively, where $`R_{\mathrm{H}_2}(r)`$, $`R_{\mathrm{HD}}(r)`$ and $`R_{\mathrm{D}_2}(r)`$ are the production rates per grain for grains of radius $`r`$.

Observations support the assumption that the grain size distribution takes the power-law form (Mathis et al., 1977; Weingartner, 2001)

$$n_{\mathrm{gr}}(r)=\frac{K}{r^\alpha }$$
(25)
in a range of sizes $`r_{\mathrm{min}}<r<r_{\mathrm{max}}`$, where

$$K=\frac{(\alpha 1)n_{\mathrm{gr}(\mathrm{total})}}{r_{\mathrm{min}}^{1\alpha }r_{\mathrm{max}}^{1\alpha }},$$
(26)
and $`n_{\mathrm{gr}(\mathrm{total})}`$ (cm-3) is the total density of grains, integrated over the whole range of grain sizes. To exemplify the role of the grain size distribution we analyze the case in which $`\alpha =3`$. This distribution has the special property that the mass of grains in the cloud is equally distributed within the range of grain sizes. The total surface area of all grains is distributed such that grains in the size range between $`(r,r+dr)`$ have a total surface area proportional to $`1/r`$. As a result, most of the surface area belongs to grains in the lower end of the size distribution. In the analysis we chose $`r_{\mathrm{min}}=4\times 10^7`$ cm and $`r_{\mathrm{max}}=2\times 10^5`$ cm. The parameter $`K`$ is chosen such that the density of large grains, of radii in the range $`5\times 10^6<r<2\times 10^5`$ cm, is $`n_{\mathrm{gr}(\mathrm{classical})}=10^{12}n_\mathrm{H}`$. Integrating the grain size distribution over the entire range we obtain $`n_{\mathrm{gr}(\mathrm{total})}=1.66\times 10^{10}n_\mathrm{H}`$. The minimum grain size is chosen to be larger than PAH molecules so that the problem of stochastic heating and desorption by photons does not dominate normal thermal evaporation. According to the calculations of Draine & Li (2001) and Li & Draine (2001), however, grains of our minimum size may experience stochastic heating of temperatures up to 40 K at time intervals ($`10^4`$ s-1) somewhat shorter than the D evaporation time at the grain temperature chosen. Such heating would effectively shorten the grain lifetime of D atoms and reduce the isotope effect. The time scale between heating events could be lengthened by extinction in the cloud.

To examine the contribution of grains of different sizes to the total production rates of H2, HD and D2 we show in Fig. 5 the differential production rates given by the integrands of Eqs. (22) and (23) as a function of the grain radius $`r`$. Two competing effects play a role in this process. The production efficiency of molecules on the grain decreases as the grain size is decreased. However, small grains account for most of the surface area and thus provide a large fraction of the molecules that are produced. The relative enhancement in the production of deuterated species arises because their production efficiency decreases more slowly than that of H2 as the grain size is reduced.

Integrating over the range of grain sizes we obtain that the production rates are $`_{\mathrm{H}_2}=3.9\times 10^{15}`$ cm-3 s-1, $`_{\mathrm{HD}}=1.2\times 10^{18}`$ cm-3 s-1 and $`_{\mathrm{D}_2}=1.5\times 10^{23}`$ cm-3 s-1. Therefore, the ratio $`_{\mathrm{HD}}/_{\mathrm{H}_2}=3.1\times 10^4`$ ($`_{\mathrm{D}_2}/_{\mathrm{H}_2}=3.9\times 10^9`$) is larger by about a factor of $`44.1`$ ($`78.3`$) than the expected value based on the ratio of the incoming fluxes, namely $`F_\mathrm{D}/F_\mathrm{H}=10^5/\sqrt{2}`$. Observations indicate that the exponent in Eq. (25) is $`\alpha =3.5`$ (Weingartner, 2001). Therefore, the contribution of very small grains may be even more dominant, leading to a larger enhancement in the formation of HD and D2 vs. H2 than in the analysis shown here.

We will now compare our results with the standard rate law used for the production of molecular hydrogen in diffuse regions. This rate law assumes that all the grains are of radius $`r=10^5`$ cm and that their density is equal to $`n_{\mathrm{gr}(\mathrm{classical})}`$. It is also assumed that the efficiency is $`\eta =1`$, namely, all the atoms that collide with grains end up in molecules. Under these assumptions the production rate of H2 molecules is given by $`_{\mathrm{H}_2}=n_\mathrm{H}v_\mathrm{H}n_{\mathrm{gr}(\mathrm{classical})}\sigma /2`$ cm-3 s-1 (Hollenbach et al., 1971), or more conveniently, by

$$_{\mathrm{H}_2}=R(\mathrm{H}_2)n_\mathrm{H}n,$$
(27)
where $`n_\mathrm{H}`$ is the density of H atoms and $`n`$ is the total density of H nuclei in all forms (although we assume here that all the hydrogen in the gas is in the atomic form, namely $`n=n_\mathrm{H}`$). The parameter $`R(\mathrm{H}_2)`$ (cm3 s-1) is the effective rate coefficient given by

$$R(\mathrm{H}_2)=5\times 10^{17}\sqrt{T_{\mathrm{gas}}/300}.$$
(28)
For the conditions considered here, $`R(\mathrm{H}_2)=2.74\times 10^{17}`$ cm3 s-1 and $`_{\mathrm{H}_2}=6.85\times 10^{14}`$ cm-3 s-1. This result is about an order of magnitude higher than our result, which is consistent with $`R(\mathrm{H}_2)=1.56\times 10^{18}`$ cm3 s-1. The reason for this apparent discrepancy is the low production efficiency of H2 molecules under these conditions, which on large grains is $`\eta =0.09`$ (given by the rate equations) and is even lower for small grains. The production rate of HD molecules (cm-3 s-1) can be expressed in a similar way, by

$$_{\mathrm{HD}}=R(\mathrm{HD})n_\mathrm{D}n,$$
(29)
where $`R(\mathrm{HD})`$ (cm3 s-1) is the effective rate coefficient for this reaction. The values of the rate coefficients $`R(\mathrm{H}_2)`$ and $`R(\mathrm{HD})`$ under the conditions studied here, are shown in Table 1, for $`\alpha =3`$ and $`3.5`$ and for several grain temperatures between 12 and 20 K.

In Fig. 6 we present the temperature dependence of the production rates of H2 and HD molecules. The production rates $`_{\mathrm{H}_2}`$, obtained from the rate equations, are shown by the solid lines for the case that the exponent of the grain size distribution is $`\alpha =3.5`$ (upper solid line) and for $`\alpha =3`$ (lower solid line). The master equation results for these two cases are also shown (by $``$ and $`\times `$, respectively). The production rates $`_{\mathrm{HD}}`$, obtained from the rate equations, are shown by the dashed lines, for $`\alpha =3.5`$ (upper dashed line) and for $`\alpha =3`$ (lower dashed line). The master equation results are also shown (by $``$ and $`\times `$, respectively). These graphs provide much insight on the enhancement of HD production. They show that for grain temperatures lower than about 16 K, the production efficiency is $`\eta 1`$ and the rate coefficient $`R(\mathrm{H}_2)`$ is consistent with Eq. (28). At this temperature range there is no enhancement in the production of HD molecules. The enhancement appears in a range of temperatures, roughly between 16 K and 20 K, where the H2 production rate quickly decreases as $`T`$ increases while the HD production rate remains steady. Above 20 K the production rates of both H2 and HD are diminished.

In Fig. 7 we present the ratio, $`_{\mathrm{HD}}/_{\mathrm{H}_2}`$, of the production rates of HD and H2 molecules as a function of the grain temperature, as obtained from the rate equations (solid line), and the master equation for $`\alpha =3`$ ($`\times `$) and $`\alpha =3.5`$ ($``$). The enhancement ratio increases sharply as the temperature is raised above 16 K and reaches a peak around 20 K. A significant enhancement is observed on large grains, for which the rate equation results apply. When the grain size distribution is taken into account, further enhancement is obtained due to the contribution of the small grains, on which the enhancement ratio is even larger.

The dependence of the ratio of the production rates $`R_{\mathrm{HD}}/R_{\mathrm{H}_2}`$ on the grain size $`S`$ and radius $`r`$ is shown in Fig. 8 for several grain temperatures. At 14 K there is no significant enhancement for any grain size. At 16 and 18 K the enhancement strongly depends on the grain size and increases sharply as the grain size decreases. At 20 K there is a large enhancement for grains of all sizes. However, at this temperature the total production of both H2 and HD is very small, thus making little contribution to the gas phase abundances.

## 7 Discussion

In this paper we have used a simple description of the grains as spherical objects. In practice the grains are expected to be amorphous with rough surfaces and may include pores. These features may increase the surface area of the grain and introduce a distribution of diffusion and desorption barriers. However, no solid data are available to quantify such effects. Moreover, the analysis of experimental results on molecular hydrogen production on olivine and amorphous carbon surfaces (Pirronello et al., 1997a, b, 1999) did not require more than one energy barrier for each process (Katz et al., 1999).

The analysis presented in this paper is restricted to the case in which the mobility of adsorbed atoms on the surface is due to thermal hopping. Tunneling of H and D atoms on the surface is not included since experimental results indicate that tunneling does not play a significant role in the formation of molecular hydrogen under astrophysically relevant conditions (Pirronello et al., 1997a, b, 1999). In a recent analysis of the results of these experiments, chemisorption processes were taken into account (Cazaux & Tielens, 2002). The study of chemisorption processes is motivated by observations of molecular hydrogen in photon-dominated regions in which the typical grain temperature is around 40 K, making the production of H2 molecules by physisorbed H atoms hopelessly inefficient.

The mechanism presented here for an enhanced production of HD and D2 molecules, vs. H2, takes place in a window of grain temperatures. For the surface parameters used here this temperature window is between 16 - 20 K. In this temperature range the production efficiency of H2 is reduced, while the efficiency of HD and D2 production remains steady. The effect is observed for grains of all sizes, but the relative enhancement is larger on small grains. It will be interesting to see if the enhancement plays a significant role in chemical models of diffuse and translucent clouds; relevant model calculations are planned by Le Petit and Roueff (private communication).

## 8 Summary

The formation of H2, HD and D2 molecules on dust grain surfaces was studied using the master equation approach. The production of these molecules was assumed to occur via the Langmuir-Hinshelwood (diffusive) mechanism. It was shown that for surface temperatures in excess of 15 K, the fractions of HD and D2 molecules in the population of hydrogen molecules that are produced are significantly enhanced compared with their expected ratio based on the gas phase abundances. The enhancement takes place on grains of all sizes, but is largest on the smallest grains considered until the temperature reaches 20 K. It arises because of the assumption that D atoms are bound slightly more strongly to grains than are H atoms. At a surface temperature of 18 K, with binding energies and rates of diffusion based on laboratory experiments, the enhancement is found to be a factor of 44.1 (78.3) for the production of HD (D2)with the assumption that the grain size distribution goes as $`r^3`$. Although the enhancement is even larger at higher temperatures, the absolute rates of production for H2 and its isotopomers become too small to be of importance. Use of the master equation method can lead to somewhat greater enhancements (factor of a few) than the simple rate equation method.

While gas-phase mechanisms (Watson, 1973) may account for the observed abundance of HD, it is possible that surface mechanisms such as the one presented here also play a role in the production of deuterated molecules. Recent observations have uncovered a variety of multiply deuterated molecular species in the interstellar medium such as D2CO and ND3 (Ceccarelli et al., 2002; van der Tak et al., 2002). Gas phase reaction networks that lead to the formation of such multiply deuterated species have been studied and their efficiency examined (Charnley et al., 1997; Rodgers & Charnley, 2001; Roberts et al., 2002, 2003). The picture in pre-stellar cores is that gas-phase fractionation can explain much of what is observed, but the picture in protostellar cores, where grain mantles tend to be evaporated back into the gas, is that surface reactions play an important role in the fractionation leading to multi-deuterated species (Stantcheva et al., 2002). It will be interesting to see if the surface mechanism presented here may be extended to produce multiply deuterated species. Preliminary work on this matter has been undertaken by Caselli et al. (2002) and Stantcheva & Herbst (2003).

## 9 acknowledgments

This work was supported by the Adler Foundation for Space Research of the Israel Science Foundation. E. H. wishes to thank the National Science Foundation (U. S.) for support of his research programme in astrochemistry.

## Appendix A Isotope effects in related experiments

Isotope effects resembling the type studied here were observed experimentally on various surfaces including metals (Hoogers et al., 1995) and silicon (Koehler et al., 1988). However, these experiments were done at much higher temperatures than those relevant for molecular hydrogen formation. Also, the surfaces on which isotope effects were studied are more ordered than the astrophysically relevant surfaces. Therefore, it may not be possible to draw conclusions from these experiments for astrophysically relevant surface processes.

A crucial feature in the recent experiments on molecular hydrogen formation is the use of two beams, of H and D atoms, and the measurement of the HD production rate (Pirronello et al., 1997a, b, 1999). The use of HD molecules improves the signal-to-noise ratio in the experiment and allows reliable measurement of the production rate. In the analysis of the experimental results it was assumed, for simplicity, that there is no isotope effect, namely, that the diffusion and desorption barriers for H and D atoms are identical (Katz et al., 1999). This assumption was made in order to limit the number of free parameters used in the fitting. Under the experimental conditions this assumption had little effect on the results of the analysis. To explain this issue, let us recall the experimental procedure. The experiment followed the standard temperature programmed desorption (TPD) procedure, which consists of two stages. In the first stage, the surface was irradiated by two beams: one of H atoms and the other of D atoms. The beams were of (approximately) equal fluxes and lasted for an equal irradiation time, giving rise to equal coverages of H and D atoms. In the second stage the surface temperature was raised at a constant rate and the production rate of HD molecules was recorded using a mass spectrometer.

The analysis of the experimental results was based on the fact that the surface was of macroscopic size, so that rate equations of the form (5) can account correctly for the atomic populations on the surface and for the production rates of molecules. In the first stage of the experiment the fluxes were $`F_\mathrm{H}=F_\mathrm{D}`$, while the diffusion and desorption rates were practically zero, due to the low surface temperature. In the second stage the two beams were turned off, while the diffusion and desorption coefficients quickly increased due to the heating of the surface.

The isotope effect studied here gives rise to a small difference in the energy barriers for diffusion of H and D atoms, $`E_0(\mathrm{H})<E_0(\mathrm{D})`$, and a larger difference in the desorption energies, where $`E_1(\mathrm{H})<E_1(\mathrm{D})`$. We will now examine how these differences would affect the results of the TPD experiment. We first consider the effect of the difference in the diffusion barriers. The rate of formation of HD molecules is described by Eq. (9), where the diffusion rates of the H and D atoms sum up. It is thus clear that even if the D atoms are immobile, this gives rise only to about a factor of 2 in the production rate compared to the case in which $`A_\mathrm{D}=A_\mathrm{H}`$. Such differences have little effect on the energy barrier for diffusion, $`E_0`$, obtained from fitting the experimental TPD curve. This is because $`E_0`$ appears in the exponent of Eq. (4). Therefore, when the experimental TPD curves are fitted by a rate equation model that includes a single energy barrier for diffusion, $`E_0`$, the result is very close to the lower diffusion barrier among the H and D, namely to $`E_0(\mathrm{H})`$.

Consider the case in which the desorption barrier of D atoms, $`E_1(\mathrm{D})`$, is larger than the desorption barrier of H atoms, $`E_1(\mathrm{H})`$. In this case, as the temperature is raised during a TPD experiment, D atoms desorb more slowly than H atoms. The formation rate of HD molecules is thus limited by the availability of H atoms. When the resulting TPD curves are fitted by a rate equation model that includes only a single energy barrier for desorption, $`E_1`$, the fitted value turns out to be a very good approximation to $`E_1(\mathrm{H})`$.

It thus turns out that slower diffusion and desorption rates of D atoms vs. H atoms, have little effect on the TPD curves that describe the rate of formation of HD molecules. When these TPD curves are fitted by a rate equation model with a single diffusion barrier, $`E_0`$, and a single desorption barrier, $`E_1`$, these barriers provide good approximations to $`E_0(\mathrm{H})`$ and $`E_1(\mathrm{H})`$. This conclusion is not limited to the case in which the experiment is done on a macroscopic surface. It is expected to be valid even if the experiment is done on a powder that consists of small grains.

The origin of the isotope effect on dust grains in the interstellar medium is the long residence time of D atoms on the grain compared with H atoms. As a result, an (H or D) atom that hits a grain is more likely to encounter a D atoms on the grain than expected based on the gas phase abundance ratio. In the TPD experiments all the H and D atoms are deposited on the surface at the initial stage of the experiment, when the surface temperature is low and the adsorbed atoms are immobile. During the heating stage, when the adsorbed atoms become mobile, the beams are turned off and there is no fresh supply of H or D atoms. Therefore, the D atoms that remain on the surface after most of the H atoms are desorbed tend to form D2 molecules. The isotope effect is thus expected to affect the TPD curve of D2 molecules. However, it does not affect the TPD curve of HD molecules, which is the one that was measured in the experiments.


text:
# Dissociation and Decay of Ultra-cold Sodium Molecules

## Abstract

The dissociation of ultracold molecules is studied by ramping an external magnetic field through a Feshbach resonance. The observed dissociation energy shows non-linear dependence on the ramp speed and directly yields the strength of the atom-molecule coupling. In addition, inelastic molecule-molecule and molecule-atom collisions are characterized.

Recently, it has become possible to create ultracold molecular gases from precooled atomic samples Wynar ; Donley ; Regal ; Jochim ; Cubizolles ; Strecker ; Durr ; Chin ; Xu ; Greiner ; Jochim2 . Extending the ultralow temperature regime from atoms to molecules is an important step towards controlling the motion of more complicated objects. The complex structure of molecules may lead to new scientific opportunities, including the search for a permanent electric dipole moment, with sensitivity much higher than for heavy atoms Hudson , and the realization of quantum fluids of bosons and fermions with anisotropic interactions Baranov . Furthermore, stable mixtures of atomic and molecular condensates are predicted to show coherent stimulation of atom-molecule or molecule-atom conversion, constituting quantum-coherent chemistry Heinzen .

To date, all realizations of ultracold molecules have bypassed the need for direct cooling of the molecules, which is difficult due to the complicated ro-vibrational structure. Rather, molecules were formed from ultracold atoms using Feshbach resonances Donley ; Regal ; Jochim ; Cubizolles ; Strecker ; Durr ; Chin ; Xu ; Greiner ; Jochim2 , where a highly-vibrational excited molecular state is magnetically tuned into resonance with a pair of colliding atoms.

In this paper, we study the dissociation and decay of such highly excited molecules. Figure 1 shows the relevant energy levels. For magnetic fields above the Feshbach resonance, the molecular state crosses the free atomic states, shown here as discrete states in a finite quantization volume. The interaction between atoms and molecules turns these crossing into anti-crossings. When the magnetic field is swept very slowly through the resonance, the molecules will follow the adiabatic curve and end up in the lowest energy state of the atoms. For faster ramps, the molecular populations will partially cross some of the low-lying states, and the dissociation products will populate several atomic states. The stronger the coupling between the molecular state and the atomic states, the faster the molecules dissociate and the smaller the energy release in the dissociation. Observing the atom-molecule coupling in one-body decay (dissociation) is a new method to experimentally determine the strength of a Feshbach resonance. Previous measurements used two- or three-body processes to characterize the Feshbach resonance and therefore required accurate knowledge of the atomic density distribution.

Collisional properties of the molecules were also studied. Inelastic collisions limit both the production of molecules and their lifetime. We observed loss of molecules by collisions both with atoms and other molecules. These two processes were studied separately because we could produce atom-molecule mixtures, as well as pure molecular samples, by separating atoms and molecules with short pulses of laser light Xu .

To generate molecules, sodium condensates in the $`|F,m_F=|1,1`$ state were prepared in an optical dipole trap. The trap frequencies of 290 Hz in the radial direction and 2.2 Hz in the axial direction yielded typical densities of $`1.7\times 10^{14}`$ cm-3 for 5 million atoms. Atoms were then spin-flipped using an adiabatic radio frequency sweep to the $`|1,1`$ state, where a 1 G wide Feshbach resonance exists at 907 G Inouye .

The magnetic field sequence used to create and detect Na2 molecules was identical to our previous work Xu . Briefly, the axial magnetic field was ramped to 903 G in 100 ms using two pairs of bias coils. In order to prepare the condensate on the negative scattering length side of the resonance, the field was stepped up to 913 G as quickly as possible ($``$ 1 $`\mu `$s) to jump through the resonance with minimal atom loss. The field was then ramped back down to 903 G in 50 $`\mu `$s to form molecules. In order to remove non-paired atoms from the trap, the sample was irradiated with a 20 $`\mu `$s pulse of resonant light. Because 903 G is far from the Feshbach resonance, the mixing between atomic and molecular states was small, and therefore molecules were transparent to this “blast” pulse. By ramping the field back to 913 G molecules were converted back to atoms. Absorption images were taken at high fields (either at 903 G or 913 G) after 10 to 17 ms ballistic expansion, with the imaging light incident along the axial direction of the condensate.

To study the momentum distribution of the back-converted atoms, the magnetic field was ramped up immediately after turning off the optical trap, or for a reference, at the end of the ballistic expansion. The difference between the energies of ballistic expansion is the released dissociation energy. Energies were obtained from the rms width of the cloud $`x^2`$ as $`E=3mx^2/2t^2`$, where $`t`$ is the ballistic expansion time, and $`m`$ is the atomic mass. Figure 2 shows that faster field ramps created hotter atoms.

An earlier theoretical treatment assumes a constant predissociation lifetime of the molecules and predicts a linear relation between dissociation energy and field ramp rate Abeelen . This theory predicts a much faster dissociation (and therefore smaller dissociation energy) than was observed. Furthermore, our data shows a non-linear dependence. Linear behavior would be expected if the lifetime of the molecules was independent of the energy $`ϵ`$ from the dissociation threshold. The fact that the slope becomes smaller for increasing ramp rate indicates that the lifetime of molecules decreases with the ramp rate. As we will show, this can be explained by an increase of the density of atomic states, leading to a $`\sqrt{ϵ}`$ dependence of the molecular decay rate (Wigner threshold law Wigner ).

The decay rate $`\mathrm{\Gamma }(ϵ)`$ follows from Fermi’s golden rule as $`\mathrm{}\mathrm{\Gamma }(ϵ)=2\pi |V_{ma}(ϵ)|^2D(ϵ)`$ Mies , where $`V_{ma}`$ is the matrix element between atomic and molecular states, which to leading order is independent of $`ϵ`$. The density of states $`D(ϵ)`$ is given by

$$D(ϵ)=\frac{V}{(2\pi )^2}\left(\frac{m}{\mathrm{}^2}\right)^{3/2}ϵ^{1/2},$$
(1)
where $`V`$ is the quantization volume for free atomic states.

An expression for the matrix element $`V_{ma}`$ is obtained by comparing the energy shift near a Feshbach resonance with second-order perturbation theory. Assuming two atoms in a volume $`V`$, the energy shift of the low-lying continuum states due to the coupling with a bound molecular state is

$$\delta (ϵ)=\frac{|V_{ma}|^2}{ϵ}=\frac{|V_{ma}|^2}{\mathrm{\Delta }\mu (BB_0)},$$
(2)
where $`\mathrm{\Delta }\mu `$ is the difference between atomic and molecular magnetic moments, $`B`$ the applied magnetic field, and $`B_0`$ the position of the Feshbach resonance.

The energy shift can also be expressed in terms of the mean field energy $`4\pi \mathrm{}^2a/mV`$, where $`a=a_{bg}\mathrm{\Delta }B/(BB_0)`$ is the scattering length near the Feshbach resonance ($`a_{bg}`$ is the background scattering length and $`\mathrm{\Delta }B`$ is the resonance width Timmermans ):

$$\delta (ϵ)=\frac{4\pi \mathrm{}^2}{mV}\frac{a_{bg}\mathrm{\Delta }B}{BB_0}.$$
(3)
Comparing eq. (2) and eq. (3) yields

$$|V_{ma}|^2=\frac{4\pi \mathrm{}^2}{mV}a_{bg}\mathrm{\Delta }\mu \mathrm{\Delta }B.$$
(4)
If the entire population is initially in the molecular state, the fraction of molecules, $`m(ϵ)`$, at energy $`ϵ`$ follows the rate equation,

$`n{\displaystyle \frac{dm(ϵ)}{dϵ}}`$ $`=`$ $`{\displaystyle \frac{dm(ϵ)}{dt}}\left({\displaystyle \frac{dϵ}{dt}}\right)^1=\mathrm{\Gamma }(ϵ)m(ϵ)\left({\displaystyle \frac{dϵ}{dt}}\right)^1`$ (5)
$`=`$ $`{\displaystyle \frac{2\pi |V_{ma}(ϵ)|^2D(ϵ)}{\mathrm{}\mathrm{\Delta }\mu |\dot{B}|}}m(ϵ).`$ (6)
Using Eqs. (1) and (4), we solve the differential equation for $`m(ϵ)`$

$`m(ϵ)`$ $`=`$ $`e^{\frac{2}{3}Cϵ^{3/2}},`$ (7)
$`C`$ $`=`$ $`{\displaystyle \frac{2\mathrm{\Delta }B}{\mathrm{}\dot{B}}}\sqrt{{\displaystyle \frac{ma_{bg}^2}{\mathrm{}^2}}}.`$
In the lab frame, the atoms have kinetic energy $`ϵ/2`$ and therefore the average energy of an atom after dissociation is

$$\delta E=_0^{\mathrm{}}\frac{ϵ}{2}\left(dm(ϵ)\right)=0.591\left(\sqrt{\frac{\mathrm{}^2}{ma_{bg}^2}}\frac{\mathrm{}\dot{B}}{2\mathrm{\Delta }B}\right)^{2/3}.$$
(8)
Using theoretical values $`a_{bg}=3.4`$ nm, $`\mathrm{\Delta }\mu /h=3.65\times 1.4`$ MHz/G, and $`\mathrm{\Delta }B=0.98`$ G Yurovsky43611 ; Yurovsky43605 , our parameter-free prediction (solid line in Fig. 2) shows good agreement with the experimental data. Alternatively, we may regard the width of the Feshbach resonance as a free parameter to obtain a fitted value of $`\mathrm{\Delta }B=1.07\pm 0.02`$ G (dotted line in Fig. 2). Compared to previous mean-field measurements Inouye ; Stenger , our determination of the resonance width is more accurate and free from systematic errors associated with the determination of atomic densities.

Further experiments with ultracold sodium molecules will critically depend on their collision properties. Therefore we also studied heating and inelastic collision processes. As shown in Fig. 3, we observed monotonic heating of the pure molecular sample over $``$ 30 ms. In addition, we observed short timescale oscillations (6 ms) in the fitted temperature (inset of Fig. 3). Such breathing oscillations were excited because the molecules were formed over the volume of the atomic condensate, which was larger than the equilibrium volume for the molecules. The absence of damping implies a collision time of at least 6 ms, or a molecular scattering length smaller than 17 nm (obtained using the expression for the collision rate $`8\pi a^2v_{th}n_m`$ where $`v_{th}`$ is the thermal velocity). It is unclear whether the oscillation disappeared due to collisions or limited signal-to-noise ratio.

The temperature of the molecular cloud saturated at $``$ 250 nK after 15 ms. A possible explanation is the balance between heating due to inelastic molecular decay and the evaporative cooling caused by the finite trap depth (1.7 $`\mu `$K). This would imply a collision time of 15 ms. However, we have no clear evidence that thermalization has occurred. Clearly, further studies of elastic collisions between ultracold molecules are necessary.

Molecules formed via Feshbach resonances are created in high vibrational states. Therefore, one expects vibrational relaxation to be a strong, inelastic decay mechanism. Vibrational energy spacings are much larger than the trap depth, leading to loss of molecules from the trap.

Figure 4(a) shows the decay of a pure molecular sample. The decay was analyzed with the rate equation

$$\frac{\dot{N}_m}{N_m}=K_{mm}n_m.$$
(9)
Here $`n_m`$ is the density of the molecules, and $`K_{mm}`$ is the molecule-molecule collision rate coefficient. Because of the changing size and temperature of the molecular cloud during the first $``$ 15 ms (Fig. 3), we only fit data points at later times, assuming a thermal equilibrium volume for the molecules. The decay at earlier times is slower, consistent with a larger molecular cloud. The fit yields a molecule-molecule collision coefficient of $`K_{mm}5.1\times 10^{11}`$ cm3/s, about 2 orders of magnitude larger than the typical values reported for fermions Jochim ; Cubizolles .

Inelastic collisions between molecules and atoms were also observed by keeping atoms in the trap (Fig. 4(b)). The decay was analyzed assuming that the loss of molecules occurred mainly due to collisions with atoms, resulting in an exponential decay:

$$\frac{\dot{N}_m}{N_m}=K_{am}n_a.$$
(10)
Here $`N_m`$ is the number of the molecules, $`n_a`$ is the density of atoms, and $`K_{am}`$ is the atom-molecule collision rate coefficient. From the fit, we extract a lifetime of 106 $`\mu `$s and a rate coefficient $`K_{am}5.5\times 10^{11}`$ cm3/s, which agrees well with theoretical predictions Yurovsky43605 ; Yurovsky43611 .

The inelastic losses determine the maximum conversion efficiency from atoms to molecules. For an adiabatic ramp, one expects close to 100% conversion efficiency. Indeed, in experiments with fermionic atoms, efficiencies up to 85% have been observed Cubizolles . Figure 5 shows the results for magnetic field ramps of different durations. The two sets of images show that applying the blast pulse dramatically improved the molecular number and temperature. Without it, a slower ramp time (4 ms) appeared to be more favorable for molecule formation (open circles in Fig. 5(c)). No molecules were observed for a 50 $`\mu `$s ramp time. However, with the blast pulse, nearly the same number of molecules was obtained for all ramp times between 50 $`\mu `$s to 4 ms (closed circles in Fig. 5(c)).

We interpret our data as the interplay of two competing processes. The adiabatic condition requires a relatively slow field ramp for efficient conversion. However, this means that the atoms and molecules spend more time near or at the Feshbach resonance, where inelastic collision rates are enhanced. In contrast to Fig.5(b), the absence of molecular signal in Fig. 5(a) for 50 $`\mu `$s ramp time reflects that the atomic density reduction due to the mean-field expansion is too slow for the molecules to survive the inelastic collisions with the atoms.

In conclusion, we observed a Wigner threshold behavior in the dissociation of ultracold molecules. We were able to characterize a Feshbach resonance using a one-body decay (dissociation) process. The rapid decay of the molecules due to collisions with atoms and other molecules imposes a severe limit to further evaporative cooling for bosons. This also explains the low conversion efficiency ($``$ 4%), in contrast to recent experiments with fermions.

The authors would like to thank D. E. Miller for experimental assistance and A. E. Leanhardt, M. Saba and D. Schneble for their critical reading of the manuscript. We also thank P. S. Julienne, B. J. Verhaar, and V. A. Yurovsky for useful discussion. This research is supported by NSF, ONR, ARO and NASA.